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

【问题标题】：Using GPU Speeding up monte carlo simulation in python在 python 中使用 GPU 加速蒙特卡罗模拟
【发布时间】：2020-03-26 08:36:16
【问题描述】：

我的代码使用蒙特卡洛技术来绘制在氢原子周围可以找到电子的点，但是在绘制这个时，程序可能需要几个小时才能生成我一直在查看的图表减少这个时间的方法并尝试使用 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)

【问题讨论】：

标签： python montecarlo

【解决方案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 开关非常有用。

【讨论】：