python - 使用GPU加速Python中的蒙特卡罗模拟

Question

我的代码使用蒙特卡罗技术来绘制在氢原子周围发现电子的点，但是在绘制该点时，程序可能需要几个小时以上才能生成我一直在查看的图表减少这个时间的方法并尝试使用 numba 我不知道我是否正确使用它，但它没有任何区别

有人可以帮我吗

from scipy.special import genlaguerre, sph_harm
from numpy import random, linspace, sqrt, pi, arccos, exp
from matplotlib.pyplot import plot, figure, show
from numba import jit

n = 2
l = 1
m = -1

Lx = 3*10**-9
Ly = Lx

a_0 = 5.29*10**-11

x = linspace(0,Lx, 100000)
fi = linspace(0,2*pi, 100000)


def fact(x):
    if x == 0:
        return 1
    else:
        return x*fact(x-1)

def row(r):
    return (2*r)/(n*a_0)

def Hydrograd(r):
    return -1*sqrt((2/(n*a_0))**3*(fact(n - l- 1))/(2*n*(fact(n+l))**3))*exp(-1*row(r)/2)*row(r)**l*genlaguerre(n-l-1, 2*l+1)(row(r))

def Hydrogsph(theta, phi):
    return sph_harm(m,l,theta,phi).real

@jit
def HydroFunc(r,theta, phi):
    return abs(Hydrograd(r))**2*abs(Hydrogsph(theta, phi))**2

def rFinder(x,y,z):
    return sqrt(x**2 + y**2 + z**2)

def phiFinder(x,y,z):
    return arccos(x/rFinder(x,y,z))

i = 0
p = 50
@jit
def monty(Lx,Ly,x,fi):
    i = 0
    p = 50
    X = []
    Y = []
    while i < p:
        xr = random.uniform(-Lx,Lx)
        yr = random.uniform(-Ly,Ly)
        if HydroFunc(rFinder(xr,yr,0),0,phiFinder(xr,yr,0))>random.uniform(0,max(HydroFunc(x,0,fi))/100):
            X.append(xr)
            Y.append(yr)
            i += 1
            print(i)
    figure()
    plot(X,Y,'.')
    show()

monty(Lx,Ly,x,fi)

Answer 1

刚刚快速玩了一下并注意到:

max(HydroFunc(x,0,fi))/100

看起来不变。

因此，您可以更改 monty 函数来执行以下操作:

def monty(Lx,Ly,x,fi):
    rmax = max(HydroFunc(x,0,fi)) / 100
    # ...
        if HydroFunc(rFinder(xr,yr,0),0,phiFinder(xr,yr,0)) > random.uniform(0, rmax):

这对我来说加快了很多速度。我注意到你正在做其他多余的计算，例如当您使用相同的参数调用外部时，phiFinder 会调用 rFinder。

您的拒绝率也非常高，您也许可以考虑缩小提案分配范围。马尔可夫链算法也可能有帮助，例如:

def mcmc(sigma, nsamples):
    rmax = max(HydroFunc(x, 0, fi)) / 100
    x1, y1 = 1e-9, 1e-9
    p1 = min(rmax, HydroFunc(rFinder(x1,y1,0),0,phiFinder(x1,y1,0)))
    for _ in range(nsamples):
        x2, y2 = random.normal([x1, y1], scale=sigma)
        p2 = min(rmax, HydroFunc(rFinder(x2,y2,0),0,phiFinder(x2,y2,0)))
        if p2 / p1 > random.uniform():
            x1, y1, p1 = x2, y2, p2
        yield x1, y1

在大约 15 秒内给我 100k(自动相关)样本，而从您的方法中获得相同的样本需要大约 350 秒(将 rmax 计算移出内部循环之后) 。据我所知，这些样本来自相同的分布，例如我运行并绘制:

out = np.array(list(mcmc(8e-10, 100000)))

# thin out the chain to reduce autocorrelation
ii = range(0, len(out), 5)
plt.scatter(out[ii,0], out[ii,1], 10, edgecolor='none', alpha=0.1)
plt.xlim(-1e-9, 1e-9)
plt.ylim(-1e-9, 1e-9)

给我:

(橙色点是您的原始实现中的 50 个样本，花费了 30 秒)

我刚刚更新到 LLVM 9，所以 numba 目前对我不起作用，但我建议查看 https://numba.pydata.org/numba-doc/dev/user/5minguide.html弄清楚它在做什么。我之前发现 nopython=True 开关非常有用。