Improve/clarify `Evoked.get_peak()`

[cbrnr]

As discussed with @drammock here, we should discuss if the current implementation of Evoked.get_peak() is sufficient or if we should improve it with a true peak detection. Currently, it only finds the min/max within an array. If we decide that the current behavior is good enough, we should at least put a note in the docstring.

[agramfort]

do you have something in mind to do something smarter for peak detection?

+1 to improve the docstring anyway

[drammock]

scipy.signal.find_peaks is the main option to consider I think

[agramfort]

ok but here we want only one peak no?

[drammock]

Yes, we want only one peak, but it is easy to use find_peaks and then select only the largest one. It also has a lot of configurable options so it would need some thought about what is optimal

[cbrnr]

👍 for scipy.signal.find_peaks, and we should probably keep the old method around by introducting a method parameter. WDYT?

[agramfort]

how should we evaluate if new suggested ways are better?

[cbrnr]

I'd make the new method the default, and maybe issue a warning that it has changed? Not sure how many people are using it.

[drammock]

how should we evaluate if new suggested ways are better?

I think results will change only in the case where the most extreme value is on an up- or down-slope at the edge of the window. If that is correct I'd view this as a bugfix. Some testing is needed though; if it turns out to give different answers in more cases, we'd need to justify the change as @agramfort implies

[agramfort]

I think results will change only in the case where the most extreme value is on an up- or down-slope at the edge of the window. If that is correct I'd view this as a bugfix.

should we just warn if this happens?

what does scipy do in this case?

Some testing is needed though; if it turns out to give different answers in more cases, we'd need to justify the change as @agramfort https://github.com/agramfort implies

yes it needs testing on real data.

[cbrnr]

What if we leave the default and just add this as a new method?

Which tests/evaluations are we talking about? Are there tests in MNE that use get_peak() to assert certain conditions? Or were you referring to tests that we don't currently have?

[agramfort]

What if we leave the default and just add this as a new method?

Which tests/evaluations are we talking about? Are there tests in MNE that use get_peak() to assert certain conditions? Or were you referring to tests that we don't currently have?

I would for example take 600 subjects from the cam-can dataset and compare the peak estimates with a scatter plot

[drammock]

I think results will change only in the case where the most extreme value is on an up- or down-slope at the edge of the window. If that is correct I'd view this as a bugfix.

should we just warn if this happens?

I guess that could be a simpler solution (check whether the "peak" is the first or last sample of the time window, and warn if so).

what does scipy do in this case?

scipy would not return a "peak" that is the first or last sample because a peak must have 2 neighbors that are lower (or for a plateau peak, the plateau neighbors must both be lower):

In [1]: from scipy.signal import find_peaks
In [2]: y = [10, 2, 1, 4, 1]
In [3]: find_peaks(y)
Out[3]: (array([3]), {})

What if we leave the default and just add this as a new method?

I think at a minimum we should add a warning for edge samples. Almost surely better would be switching to scipy.signal.find_peaks but someone needs to make the time to do some testing. I'm not even sure you need to test 600 subjects, it ought to be enough to take the sample data and test 3 different channel types, 5 channels of each, 5-10 different window widths, at 5-10 different places within the recording (and do this for both raw and evoked).

[agramfort]

I guess that could be a simpler solution (check whether the "peak" is the first or last sample of the time window, and warn if so).

let's do a first PR with this warning. Then scipy can be done in a next PR

[jona-sassenhagen]

Vaguely related - both plot_joint and joint TFR plots use automatic peak detection if users don't specify coordinates for topoplots themselves.

https://github.com/mne-tools/mne-python/blob/616b9430392d0d2bbf73063e738c7b8c378defcd/mne/time_frequency/tfr.py#L2592 https://github.com/mne-tools/mne-python/blob/44ac64b42a942846061e42a2f0064f430df56153/mne/viz/utils.py#L758