# [SciPy-user] filling array without loop...

Anne Archibald peridot.faceted@gmail....
Sun Apr 22 12:23:42 CDT 2007

On 22/04/07, fred <fredmfp@gmail.com> wrote:

> Works fine, but a bit too slow for nx,ny,nz \approx 10^2 (endless).

This probably gives you an array that needs some transpose()ing before
you ravel it, but you can more or less write your formula as is:

i = arange(nx)[newaxis,newaxis,:]
j = arange(ny)[newaxis,:,newaxis]
k = arange(nz)[:,newaxis,newaxis]

Now any elementwise arithmetic you do on i j and k will broadcast them
up to arrays of the right shape, so for example
k[:-1,:,:]*(nx*ny)+j[:,:-1,:]*nx+i[:,:,:-1]
yields a cubical matrix that is essentially the first formula you
gave. Putting them all together (I'm sure there's a tidier way to do
this, particularly if you don't care much about the order):
ca1 = array([VTK_HEXAHEDRON_NB_POINTS*ones((nz-1,ny-1,nx-1),dtype=int),
k[:-1,:,:]*(nx*ny)+j[:,:-1,:]*nx+i[:,:,:-1],
k[:-1,:,:]*(nx*ny)+j[:,:-1,:]*nx+i[:,:,1:],
k[:-1,:,:]*(nx*ny)+j[:,1:,:]*nx+i[:,:,1:],
k[:-1,:,:]*(nx*ny)+j[:,1:,:]*nx+i[:,:,:-1],
k[1:,:,:]*(nx*ny)+j[:,:-1,:]*nx+i[:,:,:-1],
k[1:,:,:]*(nx*ny)+j[:,:-1,:]*nx+i[:,:,1:],
k[1:,:,:]*(nx*ny)+j[:,1:,:]*nx+i[:,:,1:],
k[1:,:,:]*(nx*ny)+j[:,1:,:]*nx+i[:,:,:-1]], dtype=int)

This could also be improved by using one of numpy's stacking functions
to put the matrices together instead of the catch-all
too-clever-for-its-own-good "array".

The point is, elementwise operations and broadcasting make it easy to