diff --git a/doc/changes/dev/14279.bugfix.rst b/doc/changes/dev/14279.bugfix.rst new file mode 100644 index 00000000000..7617d1b15e6 --- /dev/null +++ b/doc/changes/dev/14279.bugfix.rst @@ -0,0 +1 @@ +Use a fixed normalized Kuiper-index cutoff of 0.3 for automatic CTPS ECG detection in :meth:`mne.preprocessing.ICA.find_bads_ecg`, independent of sampling frequency and trial count, by :newcontrib:`Amirali Moradniaei`. diff --git a/doc/changes/names.inc b/doc/changes/names.inc index 75fe3d0cb18..bbb0adf2be7 100644 --- a/doc/changes/names.inc +++ b/doc/changes/names.inc @@ -20,6 +20,7 @@ .. _Alexandra Corneyllie: https://github.com/Lx37 .. _Alexandre Barachant: https://alexandre.barachant.org .. _Aman Srivastava: https://github.com/aman-coder03 +.. _Amirali Moradniaei: https://github.com/AtomicGlance .. _Ana Radanovic: https://github.com/anaradanovic .. _Andrea Brovelli: https://brovelli.github.io/ .. _Andreas Højlund: https://github.com/ahoejlund diff --git a/mne/channels/_dig_montage_utils.py b/mne/channels/_dig_montage_utils.py index 31bee83648d..e574f8d27ab 100644 --- a/mne/channels/_dig_montage_utils.py +++ b/mne/channels/_dig_montage_utils.py @@ -19,7 +19,9 @@ def _read_dig_montage_egi( "hsp, hpi, elp, point_names, fif must all be None if egi is not None" ) _check_fname(fname, overwrite="read", must_exist=True) - defusedxml = _soft_import("defusedxml", "reading EGI montages") + _soft_import("defusedxml", "reading EGI montages") + import defusedxml.ElementTree + root = defusedxml.ElementTree.parse(fname).getroot() ns = root.tag[root.tag.index("{") : root.tag.index("}") + 1] sensors = root.find(f"{ns}sensorLayout/{ns}sensors") diff --git a/mne/channels/_standard_montage_utils.py b/mne/channels/_standard_montage_utils.py index b0e9aeebfb3..e77b10c78b6 100644 --- a/mne/channels/_standard_montage_utils.py +++ b/mne/channels/_standard_montage_utils.py @@ -423,7 +423,9 @@ def _read_brainvision(fname, head_size): # standard electrode positions: X-axis from T7 to T8, Y-axis from Oz to # Fpz, Z-axis orthogonal from XY-plane through Cz, fit to a sphere if # idealized (when radius=1), specified in millimeters - defusedxml = _soft_import("defusedxml", "reading BrainVision montages") + _soft_import("defusedxml", "reading BrainVision montages") + import defusedxml.ElementTree + root = defusedxml.ElementTree.parse(fname).getroot() ch_names = [s.text for s in root.findall("./Electrode/Name")] theta = [float(s.text) for s in root.findall("./Electrode/Theta")] diff --git a/mne/io/nedf/nedf.py b/mne/io/nedf/nedf.py index 63242a51bee..77815a8f7d6 100644 --- a/mne/io/nedf/nedf.py +++ b/mne/io/nedf/nedf.py @@ -55,7 +55,9 @@ def _parse_nedf_header(header): n_samples : int The number of data samples. """ - defusedxml = _soft_import("defusedxml", "reading NEDF data") + _soft_import("defusedxml", "reading NEDF data") + import defusedxml.ElementTree # ty: ignore[unresolved-import] + info = {} # nedf files have three accelerometer channels sampled at 100Hz followed # by five EEG samples + TTL trigger sampled at 500Hz diff --git a/mne/preprocessing/ica.py b/mne/preprocessing/ica.py index c8aebe7dbc5..f171e874891 100644 --- a/mne/preprocessing/ica.py +++ b/mne/preprocessing/ica.py @@ -4,7 +4,6 @@ # Copyright the MNE-Python contributors. import json -import math import warnings from collections import namedtuple from collections.abc import Sequence @@ -1617,30 +1616,6 @@ def _find_bads_ch( return labels, scores - def _get_ctps_threshold(self, pk_threshold=20): - """Automatically decide the threshold of Kuiper index for CTPS method. - - This function finds the threshold of Kuiper index based on the - threshold of pk. Kuiper statistic that minimizes the difference between - pk and the pk threshold (defaults to 20 :footcite:`DammersEtAl2008`) - is returned. It is assumed that the data are appropriately filtered and - bad data are rejected at least based on peak-to-peak amplitude - when/before running the ICA decomposition on data. - - References - ---------- - .. footbibliography:: - """ - N = self.info["sfreq"] - Vs = np.arange(1, 100) / 100 - C = math.sqrt(N) + 0.155 + 0.24 / math.sqrt(N) - # in formula (13), when k gets large, only k=1 matters for the - # summation. k*V*C thus becomes V*C - Pks = 2 * (4 * (Vs * C) ** 2 - 1) * (np.exp(-2 * (Vs * C) ** 2)) - # NOTE: the threshold of pk is transformed to Pk for comparison - # pk = -log10(Pk) - return Vs[np.argmin(np.abs(Pks - 10 ** (-pk_threshold)))] - @verbose def find_bads_ecg( self, @@ -1714,9 +1689,9 @@ def find_bads_ecg( The ``threshold``, ``method``, and ``measure`` parameters interact in the following ways: - - If ``method='ctps'``, ``threshold`` refers to the significance value - of a Kuiper statistic, and ``threshold='auto'`` will compute the - threshold automatically based on the sampling frequency. + - If ``method='ctps'``, ``threshold`` refers to the maximum normalized + Kuiper index across time, and ``threshold='auto'`` sets the threshold + to 0.3, independent of the sampling frequency and number of trials. - If ``method='correlation'`` and ``measure='correlation'``, ``threshold`` refers to the Pearson correlation value, and ``threshold='auto'`` sets the threshold to 0.9. @@ -1747,9 +1722,6 @@ def find_bads_ecg( ecg = inst.ch_names[idx_ecg] if method == "ctps": - if threshold == "auto": - threshold = self._get_ctps_threshold() - logger.info(f"Using threshold: {threshold:.2f} for CTPS ECG detection") if isinstance(inst, BaseRaw): sources = self.get_sources( create_ecg_epochs( @@ -1771,6 +1743,9 @@ def find_bads_ecg( sources = self.get_sources(inst).get_data(copy=False) else: raise ValueError("With `ctps` only Raw and Epochs input is supported") + if threshold == "auto": + threshold = 0.3 + logger.info(f"Using threshold: {threshold:.2f} for CTPS ECG detection") _, p_vals, _ = ctps(sources) scores = p_vals.max(-1) ecg_idx = np.where(scores >= threshold)[0] diff --git a/mne/preprocessing/tests/test_ica.py b/mne/preprocessing/tests/test_ica.py index d095154936c..1217de4e053 100644 --- a/mne/preprocessing/tests/test_ica.py +++ b/mne/preprocessing/tests/test_ica.py @@ -762,6 +762,24 @@ def short_raw_epochs(): return raw, epochs, epochs_eog +@pytest.mark.parametrize("n_epochs", [2, 3]) +@pytest.mark.parametrize("sfreq", [50, 100]) +def test_ctps_auto_threshold(short_raw_epochs, n_epochs, sfreq): + """Check auto uses the same normalized cutoff for different inputs.""" + _, epochs, _ = short_raw_epochs + epochs = epochs[:n_epochs].copy().resample(sfreq) + ica = _ICA(n_components=2, rng=0) + with _baseline_corrected: + ica.fit(epochs) + with catch_logging(True) as log: + auto_idx, auto_scores = ica.find_bads_ecg(epochs, threshold="auto") + manual_idx, manual_scores = ica.find_bads_ecg(epochs, threshold=0.3) + assert "Using threshold: 0.30" in log.getvalue() + assert_array_equal(auto_scores, manual_scores) + assert_array_equal(auto_idx, manual_idx) + assert_array_equal(sorted(auto_idx), np.flatnonzero(auto_scores >= 0.3)) + + @pytest.mark.slowtest @pytest.mark.parametrize("method", ["picard", "fastica"]) def test_ica_additional(method, tmp_path, short_raw_epochs): @@ -793,8 +811,6 @@ def test_ica_additional(method, tmp_path, short_raw_epochs): ica.fit(raw, np.arange(1, 6)) _assert_ica_attributes(ica, raw.get_data(np.arange(1, 6))) - # check Kuiper index threshold - assert_allclose(ica._get_ctps_threshold(), 0.5) with pytest.raises(TypeError, match="str or numeric"): ica.find_bads_ecg(raw, threshold=None) with pytest.warns(RuntimeWarning, match="is longer than the signal"): @@ -806,7 +822,11 @@ def test_ica_additional(method, tmp_path, short_raw_epochs): ) # check passing a ch_name to find_bads_ecg with pytest.warns(RuntimeWarning, match="longer"): - _, scores_1 = ica.find_bads_ecg(raw, threshold="auto") + auto_idx, scores_1 = ica.find_bads_ecg(raw, threshold="auto") + with pytest.warns(RuntimeWarning, match="longer"): + manual_idx, manual_scores = ica.find_bads_ecg(raw, threshold=0.3) + assert_array_equal(scores_1, manual_scores) + assert_array_equal(auto_idx, manual_idx) with pytest.warns(RuntimeWarning, match="longer"): _, scores_2 = ica.find_bads_ecg(raw, raw.ch_names[1], threshold="auto") assert scores_1[0] != scores_2[0] diff --git a/tutorials/preprocessing/40_artifact_correction_ica.py b/tutorials/preprocessing/40_artifact_correction_ica.py index 3d2fc6f7923..637bfe34887 100644 --- a/tutorials/preprocessing/40_artifact_correction_ica.py +++ b/tutorials/preprocessing/40_artifact_correction_ica.py @@ -444,8 +444,8 @@ # necessary to pass a specific channel name. # `~mne.preprocessing.ICA.find_bads_ecg` also has two options for its # ``method`` parameter: ``'ctps'`` (cross-trial phase statistics -# :footcite:`DammersEtAl2008`) and -# ``'correlation'`` (Pearson correlation between data and ECG channel). +# :footcite:`DammersEtAl2008`) and ``'correlation'`` (Pearson correlation +# between data and ECG channel). ica.exclude = [] # find which ICs match the ECG pattern