算法性能优化

Question

我写了以下代码。Stack溢出的一个用户建议我在这里编写代码，以便进行优化检查。

我怎么才能优化它？

我能用一个专业的，更好的形式，提高性能的代码吗？

说明:此代码计算c()函数中三个矩阵的乘法。然后，将Omo作为旧的，Omn作为新的引入c()计算L()。最后，定义两条if语句，接受大于1的Ls和随机数的R3。此代码为MCMC拟合数据。这里我适合omn和Mn。

import numpy as np
import emcee
import matplotlib.pyplot as plt
from math import *
from scipy.integrate import quad
from scipy.integrate import odeint
O_m=chi=None
xx=np.array([0.01,0.012])  
yy=np.array([32.95388698,33.87900347])
Cov=[[137,168],[28155,-2217]]     # this is a covariance matrix for simplification I chose the 2*2 one

temp=1e5
z0=0
M=2
Omo = 0.3
Odo=0.7
H0=70

def ant(z, Om, Od):
    return 1/sqrt(((1+z)**2)*(1+Om*z)-z*(2+z)*Od)
def dl(n, Om, Od, M):
    Od=1-Om
    q=quad(ant,0,xx[n],args=(Om,Od))[0]
    h=5*log10((1+xx[n])*q)
    fn=(yy[n]-M-h)                  
    return fn

def c(Om, Od, M):
    f_list = []
    for i in range(2):  # the value '2' reflects matrix size
        f_list.append(dl(i,Om,Od,M))
    rdag=[f_list]
    rmat=[[f] for f in f_list]
    a=np.dot(rdag,Cov)
    b=np.dot(a,rmat)
    Matrix=np.linalg.det(b)*0.000001
    return Matrix
N=2000
with open('txtfile.txt', 'w') as f:
    for i in range (1,N):
        R1=np.random.uniform(0,1)
        R2=np.random.uniform(0,1)
        R3=np.random.uniform(0,1)
        R4=np.random.uniform(0,1)
        def delta(r1, r2):
            sig=0.04
            d=sig*(np.sqrt(-2*np.log(r1))*np.cos(np.radians(r2)))
            return d
        Omn=Omo+delta(R1, R2)
        Odn=1-Omn
        Mn=M+delta(R3,R4)
        R3=np.random.uniform(0,1)

        def L():
            l=np.exp(-0.5*(c(Omn,Odn,Mn)-c(Omo,Odo,M)))
            return l
        if L()>1:
            O_m=Omn
            chi=c(Omn,Odn,Mn)

        elif L()>R3:
            O_m=Omn
            chi=c(Omn, Odn, Mn)


        f.write("{0}\t{1}\n".format(chi, O_m))
    print("Minimum of chi squre is")
    if chi<temp:
        temp=chi
        chimin=temp
    print(chimin)
    print(input("Press any key to exit... "))

我感谢你的帮助和关注

Answer 1

导入*

因为您只需要sqrt和来自math的log10，所以我会更改如下：

from math import *

至：

from math import sqrt, log10

PEP8

• 使用描述性方法名称。我不知道dl或c是什么意思
• Valiable (和参数)名称应该是描述性的。很难确定Omo是什么。
• 对参数和变量名使用lower_case
• 算子间距

未使用参数

def dl(n, Om, Od, M):
    Od = 1 - Om

如果随后立即覆盖它，那么为什么要传入Od呢？

实际上，Od只在ant中使用，那么为什么不直接在那里计算它，而忽略其余的呢？

全局变量

dl在参数列表中不使用xx。它也只使用n索引xx，为什么不立即传递来自xx的相关值作为参数呢？

使用此全局变量会使此方法更难用于其他值。

内函数

函数delta和L被定义为for-循环的所有2000次迭代。这是不必要的

L使用全局状态，并多次计算，为什么不直接这么做呢？

l = np.exp(-0.5 * (c(Omn, Odn, Mn) - c(Omo, Odo, M)))

