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

[fred]

Anne Archibald a écrit :
> 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):
>   
VTK does care 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)
>
>   
Ok, I think I got the trick.
Very powerful ! ;-)
> 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".
>
>   
I tried this:

           cells_array = 
vstack([hstack(VTK_HEXAHEDRON_NB_POINTS*ones((nz-1,ny-1,nx-1),dtype=int)),
                                  
hstack(k[:-1,:,:]*(nx*ny)+j[:,:-1,:]*nx+i[:,:,:-1]),
                                  
hstack(k[:-1,:,:]*(nx*ny)+j[:,:-1,:]*nx+i[:,:,1:]),
                                  
hstack(k[:-1,:,:]*(nx*ny)+j[:,1:,:]*nx+i[:,:,1:]),
                                  
hstack(k[:-1,:,:]*(nx*ny)+j[:,1:,:]*nx+i[:,:,:-1]),
                                  
hstack(k[1:,:,:]*(nx*ny)+j[:,:-1,:]*nx+i[:,:,:-1]),
                                  
hstack(k[1:,:,:]*(nx*ny)+j[:,:-1,:]*nx+i[:,:,1:]),
                                  
hstack(k[1:,:,:]*(nx*ny)+j[:,1:,:]*nx+i[:,:,1:]),
                                  
hstack(k[1:,:,:]*(nx*ny)+j[:,1:,:]*nx+i[:,:,:-1])]).transpose().ravel()

Written as this, I get a memory allocation error: for nx=ny=nz=200, it 
exceeds 2GB,
altough it works fine when array() used. Any idea ?
For this peculiar case, I get the same issue when ijk loops are used...
> 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.
>   
Ok, but the idea is to speed up the computing of the array.

Thanks.

Cheers,

-- 
http://scipy.org/FredericPetit