我有一个数组（N = 10 ^ 4），需要找到每两个条目之间的差异（根据原子的坐标计算电位）。这是我用纯python编写的代码，但效率不高，有人能告诉我如何加速它吗？（使用numpy或weave）。这里的x，y是原子坐标的数组（只是简单的一维数组）
def potential(r):
U = 4.*(np.power(r,-12) - np.power(r,-6))
return U
def total_energy(x):
E = 0.
#need to speed up this part
for i in range(N-1):
for j in range(i):
E += potential(np.sqrt((x[i]-x[j])**2))
return E
首先，您可以使用数组算术运算：
def potential(r):
return 4.*(r**(-12) - r**(-6))
def total_energy(x):
E = 0.
for i in range(N-1):
E += potential(np.sqrt((x[i]-x[:i])**2)).sum()
return E
或者您可以测试完全向量化的版本：
def total_energy(x):
b=np.diag(x).cumsum(1)-x
return potential(abs(b[np.triu_indices_from(b,1)])).sum()
return sum( potential(np.sqrt((x[i]-x[:i])**2)).sum() for i in range(N-1))。 - lejlotx 是一维的，而根据 OP 的描述，他可能在描述一个二维坐标数组。 - Danielpdist计算数组的所有成对距离。def potential(r):
U = 4.*(np.power(r,-12) - np.power(r,-6))
return U
def total_energy(x):
E = 0.
#need to speed up this part
for i in range(N): #To N here
for j in range(i):
E += potential(np.sqrt(np.sum((x[i]-x[j])**2))) #Add sum here
return E
import scipy.spatial.distance as sd
def scipy_LJ(arr, sigma=None):
"""
Computes the Lennard-Jones potential for an array (M x N) of M points
in N dimensional space. Usage of a sigma parameter is optional.
"""
if len(arr.shape)==1:
arr = arr[:,None]
r = sd.pdist(arr)
if sigma==None:
np.power(r, -6, out=r)
return np.sum(r**2 - r)*4
else:
r *= sigma
np.power(r, -6, out=r)
return np.sum(r**2 - r)*4
让我们运行一些测试：
N = 1000
points = np.random.rand(N,3)+0.1
np.allclose(total_energy(points), scipy_LJ(points))
Out[43]: True
%timeit total_energy(points)
1 loops, best of 3: 13.6 s per loop
%timeit scipy_LJ(points)
10 loops, best of 3: 24.3 ms per loop
N = 10000
points = np.random.rand(N,3)+0.1
%timeit scipy_LJ(points)
1 loops, best of 3: 3.05 s per loop
这个程序使用了2GB的内存。
我是一个有用的助手，可以翻译文本。
这里是最终答案和一些时间记录。
0) 普通版本（非常慢）
In [16]: %timeit total_energy(points)
1 loops, best of 3: 14.9 s per loop
In [9]: %timeit scipy_LJ(points)
10 loops, best of 3: 44 ms per loop
1-2) Numpy 版本
%timeit sum( potential(np.sqrt((x[i]-x[:i])**2 + (y[i]-y[:i])**2 + (z[i] - z[:i])**2)).sum() for i in range(N-1))
10 loops, best of 3: 126 ms per loop
2）极快的Fortran版本（!表示注释）
subroutine EnergyForces(Pos, PEnergy, Dim, NAtom)
implicit none
integer, intent(in) :: Dim, NAtom
real(8), intent(in), dimension(0:NAtom-1, 0:Dim-1) :: Pos
! real(8), intent(in) :: L
real(8), intent(out) :: PEnergy
real(8), dimension(Dim) :: rij, Posi
real(8) :: d2, id2, id6, id12
real(8) :: rc2, Shift
integer :: i, j
PEnergy = 0.
do i = 0, NAtom - 1
!store Pos(i,:) in a temporary array for faster access in j loop
Posi = Pos(i,:)
do j = i + 1, NAtom - 1
rij = Pos(j,:) - Posi
! rij = rij - L * dnint(rij / L)
!compute only the squared distance and compare to squared cut
d2 = sum(rij * rij)
id2 = 1. / d2 !inverse squared distance
id6 = id2 * id2 * id2 !inverse sixth distance
id12 = id6 * id6 !inverse twelvth distance
PEnergy = PEnergy + 4. * (id12 - id6)
enddo
enddo
end subroutine
调用之后
In [14]: %timeit ljlib.energyforces(points.transpose(), 3, N)
10000 loops, best of 3: 61 us per loop
3) 结论 Fortran比scipy快1000倍，比numpy快3000倍，比纯Python快数百万倍。这是因为Scipy版本创建了一个差异矩阵然后分析它，而Fortran版本则一切都在运行时完成。
感谢您的帮助。这是我找到的内容。
最短的版本
return sum( potential(np.sqrt((x[i]-x[:i])**2)).sum() for i in range(N-1))
Scipy版本也不错。
最快的版本是考虑使用f2py程序，即将瓶颈缓慢部分纯Fortran编写（非常快），编译它，然后将其作为库插入到您的Python代码中。 例如，我有： program_lj.f90
$gfortran -c program_lj.f90
如果在Fortran程序中显式地定义了所有类型，那么我们就可以继续。
$f2py -c -m program_lj program_lj.f90
编译后，唯一剩下的事情就是从Python中调用该程序。 在Python程序中：
import program_lj
result = program_lj.subroutine_in_program(parameters)
i只到N-1？在这种情况下，您将错过最后一个条目（N-1被排除在外）。x是否总是一维的，还是需要进行泛化处理？ - eickenberg