[Numpy-discussion] best way of speeding up a filtering-like algorithm

[Chris Barker]

sorry, not enough time to look closely, but a couple general comments:

On Wed, Mar 28, 2018 at 5:56 PM, Moroney, Catherine M (398E) <
Catherine.M.Moroney at jpl.nasa.gov> wrote:

> I have the following sample code (pretty simple algorithm that uses a
> rolling filter window) and am wondering what the best way is of speeding it
> up.  I tried rewriting it in Cython by pre-declaring the variables but that
> didn’t buy me a lot of time.  Then I rewrote it in Fortran (and compiled it
> with f2py) and now it’s lightning fast.
>

if done right, Cython should be almost as fast as Fortran, and just as fast
if you use the "restrict" correctly (which I hope can be done in Cython):

https://en.wikipedia.org/wiki/Pointer_aliasing


> But I would still like to know if I could rewrite it in pure
> python/numpy/scipy
>

you can use stride_tricks to make arrays "appear" to be N+1 D, to implement
windows without actually duplicating the data, and then use array
operations on them. This can buy a lot of speed, but will not be as fast
(by a factor of 10 or so) as Cython or Fortran

see:

https://github.com/PythonCHB/IRIS_Python_Class/blob/master/Numpy/code/filter_example.py
for and example in 1D



> or in Cython and get a similar speedup.
>
>
see above -- a direct port of your Fortran code to Cython should get you
within a factor of two or so of the Fortran, and then using "restrict" to
let the compiler know your pointers aren't aliased should get you the reset
of the way.

Here is an example of a Automatic Gain Control filter in 1D, iplimented in
numpy with stride_triks, and C and Cython and Fortran.

https://github.com/PythonCHB/IRIS_Python_Class/tree/master/Interfacing_C/agc_example

Note that in that example, I never got C or Cython as fast as Fortran --
but I think using "restrict" in the C would do it.

HTH,

-CHB



>
> Here is the raw Python code:
>
>
>
> def mixed_coastline_slow(nsidc, radius, count, mask=None):
>
>
>
>     nsidc_copy = numpy.copy(nsidc)
>
>
>
>     if (mask is None):
>
>         idx_coastline = numpy.where(nsidc_copy == NSIDC_COASTLINE_MIXED)
>
>     else:
>
>         idx_coastline = numpy.where(mask & (nsidc_copy ==
> NSIDC_COASTLINE_MIXED))
>
>
>
>     for (irow0, icol0) in zip(idx_coastline[0], idx_coastline[1]):
>
>
>
>         rows = ( max(irow0-radius, 0), min(irow0+radius+1,
> nsidc_copy.shape[0]) )
>
>         cols = ( max(icol0-radius, 0), min(icol0+radius+1,
> nsidc_copy.shape[1]) )
>
>         window = nsidc[rows[0]:rows[1], cols[0]:cols[1]]
>
>
>
>         npoints = numpy.where(window != NSIDC_COASTLINE_MIXED, True,
> False).sum()
>
>         nsnowice = numpy.where( (window >= NSIDC_SEAICE_LOW) & (window <=
> NSIDC_FRESHSNOW), \
>
>                                 True, False).sum()
>
>
>
>         if (100.0*nsnowice/npoints >= count):
>
>              nsidc_copy[irow0, icol0] = MISR_SEAICE_THRESHOLD
>
>
>
>     return nsidc_copy
>
>
>
> and here is my attempt at Cython-izing it:
>
>
>
> import numpy
>
> cimport numpy as cnumpy
>
> cimport cython
>
>
>
> cdef int NSIDC_SIZE  = 721
>
> cdef int NSIDC_NO_SNOW = 0
>
> cdef int NSIDC_ALL_SNOW = 100
>
> cdef int NSIDC_FRESHSNOW = 103
>
> cdef int NSIDC_PERMSNOW  = 101
>
> cdef int NSIDC_SEAICE_LOW  = 1
>
> cdef int NSIDC_SEAICE_HIGH = 100
>
> cdef int NSIDC_COASTLINE_MIXED = 252
>
> cdef int NSIDC_SUSPECT_ICE = 253
>
>
>
> cdef int MISR_SEAICE_THRESHOLD = 6
>
>
>
> def mixed_coastline(cnumpy.ndarray[cnumpy.uint8_t, ndim=2] nsidc, int
> radius, int count):
>
>
>
>      cdef int irow, icol, irow1, irow2, icol1, icol2, npoints, nsnowice
>
>      cdef cnumpy.ndarray[cnumpy.uint8_t, ndim=2] nsidc2 \
>
>         = numpy.empty(shape=(NSIDC_SIZE, NSIDC_SIZE), dtype=numpy.uint8)
>
>      cdef cnumpy.ndarray[cnumpy.uint8_t, ndim=2] window \
>
>         = numpy.empty(shape=(2*radius+1, 2*radius+1), dtype=numpy.uint8)
>
>
>
>      nsidc2 = numpy.copy(nsidc)
>
>
>
>      idx_coastline = numpy.where(nsidc2 == NSIDC_COASTLINE_MIXED)
>
>
>
>      for (irow, icol) in zip(idx_coastline[0], idx_coastline[1]):
>
>
>
>           irow1 = max(irow-radius, 0)
>
>           irow2 = min(irow+radius+1, NSIDC_SIZE)
>
>           icol1 = max(icol-radius, 0)
>
>           icol2 = min(icol+radius+1, NSIDC_SIZE)
>
>           window = nsidc[irow1:irow2, icol1:icol2]
>
>
>
>           npoints = numpy.where(window != NSIDC_COASTLINE_MIXED, True,
> False).sum()
>
>           nsnowice = numpy.where( (window >= NSIDC_SEAICE_LOW) & (window
> <= NSIDC_FRESHSNOW), \
>
>                                   True, False).sum()
>
>
>
>           if (100.0*nsnowice/npoints >= count):
>
>                nsidc2[irow, icol] = MISR_SEAICE_THRESHOLD
>
>
>
>      return nsidc2
>
>
>
> Thanks in advance for any advice!
>
>
>
> Catherine
>
>
>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at python.org
> https://mail.python.org/mailman/listinfo/numpy-discussion
>
>


-- 

Christopher Barker, Ph.D.
Oceanographer

Emergency Response Division
NOAA/NOS/OR&R            (206) 526-6959   voice
7600 Sand Point Way NE   (206) 526-6329   fax
Seattle, WA  98115       (206) 526-6317   main reception

Chris.Barker at noaa.gov
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20180329/4a4e6cd0/attachment.html>