<filename>python-code/analysis.py
# -*- coding: utf-8 -*-

''' --------- Imports ----------- '''



import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
import scipy as sp
import scipy.stats as st
from scipy.io import loadmat
from matplotlib.patches import Rectangle
import matplotlib.lines as mlines



''' --------- Upload data from .mat files ----------- '''



ornstein_t_s = loadmat('results_ornstein_sigma_tau.mat')
ornstein_c_t = loadmat('results_ornstein_ctau.mat')
ornstein_c_s = loadmat('results_ornstein_csigma.mat')
pink         = loadmat('results_pink.mat')
white        = loadmat('results_white.mat')

# keys of smaller dictionaries
list_matrices = ['DurMatrix','PercMatrix','WTAmatrix','CVmatrix']
list_sd       = ['DurSDmatrix','PercSDmatrix','WTAsdMatrix','CVsdMatrix']
new_keys      = ['DomDur','PercTime','DomIndex','CV']

# big dictionaries
main_keys     = ['OU_ts','OU_ct','OU_cs','Pink_cs','White_cs']
old_dicts     = [ornstein_t_s,ornstein_c_t,ornstein_c_s,pink,white]

data = {}
sd_data = {}
for i in range(0,5):
  old_dict = old_dicts[i]

  var = {new_keys[j]: old_dict[list_matrices[j]] for j in range(0,4)}
  data[main_keys[i]] = var

  var = {new_keys[j]: old_dict[list_sd[j]] for j in range(0,4)}
  sd_data[main_keys[i]] = var

"""Two data structures: `data` and `data_sd` --- one for the mean values 
  of the metrics and another for the standard deviations. 

*   Main keys (types of simulation):
      * 'OU_ts'
      * 'OU_ct'
      * 'OU_cs'
      * 'Pink_cs'
      * 'White_cs'

*   Second keys (metrics):
      * 'DomDur'
      * 'PercTime'
      * 'DomIndex'
      * 'CV'

**Example:** `data['OU_ts']['DomDur']` accesses the dominance duration matrix 
  of the Ornstein-Uhlenbeck simulation with varying tau and sigma. 
  `sd_data['OU_ts']['DomDur']` accesses the standard deviation of those results."""
  
  

''' --------- Data for graphs ----------- '''
  
  

# for Fig1

heat_ou_cs = np.flipud(data['OU_cs']['PercTime'])
heat_pink  = np.flipud(data['Pink_cs']['PercTime'])
heat_white = np.flipud(data['White_cs']['PercTime'])

heat_dur_ou    = np.flipud(data['OU_cs']['DomDur'])
heat_dur_pink  = np.flipud(data['Pink_cs']['DomDur'])
heat_dur_white = np.flipud(data['White_cs']['DomDur'])

# for Fig2A

hist_ou    = heat_ou_cs.ravel()
hist_pink  = heat_pink.ravel()
hist_white = heat_white.ravel()

# for Fig2B

heat_cv_ou   = np.flipud(data['OU_cs']['CV'])
heat_cv_pink = np.flipud(data['Pink_cs']['CV'])

hist_ou_cv   = heat_cv_ou.ravel()
hist_pink_cv = heat_cv_pink.ravel()

hist_ou_dur   = heat_dur_ou.ravel()
hist_pink_dur = heat_dur_pink.ravel()

mask_ou = (heat_dur_ou >= 1) & (heat_ou_cs >= 20)
mask_pi = (heat_dur_pink >= 1) & (heat_pink >= 20)

hist_ou_cv_m    = heat_cv_ou[mask_ou].ravel()
hist_pink_cv_m  = heat_cv_pink[mask_pi].ravel()

hist_ou_dur_m    = heat_dur_ou[mask_ou].ravel()
hist_pink_dur_m  = heat_dur_pink[mask_pi].ravel()

heat_cv_white     = np.flipud(data['White_cs']['CV'])
mask_white        = heat_white >=20 
hist_white_cv_m   = heat_cv_white[mask_white].ravel()
hist_white_dur_m  = heat_dur_white[mask_white].ravel()

# for Fig3

