Compare Riemann to classic in a sliding window context

Large_sliding_window.py

import numpy as np
import mne
from mne import io
from mne.datasets import sample
data_path = sample.data_path()
raw_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw.fif'
event_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw-eve.fif'
raw = io.Raw(raw_fname, preload=True)
raw.filter(1, 20, method='iir')
events = mne.read_events(event_fname)
picks = mne.pick_types(raw.info, meg='grad', eeg=False, stim=True, eog=True,
                       exclude='bads')
epochs = mne.Epochs(raw, events, dict(aud_l=1, vis_l=3), -.4, 1.2, proj=True,
                    picks=picks, baseline=None, preload=True,
                    reject=dict(grad=4000e-13, eog=150e-6))

picks = mne.pick_types(epochs.info, meg=True, exclude='bads')
X = epochs._data[:, picks, :]
y = epochs.events[:, 2]

###############################################################################
from mne.decoding import GeneralizationAcrossTime
from sklearn.pipeline import make_pipeline
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler
from pyriemann.estimation import ERPCovariances
from pyriemann.tangentspace import TangentSpace


def scorer_auc(y_true, y_pred):
    from sklearn.metrics import roc_auc_score
    from sklearn.preprocessing import LabelBinarizer
    le = LabelBinarizer()
    y_true = le.fit_transform(y_true)
    return roc_auc_score(y_true, y_pred)


class LR_proba(LogisticRegression):
    def predict(self, X):
        # continuous prediction instead of categorical one
        return super(LR_proba, self).predict_proba(X)[:, 1]


class reshape_X(object):
    def __init__(self, window=1):
        self.window = window   # in time sample

    def fit(self, X, y=None):
        return self

    def fit_transform(self, X, y=None):
        return self.transform(X)

    def transform(self, X):
        if X.shape[1] / self.window != X.shape[1] / float(self.window):
            raise ValueError('Wrong window size')
        return X.reshape([X.shape[0], X.shape[1] / self.window, self.window])


# GAT parameters
window = 20
window_s = window / epochs.info['sfreq']  # in seconds
step = 5. / epochs.info['sfreq']

kwargs = dict(test_times='diagonal', scorer=scorer_auc, cv=4,
              train_times=dict(length=window_s, step=step))

# Fit & Score on a single validation
clf_classic = make_pipeline(StandardScaler(), LR_proba())
gat = GeneralizationAcrossTime(n_jobs=-1, clf=clf_classic, **kwargs)
gat.fit(epochs)
gat.score(epochs)

clf_rieman = make_pipeline(reshape_X(window=window),
                           ERPCovariances(estimator='lwf', svd=4),
                           TangentSpace(metric='logeuclid'), LR_proba())
gat_rieman = GeneralizationAcrossTime(n_jobs=1, clf=clf_rieman, **kwargs)
gat_rieman.fit(epochs)
gat_rieman.score(epochs)

# Plot
import matplotlib.pyplot as plt
from gat.utils import subscore
from gat.plot import plot_sem
fig, ax = plt.subplots(1)
fig = gat.plot_diagonal(label='Logistic', show=False, chance=.5, color='b',
                        ax=ax)
gat_rieman.plot_diagonal(label='Riemann', ax=ax, color='r',
                              show=False, chance=False)


scores = list()
for test, train in gat.cv_:
    scores.append(subscore(gat, test))
plot_sem(gat.train_times_['times'], np.squeeze(scores), ax=ax, color='b')

scores = list()
for test, train in gat_rieman.cv_:
    scores.append(subscore(gat_rieman, test))
plot_sem(gat.train_times_['times'], np.squeeze(scores), ax=ax, color='r')

fig.show()

@dengemann