蒙特卡洛模拟
首先,我摆脱了列表,然后使用了即时编译器(numba)。如果不编译,则得到196s(您的版本),经过编译,我得到0.44s。因此,通过使用jit编译器可以改善
435 倍,这里还有一些简单的优化。
一个主要的优点是,这里也释放了GIL(全局解释器锁)。如果代码受处理器限制且不受内存带宽限制,则可以在另一个线程中运行随机数时在另一个线程中计算随机数(可以使用多个内核)。这也可以进一步提高性能,并且可以按以下方式工作:
- 计算第一个随机数块(足够小,以至于整个问题很容易放入处理器缓存中(至少为L3缓存)。
- 使用该随机数开始迭代。在运行迭代时,将计算另一块随机数。
- 重复(2),直到完成。
码
import numba as nbimport numpy as npdef 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倍。



