我已经编写了一个random_walk
模拟,使用numpy
分配数据,使用Generator执行模拟步骤。这个random_walk
只是原始代码中的MWE(与随机行走完全无关,但它是一个太大和复杂的随机数学模型,无法用作示例。然而,random_walk
MWE模拟了核心组件。
我使用生成器的原因与模拟有关。我将在无限长的时间内运行模拟,并且我只在某些角落的情况下转储数据。我可以测量角点情况发生的概率,因此我能够以高精度预先分配numpy数组(从未出错(,但这是一个上限,因此我需要计算角点情况出现的次数,然后对数据集进行切片(这在模拟中是"模拟的"(。
为了进行比较,我还编写了一种类似的天真方法,使用普通的append
列表来存储模拟数据,只有在出现拐角情况时才添加。
重要的是要知道,拐角情况将发生十亿次(将占用很大一部分内存(,但最终模拟将运行"无限时间",这是一个非常大的步骤。角落的情况就像发生了1e-10问题。
最后的代码有一个停止条件,我在这里用distance
和经典的simulation time
模拟了这个条件。
令我惊讶的是,我注意到append
方法比numpy+generators
具有更好的性能。正如我们在下面的输出中看到的那样
对于小型数据集:
%timeit random_walk_naive(max_distance=1e5, simul_time=1e4)
5.35 ms ± 190 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit random_walk_simul(max_distance=1e5, simul_time=1e4)
16.3 ms ± 567 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
对于较大的数据集
%timeit random_walk_naive(max_distance=1e12, simul_time=1e7)
12.2 s ± 760 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit random_walk_simul(max_distance=1e12, simul_time=1e7)
36 s ± 102 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
通过调用运行cProfile
,我注意到生成器调用的执行时间与天真方法相似,所有额外的时间都花在了random_walk_simul
本身上。在评估np.empty
和切片操作时,我注意到创建空数据集并对其进行切片作为回报的时间很短。根本没有贡献在行动中花费的时间。此外,代码几乎相同,只是在generator
方法中,我首先将数据分配给局部变量,然后将它们"刷新"到numpy.array
,这比直接刷新更快,因为我将使用while
循环中的值来评估停止条件。
我需要了解为什么会出现这种行为,以及是否是意料之中的事;如果不怎么解决?
完整的源代码粘贴在下面
import numpy as np
from random import random
def clip(value, lower, upper):
return lower if value < lower else upper if value > upper else value
def random_walk(s_0, a_0, pa, pb):
"""Initial position (often 0), acceleration, 0 < pa < pb < 1"""
# Time, x-position, Velocity, Acceleration
t, x, v, a = 0, s_0, 0, a_0
yield (t, x, v, a)
while True:
# Roll the dices
god_wishes = random()
if god_wishes <= pa:
# Increase acceleration
a += .005
elif god_wishes <= pb:
# Reduce acceleration
a -= .005
# Lets avoid too much acceleration
a = clip(a, -.2, .2)
# How much time has passed, since last update?
dt = random()
v += dt*a
x += dt*v
t += dt
yield (t, x, v, a)
def random_walk_simul(initial_position = 0, acceleration = 0,
prob_increase=0.005, prob_decrease=0.005,
max_distance=10000, simul_time=1000):
"""Runs a random walk simulation given parameters
Particle initial state (initial position and acceleration)
State change probability (prob_increase, prob_decrease)
Stop criteria (max_distance, simul_time)
Returns a random_walk particle data
"""
assert (0 < prob_increase+prob_decrease < 1), "Total probability should be in range [0, 1]"
rw = random_walk(initial_position, acceleration, prob_increase, prob_decrease+prob_increase)
# Over estimated given by law of large numbers expected value of a
# uniform distribution
estimated_N = int(simul_time * 2.2)
data = np.empty((estimated_N, 4))
# Runs the first iteraction
n = 0
(t, x, v, a) = rw.__next__()
data[n] = (t, x, v, a)
# While there is simulation time or not too far away
while (t < simul_time) and (np.abs(x) < max_distance):
n += 1
(t, x, v, a) = rw.__next__()
data[n] = (t, x, v, a)
return data[:n]
def random_walk_naive(initial_position = 0, acceleration = 0,
prob_increase=0.005, prob_decrease=0.005,
max_distance=10000, simul_time=1000):
"""Emulates same behavior as random_walk_simul, but use append instead numpy and generators"""
T = []
X = []
V = []
A = []
T.append(0)
X.append(initial_position)
V.append(0)
A.append(acceleration)
a = A[-1]
t = T[-1]
v = V[-1]
x = X[-1]
while (T[-1] < simul_time) and (abs(X[-1]) < max_distance):
god_wishes = random()
if god_wishes <= prob_increase:
# Increase acceleration
a += .005
elif god_wishes <= prob_increase+prob_decrease:
# Reduce acceleration
a -= .005
# Lets avoid too much acceleration
a = clip(a, -.2, .2)
dt = random()
t += dt
v += dt*a
x += dt*v
# Storing next simulation step
T.append(t)
X.append(x)
V.append(v)
A.append(a)
return np.array((T, X, V, A)).transpose()
def main():
random_walk_naive(max_distance=1e9, simul_time=1e6)
random_walk_simul(max_distance=1e9, simul_time=1e6)
if __name__ == '__main__':
main()
这可能是使用numba:的好情况
import numpy as np
from random import random
from numba import njit
# Baseline
%timeit random_walk_naive(max_distance=1e9, simul_time=1e6)
1.28 s ± 277 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
# Few adjustments for numba
@njit
def random_walk_numba(initial_position = 0, acceleration = 0,
prob_increase=0.005, prob_decrease=0.005,
max_distance=10000, simul_time=1000):
T, X, V, A = [0], [initial_position], [0], [acceleration]
t, x, v, a = T[-1], X[-1], V[-1], A[-1]
while (T[-1] < simul_time) and (abs(X[-1]) < max_distance):
god_wishes = random()
if god_wishes <= prob_increase:
# Increase acceleration
a += .005
elif god_wishes <= prob_increase+prob_decrease:
# Reduce acceleration
a -= .005
# Lets avoid too much acceleration
lower, upper = -0.2, 0.2
a = lower if a < lower else upper if a > upper else a
dt = random()
t += dt
v += dt*a
x += dt*v
# Storing next simulation step
T.append(t)
X.append(x)
V.append(v)
A.append(a)
return np.array((T, X, V, A)).transpose()
%timeit random_walk_numba(max_distance=1e9, simul_time=1e6)
172 ms ± 32.2 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
请注意,您不能调用clip
,但幸运的是,它很容易在内部重新实现。