heat_ou_ts     = np.flipud(data['OU_ts']['PercTime'])
heat_ou_ts_sd  = np.flipud(sd_data['OU_ts']['PercTime'])
heat_dur_ou_ts = np.flipud(data['OU_ts']['DomDur'])
heat_ou_ts_sd_norm = np.divide(heat_ou_ts_sd,heat_ou_ts)*100
heat_ou_ts_sd_norm = np.nan_to_num(heat_ou_ts_sd_norm)

# for Fig6B

ou_time     = loadmat('OU_time.mat')
pink_time   = loadmat('pink_time.mat')
white_time  = loadmat('white_time.mat')

ou_time_t = ou_time['t'][0,:]
ou_time_a = ou_time['FA'][0,:]
ou_time_b = ou_time['FB'][0,:]

pink_time_t = pink_time['t'][0,:]
pink_time_a = pink_time['FA'][0,:]
pink_time_b = pink_time['FB'][0,:]

white_time_t = white_time['t'][0,:]
white_time_a = white_time['FA'][0,:]
white_time_b = white_time['FB'][0,:]



''' --------- Figure 1 ----------- '''



fig, (ax1,ax2,ax3,ax4) = plt.subplots(1,4, gridspec_kw={'width_ratios': [1,1,1,0.05]})
fig.set_size_inches(35,10)

size_titles = 32
size_labels = 24
size_pad = 25

# heatmap of PDT (proportion of dominance time)

sns.heatmap(heat_ou_cs,ax=ax1,vmin=0,vmax=100,cmap=sns.color_palette("Blues"),
                 square=True,cbar=None)

sns.heatmap(heat_pink,ax=ax2,vmin=0,vmax=100,cmap=sns.color_palette("Blues"),
                 square=True,cbar=None,yticklabels=[])

sns.heatmap(heat_white,ax=ax3,vmin=0,vmax=100,cmap=sns.color_palette("Blues"),
                 square=True,yticklabels=[],cbar_ax=ax4)

# formatting and labels

cbar = ax3.collections[0].colorbar
cbar.set_ticks([0,50,100])
cbar.set_ticklabels([0,50,100])
cbar.ax.tick_params(labelsize=size_labels)
cbar.ax.set_ylabel('Dominance Time (%)',fontsize=size_labels)

# contour of mean percept duration

X = np.linspace(0,50,50)
Y = np.linspace(0,50,50)
X, Y = np.meshgrid(X,Y)

for a,h,t in zip([ax1,ax2,ax3],[heat_dur_ou,heat_dur_pink,heat_dur_white],
                 ['Ornstein-Uhlenbeck Noise','Pink Noise','White Gaussian Noise']):
  # apply gaussian filter
  c = sp.ndimage.filters.gaussian_filter(h,sigma=1.2,mode='nearest',truncate=4)
  # plot
  a.contour(X,Y,c,levels=[1],colors=['w'],linewidths=3)
  a.contourf(X,Y,c,levels=[0,1],colors=['w'],alpha=0.5)
  # formatting and labels
  a.set_xlabel(r'Noise Intensity, $\sigma$',fontsize=size_labels,labelpad=-10)
  a.set_xticks([0,50])
  a.set_xticklabels([0.005,0.25],fontsize=size_labels)
  a.tick_params(axis='x', rotation=0, pad = 10)
  a.set_title(t,{'fontsize':size_titles,'fontweight':'normal'},pad=size_pad)

# formatting and labels

ax1.annotate(r'$\bar{D}=1$',xy=(10,36),fontsize=size_labels,fontweight='bold',color='w')
ax2.annotate(r'$\bar{D}=1$',xy=(9,34),fontsize=size_labels,fontweight='bold',color='w')

ax1.set_ylabel(r'Input Contrast, c',fontsize=size_labels,labelpad=-25)
ax1.set_yticks([0,50])
ax1.set_yticklabels([1,0.02],fontsize=size_labels)

plt.savefig('fig1.pdf',bbox_inches='tight')



''' --------- Figure 2A ----------- '''



fig, ax = plt.subplots(1,1)
fig.set_size_inches(10,10)

num_bins = 500

