Многоточие вещания в numpy.einsum

Question

У меня проблема понять, почему следующий не работает:

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

Как я не знаю формы предэкспонентов, я использую Ellipsis для учета трех дополнительных измерений в предэкспонентах:

numpy.einsum('...lmn,lmno->...o', prefactor, dipoles) 

(В приведенном примере, префактор. (1, 1, 1, 160, 160, 128) и диполей. Форма (160, 160, 128, 3). При выполнении я получаю сообщение об ошибке:

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

Это работа, однако, когда я добавить многоточие на второй срок, а также:

numpy.einsum('...lmn,...lmno->...o', prefactor, dipoles) 

Просто, что я не понимаю, почему, потому что должен не нужно использовать для этого многоточие. Кто-нибудь знает, что здесь происходит?

Тот же вопрос был задан в http://comments.gmane.org/gmane.comp.python.numeric.general/53705, но пока нет удовлетворительного ответа.

Answer 1

Существует проблема GitHub для решения этой проблемы:

https://github.com/numpy/numpy/issues/2455 улучшение индекса записи в einsum (TraC# 1862)

Ошибка случай:

einsum('ij...,j->ij...',A,B) 

Текущая работа вокруг требует (пусто) эллипсис:

einsum ('ij ..., j ...-> ij ...', A, B)

Похоже, что einsum циклически перебирает аргумент string и ops, идентифицируя индексы и типы широковещания (справа, слева, в середине, ни один) и размеры op. При этом он создает numpy.nditer. При построении op_axes для nditer einsum вызывает эту ошибку. Я не знаю, слишком ли сложны критерии теста (ibroadcast >= ndim), или если для этого аргумента необходимо сделать дополнительный шаг для построения правильного op_axes.

https://github.com/numpy/numpy/issues/2619 показывает, как nditer может использоваться для репликации поведения einsum. Исходя из этого я могу повторить свой расчет следующим образом:

prefactor = np.random.random((1, 1, 1, 160, 160, 128)) 
dipoles = np.random.random((160, 160, 128, 3)) 
x = numpy.einsum('...lmn,...lmno->...o', prefactor, dipoles) 
#numpy.einsum('...lmn,lmno->...o', prefactor, dipoles) # not work 

op_axes = [[0,1,2,3,4,5,-1], [-1,-1,-1,0,1,2,3], [0,1,2,-1,-1,-1,3]] 
flags = ['reduce_ok','buffered', 'external_loop', 'delay_bufalloc', 'grow_inner'] 
op_flags = [['readonly']]*nops + [['allocate','readwrite']] 
it = np.nditer([prefactor,dipoles,None], flags, op_flags, op_axes=op_axes) 
it.operands[nops][...] = 0 
it.reset() 
#it.debug_print() 
for (x,y,w) in it: 
    w[...] += x*y 
print "\nnditer usage:" 
print it.operands[nops] # == x 
print it.operands[nops].shape # (1, 1, 1, 3) 

op_axes линия указывает на то, что einsum выводит из '...lmn,...lmno->...o'.

---

Я изучаю этот вопрос на https://github.com/hpaulj/numpy-einsum.

Там у меня есть einsum_py.py, который эмулирует np.einsum с кодом Python. Часть, относящаяся к этой проблеме, составляет parse_subscripts(), и в частности prepare_op_axes(). Похоже, что для корректного создания op_axes требуется только итерация BROADCAST_RIGHT, независимо от того, где находятся эллипсы в нижних индексах. Он также удаляет сообщение об ошибке, лежащее в основе этой проблемы.

Файл einsum.c.src в этом репозитории имеет это изменение и правильно компилируется с текущим основным дистрибутивом (просто замените файл и его сборку). Он отлично тестирует test_einsum.py, а также примеры из этой проблемы.

Я представил запрос на перенос для этого изменения.