2015-02-13 13:25 GMT+01:00 Julian Taylor <jtaylor.debian@googlemail.com>:
On 02/13/2015 01:03 PM, Francesc Alted wrote:
> 2015-02-13 12:51 GMT+01:00 Julian Taylor <jtaylor.debian@googlemail.com
> <mailto:jtaylor.debian@googlemail.com>>:
>
>     On 02/13/2015 11:51 AM, Francesc Alted wrote:
>     > Hi,
>     >
>     > I would like to vectorize the next computation:
>     >
>     > nx, ny, nz = 720, 180, 3
>     > outheight = np.arange(nz) * 3
>     > oro = np.arange(nx * ny).reshape((nx, ny))
>     >
>     > def compute1(outheight, oro):
>     >     result = np.zeros((nx, ny, nz))
>     >     for ix in range(nx):
>     >         for iz in range(nz):
>     >             result[ix, :, iz] = outheight[iz] + oro[ix, :]
>     >     return result
>     >
>     > I think this should be possible by using an advanced use of
>     broadcasting
>     > in numpy.  Anyone willing to post a solution?
>
>
>     result = outheight + oro.reshape(nx, ny, 1)
>
>
> And 4x faster for my case.  Oh my, I am afraid that my mind will never
> scratch all the amazing possibilities that broadcasting is offering :)
>
> Thank you very much for such an elegant solution!
>


if speed is a concern this is faster as it has a better data layout for
numpy during the computation, but the result may be worse layed out for
further processing

    result = outheight.reshape(nz, 1, 1) + oro
    return np.rollaxis(result, 0, 3)


Holly cow, this makes for another 4x speed improvement!  I don't think I need that much in my scenario, so I will stick with the first one (more readable and the expected data layout), but thanks a lot!

Francesc