ax.hist(hist_white,bins=num_bins,range=(0,100),color='tab:orange',alpha = 0.5,cumulative=True)
ax.hist(hist_pink,bins=num_bins,range=(0,100),color='tab:blue',alpha = 0.5,cumulative=True)
ax.hist(hist_ou,bins=num_bins,range=(0,100),color='tab:green',alpha = 0.4,cumulative=True)

size_labels = 16

ax.legend(['White Gaussian Noise','Pink Noise','Ornstein-Uhlenbeck Noise'],fontsize=size_labels,loc='lower right')

ax.set_xlabel(r'Dominance Time (%)',fontsize=size_labels,labelpad=size_pad)
ax.set_xticks(np.linspace(0,100,5))
ax.set_xticklabels(np.linspace(0,100,5,dtype=int),fontsize=size_labels)
ax.set_xlim([0,100])
ax.tick_params('y',labelsize=size_labels)
ax.set_ylabel(r'Cumulative Number of Cases',fontsize=size_labels,labelpad=size_pad)
ax.set_ylim([0,2500])

ax.spines["top"].set_visible(False)
ax.spines["right"].set_visible(False)
ax.spines['bottom'].set_color('#c4c4c4')
ax.spines['left'].set_color('#c4c4c4')

plt.savefig('fig2A.pdf',bbox_inches='tight')



''' --------- Figure 2B ----------- '''



fig, ax = plt.subplots(1,1)
fig.set_size_inches(10,10)

# mean values and standard deviation of CV and D for all noise processes

print('\nMeanDuration white: \t %.4f +- %.4f'%(np.mean(hist_white_dur_m),np.std(hist_white_dur_m)))
print('MeanDuration OU: \t %.4f +- %.4f'%(np.mean(hist_ou_dur_m),np.std(hist_ou_dur_m)))
print('MeanDuration pink: \t %.4f +- %.4f'%(np.mean(hist_pink_dur_m),np.std(hist_pink_dur_m)))

print('\nCV white: \t %.4f +- %.4f'%(np.mean(hist_white_cv_m),np.std(hist_white_cv_m)))
print('CV OU: \t\t %.4f +- %.4f'%(np.mean(hist_ou_cv_m),np.std(hist_ou_cv_m)))
print('CV pink: \t %.4f +- %.4f\n'%(np.mean(hist_pink_cv_m),np.std(hist_pink_cv_m)))

size_labels = 16

ax.spines["top"].set_visible(False)
ax.spines["right"].set_visible(False)
ax.spines['bottom'].set_color('#c4c4c4')
ax.spines['left'].set_color('#c4c4c4')
ax.tick_params(axis='both',labelsize=size_labels)


xpos_dur = [0.85,1.85]
xpos_cv  = [1.15,2.15]
wid_violin = 0.2

# violin plot durations

mins = np.array([np.min(hist_ou_dur_m),np.min(hist_pink_dur_m)])
maxes = np.array([np.max(hist_ou_dur_m),np.max(hist_pink_dur_m)])
means = np.array([np.mean(hist_ou_dur_m),np.mean(hist_pink_dur_m)])
std = np.array([np.std(hist_ou_dur_m),np.std(hist_pink_dur_m)])

parts = ax.violinplot([hist_ou_dur_m,hist_pink_dur_m],xpos_dur,widths=wid_violin,showextrema=False)
for pc in parts['bodies']:
    pc.set_facecolor('tab:blue')

ax.errorbar(np.array(xpos_dur), means, [means - mins, maxes - means], fmt='.k', ecolor='darkgray', elinewidth=1)
ax.errorbar(np.array(xpos_dur), means, std, fmt='ok', elinewidth=4, ecolor='tab:blue')

# violin plot cv

mins = np.array([np.min(hist_ou_cv_m),np.min(hist_pink_cv_m)])
maxes = np.array([np.max(hist_ou_cv_m),np.max(hist_pink_cv_m)])
means = np.array([np.mean(hist_ou_cv_m),np.mean(hist_pink_cv_m)])
std = np.array([np.std(hist_ou_cv_m),np.std(hist_pink_cv_m)])

parts = ax.violinplot([hist_ou_cv_m,hist_pink_cv_m],xpos_cv,widths=wid_violin,showextrema=False)
for pc in parts['bodies']:
    pc.set_facecolor('tab:green')

