Python/Numpy - 加速放射性衰變的蒙特卡洛方法

Question

我正試圖優化放射性同位素蒙特卡羅衰變時間的產生。 給出了nsims同位素原子，其半衰期爲t12，每個同位素衰變何時發生？ 我試圖通過使用單個numpy.random.random調用（我稱這種方法並行）爲所有未衰減的原子一次生成隨機數來優化這一點，但我希望還有更多的性能可以獲得。我還展示了一種爲每種同位素分別進行計算的方法（系列）。

import numpy as np 
import time 
import matplotlib.pyplot as plt 
import scipy.optimize 

t12 = 3.1*60. 
dt = 0.01 
ln2 = np.log(2) 
decay_exp = lambda t,A,tau: A * np.exp(-t/tau) 

def serial(nsims): 
    sim_start_time = time.clock() 
    decay_time = np.zeros(nsims) 
    for i in range(nsims): 
     t = dt 
     while decay_time[i] == 0: 
      if np.random.random() > np.exp(-ln2*dt/t12): 
       decay_time[i] = t 
      t += dt 
    sim_end_time = time.clock() 
    return (sim_end_time - sim_start_time,decay_time) 

def parallel(nsims): 
    sim_start_time = time.clock() 
    decay_time = np.zeros(nsims) 
    t = dt 
    while 0 in decay_time: 
     inot_decayed = np.where(decay_time == 0)[0] 
     idecay_check = np.random.random(len(inot_decayed)) > np.exp(-ln2*dt/t12) 
     decay_time[inot_decayed[np.where(idecay_check==True)[0]]] = t 
     t += dt 
    sim_end_time = time.clock() 
    return (sim_end_time - sim_start_time,decay_time) 

我感興趣的執行比parallel函數是純Python更好的任何建議，即不是用Cython。 這種方法已經在serial大nsims計算這個方法已經大大改善。

Answer 1

還有一些速度收益從原來的「水貨」有（矢量化是正確的單詞）執行。

備註，這是微管理，但它仍然會帶來小的性能提升。

import numpy as np 
t12 = 3.1*60. 
dt = 0.01 
ln2 = np.log(2) 

s = 98765 

def parallel(nsims): # your code, unaltered, except removed inaccurate timing method 
    decay_time = np.zeros(nsims) 
    t = dt 
    np.random.seed(s) # also had to add a seed to get comparable results 
    while 0 in decay_time: 
     inot_decayed = np.where(decay_time == 0)[0] 
     idecay_check = np.random.random(len(inot_decayed)) > np.exp(-ln2*dt/t12) 
     decay_time[inot_decayed[np.where(idecay_check==True)[0]]] = t 
     t += dt 
    return decay_time 

def parallel_micro(nsims): # micromanaged code 
    decay_time = np.zeros(nsims) 
    t = dt 
    half_time = np.exp(-ln2*dt/t12) # there was no need to calculate this again in every loop iteration 
    np.random.seed(s) # fixed seed to get comparable results 
    while 0 in decay_time: 
     inot_decayed = np.where(decay_time == 0)[0] # only here you need the call to np.where 
     # to my own surprise, len(some_array) is quicker than some_array.size (function lookup vs attribute lookup) 
     idecay_check = np.random.random(len(inot_decayed)) > half_time 
     decay_time[inot_decayed[idecay_check]] = t # no need for another np.where and certainly not for another boolean comparison 
     t += dt 
    return decay_time 

您可以使用timeit module運行計時測量。性能分析將告訴您，此處的瓶頸是致電np.where。

明知瓶頸是np.where，你可以擺脫這樣的：

def parallel_micro2(nsims): 
    decay_time = np.zeros(nsims) 
    t = dt 
    half_time = np.exp(-ln2*dt/t12) 
    np.random.seed(s) 
    indices = np.where(decay_time==0)[0] 
    u = len(indices) 
    while u: 
     decayed = np.random.random(u) > half_time 
     decay_time[indices[decayed]] = t 
     indices = indices[np.logical_not(decayed)] 
     u = len(indices) 
     t += dt 
    return decay_time 

這確實給相當大的速度增加：

In [2]: %timeit -n1 -r1 parallel_micro2(1e4) 
1 loops, best of 1: 7.81 s per loop 

In [3]: %timeit -n1 -r1 parallel_micro(1e4) 
1 loops, best of 1: 29 s per loop 

In [4]: %timeit -n1 -r1 parallel(1e4) 
1 loops, best of 1: 33.5 s per loop 

不要忘記讓當你完成優化時，擺脫np.random.seed的呼叫。