{
  "nbformat_minor": 0, 
  "nbformat": 4, 
  "cells": [
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "%matplotlib inline"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "\n# Compute mixed source space connectivity and visualize it using a circular graph\n\n\nThis example computes the all-to-all connectivity between 75 regions in\na mixed source space based on dSPM inverse solutions and a FreeSurfer cortical\nparcellation. The connectivity is visualized using a circular graph which\nis ordered based on the locations of the regions.\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "# Author: Annalisa Pascarella <a.pascarella@iac.cnr.it>\n#\n# License: BSD (3-clause)\n\nimport os.path as op\nimport numpy as np\nimport mne\n\nfrom mne.datasets import sample\nfrom mne import setup_volume_source_space, setup_source_space\nfrom mne import make_forward_solution\nfrom mne.io import read_raw_fif\nfrom mne.minimum_norm import make_inverse_operator, apply_inverse_epochs\nfrom mne.connectivity import spectral_connectivity\nfrom mne.viz import circular_layout, plot_connectivity_circle\n\n# Set dir\ndata_path = sample.data_path()\nsubject = 'sample'\ndata_dir = op.join(data_path, 'MEG', subject)\nsubjects_dir = op.join(data_path, 'subjects')\nbem_dir = op.join(subjects_dir, subject, 'bem')\n\n# Set file names\nfname_aseg = op.join(subjects_dir, subject, 'mri', 'aseg.mgz')\n\nfname_model = op.join(bem_dir, '%s-5120-bem.fif' % subject)\nfname_bem = op.join(bem_dir, '%s-5120-bem-sol.fif' % subject)\n\nfname_raw = data_dir + '/sample_audvis_filt-0-40_raw.fif'\nfname_trans = data_dir + '/sample_audvis_raw-trans.fif'\nfname_cov = data_dir + '/ernoise-cov.fif'\nfname_event = data_dir + '/sample_audvis_filt-0-40_raw-eve.fif'\n\n# List of sub structures we are interested in. We select only the\n# sub structures we want to include in the source space\nlabels_vol = ['Left-Amygdala',\n              'Left-Thalamus-Proper',\n              'Left-Cerebellum-Cortex',\n              'Brain-Stem',\n              'Right-Amygdala',\n              'Right-Thalamus-Proper',\n              'Right-Cerebellum-Cortex']\n\n# Setup a surface-based source space\nsrc = setup_source_space(subject, fname=None, subjects_dir=subjects_dir,\n                         spacing='oct6', add_dist=False)\n\n# Setup a volume source space\n# set pos=7.0 for speed issue\nvol_src = setup_volume_source_space(subject, mri=fname_aseg,\n                                    pos=7.0,\n                                    bem=fname_model,\n                                    volume_label=labels_vol,\n                                    subjects_dir=subjects_dir)\n# Generate the mixed source space\nsrc += vol_src\n\n# compute the fwd matrix\nfwd = make_forward_solution(fname_raw, fname_trans, src, fname_bem,\n                            mindist=5.0,  # ignore sources<=5mm from innerskull\n                            meg=True, eeg=False,\n                            n_jobs=1)\n\n# Load data\nraw = read_raw_fif(fname_raw, preload=True)\nnoise_cov = mne.read_cov(fname_cov)\nevents = mne.read_events(fname_event)\n\n# Add a bad channel\nraw.info['bads'] += ['MEG 2443']\n\n# Pick MEG channels\npicks = mne.pick_types(raw.info, meg=True, eeg=False, stim=False, eog=True,\n                       exclude='bads')\n\n# Define epochs for left-auditory condition\nevent_id, tmin, tmax = 1, -0.2, 0.5\nepochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,\n                    baseline=(None, 0), reject=dict(mag=4e-12, grad=4000e-13,\n                                                    eog=150e-6))\n\n# Compute inverse solution and for each epoch\nsnr = 1.0           # use smaller SNR for raw data\ninv_method = 'dSPM'  # sLORETA, MNE, dSPM\nparc = 'aparc'      # the parcellation to use, e.g., 'aparc' 'aparc.a2009s'\n\nlambda2 = 1.0 / snr ** 2\n\n# Compute inverse operator\ninverse_operator = make_inverse_operator(raw.info, fwd, noise_cov,\n                                         loose=None, depth=None,\n                                         fixed=False)\n\n\nstcs = apply_inverse_epochs(epochs, inverse_operator, lambda2, inv_method,\n                            pick_ori=None, return_generator=True)\n\n# Get labels for FreeSurfer 'aparc' cortical parcellation with 34 labels/hemi\nlabels_parc = mne.read_labels_from_annot(subject, parc=parc,\n                                         subjects_dir=subjects_dir)\n\n# Average the source estimates within each label of the cortical parcellation\n# and each sub structures contained in the src space\n# If mode = 'mean_flip' this option is used only for the cortical label\nsrc = inverse_operator['src']\nlabel_ts = mne.extract_label_time_course(stcs, labels_parc, src,\n                                         mode='mean_flip',\n                                         allow_empty=True,\n                                         return_generator=False)\n\n# We compute the connectivity in the alpha band and plot it using a circular\n# graph layout\nfmin = 8.\nfmax = 13.\nsfreq = raw.info['sfreq']  # the sampling frequency\ncon, freqs, times, n_epochs, n_tapers = spectral_connectivity(\n    label_ts, method='pli', mode='multitaper', sfreq=sfreq, fmin=fmin,\n    fmax=fmax, faverage=True, mt_adaptive=True, n_jobs=1)\n\n# We create a list of Label containing also the sub structures\nlabels_aseg = mne.get_volume_labels_from_src(src, subject, subjects_dir)\nlabels = labels_parc + labels_aseg\n\n# read colors\nnode_colors = [label.color for label in labels]\n\n# We reorder the labels based on their location in the left hemi\nlabel_names = [label.name for label in labels]\nlh_labels = [name for name in label_names if name.endswith('lh')]\nrh_labels = [name for name in label_names if name.endswith('rh')]\n\n# Get the y-location of the label\nlabel_ypos_lh = list()\nfor name in lh_labels:\n    idx = label_names.index(name)\n    ypos = np.mean(labels[idx].pos[:, 1])\n    label_ypos_lh.append(ypos)\ntry:\n    idx = label_names.index('Brain-Stem')\n    ypos = np.mean(labels[idx].pos[:, 1])\n    lh_labels.append('Brain-Stem')\n    label_ypos_lh.append(ypos)\nexcept ValueError:\n    pass\n\n\n# Reorder the labels based on their location\nlh_labels = [label for (yp, label) in sorted(zip(label_ypos_lh, lh_labels))]\n\n# For the right hemi\nrh_labels = [label[:-2] + 'rh' for label in lh_labels\n             if label != 'Brain-Stem' and label[:-2] + 'rh' in rh_labels]\n\n# Save the plot order\nnode_order = list()\nnode_order = lh_labels[::-1] + rh_labels\n\nnode_angles = circular_layout(label_names, node_order, start_pos=90,\n                              group_boundaries=[0, len(label_names) // 2])\n\n\n# Plot the graph using node colors from the FreeSurfer parcellation. We only\n# show the 300 strongest connections.\nconmat = con[:, :, 0]\nplot_connectivity_circle(conmat, label_names, n_lines=300,\n                         node_angles=node_angles, node_colors=node_colors,\n                         title='All-to-All Connectivity left-Auditory '\n                               'Condition (PLI)')\n\n# Uncomment the following line to save the figure\n'''\nimport matplotlib.pyplot as plt\nplt.savefig('circle.png', facecolor='black')\n'''"
      ], 
      "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"
      }
    }
  }
}