ax.errorbar(np.array(xpos_cv), means, [means - mins, maxes - means], fmt='.k', ecolor='darkgray', elinewidth=1)
ax.errorbar(np.array(xpos_cv), means, std, fmt='ok', elinewidth=4, ecolor='tab:green')

# formatting and labels

ax.set_xlim([0.5,2.5])
ax.set_xticks([1,2])
ax.set_xticklabels(['Ornstein-Uhlenbeck Noise','Pink Noise'])
ax.tick_params(axis='x', pad=20)
ax.set_xlabel(' ')
ax.set_ylabel(' ')
ax.set_ylim([0,5.2])
ax.set_yticks([1,2,3,4,5])
ax.set_yticklabels([1,2,3,4,5])

blue_line = mlines.Line2D([], [], color='tab:blue',lw=3)
green_line = mlines.Line2D([], [], color='tab:green',lw=3)
ax.legend(handles=[blue_line,green_line],labels=['Mean Percept Duration (s)','Coefficient of Variation of Percept Durations'],loc='upper center',markerscale=0,fontsize=size_labels)


# inset - CV
left, bottom, width, height = [0.42, 0.54, 0.2, 0.2]
ax2 = fig.add_axes([left, bottom, width, height])

ax2.errorbar(np.array([1.05,2.05]), means, [means - mins, maxes - means], fmt='.k', ecolor='lightgray', elinewidth=2)
ax2.errorbar(np.array([1.05,2.05]), means, std, fmt='ok', elinewidth=4, ecolor='tab:green')

# formatting and labels

size_labels_in = 10

ax2.set_xticks([1.05,2.05])
ax2.set_xticklabels([])
ax2.set_xlim([0.55,2.55])
ax2.set_ylim([0.4,1.2])
ax2.set_yticks([0.4,0.6,0.8,1.0,1.2])
ax2.set_yticklabels([0.4,0.6,0.8,1.0,1.2],fontsize=size_labels_in)

ax2.spines["top"].set_color('#c4c4c4')
ax2.spines["right"].set_color('#c4c4c4')
ax2.spines['bottom'].set_color('#c4c4c4')
ax2.spines['left'].set_color('#c4c4c4')


plt.savefig('fig2B.pdf',bbox_inches='tight')



''' --------- Figure 3 ----------- '''



fig, (ax1,ax2) = plt.subplots(1,2)
fig.set_size_inches(15,10)

size_labels = 16

# heatmap of PDT - OU noise with varying sigma and tau

sns.heatmap(heat_ou_ts,ax=ax1,vmin=0,vmax=100,cmap=sns.color_palette("Blues"),
                 square=True,cbar_kws={'use_gridspec':False,'location':'top',
                                       'shrink':0.8})

# formatting and labels

ax1.set_xlabel(r'Noise Intensity, $\sigma$',fontsize=size_labels,labelpad=-10)
ax1.set_xticks([0,50])
ax1.set_xticklabels([0.002,0.1],fontsize=size_labels)
ax1.tick_params(axis='x', rotation=0, pad=10)

ax1.set_ylabel(r'Noise Correlation Time, $\tau$ (ms)',fontsize=size_labels,labelpad=10)
ax1.set_yticks([0,35,40,50])
ax1.set_yticklabels([1000,300,200,20],fontsize=size_labels)
ax1.tick_params(axis='y', rotation=0)

cbar = ax1.collections[0].colorbar
cbar.set_ticks([0,50,100])
cbar.set_ticklabels([0,50,100])
cbar.ax.tick_params(labelsize=size_labels)
cbar.ax.set_xlabel('Dominance Time (%)',fontsize=size_labels,labelpad=15)

# contour of mean percept duration

contour_dur_ou = sp.ndimage.filters.gaussian_filter(heat_dur_ou_ts,sigma=1.2,mode='nearest',truncate=4)
X = np.linspace(0,50,50)
Y = np.linspace(0,50,50)
X, Y = np.meshgrid(X,Y)
ax1.contour(X,Y,contour_dur_ou,levels=[1],colors=['w'],linewidths=3)
ax1.contourf(X,Y,contour_dur_ou,levels=[0,1],colors=['w'],alpha=0.5)

