Riemannian tangent-space EEG path

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.

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
 

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:

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,