main.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## This notebook contains the necessary code to replicate results for: "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# *The connectome spectrum as a canonical basis for a sparse representation of fast brain activity*\n",
    "Article authored by: Joan Rué-Queralt, Katharina Glomb; David Pascucci; Sebastien Tourbier; Margherita Carboni; Serge Vulliémoz; Gijs Plomp; and Patric Hagmann  \n",
    "Notebook created by: Joan Rué-Queralt, Joan.Rue-Queralt@chuv.ch  \n",
    "Last update 22-September-2021"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Import necessary packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import tqdm\n",
    "import pandas as pd\n",
    "import scipy.io\n",
    "from scipy import stats\n",
    "from scipy.spatial.distance import pdist,squareform\n",
    "from scipy.stats import ranksums\n",
    "\n",
    "import matplotlib\n",
    "from matplotlib import cm\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "from sklearn.decomposition import PCA, FastICA\n",
    "import pygsp\n",
    "\n",
    "from utils.data import TsGenerator, Connectome\n",
    "from utils.data import loader_sc, loader_sc_surrogates,\\\n",
    "                       loader_sc_surrogates_geometry_preserving\n",
    "from utils.data import loader_sc_non_norm, loader_sc_surrogates_non_norm, \\\n",
    "                       loader_sc_surrogates_geometry_preserving_non_norm\n",
    "from utils.compactness import compress, compress_error, compactness_scores, \\\n",
    "                              compactness_scores_data_driven, \\\n",
    "                              compactness_scores_geometry_preserving, \\\n",
    "                              compactness_dynamics\n",
    "from utils.visualization import plot_shaded, plot_shaded_norm, plot_surface_ld, \\\n",
    "                                plot_surface_hd, fancy_fig, my_box_plot\n",
    "from utils.stats import perm_test, reject_outliers,conditional_probability\n",
    "from utils.graph import sdi, smoothness"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data directories\n",
    "\n",
    "datadir = './data'\n",
    "bidsdir = os.path.join(datadir,'bids_dir')\n",
    "scale = 3\n",
    "sc_rand_dir = os.path.join(datadir,'SC_surrogates','degree_preserving','SC_num',f'scale_{scale}')\n",
    "sc_rand_dir_geometry =  os.path.join(datadir,'SC_surrogates','geometry_preserving','SC_num',f'scale_{scale}')\n",
    "\n",
    "# define data properties\n",
    "subject_list_faces = [1,2,3,4,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20]\n",
    "subject_list_motion = [1,2,3,4,6,7,8,9,10,11,12,13,14,16,17,18,19,20]\n",
    "nsub_faces = len(subject_list_faces)\n",
    "nsub_motion = len(subject_list_motion)\n",
    "tvec = np.arange(-1500,1000,4) / 1000\n",
    "\n",
    "# define time of interest (from -100 ms to 600 ms after stimulus presentation)\n",
    "tvec_analysis = np.where((tvec>=-0.1)&(tvec<=0.6))[0] \n",
    "tvec_pre = np.where((tvec[tvec_analysis]<0))[0]\n",
    "\n",
    "# create data loader objects\n",
    "ts_gen_faces = TsGenerator(subject_list_faces,scale,'Faces',datadir,tvec_analysis,tvec_pre)\n",
    "ts_gen_motion = TsGenerator(subject_list_motion,scale,'Motion',datadir,tvec_analysis,tvec_pre)\n",
    "sc = loader_sc(scale,datadir)\n",
    "sc_surro = loader_sc_surrogates(sc_rand_dir)\n",
    "sc_surro_geometry = loader_sc_surrogates_geometry_preserving(sc_rand_dir_geometry,scale,datadir)\n",
    "\n",
    "# percentiles for compactness analyisis\n",
    "percentiles = np.arange(0,101)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load visual evoked signals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some notation for variables: \n",
    "- ts: time-series\n",
    "- roi: signal defined at the level of region of interest\n",
    "- ch: signal defined at the level of connectome harmonics\n",
    "- face: signal corresponding to the faces stimuli\n",
    "- scra: signal corresponding to the scrambled faces stimuli"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# allocate mem. for time-series\n",
    "face_roi_ts = np.zeros((sc.N,len(tvec_analysis),nsub_faces))\n",
    "face_ch_ts = np.zeros((sc.N,len(tvec_analysis),nsub_faces))\n",
    "scra_roi_ts = np.zeros((sc.N,len(tvec_analysis),nsub_faces))\n",
    "scra_ch_ts = np.zeros((sc.N,len(tvec_analysis),nsub_faces))\n",
    "\n",
    "# allocate mem. for power of time-series\n",
    "power_in_time = np.zeros((nsub_faces,len(tvec_analysis)))\n",
    "power_in_time_motion = np.zeros((nsub_motion,len(tvec_analysis)))\n",
    "\n",
    "# Compute evoked patterns and the power in time for the FACES dataset.\n",
    "\n",
    "# For each subject:\n",
    "for s in tqdm.tqdm(range(nsub_faces)):    \n",
    "    # load time-series for subject [rois x time x trials]\n",
    "    epochs,cond = ts_gen_faces.loader_ts(s)\n",
    "    # group trials by condition\n",
    "    epochs_f = epochs[:,:,cond==1] # faces \n",
    "    epochs_s = epochs[:,:,cond==0] # scrambled\n",
    "    # perform baseline (pre-stimulus) correction\n",
    "    epochs_f = epochs_f - epochs_f[:,tvec_pre].mean(1,keepdims=True)\n",
    "    epochs_s = epochs_s - epochs_s[:,tvec_pre].mean(1,keepdims=True)\n",
    "    # perform graph Fourier transform\n",
    "    epochs_graph_f = sc.gft(epochs_f)\n",
    "    epochs_graph_s = sc.gft(epochs_s)\n",
    "   \n",
    "    # compute evoked response\n",
    "    face_roi_ts[:,:,s] = epochs_f.mean(-1)    \n",
    "    face_ch_ts[:,:,s] = epochs_graph_f.mean(-1)\n",
    "    scra_roi_ts[:,:,s] = epochs_s.mean(-1)    \n",
    "    scra_ch_ts[:,:,s] = epochs_graph_s.mean(-1)\n",
    "    \n",
    "    # compute power (l2 norm across ROI or harmonics)\n",
    "    power_in_time[s] = np.linalg.norm(face_roi_ts[:,:,s],ord=2,axis=0)\n",
    "\n",
    "# Normalize power for each subject\n",
    "power_in_time /= np.sum(power_in_time,axis=1,keepdims=True)\n",
    "\n",
    "\n",
    "\n",
    "# Do the same for the MOTION dataset.\n",
    "motion_roi_ts = np.zeros((sc.N,len(tvec_analysis),nsub_motion))\n",
    "motion_ch_ts = np.zeros((sc.N,len(tvec_analysis),nsub_motion))\n",
    "random_roi_ts = np.zeros((sc.N,len(tvec_analysis),nsub_motion))\n",
    "random_ch_ts = np.zeros((sc.N,len(tvec_analysis),nsub_motion))\n",
    "\n",
    "# For each subject:\n",
    "for s in tqdm.tqdm(range(nsub_motion)):    \n",
    "    # Load time-series for subject [ROIs x Time x Trials]\n",
    "    epochs,cond = ts_gen_motion.loader_ts(s)\n",
    "    # Select only motion trials\n",
    "    epochs_m = epochs[:,:,cond==1] \n",
    "    epochs_r = epochs[:,:,cond==0] \n",
    "    # Baseline correction\n",
    "    epochs_m = epochs_m - epochs_m[:,tvec_pre].mean(1,keepdims=True)\n",
    "    epochs_r = epochs_r - epochs_r[:,tvec_pre].mean(1,keepdims=True)\n",
    "    # Graph Fourier transform\n",
    "    epochs_graph_m = sc.gft(epochs_m)\n",
    "    epochs_graph_r = sc.gft(epochs_r)\n",
    "    motion_roi_ts[:,:,s] = epochs_m.mean(-1)    \n",
    "    motion_ch_ts[:,:,s] = epochs_graph_m.mean(-1)\n",
    "    random_roi_ts[:,:,s] = epochs_r.mean(-1)    \n",
    "    random_ch_ts[:,:,s] = epochs_graph_r.mean(-1)\n",
    "    # compute power (l2 norm across ROI or harmonics)\n",
    "    power_in_time_motion[s] = np.linalg.norm(motion_roi_ts[:,:,s],ord=2,axis=0)\n",
    "# Normalize power for each subject\n",
    "\n",
    "\n",
    "power_in_time_motion /= np.sum(power_in_time_motion,axis=1,keepdims=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plots for Figure 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "font = {'size'   : 16}\n",
    "matplotlib.rc('font', **font)\n",
    "\n",
    "cmap = matplotlib.cm.get_cmap('Spectral')\n",
    "color1 = cmap(.05) # ROI \n",
    "color2 = cmap(.95) # Connectome Harmonics\n",
    "\n",
    "# Plot time-series of each ROI\n",
    "fig,axs = plt.subplots(1,2,figsize=(8 ,2))\n",
    "ax = axs[0]\n",
    "ax.pcolormesh(tvec[tvec_analysis],np.arange(sc.N),face_roi_ts.mean(-1),cmap='Greys_r',shading='auto')\n",
    "ax.axvline(0,ls='--',c='k')\n",
    "ax.set_xlabel('Time [s]')\n",
    "ax.set_ylabel('Brain area')\n",
    "\n",
    "# Plot time-series of each Connectome Harmonics\n",
    "ax = axs[1]\n",
    "ax.pcolormesh(  tvec[tvec_analysis],sc.e,face_ch_ts.mean(-1),cmap='Greys_r',shading='auto')\n",
    "ax.axvline(0,ls='--',c='k')\n",
    "ax.set_xlabel('Time [s]')\n",
    "ax.set_ylabel('Connectome \\n spectrum ($\\lambda$)')\n",
    "fig.tight_layout()\n",
    "\n",
    "# Plot the Connectome graph spectrum\n",
    "fig,ax = plt.subplots(figsize=(1.5,1))\n",
    "ax.plot(np.arange(sc.N)+1,sc.e,lw=2,c='k')\n",
    "ax.set_xlabel('Eigenvect. ordinal')\n",
    "ax.set_ylabel('Eigenval. ($\\lambda$)')\n",
    "fancy_fig(ax) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analyses for Figure 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define paths where results will be stored / loaded from"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "load_data = True\n",
    "compactness_path = os.path.join(datadir,'results_SC_num','compactness')\n",
    "dynamics_path = os.path.join(datadir,'results_SC_num','dynamics')\n",
    "corr_path = os.path.join(datadir,'results_SC_num','corr')\n",
    "topology_path = os.path.join(datadir,'results_SC_num','topology')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute signal compression in harmonics from geometry preserving surrogate graphs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if not load_data:\n",
    "    sc_surro_wwp, sc_surro_wsp, sc_surro_wssp = sc_surro_geometry\n",
    "\n",
    "    compression_error_wwp, correlation_wwp, compression_error_wsp, \\\n",
    "    correlation_wsp, compression_error_wssp, correlation_wssp = \\\n",
    "    compactness_scores_geometry_preserving(percentiles,\n",
    "                                          ts_gen_faces,\n",
    "                                          sc_surro_wwp, \n",
    "                                          [], \n",
    "                                          [],\n",
    "                                          options = ['wwp'])\n",
    "\n",
    "    np.save(os.path.join(compactness_path,'compression_wwp.npy'),compression_error_wwp)\n",
    "    np.save(os.path.join(compactness_path,'correlation_wwp.npy'),correlation_wwp)\n",
    "else:\n",
    "    compression_error_wwp = np.load(os.path.join(compactness_path,'compression_wwp.npy'))\n",
    "    correlation_wwp = np.load(os.path.join(compactness_path,'correlation_wwp.npy'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create euclidean graphs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "roifname = os.path.join(datadir, 'Lausanne2008_Yeo7RSNs.xlsx')\n",
    "roidata = pd.read_excel(roifname, sheet_name=f'SCALE {scale}')\n",
    "cort = np.where(roidata['Structure'] == 'cort')[0]\n",
    "\n",
    "euc = []\n",
    "for s,sub in enumerate(ts_gen_faces.subject_list):\n",
    "    subject = f'sub-{str(sub).zfill(2)}'\n",
    "    pos = scipy.io.loadmat(os.path.join(datadir,\n",
    "                                     'derivatives',\n",
    "                                     'cmp-v3.0.0-beta-RC1',\n",
    "                                     subject,\n",
    "                                     'connectivity',\n",
    "                                     f'{subject}_label-L2008_desc-scale{scale}_conndata-snetwork_connectivity.mat'))['nodes'][0][0][0][cort]\n",
    "\n",
    "    inv_euc = (1/squareform(pdist(pos)))\n",
    "    inv_euc = (inv_euc + inv_euc.T)*0.5\n",
    "    inv_euc[np.diag_indices(len(pos))]=0\n",
    "    euc.append( Connectome(inv_euc))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute signal compression in harmonics from euclidean graphs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils.compactness import compactness_scores_euclidean\n",
    "\n",
    "if not load_data:\n",
    "    compression_error_euc, correlation_euc = \\\n",
    "    compactness_scores_euclidean(percentiles,ts_gen_faces,euc)\n",
    "    np.save(os.path.join(compactness_path,'compression_euc.npy'),compression_error_euc)\n",
    "    np.save(os.path.join(compactness_path,'correlation_euc.npy'),correlation_euc)\n",
    "else:\n",
    "    compression_error_euc = np.load(os.path.join(compactness_path,'compression_euc.npy'))\n",
    "    correlation_euc = np.load(os.path.join(compactness_path,'correlation_euc.npy'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute signal compression in roi signal, connectome harmonics signal and harmonics from degree-preserving surrogate graphs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not load_data:\n",
    "    compression_error_roi,\\\n",
    "    correlation_roi,\\\n",
    "    compression_error_ch,\\\n",
    "    correlation_ch,\\\n",
    "    compression_error_rand_ch,\\\n",
    "    correlation_rand_ch = compactness_scores(percentiles,ts_gen_faces,sc,sc_surro)\n",
    "\n",
    "    np.save(os.path.join(compactness_path,'compression_ROI.npy'),compression_error_roi)\n",
    "    np.save(os.path.join(compactness_path,'correlation_ROI.npy'),correlation_roi)\n",
    "    np.save(os.path.join(compactness_path,'compression_GFT.npy'),compression_error_ch)\n",
    "    np.save(os.path.join(compactness_path,'correlation_GFT.npy'),correlation_ch)\n",
    "    np.save(os.path.join(compactness_path,'compression_RND.npy'),compression_error_rand_ch)\n",
    "    np.save(os.path.join(compactness_path,'correlation_RND.npy'),correlation_rand_ch)\n",
    "\n",
    "else:\n",
    "    \n",
    "    compression_error_roi = np.load(os.path.join(compactness_path,'compression_ROI.npy'))\n",
    "    compression_error_ch = np.load(os.path.join(compactness_path,'compression_GFT.npy'))\n",
    "    correlation_roi = np.load(os.path.join(compactness_path,'correlation_ROI.npy'))\n",
    "    correlation_ch = np.load(os.path.join(compactness_path,'correlation_GFT.npy'))\n",
    "    compression_error_rand_ch = np.load(os.path.join(compactness_path,'compression_RND.npy'))\n",
    "    correlation_rand_ch = np.load(os.path.join(compactness_path,'correlation_RND.npy'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Print explained variance for compressed signals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "print(f'Explained variance for compressed in ROI space: {(correlation_roi[:,95]**2).mean()}')\n",
    "print(f'Explained variance for compressed in ROI space: {(correlation_ch[:,95]**2).mean()}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute signal compression for data-driven decomposition methods (PCA and ICA)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not load_data:\n",
    "    compression_error_PCA, \\\n",
    "    correlation_PCA, \\\n",
    "    compression_error_ICA,\\\n",
    "    correlation_ICA = compactness_scores_data_driven(percentiles,ts_gen_faces)\n",
    "\n",
    "    np.save(os.path.join(compactness_path,'compression_PCA.npy'),compression_error_PCA)\n",
    "    np.save(os.path.join(compactness_path,'compression_ICA.npy'),compression_error_ICA)\n",
    "    np.save(os.path.join(compactness_path,'correlation_PCA.npy'),correlation_PCA)\n",
    "    np.save(os.path.join(compactness_path,'correlation_ICA.npy'),correlation_ICA)\n",
    "else:        \n",
    "\n",
    "    compression_error_PCA = np.load(os.path.join(compactness_path,'compression_PCA.npy'))\n",
    "    compression_error_ICA = np.load(os.path.join(compactness_path,'compression_ICA.npy'))\n",
    "    correlation_PCA = np.load(os.path.join(compactness_path,'correlation_PCA.npy'))\n",
    "    correlation_ICA = np.load(os.path.join(compactness_path,'correlation_ICA.npy'))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare statistically the compactness computed different coordinate spaces"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "data_comp = [np.mean(1-compression_error_roi,-1),np.mean(1-compression_error_ch,-1),np.mean(1-compression_error_rand_ch,-1).flatten(),np.mean(1-compression_error_wwp,-1).flatten(), np.mean(1-compression_error_PCA,-1), np.mean(1-compression_error_ICA,-1)]\n",
    "data_cor = [np.mean(correlation_roi,-1),np.mean(correlation_ch,-1),np.mean(correlation_rand_ch,-1).flatten(),np.mean(correlation_wwp,-1).flatten(), np.mean(correlation_PCA,-1), np.mean(correlation_ICA,-1)]\n",
    "\n",
    "\n",
    "pvalues_comp = np.zeros((len(data_comp),len(data_comp)))*np.nan\n",
    "pvalues_cor = np.zeros((len(data_comp),len(data_comp)))*np.nan\n",
    "\n",
    "zvalues_comp = np.zeros((len(data_comp),len(data_comp)))*np.nan\n",
    "zvalues_cor = np.zeros((len(data_comp),len(data_comp)))*np.nan\n",
    "\n",
    "ncomp = (len(data_comp)*(len(data_comp)-1))/2\n",
    "for i in range(len(data_comp)):\n",
    "    for j in range(i+1,len(data_comp)):\n",
    "        zvalues_comp[i,j],pvalues_comp[i,j] = ranksums(data_comp[i],data_comp[j])\n",
    "        zvalues_cor[i,j],pvalues_cor[i,j] = ranksums(data_cor[i],data_cor[j])        \n",
    "        \n",
    "pvalues_comp *= ncomp\n",
    "pvalues_comp[pvalues_comp>1] = 1\n",
    "pvalues_cor *= ncomp\n",
    "pvalues_cor[pvalues_cor>1] = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute compactness dynamics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not load_data:\n",
    "    compression_time_error_roi, \\\n",
    "    compression_time_error_ch = compactness_dynamics(ts_gen_faces, sc, sc_surro, percentiles)\n",
    "    \n",
    "    np.save(os.path.join(dynamics_path,'compression_ROI.npy'),compression_time_error_roi)\n",
    "    np.save(os.path.join(dynamics_path,'compression_GFT.npy'),compression_time_error_ch)\n",
    "    \n",
    "else:\n",
    "\n",
    "    compression_time_error_roi = np.load(os.path.join(dynamics_path,'compression_ROI.npy')) \n",
    "    compression_time_error_ch = np.load(os.path.join(dynamics_path,'compression_GFT.npy')) \n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare statistically compactness dynamics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Set a seed for the pseudo-random generator to replicate results\n",
    "np.random.seed(0)\n",
    "# Only compare withini the time-window of interest to reduce number of comparisons\n",
    "mask = np.array([True]*len(tvec_analysis))\n",
    "p_values,p_values_corrected,p_values_sig = perm_test(compression_time_error_roi,compression_time_error_ch,mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f'The time-points where compression is significantly different are: {tvec[tvec_analysis][p_values_sig]} seconds after stimulus presentation')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute the conditional probability  --> P( compactness | power )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "compactness_erp_roi_t = 1-compression_time_error_roi\n",
    "compactness_erp_ch_t = 1-compression_time_error_ch\n",
    "\n",
    "prob_cond1, prob_cond2, power, compactness = conditional_probability(compactness_erp_roi_t,compactness_erp_ch_t, power_in_time, tvec_analysis>=0)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Final plots Figure 2 (& Fig.S1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "font = {'size'   : 14}\n",
    "matplotlib.rc('font', **font)\n",
    "\n",
    "cmap = matplotlib.cm.get_cmap('Spectral')\n",
    "color1 = cmap(.05) # ROI \n",
    "color2 = cmap( .95) # Connectome Harmonics\n",
    "\n",
    "cmap = matplotlib.cm.get_cmap('Oranges')\n",
    "color3 = cmap(.95) # PCA \n",
    "cmap = matplotlib.cm.get_cmap('PiYG')\n",
    "color4 = cmap( .95) # ICA\n",
    "cmap = matplotlib.cm.get_cmap('cividis')\n",
    "color5 = cmap(.35) # Surrogate Harmonics\n",
    "\n",
    "colors =[color1,color2,color5,color5,color3,color4]\n",
    "labels =['$r$','$\\lambda$','$\\lambda_{dsur}$','$\\lambda_{gsur}$','PCA','ICA']\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,ax = plt.subplots(figsize=(4.2,3.15))\n",
    "plot_shaded_norm(percentiles,1-compression_error_roi,0,color1,'$r$',ax)\n",
    "plot_shaded_norm(percentiles,1-compression_error_ch,0,color2,'$\\lambda$',ax)\n",
    "plot_shaded_norm(percentiles,correlation_roi,0,color1,None,ax,ls='--')\n",
    "plot_shaded_norm(percentiles,correlation_ch,0,color2,None,ax,ls='--')\n",
    "plt.plot([-2,-1],[-2,-1],'--k',label='Correlation')\n",
    "plt.plot([-2,-1],[-2,-1],'-k',label='1 - error')\n",
    "ax.set_xlim(0,100)\n",
    "ax.set_ylim(0,1)\n",
    "ax.set_xlabel('Percentile')\n",
    "ax.legend(frameon=False)\n",
    "fancy_fig(ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,axs = plt.subplots(2,1,figsize=(5*0.8,4*0.8))\n",
    "ax = axs[0]\n",
    "ax = my_box_plot(data_cor,ax,colors,labels);\n",
    "ax.set_ylabel('Correlation')\n",
    "ax.set_ylim(0.7,1)\n",
    "ax.set_yticks((0.8,1))\n",
    "plt.setp(ax.get_xticklabels(), size=13)\n",
    "fancy_fig(ax)\n",
    "ax = axs[1]\n",
    "ax = my_box_plot(data_comp,ax,colors,labels);\n",
    "fancy_fig(ax)\n",
    "ax.set_ylim(0.3,0.9)\n",
    "ax.set_yticks((0.5,0.9))\n",
    "ax.set_ylabel('1 - error')\n",
    "ax.set_xlabel('Coordinate system')\n",
    "plt.setp(ax.get_xticklabels(), size=13)\n",
    "fig.tight_layout()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1,2,figsize=(15,6))\n",
    "logthresh = 1\n",
    "ax =axs[0]\n",
    "maxval = np.nanmax(abs(np.concatenate((zvalues_cor,zvalues_comp))))\n",
    "cmap = matplotlib.cm.get_cmap('PRGn').copy()\n",
    "cmap.set_bad(color='gray')\n",
    "im = ax.imshow(zvalues_cor,cmap=cmap,norm=matplotlib.colors.SymLogNorm(10**-logthresh,vmax=maxval,vmin=-maxval))\n",
    "ax.set_xticks(np.arange(len(data_comp)))\n",
    "ax.set_yticks(np.arange(len(data_comp)))\n",
    "ax.set_xticklabels(labels)\n",
    "ax.set_yticklabels(labels)\n",
    "plt.setp(ax.get_xticklabels(), rotation=45, ha=\"right\",\n",
    "         rotation_mode=\"anchor\")\n",
    "for i in range(len(labels)):\n",
    "    for j in range(i+1,len(labels)):\n",
    "        text = ax.text(j, i, \"{:2.0e}\".format(pvalues_cor[i, j]),\n",
    "                       {'size':'medium'},ha=\"center\", va=\"center\", \n",
    "                       color=\"w\")\n",
    "ax.tick_params(width=3, length=8)\n",
    "ax.spines['top'].set_linewidth(3)\n",
    "ax.spines['bottom'].set_linewidth(3)\n",
    "ax.spines['left'].set_linewidth(3)\n",
    "ax.spines['right'].set_linewidth(3)\n",
    "cb=plt.colorbar(im,ax = axs[0])\n",
    "\n",
    "ax =axs[1]\n",
    "im = ax.imshow(zvalues_comp,cmap=cmap,norm=matplotlib.colors.SymLogNorm(10**-logthresh,vmax=maxval,vmin=-maxval))\n",
    "ax.set_xticks(np.arange(len(data_comp)))\n",
    "ax.set_yticks(np.arange(len(data_comp)))\n",
    "ax.set_xticklabels(labels)\n",
    "ax.set_yticklabels(labels)\n",
    "plt.setp(ax.get_xticklabels(), rotation=45, ha=\"right\",\n",
    "         rotation_mode=\"anchor\")\n",
    "for i in range(len(labels)):\n",
    "    for j in range(i+1,len(labels)):\n",
    "        text = ax.text(j, i, \"{:2.0e}\".format(pvalues_comp[i, j]),\n",
    "                       {'size':'medium'},ha=\"center\", va=\"center\", \n",
    "                       color=\"w\")\n",
    "ax.tick_params(width=3, length=8)\n",
    "ax.spines['top'].set_linewidth(3)\n",
    "ax.spines['bottom'].set_linewidth(3)\n",
    "ax.spines['left'].set_linewidth(3)\n",
    "ax.spines['right'].set_linewidth(3)\n",
    "cb=plt.colorbar(im,ax = axs[1])\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "    \n",
    "fig,ax = plt.subplots(figsize=(4.2,3.15))\n",
    "ax.axvline(0,ls='-',c='k',lw=.5)\n",
    "plot_shaded_norm(tvec[tvec_analysis],1-compression_time_error_roi,0,color1,'$r$',ax)\n",
    "plot_shaded_norm(tvec[tvec_analysis],1-compression_time_error_ch,0,color2,'$\\lambda$',ax)\n",
    "ax.plot(tvec[tvec_analysis][p_values_sig],0.45*np.ones(sum(p_values_sig)),'ks',markersize=2)\n",
    "ax.set_xticks(np.arange(-.100,.601,.100))\n",
    "ax.legend(frameon=False)\n",
    "ax.set_xlim([-.100,.600])\n",
    "ax.set_xlabel('Time [s]')\n",
    "ax.set_ylabel('$c^{t}$')\n",
    "fancy_fig(ax)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=[4,3.5])\n",
    "vmax = 0.015\n",
    "vmin = 0.000\n",
    "cmap = 'magma'\n",
    "ax = plt.subplot(2,1,1)\n",
    "cax1 = plt.pcolormesh(power,compactness,prob_cond1.T, cmap=cmap,vmax=vmax,vmin=vmin,shading='auto')\n",
    "\n",
    "plt.ylabel('$c^{t}_{r}$')\n",
    "plt.xticks([])\n",
    "plt.yticks(fontsize=12)\n",
    "\n",
    "from matplotlib.ticker import FormatStrFormatter    \n",
    "ax.xaxis.set_major_formatter(FormatStrFormatter('%.3f'))\n",
    "cbar = plt.colorbar(cax1, ticks=[0, 0.005, 0.01, 0.015])\n",
    "\n",
    "ax = plt.subplot(2,1,2)\n",
    "cax2 = plt.pcolormesh(power,compactness,prob_cond2.T, cmap=cmap,vmax=vmax,vmin=vmin,shading='auto')\n",
    "plt.xlabel('Power')\n",
    "plt.ylabel('$c^{t}_{\\lambda}$')\n",
    "plt.xticks([0.003, 0.009, 0.015])\n",
    "plt.yticks(fontsize=12)\n",
    "ax.xaxis.set_major_formatter(FormatStrFormatter('%.3f'))\n",
    "cbar = plt.colorbar(cax2, ticks=[0, 0.005, 0.01, 0.015])\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,ax = plt.subplots(figsize=(6,4.5))\n",
    "data2plot = [correlation_roi,correlation_ch,correlation_rand_ch.mean(0),correlation_wwp.mean(0),correlation_PCA,correlation_ICA] \n",
    "ls = ['-','-','-','--','-','-']\n",
    "for i,label in enumerate(labels):\n",
    "    plot_shaded_norm(percentiles,data2plot[i],0,colors[i],label,ax,ls = ls[i])    \n",
    "plt.ylabel('Correlation')\n",
    "ax.set_xlabel('Percentile')\n",
    "ax.set_xlim(0,100)\n",
    "ax.set_ylim(0,1)\n",
    "\n",
    "ax.legend(frameon=False)\n",
    "fancy_fig(ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "    \n",
    "fig,ax = plt.subplots(figsize=(6,4.5))\n",
    "data2plot = [1-compression_error_roi,1-compression_error_ch,1-compression_error_rand_ch.mean(0),1-compression_error_wwp.mean(0),1-compression_error_PCA,1-compression_error_ICA] \n",
    "for i,label in enumerate(labels):\n",
    "    plot_shaded_norm(percentiles,data2plot[i],0,colors[i],label,ax,ls = ls[i])    \n",
    "plt.ylabel('Compactness')\n",
    "ax.set_xlabel('Percentile')\n",
    "ax.set_xlim(0,100)\n",
    "ax.set_ylim(0,1)\n",
    "\n",
    "ax.legend(frameon=False)\n",
    "fancy_fig(ax)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plots for Fig. S5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster_results_path = 'data/results_white_noise/'\n",
    "cluster_files = os.listdir(cluster_results_path)\n",
    "realisations = np.unique([int(file.split('[')[1].split(']')[0]) for file in cluster_files])\n",
    "\n",
    "compression_error_roi_wn  = []\n",
    "compression_error_ch_wn  = []\n",
    "compression_error_dsur_wn  = []\n",
    "compression_error_gsur_wn  = []\n",
    "compression_error_euc_wn  = []\n",
    "correlation_roi_wn  = []\n",
    "correlation_ch_wn  = []\n",
    "correlation_dsur_wn  = []\n",
    "correlation_gsur_wn  = []\n",
    "correlation_euc_wn  = []\n",
    "\n",
    "for r,realisation in enumerate(realisations):\n",
    "    compression_error_roi_wn.append(np.load(os.path.join(cluster_results_path,f'compression_err_roi_wn_r[{realisation}].npy')))\n",
    "    compression_error_ch_wn.append(np.load(os.path.join(cluster_results_path,f'compression_err_gft_wn_r[{realisation}].npy')))\n",
    "    compression_error_dsur_wn.append(np.load(os.path.join(cluster_results_path,f'compression_err_gft_sur_wn_r[{realisation}].npy')))\n",
    "    compression_error_gsur_wn.append(np.load(os.path.join(cluster_results_path,f'compression_err_gft_gsur_wn_r[{realisation}].npy')))\n",
    "    compression_error_euc_wn.append(np.load(os.path.join(cluster_results_path,f'compression_err_euc_wn_r[{realisation}].npy')))\n",
    "    correlation_roi_wn.append(np.load(os.path.join(cluster_results_path,f'correlation_roi_wn_r[{realisation}].npy')))\n",
    "    correlation_ch_wn.append(np.load(os.path.join(cluster_results_path,f'correlation_gft_wn_r[{realisation}].npy')))\n",
    "    correlation_dsur_wn.append(np.load(os.path.join(cluster_results_path,f'correlation_gft_sur_wn_r[{realisation}].npy')))\n",
    "    correlation_gsur_wn.append(np.load(os.path.join(cluster_results_path,f'correlation_gft_gsur_wn_r[{realisation}].npy')))\n",
    "    correlation_euc_wn.append(np.load(os.path.join(cluster_results_path,f'correlation_euc_wn_r[{realisation}].npy')))\n",
    "    \n",
    "compression_error_roi_wn = np.array(compression_error_roi_wn)[...,0]\n",
    "compression_error_ch_wn = np.array(compression_error_ch_wn)[...,0]\n",
    "compression_error_dsur_wn = np.array(compression_error_dsur_wn)[...,0]\n",
    "compression_error_gsur_wn = np.array(compression_error_gsur_wn)[...,0]\n",
    "compression_error_euc_wn = np.array(compression_error_euc_wn)[...,0]\n",
    "correlation_roi_wn = np.array(correlation_roi_wn)[...,0]\n",
    "correlation_ch_wn = np.array(correlation_ch_wn)[...,0]\n",
    "correlation_dsur_wn = np.array(correlation_dsur_wn)[...,0]\n",
    "correlation_gsur_wn = np.array(correlation_gsur_wn)[...,0]\n",
    "correlation_euc_wn = np.array(correlation_euc_wn)[...,0]\n",
    "\n",
    "compression_error_ch = np.load(os.path.join(compactness_path,'compression_GFT.npy'))\n",
    "correlation_ch = np.load(os.path.join(compactness_path,'correlation_GFT.npy'))\n",
    "\n",
    "compression_error_euc = np.load(os.path.join(compactness_path,'compression_euc.npy'))\n",
    "correlation_euc = np.load(os.path.join(compactness_path,'correlation_euc.npy')) \n",
    "\n",
    "data_comp = [\n",
    "    np.mean(1-compression_error_ch,-1),\n",
    "    np.mean(1-compression_error_euc,-1)\n",
    "]\n",
    "\n",
    "data_corr = [\n",
    "    np.mean(correlation_ch,-1),\n",
    "    np.mean(correlation_euc,-1)\n",
    "]\n",
    "\n",
    "\n",
    "data_comp_wn = [\n",
    "    np.mean(1-compression_error_ch_wn),\n",
    "    np.mean(1-compression_error_euc_wn),\n",
    "]\n",
    "\n",
    "data_corr_wn = [\n",
    "    np.mean(correlation_ch_wn),\n",
    "    np.mean(correlation_euc_wn),\n",
    "]\n",
    "\n",
    "data_comp_wn_flat = [\n",
    "    (1-compression_error_ch_wn).flatten(),\n",
    "    (1-compression_error_euc_wn).flatten(),\n",
    "]\n",
    "\n",
    "data_corr_wn_flat = [\n",
    "    (correlation_ch_wn).flatten(),\n",
    "    (correlation_euc_wn).flatten(),\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Figure S5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "compression_error_roi = np.load(os.path.join(compactness_path,'compression_ROI.npy'))\n",
    "correlation_roi = np.load(os.path.join(compactness_path,'correlation_ROI.npy'))\n",
    "\n",
    "compression_error_ch = np.load(os.path.join(compactness_path,'compression_GFT.npy'))\n",
    "correlation_ch = np.load(os.path.join(compactness_path,'correlation_GFT.npy'))\n",
    "\n",
    "compression_error_euc = np.load(os.path.join(compactness_path,'compression_euc.npy'))\n",
    "correlation_euc = np.load(os.path.join(compactness_path,'correlation_euc.npy')) \n",
    "\n",
    "data_comp = [\n",
    "    np.mean(1-compression_error_roi,-1),\n",
    "    np.mean(1-compression_error_ch,-1),\n",
    "    np.mean(1-compression_error_euc,-1)\n",
    "]\n",
    "\n",
    "data_corr = [\n",
    "    np.mean(correlation_roi,-1),\n",
    "    np.mean(correlation_ch,-1),\n",
    "    np.mean(correlation_euc,-1)\n",
    "]\n",
    "\n",
    "\n",
    "data_comp_wn = [\n",
    "    np.mean(1-compression_error_roi_wn),\n",
    "    np.mean(1-compression_error_ch_wn),\n",
    "    np.mean(1-compression_error_euc_wn),\n",
    "]\n",
    "\n",
    "data_corr_wn = [\n",
    "    np.mean(correlation_roi_wn),\n",
    "    np.mean(correlation_ch_wn),\n",
    "    np.mean(correlation_euc_wn),\n",
    "]\n",
    "\n",
    "data_comp_wn_flat = [\n",
    "    (1-compression_error_roi_wn).flatten(),\n",
    "    (1-compression_error_ch_wn).flatten(),\n",
    "    (1-compression_error_euc_wn).flatten(),\n",
    "]\n",
    "\n",
    "data_corr_wn_flat = [\n",
    "    (correlation_roi_wn).flatten(),\n",
    "    (correlation_ch_wn).flatten(),\n",
    "    (correlation_euc_wn).flatten(),\n",
    "]\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "font = {'size'   : 14}\n",
    "matplotlib.rc('font', **font)\n",
    "\n",
    "cmap = matplotlib.cm.get_cmap('Spectral')\n",
    "color1 = cmap(.05) # ROI \n",
    "color2 = cmap( .95) # CH\n",
    "cmap = matplotlib.cm.get_cmap('cividis')\n",
    "color5 = cmap(.35) # RND\n",
    "\n",
    "labels_delt = ['$r$','$\\lambda$','$\\mathcal{G}_{euc}$']\n",
    "colors_delt = [color1,color2,color5]\n",
    "\n",
    "data_delt_comp = [data_comp[i]-data_comp_wn[i] for i in range(len(labels_delt))]\n",
    "data_delt_corr = [data_corr[i]-data_corr_wn[i] for i in range(len(labels_delt))]\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,axs = plt.subplots(2,2,figsize=(8.5,5))\n",
    "\n",
    "ax = axs[0,0]\n",
    "ax = my_box_plot(data_corr_wn_flat,ax,colors_delt,labels_delt);\n",
    "ax.set_ylabel('Correlation')\n",
    "plt.setp(ax.get_xticklabels(), size=13)\n",
    "fancy_fig(ax)\n",
    "\n",
    "ax = axs[1,0]\n",
    "ax = my_box_plot(data_comp_wn_flat,ax,colors_delt,labels_delt);\n",
    "fancy_fig(ax)\n",
    "ax.set_ylabel('1-error')\n",
    "ax.set_xlabel('Coordinate system')\n",
    "plt.setp(ax.get_xticklabels(), size=13)\n",
    "fancy_fig(ax)\n",
    "\n",
    "ax = axs[0,1]\n",
    "ax = my_box_plot(data_delt_corr,ax,colors_delt,labels_delt);\n",
    "ax.set_ylabel('$\\delta_{correlation}$')\n",
    "plt.setp(ax.get_xticklabels(), size=13)\n",
    "fancy_fig(ax)\n",
    "\n",
    "ax = axs[1,1]\n",
    "ax = my_box_plot(data_delt_comp,ax,colors_delt,labels_delt);\n",
    "fancy_fig(ax)\n",
    "ax.set_ylabel('$\\delta_{1-error}$')\n",
    "ax.set_xlabel('Coordinate system')\n",
    "plt.setp(ax.get_xticklabels(), size=13)\n",
    "fancy_fig(ax)\n",
    "\n",
    "fig.tight_layout()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plots for Fig. S3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "roifname = os.path.join(datadir,'plotting','Lausanne2008_Yeo7RSNs.xlsx')\n",
    "roidata = pd.read_excel(roifname, sheet_name='SCALE {}'.format(scale))\n",
    "cort = np.where(roidata['Structure'] == 'cort')[0]\n",
    "ROIs = list(roidata['Label Lausanne2008'][cort])\n",
    "\n",
    "dorsal_rois_labels = ['cuneus', 'precuneus_1', 'precuneus_4', 'precuneus_5', 'inferiorparietal_1', 'inferiorparietal_3', 'superiorparietal', 'inferiorparietal_4', 'inferiorparietal_2', 'inferiorparietal_5']\n",
    "ventral_rois_labels = ['fusiform', 'inferiortemporal', 'middletemporal', 'parahippocampal', 'bankssts']\n",
    "early_rois_labels = ['occipital', 'pericalcarine', 'lingual', 'inferiorparietal_6']\n",
    "\n",
    "dorsal_rois = np.unique(np.concatenate([[i for i, elem in enumerate(ROIs) if roi in elem] for roi in dorsal_rois_labels]))\n",
    "ventral_rois = np.unique(np.concatenate([[i for i, elem in enumerate(ROIs) if roi in elem] for roi in ventral_rois_labels]))\n",
    "early_rois = np.unique(np.concatenate([[i for i, elem in enumerate(ROIs) if roi in elem] for roi in early_rois_labels]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute compression for the different visual streams"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "comp_roi_dorsal = np.zeros((len(percentiles), nsub_faces))\n",
    "comp_ch_dorsal = np.zeros((len(percentiles), nsub_faces))\n",
    "comp_roi_ventral = np.zeros((len(percentiles), nsub_faces))\n",
    "comp_ch_ventral = np.zeros((len(percentiles), nsub_faces))\n",
    "comp_roi_early = np.zeros((len(percentiles), nsub_faces))\n",
    "comp_ch_early = np.zeros((len(percentiles), nsub_faces))\n",
    "\n",
    "\n",
    "corr_roi_dorsal = np.zeros((len(percentiles), nsub_faces))\n",
    "corr_ch_dorsal = np.zeros((len(percentiles), nsub_faces))\n",
    "corr_roi_ventral = np.zeros((len(percentiles), nsub_faces))\n",
    "corr_ch_ventral = np.zeros((len(percentiles), nsub_faces))\n",
    "corr_roi_early = np.zeros((len(percentiles), nsub_faces))\n",
    "corr_ch_early = np.zeros((len(percentiles), nsub_faces))\n",
    "\n",
    "t = np.argmin(abs(tvec[tvec_analysis] - 0.170))\n",
    "for s in tqdm.tqdm(range(nsub_faces)):\n",
    "    # Load time-series for subject [ROIs x Time x Trials]\n",
    "    epochs, cond = ts_gen_faces.loader_ts(s)\n",
    "    # Select only FACES trials\n",
    "    epochs = epochs[:, :, cond==1] \n",
    "    # Baseline correction\n",
    "    epochs = epochs - epochs[:, tvec_pre].mean(1, keepdims=True)\n",
    "    # Graph Fourier transform\n",
    "    epochs_graph = sc.gft(epochs)\n",
    "    for p, percentile in enumerate(percentiles):\n",
    "\n",
    "        compressed = compress(epochs[:, t], percentile, 1)\n",
    "\n",
    "        comp_roi_dorsal[p, s] = compress_error(epochs[dorsal_rois, t].mean(1), compressed[dorsal_rois].mean(1))\n",
    "        comp_roi_ventral[p, s] = compress_error(epochs[ventral_rois, t].mean(1), compressed[ventral_rois].mean(1))\n",
    "        comp_roi_early[p, s] = compress_error(epochs[early_rois, t].mean(1), compressed[early_rois].mean(1))\n",
    "\n",
    "        corr_roi_dorsal[p, s] = np.corrcoef(epochs[dorsal_rois, t].mean(1), compressed[dorsal_rois].mean(1))[0, 1]\n",
    "        corr_roi_ventral[p, s] = np.corrcoef(epochs[ventral_rois, t].mean(1), compressed[ventral_rois].mean(1))[0, 1]\n",
    "        corr_roi_early[p, s] = np.corrcoef(epochs[early_rois, t].mean(1), compressed[early_rois].mean(1))[0, 1]\n",
    "\n",
    "        compressed = sc.igft(compress(epochs_graph[:, t], percentile, 1))\n",
    "\n",
    "        comp_ch_dorsal[p, s] = compress_error(epochs[dorsal_rois, t].mean(1), compressed[dorsal_rois].mean(1))\n",
    "        comp_ch_ventral[p, s] = compress_error(epochs[ventral_rois, t].mean(1), compressed[ventral_rois].mean(1))\n",
    "        comp_ch_early[p, s] = compress_error(epochs[early_rois, t].mean(1), compressed[early_rois].mean(1))\n",
    "\n",
    "        corr_ch_dorsal[p, s] = np.corrcoef(epochs[dorsal_rois, t].mean(1), compressed[dorsal_rois].mean(1))[0, 1]\n",
    "        corr_ch_ventral[p, s] = np.corrcoef(epochs[ventral_rois, t].mean(1), compressed[ventral_rois].mean(1))[0, 1]\n",
    "        corr_ch_early[p, s] = np.corrcoef(epochs[early_rois, t].mean(1), compressed[early_rois].mean(1))[0, 1]\n",
    "\n",
    "comp_roi_dorsal_motion = np.zeros((len(percentiles), nsub_motion))\n",
    "comp_ch_dorsal_motion = np.zeros((len(percentiles), nsub_motion))\n",
    "comp_roi_ventral_motion = np.zeros((len(percentiles), nsub_motion))\n",
    "comp_ch_ventral_motion = np.zeros((len(percentiles), nsub_motion))\n",
    "comp_roi_early_motion = np.zeros((len(percentiles), nsub_motion))\n",
    "comp_ch_early_motion = np.zeros((len(percentiles), nsub_motion))\n",
    "\n",
    "\n",
    "corr_roi_dorsal_motion = np.zeros((len(percentiles), nsub_motion))\n",
    "corr_ch_dorsal_motion = np.zeros((len(percentiles), nsub_motion))\n",
    "corr_roi_ventral_motion = np.zeros((len(percentiles), nsub_motion))\n",
    "corr_ch_ventral_motion = np.zeros((len(percentiles), nsub_motion))\n",
    "corr_roi_early_motion = np.zeros((len(percentiles), nsub_motion))\n",
    "corr_ch_early_motion = np.zeros((len(percentiles), nsub_motion))\n",
    "\n",
    "t = np.argmin(abs(tvec[tvec_analysis]-0.150))\n",
    "for s in tqdm.tqdm(range(nsub_motion)):\n",
    "    # Load time-series for subject [ROIs x Time x Trials]\n",
    "    epochs, cond = ts_gen_motion.loader_ts(s)\n",
    "    # Select only FACES trials\n",
    "    epochs = epochs[:, :, cond==1] \n",
    "    # Baseline correction\n",
    "    epochs = epochs - epochs[:, tvec_pre].mean(1, keepdims=True)\n",
    "    # Graph Fourier transform\n",
    "    epochs_graph = sc.gft(epochs)\n",
    "    for p, percentile in enumerate(percentiles):\n",
    "\n",
    "        compressed = compress(epochs[:, t], percentile, 1)\n",
    "\n",
    "        comp_roi_dorsal_motion[p, s] = compress_error(epochs[dorsal_rois, t].mean(1), compressed[dorsal_rois].mean(1))\n",
    "        comp_roi_ventral_motion[p, s] = compress_error(epochs[ventral_rois, t].mean(1), compressed[ventral_rois].mean(1))\n",
    "        comp_roi_early_motion[p, s] = compress_error(epochs[early_rois, t].mean(1), compressed[early_rois].mean(1))\n",
    "\n",
    "        corr_roi_dorsal_motion[p, s] = np.corrcoef(epochs[dorsal_rois, t].mean(1), compressed[dorsal_rois].mean(1))[0, 1]\n",
    "        corr_roi_ventral_motion[p, s] = np.corrcoef(epochs[ventral_rois, t].mean(1), compressed[ventral_rois].mean(1))[0, 1]\n",
    "        corr_roi_early_motion[p, s] = np.corrcoef(epochs[early_rois, t].mean(1), compressed[early_rois].mean(1))[0, 1]\n",
    "\n",
    "        compressed = sc.igft(compress(epochs_graph[:, t], percentile, 1))\n",
    "\n",
    "        comp_ch_dorsal_motion[p, s] = compress_error(epochs[dorsal_rois, t].mean(1), compressed[dorsal_rois].mean(1))\n",
    "        comp_ch_ventral_motion[p, s] = compress_error(epochs[ventral_rois, t].mean(1), compressed[ventral_rois].mean(1))\n",
    "        comp_ch_early_motion[p, s] = compress_error(epochs[early_rois, t].mean(1), compressed[early_rois].mean(1))\n",
    "\n",
    "        corr_ch_dorsal_motion[p, s] = np.corrcoef(epochs[dorsal_rois, t].mean(1), compressed[dorsal_rois].mean(1))[0, 1]\n",
    "        corr_ch_ventral_motion[p, s] = np.corrcoef(epochs[ventral_rois, t].mean(1), compressed[ventral_rois].mean(1))[0, 1]\n",
    "        corr_ch_early_motion[p, s] = np.corrcoef(epochs[early_rois, t].mean(1), compressed[early_rois].mean(1))[0, 1]\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plots Fig. S3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "data_comp = [1-comp_roi_dorsal.mean(0),1-comp_ch_dorsal.mean(0),1-comp_roi_ventral.mean(0),\n",
    " 1-comp_ch_ventral.mean(0),1-comp_roi_early.mean(0),1-comp_ch_early.mean(0)]\n",
    "\n",
    "data_comp_motion = [1-comp_roi_dorsal_motion.mean(0),1-comp_ch_dorsal_motion.mean(0),1-comp_roi_ventral_motion.mean(0),\n",
    " 1-comp_ch_ventral_motion.mean(0),1-comp_roi_early_motion.mean(0),1-comp_ch_early_motion.mean(0)]\n",
    "\n",
    "corr_roi_dorsal[np.isnan(corr_roi_dorsal)] = 0\n",
    "corr_roi_ventral[np.isnan(corr_roi_ventral)]   = 0\n",
    "corr_roi_early[np.isnan(corr_roi_early)] = 0\n",
    "\n",
    "corr_roi_dorsal_motion[np.isnan(corr_roi_dorsal_motion)] = 0\n",
    "corr_roi_ventral_motion[np.isnan(corr_roi_ventral_motion)] = 0\n",
    "corr_roi_early_motion[np.isnan(corr_roi_early_motion)] = 0\n",
    "\n",
    "data_cor = [corr_roi_dorsal.mean(0),corr_ch_dorsal.mean(0),corr_roi_ventral.mean(0),\n",
    " corr_ch_ventral.mean(0),corr_roi_early.mean(0),corr_ch_early.mean(0)]\n",
    "\n",
    "data_cor_motion =  [corr_roi_dorsal_motion.mean(0),corr_ch_dorsal_motion.mean(0),corr_roi_ventral_motion.mean(0),\n",
    " corr_ch_ventral_motion.mean(0),corr_roi_early_motion.mean(0),corr_ch_early_motion.mean(0)]\n",
    "\n",
    "\n",
    "cmap = matplotlib.cm.get_cmap('brg')\n",
    "colors =cmap([0, 0, 0.5, 0.5, 1, 1])\n",
    "colors[::2] = colors[::2]/1.5\n",
    "labels =['$r$', '$\\lambda$', '$r$', '$\\lambda$', '$r$', '$\\lambda$']\n",
    "\n",
    "fig,axs = plt.subplots(2,2,figsize=(14,6))\n",
    "\n",
    "\n",
    "ax = axs[0,0]\n",
    "\n",
    "ax = my_box_plot(data_cor,ax,colors,labels);\n",
    "ax.set_xticks([])\n",
    "ax.set_ylabel('Correlation')\n",
    "ax.set_title('Face 170ms')\n",
    "fancy_fig(ax)\n",
    "\n",
    "ax = axs[1,0]\n",
    "ax = my_box_plot(data_comp,ax,colors,labels);\n",
    "ax.set_ylabel('1 - error')\n",
    "fancy_fig(ax)\n",
    "ax.set_xlabel('Coordinate system')\n",
    "\n",
    "\n",
    "ax = axs[0,1]\n",
    "\n",
    "ax = my_box_plot(data_cor_motion,ax,colors,labels);\n",
    "ax.set_ylabel('Correlation')\n",
    "ax.set_title('Motion 150ms')\n",
    "ax.set_xticks([])\n",
    "fancy_fig(ax)\n",
    "\n",
    "ax = axs[1,1]\n",
    "ax = my_box_plot(data_comp_motion,ax,colors,labels);\n",
    "ax.set_ylabel('1 - error')\n",
    "fancy_fig(ax)\n",
    "ax.set_xlabel('Coordinate system')\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "strs = ['dorsal','ventral','early']\n",
    "pvalues = np.zeros((3,4))\n",
    "for i in range(3):\n",
    "    _,pvalues[i,0] = ranksums(data_comp[i*2],data_comp[i*2+1])\n",
    "    _,pvalues[i,1] = ranksums(data_cor[i*2],data_cor[i*2+1])\n",
    "    _,pvalues[i,2] = ranksums(data_comp_motion[i*2],data_comp_motion[i*2+1])\n",
    "    _,pvalues[i,3] = ranksums(data_cor_motion[i*2],data_cor_motion[i*2+1])\n",
    "\n",
    "# bonferroni correction\n",
    "pvalues = pvalues * np.multiply(*pvalues.shape)\n",
    "\n",
    "for i in range(3):\n",
    "    print('(Faces) Compression reconstruction {} ROI-vs-CH pval = {}'.format(strs[i],pvalues[i,0]))    \n",
    "    print('(Faces) Compression correlation {} ROI-vs-CH pval = {}'.format(strs[i],pvalues[i,1])) \n",
    "    \n",
    "    print('(Motion) Compression reconstruction {} ROI-vs-CH pval = {}'.format(strs[i],pvalues[i,2]))\n",
    "    print('(Motion) Compression correlation {} ROI-vs-CH pval = {}'.format(strs[i],pvalues[i,3]))     "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plots for Fig. 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute integartion / segregation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "psd, nn, x_d, x_c = sdi(sc,ts_gen_faces)\n",
    "d_norm_timeseries = np.linalg.norm(x_d,axis=0)\n",
    "c_norm_timeseries = np.linalg.norm(x_c,axis=0)\n",
    "sdi_norm_timeseries = d_norm_timeseries / c_norm_timeseries\n",
    "d_avg_timeseries = np.linalg.norm(x_d,2,0)\n",
    "c_avg_timeseries = np.linalg.norm(x_c,2,0)\n",
    "sdi_avg_timeseries = d_avg_timeseries / c_avg_timeseries\n",
    "\n",
    "psd_motion, nn_motion, x_d_motion, x_c_motion = sdi(sc,ts_gen_motion)\n",
    "d_norm_timeseries_motion = np.linalg.norm(x_d_motion,axis=0)\n",
    "c_norm_timeseries_motion = np.linalg.norm(x_c_motion,axis=0)\n",
    "sdi_norm_timeseries_motion = d_norm_timeseries_motion / c_norm_timeseries_motion\n",
    "d_avg_timeseries_motion = np.linalg.norm(x_d_motion,2,0)\n",
    "c_avg_timeseries_motion = np.linalg.norm(x_c_motion,2,0)\n",
    "sdi_avg_timeseries_motion = d_avg_timeseries_motion / c_avg_timeseries_motion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute integration / segregation for surrogate graphs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "if not load_data:\n",
    "    d_norm_timeseries_sur = []\n",
    "    c_norm_timeseries_sur = []\n",
    "    sdi_norm_timeseries_sur = []\n",
    "    d_avg_timeseries_sur = []\n",
    "    c_avg_timeseries_sur = []\n",
    "    sdi_avg_timeseries_sur = []\n",
    "\n",
    "    for s,sur in enumerate(sc_surro):\n",
    "        psd_sur, nn_sur, x_d_sur, x_c_sur = sdi(sur,ts_gen_faces)\n",
    "        d_norm_timeseries_surro = np.linalg.norm(x_d_sur,axis=0)\n",
    "        c_norm_timeseries_surro = np.linalg.norm(x_c_sur,axis=0)\n",
    "        sdi_norm_timeseries_surro = d_norm_timeseries_surro / c_norm_timeseries_surro\n",
    "        d_avg_timeseries_surro = np.linalg.norm(x_d_sur,2,0)\n",
    "        c_avg_timeseries_surro = np.linalg.norm(x_c_sur,2,0)\n",
    "        sdi_avg_timeseries_surro = d_avg_timeseries_surro / c_avg_timeseries_surro\n",
    "\n",
    "        d_norm_timeseries_sur.append(d_norm_timeseries_surro)\n",
    "        c_norm_timeseries_sur.append(c_norm_timeseries_surro)\n",
    "        sdi_norm_timeseries_sur.append(sdi_norm_timeseries_surro)\n",
    "        d_avg_timeseries_sur.append(d_avg_timeseries_surro)\n",
    "        c_avg_timeseries_sur.append(c_avg_timeseries_surro)\n",
    "        sdi_avg_timeseries_sur.append(sdi_avg_timeseries_surro)    \n",
    "\n",
    "    np.save(os.path.join(dynamics_path,'d_norm_timeseries_sur.npy'),d_norm_timeseries_sur)\n",
    "    np.save(os.path.join(dynamics_path,'c_norm_timeseries_sur.npy'),c_norm_timeseries_sur)\n",
    "    np.save(os.path.join(dynamics_path,'sdi_norm_timeseries_sur.npy'),sdi_norm_timeseries_sur)\n",
    "    np.save(os.path.join(dynamics_path,'d_avg_timeseries_sur.npy'),d_avg_timeseries_sur)\n",
    "    np.save(os.path.join(dynamics_path,'c_avg_timeseries_sur.npy'),c_avg_timeseries_sur)\n",
    "    np.save(os.path.join(dynamics_path,'sdi_avg_timeseries_sur.npy'),sdi_avg_timeseries_sur)\n",
    "\n",
    "else:\n",
    "    d_norm_timeseries_sur = np.load(os.path.join(dynamics_path,'d_norm_timeseries_sur.npy'))\n",
    "    c_norm_timeseries_sur = np.load(os.path.join(dynamics_path,'c_norm_timeseries_sur.npy'))\n",
    "    sdi_norm_timeseries_sur = np.load(os.path.join(dynamics_path,'sdi_norm_timeseries_sur.npy'))\n",
    "    d_avg_timeseries_sur = np.load(os.path.join(dynamics_path,'d_avg_timeseries_sur.npy'))\n",
    "    c_avg_timeseries_sur = np.load(os.path.join(dynamics_path,'c_avg_timeseries_sur.npy'))\n",
    "    sdi_avg_timeseries_sur = np.load(os.path.join(dynamics_path,'sdi_avg_timeseries_sur.npy'))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not load_data:\n",
    "    d_norm_timeseries_sur_motion = []\n",
    "    c_norm_timeseries_sur_motion = []\n",
    "    sdi_norm_timeseries_sur_motion = []\n",
    "    d_avg_timeseries_sur_motion = []\n",
    "    c_avg_timeseries_sur_motion = []\n",
    "    sdi_avg_timeseries_sur_motion = []\n",
    "\n",
    "    for s,sur in enumerate(sc_surro):\n",
    "        psd_sur, nn_sur, x_d_sur, x_c_sur = sdi(sur,ts_gen_motion)\n",
    "        d_norm_timeseries_surro = np.linalg.norm(x_d_sur,axis=0)\n",
    "        c_norm_timeseries_surro = np.linalg.norm(x_c_sur,axis=0)\n",
    "        sdi_norm_timeseries_surro = d_norm_timeseries_surro / c_norm_timeseries_surro\n",
    "        d_avg_timeseries_surro = np.linalg.norm(x_d_sur,2,0)\n",
    "        c_avg_timeseries_surro = np.linalg.norm(x_c_sur,2,0)\n",
    "        sdi_avg_timeseries_surro = d_avg_timeseries_surro / c_avg_timeseries_surro\n",
    "\n",
    "        d_norm_timeseries_sur_motion.append(d_norm_timeseries_surro)\n",
    "        c_norm_timeseries_sur_motion.append(c_norm_timeseries_surro)\n",
    "        sdi_norm_timeseries_sur_motion.append(sdi_norm_timeseries_surro)\n",
    "        d_avg_timeseries_sur_motion.append(d_avg_timeseries_surro)\n",
    "        c_avg_timeseries_sur_motion.append(c_avg_timeseries_surro)\n",
    "        sdi_avg_timeseries_sur_motion.append(sdi_avg_timeseries_surro)    \n",
    "\n",
    "    np.save(os.path.join(dynamics_path,'d_norm_timeseries_sur_motion.npy'),d_norm_timeseries_sur_motion)\n",
    "    np.save(os.path.join(dynamics_path,'c_norm_timeseries_sur_motion.npy'),c_norm_timeseries_sur_motion)\n",
    "    np.save(os.path.join(dynamics_path,'sdi_norm_timeseries_sur_motion.npy'),sdi_norm_timeseries_sur_motion)\n",
    "    np.save(os.path.join(dynamics_path,'d_avg_timeseries_sur_motion.npy'),d_avg_timeseries_sur_motion)\n",
    "    np.save(os.path.join(dynamics_path,'c_avg_timeseries_sur_motion.npy'),c_avg_timeseries_sur_motion)\n",
    "    np.save(os.path.join(dynamics_path,'sdi_avg_timeseries_sur_motion.npy'),sdi_avg_timeseries_sur_motion)\n",
    "\n",
    "else:\n",
    "    d_norm_timeseries_sur_motion = np.load(os.path.join(dynamics_path,'d_norm_timeseries_sur_motion.npy'))\n",
    "    c_norm_timeseries_sur_motion = np.load(os.path.join(dynamics_path,'c_norm_timeseries_sur_motion.npy'))\n",
    "    sdi_norm_timeseries_sur_motion = np.load(os.path.join(dynamics_path,'sdi_norm_timeseries_sur_motion.npy'))\n",
    "    d_avg_timeseries_sur_motion = np.load(os.path.join(dynamics_path,'d_avg_timeseries_sur_motion.npy'))\n",
    "    c_avg_timeseries_sur_motion = np.load(os.path.join(dynamics_path,'c_avg_timeseries_sur_motion.npy'))\n",
    "    sdi_avg_timeseries_sur_motion = np.load(os.path.join(dynamics_path,'sdi_avg_timeseries_sur_motion.npy'))\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute graph signal smoothness "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "signal_smoothness = np.zeros((ts_gen_faces.nsub,len(tvec_analysis)))\n",
    "for s in tqdm.tqdm(range(ts_gen_faces.nsub)):\n",
    "    # Load time-series for subject [ROIs x Time x Trials]\n",
    "    epochs,cond = ts_gen_faces.loader_ts(s)\n",
    "    # Select only FACES trials\n",
    "    epochs = epochs[:,:,cond==1] \n",
    "        \n",
    "    for t in range(len(tvec_analysis)):    \n",
    "        # Signal smoothness in Graph\n",
    "        signal_smoothness[s,t] = smoothness(epochs[:,t].mean(1),sc)\n",
    "            "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot broadcasting dynamics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "save_figs = True\n",
    "font = {'size'   : 16}\n",
    "matplotlib.rc('font', **font)\n",
    "\n",
    "cmap = matplotlib.cm.get_cmap('PRGn')\n",
    "colorL = cmap(.05) \n",
    "colorH = cmap(.95) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(4,2))\n",
    "\n",
    "plot_shaded(sc.e[:nn-1], psd[:nn-1],1,colorL,'$P_{L}$',ax)\n",
    "plot_shaded(sc.e[nn-1:], psd[nn-1:],1,colorH,'$P_{H}$',ax)\n",
    "ax.axvline(np.mean([sc.e[nn-1],sc.e[nn-2]]),color='k')\n",
    "ax.set_xlabel('$\\lambda$')\n",
    "ax.set_ylabel('Power \\n Spectral Density')\n",
    "ax.set_yscale('log')\n",
    "ax.set_xlim([np.min(sc.e),np.max(sc.e)])\n",
    "fancy_fig(ax)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,axs = plt.subplots(figsize=(4,2))\n",
    "axs.axvline(0,ls='--',c='k',lw=1)\n",
    "plot_shaded(tvec[tvec_analysis], d_avg_timeseries,1,colorH,'$P_{H}$',axs)\n",
    "plot_shaded(tvec[tvec_analysis], c_avg_timeseries,1,colorL,'$P_{L}$',axs)\n",
    "axs.set_xlim((np.min(tvec[tvec_analysis]),np.max(tvec[tvec_analysis])))\n",
    "axs.legend(loc='upper right',frameon=False)\n",
    "axs.set_xlabel('Time [s]')\n",
    "axs.set_ylabel('Power [a.u]')\n",
    "axs.set_ylim((-1,38))\n",
    "fancy_fig(axs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,axs = plt.subplots(figsize=(4,2))\n",
    "axs.axvline(0,ls='--',c='k',lw=1)\n",
    "plot_shaded(tvec[tvec_analysis], d_avg_timeseries/power_in_time.T,1,colorH,'norm-$P_{H}$',axs)\n",
    "plot_shaded(tvec[tvec_analysis], c_avg_timeseries/power_in_time.T,1,colorL,'norm-$P_{L}$',axs)\n",
    "axs.set_xlim((np.min(tvec[tvec_analysis]),np.max(tvec[tvec_analysis])))\n",
    "axs.set_xlabel('Time [s]')\n",
    "axs.set_ylabel('Power [a.u]')\n",
    "fancy_fig(axs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,axs = plt.subplots(figsize=(4,2))\n",
    "axs.axvline(0,ls='--',c='k',lw=1)\n",
    "axs_b = plt.twinx(axs)\n",
    "sdi_norm = sdi_avg_timeseries\n",
    "cmp_time_norm = compression_time_error_ch \n",
    "plot_shaded(tvec[tvec_analysis], sdi_norm,1,'gray','SDI',axs)\n",
    "plot_shaded(tvec[tvec_analysis], 1-cmp_time_norm,0,color2,'CH compactness',axs_b)\n",
    "axs.set_xlim((np.min(tvec[tvec_analysis]),np.max(tvec[tvec_analysis])))\n",
    "lns = axs.get_lines() + axs_b.get_lines()\n",
    "labs = [l.get_label() for l in lns]\n",
    "axs.legend(lns, labs, loc=0,frameon=False)\n",
    "axs.set_xlabel('Time [s]')\n",
    "axs.set_ylabel('Power [a.u]')\n",
    "axs.set_ylabel('Compactness')\n",
    "fancy_fig(axs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,ax  = plt.subplots(figsize=(4,2))\n",
    "ax.axvline(0,ls='--',c='k',lw=1)\n",
    "plot_shaded(tvec[tvec_analysis],signal_smoothness/power_in_time,axis=0,color='k',label='SC',ax=ax)\n",
    "ax.set_xlim((np.min(tvec[tvec_analysis]),np.max(tvec[tvec_analysis])))\n",
    "ax.set_xlabel('Time [s]')\n",
    "ax.set_ylabel('Dirichlet energy')\n",
    "ax.set_yscale('log')\n",
    "fancy_fig(ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,axs = plt.subplots(figsize=(4,2))\n",
    "sdi_norm = sdi_avg_timeseries\n",
    "cmp_time_norm = compression_time_error_ch \n",
    "axs.plot((d_avg_timeseries/power_in_time.T).T.flatten(),1-compression_time_error_ch.flatten(),'o',c=color2,ms=6,alpha=0.15)\n",
    "axs.set_xlabel('$P_{H}$')\n",
    "axs.set_ylabel('$c^{t}$')\n",
    "fancy_fig(axs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plots for Fig. S4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# SUPP\n",
    "fig, ax = plt.subplots(figsize=(4,2))\n",
    "\n",
    "plot_shaded(sc.e[:nn_motion-1], psd_motion[:nn_motion-1],1,colorL,'$P_{L}$',ax)\n",
    "plot_shaded(sc.e[nn_motion-1:], psd_motion[nn_motion-1:],1,colorH,'$P_{H}$',ax)\n",
    "ax.axvline(np.mean([sc.e[nn_motion-1],sc.e[nn_motion-2]]),color='k')\n",
    "ax.set_xlabel('$\\lambda$')\n",
    "ax.set_ylabel('Power \\n Spectral Density')\n",
    "ax.set_yscale('log')\n",
    "ax.set_xlim([np.min(sc.e),np.max(sc.e)])\n",
    "fancy_fig(ax)\n",
    "\n",
    "fig,axs = plt.subplots(figsize=(4,2))\n",
    "axs.axvline(0,ls='--',c='k',lw=1)\n",
    "plot_shaded(tvec[tvec_analysis], d_avg_timeseries_motion,1,colorH,'$P_{H}$',axs)\n",
    "plot_shaded(tvec[tvec_analysis], c_avg_timeseries_motion,1,colorL,'$P_{L}$',axs)\n",
    "axs.set_xlim((np.min(tvec[tvec_analysis]),np.max(tvec[tvec_analysis])))\n",
    "axs.legend(loc='upper right',frameon=False)\n",
    "axs.set_xlabel('Time [s]')\n",
    "axs.set_ylabel('Power [a.u]')\n",
    "axs.set_ylim([-1,38])\n",
    "fancy_fig(axs)\n",
    "\n",
    "fig,axs = plt.subplots(figsize=(4,2))\n",
    "axs.axvline(0,ls='--',c='k',lw=1)\n",
    "plot_shaded(tvec[tvec_analysis], d_avg_timeseries_motion/power_in_time_motion.T,1,colorH,'norm-$P_{H}$',axs)\n",
    "plot_shaded(tvec[tvec_analysis], c_avg_timeseries_motion/power_in_time_motion.T,1,colorL,'norm-$P_{L}$',axs)\n",
    "axs.set_xlim((np.min(tvec[tvec_analysis]),np.max(tvec[tvec_analysis])))\n",
    "#axs.legend()\n",
    "axs.set_xlabel('Time [s]')\n",
    "axs.set_ylabel('Power [a.u]')\n",
    "fancy_fig(axs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot broadcasting dynamics with surrogate analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute significance for broadcasting dynamics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.stats.multitest import fdrcorrection\n",
    "\n",
    "pvalue = np.zeros((ts_gen_faces.nsub,len(ts_gen_faces.tvec_analysis)))\n",
    "\n",
    "for s in range(ts_gen_faces.nsub):\n",
    "    tmp = np.sum((np.array(sdi_avg_timeseries_sur)[:,:,0]).T>sdi_avg_timeseries[:,s][:,np.newaxis],1)\n",
    "    pvalue[s] = (tmp + 1) / (len(tmp)+1)\n",
    "    \n",
    "#significant_t = (pvalue<=0.025) | (pvalue>=0.975)\n",
    "\n",
    "pvalue[pvalue>=0.5] = 1 - pvalue[pvalue>=0.5]\n",
    "# FDR\n",
    "_, pvalue = fdrcorrection(pvalue.flatten())\n",
    "pvalue = pvalue.reshape((ts_gen_faces.nsub,len(ts_gen_faces.tvec_analysis)))\n",
    "significant_t = (pvalue<=0.025)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "significant_t_vals = np.ones_like(significant_t)*np.nan\n",
    "dif = d_avg_timeseries/power_in_time.T-c_avg_timeseries/power_in_time.T\n",
    "significant_t_vals[significant_t==1] = dif.T[significant_t==1]\n",
    "order = np.argsort(np.nansum(significant_t_vals,1))[::-1]\n",
    "\n",
    "font = {'size'   : 16}\n",
    "matplotlib.rc('font', **font)\n",
    "\n",
    "fig = plt.figure(figsize=(10,5),constrained_layout=True)\n",
    "gs = fig.add_gridspec(3, 1)\n",
    "\n",
    "ax0 = fig.add_subplot(gs[:2])\n",
    "maxval = np.nanmax(abs(significant_t_vals))\n",
    "cmap = matplotlib.cm.get_cmap('PRGn').copy()\n",
    "cmap.set_bad(color='gray')\n",
    "im = ax0.pcolormesh(tvec[ts_gen_faces.tvec_analysis],np.arange(ts_gen_faces.nsub),\n",
    "                    significant_t_vals[order],cmap=cmap,vmax=maxval,vmin=-maxval,\n",
    "                    shading='auto')\n",
    "ax0.set_ylabel('Subject')\n",
    "ax0.set_xlabel('Time [s]')\n",
    "ax0.axvline(0,ls='--',c='k')\n",
    "cax = fig.colorbar(im,ax = ax0)\n",
    "ax0.tick_params(width=3, length=8)\n",
    "ax0.spines['top'].set_linewidth(3)\n",
    "ax0.spines['bottom'].set_linewidth(3)\n",
    "ax0.spines['left'].set_linewidth(3)\n",
    "ax0.spines['right'].set_linewidth(3)\n",
    "cax.ax.tick_params(width=3, length=8)\n",
    "cax.ax.spines['top'].set_linewidth(3)\n",
    "cax.ax.spines['bottom'].set_linewidth(3)\n",
    "cax.ax.spines['left'].set_linewidth(3)\n",
    "cax.ax.spines['right'].set_linewidth(3)\n",
    "\n",
    "ax1 = fig.add_subplot(gs[2],sharex=ax0)\n",
    "ax1.plot(tvec[ts_gen_faces.tvec_analysis],significant_t.mean(0)*100,'gray',lw=2)\n",
    "ax1.axvline(0,c='k',ls='--')\n",
    "ax1.set_xlim(np.min(tvec[ts_gen_faces.tvec_analysis]),np.max(tvec[ts_gen_faces.tvec_analysis]))\n",
    "ax1.set_ylim(0,100)\n",
    "ax1.set_xlabel('Time [s]')\n",
    "ax1.set_ylabel('% of subjects\\nwith significant\\n direction',color='gray')\n",
    "fancy_fig(ax1)\n",
    "ax1_t = plt.twinx(ax1)\n",
    "ax1_t.plot(tvec[ts_gen_faces.tvec_analysis],dif.mean(1),c='r',lw=2)\n",
    "ax1_t.axhline(0,c='r',ls='--',alpha=.5)\n",
    "ax1_t.set_ylabel('Broadcasting\\ndirection',color='r')\n",
    "fancy_fig(ax1_t)\n",
    "ax1_t.spines['right'].set_linewidth(3)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Broadcasting dynamics for motion dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pvalue_motion = np.zeros((ts_gen_motion.nsub,len(ts_gen_motion.tvec_analysis)))\n",
    "\n",
    "for s in range(ts_gen_motion.nsub):\n",
    "    tmp = np.sum((np.array(sdi_avg_timeseries_sur_motion)[:,:,0]).T>sdi_avg_timeseries_motion[:,s][:,np.newaxis],1)\n",
    "    pvalue_motion[s] = (tmp + 1) / (len(tmp)+1)\n",
    "    \n",
    "pvalue_motion[pvalue_motion>=0.5] = 1 - pvalue_motion[pvalue_motion>=0.5]\n",
    "# FDR\n",
    "_, pvalue_motion = fdrcorrection(pvalue_motion.flatten())\n",
    "pvalue_motion = pvalue_motion.reshape((ts_gen_motion.nsub,len(ts_gen_motion.tvec_analysis)))\n",
    "significant_t_motion = (pvalue_motion<=0.025)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "significant_t_vals_motion = np.ones_like(significant_t_motion)*np.nan\n",
    "dif_motion = d_avg_timeseries_motion/power_in_time_motion.T-c_avg_timeseries_motion/power_in_time_motion.T\n",
    "significant_t_vals_motion[significant_t_motion==1] = dif_motion.T[significant_t_motion==1]\n",
    "order = np.argsort(np.nansum(significant_t_vals_motion,1))[::-1]\n",
    "\n",
    "font = {'size'   : 16}\n",
    "matplotlib.rc('font', **font)\n",
    "\n",
    "fig = plt.figure(figsize=(13,5),constrained_layout=True)\n",
    "gs = fig.add_gridspec(3, 1)\n",
    "\n",
    "ax0 = fig.add_subplot(gs[:2])\n",
    "maxval = np.nanmax(abs(significant_t_vals_motion))\n",
    "cmap = matplotlib.cm.get_cmap('PRGn').copy()\n",
    "cmap.set_bad(color='gray')\n",
    "im = ax0.pcolormesh(tvec[ts_gen_motion.tvec_analysis],np.arange(ts_gen_motion.nsub),\n",
    "                    significant_t_vals_motion[order],cmap=cmap,vmax=maxval,vmin=-maxval,\n",
    "                    shading='auto')\n",
    "ax0.set_ylabel('Subject')\n",
    "ax0.set_xlabel('Time [s]')\n",
    "ax0.axvline(0,ls='--',c='k')\n",
    "cax = fig.colorbar(im,ax = ax0)\n",
    "ax0.tick_params(width=3, length=8)\n",
    "ax0.spines['top'].set_linewidth(3)\n",
    "ax0.spines['bottom'].set_linewidth(3)\n",
    "ax0.spines['left'].set_linewidth(3)\n",
    "ax0.spines['right'].set_linewidth(3)\n",
    "cax.ax.tick_params(width=3, length=8)\n",
    "cax.ax.spines['top'].set_linewidth(3)\n",
    "cax.ax.spines['bottom'].set_linewidth(3)\n",
    "cax.ax.spines['left'].set_linewidth(3)\n",
    "cax.ax.spines['right'].set_linewidth(3)\n",
    "\n",
    "ax1 = fig.add_subplot(gs[2],sharex=ax0)\n",
    "ax1.plot(tvec[ts_gen_motion.tvec_analysis],significant_t_motion.mean(0)*100,'gray',lw=2)\n",
    "ax1.axvline(0,c='k',ls='--')\n",
    "ax1.set_xlim(np.min(tvec[ts_gen_motion.tvec_analysis]),np.max(tvec[ts_gen_motion.tvec_analysis]))\n",
    "ax1.set_ylim(0,100)\n",
    "ax1.set_xlabel('Time [s]')\n",
    "ax1.set_ylabel('% of subjects\\nwith significant\\n direction',color='gray')\n",
    "fancy_fig(ax1)\n",
    "ax1_t = plt.twinx(ax1)\n",
    "ax1_t.plot(tvec[ts_gen_motion.tvec_analysis],dif_motion.mean(1),c='r',lw=2)\n",
    "ax1_t.axhline(0,c='r',ls='--',alpha=.5)\n",
    "ax1_t.set_ylabel('Broadcasting\\ndirection',color='r')\n",
    "fancy_fig(ax1_t)\n",
    "ax1_t.spines['right'].set_linewidth(3)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plots for figure 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute distance between trajectories"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dist_trajectories_roi = np.zeros((len(percentiles),len(tvec_analysis),nsub_faces))\n",
    "dist_trajectories_ch = np.zeros((len(percentiles),len(tvec_analysis),nsub_faces))\n",
    "dist_trajectories2_ch = np.zeros((len(percentiles),len(tvec_analysis),nsub_faces))\n",
    "\n",
    "diff_roi = face_roi_ts - scra_roi_ts\n",
    "diff_ch  =  face_ch_ts - scra_ch_ts\n",
    "\n",
    "diff_power = np.linalg.norm(diff_roi.mean(-1),2,1)    \n",
    "norm_roi_face = np.linalg.norm(face_roi_ts.mean(-1),2,1)    \n",
    "norm_ch_face = np.linalg.norm(face_ch_ts.mean(-1),2,1)    \n",
    "\n",
    "for p,percentile in enumerate(percentiles):    \n",
    "\n",
    "    compressed_epochs_diff = compress(face_roi_ts-scra_roi_ts,percentile,(1,2))    \n",
    "    roi_lowd = np.where( compressed_epochs_diff.sum((1,2)))[0]\n",
    "\n",
    "    compressed_epochs_diff = compress(face_ch_ts-scra_ch_ts,percentile,(1,2))    \n",
    "    ch_lowd = np.where( compressed_epochs_diff.sum((1,2)))[0]\n",
    "    \n",
    "    diff_ch_tmp = np.zeros_like(diff_ch)\n",
    "    diff_ch_tmp[ch_lowd] = diff_ch[ch_lowd]     \n",
    "    dist_trajectories_roi[p] = np.linalg.norm(diff_roi[roi_lowd],1,0)\n",
    "    dist_trajectories_ch[p] = np.linalg.norm(sc.igft(diff_ch_tmp),1,0)\n",
    "\n",
    "distance = (dist_trajectories_ch-dist_trajectories_roi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Select 3D-space dimensions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "compressed_epochs = compress(diff_roi, 99, (1, 2))    \n",
    "roi_lowd = np.where( compressed_epochs.sum((1, 2)))[0]\n",
    "\n",
    "diff_ch_z = diff_ch.copy()\n",
    "diff_ch_z[0] = 0\n",
    "compressed_epochs = compress(diff_ch_z, 99, (1, 2))    \n",
    "ch_lowd = np.where( compressed_epochs.sum((1, 2)))[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute correlation among dimensions in high-dimneisonal space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "if not load_data:\n",
    "\n",
    "    corr_evoked = np.zeros(nsub_faces)\n",
    "    corr_graph =  np.zeros(nsub_faces)\n",
    "    corr_pca =  np.zeros(nsub_faces)\n",
    "    corr_ica =  np.zeros(nsub_faces)\n",
    "\n",
    "    for s in tqdm.tqdm(range(nsub_faces)):\n",
    "        # Load time-series for subject [ROIs x Time x Trials]\n",
    "        epochs, cond = ts_gen_faces.loader_ts(s)\n",
    "        # Select only FACES trials\n",
    "        epochs = epochs[:, :, cond==1] \n",
    "        # Baseline correction\n",
    "        epochs = epochs - epochs[:, tvec_pre].mean(1, keepdims=True)\n",
    "        # Graph Fourier transform\n",
    "        epochs_graph = sc.gft(epochs)\n",
    "\n",
    "        evoked = epochs.mean(2)\n",
    "        evoked_graph = epochs_graph.mean(2)\n",
    "        pca = PCA()\n",
    "        ica = FastICA(max_iter=500)\n",
    "        epochs_pca = pca.fit_transform(epochs.reshape(sc.N, -1).T).T.reshape(*epochs.shape)\n",
    "        evoked_pca = epochs_pca.mean(2)\n",
    "        epochs_ica = ica.fit_transform(epochs.reshape(sc.N, -1).T).T.reshape(*epochs.shape)\n",
    "        evoked_ica = epochs_ica.mean(2)\n",
    "\n",
    "        tmp = np.corrcoef(evoked)\n",
    "        corr_evoked[s] = np.mean(abs(tmp[np.tril_indices(sc.N, -1)]))\n",
    "\n",
    "        tmp = np.corrcoef(evoked_graph)\n",
    "        corr_graph[s] = np.mean(abs(tmp[np.tril_indices(sc.N, -1)]))\n",
    "\n",
    "        tmp = np.corrcoef(evoked_pca)\n",
    "        corr_pca[s] = np.mean(abs(tmp[np.tril_indices(sc.N, -1)]))\n",
    "\n",
    "        tmp = np.corrcoef(evoked_ica)\n",
    "        corr_ica[s] = np.mean(abs(tmp[np.tril_indices(sc.N, -1)]))\n",
    "        \n",
    "    np.save(os.path.join(corr_path, 'corr_evoked.npy'), corr_evoked)\n",
    "    np.save(os.path.join(corr_path, 'corr_graph.npy'), corr_graph)\n",
    "    np.save(os.path.join(corr_path, 'corr_pca.npy'), corr_pca)\n",
    "    np.save(os.path.join(corr_path, 'corr_ica.npy'), corr_ica)\n",
    "    \n",
    "else:\n",
    "    corr_evoked = np.load(os.path.join(corr_path, 'corr_evoked.npy'))\n",
    "    corr_graph = np.load(os.path.join(corr_path, 'corr_graph.npy'))\n",
    "    corr_pca = np.load(os.path.join(corr_path, 'corr_pca.npy'))\n",
    "    corr_ica = np.load(os.path.join(corr_path, 'corr_ica.npy'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Print mean correlation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(corr_evoked.mean())\n",
    "print(corr_graph.mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare statistically"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_,pvalue = ranksums(corr_evoked, corr_graph)\n",
    "print(pvalue)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute correlation among dimensions in high-dimneisonal space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not load_data:\n",
    "    corr_evoked_lowd = np.zeros(nsub_faces)\n",
    "    corr_graph_lowd =  np.zeros(nsub_faces)\n",
    "    corr_pca_lowd =  np.zeros(nsub_faces)\n",
    "    corr_ica_lowd =  np.zeros(nsub_faces)\n",
    "\n",
    "    for s in tqdm.tqdm(range(nsub_faces)):\n",
    "        # Load time-series for subject [ROIs x Time x Trials]\n",
    "        epochs, cond = ts_gen_faces.loaderTs(s)\n",
    "        # Select only FACES trials\n",
    "        epochs = epochs[:, :, cond==1] \n",
    "        # Baseline correction\n",
    "        epochs = epochs - epochs[:, tvec_pre].mean(1, keepdims=True)\n",
    "        # Graph Fourier transform\n",
    "        epochs_graph = sc.gft(epochs)\n",
    "\n",
    "        evoked = epochs.mean(2)\n",
    "        evoked_graph = epochs_graph.mean(2)\n",
    "        pca = PCA()\n",
    "        ica = FastICA(max_iter=500)\n",
    "        epochs_pca = pca.fit_transform(epochs.reshape(sc.N, -1).T).T.reshape(*epochs.shape)\n",
    "        evoked_pca = epochs_pca.mean(2)\n",
    "        epochs_ica = ica.fit_transform(epochs.reshape(sc.N, -1).T).T.reshape(*epochs.shape)\n",
    "        evoked_ica = epochs_ica.mean(2)\n",
    "\n",
    "        tmp = np.corrcoef(evoked[roi_lowd])\n",
    "        corr_evoked_lowd[s] = np.mean(abs(tmp[np.tril_indices(len(roi_lowd), -1)]))\n",
    "\n",
    "        tmp = np.corrcoef(evoked_graph[[nh_lowd]])\n",
    "        corr_graph_lowd[s] = np.mean(abs(tmp[np.tril_indices(len(nh_lowd), -1)]))\n",
    "\n",
    "        tmp = np.corrcoef(evoked_pca[:len(roi_lowd)])\n",
    "        corr_pca_lowd[s] = np.mean(abs(tmp[np.tril_indices(len(roi_lowd), -1)]))\n",
    "\n",
    "        tmp = np.corrcoef(evoked_ica[:len(roi_lowd)])\n",
    "        corr_ica_lowd[s] = np.mean(abs(tmp[np.tril_indices(len(roi_lowd), -1)]))\n",
    "\n",
    "        np.save(os.path.join(corr_path, 'corr_evoked_lowd.npy'), corr_evoked_lowd)\n",
    "        np.save(os.path.join(corr_path, 'corr_graph_lowd.npy'), corr_graph_lowd)\n",
    "        np.save(os.path.join(corr_path, 'corr_pca_lowd.npy'), corr_pca_lowd)\n",
    "        np.save(os.path.join(corr_path, 'corr_ica_lowd.npy'), corr_ica_lowd)\n",
    "\n",
    "\n",
    "else:\n",
    "    corr_evoked_lowd = np.load(os.path.join(corr_path, 'corr_evoked_lowd.npy'))\n",
    "    corr_graph_lowd = np.load(os.path.join(corr_path, 'corr_graph_lowd.npy'))\n",
    "    corr_pca_lowd = np.load(os.path.join(corr_path, 'corr_pca_lowd.npy'))\n",
    "    corr_ica_lowd = np.load(os.path.join(corr_path, 'corr_ica_lowd.npy'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Print mean correlation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(corr_evoked_lowd.mean())\n",
    "print(corr_graph_lowd.mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare statistically"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_,pvalue = ranksums(corr_evoked_lowd, corr_graph_lowd)\n",
    "print(pvalue)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plots for Fig. 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "lw = 1\n",
    "s = 50\n",
    "alpha = 0.9\n",
    "ec = 'gray'\n",
    "\n",
    "fig = plt.figure(figsize=(3,3))\n",
    "ax = fig.add_subplot(111,projection='3d')\n",
    "ax.plot(face_roi_ts[roi_lowd[1]].mean(1),face_roi_ts[roi_lowd[0]].mean(1),face_roi_ts[roi_lowd[2]].mean(1),'.-',c='gray',lw=lw)\n",
    "ax.plot(scra_roi_ts[roi_lowd[1]].mean(1),scra_roi_ts[roi_lowd[0]].mean(1),scra_roi_ts[roi_lowd[2]].mean(1),'.-',c='gray',lw=lw)\n",
    "ax.scatter(face_roi_ts[roi_lowd[1]].mean(1),face_roi_ts[roi_lowd[0]].mean(1),face_roi_ts[roi_lowd[2]].mean(1),s=s,c=cm.Blues(tvec[tvec_analysis]),edgecolor=ec,alpha=alpha,lw=0.4)\n",
    "ax.scatter(scra_roi_ts[roi_lowd[1]].mean(1),scra_roi_ts[roi_lowd[0]].mean(1),scra_roi_ts[roi_lowd[2]].mean(1),s=s,c=cm.Oranges(tvec[tvec_analysis]),edgecolor=ec,alpha=alpha,lw=0.4)\n",
    "ax.plot(face_roi_ts[roi_lowd[1]].mean(1)[[0,-1]],face_roi_ts[roi_lowd[0]].mean(1)[[0,-1]],face_roi_ts[roi_lowd[2]].mean(1)[[0,-1]],'o',c='k',ms=10)\n",
    "ax.plot(scra_roi_ts[roi_lowd[1]].mean(1)[[0,-1]],scra_roi_ts[roi_lowd[0]].mean(1)[[0,-1]],scra_roi_ts[roi_lowd[2]].mean(1)[[0,-1]],'o',c='k',ms=10)\n",
    "\n",
    "for axis in [ax.w_xaxis, ax.w_yaxis, ax.w_zaxis]:\n",
    "    axis.line.set_linewidth(2)    \n",
    "ax.set_xticklabels([])\n",
    "ax.set_zticklabels([])\n",
    "ax.set_yticklabels([])\n",
    "fig.tight_layout()\n",
    "\n",
    "fig = plt.figure(figsize=(3,3))\n",
    "ax = fig.add_subplot(111,projection='3d')\n",
    "ax.plot(face_ch_ts[ch_lowd[1]].mean(1),face_ch_ts[ch_lowd[0]].mean(1),face_ch_ts[ch_lowd[2]].mean(1),'.-',c='gray',lw=lw)\n",
    "ax.plot(scra_ch_ts[ch_lowd[1]].mean(1),scra_ch_ts[ch_lowd[0]].mean(1),scra_ch_ts[ch_lowd[2]].mean(1),'.-',c='gray',lw=lw)\n",
    "ax.scatter(face_ch_ts[ch_lowd[1]].mean(1),face_ch_ts[ch_lowd[0]].mean(1),face_ch_ts[ch_lowd[2]].mean(1),s=s,c=cm.Blues(tvec[tvec_analysis]),edgecolor=ec,alpha=alpha,lw=0.4)\n",
    "ax.scatter(scra_ch_ts[ch_lowd[1]].mean(1),scra_ch_ts[ch_lowd[0]].mean(1),scra_ch_ts[ch_lowd[2]].mean(1),s=s,c=cm.Oranges(tvec[tvec_analysis]),edgecolor=ec,alpha=alpha,lw=0.4)\n",
    "ax.plot(scra_ch_ts[roi_lowd[1]].mean(1)[[0,-1]],scra_ch_ts[roi_lowd[0]].mean(1)[[0,-1]],scra_ch_ts[roi_lowd[2]].mean(1)[[0,-1]],'o',c='k',ms=10)\n",
    "ax.plot(scra_ch_ts[roi_lowd[1]].mean(1)[[0,-1]],scra_ch_ts[roi_lowd[0]].mean(1)[[0,-1]],scra_ch_ts[roi_lowd[2]].mean(1)[[0,-1]],'o',c='k',ms=10)\n",
    "\n",
    "for axis in [ax.w_xaxis, ax.w_yaxis, ax.w_zaxis]:\n",
    "    axis.line.set_linewidth(2)    \n",
    "ax.set_xticklabels([])\n",
    "ax.set_zticklabels([])\n",
    "ax.set_yticklabels([])\n",
    "fig.tight_layout()\n",
    "\n",
    "\n",
    "data = [corr_evoked,corr_graph,corr_pca,corr_ica]\n",
    "data_lowd = [corr_evoked_lowd,corr_graph_lowd,corr_pca_lowd,corr_ica_lowd]\n",
    "labels=['$r$','$\\lambda$','PCA','ICA']\n",
    "colors = [color1,color2,color3,color4]\n",
    "fig,axs = plt.subplots(2,1,figsize=(3,3))\n",
    "ax = axs[0]\n",
    "ax = my_box_plot(data_lowd,ax,colors,labels);\n",
    "ax.set_ylabel('Corr. (LD)')\n",
    "fancy_fig(ax)\n",
    "ax.set_ylim(-.1,1.1)\n",
    "ax.set_xticklabels([])\n",
    "\n",
    "ax = axs[1]\n",
    "ax = my_box_plot(data,ax,colors,labels);\n",
    "ax.set_ylabel('Corr. (HD)')\n",
    "fancy_fig(ax)\n",
    "ax.set_xlabel('Coordinate system')\n",
    "ax.set_ylim(-.1,1.1)\n",
    "\n",
    "\n",
    "fig,ax = plt.subplots(figsize=(4,3))\n",
    "\n",
    "ax.axvline(0,ls='--',c='k',lw=1)\n",
    "im = ax.pcolormesh(tvec[tvec_analysis],percentiles,np.median(distance,2),cmap='magma')\n",
    "ax.set_xlim((np.min(tvec[tvec_analysis]),np.max(tvec[tvec_analysis])))\n",
    "ax.set_xlabel('Time [ms]')\n",
    "ax.set_ylabel('Dimensionality [%]')\n",
    "fig.colorbar(im,ax=ax)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Map between scale 3 parcellation and Deskian atlas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "roidata_scale1 = pd.read_excel(roifname, sheet_name='SCALE 1')\n",
    "cort_scale1 = np.where(roidata_scale1['Structure'] == 'cort')[0]\n",
    "ROIs_scale1 = list(roidata_scale1['Label Lausanne2008'][cort_scale1])\n",
    "hemi_scale1 = list(roidata_scale1['Hemisphere'][cort_scale1])\n",
    "hemi_scale3 = list(roidata['Hemisphere'][cort])\n",
    "map_scale1_rh = []\n",
    "map_scale1_lh = []\n",
    "for j,roi in enumerate(ROIs_scale1):\n",
    "    if hemi_scale1[j] == 'rh':\n",
    "        map_scale1_rh.append([i for i,elem in enumerate(ROIs) if (roi == elem.split('_')[0]) & (hemi_scale3[i]=='rh')])\n",
    "    else:\n",
    "        map_scale1_lh.append([i for i,elem in enumerate(ROIs) if (roi == elem.split('_')[0]) & (hemi_scale3[i]=='lh')])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Select time point where difference is largest in the original high-dim space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "toi = np.argmax(np.linalg.norm((face_roi_ts-scra_roi_ts).mean(2),2,0))\n",
    "print(tvec[tvec_analysis][toi])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot for Fig. S6 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diff = np.zeros((sc.N,nsub_faces))\n",
    "diff[ch_lowd] = (face_ch_ts-scra_ch_ts)[ch_lowd,toi]\n",
    "diff = sc.igft(diff)\n",
    "\n",
    "diff_scale1_rh = np.array([np.sum(diff[rois],0) for rois in map_scale1_rh])\n",
    "ids_sort_rh = list(np.argsort(abs(diff_scale1_rh).mean(1))[::-1])\n",
    "labels_rh = [ROIs_scale1[i] for i in ids_sort_rh]\n",
    "boxes_rh = [abs(diff_scale1_rh[i]) for i in ids_sort_rh]\n",
    "diff_scale1_lh = np.array([np.sum(diff[rois],0) for rois in map_scale1_lh])\n",
    "ids_sort_lh = list(np.argsort(abs(diff_scale1_lh).mean(1))[::-1])\n",
    "labels_lh = [ROIs_scale1[i] for i in ids_sort_lh]\n",
    "boxes_lh = [abs(diff_scale1_lh[i]) for i in ids_sort_lh]\n",
    "\n",
    "\n",
    "\n",
    "fig,axs = plt.subplots(1,2,figsize=(3,10))\n",
    "axs[0].boxplot(boxes_lh,labels=labels_lh,vert=False,showfliers=False);\n",
    "axs[0].plot([np.mean(elem) for elem in boxes_lh],np.arange(1,len(boxes_lh)+1));\n",
    "axs[0].set_title('LH')\n",
    "axs[1].boxplot(boxes_rh,labels=labels_rh,vert=False,showfliers=False);\n",
    "axs[1].plot([np.mean(elem) for elem in boxes_rh],np.arange(1,len(boxes_rh)+1));\n",
    "axs[1].set_title('RH')\n",
    "axs[1].invert_xaxis()\n",
    "axs[1].yaxis.tick_right()\n",
    "axs[0].tick_params(width=3, length=8)\n",
    "axs[0].spines['top'].set_linewidth(0)\n",
    "axs[0].spines['bottom'].set_linewidth(3)\n",
    "axs[0].spines['left'].set_linewidth(3)\n",
    "axs[0].spines['right'].set_linewidth(0)\n",
    "\n",
    "axs[1].tick_params(width=3, length=8)\n",
    "axs[1].spines['top'].set_linewidth(0)\n",
    "axs[1].spines['bottom'].set_linewidth(3)\n",
    "axs[1].spines['left'].set_linewidth(0)\n",
    "axs[1].spines['right'].set_linewidth(3)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plots Figure 5\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import gudhi as gd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def group_agreement(norm,percentage):\n",
    "    # Initialize dims with nans\n",
    "    dims = np.zeros_like(norm)*np.nan\n",
    "    for s in range(norm.shape[0]):\n",
    "        # write down all those dims with norm larger than percentile. \n",
    "        dims_s = np.where(norm[s]>np.percentile(norm[s],percentage))[0]\n",
    "        dims[s,:len(dims_s)] = dims_s    \n",
    "    \n",
    "    # Set a counter that counts how many subjects have kept each dimensions kept before. \n",
    "    unique_dims = np.unique(dims[~np.isnan(dims)])\n",
    "    counter = np.zeros_like(unique_dims)\n",
    "    for i in range(len(unique_dims)):\n",
    "        counter[i] = len(np.where(dims == unique_dims[i])[0])\n",
    "    \n",
    "    # In dims_agreement we store only those dims kept in all subjects\n",
    "    dims_agreement = unique_dims[np.where(counter==norm.shape[0])]\n",
    "    # In agreement we count how many dims are kept for all subjects\n",
    "    agreement = sum(counter==norm.shape[0])\n",
    "    \n",
    "    return agreement, dims_agreement.astype('int')\n",
    "\n",
    "# Initialize norms and evoked matrices\n",
    "norms_roi_fro = np.zeros((nsub_faces,sc.N))\n",
    "norms_graph_fro = np.zeros((nsub_faces,sc.N))\n",
    "\n",
    "evoked_f = np.zeros((nsub_faces,sc.N,len(tvec_analysis)))\n",
    "evoked_s = np.zeros((nsub_faces,sc.N,len(tvec_analysis)))\n",
    "evoked_graph_f = np.zeros((nsub_faces,sc.N,len(tvec_analysis)))\n",
    "evoked_graph_s = np.zeros((nsub_faces,sc.N,len(tvec_analysis)))\n",
    "\n",
    "\n",
    "for s in tqdm.tqdm(range(nsub_faces)):    \n",
    "    # Load time-series for subject [ROIs x Time x Trials]\n",
    "    epochs,cond = ts_gen_faces.loader_ts(s)\n",
    "    # Baseline correction\n",
    "    epochs = epochs - epochs[:,tvec_pre].mean(1,keepdims=True)\n",
    "    \n",
    "    # Visual Evoked Potentials\n",
    "    evoked_f[s] = epochs[:,:,cond==1].mean(2)\n",
    "    evoked_s[s] = epochs[:,:,cond==0].mean(2)\n",
    "\n",
    "    # Graph Fourier transform\n",
    "    epochs_graph = sc.gft(epochs)\n",
    "    evoked_graph_f[s] = epochs_graph[:,:,cond==1].mean(2)\n",
    "    evoked_graph_s[s] = epochs_graph[:,:,cond==0].mean(2)\n",
    "    # Frobernius norm of each dimension\n",
    "    norms_roi_fro[s] = np.linalg.norm(epochs[:,:,cond==1],ord='fro',axis=(1,2))    \n",
    "    norms_graph_fro[s] = np.linalg.norm(epochs_graph[:,:,cond==1],ord='fro',axis=(1,2))\n",
    "\n",
    "# Sweep across norm space and select the dimensions with large norms present in all subjects\n",
    "percentiles3 = np.linspace(0,100,10000)\n",
    "nrois =[]\n",
    "nchs =[]\n",
    "for p,percentile in tqdm.tqdm(enumerate(percentiles3)):\n",
    "    # b includes those dimensions with norm larger than percentile that appear in all subjects.\n",
    "    _,b = group_agreement(norms_roi_fro,percentile)\n",
    "    _,b2 = group_agreement(norms_graph_fro,percentile)\n",
    "    # we store the number of dimensions -not the dimensions- that are kept in all subjects \n",
    "    nrois.append(len(b))\n",
    "    nchs.append(len(b2))\n",
    "\n",
    "    \n",
    "# See which number of kept rois and number of kept CH are equivalent, to do further analysis with the same number of dims.\n",
    "dimensions_analysis = [elem for elem in np.unique(nchs) if elem in np.unique(nrois)]\n",
    "dimensions_analysis.remove(0)\n",
    "rois_analysis = []\n",
    "ch_analysis = []\n",
    "dimensions_analysis_roi = dimensions_analysis.copy()\n",
    "dimensions_analysis_ch = dimensions_analysis.copy()\n",
    "for p,percentile in tqdm.tqdm(enumerate(percentiles3)):\n",
    "    # look again for kept dimensions\n",
    "    _,b = group_agreement(norms_roi_fro,percentile)\n",
    "    _,b2 = group_agreement(norms_graph_fro,percentile)\n",
    "    # if number of kept dimensions is the list of wanted dimensions, keep them. \n",
    "    if len(b) in dimensions_analysis_roi:\n",
    "        rois_analysis.append(b)\n",
    "        dimensions_analysis_roi.remove(len(b))\n",
    "    if len(b2) in dimensions_analysis_ch:\n",
    "        ch_analysis.append(b2)\n",
    "        dimensions_analysis_ch.remove(len(b2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute bottleneck distances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.spatial.distance import pdist,squareform\n",
    "\n",
    "if not load_data:\n",
    "    # Initialize bottleneck distances between high-dim and low-dim spaces\n",
    "    b0_ch_f = np.zeros((nsub_faces,len(dimensions_analysis)))\n",
    "    b0_ch_s = np.zeros((nsub_faces,len(dimensions_analysis)))\n",
    "    b0_rois_f = np.zeros((nsub_faces,len(dimensions_analysis)))\n",
    "    b0_rois_s = np.zeros((nsub_faces,len(dimensions_analysis)))\n",
    "\n",
    "    b1_ch_f = np.zeros((nsub_faces,len(dimensions_analysis)))\n",
    "    b1_ch_s = np.zeros((nsub_faces,len(dimensions_analysis)))\n",
    "    b1_rois_f = np.zeros((nsub_faces,len(dimensions_analysis)))\n",
    "    b1_rois_s = np.zeros((nsub_faces,len(dimensions_analysis)))\n",
    "\n",
    "    for s in tqdm.tqdm(range(nsub_faces)):    \n",
    "        # Compute distance matrix\n",
    "        dist_mat_f = squareform(pdist(face_roi_ts[:,:,s].T))\n",
    "        dist_mat_s = squareform(pdist(scra_roi_ts[:,:,s].T))\n",
    "\n",
    "        # Build the Vietoris-Rips complex        \n",
    "        # First build the skeleton-graph (MST?)\n",
    "        skeleton_f = gd.RipsComplex(distance_matrix = dist_mat_f, max_edge_length = np.max(dist_mat_f))\n",
    "        skeleton_s = gd.RipsComplex(distance_matrix = dist_mat_s, max_edge_length = np.max(dist_mat_s))\n",
    "\n",
    "        # Then build the rips simplicial complex\n",
    "        rips_simplex_tree_f = skeleton_f.create_simplex_tree(max_dimension=2)\n",
    "        rips_simplex_tree_s = skeleton_s.create_simplex_tree(max_dimension=2)\n",
    "\n",
    "        # Now compute persistence on the simplex tree\n",
    "        barcodes_rips_f = rips_simplex_tree_f.persistence()\n",
    "        barcodes_rips_s = rips_simplex_tree_s.persistence()\n",
    "        persistence_list0_f = rips_simplex_tree_f.persistence_intervals_in_dimension(0)\n",
    "        persistence_list1_f = rips_simplex_tree_f.persistence_intervals_in_dimension(1)\n",
    "        persistence_list0_s = rips_simplex_tree_s.persistence_intervals_in_dimension(0)\n",
    "        persistence_list1_s = rips_simplex_tree_s.persistence_intervals_in_dimension(1)\n",
    "\n",
    "        # Compute dim reduction and bottleneck distances\n",
    "        for d in range(len(dimensions_analysis)):    \n",
    "            ch = ch_analysis[d]\n",
    "            rois = rois_analysis[d]\n",
    "\n",
    "            # Compute distance matrix\n",
    "\n",
    "            dist_mat_ch_f = squareform(pdist(face_ch_ts[ch,:,s].T))\n",
    "            dist_mat_ch_s = squareform(pdist(scra_ch_ts[ch,:,s].T))\n",
    "            dist_mat_rois_f = squareform(pdist(face_roi_ts[rois,:,s].T))\n",
    "            dist_mat_rois_s = squareform(pdist(scra_roi_ts[rois,:,s].T))\n",
    "\n",
    "            # Build the Vietoris-Rips complex        \n",
    "            # First build the skeleton-graph (MST?)\n",
    "            skeleton_ch_f = gd.RipsComplex(distance_matrix = dist_mat_ch_f, max_edge_length = np.max(dist_mat_ch_f))\n",
    "            skeleton_ch_s = gd.RipsComplex(distance_matrix = dist_mat_ch_s, max_edge_length = np.max(dist_mat_ch_s))\n",
    "            skeleton_rois_f = gd.RipsComplex(distance_matrix = dist_mat_rois_f, max_edge_length = np.max(dist_mat_rois_f))\n",
    "            skeleton_rois_s = gd.RipsComplex(distance_matrix = dist_mat_rois_s, max_edge_length = np.max(dist_mat_rois_s))\n",
    "            # Then build the rips simplicial complex\n",
    "            rips_simplex_tree_ch_f = skeleton_ch_f.create_simplex_tree(max_dimension=2)\n",
    "            rips_simplex_tree_ch_s = skeleton_ch_s.create_simplex_tree(max_dimension=2)\n",
    "            rips_simplex_tree_rois_f = skeleton_rois_f.create_simplex_tree(max_dimension=2)\n",
    "            rips_simplex_tree_rois_s = skeleton_rois_s.create_simplex_tree(max_dimension=2)\n",
    "\n",
    "            # Now compute persistence on the simplex tree\n",
    "\n",
    "            barcodes_rips_ch_f = rips_simplex_tree_ch_f.persistence()\n",
    "            persistence_list0_ch_f = rips_simplex_tree_ch_f.persistence_intervals_in_dimension(0)\n",
    "            persistence_list1_ch_f = rips_simplex_tree_ch_f.persistence_intervals_in_dimension(1)\n",
    "\n",
    "            barcodes_rips_ch_s = rips_simplex_tree_ch_s.persistence()\n",
    "            persistence_list0_ch_s = rips_simplex_tree_ch_s.persistence_intervals_in_dimension(0)\n",
    "            persistence_list1_ch_s = rips_simplex_tree_ch_s.persistence_intervals_in_dimension(1)\n",
    "\n",
    "            barcodes_rips_rois_f = rips_simplex_tree_rois_f.persistence()\n",
    "            persistence_list0_rois_f = rips_simplex_tree_rois_f.persistence_intervals_in_dimension(0)\n",
    "            persistence_list1_rois_f = rips_simplex_tree_rois_f.persistence_intervals_in_dimension(1)\n",
    "\n",
    "            barcodes_rips_rois_s = rips_simplex_tree_rois_s.persistence()\n",
    "            persistence_list0_rois_s = rips_simplex_tree_rois_s.persistence_intervals_in_dimension(0)\n",
    "            persistence_list1_rois_s = rips_simplex_tree_rois_s.persistence_intervals_in_dimension(1)\n",
    "\n",
    "            b0_ch_f[s,d] = gd.bottleneck_distance(persistence_list0_f,persistence_list0_ch_f)\n",
    "            b0_ch_s[s,d] = gd.bottleneck_distance(persistence_list0_s,persistence_list0_ch_s)\n",
    "            b0_rois_f[s,d] = gd.bottleneck_distance(persistence_list0_f,persistence_list0_rois_f)\n",
    "            b0_rois_s[s,d] = gd.bottleneck_distance(persistence_list0_s,persistence_list0_rois_s)\n",
    "\n",
    "            b1_ch_f[s,d] = gd.bottleneck_distance(persistence_list1_f,persistence_list1_ch_f)\n",
    "            b1_ch_s[s,d] = gd.bottleneck_distance(persistence_list1_s,persistence_list1_ch_s)\n",
    "            b1_rois_f[s,d] = gd.bottleneck_distance(persistence_list1_f,persistence_list1_rois_f)\n",
    "            b1_rois_s[s,d] = gd.bottleneck_distance(persistence_list1_s,persistence_list1_rois_s)\n",
    "            \n",
    "    np.save(os.path.join(topology_path,'b0_ch_f.npy'),b0_ch_f)\n",
    "    np.save(os.path.join(topology_path,'b0_ch_s.npy'),b0_ch_s)\n",
    "    np.save(os.path.join(topology_path,'b0_rois_f.npy'),b0_rois_f)\n",
    "    np.save(os.path.join(topology_path,'b0_rois_s.npy'),b0_rois_s)\n",
    "    np.save(os.path.join(topology_path,'b1_ch_f.npy'),b1_ch_f)\n",
    "    np.save(os.path.join(topology_path,'b1_ch_s.npy'),b1_ch_s)\n",
    "    np.save(os.path.join(topology_path,'b1_rois_f.npy'),b1_rois_f)\n",
    "    np.save(os.path.join(topology_path,'b1_rois_s.npy'),b1_rois_s)\n",
    "\n",
    "\n",
    "\n",
    "else:\n",
    "    \n",
    "    b0_ch_f = np.load(os.path.join(topology_path,'b0_ch_f.npy'))\n",
    "    b0_ch_s = np.load(os.path.join(topology_path,'b0_ch_s.npy'))\n",
    "    b0_rois_f = np.load(os.path.join(topology_path,'b0_rois_f.npy'))\n",
    "    b0_rois_s = np.load(os.path.join(topology_path,'b0_rois_s.npy'))\n",
    "    b1_ch_f = np.load(os.path.join(topology_path,'b1_ch_f.npy'))\n",
    "    b1_ch_s = np.load(os.path.join(topology_path,'b1_ch_s.npy'))\n",
    "    b1_rois_f = np.load(os.path.join(topology_path,'b1_rois_f.npy'))\n",
    "    b1_rois_s = np.load(os.path.join(topology_path,'b1_rois_s.npy'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare statistically"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_,pvalue_b0f = ranksums(b0_rois_f.flatten(),b0_ch_f.flatten())\n",
    "_,pvalue_b1f = ranksums(b1_rois_f.flatten(),b1_ch_f.flatten())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(pvalue_b0f*2)\n",
    "print(pvalue_b1f*2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plots for Fig. 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fsize = (3.5*1.4,3*1.4)\n",
    "font = {'size'   : 16}\n",
    "matplotlib.rc('font', **font)\n",
    "\n",
    "nsample = len(tvec_analysis)//4\n",
    "epochs_plot = face_ch_ts.mean(2)[ch_lowd[:2]][:,::4]\n",
    "epochs_plot = epochs_plot - np.min(epochs_plot,axis=1,keepdims=True)\n",
    "epochs_plot = epochs_plot / np.max(epochs_plot,axis=1,keepdims=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(15,5))\n",
    "\n",
    "thres = [0.075,0.10,0.19]\n",
    "s = [260,470,1500]\n",
    "\n",
    "for kk in range(3):\n",
    "\n",
    "    distances = squareform(pdist(np.array([epochs_plot[0],epochs_plot[1]]).T))\n",
    "    dist_ = distances < thres[kk]\n",
    "\n",
    "    edges = []\n",
    "    for i in range(nsample):\n",
    "        for j in range(i,nsample):\n",
    "            if dist_[i,j]:\n",
    "                edges.append([[epochs_plot[k][i],epochs_plot[k][j]] for k in range(2)])\n",
    "\n",
    "\n",
    "\n",
    "    ax = fig.add_subplot(1,3,kk+1)\n",
    "    ax.plot(epochs_plot[1],epochs_plot[0],'.',c='k',ms=10)\n",
    "    ax.scatter(epochs_plot[1],epochs_plot[0],s=s[kk],c='k',edgecolor='k',alpha=0.1,lw=0.5)\n",
    "    for edge in edges:\n",
    "        ax.plot(edge[1],edge[0],c='k',lw=1,alpha=0.5)\n",
    "    ax.set_xticks([])\n",
    "    ax.set_yticks([])\n",
    "    ax.set_xlim([-0.1,1.1])\n",
    "    ax.set_ylim([-0.1,1.1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_axis = np.array([len(elem) for elem in ch_analysis])\n",
    "fig,ax = plt.subplots(figsize=fsize)\n",
    "\n",
    "plot_shaded_norm(x_axis,b0_rois_f,0,color1,label='$r$',ax = ax)\n",
    "plot_shaded_norm(x_axis,b0_ch_f,0,color2,label='$\\lambda$',ax = ax)\n",
    "ax.legend(frameon=False)\n",
    "\n",
    "ax.set_xlabel('Dimensions considered')\n",
    "ax.set_ylabel('Bottleneck distance')\n",
    "fancy_fig(ax)\n",
    "ax.set_xlim([0,np.max(x_axis)])\n",
    "ax.set_ylim([0,6])\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,ax = plt.subplots(figsize=fsize)\n",
    "\n",
    "plot_shaded_norm(x_axis,b1_rois_f,0,color1,label='$r$',ax = ax)\n",
    "plot_shaded_norm(x_axis,b1_ch_f,0,color2,label='$\\lambda$',ax = ax)\n",
    "\n",
    "ax.legend(frameon=False)\n",
    "\n",
    "ax.set_xlabel('Dimensions considered')\n",
    "ax.set_ylabel('Bottleneck distance')\n",
    "fancy_fig(ax)\n",
    "ax.set_xlim([0,np.max(x_axis)])\n",
    "ax.set_ylim([0,6])\n",
    "fig.tight_layout()\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}