mode='wrap' in find_peaks(), as in argrelextrema() – signals with periodic boundary condition

[bekemax]

Is your feature request related to a problem? Please describe.

I'm working with signals that have "periodic boundary condition" – if a signal is given with x, a numpy.array of shape (N,), then after its last, N-1-st element, x[N-1], the "next" element should be again x[0]. (So, x[N] is kind of forced to be equal to x[0]).

It would be great if I could use the find_peaks() function from scipy.signal with all its features in this setting, but if done straightforwardly, here's what that gives:

import numpy as np
from scipy.signal import find_peaks
import matplotlib.pyplot as plt

x = np.linspace(0,3*np.pi,50)
y = np.cos(x)

peaks, _ = find_peaks(y)

plt.plot(y,'.')
plt.plot(peaks, y[peaks], "x", color='red')

From the standard point of view, the 0-th element of y is not a peak, since it has no "left" neighbour to be larger of. But from the "periodic boundary condition" point of view, it is larger than both the 1-st and the -1-st element, so it is a peak.

Describe the solution you'd like
It would be great to have an argument in find_peaks() similar to that of mode='wrap' of the argrelextrema() function from scipy.signal.

Describe alternatives you've considered
If I use scipy.signal.argrelextrema() with mode='wrap', I get the desired result:

from scipy.signal import argrelextrema

xmax_clip = argrelextrema(y, comparator = np.greater, mode = 'clip')[0]
xmax_wrap = argrelextrema(y, comparator = np.greater, mode = 'wrap')[0]

fig, (ax1, ax2) = plt.subplots(1, 2)

ax1.plot(y, '.')
ax1.plot(xmax_clip, y[xmax_clip], "x", color='red')
ax1.set_title('mode=clip (default)')

ax2.plot(y, '.')
ax2.plot(xmax_wrap, y[xmax_wrap], "x", color='red')
ax2.set_title('mode=wrap')

– it identifies the 0-th element as a local maximum as well, since it is larger than both the1-st and -1-st elements.

But find_peaks() has far more functionality than argrelextrema() – the former returns the properties of the peaks, has helpful arguments like distance – to find peaks that are not closer than a given distance, etc. Actually, distance is probably the only thing to take care of if such a mode='wrap' feature is implemented for find_peaks() – since the distance is now circular, not only in the positive direction of the indices: in mode='clip' the distance between the 0-th and the N-1-st element is N-1 (for an array of shape (N,)), while in mode='wrap' it is just 1. This should be taken care of by numpy.take(mode='wrap'), as in argrelextrema().

As a temporary workaround for my needs, I can, for example, find_peaks(), do an np.roll(signal, 1) of the signal array, find_peaks() again on a roll-ed array, then np.roll(signal, -1) it back, and then merge the two peaks sets that I've found. But implementing mode='wrap' for find_peaks() could be quite helpful for others.

Additional context (e.g. screenshots)

Python 3.7.3
numpy 1.16.3
scipy 1.3.0
matplotlib 3.1.0

[tupui]

Hi, thanks for the suggestion. It seems like you would know how to do this. Would you be willing to do a PR?

[bekemax]

Yep, I guess, I can!
Looks like now find_peaks() is using _local_maxima_1d() from _peak_finding_utils, which is a Cython-function that is said to be faster than argrelextrema(), so, I guess, for such a mode='wrap' case we can just use the already-implemented argrelextrema(), sacrificing speed a little bit.
But then the distance-checking logic is using _select_by_peak_distance() from _peak_finding_utils, which doesn't have such a "circular"-distance option.
Let me think a bit what can be done there and I'll be glad to do a PR

[tupui]

Great, thanks! Don't hesitate to post on the dev mail list to enquire for help/advice about this.

[johnharveymath]

Hi there! I wonder how far you got with this. I'm on this page because I have a similar issue with find_peaks. In my case, the time series ends with a nice high peak but then drops down slightly just before the end. As a result, the right_base is just the end of the window, and the calculated prominence is very low. This is not at all desirable in my use case.

I was thinking of adding a mode where we do not calculate prominence by reference to a value at the extreme end of the window unless, of course, both the left_base and right_base are at the end. I'm not super experienced with (a) Python and (b) contributing to large projects, but it seems like these issues should be handled together. I am happy to contribute, as it seems like it's within my competence!

[lagru]

My current thoughts on this:

This request might be solved by adapting the indexing logic on the Cython level. For _local_maxima_1d this would be quite easy (if not elegant) by adding some logic to compare the first and last sample if something like mode = 'wrap' is provided. However for _select_by_peak_distance this logic would get way more complex, same for _peak_prominences and _peak_widths. Perhaps a helper function could be used to handle the indexing, wrapping and break conditions of the loops? Alternatively, enabling the directive #cython: wraparound=True (currently False) might provide a solution. However this is set at compile time and I think it comes at the cost of indexing speed.

The actual problem is not adapting the part that works similar to argrelextrema but the additional property calculation and peak filtering. If memory is not a concern I would probably advise to wrap the input signal x like @bekemax suggested.

However I can understand why this would be a nice enhancement. If I can think of a solution without large drawbacks I will follow up. @bekemax, feel free to tag me if you come up with a solution.

[nahuelalmeira]

Hi, I would like to know if there has been any progress related to this enhancement. Thanks!

[lueck1803]

Hey,
my workaround is to find a mid value of the array, split it into half and put the second half on the first position so that definitly no real local extrema is at the border of the array (which would be seen as a global extrema).

Then I apply find_peaks as usual. Afterwards I shift the gained indice list of the extrema by the length of the second half splitted list I putted in front of the first half, to get the indices of the original array.

	## get the index of the middle value of an unsorted list
	def get_index_of_midval(this_list):
		mean = np.mean(this_list)
		mid = None
		for i, el in enumerate(this_list):
			dif_el = np.abs(el-mean)
			if(mid == None):
				mid = (i,el)
			else:
				dif_mid = np.abs(mid[1] - mean)
				if(dif_el < dif_mid):
					mid = (i,el)
		return mid[0]

            splitindex = get_index_of_midval(NP)

	splited1 = NP[0:splitindex]
        splited2 = NP[splitindex:]

            ## the shifted list where the extrema arent at the borders
	shifted = np.concatenate((splited2,splited1))

            extremaIndices, _ = find_peaks(shifted)
                      
            extremaIndices = extremaIndices - len(splited2)

[lueck1803]

It's not 100% robust. It depends on your data and if you really want all local extrema.
For an periodic signal like a sinus based one, its really good enough.

For my case I have particles in a box with periodic boundaries and consider the density distribution.

For me it is better to take like 70% of max value. I changed to:

	## get the index of the value of an unsorted list which is closesd to criteria 
	def get_index_of_criteria(this_list,criteria):
		if(criteria == "mean"):
			crit = np.mean(this_list)
		elif(isinstance(criteria,tuple) and isinstance(criteria[0],str) and isinstance(criteria[1],(int,float))):
			if(criteria[0]=="max"):
				crit = max(this_list)*criteria[1]
			elif(criteria[0]=="min"):
				crit = min(this_list)*criteria[1]
			else:
				raise Exception(f"The criteria {criteria} is not available")
		else:
			raise Exception(f"The criteria {criteria} is not available")

		result = None
		for i, el in enumerate(this_list):
			dif_el = np.abs(el-crit)
			if(result == None):
				result = (i,el)
			else:
				dif_res = np.abs(result[1] - crit)
				if(dif_el < dif_res):
					result = (i,el)

		return result[0]
	#######################################################