# mark parameters used by Said and Heeger

ax1.annotate(r'$\bar{D}=1$',xy=(25,38),fontsize=size_labels,fontweight='bold',color='w')

oriX, oriY = 24, 10
lil_box = Rectangle((oriX, oriY), 1, 1,alpha=1,fc=None,ec='tab:red',lw=3)
ax1.add_patch(lil_box)

# heatmap of standard deviation of PDT

sns.heatmap(heat_ou_ts_sd_norm,vmin=0,vmax=180,ax=ax2,cmap=sns.color_palette("Blues"),
                 square=True,cbar_kws={'use_gridspec':False,'location':'top',
                                       'shrink':0.8})

# formatting and labels

ax2.set_xlabel(r'Noise Intensity, $\sigma$',fontsize=size_labels,labelpad=-10)
ax2.set_xticks([0,50])
ax2.set_xticklabels([0.002,0.1],fontsize=size_labels)
ax2.tick_params(axis='x',rotation=0,pad=10)

ax2.set_yticks([])
ax2.set_yticklabels([])

cbar = ax2.collections[0].colorbar
cbar.set_ticks([0,90,180])
cbar.set_ticklabels([0,90,180])
cbar.ax.tick_params(labelsize=size_labels)
cbar.ax.set_xlabel('Standard Deviation of Dominance Time (%)',fontsize=size_labels,labelpad=15)

plt.savefig('fig3.pdf',bbox_inches='tight')



''' --------- Figure 4 ----------- '''



fig, (ax1,ax2,ax3) = plt.subplots(1,3,gridspec_kw={'width_ratios': [1,1,1]})
fig.set_size_inches(30,10)

# define slices

low_c_ou    =    heat_dur_ou[45,:]
low_c_pink  =  heat_dur_pink[45,:]
low_c_white = heat_dur_white[45,:]

med_c_ou    =    heat_dur_ou[25,:]
med_c_pink  =  heat_dur_pink[25,:]
med_c_white = heat_dur_white[25,:]

high_c_ou    =    heat_dur_ou[5,:]
high_c_pink  =  heat_dur_pink[5,:]
high_c_white = heat_dur_white[5,:]

# formatting and labels

size_labels = 20
size_titles = 22

for a in [ax1,ax2,ax3]:
  a.set_xlim(0,0.255)
  a.set_ylim(0,3.5)
  a.set_xlabel(r'Noise intensity, $\sigma$',fontsize=size_labels,labelpad=size_pad)
  a.set_xticks(np.linspace(0,0.25,6))
  a.set_xticklabels(np.linspace(0,0.25,6).round(3),fontsize=size_labels)
  a.set_yticks(np.linspace(0,3.5,3))
  a.set_yticklabels(np.linspace(0,3.5,3).round(2),fontsize=size_labels)
  a.tick_params(axis='x',rotation=0,pad=10)

  a.spines["top"].set_visible(False)
  a.spines["right"].set_visible(False)
  a.spines['bottom'].set_color('#c4c4c4')
  a.spines['left'].set_color('#c4c4c4')

x = np.arange(0.005,0.2505,0.005)

# plot low contrast

ax1.plot(x,low_c_ou,linewidth=3,color='tab:green')
ax1.plot(x,low_c_pink,linewidth=3,color='tab:blue')
ax1.plot(x,low_c_white,linewidth=3,color='tab:orange')
ax1.set_title(r'c = 0.1',fontsize=size_titles)
ax1.set_ylabel(r'Mean Percept Duration (s)',fontsize=size_labels)

ax1.fill_between(x,low_c_ou,0,color='tab:green',alpha=0.3)
ax1.fill_between(x,low_c_pink,0,color='tab:blue',alpha=0.5)
ax1.fill_between(x,low_c_white,0,color='tab:orange',alpha=0.6)

# plot medium contrast

ax2.plot(x,med_c_ou,linewidth=3,color='tab:green')
ax2.plot(x,med_c_pink,linewidth=3,color='tab:blue')
ax2.plot(x,med_c_white,linewidth=3,color='tab:orange')
ax2.set_title(r'c = 0.5',fontsize=size_titles)
ax2.legend(['Ornstein-Uhlenbeck Noise','Pink Noise','White Gaussian Noise'],fontsize=size_titles,loc='upper center',
           bbox_to_anchor=(0.5,-0.15),ncol=3)

