<div dir="ltr">Well, one tip to start with:<div><br></div><div><span style="color:rgb(33,33,33);font-family:Courier;font-size:14.6667px">numpy.where(some_comparison, True, False)</span> <br></div><div><br></div><div>is the same as but slower than</div><div><br></div><div><span style="color:rgb(33,33,33);font-family:Courier;font-size:14.6667px">some_comparison</span> <br></div><div><br></div><div>Eric</div></div><br><div class="gmail_quote"><div dir="ltr">On Wed, 28 Mar 2018 at 18:36 Moroney, Catherine M (398E) <<a href="mailto:Catherine.M.Moroney@jpl.nasa.gov">Catherine.M.Moroney@jpl.nasa.gov</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div lang="EN-US" link="#0563C1" vlink="#954F72">
<div class="m_-3945452794008495637WordSection1">
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Helvetica">Hello,<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Helvetica"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Helvetica">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.<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Helvetica"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Helvetica">Here is the raw Python code:
<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Helvetica"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier">def mixed_coastline_slow(nsidc, radius, count, mask=None):<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"> nsidc_copy = numpy.copy(nsidc)<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"> if (mask is None):<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"> idx_coastline = numpy.where(nsidc_copy == NSIDC_COASTLINE_MIXED)<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"> else:<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"> idx_coastline = numpy.where(mask & (nsidc_copy == NSIDC_COASTLINE_MIXED))<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"> <u></u>
<u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"> for (irow0, icol0) in zip(idx_coastline[0], idx_coastline[1]):<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"> rows = ( max(irow0-radius, 0), min(irow0+radius+1, nsidc_copy.shape[0]) )<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"> cols = ( max(icol0-radius, 0), min(icol0+radius+1, nsidc_copy.shape[1]) )<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"> window = nsidc[rows[0]:rows[1], cols[0]:cols[1]]<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"> npoints = numpy.where(window != NSIDC_COASTLINE_MIXED, True, False).sum()<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"> nsnowice = numpy.where( (window >= NSIDC_SEAICE_LOW) & (window <= NSIDC_FRESHSNOW), \<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"> True, False).sum()<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"> if (100.0*nsnowice/npoints >= count):<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"> nsidc_copy[irow0, icol0] = MISR_SEAICE_THRESHOLD<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"> return nsidc_copy<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Helvetica">and here is my attempt at Cython-izing it:<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Helvetica"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier">import numpy<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier">cimport numpy as cnumpy<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier">cimport cython<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier">cdef int NSIDC_SIZE = 721<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier">cdef int NSIDC_NO_SNOW = 0<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier">cdef int NSIDC_ALL_SNOW = 100<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier">cdef int NSIDC_FRESHSNOW = 103<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier">cdef int NSIDC_PERMSNOW = 101<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier">cdef int NSIDC_SEAICE_LOW = 1<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier">cdef int NSIDC_SEAICE_HIGH = 100<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier">cdef int NSIDC_COASTLINE_MIXED = 252<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier">cdef int NSIDC_SUSPECT_ICE = 253<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier">cdef int MISR_SEAICE_THRESHOLD = 6<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier">def mixed_coastline(cnumpy.ndarray[cnumpy.uint8_t, ndim=2] nsidc, int radius, int count):<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"> cdef int irow, icol, irow1, irow2, icol1, icol2, npoints, nsnowice<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"> cdef cnumpy.ndarray[cnumpy.uint8_t, ndim=2] nsidc2 \<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"> = numpy.empty(shape=(NSIDC_SIZE, NSIDC_SIZE), dtype=numpy.uint8)<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"> cdef cnumpy.ndarray[cnumpy.uint8_t, ndim=2] window \<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"> = numpy.empty(shape=(2*radius+1, 2*radius+1), dtype=numpy.uint8)<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"> nsidc2 = numpy.copy(nsidc)<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"> idx_coastline = numpy.where(nsidc2 == NSIDC_COASTLINE_MIXED)<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"> <u></u>
<u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"> for (irow, icol) in zip(idx_coastline[0], idx_coastline[1]):<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"> irow1 = max(irow-radius, 0)<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"> irow2 = min(irow+radius+1, NSIDC_SIZE)<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"> icol1 = max(icol-radius, 0)<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"> icol2 = min(icol+radius+1, NSIDC_SIZE)<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"> window = nsidc[irow1:irow2, icol1:icol2]<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"> npoints = numpy.where(window != NSIDC_COASTLINE_MIXED, True, False).sum()<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"> nsnowice = numpy.where( (window >= NSIDC_SEAICE_LOW) & (window <= NSIDC_FRESHSNOW), \<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"> True, False).sum()<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"> if (100.0*nsnowice/npoints >= count):<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"> nsidc2[irow, icol] = MISR_SEAICE_THRESHOLD<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"> return nsidc2<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier">Thanks in advance for any advice!<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier">Catherine<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Courier"><u></u> <u></u></span></p>
</div>
</div>
_______________________________________________<br>
NumPy-Discussion mailing list<br>
<a href="mailto:NumPy-Discussion@python.org" target="_blank">NumPy-Discussion@python.org</a><br>
<a href="https://mail.python.org/mailman/listinfo/numpy-discussion" rel="noreferrer" target="_blank">https://mail.python.org/mailman/listinfo/numpy-discussion</a><br>
</blockquote></div>