[SciPy-user] Finding local minima of greater than a given depth

[Zachary Pincus]

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 at scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-user