Обновление узловых значений без цикла for, используя numpy

Question

Я пытаюсь обновить узловые значения в сетке из значений элементов.

В массиве faces я определяю идентификаторы узлов элемента (предположим, что Iтолько два элемента):

faces = np.array([[0, 1, 2], [1, 3, 2]])

Массив force_el содержит, скажем, силу, действующую одинаково на каждый узел элемента:

force_el = np.array([[0.7, 1.1], [1.2, 0.3]])

Теперь я хотел быобновите узловые силы force_node:

force_node = np.zeros((4, force_el.shape[1]))
for face, fel in zip(faces, force_el):
    force_node[face.ravel(), :] += fel

, чтобы получить результат:

>>> force_node
array([[0.7, 1.1],
       [1.9, 1.4],
       [1.9, 1.4],
       [1.2, 0.3]])

Поскольку это обновление должно выполняться много раз (порядка 100k-1m раз), Я пытаюсь оптимизировать его, но не вижу хорошего решения.

Answer 1

Вы можете использовать некоторые matrix-multiplication force -

out_nrows = 4 # number of nodes
mask = np.zeros((len(faces),out_nrows),dtype=bool)
np.put_along_axis(mask,faces,True,axis=1)
force_node_out = mask.T.dot(force_el)

С небольшим количеством столбцов в force_el мы также можем использовать np.bincount для еще лучшей производительности -

out_nrows = 4 # number of nodes
out = np.zeros((out_nrows, force_el.shape[1]))
n = faces.shape[1]
l = force_el.shape[1]
for i in range(n):
    for j in range(l):
        out[:,j] += np.bincount(faces[:,i],force_el[:,j],minlength=out_nrows)

Сроки -

In [35]: # Setup data (from OP's comments)
    ...: np.random.seed(0)
    ...: faces=np.array([np.random.choice(1800,3,replace=0) for i in range(3500)])
    ...: force_el = np.random.rand(len(faces),3)

In [36]: %%timeit # Original loopy soln
    ...: out_nrows = 1800
    ...: force_node = np.zeros((out_nrows, force_el.shape[1]))
    ...: for face, fel in zip(faces, force_el):
    ...:     force_node[face.ravel(), :] += fel
100 loops, best of 3: 16.1 ms per loop

In [37]: %%timeit # @RafaelC's soln with np.add.at
    ...: force_node = np.zeros((1800, force_el.shape[1]))
    ...: np.add.at(force_node, faces, force_el[:,None])
100 loops, best of 3: 2.45 ms per loop

In [38]: %%timeit # Posted in this post that uses matrix-multiplication
    ...: out_nrows = 1800
    ...: mask = np.zeros((len(faces),out_nrows),dtype=bool)
    ...: np.put_along_axis(mask,faces,True,axis=1)
    ...: force_node_out = mask.T.dot(force_el)
10 loops, best of 3: 38.4 ms per loop

In [39]: %%timeit # Posted in this post that uses bincount
    ...: out_nrows = 1800
    ...: out = np.zeros((out_nrows, force_el.shape[1]))
    ...: n = faces.shape[1]
    ...: l = force_el.shape[1]
    ...: for i in range(n):
    ...:     for j in range(l):
    ...:         out[:,j]+=np.bincount(faces[:,i],force_el[:,j],minlength=out_nrows)
10000 loops, best of 3: 149 µs per loop

Answer 2

Вы должны использовать преимущества numpy вещания, когда это возможно.

Использование np.add.at

np.add.at(force_node, faces, force_el[:,None])