[SciPy-user] Extrema finding

[iCy-fLaME]

On Thu, Oct 9, 2008 at 1:21 PM, John Hunter <jdh2358 at gmail.com> wrote:
> On Thu, Oct 9, 2008 at 6:44 AM, iCy-fLaME <icy.flame.gm at gmail.com> wrote:
>> I am trying to find a list of all maxima and minima, each for a given
>> (1D) numpy array. Anyone know of a quick way to do it?
>
>
>>
>> Ideally the function will return the extrema values and their
>> positions. Relatively simple function to implement in Python, but that
>> would be painfully slow. The typical data array I am looking at, has
>> approximately 500k elements of double precision float.
>
> I use the following to find the indices of extrema:
>
>
> def extrema_sampling(x, withend=False):
>    """
>    return the indices into the local maxima or minima of x
>    if withend, include the 0 and end points in the sampling
>    """
>    d = np.diff(x)
>    d2 = d[:-1]*d[1:]
>    ind = []
>    if withend:
>        ind.append(0)
>    ind.extend(np.nonzero(d2<0)[0]+1)
>    if withend and ind[-1]!=len(x)-1:
>        ind.append(len(x)-1)
>    return np.array(ind)

This is an excellent way of indexing extrema! Although it doesnt
return exactly what I wanted and it wont catch gradient changes to
zero.

Inspired by the above, I have come up with the following, hope it can
be of some use to others:


def extrema(x, max = True, min = True, strict = False, withend = False):
	"""
	This function will index the extrema of a given array x.
	
	Options:
		max		If true, will index maxima
		min		If true, will index minima
		strict		If true, will not index changes to zero gradient
		withend	If true, always include x[0] and x[-1]
	
	This function will return a tuple of extrema indexies and values
	"""
	
	# This is the gradient
	from numpy import zeros
	dx = zeros(len(x))
	from numpy import diff
	dx[1:] = diff(x)
	dx[0] = dx[1]
	
	# Clean up the gradient in order to pick out any change of sign
	from numpy import sign
	dx = sign(dx)
	
	# define the threshold for whether to pick out changes to zero gradient
	threshold = 0
	if strict:
		threshold = 1
		
	# Second order diff to pick out the spikes
	d2x = diff(dx)
	
	if max and min:
		d2x = abs(d2x)
	elif max:
		d2x = -d2x
	
	# Take care of the two ends
	if withend:
		d2x[0] = 2
		d2x[-1] = 2
	
	# Sift out the list of extremas
	from numpy import nonzero
	ind = nonzero(d2x > threshold)[0]
	
	return ind, x[ind]