{
  "nbformat_minor": 0, 
  "nbformat": 4, 
  "cells": [
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "%matplotlib inline"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "\n# Decoding source space data\n\n\nDecoding, a.k.a MVPA or supervised machine learning applied to MEG\ndata in source space on the left cortical surface. Here f-test feature\nselection is employed to confine the classification to the potentially\nrelevant features. The classifier then is trained to selected features of\nepochs in source space.\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "# Author: Denis A. Engemann <denis.engemann@gmail.com>\n#         Alexandre Gramfort <alexandre.gramfort@telecom-paristech.fr>\n#\n# License: BSD (3-clause)\n\nimport mne\nimport os\nimport numpy as np\nfrom mne import io\nfrom mne.datasets import sample\nfrom mne.minimum_norm import apply_inverse_epochs, read_inverse_operator\n\nprint(__doc__)\n\ndata_path = sample.data_path()\nfname_fwd = data_path + 'MEG/sample/sample_audvis-meg-oct-6-fwd.fif'\nfname_evoked = data_path + '/MEG/sample/sample_audvis-ave.fif'\nsubjects_dir = data_path + '/subjects'\nsubject = os.environ['SUBJECT'] = subjects_dir + '/sample'\nos.environ['SUBJECTS_DIR'] = subjects_dir"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Set parameters\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "raw_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw.fif'\nevent_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw-eve.fif'\nfname_cov = data_path + '/MEG/sample/sample_audvis-cov.fif'\nlabel_names = 'Aud-rh', 'Vis-rh'\nfname_inv = data_path + '/MEG/sample/sample_audvis-meg-oct-6-meg-inv.fif'\n\ntmin, tmax = -0.2, 0.5\nevent_id = dict(aud_r=2, vis_r=4)  # load contra-lateral conditions\n\n# Setup for reading the raw data\nraw = io.read_raw_fif(raw_fname, preload=True)\nraw.filter(2, None)  # replace baselining with high-pass\nevents = mne.read_events(event_fname)\n\n# Set up pick list: MEG - bad channels (modify to your needs)\nraw.info['bads'] += ['MEG 2443']  # mark bads\npicks = mne.pick_types(raw.info, meg=True, eeg=False, stim=True, eog=True,\n                       exclude='bads')\n\n# Read epochs\nepochs = mne.Epochs(raw, events, event_id, tmin, tmax, proj=True,\n                    picks=picks, baseline=None, preload=True,\n                    reject=dict(grad=4000e-13, eog=150e-6),\n                    decim=5)  # decimate to save memory and increase speed\n\nepochs.equalize_event_counts(list(event_id.keys()))\nepochs_list = [epochs[k] for k in event_id]\n\n# Compute inverse solution\nsnr = 3.0\nlambda2 = 1.0 / snr ** 2\nmethod = \"dSPM\"  # use dSPM method (could also be MNE or sLORETA)\nn_times = len(epochs.times)\nn_vertices = 3732\nn_epochs = len(epochs.events)\n\n# Load data and compute inverse solution and stcs for each epoch.\n\nnoise_cov = mne.read_cov(fname_cov)\ninverse_operator = read_inverse_operator(fname_inv)\nX = np.zeros([n_epochs, n_vertices, n_times])\n\n# to save memory, we'll load and transform our epochs step by step.\nfor condition_count, ep in zip([0, n_epochs // 2], epochs_list):\n    stcs = apply_inverse_epochs(ep, inverse_operator, lambda2,\n                                method, pick_ori=\"normal\",  # saves us memory\n                                return_generator=True)\n    for jj, stc in enumerate(stcs):\n        X[condition_count + jj] = stc.lh_data"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Decoding in sensor space using a linear SVM\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "# Make arrays X and y such that :\n# X is 3d with X.shape[0] is the total number of epochs to classify\n# y is filled with integers coding for the class to predict\n# We must have X.shape[0] equal to y.shape[0]\n\n# we know the first half belongs to the first class, the second one\ny = np.repeat([0, 1], len(X) / 2)   # belongs to the second class\nX = X.reshape(n_epochs, n_vertices * n_times)\n# we have to normalize the data before supplying them to our classifier\nX -= X.mean(axis=0)\nX /= X.std(axis=0)\n\n# prepare classifier\nfrom sklearn.svm import SVC  # noqa\nfrom sklearn.cross_validation import ShuffleSplit  # noqa\n\n# Define a monte-carlo cross-validation generator (reduce variance):\nn_splits = 10\nclf = SVC(C=1, kernel='linear')\ncv = ShuffleSplit(len(X), n_splits, test_size=0.2)\n\n# setup feature selection and classification pipeline\nfrom sklearn.feature_selection import SelectKBest, f_classif  # noqa\nfrom sklearn.pipeline import Pipeline  # noqa\n\n# we will use an ANOVA f-test to preselect relevant spatio-temporal units\nfeature_selection = SelectKBest(f_classif, k=500)  # take the best 500\n# to make life easier we will create a pipeline object\nanova_svc = Pipeline([('anova', feature_selection), ('svc', clf)])\n\n# initialize score and feature weights result arrays\nscores = np.zeros(n_splits)\nfeature_weights = np.zeros([n_vertices, n_times])\n\n# hold on, this may take a moment\nfor ii, (train, test) in enumerate(cv):\n    anova_svc.fit(X[train], y[train])\n    y_pred = anova_svc.predict(X[test])\n    y_test = y[test]\n    scores[ii] = np.sum(y_pred == y_test) / float(len(y_test))\n    feature_weights += feature_selection.inverse_transform(clf.coef_) \\\n        .reshape(n_vertices, n_times)\n\nprint('Average prediction accuracy: %0.3f | standard deviation:  %0.3f'\n      % (scores.mean(), scores.std()))\n\n# prepare feature weights for visualization\nfeature_weights /= (ii + 1)  # create average weights\n# create mask to avoid division error\nfeature_weights = np.ma.masked_array(feature_weights, feature_weights == 0)\n# normalize scores for visualization purposes\nfeature_weights /= feature_weights.std(axis=1)[:, None]\nfeature_weights -= feature_weights.mean(axis=1)[:, None]\n\n# unmask, take absolute values, emulate f-value scale\nfeature_weights = np.abs(feature_weights.data) * 10\n\nvertices = [stc.lh_vertno, np.array([], int)]  # empty array for right hemi\nstc_feat = mne.SourceEstimate(feature_weights, vertices=vertices,\n                              tmin=stc.tmin, tstep=stc.tstep,\n                              subject='sample')\n\nbrain = stc_feat.plot(views=['lat'], transparent=True,\n                      initial_time=0.1, time_unit='s')"
      ], 
      "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"
      }
    }
  }
}