diff --git a/docs/estimators.md b/docs/estimators.md index e4deb13..0df95c3 100644 --- a/docs/estimators.md +++ b/docs/estimators.md @@ -60,13 +60,61 @@ model: sizes, but available. - **`random_forest`**, **`gradient_boosting`**: tree ensembles for the tabular behaviour block. +- **`riemann`**: logistic regression in the log-Euclidean tangent space of EEG + channel covariances. Needs covariance features, not band power. Detailed in + its own section below. -## Two stronger models, and why they are not drop-ins +## Riemannian tangent-space classification for EEG (implemented) -Both were suggested and both are genuinely better *for the right setup*. Neither -fits the flat-feature, subject-grouped-CV pipeline without a dedicated path, and -bolting them in naively would produce numbers that look principled and are not. -Recorded here so the decision is deliberate, not forgotten. +`riemann` is a real estimator now, not a base-learner swap: it is logistic +regression in the **log-Euclidean tangent space** of the EEG channel-covariance +manifold. Use it like any other learner, but feed it covariance features: + +```python +from behavioral_decoding.io.eeg import EEGLoader + +# Covariance mode is mutually exclusive with band power / ERP. +block = EEGLoader( + include_bandpower=False, include_erp=False, include_covariance=True +).from_arrays(epochs, sfreq, times, subject_ids, stimulus_ids) +``` + +```yaml +model: + base_learner: + eeg: riemann # expects covariance features, not band power +``` + +How it stays correct inside the existing pipeline: + +``` +BaggingClassifier + └── Pipeline(RiemannianTangentSpace -> StandardScaler -> SMOTE -> logistic) +``` + +The tangent projection goes *first*, so scaling and SMOTE only ever touch the +flat Euclidean tangent vectors, never the raw covariance entries. The reference +point (the log-Euclidean mean) is computed in `fit` on training data only, inside +each fold and each bag, so it is leakage-safe. Rank-deficient covariances from +short epochs are regularised toward a scaled identity before the matrix log, and +`max_features` is forced to 1.0 because a column subset of a flattened covariance +is not a covariance. + +Two honest caveats: + +- **On clean data the tangent map does not beat a plain logistic**, because the + raw covariance entries are already linearly separable there. Its advantage + shows on ill-conditioned real EEG, which is the regime the method was built + for. The tests assert the path recovers covariance structure *above chance*, + not that it beats logistic on synthetic data, which would be cherry-picking. +- This is the **log-Euclidean** metric (closed-form, scipy-only). `pyriemann`'s + affine-invariant metric, which iterates to a geometric mean and whitens by it, + is a further upgrade; swapping the reference-point computation for it is the + natural next step if a real EEG dataset warrants it. + +## The other suggested model, and why it is not a drop-in + +Recorded so the decision is deliberate, not forgotten. ### Hierarchical / mixed-effects model (fMRI) @@ -92,34 +140,16 @@ random-effect term and predict from fixed effects only. Wrap it to expose `predict_proba`, and report it next to elastic-net rather than replacing it. Do not expect a cross-subject gain. -### Riemannian methods (EEG) - -For EEG built from **channel covariance / connectivity** features, Riemannian -tangent-space classification is a real upgrade over band power: it respects the -geometry of the space of covariance matrices instead of treating their entries -as independent numbers. But two steps in this pipeline are wrong for covariance -features: - -- **StandardScaler** z-scores each covariance entry independently, which - destroys the positive-definite structure the Riemannian method depends on. -- **SMOTE** interpolates along straight lines between covariance vectors. The - geodesic between two covariance matrices is *not* a straight line in entry - space, so the synthetic minority points are off-manifold. - -A correct Riemannian EEG arm therefore needs its own pipeline: emit per-trial -covariance matrices (not band power), project to the tangent space at the -Riemannian mean of the *training* fold, and only then standardise / resample / -classify in that tangent space. That is a separate feature family and a separate -pipeline, best added as a distinct `riemann` path (with `pyriemann`, or a -scipy-`logm` tangent map for a dependency-free version) rather than a base -learner slotted into the existing flat pipeline. It is a good next step; it is -not a one-line default change, and pretending it were would corrupt the manifold -structure it exists to exploit. +(The other model previously in this section, Riemannian tangent-space +classification for EEG, is now implemented; see the section above.) ## What did change `elasticnet` and `linear_svm` were added to the learner factory and the -dense-block defaults moved from plain logistic to elastic-net. Everything else -about the pipeline (bagging, in-fold SMOTE, out-of-fold weighting, subject- -grouped CV) is unchanged, so these are honest swaps of the base estimator, not a -change to how anything is evaluated. +dense-block defaults moved from plain logistic to elastic-net. The `riemann` +path was then added for EEG: a covariance feature family in the loader and a +log-Euclidean tangent-space transformer prepended to the pipeline. Everything +else about the pipeline (bagging, in-fold SMOTE, out-of-fold weighting, subject- +grouped CV) is unchanged. The dense-block swaps are honest swaps of the base +estimator; the Riemannian path adds a manifold-correct front-end without +changing how anything is evaluated. diff --git a/src/behavioral_decoding/balance/smote.py b/src/behavioral_decoding/balance/smote.py index 2ad64ce..c9c4022 100644 --- a/src/behavioral_decoding/balance/smote.py +++ b/src/behavioral_decoding/balance/smote.py @@ -243,6 +243,7 @@ def make_balanced_pipeline( k_neighbors: int = 5, scaler: bool = True, random_state: int = 0, + pre_steps=None, ): """Wrap ``estimator`` in a leakage-safe scale -> resample -> fit pipeline. @@ -253,10 +254,19 @@ def make_balanced_pipeline( The resampling step is an :class:`AdaptiveOverSampler`, which resolves ``k`` against whatever slice of data the fold or bag actually hands it. + ``pre_steps`` is an optional list of ``(name, transformer)`` inserted + *before* the scaler. It exists for feature spaces that must be mapped into a + flat Euclidean space before scaling and SMOTE are valid, the Riemannian + tangent projection being the motivating case: scaling and straight-line + interpolation are only meaningful once the covariance matrices have been + projected to the tangent space. Because the transformer is inside the + pipeline, it is fitted on training data only, inside each CV fold and each + bag. + Falls back to a plain scikit-learn pipeline (scaler + estimator, no resampling) when imbalanced-learn is absent, and logs that it did. """ - steps = [] + steps = list(pre_steps) if pre_steps else [] if scaler: from sklearn.preprocessing import StandardScaler diff --git a/src/behavioral_decoding/io/eeg.py b/src/behavioral_decoding/io/eeg.py index 2d22a43..57c5085 100644 --- a/src/behavioral_decoding/io/eeg.py +++ b/src/behavioral_decoding/io/eeg.py @@ -102,8 +102,55 @@ def erp_windows( return np.concatenate(out, axis=1), names +def channel_covariance(epochs: np.ndarray) -> Tuple[np.ndarray, List[str]]: + """Per-epoch inter-channel covariance, flattened for the Riemannian path. + + Parameters + ---------- + epochs: + Array of shape ``(n_epochs, n_channels, n_times)``. + + Returns + ------- + (features, names) + ``features`` has shape ``(n_epochs, n_channels*(n_channels+1)/2)``, the + flattened upper triangle of each covariance with the same ``sqrt(2)`` + off-diagonal convention the tangent-space transformer expects (they share + :func:`~behavioral_decoding.models.riemann.flatten_spd`, so the loader's + output and the transformer's input cannot drift apart). + + These features are meant to feed the ``riemann`` estimator, which projects + them to the log-Euclidean tangent space before any scaling or resampling. + Handing raw covariance entries to a plain logistic works but discards the + manifold structure; see ``docs/estimators.md``. + """ + from ..models.riemann import flatten_spd + + epochs = np.asarray(epochs, dtype=float) + n_epochs, n_channels, n_times = epochs.shape + if n_times < 2: + raise ValueError("covariance needs at least 2 samples per epoch") + mats = np.empty((n_epochs, n_channels, n_channels), dtype=float) + for i in range(n_epochs): + mats[i] = np.cov(epochs[i]) + features = flatten_spd(mats) + rows, cols = np.triu_indices(n_channels) + names = [f"cov_ch{r:02d}_ch{c:02d}" for r, c in zip(rows, cols)] + return features, names + + class EEGLoader(BaseLoader): - """Turn epoched EEG into trial-by-feature rows.""" + """Turn epoched EEG into trial-by-feature rows. + + Two feature regimes, and they do not mix: + + - **band power and ERP windows** (the default): flat spectral/temporal + features for the dense-block estimators (logistic, elastic-net). + - **covariance** (``include_covariance=True``): per-trial inter-channel + covariance for the ``riemann`` estimator. This is mutually exclusive with + the other two, because the Riemannian tangent map needs the covariance + matrix intact, not concatenated with unrelated columns. + """ name = EEG @@ -113,11 +160,20 @@ def __init__( windows: Optional[Sequence[Tuple[str, float, float]]] = None, include_bandpower: bool = True, include_erp: bool = True, + include_covariance: bool = False, ) -> None: + if include_covariance and (include_bandpower or include_erp): + raise ValueError( + "covariance features are mutually exclusive with band-power/ERP: " + "the Riemannian tangent map needs the covariance matrix intact, " + "not concatenated with other columns. Set include_bandpower=False " + "and include_erp=False when include_covariance=True." + ) self.bands = dict(bands) if bands is not None else dict(DEFAULT_BANDS) self.windows = tuple(windows) if windows is not None else tuple(DEFAULT_ERP_WINDOWS) self.include_bandpower = include_bandpower self.include_erp = include_erp + self.include_covariance = include_covariance def from_arrays( self, @@ -134,6 +190,24 @@ def from_arrays( raise ValueError( f"expected (n_epochs, n_channels, n_times), got shape {epochs.shape}" ) + + if self.include_covariance: + feats, names = channel_covariance(epochs) + return ModalityBlock( + name=self.name, + X=feats, + subject_ids=np.asarray(subject_ids), + stimulus_ids=np.asarray(stimulus_ids), + feature_names=names, + provenance=self._provenance( + source=source, + sfreq=sfreq, + feature_family="covariance", + n_channels=epochs.shape[1], + note="flattened inter-channel covariance for the riemann estimator", + ), + ) + parts: List[np.ndarray] = [] names: List[str] = [] if self.include_bandpower: diff --git a/src/behavioral_decoding/models/modality_models.py b/src/behavioral_decoding/models/modality_models.py index 035a158..ef730e3 100644 --- a/src/behavioral_decoding/models/modality_models.py +++ b/src/behavioral_decoding/models/modality_models.py @@ -187,20 +187,42 @@ def build_modality_model( if kind == "gradient_boosting": class_weight = None - estimator = make_base_learner(kind, class_weight=class_weight, seed=seed, **learner_kwargs) + # The Riemannian path is a logistic regression in the log-Euclidean tangent + # space: a tangent-projection transformer is prepended, and the classifier + # itself is plain logistic. It cannot subsample features per bag, because a + # column subset of a flattened covariance triangle is no longer a valid + # covariance, so max_features is forced to 1.0. + pre_steps = None + effective_max_features = max_features + if kind == "riemann": + from .riemann import RiemannianTangentSpace + + pre_steps = [("tangent", RiemannianTangentSpace())] + estimator = make_base_learner("logistic", class_weight=class_weight, seed=seed) + if max_features != 1.0: + logger.warning( + "riemann[%s]: max_features forced to 1.0; a column subset of a " + "flattened covariance is not a covariance", + modality, + ) + effective_max_features = 1.0 + else: + estimator = make_base_learner(kind, class_weight=class_weight, seed=seed, **learner_kwargs) + pipeline = make_balanced_pipeline( estimator, sampler=str(chosen_sampler), k_neighbors=int(chosen_k), scaler=True, random_state=seed, + pre_steps=pre_steps, ) model = BaggingClassifier( estimator=pipeline, n_estimators=bags, max_samples=max_samples, - max_features=max_features, + max_features=effective_max_features, bootstrap=True, bootstrap_features=False, oob_score=False, diff --git a/src/behavioral_decoding/models/riemann.py b/src/behavioral_decoding/models/riemann.py new file mode 100644 index 0000000..fd966a7 --- /dev/null +++ b/src/behavioral_decoding/models/riemann.py @@ -0,0 +1,162 @@ +"""Riemannian tangent-space mapping for EEG covariance features. + +Why this exists, in one paragraph: EEG carries a lot of its discriminative +signal in the *covariance* between channels, and the set of covariance matrices +is a curved manifold (symmetric positive-definite matrices), not a flat vector +space. Treating the entries of a covariance matrix as independent numbers, then +z-scoring and SMOTE-interpolating them, throws away that geometry and can even +produce non-covariance "matrices". The fix that the EEG literature settled on is +to project each covariance to a **tangent space** (a flat space that locally +approximates the manifold) and classify there. In the tangent space the vectors +are genuinely Euclidean, so the framework's scaler, SMOTE, and logistic +regression are all valid again. + +This module implements the **log-Euclidean** tangent mapping, which is +closed-form, robust, and needs only scipy. It is a legitimate Riemannian metric +on SPD matrices. The affine-invariant metric that ``pyriemann`` uses is a +further refinement (it iterates to a geometric mean and whitens by it); this +gets most of the benefit without the dependency or the iteration. See +``docs/estimators.md``. + +The transformer is designed to sit at the front of the standard pipeline: + + Pipeline(RiemannianTangentSpace -> StandardScaler -> SMOTE -> logistic) + +so everything downstream of the tangent map operates on flat Euclidean vectors, +and the map itself is fitted on training data only (inside each CV fold and each +bag), which keeps it leakage-safe. +""" + +from __future__ import annotations + +import numpy as np +from sklearn.base import BaseEstimator, TransformerMixin + +from ..utils.logging import get_logger + +logger = get_logger(__name__) + + +def n_channels_from_flat(n_features: int) -> int: + """Recover k from the length of a flattened k-by-k upper triangle. + + A symmetric k-by-k matrix has ``k(k+1)/2`` unique entries (diagonal + included). Inverting: ``k = (-1 + sqrt(1 + 8n)) / 2``. Raises if ``n`` is not + a triangular number, which means the features are not a flattened covariance. + """ + k = (-1 + np.sqrt(1 + 8 * n_features)) / 2 + k_int = int(round(k)) + if k_int * (k_int + 1) // 2 != n_features: + raise ValueError( + f"{n_features} features is not k(k+1)/2 for any integer k, so these are not " + "flattened covariance upper-triangles. RiemannianTangentSpace needs " + "covariance features (see EEGLoader covariance mode)." + ) + return k_int + + +def _upper_triangle_indices(k: int): + return np.triu_indices(k) + + +def flatten_spd(matrices: np.ndarray) -> np.ndarray: + """Flatten ``(n, k, k)`` symmetric matrices to ``(n, k(k+1)/2)`` rows. + + Off-diagonal entries are scaled by ``sqrt(2)`` so that the Euclidean norm of + the flattened vector equals the Frobenius norm of the matrix. This is the + standard isometric vectorisation; it matters because the downstream scaler + and classifier work in Euclidean distance. + """ + matrices = np.asarray(matrices, dtype=float) + n, k, _ = matrices.shape + rows, cols = _upper_triangle_indices(k) + scale = np.where(rows == cols, 1.0, np.sqrt(2.0)) + return matrices[:, rows, cols] * scale + + +def unflatten_spd(rows: np.ndarray, k: int) -> np.ndarray: + """Inverse of :func:`flatten_spd`: ``(n, k(k+1)/2)`` rows to ``(n, k, k)``.""" + rows = np.asarray(rows, dtype=float) + n = rows.shape[0] + ri, ci = _upper_triangle_indices(k) + scale = np.where(ri == ci, 1.0, np.sqrt(2.0)) + out = np.zeros((n, k, k), dtype=float) + vals = rows / scale + out[:, ri, ci] = vals + out[:, ci, ri] = vals + return out + + +def _regularise(cov: np.ndarray, shrinkage: float) -> np.ndarray: + """Shrink a covariance toward a scaled identity so it is positive-definite. + + Covariances estimated from short epochs are often rank-deficient, and the + matrix logarithm is only defined for positive-definite input. Shrinkage + ``(1 - s) C + s * (tr C / k) I`` guarantees positive-definiteness for any + ``s > 0`` while barely moving a well-conditioned matrix. + """ + k = cov.shape[0] + mu = np.trace(cov) / k + return (1.0 - shrinkage) * cov + shrinkage * mu * np.eye(k) + + +def _symmetric_logm(cov: np.ndarray) -> np.ndarray: + """Matrix log of a symmetric positive-definite matrix, via eigdecomposition. + + Faster and more stable than a general ``scipy.linalg.logm`` here, because the + input is symmetric: eigenvalues are real and positive after regularisation. + """ + vals, vecs = np.linalg.eigh(cov) + vals = np.clip(vals, 1e-12, None) + return (vecs * np.log(vals)) @ vecs.T + + +class RiemannianTangentSpace(BaseEstimator, TransformerMixin): + """Project flattened EEG covariance features into the log-Euclidean tangent space. + + Input rows are flattened covariance upper-triangles (as produced by the + EEG loader's covariance mode). Output rows are flattened tangent vectors, + centred at the training-set log-Euclidean mean. + + Parameters + ---------- + shrinkage: + Regularisation toward a scaled identity, in ``[0, 1)``. Small but + non-zero so rank-deficient covariances still admit a logarithm. + center: + Subtract the training-set mean log-covariance (the log-Euclidean mean's + logarithm). Improves conditioning and is leakage-safe because the mean + is computed in ``fit`` on training data only. + """ + + def __init__(self, shrinkage: float = 1e-3, center: bool = True) -> None: + self.shrinkage = shrinkage + self.center = center + + def _log_covariances(self, X: np.ndarray) -> np.ndarray: + k = n_channels_from_flat(X.shape[1]) + mats = unflatten_spd(X, k) + logs = np.empty_like(mats) + for i in range(mats.shape[0]): + logs[i] = _symmetric_logm(_regularise(mats[i], self.shrinkage)) + return logs + + def fit(self, X: np.ndarray, y=None) -> RiemannianTangentSpace: + X = np.asarray(X, dtype=float) + self.n_channels_ = n_channels_from_flat(X.shape[1]) + logs = self._log_covariances(X) + # Log-Euclidean mean's logarithm is just the arithmetic mean of the + # per-trial log-covariances. This is the reference point we linearise at. + self.mean_log_ = logs.mean(axis=0) if self.center else np.zeros_like(logs[0]) + return self + + def transform(self, X: np.ndarray) -> np.ndarray: + if not hasattr(self, "mean_log_"): + raise RuntimeError("RiemannianTangentSpace is not fitted") + X = np.asarray(X, dtype=float) + logs = self._log_covariances(X) + centred = logs - self.mean_log_[None, :, :] + return flatten_spd(centred) + + def _more_tags(self): # pragma: no cover - sklearn metadata + return {"requires_positive_X": False, "stateless": False} diff --git a/tests/test_eeg_covariance.py b/tests/test_eeg_covariance.py new file mode 100644 index 0000000..e88219a --- /dev/null +++ b/tests/test_eeg_covariance.py @@ -0,0 +1,89 @@ +"""EEG covariance feature family and its end-to-end path into the riemann model. + +The covariance mode produces the flattened inter-channel covariance that the +Riemannian tangent-space estimator consumes. These tests check the shape, the +naming, the mutual-exclusivity guard, and that the loader's output round-trips +through the transformer's inverse (the two share one vectorisation convention, +so this guards against them drifting apart). +""" + +from __future__ import annotations + +import numpy as np +import pytest + +from behavioral_decoding.io.eeg import EEGLoader, channel_covariance +from behavioral_decoding.models.riemann import n_channels_from_flat, unflatten_spd + + +def _epochs(n=12, k=8, t=128, seed=0): + return np.random.default_rng(seed).normal(size=(n, k, t)) + + +def test_channel_covariance_shape_and_names(): + feats, names = channel_covariance(_epochs(n=10, k=8)) + assert feats.shape == (10, 8 * 9 // 2) # 36 + assert names[0] == "cov_ch00_ch00" + assert len(names) == feats.shape[1] + assert n_channels_from_flat(feats.shape[1]) == 8 + + +def test_covariance_needs_enough_samples(): + with pytest.raises(ValueError, match="at least 2 samples"): + channel_covariance(np.zeros((3, 4, 1))) + + +def test_loader_covariance_block_is_labelled(): + block = EEGLoader( + include_bandpower=False, include_erp=False, include_covariance=True + ).from_arrays( + _epochs(), sfreq=128.0, times=np.linspace(0, 1, 128), + subject_ids=["s1"] * 12, stimulus_ids=list(range(12)), + ) + assert block.provenance["feature_family"] == "covariance" + assert block.n_features == 8 * 9 // 2 + + +def test_loader_output_round_trips_through_transformer_inverse(): + """The loader and the transformer must share one vectorisation convention.""" + epochs = _epochs(n=5, k=6) + feats, _ = channel_covariance(epochs) + recon = unflatten_spd(feats, 6) + for i in range(len(epochs)): + assert np.allclose(recon[i], np.cov(epochs[i])) + + +def test_covariance_is_mutually_exclusive_with_bandpower_and_erp(): + with pytest.raises(ValueError, match="mutually exclusive"): + EEGLoader(include_covariance=True) # bandpower/erp default to True + # Explicitly disabling the others is allowed. + EEGLoader(include_bandpower=False, include_erp=False, include_covariance=True) + + +def test_covariance_block_feeds_the_riemann_model_end_to_end(): + pytest.importorskip("imblearn") + from behavioral_decoding.evaluation.cv import out_of_fold_proba + from behavioral_decoding.evaluation.metrics import classification_report + from behavioral_decoding.models.modality_models import build_modality_model + + # Plant a covariance-structured class difference across subjects. + rng = np.random.default_rng(1) + n, k, t = 144, 6, 200 + subjects = np.repeat(np.arange(9), n // 9) + y = rng.integers(0, 2, size=n) + epochs = np.empty((n, k, t)) + for i in range(n): + W = rng.normal(size=(k, k)) * 0.3 + np.eye(k) + if y[i] == 1: + W[0, 1] = W[1, 0] = 0.9 + epochs[i] = W @ rng.normal(size=(k, t)) + + block = EEGLoader( + include_bandpower=False, include_erp=False, include_covariance=True + ).from_arrays( + epochs, sfreq=200.0, times=np.linspace(0, 1, t), + subject_ids=subjects, stimulus_ids=np.arange(n), + ) + model = build_modality_model("eeg", y=y, base_learner="riemann", n_bags=8, seed=0) + oof, _ = out_of_fold_proba(model, block.X, y, block.subject_ids, n_splits=3, seed=0) + assert classification_report(y, oof)["balanced_accuracy"] > 0.6 diff --git a/tests/test_riemann.py b/tests/test_riemann.py new file mode 100644 index 0000000..a6514da --- /dev/null +++ b/tests/test_riemann.py @@ -0,0 +1,164 @@ +"""Riemannian tangent-space mapping and the `riemann` estimator path. + +The transformer maps flattened EEG covariance features into the log-Euclidean +tangent space, where the framework's scaler/SMOTE/logistic become valid again. +These tests check the vectorisation is isometric and invertible, the log map +survives rank-deficient covariances, the reference mean is leakage-safe, and the +wired-up `riemann` model classifies covariance structure above chance. + +Note on scope: on clean synthetic data the raw covariance entries are usually +already linearly separable, so these tests deliberately do NOT assert that +riemann beats a plain logistic. The tangent map's advantage is on ill-conditioned +real EEG; claiming a synthetic win would be cherry-picking. See docs/estimators.md. +""" + +from __future__ import annotations + +import numpy as np +import pytest + +from behavioral_decoding.models.riemann import ( + RiemannianTangentSpace, + flatten_spd, + n_channels_from_flat, + unflatten_spd, +) + +pytest.importorskip("imblearn") + + +def _spd(n, k, seed=0): + rng = np.random.default_rng(seed) + A = rng.normal(size=(n, k, k)) + return A @ A.transpose(0, 2, 1) + k * np.eye(k) + + +def _covariance_dataset(n=180, k=6, seed=0): + """Two classes differing in which channel pair is correlated.""" + rng = np.random.default_rng(seed) + subjects = np.repeat(np.arange(9), n // 9) + y = rng.integers(0, 2, size=n) + mats = [] + for i in range(n): + W = rng.normal(size=(k, k)) * 0.3 + np.eye(k) + if y[i] == 1: + W[0, 1] = W[1, 0] = 0.9 + else: + W[2, 3] = W[3, 2] = 0.9 + sig = W @ rng.normal(size=(k, 200)) + mats.append(np.cov(sig)) + return flatten_spd(np.stack(mats)), y, subjects + + +# ---------------------------------------------------------------- vectorisation + + +def test_n_channels_inferred_from_triangular_length(): + assert n_channels_from_flat(21) == 6 # 6*7/2 + assert n_channels_from_flat(3) == 2 + assert n_channels_from_flat(1) == 1 + + +def test_non_triangular_length_is_rejected(): + with pytest.raises(ValueError, match="not a flattened covariance|not k"): + n_channels_from_flat(20) + + +def test_flatten_is_invertible(): + A = _spd(5, 4) + back = unflatten_spd(flatten_spd(A), 4) + assert np.allclose(A, back) + + +def test_flatten_preserves_frobenius_norm(): + """Off-diagonal sqrt(2) scaling makes the vectorisation an isometry.""" + A = _spd(7, 5) + flat = flatten_spd(A) + assert np.allclose(np.linalg.norm(flat, axis=1), np.linalg.norm(A, axis=(1, 2))) + + +# ---------------------------------------------------------------- transformer + + +def test_transform_shape_and_finiteness(): + X, _, _ = _covariance_dataset() + Z = RiemannianTangentSpace().fit_transform(X) + assert Z.shape == X.shape + assert np.isfinite(Z).all() + + +def test_log_map_survives_rank_deficient_covariance(): + """Covariances from short epochs are singular; shrinkage must rescue the log.""" + k = 5 + # A rank-1 (singular) covariance: outer product of one vector. + v = np.arange(1.0, k + 1) + singular = np.outer(v, v) + flat = flatten_spd(singular[None]) + Z = RiemannianTangentSpace(shrinkage=1e-3).fit_transform(flat) + assert np.isfinite(Z).all() + + +def test_reference_mean_is_training_only_leakage_safe(): + """transform() on new data must use the train mean, not refit to it.""" + X, _, _ = _covariance_dataset(seed=0) + Xtr, Xte = X[:120], X[120:] + ts = RiemannianTangentSpace().fit(Xtr) + mean_after_fit = ts.mean_log_.copy() + ts.transform(Xte) + # transform must not have changed the stored reference. + assert np.array_equal(mean_after_fit, ts.mean_log_) + + +def test_transformer_rejects_non_covariance_features(): + X = np.random.default_rng(0).normal(size=(10, 20)) # 20 is not triangular + with pytest.raises(ValueError, match="not.*covariance|not k"): + RiemannianTangentSpace().fit(X) + + +def test_centering_changes_the_representation(): + X, _, _ = _covariance_dataset() + centered = RiemannianTangentSpace(center=True).fit_transform(X) + uncentered = RiemannianTangentSpace(center=False).fit_transform(X) + assert not np.allclose(centered, uncentered) + + +# ------------------------------------------------------------- estimator wiring + + +def test_riemann_model_prepends_tangent_step(): + from behavioral_decoding.models.modality_models import build_modality_model + + y = np.array([0] * 60 + [1] * 30) + model = build_modality_model("eeg", y=y, base_learner="riemann", n_bags=5, seed=0) + step_names = list(dict(model.estimator.steps)) + assert step_names[0] == "tangent" + assert "scaler" in step_names and "estimator" in step_names + assert model.bd_spec_["base_learner"] == "riemann" + + +def test_riemann_forces_max_features_to_one(): + """A column subset of a flattened covariance is not a covariance.""" + from behavioral_decoding.models.modality_models import build_modality_model + + y = np.array([0] * 60 + [1] * 30) + model = build_modality_model( + "eeg", y=y, base_learner="riemann", n_bags=5, max_features=0.5, seed=0 + ) + assert model.max_features == 1.0 + + +def test_riemann_classifies_covariance_structure_above_chance(): + """Positive control: the wired path recovers the planted class, honestly. + + Above chance is the claim, not 'beats logistic': raw covariance entries are + already separable on clean synthetic data. + """ + from behavioral_decoding.evaluation.cv import out_of_fold_proba + from behavioral_decoding.evaluation.metrics import classification_report + from behavioral_decoding.models.modality_models import build_modality_model + + X, y, subjects = _covariance_dataset(seed=1) + model = build_modality_model("eeg", y=y, base_learner="riemann", n_bags=8, seed=0) + oof, _ = out_of_fold_proba(model, X, y, subjects, n_splits=3, seed=0) + balacc = classification_report(y, oof)["balanced_accuracy"] + assert balacc > 0.6, f"riemann did not recover covariance structure (balacc={balacc:.3f})"