Numpy 中的矢量化 - 广播

Question

我在 python 中有一个包含以下元素的代码：

我有一个类似这样的intensities 向量：

array([ 1142.,  1192.,  1048., ...,    29.,    18.,    35.])

我还有一个x 向量，看起来像这样：

array([   0,    1,    1, ..., 1060, 1060, 1061])

然后，我使用 for 循环填充另一个向量 radialDistribution，如下所示：

for i in range(1000):
    radialDistribution[i] = sum(intensities[np.where(x == i)]) / len(np.where(x == i)[0])

问题是完成它需要 20 秒...因此我想对其进行矢量化。但是我对在 Numpy 中进行广播还是很陌生，并且在那里没有找到太多东西......因此我需要你的帮助。

我试过了，但没有用：

i= np.ogrid[:1000]
intensities[i] = sum(sortedIntensities1D[np.where(sortedDists1D == i)]) / len(np.where(sortedDists1D == i)[0])

你能帮我告诉我应该去哪里学习 Numpy 的向量化过程吗？

提前感谢您的宝贵帮助！

Answer 1

如果您的x 向量具有从 0 开始的连续整数，那么您可以这样做：

radialDistribution = np.bincount(x, weights=intensities) / np.bincount(x)

Answer 2

Here 是我在 numpy 中实现 group_by 功能。它在概念上类似于 pandas 解决方案；除了这不需要 pandas，而且在我看来应该成为 numpy 核心的一部分。

使用此功能，您的代码将如下所示：

radialDistribution = group_by(x).mean(intensities)

很快就会完成。

还要查看最后定义的 test_radial 函数，它可能更接近您的最终目标。

Answer 3

这是一个使用广播的方法：

# arrays need to be at least 2D for broadcasting
x = np.atleast_2d(x)

# create vector of indices
i = np.atleast_2d(np.arange(x.size))

# do the vectorized calculation
bool_eq = (x == i.T)
totals = np.sum(np.where(bool_eq, intensities, 0), axis=1)
rD = totals / np.sum(bool_eq, axis=1)

这使用了两次广播：在操作x == i.T 和调用np.where。不幸的是，上面的代码非常慢，甚至比原来的还要慢。这里的主要瓶颈是np.where，在这种情况下，我们可以通过获取布尔数组和强度的乘积（也可以通过广播）来加速：

totals = np.sum(bool_eq*intensities, axis=1)

这本质上与矩阵向量积相同，所以我们可以这样写：

totals = np.dot(intensities, bool_eq.T)

最终结果是比原始代码更快的代码（至少在中间数组的内存使用成为限制因素之前），但正如其他答案之一所建议的那样，使用迭代方法可能会更好.

编辑：使用np.einsum 更快（在我的试用版中）：

totals = np.einsum('ij,j', bool_eq, intensities)

Answer 4

基于我在https://stackoverflow.com/a/22265803/901925 中的itertools.groupby 解决方案，这是一个适用于2 个小型阵列的解决方案。

import numpy as np
import itertools
intensities = np.arange(12,dtype=float)
x=np.array([1,0,1,2,2,1,0,0,1,2,1,0]) # general, not sorted or consecutive

首先是一个 bincount 解决方案，针对不连续的值进行调整

# using bincount
# if 'x' are not consecutive
J=np.bincount(x)>0
print np.bincount(x,weights=intensities)[J]/np.bincount(x)[J]

现在是groupby 解决方案

# using groupby;
# sort if need
I=np.argsort(x)
x=x[I]
intensities=intensities[I]

# make a record array for use by groupby
xi=np.zeros(shape=x.shape, dtype=[('intensities',float),('x',int)])
xi['intensities']=intensities
xi['x']=x

g=itertools.groupby(xi, lambda z:z['x'])
xx=np.array([np.array([z[0] for z in y[1]]).mean() for y in g])
print xx

这是一个紧凑的numpy 解决方案，使用np.unique 和np.split 的return_index 选项。 x 应该被排序。我对大型数组的速度并不乐观，因为除了理解之外，unique 和 split 还会有迭代。

[values, index] = np.unique(x, return_index=True)
[y.mean() for y in np.split(intensities, index[1:])]