{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\nCreated on Mon Nov 22 21:19:37 2021\n\n@author: maout\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\n\n\n#from matplotlib import pyplot as plt\n\nfrom scipy.spatial.distance import cdist\nimport time\nimport ot\nimport numba\nimport math\nimport random\nimport sys\nimport pickle\n\ndim = 6\nh = 0.001\nt1 = 0\nt2 = 0.5\nT = t2-t1\ntimegrid = np.arange(0,T+h/2,h)\nN = 3000#0#0#0#2000\n#g = 1\nk = timegrid.size\nM = 300\n\ny2 = np.ones(dim)  \nsigmas = np.array([1.0,1.5])\nrep_bridge = 1   #10 different bridge instances for every setting\nreps = 20 ##instanses for stochastic path evaluation of each bridge\n# dy0 = np.array([np.pi/4, np.pi/2, np.pi, 3*np.pi/2])    \n# y0s = np.zeros((dim,rep_bridge, dy0.size))\nrandom.seed(22)\n# for i in range(rep_bridge):\n#     y0s[0, i, :] = np.random.uniform( low=0, high=2*np.pi,size=1 ) \n#     y0s[1, i, :] = (y0s[0, i, :] + dy0) %(2* np.pi)\n\n\n\nKs   = np.linspace(0,6,7) \nrepetition = 0\nRttcont = np.zeros((Ks.size,  sigmas.size, timegrid.size,reps   ))*np.nan\nRttnon = np.zeros((  Ks.size,  sigmas.size, timegrid.size,reps   ))*np.nan\nused_us = np.zeros((dim,  Ks.size,  sigmas.size,timegrid.size,  reps   ))*np.nan    \nfor gii in range(0,2):\n    noise = gii+1\n    for ki in range(Ks.size):   \n        K = ki#Ks[ki]\n        naming = 'Kuramoto_6N\\\\dim6\\\\Latest%d_N_systematic_Kuramoto_inhomogeneous_k_%d_gi_%d_N_%d_M_%d_repetition_%d'%(dim, K,noise,N, M,repetition)\n          \n        try:\n            file = open(naming+'.dat','rb')\n            to_save = pickle.load(file)   \n            # naming2 = 'kuramot\\\\2N_systematic_Kuramoto_homogeneous_repb_%d_ki_%d_yi_%d__gi_%d_N_%d_M_%d_UNCONTROLLED'%(bi, ki,yi,gii,N, M)\n            # file2 = open(naming2+'.dat','rb')\n            # to_save2 = pickle.load(file2)\n            \n            #Fcont = to_save['Fcont'] \n            Rttcont[ki,gii,:,:] = to_save['Rttcont']  \n            #Fnon = to_save2['Fnon'] \n            Rttnon[ki, gii,:,:] = to_save['Rttnon']  \n            #Kin = to_save['K']                     \n                            \n            used_us[:,ki, gii,:] = to_save['used_us']   \n            \n        except FileNotFoundError:\n            i=0\n            print(naming)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from matplotlib import pyplot as plt                     \n\n#figure and plotting settings\nimport seaborn as sns\nplt.rcParams['figure.dpi'] = 300\nplt.rcParams['savefig.dpi'] = 300\nplt.rcParams['savefig.facecolor'] = (1,1,1,0)\nplt.rcParams['savefig.bbox'] = \"tight\"\nplt.rcParams['font.size'] = 10*1.2\n# plt.rcParams['font.family'] = 'sans-serif'     # not available in Colab\n# plt.rcParams['font.sans-serif'] = 'Helvetica'  # not available in Colab\nplt.rcParams['pdf.fonttype'] = 42\nplt.rcParams['xtick.labelsize'] = 10*1.2\nplt.rcParams['ytick.labelsize'] = 10*1.2\nplt.rcParams['axes.labelsize'] = 12*1.2\nplt.rcParams['axes.linewidth'] = 2\nplt.rcParams['axes.spines.top'] = False\nplt.rcParams['axes.spines.right'] = False\n#plt.rcParams['axes.Axes.tick_params'] = True\n#plt.rcParams['axes.solid_capstyle'] = 'round'\nplt.rc('axes',edgecolor='#4f4949')\nplt.rcParams['figure.frameon'] = False\nplt.rcParams['figure.subplot.hspace'] = 0.51\nplt.rcParams['figure.subplot.wspace'] = 0.51\nplt.rcParams['figure.subplot.left'] = 0.1\nplt.rcParams['figure.subplot.right'] = 0.9\nplt.rcParams['figure.subplot.top'] = 0.9\nplt.rcParams['figure.subplot.bottom'] = 0.1\nplt.rcParams['lines.solid_capstyle'] = 'round'\nplt.rcParams['lines.solid_joinstyle'] = 'round'\n#plt.rcParams['xtick.major.size'] = 20\n#plt.rcParams['xtick.major.width'] = 4\nplt.rcParams['text.usetex'] = True"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "onset_cont = np.zeros((Ks.size,2,20)) *np.nan\nonset_uncont = np.zeros((Ks.size,2,20)) *np.nan\n\ncont_synched_durations_m = np.zeros((Ks.size,2,20)) *np.nan\ncont_synched_durations_std  = np.zeros((Ks.size,2,20)) *np.nan\ncont_synched_durations_max = np.zeros((Ks.size,2,20)) *np.nan\ncont_synched_durations_sum  = np.zeros((Ks.size,2,20)) *np.nan\n\nuncont_synched_durations_m = np.zeros((Ks.size,2,20)) *np.nan\nuncont_synched_durations_std  = np.zeros((Ks.size,2,20)) *np.nan\nuncont_synched_durations_max = np.zeros((Ks.size,2,20)) *np.nan\nuncont_synched_durations_sum  = np.zeros((Ks.size,2,20)) *np.nan\n\ncounter_cont = np.zeros((Ks.size,2))\ncounter_uncont =  np.zeros((Ks.size,2))\n\nthreshold = 0.95\nmin_duri = 20\n\nfor ki in range(Ks.size):\n    \n    for gi in range(2):#3                \n        for repi in range(20):                    \n            positions = np.where( Rttcont[ ki,gi,1:-2, repi] >=threshold)[0]  ###this one detects if there is at all any synchronisation\n            synch_or_not = Rttcont[ ki, gi,1:-2, repi] >=threshold ##this is boolean indicating the synchronised positions\n            \n            if positions.size==0: \n                onset_cont[ki,gi,repi] = np.nan\n            else:\n                diffpos = np.diff(synch_or_not.astype(int)   ) #difference between consecutive steps\n                \n                diffpos2 = np.zeros(diffpos.size+1) # extend by one step at the beginning to detect \n                diffpos2[1:] = diffpos\n                starts = np.argwhere(diffpos2 == 1)\n                stops = np.argwhere(diffpos2 == -1)                        \n                if starts.size == 0:\n                    onset_cont[ki,gi,repi] = np.nan\n                elif stops.size == 0:\n                    onset_cont[ki,gi,repi] = timegrid[ starts[0] ]\n                    counter_cont[ki,gi] += 1\n                    dur = timegrid[ -1 ] - timegrid[ starts[0] ]\n                    cont_synched_durations_m[ki,gi,repi]  = dur\n                    cont_synched_durations_std[ki,gi,repi]  = 0\n                    cont_synched_durations_max[ki,gi,repi]  = dur\n                    cont_synched_durations_sum[ki,gi,repi]  = np.sum(  synch_or_not )/ (timegrid.size-1 - starts[0] )\n                else:\n                    if stops[0,0] < starts[0,0]:\n                        stops = stops[1:,:]\n                    durations = stops[:,0] - starts[:stops.size, 0]                            \n                    synchned_more_than_50 = np.where( durations>=min_duri )[0]\n                    if synchned_more_than_50.size == 0:\n                        onset_cont[ki,gi,repi] = np.nan\n                    else:                                \n                        onset_cont[ki,gi,repi] = timegrid[ starts[synchned_more_than_50[0]] ]\n                        counter_cont[ki,gi] += 1                                \n                        cont_synched_durations_m[ki,gi,repi]  = np.mean(  durations  )\n                        cont_synched_durations_std[ki,gi,repi]  = np.std(  durations )\n                        cont_synched_durations_max[ki,gi,repi]  = np.max(  durations )\n                        cont_synched_durations_sum[ki,gi,repi]  = np.sum(  synch_or_not[ starts[synchned_more_than_50[0]][0] :   ] )/(timegrid.size-1 -starts[synchned_more_than_50[0]] ) ##timesteps in synchronied\n\n               \n                \n            positions = np.where( Rttnon[ ki,  gi, 1:-2,repi] >=threshold)[0]\n            synch_or_not = Rttnon[ ki,  gi,1:-2, repi] >=threshold                     \n            if positions.size ==0: \n                onset_uncont[ki,gi,repi] = np.nan\n            else:\n                diffpos = np.diff(synch_or_not.astype(int)   ) #difference between consecutive steps\n                \n                diffpos2 = np.zeros(diffpos.size+1) # extend by one step at the beginning to detect \n                diffpos2[1:] = diffpos\n                starts = np.argwhere(diffpos2 == 1)\n                stops = np.argwhere(diffpos2 == -1)\n                #start_stop =[starts, stops - starts]\n                if starts.size == 0:\n                    onset_uncont[ki,gi,repi] = np.nan\n                elif stops.size == 0:\n                    onset_uncont[ki,gi,repi] = timegrid[ starts[0] ]\n                    counter_uncont[ki,gi] += 1\n                    dur = timegrid[ -1 ] - timegrid[ starts[0] ]\n                    uncont_synched_durations_m[ki,gi,repi]  = dur\n                    uncont_synched_durations_std[ki,gi,repi]  = 0\n                    uncont_synched_durations_max[ki,gi,repi]  = dur\n                    uncont_synched_durations_sum[ki,gi,repi]  = np.sum(  synch_or_not )/ (timegrid.size-1 - starts[0] )\n                else:\n                    if stops[0,0] < starts[0,0]:\n                        stops = stops[1:,:]\n                    durations = stops[:,0] - starts[:stops.size, 0]                            \n                    synchned_more_than_50 = np.where( durations>=min_duri )[0]\n                    \n                    if synchned_more_than_50.size == 0:\n                        onset_uncont[ki,gi,repi] = np.nan\n                    else:                                \n                        onset_uncont[ki,gi,repi] = timegrid[starts[synchned_more_than_50[0]] ]\n                        counter_uncont[ki,gi] += 1\n                        uncont_synched_durations_m[ki,gi,repi]  = np.mean(  durations )\n                        uncont_synched_durations_std[ki,gi,repi]  = np.std(  durations )\n                        uncont_synched_durations_max[ki,gi,repi]  = np.max(  durations  )\n                        uncont_synched_durations_sum[ki,gi,repi]  = np.sum(  synch_or_not[ starts[synchned_more_than_50[0]][0] :   ] ) /(timegrid.size-1 -starts[synchned_more_than_50[0]] )"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import seaborn as sns\nfrom matplotlib import pyplot as plt\n\n\npali = sns.diverging_palette(145, 300, s=80, as_cmap=True)\ncols = [pali(0.99), pali(0.95), pali(0.85), pali(0.80), pali(0.75), pali(0.70) ]\ncols2 = [pali(0.),pali(0.05),pali(0.1),pali(0.15),pali(0.2),pali(0.25)]\n\nfrecmap = plt.get_cmap( 'plasma')\n# minw = np.abs(np.min(ws))+0.3\n# maxw = np.max(ws)\n# intervalw = maxw +minw +0.4\n# print([ ( (wsi + np.pi/2)/(np.pi)   ) for wsi in ws       ])\n# wscols = [ frecmap( (wsi + minw)/(intervalw)   ) for wsi in ws       ]\n\nfig9 = plt.figure(constrained_layout=False, figsize=(6,3))\ngs1 = fig9.add_gridspec(nrows=4, ncols=8, wspace=1.5, hspace=1.2)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "ax01 = fig9.add_subplot(gs1[0:2,0:2 ]) \n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig9.text(0.7, -0.05, 'time', ha='center',fontsize=16)\nfig9.text(0.45, 0.5, r'order  param. $R$', va='center', rotation='vertical',fontsize=14)\n\nfig9.text(0.7, 0.95, r'coupling  $J= 3.0$', ha='center',fontsize=10)\nfig9.text(0.95, 0.95, r'noise', ha='center',fontsize=10)\nfig9.text(0.95, 0.73, r'$\\sigma= 1.0$', ha='center',fontsize=10)\nfig9.text(0.95, 0.23, r'$\\sigma= 1.5$', ha='center',fontsize=10)\n\n\n\nax02 = fig9.add_subplot(gs1[0:4,0:4 ]) \n\n#color = next(ax02._get_lines.prop_cycler)['color']\nplt.plot(Ks[1:], np.nanmean(np.nanmean(Rttcont[1 :, 0,:-2 ], axis=-2), axis=(-1)), c=cols[5],marker='^',label=r'$\\sigma=1.0$',zorder=5 ,lw=2,markersize=6,markeredgecolor='#4f4949')  \n\n#color = next(ax02._get_lines.prop_cycler)['color']\nplt.plot(Ks[1:], np.nanmean(np.nanmean(Rttcont[1 :, 1 ,:-2], axis=-2), axis=(-1)) , c=cols[1],marker='.',label=r'$\\sigma=1.5$' ,zorder=5,lw=2,markersize=10,markeredgecolor='#4f4949')    \nplt.plot(Ks[1:], np.nanmean(np.nanmean(Rttnon[1:, 0,:-2 ], axis=-2), axis=(-1)) ,'--', c=cols2[5],marker='^',label=r'$\\sigma=1.0$',lw=2 ,markersize=6,markeredgecolor='#4f4949')                                   \nplt.plot(Ks[1:], np.nanmean(np.nanmean(Rttnon[1 :, 1,:-2 ], axis=-2), axis=(-1)) ,'--', c=cols2[1],marker='.' ,label=r'$\\sigma=1.5$',lw=2,markersize=10,markeredgecolor='#4f4949')\nplt.ylim(0.5,1)    \n#plt.xlim(0,4.01)\n\nlegend = ax02.legend()        \nlegend.get_frame().set_linewidth(1.8)\nlegend.get_frame().set_facecolor('white')\nlegend.get_frame().set_edgecolor('white')\nhandles, labels = ax02.get_legend_handles_labels()\n#handles = [ handles[-2],handles[0] , handles[-1]]\n#labels= [ labels[-2],labels[0] , labels[-1]]\nif True:\n    ax02.legend(handles, labels, title=r'controlled $\\,$    uncontrolled',\n              handletextpad=0.5, columnspacing=0.5,handlelength=1.3, bbox_to_anchor=[0.10, 0.0],\n              loc=3, ncol=2, frameon=True,fontsize = 10,shadow=None,framealpha =0,edgecolor ='#0a0a0a')    \n\nplt.setp(plt.gca().get_legend().get_title(), fontsize='10') \nplt.ylabel(r'order  param. $R$') \nplt.xlabel(r'coupling $J$')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ax1 = fig9.add_subplot(gs1[0:2,4:6 ] ) \n\nplt.plot(timegrid[:-2],Rttcont[3,0,:-2,:],c=cols[5],alpha=1),\nplt.plot(timegrid[:-2],np.mean(Rttcont[3,0,:-2,:],axis=1), '--',c='#4f4949')\n\nax1.set_ylim(0.0,1.01)\n#plt.yticks([0.5,1.0])\nplt.xticks([0.5,1.5])\nax1b = fig9.add_subplot(gs1[0:2,6:8 ] ,  sharey = ax1) \n\nplt.plot(timegrid[:-2],Rttnon[3,0,:-2,:],c=cols2[5],alpha=0.5)\nax1b.set_ylim(0.0,1.01)\nplt.plot(timegrid[:-2],np.mean(Rttnon[3,0,:-2,:],axis=1), '--',c='#4f4949')\nplt.plot([0,timegrid[-2]],[1,1], '--',c='silver',zorder=0)\n#plt.ylabel(r'order  param. $R$')   \nplt.xticks([0.5,1.5])\n#plt.yticks([0.5,1.0])\nax2 = fig9.add_subplot(gs1[ 2:4,4:6], sharex = ax1)\nax2.set_ylim(0.0,1.01)\nplt.plot(timegrid[:-2],Rttcont[3,1,:-2,:],c=cols[1],alpha=1),\nplt.plot(timegrid[:-2],np.mean(Rttcont[3,1,:-2,:],axis=1), '--',c='#4f4949')\n#plt.xlabel('time')\n#plt.ylabel(r'order  param. $R$')   \n#plt.xticks([0.5,1.5])\n#plt.yticks([0.5,1.0])\nax2b = fig9.add_subplot(gs1[ 2:4,6:8], sharex = ax1b, sharey = ax2) \nplt.plot(timegrid[:-2],Rttnon[3,1,:-2,:],c=cols2[1],alpha=0.5)\nplt.plot(timegrid[:-2],np.mean(Rttnon[3,1,:-2,:],axis=1), '--',c='#4f4949')\nplt.plot([0,timegrid[-2]],[1,1], '--',c='silver',zorder=0)\nax2b.set_ylim(0.0,1.01)\n#plt.xticks([0.5,1.5])\n\n\n\nplt.savefig('systematic_Kuramoto_6N.png', bbox_inches='tight',dpi=300 , pad_inches = 0.)\nplt.savefig('systematic_Kuramoto_6N.pdf', bbox_inches='tight',dpi=300,  pad_inches = 0.)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ax6 = fig9.add_subplot(gs1[0:4,4:8 ]) \n\nplt.plot(Ks, np.nanmean(onset_cont[:, 0], axis=(-1)) , c=cols[5],marker='^',label=r'$\\sigma=0.5$',zorder=5 ,lw=2.8,markersize=6,markeredgecolor='#4f4949') \nplt.plot(Ks, np.nanmean(onset_uncont[:, 0], axis=(-1)),'--', c=cols2[5],marker='^',label=r'$\\sigma=0.5$' ,zorder=5,lw=2.8,markersize=6,markeredgecolor='#4f4949')  \n\nplt.plot(Ks, np.nanmean(onset_cont[:, 1], axis=(-1)) , c=cols[1],marker='.',label=r'$\\sigma=1.0$',lw=2.8 ,markersize=12,markeredgecolor='#4f4949') \nplt.plot(Ks, np.nanmean(onset_uncont[:, 1], axis=(-1)),'--', c=cols2[0],marker='.' ,label=r'$\\sigma=1.0$',lw=2.8,markersize=12,markeredgecolor='#4f4949')   \n\n                   \nplt.ylabel(r'onset of synchrony $t^{syn}$') \nplt.xlabel(r'coupling $J$')                \n\n\n\nax6.tick_params(axis='both',which='both',direction='in', length=3, width=1,colors='#4f4949',zorder=3)\nax6.tick_params(bottom=True, top=True, left=True, right=True)\nax6.spines['top'].set_visible(True)\nax6.spines['right'].set_visible(True)\nax6.minorticks_on()\nax6.tick_params(axis='both',which='major',direction='in', length=3.5, width=1,colors='#4f4949',zorder=3)\nax6.tick_params(axis='both',which='minor',direction='in', length=2.5, width=0.5,colors='#4f4949',zorder=3)\nax6.tick_params(bottom=True, top=True, left=True, right=True)\nax6.tick_params(axis='both', which='minor', bottom=True, top=True, left=True, right=True)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ax7 = fig9.add_subplot(gs1[0:4,8:12 ]) \n\n\n\nplt.plot(Ks, np.nanmean(cont_synched_durations_sum[:, 0], axis=(-1)) , c=cols[5],marker='^',label=r'$\\sigma=1.0$',zorder=2 ,lw=2.8,markersize=6,markeredgecolor='#4f4949') \nplt.plot(Ks, np.nanmean(uncont_synched_durations_sum[:, 0], axis=(-1)),'--', c=cols2[5],marker='^',label=r'$\\sigma=1.0$' ,zorder=2,lw=2.8,markersize=6,markeredgecolor='#4f4949')  \n\n\nplt.plot(Ks, np.nanmean(cont_synched_durations_sum[:, 1], axis=(-1)) , c=cols[1],marker='.',label=r'$\\sigma=1.5$',lw=2.8 ,markersize=12,markeredgecolor='#4f4949') \nplt.plot(Ks, np.nanmean(uncont_synched_durations_sum[:, 1], axis=(-1)),'--', c=cols2[0],marker='.' ,label=r'$\\sigma=1.5$',lw=2.8,markersize=12,markeredgecolor='#4f4949')   \n\nplt.ylabel('$\\%$ time synchronised\\nafter $t^{syn}$', multialignment='center',labelpad=-0.1) # \"Mat\\nTTp\\n123\"\nplt.xlabel(r'coupling $J$') \nplt.yticks([0, 0.5,1])\nax7 = plt.gca()\nax7.tick_params(axis='both',which='both',direction='in', length=3, width=1,colors='#4f4949',zorder=3)\nax7.tick_params(bottom=True, top=True, left=True, right=True)\nax7.spines['top'].set_visible(True)\nax7.spines['right'].set_visible(True)\nax7.minorticks_on()\nax7.tick_params(axis='both',which='major',direction='in', length=3.5, width=1,colors='#4f4949',zorder=3)\nax7.tick_params(axis='both',which='minor',direction='in', length=2.5, width=0.5,colors='#4f4949',zorder=3)\nax7.tick_params(bottom=True, top=True, left=True, right=True)\nax7.tick_params(axis='both', which='minor', bottom=True, top=True, left=True, right=True)\n\nlegend = ax7.legend()        \nlegend.get_frame().set_linewidth(1.8)\nlegend.get_frame().set_facecolor('white')\nlegend.get_frame().set_edgecolor('white')\nhandles, labels = ax7.get_legend_handles_labels()\nhandles = [ handles[0],handles[2] , handles[1], handles[3]]\nlabels= [ labels[0],labels[2] , labels[1], labels[3]]\nif True:\n    ax7.legend(handles, labels, title=r'$\\,$controlled $\\,\\,$ uncontrolled',\n              handletextpad=0.5, columnspacing=1,handlelength=0.45, bbox_to_anchor=[-0.65, 0.99],\n              loc=3, ncol=2, frameon=True,fontsize = 'small',shadow=None,framealpha =0,edgecolor ='#0a0a0a')    \nplt.setp(plt.gca().get_legend().get_title(), fontsize='10') \nplt.subplots_adjust(wspace=0.325)#, hspace=0)"
      ]
    }
  ],
  "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.9.15"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}