如何在numpy中进行分散/聚集操作

Question

假设我有数组：

a = array((1,2,3,4,5))
indices = array((1,1,1,1))

我做手术：

^{pr2}$

结果是

array([1, 3, 3, 4, 5])

换句话说，indices中的重复项被忽略

如果我不想忽略重复项，结果是：

array([1, 6, 3, 4, 5])

我该怎么做？在

上面的例子有些琐碎，下面就是我要做的：

def inflate(self,pressure):
    faceforces = pressure * cross(self.verts[self.faces[:,1]]-self.verts[self.faces[:,0]], self.verts[self.faces[:,2]]-self.verts[self.faces[:,0]])
    self.verts[self.faces[:,0]] += faceforces
    self.verts[self.faces[:,1]] += faceforces
    self.verts[self.faces[:,2]] += faceforces

def constrain_lengths(self):
    vectors = self.verts[self.constraints[:,1]] - self.verts[self.constraints[:,0]]
    lengths = sqrt(sum(square(vectors), axis=1))
    correction = 0.5 * (vectors.T * (1 - (self.restlengths / lengths))).T
    self.verts[self.constraints[:,0]] += correction
    self.verts[self.constraints[:,1]] -= correction

def compute_normals(self):
    facenormals = cross(self.verts[self.faces[:,1]]-self.verts[self.faces[:,0]], self.verts[self.faces[:,2]]-self.verts[self.faces[:,0]])
    self.normals.fill(0)
    self.normals[self.faces[:,0]] += facenormals
    self.normals[self.faces[:,1]] += facenormals
    self.normals[self.faces[:,2]] += facenormals
    lengths = sqrt(sum(square(self.normals), axis=1))
    self.normals = (self.normals.T / lengths).T

由于在索引赋值操作中忽略了重复项，我得到了一些非常错误的结果。在

Answer 1

我不知道有什么方法比：

for face in self.faces[:,0]:
    self.verts[face] += faceforces

你也可以自我面孔放入一个由3个字典组成的数组中，其中的键对应于面，值对应于需要添加的次数。你会得到这样的代码：

^{pr2}$

可能更快。我真的希望有人能想出一个更好的方法，因为我想在今天早些时候帮助别人加速他们的代码。在

Answer 2

numpy的histogram函数是一个分散操作。在

a += histogram(indices, bins=a.size, range=(0, a.size))[0]

您可能需要小心，因为如果indices包含整数，那么小的舍入错误可能会导致值在错误的bucket中结束。在这种情况下使用：

a += histogram(indices, bins=a.size, range=(-0.5, a.size-0.5))[0]

把每个索引放到每个箱子的中心。在

更新：这是有效的。但我建议使用@Eelco Hoogendoorn基于numpy.add.at的答案。在

Answer 3

派对稍微晚了一点，但看到这种手术是多么普遍，而且它似乎仍然不是标准纽姆的一部分，我把我的解决方案放在这里作为参考：

def scatter(rowidx, vals, target):
    """compute target[rowidx] += vals, allowing for repeated values in rowidx"""
    rowidx = np.ravel(rowidx)
    vals   = np.ravel(vals)
    cols   = len(vals)
    data   = np.ones(cols)
    colidx = np.arange(cols)
    rows   = len(target)
    from scipy.sparse import coo_matrix
    M = coo_matrix((data,(rowidx,colidx)), shape=(rows, cols))
    target += M*vals
def gather(idx, vals):
    """for symmetry with scatter"""
    return vals[idx]

对于初学者来说，numpy中的自定义C例程的速度可能会快一倍，消除了多余的1的分配和乘法，但与python中的循环相比，它在性能上有着天壤之别。在

除了性能方面的考虑，它在风格上更符合其他numpy矢量化代码，使用分散操作，而不是在代码中混搭一些for循环。在

编辑：

好吧，忘了上面的事。从最新的1.8版本开始，numpy以最佳效率直接支持分散操作。在

^{pr2}$