在蒙特卡洛模拟中改进python代码
危彬彬
问题内容：
问题答案：
我已经为“ 2d有效ising模型”编写了蒙特卡洛模拟，并且试图改善运行时间。
我的代码做什么：我为每个粒子（rgrid和mgrid）创建一个用于粒子数量（r）的矩阵和一个用于磁化强度的矩阵。粒子的自旋可以是-1/1，因此磁化强度范围为[-r，r]，步长为2。
然后选择一个随机点和一个随机粒子（+1或-1）。由于概率取决于每个位置的正/负粒子数量，因此我创建了2个数组并将其压缩，这样我就可以得到正粒子的合适数量，即[[-3，0），（-1，1），（
1，2），（3，3）]。使用3个粒子，我的磁化强度为（-3，-1、1、3），其中具有（0、1、2、3）+1个粒子。
之后，我计算该点的概率并选择一个动作：spinflip，向上/向下跳，向左/向右跳，什么也不做。现在，我必须移动粒子（或不移动粒子）并更改两个点的磁体/密度（并检查周期性边界条件）。
这是我的代码：
from __future__ import print_function
from __future__ import division
from datetime import datetime
import numpy as np
import math
import matplotlib.pyplot as plt
import cProfile
pr = cProfile.Profile()
pr.enable()
m = 10 # zeilen, spalten
j = 1000 # finale zeit
r = 3 # platzdichte
b = 1.6 # beta
e = 0.9 # epsilon
M = m * m # platzanzahl
N = M * r # teilchenanzahl
dt = 1 / (4 * np.exp(b)) # delta-t
i = 0
rgrid = r * np.ones((m, m)).astype(int) # dichte-matrix, rho = n(+) + n(-)
magrange = np.arange(-r, r + 1, 2) # mögliche magnetisierungen [a, b), schrittweite
mgrid = np.random.choice(magrange, (m, m)) # magnetisierungs-matrix m = n(+) - (n-)
def flip():
mgrid[math.trunc(x / m), x % m] -= 2 * spin
def up():
y = x - m
if y < 0: # periodische randbedingung hoch
y += m * m
x1 = math.trunc(x / m)
x2 = x % m
y1 = math.trunc(y / m)
y2 = y % m
rgrid[x1, x2] -= 1 # [zeile, spalte] masse -1
rgrid[y1, y2] += 1 # [zeile, spalte] masse +1
mgrid[x1, x2] -= spin # [zeile, spalte] spinänderung alter platz
mgrid[y1, y2] += spin # [zeile, spalte] spinänderung neuer platz
def down():
y = x + m
if y >= m * m: # periodische randbedingung unten
y -= m * m
x1 = math.trunc(x / m)
x2 = x % m
y1 = math.trunc(y / m)
y2 = y % m
rgrid[x1, x2] -= 1
rgrid[y1, y2] += 1
mgrid[x1, x2] -= spin
mgrid[y1, y2] += spin
def left():
y = x - 1
if math.trunc(y / m) < math.trunc(x / m): # periodische randbedingung links
y += m
x1 = math.trunc(x / m)
x2 = x % m
y1 = math.trunc(y / m)
y2 = y % m
rgrid[x1, x2] -= 1
rgrid[y1, y2] += 1
mgrid[x1, x2] -= spin
mgrid[y1, y2] += spin
def right():
y = x + 1
if math.trunc(y / m) > math.trunc(x / m): # periodische randbedingung rechts
y -= m
x1 = math.trunc(x / m)
x2 = x % m
y1 = math.trunc(y / m)
y2 = y % m
rgrid[x1, x2] -= 1
rgrid[y1, y2] += 1
mgrid[x1, x2] -= spin
mgrid[y1, y2] += spin
while i < j:
# 1. platz aussuchen
x = np.random.randint(M) # wähle zufälligen platz aus
if rgrid.item(x) != 0:
i += dt / N
# 2. teilchen aussuchen
li1 = np.arange(-abs(rgrid.item(x)), abs(rgrid.item(x)) + 1, 2)
li2 = np.arange(0, abs(rgrid.item(x)) + 1)
li3 = zip(li1, li2) # list1 und list2 als tupel in list3
results = [item[1] for item in li3 if item[0] == mgrid.item(x)] # gebe 2. element von tupel aus für passende magnetisierung
num = int(''.join(map(str, results))) # wandle listeneintrag in int um
spin = 1.0 if np.random.random() < num / rgrid.item(x) else -1.0
# 3. ereignis aussuchen
p = np.random.random()
p1 = np.exp(- spin * b * mgrid.item(x) / rgrid.item(x)) * dt # flip
p2 = dt # hoch
p3 = dt # runter
p4 = (1 - spin * e) * dt # links
p5 = (1 + spin * e) * dt # rechts
p6 = 1 - (4 + p1) * dt # nichts
if p < p6:
continue
elif p < p6 + p1:
flip()
continue
elif p < p6 + p1 + p2:
up()
continue
elif p < p6 + p1 + p2 + p3:
down()
continue
elif p < p6 + p1 + p2 + p3 + p4:
left()
continue
else:
right()
continue
pr.disable()
pr.print_stats(sort='cumtime')
我已经做了些什么来加快速度：
- 使用cProfile来查看我使用的random.choice非常慢，并将其更改为旋转和概率p的if条件
- 安装了cython，并用于
import pyximport; pyximport.install()创建编译的cython文件。这没有任何改善，经过检查，cython -a script.py我发现我需要更多的静态变量来获得一些改善。
cProfile现在运行向我显示此：
188939207 function calls in 151.834 seconds
Ordered by: cumulative time
ncalls tottime percall cumtime percall filename:lineno(function)
5943639 10.582 0.000 40.678 0.000 {method 'join' of 'str' objects}
5943639 32.543 0.000 37.503 0.000 script.py:107(<listcomp>)
5943639 4.722 0.000 30.096 0.000 numeric.py:1905(array_str)
8276949 25.852 0.000 25.852 0.000 {method 'randint' of 'mtrand.RandomState' objects}
5943639 11.855 0.000 25.374 0.000 arrayprint.py:381(wrapper)
11887279 14.403 0.000 14.403 0.000 {built-in method numpy.core.multiarray.arange}
80651998 13.559 0.000 13.559 0.000 {method 'item' of 'numpy.ndarray' objects}
5943639 8.427 0.000 9.364 0.000 arrayprint.py:399(array2string)
11887278 8.817 0.000 8.817 0.000 {method 'random_sample' of 'mtrand.RandomState' objects}
579016 7.351 0.000 7.866 0.000 script.py:79(right)
300021 3.669 0.000 3.840 0.000 script.py:40(up)
152838 1.950 0.000 2.086 0.000 script.py:66(left)
17830917 1.910 0.000 1.910 0.000 {built-in method builtins.abs}
176346 1.147 0.000 1.217 0.000 script.py:37(flip)
5943639 1.131 0.000 1.131 0.000 {method 'discard' of 'set' objects}
5943639 1.054 0.000 1.054 0.000 {built-in method _thread.get_ident}
5943639 1.010 0.000 1.010 0.000 {method 'add' of 'set' objects}
5943639 0.961 0.000 0.961 0.000 {built-in method builtins.id}
3703804 0.892 0.000 0.892 0.000 {built-in method math.trunc}
我用join得到关于点+1的粒子数的整数值，因为我需要为我的概率。
如果我要运行什么大不了的事情一样m=400，r=3，j=300000（记者：最后一次），我期待在我目前的速度大约4年运行时间。
任何帮助是极大的赞赏。
问题答案：
蒙特卡洛模拟
首先，我摆脱了列表，然后使用了即时编译器（numba）。如果不编译，则得到196s（您的版本），经过编译，我得到0.44s。因此，通过使用jit编译器可以改善
435 倍，这里还有一些简单的优化。
一个主要的优点是，这里也释放了GIL（全局解释器锁）。如果代码受处理器限制且不受内存带宽限制，则可以在另一个线程中运行随机数时在另一个线程中计算随机数（可以使用多个内核）。这也可以进一步提高性能，并且可以按以下方式工作：
- 计算第一个随机数块（足够小，以至于整个问题很容易放入处理器缓存中（至少为L3缓存）。
- 使用该随机数开始迭代。在运行迭代时，将计算另一块随机数。
- 重复（2），直到完成。
码
import numba as nb
import numpy as np
def calc (m,j,e,r,dt,b,rgrid,mgrid):
M=m*m
N = M * r
i=0
while i < j:
# 1. platz aussuchen
x = np.random.randint(M) # wähle zufälligen platz aus
if rgrid[x] != 0:
i += dt / N
########
# 2. teilchen aussuchen
#li1 = np.arange(-abs(rgrid[x]), abs(rgrid[x]) + 1, 2)
#li2 = np.arange(0, abs(rgrid[x]) + 1)
#li3 = zip(li1, li2) # list1 und list2 als tupel in list3
#results = [item[1] for item in li3 if item[0] == mgrid[x]] # gebe 2. element von tupel aus für passende magnetisierung
#num = int(''.join(map(str, results))) # wandle listeneintrag in int um
#######
# This should be equivalent
jj=0 #li2 starts with 0 and has a increment of 1
for ii in range(-abs(rgrid[x]),abs(rgrid[x])+1, 2):
if (ii==mgrid[x]):
num=jj
break
jj+=1
spin = 1.0 if np.random.random() < num / rgrid[x] else -1.0
# 3. ereignis aussuchen
p = np.random.random()
p1 = np.exp(- spin * b * mgrid[x] / rgrid[x]) * dt # flip
p2 = dt # hoch
p3 = dt # runter
p4 = (1 - spin * e) * dt # links
p5 = (1 + spin * e) * dt # rechts
p6 = 1 - (4 + p1) * dt # nichts
if p < p6:
continue
elif p < p6 + p1:
#flip()
mgrid[x] -= 2 * spin
continue
elif p < p6 + p1 + p2:
#up()
y = x - m
if y < 0: # periodische randbedingung hoch
y += m * m
rgrid[x] -= 1 # [zeile, spalte] masse -1
rgrid[y] += 1 # [zeile, spalte] masse +1
mgrid[x] -= spin # [zeile, spalte] spinänderung alter platz
mgrid[y] += spin # [zeile, spalte] spinänderung neuer platz
continue
elif p < p6 + p1 + p2:
#down()
y = x + m
if y >= m * m: # periodische randbedingung unten
y -= m * m
rgrid[x] -= 1
rgrid[y] += 1
mgrid[x] -= spin
mgrid[y] += spin
continue
elif p < p6 + p2 + p3:
#left()
y = x - 1
if (y // m) < (x // m): # periodische randbedingung links
y += m
rgrid[x] -= 1
rgrid[y] += 1
mgrid[x] -= spin
mgrid[y] += spin
continue
else:
#right()
y = x + 1
if (y // m) > (x // m): # periodische randbedingung rechts
y -= m
rgrid[x] -= 1
rgrid[y] += 1
mgrid[x] -= spin
mgrid[y] += spin
continue
return (mgrid,rgrid)
现在测试的主要功能是：
def main():
m = 10 # zeilen, spalten
j = 1000 # finale zeit
r = 3 # platzdichte
b = 1.6 # beta
e = 0.9 # epsilon
M = m * m # platzanzahl
N = M * r # teilchenanzahl
dt = 1 / (4 * np.exp(b)) # delta-t
i = 0
rgrid = r * np.ones((m, m),dtype=np.int) #don't convert the array build it up with the right datatype # dichte-matrix, rho = n(+) + n(-)
magrange = np.arange(-r, r + 1, 2) # mögliche magnetisierungen [a, b), schrittweite
mgrid = np.random.choice(magrange, (m, m)) # magnetisierungs-matrix m = n(+) - (n-)
#Compile the function
nb_calc = nb.njit(nb.types.Tuple((nb.int32[:], nb.int32[:]))(nb.int32, nb.int32,nb.float32,nb.int32,nb.float32,nb.float32,nb.int32[:], nb.int32[:]),nogil=True)(calc)
Results=nb_calc(m,j,e,r,dt,b,rgrid.flatten(),mgrid.flatten())
#Get the results
mgrid_new=Results[0].reshape(mgrid.shape)
rgrid_new=Results[1].reshape(rgrid.shape)
编辑 我已经重写了代码“ 2.Teilchen aussuchen”，并对代码进行了重新设计，使其可用于标量索引。这样可以使速度额外提高4倍。
类似资料：
-
1 前言 在上一篇文章中，我们介绍了基于Bellman方程而得到的Policy Iteration和Value Iteration两种基本的算法，但是这两种算法实际上很难直接应用，原因在于依然是偏于理想化的两个算法，需要知道状态转移概率，也需要遍历所有的状态。对于遍历状态这个事，我们当然可以不用做到完全遍历，而只需要尽可能的通过探索来遍及各种状态即可。而对于状态转移概率，也就是依赖于模型Model
-
个人觉得，整个 AplphaGo 对于机器学习来说，最核心的算法就是深度学习（Deep Learning）和增强学习（Reinforcement Learning）。蒙特卡洛树搜索 MCTS 是一个搜索框架，将这些机器学习的技术融合在了一起。今天这篇文章的重点在深度学习，增强学习以后再说。 蒙特卡洛树搜索 每个博弈类的人工智能算法的基础都是一个搜索算法。比如我们上学时学习的 A-star 算法，a
-
我试图模拟来自rho=0.7的AR（1）模型的数据（Y）。然后我将使用这些数据在截距上运行Y的回归（通过这样做，参数估计成为Y的平均值），然后使用鲁棒的标准错误。我想对这个假设运行一个蒙特卡罗模拟，使用2000次重复不同的滞后值。目的是显示当滞后变化时Newey West估计器的有限样本性能 我的问题是：上面的代码是进行这种模拟的正确方法吗？如果是，我如何得到一个代码来重复这个过程在HAC测试中的
-
我正在写一个蒙特卡罗模拟来检查有多少次y不是紧挨着另一个y。我变出了一个40 x和10 y的向量，放置在向量中的随机位置。我的目标是计算向量中没有任何相邻y的概率。以下是我尝试过的： 结果是一个非常小的数字，这对我来说似乎没有意义。
-
TLDR MCTS代理实现在本地无错误运行，实现了 嗨，我目前正忙于这个项目，我需要在两周内完成我注册的项目之前完成这个项目。我的任务，我已经完成了，是实现一个代理，在两个国际象棋骑士之间的隔离游戏中，与启发式驱动的minimax代理对抗。游戏的完整实现细节可以在这里找到。在我的项目中，游戏将在一个9 x 11的棋盘上进行，使用位棋盘编码。我的MCTS实现非常简单，紧跟本文（第6页）提供的伪代码。