BUG: Raise early on non-finite values in PSD (Welch) and ICA.fit (Fix… (mne-tools#13486)

Co-authored-by: Daniel McCloy <[email protected]>
Co-authored-by: Eric Larson <[email protected]>

Author: 3 people

doc/changes/dev/13486.bugfix.rst

@@ -0,0 +1,3 @@
+Improved error message when non-finite values (NaN/Inf) are detected in calls to
+:meth:`inst.compute_psd(method="welch") <mne.io.Raw.compute_psd>` or
+:meth:`ICA.fit() <mne.preprocessing.ICA.fit>`, by :newcontrib:`Emma Zhang`.

doc/changes/names.inc

@@ -79,6 +79,7 @@
 .. _Eduard Ort: https://github.com/eort
 .. _Emily Stephen: https://github.com/emilyps14
 .. _Emma Bailey: https://www.cbs.mpg.de/employees/bailey
+.. _Emma Zhang: https://portfolio-production-ed03.up.railway.app/
 .. _Emmanuel Ferdman: https://github.com/emmanuel-ferdman
 .. _Emrecan Çelik: https://github.com/emrecncelik
 .. _Enrico Varano: https://github.com/enricovara/

doc/sphinxext/related_software.py

@@ -64,10 +64,6 @@
         "Home-page": "https://github.com/mind-inria/picard",
         "Summary": "Preconditioned ICA for Real Data",
     },
-    "mne-features": {
-        "Home-page": "https://mne.tools/mne-features",
-        "Summary": "MNE-Features software for extracting features from multivariate time series",  # noqa: E501
-    },
     "mffpy": {
         "Home-page": "https://github.com/BEL-Public/mffpy",
         "Summary": "Reader and Writer for Philips' MFF file format.",

mne/preprocessing/ica.py

@@ -891,6 +891,9 @@ def _pre_whiten(self, data):
 
     def _fit(self, data, fit_type):
         """Aux function."""
+        if not np.isfinite(data).all():
+            raise ValueError("Input data contains non-finite values (NaN/Inf). ")
+
         random_state = check_random_state(self.random_state)
         n_channels, n_samples = data.shape
         self._compute_pre_whitener(data)

mne/preprocessing/tests/test_ica.py

@@ -1746,3 +1746,30 @@ def test_ica_get_sources_concatenated():
     # get sources
     raw_sources = ica.get_sources(raw_concat)  # but this only has 3 seconds of data
     assert raw_concat.n_times == raw_sources.n_times  # this will fail
+
+
+@pytest.mark.filterwarnings(
+    "ignore:The data has not been high-pass filtered.:RuntimeWarning"
+)
+@pytest.mark.filterwarnings(
+    "ignore:invalid value encountered in subtract:RuntimeWarning"
+)
+def test_ica_rejects_nonfinite():
+    """ICA.fit should fail early on NaN/Inf in the input data."""
+    info = create_info(["Fz", "Cz", "Pz", "Oz"], sfreq=100.0, ch_types="eeg")
+    rng = np.random.RandomState(1)
+    data = rng.standard_normal(size=(4, 1000))
+
+    # Case 1: NaN
+    raw = RawArray(data.copy(), info)
+    raw._data[0, 25] = np.nan
+    ica = ICA(n_components=2, random_state=0, method="fastica", max_iter="auto")
+    with pytest.raises(ValueError, match=r"Input data contains non-finite values"):
+        ica.fit(raw)
+
+    # Case 2: Inf
+    raw = RawArray(data.copy(), info)
+    raw._data[1, 50] = np.inf
+    ica = ICA(n_components=2, random_state=0, method="fastica", max_iter="auto")
+    with pytest.raises(ValueError, match=r"Input data contains non-finite values"):
+        ica.fit(raw)

mne/time_frequency/psd.py

@@ -9,7 +9,7 @@
 from scipy.signal import spectrogram
 
 from ..parallel import parallel_func
-from ..utils import _check_option, _ensure_int, logger, verbose
+from ..utils import _check_option, _ensure_int, logger, verbose, warn
 from ..utils.numerics import _mask_to_onsets_offsets
 
 
@@ -200,6 +200,42 @@ def psd_array_welch(
     del freq_mask
     freqs = freqs[freq_sl]
 
+    step = max(int(n_per_seg) - int(n_overlap), 1)
+    if n_times >= n_per_seg:
+        n_segments = 1 + (n_times - n_per_seg) // step
+        analyzed_end = step * (n_segments - 1) + n_per_seg
+    else:
+        n_segments = 0
+        analyzed_end = 0
+
+    nan_mask_full = np.isnan(x)
+    nan_present = nan_mask_full.any()
+    if nan_present:
+        good_mask_full = ~nan_mask_full
+        aligned_nan = np.allclose(good_mask_full, good_mask_full[[0]], equal_nan=True)
+    else:
+        aligned_nan = False
+
+    if analyzed_end > 0:
+        # Inf always counts as non-finite per-channel
+        nonfinite_mask = np.isinf(x[..., :analyzed_end])
+        # NaNs count per-channel only if NOT aligned (i.e., not annotations)
+        if nan_present and not aligned_nan:
+            nonfinite_mask |= nan_mask_full[..., :analyzed_end]
+        bad_ch = nonfinite_mask.any(axis=-1)
+    else:
+        bad_ch = np.zeros(x.shape[0], dtype=bool)
+
+    if bad_ch.any():
+        warn(
+            "Non-finite values (NaN/Inf) detected in some channels; PSD for "
+            "those channels will be NaN.",
+        )
+        # avoid downstream NumPy warnings by zeroing bad channels;
+        # will overwrite their PSD rows with NaN at the end
+        x = x.copy()
+        x[bad_ch] = 0.0
+
     # Parallelize across first N-1 dimensions
     logger.debug(
         f"Spectogram using {n_fft}-point FFT on {n_per_seg} samples with "
@@ -217,11 +253,9 @@ def psd_array_welch(
         window=window,
         mode=mode,
     )
-    if np.any(np.isnan(x)):
-        good_mask = ~np.isnan(x)
-        # NaNs originate from annot, so must match for all channels. Note that we CANNOT
-        # use np.testing.assert_allclose() here; it is strict about shapes/broadcasting
-        assert np.allclose(good_mask, good_mask[[0]], equal_nan=True)
+    if nan_present and aligned_nan:
+        # Aligned NaNs across channels → treat as bad annotations.
+        good_mask = ~nan_mask_full
         t_onsets, t_offsets = _mask_to_onsets_offsets(good_mask[0])
         x_splits = [x[..., t_ons:t_off] for t_ons, t_off in zip(t_onsets, t_offsets)]
         # weights reflect the number of samples used from each span. For spans longer
@@ -257,6 +291,12 @@ def func(*args, **kwargs):
                 return _func(*args, **kwargs)
 
     else:
+        # Either no NaNs, or NaNs are not aligned across channels.
+        if nan_present and not aligned_nan:
+            logger.info(
+                "NaN masks are not aligned across channels; treating NaNs as "
+                "per-channel contamination."
+            )
         x_splits = [arr for arr in np.array_split(x, n_jobs) if arr.size != 0]
         agg_func = np.concatenate
         func = _func
@@ -268,5 +308,9 @@ def func(*args, **kwargs):
     shape = dshape + (len(freqs),)
     if average is None:
         shape = shape + (-1,)
+
+    if bad_ch.any():
+        psds[bad_ch] = np.nan
+
     psds.shape = shape
     return psds, freqs

mne/time_frequency/tests/test_psd.py

@@ -226,3 +226,36 @@ def test_psd_array_welch_n_jobs():
     data = np.zeros((1, 2048))
     psd_array_welch(data, 1024, n_jobs=1)
     psd_array_welch(data, 1024, n_jobs=2)
+
+
+def test_psd_nan_in_data():
+    """psd_array_welch should fail if +Inf lies inside analyzed samples."""
+    n_samples, n_fft, n_overlap = 2048, 256, 128
+    rng = np.random.RandomState(0)
+    x = rng.standard_normal(size=(2, n_samples))
+    # Put +Inf inside the series; this falls within Welch windows
+    x[0, 800] = np.inf  # Channel 0 has Inf → bad channel
+    with pytest.warns(RuntimeWarning, match="Non-finite values"):
+        psds, freqs = psd_array_welch(x, float(n_fft), n_fft=n_fft, n_overlap=n_overlap)
+
+    # Channel 0 is contaminated → NaN PSD
+    assert np.isnan(psds[0]).all()
+
+    # Channel 1 is clean → has finite PSD values
+    assert np.isfinite(psds[1]).any()
+
+
+def test_psd_misaligned_nan_across_channels():
+    """If NaNs are present but masks are NOT aligned across channels."""
+    n_samples, n_fft, n_overlap = 2048, 256, 128
+    rng = np.random.RandomState(42)
+    x = rng.standard_normal(size=(2, n_samples))
+    # NaN only in ch0; ch1 has no NaN => masks not aligned -> should raise
+    x[0, 500] = np.nan
+    with pytest.warns(RuntimeWarning, match="Non-finite values"):
+        psds, freqs = psd_array_welch(x, float(n_fft), n_fft=n_fft, n_overlap=n_overlap)
+    # Bad channel gets NaN PSD
+    assert np.isnan(psds[0]).all()
+
+    # Good channel retains finite values
+    assert np.isfinite(psds[1]).any()