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.moro...@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
>
--
(\__/)
( O.o)
( > <) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes
de dominación mundial.
_______________________________________________
NumPy-Discussion mailing list
NumPy-Discussion@python.org
https://mail.python.org/mailman/listinfo/numpy-discussion
