как сделать разброс / собрать операции в NumPy

Question

допустим, у меня есть массивы:

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

и я выполняю операцию:

a[indices] += 1

результат

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

numpy Функция histogram - это операция разброса.

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

Возможно, вам придется позаботиться, потому что, если indices содержит целые числа, небольшие ошибки округленияможет привести к тому, что значения окажутся в неправильном ведре.В этом случае используйте:

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

, чтобы поместить каждый указатель в центр каждой ячейки.

Обновление: это работает.Но я рекомендую использовать ответ @Eelco Hoogendoorn, основанный на numpy.add.at.

Answer 2

Немного опоздал на вечеринку, но, видя, как часто эта операция требуется, и тот факт, что она все еще не является частью стандартного numpy, я поставлю здесь мое решение для справки:

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]

Пользовательская подпрограмма C в numpy может быть в два раза быстрее, исключая излишнее распределение и умножение на единицу, для начинающих, но она сильно отличается от производительности по сравнению с циклом в python.

Помимо соображений производительности, стилистически гораздо больше соответствует другой код с векторизованным вектором для использования операции разброса, а не для перемешивания некоторых циклов в вашем коде.

Редактировать:

Хорошо,забудь о вышесказанном.Начиная с последней версии 1.8, операции разброса теперь напрямую поддерживаются в numpy с оптимальной эффективностью.

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

Answer 3

Я не знаю, как сделать это быстрее, чем:

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

Вы также можете превратить self.faces в массив из 3 словарей, где ключи соответствуют лицу изначение к числу раз, которое это должно быть добавлено.Затем вы получите код:

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

, который может быть быстрее.Я действительно надеюсь, что кто-то придумает лучший способ, потому что я хотел сделать это, пытаясь помочь кому-то ускорить его код сегодня.