montecarlo和Metropolis算法在Python中非常慢

Question

我正在尝试用Python实现一个简单的montecarlo（我对它还比较陌生）。来自C语言的我可能走的是最错误的路径，因为我的代码对于我的要求来说太慢了：我有一个潜在的硬球（见代码中的V_pot(r)），可以容纳60个3d粒子和周期性边界条件（PBC），所以我定义了以下函数

import timeit
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
from numpy import inf
#
L, kb, d, eps, DIM = 100, 1, 1, 1, 3
r_c, T = L/2, eps/(.5*kb)
beta = 1/(kb*T)
#
def dist(A, B):
    d = A - B
    d -= L*np.around(d/L)
    return np.sqrt(np.sum(d**2))

#
def V_pot(r):
    V = -eps*(d**6/r**6 - d**6/r_c**6)
    if r > r_c:
        V = 0
    elif r < d:
        V = inf

    return V

#
def ener(config):
    V_jk_val, j = 0, N
    #
    while (j > 0):
        j -= 1
        i = 0
        while (i < j):
            V_jk_val += V_pot(dist(config[j,:], config[i,:]))
            i += 1

    #
    return V_jk_val

#
def acc(en_n, en_o):
    d_en = en_n-en_o
    if (d_en <= 0):
        acc_val = 1
    else:
        acc_val = np.exp(-beta*(d_en))

    return acc_val

#

然后，从配置开始（数组的每一行代表一个三维粒子的坐标）

config = np.array([[16.24155657, 57.41672173, 94.39565792],
       [76.38121764, 55.88334066,  5.72255163],
       [38.41393783, 58.09432145,  6.26448054],
       [86.44286438, 61.37100899, 91.97737383],
       [37.7315366 , 44.52697269, 23.86320444],
       [ 0.59231801, 39.20183376, 89.63974115],
       [38.00998141,  3.84363202, 52.74021401],
       [99.53480756, 69.97688928, 21.43528924],
       [49.62030291, 93.60889503, 15.73723259],
       [54.49195524,  0.6431965 , 25.37401196],
       [33.82527814, 25.37776021, 67.4320553 ],
       [64.61952893, 46.8407798 ,  4.93960443],
       [60.47322732, 16.48140136, 33.26481306],
       [19.71667792, 46.56999616, 35.61044526],
       [ 5.33252557,  4.44393836, 60.55759256],
       [44.95897856,  7.81728046, 10.26000715],
       [86.5548395 , 49.74079452,  4.80480133],
       [52.47965686, 42.831448  , 22.03890639],
       [ 2.88752006, 59.84605062, 22.75760029],
       [ 9.49231045, 42.08653603, 40.63380097],
       [13.90093641, 74.40377984, 32.62917915],
       [97.44839233, 90.47695772, 91.60794836],
       [51.29501624, 27.03796277, 57.09525454],
       [10.30180295, 21.977336  , 69.54173272],
       [59.61327648, 14.29582325, 11.70942289],
       [89.52722796, 26.87758644, 76.34934637],
       [82.03736088, 78.5665713 , 23.23587395],
       [79.77571695, 66.140968  , 53.6784269 ],
       [82.86070472, 40.82189833, 51.48739072],
       [99.05647523, 98.63386809,  6.33888993],
       [31.02997123, 66.99709163, 95.88332332],
       [97.71654767, 59.24793618,  5.20183793],
       [ 6.79964473, 45.01258652, 48.69477807],
       [93.34977049, 55.20537774, 82.35693526],
       [17.35577815, 20.45936211, 29.27981422],
       [55.51942207, 52.22875901,  3.6616131 ],
       [61.45612224, 36.50170405, 62.89796773],
       [23.55822368,  7.09069623, 37.38274914],
       [39.57082799, 58.95457592, 48.0304924 ],
       [93.94997617, 64.34383203, 77.63346308],
       [17.47989107, 90.01113402, 81.00648645],
       [86.79068539, 66.35768515, 56.64402907],
       [98.71924121, 38.33749023, 73.4715132 ],
       [ 0.42356139, 78.32172925, 15.19883322],
       [77.75572529,  2.60088767, 56.4683935 ],
       [49.76486142,  3.01800153, 93.48019286],
       [42.54483899,  4.27174457,  4.38942325],
       [66.75777178, 41.1220603 , 19.64484167],
       [19.69520773, 41.09230171,  2.51986091],
       [73.20493772, 73.16590392, 99.19174281],
       [94.16756184, 72.77653334, 10.32128552],
       [29.95281655, 27.58596604, 85.12791195],
       [ 2.44803886, 32.82333962, 41.6654683 ],
       [23.9665915 , 49.94906612, 37.42701059],
       [30.40282934, 39.63854309, 47.16572743],
       [56.04809276, 30.19705527, 29.15729635],
       [ 2.50566522, 70.37965564, 16.78016719],
       [28.39713572,  4.04948368, 27.72615789],
       [26.11873563, 41.49557167, 14.38703697],
       [81.91731981, 12.10514972, 12.03083427]])

