From fe2a44b82dee68b4854d2f13e802564d9bfa0846 Mon Sep 17 00:00:00 2001 From: qbarthelemy Date: Fri, 16 Oct 2026 21:34:04 +0200 Subject: [PATCH 1/3] Add CrossSystemEvaluation --- docs/source/api.rst | 1 + docs/source/whats_new.rst | 1 + .../plot_cross_system_transfer_maskedmean.py | 577 ++++++++++++++++++ moabb/evaluations/__init__.py | 1 + moabb/evaluations/base.py | 16 +- moabb/evaluations/evaluations.py | 101 +++ moabb/evaluations/utils.py | 7 +- moabb/tests/test_cross_system_evaluation.py | 225 +++++++ moabb/tests/test_evaluations.py | 123 ++++ 9 files changed, 1046 insertions(+), 6 deletions(-) create mode 100644 examples/how_to_benchmark/plot_cross_system_transfer_maskedmean.py create mode 100644 moabb/tests/test_cross_system_evaluation.py diff --git a/docs/source/api.rst b/docs/source/api.rst index e29fa40be..f6cd274eb 100644 --- a/docs/source/api.rst +++ b/docs/source/api.rst @@ -535,6 +535,7 @@ accuracy, across-subject accuracy, or other transfer learning settings. WithinSessionEvaluation CrossSessionEvaluation CrossSubjectEvaluation + CrossSystemEvaluation .. autosummary:: :toctree: generated/ diff --git a/docs/source/whats_new.rst b/docs/source/whats_new.rst index e0bfa9e8a..31d1c817c 100644 --- a/docs/source/whats_new.rst +++ b/docs/source/whats_new.rst @@ -23,6 +23,7 @@ Version 1.8.0 (2026-10-06) Enhancements ~~~~~~~~~~~~ +- Add :class:`moabb.evaluations.CrossSystemEvaluation`, a cross-system evaluation where each dataset is a system (an EEG device with its own montage and sampling rate): for each target system, chosen with ``target_datasets``, pipelines are trained on all other systems. Trials are kept as recorded, with the channels and sampling rate of their system (an object array of trials when systems differ), handling heterogeneous systems being left to the pipelines; pipelines can request the ``systems`` fit metadata through scikit-learn's metadata routing, and target calibration presets are the ones of :class:`moabb.evaluations.CrossSubjectMode`. A new example stacks the systems on the union of their channels, and compares a channel-intersection baseline with a minimum distance to Riemannian masked mean classifier (by `Quentin Barthelemy`_). - Add ``n_jobs_data`` to all evaluations (default ``1``, unchanged behavior): it is forwarded as ``n_jobs`` to :meth:`~moabb.paradigms.base.BaseParadigm.get_data`, so subjects are loaded in parallel. :meth:`moabb.datasets.base.BaseDataset.get_data` already supported ``n_jobs`` but evaluations never passed it, so loading was always sequential; this matters when reads are latency-bound, e.g. raw EDF files on a network filesystem (:gh:`1250` by `Bruno Aristimunha`_). - Install a ``moabb`` console script (``[project.scripts] moabb = "moabb.run:main"``), so benchmarks run with ``moabb --pipelines=./pipelines/ --evaluations WithinSession`` instead of ``python -m moabb.run ...``; the command-line logic of :mod:`moabb.run` moved into a reusable ``main(argv=None)`` (:gh:`175` by `Bruno Aristimunha`_). - Test MOABB on free-threaded CPython 3.14 (``3.14t``) with an allowed-to-fail CI job that also checks importing MOABB keeps the GIL disabled, and require ``h5py>=3.16.0`` on Python 3.14+, the first h5py release with free-threaded (``cp314t``) wheels (:gh:`829` by `Bruno Aristimunha`_). diff --git a/examples/how_to_benchmark/plot_cross_system_transfer_maskedmean.py b/examples/how_to_benchmark/plot_cross_system_transfer_maskedmean.py new file mode 100644 index 000000000..6b486e902 --- /dev/null +++ b/examples/how_to_benchmark/plot_cross_system_transfer_maskedmean.py @@ -0,0 +1,577 @@ +""" +========================================================== +Cross-system transfer with minimum distance to masked mean +========================================================== + +A cross-system benchmark asks whether a model trained on several systems can +generalize to a system it has never seen. A *system* is an acquisition setup: +an EEG device with its own montage and sampling rate. The source systems do +not necessarily share the same channels, neither with one another nor with the +target system, and they usually record different subjects. + +This tutorial introduces the Minimum Distance to Mean (MDM) [1]_ classifier +used with the Riemannian masked mean [2]_, which averages covariance matrices +that are only partially observed. We will: + +1. stack every system on the union of their channels; +2. describe the channels recorded by each system with a mask matrix; +3. estimate covariance matrices at the sampling rate of each system; +4. estimate class centroids with a masked mean of the covariance matrices; +5. compare this classifier with two simpler ways of handling missing channels. +""" + +# Author: Quentin Barthélemy +# +# License: BSD (3-clause) + +import matplotlib.pyplot as plt +import mne +import numpy as np +from matplotlib.colors import ListedColormap +from matplotlib.lines import Line2D +from pyriemann.classification import MDM +from pyriemann.estimation import Covariances +from pyriemann.geometry.distance import distance_riemann +from pyriemann.geometry.mean import maskedmean_riemann +from sklearn import config_context +from sklearn.base import BaseEstimator, TransformerMixin +from sklearn.pipeline import make_pipeline + +from moabb.datasets.fake import FakeDataset +from moabb.evaluations import CrossSystemEvaluation +from moabb.paradigms import LeftRightImagery + + +############################################################################### +# The cross-system transfer protocol +# ---------------------------------- +# +# We simulate three systems. The two source systems record seven channels +# each, but not the same ones. The target system records six channels: only +# Cz and Pz are recorded by both sources, while C3 and P3 are recorded by +# Source 1 only, and C4 and P4 by Source 2 only. The three systems also have +# different sampling rates. +# +# :class:`~moabb.evaluations.CrossSystemEvaluation` treats each dataset as a +# system: the pipeline is trained on all trials of the source systems and +# scored on each subject/session of the target system. It gives every trial as +# recorded, with the channels and the sampling rate of its system: handling +# heterogeneous systems is the job of the pipeline. + +systems = { + "Source 1": {"channels": ("F3", "Fz", "F4", "Cz", "C4", "P3", "Pz"), "sfreq": 128}, + "Source 2": {"channels": ("F3", "Fz", "F4", "C3", "Cz", "Pz", "P4"), "sfreq": 256}, + "Target": {"channels": ("P3", "Pz", "P4", "C3", "Cz", "C4"), "sfreq": 200}, +} +union = list( + dict.fromkeys(ch for system in systems.values() for ch in system["channels"]) +) +datasets = { + name: FakeDataset( + ["left_hand", "right_hand"], + channels=system["channels"], + sfreq=system["sfreq"], + n_subjects=n_subjects, + n_sessions=n_sessions, + seed=seed, + ) + for (name, system), n_subjects, n_sessions, seed in zip( + systems.items(), (3, 2, 3), (2, 3, 2), (42, 43, 44) + ) +} + +info = mne.create_info(union, sfreq=128, ch_types="eeg") +info.set_montage("standard_1020") + +PRESENT, ABSENT = "#2CA02C", "white" + + +def plot_channels(ax, groups, title): + """Plot the union channels on a head, one color per group of channels. + + ``groups`` is a list of ``(channel names, color, label)`` tuples. + """ + groups = [group for group in groups if len(group[0])] + mne.viz.plot_sensors( + info, + kind="topomap", + ch_groups=[[union.index(ch) for ch in channels] for channels, _, _ in groups], + cmap=ListedColormap([color for _, color, _ in groups]), + show_names=True, + axes=ax, + show=False, + pointsize=80, + linewidth=1, + ) + ax.set_title(title) + + +def channel_legend(fig, entries): + handles = [ + Line2D( + [], [], marker="o", linestyle="", markersize=9, markeredgecolor="k", color=c + ) + for _, c in entries + ] + fig.legend(handles, [label for label, _ in entries], loc="lower center", ncols=3) + + +fig, axes = plt.subplots(1, len(systems), figsize=(10, 4.2)) +fig.subplots_adjust(left=0.02, right=0.98, top=0.86, bottom=0.12, wspace=0.05) +for ax, (name, system) in zip(axes, systems.items()): + recorded = system["channels"] + missing = [ch for ch in union if ch not in recorded] + plot_channels( + ax, + [(recorded, PRESENT, "recorded"), (missing, ABSENT, "missing")], + f"{name}\n{len(recorded)} channels, {system['sfreq']} Hz", + ) +channel_legend(fig, [("recorded", PRESENT), ("missing", ABSENT)]) +plt.show() + + +############################################################################### +# Channel union and mask matrices +# ------------------------------- +# +# Masked covariance methods need a common channel space. Let :math:`p` be the +# number of channels in the union of all systems. The first step of the +# pipelines zero-fills a trial recorded by a system with :math:`p_s` channels +# on the missing channels, so that every trial is an array of :math:`p` +# channels. Its covariance matrix :math:`C` is therefore only meaningful on the +# recorded channels. +# +# The recorded channels of system :math:`s` are described by a semi-orthogonal +# mask matrix :math:`M_s \in \mathbb{R}^{p \times p_s}`, built from the columns +# of the identity matrix of the recorded channels. The masked covariance +# matrix :math:`M_s^\top C M_s` is the :math:`p_s \times p_s` covariance matrix +# of the recorded channels, which is SPD. +# +# The masks are metadata, known from the montages of the devices: starting +# from the channel union, the mask of a system is 1 on the channels it records +# and 0 on the other channels of the union. + +montages = {datasets[name].code: system["channels"] for name, system in systems.items()} +system_masks = {code: np.isin(union, channels) for code, channels in montages.items()} +target = datasets["Target"].code + +############################################################################### +# To place each training trial on the channel union, the first step of the +# pipelines needs the montage of its system: it requests the ``systems`` fit +# metadata, routed by MOABB through scikit-learn's metadata routing. At +# prediction, all trials come from the target system. + + +class ChannelUnion(TransformerMixin, BaseEstimator): + """Zero-fill the trials of each system on the channel union. + + ``montages`` maps each system to its channel names, and ``target`` is the + target system, whose montage is applied at prediction. + """ + + def __init__(self, montages=None, target=None): + self.montages = montages + self.target = target + + def fit(self, X, y=None, *, systems=None): + self.channels_ = list( + dict.fromkeys(ch for channels in self.montages.values() for ch in channels) + ) + return self + + def fit_transform(self, X, y=None, *, systems=None): + if systems is None: + raise ValueError("ChannelUnion needs `systems` metadata.") + return self.fit(X, y)._stack(X, systems) + + def transform(self, X): + return self._stack(X, [self.target] * len(X)) + + def _stack(self, X, systems): + X_union = np.empty(len(X), dtype=object) + for i, (trial, system) in enumerate(zip(X, systems)): + X_union[i] = np.zeros((len(self.channels_), trial.shape[-1])) + X_union[i][[self.channels_.index(ch) for ch in self.montages[system]]] = trial + return X_union + + +def channel_union(): + with config_context(enable_metadata_routing=True): + return ChannelUnion(montages, target).set_fit_request(systems=True) + + +############################################################################### +# Different sampling rates +# ------------------------ +# +# MOABB does not resample the systems: each trial is kept at the sampling rate +# of its system. Resampling is not needed either: covariance matrices only +# require that all trials cover windows of the same duration in seconds, and +# they have the same size whatever the number of time samples. +# +# As the systems have different channels and numbers of time samples, MOABB +# gives the trials as an object array, one array of shape +# :math:`(p_s, n_{times})` per trial, and so does the channel-union step, with +# arrays of shape :math:`(p, n_{times})`. pyRiemann's +# :class:`~pyriemann.estimation.Covariances` expects trials of the same length, +# so the estimator below applies it to the groups of trials of the same +# length. + + +class TrialCovariances(Covariances): + """Covariance matrices of trials with different numbers of time samples.""" + + def transform(self, X): + n_times = np.array([trial.shape[-1] for trial in X]) + covs = np.empty((len(X), X[0].shape[0], X[0].shape[0])) + for n in np.unique(n_times): + covs[n_times == n] = super().transform(np.stack(X[n_times == n])) + return covs + + +############################################################################### +# Minimum distance to masked mean +# ------------------------------- +# +# MDM represents each class :math:`k` by the Riemannian mean :math:`G_k` of +# its training covariance matrices. With partially observed matrices, the +# masked Riemannian mean [2]_ solves +# +# .. math:: +# +# G_k = \arg\min_{G} \sum_{i \in \mathcal{I}_k} +# d_R^2 \left( M_i^\top C_i M_i, M_i^\top G M_i \right), +# +# where :math:`d_R` is the affine-invariant Riemannian distance and +# :math:`M_i` is the mask of the system of trial :math:`i`. Every source +# system thus contributes to the prototype on the channels it records, and +# :math:`G_k` is defined on the complete channel union. +# +# A target trial :math:`C` is compared with the centroids on the target +# channels only, with the target mask :math:`M_t`: +# +# .. math:: +# +# \hat{k} = \arg\min_k d_R \left( M_t^\top C M_t, M_t^\top G_k M_t \right). +# +# The classifier receives the masks of the systems and the target system as +# parameters. To know the system of each training trial, it also requests the +# ``systems`` fit metadata. +# +# The covariance matrices are estimated with the sample covariance matrix: +# a masked block of the covariance of zero-filled data is then exactly the +# covariance matrix of the recorded channels. A shrinkage estimator would mix +# in the zero-filled channels. + + +class MDMM(MDM): + """Minimum distance to masked Riemannian mean. + + ``masks`` maps each system to its boolean mask on the channel union, and + ``target`` is the target system, whose mask is applied at prediction. + """ + + def __init__(self, masks=None, target=None, metric="riemann", n_jobs=1): + super().__init__(metric=metric, n_jobs=n_jobs) + self.masks = masks + self.target = target + + def _mask(self, system): + """Semi-orthogonal mask matrix of a system.""" + return np.eye(len(self.masks[system]))[:, self.masks[system]] + + def fit(self, X, y, *, systems=None): + if systems is None: + raise ValueError("MDMM needs `systems` metadata.") + X, y = np.asarray(X), np.asarray(y) + self.classes_ = np.unique(y) + masks = [self._mask(system) for system in systems] + self.covmeans_ = np.stack( + [ + maskedmean_riemann(X[y == c], [m for m, yi in zip(masks, y) if yi == c]) + for c in self.classes_ + ] + ) + return self + + def _predict_distances(self, X): + mask = self._mask(self.target) + maskedX = mask.T @ np.asarray(X) @ mask + return np.stack( + [distance_riemann(maskedX, mask.T @ G @ mask) for G in self.covmeans_], axis=1 + ) + + +############################################################################### +# Two simpler baselines +# --------------------- +# +# The first baseline ignores the masks: the MDM is trained on the covariance +# matrices of zero-filled data, estimated with a shrinkage estimator to keep +# them SPD. The centroids then mix systems whose zero-filled channels differ. +# +# The second baseline keeps only the channels recorded by every source and by +# the target system: here, only Cz and Pz remain. It discards most of the +# spatial information, in particular the lateral C3 and C4 channels. The +# montages of the systems are known in advance, so the common channels are +# given to the pipeline as a parameter. +# +# Note that none of these three methods applies a recentering of the +# covariances by subject. + + +class ChannelSelector(TransformerMixin, BaseEstimator): + """Keep a subset of the channels of covariance matrices.""" + + def __init__(self, picks=None): + self.picks = picks + + def fit(self, X, y=None): + return self + + def transform(self, X): + return np.asarray(X)[:, self.picks][:, :, self.picks] + + +common = np.flatnonzero(np.all(list(system_masks.values()), axis=0)) + +with config_context(enable_metadata_routing=True): + mdmm = MDMM(masks=system_masks, target=target).set_fit_request(systems=True) + +pipelines = { + "Zero-filled MDM": make_pipeline(channel_union(), TrialCovariances("oas"), MDM()), + "Intersection MDM": make_pipeline( + channel_union(), TrialCovariances("scm"), ChannelSelector(common), MDM() + ), + "Masked-mean MDM": make_pipeline(channel_union(), TrialCovariances("scm"), mdmm), +} + + +############################################################################### +# Channels used by each cross-system classifier +# --------------------------------------------- +# +# When the target system is held out, the three classifiers do not use the +# same channels, neither to train on the source systems nor to predict the +# target trials: +# +# * the zero-filled MDM uses every channel of the union, including channels +# that are zero-filled because a source, or the target, does not record them; +# * the intersection MDM uses only the channels recorded by all systems; +# * the masked-mean MDM estimates its centroids on the whole union, each +# source contributing on the channels it records, and predicts on the +# channels recorded by the target system. It never uses zero-filled values. + +ZERO_FILLED = "#F2CF5B" +source_masks = np.array( + [system_masks[datasets[name].code] for name in ("Source 1", "Source 2")] +) +target_mask = system_masks[target] + + +def channel_groups(used, zero_filled): + """Split the union channels into recorded, zero-filled and unused ones.""" + used, zero_filled = np.asarray(used), np.asarray(zero_filled) + names = np.asarray(union) + return [ + (names[used & ~zero_filled], PRESENT, "used"), + (names[used & zero_filled], ZERO_FILLED, "used, zero-filled"), + (names[~used], ABSENT, "unused"), + ] + + +everywhere = np.ones(len(union), dtype=bool) +nowhere = np.zeros(len(union), dtype=bool) +channel_usage = { + "Zero-filled MDM": ( + channel_groups(everywhere, ~source_masks.all(axis=0)), + channel_groups(everywhere, ~target_mask), + ), + "Intersection MDM": ( + channel_groups(source_masks.all(axis=0) & target_mask, nowhere), + channel_groups(source_masks.all(axis=0) & target_mask, nowhere), + ), + "Masked-mean MDM": ( + channel_groups(source_masks.any(axis=0), nowhere), + channel_groups(target_mask, nowhere), + ), +} + +fig, axes = plt.subplots(2, len(channel_usage), figsize=(10, 8)) +fig.subplots_adjust( + left=0.02, right=0.98, top=0.88, bottom=0.07, wspace=0.05, hspace=0.25 +) +for col, (pipeline, (train_groups, test_groups)) in enumerate(channel_usage.items()): + plot_channels(axes[0, col], train_groups, f"{pipeline}\ntraining on sources") + plot_channels(axes[1, col], test_groups, "prediction on target") +fig.suptitle("Channels used by each classifier when the target system is held out") +channel_legend( + fig, [("used", PRESENT), ("used, zero-filled", ZERO_FILLED), ("unused", ABSENT)] +) +plt.show() + + +############################################################################### +# Run the cross-system benchmark +# ------------------------------ +# +# The systems are deterministic :class:`~moabb.datasets.fake.FakeDataset` +# instances, which keep the tutorial fast and download-free. Their scores have +# no scientific meaning; they are used to expose the complete evaluation path. +# +# The three systems share the same epoch duration of 3 seconds. +# ``target_datasets`` sets the target system, the other systems being the +# sources; by default, each system is used as target in turn. ``cs_mode`` sets +# the target data available during training: by default, the pipeline sees no +# target trial. + +paradigm = LeftRightImagery() + +evaluation = CrossSystemEvaluation( + paradigm=paradigm, + datasets=list(datasets.values()), + target_datasets=[datasets["Target"]], + overwrite=True, + suffix="cross_system_maskedmean", +) +results = evaluation.process(pipelines) + +print( + results.groupby("pipeline")["score"] + .agg(["mean", "std", "count"]) + .sort_values("mean", ascending=False) + .round(3) +) + + +############################################################################### +# Inspect paired target results +# ----------------------------- +# +# Each thin line below connects the three pipelines on the same target subject +# and session. The diamonds show pipeline means. +# +# Do not interpret the ranking produced by ``FakeDataset``: its labels contain +# no physiological effect, so all methods fluctuate around chance. A proper +# study would repeat the benchmark over real datasets recorded with different +# devices, and use MOABB's statistical-analysis utilities. + +paired = results.pivot(index=["subject", "session"], columns="pipeline", values="score") +order = list(pipelines) +scores = paired[order].dropna().to_numpy() +colors = ["#4C78A8", "#F58518", "#54A24B"] +positions = np.arange(len(order)) +offsets = np.linspace(-0.06, 0.06, len(scores)) + +fig, ax = plt.subplots(figsize=(7.5, 5)) +for row, offset in zip(scores, offsets): + ax.plot(positions + offset, row, color="0.75", linewidth=1, zorder=1) + +for position, color, values in zip(positions, colors, scores.T): + 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=positions, + xticklabels=order, + ylabel="ROC AUC", + title="Paired target-subject/session scores on the unseen system", +) +ax.grid(axis="y", alpha=0.25) +ax.legend() +fig.tight_layout() +plt.show() + + +############################################################################### +# Masked-mean centroids and system masks +# -------------------------------------- +# +# To look at the class centroids :math:`G_k` of the masked-mean MDM, we fit it +# on the two source systems, as in the benchmark: the trials of each system +# are zero-filled on the channel union, and their covariance matrices are +# averaged with the mask matrix :math:`M_s` of their system. +# +# The centroids are defined on the complete channel union, although no source +# records all of its channels. Each system only sees the block +# :math:`M_s^\top G_k M_s` selected by its mask: the target system is +# classified on its own channels. + +trials, labels, trial_systems = [], [], [] +for name in ("Source 1", "Source 2"): + epochs, y, _ = paradigm.get_data(datasets[name], return_epochs=True) + trials += list(epochs.get_data() * datasets[name].unit_factor) + labels += list(y) + trial_systems += [datasets[name].code] * len(y) +X = np.empty(len(trials), dtype=object) +for i, trial in enumerate(trials): + X[i] = trial + +X_union = ChannelUnion(montages, target).fit_transform(X, systems=trial_systems) +covs = TrialCovariances("scm").transform(X_union) +centroids = ( + MDMM(masks=system_masks, target=target) + .fit(covs, labels, systems=trial_systems) + .covmeans_ +) + +fig = plt.figure(figsize=(10, 9)) +grid = fig.add_gridspec(2, 6, height_ratios=(1.1, 1), hspace=0.45, wspace=0.9) +vmax = np.abs(centroids).max() +for col, (klass, centroid) in enumerate(zip(np.unique(labels), centroids)): + ax = fig.add_subplot(grid[1, 3 * col : 3 * col + 3]) + image = ax.imshow(centroid, cmap="RdBu_r", vmin=-vmax, vmax=vmax) + ax.set( + xticks=range(len(union)), + yticks=range(len(union)), + xticklabels=union, + yticklabels=union, + title=f"Centroid of class {klass}", + ) +fig.colorbar(image, ax=fig.axes, location="right", shrink=0.8, label="Covariance") + +for col, name in enumerate(systems): + mask = system_masks[datasets[name].code] + ax = fig.add_subplot(grid[0, 2 * col : 2 * col + 2]) + ax.imshow( + np.eye(len(union))[:, mask], + cmap=ListedColormap([ABSENT, PRESENT]), + vmin=0, + vmax=1, + aspect="auto", + ) + ax.set( + xticks=range(mask.sum()), + yticks=range(len(union)), + xticklabels=np.asarray(union)[mask], + yticklabels=union, + title=f"{name}: mask $M_s$ ({len(union)} x {mask.sum()})", + ) + ax.tick_params(axis="x", labelrotation=90) + ax.set_xticks(np.arange(mask.sum() + 1) - 0.5, minor=True) + ax.set_yticks(np.arange(len(union) + 1) - 0.5, minor=True) + ax.grid(which="minor", color="0.75", linewidth=0.8) + ax.tick_params(which="minor", length=0) +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 +# .. [2] Yger, F., Chevallier, S., Barthélemy, Q., & Sra, S. (2020). +# Geodesically-convex optimization for averaging partially observed +# covariance matrices. *Asian Conference on Machine Learning (ACML)*, +# 417-432. https://hal.archives-ouvertes.fr/hal-02984423 diff --git a/moabb/evaluations/__init__.py b/moabb/evaluations/__init__.py index 0bb2b8c1e..274861ba3 100644 --- a/moabb/evaluations/__init__.py +++ b/moabb/evaluations/__init__.py @@ -7,6 +7,7 @@ from .evaluations import ( CrossSessionEvaluation, CrossSubjectEvaluation, + CrossSystemEvaluation, WithinSessionEvaluation, WithinSubjectEvaluation, ) diff --git a/moabb/evaluations/base.py b/moabb/evaluations/base.py index 968ddebb9..ad159657f 100644 --- a/moabb/evaluations/base.py +++ b/moabb/evaluations/base.py @@ -163,6 +163,7 @@ def _evaluate_fold( split_metadata=None, calib_idx=None, score_subjects=None, + extra_fit_params=None, ): """Evaluate a single CV fold. Pure function, no shared mutable state. @@ -194,6 +195,10 @@ def _evaluate_fold( Held-out subjects whose results are still missing for this pipeline. When provided, the fold is fitted once and scored separately for these subjects. + extra_fit_params : dict | None + Evaluation-specific fit metadata (e.g. the systems of the training + trials of a cross-system fold). It is merged into the transfer fit + metadata and routed only to the pipeline steps that request it. Returns ------- @@ -262,6 +267,9 @@ def _evaluate_fold( else: fit_params["X_target_labeled"] = X_calib fit_params["y_target_labeled"] = y_calib + if extra_fit_params: + fit_params.update(extra_fit_params) + if fit_params: with config_context(enable_metadata_routing=True): requested = get_routing_for_object(cvclf).consumes("fit", set(fit_params)) fit_params = { @@ -1397,11 +1405,11 @@ def _process_parallel(self, pipelines, param_grid, postprocess_pipeline): self._log_result(res) self._push_results_batch(batch, pipelines, process_pipeline) - res_per_db.append( - self.results.to_dataframe( - pipelines=pipelines, process_pipeline=process_pipeline - ) + # Datasets sharing a pipeline digest share their results group. + df = self.results.to_dataframe( + pipelines=pipelines, process_pipeline=process_pipeline ) + res_per_db.append(df[df["dataset"] == dataset.code]) return pd.concat(res_per_db, ignore_index=True) diff --git a/moabb/evaluations/evaluations.py b/moabb/evaluations/evaluations.py index 344cab344..c43603809 100644 --- a/moabb/evaluations/evaluations.py +++ b/moabb/evaluations/evaluations.py @@ -2,6 +2,7 @@ from typing import TYPE_CHECKING, Optional import numpy as np +import pandas as pd from sklearn.base import clone from sklearn.model_selection import ( BaseCrossValidator, @@ -750,6 +751,106 @@ def _get_incompatibility_reason(self, dataset): return "requirements not met" +class CrossSystemEvaluation(CrossSubjectEvaluation): + """Cross-system evaluation performance. + + A *system* is an acquisition setup: an EEG device with its own montage and + sampling rate, which strongly influences EEG data [1]_. Each dataset is one + system, and systems do not need to share their channels or sampling rate. + For each target system, the pipeline is trained on all other (source) + systems and scored on each subject/session of the target system, a + cross-device transfer learning setting [2]_. + + Each trial is kept as recorded by its system, with its own channels and + sampling rate: handling heterogeneous systems is left to the pipelines. + When all systems share the same channels and number of time samples, ``X`` + is an array of shape ``(n_trials, n_channels, n_times)``; otherwise, it is + an object array of shape ``(n_trials,)`` whose elements are arrays of shape + ``(n_channels_s, n_times_s)`` of their system :math:`s`. + + In addition to the fit metadata of + :class:`~moabb.evaluations.CrossSubjectEvaluation`, pipelines can request, + through scikit-learn's metadata routing (``set_fit_request``), ``systems``: + the system (dataset code) of each training trial. + + Parameters + ---------- + target_datasets : list of :class:`~moabb.datasets.base.BaseDataset` or None + Systems used as targets, each one in turn, among ``datasets``. If None, + every system is used as target (leave-one-system-out). + Defaults to ``None``. + **kwargs + Parameters of :class:`~moabb.evaluations.CrossSubjectEvaluation`. + + Notes + ----- + .. versionadded:: 1.8.0 + + References + ---------- + .. [1] Melnik, A., Legkov, P., Izdebski, K., Kärcher, S. M., Hairston, + W. D., Ferris, D. P., & König, P. (2017). Systems, Subjects, + Sessions: To What Extent Do These Factors Influence EEG Data? + *Frontiers in Human Neuroscience*, 11, 150. + https://doi.org/10.3389/fnhum.2017.00150 + .. [2] Wu, D., Xu, Y., & Lu, B.-L. (2020). Transfer Learning for EEG-Based + Brain-Computer Interfaces: A Review of Progress Made Since 2016. + *IEEE Transactions on Cognitive and Developmental Systems*, 14(1), + 4-19. https://doi.org/10.1109/TCDS.2020.3007453 + """ + + _eval_type = "CrossSystem" + + def __init__(self, *args, target_datasets=None, **kwargs): + super().__init__(*args, **kwargs) + self.groups = "system" + # All systems are loaded as sources, only targets are evaluated. + self.systems = self.datasets + if target_datasets is not None: + self.datasets = list(target_datasets) + + def is_valid(self, dataset): + return True + + def _load_data( + self, dataset, run_pipes, process_pipeline, postprocess_pipeline, subjects=None + ): + """Concatenate the trials of all systems, as recorded.""" + if postprocess_pipeline is not None: + raise ValueError( + "CrossSystemEvaluation does not support postprocess_pipeline." + ) + systems = [] + for system in self.systems: + epochs, y, metadata = self.paradigm.get_data( + system, return_epochs=True, cache_config=self.cache_config + ) + X = epochs.get_data() * system.unit_factor + systems.append((epochs.ch_names, X, y, metadata.assign(system=system.code))) + + if len({(tuple(ch_names), X.shape[-1]) for ch_names, X, _, _ in systems}) == 1: + X_all = np.concatenate([X for _, X, _, _ in systems]) + else: + trials = [trial for _, X, _, _ in systems for trial in X] + X_all = np.empty(len(trials), dtype=object) + for i, trial in enumerate(trials): + X_all[i] = trial + y_all = np.concatenate([y for _, _, y, _ in systems]) + metadata = pd.concat([metadata for *_, metadata in systems], ignore_index=True) + return X_all, y_all, metadata + + def _build_task_list(self, dataset, X, y, metadata, *args): + """Keep the fold holding out ``dataset`` and route the systems.""" + tasks = super()._build_task_list(dataset, X, y, metadata, *args) + systems = metadata["system"].to_numpy() + tasks = [ + task for task in tasks if np.all(systems[task["test_idx"]] == dataset.code) + ] + for task in tasks: + task["extra_fit_params"] = {"systems": systems[task["train_idx"]]} + return tasks + + class WithinSubjectEvaluation(BaseEvaluation): """Within-subject k-fold cross-validation pooling all sessions. diff --git a/moabb/evaluations/utils.py b/moabb/evaluations/utils.py index 7791c4e53..ca4fd7492 100644 --- a/moabb/evaluations/utils.py +++ b/moabb/evaluations/utils.py @@ -547,10 +547,13 @@ def _pipeline_requires_epochs(pipeline): def _get_nchan(X): - """Extract number of channels from data (Epochs or ndarray).""" + """Extract number of channels from data (Epochs, ndarray or object ndarray).""" from mne.epochs import BaseEpochs - return X.info["nchan"] if isinstance(X, BaseEpochs) else X.shape[1] + if isinstance(X, BaseEpochs): + return X.info["nchan"] + # An object array holds one (n_channels, n_times) array per trial. + return X[0].shape[0] if X.dtype == object else X.shape[1] class Emissions: diff --git a/moabb/tests/test_cross_system_evaluation.py b/moabb/tests/test_cross_system_evaluation.py new file mode 100644 index 000000000..9bfe63c17 --- /dev/null +++ b/moabb/tests/test_cross_system_evaluation.py @@ -0,0 +1,225 @@ +"""Tests for the cross-system CrossSystemEvaluation.""" + +import numpy as np +import pandas as pd +import pytest +import sklearn.base +from pyriemann.estimation import Covariances +from pyriemann.spatialfilters import CSP +from sklearn import config_context +from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA +from sklearn.pipeline import FunctionTransformer, make_pipeline + +from moabb.datasets.fake import FakeDataset +from moabb.evaluations import CrossSubjectMode, CrossSystemEvaluation +from moabb.paradigms import LeftRightImagery + + +_CAPTURE = [] + + +class _SystemRecorder(sklearn.base.TransformerMixin, sklearn.base.BaseEstimator): + """A system-aware step that records the cross-system metadata it is routed.""" + + def fit(self, X, y=None, systems=None, X_target_unlabeled=None): + _CAPTURE.append( + {"X": X, "systems": systems, "X_target_unlabeled": X_target_unlabeled} + ) + return self + + def transform(self, X): + return X + + +def _trial_variances(X): + """One feature per trial, whatever its channels and time samples.""" + return np.array([[np.var(trial)] for trial in X]) + + +@pytest.fixture(autouse=True) +def _clear_capture(): + _CAPTURE.clear() + yield + _CAPTURE.clear() + + +def _datasets(sfreqs=(128, 128), channels=(("C3", "Cz", "C4"), ("Cz", "C4", "Pz", "P4"))): + """Two fake systems, by default with different channels.""" + return [ + FakeDataset( + ["left_hand", "right_hand"], + channels=channels[0], + n_subjects=2, + n_sessions=1, + sfreq=sfreqs[0], + seed=1, + ), + FakeDataset( + ["left_hand", "right_hand"], + channels=channels[1], + n_subjects=3, + n_sessions=2, + sfreq=sfreqs[1], + seed=2, + ), + ] + + +def _recorder_pipeline(): + with config_context(enable_metadata_routing=True): + step = _SystemRecorder().set_fit_request(systems=True, X_target_unlabeled=True) + return make_pipeline(step, FunctionTransformer(_trial_variances), LDA()) + + +def _evaluation(tmp_path, datasets, **kwargs): + kwargs.setdefault("paradigm", LeftRightImagery()) + kwargs.setdefault("overwrite", True) + return CrossSystemEvaluation(datasets=datasets, hdf5_path=tmp_path, **kwargs) + + +def _source_and_target(capture, datasets): + (source,) = np.unique(capture["systems"]) + (target,) = {dataset.code for dataset in datasets} - {source} + return source, target + + +def _identity(X): + return X + + +def test_trials_keep_the_channels_of_their_system(tmp_path): + ds_a, ds_b = _datasets() + evaluation = _evaluation(tmp_path, [ds_a, ds_b]) + X, y, metadata = evaluation._load_data(ds_a, {}, None, None) + + # Systems with different channels give an object array of trials. + assert X.dtype == object and X.shape == (len(y),) + n_channels = {ds_a.code: 3, ds_b.code: 4} + for trial, system in zip(X, metadata["system"]): + assert trial.shape == (n_channels[system], 3 * 128 + 1) + + results = evaluation.process({"rec": _recorder_pipeline()}) + # One row per target subject/session: 2 x 1 for A, 3 x 2 for B. + assert results.groupby("dataset").size().to_dict() == {ds_a.code: 2, ds_b.code: 6} + assert len(_CAPTURE) == 2 + for capture in _CAPTURE: + source, _ = _source_and_target(capture, [ds_a, ds_b]) + # Sources are trained on their own trials, routed with their system. + assert {trial.shape[0] for trial in capture["X"]} == {n_channels[source]} + assert len(capture["systems"]) == len(capture["X"]) + # Default TRAIN mode: no target trial is routed. + assert capture["X_target_unlabeled"] is None + + +def test_trials_keep_the_sampling_rate_of_their_system(tmp_path): + datasets = _datasets(sfreqs=(128, 256), channels=(("C3", "Cz", "C4"),) * 2) + evaluation = _evaluation(tmp_path, datasets) + X, _, metadata = evaluation._load_data(datasets[0], {}, None, None) + + assert X.dtype == object + n_times = {datasets[0].code: 3 * 128 + 1, datasets[1].code: 3 * 256 + 1} + for trial, system in zip(X, metadata["system"]): + assert trial.shape == (3, n_times[system]) + assert len(evaluation.process({"rec": _recorder_pipeline()})) == 8 + + +def test_same_channels_and_sampling_rate_give_regular_array(tmp_path): + datasets = _datasets(channels=(("C3", "Cz", "C4"),) * 2) + evaluation = _evaluation(tmp_path, datasets) + X, y, _ = evaluation._load_data(datasets[0], {}, None, None) + assert X.dtype != object and X.shape == (len(y), 3, 3 * 128 + 1) + + # A standard pipeline then runs as is. + results = evaluation.process( + {"csp": make_pipeline(Covariances("oas"), CSP(2), LDA())} + ) + assert len(results) == 8 + assert results["score"].between(0, 1).all() + + +def test_calibration_routes_target_slice(tmp_path): + ds_a, ds_b = _datasets() + results = _evaluation( + tmp_path, [ds_a, ds_b], cs_mode=CrossSubjectMode.TRAIN_AND_TARGET_UNLABELED_20P + ).process({"rec": _recorder_pipeline()}) + + (capture,) = [ + c for c in _CAPTURE if _source_and_target(c, [ds_a, ds_b])[1] == ds_a.code + ] + X_target = capture["X_target_unlabeled"] + assert X_target is not None and len(X_target) > 0 + # Target trials keep the channels of the target system. + assert {trial.shape[0] for trial in X_target} == {3} + # Calibration and scored target trials partition the target system. + _, y_a, _ = LeftRightImagery().get_data(ds_a) + scored = results.loc[results["dataset"] == ds_a.code, "samples_test"].sum() + assert scored + len(X_target) == len(y_a) + + +def test_target_datasets_restricts_targets(tmp_path): + ds_a, ds_b = _datasets() + ds_c = FakeDataset( + ["left_hand", "right_hand"], + channels=("C3", "Pz"), + n_subjects=2, + n_sessions=1, + seed=3, + ) + results = _evaluation(tmp_path, [ds_a, ds_b, ds_c], target_datasets=[ds_c]).process( + {"rec": _recorder_pipeline()} + ) + + assert set(results["dataset"]) == {ds_c.code} + assert len(results) == 2 + (capture,) = _CAPTURE + # The non-target systems are the sources. + assert set(capture["systems"]) == {ds_a.code, ds_b.code} + + +def test_trialwise_mode(tmp_path): + evaluation = _evaluation( + tmp_path, _datasets(), cs_mode=CrossSubjectMode.TRAIN_TRIALWISE + ) + assert len(evaluation.process({"rec": _recorder_pipeline()})) == 8 + + +def test_results_are_cached(tmp_path): + datasets = _datasets() + pipelines = {"rec": _recorder_pipeline()} + first = _evaluation(tmp_path, datasets).process(pipelines) + n_fits = len(_CAPTURE) + + cached = _evaluation(tmp_path, datasets, overwrite=False).process(pipelines) + + assert len(_CAPTURE) == n_fits + keys = ["dataset", "subject", "session"] + pd.testing.assert_frame_equal( + first.sort_values(keys).reset_index(drop=True), + cached.sort_values(keys).reset_index(drop=True), + check_like=True, + ) + + +def test_evaluate_single_target(tmp_path): + ds_a, ds_b = _datasets() + evaluation = _evaluation(tmp_path, [ds_a, ds_b]) + process_pipeline = evaluation.paradigm.make_process_pipelines(ds_b)[0] + results = list( + evaluation.evaluate( + ds_b, + {"rec": _recorder_pipeline()}, + param_grid=None, + process_pipeline=process_pipeline, + ) + ) + assert len(results) == 6 + assert {res["dataset"].code for res in results} == {ds_b.code} + + +def test_rejects_postprocess_pipeline(tmp_path): + evaluation = _evaluation(tmp_path, _datasets()) + with pytest.raises(ValueError, match="postprocess_pipeline"): + evaluation.process( + {"rec": _recorder_pipeline()}, + postprocess_pipeline=FunctionTransformer(_identity), + ) diff --git a/moabb/tests/test_evaluations.py b/moabb/tests/test_evaluations.py index 3fc6407bb..08c5440be 100644 --- a/moabb/tests/test_evaluations.py +++ b/moabb/tests/test_evaluations.py @@ -44,6 +44,11 @@ def _identity(x): return x +def _trial_variances(X): + """One feature per trial, whatever its channels and time samples.""" + return np.array([[np.var(trial)] for trial in X]) + + class OpaqueLeaveOneGroupOut(LeaveOneGroupOut): """Real group CV whose constructor signature cannot be inspected.""" @@ -1595,6 +1600,15 @@ def test_incompatibility_error_message(self): class TestUtilEvaluation: + def test_get_nchan(self): + from moabb.evaluations.utils import _get_nchan + + assert _get_nchan(np.zeros((4, 3, 10))) == 3 + # Trials of different lengths, as given by CrossSystemEvaluation. + trials = np.empty(2, dtype=object) + trials[0], trials[1] = np.zeros((3, 10)), np.zeros((3, 20)) + assert _get_nchan(trials) == 3 + def test_save_model_cv(self): model = Dummy() save_path = "test_save_path" @@ -1997,6 +2011,52 @@ def test_cross_subject_process_structure(self, tmp_path): # 2 subjects × 2 sessions = 4 results assert len(results) == 4 + def test_cross_system_process_structure(self, tmp_path): + """CrossSystem process() returns one result per target subject/session.""" + systems = [ + FakeDataset( + ["left_hand", "right_hand"], channels=channels, n_subjects=2, seed=seed + ) + for channels, seed in ((("C3", "Cz", "C4"), 1), (("Cz", "C4", "Pz"), 2)) + ] + evaluation = ev.CrossSystemEvaluation( + paradigm=FakeImageryParadigm(), + datasets=systems, + overwrite=True, + hdf5_path=str(tmp_path / "parallel_test"), + ) + # Trials keep the channels of their system, so the pipeline must not + # depend on them. + results = evaluation.process( + {"V": make_pipeline(FunctionTransformer(_trial_variances), LDA())} + ) + # 2 systems as targets x 2 subjects x 2 sessions = 8 results + assert len(results) == 8 + assert set(results["dataset"]) == {system.code for system in systems} + + def test_process_keeps_rows_of_each_dataset_once(self, tmp_path): + """Datasets sharing a pipeline digest do not duplicate their results.""" + # Same parameters except the channels: same process pipeline digest. + datasets = [ + FakeDataset( + ["left_hand", "right_hand"], channels=channels, n_subjects=2, seed=3 + ) + for channels in (("C3", "Cz", "C4"), ("C3", "Cz", "C4", "Pz")) + ] + evaluation = ev.CrossSubjectEvaluation( + paradigm=FakeImageryParadigm(), + datasets=datasets, + overwrite=True, + hdf5_path=str(tmp_path / "parallel_test"), + ) + results = evaluation.process( + {"C": make_pipeline(Covariances("oas"), CSP(2), LDA())} + ) + # 2 datasets x 2 subjects x 2 sessions = 8 results, each once + keys = ["dataset", "subject", "session", "pipeline"] + assert len(results) == 8 + assert not results.duplicated(keys).any() + def test_cross_subject_multisubject_fold_preserves_provenance(self, tmp_path): """A fold holding out multiple subjects is scored per subject/session.""" ds = FakeDataset(["left_hand", "right_hand"], n_subjects=4, n_sessions=2, seed=23) @@ -3004,3 +3064,66 @@ def test_cross_subject_cv_class_with_reserved_metadata_columns_saves(tmp_path): assert set(results["subject"]) == {"2"} assert set(results["session"]) == {"1"} + + +# --------------------------------------------------------------------------- +# Evaluation-specific fit metadata routed by _evaluate_fold. +# --------------------------------------------------------------------------- + +_EXTRA_CAPTURE = [] + + +class _ExtraRecorder(sklearn.base.TransformerMixin, sklearn.base.BaseEstimator): + """Records the evaluation-specific metadata it is routed.""" + + def fit(self, X, y=None, systems=None): + _EXTRA_CAPTURE.append(systems) + return self + + def transform(self, X): + return X + + +@pytest.mark.parametrize("requested", [True, False]) +def test_evaluate_fold_routes_extra_fit_params(requested): + """extra_fit_params reach only the steps requesting them, without a + transfer split, and are ignored by the other pipelines.""" + from sklearn import config_context + + from moabb.evaluations.base import _evaluate_fold + + _EXTRA_CAPTURE.clear() + with config_context(enable_metadata_routing=True): + step = _ExtraRecorder().set_fit_request(systems=requested) + pipe = make_pipeline(step, Covariances("oas"), CSP(2), LDA()) + + evaluation = ev.CrossSessionEvaluation( + paradigm=FakeImageryParadigm(), datasets=[dataset], overwrite=True + ) + X, y, metadata = FakeImageryParadigm().get_data(dataset, subjects=[1]) + train = np.flatnonzero(metadata["session"] == metadata["session"].iloc[0]) + test = np.flatnonzero(metadata["session"] != metadata["session"].iloc[0]) + systems = np.array(["system"] * len(y)) + + results = _evaluate_fold( + X, + y, + metadata, + evaluation._build_eval_config(None), + dataset, + "P", + pipe, + train, + test, + 1, + metadata["session"].iloc[-1], + 0, + extra_fit_params={"systems": systems[train]}, + ) + + assert len(results) == 1 + (routed,) = _EXTRA_CAPTURE + if requested: + np.testing.assert_array_equal(routed, systems[train]) + else: + assert routed is None From 9f0483186e707ca06ce46c1c0f6d67bd0ddf1526 Mon Sep 17 00:00:00 2001 From: Bru Date: Sat, 24 Oct 2026 09:06:41 -0700 Subject: [PATCH 2/3] CrossSystemEvaluation: load the systems once, reject unknown targets and epoch pipelines Review follow-ups: the systems were reloaded in full for every target (k^2 loads), a target_datasets entry outside datasets gave an empty result, and pipelines needing mne.Epochs silently received arrays. --- moabb/evaluations/evaluations.py | 20 ++++++++++++++++++++ moabb/tests/test_cross_system_evaluation.py | 7 +++++++ 2 files changed, 27 insertions(+) diff --git a/moabb/evaluations/evaluations.py b/moabb/evaluations/evaluations.py index c43603809..01c158c5c 100644 --- a/moabb/evaluations/evaluations.py +++ b/moabb/evaluations/evaluations.py @@ -33,6 +33,7 @@ _average_scores, _carbonfootprint, _create_scorer, + _pipeline_requires_epochs, _update_result_with_scores, ) @@ -806,7 +807,15 @@ def __init__(self, *args, target_datasets=None, **kwargs): self.groups = "system" # All systems are loaded as sources, only targets are evaluated. self.systems = self.datasets + self._trials = None if target_datasets is not None: + codes = {system.code for system in self.systems} + unknown = [d.code for d in target_datasets if d.code not in codes] + if unknown: + raise ValueError( + f"target_datasets {unknown} are not among the valid datasets " + f"{sorted(codes)}." + ) self.datasets = list(target_datasets) def is_valid(self, dataset): @@ -820,6 +829,17 @@ def _load_data( raise ValueError( "CrossSystemEvaluation does not support postprocess_pipeline." ) + if any(_pipeline_requires_epochs(clf) for clf in run_pipes.values()): + raise ValueError( + "CrossSystemEvaluation gives pipelines arrays, not mne.Epochs." + ) + # ponytail: every target needs the same trials of all systems, so load + # them once per evaluation (held in memory for the whole run). + if self._trials is None: + self._trials = self._load_systems() + return self._trials + + def _load_systems(self): systems = [] for system in self.systems: epochs, y, metadata = self.paradigm.get_data( diff --git a/moabb/tests/test_cross_system_evaluation.py b/moabb/tests/test_cross_system_evaluation.py index 9bfe63c17..4e13b5755 100644 --- a/moabb/tests/test_cross_system_evaluation.py +++ b/moabb/tests/test_cross_system_evaluation.py @@ -223,3 +223,10 @@ def test_rejects_postprocess_pipeline(tmp_path): {"rec": _recorder_pipeline()}, postprocess_pipeline=FunctionTransformer(_identity), ) + + +def test_unknown_target_dataset_raises(tmp_path): + ds_a, ds_b = _datasets() + other = FakeDataset(["left_hand", "right_hand"], channels=("C3", "Pz"), n_subjects=2) + with pytest.raises(ValueError, match="not among the valid datasets"): + _evaluation(tmp_path, [ds_a, ds_b], target_datasets=[other]) From bf93c3d183420d9a293d959d9f21682950819f44 Mon Sep 17 00:00:00 2001 From: Bru Date: Mon, 26 Oct 2026 00:47:49 +0100 Subject: [PATCH 3/3] Turn CrossSystemEvaluation into a tutorial on the existing API Review of #1305 (Bruno): no new evaluation. The example now pools the devices in a CompoundDataset and holds out one device at a time with CrossSubjectEvaluation(cv_class=LeaveOneGroupOut, groups='device'), on the nested trials of #1327 and the device column of #1328. - remove CrossSystemEvaluation, its tests and extra_fit_params; - _evaluate_fold routes any metadata column of the training trials to the steps requesting it (set_fit_request(device=True)), replacing the cross-system-only 'systems' routing; - FakeDataset takes a device; - the example's ChannelUnion and MDMM read the device of each training trial and infer the held-out device, so every device can be the target; - keep the duplicated-results fix of process() and its test. --- docs/source/api.rst | 1 - docs/source/whats_new.rst | 3 +- .../plot_cross_system_transfer_maskedmean.py | 145 ++++++----- moabb/datasets/fake.py | 8 + moabb/evaluations/__init__.py | 1 - moabb/evaluations/base.py | 12 +- moabb/evaluations/evaluations.py | 121 --------- moabb/tests/test_cross_system_evaluation.py | 232 ------------------ moabb/tests/test_evaluations.py | 88 ++++--- 9 files changed, 149 insertions(+), 462 deletions(-) delete mode 100644 moabb/tests/test_cross_system_evaluation.py diff --git a/docs/source/api.rst b/docs/source/api.rst index 2fee831f3..5cf1778d1 100644 --- a/docs/source/api.rst +++ b/docs/source/api.rst @@ -535,7 +535,6 @@ accuracy, across-subject accuracy, or other transfer learning settings. WithinSessionEvaluation CrossSessionEvaluation CrossSubjectEvaluation - CrossSystemEvaluation .. autosummary:: :toctree: generated/ diff --git a/docs/source/whats_new.rst b/docs/source/whats_new.rst index c63e4f089..218c9d3f0 100644 --- a/docs/source/whats_new.rst +++ b/docs/source/whats_new.rst @@ -22,7 +22,7 @@ Version 1.8.0 (2026-10-06) Enhancements ~~~~~~~~~~~~ -- Add :class:`moabb.evaluations.CrossSystemEvaluation`, a cross-system evaluation where each dataset is a system (an EEG device with its own montage and sampling rate): for each target system, chosen with ``target_datasets``, pipelines are trained on all other systems. Trials are kept as recorded, with the channels and sampling rate of their system (an object array of trials when systems differ), handling heterogeneous systems being left to the pipelines; pipelines can request the ``systems`` fit metadata through scikit-learn's metadata routing, and target calibration presets are the ones of :class:`moabb.evaluations.CrossSubjectMode`. A new example stacks the systems on the union of their channels, and compares a channel-intersection baseline with a minimum distance to Riemannian masked mean classifier (by `Quentin Barthelemy`_). +- Add the tutorial :ref:`sphx_glr_auto_examples_how_to_benchmark_plot_cross_system_transfer_maskedmean.py`: cross-system transfer with a minimum distance to masked Riemannian mean, on datasets recorded with different montages and sampling rates, pooled in a :class:`~moabb.datasets.compound_dataset.CompoundDataset` and evaluated with :class:`~moabb.evaluations.CrossSubjectEvaluation` (``cv_class=LeaveOneGroupOut, groups="device"``). To support it, the evaluations route any metadata column of the training trials (``device``, ``dataset``, ``subject``...) to the pipeline steps that request it with scikit-learn metadata routing (``set_fit_request(device=True)``), and :class:`~moabb.datasets.fake.FakeDataset` takes a ``device`` (:gh:`1305` by `Quentin Barthelemy`_). - The metadata returned by :meth:`moabb.paradigms.base.BaseParadigm.get_data` has two new columns: ``dataset``, the code of the dataset each trial comes from (for a :class:`~moabb.datasets.compound_dataset.CompoundDataset`, the dataset the subject comes from), and ``device``, the recording device. ``device`` is the canonical name of ``METADATA.acquisition.hardware`` given by :func:`moabb.datasets.metadata.devices.canonical_device` (the 107 free-text descriptions map to 30 devices, e.g. every spelling of the g.tec g.USBamp; ``"unknown"`` when the dataset gives none); :class:`moabb.datasets.Peters2026` returns the DSI-24 or DSI-Flex of each visit. ``CrossSubjectEvaluation(cv_class=LeaveOneGroupOut, groups="dataset")`` holds out one dataset at a time and ``groups="device"`` one device at a time, and :class:`~moabb.evaluations.CrossDatasetSplitter` works on paradigm output (:gh:`1328` by `Bruno Aristimunha`_). - :meth:`moabb.paradigms.base.BaseParadigm.get_data` no longer fails when the recordings it stacks have different channels or numbers of samples, e.g. datasets with different montages or sampling rates pooled in a :class:`~moabb.datasets.compound_dataset.CompoundDataset`: it returns a nested object array of shape ``(n_trials,)`` holding each trial's ``(n_channels, n_times)`` array, and the evaluations accept it, so leave-one-dataset-out runs with :class:`~moabb.evaluations.CrossSubjectEvaluation` and ``cv_class=LeaveOneGroupOut``. Recordings of one shape still give a regular array; ``return_epochs=True`` raises a clear error for such recordings, as ``mne.Epochs`` cannot hold them (:gh:`1327` by `Bruno Aristimunha`_). - Versioned documentation: release tags (``v*``) deploy the "stable" docs to https://moabb.neurotechx.com/docs/ and ``develop`` deploys the "dev" docs to https://moabb.neurotechx.com/dev/, so the default documentation matches the PyPI release; a version switcher in the navbar (``_static/versions.json``) moves between them (:gh:`445`, :gh:`983` by `Bruno Aristimunha`_). @@ -95,6 +95,7 @@ Requirements Bugs ~~~~ +- Fix the evaluations returning each result several times when datasets share a processing pipeline (e.g. :class:`~moabb.datasets.fake.FakeDataset` instances that differ only in their channels): ``process`` collected the whole results group of the pipeline for every dataset, instead of the rows of that dataset (:gh:`1305` by `Quentin Barthelemy`_). - Fix the evaluations training and scoring filter-bank pipelines on the first band only: since 1.4.0 (:gh:`774`), :class:`moabb.evaluations.base.BaseEvaluation` passed ``get_data`` the processing pipeline of the first filter alone, so with :class:`moabb.paradigms.FilterBankLeftRightImagery`, :class:`moabb.paradigms.FilterBankMotorImagery` or :class:`moabb.paradigms.FilterBankSSVEP` the data reaching an array pipeline had one band instead of one per filter (shape ``(n_trials, n_channels, n_times, 1)``). FBCSP and the SSVEP MDM, TS + LR and TS + SVM pipelines scored well below their published values (Nakanishi2015 SSVEP TS + LR 0.26 instead of 0.87; with every band, 0.89). Epoch-based pipelines (TRCA, TDCA, CCA, deep learning) were not affected. Results of filter-bank pipelines computed with MOABB 1.4.0 to 1.7.2 must be recomputed (:gh:`1291` by `Bruno Aristimunha`_). - :class:`moabb.datasets.BNCI2015_004` names its five mental tasks again (word association, mental subtraction, spatial navigation, right hand and feet imagery, coded 1-5 as in the .mat files) and epochs the imagery period [3, 10] s; since v1.5 the events were those of another dataset, the window covered the fixation cross, and :class:`moabb.paradigms.MotorImagery` rejected the dataset (:gh:`1323` by `Bruno Aristimunha`_) - :class:`moabb.datasets.Wang2016` lists all 35 subjects (S35.mat was never loaded) (:gh:`1322` by `Rızgar Ozan`_, merged in :gh:`1323`) diff --git a/examples/how_to_benchmark/plot_cross_system_transfer_maskedmean.py b/examples/how_to_benchmark/plot_cross_system_transfer_maskedmean.py index 6b486e902..35b13c932 100644 --- a/examples/how_to_benchmark/plot_cross_system_transfer_maskedmean.py +++ b/examples/how_to_benchmark/plot_cross_system_transfer_maskedmean.py @@ -9,6 +9,14 @@ not necessarily share the same channels, neither with one another nor with the target system, and they usually record different subjects. +MOABB needs no dedicated evaluation for this. The datasets are pooled in a +:class:`~moabb.datasets.compound_dataset.CompoundDataset`, the paradigm keeps +each trial with the channels and sampling rate of its device, and +:class:`~moabb.evaluations.CrossSubjectEvaluation` holds out one device at a +time with ``cv_class=LeaveOneGroupOut, groups="device"``. The ``device`` +column of the metadata is the recording device, not the dataset: one dataset +can use several devices, and several datasets can share one. + This tutorial introduces the Minimum Distance to Mean (MDM) [1]_ classifier used with the Riemannian masked mean [2]_, which averages covariance matrices that are only partially observed. We will: @@ -35,10 +43,12 @@ from pyriemann.geometry.mean import maskedmean_riemann from sklearn import config_context from sklearn.base import BaseEstimator, TransformerMixin +from sklearn.model_selection import LeaveOneGroupOut from sklearn.pipeline import make_pipeline +from moabb.datasets.compound_dataset import CompoundDataset from moabb.datasets.fake import FakeDataset -from moabb.evaluations import CrossSystemEvaluation +from moabb.evaluations import CrossSubjectEvaluation from moabb.paradigms import LeftRightImagery @@ -52,11 +62,13 @@ # Source 1 only, and C4 and P4 by Source 2 only. The three systems also have # different sampling rates. # -# :class:`~moabb.evaluations.CrossSystemEvaluation` treats each dataset as a -# system: the pipeline is trained on all trials of the source systems and -# scored on each subject/session of the target system. It gives every trial as -# recorded, with the channels and the sampling rate of its system: handling -# heterogeneous systems is the job of the pipeline. +# Each simulated dataset is recorded with its own device, given by +# ``device=``; with real datasets the ``device`` column comes from their +# ``METADATA``. When the held-out device is the target, the pipeline is trained +# on all trials of the other devices and scored on each subject/session of the +# target device. Every trial is given as recorded, with the channels and the +# sampling rate of its device: handling heterogeneous systems is the job of the +# pipeline. systems = { "Source 1": {"channels": ("F3", "Fz", "F4", "Cz", "C4", "P3", "Pz"), "sfreq": 128}, @@ -74,6 +86,7 @@ n_subjects=n_subjects, n_sessions=n_sessions, seed=seed, + device=name, ) for (name, system), n_subjects, n_sessions, seed in zip( systems.items(), (3, 2, 3), (2, 3, 2), (42, 43, 44) @@ -151,41 +164,43 @@ def channel_legend(fig, entries): # from the channel union, the mask of a system is 1 on the channels it records # and 0 on the other channels of the union. -montages = {datasets[name].code: system["channels"] for name, system in systems.items()} -system_masks = {code: np.isin(union, channels) for code, channels in montages.items()} -target = datasets["Target"].code +montages = {name: system["channels"] for name, system in systems.items()} +system_masks = {name: np.isin(union, channels) for name, channels in montages.items()} +target = "Target" ############################################################################### # To place each training trial on the channel union, the first step of the -# pipelines needs the montage of its system: it requests the ``systems`` fit -# metadata, routed by MOABB through scikit-learn's metadata routing. At -# prediction, all trials come from the target system. +# pipelines needs the device of each trial: it requests the ``device`` metadata +# column, which MOABB routes to the steps asking for it through +# scikit-learn's metadata routing (``set_fit_request(device=True)``; any +# metadata column can be requested this way). At prediction, all trials come +# from the held-out device: the one missing from the training devices. class ChannelUnion(TransformerMixin, BaseEstimator): """Zero-fill the trials of each system on the channel union. - ``montages`` maps each system to its channel names, and ``target`` is the - target system, whose montage is applied at prediction. + ``montages`` maps each device to its channel names. The montage applied at + prediction is the one of the held-out device. """ - def __init__(self, montages=None, target=None): + def __init__(self, montages=None): self.montages = montages - self.target = target - def fit(self, X, y=None, *, systems=None): + def fit(self, X, y=None, *, device=None): + if device is None: + raise ValueError("ChannelUnion needs the `device` metadata.") self.channels_ = list( dict.fromkeys(ch for channels in self.montages.values() for ch in channels) ) + self.target_ = held_out(self.montages, device) return self - def fit_transform(self, X, y=None, *, systems=None): - if systems is None: - raise ValueError("ChannelUnion needs `systems` metadata.") - return self.fit(X, y)._stack(X, systems) + def fit_transform(self, X, y=None, *, device=None): + return self.fit(X, y, device=device)._stack(X, device) def transform(self, X): - return self._stack(X, [self.target] * len(X)) + return self._stack(X, [self.target_] * len(X)) def _stack(self, X, systems): X_union = np.empty(len(X), dtype=object) @@ -195,9 +210,15 @@ def _stack(self, X, systems): return X_union +def held_out(devices, training_devices): + """The device missing from the training devices.""" + (device,) = set(devices) - set(training_devices) + return device + + def channel_union(): with config_context(enable_metadata_routing=True): - return ChannelUnion(montages, target).set_fit_request(systems=True) + return ChannelUnion(montages).set_fit_request(device=True) ############################################################################### @@ -254,9 +275,9 @@ def transform(self, X): # # \hat{k} = \arg\min_k d_R \left( M_t^\top C M_t, M_t^\top G_k M_t \right). # -# The classifier receives the masks of the systems and the target system as -# parameters. To know the system of each training trial, it also requests the -# ``systems`` fit metadata. +# The classifier receives the masks of the devices as a parameter. To know the +# device of each training trial, and so the held-out device, it also requests +# the ``device`` metadata. # # The covariance matrices are estimated with the sample covariance matrix: # a masked block of the covariance of zero-filled data is then exactly the @@ -267,25 +288,25 @@ def transform(self, X): class MDMM(MDM): """Minimum distance to masked Riemannian mean. - ``masks`` maps each system to its boolean mask on the channel union, and - ``target`` is the target system, whose mask is applied at prediction. + ``masks`` maps each device to its boolean mask on the channel union. The + mask applied at prediction is the one of the held-out device. """ - def __init__(self, masks=None, target=None, metric="riemann", n_jobs=1): + def __init__(self, masks=None, metric="riemann", n_jobs=1): super().__init__(metric=metric, n_jobs=n_jobs) self.masks = masks - self.target = target def _mask(self, system): """Semi-orthogonal mask matrix of a system.""" return np.eye(len(self.masks[system]))[:, self.masks[system]] - def fit(self, X, y, *, systems=None): - if systems is None: - raise ValueError("MDMM needs `systems` metadata.") + def fit(self, X, y, *, device=None): + if device is None: + raise ValueError("MDMM needs the `device` metadata.") X, y = np.asarray(X), np.asarray(y) self.classes_ = np.unique(y) - masks = [self._mask(system) for system in systems] + self.target_ = held_out(self.masks, device) + masks = [self._mask(d) for d in device] self.covmeans_ = np.stack( [ maskedmean_riemann(X[y == c], [m for m, yi in zip(masks, y) if yi == c]) @@ -295,7 +316,7 @@ def fit(self, X, y, *, systems=None): return self def _predict_distances(self, X): - mask = self._mask(self.target) + mask = self._mask(self.target_) maskedX = mask.T @ np.asarray(X) @ mask return np.stack( [distance_riemann(maskedX, mask.T @ G @ mask) for G in self.covmeans_], axis=1 @@ -336,7 +357,7 @@ def transform(self, X): common = np.flatnonzero(np.all(list(system_masks.values()), axis=0)) with config_context(enable_metadata_routing=True): - mdmm = MDMM(masks=system_masks, target=target).set_fit_request(systems=True) + mdmm = MDMM(masks=system_masks).set_fit_request(device=True) pipelines = { "Zero-filled MDM": make_pipeline(channel_union(), TrialCovariances("oas"), MDM()), @@ -363,9 +384,7 @@ def transform(self, X): # channels recorded by the target system. It never uses zero-filled values. ZERO_FILLED = "#F2CF5B" -source_masks = np.array( - [system_masks[datasets[name].code] for name in ("Source 1", "Source 2")] -) +source_masks = np.array([system_masks[name] for name in ("Source 1", "Source 2")]) target_mask = system_masks[target] @@ -419,27 +438,38 @@ def channel_groups(used, zero_filled): # instances, which keep the tutorial fast and download-free. Their scores have # no scientific meaning; they are used to expose the complete evaluation path. # -# The three systems share the same epoch duration of 3 seconds. -# ``target_datasets`` sets the target system, the other systems being the -# sources; by default, each system is used as target in turn. ``cs_mode`` sets -# the target data available during training: by default, the pipeline sees no -# target trial. +# The three systems share the same epoch duration of 3 seconds. Their subjects +# are pooled in one :class:`~moabb.datasets.compound_dataset.CompoundDataset`, +# and ``groups="device"`` holds out one device at a time: each device is the +# target in turn, the pipeline seeing no trial of it during training. paradigm = LeftRightImagery() +pool = CompoundDataset( + [ + (ds, subject, None, None) + for ds in datasets.values() + for subject in ds.subject_list + ], + code="CompoundDatasetCrossSystem", + interval=datasets["Target"].interval, +) -evaluation = CrossSystemEvaluation( +evaluation = CrossSubjectEvaluation( paradigm=paradigm, - datasets=list(datasets.values()), - target_datasets=[datasets["Target"]], + datasets=[pool], + cv_class=LeaveOneGroupOut, + groups="device", overwrite=True, suffix="cross_system_maskedmean", ) results = evaluation.process(pipelines) +# The device of each pooled subject, to read the results per held-out device. +device_of = {i + 1: ds.device for i, (ds, *_) in enumerate(pool.subjects_list)} +results["device"] = results["subject"].astype(int).map(device_of) print( - results.groupby("pipeline")["score"] + results.groupby(["device", "pipeline"])["score"] .agg(["mean", "std", "count"]) - .sort_values("mean", ascending=False) .round(3) ) @@ -456,7 +486,10 @@ def channel_groups(used, zero_filled): # study would repeat the benchmark over real datasets recorded with different # devices, and use MOABB's statistical-analysis utilities. -paired = results.pivot(index=["subject", "session"], columns="pipeline", values="score") +target_results = results[results["device"] == target] +paired = target_results.pivot( + index=["subject", "session"], columns="pipeline", values="score" +) order = list(pipelines) scores = paired[order].dropna().to_numpy() colors = ["#4C78A8", "#F58518", "#54A24B"] @@ -511,18 +544,14 @@ def channel_groups(used, zero_filled): epochs, y, _ = paradigm.get_data(datasets[name], return_epochs=True) trials += list(epochs.get_data() * datasets[name].unit_factor) labels += list(y) - trial_systems += [datasets[name].code] * len(y) + trial_systems += [name] * len(y) X = np.empty(len(trials), dtype=object) for i, trial in enumerate(trials): X[i] = trial -X_union = ChannelUnion(montages, target).fit_transform(X, systems=trial_systems) +X_union = ChannelUnion(montages).fit_transform(X, device=trial_systems) covs = TrialCovariances("scm").transform(X_union) -centroids = ( - MDMM(masks=system_masks, target=target) - .fit(covs, labels, systems=trial_systems) - .covmeans_ -) +centroids = MDMM(masks=system_masks).fit(covs, labels, device=trial_systems).covmeans_ fig = plt.figure(figsize=(10, 9)) grid = fig.add_gridspec(2, 6, height_ratios=(1.1, 1), hspace=0.45, wspace=0.9) @@ -540,7 +569,7 @@ def channel_groups(used, zero_filled): fig.colorbar(image, ax=fig.axes, location="right", shrink=0.8, label="Covariance") for col, name in enumerate(systems): - mask = system_masks[datasets[name].code] + mask = system_masks[name] ax = fig.add_subplot(grid[0, 2 * col : 2 * col + 2]) ax.imshow( np.eye(len(union))[:, mask], diff --git a/moabb/datasets/fake.py b/moabb/datasets/fake.py index fe11070dd..f3fda277b 100644 --- a/moabb/datasets/fake.py +++ b/moabb/datasets/fake.py @@ -46,6 +46,9 @@ class FakeDataset(BaseDataset): If True, pass events through stim channel. annotations: bool If True, pass events through Annotations. + device: str, default "unknown" + Recording device reported in the ``device`` metadata column, e.g. to + simulate datasets recorded with different devices. """ def __init__( @@ -67,8 +70,10 @@ def __init__( sessions=None, *, return_all_modalities=False, + device="unknown", **kwargs, ): + self.device = device deprecated_renames = { "EventList": "event_list", "NSessions": "n_sessions", @@ -142,6 +147,9 @@ def __init__( # MNE's get_config() checks environment variables first. os.environ[key] = temp_dir + def _device(self, subject, session): + return self.device + def _get_single_subject_data(self, subject): if self.seed is not None: np.random.seed(self.seed + subject) diff --git a/moabb/evaluations/__init__.py b/moabb/evaluations/__init__.py index 274861ba3..0bb2b8c1e 100644 --- a/moabb/evaluations/__init__.py +++ b/moabb/evaluations/__init__.py @@ -7,7 +7,6 @@ from .evaluations import ( CrossSessionEvaluation, CrossSubjectEvaluation, - CrossSystemEvaluation, WithinSessionEvaluation, WithinSubjectEvaluation, ) diff --git a/moabb/evaluations/base.py b/moabb/evaluations/base.py index 226461bba..d62f536ff 100644 --- a/moabb/evaluations/base.py +++ b/moabb/evaluations/base.py @@ -163,7 +163,6 @@ def _evaluate_fold( split_metadata=None, calib_idx=None, score_subjects=None, - extra_fit_params=None, ): """Evaluate a single CV fold. Pure function, no shared mutable state. @@ -195,10 +194,6 @@ def _evaluate_fold( Held-out subjects whose results are still missing for this pipeline. When provided, the fold is fitted once and scored separately for these subjects. - extra_fit_params : dict | None - Evaluation-specific fit metadata (e.g. the systems of the training - trials of a cross-system fold). It is merged into the transfer fit - metadata and routed only to the pipeline steps that request it. Returns ------- @@ -258,7 +253,10 @@ def _evaluate_fold( if calibration_labeled: y_calib = y[calib_idx] if mne_labels else le.transform(y[calib_idx]) - fit_params = {} + # Any metadata column (``device``, ``dataset``, ``subject``...) of the + # training trials can be requested as fit metadata, e.g. + # ``step.set_fit_request(device=True)``; only requested keys are passed. + fit_params = {col: metadata[col].to_numpy()[train_idx] for col in metadata.columns} 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: @@ -267,8 +265,6 @@ def _evaluate_fold( else: fit_params["X_target_labeled"] = X_calib fit_params["y_target_labeled"] = y_calib - if extra_fit_params: - fit_params.update(extra_fit_params) if fit_params: with config_context(enable_metadata_routing=True): requested = get_routing_for_object(cvclf).consumes("fit", set(fit_params)) diff --git a/moabb/evaluations/evaluations.py b/moabb/evaluations/evaluations.py index 01c158c5c..344cab344 100644 --- a/moabb/evaluations/evaluations.py +++ b/moabb/evaluations/evaluations.py @@ -2,7 +2,6 @@ from typing import TYPE_CHECKING, Optional import numpy as np -import pandas as pd from sklearn.base import clone from sklearn.model_selection import ( BaseCrossValidator, @@ -33,7 +32,6 @@ _average_scores, _carbonfootprint, _create_scorer, - _pipeline_requires_epochs, _update_result_with_scores, ) @@ -752,125 +750,6 @@ def _get_incompatibility_reason(self, dataset): return "requirements not met" -class CrossSystemEvaluation(CrossSubjectEvaluation): - """Cross-system evaluation performance. - - A *system* is an acquisition setup: an EEG device with its own montage and - sampling rate, which strongly influences EEG data [1]_. Each dataset is one - system, and systems do not need to share their channels or sampling rate. - For each target system, the pipeline is trained on all other (source) - systems and scored on each subject/session of the target system, a - cross-device transfer learning setting [2]_. - - Each trial is kept as recorded by its system, with its own channels and - sampling rate: handling heterogeneous systems is left to the pipelines. - When all systems share the same channels and number of time samples, ``X`` - is an array of shape ``(n_trials, n_channels, n_times)``; otherwise, it is - an object array of shape ``(n_trials,)`` whose elements are arrays of shape - ``(n_channels_s, n_times_s)`` of their system :math:`s`. - - In addition to the fit metadata of - :class:`~moabb.evaluations.CrossSubjectEvaluation`, pipelines can request, - through scikit-learn's metadata routing (``set_fit_request``), ``systems``: - the system (dataset code) of each training trial. - - Parameters - ---------- - target_datasets : list of :class:`~moabb.datasets.base.BaseDataset` or None - Systems used as targets, each one in turn, among ``datasets``. If None, - every system is used as target (leave-one-system-out). - Defaults to ``None``. - **kwargs - Parameters of :class:`~moabb.evaluations.CrossSubjectEvaluation`. - - Notes - ----- - .. versionadded:: 1.8.0 - - References - ---------- - .. [1] Melnik, A., Legkov, P., Izdebski, K., Kärcher, S. M., Hairston, - W. D., Ferris, D. P., & König, P. (2017). Systems, Subjects, - Sessions: To What Extent Do These Factors Influence EEG Data? - *Frontiers in Human Neuroscience*, 11, 150. - https://doi.org/10.3389/fnhum.2017.00150 - .. [2] Wu, D., Xu, Y., & Lu, B.-L. (2020). Transfer Learning for EEG-Based - Brain-Computer Interfaces: A Review of Progress Made Since 2016. - *IEEE Transactions on Cognitive and Developmental Systems*, 14(1), - 4-19. https://doi.org/10.1109/TCDS.2020.3007453 - """ - - _eval_type = "CrossSystem" - - def __init__(self, *args, target_datasets=None, **kwargs): - super().__init__(*args, **kwargs) - self.groups = "system" - # All systems are loaded as sources, only targets are evaluated. - self.systems = self.datasets - self._trials = None - if target_datasets is not None: - codes = {system.code for system in self.systems} - unknown = [d.code for d in target_datasets if d.code not in codes] - if unknown: - raise ValueError( - f"target_datasets {unknown} are not among the valid datasets " - f"{sorted(codes)}." - ) - self.datasets = list(target_datasets) - - def is_valid(self, dataset): - return True - - def _load_data( - self, dataset, run_pipes, process_pipeline, postprocess_pipeline, subjects=None - ): - """Concatenate the trials of all systems, as recorded.""" - if postprocess_pipeline is not None: - raise ValueError( - "CrossSystemEvaluation does not support postprocess_pipeline." - ) - if any(_pipeline_requires_epochs(clf) for clf in run_pipes.values()): - raise ValueError( - "CrossSystemEvaluation gives pipelines arrays, not mne.Epochs." - ) - # ponytail: every target needs the same trials of all systems, so load - # them once per evaluation (held in memory for the whole run). - if self._trials is None: - self._trials = self._load_systems() - return self._trials - - def _load_systems(self): - systems = [] - for system in self.systems: - epochs, y, metadata = self.paradigm.get_data( - system, return_epochs=True, cache_config=self.cache_config - ) - X = epochs.get_data() * system.unit_factor - systems.append((epochs.ch_names, X, y, metadata.assign(system=system.code))) - - if len({(tuple(ch_names), X.shape[-1]) for ch_names, X, _, _ in systems}) == 1: - X_all = np.concatenate([X for _, X, _, _ in systems]) - else: - trials = [trial for _, X, _, _ in systems for trial in X] - X_all = np.empty(len(trials), dtype=object) - for i, trial in enumerate(trials): - X_all[i] = trial - y_all = np.concatenate([y for _, _, y, _ in systems]) - metadata = pd.concat([metadata for *_, metadata in systems], ignore_index=True) - return X_all, y_all, metadata - - def _build_task_list(self, dataset, X, y, metadata, *args): - """Keep the fold holding out ``dataset`` and route the systems.""" - tasks = super()._build_task_list(dataset, X, y, metadata, *args) - systems = metadata["system"].to_numpy() - tasks = [ - task for task in tasks if np.all(systems[task["test_idx"]] == dataset.code) - ] - for task in tasks: - task["extra_fit_params"] = {"systems": systems[task["train_idx"]]} - return tasks - - class WithinSubjectEvaluation(BaseEvaluation): """Within-subject k-fold cross-validation pooling all sessions. diff --git a/moabb/tests/test_cross_system_evaluation.py b/moabb/tests/test_cross_system_evaluation.py deleted file mode 100644 index 4e13b5755..000000000 --- a/moabb/tests/test_cross_system_evaluation.py +++ /dev/null @@ -1,232 +0,0 @@ -"""Tests for the cross-system CrossSystemEvaluation.""" - -import numpy as np -import pandas as pd -import pytest -import sklearn.base -from pyriemann.estimation import Covariances -from pyriemann.spatialfilters import CSP -from sklearn import config_context -from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA -from sklearn.pipeline import FunctionTransformer, make_pipeline - -from moabb.datasets.fake import FakeDataset -from moabb.evaluations import CrossSubjectMode, CrossSystemEvaluation -from moabb.paradigms import LeftRightImagery - - -_CAPTURE = [] - - -class _SystemRecorder(sklearn.base.TransformerMixin, sklearn.base.BaseEstimator): - """A system-aware step that records the cross-system metadata it is routed.""" - - def fit(self, X, y=None, systems=None, X_target_unlabeled=None): - _CAPTURE.append( - {"X": X, "systems": systems, "X_target_unlabeled": X_target_unlabeled} - ) - return self - - def transform(self, X): - return X - - -def _trial_variances(X): - """One feature per trial, whatever its channels and time samples.""" - return np.array([[np.var(trial)] for trial in X]) - - -@pytest.fixture(autouse=True) -def _clear_capture(): - _CAPTURE.clear() - yield - _CAPTURE.clear() - - -def _datasets(sfreqs=(128, 128), channels=(("C3", "Cz", "C4"), ("Cz", "C4", "Pz", "P4"))): - """Two fake systems, by default with different channels.""" - return [ - FakeDataset( - ["left_hand", "right_hand"], - channels=channels[0], - n_subjects=2, - n_sessions=1, - sfreq=sfreqs[0], - seed=1, - ), - FakeDataset( - ["left_hand", "right_hand"], - channels=channels[1], - n_subjects=3, - n_sessions=2, - sfreq=sfreqs[1], - seed=2, - ), - ] - - -def _recorder_pipeline(): - with config_context(enable_metadata_routing=True): - step = _SystemRecorder().set_fit_request(systems=True, X_target_unlabeled=True) - return make_pipeline(step, FunctionTransformer(_trial_variances), LDA()) - - -def _evaluation(tmp_path, datasets, **kwargs): - kwargs.setdefault("paradigm", LeftRightImagery()) - kwargs.setdefault("overwrite", True) - return CrossSystemEvaluation(datasets=datasets, hdf5_path=tmp_path, **kwargs) - - -def _source_and_target(capture, datasets): - (source,) = np.unique(capture["systems"]) - (target,) = {dataset.code for dataset in datasets} - {source} - return source, target - - -def _identity(X): - return X - - -def test_trials_keep_the_channels_of_their_system(tmp_path): - ds_a, ds_b = _datasets() - evaluation = _evaluation(tmp_path, [ds_a, ds_b]) - X, y, metadata = evaluation._load_data(ds_a, {}, None, None) - - # Systems with different channels give an object array of trials. - assert X.dtype == object and X.shape == (len(y),) - n_channels = {ds_a.code: 3, ds_b.code: 4} - for trial, system in zip(X, metadata["system"]): - assert trial.shape == (n_channels[system], 3 * 128 + 1) - - results = evaluation.process({"rec": _recorder_pipeline()}) - # One row per target subject/session: 2 x 1 for A, 3 x 2 for B. - assert results.groupby("dataset").size().to_dict() == {ds_a.code: 2, ds_b.code: 6} - assert len(_CAPTURE) == 2 - for capture in _CAPTURE: - source, _ = _source_and_target(capture, [ds_a, ds_b]) - # Sources are trained on their own trials, routed with their system. - assert {trial.shape[0] for trial in capture["X"]} == {n_channels[source]} - assert len(capture["systems"]) == len(capture["X"]) - # Default TRAIN mode: no target trial is routed. - assert capture["X_target_unlabeled"] is None - - -def test_trials_keep_the_sampling_rate_of_their_system(tmp_path): - datasets = _datasets(sfreqs=(128, 256), channels=(("C3", "Cz", "C4"),) * 2) - evaluation = _evaluation(tmp_path, datasets) - X, _, metadata = evaluation._load_data(datasets[0], {}, None, None) - - assert X.dtype == object - n_times = {datasets[0].code: 3 * 128 + 1, datasets[1].code: 3 * 256 + 1} - for trial, system in zip(X, metadata["system"]): - assert trial.shape == (3, n_times[system]) - assert len(evaluation.process({"rec": _recorder_pipeline()})) == 8 - - -def test_same_channels_and_sampling_rate_give_regular_array(tmp_path): - datasets = _datasets(channels=(("C3", "Cz", "C4"),) * 2) - evaluation = _evaluation(tmp_path, datasets) - X, y, _ = evaluation._load_data(datasets[0], {}, None, None) - assert X.dtype != object and X.shape == (len(y), 3, 3 * 128 + 1) - - # A standard pipeline then runs as is. - results = evaluation.process( - {"csp": make_pipeline(Covariances("oas"), CSP(2), LDA())} - ) - assert len(results) == 8 - assert results["score"].between(0, 1).all() - - -def test_calibration_routes_target_slice(tmp_path): - ds_a, ds_b = _datasets() - results = _evaluation( - tmp_path, [ds_a, ds_b], cs_mode=CrossSubjectMode.TRAIN_AND_TARGET_UNLABELED_20P - ).process({"rec": _recorder_pipeline()}) - - (capture,) = [ - c for c in _CAPTURE if _source_and_target(c, [ds_a, ds_b])[1] == ds_a.code - ] - X_target = capture["X_target_unlabeled"] - assert X_target is not None and len(X_target) > 0 - # Target trials keep the channels of the target system. - assert {trial.shape[0] for trial in X_target} == {3} - # Calibration and scored target trials partition the target system. - _, y_a, _ = LeftRightImagery().get_data(ds_a) - scored = results.loc[results["dataset"] == ds_a.code, "samples_test"].sum() - assert scored + len(X_target) == len(y_a) - - -def test_target_datasets_restricts_targets(tmp_path): - ds_a, ds_b = _datasets() - ds_c = FakeDataset( - ["left_hand", "right_hand"], - channels=("C3", "Pz"), - n_subjects=2, - n_sessions=1, - seed=3, - ) - results = _evaluation(tmp_path, [ds_a, ds_b, ds_c], target_datasets=[ds_c]).process( - {"rec": _recorder_pipeline()} - ) - - assert set(results["dataset"]) == {ds_c.code} - assert len(results) == 2 - (capture,) = _CAPTURE - # The non-target systems are the sources. - assert set(capture["systems"]) == {ds_a.code, ds_b.code} - - -def test_trialwise_mode(tmp_path): - evaluation = _evaluation( - tmp_path, _datasets(), cs_mode=CrossSubjectMode.TRAIN_TRIALWISE - ) - assert len(evaluation.process({"rec": _recorder_pipeline()})) == 8 - - -def test_results_are_cached(tmp_path): - datasets = _datasets() - pipelines = {"rec": _recorder_pipeline()} - first = _evaluation(tmp_path, datasets).process(pipelines) - n_fits = len(_CAPTURE) - - cached = _evaluation(tmp_path, datasets, overwrite=False).process(pipelines) - - assert len(_CAPTURE) == n_fits - keys = ["dataset", "subject", "session"] - pd.testing.assert_frame_equal( - first.sort_values(keys).reset_index(drop=True), - cached.sort_values(keys).reset_index(drop=True), - check_like=True, - ) - - -def test_evaluate_single_target(tmp_path): - ds_a, ds_b = _datasets() - evaluation = _evaluation(tmp_path, [ds_a, ds_b]) - process_pipeline = evaluation.paradigm.make_process_pipelines(ds_b)[0] - results = list( - evaluation.evaluate( - ds_b, - {"rec": _recorder_pipeline()}, - param_grid=None, - process_pipeline=process_pipeline, - ) - ) - assert len(results) == 6 - assert {res["dataset"].code for res in results} == {ds_b.code} - - -def test_rejects_postprocess_pipeline(tmp_path): - evaluation = _evaluation(tmp_path, _datasets()) - with pytest.raises(ValueError, match="postprocess_pipeline"): - evaluation.process( - {"rec": _recorder_pipeline()}, - postprocess_pipeline=FunctionTransformer(_identity), - ) - - -def test_unknown_target_dataset_raises(tmp_path): - ds_a, ds_b = _datasets() - other = FakeDataset(["left_hand", "right_hand"], channels=("C3", "Pz"), n_subjects=2) - with pytest.raises(ValueError, match="not among the valid datasets"): - _evaluation(tmp_path, [ds_a, ds_b], target_datasets=[other]) diff --git a/moabb/tests/test_evaluations.py b/moabb/tests/test_evaluations.py index 32adbc88a..f5ad5d7a0 100644 --- a/moabb/tests/test_evaluations.py +++ b/moabb/tests/test_evaluations.py @@ -45,11 +45,6 @@ def _identity(x): return x -def _trial_variances(X): - """One feature per trial, whatever its channels and time samples.""" - return np.array([[np.var(trial)] for trial in X]) - - class OpaqueLeaveOneGroupOut(LeaveOneGroupOut): """Real group CV whose constructor signature cannot be inspected.""" @@ -1605,7 +1600,7 @@ def test_get_nchan(self): from moabb.evaluations.utils import _get_nchan assert _get_nchan(np.zeros((4, 3, 10))) == 3 - # Trials of different lengths, as given by CrossSystemEvaluation. + # Trials of different lengths, as given for recordings of different rates. trials = np.empty(2, dtype=object) trials[0], trials[1] = np.zeros((3, 10)), np.zeros((3, 20)) assert _get_nchan(trials) == 3 @@ -2012,29 +2007,6 @@ def test_cross_subject_process_structure(self, tmp_path): # 2 subjects × 2 sessions = 4 results assert len(results) == 4 - def test_cross_system_process_structure(self, tmp_path): - """CrossSystem process() returns one result per target subject/session.""" - systems = [ - FakeDataset( - ["left_hand", "right_hand"], channels=channels, n_subjects=2, seed=seed - ) - for channels, seed in ((("C3", "Cz", "C4"), 1), (("Cz", "C4", "Pz"), 2)) - ] - evaluation = ev.CrossSystemEvaluation( - paradigm=FakeImageryParadigm(), - datasets=systems, - overwrite=True, - hdf5_path=str(tmp_path / "parallel_test"), - ) - # Trials keep the channels of their system, so the pipeline must not - # depend on them. - results = evaluation.process( - {"V": make_pipeline(FunctionTransformer(_trial_variances), LDA())} - ) - # 2 systems as targets x 2 subjects x 2 sessions = 8 results - assert len(results) == 8 - assert set(results["dataset"]) == {system.code for system in systems} - def test_process_keeps_rows_of_each_dataset_once(self, tmp_path): """Datasets sharing a pipeline digest do not duplicate their results.""" # Same parameters except the channels: same process pipeline digest. @@ -3068,17 +3040,17 @@ def test_cross_subject_cv_class_with_reserved_metadata_columns_saves(tmp_path): # --------------------------------------------------------------------------- -# Evaluation-specific fit metadata routed by _evaluate_fold. +# Metadata columns routed as fit metadata by _evaluate_fold. # --------------------------------------------------------------------------- _EXTRA_CAPTURE = [] class _ExtraRecorder(sklearn.base.TransformerMixin, sklearn.base.BaseEstimator): - """Records the evaluation-specific metadata it is routed.""" + """Records the ``device`` metadata it is routed.""" - def fit(self, X, y=None, systems=None): - _EXTRA_CAPTURE.append(systems) + def fit(self, X, y=None, device=None): + _EXTRA_CAPTURE.append(device) return self def transform(self, X): @@ -3086,16 +3058,16 @@ def transform(self, X): @pytest.mark.parametrize("requested", [True, False]) -def test_evaluate_fold_routes_extra_fit_params(requested): - """extra_fit_params reach only the steps requesting them, without a - transfer split, and are ignored by the other pipelines.""" +def test_evaluate_fold_routes_requested_metadata_columns(requested): + """A metadata column reaches only the steps requesting it, without a + transfer split, and is ignored by the other pipelines.""" from sklearn import config_context from moabb.evaluations.base import _evaluate_fold _EXTRA_CAPTURE.clear() with config_context(enable_metadata_routing=True): - step = _ExtraRecorder().set_fit_request(systems=requested) + step = _ExtraRecorder().set_fit_request(device=requested) pipe = make_pipeline(step, Covariances("oas"), CSP(2), LDA()) evaluation = ev.CrossSessionEvaluation( @@ -3104,7 +3076,7 @@ def test_evaluate_fold_routes_extra_fit_params(requested): X, y, metadata = FakeImageryParadigm().get_data(dataset, subjects=[1]) train = np.flatnonzero(metadata["session"] == metadata["session"].iloc[0]) test = np.flatnonzero(metadata["session"] != metadata["session"].iloc[0]) - systems = np.array(["system"] * len(y)) + metadata["device"] = np.where(np.arange(len(y)) % 2, "Amp X", "Amp Y") results = _evaluate_fold( X, @@ -3119,13 +3091,12 @@ def test_evaluate_fold_routes_extra_fit_params(requested): 1, metadata["session"].iloc[-1], 0, - extra_fit_params={"systems": systems[train]}, ) assert len(results) == 1 (routed,) = _EXTRA_CAPTURE if requested: - np.testing.assert_array_equal(routed, systems[train]) + np.testing.assert_array_equal(routed, metadata["device"].to_numpy()[train]) else: assert routed is None @@ -3294,3 +3265,40 @@ def test_cross_subject_leave_one_device_out(tmp_path): per_subject = len(y) / 6 # Amp X subjects (1-4) train on Amp Y's two subjects, and the reverse. assert list(results["samples"]) == [2 * per_subject] * 4 + [4 * per_subject] * 2 + + +def test_leave_one_device_out_routes_the_training_devices(tmp_path): + from sklearn import config_context + + from moabb.datasets.compound_dataset import CompoundDataset + + a, b, c = ( + FakeDataset( + ["left_hand", "right_hand"], n_subjects=1, n_sessions=1, seed=s, device=d + ) + for s, d in ((1, "Amp X"), (2, "Amp X"), (3, "Amp Y")) + ) + b.code, c.code = "FakeDatasetB", "FakeDatasetC" + pool = CompoundDataset( + [(d, 1, None, None) for d in (a, b, c)], + code="CompoundDatasetABC", + interval=a.interval, + ) + _EXTRA_CAPTURE.clear() + with config_context(enable_metadata_routing=True): + step = _ExtraRecorder().set_fit_request(device=True) + evaluation = ev.CrossSubjectEvaluation( + paradigm=FakeImageryParadigm(), + datasets=[pool], + cv_class=LeaveOneGroupOut, + groups="device", + hdf5_path=str(tmp_path), + overwrite=True, + ) + results = evaluation.process( + {"P": make_pipeline(step, Covariances("oas"), CSP(2), LDA())} + ) + assert len(results) == 3 + # One fold per device; each model only saw the other device. + seen = sorted(tuple(sorted(set(devices))) for devices in _EXTRA_CAPTURE) + assert seen == [("Amp X",), ("Amp Y",)]