[SciPy-user] Finding local minima of greater than a given depth
Zachary Pincus
zachary.pincus@yale....
Thu Aug 14 22:49:04 CDT 2008
Hi,
Not that it likely matters in this case, but below is a version of the
local_minima function that uses some more advanced numpy features,
like "fancy indexing". These features really make life easier in many
cases (and are usually faster than explicit loops), so it's really
worth learning how they work:
def local_minima_fancy(fits, window=15):
from scipy.ndimage import minimum_filter
fits = numpy.asarray(fits)
minfits = minimum_filter(fits, size=window, mode="wrap")
minima_mask = fits == minfits
good_indices = numpy.arange(len(fits))[minima_mask]
good_fits = fits[minima_mask]
order = good_fits.argsort()
return good_indices[order], good_fits[order]
We have two types of fancy indexing here. The first takes a boolean
"mask":
good_fits = fits[minima_mask]
returns only the elements of fits where the minima_mask array is true.
It's equivalent to:
good_fits = numpy.array([fits[i] for i, m in enumerate(minima_mask)
if m])
The second takes a list/array of indices:
return good_indices[order]
is equivalent to:
return numpy.array([good_indices[i] for i in order])
Also note that the original function you sent has a slight bug when
there are multiple minima with the same value: list.index returns the
index of the *first* entry in the list with that value, so the indices
of later minima will be incorrect. Plus, the function does some
unnecessary work:
good_fits = [ fit for fit in minima ]
makes an unneeded copy of minima, and
good_indices = [ fits.index(fit) for fit in minima ]
requires traversing the fits list once for each minima (the "index"
method runs a linear search), when only one traversal should be
required.
Here's a fixed and tuned-up version with just list processing:
def local_minima_fixed(fits, window=15):
from scipy.ndimage.filters import minimum_filter as min_filter
minfits = min_filter(fits, size=window, mode="wrap")
minima_and_indices = []
for i, (fit, minfit) in enumerate(zip(fits, minfits)):
if fit == minfit:
minima_and_indices.append([fit, i])
minima_and_indices.sort()
good_fits, good_indices = zip(*minima_and_indices)
return good_indices, good_fits
For reference, here are some timings (made with ipython's excellent
timeit magic command):
In [60]: a = numpy.random.randint(400, size=2000)
In [61]: timeit local_minima(list(a), 5)
100 loops, best of 3: 7.38 ms per loop
In [62]: timeit local_minima_fixed(list(a), 5)
100 loops, best of 3: 2.69 ms per loop
In [63]: timeit local_minima_fancy(list(a), 5)
1000 loops, best of 3: 973 [micro]s per loop
As above, the speed of this routine probably doesn't matter too much,
but it's a useful exercise to understand how and why these other
versions work (if some step doesn't make sense, try working through
the steps in the interpreter to see what's happening -- this is a good
way to learn some nice features of python and numpy), and how the
"fancy" version uses numpy's advanced features to good effect.
Zach
> def local_minima(fits, window=15): #{{{
> """
> Find the local minima within fits, and return them and their
> indices.
>
> Returns a list of indices at which the minima were found, and a
> list of the
> minima, sorted in order of increasing minimum. The keyword
> argument window
> determines how close two local minima are allowed to be to one
> another. If
> two local minima are found closer together than that, then the
> lowest of
> them is taken as the real minimum. window=1 will return all
> local minima.
>
> """
> from scipy.ndimage.filters import minimum_filter as min_filter
>
> minfits = min_filter(fits, size=window, mode="wrap")
>
> minima = []
> for i in range(len(fits)):
> if fits[i] == minfits[i]:
> minima.append(fits[i])
>
> minima.sort()
>
> good_indices = [ fits.index(fit) for fit in minima ]
> good_fits = [ fit for fit in minima ]
>
> return(good_indices, good_fits)
>
> _______________________________________________
> SciPy-user mailing list
> SciPy-user@scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-user
More information about the SciPy-user
mailing list