ax2.fill_between(x,med_c_ou,0,color='tab:green',alpha=0.3)
ax2.fill_between(x,med_c_pink,0,color='tab:blue',alpha=0.5)
ax2.fill_between(x,med_c_white,0,color='tab:orange',alpha=0.6)

# plot high contrast

ax3.plot(x,high_c_ou,linewidth=3,color='tab:green')
ax3.plot(x,high_c_pink,linewidth=3,color='tab:blue')
ax3.plot(x,high_c_white,linewidth=3,color='tab:orange')
ax3.set_title(r'c = 0.9',fontsize=size_titles)

ax3.fill_between(x,high_c_ou,0,color='tab:green',alpha=0.3)
ax3.fill_between(x,high_c_pink,0,color='tab:blue',alpha=0.5)
ax3.fill_between(x,high_c_white,0,color='tab:orange',alpha=0.6)

plt.savefig('fig4.pdf',bbox_inches='tight')



''' --------- Figure 5 ----------- '''



fig, (ax1,ax2,ax3) = plt.subplots(1,3,gridspec_kw={'width_ratios': [1,1,1]})
fig.set_size_inches(30,10)

# define slices

low_s_ou    = np.flipud(heat_ou_cs[:,2])
low_s_pink  =  np.flipud(heat_pink[:,2])
low_s_white = np.flipud(heat_white[:,2])

med_s_ou    = np.flipud(heat_ou_cs[:,24])
med_s_pink  =  np.flipud(heat_pink[:,24])
med_s_white = np.flipud(heat_white[:,24])

high_s_ou    = np.flipud(heat_ou_cs[:,49])
high_s_pink  =  np.flipud(heat_pink[:,49])
high_s_white = np.flipud(heat_white[:,49])

# formatting and labels

size_labels = 20
size_titles = 22

for a in [ax1,ax2,ax3]:
  a.set_xlim(0,1.02)
  a.set_ylim(0,100)
  a.set_xlabel(r'Input Contrast, c',fontsize=size_labels,labelpad=size_pad)
  a.set_xticks(np.linspace(0,1,5))
  a.set_xticklabels(np.linspace(0,1,5).round(3),fontsize=size_labels)
  a.set_yticks(np.linspace(0,100,5))
  a.set_yticklabels(np.linspace(0,100,5,dtype=int),fontsize=size_labels)
  a.tick_params(axis='x',rotation=0,pad=10)

  a.spines["top"].set_visible(False)
  a.spines["right"].set_visible(False)
  a.spines['bottom'].set_color('#c4c4c4')
  a.spines['left'].set_color('#c4c4c4')

x = np.arange(0.02,1.02,0.02)

# plot low noise

ax1.plot(x,low_s_ou,linewidth=3,color='tab:green')
ax1.plot(x,low_s_pink,linewidth=3,color='tab:blue')
ax1.plot(x,low_s_white,linewidth=3,color='tab:orange')
ax1.set_title(r'$\sigma$ = 0.06',fontsize=size_titles)
ax1.set_ylabel(r'Dominance Time (%)',fontsize=size_labels)

ax1.fill_between(x,low_s_ou,0,color='tab:green',alpha=0.3)
ax1.fill_between(x,low_s_pink,0,color='tab:blue',alpha=0.5)
ax1.fill_between(x,low_s_white,0,color='tab:orange',alpha=0.6)

# plot medium noise

ax2.plot(x,med_s_ou,linewidth=3,color='tab:green')
ax2.plot(x,med_s_pink,linewidth=3,color='tab:blue')
ax2.plot(x,med_s_white,linewidth=3,color='tab:orange')
ax2.set_title(r'$\sigma$ = 0.125',fontsize=size_titles)
ax2.legend(['Ornstein-Uhlenbeck Noise','Pink Noise','White Gaussian Noise'],fontsize=size_titles,loc='upper center',
           bbox_to_anchor=(0.5,-0.15),ncol=3)

