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
_______________________________________________
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
