ENH: Viz for spatial filters (#13332)

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Eric Larson <[email protected]>

Author: 3 people

doc/_includes/ged.rst

@@ -14,7 +14,7 @@ This section describes the mathematical formulation and application of
 Generalized Eigendecomposition (GED), often used in spatial filtering
 and source separation algorithms, such as :class:`mne.decoding.CSP`, 
 :class:`mne.decoding.SPoC`, :class:`mne.decoding.SSD` and 
-:class:`mne.preprocessing.Xdawn`.
+:class:`mne.decoding.XdawnTransformer`.
 
 The core principle of GED is to find a set of channel weights (spatial filter) 
 that maximizes the ratio of signal power between two data features. 

doc/api/decoding.rst

@@ -30,6 +30,7 @@ Decoding
    SPoC
    SSD
    XdawnTransformer
+   SpatialFilter
 
 Functions that assist with decoding and model fitting:
 
@@ -39,3 +40,4 @@ Functions that assist with decoding and model fitting:
    compute_ems
    cross_val_multiscore
    get_coef
+   get_spatial_filter_from_estimator

doc/changes/dev/13332.newfeature.rst

@@ -0,0 +1,4 @@
+Implement :class:`mne.decoding.SpatialFilter` class returned by :func:`mne.decoding.get_spatial_filter_from_estimator` for
+visualisation of filters and patterns for :class:`mne.decoding.LinearModel` 
+and additionally eigenvalues for GED-based transformers such as 
+:class:`mne.decoding.XdawnTransformer`, :class:`mne.decoding.CSP`, by `Gennadiy Belonosov`_.

examples/decoding/decoding_csp_eeg.py

@@ -14,6 +14,7 @@
 `PhysioNet documentation page <https://physionet.org/content/eegmmidb/1.0.0/>`_.
 The dataset is available at PhysioNet :footcite:`GoldbergerEtAl2000`.
 """
+
 # Authors: Martin Billinger <[email protected]>
 #
 # License: BSD-3-Clause
@@ -30,7 +31,7 @@
 from mne import Epochs, pick_types
 from mne.channels import make_standard_montage
 from mne.datasets import eegbci
-from mne.decoding import CSP
+from mne.decoding import CSP, get_spatial_filter_from_estimator
 from mne.io import concatenate_raws, read_raw_edf
 
 print(__doc__)
@@ -95,10 +96,11 @@
 class_balance = max(class_balance, 1.0 - class_balance)
 print(f"Classification accuracy: {np.mean(scores)} / Chance level: {class_balance}")
 
-# plot CSP patterns estimated on full data for visualization
+# plot eigenvalues and patterns estimated on full data for visualization
 csp.fit_transform(epochs_data, labels)
-
-csp.plot_patterns(epochs.info, ch_type="eeg", units="Patterns (AU)", size=1.5)
+spf = get_spatial_filter_from_estimator(csp, info=epochs.info)
+spf.plot_scree()
+spf.plot_patterns(components=np.arange(4))
 
 # %%
 # Look at performance over time

examples/decoding/decoding_spoc_CMC.py

@@ -32,7 +32,7 @@
 import mne
 from mne import Epochs
 from mne.datasets.fieldtrip_cmc import data_path
-from mne.decoding import SPoC
+from mne.decoding import SPoC, get_spatial_filter_from_estimator
 
 # Define parameters
 fname = data_path() / "SubjectCMC.ds"
@@ -82,9 +82,18 @@
 # Plot the contributions to the detected components (i.e., the forward model)
 
 spoc.fit(X, y)
-spoc.plot_patterns(meg_epochs.info)
+spf = get_spatial_filter_from_estimator(spoc, info=meg_epochs.info)
+spf.plot_scree()
+
+# Plot patterns for the first three components
+# with largest absolute generalized eigenvalues,
+# as we can see on the scree plot
+spf.plot_patterns(components=[0, 1, 2])
+
 
 ##############################################################################
 # References
 # ----------
 # .. footbibliography::
+
+# %%

examples/decoding/decoding_xdawn_eeg.py

@@ -10,6 +10,7 @@
 Channels are concatenated and rescaled to create features vectors that will be
 fed into a logistic regression.
 """
+
 # Authors: Alexandre Barachant <[email protected]>
 #
 # License: BSD-3-Clause
@@ -26,10 +27,9 @@
 from sklearn.pipeline import make_pipeline
 from sklearn.preprocessing import MinMaxScaler
 
-from mne import Epochs, EvokedArray, create_info, io, pick_types, read_events
+from mne import Epochs, io, pick_types, read_events
 from mne.datasets import sample
-from mne.decoding import Vectorizer
-from mne.preprocessing import Xdawn
+from mne.decoding import Vectorizer, XdawnTransformer, get_spatial_filter_from_estimator
 
 print(__doc__)
 
@@ -71,31 +71,33 @@
 
 # Create classification pipeline
 clf = make_pipeline(
-    Xdawn(n_components=n_filter),
+    XdawnTransformer(n_components=n_filter),
     Vectorizer(),
     MinMaxScaler(),
     OneVsRestClassifier(LogisticRegression(penalty="l1", solver="liblinear")),
 )
 
-# Get the labels
-labels = epochs.events[:, -1]
+# Get the data and labels
+# X is of shape (n_epochs, n_channels, n_times)
+X = epochs.get_data(copy=False)
+y = epochs.events[:, -1]
 
 # Cross validator
 cv = StratifiedKFold(n_splits=10, shuffle=True, random_state=42)
 
 # Do cross-validation
-preds = np.empty(len(labels))
-for train, test in cv.split(epochs, labels):
-    clf.fit(epochs[train], labels[train])
-    preds[test] = clf.predict(epochs[test])
+preds = np.empty(len(y))
+for train, test in cv.split(epochs, y):
+    clf.fit(X[train], y[train])
+    preds[test] = clf.predict(X[test])
 
 # Classification report
 target_names = ["aud_l", "aud_r", "vis_l", "vis_r"]
-report = classification_report(labels, preds, target_names=target_names)
+report = classification_report(y, preds, target_names=target_names)
 print(report)
 
 # Normalized confusion matrix
-cm = confusion_matrix(labels, preds)
+cm = confusion_matrix(y, preds)
 cm_normalized = cm.astype(float) / cm.sum(axis=1)[:, np.newaxis]
 
 # Plot confusion matrix
@@ -109,30 +111,35 @@
 ax.set(ylabel="True label", xlabel="Predicted label")
 
 # %%
-# The ``patterns_`` attribute of a fitted Xdawn instance (here from the last
-# cross-validation fold) can be used for visualization.
-
-fig, axes = plt.subplots(
-    nrows=len(event_id),
-    ncols=n_filter,
-    figsize=(n_filter, len(event_id) * 2),
-    layout="constrained",
+# Patterns of a fitted XdawnTransformer instance (here from the last
+# cross-validation fold) can be visualized using SpatialFilter container.
+
+# Instantiate SpatialFilter
+spf = get_spatial_filter_from_estimator(
+    clf, info=epochs.info, step_name="xdawntransformer"
+)
+
+# Let's first examine the scree plot of generalized eigenvalues
+# for each class.
+spf.plot_scree(title="")
+
+# We can see that for all four classes ~five largest components
+# capture most of the variance, let's plot their patterns.
+# Each class will now return its own figure
+components_to_plot = np.arange(5)
+figs = spf.plot_patterns(
+    # Indices of patterns to plot,
+    # we will plot the first three for each class
+    components=components_to_plot,
+    show=False,  # to set the titles below
 )
-fitted_xdawn = clf.steps[0][1]
-info = create_info(epochs.ch_names, 1, epochs.get_channel_types())
-info.set_montage(epochs.get_montage())
-for ii, cur_class in enumerate(sorted(event_id)):
-    cur_patterns = fitted_xdawn.patterns_[cur_class]
-    pattern_evoked = EvokedArray(cur_patterns[:n_filter].T, info, tmin=0)
-    pattern_evoked.plot_topomap(
-        times=np.arange(n_filter),
-        time_format="Component %d" if ii == 0 else "",
-        colorbar=False,
-        show_names=False,
-        axes=axes[ii],
-        show=False,
-    )
-    axes[ii, 0].set(ylabel=cur_class)
+
+# Set the class titles
+event_id_reversed = {v: k for k, v in event_id.items()}
+for fig, class_idx in zip(figs, clf[0].classes_):
+    class_name = event_id_reversed[class_idx]
+    fig.suptitle(class_name, fontsize=16)
+
 
 # %%
 # References

examples/decoding/linear_model_patterns.py

@@ -14,6 +14,7 @@
 Note patterns/filters in MEG data are more similar than EEG data
 because the noise is less spatially correlated in MEG than EEG.
 """
+
 # Authors: Alexandre Gramfort <[email protected]>
 #          Romain Trachel <[email protected]>
 #          Jean-Rémi King <[email protected]>
@@ -28,11 +29,16 @@
 from sklearn.preprocessing import StandardScaler
 
 import mne
-from mne import EvokedArray, io
+from mne import io
 from mne.datasets import sample
 
 # import a linear classifier from mne.decoding
-from mne.decoding import LinearModel, Vectorizer, get_coef
+from mne.decoding import (
+    LinearModel,
+    SpatialFilter,
+    Vectorizer,
+    get_spatial_filter_from_estimator,
+)
 
 print(__doc__)
 
@@ -77,20 +83,38 @@
 X = scaler.fit_transform(meg_data)
 model.fit(X, labels)
 
-# Extract and plot spatial filters and spatial patterns
+coefs = dict()
 for name, coef in (("patterns", model.patterns_), ("filters", model.filters_)):
     # We fit the linear model on Z-scored data. To make the filters
     # interpretable, we must reverse this normalization step
     coef = scaler.inverse_transform([coef])[0]
 
     # The data was vectorized to fit a single model across all time points and
     # all channels. We thus reshape it:
-    coef = coef.reshape(len(meg_epochs.ch_names), -1)
-
-    # Plot
-    evoked = EvokedArray(coef, meg_epochs.info, tmin=epochs.tmin)
-    fig = evoked.plot_topomap()
-    fig.suptitle(f"MEG {name}")
+    coefs[name] = coef.reshape(len(meg_epochs.ch_names), -1).T
+
+# Now we can instantiate the visualization container
+spf = SpatialFilter(info=meg_epochs.info, **coefs)
+fig = spf.plot_patterns(
+    # we will automatically select patterns
+    components="auto",
+    # as our filters and patterns correspond to actual times
+    # we can align them
+    tmin=epochs.tmin,
+    units="fT",  # it's physical - we inversed the scaling
+    show=False,  # to set the title below
+    name_format=None,  # to plot actual times
+)
+fig.suptitle("MEG patterns")
+# Same for filters
+fig = spf.plot_filters(
+    components="auto",
+    tmin=epochs.tmin,
+    units="fT",
+    show=False,
+    name_format=None,
+)
+fig.suptitle("MEG filters")
 
 # %%
 # Let's do the same on EEG data using a scikit-learn pipeline
@@ -107,15 +131,26 @@
     ),
 )
 clf.fit(X, y)
-
-# Extract and plot patterns and filters
-for name in ("patterns_", "filters_"):
-    # The `inverse_transform` parameter will call this method on any estimator
-    # contained in the pipeline, in reverse order.
-    coef = get_coef(clf, name, inverse_transform=True)
-    evoked = EvokedArray(coef, epochs.info, tmin=epochs.tmin)
-    fig = evoked.plot_topomap()
-    fig.suptitle(f"EEG {name[:-1]}")
+spf = get_spatial_filter_from_estimator(
+    clf, info=epochs.info, inverse_transform=True, step_name="linearmodel"
+)
+fig = spf.plot_patterns(
+    components="auto",
+    tmin=epochs.tmin,
+    units="uV",
+    show=False,
+    name_format=None,
+)
+fig.suptitle("EEG patterns")
+# Same for filters
+fig = spf.plot_filters(
+    components="auto",
+    tmin=epochs.tmin,
+    units="uV",
+    show=False,
+    name_format=None,
+)
+fig.suptitle("EEG filters")
 
 # %%
 # References

examples/decoding/ssd_spatial_filters.py

@@ -26,7 +26,7 @@
 import mne
 from mne import Epochs
 from mne.datasets.fieldtrip_cmc import data_path
-from mne.decoding import SSD
+from mne.decoding import SSD, get_spatial_filter_from_estimator
 
 # %%
 # Define parameters
@@ -70,8 +70,8 @@
 # (W^{-1}) or by multiplying the noise cov with the filters Eq. (22) (C_n W)^t.
 # We rely on the inversion approach here.
 
-pattern = mne.EvokedArray(data=ssd.patterns_[:4].T, info=ssd.info)
-pattern.plot_topomap(units=dict(mag="A.U."), time_format="")
+spf = get_spatial_filter_from_estimator(ssd, info=ssd.info)
+spf.plot_patterns(components=list(range(4)))
 
 # The topographies suggest that we picked up a parietal alpha generator.
 
@@ -150,8 +150,8 @@
 ssd_epochs.fit(X=epochs.get_data(copy=False))
 
 # Plot topographies.
-pattern_epochs = mne.EvokedArray(data=ssd_epochs.patterns_[:4].T, info=ssd_epochs.info)
-pattern_epochs.plot_topomap(units=dict(mag="A.U."), time_format="")
+spf = get_spatial_filter_from_estimator(ssd_epochs, info=ssd_epochs.info)
+spf.plot_patterns(components=list(range(4)))
 # %%
 # References
 # ----------

mne/decoding/__init__.pyi

@@ -11,6 +11,7 @@ __all__ = [
     "SSD",
     "Scaler",
     "SlidingEstimator",
+    "SpatialFilter",
     "TemporalFilter",
     "TimeDelayingRidge",
     "TimeFrequency",
@@ -21,6 +22,7 @@ __all__ = [
     "compute_ems",
     "cross_val_multiscore",
     "get_coef",
+    "get_spatial_filter_from_estimator",
 ]
 from .base import (
     BaseEstimator,
@@ -33,6 +35,7 @@ from .csp import CSP, SPoC
 from .ems import EMS, compute_ems
 from .receptive_field import ReceptiveField
 from .search_light import GeneralizingEstimator, SlidingEstimator
+from .spatial_filter import SpatialFilter, get_spatial_filter_from_estimator
 from .ssd import SSD
 from .time_delaying_ridge import TimeDelayingRidge
 from .time_frequency import TimeFrequency