{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import scipy.stats as sps\n",
    "from scipy.io import loadmat\n",
    "import scipy.signal as sig\n",
    "import matplotlib.pyplot as pl\n",
    "from mpl_toolkits.axes_grid1.inset_locator import inset_axes\n",
    "%matplotlib tk\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from matplotlib import cm\n",
    "from matplotlib.patches import Circle\n",
    "from matplotlib.gridspec import GridSpec\n",
    "import copy\n",
    "%matplotlib tk\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Resting state dynamics, SI fig 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pathFile = './data/rest/Globals_v5ms_grid_globalC_'\n",
    "globalC = np.concatenate((np.array([0.01, 0.05, 0.1, 0.15,0.2,0.25]), np.arange(0.3,2.,0.05)))\n",
    "\n",
    "nL = 13\n",
    "nG = globalC.shape[0]\n",
    "nMo = 25\n",
    "T = 60000\n",
    "\n",
    "mG = np.empty((13,globalC.shape[0]))\n",
    "mL = np.empty((13,globalC.shape[0]))\n",
    "mGt = np.empty((13,globalC.shape[0], T, nMo), dtype=np.float32)\n",
    "mL[:] = np.nan\n",
    "mG[:] = np.nan\n",
    "\n",
    "for iX, lC in enumerate(globalC):\n",
    "    file = pathFile + str(np.round(lC,3)) + '_co.npz'\n",
    "    with np.load(file) as data:\n",
    "        mG[:,iX] = data['gMet'].ravel()\n",
    "        mL[:,iX] = data['lMet'].ravel()\n",
    "        mGt[:,iX,:,:] = np.squeeze(data['gMetT'])\n",
    "        localC  = data['localC']\n",
    "\n",
    "mGc = np.mean(np.std(mGt,axis=2),axis=2)\n",
    "\n",
    "mL[np.isnan(mL)] = np.nanmin(mL)\n",
    "mG[np.isnan(mG)] = np.nanmin(mG)\n",
    "mGc[np.isnan(mGc)] = np.nanmin(mGc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "autocorrelation of global metatability"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "acMet = np.zeros((nL, nG))\n",
    "\n",
    "for i, iXl in enumerate(range(localC.shape[0])):\n",
    "    for iXg in range(globalC.shape[0]):\n",
    "        acM = 0\n",
    "        for iXm in range(nMo): # number of models\n",
    "#             autocorrelation\n",
    "            y = mGt[iXl, iXg, :, iXm] - np.mean(mGt[iXl, iXg, :, iXm])\n",
    "            ac_ = sig.correlate(y, y, mode='full')\n",
    "            ac_ /= ac_.max()\n",
    "            acM += np.mean(ac_[:int(ac_.shape[0]/2)]**2)\n",
    "        acMet[i, iXg] = acM / nMo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "SI fig 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lC = 5\n",
    "gC = 3\n",
    "fontS=9\n",
    "figlabels = ['Local Metasta.', 'Global Metasta.', 'Autocorrelation\\nGlobal Order']\n",
    "im = [None]*3\n",
    "fgrid = dict(hspace=0.1,top=0.9,bottom=0.2, left=0.1, right=0.9)\n",
    "fig, axs = pl.subplots(3,1,gridspec_kw=fgrid)\n",
    "\n",
    "fig.set_figheight(7)\n",
    "fig.set_figwidth(7)\n",
    "fig.set_dpi(900)\n",
    "\n",
    "im[0] = axs[0].pcolor(mL[1:,1:], cmap=cm.magma, vmin=0, vmax=np.nanmax(mLm[1:,1:]))\n",
    "im[1] = axs[1].pcolor(mGc[1:,1:], cmap=cm.magma, vmin=0, vmax=np.nanmax(mGcm[1:,1:]))\n",
    "im[2] = axs[2].pcolor(acMet[1:,1:], cmap=cm.magma_r,  vmin=0, vmax=np.nanmax(acMet[1:,1:]))\n",
    "\n",
    "axs[0].set_title('Resting state coordination dynamics')\n",
    "for i in range(3):\n",
    "    axc = pl.colorbar(im[i],ax=axs[i],ticks=[0,np.nanmax(im[i].get_array())])\n",
    "    axc.set_label(figlabels[i])\n",
    "    axc.ax.yaxis.tick_left()\n",
    "    axc.ax.set_yticklabels(['0',\n",
    "                            str(np.round(np.nanmax(im[i].get_array()),4))],rotation=90)\n",
    "    axc.ax.get_ymajorticklabels()[0].set_verticalalignment('bottom')  \n",
    "    axc.ax.get_ymajorticklabels()[1].set_verticalalignment('top') \n",
    "    axs[i].yaxis.set_ticks_position('both')\n",
    "    axs[i].set_yticks(np.arange(localC.shape[0]-1,None,2)+0.5)\n",
    "    axs[i].set_yticklabels(np.round(localC[1::2],2))\n",
    "    axs[i].set_aspect('equal')\n",
    "    axs[i].set_ylabel('Local\\nCoupling',fontsize=fontS)\n",
    "    circ = Circle((lC+0.5,gC+0.5), radius=0.5, facecolor='None', edgecolor='g', lw=2)\n",
    "    axs[i].add_patch(circ)\n",
    "    axs[i].set_xticklabels([''])\n",
    "    axs[i].set_xticks(np.arange(globalC.shape[0]-1,None,2)+0.5)\n",
    "    axs[i].xaxis.set_ticks_position('both')\n",
    "axs[2].set_xticklabels(np.round(globalC[1::2],2),rotation='vertical')\n",
    "    \n",
    "axs[-1].set_xlabel('Global Coupling',fontsize=fontS)\n",
    "axs[-1].patch.set(hatch='xx', edgecolor='k')\n",
    "axs[-1].patch.set(color='k', fill=True,hatch='xx',edgecolor='w')\n",
    "\n",
    "# fname  = './dataSave/FigureSI_Rest.png'  \n",
    "# pl.savefig(fname, dpi=None, facecolor='w', edgecolor=None,\n",
    "#         orientation='portrait', papertype=None, format='png',\n",
    "#         transparent=False, bbox_inches=None, pad_inches=0.,\n",
    "#         frameon=None)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Evolution of optimizers, SI Fig 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with np.load('./data/evolutionIter.npz') as data:\n",
    "    scor = data['sCor']\n",
    "    fit  = data['fit']\n",
    "    ener = data['ener']\n",
    "    puls = data['puls']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stg = 3 # Stage to be plotted\n",
    "isl = 0 # I sland to be plotted\n",
    "\n",
    "fitMask= fit[isl][stg] > 0\n",
    "fit0   = fit[isl][stg][fitMask]\n",
    "scor0  = scor[isl][stg][fitMask]\n",
    "ener0  = ener[isl][stg][fitMask]\n",
    "puls0  = puls[isl][stg][fitMask,:]\n",
    "iXsort = np.argsort(fit0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "SI Fig 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ylabels = ['Fit\\nIndex', 'Spear.\\n Corr.','Total \\n input' , 'Local \\n Input','Local \\n Input Sorted']\n",
    "gridHratios = [0.13,0.13,0.13,0.3,0.3]\n",
    "fgrid = dict(height_ratios=gridHratios,\n",
    "             hspace=0.05,right=0.88)\n",
    "fig, axs = pl.subplots(5, 1, gridspec_kw= fgrid)\n",
    "axs[0].set_title('Evolution: Island 1, Stage 4')\n",
    "axs[0].plot(fit0,'.', markersize=0.1)\n",
    "axs[0].plot(fit0[iXsort],'.', markersize=0.1)\n",
    "# axs[0].legend(['iterations', 'sorted'])\n",
    "axs[1].plot(scor0,'.', markersize=0.1)\n",
    "axs[1].plot(scor0[iXsort],'.', markersize=0.05)\n",
    "axs[2].plot(ener0,'.', markersize=0.1)\n",
    "axs[2].plot(ener0[iXsort],'.', markersize=0.05)\n",
    "axs[2].set_ylim([0,60000])\n",
    "axs[2].set_yticks([0,30000, 60000])\n",
    "axs[3].matshow(puls0.T,aspect='auto',cmap='bwr',vmin=-1500,vmax=1500)\n",
    "im = axs[4].matshow(puls0[iXsort,:].T,aspect='auto',cmap='bwr',vmin=-1500,vmax=1500)\n",
    "axs[4].xaxis.tick_bottom()\n",
    "axs[4].set_xlabel('iterations')\n",
    "\n",
    "\n",
    "\n",
    "axc = fig.add_axes((.89,0.14,0.01,0.4))\n",
    "pl.colorbar(im, cax=axc)\n",
    "\n",
    "yLoc = (np.cumsum(np.flipud(gridHratios))) *0.8 +0.06# * (1-bottOff) + bottOff\n",
    "ids  = ['(E)','(D)','(C)','(B)','(A)']\n",
    "for i in range(5):\n",
    "    axs[i].set_ylabel(ylabels[i])\n",
    "    axs[i].yaxis.tick_right()\n",
    "    axs[i].set_xlim([0,fit0.shape[0]])\n",
    "    fig.text(0.07, yLoc[i],ids[i],ha='left',va='center')\n",
    "    if i < 2:\n",
    "        axs[i].set_ylim([0,1])\n",
    "        axs[i].yaxis.get_major_ticks()[1].label2.set_verticalalignment('top')  \n",
    "        axs[i].yaxis.get_major_ticks()[0].label2.set_verticalalignment('bottom')  \n",
    "    if i == 2:\n",
    "        axs[i].yaxis.get_major_ticks()[-1].label2.set_verticalalignment('top')  \n",
    "        axs[i].yaxis.get_major_ticks()[0].label2.set_verticalalignment('bottom')  \n",
    "    if i < 4:\n",
    "        axs[i].set_xticks([])\n",
    "    if i > 2:\n",
    "        axs[i].set_yticks([])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Return to metastability, SI Fig 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file = './data/metaFC1000afterPulse_G0p3L0p7.npz'\n",
    "\n",
    "with np.load(file) as data:\n",
    "    siMglo = data['metaG']\n",
    "    siMloc = data['metaL']\n",
    "    si2fc  = data['sim2eFC']\n",
    "\n",
    "\n",
    "fs = 100\n",
    "t  = np.linspace(1/fs,siMglo.shape[1]/fs,siMglo.shape[1])\n",
    "tR = np.arange(1/fs,1.5+1/fs,1/fs) # trial length\n",
    "nR = tR.shape[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "SI Fig 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ylabels = ['Global Meta.','Local Meta.','fit to FC']\n",
    "fgrid = dict(height_ratios=[0.33,0.33,0.33],\n",
    "             hspace=0.1)\n",
    "fig, axs = pl.subplots(3, 1, gridspec_kw= fgrid)\n",
    "\n",
    "\n",
    "axs[0].plot(t, siMglo[-1,:])\n",
    "p1 = axs[1].plot(t, siMloc[-1,:,:,0].T, linewidth=0.5,alpha=0.7)\n",
    "axs[2].plot(t, si2fc.T)\n",
    "for i in range(3):\n",
    "    axs[i].set_ylabel(ylabels[i])\n",
    "    if i < 2:\n",
    "        axs[i].set_xticks([])\n",
    "        axs[i].set_ylim([0,1])\n",
    "    else:\n",
    "        axs[i].set_xlabel('Time [seconds]')\n",
    "        axs[i].set_ylim([-0.5,1])\n",
    "        axs[i].legend(loc=1,labels=['stage 1','stage 2','stage 3','stage 4'],\n",
    "           labelspacing=0.1,framealpha=0.6,edgecolor='w',handletextpad=0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}