cómo hacer de dispersión/agrupación operaciones en numpy

Question

Digamos que tengo matrices:

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

y llevo a cabo la operación:

a[indices] += 1 

el resultado es

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

en otras palabras, , los duplicados en indices se ignoran

si quería los duplicados no debe ser ignorado, lo que resulta en:

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

¿cómo iba a ir sobre esto?

el ejemplo anterior es algo trivial, lo que sigue es exactamente lo que estoy tratando de hacer:

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 

he estado consiguiendo unos resultados muy con errores como resultado de duplicados ser ignorado en mis operaciones de asignación indexados.

Answer 1

no sé de una manera de hacerlo que es más rápido que:

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

También puede hacer self.faces en una matriz de 3 diccionarios donde las claves se corresponden con la cara y el valor al número de veces que necesita ser agregado. Luego obtendría código como:

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

que podría ser más rápido. Espero que a alguien se le ocurra una mejor manera porque quería hacer esto cuando trato de ayudar a alguien a acelerar su código el día de hoy.

Answer 2

numpy La función histogram es una operación de dispersión.

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

Es posible que tenga que tomar un cierto cuidado porque si indices contiene números enteros, los pequeños errores de redondeo puede resultar en valores de acabar en el cubo mal. En el que caso de uso:

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

para obtener cada índice en el centro de cada bin.

Actualización: esto funciona. Pero recomiendo usar la respuesta de @Eelco Hoogendoorn basada en numpy.add.at.

Answer 3

poco tarde a la fiesta, pero al ver la forma comúnmente se requiere esta operación, y el hecho de que todavía no parece ser parte de numpy estándar, mal poner mi solución aquí por referencia:

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] 

Una rutina C personalizada en numpy fácilmente podría ser dos veces más rápida, eliminando la asignación superflua y la multiplicación con unos, para empezar, pero hace una gran diferencia en el rendimiento frente a un bucle en python.

Aparte de las consideraciones de rendimiento, es estilísticamente mucho más en línea con otros códigos numpy-vectorized para usar una operación de dispersión, en lugar de combinar algunos bucles para su código.

Edit:

Ok, olvídate de lo anterior. A partir de la versión 1.8 más reciente, realizar operaciones de dispersión ahora se admite directamente en numpy con una eficiencia óptima.

def scatter(idx, vals, target): 
    """target[idx] += vals, but allowing for repeats in idx""" 
    np.add.at(target, idx.ravel(), vals.ravel())