[FIX] ASR: uncentered scm block covariance + overlap/bound/jump fixes (#110)

Fixes a coordinated cluster of issues in `block_covariance`:

- **Uncentered scm covariance:** for `estimator="scm"` compute the
uncentered per-block second moment `E[x x.T] = B @ B.T / window`
directly, matching `ASR.transform` (`cov = 1/N * X @ X.T`) and the ASR
spec (the calibration covariance is not mean-subtracted). pyriemann's
`scm` mean-centers each block, which diverges the calibration matrix `M`
from the correct value on real data and makes calibrate inconsistent
with transform. Other estimators still route through pyriemann
`covariances`.
- **Loop bound:** the loop bound is now computed against the padded
length. Previously `n_samples` was captured before padding, so padding
added zero iterations and trailing samples were dropped. A window too
large to form any complete block now raises a clear `ValueError`.
- **Overlap convention:** `jump` now uses `round(window * (1 -
overlap))` so `overlap` is the fraction of overlap (higher = more
overlap), matching `clean_windows` and `asr_calibrate`. Previously a
larger overlap produced a larger jump (less actual overlap).
- **Robustness:** `jump` is clamped to `>= 1` so zero-overlap no longer
loops forever, and `window` is cast to `int` so fractional windows don't
break slicing.

### Testing
Adds one regression test per fix. Full `tests/test_asr.py` passes;
`ruff` clean.

---------

Co-authored-by: nbara <10333715+nbara@users.noreply.github.com>

Author: sappelhoff

meegkit/utils/covariances.py

@@ -24,7 +24,7 @@ def block_covariance(data, window=128, overlap=0.5, padding=True, estimator="cov
     window : int
         Window size.
     overlap : float
-        Overlap between successive windows.
+        Fraction of overlap between successive windows (higher = more overlap).
 
     Returns
     -------
@@ -33,19 +33,30 @@ def block_covariance(data, window=128, overlap=0.5, padding=True, estimator="cov
 
     """
     assert 0 <= overlap < 1, "overlap must be < 1"
+    window = int(window)  # window may be fractional
     blocks = []
     n_chans, n_samples = data.shape
     if padding:  # pad data with zeros
-        pad = np.zeros((n_chans, int(window / 2)))
+        pad = np.zeros((n_chans, window // 2))
         data = np.concatenate((pad, data, pad), axis=1)
+        n_samples = data.shape[1]  # bound against the padded length
 
-    jump = int(window * overlap)
+    jump = max(int(round(window * (1 - overlap))), 1)  # >= 1 so overlap=0 advances
     ix = 0
-    while (ix + window < n_samples):
+    while ix + window <= n_samples:
         blocks.append(data[:, ix:ix + window])
-        ix = ix + jump
-
-    return covariances(np.array(blocks), estimator=estimator)
+        ix += jump
+
+    if len(blocks) == 0:
+        raise ValueError(
+            "block_covariance: window is too large for the given data "
+            "(no complete blocks).")
+
+    blocks = np.array(blocks)
+    if estimator == "scm":
+        # uncentered second moment E[x x.T] per block (not mean-subtracted)
+        return blocks @ blocks.transpose(0, 2, 1) / window
+    return covariances(blocks, estimator=estimator)
 
 
 def cov_lags(X, Y, shifts=None):

tests/test_asr.py

@@ -10,6 +10,7 @@
 
 from meegkit.asr import ASR, asr_calibrate, asr_process, clean_windows
 from meegkit.utils.asr import SHAPE_RANGE, fit_eeg_distribution, yulewalk, yulewalk_filter
+from meegkit.utils.covariances import block_covariance
 from meegkit.utils.matrix import sliding_window
 
 # Data files
@@ -359,7 +360,7 @@ def test_asr_calibrate_too_short():
     with pytest.raises(ValueError, match="shorter than one analysis window"):
         asr_calibrate(X_short, 250)
 
-        
+
 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
@@ -451,6 +452,91 @@ def test_fit_eeg_distribution_default_step_sizes():
     assert not np.allclose(sig_old, sig_explicit, rtol=1e-9)
 
 
+def test_block_covariance_uncentered_scm():
+    """scm blocks use the uncentered second moment, not the mean-subtracted cov.
+
+    A large per-channel DC offset makes the two estimates diverge strongly.
+    """
+    W = 50
+    data = rng.standard_normal((4, 250)) + np.arange(1, 5)[:, None] * 100.0
+
+    cov = block_covariance(data, window=W, overlap=0.5, padding=False,
+                           estimator="scm")
+
+    # Reference: same block start indices the function uses.
+    jump = max(int(round(W * (1 - 0.5))), 1)
+    n_samples = data.shape[1]
+    ref, centered = [], []
+    ix = 0
+    while ix + W <= n_samples:
+        B = data[:, ix:ix + W]
+        ref.append(B @ B.T / W)          # uncentered second moment
+        centered.append(np.cov(B, bias=True))  # mean-subtracted
+        ix += jump
+    ref = np.array(ref)
+    centered = np.array(centered)
+
+    assert cov.shape == ref.shape
+    assert np.allclose(cov, ref, rtol=1e-10)
+    # Must NOT be the centered (mean-subtracted) version.
+    assert not np.allclose(cov, centered)
+
+
+def test_block_covariance_overlap_semantics():
+    """Higher overlap yields more blocks."""
+    W = 50
+    data = rng.standard_normal((4, 500))
+    n_hi = block_covariance(data, window=W, overlap=0.8).shape[0]
+    n_lo = block_covariance(data, window=W, overlap=0.2).shape[0]
+    assert n_hi > n_lo
+
+
+def test_block_covariance_padding_bound():
+    """Padding adds blocks (loop bound uses the padded length)."""
+    W = 100
+    data = rng.standard_normal((8, 300))
+    n_pad = block_covariance(data, window=W, overlap=0.5, padding=True).shape[0]
+    n_nopad = block_covariance(data, window=W, overlap=0.5,
+                               padding=False).shape[0]
+    assert n_pad > n_nopad
+
+
+def test_block_covariance_empty_guard():
+    """A window too large for the data raises a clear error."""
+    data = rng.standard_normal((4, 30))
+    with pytest.raises(ValueError, match="too large"):
+        block_covariance(data, window=100, padding=False)
+
+
+def test_block_covariance_exact_one_window():
+    """When n_samples equals window, exactly one block is returned."""
+    data = rng.standard_normal((4, 100))
+    cov = block_covariance(data, window=100, overlap=0.5, padding=False,
+                           estimator="scm")
+
+    assert cov.shape == (1, 4, 4)
+
+
+def test_block_covariance_no_hang_and_float_window():
+    """Zero-overlap does not hang and float windows are accepted."""
+    raw = np.load(os.path.join(THIS_FOLDER, "data", "eeg_raw.npy"))
+    sfreq = 250
+    X_small = raw[:, 5 * sfreq:15 * sfreq]  # 8 chans x 10 s
+
+    # win_overlap=0.0 previously looped forever (jump=0)
+    asr = ASR(sfreq=sfreq, win_overlap=0.0)
+    asr.fit(X_small)
+    M = asr.state_["M"]
+    T = asr.state_["T"]
+    assert np.isfinite(M).all()
+    assert np.isfinite(T).all()
+
+    # Float window must not break slicing (int-cast).
+    cov = block_covariance(X_small, window=100.0, overlap=0.5,
+                           estimator="scm")
+    assert np.isfinite(cov).all()
+
+
 if __name__ == "__main__":
     pytest.main([__file__])
     # test_yulewalk(250, True)

tests/test_pyriemann_imports.py

@@ -45,10 +45,14 @@ def test_block_covariance_uses_pyriemann_covariances():
     from pyriemann.geometry.covariance import covariances
 
     data = 0.1 + 0.01 * (1 + np.arange(24, dtype=float).reshape(3, 8))
+    window = 4
+    overlap = 0.5
+    jump = max(int(round(window * (1 - overlap))), 1)
+    n_samples = data.shape[1]
     expected_blocks = np.array(
-        [data[:, start:start + 4] for start in range(0, 4, 2)]
+        [data[:, start:start + window] for start in range(0, n_samples - window + 1, jump)]
     )
 
-    actual = block_covariance(data, window=4, overlap=0.5, padding=False)
+    actual = block_covariance(data, window=window, overlap=overlap, padding=False)
 
     np.testing.assert_allclose(actual, covariances(expected_blocks, estimator="cov"))