我觉得应该有一种快速的方法来加速这段代码。我想答案是here ,但我似乎无法以那种格式解决我的问题。我试图解决的根本问题是找到平行和垂直分量的逐点差异,并创建这些差异的二维直方图。
out = np.zeros((len(rpbins)-1,len(pibins)-1))
tmp = np.zeros((len(x),2))
for i in xrange(len(x)):
tmp[:,0] = x - x[i]
tmp[:,1] = y - y[i]
para = np.sum(tmp**2,axis=-1)**(1./2)
perp = np.abs(z - z[i])
H, _, _ = np.histogram2d(para, perp, bins=[rpbins, pibins])
out += H
最佳答案
向量化这样的事情很棘手,因为要摆脱 n
元素的循环,你必须构造一个 (n, n)
数组,所以对于大输入您可能会获得比 Python 循环更差的性能。但这是可以做到的:
mask = np.triu_indices(x.shape[0], 1)
para = np.sqrt((x[:, None] - x)**2 + (y[:, None] - y)**2)
perp = np.abs(z[:, None] - z)
hist, _, _ = np.histogram2d(para[mask], perp[mask], bins=[rpbins, pibins])
mask
是为了避免对每个距离计算两次。我还将对角线偏移设置为 1
,以避免在直方图中包含每个点到自身的 0
距离。但是,如果您不使用它索引 para
和 perp
,您将获得与您的代码完全相同的结果。
使用此示例数据:
items = 100
rpbins, pibins = np.linspace(0, 1, 3), np.linspace(0, 1, 3)
x = np.random.rand(items)
y = np.random.rand(items)
z = np.random.rand(items)
我为我的hist
和你的out
得到了这个:
>>> hist
array([[ 1795., 651.],
[ 1632., 740.]])
>>> out
array([[ 3690., 1302.],
[ 3264., 1480.]])
和out[i, j] = 2 * hist[i, j]
除了i = j = 0
,其中out[0, 0 ] = 2 * hist[0, 0] + items
因为每个项目与其自身的距离为 0
。
编辑 在 tcaswell 发表评论后尝试了以下操作:
items = 1000
rpbins, pibins = np.linspace(0, 1, 3), np.linspace(0, 1, 3)
x, y, z = np.random.rand(3, items)
def hist1(x, y, z, rpbins, pibins) :
mask = np.triu_indices(x.shape[0], 1)
para = np.sqrt((x[:, None] - x)**2 + (y[:, None] - y)**2)
perp = np.abs(z[:, None] - z)
hist, _, _ = np.histogram2d(para[mask], perp[mask], bins=[rpbins, pibins])
return hist
def hist2(x, y, z, rpbins, pibins) :
mask = np.triu_indices(x.shape[0], 1)
para = np.sqrt((x[:, None] - x)[mask]**2 + (y[:, None] - y)[mask]**2)
perp = np.abs((z[:, None] - z)[mask])
hist, _, _ = np.histogram2d(para, perp, bins=[rpbins, pibins])
return hist
def hist3(x, y, z, rpbins, pibins) :
mask = np.triu_indices(x.shape[0], 1)
para = np.sqrt(((x[:, None] - x)**2 + (y[:, None] - y)**2)[mask])
perp = np.abs((z[:, None] - z)[mask])
hist, _, _ = np.histogram2d(para, perp, bins=[rpbins, pibins])
return hist
In [10]: %timeit -n1 -r10 hist1(x, y, z, rpbins, pibins)
1 loops, best of 10: 289 ms per loop
In [11]: %timeit -n1 -r10 hist2(x, y, z, rpbins, pibins)
1 loops, best of 10: 294 ms per loop
In [12]: %timeit -n1 -r10 hist3(x, y, z, rpbins, pibins)
1 loops, best of 10: 278 ms per loop
似乎大部分时间都花在实例化新数组上,而不是进行实际计算,因此虽然有一些效率可以刮掉,但实际上并不多。
关于python - 如何在 numpy 中向量化这个循环差异?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/14951676/