Robert Kern wrote: > On Wed, Mar 25, 2009 at 00:09, Brennan Williams > <brennan.willi...@visualreservoir.com> wrote: > >> Robert Kern wrote: >> >>> On Tue, Mar 24, 2009 at 18:29, Brennan Williams >>> <brennan.willi...@visualreservoir.com> wrote: >>> >>> >>>> I have an array (porvatt.yarray) of ni*nj*nk values. >>>> I want to create two further arrays. >>>> >>>> activeatt.yarray is of size ni*nj*nk and is a pointer array to an active >>>> cell number. If a cell is inactive then its activeatt.yarray value will be >>>> 0 >>>> >>>> ijkatt.yarray is of size nactive, the number of active cells (which I >>>> already know). ijkatt.yarray holds the ijk cell number for each active >>>> cell. >>>> >>>> >>>> My code looks something like... >>>> >>>> activeatt.yarray=zeros(ncells,dtype=int) >>>> ijkatt.yarray=zeros(nactivecells,dtype=int) >>>> >>>> iactive=-1 >>>> ni=currentgrid.ni >>>> nj=currentgrid.nj >>>> nk=currentgrid.nk >>>> for ijk in range(0,ni*nj*nk): >>>> if porvatt.yarray[ijk]>0: >>>> iactive+=1 >>>> activeatt.yarray[ijk]=iactive >>>> ijkatt.yarray[iactive]=ijk >>>> >>>> I may often have a million+ cells. >>>> So the code above is slow. >>>> How can I speed it up? >>>> >>>> >>> mask = (porvatt.yarray.flat > 0) >>> ijkatt.yarray = np.nonzero(mask) >>> >>> # This is not what your code does, but what I think you want. >>> # Where porvatt.yarray is inactive, activeatt.yarray is -1. >>> # 0 might be an active cell. >>> activeatt.yarray = np.empty(ncells, dtype=int) >>> activeatt.yarray.fill(-1) >>> activeatt.yarray[mask] = ijkatt.yarray >>> >>> >>> >>> >> Thanks. Concise & fast. This is what I've got so far (minor mods from >> the above).... >> >> from numpy import * >> ... >> mask=porvatt.yarray>0.0 >> ijkatt.yarray=nonzero(mask)[0] >> activeindices=arange(0,ijkatt.yarray.size) >> activeatt.yarray = empty(ncells, dtype=int) >> activeatt.yarray.fill(-1) >> activeatt.yarray[mask] = activeindices >> >> I have... >> >> ijkatt.yarray=nonzero(mask)[0] >> >> because it looks like nonzero returns a tuple of arrays rather than an >> array. >> > > Yes. Apologies. > > Apology accepted. Don't do it again.
On a more serious note, it is clear that, as expected, operating on elements of an array inside a Python for loop is slow for large arrays. Soon I will be writing an import interface to read corner point grid geometries and I'm currently looking at vtk unstructured grids etc. Most of the numpy vectorization is aimed at relatively simply structured arrays on the basis that you'll never meet everyone's needs/data structures. So I presume that if I find I have a bottleneck in my code which looks specific to my data structures I should then look at offloading that to C or Fortran? (assuming I can't find it in numpy or scipy). I'm already doing this to read in data from Fortran binary files although I actually decided to code it in C rather than use Fortran. >> I used >> >> activeindices=arange(0,ijkatt.yarray.size) >> >> and >> >> activeatt.yarray[mask] = activeindices >> > > Yes. You are correct. > > _______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion