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

Thu Mar 29 11:14:16 EDT 2018

Install snakeviz to visualize what’s taking all the time.

You might want to check out numba.jit(nopython) for optimizing specific
sections.

On Wed, Mar 28, 2018 at 9:10 PM Joseph Fox-Rabinovitz <
jfoxrabinovitz at gmail.com> wrote:

> It looks like you are creating a coastline mask (or a coastline mask +
> some other mask), and computing the ratio of two quantities in a
> particular window around each point. If your coastline covers a
> sufficiently large portion of the image, you may get quite a bit of
> mileage using an efficient convolution instead of summing the windows
> directly. For example, you could use scipy.signal.convolve2d with
> inputs being (nsidc_copy != NSIDC_COASTLINE_MIXED), (nsidc_copy ==
> NSIDC_SEAICE_LOW & nsdic_copy == NSIDC_FRESHSNOW) for the frst array,
> and a (2*radius x 2*radius) array of ones for the second. You may have
> to center the block of ones in an array of zeros the same size as
> nsdic_copy, but I am not sure about that.
>
> Another option you may want to try is implementing your window
> movement more efficiently. If you step your window center along using
> an algorithm like flood-fill, you can insure that there will be very
> large overlap between successive steps (even if there is a break in
> the coastline). That means that you can reuse most of the data you've
> extracted. You will only need to subtract off the non-overlapping
> portion of the previous window and add in the non-overlapping portion
> of the updated window. If radius is 16, giving you a 32x32 window, you
> go from summing ~1000 pixels per quantity of interest, to summing only
> ~120 if the window moves along a diagonal, and only 64 if it moves
> vertically or horizontally. While an algorithm like this will probably
> give you the greatest boost, it is a pain to implement.
>
> If I had to guess, this looks like L2 processing for a multi-spectral
> instrument. If you don't mind me asking, what mission is this for? I'm
> working on space-looking detectors at the moment, but have spent many
> years on the L0, L1b and L1 portions of the GOES-R ground system.
>
> - Joe
>
> On Wed, Mar 28, 2018 at 9:43 PM, Eric Wieser
> <wieser.eric+numpy at gmail.com> wrote:
> > Well, one tip to start with:
> >
> > numpy.where(some_comparison, True, False)
> >
> > is the same as but slower than
> >
> > some_comparison
> >
> > Eric
> >
> > On Wed, 28 Mar 2018 at 18:36 Moroney, Catherine M (398E)
> > <Catherine.M.Moroney at jpl.nasa.gov> wrote:
> >>
> >> Hello,
> >>
> >>
> >>
> >> 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.  But I would still like to know
> if I
> >> could rewrite it in pure python/numpy/scipy or in Cython and get a
> similar
> >> speedup.
> >>
> >>
> >>
> >> 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
> >
> >
> > _______________________________________________
> > NumPy-Discussion mailing list
> > NumPy-Discussion at python.org
> > https://mail.python.org/mailman/listinfo/numpy-discussion
> >
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at python.org
> https://mail.python.org/mailman/listinfo/numpy-discussion
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20180329/9e88c0b9/attachment-0001.html>