在 python 中使用 GPU 加速 monte carlo 模拟
Using GPU Speeding up monte carlo simulation in python
我的代码使用 monte carlo 技术来绘制氢原子周围有电子的点,但是当涉及到绘制这个时,程序可能需要几个小时以上才能生成图表 I一直在寻找减少这段时间的方法并尝试使用 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)
刚刚快速玩了一下,发现:
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
开关非常有用。
我的代码使用 monte carlo 技术来绘制氢原子周围有电子的点,但是当涉及到绘制这个时,程序可能需要几个小时以上才能生成图表 I一直在寻找减少这段时间的方法并尝试使用 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)
刚刚快速玩了一下,发现:
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
开关非常有用。