From 389a1844b022f0a13c631062ab9dfac96d5ce773 Mon Sep 17 00:00:00 2001 From: Stefan Appelhoff Date: Mon, 13 Jul 2026 11:23:44 +0200 Subject: [PATCH] [MAINT] ASR: defensive hardening and docstring fixes - asr_process: read state by key instead of relying on dict order - fit_eeg_distribution: nudge gridsearch arange endpoint (mirror cols grid) - clean_windows: build window offsets from non-truncated N (phase drift) - ASR: expose max_bad_chans as a constructor parameter (default 0.3) - docs: fix clean_windows win_len default (0.5); document the mask3 term --- meegkit/asr.py | 18 ++++++++++---- meegkit/utils/asr.py | 4 ++- tests/test_asr.py | 58 ++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 74 insertions(+), 6 deletions(-) diff --git a/meegkit/asr.py b/meegkit/asr.py index 75ea41aa..be90deab 100755 --- a/meegkit/asr.py +++ b/meegkit/asr.py @@ -44,6 +44,10 @@ class ASR: Window overlap fraction. The fraction of two successive windows that overlaps. Higher overlap ensures that fewer artifact portions are going to be missed, but is slower (default=0.66). + max_bad_chans : float + The maximum number or fraction of bad channels that a retained window + may still contain (more than this and it is removed) during + calibration window selection (default=0.3). max_dropout_fraction : float Maximum fraction of windows that can be subject to signal dropouts (e.g., sensor unplugged), used for threshold estimation (default=0.1). @@ -94,7 +98,8 @@ class ASR: """ def __init__(self, *, sfreq=250, cutoff=5, blocksize=100, win_len=0.5, - win_overlap=0.66, max_dropout_fraction=0.1, + win_overlap=0.66, max_bad_chans=0.3, + max_dropout_fraction=0.1, min_clean_fraction=0.25, method="euclid", memory=None, estimator="scm", **kwargs): @@ -104,7 +109,7 @@ def __init__(self, *, sfreq=250, cutoff=5, blocksize=100, win_len=0.5, self.win_overlap = win_overlap self.max_dropout_fraction = max_dropout_fraction self.min_clean_fraction = min_clean_fraction - self.max_bad_chans = 0.3 + self.max_bad_chans = max_bad_chans self.method = method if memory is None: self.memory = int(2 * sfreq) # smoothing window for covariances @@ -290,7 +295,7 @@ def clean_windows(X, sfreq, max_bad_chans=0.2, zthresholds=[-3.5, 5], Window length that is used to check the data for artifact content. This is ideally as long as the expected time scale of the artifacts but not shorter than half a cycle of the high-pass filter that was - used. Default: 1. + used. Default: 0.5. win_overlap : float Window overlap fraction. The fraction of two successive windows that overlaps. Higher overlap ensures that fewer artifact portions are @@ -330,7 +335,8 @@ def clean_windows(X, sfreq, max_bad_chans=0.2, zthresholds=[-3.5, 5], # set data indices [nc, ns] = X.shape N = int(win_len * sfreq) - offsets = np.round(np.arange(0, ns - N, (N * (1 - win_overlap)))) + N_raw = win_len * sfreq # non-truncated N, avoids step-size phase drift + offsets = np.round(np.arange(0, ns - N, N_raw * (1 - win_overlap))) offsets = offsets.astype(int) logging.debug("[ASR] Determining channel-wise rejection thresholds") @@ -360,6 +366,8 @@ def clean_windows(X, sfreq, max_bad_chans=0.2, zthresholds=[-3.5, 5], if np.min(zthresholds) < 0: mask2 = (swz[1 + int(max_bad_chans - 1), :] < np.min(zthresholds)) + # extra meegkit-specific criterion: drop windows whose across-channel + # z-scores are nearly flat (very low MAD or std) bad_by_mad = mad(wz, c=1, axis=0) < .1 bad_by_std = np.std(wz, axis=0) < .1 mask3 = np.logical_or(bad_by_mad, bad_by_std) @@ -558,7 +566,7 @@ def asr_process(X, X_filt, state, cov=None, detrend=False, method="riemann", Output ASR parameters. """ - M, T, R = state.values() + M, T, R = state["M"], state["T"], state["R"] [nc, ns] = X.shape if cov is None: diff --git a/meegkit/utils/asr.py b/meegkit/utils/asr.py index 3c5ced93..0cc11b66 100755 --- a/meegkit/utils/asr.py +++ b/meegkit/utils/asr.py @@ -114,7 +114,9 @@ def fit_eeg_distribution(X, min_clean_fraction=0.25, max_dropout_fraction=0.1, opt_lu = np.inf opt_bounds = np.inf opt_beta = np.inf - gridsearch = np.round(n * np.arange(max_width, min_width, -step_sizes[1])) + # nudge the stop so min_width isn't dropped by float error (as for cols) + gridsearch = np.round( + n * np.arange(max_width, min_width - step_sizes[1] * 1e-9, -step_sizes[1])) for m in gridsearch.astype(int): mcurr = m - 1 nbins = int(np.round(3 * np.log2(1 + m / 2))) diff --git a/tests/test_asr.py b/tests/test_asr.py index 811bcb38..6275857f 100644 --- a/tests/test_asr.py +++ b/tests/test_asr.py @@ -276,6 +276,64 @@ def test_asr_class(method, reref, show=False): ASR(sfreq=125) ASR(Sfreq=150) +def test_asr_max_bad_chans_param(): + """max_bad_chans is exposed on ASR and defaults to 0.3.""" + assert ASR().max_bad_chans == 0.3 + assert ASR(max_bad_chans=0.2).max_bad_chans == 0.2 + + +def test_asr_process_state_key_order(): + """asr_process reads state by key, not by dict insertion order. + + M and T are distinct and the keep mask is non-trivial, so a positional + unpack of a reordered dict would change the output. + """ + nc, ns = 4, 8 + X = np.arange(nc * ns, dtype=float).reshape(nc, ns) + 1.0 + X_filt = X.copy() + # T's large column-2 norm keeps component 2 while M=eye rejects it, so the + # two dict orderings give different keep masks (hence different R) under a + # positional unpack. + cov = np.diag([1.0, 1.0, 1000.0, 1000.0]) + T = np.eye(nc) + T[2, 2] = 40.0 # column-2 sum of squares = 1600 > 1000 + M = np.eye(nc) + + state_ordered = dict(M=M, T=T, R=None) + out_ordered, _ = asr_process( + X, X_filt, state_ordered, cov=cov.copy(), method="euclid") + + # Same content, different insertion order + state_reordered = dict(T=T, M=M, R=None) + out_reordered, _ = asr_process( + X, X_filt, state_reordered, cov=cov.copy(), method="euclid") + + np.testing.assert_allclose(out_ordered, out_reordered) + + +def test_clean_windows_offset_phase_drift(): + """Offsets use the non-truncated win_len*sfreq to avoid phase drift. + + Diverges from truncated-N spacing only when win_len*sfreq is non-integer. + """ + sfreq = 251 + win_len = 0.5 # win_len * sfreq = 125.5, non-integer + win_overlap = 0.66 + ns = 5000 + + N = int(win_len * sfreq) + N_raw = win_len * sfreq + offsets_raw = np.round( + np.arange(0, ns - N, N_raw * (1 - win_overlap))).astype(int) + offsets_truncated = np.round( + np.arange(0, ns - N, N * (1 - win_overlap))).astype(int) + + assert len(offsets_raw) == len(offsets_truncated) + # The two stepping schemes diverge (truncated-N accumulates drift) + assert not np.array_equal(offsets_raw, offsets_truncated) + assert abs(int(offsets_raw[-1]) - int(offsets_truncated[-1])) >= 10 + + if __name__ == "__main__": pytest.main([__file__]) # test_yulewalk(250, True)