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

[Anne Archibald]

On 22/04/07, fred <fredmfp at 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
transcribe your generation routine.

You could also, and again this may take some transposing, use something like
indices = arange(nx*ny*nz).reshape((nz,ny,nx))
so that indices[k,j,i] replaces k*(nx*ny)+j*nx+i. It'll be slower,
though hopefully more maintainable.


Anne M. Archibald