加速python中的numpy循环？

Question

考虑以下使用numpy数组的代码，它非常慢：

# Intersection of an octree and a trajectory
def intersection(octree, trajectory):
    # Initialize numpy arrays
    ox = octree.get("x")
    oy = octree.get("y")
    oz = octree.get("z")
    oe = octree.get("extent")/2
    tx = trajectory.get("x")
    ty = trajectory.get("y")
    tz = trajectory.get("z")
    result = np.zeros(np.size(ox))
    # Loop over elements
    for i in range(0, np.size(tx)):
        for j in range(0, np.size(ox)):
            if (tx[i] > ox[j]-oe[j] and 
                tx[i] < ox[j]+oe[j] and 
                ty[i] > oy[j]-oe[j] and 
                ty[i] < oy[j]+oe[j] and 
                tz[i] > oz[j]-oe[j] and 
                tz[i] < oz[j]+oe[j]):
                result[j] += 1
    # Finalize
    return result

如何重写函数以加快计算速度？（np.size(tx) == 10000和np.size(ox) == 100000）

Answer 1

你分配了10000张100000大小的列表。要做的第一件事是停止对嵌套的j循环使用range，而使用生成器版本xrange。这将节省您分配所有这些列表的时间和空间。在

下一个方法是使用矢量化操作：

for i in xrange(0, np.size(tx)):
    index = (ox-oe < tx[i]) & (ox+oe > tx[i]) & (oy-oe < ty[i]) & (oy+oe > ty[i]) & (oz-oe < tz[i]) & (oz+oe > tz[i])
    result[index] += 1  

Answer 2

通过在PyPy:http://pypy.org/（在https://bitbucket.org/pypy/numpy上对NumPy集成的说明）下运行这段代码，您可能会得到很好的结果

Answer 3

我认为这将给双循环同样的结果，并且更快：

for j in xrange(np.size(ox)):
    result[j] += sum( abs(tx-ox[j])<oe[j] & abs(ty-oy[j])<oe[j] & abs(tz-oz[j])<oe[j] )

要实现这一点：1）重新排序循环（即交换它们），因为循环内没有任何更改；2）将result[j]拉到i循环之外；3）将所有t>ox-oe and t<ox+oe转换为abs(t-ox)<oe（虽然这可能不是一个巨大的加速，但更容易阅读）。在

既然你没有可运行的代码，我也不想为此建立一个测试，我不能百分之百地确定这是正确的。在