[NumPy-Tickets] [NumPy] #2027: einsum/iterator bug

[NumPy Trac]

#2027: einsum/iterator bug
------------------------+---------------------------------------------------
 Reporter:  mwiebe      |       Owner:  somebody
     Type:  defect      |      Status:  new     
 Priority:  high        |   Milestone:  1.7.0   
Component:  numpy.core  |     Version:  1.6.1   
 Keywords:              |  
------------------------+---------------------------------------------------

Comment(by mwiebe):

 Here's some code which replicates the nditer einsum uses. The output from
 the debug_print here matches the output of a NpyIter_DebugPrint inserted
 in the einsum C code (einsum.c.src line 3083) in all the salient details.
 This suggests it may not be an iterator bug, but rather something in the
 einsum looping code.

 {{{
 import numpy;
 import numpy as np;
 print(numpy.version.full_version);
 a = numpy.arange(1, 3);
 b = numpy.arange(1, 5).reshape(2, 2);
 c = numpy.arange(1, 9).reshape(4, 2);

 print "\nground truth:"
 print numpy.array([[[a[x]*b[y, x]*c[z, x] for y in range(2)] for z in
 range(4)] for x in range(2)])

 it = np.nditer([a,b,c,None],
             ['reduce_ok','buffered', 'external_loop',
              'delay_bufalloc', 'grow_inner'],
             [['readonly']]*3 + [['allocate','readwrite']],
             op_axes=[[0,-1,-1],[1,-1,0],[1,0,-1],[0,1,2]])
 it.operands[3][...] = 0
 it.reset()
 it.debug_print()
 for (x,y,z,w) in it:
     w[...] += x*y*z
 print "\nnditer usage:"
 print it.operands[3]

 print "\neinsum:"
 print numpy.einsum('x,yx,zx->xzy', a, b, c)
 }}}

-- 
Ticket URL: <http://projects.scipy.org/numpy/ticket/2027#comment:1>
NumPy <http://projects.scipy.org/numpy>
My example project