Hi Catherine,

One problem with sliding window algorithms is that the straightforward approach can be very inefficient. Ideally you would want to not recompute your windowed quantity from all points in the window, but to reuse the result from an overlapping window and only take into account the points that have changed in the sliding of the window. In your case this can be efficiently done using a summed area table. Consider these two auxiliary functions:

def summed_area_table(array):

rows, cols = array.shape

out = np.zeros((rows + 1, cols + 1), np.intp)

np.cumsum(array, axis=0, out=out[1:, 1:])

np.cumsum(out[1:, 1:], axis=1, out=out[1:, 1:])

return out

def windowed_sum_from_summed_area_table(array, size):

sat = summed_area_table(array)

return (sat[:-size, :-size] + sat[size:, size:] - sat[:-size, size:] -

sat[size:, -size:])

Using these, you can compute npoints and nsnowice for all points in your input nsidc array as:

mask_coastline = nsidc == NSIDC_COASTLINE_MIXED

mask_not_coastline = ~mask_coastline

mask_snowice = (nsidc >= NSIDC_SEAICE_LOW) & (nsidc <= NSIDC_FRESHSNOW)

nsnowice = windowed_sum_from_summed_area_table(mask_snowice, 2*radius + 1)

npoints = windowed_sum_from_summed_area_table(mask_not_coastline, 2*radius + 1)

From here it should be more or less straightforward to reproduce the rest of your calculations. As written this code only handles points a distance of at least radius from an array edge. If the edges are important to you, they can also be extracted from the summed area table, but the expressions get ugly: it may be cleaner, even if slower, to pad the masks with zeros before summing them up. Also, if the fraction of points that are in mask_coastline is very small, you may be doing way too many unnecessary calculations.

Good luck!

Jaime

On Thu, Mar 29, 2018 at 3:36 AM Moroney, Catherine M (398E) <Catherine.M.Moroney@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@python.org
https://mail.python.org/mailman/listinfo/numpy-discussion