{
  "nbformat_minor": 0, 
  "nbformat": 4, 
  "cells": [
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "%matplotlib inline"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "\n# Compute all-to-all connectivity in sensor space\n\n\nComputes the Phase Lag Index (PLI) between all gradiometers and shows the\nconnectivity in 3D using the helmet geometry. The left visual stimulation data\nare used which produces strong connectvitiy in the right occipital sensors.\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\nfrom scipy import linalg\n\nimport mne\nfrom mne import io\nfrom mne.connectivity import spectral_connectivity\nfrom mne.datasets import sample\n\nprint(__doc__)"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Set parameters\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "data_path = sample.data_path()\nraw_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'\n\n# Setup for reading the raw data\nraw = io.read_raw_fif(raw_fname)\nevents = mne.read_events(event_fname)\n\n# Add a bad channel\nraw.info['bads'] += ['MEG 2443']\n\n# Pick MEG gradiometers\npicks = mne.pick_types(raw.info, meg='grad', eeg=False, stim=False, eog=True,\n                       exclude='bads')\n\n# Create epochs for the visual condition\nevent_id, tmin, tmax = 3, -0.2, 0.5\nepochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,\n                    baseline=(None, 0), reject=dict(grad=4000e-13, eog=150e-6))\n\n# Compute connectivity for band containing the evoked response.\n# We exclude the baseline period\nfmin, fmax = 3., 9.\nsfreq = raw.info['sfreq']  # the sampling frequency\ntmin = 0.0  # exclude the baseline period\ncon, freqs, times, n_epochs, n_tapers = spectral_connectivity(\n    epochs, method='pli', mode='multitaper', sfreq=sfreq, fmin=fmin, fmax=fmax,\n    faverage=True, tmin=tmin, mt_adaptive=False, n_jobs=1)\n\n# the epochs contain an EOG channel, which we remove now\nch_names = epochs.ch_names\nidx = [ch_names.index(name) for name in ch_names if name.startswith('MEG')]\ncon = con[idx][:, idx]\n\n# con is a 3D array where the last dimension is size one since we averaged\n# over frequencies in a single band. Here we make it 2D\ncon = con[:, :, 0]\n\n# Now, visualize the connectivity in 3D\nfrom mayavi import mlab  # noqa\n\nmlab.figure(size=(600, 600), bgcolor=(0.5, 0.5, 0.5))\n\n# Plot the sensor locations\nsens_loc = [raw.info['chs'][picks[i]]['loc'][:3] for i in idx]\nsens_loc = np.array(sens_loc)\n\npts = mlab.points3d(sens_loc[:, 0], sens_loc[:, 1], sens_loc[:, 2],\n                    color=(1, 1, 1), opacity=1, scale_factor=0.005)\n\n# Get the strongest connections\nn_con = 20  # show up to 20 connections\nmin_dist = 0.05  # exclude sensors that are less than 5cm apart\nthreshold = np.sort(con, axis=None)[-n_con]\nii, jj = np.where(con >= threshold)\n\n# Remove close connections\ncon_nodes = list()\ncon_val = list()\nfor i, j in zip(ii, jj):\n    if linalg.norm(sens_loc[i] - sens_loc[j]) > min_dist:\n        con_nodes.append((i, j))\n        con_val.append(con[i, j])\n\ncon_val = np.array(con_val)\n\n# Show the connections as tubes between sensors\nvmax = np.max(con_val)\nvmin = np.min(con_val)\nfor val, nodes in zip(con_val, con_nodes):\n    x1, y1, z1 = sens_loc[nodes[0]]\n    x2, y2, z2 = sens_loc[nodes[1]]\n    points = mlab.plot3d([x1, x2], [y1, y2], [z1, z2], [val, val],\n                         vmin=vmin, vmax=vmax, tube_radius=0.001,\n                         colormap='RdBu')\n    points.module_manager.scalar_lut_manager.reverse_lut = True\n\n\nmlab.scalarbar(title='Phase Lag Index (PLI)', nb_labels=4)\n\n# Add the sensor names for the connections shown\nnodes_shown = list(set([n[0] for n in con_nodes] +\n                       [n[1] for n in con_nodes]))\n\nfor node in nodes_shown:\n    x, y, z = sens_loc[node]\n    mlab.text3d(x, y, z, raw.ch_names[picks[node]], scale=0.005,\n                color=(0, 0, 0))\n\nview = (-88.7, 40.8, 0.76, np.array([-3.9e-4, -8.5e-3, -1e-2]))\nmlab.view(*view)"
      ], 
      "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"
      }
    }
  }
}