[Numpy-discussion] trying to speed up the following....

[Brennan Williams]

Robert Kern wrote:
> On Wed, Mar 25, 2009 at 00:09, Brennan Williams
> <brennan.williams at visualreservoir.com> wrote:
>   
>> Robert Kern wrote:
>>     
>>> On Tue, Mar 24, 2009 at 18:29, Brennan Williams
>>> <brennan.williams at 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.
>
>