Finding local minima of greater than a given depth
![](https://secure.gravatar.com/avatar/78bed16eb45e28e567b06e56a1407acc.jpg?s=120&d=mm&r=g)
Is there a function within scipy somewhere which will, given an array representing values of a function, find all the local minima having a depth greater than some specified minimum? The following works great for smooth functions, but when the data has noise in it, it also returns all of the (very) local minima, which I don't want. The functions I'm working with are periodic (hence the modulo in the indices for endpoint cases). Or, if there isn't such a built in functionality, what's the right way to measure the depth of a local minimum? def local_minima(fitlist): minima = [] for i in range(len(fitlist)): if fitlist[i] < fitlist[mod(i+1,len(fitlist))] and fitlist[i] < fitlist[mod(i-1,len(fitlist))]: minima.append(fitlist[i]) minima.sort() good_indices = [ fitlist.index(fit) for fit in minima ] good_fits = [ fit for fit in minima ] return(good_indices, good_fits) -- Zane Selvans Amateur Earthling http://zaneselvans.org zane@ideotrope.org 303/815-6866 PGP Key: 55E0815F
![](https://secure.gravatar.com/avatar/790fd6e3ccfeaf2865c9294d2df0fa3e.jpg?s=120&d=mm&r=g)
Is there a function within scipy somewhere which will, given an array representing values of a function, find all the local minima having a depth greater than some specified minimum? The following works great for smooth functions, but when the data has noise in it, it also returns all of the (very) local minima, which I don't want.
You could low-pass filter out the (presumably high frequency) noise, which might introduce a slight phase change to your data. But you can use the local minima of the filtered data as good starting points for a better search in your original data. You might need to fit local polynomials to your data near these to find the minimum without introducing too much statistical bias (e.g. by just taking the smallest data point, which might only be the smallest because of the noise). There isn't a scipy function to do this in one go, but there are filter functions in scipy.signal and some polynomial fitting functionality recently added by Anne Archibald which I use for this kind of problem (search these archives for reference to that). The depth of the minimum is better defined once you fit a function to the region, but the appropriate size of that region is context dependent. It could be meaningfully measured w.r.t. the next nearest local maximum, or to a global average or maximum level derived from the data. There's no one way to do it. -Rob
![](https://secure.gravatar.com/avatar/afdaaab755ef79ac9e1374882d60ae9f.jpg?s=120&d=mm&r=g)
Is there a function within scipy somewhere which will, given an array representing values of a function, find all the local minima having a depth greater than some specified minimum? The following works great for smooth functions, but when the data has noise in it, it also returns all of the (very) local minima, which I don't want. The functions I'm working with are periodic (hence the modulo in the indices for endpoint cases). Or, if there isn't such a built in functionality, what's the right way to measure the depth of a local minimum?
You could measure "depth" of minima (of a 1D array) by also finding the flanking maxima and looking at the distance between them. Or any of the other methods Rob suggested. Another way to find local minima in a noise-robust manner that I've often seen is to not look for a minimum "depth", but for a minimum distance between minima. This is easy to implement using scipy.ndimage's minimum filter, which sets each element of an array to the minimum value seen over a specified neighborhood of that element. Then you just check for array elements where the element is equal to the minimum in the neighborhood... I'd also suggest smoothing the data a bit with a gaussian to get rid or some of the noise. Scipy.ndimage also provides these filters. Zach PS. Here's my implementation... it returns the indices of the local maxima in a list. Also, the min-distance is in terms of manhattan distance, not euclidian, so be warned. For a 2D array, the returned list will have two elements -- the row- indices of the maxima and the column-indices of the maxima. There's probably a better way to do that, but this is what I have. def local_maxima(array, min_distance = 1, periodic=False, edges_allowed=True): """Find all local maxima of the array, separated by at least min_distance.""" import scipy.ndimage as ndimage array = numpy.asarray(array) cval = 0 if periodic: mode = 'wrap' elif edges_allowed: mode = 'nearest' else: mode = 'constant' cval = array.max()+1 max_points = array == ndimage.maximum_filter(array, 1+2*min_distance, mode=mode, cval=cval) return [indices[max_points] for indices in numpy.indices(array.shape)] On Aug 14, 2008, at 3:40 PM, Zane Selvans wrote:
Is there a function within scipy somewhere which will, given an array representing values of a function, find all the local minima having a depth greater than some specified minimum? The following works great for smooth functions, but when the data has noise in it, it also returns all of the (very) local minima, which I don't want. The functions I'm working with are periodic (hence the modulo in the indices for endpoint cases). Or, if there isn't such a built in functionality, what's the right way to measure the depth of a local minimum?
def local_minima(fitlist): minima = []
for i in range(len(fitlist)): if fitlist[i] < fitlist[mod(i+1,len(fitlist))] and fitlist[i] < fitlist[mod(i-1,len(fitlist))]: minima.append(fitlist[i])
minima.sort()
good_indices = [ fitlist.index(fit) for fit in minima ] good_fits = [ fit for fit in minima ]
return(good_indices, good_fits)
-- Zane Selvans Amateur Earthling http://zaneselvans.org zane@ideotrope.org 303/815-6866 PGP Key: 55E0815F
_______________________________________________ SciPy-user mailing list SciPy-user@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-user
![](https://secure.gravatar.com/avatar/78bed16eb45e28e567b06e56a1407acc.jpg?s=120&d=mm&r=g)
Zachary Pincus <zachary.pincus <at> yale.edu> writes:
Another way to find local minima in a noise-robust manner that I've often seen is to not look for a minimum "depth", but for a minimum distance between minima. This is easy to implement using scipy.ndimage's minimum filter, which sets each element of an array to the minimum value seen over a specified neighborhood of that element. Then you just check for array elements where the element is equal to the minimum in the neighborhood...
I'd also suggest smoothing the data a bit with a gaussian to get rid or some of the noise. Scipy.ndimage also provides these filters.
Great! This works well: 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)
![](https://secure.gravatar.com/avatar/afdaaab755ef79ac9e1374882d60ae9f.jpg?s=120&d=mm&r=g)
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
participants (3)
-
Rob Clewley
-
Zachary Pincus
-
Zane Selvans