[MAINT] ASR: defensive hardening and docstring fixes

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:

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)))

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)