{
  "nbformat_minor": 0, 
  "nbformat": 4, 
  "cells": [
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "%matplotlib inline"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "\n===========================================================================\nMotor imagery decoding from EEG data using the Common Spatial Pattern (CSP)\n===========================================================================\n\nDecoding of motor imagery applied to EEG data decomposed using CSP.\nHere the classifier is applied to features extracted on CSP filtered signals.\n\nSee http://en.wikipedia.org/wiki/Common_spatial_pattern and [1]_. The EEGBCI\ndataset is documented in [2]_. The data set is available at PhysioNet [3]_.\n\nReferences\n----------\n\n.. [1] Zoltan J. Koles. The quantitative extraction and topographic mapping\n       of the abnormal components in the clinical EEG. Electroencephalography\n       and Clinical Neurophysiology, 79(6):440--447, December 1991.\n.. [2] Schalk, G., McFarland, D.J., Hinterberger, T., Birbaumer, N.,\n       Wolpaw, J.R. (2004) BCI2000: A General-Purpose Brain-Computer Interface\n       (BCI) System. IEEE TBME 51(6):1034-1043.\n.. [3] Goldberger AL, Amaral LAN, Glass L, Hausdorff JM, Ivanov PCh, Mark RG,\n       Mietus JE, Moody GB, Peng C-K, Stanley HE. (2000) PhysioBank,\n       PhysioToolkit, and PhysioNet: Components of a New Research Resource for\n       Complex Physiologic Signals. Circulation 101(23):e215-e220.\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "# Authors: Martin Billinger <martin.billinger@tugraz.at>\n#\n# License: BSD (3-clause)\n\nimport numpy as np\nimport matplotlib.pyplot as plt\n\nfrom mne import Epochs, pick_types, find_events\nfrom mne.channels import read_layout\nfrom mne.io import concatenate_raws, read_raw_edf\nfrom mne.datasets import eegbci\nfrom mne.decoding import CSP\n\nprint(__doc__)\n\n# #############################################################################\n# # Set parameters and read data\n\n# avoid classification of evoked responses by using epochs that start 1s after\n# cue onset.\ntmin, tmax = -1., 4.\nevent_id = dict(hands=2, feet=3)\nsubject = 1\nruns = [6, 10, 14]  # motor imagery: hands vs feet\n\nraw_fnames = eegbci.load_data(subject, runs)\nraw_files = [read_raw_edf(f, preload=True) for f in raw_fnames]\nraw = concatenate_raws(raw_files)\n\n# strip channel names of \".\" characters\nraw.rename_channels(lambda x: x.strip('.'))\n\n# Apply band-pass filter\nraw.filter(7., 30.)\n\nevents = find_events(raw, shortest_event=0, stim_channel='STI 014')\n\npicks = pick_types(raw.info, meg=False, eeg=True, stim=False, eog=False,\n                   exclude='bads')\n\n# Read epochs (train will be done only between 1 and 2s)\n# Testing will be done with a running classifier\nepochs = Epochs(raw, events, event_id, tmin, tmax, proj=True, picks=picks,\n                baseline=None, preload=True)\nepochs_train = epochs.copy().crop(tmin=1., tmax=2.)\nlabels = epochs.events[:, -1] - 2"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Classification with linear discrimant analysis\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "from sklearn.lda import LDA  # noqa\nfrom sklearn.cross_validation import ShuffleSplit  # noqa\n\n# Assemble a classifier\nlda = LDA()\ncsp = CSP(n_components=4, reg=None, log=True)\n\n# Define a monte-carlo cross-validation generator (reduce variance):\ncv = ShuffleSplit(len(labels), 10, test_size=0.2, random_state=42)\nscores = []\nepochs_data = epochs.get_data()\nepochs_data_train = epochs_train.get_data()\n\n# Use scikit-learn Pipeline with cross_val_score function\nfrom sklearn.pipeline import Pipeline  # noqa\nfrom sklearn.cross_validation import cross_val_score  # noqa\nclf = Pipeline([('CSP', csp), ('LDA', lda)])\nscores = cross_val_score(clf, epochs_data_train, labels, cv=cv, n_jobs=1)\n\n# Printing the results\nclass_balance = np.mean(labels == labels[0])\nclass_balance = max(class_balance, 1. - class_balance)\nprint(\"Classification accuracy: %f / Chance level: %f\" % (np.mean(scores),\n                                                          class_balance))\n\n# plot CSP patterns estimated on full data for visualization\ncsp.fit_transform(epochs_data, labels)\n\nevoked = epochs.average()\nevoked.data = csp.patterns_.T\nevoked.times = np.arange(evoked.data.shape[0])\n\nlayout = read_layout('EEG1005')\nevoked.plot_topomap(times=[0, 1, 2, 3, 4, 5], ch_type='eeg', layout=layout,\n                    scale_time=1, time_format='%i', scale=1,\n                    unit='Patterns (AU)', size=1.5)"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Look at performance over time\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "sfreq = raw.info['sfreq']\nw_length = int(sfreq * 0.5)   # running classifier: window length\nw_step = int(sfreq * 0.1)  # running classifier: window step size\nw_start = np.arange(0, epochs_data.shape[2] - w_length, w_step)\n\nscores_windows = []\n\nfor train_idx, test_idx in cv:\n    y_train, y_test = labels[train_idx], labels[test_idx]\n\n    X_train = csp.fit_transform(epochs_data_train[train_idx], y_train)\n    X_test = csp.transform(epochs_data_train[test_idx])\n\n    # fit classifier\n    lda.fit(X_train, y_train)\n\n    # running classifier: test classifier on sliding window\n    score_this_window = []\n    for n in w_start:\n        X_test = csp.transform(epochs_data[test_idx][:, :, n:(n + w_length)])\n        score_this_window.append(lda.score(X_test, y_test))\n    scores_windows.append(score_this_window)\n\n# Plot scores over time\nw_times = (w_start + w_length / 2.) / sfreq + epochs.tmin\n\nplt.figure()\nplt.plot(w_times, np.mean(scores_windows, 0), label='Score')\nplt.axvline(0, linestyle='--', color='k', label='Onset')\nplt.axhline(0.5, linestyle='-', color='k', label='Chance')\nplt.xlabel('time (s)')\nplt.ylabel('classification accuracy')\nplt.title('Classification score over time')\nplt.legend(loc='lower right')\nplt.show()"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }
  ], 
  "metadata": {
    "kernelspec": {
      "display_name": "Python 2", 
      "name": "python2", 
      "language": "python"
    }, 
    "language_info": {
      "mimetype": "text/x-python", 
      "nbconvert_exporter": "python", 
      "name": "python", 
      "file_extension": ".py", 
      "version": "2.7.13", 
      "pygments_lexer": "ipython2", 
      "codemirror_mode": {
        "version": 2, 
        "name": "ipython"
      }
    }
  }
}