diff --git a/meegkit/utils/covariances.py b/meegkit/utils/covariances.py index 71ca5b1f..bf31093a 100644 --- a/meegkit/utils/covariances.py +++ b/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): diff --git a/tests/test_asr.py b/tests/test_asr.py index d17db8c1..635a3055 100644 --- a/tests/test_asr.py +++ b/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) diff --git a/tests/test_pyriemann_imports.py b/tests/test_pyriemann_imports.py index 9b82e089..eade9711 100644 --- a/tests/test_pyriemann_imports.py +++ b/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"))