{
  "nbformat_minor": 0, 
  "nbformat": 4, 
  "cells": [
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "%matplotlib inline"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "\n\n# 2 samples permutation test on source data with spatio-temporal clustering\n\n\nTests if the source space data are significantly different between\n2 groups of subjects (simulated here using one subject's data).\nThe multiple comparisons problem is addressed with a cluster-level\npermutation test across space and time.\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "# Authors: Alexandre Gramfort <alexandre.gramfort@telecom-paristech.fr>\n#          Eric Larson <larson.eric.d@gmail.com>\n# License: BSD (3-clause)\n\nimport os.path as op\n\nimport numpy as np\nfrom scipy import stats as stats\n\nimport mne\nfrom mne import spatial_tris_connectivity, grade_to_tris\nfrom mne.stats import spatio_temporal_cluster_test, summarize_clusters_stc\nfrom mne.datasets import sample\n\nprint(__doc__)"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Set parameters\n--------------\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "data_path = sample.data_path()\nstc_fname = data_path + '/MEG/sample/sample_audvis-meg-lh.stc'\nsubjects_dir = data_path + '/subjects'\n\n# Load stc to in common cortical space (fsaverage)\nstc = mne.read_source_estimate(stc_fname)\nstc.resample(50, npad='auto')\n\nstc = mne.morph_data('sample', 'fsaverage', stc, grade=5, smooth=20,\n                     subjects_dir=subjects_dir)\nn_vertices_fsave, n_times = stc.data.shape\ntstep = stc.tstep\n\nn_subjects1, n_subjects2 = 7, 9\nprint('Simulating data for %d and %d subjects.' % (n_subjects1, n_subjects2))\n\n#    Let's make sure our results replicate, so set the seed.\nnp.random.seed(0)\nX1 = np.random.randn(n_vertices_fsave, n_times, n_subjects1) * 10\nX2 = np.random.randn(n_vertices_fsave, n_times, n_subjects2) * 10\nX1[:, :, :] += stc.data[:, :, np.newaxis]\n# make the activity bigger for the second set of subjects\nX2[:, :, :] += 3 * stc.data[:, :, np.newaxis]\n\n#    We want to compare the overall activity levels for each subject\nX1 = np.abs(X1)  # only magnitude\nX2 = np.abs(X2)  # only magnitude"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Compute statistic\n-----------------\n\nTo use an algorithm optimized for spatio-temporal clustering, we\njust pass the spatial connectivity matrix (instead of spatio-temporal)\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "print('Computing connectivity.')\nconnectivity = spatial_tris_connectivity(grade_to_tris(5))\n\n#    Note that X needs to be a list of multi-dimensional array of shape\n#    samples (subjects_k) x time x space, so we permute dimensions\nX1 = np.transpose(X1, [2, 1, 0])\nX2 = np.transpose(X2, [2, 1, 0])\nX = [X1, X2]\n\n#    Now let's actually do the clustering. This can take a long time...\n#    Here we set the threshold quite high to reduce computation.\np_threshold = 0.0001\nf_threshold = stats.distributions.f.ppf(1. - p_threshold / 2.,\n                                        n_subjects1 - 1, n_subjects2 - 1)\nprint('Clustering.')\nT_obs, clusters, cluster_p_values, H0 = clu =\\\n    spatio_temporal_cluster_test(X, connectivity=connectivity, n_jobs=1,\n                                 threshold=f_threshold)\n#    Now select the clusters that are sig. at p < 0.05 (note that this value\n#    is multiple-comparisons corrected).\ngood_cluster_inds = np.where(cluster_p_values < 0.05)[0]"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Visualize the clusters\n----------------------\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "print('Visualizing clusters.')\n\n#    Now let's build a convenient representation of each cluster, where each\n#    cluster becomes a \"time point\" in the SourceEstimate\nfsave_vertices = [np.arange(10242), np.arange(10242)]\nstc_all_cluster_vis = summarize_clusters_stc(clu, tstep=tstep,\n                                             vertices=fsave_vertices,\n                                             subject='fsaverage')\n\n#    Let's actually plot the first \"time point\" in the SourceEstimate, which\n#    shows all the clusters, weighted by duration\nsubjects_dir = op.join(data_path, 'subjects')\n# blue blobs are for condition A != condition B\nbrain = stc_all_cluster_vis.plot('fsaverage', hemi='both', colormap='mne',\n                                 views='lateral', subjects_dir=subjects_dir,\n                                 time_label='Duration significant (ms)')\nbrain.save_image('clusters.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"
      }
    }
  }
}