diff --git a/docs/source/api.rst b/docs/source/api.rst index 90d2b6b8c7..4c8a5999a9 100644 --- a/docs/source/api.rst +++ b/docs/source/api.rst @@ -467,6 +467,16 @@ accuracy, across-subject accuracy, or other transfer learning settings. CrossSessionSplitter CrossSubjectSplitter +A cross-subject transfer protocol additionally states what the estimator is +allowed to see of the held-out target subject, and how the rest of that +subject is scored. + +.. autosummary:: + :toctree: generated/ + :template: class.rst + + CrossSubjectMode + --------- Utilities --------- diff --git a/docs/source/whats_new.rst b/docs/source/whats_new.rst index d8a04bf21c..f3160f032d 100644 --- a/docs/source/whats_new.rst +++ b/docs/source/whats_new.rst @@ -23,6 +23,7 @@ Version 1.6 (Source - GitHub) Enhancements ~~~~~~~~~~~~ +- Add cross-subject transfer learning to :class:`moabb.evaluations.CrossSubjectEvaluation` through an optional target-calibration slice. :class:`moabb.evaluations.splitters.CrossSubjectSplitter` gains ``calibration_size`` (the fraction of each held-out subject/session pair set aside for adaptation, in ``[0, 1]``) and ``calibration_labeled``; when ``calibration_size > 0`` each fold becomes ``(train, calibration, test)``. The calibration trials never enter the training fold - they are handed to the pipeline steps that opt in via ``set_fit_request`` (``X_target_unlabeled``, or ``X_target_labeled`` and ``y_target_labeled`` when ``calibration_labeled=True``, plus the per-trial ``subjects`` array), so train, calibration and test stay trial-disjoint. Named presets are available through the new :class:`moabb.evaluations.CrossSubjectMode` enum and the ``cs_mode`` argument, covering train-only, unlabeled target adaptation at 20 / 50 / 100 percent, and labeled target calibration at 20 / 50 percent; ``TRAIN_TRIALWISE`` additionally scores one target trial at a time, so a method cannot exploit statistics of the whole test block. Two examples are added under ``examples/how_to_benchmark/`` (:gh:`1093` by `Bruno Aristimunha`_ and `Anton Andreev`_) - Add :class:`moabb.datasets.Wang2026`, one motor-imagery dataset from the Wang et al. 2026 sensory-guided joint-learning study. The release contains 39 globally unique participants: 31 in the primary randomized experiment (joint learning, n=15; BCI2000 control, n=8; tactile control, n=8) plus an independently recruited EEGNet-control cohort (n=8). A ``group`` filter preserves stable global IDs across the four cohort archives. The longitudinal protocol used online 1D and 2D cursor control with 62 EEG channels at 1000 Hz and four classes over at least four regular sessions plus a baseline. Bounded 8-MiB HTTP range read-ahead extracts one subject without downloading a complete 8.7--25.3 GB archive. Raw construction is explicitly disabled until the authors supply the physical units or gain/offset needed to convert released values to volts (:gh:`1126` by `Paul-Adrien Graignic`_ and `Bruno Aristimunha`_) - Add :class:`moabb.datasets.Lenaig2026` - SSAEP-BCI data of 48 participants in response to a set of four auditory stimuli: a pure tone (used as a reference), cicada song and cat's purr, and brownian noise. EEG acquisition is performed using a 24-channel (international 10-20 system, passive electrodes, impedance maintained below 10 kΩ) at a sampling rate of 500 Hz. The stimuli are amplitude-modulated by a 40 Hz sinusoid and have a duration of 10 seconds. The experiment is conducted at two loudness levels (60 and 66 phons, diotic presentation), with 24 participants each. The measurement consists in one session of two 10-minute runs (separated by a 5-minute break), each including 50 trials (10 repetitions per condition). (:gh:`1121` by `Henrique Lefundes`_) - Add :class:`moabb.datasets.Schrag2026Pediatric` — open-access pediatric SSVEP-BCI dataset (47 children aged 5-18, g.tec g.GAMMAsys + g.USBamp at 256 Hz, 16 channels) covering both an online 4-target SSVEP game (6.25 / 10 / 11.11 / 14.28 Hz) and an opt-in 12-stimulus personalization recording (4 contrasts x 3 sizes at 10 Hz). XDF + Unity markers; trial labels are derived from the matching ``Movements/`` CSV (live fbCCA classifier output). Single 1.2 GB zip on Zenodo (``10.5281/zenodo.19440997``) extracted per-subject on first use; the SSVEP game is exposed as two runs (standard and personal stimulus) of a single session (by `Bruno Aristimunha`_ and `Emily Schrag`_) diff --git a/examples/advanced_examples/plot_hinss2021_classification.py b/examples/advanced_examples/plot_hinss2021_classification.py index 7d1ab94af9..f470cd6f9d 100644 --- a/examples/advanced_examples/plot_hinss2021_classification.py +++ b/examples/advanced_examples/plot_hinss2021_classification.py @@ -27,7 +27,7 @@ from pyriemann.estimation import Covariances from pyriemann.spatialfilters import Xdawn from pyriemann.tangentspace import TangentSpace -from sklearn.base import TransformerMixin +from sklearn.base import BaseEstimator, TransformerMixin from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA from sklearn.pipeline import make_pipeline @@ -52,7 +52,7 @@ # select electrodes based on the covariance information -class EpochSelectChannel(TransformerMixin): +class EpochSelectChannel(TransformerMixin, BaseEstimator): """Select channels based on covariance information.""" def __init__(self, n_chan, cov_est): diff --git a/examples/how_to_benchmark/plot_cross_subject_transfer_rpa.py b/examples/how_to_benchmark/plot_cross_subject_transfer_rpa.py new file mode 100644 index 0000000000..8dc55ee7fe --- /dev/null +++ b/examples/how_to_benchmark/plot_cross_subject_transfer_rpa.py @@ -0,0 +1,326 @@ +""" +================================================= +Cross-subject transfer with Riemannian alignment +================================================= + +A cross-subject benchmark asks whether a model trained on several people can +generalize to a person it has never seen. This is harder than a random +train/test split because EEG covariance matrices vary substantially between +people, even when they perform the same task. + +This tutorial introduces a target-aware alternative inspired by Riemannian +Procrustes Analysis (RPA) [1]_. We will: + +1. separate source, target-calibration, and target-test trials; +2. recenter each subject's covariance matrices on the SPD manifold; +3. route an unlabeled target slice through a scikit-learn pipeline; and +4. compare the aligned pipeline with a source-only baseline. + +The example focuses on the *recentering* step of RPA. Full RPA also includes +scaling and rotation. Keeping one operation here makes both the geometry and +MOABB's transfer-learning interface visible. +""" + +# Authors: Anton Andonov +# Bruno Aristimunha +# +# License: BSD (3-clause) + +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +from matplotlib.patches import Rectangle +from pyriemann.estimation import Covariances +from pyriemann.preprocessing import Whitening +from pyriemann.tangentspace import TangentSpace +from sklearn import config_context +from sklearn.base import BaseEstimator, TransformerMixin +from sklearn.linear_model import LogisticRegression +from sklearn.pipeline import make_pipeline + +from moabb.datasets.fake import FakeDataset +from moabb.evaluations import CrossSubjectEvaluation +from moabb.evaluations.protocols import CrossSubjectMode +from moabb.paradigms import LeftRightImagery + + +############################################################################### +# The cross-subject transfer protocol +# ----------------------------------- +# +# Leave-one-subject-out evaluation repeats the following experiment: choose one +# person as the target and train on all remaining source subjects. A standard +# ``TRAIN`` fold exposes no target trial during fitting. +# +# A target-aware mode divides the held-out subject into two parts. For +# ``TRAIN_AND_TARGET_UNLABELED_20P``, the first 20% of every target session may +# be used without labels to estimate an alignment reference. The remaining 80% +# are untouched until scoring. Thus the calibration trials are neither source +# training samples nor test samples. + +fig, ax = plt.subplots(figsize=(10, 3.4)) +colors = {"source": "#4C78A8", "calibration": "#F2CF5B", "test": "#E45756"} + +for row, subject in enumerate(["Source 1", "Source 2", "Source 3", "Target"]): + if subject == "Target": + pieces = [ + (0.0, 0.2, "calibration\nX only", "calibration"), + (0.2, 0.8, "scored target trials", "test"), + ] + else: + pieces = [(0.0, 1.0, "source training trials", "source")] + + for left, width, label, role in pieces: + ax.add_patch( + Rectangle( + (left, row - 0.32), width, 0.64, facecolor=colors[role], edgecolor="white" + ) + ) + ax.text(left + width / 2, row, label, ha="center", va="center", fontsize=9) + +ax.set( + xlim=(0, 1), + ylim=(-0.7, 3.7), + yticks=range(4), + yticklabels=["Source 1", "Source 2", "Source 3", "Target"], + xlabel="Fraction of each subject's trials", + title="One TRAIN_AND_TARGET_UNLABELED_20P fold", +) +ax.invert_yaxis() +for spine in ("top", "right", "left"): + ax.spines[spine].set_visible(False) +ax.tick_params(axis="y", length=0) +fig.tight_layout() +plt.show() + + +############################################################################### +# Why align covariance matrices? +# ------------------------------ +# +# A trial with :math:`p` EEG channels is summarized by a +# :math:`p \\times p` symmetric positive-definite (SPD) covariance matrix +# :math:`C`. SPD matrices do not form a flat Euclidean space, so their average +# is represented by a Riemannian mean. +# +# For a domain :math:`d` (one source subject or the target subject), let +# :math:`G_d` be its Riemannian mean. Recentering applies +# +# .. math:: +# +# C' = G_d^{-1/2} C G_d^{-1/2}. +# +# The Riemannian mean of the transformed domain is the identity matrix. Each +# subject therefore keeps its trial-to-trial structure while a large part of +# its subject-specific covariance offset is removed. +# +# At training time, every source trial must use the reference of its own +# subject. At prediction time, every held-out trial must use the target +# reference estimated from the permitted *unlabeled* calibration slice. These +# two paths explain why the transformer defines both ``fit_transform`` and +# ``transform``. + + +class RiemannianAlignment(TransformerMixin, BaseEstimator): + """Recenter source and target covariance matrices by domain.""" + + def fit(self, X, y=None, *, subjects=None, X_target_unlabeled=None): + if subjects is None or X_target_unlabeled is None: + raise ValueError( + "RiemannianAlignment needs `subjects` and `X_target_unlabeled` metadata." + ) + + X = np.asarray(X) + subjects = np.asarray(subjects) + self.source_whiteners_ = { + subject: Whitening(metric="riemann").fit(X[subjects == subject]) + for subject in np.unique(subjects) + } + self.target_whitener_ = Whitening(metric="riemann").fit( + np.asarray(X_target_unlabeled) + ) + return self + + def fit_transform(self, X, y=None, *, subjects=None, X_target_unlabeled=None): + """Fit domain references and align the source training trials.""" + self.fit(X, y, subjects=subjects, X_target_unlabeled=X_target_unlabeled) + subjects = np.asarray(subjects) + X_aligned = np.empty_like(X) + for subject, whitener in self.source_whiteners_.items(): + mask = subjects == subject + X_aligned[mask] = whitener.transform(X[mask]) + return X_aligned + + def transform(self, X): + """Align unseen trials with the target reference.""" + return self.target_whitener_.transform(X) + + +############################################################################### +# Route the target slice through the pipeline +# ------------------------------------------- +# +# MOABB initially owns raw EEG epochs for both source and target trials. The +# alignment step, however, sits after ``Covariances`` and therefore expects SPD +# matrices. Scikit-learn's metadata routing handles this representation change: +# +# * ``set_fit_request`` declares the two fields consumed by the alignment step; +# * ``transform_input=["X_target_unlabeled"]`` sends the target slice through +# the already-fitted pipeline prefix before delivering it to that step. +# +# Consequently, ``RiemannianAlignment.fit`` receives source covariances as +# ``X`` and target covariances as ``X_target_unlabeled``. No MOABB-specific +# pipeline class is needed. pyRiemann 0.12 marks the stateless +# ``Covariances`` transformer as fitted, so it works directly with +# ``transform_input``. + +with config_context(enable_metadata_routing=True): + alignment = RiemannianAlignment().set_fit_request( + subjects=True, X_target_unlabeled=True + ) + +aligned_pipeline = make_pipeline( + Covariances("oas"), + alignment, + TangentSpace(metric="riemann"), + LogisticRegression(max_iter=500), + transform_input=["X_target_unlabeled"], +) + +source_only_pipeline = make_pipeline( + Covariances("oas"), TangentSpace(metric="riemann"), LogisticRegression(max_iter=500) +) + +############################################################################### +# The resulting data flow is compact: metadata is transformed only until the +# step that requests it, while the ordinary source ``X`` continues through the +# complete pipeline. + +fig, ax = plt.subplots(figsize=(10, 3.2)) +ax.axis("off") +nodes = [ + (0.03, 0.68, "Source epochs\n+ subject IDs", "#D8ECFF"), + (0.03, 0.20, "Unlabeled target\nepochs", "#FFF0CC"), + (0.35, 0.44, "Covariances", "#E8E8E8"), + (0.58, 0.44, "Riemannian\nalignment", "#E7DDFF"), + (0.82, 0.44, "Tangent space\n+ classifier", "#DFF4DF"), +] +for x, y, label, color in nodes: + ax.text( + x, + y, + label, + ha="left", + va="center", + bbox={"boxstyle": "round,pad=0.5", "facecolor": color, "edgecolor": "0.35"}, + ) +for start, end in [ + ((0.23, 0.68), (0.34, 0.52)), + ((0.23, 0.20), (0.34, 0.38)), + ((0.49, 0.44), (0.57, 0.44)), + ((0.73, 0.44), (0.81, 0.44)), +]: + ax.annotate("", xy=end, xytext=start, arrowprops={"arrowstyle": "->", "lw": 1.8}) +ax.set_title("Source data and routed target metadata share the fitted prefix") +fig.tight_layout() +plt.show() + + +############################################################################### +# Run the two benchmark contracts +# -------------------------------- +# +# A deterministic :class:`~moabb.datasets.fake.FakeDataset` keeps the tutorial +# fast and download-free. Its scores have no scientific meaning; it is used to +# expose the complete evaluation path. +# +# The baseline uses ``TRAIN`` and therefore sees no target data during fitting. +# The aligned pipeline uses ``TRAIN_AND_TARGET_UNLABELED_20P``. Recording the +# mode next to every result is essential because the two scores answer different +# questions and use different numbers of scored target trials. + +dataset = FakeDataset(["left_hand", "right_hand"], n_subjects=4, n_sessions=2, seed=42) +paradigm = LeftRightImagery() + +baseline = CrossSubjectEvaluation( + paradigm=paradigm, + datasets=[dataset], + cs_mode=CrossSubjectMode.TRAIN, + overwrite=True, + suffix="rpa_source_only", +).process({"Source only": source_only_pipeline}) +baseline["protocol"] = "Source only" + +aligned = CrossSubjectEvaluation( + paradigm=paradigm, + datasets=[dataset], + cs_mode=CrossSubjectMode.TRAIN_AND_TARGET_UNLABELED_20P, + overwrite=True, + suffix="rpa_unlabeled_20p", +).process({"Riemannian alignment": aligned_pipeline}) +aligned["protocol"] = "20% unlabeled target" + +results = pd.concat([baseline, aligned], ignore_index=True) +print( + results.groupby(["protocol", "pipeline"])["score"] + .agg(["mean", "std", "count"]) + .round(3) +) + + +############################################################################### +# Inspect paired target results +# ----------------------------- +# +# Each thin line below connects the same held-out subject and session. The +# diamonds show pipeline means. With a real dataset, this pairing is useful +# because subject difficulty often dominates the score variation. +# +# Do not interpret the random ranking produced by ``FakeDataset``. A proper +# study would repeat the benchmark over real datasets and compare methods under +# the same target-access protocol. In particular, a source-only score and a +# 20%-calibration score should never be presented as if they had identical +# information budgets. + +paired = results.pivot(index=["subject", "session"], columns="protocol", values="score") +order = ["Source only", "20% unlabeled target"] +colors = ["#4C78A8", "#E45756"] + +fig, ax = plt.subplots(figsize=(7.5, 5)) +for row in paired[order].dropna().to_numpy(): + ax.plot([0, 1], row, color="0.75", linewidth=1, zorder=1) + +for position, (protocol, color) in enumerate(zip(order, colors, strict=True)): + values = paired[protocol].dropna().to_numpy() + offsets = np.linspace(-0.06, 0.06, len(values)) + ax.scatter(position + offsets, values, color=color, alpha=0.7, zorder=2) + ax.scatter( + position, + values.mean(), + marker="D", + s=95, + color=color, + edgecolor="black", + zorder=3, + ) + +ax.axhline(0.5, color="black", linestyle="--", linewidth=1, label="Chance") +ax.set( + xticks=[0, 1], + xticklabels=order, + ylabel="ROC AUC", + title="Cross-subject scores under two target-access contracts", +) +ax.grid(axis="y", alpha=0.25) +ax.legend() +fig.tight_layout() +plt.show() + + +############################################################################### +# References +# ---------- +# .. [1] Rodrigues, P. L. C., Jutten, C., & Congedo, M. (2019). +# Riemannian Procrustes Analysis: Transfer Learning for Brain-Computer +# Interfaces. *IEEE Transactions on Biomedical Engineering*, 66(8), +# 2390-2401. https://doi.org/10.1109/TBME.2018.2889705 diff --git a/examples/how_to_benchmark/plot_cross_subject_transfer_sp_mdm.py b/examples/how_to_benchmark/plot_cross_subject_transfer_sp_mdm.py new file mode 100644 index 0000000000..87ad348755 --- /dev/null +++ b/examples/how_to_benchmark/plot_cross_subject_transfer_sp_mdm.py @@ -0,0 +1,314 @@ +""" +========================================================= +Trialwise cross-subject subject-prototype classification +========================================================= + +Minimum Distance to Mean (MDM) [1]_ represents each class by one Riemannian +mean covariance matrix. In a cross-subject setting, pooling every source +subject into one mean can hide useful subject structure. + +This tutorial builds a subject-prototype MDM (SP-MDM) classifier. It learns one +class prototype *per source subject*, routes subject identifiers to the +classifier during fitting, and predicts a held-out target trial from its mean +distance to those prototypes. + +We also use MOABB's strict ``TRAIN_TRIALWISE`` protocol. The fitted model sees +one target trial per prediction call, so it cannot estimate a normalization, +alignment, or other statistic from the complete target test block. +""" + +# Authors: Anton Andonov +# Bruno Aristimunha +# +# License: BSD (3-clause) + +import matplotlib.pyplot as plt +import numpy as np +from pyriemann.classification import MDM +from pyriemann.estimation import Covariances +from pyriemann.geometry.distance import distance_riemann +from pyriemann.geometry.mean import mean_riemann +from sklearn import config_context +from sklearn.base import BaseEstimator, ClassifierMixin +from sklearn.pipeline import make_pipeline + +from moabb.datasets.fake import FakeDataset +from moabb.evaluations import CrossSubjectEvaluation +from moabb.evaluations.protocols import CrossSubjectMode +from moabb.paradigms import LeftRightImagery + + +############################################################################### +# From one class mean to subject prototypes +# ----------------------------------------- +# +# Standard MDM estimates one prototype :math:`G_k` for every class :math:`k` +# after pooling all training trials. SP-MDM instead estimates +# :math:`G_{s,k}` for every source subject :math:`s` and class :math:`k`. +# +# For a new covariance matrix :math:`C`, its distance to class :math:`k` is +# +# .. math:: +# +# D_k(C) = \\frac{1}{|S_k|} +# \\sum_{s \\in S_k} d_R(C, G_{s,k}), +# +# where :math:`d_R` is the affine-invariant Riemannian distance and :math:`S_k` +# contains the source subjects for which class :math:`k` is available. The +# predicted class minimizes :math:`D_k(C)`. +# +# The sketch below shows why subject IDs matter. Blue and orange denote the two +# classes; marker shapes denote source subjects. Pooling would replace each +# colored set with one mean, whereas SP-MDM retains all six reference points. + +prototype_xy = np.array( + [[-1.8, 0.8], [-1.2, -0.7], [-0.5, 1.2], [1.5, 0.9], [0.8, -0.8], [1.9, -0.4]] +) +target_xy = np.array([0.35, 0.15]) +prototype_class = np.array([0, 0, 0, 1, 1, 1]) +markers = ["o", "s", "^"] +colors = ["#4C78A8", "#F58518"] + +fig, ax = plt.subplots(figsize=(7.5, 4.8)) +for index, point in enumerate(prototype_xy): + klass = prototype_class[index] + subject = index % 3 + ax.scatter( + *point, + s=110, + marker=markers[subject], + color=colors[klass], + edgecolor="black", + zorder=3, + ) + ax.plot( + [target_xy[0], point[0]], + [target_xy[1], point[1]], + color=colors[klass], + alpha=0.28, + linewidth=1.5, + ) + +ax.scatter( + *target_xy, + marker="*", + s=260, + color="#54A24B", + edgecolor="black", + label="Target trial", + zorder=4, +) +for subject, marker in enumerate(markers, start=1): + ax.scatter([], [], marker=marker, color="0.65", label=f"Source subject {subject}") +ax.set( + xlabel="Schematic manifold coordinate 1", + ylabel="Schematic manifold coordinate 2", + title="A target trial is compared with every subject/class prototype", +) +ax.legend(ncols=2, frameon=False) +ax.grid(alpha=0.2) +fig.tight_layout() +plt.show() + + +############################################################################### +# Implement SP-MDM as a small scikit-learn classifier +# --------------------------------------------------- +# +# The covariance matrices and task labels are ordinary ``X`` and ``y``. +# ``subjects`` is metadata: it identifies the domain of each training row but +# is not a prediction feature. The estimator therefore exposes it only as an +# optional keyword in ``fit`` and later requests it through scikit-learn's +# metadata-routing API. +# +# ``decision_function`` returns a signed score for binary ROC AUC. A positive +# value means that the covariance is closer to the second class than to the +# first. + + +class SubjectPrototypeMDM(ClassifierMixin, BaseEstimator): + """MDM using one class prototype per source subject.""" + + def fit(self, X, y, *, subjects=None): + if subjects is None: + raise ValueError("SubjectPrototypeMDM needs `subjects` metadata.") + + X = np.asarray(X) + y = np.asarray(y) + subjects = np.asarray(subjects) + self.classes_ = np.unique(y) + self.prototypes_ = { + klass: [ + mean_riemann(X[(subjects == subject) & (y == klass)]) + for subject in np.unique(subjects) + if np.any((subjects == subject) & (y == klass)) + ] + for klass in self.classes_ + } + return self + + def _distances(self, X): + return np.asarray( + [ + [ + np.mean([distance_riemann(cov, ref) for ref in self.prototypes_[k]]) + for k in self.classes_ + ] + for cov in X + ] + ) + + def predict(self, X): + return self.classes_[np.argmin(self._distances(X), axis=1)] + + def decision_function(self, X): + distances = self._distances(X) + if len(self.classes_) == 2: + return distances[:, 0] - distances[:, 1] + return -distances + + +############################################################################### +# Build a subject-aware pipeline +# ------------------------------ +# +# Only ``SubjectPrototypeMDM`` requests ``subjects``. ``Covariances`` remains a +# standard pyRiemann transformer, and a baseline MDM pipeline requests no +# metadata at all. + +with config_context(enable_metadata_routing=True): + sp_mdm = SubjectPrototypeMDM().set_fit_request(subjects=True) + +pipelines = { + "Pooled MDM": make_pipeline(Covariances("oas"), MDM(metric="riemann")), + "Subject-prototype MDM": make_pipeline(Covariances("oas"), sp_mdm), +} + + +############################################################################### +# What does trialwise scoring guarantee? +# -------------------------------------- +# +# ``TRAIN`` and ``TRAIN_TRIALWISE`` expose exactly the same source information +# during fitting: neither provides target calibration trials. Their difference +# is the batch presented during scoring. +# +# MOABB implements the strict mode with scikit-learn's +# :class:`~sklearn.frozen.FrozenEstimator`, +# :func:`~sklearn.model_selection.cross_val_predict`, and +# :class:`~sklearn.model_selection.LeaveOneOut`. ``FrozenEstimator`` makes every +# target-side ``fit`` a no-op, while ``LeaveOneOut`` makes every prediction fold +# contain one trial. The individual responses are collected before the +# session-level ROC AUC is computed. +# +# This restriction matters for a method that might otherwise compute target +# batch statistics inside ``predict``. It is redundant for a strictly inductive +# estimator such as the SP-MDM implementation above, but it makes the benchmark +# contract explicit and enforceable. + +fig, ax = plt.subplots(figsize=(9, 3.4)) +trial_x = np.arange(6) +for row, (label, color) in enumerate( + [("Blockwise: one call", "#D8ECFF"), ("Trialwise: six calls", "#FFF0CC")] +): + ax.scatter( + trial_x, np.full(6, 1 - row), marker="s", s=720, color=color, edgecolor="0.35" + ) + for trial in trial_x: + ax.text(trial, 1 - row, f"T{trial + 1}", ha="center", va="center") + ax.text(6.0, 1 - row, label, va="center", fontsize=10) + +ax.plot([-0.45, 5.45], [1.42, 1.42], color="#4C78A8", linewidth=2) +ax.plot([-0.45, -0.45], [1.28, 1.42], color="#4C78A8", linewidth=2) +ax.plot([5.45, 5.45], [1.28, 1.42], color="#4C78A8", linewidth=2) +ax.set(xlim=(-0.8, 7.6), ylim=(-0.55, 1.75), yticks=[], xticks=[]) +ax.set_title("Information available to each prediction call") +for spine in ax.spines.values(): + spine.set_visible(False) +fig.tight_layout() +plt.show() + + +############################################################################### +# Run a fair trialwise comparison +# ------------------------------- +# +# Both pipelines are evaluated under the same ``TRAIN_TRIALWISE`` protocol, so +# any difference comes from the classifier rather than a different target-data +# budget. The deterministic fake dataset keeps the example fast and +# download-free; its labels contain no meaningful physiological effect. + +dataset = FakeDataset(["left_hand", "right_hand"], n_subjects=4, n_sessions=2, seed=42) +paradigm = LeftRightImagery() +evaluation = CrossSubjectEvaluation( + paradigm=paradigm, + datasets=[dataset], + cs_mode=CrossSubjectMode.TRAIN_TRIALWISE, + overwrite=True, + suffix="subject_prototype_mdm", +) +results = evaluation.process(pipelines) + +print( + results.groupby("pipeline")["score"] + .agg(["mean", "std", "count"]) + .sort_values("mean", ascending=False) + .round(3) +) + + +############################################################################### +# Inspect paired held-out subjects +# -------------------------------- +# +# Every gray line connects the two classifiers on the same target subject and +# session. Diamonds show means. Because this is fake data, both methods should +# fluctuate around chance and the ranking must not be interpreted as evidence. +# +# On real datasets, use the same paired structure across many subjects and +# datasets, then apply MOABB's statistical-analysis utilities. Trialwise +# prediction is intentionally more expensive than blockwise prediction because +# the frozen estimator is invoked separately for every target trial. + +paired = results.pivot(index=["subject", "session"], columns="pipeline", values="score") +order = ["Pooled MDM", "Subject-prototype MDM"] +colors = ["#4C78A8", "#F58518"] + +fig, ax = plt.subplots(figsize=(7.5, 5)) +for row in paired[order].dropna().to_numpy(): + ax.plot([0, 1], row, color="0.75", linewidth=1, zorder=1) + +for position, (pipeline, color) in enumerate(zip(order, colors, strict=True)): + values = paired[pipeline].dropna().to_numpy() + offsets = np.linspace(-0.06, 0.06, len(values)) + ax.scatter(position + offsets, values, color=color, alpha=0.7, zorder=2) + ax.scatter( + position, + values.mean(), + marker="D", + s=95, + color=color, + edgecolor="black", + zorder=3, + ) + +ax.axhline(0.5, color="black", linestyle="--", linewidth=1, label="Chance") +ax.set( + xticks=[0, 1], + xticklabels=order, + ylabel="ROC AUC", + title="Paired target-subject/session scores", +) +ax.grid(axis="y", alpha=0.25) +ax.legend() +fig.tight_layout() +plt.show() + + +############################################################################### +# References +# ---------- +# .. [1] Barachant, A., Bonnet, S., Congedo, M., & Jutten, C. (2012). +# Multiclass Brain-Computer Interface Classification by Riemannian +# Geometry. *IEEE Transactions on Biomedical Engineering*, 59(4), +# 920-928. https://doi.org/10.1109/TBME.2011.2172210 diff --git a/moabb/evaluations/__init__.py b/moabb/evaluations/__init__.py index ee17fa525b..0bb2b8c1e1 100644 --- a/moabb/evaluations/__init__.py +++ b/moabb/evaluations/__init__.py @@ -10,6 +10,7 @@ WithinSessionEvaluation, WithinSubjectEvaluation, ) +from .protocols import CrossSubjectMode from .splitters import ( CrossDatasetSplitter, CrossSessionSplitter, diff --git a/moabb/evaluations/base.py b/moabb/evaluations/base.py index b5747ea43b..dffb908315 100644 --- a/moabb/evaluations/base.py +++ b/moabb/evaluations/base.py @@ -12,9 +12,13 @@ import numpy as np import pandas as pd from joblib import Parallel, delayed +from sklearn import config_context from sklearn.base import BaseEstimator, clone -from sklearn.model_selection import StratifiedKFold +from sklearn.frozen import FrozenEstimator +from sklearn.metrics import accuracy_score, roc_auc_score +from sklearn.model_selection import LeaveOneOut, StratifiedKFold, cross_val_predict from sklearn.preprocessing import LabelEncoder +from sklearn.utils.metadata_routing import get_routing_for_object from moabb.analysis import Results from moabb.datasets.base import ( # noqa: F401 - CacheConfig used in type hints @@ -44,6 +48,39 @@ log = logging.getLogger(__name__) + +def _score_trialwise(estimator, X, y, scoring): + """Score a fitted estimator from isolated, one-trial predictions. + + ``FrozenEstimator`` prevents refitting on target data and ``LeaveOneOut`` + makes every prediction fold contain exactly one target trial. + """ + if scoring not in ("accuracy", "roc_auc"): + raise TypeError( + "Trialwise scoring supports the built-in 'accuracy' and 'roc_auc' " + "metrics. Set one of them on the paradigm or use a blockwise " + "CrossSubjectMode for a custom scorer." + ) + + if scoring == "accuracy": + method = "predict" + elif hasattr(estimator, "decision_function"): + method = "decision_function" + else: + method = "predict_proba" + + response = cross_val_predict( + FrozenEstimator(estimator), X, y, cv=LeaveOneOut(), method=method + ) + if scoring == "accuracy": + score = accuracy_score(y, response) + else: + if response.ndim == 2: + response = response[:, 1] + score = roc_auc_score(y, response) + return {"score": score} + + # Making the optuna soft dependency @@ -100,6 +137,7 @@ def _evaluate_fold( session, cv_ind, split_metadata=None, + calib_idx=None, ): """Evaluate a single CV fold. Pure function, no shared mutable state. @@ -140,12 +178,12 @@ def _evaluate_fold( score_per_session = config["score_per_session"] mne_labels = config["mne_labels"] codecarbon_config = config["codecarbon_config"] + trialwise = config.get("trialwise", False) - # Label encode per fold (matching old per-session/per-subject scoping) + # Fit the encoder on source labels only. Target labels never influence the + # training representation. if not mne_labels: - le = LabelEncoder() - combined_y = np.concatenate([y[train_idx], y[test_idx]]) - le.fit(combined_y) + le = LabelEncoder().fit(y[train_idx]) y_train = le.transform(y[train_idx]) y_test = le.transform(y[test_idx]) else: @@ -174,6 +212,32 @@ def _evaluate_fold( tracker = emissions_obj.create_tracker() tracker.start() + # Optional transfer-learning calibration slice (raw). The protocol decides + # whether it is offered as unlabeled or labeled target metadata. + X_calib = y_calib = None + if calib_idx is not None and len(calib_idx): + X_calib = X[calib_idx] + calibration_labeled = split_metadata is not None and split_metadata.get( + "calibration_labeled", False + ) + if calibration_labeled: + y_calib = y[calib_idx] if mne_labels else le.transform(y[calib_idx]) + + fit_params = {} + if split_metadata is not None and "calibration_size" in split_metadata: + fit_params["subjects"] = metadata["subject"].to_numpy()[train_idx] + if X_calib is not None: + if y_calib is None: + fit_params["X_target_unlabeled"] = X_calib + else: + fit_params["X_target_labeled"] = X_calib + fit_params["y_target_labeled"] = y_calib + with config_context(enable_metadata_routing=True): + requested = get_routing_for_object(cvclf).consumes("fit", set(fit_params)) + fit_params = { + name: value for name, value in fit_params.items() if name in requested + } + # Fit model task_name = None emissions = math.nan @@ -181,7 +245,8 @@ def _evaluate_fold( task_name = str(uuid4()) tracker.start_task(task_name) t_start = perf_counter() - cvclf.fit(X[train_idx], y_train) + with config_context(enable_metadata_routing=True): + cvclf.fit(X[train_idx], y_train, **fit_params) duration = perf_counter() - t_start if tracker is not None: emissions_data = tracker.stop_task() @@ -206,7 +271,7 @@ def _evaluate_fold( ) _save_model_cv(model=cvclf, save_path=model_save_path, cv_index=str(cv_ind)) - scorer = _create_scorer(cvclf, scoring) + scorer = None if trialwise else _create_scorer(cvclf, scoring) # Build score groups: per-session or full test set if score_per_session: @@ -222,7 +287,10 @@ def _evaluate_fold( for group_idx, group_y, group_session in score_groups: is_error = False try: - score = scorer(cvclf, X[group_idx], group_y) + if trialwise: + score = _score_trialwise(cvclf, X[group_idx], group_y, scoring) + else: + score = scorer(cvclf, X[group_idx], group_y) except ValueError as err: if error_score == "raise": raise err @@ -676,6 +744,7 @@ def _build_eval_config(self, param_grid): self.emissions.codecarbon_config if _carbonfootprint else None ), "score_per_session": self._score_per_session, + "trialwise": getattr(self, "trialwise", False), "param_grid": None, # overridden per-task below if needed } @@ -683,13 +752,16 @@ def _build_eval_config(self, param_grid): def _preview_splits(splitter, y, metadata): """Materialize folds up front with optional splitter metadata.""" preview = [] - for cv_ind, (train_idx, test_idx) in enumerate(splitter.split(y, metadata)): + # ``*cal`` absorbs the optional calibration slice from a transfer + # splitter; a plain 2-tuple splitter gives cal == [] (no calibration). + for cv_ind, (train_idx, *cal, test_idx) in enumerate(splitter.split(y, metadata)): + calib_idx = cal[0] if cal else train_idx[:0] split_metadata = None if hasattr(splitter, "get_metadata"): split_metadata = splitter.get_metadata() if split_metadata is not None: split_metadata = dict(split_metadata) - preview.append((cv_ind, train_idx, test_idx, split_metadata)) + preview.append((cv_ind, train_idx, calib_idx, test_idx, split_metadata)) return preview def _build_task_list( @@ -700,7 +772,7 @@ def _build_task_list( config = self._build_eval_config(param_grid) fold_preview = self._preview_splits(splitter, y, metadata) - for cv_ind, train_idx, test_idx, split_meta in fold_preview: + for cv_ind, train_idx, calib_idx, test_idx, split_meta in fold_preview: test_meta = metadata.iloc[test_idx] subject = test_meta["subject"].iloc[0] @@ -728,6 +800,7 @@ def _build_task_list( "session": session, "cv_ind": cv_ind, "split_metadata": split_meta, + "calib_idx": calib_idx, } ) return tasks diff --git a/moabb/evaluations/evaluations.py b/moabb/evaluations/evaluations.py index 2d4221973b..cbc1559f93 100644 --- a/moabb/evaluations/evaluations.py +++ b/moabb/evaluations/evaluations.py @@ -8,6 +8,7 @@ from tqdm import tqdm from moabb.evaluations.base import BaseEvaluation +from moabb.evaluations.protocols import CrossSubjectMode, validate_transfer_protocol from moabb.evaluations.splitters import ( CrossSessionSplitter, CrossSubjectSplitter, @@ -459,6 +460,28 @@ class CrossSubjectEvaluation(BaseEvaluation): n_splits : int or None Number of splits for cross-validation. If None, the number of splits is equal to the number of subjects. Defaults to ``None``. + cv_class : type or None + Cross-validation strategy used to hold out subjects (e.g. + ``LeaveOneGroupOut``, ``GroupShuffleSplit``, ``GroupKFold``). Defaults to + ``None`` (``LeaveOneGroupOut``, or ``GroupKFold`` when ``n_splits`` is set). + cv_kwargs : dict + Keyword arguments for ``cv_class``. ``calibration_size`` (float in + ``[0, 1]``, default ``0.0``) enables transfer learning: when ``> 0`` each + fold becomes ``(train, calibration, test)``. The fraction is taken + within every held-out subject/session pair so each remains scorable, + and the calibration slice is routed (raw) to pipeline steps via + ``set_fit_request``. With ``calibration_labeled=False``, only + ``X_target_unlabeled`` may be routed. With ``calibration_labeled=True``, + ``X_target_labeled`` and ``y_target_labeled`` may be routed. + Labeled calibration is only allowed with ``calibration_size <= 0.5``. + cs_mode : CrossSubjectMode or str, default=CrossSubjectMode.TRAIN + Named cross-subject protocol preset. By default, this is the standard + train-only cross-subject evaluation with no target calibration. The + ``TRAIN_TRIALWISE`` mode additionally enforces one-trial-at-a-time + prediction during scoring and supports the built-in ``"accuracy"`` and + ``"roc_auc"`` metrics. Cannot be combined with manual + ``calibration_size`` or ``calibration_labeled`` in ``cv_kwargs``, except + for the default ``TRAIN`` mode. Notes ----- @@ -470,8 +493,48 @@ class CrossSubjectEvaluation(BaseEvaluation): _score_per_session = True _needs_all_subjects = True + def __init__(self, *args, cs_mode=CrossSubjectMode.TRAIN, **kwargs): + cv_kwargs = dict(kwargs.get("cv_kwargs") or {}) + + if cs_mode is None: + cs_mode = CrossSubjectMode.TRAIN + + cs_mode = CrossSubjectMode(cs_mode) + self.cs_mode = cs_mode + + # Manual cv_kwargs still work when the default train-only blockwise + # mode is used. + has_manual_calibration = ( + "calibration_size" in cv_kwargs or "calibration_labeled" in cv_kwargs + ) + + if has_manual_calibration and cs_mode != CrossSubjectMode.TRAIN: + raise ValueError( + "Pass either cs_mode or calibration_size/calibration_labeled, not both." + ) + + if not has_manual_calibration: + cv_kwargs["calibration_size"] = cs_mode.calibration_size + cv_kwargs["calibration_labeled"] = cs_mode.calibration_labeled + + self.trialwise = cs_mode.trialwise + + validate_transfer_protocol( + cv_kwargs.get("calibration_size", 0.0), + cv_kwargs.get("calibration_labeled", False), + ) + + kwargs["cv_kwargs"] = cv_kwargs + super().__init__(*args, **kwargs) + def _create_splitter(self): - """Create the CrossSubjectSplitter for parallel evaluation.""" + """Create the CrossSubjectSplitter for parallel evaluation. + + ``calibration_size`` and ``calibration_labeled`` passed via + ``cv_kwargs`` turn each fold into a transfer split: + + ``(train, calibration, test)``. + """ if self.n_splits is None: default_class = LeaveOneGroupOut default_kwargs = {} @@ -479,6 +542,7 @@ def _create_splitter(self): default_class = GroupKFold default_kwargs = {"n_splits": self.n_splits} + default_kwargs.update(self.cv_kwargs) cv_class, cv_kwargs = self._resolve_cv(default_class, default_kwargs) if self.groups is not None: cv_kwargs = {**cv_kwargs, "groups": self.groups} @@ -486,7 +550,6 @@ def _create_splitter(self): cv_class=cv_class, random_state=self.random_state, **cv_kwargs ) - # flake8: noqa: C901 def evaluate( self, dataset: "BaseDataset", @@ -498,101 +561,16 @@ def evaluate( if not self.is_valid(dataset): reason = self._get_incompatibility_reason(dataset) raise AssertionError( - f"Dataset '{dataset.code}' is not appropriate for {self.__class__.__name__}: {reason}" - ) - # this is a bit awkward, but we need to check if at least one pipe - # have to be run before loading the data. If at least one pipeline - # need to be run, we have to load all the data. - # we might need a better granularity, if we query the DB - run_pipes = {} - for subject in dataset.subject_list: - run_pipes.update( - self.results.not_yet_computed( - pipelines, dataset, subject, process_pipeline - ) + f"Dataset '{dataset.code}' is not appropriate for " + f"{self.__class__.__name__}: {reason}" ) - if len(run_pipes) == 0: - return - - X, y, metadata = self._load_data( - dataset, run_pipes, process_pipeline, postprocess_pipeline + yield from self._evaluate_parallel_dataset( + dataset=dataset, + pipelines=pipelines, + param_grid=param_grid, + process_pipeline=process_pipeline, + postprocess_pipeline=postprocess_pipeline, ) - le = LabelEncoder() - y = y if self.mne_labels else le.fit_transform(y) - - # extract metadata - groups = metadata.subject.values - sessions = metadata.session.values - n_subjects = len(dataset.subject_list) - nchan = self._get_nchan(X) - - # perform leave one subject out CV - self.cv = self._create_splitter() - if self.n_splits is not None and self.cv_class is None: - n_subjects = self.n_splits - - inner_cv = StratifiedKFold(3, shuffle=True, random_state=self.random_state) - - if _carbonfootprint: - # Initialise CodeCarbon per cross-validation - tracker = self.emissions.create_tracker() - tracker.start() - - # Progressbar at subject level - for cv_ind, (train, test) in enumerate( - tqdm( - self.cv.split(y, metadata), - total=n_subjects, - desc=f"{dataset.code}-CrossSubject", - ) - ): - subject = groups[test[0]] - # now we can check if this subject has results - run_pipes = self.results.not_yet_computed( - pipelines, dataset, subject, process_pipeline - ) - # iterate over pipelines - for name, clf in run_pipes.items(): - clf = self._grid_search( - param_grid=param_grid, name=name, grid_clf=clf, inner_cv=inner_cv - ) - cvclf = clone(clf) - - duration, emissions, task_name = self._fit_cv( - cvclf, X[train], y[train], tracker if _carbonfootprint else None - ) - self._maybe_save_model_cv( - cvclf, dataset, subject, "", name, cv_ind, eval_type="CrossSubject" - ) - - # Create scorer once per pipeline - scorer = _create_scorer(cvclf, self.paradigm.scoring) - - # Evaluate on each session - for session in np.unique(sessions[test]): - ix = sessions[test] == session - - res = self._build_scored_result( - dataset, - subject, - session, - name, - len(train), - nchan, - duration, - scorer, - cvclf, - X[test[ix]], - y[test[ix]], - ) - - if _carbonfootprint: - self._attach_emissions(res, emissions, task_name) - - yield res - - if _carbonfootprint: - tracker.stop() def is_valid(self, dataset: "BaseDataset") -> bool: return len(dataset.subject_list) > 1 @@ -600,11 +578,13 @@ def is_valid(self, dataset: "BaseDataset") -> bool: def _get_incompatibility_reason(self, dataset): """Get specific reason for dataset incompatibility.""" n_subjects = len(dataset.subject_list) + if n_subjects <= 1: return ( f"dataset has only {n_subjects} subject(s), " f"but {self.__class__.__name__} requires at least 2 subjects" ) + return "requirements not met" diff --git a/moabb/evaluations/protocols.py b/moabb/evaluations/protocols.py new file mode 100644 index 0000000000..7821ce9c85 --- /dev/null +++ b/moabb/evaluations/protocols.py @@ -0,0 +1,64 @@ +"""Named target-access protocols for cross-subject evaluation.""" + +from enum import Enum + + +class CrossSubjectMode(str, Enum): + """Target data made available in a cross-subject benchmark. + + A mode fixes the target calibration fraction, whether calibration labels + are routed, and whether prediction is blockwise or trialwise. + + Pass a member to the ``cs_mode`` parameter of + :class:`moabb.evaluations.CrossSubjectEvaluation`. + """ + + def __new__(cls, value, calibration_size, calibration_labeled, trialwise=False): + member = str.__new__(cls, value) + member._value_ = value + member.calibration_size = calibration_size + member.calibration_labeled = calibration_labeled + member.trialwise = trialwise + return member + + # Train only on source subjects and predict the target block normally. + TRAIN = ("train", 0.0, False) + + # Train only on source subjects and predict one target trial at a time. + TRAIN_TRIALWISE = ("train_trialwise", 0.0, False, True) + + # Use an unlabeled target slice for adaptation. + TRAIN_AND_TARGET_UNLABELED_20P = ("train_and_target_unlabeled_20p", 0.2, False) + TRAIN_AND_TARGET_UNLABELED_50P = ("train_and_target_unlabeled_50p", 0.5, False) + + # Transductive: adapt on the same unlabeled target block that is scored. + TRAIN_AND_TARGET_UNLABELED_FULL = ("train_and_target_unlabeled_full", 1.0, False) + + # Use a labeled target slice for calibration. + TRAIN_AND_TARGET_LABELED_20P = ("train_and_target_labeled_20p", 0.2, True) + TRAIN_AND_TARGET_LABELED_50P = ("train_and_target_labeled_50p", 0.5, True) + + +def validate_transfer_protocol(calibration_size, calibration_labeled): + if isinstance(calibration_size, bool) or not isinstance( + calibration_size, (int, float) + ): + raise TypeError( + f"calibration_size must be a number. Got {type(calibration_size).__name__}." + ) + + calibration_size = float(calibration_size) + + if not 0.0 <= calibration_size <= 1.0: + raise ValueError(f"calibration_size must be in [0, 1]. Got {calibration_size!r}.") + + if not isinstance(calibration_labeled, bool): + raise TypeError( + "calibration_labeled must be a bool. " + f"Got {type(calibration_labeled).__name__}." + ) + + if calibration_labeled and calibration_size > 0.5: + raise ValueError( + "calibration_labeled=True is only allowed with calibration_size <= 0.5." + ) diff --git a/moabb/evaluations/splitters.py b/moabb/evaluations/splitters.py index 2686a952c5..a94f24a7cc 100644 --- a/moabb/evaluations/splitters.py +++ b/moabb/evaluations/splitters.py @@ -14,6 +14,8 @@ from sklearn.model_selection._split import GroupsConsumerMixin from sklearn.utils import check_random_state +from moabb.evaluations.protocols import validate_transfer_protocol + log = logging.getLogger(__name__) @@ -626,6 +628,15 @@ class CrossSubjectSplitter(BaseCrossValidator): Controls the randomness of the cross-validation. Pass an int for reproducible output across multiple calls. Defaults to ``None``. + calibration_size : float, default=0.0 + Fraction of each held-out subject/session pair reserved for + calibration. Values in ``(0, 1)`` yield + ``(train, calibration, test)``; ``0`` keeps the ordinary + ``(train, test)`` split, and ``1`` uses the full target block both for + unlabeled adaptation and scoring. + calibration_labeled : bool, default=False + If True, route labels with the calibration slice. Labeled calibration + is limited to at most half of the target trials. cv_kwargs : dict Additional arguments to pass to the inner cross-validation strategy. A callable value is resolved against the metadata at ``split`` time @@ -636,6 +647,10 @@ class CrossSubjectSplitter(BaseCrossValidator): train : ndarray The training set indices for that split. + calibration : ndarray + The held-out calibration indices, yielded only when + ``calibration_size`` is greater than zero. + test : ndarray The testing set indices for that split. """ @@ -645,13 +660,20 @@ def __init__( cv_class: type[BaseCrossValidator] = LeaveOneGroupOut, groups="subject", random_state: int = None, + calibration_size: float = 0.0, + calibration_labeled: bool = False, **cv_kwargs, ): + validate_transfer_protocol(calibration_size, calibration_labeled) + self.cv_class = cv_class # ``groups`` selects what defines a fold: a metadata column ("subject"), # a list of columns (["subject", "session"]), or a callable # ``metadata -> array``. Fed straight to the stock sklearn cv_class. self.groups = groups + self.calibration_size = calibration_size + self.calibration_labeled = calibration_labeled + self.random_state = random_state self.cv_kwargs = cv_kwargs self._cv_kwargs = dict(**cv_kwargs) @@ -705,8 +727,42 @@ def split(self, y, metadata): split_kwargs["groups"] = _resolve_groups(self.groups, metadata) for train_session_idx, test_session_idx in splitter.split(**split_kwargs): - self._last_split_metadata = _splitter_metadata(splitter) - yield all_index[train_session_idx], all_index[test_session_idx] + self._last_split_metadata = _splitter_metadata(splitter) or {} + self._last_split_metadata.update( + { + "calibration_size": self.calibration_size, + "calibration_labeled": self.calibration_labeled, + } + ) + train_idx = all_index[train_session_idx] + target_idx = all_index[test_session_idx] + + if self.calibration_size > 0: + if self.calibration_size == 1: + calib_idx = test_idx = target_idx + else: + target_metadata = metadata.loc[target_idx, ["subject", "session"]] + target_subjects = target_metadata["subject"].to_numpy() + target_sessions = target_metadata["session"].to_numpy() + calibration_mask = np.zeros(len(target_idx), dtype=bool) + for subject in np.unique(target_subjects): + for session in np.unique( + target_sessions[target_subjects == subject] + ): + positions = np.flatnonzero( + (target_subjects == subject) + & (target_sessions == session) + ) + n_calib = int( + np.floor(self.calibration_size * len(positions)) + ) + n_calib = min(max(n_calib, 1), len(positions) - 1) + calibration_mask[positions[:n_calib]] = True + calib_idx = target_idx[calibration_mask] + test_idx = target_idx[~calibration_mask] + yield train_idx, calib_idx, test_idx + else: + yield train_idx, target_idx def get_metadata(self): """Return metadata for the most recent split.""" diff --git a/moabb/tests/test_evaluations.py b/moabb/tests/test_evaluations.py index a74db9b54b..028c7f9e47 100644 --- a/moabb/tests/test_evaluations.py +++ b/moabb/tests/test_evaluations.py @@ -455,6 +455,38 @@ def test_within_n_splits_drives_n_folds(klass): os.remove(e.results.filepath) +def test_resolve_cv_honours_cv_kwargs_without_forcing_defaults(): + """``cv_kwargs`` always reaches the splitter; ``default_kwargs`` never + reaches a user-supplied ``cv_class`` that cannot accept it.""" + from sklearn.model_selection import LeaveOneGroupOut + + kw = { + "paradigm": FakeImageryParadigm(), + "datasets": [dataset], + "hdf5_path": "res_test", + } + evals = [] + try: + # Default cv_class: cv_kwargs must survive (this is how + # calibration_size is passed). + e = ev.CrossSubjectEvaluation(cv_kwargs={"calibration_size": 0.5}, **kw) + evals.append(e) + assert e._create_splitter().calibration_size == 0.5 + + # User-supplied cv_class: n_splits is a default of GroupKFold, not of + # LeaveOneGroupOut, which takes no arguments at all. + e = ev.CrossSubjectEvaluation(cv_class=LeaveOneGroupOut, n_splits=3, **kw) + evals.append(e) + splitter = e._create_splitter() + assert splitter._cv_kwargs == {} + _, y, metadata = FakeImageryParadigm().get_data(dataset) + assert splitter.get_n_splits(metadata) == metadata["subject"].nunique() + finally: + for e in evals: + if os.path.isfile(e.results.filepath): + os.remove(e.results.filepath) + + class Test_CrossSubj(TestWithinSess): def setup_method(self): self.eval = ev.CrossSubjectEvaluation( @@ -1068,3 +1100,266 @@ def test_non_numeric_score_fold_does_not_abort(self): assert len(agg) == 1 # Non-numeric fold is coerced to NaN and skipped by mean -> only 0.7 left. np.testing.assert_almost_equal(agg[0]["score"], 0.7) + + +# Transfer-learning calibration through the stock CrossSubjectEvaluation. +_TRANSFER_CAPTURE = [] + + +class _TransferRecorder(sklearn.base.TransformerMixin, sklearn.base.BaseEstimator): + """A target-aware step that records the transfer metadata it is routed.""" + + def fit(self, X, y=None, subjects=None, X_target_unlabeled=None): + _TRANSFER_CAPTURE.append( + { + "n_subjects": 0 if subjects is None else len(subjects), + "n_target": 0 if X_target_unlabeled is None else len(X_target_unlabeled), + } + ) + return self + + def transform(self, X): + return X + + +def test_cross_subject_calibration_routes_to_estimator(): + """calibration_size>0 routes subjects + the (raw) calibration slice to the + steps that request it, via the unchanged CrossSubjectEvaluation.""" + from sklearn import config_context + + _TRANSFER_CAPTURE.clear() + with config_context(enable_metadata_routing=True): + step = _TransferRecorder().set_fit_request(subjects=True, X_target_unlabeled=True) + pipe = make_pipeline(Covariances("oas"), step, CSP(8), LDA()) + ds = FakeDataset(["left_hand", "right_hand"], n_subjects=3, n_sessions=2, seed=9) + + evaluation = ev.CrossSubjectEvaluation( + paradigm=FakeImageryParadigm(), + datasets=[ds], + cv_kwargs={"calibration_size": 0.5}, + overwrite=True, + n_jobs=1, + suffix="calibroute", + ) + results = evaluation.process(pipelines=OrderedDict([("T", pipe)])) + + assert len(results) > 0 + assert _TRANSFER_CAPTURE, "transfer step was never fitted" + assert all(c["n_subjects"] > 0 for c in _TRANSFER_CAPTURE) + assert all(c["n_target"] > 0 for c in _TRANSFER_CAPTURE) + + +def test_cross_subject_calibration_leaves_plain_pipeline_unaffected(): + """A plain pipeline runs through calibration_size>0 unchanged (calib ignored).""" + pipe = make_pipeline(Covariances("oas"), CSP(8), LDA()) + ds = FakeDataset(["left_hand", "right_hand"], n_subjects=3, n_sessions=2, seed=9) + + evaluation = ev.CrossSubjectEvaluation( + paradigm=FakeImageryParadigm(), + datasets=[ds], + cv_kwargs={"calibration_size": 0.5}, + overwrite=True, + n_jobs=1, + suffix="calibplain", + ) + results = evaluation.process(pipelines=OrderedDict([("P", pipe)])) + assert len(results) > 0 + + +def test_cross_subject_calibration_with_custom_cv(): + """cv_class is exposed (like WithinSession) and composes with calibration.""" + from sklearn import config_context + from sklearn.model_selection import GroupShuffleSplit + + _TRANSFER_CAPTURE.clear() + with config_context(enable_metadata_routing=True): + step = _TransferRecorder().set_fit_request(subjects=True, X_target_unlabeled=True) + pipe = make_pipeline(Covariances("oas"), step, CSP(8), LDA()) + ds = FakeDataset(["left_hand", "right_hand"], n_subjects=5, n_sessions=2, seed=9) + + evaluation = ev.CrossSubjectEvaluation( + paradigm=FakeImageryParadigm(), + datasets=[ds], + cv_class=GroupShuffleSplit, + cv_kwargs={"calibration_size": 0.5, "n_splits": 2}, + random_state=0, + overwrite=True, + n_jobs=1, + suffix="calibcv", + ) + results = evaluation.process(pipelines=OrderedDict([("T", pipe)])) + + assert len(results) > 0 + assert _TRANSFER_CAPTURE, "transfer step was never fitted" + assert all(c["n_subjects"] > 0 and c["n_target"] > 0 for c in _TRANSFER_CAPTURE) + + +# --------------------------------------------------------------------------- +# CrossSubjectMode presets and the trialwise scoring guarantee. +# --------------------------------------------------------------------------- + + +class _PredictSpy(sklearn.base.ClassifierMixin, sklearn.base.BaseEstimator): + """Records the batch size of every predict/score call.""" + + def __init__(self): + self.calls = [] + + def fit(self, X, y): + self.classes_ = np.unique(y) + return self + + def predict(self, X): + self.calls.append(("predict", len(X))) + return np.zeros(len(X), dtype=int) + + def score(self, X, y): + self.calls.append(("score", len(X))) + return 1.0 + + +def test_trialwise_score_predicts_one_trial_at_a_time(): + """FrozenEstimator + LeaveOneOut isolates every target prediction.""" + from moabb.evaluations.base import _score_trialwise + + X = np.random.RandomState(0).randn(8, 3) + y = np.array([0, 1] * 4) + + spy = _PredictSpy().fit(X, y) + score = _score_trialwise(spy, X, y, "accuracy") + + assert score == {"score": 0.5} + assert spy.calls, "the frozen estimator was never called" + assert {n for _, n in spy.calls} == {1} + + +def test_trialwise_score_refuses_passthrough_scoring(): + """scoring=None must not silently call the estimator on the whole block. + + ``check_scoring(estimator, scoring=None)`` returns a passthrough scorer that + calls ``estimator.score(X, y)``. If that is delegated to the frozen + estimator, the whole target test block goes through in one call and the + trialwise guarantee is void. + """ + from moabb.evaluations.base import _score_trialwise + + X = np.random.RandomState(0).randn(8, 3) + y = np.array([0, 1] * 4) + + spy = _PredictSpy().fit(X, y) + with pytest.raises(TypeError, match="Trialwise scoring supports"): + _score_trialwise(spy, X, y, None) + assert spy.calls == [] + + +def test_trialwise_score_supports_probability_scorers(): + """The official CV recipe also supports metrics such as ROC AUC.""" + from moabb.evaluations.base import _score_trialwise + + X = np.random.RandomState(0).randn(8, 3) + y = np.array([0, 1] * 4) + + score = _score_trialwise(LDA().fit(X, y), X, y, "roc_auc") + assert 0.0 <= score["score"] <= 1.0 + + +@pytest.mark.parametrize( + "mode,calibration_size,calibration_labeled,trialwise", + [ + ("TRAIN", 0.0, False, False), + ("TRAIN_TRIALWISE", 0.0, False, True), + ("TRAIN_AND_TARGET_UNLABELED_20P", 0.2, False, False), + ("TRAIN_AND_TARGET_UNLABELED_50P", 0.5, False, False), + ("TRAIN_AND_TARGET_UNLABELED_FULL", 1.0, False, False), + ("TRAIN_AND_TARGET_LABELED_20P", 0.2, True, False), + ("TRAIN_AND_TARGET_LABELED_50P", 0.5, True, False), + ], +) +def test_cs_mode_resolves_to_splitter_kwargs( + mode, calibration_size, calibration_labeled, trialwise +): + """Every CrossSubjectMode maps to the documented calibration settings.""" + from moabb.evaluations import CrossSubjectMode + + evaluation = ev.CrossSubjectEvaluation( + paradigm=FakeImageryParadigm(), + datasets=[dataset], + hdf5_path="res_test", + cs_mode=getattr(CrossSubjectMode, mode), + ) + try: + splitter = evaluation._create_splitter() + assert splitter.calibration_size == calibration_size + assert splitter.calibration_labeled is calibration_labeled + assert evaluation.trialwise is trialwise + # A plain string is accepted and normalised to the enum member. + by_value = ev.CrossSubjectEvaluation( + paradigm=FakeImageryParadigm(), + datasets=[dataset], + hdf5_path="res_test", + cs_mode=getattr(CrossSubjectMode, mode).value, + ) + assert by_value.cs_mode is getattr(CrossSubjectMode, mode) + finally: + if os.path.isfile(evaluation.results.filepath): + os.remove(evaluation.results.filepath) + + +def test_cs_mode_rejects_manual_calibration_kwargs(): + """cs_mode and manual calibration kwargs are mutually exclusive.""" + from moabb.evaluations import CrossSubjectMode + + with pytest.raises(ValueError, match="not both"): + ev.CrossSubjectEvaluation( + paradigm=FakeImageryParadigm(), + datasets=[dataset], + hdf5_path="res_test", + cs_mode=CrossSubjectMode.TRAIN_AND_TARGET_UNLABELED_20P, + cv_kwargs={"calibration_size": 0.5}, + ) + + +def test_cs_mode_rejects_labeled_calibration_above_half(): + """calibration_labeled=True is capped at calibration_size <= 0.5.""" + with pytest.raises(ValueError, match="calibration_labeled"): + ev.CrossSubjectEvaluation( + paradigm=FakeImageryParadigm(), + datasets=[dataset], + hdf5_path="res_test", + cv_kwargs={"calibration_size": 1.0, "calibration_labeled": True}, + ) + + +def test_cs_mode_trialwise_runs_end_to_end(): + """TRAIN_TRIALWISE uses frozen leave-one-out predictions and still reports + per-session rows, matching the blockwise TRAIN baseline.""" + from moabb.evaluations import CrossSubjectMode + + pipe = make_pipeline(Covariances("oas"), CSP(8), LDA()) + ds = FakeDataset(["left_hand", "right_hand"], n_subjects=3, n_sessions=2, seed=9) + + scores = {} + for mode, suffix in [ + (CrossSubjectMode.TRAIN, "csmodetrain"), + (CrossSubjectMode.TRAIN_TRIALWISE, "csmodetrialwise"), + ]: + evaluation = ev.CrossSubjectEvaluation( + paradigm=FakeImageryParadigm(), + datasets=[ds], + cs_mode=mode, + overwrite=True, + n_jobs=1, + suffix=suffix, + ) + results = evaluation.process(pipelines=OrderedDict([("P", pipe)])) + scores[mode] = results + + assert len(scores[CrossSubjectMode.TRAIN_TRIALWISE]) == len( + scores[CrossSubjectMode.TRAIN] + ) + # Trialwise vs blockwise only changes the batching, not the predictions, + # for a pipeline that treats trials independently. + np.testing.assert_allclose( + scores[CrossSubjectMode.TRAIN_TRIALWISE]["score"].to_numpy(), + scores[CrossSubjectMode.TRAIN]["score"].to_numpy(), + ) diff --git a/moabb/tests/test_splits.py b/moabb/tests/test_splits.py index 879c3c7812..e411750b8e 100644 --- a/moabb/tests/test_splits.py +++ b/moabb/tests/test_splits.py @@ -855,3 +855,101 @@ def test_within_subject_groups_routes_through(data): same_subject_train = metadata.loc[train] same_subject_train = same_subject_train[same_subject_train["subject"] == subj] assert held_out_session not in set(same_subject_train["session"]) + + +# --------------------------------------------------------------------------- +# Cross-subject transfer learning: the target-calibration slice. +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize("calibration_size", [0.0, 0.3, 1.0]) +def test_cross_subject_calibration(calibration_size, data): + """Calibration folds are consumed uniformly with ``train, *cal, test``.""" + _, y, metadata = data + base = CrossSubjectSplitter() + cal_split = CrossSubjectSplitter(calibration_size=calibration_size) + + base_folds = list(base.split(y, metadata)) + cal_folds = list(cal_split.split(y, metadata)) + assert len(cal_folds) == len(base_folds) == cal_split.get_n_splits(metadata) + + for (b_train, b_test), fold in zip(base_folds, cal_folds): + train, *cal, test = fold # generic consumption: 2- or 3-tuple + calib = cal[0] if cal else test[:0] + + assert np.array_equal(train, b_train) + assert np.intersect1d(train, b_test).size == 0 + + if calibration_size == 0.0: + assert calib.size == 0 + assert np.array_equal(test, b_test) + elif calibration_size == 1.0: + assert np.array_equal(calib, b_test) + assert np.array_equal(test, b_test) + else: + assert calib.size >= 1 and test.size >= 1 + assert np.array_equal(np.union1d(calib, test), b_test) + + +def test_cross_subject_calibration_invalid_size(): + with pytest.raises(ValueError): + CrossSubjectSplitter(calibration_size=1.5) + + +@pytest.mark.parametrize("calibration_labeled", [False, True]) +def test_cross_subject_calibration_leakage_boundary(calibration_labeled, data): + """Pin the leakage contract of the transfer split. + + The held-out fold is a single target subject. The calibration slice is + carved out of that target subject and reaches the estimator only through + ``fit`` metadata routing -- never through ``train_idx`` -- and it is + removed from the scored test set. Train, calibration and test are + pairwise trial-disjoint. + """ + _, y, metadata = data + splitter = CrossSubjectSplitter( + calibration_size=0.2, calibration_labeled=calibration_labeled + ) + folds = list(splitter.split(y, metadata)) + assert len(folds) == metadata["subject"].nunique() + + for train_idx, calib_idx, test_idx in folds: + target = set(metadata.loc[test_idx, "subject"]) + # Exactly one held-out target subject per fold. + assert len(target) == 1 + # The calibration slice belongs to that same target subject... + assert set(metadata.loc[calib_idx, "subject"]) == target + assert calib_idx.size > 0 + # ...and no trial of the target subject is in the training fold. + assert target.isdisjoint(set(metadata.loc[train_idx, "subject"])) + # Pairwise trial-disjoint: nothing is scored that was also fitted. + assert np.intersect1d(train_idx, test_idx).size == 0 + assert np.intersect1d(train_idx, calib_idx).size == 0 + assert np.intersect1d(calib_idx, test_idx).size == 0 + # Calibration + test partition the target subject exactly. + target_idx = metadata.index[metadata["subject"].isin(target)].to_numpy() + assert np.array_equal(np.union1d(calib_idx, test_idx), target_idx) + + +def test_cross_subject_calibration_keeps_every_target_session(data): + """Calibration must not consume an entire target session.""" + _, y, metadata = data + cv_kwargs = { + "cv_class": GroupShuffleSplit, + "n_splits": 2, + "test_size": 2, + "random_state": 0, + } + baseline = list(CrossSubjectSplitter(**cv_kwargs).split(y, metadata)) + calibrated = list( + CrossSubjectSplitter(calibration_size=0.5, **cv_kwargs).split(y, metadata) + ) + + for (_, base_test), (_, _calib, test) in zip(baseline, calibrated): + assert set( + metadata.loc[test, ["subject", "session"]].itertuples(index=False, name=None) + ) == set( + metadata.loc[base_test, ["subject", "session"]].itertuples( + index=False, name=None + ) + )