ax2.fill_between(x,med_s_ou,0,color='tab:green',alpha=0.3)
ax2.fill_between(x,med_s_pink,0,color='tab:blue',alpha=0.5)
ax2.fill_between(x,med_s_white,0,color='tab:orange',alpha=0.6)

# plot high noise

ax3.plot(x,high_s_ou,linewidth=3,color='tab:green')
ax3.plot(x,high_s_pink,linewidth=3,color='tab:blue')
ax3.plot(x,high_s_white,linewidth=3,color='tab:orange')
ax3.set_title(r'$\sigma$ = 0.25',fontsize=size_titles)

ax3.fill_between(x,high_s_ou,0,color='tab:green',alpha=0.3)
ax3.fill_between(x,high_s_pink,0,color='tab:blue',alpha=0.5)
ax3.fill_between(x,high_s_white,0,color='tab:orange',alpha=0.6)

plt.savefig('fig5.pdf',bbox_inches='tight')



''' --------- Figure 6B ----------- '''



fig, (ax1,ax2,ax3) = plt.subplots(3,1,gridspec_kw={'height_ratios': [1,1,1]})
fig.set_size_inches(20,15)

size_labels = 16
size_titles = 24

for a in [ax1,ax2,ax3]:
  a.set_xlim(0,20000)
  a.set_ylim(0,1)
  a.set_ylabel(r'Firing Rate',fontsize=size_labels)
  a.set_xticks(np.linspace(0,20000,6))
  a.set_xticklabels([])
  a.set_yticks(np.linspace(0,1,3))
  a.set_yticklabels(np.linspace(0,1,3).round(3),fontsize=size_labels)

  a.spines["top"].set_visible(False)
  a.spines["right"].set_visible(False)
  a.spines['bottom'].set_color('#c4c4c4')
  a.spines['left'].set_color('#c4c4c4')

# plot time course Ornstein-Uhlenbeck noise

ax1.plot(ou_time_t,ou_time_a,linewidth=3,color='tab:orange')
ax1.plot(ou_time_t,ou_time_b,linewidth=3,color='tab:blue')
ax1.set_title(r'Ornstein-Uhlenbeck Noise',fontsize=size_titles)
ax1.legend(['Binocular Layer A','Binocular Layer B'],fontsize=size_labels)
ax1.fill_between(ou_time_t, ou_time_a, ou_time_b, where=(ou_time_a > ou_time_b), color='tab:orange', alpha=0.3)
ax1.fill_between(ou_time_t, ou_time_a, ou_time_b, where=(ou_time_a < ou_time_b), color='tab:blue', alpha=0.3)

# plot time course pink noise

ax2.plot(pink_time_t,pink_time_a,linewidth=3,color='tab:orange')
ax2.plot(pink_time_t,pink_time_b,linewidth=3,color='tab:blue')
ax2.set_title(r'Pink Noise',fontsize=size_titles)
ax2.legend(['Binocular Layer A','Binocular Layer B'],fontsize=size_labels)
ax2.fill_between(pink_time_t, pink_time_a, pink_time_b, where=(pink_time_a > pink_time_b), color='tab:orange', alpha=0.3)
ax2.fill_between(pink_time_t, pink_time_a, pink_time_b, where=(pink_time_a < pink_time_b), color='tab:blue', alpha=0.3)

# plot time course white noise

ax3.plot(white_time_t,white_time_a,linewidth=3,color='tab:orange')
ax3.plot(white_time_t,white_time_b,linewidth=3,color='tab:blue')
ax3.legend(['Binocular Layer A','Binocular Layer B'],fontsize=size_labels)
ax3.set_title(r'White Gaussian Noise',fontsize=size_titles)
ax3.set_xlabel(r'Time (s)',fontsize=size_labels)
ax3.set_xticklabels(np.linspace(0,20,6,dtype=int),fontsize=size_labels)
ax3.fill_between(white_time_t, white_time_a, white_time_b, where=(white_time_a > white_time_b), color='tab:orange', alpha=0.3)
ax3.fill_between(white_time_t, white_time_a, white_time_b, where=(white_time_a < white_time_b), color='tab:blue', alpha=0.3)

plt.savefig('fig6B.pdf',bbox_inches='tight')

