[Numpy-discussion] broacasting question

Derek Homeier derek@astro.physik.uni-goettingen...
Thu Jun 30 18:31:31 CDT 2011

```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() 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

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

so YMMV...

Cheers,
Derek

```