{
  "nbformat_minor": 0, 
  "nbformat": 4, 
  "cells": [
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "%matplotlib inline"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "\n# Compute a sparse inverse solution using the Gamma-Map empirical Bayesian method\n\n\nSee Wipf et al. \"A unified Bayesian framework for MEG/EEG source imaging.\"\nNeuroImage, vol. 44, no. 3, pp. 947?66, Mar. 2009.\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "# Author: Martin Luessi <mluessi@nmr.mgh.harvard.edu>\n#\n# License: BSD (3-clause)\n\nimport numpy as np\n\nimport mne\nfrom mne.datasets import sample\nfrom mne.inverse_sparse import gamma_map\nfrom mne.viz import plot_sparse_source_estimates\n\nprint(__doc__)\n\ndata_path = sample.data_path()\nsubjects_dir = data_path + '/subjects'\nfwd_fname = data_path + '/MEG/sample/sample_audvis-meg-eeg-oct-6-fwd.fif'\nevoked_fname = data_path + '/MEG/sample/sample_audvis-ave.fif'\ncov_fname = data_path + '/MEG/sample/sample_audvis-cov.fif'\n\n# Read the evoked response and crop it\ncondition = 'Left visual'\nevoked = mne.read_evokeds(evoked_fname, condition=condition,\n                          baseline=(None, 0))\nevoked.crop(tmin=-50e-3, tmax=300e-3)\n\n# Read the forward solution\nforward = mne.read_forward_solution(fwd_fname, surf_ori=True,\n                                    force_fixed=False)\n\n# Read noise noise covariance matrix and regularize it\ncov = mne.read_cov(cov_fname)\ncov = mne.cov.regularize(cov, evoked.info)\n\n# Run the Gamma-MAP method\nalpha = 0.5\nstc, residual = gamma_map(evoked, forward, cov, alpha, xyz_same_gamma=True,\n                          return_residual=True)\n\n# View in 2D and 3D (\"glass\" brain like 3D plot)\n\n# Show the sources as spheres scaled by their strength\nscale_factors = np.max(np.abs(stc.data), axis=1)\nscale_factors = 0.5 * (1 + scale_factors / np.max(scale_factors))\n\nplot_sparse_source_estimates(\n    forward['src'], stc, bgcolor=(1, 1, 1),\n    modes=['sphere'], opacity=0.1, scale_factors=(scale_factors, None),\n    fig_name=\"Gamma-MAP\")\n\n# Show the evoked response and the residual for gradiometers\nylim = dict(grad=[-120, 120])\nevoked.pick_types(meg='grad', exclude='bads')\nevoked.plot(titles=dict(grad='Evoked Response Gradiometers'), ylim=ylim,\n            proj=True)\n\nresidual.pick_types(meg='grad', exclude='bads')\nresidual.plot(titles=dict(grad='Residuals Gradiometers'), ylim=ylim,\n              proj=True)"
      ], 
      "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"
      }
    }
  }
}