Elipsis emitiendo en numpy.einsum

Tengo problemas para entender por qué lo siguiente no funciona:

Tengo un prefactor de matriz que puede ser tridimensional o de seis dimensiones. Tengo una serie de dipolos que tiene cuatro dimensiones. Las tres primeras dimensiones de los dipolos coinciden con las tres últimas dimensiones del prefactor .

Como no conozco la forma del prefactor , estoy usando una elipsis para tener en cuenta las tres dimensiones opcionales en el prefactor :

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

(En el ejemplo aquí, prefactor.shape es (1, 1, 1, 160, 160, 128) y dipoles.shape es (160, 160, 128, 3). Al ejecutarlo, aparece el error:

El operando 1 no tenía suficientes dimensiones para coincidir con la transmisión, y no pudo extenderse porque los subíndices de sum de einstein se especificaron tanto al inicio como al final

Sin embargo, funciona cuando agrego puntos suspensivos al segundo término también:

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

Solo que no entiendo por qué, porque no debería haber necesidad de una elipsis allí. ¿Alguien sabe lo que está pasando aquí?

La misma pregunta se ha formulado en http://comments.gmane.org/gmane.comp.python.numeric.general/53705 pero todavía no hay una respuesta satisfactoria.

Hay un problema con github para este problema:

https://github.com/numpy/numpy/issues/2455 mejora de la notación de índice en einsum (Trac # 1862)

Caso de error:

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

El trabajo actual requiere (vacío) puntos suspensivos:

einsum (‘ij …, j …-> ij …’, A, B)

Parece que einsum recorre el argumento de cadena y las operaciones varias veces, identificando los índices y los tipos de transmisión (derecha, izquierda, centro, ninguno) y las dimensiones op. Con esto construye un numpy.nditer . Es mientras se construyen op_axes para el nditer que einsum genera este error. No sé si los criterios de prueba son demasiado ajustados ( ibroadcast >= ndim ), o si es necesario dar un paso adicional para construir los op_axes correctos para este argumento.

https://github.com/numpy/numpy/issues/2619 muestra cómo se puede usar einsum para replicar el comportamiento de einsum . Trabajando de esto puedo replicar tu cálculo así:

 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) 

La línea op_axes es indicativa de lo que einsum deduce de '...lmn,...lmno->...o' .


Estoy explorando este problema en https://github.com/hpaulj/numpy-einsum .

Ahí tengo un einsum_py.py que emula np.einsum con el código Python. La parte que es relevante para este problema es parse_subscripts() , y en particular prepare_op_axes() . Parece que solo se necesita la iteración BROADCAST_RIGHT (comenzando desde el final) para crear correctamente op_axes , independientemente de dónde estén los puntos suspensivos en los subíndices. También elimina el mensaje de error que se encuentra en el núcleo de este problema.

El archivo einsum.c.src en ese repository tiene este cambio y se comstack correctamente con la distribución maestra actual (solo reemplace el archivo y compile). Prueba bien contra test_einsum.py , así como ejemplos de este problema.

He enviado una solicitud de extracción para este cambio.