From 142fbdd7ec3952109cb17453246d2e6450c43926 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Richard=20H=C3=B6chenberger?= Date: Sat, 17 Apr 2021 12:11:17 +0200 Subject: [PATCH 1/9] Introduce lower bounds for ExG peak detection thresholds Fixes #9273 --- doc/changes/latest.inc | 2 ++ mne/preprocessing/_peak_finder.py | 21 ++++++++++++++------- mne/preprocessing/ecg.py | 20 +++++++++++++++++--- mne/preprocessing/ica.py | 2 +- mne/preprocessing/tests/test_ecg.py | 17 +++++++++++++---- mne/preprocessing/tests/test_eog.py | 12 ++++++++++++ mne/preprocessing/tests/test_ssp.py | 6 +----- 7 files changed, 60 insertions(+), 20 deletions(-) diff --git a/doc/changes/latest.inc b/doc/changes/latest.inc index 173f1b11679..c72ce384acd 100644 --- a/doc/changes/latest.inc +++ b/doc/changes/latest.inc @@ -297,6 +297,8 @@ Bugs - Fix bug where ``backend='notebook'`` could not be used in :meth:`mne.SourceEstimate.plot` (:gh:`9305` by `Jean-Remi King`_) +- `mne.preprocessing.find_eog_events`, `mne.preprocessing.find_ecg_events`, `mne.preprocessing.create_eog_epochs`, and `mne.preprocessing.create_ecg_epochs` now have better handling for situations where the signal is flat or the provided threshold value is too small due to floating point imprecisions (:gh:`xxx` by `Richard Höchenberger`_) + - `mne.preprocessing.compute_proj_eog` and `mne.preprocessing.compute_proj_ecg` now return empty lists if no EOG or ECG events, respectively, could be found. Previously, we'd return ``None`` in these situations, which does not match the documented behavior of returning a list of projectors (:gh:`9277` by `Richard Höchenberger`_) API changes diff --git a/mne/preprocessing/_peak_finder.py b/mne/preprocessing/_peak_finder.py index 9bac177988a..894be1f566c 100644 --- a/mne/preprocessing/_peak_finder.py +++ b/mne/preprocessing/_peak_finder.py @@ -50,11 +50,19 @@ def peak_finder(x0, thresh=None, extrema=1, verbose=None): s = x0.size if x0.ndim >= 2 or s == 0: - raise ValueError('The input data must be a non empty 1D vector') + raise ValueError('The input data must be a non-empty 1D vector') + eps = 1e2 * np.finfo(x0.dtype).eps if thresh is None: thresh = (np.max(x0) - np.min(x0)) / 4 + if thresh < eps: + thresh = eps logger.debug('Peak finder automatic threshold: %0.2g' % (thresh,)) + elif thresh < eps: + raise ValueError( + f'The selected threshold value is too small for peak detection ' + f'due to inevitable floating-point imprecisions.\nPlease pick a ' + f'value >= {eps:g} (You passed thresh={thresh})') assert extrema in [-1, 1] @@ -160,10 +168,6 @@ def peak_finder(x0, thresh=None, extrema=1, verbose=None): peak_mags = [] peak_inds = [] - # Change sign of data if was finding minima - if extrema < 0: - peak_mags *= -1.0 - # ensure output type array if not isinstance(peak_inds, np.ndarray): peak_inds = np.atleast_1d(peak_inds).astype('int64') @@ -171,8 +175,11 @@ def peak_finder(x0, thresh=None, extrema=1, verbose=None): if not isinstance(peak_mags, np.ndarray): peak_mags = np.atleast_1d(peak_mags).astype('float64') - # Plot if no output desired - if len(peak_inds) == 0: + # Change sign of data if was finding minima + if extrema < 0 and peak_mags.size > 0: + peak_mags *= -1.0 + + if peak_inds.size == 0: logger.info('No significant peaks found') else: logger.info('Found %d significant peak%s' diff --git a/mne/preprocessing/ecg.py b/mne/preprocessing/ecg.py index 0ed1e0e52e7..3559c990f72 100644 --- a/mne/preprocessing/ecg.py +++ b/mne/preprocessing/ecg.py @@ -31,7 +31,7 @@ def qrs_detector(sfreq, ecg, thresh_value=0.6, levels=2.5, n_thresh=3, Sampling rate ecg : array ECG signal - thresh_value : float | str + thresh_value : float | 'auto' qrs detection threshold. Can also be "auto" for automatic selection of threshold. levels : float @@ -65,18 +65,32 @@ def qrs_detector(sfreq, ecg, thresh_value=0.6, levels=2.5, n_thresh=3, n_points = len(ecg_abs) + # get peak amplitudes in the 1st, 2nd, and 3rd second of the recording maxpt = np.empty(3) maxpt[0] = np.max(ecg_abs[:init]) maxpt[1] = np.max(ecg_abs[init:init * 2]) maxpt[2] = np.max(ecg_abs[init * 2:init * 3]) init_max = np.mean(maxpt) + eps = 1e2 * np.finfo(ecg_abs.dtype).eps + + if init_max < eps: + warn('The first 3 seconds of your ECG data appear to be flat. Please ' + 'inspect your data to ensure this is intentional.') + init_max = eps if thresh_value == 'auto': thresh_runs = np.arange(0.3, 1.1, 0.05) elif isinstance(thresh_value, str): raise ValueError('threshold value must be "auto" or a float') else: + if init_max * thresh_value < eps: + raise ValueError( + f'The selected threshold value is too small for peak ' + f'detection with the provided data, due to inevitable ' + f'floating-point imprecisions.\n' + f'Please pick a value >= {eps/init_max:g} ' + f'(You passed thresh={thresh_value})') thresh_runs = [thresh_value] # Try a few thresholds (or just one) @@ -100,7 +114,7 @@ def qrs_detector(sfreq, ecg, thresh_value=0.6, levels=2.5, n_thresh=3, else: ii += 1 - if len(rms) == 0: + if not rms: rms.append(0.0) time.append(0.0) time = np.array(time) @@ -151,7 +165,7 @@ def find_ecg_events(raw, event_id=999, ch_name=None, tstart=0.0, %(ecg_ch_name)s %(ecg_tstart)s %(ecg_filter_freqs)s - qrs_threshold : float | str + qrs_threshold : float | 'auto' Between 0 and 1. qrs detection threshold. Can also be "auto" to automatically choose the threshold that generates a reasonable number of heartbeats (40-160 beats / min). diff --git a/mne/preprocessing/ica.py b/mne/preprocessing/ica.py index 07a68e2dc62..e0ada5febb9 100644 --- a/mne/preprocessing/ica.py +++ b/mne/preprocessing/ica.py @@ -1246,7 +1246,7 @@ def find_bads_ecg(self, inst, ch_name=None, threshold='auto', start=None, The name of the channel to use for ECG peak detection. The argument is mandatory if the dataset contains no ECG channels. - threshold : float | str + threshold : float | 'auto' The value above which a feature is classified as outlier. If 'auto' and method is 'ctps', automatically compute the threshold. If 'auto' and method is 'correlation', defaults to 3.0. The default diff --git a/mne/preprocessing/tests/test_ecg.py b/mne/preprocessing/tests/test_ecg.py index 75148098b16..9569927bca8 100644 --- a/mne/preprocessing/tests/test_ecg.py +++ b/mne/preprocessing/tests/test_ecg.py @@ -86,12 +86,21 @@ def test_find_ecg(): assert 'ECG-SYN' not in raw.ch_names assert 'ECG-SYN' not in ecg_epochs.ch_names - # Test behavior if no peaks can be found -> achieve this by providing - # all-zero'd data - raw._data[ecg_idx] = 0. + # test behavior if no peaks can be found -> achieve this by setting a high + # threshold ecg_events, _, average_pulse, ecg = find_ecg_events( - raw, ch_name=raw.ch_names[ecg_idx], return_ecg=True + raw, ch_name=raw.ch_names[ecg_idx], return_ecg=True, qrs_threshold=1e5 ) assert ecg_events.size == 0 assert average_pulse == 0 assert np.allclose(ecg, np.zeros_like(ecg)) + + # test behavior for an entirely flat channel + raw._data[ecg_idx] = 0. + with pytest.warns(RuntimeWarning, match='first 3 seconds.*flat'): + ecg_events, _, average_pulse, ecg = find_ecg_events( + raw, ch_name=raw.ch_names[ecg_idx], return_ecg=True + ) + assert ecg_events.size == 0 + assert average_pulse == 0 + assert np.allclose(ecg, np.zeros_like(ecg)) diff --git a/mne/preprocessing/tests/test_eog.py b/mne/preprocessing/tests/test_eog.py index f0435c42878..8eb05ae4300 100644 --- a/mne/preprocessing/tests/test_eog.py +++ b/mne/preprocessing/tests/test_eog.py @@ -1,5 +1,6 @@ import os.path as op +import numpy as np import pytest from mne import Annotations @@ -39,3 +40,14 @@ def test_find_eog(): events = find_eog_events(raw, ch_name=['EEG 060', 'EOG 061']) assert len(events) == 4 + + # test too low threshold + eps = np.finfo(raw.get_data().dtype).eps + with pytest.raises(ValueError, match='threshold value is too small'): + find_eog_events(raw, ch_name=None, thresh=eps) + + # test with an entirely flat signal + eog = raw.copy().pick_types(meg=False, eog=True).load_data() + eog._data[:] = 0. + events = find_eog_events(eog, ch_name=None) + assert events.size == 0 diff --git a/mne/preprocessing/tests/test_ssp.py b/mne/preprocessing/tests/test_ssp.py index 4307dc9f6ea..934084717a4 100644 --- a/mne/preprocessing/tests/test_ssp.py +++ b/mne/preprocessing/tests/test_ssp.py @@ -63,11 +63,7 @@ def test_compute_proj_ecg(short_raw, average): projs, events, drop_log = compute_proj_ecg( raw, n_mag=2, n_grad=2, n_eeg=2, ch_name='MEG 1531', bads=[], average=average, avg_ref=True, no_proj=True, l_freq=None, - h_freq=None, tmax=dur_use, return_drop_log=True, - # XXX can be removed once - # XXX https://github.com/mne-tools/mne-python/issues/9273 - # XXX has been resolved: - qrs_threshold=1e-15) + h_freq=None, tmax=dur_use, return_drop_log=True) assert projs == [] assert len(events) == len(drop_log) From 1dffa5562317195edd0109c652996caa529fa9d3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Richard=20H=C3=B6chenberger?= Date: Sun, 18 Apr 2021 11:19:49 +0200 Subject: [PATCH 2/9] Remove unnecessary test case --- mne/preprocessing/tests/test_ecg.py | 1 - 1 file changed, 1 deletion(-) diff --git a/mne/preprocessing/tests/test_ecg.py b/mne/preprocessing/tests/test_ecg.py index 9569927bca8..fa18e66f0fc 100644 --- a/mne/preprocessing/tests/test_ecg.py +++ b/mne/preprocessing/tests/test_ecg.py @@ -93,7 +93,6 @@ def test_find_ecg(): ) assert ecg_events.size == 0 assert average_pulse == 0 - assert np.allclose(ecg, np.zeros_like(ecg)) # test behavior for an entirely flat channel raw._data[ecg_idx] = 0. From 698653de831a21991b85204be8e8101c2e008905 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Richard=20H=C3=B6chenberger?= Date: Sun, 18 Apr 2021 11:48:47 +0200 Subject: [PATCH 3/9] Handle eps for exact dtypes --- mne/preprocessing/_peak_finder.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/mne/preprocessing/_peak_finder.py b/mne/preprocessing/_peak_finder.py index 894be1f566c..cc7cd941769 100644 --- a/mne/preprocessing/_peak_finder.py +++ b/mne/preprocessing/_peak_finder.py @@ -52,7 +52,11 @@ def peak_finder(x0, thresh=None, extrema=1, verbose=None): if x0.ndim >= 2 or s == 0: raise ValueError('The input data must be a non-empty 1D vector') - eps = 1e2 * np.finfo(x0.dtype).eps + try: + eps = 1e2 * np.finfo(x0.dtype).eps + except ValueError: # int dtypes are exact + eps = 0. + if thresh is None: thresh = (np.max(x0) - np.min(x0)) / 4 if thresh < eps: From aa587a1a0c8ab5e2779406a66efe4fd569bd4402 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Richard=20H=C3=B6chenberger?= Date: Sun, 18 Apr 2021 11:53:16 +0200 Subject: [PATCH 4/9] Add a sanity check for ECG thresholds --- mne/preprocessing/_peak_finder.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/mne/preprocessing/_peak_finder.py b/mne/preprocessing/_peak_finder.py index cc7cd941769..d648afeaa2d 100644 --- a/mne/preprocessing/_peak_finder.py +++ b/mne/preprocessing/_peak_finder.py @@ -57,6 +57,9 @@ def peak_finder(x0, thresh=None, extrema=1, verbose=None): except ValueError: # int dtypes are exact eps = 0. + if thresh is not None: + assert thresh >= 0. + if thresh is None: thresh = (np.max(x0) - np.min(x0)) / 4 if thresh < eps: From 20361b451f83f427237fca692758f29ed4d04c6f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Richard=20H=C3=B6chenberger?= Date: Mon, 19 Apr 2021 10:21:11 +0200 Subject: [PATCH 5/9] Cast to float if need be --- mne/preprocessing/_peak_finder.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/mne/preprocessing/_peak_finder.py b/mne/preprocessing/_peak_finder.py index d648afeaa2d..69f02081f93 100644 --- a/mne/preprocessing/_peak_finder.py +++ b/mne/preprocessing/_peak_finder.py @@ -52,10 +52,11 @@ def peak_finder(x0, thresh=None, extrema=1, verbose=None): if x0.ndim >= 2 or s == 0: raise ValueError('The input data must be a non-empty 1D vector') - try: - eps = 1e2 * np.finfo(x0.dtype).eps - except ValueError: # int dtypes are exact - eps = 0. + # cast to float, but only if we must + if not isinstance(x0, (np.single, np.double, np.longdouble)): + x0 = x0.astype(float) + + eps = 1e2 * np.finfo(x0.dtype).eps if thresh is not None: assert thresh >= 0. From 640df6469934826b7910039f44ba69c8943276ab Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Richard=20H=C3=B6chenberger?= Date: Mon, 19 Apr 2021 10:39:52 +0200 Subject: [PATCH 6/9] Filter deprecation warning --- doc/conf.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/doc/conf.py b/doc/conf.py index de86ec90dfe..c69dfe80326 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -791,7 +791,15 @@ def reset_warnings(gallery_conf, fname): 'ignore', message="can't resolve package from", category=ImportWarning) warnings.filterwarnings( 'ignore', message='.*mne-realtime.*', category=DeprecationWarning) + warnings.filterwarnings( + 'ignore', message='.*mne-realtime.*', category=DeprecationWarning) + # https://github.com/mne-tools/mne-python/pull/9311#issuecomment-822042591 + warnings.filterwarnings( + 'ignore', module=r'sphinx\.jinja2glue', + message=".*the old name will be removed.*", + category=DeprecationWarning + ) reset_warnings(None, None) From 32aac0c14346cbc5fef730e02f9a7c6453410f12 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Richard=20H=C3=B6chenberger?= Date: Mon, 19 Apr 2021 10:41:19 +0200 Subject: [PATCH 7/9] Review suggestion --- mne/preprocessing/_peak_finder.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/mne/preprocessing/_peak_finder.py b/mne/preprocessing/_peak_finder.py index 69f02081f93..3335b075a9f 100644 --- a/mne/preprocessing/_peak_finder.py +++ b/mne/preprocessing/_peak_finder.py @@ -52,8 +52,7 @@ def peak_finder(x0, thresh=None, extrema=1, verbose=None): if x0.ndim >= 2 or s == 0: raise ValueError('The input data must be a non-empty 1D vector') - # cast to float, but only if we must - if not isinstance(x0, (np.single, np.double, np.longdouble)): + if x0.dtype.kind in 'uib': # integer dtype x0 = x0.astype(float) eps = 1e2 * np.finfo(x0.dtype).eps From 7eb97b69fad06e472d82900d4779a465815963d0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Richard=20H=C3=B6chenberger?= Date: Mon, 19 Apr 2021 10:44:06 +0200 Subject: [PATCH 8/9] flake --- doc/conf.py | 1 + 1 file changed, 1 insertion(+) diff --git a/doc/conf.py b/doc/conf.py index c69dfe80326..b22a0a3fe1f 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -801,6 +801,7 @@ def reset_warnings(gallery_conf, fname): category=DeprecationWarning ) + reset_warnings(None, None) From 6164e2f18815aa138746e519d60232c6637bcdfc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Richard=20H=C3=B6chenberger?= Date: Mon, 19 Apr 2021 16:51:00 +0200 Subject: [PATCH 9/9] Revert "Filter deprecation warning" This reverts commit 640df6469934826b7910039f44ba69c8943276ab. --- doc/conf.py | 8 -------- 1 file changed, 8 deletions(-) diff --git a/doc/conf.py b/doc/conf.py index b22a0a3fe1f..18422dbcf04 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -791,15 +791,7 @@ def reset_warnings(gallery_conf, fname): 'ignore', message="can't resolve package from", category=ImportWarning) warnings.filterwarnings( 'ignore', message='.*mne-realtime.*', category=DeprecationWarning) - warnings.filterwarnings( - 'ignore', message='.*mne-realtime.*', category=DeprecationWarning) - # https://github.com/mne-tools/mne-python/pull/9311#issuecomment-822042591 - warnings.filterwarnings( - 'ignore', module=r'sphinx\.jinja2glue', - message=".*the old name will be removed.*", - category=DeprecationWarning - ) reset_warnings(None, None)