numpy : efficient sum computations

[Robert Kern]

TG wrote:
> Okay, another one which I don't have answer for.
> 
> it is the reverse case, sort of :
> 
>> phi.shape (x,y)
>> d.shape (a,b)
> 
> I want to return m :
>> m.shape = (x,y,a,b)
> with
> m[x,y] = d * phi[x,y]
> 
> currently, my code is :
>> m = empty(phi.shape + d.shape)
>> m[:,:] = d
> 
> this repeats the matrix d x*y times, which is what I want. But now, if
> I do :
>> m[:,:] = d * phi[:,:]
> in order to multiply each "d" by a specific value of "phi", it doesn't
> work.
> 
> <type 'exceptions.ValueError'>: shape mismatch: objects cannot be
> broadcast to a single shape
> 
> Right now I'm stuck here. If anyone can help, I would be glad.

Use broadcasting to get the desired behavior. It's a widely used feature of
numpy that takes some getting used to, but it is quite powerful. Here is some
documentation:

  http://www.scipy.org/EricsBroadcastingDoc

It still references Numeric, numpy's predecessor, but the concepts are still
valid. Use the numpy.newaxis object to add new axes to phi and d in the right
places, then multiply them together:


In [8]: from numpy import arange, newaxis

In [9]: phi = arange(5*4).reshape((5,4))

In [10]: phi.shape
Out[10]: (5, 4)

In [11]: d = arange(3*2).reshape((3,2))

In [12]: d.shape
Out[12]: (3, 2)

In [13]: m = phi[:,:,newaxis,newaxis] * d[newaxis,newaxis,:,:]

In [14]: m.shape
Out[14]: (5, 4, 3, 2)

In [17]: m[1, 2]
Out[17]:
array([[ 0,  6],
       [12, 18],
       [24, 30]])

In [18]: d * phi[1, 2]
Out[18]:
array([[ 0,  6],
       [12, 18],
       [24, 30]])

In [19]: for i in range(phi.shape[0]):
   ....:     for j in range(phi.shape[1]):
   ....:         assert (m[i,j] == d*phi[i,j]).all()
   ....:
   ....:

In [20]:

-- 
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless enigma
 that is made terrible by our own mad attempt to interpret it as though it had
 an underlying truth."
  -- Umberto Eco