时间戳

首页 > 解决方案 > Basin Hopping方法的编程和优化代码

python - Basin Hopping方法的编程和优化代码

问题描述

所以我正在编写一个从零开始的 Basin Hoping 以找到相互作用粒子系统的潜在最小值，到目前为止我已经获得了很好的结果，但是由于我增加了系统的粒子数量，所以代码需要更多时间来运行. 我正在使用 scipy 共轭梯度来查找局部最小值。我没有收到任何错误消息，程序似乎运行良好，但我想知道如何优化代码以减少计算时间。

我将盆地跳跃定义为一个函数，给定：

def bhmet(n,itmax, pot,derpot, T):
i = 1 
phi = np.random.rand(n)*2*np.pi
theta = np.random.rand(n)*np.pi
x = 3*np.sin(theta)*np.cos(phi)
y = 3*np.sin(theta)*np.sin(phi)
z = 3*np.cos(theta)
xyzat = np.hstack((x,y,z))
vmintot = 0 
while i <= itmax: 
    print(i)   
    plmin = optimize.fmin_cg(pot,xyzat,derpot,gtol = 1e-5,) #posiciones para el mínimo local.4
    if pot(plmin) < vmintot or np.exp((1/T)*(vmintot-pot(plmin))) > np.random.rand():
        vmintot = pot(plmin)
    xyzat = plmin + 2*0.3*(np.random.rand(len(plmin))-0.5)
    i = i + 1  
return plmin,vmintot

我尝试将初始条件（第一个'xyzat'）定义为矩阵，但 scipy.optimize.fmin_cg 的参数被请求为数组（这就是为什么在函数中我将数组重塑为矩阵）。

我正在搜索全局最小值的函数是：

def ljpot(posiciones,):
r = np.array([])
matpos = np.zeros((int((1/3)*len(posiciones)),3))
matpos[:,0] = posiciones[0:int((1/3)*len(posiciones))]
matpos[:,1] = posiciones[int((1/3)*len(posiciones)):int((2/3)*len(posiciones))]
matpos[:,2] = posiciones[int((2/3)*len(posiciones)):]    
for j in range(0,np.shape(matpos)[0]):
    for k in range(j+1,np.shape(matpos)[0]):
        ri = np.sqrt(sum((matpos[k,:]-matpos[j,:])**2))
        r = np.append(r,ri) 
V = 4*((1/r)**12-(1/r)**6)
vt = sum(V)
return vt

它的梯度是：

def gradpot(posiciones,):
gradv = np.array([])
matposg = np.zeros((int((1/3)*len(posiciones)),3))
matposg[:,0] = posiciones[:int((1/3)*len(posiciones))]
matposg[:,1] = posiciones[int((1/3)*len(posiciones)):int((2/3)*len(posiciones))]
matposg[:,2] = posiciones[int((2/3)*len(posiciones)):]   
for w in range(0,np.shape(matposg)[1]): #índice que cambia de columna. 
    for k in range(0,np.shape(matposg)[0]): #índice que cambia de fila.
        rkj = np.array([])
        xkj = np.array([])
        for j in range(0,np.shape(matposg)[0]): #también este cambia de fila. 
            if j != k:
                r = np.sqrt(sum((matposg[j,:]-matposg[k,:])**2))
                rkj = np.append(rkj,r) 
                xkj = np.append(xkj,matposg[j,w])
        dEdxj = sum(4*(6*(1/rkj)**8- 12*(1/rkj)**14)*(matposg[k,w]-xkj))
        gradv = np.append(gradv,dEdxj)
return gradv

我将数组输入转换为矩阵的原因是，对于每个粒子，都有三个坐标 x、y、z，因此矩阵的列是每个粒子的 x、y、z。我尝试使用 np.reshape() 来执行此操作，但对于程序已经获得正确结果的系统，它似乎给了我错误的结果。

代码似乎运行良好，但随着粒子数量的增加，运行时间呈指数增长。我知道全局优化可能需要很长时间，但也许我对代码有点混乱。我不知道减少运行时间的方法是否很明显，我对优化代码有点陌生，如果是这样，对不起。当然，欢迎任何建议。非常感谢大家！

标签： pythonoptimizationspyder

解决方案

快速浏览后我注意到两件事，您绝对可以节省一些时间。两者之前都需要更多思考，但之后您将获得优化且更清晰的代码。

1.尽量避免使用append。 append是非常低效的。你从一个空数组开始，之后每次都需要分配更多的内存。这会导致内存处理效率低下，因为每次附加数字时都会复制数组。数组越长，效率就越低append。

替代方案：用于np.zeros((m,n))初始化数组，使用m和n是它最终的大小。然后您需要一个将新值放在相应位置的计数器。如果在计算之前没有定义数组的大小，可以将其初始化为一个大数，然后将其剪切。

2.尽量避免使用for循环。它们通常非常慢，尤其是在遍历大型矩阵时，因为您需要单独索引每个条目。

替代方案：矩阵运算通常要快得多。例如，您可以先定义两个矩阵，而不是在两个循环中，它们r = np.sqrt(sum((matposg[j,:]-matposg[k,:])**2))对应于and （应该可以不使用循环！）然后简单地使用forABmatposg[j,:]matposg[k,:]r = np.linalg.norm(A-B)

推荐阅读