为 Basin Hopping 方法编程和优化代码
Programing and optimizing code for a Basin Hopping method
所以我正在编写一个从零开始的 Basin Hoping 程序来找到相互作用粒子系统的潜在最小值,到目前为止我已经获得了很好的结果,但是自从我增加了系统的粒子数量后,代码运行 需要更多时间。我正在使用 scipy 共轭梯度来寻找局部最小值。我没有收到任何错误消息,程序似乎运行良好,但我想知道如何优化代码以减少计算时间。
我将 Basin Hopping 定义为一个函数,给定:
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() 来做到这一点,但对于程序已经获得正确结果的系统,它似乎给了我错误的结果。
代码似乎运行良好,但随着我增加粒子数量,运行ning 时间呈指数增长。我知道全局优化可能需要很长时间,但也许我在处理代码时有点乱。我不知道减少 运行ning 时间的方法是否很明显,我对优化代码有点陌生,如果是这样的话,我很抱歉。当然,欢迎任何建议。非常感谢大家!
快速浏览后我注意到两件事,您绝对可以节省一些时间。之前两者都需要更多思考,但之后您将获得优化和更清晰的代码。
1.尽量避免使用 append
。 append
效率非常低。您从一个空数组开始,之后每次都需要分配更多内存。这会导致内存处理效率低下,因为每次附加数字时都会复制数组。数组越长,append
变得越低效。
备选方案: 使用 np.zeros((m,n))
初始化数组,m
和 n
是数组最终的大小。然后您需要一个计数器,将新值放在相应的位置。如果你在计算前没有定义数组的大小,你可以将它初始化为一个大数,然后切割它。
2。尽量避免使用 for
循环。 它们通常很慢,尤其是在遍历大型矩阵时,因为您需要单独索引每个条目。
备选方案: 矩阵运算通常要快得多。例如,不是在两个 for
循环中使用 r = np.sqrt(sum((matposg[j,:]-matposg[k,:])**2))
,而是可以先定义两个矩阵 A
和 B
,它们对应于 matposg[j,:]
和 matposg[k,:]
(应该可以不使用循环!)然后简单地使用 r = np.linalg.norm(A-B)
所以我正在编写一个从零开始的 Basin Hoping 程序来找到相互作用粒子系统的潜在最小值,到目前为止我已经获得了很好的结果,但是自从我增加了系统的粒子数量后,代码运行 需要更多时间。我正在使用 scipy 共轭梯度来寻找局部最小值。我没有收到任何错误消息,程序似乎运行良好,但我想知道如何优化代码以减少计算时间。
我将 Basin Hopping 定义为一个函数,给定:
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() 来做到这一点,但对于程序已经获得正确结果的系统,它似乎给了我错误的结果。
代码似乎运行良好,但随着我增加粒子数量,运行ning 时间呈指数增长。我知道全局优化可能需要很长时间,但也许我在处理代码时有点乱。我不知道减少 运行ning 时间的方法是否很明显,我对优化代码有点陌生,如果是这样的话,我很抱歉。当然,欢迎任何建议。非常感谢大家!
快速浏览后我注意到两件事,您绝对可以节省一些时间。之前两者都需要更多思考,但之后您将获得优化和更清晰的代码。
1.尽量避免使用 append
。 append
效率非常低。您从一个空数组开始,之后每次都需要分配更多内存。这会导致内存处理效率低下,因为每次附加数字时都会复制数组。数组越长,append
变得越低效。
备选方案: 使用 np.zeros((m,n))
初始化数组,m
和 n
是数组最终的大小。然后您需要一个计数器,将新值放在相应的位置。如果你在计算前没有定义数组的大小,你可以将它初始化为一个大数,然后切割它。
2。尽量避免使用 for
循环。 它们通常很慢,尤其是在遍历大型矩阵时,因为您需要单独索引每个条目。
备选方案: 矩阵运算通常要快得多。例如,不是在两个 for
循环中使用 r = np.sqrt(sum((matposg[j,:]-matposg[k,:])**2))
,而是可以先定义两个矩阵 A
和 B
,它们对应于 matposg[j,:]
和 matposg[k,:]
(应该可以不使用循环!)然后简单地使用 r = np.linalg.norm(A-B)