{
  "nbformat_minor": 0, 
  "nbformat": 4, 
  "cells": [
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "%matplotlib inline"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "\n# From raw data to dSPM on SPM Faces dataset\n\n\nRuns a full pipeline using MNE-Python:\n\n    - artifact removal\n    - averaging Epochs\n    - forward model computation\n    - source reconstruction using dSPM on the contrast : \"faces - scrambled\"\n\n<div class=\"alert alert-info\"><h4>Note</h4><p>This example does quite a bit of processing, so even on a\n          fast machine it can take several minutes to complete.</p></div>\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "# Authors: Alexandre Gramfort <alexandre.gramfort@telecom-paristech.fr>\n#          Denis Engemann <denis.engemann@gmail.com>\n#\n# License: BSD (3-clause)\n\nimport os.path as op\nimport matplotlib.pyplot as plt\n\nimport mne\nfrom mne.datasets import spm_face\nfrom mne.preprocessing import ICA, create_eog_epochs\nfrom mne import io, combine_evoked\nfrom mne.minimum_norm import make_inverse_operator, apply_inverse\n\nprint(__doc__)\n\ndata_path = spm_face.data_path()\nsubjects_dir = data_path + '/subjects'"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Load and filter data, set up epochs\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "raw_fname = data_path + '/MEG/spm/SPM_CTF_MEG_example_faces%d_3D.ds'\n\nraw = io.read_raw_ctf(raw_fname % 1, preload=True)  # Take first run\n# Here to save memory and time we'll downsample heavily -- this is not\n# advised for real data as it can effectively jitter events!\nraw.resample(120., npad='auto')\n\npicks = mne.pick_types(raw.info, meg=True, exclude='bads')\nraw.filter(1, 30, method='iir')\n\nevents = mne.find_events(raw, stim_channel='UPPT001')\n\n# plot the events to get an idea of the paradigm\nmne.viz.plot_events(events, raw.info['sfreq'])\n\nevent_ids = {\"faces\": 1, \"scrambled\": 2}\n\ntmin, tmax = -0.2, 0.6\nbaseline = None  # no baseline as high-pass is applied\nreject = dict(mag=5e-12)\n\nepochs = mne.Epochs(raw, events, event_ids, tmin, tmax,  picks=picks,\n                    baseline=baseline, preload=True, reject=reject)\n\n# Fit ICA, find and remove major artifacts\nica = ICA(n_components=0.95, random_state=0).fit(raw, decim=1, reject=reject)\n\n# compute correlation scores, get bad indices sorted by score\neog_epochs = create_eog_epochs(raw, ch_name='MRT31-2908', reject=reject)\neog_inds, eog_scores = ica.find_bads_eog(eog_epochs, ch_name='MRT31-2908')\nica.plot_scores(eog_scores, eog_inds)  # see scores the selection is based on\nica.plot_components(eog_inds)  # view topographic sensitivity of components\nica.exclude += eog_inds[:1]  # we saw the 2nd ECG component looked too dipolar\nica.plot_overlay(eog_epochs.average())  # inspect artifact removal\nica.apply(epochs)  # clean data, default in place\n\nevoked = [epochs[k].average() for k in event_ids]\n\ncontrast = combine_evoked(evoked, weights=[-1, 1])  # Faces - scrambled\n\nevoked.append(contrast)\n\nfor e in evoked:\n    e.plot(ylim=dict(mag=[-400, 400]))\n\nplt.show()\n\n# estimate noise covarariance\nnoise_cov = mne.compute_covariance(epochs, tmax=0, method='shrunk')"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Visualize fields on MEG helmet\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "trans_fname = data_path + ('/MEG/spm/SPM_CTF_MEG_example_faces1_3D_'\n                           'raw-trans.fif')\n\nmaps = mne.make_field_map(evoked[0], trans_fname, subject='spm',\n                          subjects_dir=subjects_dir, n_jobs=1)\n\nevoked[0].plot_field(maps, time=0.170)"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Compute forward model\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "# Make source space\nsrc_fname = data_path + '/subjects/spm/bem/spm-oct-6-src.fif'\nif not op.isfile(src_fname):\n    src = mne.setup_source_space('spm', src_fname, spacing='oct6',\n                                 subjects_dir=subjects_dir, overwrite=True)\nelse:\n    src = mne.read_source_spaces(src_fname)\n\nbem = data_path + '/subjects/spm/bem/spm-5120-5120-5120-bem-sol.fif'\nforward = mne.make_forward_solution(contrast.info, trans_fname, src, bem)\nforward = mne.convert_forward_solution(forward, surf_ori=True)"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Compute inverse solution\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "snr = 3.0\nlambda2 = 1.0 / snr ** 2\nmethod = 'dSPM'\n\ninverse_operator = make_inverse_operator(contrast.info, forward, noise_cov,\n                                         loose=0.2, depth=0.8)\n\n# Compute inverse solution on contrast\nstc = apply_inverse(contrast, inverse_operator, lambda2, method, pick_ori=None)\n# stc.save('spm_%s_dSPM_inverse' % constrast.comment)\n\n# Plot contrast in 3D with PySurfer if available\nbrain = stc.plot(hemi='both', subjects_dir=subjects_dir, initial_time=0.170,\n                 views=['ven'])\n# brain.save_image('dSPM_map.png')"
      ], 
      "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"
      }
    }
  }
}