加速Ising/PottsModel MonteC的Python/numy代码

Question

对于一个小型项目，我用Python/Numpy编写了一个2D Ising-/Potts模型montecarlo模拟，下面是（简化的）代码。

基本上，该代码执行以下操作：

• 在[1，orientations]中设置一个用随机整数（orientations）填充的NXM数组
• 对于每个时间步（MCS），数组的每个像素都将以伪随机顺序访问一次（因此，在（）和 index = (a*i + rand) % (N**2) x = index % N y = index // N）
• 检查8个相邻数组条目（周期性边界条件），并将条目更改为相邻值之一
• 如果新配置的能量变低，则接受更改，否则将拒绝更改，除非满足某个条件

我试图尽我所能加快它的速度，但是对于大型数组（N，M>；500），代码并不是很快。因为我需要大约10e5 MCS的阵列才能看到明显的趋势，实现了

1 loop, best of 3: 276 ms per loop

对于100x100阵列上的1个MCS来说，这还不够。不幸的是，由于缺乏经验，我不知道如何提高性能。

我假设Neighbors（）和calc_dE（）函数是瓶颈，尤其是嵌套循环，但我找不到加快速度的方法。我的cython尝试并不是很成功，因为我以前从未对cython做过任何事情，所以我愿意接受任何建议。

代码：

（pyplot命令仅用于可视化，通常会添加注释。）

^{pr2}$

Answer 1

不确定在依赖关系方面是否有任何限制，但我肯定会研究Numba。它提供了一组修饰符（njit），如果您使用兼容的类型（例如numpy数组，而不是pandas数据帧），可以将代码编译为机器代码，并使其速度更快。在

另外，不知道你在看什么样的规模，但我相信你可以找到比手动实现for循环更好的优化素数搜索算法的例子。在

否则，您总是可以依靠Cython，但这需要重新编写代码。在

---

编辑：好的，我用numba试一试。在

一些注意事项：

1. 在函数中移动了整个for循环，以便可以用njit来修饰它
2. 在Neighbors中，我不得不将rows和{}从list改为{}，因为numba不接受通过列表建立索引
3. 我将np.random.random_integers替换为np.random.randint，因为前者已被弃用
4. 我将math.exp替换为np.exp，这将稍微提高性能（除了节省导入）

import numpy as np
from numba import njit

def largest_primes_under(N):
    n = N - 1
    while n >= 2:
        if all(n % d for d in range(2, int(n ** 0.5 + 1))):
            return n
        n -= 1

@njit
def Neighbors(Lattice,i,j,n=1):
    ''' Returns an flat array of all neighboring sites in the n-th coordination sphere including the center'''
    N, M = Lattice.shape
    rows = np.array([(i-1) % N, i, (i+1) % N])
    cols = np.array([(j-1) % N, j, (j+1) % M])
    return Lattice[rows,:][:,cols].flatten()

@njit
def calc_dE(Lattice, x, y, z):
    N, M = Lattice.shape
    old_energy = 0
    new_energy = 0
    for i in [0,1,-1]:
        for j in [0,1,-1]:
            if i == 0 and j == 0: 
                continue
            if Lattice[x%N,y%M] == Lattice[(x+i)%N,(y+j)%M]:
                old_energy += 1
            elif z == Lattice[(x+i)%N,(y+j)%M]: 
                new_energy += 1 
    return old_energy-new_energy

@njit
def fun(L, MCS, a):
    N, M = L.shape

    for t in range(1,MCS+1):
        rand = np.random.randint(N*M)
        for i in range(0,N**2):
            index = (a*i + rand) % (N**2)
            x = index % N
            y = index // N
            n = Neighbors(L,x,y)
            if len(n)-1 == 0: continue
            else: z = np.random.choice(n)
            dE = calc_dE(L,x,y,z)
            if  (dE < 0): L[x%N,y%N] = z      
            elif np.random.sample() < np.exp(-dE*2.5): L[x%N,y%N] = z    
    return L

运行相同的例子

^{pr2}$

通过%timeit fun(L, MCS, a)（在Jupyter）给了我

16.9 ms ± 853 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

比你现在的速度快15倍。可能还有更多的优化可以做，numba的好处是我在没有深入研究或显著改变代码实现方式的情况下获得了x15的加速。在

一些一般性的观察：

1. 在Neighbors中，参数/参数n没有使用，因此为了清晰起见，应该删除它（或更新代码）
2. 在Python中，通常希望对函数名和变量使用小写。大写它通常保留给类（不是对象）和“全局”变量
3. 您上面关于largest_primes_under只被调用一次的评论绝对是正确的，我应该更仔细地看一下代码。在

premature optimization is the root of all evil