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@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... 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...

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@python.org https://mail.python.org/mailman/listinfo/numpy-discussion

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 https://en.wikipedia.org/wiki/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