# sphinx_gallery_thumbnail_number = 6

import numpy as np
import matplotlib.pyplot as plt

from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression

import mne
from mne.datasets import sample
from mne.decoding import (SlidingEstimator, GeneralizingEstimator, Scaler,
                          cross_val_multiscore, LinearModel, get_coef,
                          Vectorizer, CSP)

data_path = sample.data_path()

raw_fname = data_path + '/MEG/sample/sample_audvis_raw.fif'
tmin, tmax = -0.200, 0.500
event_id = {'Auditory/Left': 1, 'Visual/Left': 3}  # just use two
raw = mne.io.read_raw_fif(raw_fname, preload=True)

# The subsequent decoding analyses only capture evoked responses, so we can
# low-pass the MEG data. Usually a value more like 40 Hz would be used,
# but here low-pass at 20 so we can more heavily decimate, and allow
# the examlpe to run faster. The 2 Hz high-pass helps improve CSP.
raw.filter(2, 20)
events = mne.find_events(raw, 'STI 014')

# Set up pick list: EEG + MEG - bad channels (modify to your needs)
raw.info['bads'] += ['MEG 2443', 'EEG 053']  # bads + 2 more
picks = mne.pick_types(raw.info, meg='grad', eeg=False, stim=True, eog=True,
                       exclude='bads')

# Read epochs
epochs = mne.Epochs(raw, events, event_id, tmin, tmax, proj=True,
                    picks=picks, baseline=(None, 0.), preload=True,
                    reject=dict(grad=4000e-13, eog=150e-6), decim=10)
epochs.pick_types(meg=True, exclude='bads')  # remove stim and EOG

X = epochs.get_data()  # MEG signals: n_epochs, n_meg_channels, n_times
y = epochs.events[:, 2]  # target: Audio left or right

"""
Explanation: ===============
Decoding (MVPA)
===============
   :depth: 3
.. include:: ../git_links.inc
Design philosophy
Decoding (a.k.a. MVPA) in MNE largely follows the machine
learning API of the scikit-learn package.
Each estimator implements fit, transform, fit_transform, and
(optionally) inverse_transform methods. For more details on this design,
visit scikit-learn_. For additional theoretical insights into the decoding
framework in MNE, see [1]_.
For ease of comprehension, we will denote instantiations of the class using
the same name as the class but in small caps instead of camel cases.
Let's start by loading data for a simple two-class problem:
End of explanation
"""


# Uses all MEG sensors and time points as separate classification
# features, so the resulting filters used are spatio-temporal
clf = make_pipeline(Scaler(epochs.info),
                    Vectorizer(),
                    LogisticRegression(solver='lbfgs'))

scores = cross_val_multiscore(clf, X, y, cv=5, n_jobs=1)

# Mean scores across cross-validation splits
score = np.mean(scores, axis=0)
print('Spatio-temporal: %0.1f%%' % (100 * score,))

"""
Explanation: Transformation classes
Scaler
^^^^^^
The :class:mne.decoding.Scaler will standardize the data based on channel
scales. In the simplest modes scalings=None or scalings=dict(...),
each data channel type (e.g., mag, grad, eeg) is treated separately and
scaled by a constant. This is the approach used by e.g.,
:func:mne.compute_covariance to standardize channel scales.
If scalings='mean' or scalings='median', each channel is scaled using
empirical measures. Each channel is scaled independently by the mean and
standand deviation, or median and interquartile range, respectively, across
all epochs and time points during :class:~mne.decoding.Scaler.fit
(during training). The :meth:~mne.decoding.Scaler.transform method is
called to transform data (training or test set) by scaling all time points
and epochs on a channel-by-channel basis. To perform both the fit and
transform operations in a single call, the
:meth:~mne.decoding.Scaler.fit_transform method may be used. To invert the
transform, :meth:~mne.decoding.Scaler.inverse_transform can be used. For
scalings='median', scikit-learn_ version 0.17+ is required.
<div class="alert alert-info"><h4>Note</h4><p>Using this class is different from directly applying
          :class:`sklearn.preprocessing.StandardScaler` or
          :class:`sklearn.preprocessing.RobustScaler` offered by
          scikit-learn_. These scale each *classification feature*, e.g.
          each time point for each channel, with mean and standard
          deviation computed across epochs, whereas
          :class:`mne.decoding.Scaler` scales each *channel* using mean and
          standard deviation computed across all of its time points
          and epochs.</p></div>

Vectorizer
^^^^^^^^^^
Scikit-learn API provides functionality to chain transformers and estimators
by using :class:sklearn.pipeline.Pipeline. We can construct decoding
pipelines and perform cross-validation and grid-search. However scikit-learn
transformers and estimators generally expect 2D data
(n_samples * n_features), whereas MNE transformers typically output data
with a higher dimensionality
(e.g. n_samples * n_channels * n_frequencies * n_times). A Vectorizer
therefore needs to be applied between the MNE and the scikit-learn steps
like:
End of explanation
"""


