langouQuantBlind.ipynb

{
 "metadata": {
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.5"
  },
  "orig_nbformat": 4,
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3.9.5 64-bit ('mne': conda)"
  },
  "interpreter": {
   "hash": "e3970cf2a056ed71c6703e76baeca137440641ad1b48b7e7a1b0e3e35696df1b"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2,
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Initialize"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# Authors: Martin Luessi <mluessi@nmr.mgh.harvard.edu>\n",
    "#          Alexandre Gramfort <alexandre.gramfort@inria.fr>\n",
    "#          Nicolas P. Rougier (graph code borrowed from his matplotlib gallery)\n",
    "#\n",
    "# License: BSD (3-clause)\n",
    "\n",
    "import numpy as np\n",
    "import os.path as op\n",
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "from langouEEG import *\n",
    "from scipy.stats import *\n",
    "\n",
    "import mne\n",
    "import pickle\n",
    "from mne.datasets import sample\n",
    "from mne.minimum_norm import apply_inverse_epochs, read_inverse_operator\n",
    "from mne.connectivity import spectral_connectivity\n",
    "from mne.viz import circular_layout, plot_connectivity_circle\n",
    "import mne\n",
    "from mne.datasets import eegbci\n",
    "from mne.datasets import fetch_fsaverage\n",
    "from mne.datasets import sample\n",
    "from mne.minimum_norm import make_inverse_operator, apply_inverse\n",
    "from mne.minimum_norm import write_inverse_operator\n",
    "\n",
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D  # noqa\n",
    "import mne\n",
    "\n",
    "sample_data_folder = mne.datasets.sample.data_path()\n",
    "dataRoot = \"/data/home/viscent/Light\"\n",
    "# Download fsaverage files\n",
    "fs_dir = fetch_fsaverage(verbose=True)\n",
    "isMale = False\n",
    "isAll = True\n",
    "isBlind = True\n",
    "if not isAll:\n",
    "    result_dir = op.join(dataRoot,'result','male' if isMale else 'female')\n",
    "else:\n",
    "    result_dir = op.join(dataRoot,'result','all')\n",
    "if isBlind:\n",
    "    result_dir = op.join(result_dir,'Blind')\n",
    "subjects_dir = op.dirname(fs_dir)\n",
    "if not op.exists(result_dir):\n",
    "    os.mkdir(result_dir)\n",
    "# The files live in:\n",
    "subject = 'fsaverage'\n",
    "trans = 'fsaverage'  # MNE has a built-in fsaverage transformation\n",
    "src = op.join(fs_dir, 'bem', 'fsaverage-ico-5-src.fif')\n",
    "bem = op.join(fs_dir, 'bem', 'fsaverage-5120-5120-5120-bem-sol.fif')\n",
    "print(__doc__)"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "\n",
    "import os\n",
    "os.system('/usr/bin/Xvfb :99 -screen 0 1920x1080x24 &')\n",
    "os.environ['DISPLAY'] = ':99'\n",
    "\n",
    "import panel as pn\n",
    "pn.extension('vtk')"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Load data"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "epochs_4F = dict()\n",
    "epochs_RF = dict()\n",
    "epochs_4R = dict()\n",
    "epochs_RR = dict()\n",
    "for subject_name in range(21,24):\n",
    "    if not isAll:\n",
    "        if not (isMale ^ (subject_name in [7,8,11,17])):\n",
    "            continue\n",
    "    if subject_name<10:\n",
    "        subject_name='S0'+str(subject_name)\n",
    "    else:\n",
    "        subject_name='S'+str(subject_name)\n",
    "    with open(dataRoot+'/clean_data_av/'+subject_name+'_clean.lgeeg','rb') as f:\n",
    "        raw=pickle.load(f)\n",
    "    events, event_dict=extractEvents(raw)\n",
    "    picks = mne.pick_types(raw.info, meg=False, eeg=True, stim=False, eog=True,\n",
    "                        exclude='bads')\n",
    "    if isBlind:\n",
    "        epoch_RR,epoch_RF,epoch_4R,epoch_4F = extractEpochsBlind(raw,events,picks)\n",
    "    else:\n",
    "        epoch_RR,epoch_RF,epoch_4R,epoch_4F = extractEpochs(raw,events,picks)\n",
    "    epochs_4F[subject_name]=epoch_4F\n",
    "    epochs_RF[subject_name]=epoch_RF\n",
    "    epochs_RR[subject_name]=epoch_RR\n",
    "    epochs_4R[subject_name]=epoch_4R\n",
    "    if not op.exists(os.path.join(dataRoot,'fwd_solution.lgeeg')):\n",
    "        fwd = mne.make_forward_solution(raw.info, trans=trans, src=src,\n",
    "                                        bem=bem, eeg=True, mindist=5.0, n_jobs=1)\n",
    "        print(fwd)\n",
    "        mne.write_forward_solution(os.path.join(dataRoot,'fwd_solution.lgeeg'),fwd,overwrite=True)\n",
    "    else:\n",
    "        fwd = mne.read_forward_solution(os.path.join(dataRoot,'fwd_solution.lgeeg'))\n",
    "    data_path = sample.data_path()\n",
    "\n",
    "    picks = mne.pick_types(raw.info, meg=False, eeg=True, stim=False, eog=True,\n",
    "                        exclude='bads')\n",
    "\n",
    "\n",
    " \n",
    "\n",
    "\n"
   ],
   "outputs": [],
   "metadata": {
    "tags": []
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Source Estimation\n"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "cons_40 = []\n",
    "cons_rand = []\n",
    "act_40 = []\n",
    "act_40_paired = []\n",
    "act_rand = []\n",
    "act_rand_paired = []\n",
    "for subject_name,_ in epochs_4F.items():\n",
    "    # Compute inverse solution and for each epoch. By using \"return_generator=True\"\n",
    "    # stcs will be a generator object instead of a list.\n",
    "    epoch_4F = epochs_4F[subject_name]\n",
    "    epoch_RF = epochs_RF[subject_name]\n",
    "    result_dir = op.join(dataRoot,'result','male' if isMale else 'female',subject_name)\n",
    "    fname_stc_4F = op.join(result_dir,subject_name+'_stc_4F.lgeeg')\n",
    "    fname_stc_RF = op.join(result_dir,subject_name+'_stc_RF.lgeeg')\n",
    "    fname_tl_4F = op.join(result_dir,subject_name+'_tl_4F.lgeeg')\n",
    "    fname_tl_RF = op.join(result_dir,subject_name+'_tl_RF.lgeeg')\n",
    "    fname_inv = os.path.join(dataRoot,'inv_operators',subject_name+'_inv.lgeeg')\n",
    "    if not os.path.exists(result_dir):\n",
    "        os.mkdir(result_dir)\n",
    "\n",
    "\n",
    "\n",
    "    fname_inv_r = os.path.join(dataRoot,'inv_operators',subject_name+'_r_inv.lgeeg')\n",
    "    fname_cov_r = os.path.join(dataRoot,'noise_covariance',subject_name+'_r_cov.lgeeg')\n",
    "    fname_inv_4 = os.path.join(dataRoot,'inv_operators',subject_name+'_4_inv.lgeeg')\n",
    "    fname_cov_4 = os.path.join(dataRoot,'noise_covariance',subject_name+'_4_cov.lgeeg')\n",
    "\n",
    "\n",
    "\n",
    "    if not os.path.exists(fname_cov_4):  \n",
    "        noise_cov = mne.compute_covariance(\n",
    "            epochs_4R[subject_name], tmax=80., method=['shrunk', 'empirical'], rank=None, verbose=True)\n",
    "        mne.write_cov(fname_cov_4,noise_cov)\n",
    "    else:   \n",
    "    # Load data\n",
    "        noise_cov = mne.read_cov(fname_cov_4)\n",
    "    \n",
    "    if not os.path.exists(fname_inv_4):  \n",
    "        inverse_operator = make_inverse_operator(\n",
    "            raw.info, fwd, noise_cov, loose=0.2, depth=0.8)\n",
    "        write_inverse_operator(fname_inv_4,inverse_operator)\n",
    "    inverse_operator = read_inverse_operator(fname_inv_4)\n",
    "\n",
    "\n",
    "\n",
    "    snr = 1.0  # use lower SNR for single epochs\n",
    "    lambda2 = 1.0 / snr ** 2\n",
    "    method = \"dSPM\"  # use dSPM method (could also be MNE or sLORETA)\n",
    "    # Get labels for FreeSurfer 'aparc' cortical parcellation with 34 labels/hemi\n",
    "    labels = mne.read_labels_from_annot('fsaverage', parc='aparc',\n",
    "                                        subjects_dir=subjects_dir)[:68]\n",
    "    label_colors = [label.color for label in labels]\n",
    "    # Average the source estimates within each label using sign-flips to reduce\n",
    "    # signal cancellations, also here we return a generator\n",
    "    src = inverse_operator['src']\n",
    "    if op.exists(fname_tl_4F):\n",
    "        with open(fname_tl_4F,'rb') as f:\n",
    "            label_ts = pickle.load(f)\n",
    "    else:\n",
    "        stcs = apply_inverse_epochs(epoch_4F, inverse_operator, lambda2, method,\n",
    "                                    pick_ori=\"normal\", return_generator=False)\n",
    "        label_ts = mne.extract_label_time_course(\n",
    "            stcs, labels,  src,allow_empty=False, mode='mean_flip', return_generator=False)\n",
    "        with open(fname_tl_4F,'wb') as f:\n",
    "            pickle.dump(label_ts,f)\n",
    "    fmin = 8.\n",
    "    fmax = 13.\n",
    "    sfreq = raw.info['sfreq']  # the sampling frequency\n",
    "    con_methods = ['pli', 'wpli2_debiased', 'ciplv']\n",
    "    if not os.path.exists(op.join(result_dir,'cons')):\n",
    "        os.mkdir(op.join(result_dir,'cons'))\n",
    "    if not op.exists(op.join(result_dir,'cons','40_con.lgeeg')):\n",
    "        con, freqs, times, n_epochs, n_tapers = spectral_connectivity(\n",
    "            label_ts, method=con_methods, mode='multitaper', sfreq=sfreq, fmin=fmin,\n",
    "            fmax=fmax, faverage=True, mt_adaptive=True, n_jobs=1)\n",
    "        with open(op.join(result_dir,'cons','40_con.lgeeg'),'wb') as f:\n",
    "            pickle.dump(con,f)\n",
    "    else:\n",
    "        with open(op.join(result_dir,'cons','40_con.lgeeg'),'rb') as f:\n",
    "            con = pickle.load(f)\n",
    "    label_names = [label.name for label in labels]\n",
    "    con = np.squeeze(con)\n",
    "    cons_40.append(con)\n",
    "    pd_40=pd.DataFrame(con[1])\n",
    "    pd_40.columns = label_names\n",
    "    pd_40.index = label_names\n",
    "    pd_40.to_excel(op.join(result_dir,'40_conn.xlsx'))\n",
    "    con_40 = con.copy()\n",
    "    pd_40=pd.DataFrame(np.mean(label_ts,axis=2))\n",
    "    pd_40.columns = label_names\n",
    "    pd_40.to_excel(op.join(result_dir,'40_activation.xlsx'))\n",
    "    act_40.append(np.mean(np.array(label_ts),axis=2))\n",
    "    act_40_paired.append(np.mean(np.mean(np.array(label_ts),axis=2),axis=0))\n",
    "\n",
    "    if not os.path.exists(fname_cov_r):  \n",
    "        noise_cov = mne.compute_covariance(\n",
    "            epochs_RR[subject_name], tmax=80., method=['shrunk', 'empirical'], rank=None, verbose=True)\n",
    "        mne.write_cov(fname_cov_r,noise_cov)\n",
    "    else:   \n",
    "    # Load data\n",
    "        noise_cov = mne.read_cov(fname_cov_r)\n",
    "    \n",
    "    if not os.path.exists(fname_inv_r):  \n",
    "        inverse_operator = make_inverse_operator(\n",
    "            raw.info, fwd, noise_cov, loose=0.2, depth=0.8)\n",
    "        write_inverse_operator(fname_inv_r,inverse_operator)\n",
    "    inverse_operator = read_inverse_operator(fname_inv_r)\n",
    "\n",
    "\n",
    "    lambda2 = 1.0 / snr ** 2\n",
    "    method = \"dSPM\"  # use dSPM method (could also be MNE or sLORETA)\n",
    "\n",
    "\n",
    "\n",
    "    \n",
    "    if op.exists(fname_tl_RF):\n",
    "        with open(fname_tl_RF,'rb') as f:\n",
    "            label_ts = pickle.load(f)\n",
    "    else:\n",
    "        stcs = apply_inverse_epochs(epoch_RF, inverse_operator, lambda2, method,\n",
    "                                    pick_ori=\"normal\", return_generator=False)\n",
    "        label_ts = mne.extract_label_time_course(\n",
    "            stcs, labels,  src,allow_empty=False, mode='mean_flip', return_generator=False)\n",
    "        with open(fname_tl_RF,'wb') as f:\n",
    "            pickle.dump(label_ts,f)\n",
    "\n",
    "    \n",
    "    if not op.exists(op.join(result_dir,'cons','rand_con.lgeeg')):\n",
    "        con, freqs, times, n_epochs, n_tapers = spectral_connectivity(\n",
    "            label_ts, method=con_methods, mode='multitaper', sfreq=sfreq, fmin=fmin,\n",
    "            fmax=fmax, faverage=True, mt_adaptive=True, n_jobs=1)\n",
    "        with open(op.join(result_dir,'cons','rand_con.lgeeg'),'wb') as f:\n",
    "            pickle.dump(con,f)\n",
    "    else:\n",
    "        with open(op.join(result_dir,'cons','rand_con.lgeeg'),'rb') as f:\n",
    "            con = pickle.load(f)\n",
    "    # con is a 3D array, get the connectivity for the first (and only) freq. band\n",
    "    # for each method\n",
    "\n",
    "    con = np.squeeze(con)\n",
    "    cons_rand.append(con)\n",
    "    con_rand = con.copy()\n",
    "    pd_rand=pd.DataFrame(con[1])\n",
    "    pd_rand.columns = label_names\n",
    "    pd_rand.index = label_names\n",
    "    pd_rand.to_excel(op.join(result_dir,'rand_conn.xlsx'))\n",
    "    pd_delta=pd.DataFrame(con_40[1]-con_rand[1])\n",
    "    pd_delta.columns = label_names\n",
    "    pd_delta.index = label_names\n",
    "    pd_delta.to_excel(op.join(result_dir,'delta_conn.xlsx'))\n",
    "    pd_rand=pd.DataFrame(np.mean(label_ts,axis=2))\n",
    "    pd_rand.columns = label_names\n",
    "    pd_rand.to_excel(op.join(result_dir,'rand_activation.xlsx'))\n",
    "    act_rand.append(np.mean(np.array(label_ts),axis=2))\n",
    "    act_rand_paired.append(np.mean(np.mean(np.array(label_ts),axis=2),axis=0))\n",
    "cons_40 = np.array(cons_40)[:,1,...]\n",
    "cons_rand = np.array(cons_rand)[:,1,...]\n",
    "act_40 = np.vstack(act_40)\n",
    "act_rand = np.vstack(act_rand)\n",
    "act_40_paired = np.array(act_40_paired)\n",
    "act_rand_paired = np.array(act_rand_paired)"
   ],
   "outputs": [],
   "metadata": {
    "tags": []
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "act_40_norm = act_40.copy()\n",
    "for i in range(act_40.shape[0]):\n",
    "    act_40_norm[i] = (act_40[i]-act_40[i].mean())/act_40[i].std()\n",
    "act_40_lh = act_40_norm[:,::2]\n",
    "act_40_rh = act_40_norm[:,1::2]\n",
    "act_40_lh_mean = act_40_lh.mean(axis=1)\n",
    "act_40_rh_mean = act_40_rh.mean(axis=1)"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "plt.plot(act_40_lh_mean-act_40_rh_mean)"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# print('The p-value of L/R difference is: '+str(ttest_rel(act_40_lh_mean,act_40_rh_mean,axis=0).pvalue))\n",
    "pval = ttest_ind(act_40,act_rand,axis=0).pvalue\n",
    "plt.plot(np.where(pval<0.05,1,0))\n",
    "for i in range(68):\n",
    "    if pval[i] <0.05 :\n",
    "        print(label_names[i])\n",
    "# fig, ax = plt.subplots()\n",
    "# ax.pcolormesh(act_40_lh-act_40_rh)"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "stcs = []\n",
    "for subject_name,epoch in epochs_4F.items():\n",
    "    fname_inv_r = os.path.join(dataRoot,'inv_operators',subject_name+'_r_inv.lgeeg')\n",
    "    fname_cov_r = os.path.join(dataRoot,'noise_covariance',subject_name+'_r_cov.lgeeg')\n",
    "    fname_inv_4 = os.path.join(dataRoot,'inv_operators',subject_name+'_4_inv.lgeeg')\n",
    "    fname_cov_4 = os.path.join(dataRoot,'noise_covariance',subject_name+'_4_cov.lgeeg')\n",
    "\n",
    "\n",
    "\n",
    "    if not os.path.exists(fname_cov_4):  \n",
    "        noise_cov = mne.compute_covariance(\n",
    "            epochs_4R[subject_name], tmax=80., method=['shrunk', 'empirical'], rank=None, verbose=True)\n",
    "        mne.write_cov(fname_cov_4,noise_cov)\n",
    "    else:   \n",
    "    # Load data\n",
    "        noise_cov = mne.read_cov(fname_cov_4)\n",
    "    \n",
    "    if not os.path.exists(fname_inv_4):  \n",
    "        inverse_operator = make_inverse_operator(\n",
    "            raw.info, fwd, noise_cov, loose=0.2, depth=0.8)\n",
    "        write_inverse_operator(fname_inv_4,inverse_operator)\n",
    "    inverse_operator = read_inverse_operator(fname_inv_4)\n",
    "    stc = apply_inverse_epochs(epoch, inverse_operator, lambda2, method,\n",
    "                                pick_ori=\"normal\", return_generator=False)\n",
    "    stcs.append(stc)\n",
    "stc_mean = np.mean(stcs)\n",
    "    \n"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "pn.extension('vtk') \n",
    "vertno_max, time_max = stc_mean.get_peak(hemi='rh')\n",
    "\n",
    "subjects_dir = data_path + '/subjects'\n",
    "surfer_kwargs = dict(\n",
    "    hemi='rh', subjects_dir=subjects_dir,\n",
    "    clim=dict(kind='value', lims=[0, 2, 4]),backend= 'matplotlib', views='dor',\n",
    "    initial_time=time_max, time_unit='s', size=(800, 800), smoothing_steps=10)\n",
    "brain = stc_mean.plot(**surfer_kwargs)\n",
    "vertno_max, time_max = stc_mean.get_peak(hemi='lh')\n",
    "\n",
    "subjects_dir = data_path + '/subjects'\n",
    "surfer_kwargs = dict(\n",
    "    hemi='lh', subjects_dir=subjects_dir,\n",
    "    clim=dict(kind='value', lims=[0, 2, 4]),backend= 'matplotlib', views='dor',\n",
    "    initial_time=time_max, time_unit='s', size=(800, 800), smoothing_steps=10)\n",
    "brain = stc_mean.plot(**surfer_kwargs)\n",
    "\n",
    "\n"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "stcs = []\n",
    "for subject_name,epoch in epochs_RF.items():\n",
    "    fname_inv_r = os.path.join(dataRoot,'inv_operators',subject_name+'_r_inv.lgeeg')\n",
    "    fname_cov_r = os.path.join(dataRoot,'noise_covariance',subject_name+'_r_cov.lgeeg')\n",
    "    fname_inv_4 = os.path.join(dataRoot,'inv_operators',subject_name+'_4_inv.lgeeg')\n",
    "    fname_cov_4 = os.path.join(dataRoot,'noise_covariance',subject_name+'_4_cov.lgeeg')\n",
    "\n",
    "\n",
    "\n",
    "    if not os.path.exists(fname_cov_r):  \n",
    "        noise_cov = mne.compute_covariance(\n",
    "            epochs_RR[subject_name], tmax=80., method=['shrunk', 'empirical'], rank=None, verbose=True)\n",
    "        mne.write_cov(fname_cov_r,noise_cov)\n",
    "    else:   \n",
    "    # Load data\n",
    "        noise_cov = mne.read_cov(fname_cov_r)\n",
    "    \n",
    "    if not os.path.exists(fname_inv_r):  \n",
    "        inverse_operator = make_inverse_operator(\n",
    "            raw.info, fwd, noise_cov, loose=0.2, depth=0.8)\n",
    "        write_inverse_operator(fname_inv_r,inverse_operator)\n",
    "    inverse_operator = read_inverse_operator(fname_inv_r)\n",
    "    stc = apply_inverse_epochs(epoch, inverse_operator, lambda2, method,\n",
    "                                pick_ori=\"normal\", return_generator=False)\n",
    "    stcs.append(stc)\n",
    "stc_mean = np.mean(stcs)"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "pn.extension('vtk') \n",
    "vertno_max, time_max = stc_mean.get_peak(hemi='rh')\n",
    "\n",
    "subjects_dir = data_path + '/subjects'\n",
    "surfer_kwargs = dict(\n",
    "    hemi='rh', subjects_dir=subjects_dir,\n",
    "    clim=dict(kind='value', lims=[0, 2, 4]),backend= 'matplotlib', views='dor',\n",
    "    initial_time=time_max, time_unit='s', size=(800, 800), smoothing_steps=10)\n",
    "plt.figure(dpi=600)\n",
    "brain = stc_mean.plot(**surfer_kwargs)\n",
    "vertno_max, time_max = stc_mean.get_peak(hemi='lh')\n",
    "\n",
    "subjects_dir = data_path + '/subjects'\n",
    "surfer_kwargs = dict(\n",
    "    hemi='lh', subjects_dir=subjects_dir,\n",
    "    clim=dict(kind='value', lims=[0, 2, 4]),backend= 'matplotlib', views='dor',\n",
    "    initial_time=time_max, time_unit='s', size=(800, 800), smoothing_steps=10)\n",
    "brain = stc_mean.plot(**surfer_kwargs)\n"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "type(brain)"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "pval = 1-ttest_rel(cons_40,cons_rand,axis=0).pvalue\n",
    "pval_mat=np.where(pval>0.95,pval-0.95,0)\n",
    "fig, ax = plt.subplots()\n",
    "ax.pcolormesh(pval_mat)"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "pval = 1-ttest_rel(cons_40,cons_rand,axis=0).pvalue\n",
    "fig, ax = plt.subplots()\n",
    "pmin = 0.99\n",
    "ax.pcolormesh(np.where(pval>pmin,pval-pmin,0))\n",
    "fig, ax = plt.subplots()\n",
    "ax.pcolormesh(np.where(pval>0.999,pval-pmin,0))"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "result_dir = op.join(dataRoot, 'results','overall')\n",
    "if not op.exists(result_dir):\n",
    "    os.makedirs(result_dir)\n",
    "# First, we reorder the labels based on their location in the left hemi\n",
    "label_names = [label.name for label in labels]\n",
    "\n",
    "lh_labels = [name for name in label_names if name.endswith('lh')]\n",
    "\n",
    "# Get the y-location of the label\n",
    "label_ypos = list()\n",
    "for name in lh_labels:\n",
    "    idx = label_names.index(name)\n",
    "    ypos = np.mean(labels[idx].pos[:, 1])\n",
    "    label_ypos.append(ypos)\n",
    "\n",
    "# Reorder the labels based on their location\n",
    "lh_labels = [label for (yp, label) in sorted(zip(label_ypos, lh_labels))]\n",
    "\n",
    "# For the right hemi\n",
    "rh_labels = [label[:-2] + 'rh' for label in lh_labels]\n",
    "\n",
    "# Save the plot order and create a circular layout\n",
    "node_order = list()\n",
    "node_order.extend(lh_labels[::-1])  # reverse the order\n",
    "node_order.extend(rh_labels)\n",
    "# node_order = node_order[:69]\n",
    "\n",
    "node_angles = circular_layout(label_names, node_order, start_pos=90,\n",
    "                              group_boundaries=[0, len(label_names) / 2])\n",
    "\n",
    "# Plot the graph using node colors from the FreeSurfer parcellation. We only\n",
    "# show the 300 strongest connections.\n",
    "fig,ax=plot_connectivity_circle(cons_40.mean(axis=0), label_names, n_lines=300,\n",
    "                         node_angles=node_angles, node_colors=label_colors, \n",
    "                         title='All-to-All Connectivity 40 Hz '\n",
    "                               'Condition (PLI)')\n",
    "fig.savefig(op.join(result_dir,'40-connectivity.png'),facecolor='black')\n",
    "fig,ax=plot_connectivity_circle(cons_rand.mean(axis=0), label_names, n_lines=300,\n",
    "                         node_angles=node_angles, node_colors=label_colors, \n",
    "                         title='All-to-All Connectivity Random '\n",
    "                               'Condition (PLI)')\n",
    "fig.savefig(op.join(result_dir,'rand-connectivity.png'),facecolor='black')\n",
    "dcon_res = cons_40.mean(axis=0) - cons_rand.mean(axis=0)\n",
    "fig,ax=plot_connectivity_circle(np.where(dcon_res<0,0,dcon_res), label_names, n_lines=300,\n",
    "                         node_angles=node_angles, node_colors=label_colors, colormap='hot', \n",
    "                         title='All-to-All Connectivity Upregulated '\n",
    "                               '(PLI)')\n",
    "fig.savefig(op.join(result_dir,'up-connectivity.png'),facecolor='black')\n",
    "fig,ax=plot_connectivity_circle(np.where(dcon_res>0,0,-1*dcon_res), label_names, n_lines=300,\n",
    "                         node_angles=node_angles, node_colors=label_colors, colormap='gray',\n",
    "                         title='All-to-All Connectivity Downregulated '\n",
    "                               '(PLI)')\n",
    "fig.savefig(op.join(result_dir,'down-connectivity.png'),facecolor='black')\n"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "pval = 1-ttest_rel(cons_40,cons_rand,axis=0).pvalue\n",
    "pval_mat = np.where(pval>0.95,1-pval,0)\n",
    "fig,ax=plot_connectivity_circle(pval_mat, label_names, n_lines=300,\n",
    "                         node_angles=node_angles, node_colors=label_colors, colormap='hot', \n",
    "                         title='Robustly Altered Connectivity'\n",
    "                               '(p<0.05)')\n",
    "fig.savefig(op.join(result_dir,'t-test.png'),facecolor='black')\n",
    "\n",
    "pval_mat = np.where(pval>0.99,1-pval,0)\n",
    "fig,ax=plot_connectivity_circle(pval_mat, label_names, n_lines=300,\n",
    "                         node_angles=node_angles, node_colors=label_colors, colormap='hot', \n",
    "                         title='Robustly Altered Connectivity'\n",
    "                               '(p<0.01)')\n",
    "fig.savefig(op.join(result_dir,'t-test.png'),facecolor='black')\n",
    "\n",
    "pval_mat = np.where(pval>0.999,1-pval,0)\n",
    "fig,ax=plot_connectivity_circle(pval_mat, label_names, n_lines=300,\n",
    "                         node_angles=node_angles, node_colors=label_colors, colormap='hot', \n",
    "                         title='Robustly Altered Connectivity'\n",
    "                               '(p<0.001)')\n",
    "fig.savefig(op.join(result_dir,'t-test.png'),facecolor='black')"
   ],
   "outputs": [],
   "metadata": {}
  }
 ]
}