python - 在 python 中加速 monte carlo 代码

Question

考虑从 [0,T) 开始按递增顺序给出的点 Y。我们要将这些点视为位于圆周 T 的圆上。现在考虑来自 [0,T) 的点 X 也位于圆周 T 的圆上。

我们说 X 和 Y 之间的距离是 X 中的每个点与其在 Y 中最近的点之间的绝对距离之和，这两个点都被认为位于一个圆中。将此距离写为 Delta(X, Y)。

我正在尝试找到一种快速方法来近似计算 X 的所有可能旋转的圆之间的距离分布。我目前正在通过蒙特卡罗模拟来完成这项工作。首先是我制作一些假数据的代码。

import random
import numpy as np
from bisect import bisect_left

def simul(rate, T):
    time = np.random.exponential(rate)
    times = [0]
    newtime = times[-1]+time
    while (newtime < T):
        times.append(newtime)
        newtime = newtime+np.random.exponential(rate)
    return times[1:]

现在的代码是求两个圆之间的距离。

def takeClosest(myList, myNumber, T):
    """
    Assumes myList is sorted. Returns closest value to myNumber in a circle of circumference T.

    If two numbers are equally close, return the smallest number.
    """
    pos = bisect_left(myList, myNumber)
    if (pos == 0 and myList[pos] != myNumber):
        before = myList[pos - 1] - T
        after = myList[0]
    elif (pos == len(myList)):
        before = myList[pos-1]
        after = myList[0] + T
    else:
        before = myList[pos - 1]
        after = myList[pos]
    if after - myNumber < myNumber - before:
        return after
    else:
        return before

def circle_dist(timesY, timesX):
    dist = 0
    for t in timesX:
        closest_number = takeClosest(timesY, t, T)
        dist += np.abs(closest_number - t)
    return dist

现在是制作数据和尝试 1000 次不同随机旋转的主要代码。

T = 50000
timesX = simul(1, T)
timesY = simul(10, T)

dists=[]
iters = 100
for i in xrange(iters):
       offset = np.random.randint(0,T)
       timesX = [(t+offset) % T for t in timesX]
       dists.append(circle_dist(timesY, timesX))

我们现在可以打印出任何我们喜欢的距离统计数据。我对方差特别感兴趣。

print "Variance is ", np.var(dists)

不幸的是，我经常需要这样做，目前大约需要 16 秒。我觉得这有点令人惊讶，它是如此之慢。非常感谢收到有关如何加快速度的任何建议。

---

编辑 1. 将迭代次数减少到 100(之前的值与我的计时不正确)。这在我的电脑上大约需要 16 秒。

编辑 2. 修复了 takeClosest 中的错误

Answer 1

编辑:我刚刚注意到性能优化有点过早，因为表达式 closest_number - t 不是“圆”上 distance 的任何定义的有效实现 - 这只是一个开放式直线上的距离

示例测试用例(伪代码):

T = 10
X = [1, 2]
Y = [9]
dist(X, Y) = dist(1, 9) + dist(2, 9)
dist_on_line   = 8 + 7 = 15
dist_on_circle = 2 + 3 = 5

注意圆圈[0,10)的定义意味着dist(0, 10)没有定义，但在极限上它趋近于0:lim(dist(0, t), t->10) = 0

圆上距离的正确实现是:

dist_of_t = min(t - closest_number_before_t,
                closes_number_after_t - t,
                T - t + closes_number_before_t,
                T - closest_number_after_t + t)

---

原答案:

您可以旋转和迭代 timesY 而不是 timesX，因为该数组小一个数量级 - 与迭代所有元素 (bisect_left) 相比，执行 timeX 的 O(logn) 可以忽略不计 (O(n))

但是恕我直言，真正的减速是因为 Python 动态类型(每次尝试将 timesX 中的 ~50000 项与其他值进行比较时都必须检查类型兼容性)=> 将 timesX 和 timesY 转换为numpy 数组应该有帮助，如果这还不够，CPU 加速(cython、numba、...)是您认为需要的