csp = CSP(n_components=3, norm_trace=False)
clf = make_pipeline(csp, LogisticRegression(solver='lbfgs'))
scores = cross_val_multiscore(clf, X, y, cv=5, n_jobs=1)
print('CSP: %0.1f%%' % (100 * scores.mean(),))

"""
Explanation: PSDEstimator
^^^^^^^^^^^^
The :class:mne.decoding.PSDEstimator
computes the power spectral density (PSD) using the multitaper
method. It takes a 3D array as input, converts it into 2D and computes the
PSD.
FilterEstimator
^^^^^^^^^^^^^^^
The :class:mne.decoding.FilterEstimator filters the 3D epochs data.
Spatial filters
Just like temporal filters, spatial filters provide weights to modify the
data along the sensor dimension. They are popular in the BCI community
because of their simplicity and ability to distinguish spatially-separated
neural activity.
Common spatial pattern
^^^^^^^^^^^^^^^^^^^^^^
:class:mne.decoding.CSP is a technique to analyze multichannel data based
on recordings from two classes [2]_ (see also
http://en.wikipedia.org/wiki/Common_spatial_pattern).
Let $X \in R^{C\times T}$ be a segment of data with
$C$ channels and $T$ time points. The data at a single time point
is denoted by $x(t)$ such that $X=[x(t), x(t+1), ..., x(t+T-1)]$.
Common spatial pattern (CSP) finds a decomposition that projects the signal
in the original sensor space to CSP space using the following transformation:
\begin{align}x_{CSP}(t) = W^{T}x(t)
   :label: csp\end{align}
where each column of $W \in R^{C\times C}$ is a spatial filter and each
row of $x_{CSP}$ is a CSP component. The matrix $W$ is also
called the de-mixing matrix in other contexts. Let
$\Sigma^{+} \in R^{C\times C}$ and $\Sigma^{-} \in R^{C\times C}$
be the estimates of the covariance matrices of the two conditions.
CSP analysis is given by the simultaneous diagonalization of the two
covariance matrices
\begin{align}W^{T}\Sigma^{+}W = \lambda^{+}
   :label: diagonalize_p\end{align}
\begin{align}W^{T}\Sigma^{-}W = \lambda^{-}
   :label: diagonalize_n\end{align}
where $\lambda^{C}$ is a diagonal matrix whose entries are the
eigenvalues of the following generalized eigenvalue problem
\begin{align}\Sigma^{+}w = \lambda \Sigma^{-}w
   :label: eigen_problem\end{align}
Large entries in the diagonal matrix corresponds to a spatial filter which
gives high variance in one class but low variance in the other. Thus, the
filter facilitates discrimination between the two classes.
.. topic:: Examples
* `sphx_glr_auto_examples_decoding_plot_decoding_csp_eeg.py`
* `sphx_glr_auto_examples_decoding_plot_decoding_csp_timefreq.py`

<div class="alert alert-info"><h4>Note</h4><p>The winning entry of the Grasp-and-lift EEG competition in Kaggle used
    the :class:`~mne.decoding.CSP` implementation in MNE and was featured as
    a `script of the week`_.</p></div>

We can use CSP with these data with:
End of explanation
"""


# Fit CSP on full data and plot
csp.fit(X, y)
csp.plot_patterns(epochs.info)
csp.plot_filters(epochs.info, scalings=1e-9)

"""
Explanation: Source power comodulation (SPoC)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Source Power Comodulation (:class:mne.decoding.SPoC) [3]_
identifies the composition of
orthogonal spatial filters that maximally correlate with a continuous target.
SPoC can be seen as an extension of the CSP where the target is driven by a
continuous variable rather than a discrete variable. Typical applications
include extraction of motor patterns using EMG power or audio patterns using
sound envelope.
.. topic:: Examples
* `sphx_glr_auto_examples_decoding_plot_decoding_spoc_CMC.py`

xDAWN
^^^^^
:class:mne.preprocessing.Xdawn is a spatial filtering method designed to
improve the signal to signal + noise ratio (SSNR) of the ERP responses [4]_.
Xdawn was originally
designed for P300 evoked potential by enhancing the target response with
respect to the non-target response. The implementation in MNE-Python is a
generalization to any type of ERP.
.. topic:: Examples
* `sphx_glr_auto_examples_preprocessing_plot_xdawn_denoising.py`
* `sphx_glr_auto_examples_decoding_plot_decoding_xdawn_eeg.py`

Effect-matched spatial filtering
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
The result of :class:mne.decoding.EMS is a spatial filter at each time
point and a corresponding time course [5]_.
Intuitively, the result gives the similarity between the filter at
each time point and the data vector (sensors) at that time point.
.. topic:: Examples
* `sphx_glr_auto_examples_decoding_plot_ems_filtering.py`

Patterns vs. filters
^^^^^^^^^^^^^^^^^^^^
When interpreting the components of the CSP (or spatial filters in general),
it is often more intuitive to think about how $x(t)$ is composed of
the different CSP components $x_{CSP}(t)$. In other words, we can
rewrite Equation :eq:csp as follows:
\begin{align}x(t) = (W^{-1})^{T}x_{CSP}(t)
   :label: patterns\end{align}
The columns of the matrix $(W^{-1})^T$ are called spatial patterns.
This is also called the mixing matrix. The example
sphx_glr_auto_examples_decoding_plot_linear_model_patterns.py
discusses the difference between patterns and filters.
These can be plotted with:
End of explanation
"""


