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

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)

标签： pythonmontecarlo

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开关非常有用。