Note
Click here to download the full example code
Created on Tue Nov 9 19:03:18 2021
@author: maout
import numpy as np
#from matplotlib import pyplot as plt
from scipy.spatial.distance import cdist
import time
import ot
import numba
import math
import random
import sys
import pickle
dim = 2
h = 0.001
t1 = 0
t2 = 1.5
T = t2-t1
timegrid = np.arange(0,T+h/2,h)
N = 2000#0#0#0#2000
#g = 1
k = timegrid.size
M = 80
y2 = np.ones(dim)
sigmas = np.array([0.1,0.5,1])
rep_bridge = 5 #10 different bridge instances for every setting
reps = 20 ##instanses for stochastic path evaluation of each bridge
# dy0 = np.array([np.pi/4, np.pi/2, np.pi, 3*np.pi/2])
# y0s = np.zeros((dim,rep_bridge, dy0.size))
random.seed(22)
# for i in range(rep_bridge):
# y0s[0, i, :] = np.random.uniform( low=0, high=2*np.pi,size=1 )
# y0s[1, i, :] = (y0s[0, i, :] + dy0) %(2* np.pi)
Ks = np.linspace(0,4,11)
dy0 = np.array([np.pi/4, np.pi/2])
y0s = np.zeros((dim,rep_bridge, dy0.size))
Rttcont = np.zeros(( rep_bridge, Ks.size, dy0.size, sigmas.size, timegrid.size,reps ))*np.nan
Rttnon = np.zeros(( rep_bridge, Ks.size, dy0.size, sigmas.size, timegrid.size,reps ))*np.nan
used_us = np.zeros((dim, rep_bridge, Ks.size, dy0.size, sigmas.size,timegrid.size, reps ))*np.nan
bis = [2,3,4]
for bi in bis:#range(rep_bridge):
###loop over bridge iterations
#for ki,K in enumerate(Ks): ##loop over interaction strengths
for ki,K in enumerate(Ks):
for yi in range(dy0.size): ###loop over initial distances
for gii,g in enumerate(sigmas):
#naming = 'kuramot\\2N_systematic_Kuramoto_homogeneous_repb_%d_ki_%d_yi_%d__gi_%d_N_%d_M_%d'%(bi, ki,yi,gii,N, M)
naming = 'kurgamot\\New_2N_systematic_Kuramoto_homogeneous_repb_%d_ki_%d_yi_%d__gi_%d_N_%d_M_%d'%(bi, ki,yi,gii,N, M)
try:
file = open(naming+'.dat','rb')
to_save = pickle.load(file)
# 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)
# file2 = open(naming2+'.dat','rb')
# to_save2 = pickle.load(file2)
#Fcont = to_save['Fcont']
Rttcont[bi,ki,yi, gii,:,:] = to_save['Rttcont']
#Fnon = to_save2['Fnon']
Rttnon[bi,ki,yi, gii,:,:] = to_save['Rttnon']
Kin = to_save['K']
y0s[:, bi, yi] = to_save['y0']
used_us[:,bi,ki,yi, gii,:] = to_save['used_us']
except FileNotFoundError:
i=0
print(naming)
from scipy.spatial.distance import cdist
import time
import ot
import numba
import math
import random
import sys
import pickle
dim = 2
h = 0.001
t1 = 0
t2 = 1.5
T = t2-t1
timegrid = np.arange(0,T+h/2,h)
N = 2000#0#0#0#2000
#g = 1
k = timegrid.size
M = 80
bi = 2
ki = 3
gii = 2
yi = 1
g=1
naming = 'kurgamot\\New_2N_systematic_Kuramoto_homogeneous_repb_%d_ki_%d_yi_%d__gi_%d_N_%d_M_%d'%(bi, ki,yi,gii,N, M)
file = open(naming+'.dat','rb')
to_save = pickle.load(file)
Rttconti = to_save['Rttcont']
Fcont = to_save['Fcont']
Fnon = to_save['Fnon']
Rttnoni = to_save['Rttnon']
used_u = to_save['used_us']
import seaborn as sns
from matplotlib import pyplot as plt
repi = 0
pali = sns.diverging_palette(145, 300, s=60, as_cmap=True)
from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes, mark_inset
import matplotlib.gridspec as gridspec
plt.figure(figsize= (8,2.)),
plt.subplot(1,3,1),
#cols = [pali(0), pali(0.99) ]
#cols2 = [pali(0.3),pali(0.7)]
cols = [pali(0.99), pali(0.8) ]
cols2 = [pali(0.),pali(0.2)]
for di in range(dim):
plt.plot(timegrid[:-1],Fcont[di,:-1,repi],'-',lw=4,c=cols[di],label=r':math:`\theta_{%d}`'%(di+1))
plt.plot(timegrid[:-1],Fnon[di,:-1,repi],'.',lw=3,c=cols2[di],alpha=0.7,zorder=0)
for di in range(dim):
plt.plot(timegrid[:-1],Fnon[di,:-1, repi],'-',lw=3,c=cols2[di],zorder=0,label=r':math:`\theta_{%d}`'%(di+1)) #this line is only for the legend entry
plt.xlabel('time')
plt.ylabel(r'phase :math:`\theta`')
ax = plt.gca()
ax.tick_params(axis='both',which='both',direction='out', length=3, width=1,colors='#4f4949',zorder=3)
ax.tick_params(bottom=True, top=False, left=True, right=False)
ax.set_yticks([0,np.pi/2,np.pi, 3*np.pi/2, 2*np.pi])
ax.set_yticklabels([r':math:`0`',r':math:`\pi/2`',r':math:`\pi`', r':math:`{3 \pi}/{2}`', r':math:`2 \pi`'])
# legend = plt.legend(frameon = 1,prop={'size': 15})
# frame = legend.get_frame()
# frame.set_facecolor('white')
# frame.set_edgecolor('white')
plt.xticks([0,0.5,1,1.5])
legend = ax.legend()
legend.get_frame().set_linewidth(1.8)
legend.get_frame().set_facecolor('white')
legend.get_frame().set_edgecolor('white')
handles, labels = ax.get_legend_handles_labels()
if True:#g==0.5:
ax.legend(handles, labels, title='controlled uncontrolled',
handletextpad=0.5, columnspacing=3.2,handlelength=0.8, bbox_to_anchor=[0.0, 0.785],
loc=3, ncol=2, frameon=True,fontsize = 'small',shadow=None,framealpha =0,edgecolor ='#0a0a0a')
else:
ax.legend(handles, labels, title='controlled uncontrolled',
handletextpad=0.5, columnspacing=3.2,handlelength=0.8,
loc=3, ncol=2, frameon=True,fontsize = 8,shadow=None,framealpha =0,edgecolor ='#0a0a0a')
plt.setp(plt.gca().get_legend().get_title(), fontsize='10')
plt.subplot(1,3,2),
plt.plot(timegrid[:-1],Rttnoni[:-1,repi],lw=3,c= cols2[0],zorder=3,label='uncontrolled'),
plt.plot(timegrid[:-1],Rttconti[:-1,repi],lw=3, c=cols[0],zorder=4,label='controlled'),
plt.plot([0,timegrid[-1]], [1/np.sqrt(dim),1/np.sqrt(dim) ],linestyle=(0,(4,3)),lw=3, c='#ee9222',alpha=0.5,label='independent'),
plt.plot([0,timegrid[-1]], [1,1 ],linestyle=(0,(4,3)),c='grey',lw=3,dash_capstyle = "round"),
plt.xlabel('time')
plt.xticks([0,0.5,1,1.5])
plt.ylabel(r'order param. :math:`R`')
ax7 = plt.gca()
ax7.tick_params(axis='both',which='both',direction='out', length=3, width=1,colors='#4f4949',zorder=3)
ax7.tick_params(bottom=True, top=False, left=True, right=False)
legend = ax7.legend()
legend.get_frame().set_linewidth(1.8)
legend.get_frame().set_facecolor('white')
legend.get_frame().set_edgecolor('white')
handles, labels = ax7.get_legend_handles_labels()
handles = [ handles[-2],handles[0] , handles[-1]]
labels= [ labels[-2],labels[0] , labels[-1]]
if True:
ax7.legend(handles, labels,
handletextpad=0.5, columnspacing=0.3,handlelength=0.5, bbox_to_anchor=[-0.4, 0.99],
loc=3, ncol=3, frameon=True,fontsize = 'small',shadow=None,framealpha =0,edgecolor ='#0a0a0a')
###########
# axins = zoomed_inset_axes(ax,2, loc=1, bbox_to_anchor=(2400,800))
# mark_inset(ax, axins, loc1=2, loc2=1, fc="none", ec="0.5")
# axins.set_xlim([0.5,1.2])
# axins.set_ylim([0.95,1.01])
# # # Plot zoom window
# axins.plot(timegrid[:tti],Rtt[:tti],lw=3,c= cols2[0],zorder=3),
# axins.plot(timegrid[:tti],Rttcont[:tti],lw=3, c=cols[0],zorder=4),
# axins.plot([0,timegrid[tti]], [1,1 ],linestyle=(0,(4,3)),c='grey',lw=3,dash_capstyle = "round"),
# axins.tick_params(bottom=False, top=False, left=True, right=False)
plt.subplot(1,3,3),
#for ti in range(2):
for di in range(dim):
plt.plot(timegrid[:],h*g**2*used_u[di,:,repi],lw=2.5, c=cols[di],label=r':math:`u_{%d}`'%(di+1))
plt.xticks([0,0.5,1,1.5])
plt.xlabel('time')
plt.ylabel(r'control :math:`u`')
ax = plt.gca()
ax.tick_params(axis='both',which='both',direction='out', length=3, width=1,colors='#4f4949',zorder=3)
ax.tick_params(bottom=True, top=False, left=True, right=False)
legend = ax.legend()
legend.get_frame().set_linewidth(1.8)
legend.get_frame().set_facecolor('white')
legend.get_frame().set_edgecolor('white')
handles, labels = ax.get_legend_handles_labels()
ax.legend(handles, labels, title=None,
handletextpad=0.5, columnspacing=0,handlelength=1, bbox_to_anchor=[0.50, 0.35],
loc=1, ncol=1, frameon=True,fontsize = 'small',shadow=None,framealpha =0,edgecolor ='#0a0a0a')
naming2 = 'Kuramoto_2N_trajectories'
##%%
#plt.savefig( 'naming.png', bbox_inches='tight')
#plt.savefig( 'naming.pdf', bbox_inches='tight')
plt.savefig( naming2+'.png', bbox_inches='tight')
plt.savefig(naming2 +'.pdf', bbox_inches='tight')
from matplotlib import pyplot as plt
#figure and plotting settings
import seaborn as sns
plt.rcParams['figure.dpi'] = 300
plt.rcParams['savefig.dpi'] = 300
plt.rcParams['savefig.facecolor'] = (1,1,1,0)
plt.rcParams['savefig.bbox'] = "tight"
plt.rcParams['font.size'] = 10*1.2
# plt.rcParams['font.family'] = 'sans-serif' # not available in Colab
# plt.rcParams['font.sans-serif'] = 'Helvetica' # not available in Colab
plt.rcParams['pdf.fonttype'] = 42
plt.rcParams['xtick.labelsize'] = 10*1.2
plt.rcParams['ytick.labelsize'] = 10*1.2
plt.rcParams['axes.labelsize'] = 12*1.2
plt.rcParams['axes.linewidth'] = 2
plt.rcParams['axes.spines.top'] = False
plt.rcParams['axes.spines.right'] = False
#plt.rcParams['axes.Axes.tick_params'] = True
#plt.rcParams['axes.solid_capstyle'] = 'round'
plt.rc('axes',edgecolor='#4f4949')
plt.rcParams['figure.frameon'] = False
plt.rcParams['figure.subplot.hspace'] = 0.51
plt.rcParams['figure.subplot.wspace'] = 0.51
plt.rcParams['figure.subplot.left'] = 0.1
plt.rcParams['figure.subplot.right'] = 0.9
plt.rcParams['figure.subplot.top'] = 0.9
plt.rcParams['figure.subplot.bottom'] = 0.1
plt.rcParams['lines.solid_capstyle'] = 'round'
plt.rcParams['lines.solid_joinstyle'] = 'round'
#plt.rcParams['xtick.major.size'] = 20
#plt.rcParams['xtick.major.width'] = 4
plt.rcParams['text.usetex'] = True
bi=2
###bi=22 , k=4 and k=5 is good for plotting
# R [ bi, ki, yi, gi, :,:]
plt.figure(),
for ki in range(Ks.size):
plt.subplot(2,6,ki+1)
plt.title(ki)
yi=0
#for yi in range(dy0.size):
#plt.subplot(1,dy0.size,yi+1)
plt.plot(timegrid[:-2], np.nanmean(Rttcont[2:4,ki, :, :2,:-2 ], axis=(0,2,-1)).T )
plt.plot(timegrid[:-2], np.nanmean(Rttnon[2:4,ki, :, :2,:-2 ], axis=(0,2,-1)).T ,'--',alpha=0.4 )
plt.ylim(0.5,1.1)
# R [ bi, ki, yi, gi, :,:]
bi=3
plt.figure(),
for yi in range(dy0.size):
plt.subplot(1,dy0.size,yi+1)
ax = plt.gca()
# color = next(ax._get_lines.prop_cycler)['color']
# plt.plot(Ks, np.nanmean(np.nanmean(Rttcont[bi,:, yi, 0 ], axis=-2), axis=(-1)) , c=color)
# plt.plot(Ks, np.nanmean(np.nanmean(Rttnon[bi,:, yi, 0 ], axis=-2), axis=(-1)) ,'--', c=color)
color = next(ax._get_lines.prop_cycler)['color']
plt.plot(Ks[0::2], np.nanmean(np.nanmean(Rttcont[2:5,0::2, yi, 1,:-2 ], axis=-2), axis=(0,-1)), c=color,marker='.' )
plt.plot(Ks[0::2], np.nanmean(np.nanmean(Rttnon[2:5,0::2, yi, 1,:-2 ], axis=-2), axis=(0,-1)) ,'--', c=color,marker='.' )
color = next(ax._get_lines.prop_cycler)['color']
plt.plot(Ks[0::2], np.nanmean(np.nanmean(Rttcont[2:5,0::2, yi, 2 ,:-2], axis=-2), axis=(0,-1)) , c=color,marker='.' )
plt.plot(Ks[0::2], np.nanmean(np.nanmean(Rttnon[2:5,0::2, yi, 2,:-2 ], axis=-2), axis=(0,-1)) ,'--', c=color,marker='.' )
plt.ylim(0.5,1)
plt.figure()
plt.plot(Rttcont[2,3,1,2,:-2]),
plt.plot(Rttnon[2,3,1,2,:-2],'grey')
plt.figure()
plt.plot(Rttcont[2,5,1,1,:-2]),
plt.plot(Rttnon[2,5,1,1,:-2],'grey',alpha=0.5)
plt.figure()
plt.plot(Rttcont[2,5,1,2,:-2]),
plt.plot(Rttnon[2,5,1,2,:-2],'grey',alpha=0.5)
bis= 2
bie = 5
plt.figure(),
for yi in range(dy0.size):
plt.subplot(1,dy0.size,yi+1)
ax = plt.gca()
# color = next(ax._get_lines.prop_cycler)['color']
# plt.plot(Ks, np.nanmean(np.nanmean(Rttcont[bi,:, yi, 0 ], axis=-2), axis=(-1)) , c=color)
# plt.plot(Ks, np.nanmean(np.nanmean(Rttnon[bi,:, yi, 0 ], axis=-2), axis=(-1)) ,'--', c=color)
color = next(ax._get_lines.prop_cycler)['color']
plt.plot(Ks[1::2], np.nanmean(np.nanmean(Rttcont[bis:bie,1::2, yi, 1,:-2 ], axis=-2), axis=(0,-1)), c=color,marker='.' )
plt.plot(Ks[1::2], np.nanmean(np.nanmean(Rttnon[bis:bie,1::2, yi, 1,:-2 ], axis=-2), axis=(0,-1)) ,'--', c=color,marker='.' )
color = next(ax._get_lines.prop_cycler)['color']
plt.plot(Ks[1::2], np.nanmean(np.nanmean(Rttcont[bis:bie,1::2, yi, 2 ,:-2], axis=-2), axis=(0,-1)) , c=color,marker='.' )
plt.plot(Ks[1::2], np.nanmean(np.nanmean(Rttnon[bis:bie,1::2, yi, 2,:-2 ], axis=-2), axis=(0,-1)) ,'--', c=color,marker='.' )
plt.ylim(0.5,1)
bi=3
plt.figure(),
for yi in range(0,4):
plt.subplot(1,4,yi+1)
ax = plt.gca()
# color = next(ax._get_lines.prop_cycler)['color']
# plt.plot(Ks[1::2], np.nanmean(np.nanmean(Rttcont[yi+0:5,1::2, :, 0,:-2 ], axis=-2), axis=(0,2,-1)), c=color,marker='.' )
# plt.plot(Ks[1::2], np.nanmean(np.nanmean(Rttnon[yi+0:5,1::2, :, 0,:-2 ], axis=-2), axis=(0,2,-1)) ,'--', c=color,marker='.' )
# plt.plot(Ks, np.nanmean(np.nanmean(Rttnon[bi,:, yi, 0 ], axis=-2), axis=(-1)) ,'--', c=color)
color = next(ax._get_lines.prop_cycler)['color']
plt.plot(Ks[1:-1:2], np.nanmean(np.nanmean(Rttcont[yi+0:5,1::2, :, 1,:-2 ], axis=-2), axis=(0,2,-1)), c=color,marker='.' )
plt.plot(Ks[1:-1:2], np.nanmean(np.nanmean(Rttnon[yi+0:5,1::2, :, 1,:-2 ], axis=-2), axis=(0,2,-1)) ,'--', c=color,marker='.' )
color = next(ax._get_lines.prop_cycler)['color']
plt.plot(Ks[1:-1:2], np.nanmean(np.nanmean(Rttcont[yi+0:5,1::2, :, 2 ,:-2], axis=-2), axis=(0,2,-1)) , c=color,marker='.' )
plt.plot(Ks[1:-1:2], np.nanmean(np.nanmean(Rttnon[yi+0:5, 1::2, :, 2,:-2 ], axis=-2), axis=(0,2,-1)) ,'--', c=color,marker='.' )
plt.ylim(0.5,1)
#figure and plotting settings
import seaborn as sns
plt.rcParams['figure.dpi'] = 300
plt.rcParams['savefig.dpi'] = 300
plt.rcParams['savefig.facecolor'] = (1,1,1,0)
plt.rcParams['savefig.bbox'] = "tight"
plt.rcParams['font.size'] = 10*1.2
# plt.rcParams['font.family'] = 'sans-serif' # not available in Colab
# plt.rcParams['font.sans-serif'] = 'Helvetica' # not available in Colab
plt.rcParams['pdf.fonttype'] = 42
plt.rcParams['xtick.labelsize'] = 10*1.2
plt.rcParams['ytick.labelsize'] = 10*1.2
plt.rcParams['axes.labelsize'] = 12*1.2
plt.rcParams['axes.linewidth'] = 2
plt.rcParams['axes.spines.top'] = False
plt.rcParams['axes.spines.right'] = False
#plt.rcParams['axes.Axes.tick_params'] = True
#plt.rcParams['axes.solid_capstyle'] = 'round'
plt.rc('axes',edgecolor='#4f4949')
plt.rcParams['figure.frameon'] = False
plt.rcParams['figure.subplot.hspace'] = 0.51
plt.rcParams['figure.subplot.wspace'] = 0.51
plt.rcParams['figure.subplot.left'] = 0.1
plt.rcParams['figure.subplot.right'] = 0.9
plt.rcParams['figure.subplot.top'] = 0.9
plt.rcParams['figure.subplot.bottom'] = 0.1
plt.rcParams['lines.solid_capstyle'] = 'round'
plt.rcParams['lines.solid_joinstyle'] = 'round'
plt.rcParams['xtick.minor.size'] = 10
#plt.rcParams['xtick.major.width'] = 4
plt.rcParams['text.usetex'] = True
pali = sns.diverging_palette(145, 300, s=80, as_cmap=True)
cols = [pali(0.99), pali(0.95), pali(0.85), pali(0.80), pali(0.75), pali(0.70) ]
cols2 = [pali(0.),pali(0.05),pali(0.1),pali(0.15),pali(0.2),pali(0.25)]
frecmap = plt.get_cmap( 'plasma')
# minw = np.abs(np.min(ws))+0.3
# maxw = np.max(ws)
# intervalw = maxw +minw +0.4
# print([ ( (wsi + np.pi/2)/(np.pi) ) for wsi in ws ])
# wscols = [ frecmap( (wsi + minw)/(intervalw) ) for wsi in ws ]
fig9 = plt.figure(constrained_layout=False, figsize=(6,3))
gs1 = fig9.add_gridspec(nrows=4, ncols=8, wspace=1.5, hspace=1.2)
ax01 = fig9.add_subplot(gs1[0:2,0:2 ])
fig9.text(0.7, -0.05, 'time', ha='center',fontsize=16)
fig9.text(0.45, 0.5, r'order param. :math:`R`', va='center', rotation='vertical',fontsize=14)
fig9.text(0.7, 0.95, r'coupling :math:`J= 2.0`', ha='center',fontsize=10)
fig9.text(0.95, 0.95, r'noise', ha='center',fontsize=10)
fig9.text(0.95, 0.73, r':math:`\sigma= 0.5`', ha='center',fontsize=10)
fig9.text(0.95, 0.23, r':math:`\sigma= 1.0`', ha='center',fontsize=10)
# plt.ylabel(r'order param. :math:`R`')
# plt.xlabel('time')
# plt.ylabel(r'order param. :math:`R`')
ax02 = fig9.add_subplot(gs1[0:4,0:4 ])
#color = next(ax02._get_lines.prop_cycler)['color']
plt.plot(Ks[1:-1:2], np.nanmean(np.nanmean(Rttcont[2:5,1::2, :, 1,:-2 ], axis=-2), axis=(0,2,-1)), c=cols[5],marker='^',label=r':math:`\sigma=0.5`',zorder=5 ,lw=2,markersize=6,markeredgecolor='#4f4949')
#color = next(ax02._get_lines.prop_cycler)['color']
plt.plot(Ks[1:-1:2], np.nanmean(np.nanmean(Rttcont[2:5,1::2, :, 2 ,:-2], axis=-2), axis=(0,2,-1)) , c=cols[1],marker='.',label=r':math:`\sigma=1.0`' ,zorder=5,lw=2,markersize=10,markeredgecolor='#4f4949')
plt.plot(Ks[1:-1:2], np.nanmean(np.nanmean(Rttnon[2:5,1::2, :, 1,:-2 ], axis=-2), axis=(0,2,-1)) ,'--', c=cols2[5],marker='^',label=r':math:`\sigma=0.5`',lw=2 ,markersize=6,markeredgecolor='#4f4949')
plt.plot(Ks[1:-1:2], np.nanmean(np.nanmean(Rttnon[2:5, 1::2, :, 2,:-2 ], axis=-2), axis=(0,2,-1)) ,'--', c=cols2[1],marker='.' ,label=r':math:`\sigma=1.0`',lw=2,markersize=10,markeredgecolor='#4f4949')
plt.ylim(0.5,1)
#plt.xlim(0,4.01)
legend = ax02.legend()
legend.get_frame().set_linewidth(1.8)
legend.get_frame().set_facecolor('white')
legend.get_frame().set_edgecolor('white')
handles, labels = ax02.get_legend_handles_labels()
#handles = [ handles[-2],handles[0] , handles[-1]]
#labels= [ labels[-2],labels[0] , labels[-1]]
if True:
ax02.legend(handles, labels, title=r'controlled :math:`\,` uncontrolled',
handletextpad=0.5, columnspacing=0.5,handlelength=1.3, bbox_to_anchor=[-0.0, 0.0],
loc=3, ncol=2, frameon=True,fontsize = 10,shadow=None,framealpha =0,edgecolor ='#0a0a0a')
plt.setp(plt.gca().get_legend().get_title(), fontsize='10')
plt.ylabel(r'order param. :math:`R`')
plt.xlabel(r'coupling :math:`J`')
ax1 = fig9.add_subplot(gs1[0:2,4:6 ] )
plt.plot(timegrid[:-2],Rttcont[2,5,1,1,:-2],c=cols[5],alpha=1),
plt.plot(timegrid[:-2],np.mean(Rttcont[2,5,1,1,:-2],axis=1), '--',c='#4f4949')
ax1.set_ylim(0.0,1.01)
#plt.yticks([0.5,1.0])
plt.xticks([0.5,1.5])
ax1b = fig9.add_subplot(gs1[0:2,6:8 ] , sharey = ax1)
plt.plot(timegrid[:-2],Rttnon[2,5,1,1,:-2],c=cols2[5],alpha=0.5)
ax1b.set_ylim(0.0,1.01)
plt.plot(timegrid[:-2],np.mean(Rttnon[2,5,1,1,:-2],axis=1), '--',c='#4f4949')
#plt.ylabel(r'order param. :math:`R`')
plt.xticks([0.5,1.5])
#plt.yticks([0.5,1.0])
ax2 = fig9.add_subplot(gs1[ 2:4,4:6], sharex = ax1)
ax2.set_ylim(0.0,1.01)
plt.plot(timegrid[:-2],Rttcont[2,5,1,2,:-2],c=cols[1],alpha=1),
plt.plot(timegrid[:-2],np.mean(Rttcont[2,5,1,2,:-2],axis=1), '--',c='#4f4949')
#plt.xlabel('time')
#plt.ylabel(r'order param. :math:`R`')
plt.xticks([0.5,1.5])
#plt.yticks([0.5,1.0])
ax2b = fig9.add_subplot(gs1[ 2:4,6:8], sharex = ax1b, sharey = ax2)
plt.plot(timegrid[:-2],Rttnon[2,5,1,2,:-2],c=cols2[1],alpha=0.5)
plt.plot(timegrid[:-2],np.mean(Rttnon[2,5,1,2,:-2],axis=1), '--',c='#4f4949')
ax2b.set_ylim(0.0,1.01)
plt.xticks([0.5,1.5])
#plt.savefig("Kuramoto_2.png", bbox_inches='tight')
#plt.savefig("Kuramoto_2.pdf", bbox_inches='tight')
# R [ bi, ki, yi, gi, :,:]
onset_cont = np.zeros((3,Ks.size,2,3,20)) *np.nan
onset_uncont = np.zeros((3,Ks.size,2,3,20)) *np.nan
cont_synched_durations_m = np.zeros((3,Ks.size,2,3,20)) *np.nan
cont_synched_durations_std = np.zeros((3,Ks.size,2,3,20)) *np.nan
cont_synched_durations_max = np.zeros((3,Ks.size,2,3,20)) *np.nan
cont_synched_durations_sum = np.zeros((3,Ks.size,2,3,20)) *np.nan
uncont_synched_durations_m = np.zeros((3,Ks.size,2,3,20)) *np.nan
uncont_synched_durations_std = np.zeros((3,Ks.size,2,3,20)) *np.nan
uncont_synched_durations_max = np.zeros((3,Ks.size,2,3,20)) *np.nan
uncont_synched_durations_sum = np.zeros((3,Ks.size,2,3,20)) *np.nan
counter_cont = np.zeros((3,Ks.size,2,3))
counter_uncont = np.zeros((3,Ks.size,2,3))
threshold = 0.99
min_duri = 20
for bi in range(3):
bii = bi+2 ###bi should start from 2 for the data matrices since I don't have simulation results for bi:0-1
for ki in range(Ks.size):
for yi in range(2):#
for gi in range(1,3):#3
for repi in range(20):
positions = np.where( Rttcont[bii, ki, yi, gi,1:-2, repi] >=threshold)[0] ###this one detects if there is at all any synchronisation
synch_or_not = Rttcont[bii, ki, yi, gi,1:-2, repi] >=threshold ##this is boolean indicating the synchronised positions
if positions.size==0:
onset_cont[bi,ki,yi,gi,repi] = np.nan
else:
diffpos = np.diff(synch_or_not.astype(int) ) #difference between consecutive steps
diffpos2 = np.zeros(diffpos.size+1) # extend by one step at the beginning to detect
diffpos2[1:] = diffpos
starts = np.argwhere(diffpos2 == 1)
stops = np.argwhere(diffpos2 == -1)
if starts.size == 0:
onset_cont[bi,ki,yi,gi,repi] = np.nan
elif stops.size == 0:
onset_cont[bi,ki,yi,gi,repi] = timegrid[ starts[0] ]
counter_cont[bi,ki,yi,gi] += 1
dur = timegrid[ -1 ] - timegrid[ starts[0] ]
cont_synched_durations_m[bi,ki,yi,gi,repi] = dur
cont_synched_durations_std[bi,ki,yi,gi,repi] = 0
cont_synched_durations_max[bi,ki,yi,gi,repi] = dur
cont_synched_durations_sum[bi,ki,yi,gi,repi] = np.sum( synch_or_not )/ (timegrid.size-1 - starts[0] )
else:
durations = stops[:,0] - starts[:stops.size, 0]
synchned_more_than_50 = np.where( durations>=min_duri )[0]
if synchned_more_than_50.size == 0:
onset_cont[bi,ki,yi,gi,repi] = np.nan
else:
onset_cont[bi,ki,yi,gi,repi] = timegrid[ starts[synchned_more_than_50[0]] ]
counter_cont[bi,ki,yi,gi] += 1
cont_synched_durations_m[bi,ki,yi,gi,repi] = np.mean( durations )
cont_synched_durations_std[bi,ki,yi,gi,repi] = np.std( durations )
cont_synched_durations_max[bi,ki,yi,gi,repi] = np.max( durations )
cont_synched_durations_sum[bi,ki,yi,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
positions = np.where( Rttnon[bii, ki, yi, gi, 1:-2,repi] >=threshold)[0]
synch_or_not = Rttnon[bii, ki, yi, gi,1:-2, repi] >=threshold
if positions.size ==0:
onset_uncont[bi,ki,yi,gi,repi] = np.nan
else:
diffpos = np.diff(synch_or_not.astype(int) ) #difference between consecutive steps
diffpos2 = np.zeros(diffpos.size+1) # extend by one step at the beginning to detect
diffpos2[1:] = diffpos
starts = np.argwhere(diffpos2 == 1)
stops = np.argwhere(diffpos2 == -1)
#start_stop =[starts, stops - starts]
if starts.size == 0:
onset_uncont[bi,ki,yi,gi,repi] = np.nan
elif stops.size == 0:
onset_uncont[bi,ki,yi,gi,repi] = timegrid[ starts[0] ]
counter_uncont[bi,ki,yi,gi] += 1
dur = timegrid[ -1 ] - timegrid[ starts[0] ]
uncont_synched_durations_m[bi,ki,yi,gi,repi] = dur
uncont_synched_durations_std[bi,ki,yi,gi,repi] = 0
uncont_synched_durations_max[bi,ki,yi,gi,repi] = dur
uncont_synched_durations_sum[bi,ki,yi,gi,repi] = np.sum( synch_or_not )/ (timegrid.size-1 - starts[0] )
else:
durations = stops[:,0] - starts[:stops.size, 0]
synchned_more_than_50 = np.where( durations>=min_duri )[0]
if synchned_more_than_50.size == 0:
onset_uncont[bi,ki,yi,gi,repi] = np.nan
else:
onset_uncont[bi,ki,yi,gi,repi] = timegrid[starts[synchned_more_than_50[0]] ]
counter_uncont[bi,ki,yi,gi] += 1
uncont_synched_durations_m[bi,ki,yi,gi,repi] = np.mean( durations )
uncont_synched_durations_std[bi,ki,yi,gi,repi] = np.std( durations )
uncont_synched_durations_max[bi,ki,yi,gi,repi] = np.max( durations )
uncont_synched_durations_sum[bi,ki,yi,gi,repi] = np.sum( synch_or_not[ starts[synchned_more_than_50[0]][0] : ] ) /(timegrid.size-1 -starts[synchned_more_than_50[0]] )
count_perc = np.sum(counter_uncont[0:3,1::2,:1, :],axis=0)/60 #counts the percentage of synchronised uncontrolled network for the annotations
plt.figure(figsize=(6.5,2))
yi = 0
gi = 1
plt.subplot(1,2,1)
ax1 = plt.gca()
#er1 = plt.errorbar(Ks[1:-1:2], np.nanmean(onset_cont[0:3,1::2,yi, 1], axis=(0,-1)), yerr=np.nanstd(onset_cont[0:3,1::2,yi, 1], axis=(0,-1)), marker="o", linestyle="-",uplims=True,lolims=True,)
#er2 = plt.errorbar(Ks[1:-1:2], np.nanmean(onset_uncont[0:3,1::2,yi, 1], axis=(0,-1)), yerr=np.nanstd(onset_uncont[0:3,1::2,yi, 1], axis=(0,-1)), marker="o", linestyle="-",uplims=True,lolims=True)
plt.plot(Ks[1:-1:2], np.nanmean(onset_cont[0:3,1::2,:1, 1], axis=(0,2,-1)) , c=cols[5],marker='^',label=r':math:`\sigma=0.5`',zorder=5 ,lw=2.8,markersize=6,markeredgecolor='#4f4949')
plt.plot(Ks[1:-1:2], np.nanmean(onset_uncont[0:3,1::2,:1, 1], axis=(0,2,-1)),'--', c=cols2[5],marker='^',label=r':math:`\sigma=0.5`' ,zorder=5,lw=2.8,markersize=6,markeredgecolor='#4f4949')
plt.text(0.29, 0.49,"%.2f"%count_perc[0, 0 ,1], fontsize =8 , color = "grey")
plt.text(0.97, 0.57,"%.2f"%count_perc[1, 0 ,1], fontsize =8 , color = "grey")
plt.text(2.05, 0.44,"%.2f"%count_perc[2, 0 ,1], fontsize =8 , color = "grey")
plt.text(0.29, 0.28,"%.2f"%count_perc[0, 0 ,2], fontsize =8 , color = "grey")
plt.text(1.1, 0.37,"%.2f"%count_perc[1, 0 ,2], fontsize =8 , color = "grey")
plt.text(1.89, 0.28,"%.2f"%count_perc[2, 0 ,2], fontsize =8 , color = "grey")
plt.plot(Ks[1:-1:2], np.nanmean(onset_cont[0:3,1::2,:1, 2], axis=(0,2,-1)) , c=cols[1],marker='.',label=r':math:`\sigma=1.0`',lw=2.8 ,markersize=12,markeredgecolor='#4f4949')
plt.plot(Ks[1:-1:2], np.nanmean(onset_uncont[0:3,1::2,:1, 2], axis=(0,2,-1)),'--', c=cols2[0],marker='.' ,label=r':math:`\sigma=1.0`',lw=2.8,markersize=12,markeredgecolor='#4f4949')
#plt.yscale('log')
plt.ylabel(r'onset of synchrony :math:`t^{syn}`')
plt.xlabel(r'coupling :math:`J`')
ax6 = plt.gca()
ax6.tick_params(axis='both',which='both',direction='in', length=3, width=1,colors='#4f4949',zorder=3)
ax6.tick_params(bottom=True, top=True, left=True, right=True)
ax6.spines['top'].set_visible(True)
ax6.spines['right'].set_visible(True)
ax6.minorticks_on()
ax6.tick_params(axis='both',which='major',direction='in', length=3.5, width=1,colors='#4f4949',zorder=3)
ax6.tick_params(axis='both',which='minor',direction='in', length=2.5, width=0.5,colors='#4f4949',zorder=3)
ax6.tick_params(bottom=True, top=True, left=True, right=True)
ax6.tick_params(axis='both', which='minor', bottom=True, top=True, left=True, right=True)
plt.subplot(1,2,2)
plt.plot(Ks[1:-1:2], np.nanmean(cont_synched_durations_sum[0:3,1::2,:1, 1], axis=(0,2,-1)) , c=cols[5],marker='^',label=r':math:`\sigma=0.5`',zorder=2 ,lw=2.8,markersize=6,markeredgecolor='#4f4949')
plt.plot(Ks[1:-1:2], np.nanmean(uncont_synched_durations_sum[0:3,1::2,:1, 1], axis=(0,2,-1)),'--', c=cols2[5],marker='^',label=r':math:`\sigma=0.5`' ,zorder=2,lw=2.8,markersize=6,markeredgecolor='#4f4949')
plt.plot(Ks[1:-1:2], np.nanmean(cont_synched_durations_sum[0:3,1::2,:1, 2], axis=(0,2,-1)) , c=cols[1],marker='.',label=r':math:`\sigma=1.0`',lw=2.8 ,markersize=12,markeredgecolor='#4f4949')
plt.plot(Ks[1:-1:2], np.nanmean(uncont_synched_durations_sum[0:3,1::2,:1, 2], axis=(0,2,-1)),'--', c=cols2[0],marker='.' ,label=r':math:`\sigma=1.0`',lw=2.8,markersize=12,markeredgecolor='#4f4949')
plt.ylabel(':math:`\%` time synchronised\nafter :math:`t^{syn}`', multialignment='center',labelpad=-0.1) # "Mat\nTTp\n123"
plt.xlabel(r'coupling :math:`J`')
plt.yticks([0, 0.5,1])
ax7 = plt.gca()
ax7.tick_params(axis='both',which='both',direction='in', length=3, width=1,colors='#4f4949',zorder=3)
ax7.tick_params(bottom=True, top=True, left=True, right=True)
ax7.spines['top'].set_visible(True)
ax7.spines['right'].set_visible(True)
ax7.minorticks_on()
ax7.tick_params(axis='both',which='major',direction='in', length=3.5, width=1,colors='#4f4949',zorder=3)
ax7.tick_params(axis='both',which='minor',direction='in', length=2.5, width=0.5,colors='#4f4949',zorder=3)
ax7.tick_params(bottom=True, top=True, left=True, right=True)
ax7.tick_params(axis='both', which='minor', bottom=True, top=True, left=True, right=True)
legend = ax7.legend()
legend.get_frame().set_linewidth(1.8)
legend.get_frame().set_facecolor('white')
legend.get_frame().set_edgecolor('white')
handles, labels = ax7.get_legend_handles_labels()
handles = [ handles[0],handles[2] , handles[1], handles[3]]
labels= [ labels[0],labels[2] , labels[1], labels[3]]
if True:
ax7.legend(handles, labels, title=r':math:`\,`controlled :math:`\,\,` uncontrolled',
handletextpad=0.5, columnspacing=1,handlelength=0.45, bbox_to_anchor=[-0.65, 0.99],
loc=3, ncol=2, frameon=True,fontsize = 'small',shadow=None,framealpha =0,edgecolor ='#0a0a0a')
plt.setp(plt.gca().get_legend().get_title(), fontsize='10')
plt.subplots_adjust(wspace=0.325)#, hspace=0)
plt.savefig("Kuramoto_2_onset.png", bbox_inches='tight')
plt.savefig("Kuramoto_2_onset.pdf", bbox_inches='tight')
"""
plt.subplot(1,2,2)#,sharey=ax1)
# plt.plot(Ks[1:-1:2], np.nanmean(onset_cont[:,1::2,yi, 1], axis=(0,-1)) , c=cols[5],marker='^',label=r':math:`\sigma=0.5`',zorder=5 ,lw=2.8,markersize=6,markeredgecolor='#4f4949')
# plt.plot(Ks[1:-1:2], np.nanmean(onset_uncont[:,1::2,yi, 1], axis=(0,-1)),'--', c=cols2[5],marker='.',label=r':math:`\sigma=1.0`' ,zorder=5,lw=2.8,markersize=10,markeredgecolor='#4f4949')
plt.plot(Ks[1:-1:2], np.nanmean(onset_cont[0:3,1::2,0, 2], axis=(0,-1)) , c=cols[1],marker='^',label=r':math:`\sigma=0.5`',lw=2.8 ,markersize=6,markeredgecolor='#4f4949')
plt.plot(Ks[1:-1:2], np.nanmean(onset_uncont[0:3,1::2,0, 2], axis=(0,-1)),'--', c=cols2[2],marker='.' ,label=r':math:`\sigma=1.0`',lw=2.8,markersize=12,markeredgecolor='#4f4949')
plt.plot(Ks[1:-1:2], np.nanmean(onset_cont[0:3,1::2,1, 2], axis=(0,-1)) , c=cols[1],marker='d',label=r':math:`\sigma=0.5`',lw=2.8 ,markersize=6,markeredgecolor='#4f4949')
plt.plot(Ks[1:-1:2], np.nanmean(onset_uncont[0:3,1::2,1, 2], axis=(0,-1)),'--', c=cols2[2],marker='p' ,label=r':math:`\sigma=1.0`',lw=2.8,markersize=7,markeredgecolor='#4f4949')
#plt.yscale('log')
plt.ylabel(r'onset of synchrony :math:`t^{syn}`')
plt.xlabel(r'coupling :math:`J`')
"""
ki=1
yi=0
gi=1
bi=1
for repi in range(1):
positions = np.where( Rttnon[bii, ki, yi, gi, 1:-2,repi] >=threshold)[0]
print(positions)
if positions.size ==0:
onset_uncont[bi,ki,yi,gi,repi] = np.nan
else:
diffpos = np.diff(positions)
positions_consecutive = np.where(diffpos < 5)[0]
if positions_consecutive.size ==0:
onset_uncont[bi,ki,yi,gi,repi] = np.nan
else:
onset_uncont[bi,ki,yi,gi,repi] = timegrid[positions[positions_consecutive[0]] ]
print(onset_uncont[bi,ki,yi,gi,repi])
# print()
# print()
Total running time of the script: ( 0 minutes 0.000 seconds)