我用下面的代码进行了5000个时间步的模拟

N = 60
TIME_MC = 5000
DELTA_LIST = [d]
#d/6, d/3, d, 2*d, 3*d 
np.random.seed(19680801)
en_mc_delta = np.zeros((TIME_MC, len(DELTA_LIST)))
start = timeit.default_timer()
config_tmp = config
#
for iD, Delta in enumerate(DELTA_LIST):
    t=0

    while (t < TIME_MC):
        for k in range(N):
            RND = np.random.rand()
            config_tmp[k,:] = config[k,:] + Delta*(np.random.random_sample((1,3))-.5)
            en_o, en_n = ener(config), ener(config_tmp)
            ACC = acc(en_n, en_o)
            if (RND < ACC):
                config[k,:] = config_tmp[k,:]
                en_o = en_n


        en_mc_delta[t][iD] = en_o
        t += 1


stop = timeit.default_timer()
print('Time: ', stop-start)

遵循Metropolis算法的规则，接受用config_tmp[k,:] = config[k,:] + Delta*(np.random.random_sample((1,3))-.5)提取的提议的移动。你知道吗

我做了一些尝试来检查代码被卡住的地方，我发现函数ener（也因为函数dist）非常慢：它需要~0.02s之类的东西来计算一个配置的能量，这意味着~6000s周围的东西来运行完整的模拟（60个粒子，5000个提议的移动）。你知道吗

它的外层只是计算Delta的不同值的结果。你知道吗

用TIME_MC=60运行这段代码可以让您了解这段代码（~218s）有多慢，如果用C实现的话，这只需要几秒钟。我读了一些关于如何加速Python代码的其他问题，但我不明白如何在这里做。你知道吗

编辑：

我现在几乎可以肯定问题出在函数dist上，因为仅仅计算两个3D向量之间的PBC距离就需要~0.0012s，当你计算5000*60次的时候，它会给出很长的时间。你知道吗

Answer 1

请注意，这是对原问题的部分回答。你知道吗

下面是一个例子，“展开”numpy的函数如何在替换为更直接的距离计算时提高性能。请注意，这并没有被证实是等效的，特别是在四舍五入方面。我认为这个原则仍然适用。你知道吗

import random
import time
import numpy as np

L = 100
inv_L = 0.01
vec_length = 10
repetitions = 100000

def dist_np(A, B):
    d = A - B
    d -= L*np.around(d/L)
    return np.sqrt(np.sum(d**2))

def dist_direct(A, B):
    sum = 0
    for i in range(0, len(A)):
        diff = (A[0,i] - B[0,i])
        diff -= L * int(diff * inv_L)
        sum += diff * diff
    return np.sqrt(sum)

vec1 = np.zeros((1,vec_length))
vec2 = np.zeros((1,vec_length))

for i in range(0, vec_length):
    vec1[0,i] = random.random()
    vec2[0,i] = random.random()

print("with numpy method:")
start = time.time()
for i in range(0, repetitions):
    dist_np(vec1, vec2)
print("done in {}".format(time.time() - start))

print("with direct method:")
start = time.time()
for i in range(0, repetitions):
    dist_direct(vec1, vec2)
print("done in {}".format(time.time() - start))

输出：

with numpy method:
done in 6.332799911499023
with direct method:
done in 1.0938000679016113

利用平均向量长度和重复次数来确定最佳点在哪里。我希望在改变这些元参数时，性能增益不是恒定的。你知道吗