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

Question

我正在尝试优化放射性同位素 Monte Carlo 衰变时间的生成。给定 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)

我对任何比纯 python 的 parallel 函数表现更好的建议感兴趣，即不是 cython。此方法已经大大改进了为大型 nsims 计算此值的 serial 方法。

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 的调用。