[Numpy-discussion] broacasting question

[Thomas K Gamble]

> On 30.06.2011, at 11:57PM, Thomas K Gamble wrote:
> >> np.add(b.reshape(2048,3136) * c, d, out=a[:,:3136])
> >> 
> >> But to say whether this is really the equivalent result to what IDL
> >> does, one would have to study the IDL manual in detail or directly
> >> compare the output (e.g. check what happens to the values in
> >> a[:,3136:]...)
> >> 
> >> Cheers,
> >> 
> >> 							Derek
> > 
> > Your post gave me the cluse I needed.
> > 
> > I had my shapes slightly off in the example I gave, but if I try:
> > 
> > a = reshape(b.flatten('F') * c.flatten('F') + d.flatten('F'), b.shape,
> > order='F')
> > 
> > I get a result in line with the IDL result.
> > 
> > Another example with different total size arrays:
> > 
> > b = np.ndarray((2048,3577), dtype=float)
> > c = np.ndarray((256,25088), dtype=float)
> > 
> > a= reshape(b.flatten('F')[:size(c)]/c.flatten('F'), c.shape, order='F')
> > 
> > This also gives a result like that of IDL.
> 
> Right, I forgot to point out that there are at least 2 ways to bring the
> arrays into compatible shapes (that's the reason broadcasting does not
> work here, because numpy only does automatic broadcasting if there is an
> unambiguous way to do so). So the IDL arrays being Fortran-ordered is the
> essential bit of information here. Just two remarks:
> I. Assigning a = reshape(b.flatten('F')[:size(c)]/c.flatten('F'), c.shape,
> order='F') as above will create a new array of shape c.shape - if you
> wanted to put your results into an existing array of shape(2048,3577),
> you'd still have to explicitly say a[:,:3136] = ... II. The flatten()

That was the error in my first example I refered to.  I confused it with the 
second 'divide' example and probably should have used different variable names 
to avoid confusing things further.  Sorry.

> operations and the assignment above all create full copies of the arrays,
> thus the np.add ufunc above together with simple reshape operations might
> improve performance somewhat - however keeping the Fortran order also
> requires some costly transpositions, as for your last example

Right now, I'm just interested in getting the right answer.  Once I have that, 
I'll work on performance.  Unfortunately the order does seem to make a 
difference.

> 
> a = np.divide(b.T[:3136].reshape(c.T.shape).T, c, out=a)

Interesting.  

> 
> so YMMV...
> 
> Cheers,
> 						Derek

Thanks for your help.

--
Thomas K. Gamble tkgamble at windstream.net
LANL employee waiting out the Las Conchas fire.