# We will train the classifier on all left visual vs auditory trials on MEG

clf = make_pipeline(StandardScaler(), LogisticRegression(solver='lbfgs'))

time_decod = SlidingEstimator(clf, n_jobs=1, scoring='roc_auc', verbose=True)
scores = cross_val_multiscore(time_decod, X, y, cv=5, n_jobs=1)

# Mean scores across cross-validation splits
scores = np.mean(scores, axis=0)

# Plot
fig, ax = plt.subplots()
ax.plot(epochs.times, scores, label='score')
ax.axhline(.5, color='k', linestyle='--', label='chance')
ax.set_xlabel('Times')
ax.set_ylabel('AUC')  # Area Under the Curve
ax.legend()
ax.axvline(.0, color='k', linestyle='-')
ax.set_title('Sensor space decoding')

"""
Explanation: Decoding over time
This strategy consists in fitting a multivariate predictive model on each
time instant and evaluating its performance at the same instant on new
epochs. The :class:mne.decoding.SlidingEstimator will take as input a
pair of features $X$ and targets $y$, where $X$ has
more than 2 dimensions. For decoding over time the data $X$
is the epochs data of shape n_epochs x n_channels x n_times. As the
last dimension of $X$ is the time, an estimator will be fit
on every time instant.
This approach is analogous to SlidingEstimator-based approaches in fMRI,
where here we are interested in when one can discriminate experimental
conditions and therefore figure out when the effect of interest happens.
When working with linear models as estimators, this approach boils
down to estimating a discriminative spatial filter for each time instant.
Temporal decoding
^^^^^^^^^^^^^^^^^
We'll use a Logistic Regression for a binary classification as machine
learning model.
End of explanation
"""


clf = make_pipeline(StandardScaler(),
                    LinearModel(LogisticRegression(solver='lbfgs')))
time_decod = SlidingEstimator(clf, n_jobs=1, scoring='roc_auc', verbose=True)
time_decod.fit(X, y)

coef = get_coef(time_decod, 'patterns_', inverse_transform=True)
evoked = mne.EvokedArray(coef, epochs.info, tmin=epochs.times[0])
joint_kwargs = dict(ts_args=dict(time_unit='s'),
                    topomap_args=dict(time_unit='s'))
evoked.plot_joint(times=np.arange(0., .500, .100), title='patterns',
                  **joint_kwargs)

"""
Explanation: You can retrieve the spatial filters and spatial patterns if you explicitly
use a LinearModel
End of explanation
"""


# define the Temporal generalization object
time_gen = GeneralizingEstimator(clf, n_jobs=1, scoring='roc_auc',
                                 verbose=True)

scores = cross_val_multiscore(time_gen, X, y, cv=5, n_jobs=1)

# Mean scores across cross-validation splits
scores = np.mean(scores, axis=0)

# Plot the diagonal (it's exactly the same as the time-by-time decoding above)
fig, ax = plt.subplots()
ax.plot(epochs.times, np.diag(scores), label='score')
ax.axhline(.5, color='k', linestyle='--', label='chance')
ax.set_xlabel('Times')
ax.set_ylabel('AUC')
ax.legend()
ax.axvline(.0, color='k', linestyle='-')
ax.set_title('Decoding MEG sensors over time')

"""
Explanation: Temporal generalization
^^^^^^^^^^^^^^^^^^^^^^^
Temporal generalization is an extension of the decoding over time approach.
It consists in evaluating whether the model estimated at a particular
time instant accurately predicts any other time instant. It is analogous to
transferring a trained model to a distinct learning problem, where the
problems correspond to decoding the patterns of brain activity recorded at
distinct time instants.
The object to for Temporal generalization is
:class:mne.decoding.GeneralizingEstimator. It expects as input $X$
and $y$ (similarly to :class:~mne.decoding.SlidingEstimator) but
generates predictions from each model for all time instants. The class
:class:~mne.decoding.GeneralizingEstimator is generic and will treat the
last dimension as the one to be used for generalization testing. For
convenience, here, we refer to it as different tasks. If $X$
corresponds to epochs data then the last dimension is time.
This runs the analysis used in [6] and further detailed in [7]:
End of explanation
"""


fig, ax = plt.subplots(1, 1)
im = ax.imshow(scores, interpolation='lanczos', origin='lower', cmap='RdBu_r',
               extent=epochs.times[[0, -1, 0, -1]], vmin=0., vmax=1.)
ax.set_xlabel('Testing Time (s)')
ax.set_ylabel('Training Time (s)')
ax.set_title('Temporal generalization')
ax.axvline(0, color='k')
ax.axhline(0, color='k')
plt.colorbar(im, ax=ax)

"""
Explanation: Plot the full (generalization) matrix:
End of explanation
"""