把它变成变量而不是函数？

分离函数

即使在一个简单的脚本中，您也可以将结果的生成与输入的表示分开。

使用if __name__ == '__main'__:

因此，您的程序可以运行脚本，但也可以作为模块导入。

R3

为什么要生成一个同名的新随机数？要么是相同的数字，所以不需要生成它，要么是不同的数字，所以给它一个新的名称

此外，如果L()大于0到1之间的一个数字，它将大于1，因此if .. elif是不必要的。

幻数

c中的0.000001和delta中的0.04从何而来？如果这些常量，最好是这样命名。

避免重复

程序的设置方式是，Cov经常转换为numpy数组，而c(Omn、Odn、Mn)和c(Omo, Odo, M)是一次又一次的计算，参数没有任何变化，最好只做一次。

f_list

而不是附加，我会使用一个列表理解。

更进一步，既然您仍然需要它作为np.array，为什么不使用np.fromiter和生成器表达式呢？为什么要从二维数组中提取行列式，而不是对一维数组进行计算呢？

总之，我最终会有这样的结果：

import numpy as np
import emcee
import matplotlib.pyplot as plt
from math import sqrt, log10
from scipy.integrate import quad
from scipy.integrate import odeint
from io import StringIO


def ant(z, Om):
    Od = 1 - Om
    return 1 / sqrt(((1 + z) ** 2) * (1 + Om * z) - z * (2 + z) * Od)


def dl(x, y, Om, M):
    q = quad(ant, 0, x, args=(Om,))[0]
    h = 5 * log10((1 + x) * q)
    return y - M - h


def c(xx, yy, cov, Om, M):
    f_list = np.fromiter(
    (dl(x, y, Om, M) for x, y in zip(xx, yy)), 
    dtype=float,
    count=len(xx),
    )
    a = np.dot(f_list, cov)
    b = np.dot(a, f_list.T)
    return b * 0.000001

def delta(r1, r2):
    sig = 0.04
    return sig * (np.sqrt(-2 * np.log(r1)) * np.cos(np.radians(r2)))

    
def calculation(xx, yy, cov, M, Omo, N):

    c0 = c(xx, yy, cov, Omo, M)
    for i in range(1, N):
        R1 = np.random.uniform(0, 1)
        R2 = np.random.uniform(0, 1)
        R3 = np.random.uniform(0, 1)
        R4 = np.random.uniform(0, 1)
        R5 = np.random.uniform(0, 1)

        Omn = Omo + delta(R1, R2)
        Mn = M + delta(R3, R4)

        c_ = c(xx, yy, cov, Omn, Mn)
        l = np.exp(-0.5 * (c_ - c0))

        if l > R5:
            O_m = Omn
            chi = c_

        yield chi, O_m
        
        
def main(filehandle):
    temp = 1e5
    z0 = 0
    M = 2
    Omo = 0.3
    H0 = 70

    N = 2000
    xx = np.array([0.01,0.012])  
    yy = np.array([32.95388698,33.87900347])
    cov = np.array([[137,168],[28155,-2217]])
    
    
    
    for chi, O_m in calculation(xx, yy, cov, M, Omo, N):
        filehandle.write("{0}\t{1}\n".format(chi, O_m))
        
    if chi < temp:
        temp = chi
        chimin = temp
    return chimin
    
if __name__ == '__main__':
    # with StringIO() as filehandle:
    with open('txtfile.txt', 'w') as filehandle:
        chimin = main(filehandle)
        # print(filehandle.getvalue())
        print("Minimum of chi squre is")

        print(chimin)
        print(input("Press any key to exit... "))

我不知道这是否会显着地加速您的算法，但它会有帮助，特别是更简单的c。对于这个简单的输入，它在我的机器上所需的时间减少了33%。

更大的收益将通过矢量化得到更多，但为此，我需要对正在发生的事情有更好的数学理解