Use direct NumPy MAD in the windowed correlation loop (#197)

* Use direct numpy MAD in windowed correlation loop

* Add changelog entry

* Unify MAD computations behind a single _mad helper

Replace the mix of scipy.stats.median_abs_deviation calls and inline
NumPy MAD implementations in NoisyChannels with one direct-NumPy helper,
pyprep.utils._mad, used everywhere a median absolute deviation is needed
(find_bad_by_nan_flat, find_bad_by_hfnoise, the find_bad_by_correlation
windowed loop, and find_bad_by_PSD's robust_zscore).

_mad is numerically identical to scipy's median_abs_deviation with its
default scale=1 and nan_policy="propagate", so output is bit-identical;
it avoids scipy's per-call input validation and NaN scan, preserving the
windowed-loop speedup from the prior commit.

The NaN-robust spread in find_bad_by_hfnoise stays inline (it centers on
nanmedian, not median) with a comment explaining why. Re-add the _mad
unit test (now also asserting scipy equivalence and NaN propagation),
drop the now-unused scipy import, and update the changelog.

Author: sappelhoff

docs/changelog.rst

@@ -26,6 +26,7 @@ Version 0.8.0 (unreleased)
 
 Changelog
 ~~~~~~~~~
+- Unified all median absolute deviation (MAD) computations in :class:`~pyprep.NoisyChannels` behind a single direct-NumPy helper (``pyprep.utils._mad``), replacing the previous mix of :func:`scipy.stats.median_abs_deviation` calls and inline implementations. This keeps the speedup of the windowed loop in :meth:`~pyprep.NoisyChannels.find_bad_by_correlation` while producing bit-identical output, by `Stefan Appelhoff`_ (:gh:`197`)
 - :class:`~pyprep.NoisyChannels` now detrends the EEG signal in place during initialization, avoiding two transient full-size copies of the recording; :func:`pyprep.removeTrend.removeTrend` gained a ``copy`` keyword argument (default ``True``, preserving the previous behavior) to support this, by `Stefan Appelhoff`_ (:gh:`196`)
 - :meth:`~pyprep.NoisyChannels.find_bad_by_correlation` (and any other method that uses the internally band-pass-filtered signal) is now substantially faster: ``NoisyChannels._get_filtered_data`` filters all channels at once instead of channel-by-channel, with bit-identical output, by `Stefan Appelhoff`_ (:gh:`195`)
 

pyprep/find_noisy_channels.py

@@ -7,11 +7,10 @@
 import numpy as np
 from mne.utils import check_random_state, logger
 from scipy import signal
-from scipy.stats import median_abs_deviation
 
 from pyprep.ransac import find_bad_by_ransac
 from pyprep.removeTrend import removeTrend
-from pyprep.utils import _filter_design, _mat_iqr, _mat_quantile
+from pyprep.utils import _filter_design, _mad, _mat_iqr, _mat_quantile
 
 
 class NoisyChannels:
@@ -396,7 +395,7 @@ def find_bad_by_nan_flat(self, flat_threshold=1e-15):
         nan_channels = self.ch_names_original[nan_channel_mask]
 
         # Detect channels with flat or extremely weak signals
-        flat_by_mad = median_abs_deviation(EEGData, axis=1) < flat_threshold
+        flat_by_mad = _mad(EEGData, axis=1) < flat_threshold
         flat_by_stdev = np.std(EEGData, axis=1) < flat_threshold
         flat_channel_mask = flat_by_mad | flat_by_stdev
         flat_channels = self.ch_names_original[flat_channel_mask]
@@ -485,9 +484,12 @@ def find_bad_by_hfnoise(self, HF_zscore_threshold=5.0):
         # < 50 Hz amplitude for each channel and get robust z-scores of values
         if self.sample_rate > 100:
             noisiness = np.divide(
-                median_abs_deviation(self.EEGData - self.EEGFiltered, axis=1),
-                median_abs_deviation(self.EEGFiltered, axis=1),
+                _mad(self.EEGData - self.EEGFiltered, axis=1),
+                _mad(self.EEGFiltered, axis=1),
             )
+            # NaN-robust z-score: the center uses nanmedian (noisiness may hold
+            # NaNs), so the MAD-style spread is computed inline rather than with
+            # ``_mad``, which centers on the plain median.
             noise_median = np.nanmedian(noisiness)
             noise_sd = np.median(np.abs(noisiness - noise_median)) * MAD_TO_SD
             noise_zscore[self.usable_idx] = (noisiness - noise_median) / noise_sd
@@ -570,7 +572,7 @@ def find_bad_by_correlation(
             channel_amplitudes[w, usable] = _mat_iqr(eeg_raw, axis=1) * IQR_TO_SD
 
             # Check for any channel dropouts (flat signal) within the window
-            eeg_amplitude = median_abs_deviation(eeg_filtered, axis=1)
+            eeg_amplitude = _mad(eeg_filtered, axis=1)
             dropout[w, usable] = eeg_amplitude == 0
 
             # Exclude any dropout chans from further calculations (avoids div-by-zero)
@@ -580,7 +582,7 @@ def find_bad_by_correlation(
             eeg_amplitude = eeg_amplitude[eeg_amplitude > 0]
 
             # Get high-frequency noise ratios for the window
-            high_freq_amplitude = median_abs_deviation(eeg_raw - eeg_filtered, axis=1)
+            high_freq_amplitude = _mad(eeg_raw - eeg_filtered, axis=1)
             noiselevels[w, usable] = high_freq_amplitude / eeg_amplitude
 
             # Get inter-channel correlations for the window
@@ -704,8 +706,7 @@ def find_bad_by_PSD(self, zscore_threshold=3.0, fmin=1.0, fmax=45.0):
         def robust_zscore(values):
             """Compute robust z-scores using MAD."""
             median = np.median(values)
-            mad = np.median(np.abs(values - median))
-            sd = mad * MAD_TO_SD
+            sd = _mad(values) * MAD_TO_SD
             if sd > 0:
                 return (values - median) / sd
             return np.zeros_like(values)

pyprep/utils.py

@@ -474,6 +474,39 @@ def _correlate_arrays(a, b, matlab_strict=False):
         return np.diag(np.corrcoef(a, b)[:n_chan, n_chan:])
 
 
+def _mad(x, axis=None):
+    """Calculate the median absolute deviation from the median (MAD).
+
+    This is a direct NumPy implementation that is numerically identical to
+    :func:`scipy.stats.median_abs_deviation` with its default ``scale=1`` and
+    ``nan_policy="propagate"`` (i.e. any ``NaN`` along ``axis`` propagates to
+    the output). It is used in place of the SciPy function because the latter
+    performs per-call input validation and a full ``NaN`` scan of the input,
+    which is a measurable overhead when the MAD is computed repeatedly inside
+    the windowed loop of
+    :meth:`~pyprep.NoisyChannels.find_bad_by_correlation`.
+
+    Parameters
+    ----------
+    x : np.ndarray
+        A numeric array to summarize.
+    axis : int | None
+        Axis along which the MAD is calculated. If ``None``, the MAD is
+        calculated over the flattened array. Defaults to ``None``.
+
+    Returns
+    -------
+    mad : scalar | np.ndarray
+        If ``axis`` is ``None``, the MAD of the full array as a single value.
+        Otherwise, an :class:`~numpy.ndarray` with the MAD for each slice along
+        the specified axis.
+
+    """
+    x = np.asarray(x)
+    median = np.median(x, axis=axis, keepdims=True)
+    return np.median(np.abs(x - median), axis=axis)
+
+
 def _filter_design(N_order, amp, freq):
     """Create FIR low-pass filter for EEG data using frequency sampling method.
 

tests/test_utils.py

@@ -5,11 +5,13 @@
 
 import numpy as np
 import pytest
+from scipy.stats import median_abs_deviation
 
 from pyprep.utils import (
     _correlate_arrays,
     _eeglab_create_highpass,
     _get_random_subset,
+    _mad,
     _mat_iqr,
     _mat_quantile,
     _mat_round,
@@ -153,3 +155,29 @@ def test_eeglab_create_highpass():
     expected_val = 0.9961
     actual_val = vals[len(vals) // 2]
     assert np.isclose(expected_val, actual_val, atol=0.001)
+
+
+def test_mad():
+    """Test the median absolute deviation from the median (MAD) function."""
+    # Generate test data
+    tst = np.array([[1, 2, 3, 4, 8], [80, 10, 20, 30, 40], [100, 200, 800, 300, 400]])
+    expected = np.asarray([1, 10, 100])
+
+    # Compare output to expected results
+    assert all(np.equal(_mad(tst, axis=1), expected))
+    assert all(np.equal(_mad(tst.T, axis=0), expected))
+    assert _mad(tst) == 28  # Matches robust.mad from statsmodels
+
+    # _mad must stay numerically identical to the scipy function it replaces
+    # (scale=1, nan_policy="propagate"), which is what keeps its output
+    # bit-identical while avoiding scipy's per-call overhead.
+    rng = np.random.default_rng(42)
+    data = rng.standard_normal((8, 200))
+    assert np.array_equal(_mad(data, axis=1), median_abs_deviation(data, axis=1))
+    assert np.array_equal(_mad(data, axis=0), median_abs_deviation(data, axis=0))
+    assert _mad(data) == median_abs_deviation(data, axis=None)
+
+    # A NaN propagates along the reduction axis, matching scipy's default policy
+    data[0, 0] = np.nan
+    assert np.isnan(_mad(data, axis=1)[0])
+    assert not np.isnan(_mad(data, axis=1)[1])