# -*- coding: utf-8 -*-
"""
Created on Fri Mar 12 11:35:33 2021

@author: User
"""

import numpy as np
import matplotlib.pyplot as plt 
import scipy.io as sio
from scipy.stats import ranksums, pearsonr, sem, wilcoxon
import addcopyfighandler
import matplotlib.cm as cm
from brian2 import ms, psiemens
import matplotlib.colors as nc
import matplotlib as mpl
from mpl_toolkits.axes_grid1 import make_axes_locatable
from scipy.interpolate import interp2d

mpl.rcParams['font.sans-serif'] = "Arial"
mpl.rcParams['font.family'] = "sans-serif"
mpl.rcParams['font.size'] = 8
mpl.rcParams["xtick.direction"] = "in"
mpl.rcParams["ytick.direction"] = "in"
mpl.rcParams["legend.markerscale"] = 0.9

def simulations(N_units,coefC=0,d_LF=0):
    # run 6 simulations varying context and probes
    # outputs:
    # spike counts of probes
    # firing rate during maskers  
    
    probe=np.repeat(range(1,N_units+1),20)
    masker=np.repeat(range(1,N_units+1),20)
    
    # SIMULATION 1: navigation - after silence
    context=0
    exp=1
    exec(open("./input.py").read())  
    exec(open("./model.py").read()) 
    col=spike_counts(ac,onset,context)
    probe=np.column_stack((probe, col))
    
    # SIMULATION 2: communication - after silence
    exp=2
    exec(open("./input.py").read())  
    exec(open("./model.py").read()) 
    col=spike_counts(ac,onset,context)
    probe=np.column_stack((probe, col))
    
    # SIMULATION 3: navigation - navigation context
    context=1
    exp=1
    exec(open("./input.py").read())  
    exec(open("./model.py").read()) 
    col=spike_counts(ac,onset,context)
    probe=np.column_stack((probe, col))
    col=firing_context(ac,onset,context)
    colLF=firing_context(sp_low,onset,context)
    colHF=firing_context(sp_high,onset,context)
    masker=np.column_stack((masker, col))
    masker=np.column_stack((masker, colLF))
    masker=np.column_stack((masker, colHF))
    
    # SIMULATION 4: communication - navigation context
    exp=2
    exec(open("./input.py").read())  
    exec(open("./model.py").read()) 
    col=spike_counts(ac,onset,context)
    probe=np.column_stack((probe, col))
    
    # SIMULATION 5: navigation - communication context
    context=2
    exp=1
    exec(open("./input.py").read())  
    exec(open("./model.py").read()) 
    col=spike_counts(ac,onset,context)
    probe=np.column_stack((probe, col))
    col=firing_context(ac,onset,context)
    colLF=firing_context(sp_low,onset,context)
    colHF=firing_context(sp_high,onset,context)
    masker=np.column_stack((masker, col))
    masker=np.column_stack((masker, colLF))
    masker=np.column_stack((masker, colHF))
    
    # SIMULATION 6: communication - communication context
    exp=2
    exec(open("./input.py").read())  
    exec(open("./model.py").read()) 
    col=spike_counts(ac,onset,context)
    probe=np.column_stack((probe, col))
    
    return probe,masker


def spike_counts(ac,onset,context):
    # spike counts 50 ms after probe onset per trial
    spike_trains = ac.spike_trains()
    L=len(spike_trains)
    
    output=[]
    for u in range(L):
        spikes_i = np.asarray(spike_trains[u])*1000
        spikes_i = spikes_i[spikes_i>onset[context]] 
        output.append(len(spikes_i))

    return output
    

def firing_context(ac,onset,context):
    # firing rate (#sp/sec) during the context sequences per trial
    preT=onset[0]
    dur=[1.3311*1000, 1.9266*1000]
    spike_trains = ac.spike_trains()
    L=len(spike_trains)
    
    output=[]
    for u in range(L):
        spikes_i = np.asarray(spike_trains[u])*1000
        spikes_i = spikes_i[spikes_i>preT] 
        spikes_i = spikes_i[spikes_i<preT+dur[context-1]]
        output.append(len(spikes_i))   #/(dur[context-1]/1000)

    return output

def context_effect(sp):
    # calculate context effect (c.e.) for each combination per unit 
    trials=20
    N_units=int(len(sp)/trials)
    # separate spike counts per protocol
    nav_iso=sp[:,1]
    com_iso=sp[:,2]
    nav_c1=sp[:,3]
    com_c1=sp[:,4]
    nav_c2=sp[:,5]
    com_c2=sp[:,6]
    # calculate mean across trials 
    nav_iso=np.mean(nav_iso.reshape(N_units,trials),axis=1)
    com_iso=np.mean(com_iso.reshape(N_units,trials),axis=1)
    nav_c1=np.mean(nav_c1.reshape(N_units,trials),axis=1)
    com_c1=np.mean(com_c1.reshape(N_units,trials),axis=1)
    nav_c2=np.mean(nav_c2.reshape(N_units,trials),axis=1)
    com_c2=np.mean(com_c2.reshape(N_units,trials),axis=1)
    # calculate context with mean spikes counts 
    ce_nav_exp=(nav_c1-nav_iso)/(nav_c1+nav_iso)
    ce_nav_unexp=(com_c1-com_iso)/(com_c1+com_iso)
    ce_com_exp=(com_c2-com_iso)/(com_c2+com_iso)
    ce_com_unexp=(nav_c2-nav_iso)/(nav_c2+nav_iso)
    
    output=np.column_stack((ce_nav_exp, ce_nav_unexp, ce_com_exp,ce_com_unexp))
    
    return output

def firing_rate_masker(sp):
    # average firing rate (#sp/sec) during the context sequences per unit
    trials=20
    N_units=int(len(sp)/trials)
    
    if len(sp[0])==3:
        # separate rate per masker
        nav=sp[:,1]
        com=sp[:,2]
        # separate rate counts per unit
        nav=np.mean(nav.reshape(N_units,trials),axis=1)
        com=np.mean(com.reshape(N_units,trials),axis=1)
        
        output=np.column_stack((nav, com))
        
    elif len(sp[0])==7:
        # separate rate per masker
        nav=sp[:,1]
        lf_nav=sp[:,2]
        hf_nav=sp[:,3]
        com=sp[:,4]
        lf_com=sp[:,5]
        hf_com=sp[:,6]
        
        # separate rate counts per unit
        nav=np.mean(nav.reshape(N_units,trials),axis=1)
        lf_nav=np.mean(lf_nav.reshape(N_units,trials),axis=1)
        hf_nav=np.mean(hf_nav.reshape(N_units,trials),axis=1)
        com=np.mean(com.reshape(N_units,trials),axis=1)
        lf_com=np.mean(lf_com.reshape(N_units,trials),axis=1)
        hf_com=np.mean(hf_com.reshape(N_units,trials),axis=1)
        
        output=np.column_stack((nav,lf_nav,hf_nav,com,lf_com,hf_com))
    
    return output

def one_simulation(N_units,context=0):
    # run 2 simulations varying context
    # outputs: raster plot and synaptic depression

    # SIMULATION 1: navigation -  context
    exp=1
    exec(open("./input.py").read())  
    exec(open("./model.py").read()) 
    
    mpl.rcParams["xtick.direction"] = "out"
    mpl.rcParams["ytick.direction"] = "out"

    fig, ax = plt.subplots(3,2,sharex=True,sharey=True,figsize=(4,2.5))
    
    ax[0,0].scatter(sp_low.t, sp_low.i,s=0.5,c='r',edgecolors='none')  
    ax[1,0].scatter(sp_high.t, sp_high.i,s=0.5,c='b',edgecolors='none')  
    ax[2,0].scatter(ac.t, ac.i,s=0.5,c='gray',edgecolors='none')  
    
    ax[0,0].set_ylabel('LF trials')
    ax[1,0].set_ylabel('HF trials')
    ax[2,0].set_ylabel('cortical trials')
    ax[2,0].set_xlabel('Time (s)')


    # average presynaptic adaptation across trials
    h=N_units*20
    dt_=0.1
    nav_presyn=mon.x_S
    nav_high=nav_presyn[h:]
    nav_low=nav_presyn[:h] 
    nav_pre_high=np.mean(nav_high,0)
    nav_pre_low=np.mean(nav_low,0)
    
    # average postsynaptic adaptation across trials
    nav_postsyn=pos.V_th*1000   # to mV 
    nav_post=np.mean(nav_postsyn,0)
    
    # plotting synaptic depression
    ax2 = ax[0,0].twinx()
    ax2.plot(mon.t, nav_pre_low,color='r',linewidth=0.75)
    ax2.set_ylim(0.1,1.01)

    for tl in ax2.get_yticklabels():
        tl.set_color('r')
        
    ax3 = ax[1,0].twinx()
    ax3.plot(mon.t, nav_pre_high,color='b',linewidth=0.75)
    ax3.set_ylim(0.1,1.01)

    for tl in ax3.get_yticklabels():
        tl.set_color('b')
    
    # plotting postsynaptic adaptation
    ax4 = ax[2,0].twinx()
    ax4.plot(pos.t, nav_post,color='k',linewidth=0.75)
    ax4.set_ylim(-50,-40)
    
    
    fig2, axz = plt.subplots(1,1,figsize=(1,1))
    axz.hist(ac.t/ms,bins=50,range=(onset[context],onset[context]++50),rwidth=0.8,color='k')
    axz.xaxis.set_visible(False)
    axz.yaxis.set_visible(False)
    axz.patch.set_facecolor('blue')
    axz.patch.set_alpha(0.2)
    axz.set_ylim(0,15)
    fig2.tight_layout()
    
    
    # SIMULATION 2: navigation - context
    exp=2
    exec(open("./input.py").read())  
    exec(open("./model.py").read()) 

    ax[0,1].scatter(sp_low.t, sp_low.i,s=0.5,c='r',edgecolors='none')  
    ax[1,1].scatter(sp_high.t, sp_high.i,s=0.5,c='b',edgecolors='none')  
    ax[2,1].scatter(ac.t, ac.i,s=0.5,c='gray',edgecolors='none')  
    ax[2,1].set_xlabel('Time (s)')
    
    # average presynaptic adaptation across trials
    h=N_units*20
    nav_presyn=mon.x_S
    nav_high=nav_presyn[h:]
    nav_low=nav_presyn[:h] 
    nav_pre_high=np.mean(nav_high,0)
    nav_pre_low=np.mean(nav_low,0)
    
    # average postsynaptic adaptation across trials
    com_postsyn=pos.V_th*1000   # to mV 
    com_post=np.mean(com_postsyn,0)
    
    # plotting synaptic depression
    ax4 = ax[0,1].twinx()
    ax4.plot(mon.t, nav_pre_low,color='r',linewidth=0.75)
    ax4.set_ylim(0.1,1.01)
    ax4.set_ylabel('strength\n' + r'LF$_{\rm syn}$', color='r')
    for tl in ax4.get_yticklabels():
        tl.set_color('r')
        
    ax5 = ax[1,1].twinx()
    ax5.plot(mon.t,nav_pre_high,color='b',linewidth=0.75)
    ax5.set_ylim(0.1,1.01)
    ax5.set_ylabel('strength\n' + r'HF$_{\rm syn}$', color='b')
   
    for tl in ax5.get_yticklabels():
        tl.set_color('b')
        
    # plotting postsynaptic adaptation
    ax6 = ax[2,1].twinx()
    ax6.plot(pos.t,com_post,color='k',linewidth=0.75)
    ax6.set_ylim(-50,-40) 
    
    ax6.set_ylabel('spiking\n' + 'threshold', color='k')
    
    ax = ax.flatten() 
    for axi in ax:
        axi.spines['top'].set_visible(False)
        axi.spines['right'].set_visible(False)
        
    ax2.spines['top'].set_visible(False)     
    ax3.spines['top'].set_visible(False)   
    ax4.spines['top'].set_visible(False)     
    ax5.spines['top'].set_visible(False)     
    
    fig.tight_layout(w_pad=3,h_pad=0.1)

    fig3, axe = plt.subplots(1,1,figsize=(1,1))
    axe.hist(ac.t/ms,bins=50,range=(onset[context],onset[context]++50),rwidth=0.8,color='k')
    axe.xaxis.set_visible(False)
    axe.yaxis.set_visible(False)
    axe.patch.set_facecolor('red')
    axe.patch.set_alpha(0.2)
    axe.set_ylim(0,15)
    fig3.tight_layout()
    
    return fig, fig2, fig3


def di_2parameters(sims):

    # import experimental data - cliff delta values for 45 units 
    mat_contents = sio.loadmat(r'E:\Users\User\Documents\model-context\cliff_data-416.mat')
    a=mat_contents['cliff_data']
    
    exp_data_iso = a[:,0]
    exp_data_nav = a[:,1]
    exp_data_com = a[:,2]
    
    del mat_contents
    del a

    parameters1=sims[0]
    parameters2=sims[1]
    sims = np.delete(sims,[0,1], 0)
    ntrials=len(sims[0])
    nx=len(parameters1)
    ny=len(parameters2)
    ntot=nx*ny
    N_units=ntrials/20
    
    # calculate cliff delta from simulations per unit 
    sims_iso=[]
    sims_nav=[]
    sims_com=[]
    stats_iso=[]
    stats_nav=[]
    stats_com=[]
    for s in range(ntot):
        s_i=sims[s]
        units_iso=[]
        units_nav=[]
        units_com=[]
        for u in range(int(N_units)): 
            sp_e=s_i[20*u:20*u+20,1]
            sp_d=s_i[20*u:20*u+20,2]
            c,size = cliffsDelta(sp_e, sp_d)  
            units_iso.append(c)
            sp_e=s_i[20*u:20*u+20,3]
            sp_d=s_i[20*u:20*u+20,4]
            c,size = cliffsDelta(sp_e, sp_d)  
            units_nav.append(c)
            sp_e=s_i[20*u:20*u+20,5]
            sp_d=s_i[20*u:20*u+20,6]
            c,size = cliffsDelta(sp_e, sp_d)  
            units_com.append(c)
        sims_iso.append(units_iso)    
        sims_nav.append(units_nav)   
        sims_com.append(units_com) 
        
        # check for fitting experimental data
        H,p=ranksums(units_iso, exp_data_iso)
        if p>0.05:   # equal distributions
            stats_iso.append(1)
        else:
            stats_iso.append(0) 
            
        H,p=ranksums(units_nav, exp_data_nav)
        if p>0.05:   # equal distributions
            stats_nav.append(1)
        else:
            stats_nav.append(0) 
            
        H,p=ranksums(units_com, exp_data_com)
        if p>0.05:   # equal distributions
            stats_com.append(1)
        else:
            stats_com.append(0) 

    # avearage cliff delta across units 
    av_iso=np.mean(sims_iso,1)
    av_nav=np.mean(sims_nav,1)
    av_com=np.mean(sims_com,1)
    
    # convert into 2D matrix results and stats 
    av_iso=np.reshape(av_iso,(nx,ny)).T
    av_nav=np.reshape(av_nav,(nx,ny)).T
    av_com=np.reshape(av_com,(nx,ny)).T
    stats_iso=np.reshape(stats_iso,(nx,ny)).T
    stats_nav=np.reshape(stats_nav,(nx,ny)).T
    stats_com=np.reshape(stats_com,(nx,ny)).T
    
    # plotting
    fig, ax = plt.subplots(1,3,sharex=True,sharey=True,figsize=(6.5,2.5))
    li=0.5      # min and max of colorbar

    min_x=min(parameters1)-((parameters1[1]-parameters1[0])/2)
    max_x=max(parameters1)+((parameters1[1]-parameters1[0])/2)
    min_y=min(parameters2)-((parameters2[1]-parameters2[0])/2)
    max_y=max(parameters2)+((parameters2[1]-parameters2[0])/2)
    
    x_iso=np.nonzero(stats_iso)[1]
    y_iso=np.nonzero(stats_iso)[0]   
    ax[0].imshow(av_iso,cmap='RdBu',origin='lower',vmin=-li,vmax=li,extent=[min_x,max_x,min_y,max_y])  
    ax[0].scatter(parameters1[x_iso],parameters2[y_iso], s=5, c='k', marker="*")
    
    x_nav=np.nonzero(stats_nav)[1]
    y_nav=np.nonzero(stats_nav)[0]   
    ax[1].imshow(av_nav,cmap='RdBu',origin='lower',vmin=-li,vmax=li,extent=[min_x,max_x,min_y,max_y])
    ax[1].scatter(parameters1[x_nav],parameters2[y_nav], s=5, c='k', marker="*")
    
    x_com=np.nonzero(stats_com)[1]
    y_com=np.nonzero(stats_com)[0]   
    ax[2].imshow(av_com,cmap='RdBu',origin='lower',vmin=-li,vmax=li,extent=[min_x,max_x,min_y,max_y])
    ax[2].scatter(parameters1[x_com],parameters2[y_com], s=5, c='k', marker="*")


    
    # set which combination satisfy all the contexts
    xy_iso=np.array(list(zip(x_iso,y_iso)))
    xy_nav=np.array(list(zip(x_nav,y_nav)))
    xy_com=np.array(list(zip(x_com,y_com)))
    combi = np.array([x for x in set(tuple(x) for x in xy_nav) & set(tuple(x) for x in xy_com)])
    combi = np.array([x for x in set(tuple(x) for x in combi) & set(tuple(x) for x in xy_iso)])

    ax[0].scatter(parameters1[combi[:,0]],parameters2[combi[:,1]], s=50, facecolors='none', edgecolors='k')
    ax[1].scatter(parameters1[combi[:,0]],parameters2[combi[:,1]], s=50, facecolors='none', edgecolors='k')
    ax[2].scatter(parameters1[combi[:,0]],parameters2[combi[:,1]], s=50, facecolors='none', edgecolors='k')
    
    
    return fig


def cliffsDelta(lst1, lst2, **dull):

    """Returns delta and true if there are more than 'dull' differences"""
    if not dull:
        dull = {'small': 0.147, 'medium': 0.33, 'large': 0.474} # effect sizes from (Hess and Kromrey, 2004)
    m, n = len(lst1), len(lst2)
    lst2 = sorted(lst2)
    j = more = less = 0
    for repeats, x in runs(sorted(lst1)):
        while j <= (n - 1) and lst2[j] < x:
            j += 1
        more += j*repeats
        while j <= (n - 1) and lst2[j] == x:
            j += 1
        less += (n - j)*repeats
    d = (more - less) / (m*n)
    size = lookup_size(d, dull)
    return d, size


def lookup_size(delta: float, dull: dict) -> str:
    """
    :type delta: float
    :type dull: dict, a dictionary of small, medium, large thresholds.
    """
    delta = abs(delta)
    if delta < dull['small']:
        return 'negligible'
    if dull['small'] <= delta < dull['medium']:
        return 'small'
    if dull['medium'] <= delta < dull['large']:
        return 'medium'
    if delta >= dull['large']:
        return 'large'


def runs(lst):
    """Iterator, chunks repeated values"""
    for j, two in enumerate(lst):
        if j == 0:
            one, i = two, 0
        if one != two:
            yield j - i, one
            i = j
        one = two
    yield j - i + 1, two
    
    
def context_effect_trial(sp):
    # calculate context effect (c.e.) for each combination per trial 
    trials=20
    N_units=int(len(sp)/trials)
    # separate spike counts per protocol
    nav_iso=sp[:,1]
    com_iso=sp[:,2]
    nav_c1=sp[:,3]
    com_c1=sp[:,4]
    nav_c2=sp[:,5]
    com_c2=sp[:,6]
    # calculate mean across trials for isolation only 
    nav_iso=np.mean(nav_iso.reshape(N_units,trials),axis=1)
    com_iso=np.mean(com_iso.reshape(N_units,trials),axis=1)
    nav_iso=np.repeat(nav_iso,trials)
    com_iso=np.repeat(com_iso,trials)
    # calculate context with spikes counts per trial 
    ce_nav_exp=(nav_c1-nav_iso)/(nav_c1+nav_iso)
    ce_nav_unexp=(com_c1-com_iso)/(com_c1+com_iso)
    ce_com_exp=(com_c2-com_iso)/(com_c2+com_iso)
    ce_com_unexp=(nav_c2-nav_iso)/(nav_c2+nav_iso)
    
    output=np.column_stack((ce_nav_exp, ce_nav_unexp, ce_com_exp,ce_com_unexp))
    
    return output    

def plot_sp_vs_supp(sims):
    
    parameters=sims[0]
    sims = np.delete(sims,[0], 0)
    ntot=len(parameters)
    colors = cm.rainbow(np.linspace(0, 1, ntot))
    
    all_simulations=np.empty((0,6), int)
    fig=plt.figure(1,figsize=(6,3.5))
    ax = plt.gca()
    size=2
    A=-100
    Ny=80
    ny=-50
    Cy=100
    cy=-50
    nx=10
    Nx=150
    cx=280 #510
    Cx=545
    for sim in range(ntot): 
        temp=sims[sim]
        all_simulations=np.row_stack((all_simulations,temp))
        ce_nav_exp=temp[:,0]
        ce_nav_unexp=temp[:,1]
        ce_com_exp=temp[:,2]
        ce_com_unexp=temp[:,3]
        rate_nav=temp[:,4]
        rate_com=temp[:,5]
        
        
        ax=plt.subplot(2,2,1)
        ax.scatter(rate_nav,A*ce_nav_exp,size,color=colors[sim])
        ax.set_ylim(ny, Ny)
        ax.set_xlim(nx, Nx)

        
        ax=plt.subplot(2,2,2)
        ax.scatter(rate_nav,A*ce_nav_unexp,size,color=colors[sim])
        ax.set_ylim(ny, Ny)
        ax.set_xlim(nx, Nx)
        ax.set_yticks([])

        ax=plt.subplot(2,2,3)
        ax.scatter(rate_com,A*ce_com_exp,size,color=colors[sim])
        ax.set_ylim(cy, Cy)
        ax.set_xlim(cx, Cx)

        ax=plt.subplot(2,2,4)
        ax.scatter(rate_com,A*ce_com_unexp,size,color=colors[sim])
        ax.set_ylim(cy, Cy)
        ax.set_xlim(cx, Cx)
        ax.set_yticks([])
        
    R1,p1=pearsonr(all_simulations[:,4],A*all_simulations[:,0])
    R2,p2=pearsonr(all_simulations[:,4],A*all_simulations[:,1])
    R3,p3=pearsonr(all_simulations[:,5],A*all_simulations[:,2])
    R4,p4=pearsonr(all_simulations[:,5],A*all_simulations[:,3])
    
    ax=plt.subplot(2,2,1)
    ax.text(100,-40,str(round(p1, 2)))  #9
    print(R1)
    fit_lin=np.poly1d(np.polyfit(all_simulations[:,4],A*all_simulations[:,0], 1))
    ax.plot(all_simulations[:,4],fit_lin(all_simulations[:,4]),'k',linewidth=0.75)
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    ax.yaxis.set_ticks_position('left')
    ax.xaxis.set_ticks_position('bottom')
    ax.spines['bottom'].set_linewidth(1.5)
    ax.spines['left'].set_linewidth(1.5)
    
    ax=plt.subplot(2,2,2)
    ax.text(100,-40,str(round(p2, 2)))
    fit_lin=np.poly1d(np.polyfit(all_simulations[:,4],A*all_simulations[:,1], 1))
    ax.plot(all_simulations[:,4],fit_lin(all_simulations[:,4]),'k',linewidth=0.75)
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    ax.yaxis.set_ticks_position('left')
    ax.xaxis.set_ticks_position('bottom')
    ax.spines['bottom'].set_linewidth(1.5)
    ax.spines['left'].set_linewidth(1.5)
    
    ax=plt.subplot(2,2,3)
    ax.text(530,-40,str(round(p3, 105)))
    fit_lin=np.poly1d(np.polyfit(all_simulations[:,5],A*all_simulations[:,2], 1))
    ax.plot(all_simulations[:,5],fit_lin(all_simulations[:,5]),'k',linewidth=0.75)
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    ax.yaxis.set_ticks_position('left')
    ax.xaxis.set_ticks_position('bottom')
    ax.spines['bottom'].set_linewidth(1.5)
    ax.spines['left'].set_linewidth(1.5)
    
    ax=plt.subplot(2,2,4)
    ax.text(530,-40,str(round(p4, 48)))
    fit_lin=np.poly1d(np.polyfit(all_simulations[:,5],A*all_simulations[:,3], 1))
    ax.plot(all_simulations[:,5],fit_lin(all_simulations[:,5]),'k',linewidth=0.75)
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    ax.yaxis.set_ticks_position('left')
    ax.xaxis.set_ticks_position('bottom')
    ax.spines['bottom'].set_linewidth(1.5)
    ax.spines['left'].set_linewidth(1.5)

    return fig

def plot_delays_di(sims):
    
    parameters1=sims[0]
    sims = np.delete(sims,0, 0)
    ntrials=len(sims[0])
    ntot=len(parameters1)
    N_units=ntrials/20
    
    # calculate cliff delta from simulations per unit 
    sims_iso=[]
    sims_nav=[]
    sims_com=[]

    for s in range(ntot):
        s_i=sims[s]
        units_iso=[]
        units_nav=[]
        units_com=[]
        for u in range(int(N_units)): 
            sp_e=s_i[20*u:20*u+20,1]
            sp_d=s_i[20*u:20*u+20,2]
            c,size = cliffsDelta(sp_e, sp_d)  
            units_iso.append(c)
            sp_e=s_i[20*u:20*u+20,3]
            sp_d=s_i[20*u:20*u+20,4]
            c,size = cliffsDelta(sp_e, sp_d)  
            units_nav.append(c)
            sp_e=s_i[20*u:20*u+20,5]
            sp_d=s_i[20*u:20*u+20,6]
            c,size = cliffsDelta(sp_e, sp_d)  
            units_com.append(c)
        sims_iso.append(units_iso)    
        sims_nav.append(units_nav)   
        sims_com.append(units_com) 
        

    # plotting mean and error
    iso=np.mean(sims_iso,1)
    nav=np.mean(sims_nav,1)
    com=np.mean(sims_com,1)
    yiso=np.std(sims_iso,1)
    ynav=np.std(sims_nav,1)
    ycom=np.std(sims_com,1)
    
    pvalues_n=[]
    pvalues_c=[]
    for i in range(ntot):
        w_n, p_n = ranksums(np.asarray(sims_iso).flatten(),sims_nav[i])
        w_c, p_c = ranksums(np.asarray(sims_iso).flatten(),sims_com[i])
        pvalues_n.append(p_n)
        pvalues_c.append(p_c)
        print(p_c)
    pvalues=np.column_stack((pvalues_n, pvalues_c))   
    np.save('pvalues_di',pvalues)   
        
    fig, ax = plt.subplots(2,1,sharex=True,sharey=True,figsize=(2.3,2.2))

    
    ax[0].plot(parameters1,iso,color='k',lw=0.75)
    ax[0].plot(parameters1,iso,'.',markersize=3,color='k')
    ax[0].fill_between(parameters1, iso-yiso, iso+yiso ,alpha=0.3, facecolor='k')
    
    ax[1].plot(parameters1,nav,color='b',lw=0.75)
    ax[1].plot(parameters1,nav,'.',markersize=3,color='b')
    ax[1].fill_between(parameters1, nav-ynav, nav+ynav ,alpha=0.3, facecolor='b')

    ax[1].plot(parameters1,com,color='r',lw=0.75)
    ax[1].plot(parameters1,com,'.',markersize=3,color='r')
    ax[1].fill_between(parameters1, com-yiso, com+ycom ,alpha=0.3, facecolor='r')
    
    ax[1].axhline(y=0,color='k',lw=0.5,ls='--')
    
    ## add the real data  
    # import experimental data - cliff delta values for 45 units 
    mat_contents = sio.loadmat(r'E:\Users\User\Documents\model-context\cliff_data.mat')
    a=mat_contents['cliff_data']
    
    data_iso = np.mean(a[:,0])
    data_nav = np.mean(a[:,1])
    data_com = np.mean(a[:,2])
    ydata_iso = np.std(a[:,0])
    ydata_nav = np.std(a[:,1])
    ydata_com = np.std(a[:,2])
    
    del mat_contents
    del a
    
    # import experimental data - cliff delta values for 45 units 
    mat_contents = sio.loadmat(r'E:\Users\User\Documents\model-context\cliff_data-416.mat')
    a=mat_contents['cliff_data']
    
    data_iso_416 = np.mean(a[:,0])
    data_nav_416 = np.mean(a[:,1])
    data_com_416 = np.mean(a[:,2])
    ydata_iso_416 = np.std(a[:,0])
    ydata_nav_416 = np.std(a[:,1])
    ydata_com_416 = np.std(a[:,2])
    
    del mat_contents
    del a
    
 
    ax[0].errorbar([60,416],[data_iso, data_iso_416], yerr=[ydata_iso, ydata_iso_416], fmt='o',color='k', mfc='w',ms=4)
    ax[0].set_ylim(-1,0.5)

    ax[1].errorbar([60,416],[data_nav, data_nav_416], yerr=[ydata_nav, ydata_nav_416], fmt='o',color='b', mfc='w',ms=4)   
    ax[1].errorbar([60,416],[data_com, data_com_416], yerr=[ydata_com, ydata_com_416], fmt='o',color='r', mfc='w',ms=4)
    ax[1].set_ylim(-1,0.5)
    ax[1].set_xlim(0,parameters1[-1])
    ax[1].set_xlabel('gaps (ms)')
    
    ax = ax.flatten() 
    for axi in ax:
        axi.spines['top'].set_visible(False)
        axi.spines['right'].set_visible(False)
        axi.spines['bottom'].set_linewidth(1)
        axi.spines['left'].set_linewidth(1)
#        axi.axhline(y=0,color='gray',lw=2,ls='--')

    plt.figtext(0,0.3, 'discriminability index', rotation=90)
    
    fig.tight_layout()    #w_pad=0.6,h_pad=0.6
    
    return fig


def plot_pvalues(sims):
    
    # import experimental data - cliff delta values for 45 units 
    mat_contents = sio.loadmat(r'E:\Users\User\Documents\model-context\cliff_data.mat')
    a=mat_contents['cliff_data']
    
    exp_data_iso = a[:,0]
    exp_data_nav = a[:,1]
    exp_data_com = a[:,2]
    
    del mat_contents
    del a

    parameters1=sims[0]
    parameters2=sims[1]
    sims = np.delete(sims,[0,1], 0)
    ntrials=len(sims[0])
    nx=len(parameters1)
    ny=len(parameters2)
    ntot=nx*ny
    N_units=ntrials/20
    
    # calculate cliff delta from simulations per unit 
    sims_iso=[]
    sims_nav=[]
    sims_com=[]
    stats_iso=[]
    stats_nav=[]
    stats_com=[]
    for s in range(ntot):
        s_i=sims[s]
        units_iso=[]
        units_nav=[]
        units_com=[]
        for u in range(int(N_units)): 
            sp_e=s_i[20*u:20*u+20,1]
            sp_d=s_i[20*u:20*u+20,2]
            c,size = cliffsDelta(sp_e, sp_d)  
            units_iso.append(c)
            sp_e=s_i[20*u:20*u+20,3]
            sp_d=s_i[20*u:20*u+20,4]
            c,size = cliffsDelta(sp_e, sp_d)  
            units_nav.append(c)
            sp_e=s_i[20*u:20*u+20,5]
            sp_d=s_i[20*u:20*u+20,6]
            c,size = cliffsDelta(sp_e, sp_d)  
            units_com.append(c)
        sims_iso.append(units_iso)    
        sims_nav.append(units_nav)   
        sims_com.append(units_com) 
        
        # check for fitting experimental data
        H,p=ranksums(units_iso, exp_data_iso)
        stats_iso.append(p)
    
        H,p=ranksums(units_nav, exp_data_nav)
        stats_nav.append(p)

        H,p=ranksums(units_com, exp_data_com)
        stats_com.append(p)

    
    # convert into 2D matrix stats 
    stats_iso=np.reshape(stats_iso,(nx,ny)).T
    stats_nav=np.reshape(stats_nav,(nx,ny)).T
    stats_com=np.reshape(stats_com,(nx,ny)).T
    
    # plotting
    fig=plt.figure(1,figsize=(10,3))
    ax = plt.gca()
    xx=1        # factor to change the scale of the units x and y 
    dc=2        # number of decimals to plot in ticks
    
    ax=plt.subplot(131)
 
    ax.imshow(stats_iso,cmap='RdBu',origin='lower')

    
    ax.set_xticks(np.arange(nx))
    ax.set_xticklabels([str(round(float(r), dc)) for r in parameters1*xx])
    ax.set_yticks(np.arange(ny))
    ax.set_yticklabels([str(round(float(r), dc)) for r in parameters2*xx])

    ax=plt.subplot(132)

    ax.imshow(stats_nav,cmap='RdBu',origin='lower')

    
    ax.set_xticks(np.arange(nx))
    ax.set_xticklabels([str(round(float(r), dc)) for r in parameters1*xx])
    ax.set_yticks(np.arange(ny))
    ax.set_yticklabels([str(round(float(r), dc)) for r in parameters2*xx])


    ax=plt.subplot(133)

    ax.imshow(stats_com,cmap='RdBu',origin='lower')


    ax.set_xticks(np.arange(nx))
    ax.set_xticklabels([str(round(float(r), dc)) for r in parameters1*xx])
    ax.set_yticks(np.arange(ny))
    ax.set_yticklabels([str(round(float(r), dc)) for r in parameters2*xx])
    
    
    
    return fig 


def ce_2parameters(sims):

    # import experimental data - context effect values for 45 units 
    mat_contents = sio.loadmat(r'E:\Users\User\Documents\model-context\ei_data-an.mat')
    a=mat_contents['ei']
    
    data_exp_nav = a[:,0]
    data_unexp_nav = a[:,1]
    data_exp_com = a[:,2]
    data_unexp_com = a[:,3]
    
    del mat_contents
    del a

    parameters1=sims[0]
    parameters2=sims[1]
    sims = np.delete(sims,[0,1], 0)
    ntrials=len(sims[0])
    nx=len(parameters1)
    ny=len(parameters2)
    ntot=nx*ny
    N_units=ntrials/20
    
    # calculate context effect from simulations per unit 
    sims_navE=[]
    sims_comE=[]
    sims_navU=[]
    sims_comU=[]
    sims_diff_nav=[]
    sims_diff_com=[]
    
    stats_diff_nav=[]
    stats_diff_com=[]
    stats_navE=[]
    stats_navU=[]
    stats_comE=[]
    stats_comU=[]
    for s in range(ntot):
        s_i=sims[s]
        
        units_navE=[]
        units_navU=[]
        units_comE=[]
        units_comU=[]
        units_diff_nav=[]
        units_diff_com=[]
        
        for u in range(int(N_units)): 
            sp_e_iso=np.mean(s_i[20*u:20*u+20,1])
            sp_d_iso=np.mean(s_i[20*u:20*u+20,2])
            sp_e_nav=np.mean(s_i[20*u:20*u+20,3])
            sp_d_nav=np.mean(s_i[20*u:20*u+20,4])
            sp_e_com=np.mean(s_i[20*u:20*u+20,5])
            sp_d_com=np.mean(s_i[20*u:20*u+20,6])
            ce_exp_nav=(sp_e_nav-sp_e_iso)/(sp_e_nav+sp_e_iso)
            ce_unexp_nav=(sp_d_nav-sp_d_iso)/(sp_d_nav+sp_d_iso)
            ce_exp_com=(sp_d_com-sp_d_iso)/(sp_d_com+sp_d_iso)
            ce_unexp_com=(sp_e_com-sp_e_iso)/(sp_e_com+sp_e_iso)
            
            units_navE.append(ce_exp_nav)
            units_navU.append(ce_unexp_nav)
            units_comE.append(ce_exp_com)
            units_comU.append(ce_unexp_com)
            
            units_diff_nav.append(ce_unexp_nav-ce_exp_nav)
            units_diff_com.append(ce_unexp_com-ce_exp_com)
            
        sims_diff_nav.append(units_diff_nav)  
        sims_diff_com.append(units_diff_com)   
        
        sims_navE.append(units_navE)   
        sims_navU.append(units_navU) 
        
        sims_comE.append(units_comE) 
        sims_comU.append(units_comU)   
        
        # check for fitting experimental data
        H,p=ranksums(units_navE, data_exp_nav)
        if p>0.05:   # equal distributions
            stats_navE.append(1)
        else:
            stats_navE.append(0) 
            
        H,p=ranksums(units_navU, data_unexp_nav)
        if p>0.05:   # equal distributions
            stats_navU.append(1)
        else:
            stats_navU.append(0) 
            
        H,p=ranksums(units_comE, data_exp_com)
        if p>0.05:   # equal distributions
            stats_comE.append(1)
        else:
            stats_comE.append(0) 
            
        H,p=ranksums(units_comU, data_unexp_com)
        if p>0.05:   # equal distributions
            stats_comU.append(1)
        else:
            stats_comU.append(0) 
            
        H,p=ranksums(units_diff_nav, data_unexp_nav-data_exp_nav)
        if p>0.05:   # equal distributions
            stats_diff_nav.append(1)
        else:
            stats_diff_nav.append(0) 
            
        H,p=ranksums(units_diff_com, data_unexp_com-data_exp_com)
        if p>0.05:   # equal distributions
            stats_diff_com.append(1)
        else:
            stats_diff_com.append(0)     

    # avearage context effect across units 
    av_navE=np.mean(sims_navE,1)
    av_navU=np.mean(sims_navU,1)
    av_comE=np.mean(sims_comE,1)
    av_comU=np.mean(sims_comU,1)
    av_diff_nav=np.mean(sims_diff_nav,1)
    av_diff_com=np.mean(sims_diff_com,1)
    
    # convert into 2D matrix results and stats 
    av_navE=np.reshape(av_navE,(nx,ny)).T
    av_navU=np.reshape(av_navU,(nx,ny)).T
    av_comE=np.reshape(av_comE,(nx,ny)).T
    av_comU=np.reshape(av_comU,(nx,ny)).T
    av_diff_nav=np.reshape(av_diff_nav,(nx,ny)).T
    av_diff_com=np.reshape(av_diff_com,(nx,ny)).T
    
    stats_navE=np.reshape(stats_navE,(nx,ny)).T
    stats_navU=np.reshape(stats_navU,(nx,ny)).T
    stats_comE=np.reshape(stats_comE,(nx,ny)).T
    stats_comU=np.reshape(stats_comU,(nx,ny)).T
    stats_diff_nav=np.reshape(stats_diff_nav,(nx,ny)).T
    stats_diff_com=np.reshape(stats_diff_com,(nx,ny)).T
    
    # plotting
    fig, ax = plt.subplots(2,3,sharex=True,sharey=True,figsize=(7,4.5))
    li=1      # min and max of colorbar

    min_x=min(parameters1)-((parameters1[1]-parameters1[0])/2)
    max_x=max(parameters1)+((parameters1[1]-parameters1[0])/2)
    min_y=min(parameters2)-((parameters2[1]-parameters2[0])/2)
    max_y=max(parameters2)+((parameters2[1]-parameters2[0])/2)


    x_isoN=np.nonzero(stats_navE)[1]
    y_isoN=np.nonzero(stats_navE)[0]   
    ax[0,0].imshow(av_navE,cmap='RdBu',origin='lower',vmin=-li,vmax=li,extent=[min_x,max_x,min_y,max_y])
    ax[0,0].scatter(parameters1[x_isoN],parameters2[y_isoN], s=5, c='k', marker="*")
    
    x_navN=np.nonzero(stats_navU)[1]
    y_navN=np.nonzero(stats_navU)[0]   
    ax[0,1].imshow(av_navU,cmap='RdBu',origin='lower',vmin=-li,vmax=li,extent=[min_x,max_x,min_y,max_y])
    ax[0,1].scatter(parameters1[x_navN],parameters2[y_navN], s=5, c='k', marker="*")
    
    x_comN=np.nonzero(stats_diff_nav)[1]
    y_comN=np.nonzero(stats_diff_nav)[0]   
    ax[0,2].imshow(av_diff_nav,cmap='RdBu',origin='lower',vmin=-li,vmax=li,extent=[min_x,max_x,min_y,max_y])
    ax[0,2].scatter(parameters1[x_comN],parameters2[y_comN], s=5, c='k', marker="*")

    x_isoC=np.nonzero(stats_comE)[1]
    y_isoC=np.nonzero(stats_comE)[0]   
    ax[1,0].imshow(av_comE,cmap='RdBu',origin='lower',vmin=-li,vmax=li,extent=[min_x,max_x,min_y,max_y])
    ax[1,0].scatter(parameters1[x_isoC],parameters2[y_isoC], s=5, c='k', marker="*")
    
    x_navC=np.nonzero(stats_comU)[1]
    y_navC=np.nonzero(stats_comU)[0]   
    ax[1,1].imshow(av_comU,cmap='RdBu',origin='lower',vmin=-li,vmax=li,extent=[min_x,max_x,min_y,max_y])
    ax[1,1].scatter(parameters1[x_navC],parameters2[y_navC], s=5, c='k', marker="*")

    x_comC=np.nonzero(stats_diff_com)[1]
    y_comC=np.nonzero(stats_diff_com)[0]   
    ax[1,2].imshow(av_diff_com,cmap='RdBu',origin='lower',vmin=-li,vmax=li,extent=[min_x,max_x,min_y,max_y])
    ax[1,2].scatter(parameters1[x_comC],parameters2[y_comC], s=5, c='k', marker="*")
    
    
    # set which combination satisfy all the contexts
    xy_isoN=np.array(list(zip(x_isoN,y_isoN)))
    xy_navN=np.array(list(zip(x_navN,y_navN)))
    xy_comN=np.array(list(zip(x_comN,y_comN)))
    xy_isoC=np.array(list(zip(x_isoC,y_isoC)))
    xy_navC=np.array(list(zip(x_navC,y_navC)))
    xy_comC=np.array(list(zip(x_comC,y_comC)))
    
    combiN = np.array([x for x in set(tuple(x) for x in xy_navN) & set(tuple(x) for x in xy_comN)])
    combiN = np.array([x for x in set(tuple(x) for x in combiN) & set(tuple(x) for x in xy_isoN)])
    combiC = np.array([x for x in set(tuple(x) for x in xy_navC) & set(tuple(x) for x in xy_comC)])
    combiC = np.array([x for x in set(tuple(x) for x in combiC) & set(tuple(x) for x in xy_isoC)])
    combi = np.array([x for x in set(tuple(x) for x in combiC) & set(tuple(x) for x in combiN)])
 
    
    ax[0,0].scatter(parameters1[combi[:,0]],parameters2[combi[:,1]], s=50, facecolors='none', edgecolors='k')
    ax[0,1].scatter(parameters1[combi[:,0]],parameters2[combi[:,1]], s=50, facecolors='none', edgecolors='k')
    ax[0,2].scatter(parameters1[combi[:,0]],parameters2[combi[:,1]], s=50, facecolors='none', edgecolors='k')
    ax[1,0].scatter(parameters1[combi[:,0]],parameters2[combi[:,1]], s=50, facecolors='none', edgecolors='k')
    ax[1,1].scatter(parameters1[combi[:,0]],parameters2[combi[:,1]], s=50, facecolors='none', edgecolors='k')
    ax[1,2].scatter(parameters1[combi[:,0]],parameters2[combi[:,1]], s=50, facecolors='none', edgecolors='k')
    
    return fig

def plot_delays_ce(sims):
    
    parameters1=sims[0]
    sims = np.delete(sims,0, 0)
    ntrials=len(sims[0])
    ntot=len(parameters1)
    N_units=ntrials/20
    
    # calculate cliff delta from simulations per unit 
    sims_navE=[]
    sims_navU=[]
    sims_comE=[]
    sims_comU=[]

    for s in range(ntot):
        s_i=sims[s]
        units_navE=[]
        units_navU=[]
        units_comE=[]
        units_comU=[]
        for u in range(int(N_units)): 
            sp_e_iso=np.mean(s_i[20*u:20*u+20,1])
            sp_d_iso=np.mean(s_i[20*u:20*u+20,2])
            sp_e_nav=np.mean(s_i[20*u:20*u+20,3])
            sp_d_nav=np.mean(s_i[20*u:20*u+20,4])
            sp_e_com=np.mean(s_i[20*u:20*u+20,5])
            sp_d_com=np.mean(s_i[20*u:20*u+20,6])
            ce_exp_nav=(sp_e_nav-sp_e_iso)/(sp_e_nav+sp_e_iso)
            ce_unexp_nav=(sp_d_nav-sp_d_iso)/(sp_d_nav+sp_d_iso)
            ce_exp_com=(sp_d_com-sp_d_iso)/(sp_d_com+sp_d_iso)
            ce_unexp_com=(sp_e_com-sp_e_iso)/(sp_e_com+sp_e_iso)
            
            units_navE.append(ce_exp_nav)
            units_navU.append(ce_unexp_nav)
            units_comE.append(ce_exp_com)
            units_comU.append(ce_unexp_com)
            
            
        sims_navE.append(units_navE)    
        sims_navU.append(units_navU)   
        sims_comE.append(units_comE) 
        sims_comU.append(units_comU) 

    # plotting mean and error
    navE=np.mean(sims_navE,1)
    navU=np.mean(sims_navU,1)
    comE=np.mean(sims_comE,1)
    comU=np.mean(sims_comU,1)
    ynavE=np.std(sims_navE,1)
    ynavU=np.std(sims_navU,1)
    ycomE=np.std(sims_comE,1)
    ycomU=np.std(sims_comU,1)
    
    pvalues_nexp=[]
    pvalues_nunexp=[]
    pvalues_cexp=[]
    pvalues_cunexp=[]
    for i in range(ntot):
        w_ne, p_ne = wilcoxon(sims_navE[i])
        w_ne, p_nu = wilcoxon(sims_navU[i])
        w_ne, p_ce = wilcoxon(sims_comE[i])
        w_ne, p_cu = wilcoxon(sims_comU[i])
        pvalues_nexp.append(p_ne)
        pvalues_nunexp.append(p_nu)
        pvalues_cexp.append(p_ce)
        pvalues_cunexp.append(p_cu)
        print(p_ne)
    pvalues=np.column_stack((pvalues_nexp,pvalues_nunexp,pvalues_cexp,pvalues_cunexp))   
    np.save('pvalues_ce',pvalues) 

    fig, ax = plt.subplots(2,2,sharex=True,sharey=True,figsize=(2.6,1.9))


    ax[0,0].plot(parameters1,navE,color='b',lw=0.75)
    ax[0,0].plot(parameters1,navE,'.',markersize=3,color='b')
    ax[0,0].fill_between(parameters1, navE-ynavE, navE+ynavE ,alpha=0.3, facecolor='b')

    ax[1,0].plot(parameters1,navU,color='b',lw=0.75) 
    ax[1,0].plot(parameters1,navU,'.',markersize=3,color='b')     
    ax[1,0].fill_between(parameters1, navU-ynavU, navU+ynavU ,alpha=0.3, facecolor='b')
    
    ax[0,1].plot(parameters1,comE,color='r',lw=0.75)
    ax[0,1].plot(parameters1,comE,'.',markersize=3,color='r')
    ax[0,1].fill_between(parameters1, comE-ycomE, comE+ycomE ,alpha=0.3, facecolor='r')

    ax[1,1].plot(parameters1,comU,color='r',lw=0.75)
    ax[1,1].plot(parameters1,comU,'.',markersize=3,color='r')
    ax[1,1].fill_between(parameters1, comU-ycomU, comU+ycomU ,alpha=0.3, facecolor='r')
 
    ## add the real data  
    # import experimental data - context effect values for 45 units after 60 ms
    mat_contents = sio.loadmat(r'E:\Users\User\Documents\model-context\ei_data.mat')
    a=mat_contents['ei']
    
    data_exp_nav = np.mean(a[:,0])
    data_unexp_nav = np.mean(a[:,1])
    data_exp_com = np.mean(a[:,2])
    data_unexp_com = np.mean(a[:,3])

    ydata_exp_nav = np.std(a[:,0])
    ydata_unexp_nav = np.std(a[:,1])
    ydata_exp_com = np.std(a[:,2])
    ydata_unexp_com = np.std(a[:,3])
    
    del mat_contents
    del a
    
    # import experimental data - context effect values for 45 units after 416
    mat_contents = sio.loadmat(r'E:\Users\User\Documents\model-context\ei_data-416.mat')
    a=mat_contents['ei']
    
    data_exp_nav_416 = np.mean(a[:,0])
    data_unexp_nav_416 = np.mean(a[:,1])
    data_exp_com_416 = np.mean(a[:,2])
    data_unexp_com_416 = np.mean(a[:,3])

    ydata_exp_nav_416 = np.std(a[:,0])
    ydata_unexp_nav_416 = np.std(a[:,1])
    ydata_exp_com_416 = np.std(a[:,2])
    ydata_unexp_com_416 = np.std(a[:,3])
    
    del mat_contents
    del a
    

    ax[0,0].errorbar([60,416],[data_exp_nav, data_exp_nav_416], yerr=[ydata_exp_nav, ydata_exp_nav_416], fmt='o',color='b', mfc='w',label='nav exp',ms=4,lw=0.75)
    ax[1,0].errorbar([60,416],[data_unexp_nav, data_unexp_nav_416], yerr=[ydata_unexp_nav, ydata_unexp_nav_416], fmt='o',color='b', mfc='w',label='nav unexp',ms=4,lw=0.75)
    ax[0,1].errorbar([60,416],[data_exp_com, data_exp_com_416], yerr=[ydata_exp_com, ydata_exp_com_416], fmt='o',color='r', mfc='w',label='com exp',ms=4,lw=0.75)
    ax[1,1].errorbar([60,416],[data_unexp_com, data_unexp_com_416], yerr=[ydata_unexp_com, ydata_unexp_com_416], fmt='o',color='r', mfc='w',label='com unexp',ms=4,lw=0.75)
    
#    ax[0].legend(loc='lower right',frameon=False)
#    ax[1].legend(loc='lower right',frameon=False)
#    ax[2].legend(loc='lower right',frameon=False)
#    ax[3].legend(loc='lower right',frameon=False)
    
    ax = ax.flatten() 
    for axi in ax:
        axi.spines['top'].set_visible(False)
        axi.spines['right'].set_visible(False)
        axi.spines['bottom'].set_linewidth(0.75)
        axi.spines['left'].set_linewidth(0.75)
        axi.axhline(y=0,color='k',lw=0.5,ls='--')

    plt.figtext(0.525,0.3, 'context effect (com)', rotation=90)
    plt.figtext(0,0.3, 'context effect (nav)', rotation=90)
    plt.figtext(0.5,0, 'gaps(ms)', rotation=0, verticalalignment='bottom')
    fig.tight_layout(h_pad=2)  #w_pad=0.6,h_pad=0.6    

    return fig

def variables_delay(N_units):
    delay=3000
    trials=20
    dt_=0.1
 
    # SIMULATION NAVIGATION CONTEXT
    context=1
    exp=1
    exec(open("./input.py").read())  
    exec(open("./model.py").read()) 
    nav_postsyn=pos.V_th*1000   # to mV 
    nav_presyn=mon.x_S
    nav_high=nav_presyn[trials*N_units:]
    nav_low=nav_presyn[:trials*N_units]
    nav_tot=np.shape(nav_presyn)[1]*dt_
    nav_time=np.arange(0,nav_tot,dt_)
    # average acros trials
    nav_post=np.mean(nav_postsyn,0)
    nav_pre_high=np.mean(nav_high,0)
    nav_pre_low=np.mean(nav_low,0)
    
    # SIMULATION COMMUNICATION CONTEXT
    context=2
    exec(open("./input.py").read())  
    exec(open("./model.py").read()) 
    com_postsyn=pos.V_th*1000
    com_presyn=mon.x_S
    com_high=com_presyn[trials*N_units:]
    com_low=com_presyn[:trials*N_units]
    com_tot=np.shape(com_presyn)[1]*dt_
    com_time=np.arange(0,com_tot,dt_)
    # average across trials
    com_post=np.mean(com_postsyn,0)
    com_pre_high=np.mean(com_high,0)
    com_pre_low=np.mean(com_low,0)

    # plotting
    vm_i=-50
    vm_s=-44
    xs_i=0.3
    xs_s=1

    fig, ax = plt.subplots(2,2,sharex=True,sharey='row',figsize=(3.54,2))
    
    # navigation postsyn
    width_context=(onset[1]-delay)-onset[0]
    rect_context=mpl.patches.Rectangle((onset[0],vm_i),width_context,(vm_s-vm_i),facecolor='b',alpha=0.2)
    ax[0,0].add_patch(rect_context)
    ax[0,0].plot(nav_time,nav_post,c='k',lw=0.75)
    ax[0,0].vlines(onset[1]-delay+60,vm_i,vm_s,lw=0.75)
    ax[0,0].vlines(onset[1]-delay+416,vm_i,vm_s,lw=0.75)
    ax[0,0].set_ylim(vm_i,vm_s)
    ax[0,0].set_xlim(500,6000)
    ax[0,0].set_ylabel('spiking\n threshold')
    
    # communication postsyn
    width_context=(onset[2]-delay)-onset[0]
    rect_context=mpl.patches.Rectangle((onset[0],vm_i),width_context,(vm_s-vm_i),facecolor='r',alpha=0.2)
    ax[0,1].add_patch(rect_context)
    ax[0,1].plot(com_time,com_post,c='k',lw=0.75)
    ax[0,1].vlines(onset[2]-delay+60,vm_i,vm_s,lw=0.75)
    ax[0,1].vlines(onset[2]-delay+416,vm_i,vm_s,lw=0.75)
    ax[0,1].set_ylim(vm_i,vm_s)
    ax[0,1].set_xlim(500,6000)
    
    # navigation presyn
    width_context=(onset[1]-delay)-onset[0]
    rect_context=mpl.patches.Rectangle((onset[0],xs_i),width_context,(xs_s-xs_i),facecolor='b',alpha=0.2)
    ax[1,0].add_patch(rect_context)
    ax[1,0].plot(nav_time,nav_pre_high,c='b',lw=0.75)
    ax[1,0].plot(nav_time,nav_pre_low,c='r',lw=0.75)
    ax[1,0].vlines(onset[1]-delay+60,xs_i,xs_s,lw=0.75)
    ax[1,0].vlines(onset[1]-delay+416,xs_i,xs_s,lw=0.75)
    ax[1,0].set_ylim(xs_i,xs_s)
    ax[1,0].set_xlim(500,6000)
    ax[1,0].set_ylabel('strength\n synapse')
    
    # communication presyn
    width_context=(onset[2]-delay)-onset[0]
    rect_context=mpl.patches.Rectangle((onset[0],xs_i),width_context,(xs_s-xs_i),facecolor='r',alpha=0.2)
    ax[1,1].add_patch(rect_context)
    ax[1,1].plot(com_time,com_pre_high,c='b',lw=0.75)
    ax[1,1].plot(com_time,com_pre_low,c='r',lw=0.75)
    ax[1,1].vlines(onset[2]-delay+60,xs_i,xs_s,lw=0.75)
    ax[1,1].vlines(onset[2]-delay+416,xs_i,xs_s,lw=0.75)
    ax[1,1].set_ylim(xs_i,xs_s)
    ax[1,1].set_xlim(500,6000)
    
    
    ax = ax.flatten() 
    for axi in ax:
        axi.spines['top'].set_visible(False)
        axi.spines['right'].set_visible(False)
        axi.spines['bottom'].set_linewidth(0.75)
        axi.spines['left'].set_linewidth(0.75)
   
    plt.figtext(0.5, 0.01, 'Time (ms)', rotation=0, verticalalignment='bottom')
    
    fig.tight_layout() 
    
    return fig 
    
def spikes_input():

    N_units=1
    context=1
    exp=1
    exec(open("./input.py").read())  
    exec(open("./model.py").read())
    HF=firing_context(sp_high,onset,context)
    LF=firing_context(sp_low,onset,context)
    
    output=[x + y for (x, y) in zip(LF, HF)]
    av_spikes_c1=sum(output) / len(output)
    
    context=2
    exec(open("./input.py").read())  
    exec(open("./model.py").read())
    HF=firing_context(sp_high,onset,context)
    LF=firing_context(sp_low,onset,context)
    
    output=[x + y for (x, y) in zip(LF, HF)]
    av_spikes_c2=sum(output) / len(output)
    
    return av_spikes_c1,av_spikes_c2

def sim_isolation(N_units,we_LF=6,we_HF=6):
    
    probe=np.repeat(range(1,N_units+1),20)
    # SIMULATION 1: navigation - after silence
    context=0
    exp=1
    exec(open("./input.py").read())  
    exec(open("./model.py").read()) 
    col=spike_counts(ac,onset,context)
    probe=np.column_stack((probe, col))
    
    # SIMULATION 2: communication - after silence
    exp=2
    exec(open("./input.py").read())  
    exec(open("./model.py").read()) 
    col=spike_counts(ac,onset,context)
    probe=np.column_stack((probe, col))
    
    return probe
    
def plot_isolation(sims):
    # import experimental data - cliff delta values for 45 units 
    mat_contents = sio.loadmat(r'E:\Users\User\Documents\model-context\cliff_data-pref.mat')
    a=mat_contents['cliff_data']
    
    exp_data_iso = a[:,0]

    del mat_contents
    del a
    
    # plot cliff delta and spikecounts for probes in silence
    parameters1=sims[0]
    parameters2=sims[1]
    sims = np.delete(sims,[0,1], 0)
    ntrials=len(sims[0])
    nx=len(parameters1)
    ny=len(parameters2)
    ntot=nx*ny
    N_units=ntrials/20
    
    # calculate cliff delta and spikecounts from simulations per unit 
    sims_iso=[]
    sp_nav=[]
    sp_com=[]
    stats_iso=[]
    for s in range(ntot):
        s_i=sims[s]
        units_iso=[]
        units_sp_nav=[]
        units_sp_com=[]
        for u in range(int(N_units)): 
            sp_e=s_i[20*u:20*u+20,1]
            sp_d=s_i[20*u:20*u+20,2]
            c,size = cliffsDelta(sp_e, sp_d)  
            units_iso.append(c)
            units_sp_nav.append(np.mean(sp_e)) # average across trials 
            units_sp_com.append(np.mean(sp_d))
            
        sims_iso.append(units_iso)    
        sp_nav.append(units_sp_nav)   
        sp_com.append(units_sp_com) 
        
        # check for fitting experimental data
        H,p=ranksums(units_iso, exp_data_iso)
        if p>0.05:   # equal distributions
            stats_iso.append(1)
        else:
            stats_iso.append(0) 
             
            
            
    # avearage cliff delta and spikecounts across units 
    av_iso=np.mean(sims_iso,1)
    av_nav=np.mean(sp_nav,1)
    av_com=np.mean(sp_com,1)
    
    # convert into 2D matrix results and stats 
    av_iso=np.reshape(av_iso,(nx,ny)).T
    av_nav=np.reshape(av_nav,(nx,ny)).T
    av_com=np.reshape(av_com,(nx,ny)).T
    stats_iso=np.reshape(stats_iso,(nx,ny)).T
    
    # plotting
    fig, ax = plt.subplots(1,3,sharex=True,sharey=True,figsize=(6.5,2.5))
    li=1      # min and max of colorbar
    xx=1        # factor to change the scale of the units x and y 
    dc=0      # number of decimals to plot in ticks
    mi=10     # max spikecounts to plot in colorbar 

    min_x=min(parameters1)-((parameters1[1]-parameters1[0])/2)
    max_x=max(parameters1)+((parameters1[1]-parameters1[0])/2)
    min_y=min(parameters2)-((parameters2[1]-parameters2[0])/2)
    max_y=max(parameters2)+((parameters2[1]-parameters2[0])/2)
    

    ax[0].imshow(av_iso,cmap='RdBu',origin='lower',vmin=-li,vmax=li,extent=[min_x,max_x,min_y,max_y])  
    ax[1].imshow(av_nav,cmap='rainbow',origin='lower',vmin=0,vmax=mi,extent=[min_x,max_x,min_y,max_y])
    ax[2].imshow(av_com,cmap='rainbow',origin='lower',vmin=0,vmax=mi,extent=[min_x,max_x,min_y,max_y])

    x_iso=np.nonzero(stats_iso)[1]
    y_iso=np.nonzero(stats_iso)[0]    
    ax[0].scatter(parameters1[x_iso],parameters2[y_iso], s=5, c='k', marker="*")
 

    
    return fig

def discriminability_index(probe):
    # calculate context effect (c.e.) for each combination per unit 
    trials=20
    N_units=int(len(probe)/trials)
    
    output=[]
    for u in range(int(N_units)): 
        cliff_iso,size = cliffsDelta(probe[20*u:20*u+20,1], probe[20*u:20*u+20,2]) 
        cliff_nav,size = cliffsDelta(probe[20*u:20*u+20,3], probe[20*u:20*u+20,4]) 
        cliff_com,size = cliffsDelta(probe[20*u:20*u+20,5], probe[20*u:20*u+20,6]) 
        output.append([cliff_nav, cliff_iso, cliff_com])

    output=np.asarray(output)        
    return output

def plot_ce_di_multiple(sims1, sims2, sims3, sims4=None):
    
    # c.e.
    ce1=context_effect(sims1)
    ce2=context_effect(sims2)
    ce3=context_effect(sims3)

    # d.i.
    di1=discriminability_index(sims1)
    di2=discriminability_index(sims2)
    di3=discriminability_index(sims3)
  

    fig, ax = plt.subplots(1,1,figsize=(2.85,1.5))  #2.5
    
    aa=0.5
    colors=[nc.to_rgba('grey', alpha=aa),nc.to_rgba('green', alpha=aa),nc.to_rgba('orange', alpha=aa)]
    
    ax.axhline(y=0,color='k',lw=0.5,ls='--')
    b1=ax.boxplot(di1,positions=[1,2,3], patch_artist=True,widths=0.8)
    for item in ['whiskers', 'fliers', 'medians', 'caps','boxes']:
        for c in b1[item]:
            c.set(color='k', linewidth=0.5)
    for i in b1['boxes']:
        i.set_facecolor(colors[0])
    for i in b1['fliers']:
        i.set_markersize(2)
        i.set_markeredgewidth(0.5)
        
    b2=ax.boxplot(di2,positions=[5,6,7], patch_artist=True,widths=0.8)
    for item in ['whiskers', 'fliers', 'medians', 'caps','boxes']:
        for c in b2[item]:
            c.set(color='k', linewidth=0.5)
    for i in b2['boxes']:
        i.set_facecolor(colors[1])  
    for i in b2['fliers']:
        i.set_markersize(2)
        i.set_markeredgewidth(0.5)
        
    b3=ax.boxplot(di3,positions=[9,10,11], patch_artist=True,widths=0.8)
    for item in ['whiskers', 'fliers', 'medians', 'caps','boxes']:
        for c in b3[item]:
            c.set(color='k', linewidth=0.5)
    for i in b3['boxes']:
        i.set_facecolor(colors[2])
    for i in b3['fliers']:
        i.set_markersize(2)   
        i.set_markeredgewidth(0.5)
        
    ax.set_ylim(-1,1)
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.yaxis.set_ticks_position('left')
    ax.xaxis.set_ticks_position('bottom')

    for axis in ['top','bottom','left','right']:
      ax.spines[axis].set_linewidth(0.75)
    plt.locator_params(axis='y',nbins=5)
  
        
    ax.tick_params(axis='x', rotation=45)
    
    ax.set_ylabel('discriminability index')
    
    if sims4 is not None:
        ce4=context_effect(sims4)
        di4=discriminability_index(sims4)
        b4=ax.boxplot(di4,positions=[13,14,15], patch_artist=True,widths=0.8)
        for item in ['whiskers', 'fliers', 'medians', 'caps','boxes']:
            for c in b4[item]:
                c.set(color='k', linewidth=0.5)
        for i in b4['boxes']:
            i.set_facecolor(nc.to_rgba('blue', alpha=aa))  
        for i in b4['fliers']:
            i.set_markersize(2)
            i.set_markeredgewidth(0.5)
        ax.set_xticklabels(['ech','silence','com','ech','silence','com','ech','silence','com','ech','silence','com'],ha='right')
    else:
        ax.set_xticklabels(['ech','silence','com','ech','silence','com','ech','silence','com'],ha='right')
    
    
    
    fig.tight_layout()
    
    
    # context effect plots
    fig2, ax = plt.subplots(2,1,sharex=True,sharey=True,figsize=(2.5,3))
    
    aa=0.5
    colors=[nc.to_rgba('grey', alpha=aa),nc.to_rgba('green', alpha=aa),nc.to_rgba('orange', alpha=aa)]
    
    ax[0].axhline(y=0,color='k',lw=0.5,ls='--')
    ax[1].axhline(y=0,color='k',lw=0.5,ls='--')
    
    nav1=ce1[:,0:2]
    com1=ce1[:,2:4]
    b1=ax[0].boxplot(nav1,positions=[1,2], patch_artist=True,widths=0.8)
    for item in ['whiskers', 'fliers', 'medians', 'caps','boxes']:
        for c in b1[item]:
            c.set(color='k', linewidth=0.5)
    for i in b1['boxes']:
        i.set_facecolor(colors[0])
    for i in b1['fliers']:
        i.set_markersize(2)
        i.set_markeredgewidth(0.5)        
        
    b1=ax[1].boxplot(com1,positions=[1,2], patch_artist=True,widths=0.8)
    for item in ['whiskers', 'fliers', 'medians', 'caps','boxes']:
        for c in b1[item]:
            c.set(color='k', linewidth=0.5)
    for i in b1['boxes']:
        i.set_facecolor(colors[0])
    for i in b1['fliers']:
        i.set_markersize(2)
        i.set_markeredgewidth(0.5)
        
    nav2=ce2[:,0:2]
    com2=ce2[:,2:4]    
    b2=ax[0].boxplot(nav2,positions=[4,5], patch_artist=True,widths=0.8)
    for item in ['whiskers', 'fliers', 'medians', 'caps','boxes']:
        for c in b2[item]:
            c.set(color='k', linewidth=0.5)
    for i in b2['boxes']:
        i.set_facecolor(colors[1])
    for i in b2['fliers']:
        i.set_markersize(2)
        i.set_markeredgewidth(0.5)
        
    b2=ax[1].boxplot(com2,positions=[4,5], patch_artist=True,widths=0.8)
    for item in ['whiskers', 'fliers', 'medians', 'caps','boxes']:
        for c in b2[item]:
            c.set(color='k', linewidth=0.5)
    for i in b2['boxes']:
        i.set_facecolor(colors[1])  
    for i in b2['fliers']:
        i.set_markersize(2)
        i.set_markeredgewidth(0.5)
        
    nav3=ce3[:,0:2]
    com3=ce3[:,2:4]    
    b3=ax[0].boxplot(nav3,positions=[7,8], patch_artist=True,widths=0.8)
    for item in ['whiskers', 'fliers', 'medians', 'caps','boxes']:
        for c in b3[item]:
            c.set(color='k', linewidth=0.5)
    for i in b3['boxes']:
        i.set_facecolor(colors[2])
    for i in b3['fliers']:
        i.set_markersize(2)
        i.set_markeredgewidth(0.5)
        
    b3=ax[1].boxplot(com3,positions=[7,8], patch_artist=True,widths=0.8)
    for item in ['whiskers', 'fliers', 'medians', 'caps','boxes']:
        for c in b3[item]:
            c.set(color='k', linewidth=0.5)
    for i in b3['boxes']:
        i.set_facecolor(colors[2])
    for i in b3['fliers']:
        i.set_markersize(2)
        i.set_markeredgewidth(0.5)
        
    if sims4 is not None:
        nav4=ce4[:,0:2]
        com4=ce4[:,2:4]    
        b4=ax[0].boxplot(nav4,positions=[10,11], patch_artist=True,widths=0.8)
        for item in ['whiskers', 'fliers', 'medians', 'caps','boxes']:
            for c in b4[item]:
                c.set(color='k', linewidth=0.5)
        for i in b4['boxes']:
            i.set_facecolor(nc.to_rgba('blue', alpha=aa))
        for i in b4['fliers']:
            i.set_markersize(2)
            i.set_markeredgewidth(0.5)
            
        b4=ax[1].boxplot(com4,positions=[10,11], patch_artist=True,widths=0.8)
        for item in ['whiskers', 'fliers', 'medians', 'caps','boxes']:
            for c in b4[item]:
                c.set(color='k', linewidth=0.5)
        for i in b4['boxes']:
            i.set_facecolor(nc.to_rgba('blue', alpha=aa))
        for i in b4['fliers']:
            i.set_markersize(2)
            i.set_markeredgewidth(0.5)

        ax[1].set_xticks([1,2,4,5,7,8,10,11])
        ax[1].set_xticklabels(['match','mismatch','match','mismatch','match','mismatch','match','mismatch'],ha='right')
    else:
        ax[1].set_xticks([1,2,4,5,7,8])
        ax[1].set_xticklabels(['match','mismatch','match','mismatch','match','mismatch'],ha='right')
        
        
    for axi in ax: 
        axi.set_ylim(-1,0.5)
        axi.spines['top'].set_visible(False)
        axi.spines['right'].set_visible(False)
        axi.yaxis.set_ticks_position('left')
        axi.xaxis.set_ticks_position('bottom')
    
        for axis in ['top','bottom','left','right']:
          axi.spines[axis].set_linewidth(0.75)
        plt.locator_params(axis='y',nbins=5)
   
        axi.tick_params(axis='x', rotation=45)
        axi.set_ylabel('context effect')
        
        
    fig2.tight_layout()
    
    
    return fig, fig2


def plot_experimental():
    
    # import experimental data - context effect values for 45 units after 60 ms
    mat_contents = sio.loadmat(r'E:\Users\User\Documents\model-context\ei_data.mat')
    a=mat_contents['ei']

    ce_exp_nav = a[:,0]
    ce_unexp_nav = a[:,1]
    ce_exp_com = a[:,2]
    ce_unexp_com = a[:,3]
    ce_nav=[ce_exp_nav, ce_unexp_nav]
    ce_com=[ce_exp_com, ce_unexp_com]
    
    del mat_contents
    del a
    
    # import experimental data - cliff delta values for 45 units 
    mat_contents = sio.loadmat(r'E:\Users\User\Documents\model-context\cliff_data-pref.mat')
    a=mat_contents['cliff_data']
    
    di_iso = a[:,0]
    di_nav = a[:,1]
    di_com = a[:,2]
    di=[di_nav, di_iso, di_com]
    
    del mat_contents
    del a

    # plot discriminability index
    
    fig, ax = plt.subplots(1,1,figsize=(1.1,1.5))
    
    aa=0.5
    colors=[nc.to_rgba('white', alpha=aa),nc.to_rgba('green', alpha=aa),nc.to_rgba('orange', alpha=aa)]
    
    ax.axhline(y=0,color='k',lw=0.5,ls='--')
    b1=ax.boxplot(di,positions=[1,2,3], patch_artist=True,widths=0.7)
    for item in ['whiskers', 'fliers', 'medians', 'caps','boxes']:
        for c in b1[item]:
            c.set(color='k', linewidth=0.5)
    for i in b1['boxes']:
        i.set_facecolor(colors[0])
    for i in b1['fliers']:
        i.set_markersize(2)
        i.set_markeredgewidth(0.5)
        
        
    ax.set_ylim(-1,1)
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.yaxis.set_ticks_position('left')
    ax.xaxis.set_ticks_position('bottom')

    for axis in ['top','bottom','left','right']:
      ax.spines[axis].set_linewidth(0.75)
    plt.locator_params(axis='y',nbins=5)
  
    ax.set_xticklabels(['ech','silence','com'],ha='right')
    ax.tick_params(axis='x', rotation=45)
    
    ax.set_ylabel('discriminability index')
    
    fig.tight_layout()
    
    # plot context effect 
    fig2, ax = plt.subplots(1,1,figsize=(1.5,1.5))
    
    aa=0.5
    colors=[nc.to_rgba('white', alpha=aa),nc.to_rgba('green', alpha=aa),nc.to_rgba('orange', alpha=aa)]
    
    ax.axhline(y=0,color='k',lw=0.5,ls='--')
    
    b1=ax.boxplot(ce_nav,positions=[1,2], patch_artist=True,widths=0.5)
    for item in ['whiskers', 'fliers', 'medians', 'caps','boxes']:
        for c in b1[item]:
            c.set(color='k', linewidth=0.5)
    for i in b1['boxes']:
        i.set_facecolor(colors[0])
    for i in b1['fliers']:
        i.set_markersize(2)
        i.set_markeredgewidth(0.5)
        
    b1=ax.boxplot(ce_com,positions=[3,4], patch_artist=True,widths=0.5)
    for item in ['whiskers', 'fliers', 'medians', 'caps','boxes']:
        for c in b1[item]:
            c.set(color='k', linewidth=0.5)
    for i in b1['boxes']:
        i.set_facecolor(colors[0])
    for i in b1['fliers']:
        i.set_markersize(2)
        i.set_markeredgewidth(0.5)

    ax.set_xticks([1,2,3,4])
    ax.set_xticklabels(['match','mismatch','match','mismatch'],ha='right')
        
        
  
    ax.set_ylim(-1,0.5)
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.yaxis.set_ticks_position('left')
    ax.xaxis.set_ticks_position('bottom')
    
    for axis in ['top','bottom','left','right']:
        ax.spines[axis].set_linewidth(0.75)
    plt.locator_params(axis='y',nbins=5)
   
    ax.tick_params(axis='x', rotation=45)
    ax.set_ylabel('context effect')
        
        
    fig2.tight_layout()
    
    return fig, fig2


def plot_pref_com(sims):
    
    # discriminability index from simulations
    di=discriminability_index(sims)
    
    fig, ax = plt.subplots(1,1,figsize=(1.1,1.5))
    
    aa=0.5
    colors=[nc.to_rgba('blue', alpha=aa),nc.to_rgba('green', alpha=aa),nc.to_rgba('white', alpha=aa)]
    
    ax.axhline(y=0,color='k',lw=0.5,ls='--')
    b1=ax.boxplot(di,positions=[1,2,3], patch_artist=True,widths=0.7)
    for item in ['whiskers', 'fliers', 'medians', 'caps','boxes']:
        for c in b1[item]:
            c.set(color='k', linewidth=0.5)
    for i in b1['boxes']:
        i.set_facecolor(colors[0])
    for i in b1['fliers']:
        i.set_markersize(2)
        i.set_markeredgewidth(0.5)
        
        
    ax.set_ylim(-1,1)
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.yaxis.set_ticks_position('left')
    ax.xaxis.set_ticks_position('bottom')

    for axis in ['top','bottom','left','right']:
      ax.spines[axis].set_linewidth(0.75)
    plt.locator_params(axis='y',nbins=5)
  
    ax.set_xticklabels(['ech','silence','com'],ha='right')
    ax.tick_params(axis='x', rotation=45)
    
    ax.set_ylabel('discriminability index')
    
    fig.tight_layout()
    
    return fig 

def plot_pvalues_di(sims):
    delays=np.arange(50,2200,150)
    
    fig, ax = plt.subplots(1,1,figsize=(2.3,1))
    
    ax.hlines(-np.log(0.05),0,delays[-1],linewidth=0.5)
    ax.vlines(delays[4],-15,85,linewidth=0.5,linestyle=(0,(5,5)))
    ax.vlines(delays[10],-15,85,linewidth=0.5,linestyle=(0,(5,5)))
    ax.plot(delays,-np.log(sims[:,0]),'b.',markersize=3)
    ax.plot(delays,-np.log(sims[:,1]),'r.',markersize=3)


    ax.set_ylabel('-log(p-value)')
    ax.set_xticklabels([])
#    ax.set_xlabel('gaps (ms)')
    
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.yaxis.set_ticks_position('left')
    ax.xaxis.set_ticks_position('bottom')

    for axis in ['top','bottom','left','right']:
        ax.spines[axis].set_linewidth(0.75)
        
    plt.locator_params(axis='y',nbins=5)
    ax.set_xlim(0,delays[-1])
    ax.set_ylim(-15,85)
    
    fig.tight_layout()
    
    return fig

def plot_pvalues_ce(sims):
    delays=np.arange(50,2200,150)
    
    fig, ax = plt.subplots(1,2,sharex=True,sharey=True,figsize=(2.8,1))
    

    ax[0].axhline(y=-np.log(0.05),color='k',lw=0.5,ls='-')
    ax[0].vlines(delays[10],-15,60,linewidth=0.5,linestyle=(0,(5,5)))
    ax[0].plot(delays,-np.log(sims[:,0]),'b.',markersize=3)
    ax[0].plot(delays,-np.log(sims[:,1]),'r.',markersize=3)


    ax[0].set_ylabel('-log(p-value)')
    ax[0].set_xticklabels([])
    
    ax[0].spines['top'].set_visible(False)
    ax[0].spines['right'].set_visible(False)
    ax[0].yaxis.set_ticks_position('left')
    ax[0].xaxis.set_ticks_position('bottom')

    for axis in ['top','bottom','left','right']:
        ax[0].spines[axis].set_linewidth(0.75)
        
    plt.locator_params(axis='y',nbins=5)

    ax[0].set_ylim(-5,50)

    ax[1].axhline(y=-np.log(0.05),color='k',lw=0.5,ls='-')
    ax[1].vlines(delays[4],-15,60,linewidth=0.5,linestyle=(0,(5,5)))
    ax[1].plot(delays,-np.log(sims[:,2]),'r.',markersize=3)
    ax[1].plot(delays,-np.log(sims[:,3]),'b.',markersize=3)

    ax[1].spines['top'].set_visible(False)
    ax[1].spines['right'].set_visible(False)
    ax[1].yaxis.set_ticks_position('left')
    ax[1].xaxis.set_ticks_position('bottom')


    for axis in ['top','bottom','left','right']:
      ax[1].spines[axis].set_linewidth(0.75)
    plt.locator_params(axis='y',nbins=5)

    fig.tight_layout()
    
    return fig

def dd_2parameters(sims):

    # import experimental data - context effect values for 45 units awake
    mat_contents = sio.loadmat(r'E:\Users\User\Documents\model-context\ei_data-an.mat')
    a=mat_contents['ei']
    
    data_exp_nav = a[:,0]
    data_unexp_nav = a[:,1]
    data_exp_com = a[:,2]
    data_unexp_com = a[:,3]
    
    dd_nav=(data_unexp_nav-data_exp_nav)/2
    dd_com=(data_unexp_com-data_exp_com)/2
    
    del mat_contents
    del a

    parameters1=sims[0]
    parameters2=sims[1]
    sims = np.delete(sims,[0,1], 0)
    ntrials=len(sims[0])
    nx=len(parameters1)
    ny=len(parameters2)
    ntot=nx*ny
    N_units=ntrials/20
    
    # calculate context effect from simulations per unit 

    sims_dd_nav=[]
    sims_dd_com=[]
    
    stats_dd_nav=[]
    stats_dd_com=[]

    for s in range(ntot):
        s_i=sims[s]
        
        units_dd_nav=[]
        units_dd_com=[]
        
        for u in range(int(N_units)): 
            sp_e_iso=np.mean(s_i[20*u:20*u+20,1])
            sp_d_iso=np.mean(s_i[20*u:20*u+20,2])
            sp_e_nav=np.mean(s_i[20*u:20*u+20,3])
            sp_d_nav=np.mean(s_i[20*u:20*u+20,4])
            sp_e_com=np.mean(s_i[20*u:20*u+20,5])
            sp_d_com=np.mean(s_i[20*u:20*u+20,6])
            ce_exp_nav=(sp_e_nav-sp_e_iso)/(sp_e_nav+sp_e_iso)
            ce_unexp_nav=(sp_d_nav-sp_d_iso)/(sp_d_nav+sp_d_iso)
            ce_exp_com=(sp_d_com-sp_d_iso)/(sp_d_com+sp_d_iso)
            ce_unexp_com=(sp_e_com-sp_e_iso)/(sp_e_com+sp_e_iso)
            
            units_dd_nav.append((ce_unexp_nav-ce_exp_nav)/2)
            units_dd_com.append((ce_unexp_com-ce_exp_com)/2)
            
        sims_dd_nav.append(units_dd_nav)  
        sims_dd_com.append(units_dd_com)   
        
        # check for fitting experimental data
        H,p=ranksums(units_dd_nav, dd_nav)
        if p>0.05:   # equal distributions
            stats_dd_nav.append(1)
        else:
            stats_dd_nav.append(0) 
            
        H,p=ranksums(units_dd_com, dd_com)
        if p>0.05:   # equal distributions
            stats_dd_com.append(1)
        else:
            stats_dd_com.append(0)     

#    # avearage context effect across units 
#    av_dd_nav=np.mean(sims_dd_nav,1)
#    av_dd_com=np.mean(sims_dd_com,1)
#    
#    # convert into 2D matrix results and stats 
#    av_dd_nav=np.reshape(av_dd_nav,(nx,ny)).T
#    av_dd_com=np.reshape(av_dd_com,(nx,ny)).T
#    stats_dd_nav=np.reshape(stats_dd_nav,(nx,ny)).T
#    stats_dd_com=np.reshape(stats_dd_com,(nx,ny)).T
#    
#    # subsets of simulations depending on parameters combinations 
    sims_dd_nav=np.array(sims_dd_nav)
    sims_dd_com=np.array(sims_dd_com)
#    # only those that are the same for LF and HF (diagonal)
#    diag_mask=np.eye(nx,dtype=bool)
#    diag=diag_mask.flatten()
#    same_nav=sims_dd_nav[diag]
#    same_com=sims_dd_com[diag]
#    # those that LF > HF (at the bottom of the diag)
#    bottom_mask=np.tril(~diag_mask)
#    bottom=bottom_mask.flatten()
#    lower_nav=sims_dd_nav[bottom]
#    lower_com=sims_dd_com[bottom]
#    # those that HF > LF (on top of the diag)
#    top_mask=np.triu(~diag_mask)
#    top=top_mask.flatten()
#    higher_nav=sims_dd_nav[top]
#    higher_com=sims_dd_com[top]
    
    # plotting
    A_nav=sims_dd_nav   #same_nav  #sims_dd_nav  # same_nav   # 
    A_com=sims_dd_com  #same_com  #sims_dd_com  # same_com
    
    fig, ax = plt.subplots(1,1,sharex=True,sharey=True,figsize=(4,4))
    nsims=A_nav.shape[0]
    colors = cm.rainbow(np.linspace(0, 1, nsims))
    for i in range(nsims):
        ax.scatter(A_nav[i,:],A_com[i,:],s=5,color=colors[i])
    minX=-0.5
    maxX=0.5
    ax.set_xlim(minX,maxX)
    ax.set_ylim(minX,maxX)
    ax.hlines(y=0,xmin=minX,xmax=maxX)
    ax.vlines(x=0,ymin=minX,ymax=maxX)
   
    return fig


def changes_dd(sims,ind=2):  # ind=2: selectivity; ind=1: selectivity - only one ; ind=3: anesthesia
    # plot d.d and spike counts during context changing either selectivity or synapses
    parameters1=sims[0]
    parameters2=sims[1]
    n_com=len(parameters1)
    
    if ind==2:   # when changing selectivity - includes navigation     
        n_nav=len(parameters2)
        sims = np.delete(sims,[0,1], 0)
    elif ind==1:       # when changing synaptic properties - do not include navigation
        n_nav=0
        sims = np.delete(sims,[0,1], 0)
    else:
        n_nav=0
        sims = np.delete(sims,[0], 0)
    # simulations changing selectivity of communication context 
    sp_nav1=[]
    sp_com1=[]
    dd_nav1=[]
    dd_com1=[]
    for s in range(n_com):
        s_i=sims[s]
        dd_nav_p=((s_i[:,1]-s_i[:,0])/2)
        dd_com_p=((s_i[:,3]-s_i[:,2])/2)
        sp_nav_p=(s_i[:,4])
        sp_com_p=(s_i[:,5])
        dd_nav1.append(dd_nav_p)
        dd_com1.append(dd_com_p)
        sp_nav1.append(sp_nav_p)
        sp_com1.append(sp_com_p)
    # simulations changing selectivity of navigation context 
    sp_nav2=[]
    sp_com2=[]
    dd_nav2=[]
    dd_com2=[]
    for j in range(n_nav):
        s_i=sims[s+1+j]
        dd_nav_p=((s_i[:,1]-s_i[:,0])/2)
        dd_com_p=((s_i[:,3]-s_i[:,2])/2)
        sp_nav_p=(s_i[:,4])
        sp_com_p=(s_i[:,5])
        dd_nav2.append(dd_nav_p)
        dd_com2.append(dd_com_p)
        sp_nav2.append(sp_nav_p)
        sp_com2.append(sp_com_p)    
    
    # plotting
    fig, ax = plt.subplots(1,2,sharex=True,figsize=(4,2))
    parameters1=np.linspace(0,1,n_com)
    parameters2=np.linspace(0,1,n_nav)
    
    ax[0].axhline(y=0,color='k',lw=0.5,ls='--')
#    ax[0].set_ylim(-0.05, 0.20)  
    ax[0].set_ylabel('d.d.')
    ax[0].axhline(y=0.09,color='b',lw=0.5,ls='--')
    ax[0].axhline(y=0.05,color='r',lw=0.5,ls='--')
    
    ax[0].plot(parameters1,np.mean(dd_com1,1),color='r',lw=0.75)
    ax[0].fill_between(parameters1, np.mean(dd_com1,1)-sem(dd_com1,1), np.mean(dd_com1,1)+sem(dd_com1,1), alpha=0.3, facecolor='r')
    if ind==2:
        ax[0].plot(parameters2,np.mean(dd_nav2,1),color='b',lw=0.75)
        ax[0].fill_between(parameters2, np.mean(dd_nav2,1)-sem(dd_nav2,1), np.mean(dd_nav2,1)+sem(dd_nav2,1) ,alpha=0.3, facecolor='b')
    else:
        ax[0].plot(parameters1,np.mean(dd_nav1,1),color='b',lw=0.75)
        ax[0].fill_between(parameters1, np.mean(dd_nav1,1)-sem(dd_nav1,1), np.mean(dd_nav1,1)+sem(dd_nav1,1) ,alpha=0.3, facecolor='b')

    ax[1].set_ylabel('# spikes context')
#    ax[1].set_ylim(80, 200)  
    ax[1].axhline(y=59,color='b',lw=0.5,ls='--')
    ax[1].axhline(y=81,color='r',lw=0.5,ls='--')
    ax[1].axhline(y=44,color='b',lw=0.5,ls='--')
    ax[1].axhline(y=72,color='r',lw=0.5,ls='--')
 
    ax[1].plot(parameters1,np.mean(sp_com1,1),color='r',lw=0.75) 
    ax[1].fill_between(parameters1, np.mean(sp_com1,1)-sem(sp_com1,1), np.mean(sp_com1,1)+sem(sp_com1,1) ,alpha=0.3, facecolor='r')
    if ind==2:
        ax[1].plot(parameters2,np.mean(sp_nav2,1),color='b',lw=0.75) 
        ax[1].fill_between(parameters2, np.mean(sp_nav2,1)-sem(sp_nav2,1), np.mean(sp_nav2,1)+sem(sp_nav2,1) ,alpha=0.3, facecolor='b')
    else:
        ax[1].plot(parameters1,np.mean(sp_nav1,1),color='b',lw=0.75) 
        ax[1].fill_between(parameters1, np.mean(sp_nav1,1)-sem(sp_nav1,1), np.mean(sp_nav1,1)+sem(sp_nav1,1) ,alpha=0.3, facecolor='b')
 
    ax = ax.flatten() 
    for axi in ax:
        axi.spines['top'].set_visible(False)
        axi.spines['right'].set_visible(False)
        axi.spines['bottom'].set_linewidth(0.75)
        axi.spines['left'].set_linewidth(0.75)
    
    if ind==2 or ind==1:
        plt.figtext(0.5,0, 'input selectivity', rotation=0, verticalalignment='bottom')
    else:
        plt.figtext(0.5,0, 'synaptic stregth', rotation=0, verticalalignment='bottom')
    
    fig.tight_layout(h_pad=3)  #w_pad=0.6,h_pad=0.6    
    
    return fig


def marbels_dd(sims):
    # plot spike counts during context in: cortex, LF input and HF input
    parameters1=sims[0]
    parameters2=sims[1]
    n_com=len(parameters1)
    n_nav=len(parameters2)
    sims = np.delete(sims,[0,1], 0)

    # simulations changing selectivity of communication context 
    dd_com_com=[]
    sp_com_com=[]
    sp_lf_com=[]
    sp_hf_com=[]
    for s in range(n_com):
        s_i=sims[s]
#        dd_nav=((s_i[:,1]-s_i[:,0])/2)
        dd_com=((s_i[:,3]-s_i[:,2])/2)
#        sp_nav=(s_i[:,4])
#        sp_LF=(s_i[:,5])
#        sp_HF=(s_i[:,6])
        sp_com=(s_i[:,7])
        sp_LF=(s_i[:,8])
        sp_HF=(s_i[:,9])
        dd_com_com.append(dd_com)
        sp_com_com.append(sp_com)
        sp_lf_com.append(sp_LF)
        sp_hf_com.append(sp_HF)

    # simulations changing selectivity of navigation context 
    dd_nav_nav=[]
    sp_nav_nav=[]
    sp_lf_nav=[]
    sp_hf_nav=[]
    for j in range(n_nav):
        s_i=sims[s+1+j]
        dd_nav=((s_i[:,1]-s_i[:,0])/2)
#        dd_com=((s_i[:,3]-s_i[:,2])/2)
        sp_nav=(s_i[:,4])
        sp_LF=(s_i[:,5])
        sp_HF=(s_i[:,6])
#        sp_com=(s_i[:,7])
#        sp_LF=(s_i[:,8])
#        sp_HF=(s_i[:,9])
        dd_nav_nav.append(dd_nav)
        sp_nav_nav.append(sp_nav)
        sp_lf_nav.append(sp_LF)
        sp_hf_nav.append(sp_HF)
    
    # plotting
    fig, ax = plt.subplots(1,2,sharex=True,figsize=(4,2))
    parameters1=np.linspace(0,1,n_com)
    parameters2=np.linspace(0,1,n_nav)
    
    # communication change selectivity
    a=sp_com_com
    color_i='k'
    ax2 = ax[0].twinx()
    ax2.plot(parameters1,np.mean(a,1),color=color_i,lw=0.75)
    ax2.fill_between(parameters1, np.mean(a,1)-sem(a,1), np.mean(a,1)+sem(a,1), alpha=0.3, facecolor=color_i)
    ax2.set_ylim(40,99)
    ax2.spines['top'].set_visible(False)
    ax2.set_ylabel('# spikes cortex')
    
    a=sp_lf_com
    color_i='r'
    ax[0].plot(parameters1,np.mean(a,1),color=color_i,lw=0.75)
    ax[0].fill_between(parameters1, np.mean(a,1)-sem(a,1), np.mean(a,1)+sem(a,1), alpha=0.3, facecolor=color_i)
    a=sp_hf_com
    color_i='b'
    ax[0].plot(parameters1,np.mean(a,1),color=color_i,lw=0.75)
    ax[0].fill_between(parameters1, np.mean(a,1)-sem(a,1), np.mean(a,1)+sem(a,1), alpha=0.3, facecolor=color_i)
    
#    ax[0].set_ylim(0,50)
    ax[0].set_ylabel('# spikes inputs')
    ax[0].set_xlabel('input selectivity to com')

    # navigation change selectivity
    a=sp_nav_nav
    color_i='k'
    ax2 = ax[1].twinx()
    ax2.plot(parameters1,np.mean(a,1),color=color_i,lw=0.75)
    ax2.fill_between(parameters1, np.mean(a,1)-sem(a,1), np.mean(a,1)+sem(a,1), alpha=0.3, facecolor=color_i)
    ax2.set_ylim(40,99)
    ax2.spines['top'].set_visible(False)
    ax2.set_ylabel('# spikes cortex')
    
    
    a=sp_lf_nav
    color_i='r'
    ax[1].plot(parameters1,np.mean(a,1),color=color_i,lw=0.75)
    ax[1].fill_between(parameters1, np.mean(a,1)-sem(a,1), np.mean(a,1)+sem(a,1), alpha=0.3, facecolor=color_i)
    a=sp_hf_nav
    color_i='b'
    ax[1].plot(parameters1,np.mean(a,1),color=color_i,lw=0.75)
    ax[1].fill_between(parameters1, np.mean(a,1)-sem(a,1), np.mean(a,1)+sem(a,1), alpha=0.3, facecolor=color_i)

#    ax[1].set_ylabel('# spikes context')
    ax[1].set_xlabel('input selectivity to nav')
      
    ax = ax.flatten() 
    for axi in ax:
        axi.spines['top'].set_visible(False)
        axi.spines['right'].set_visible(False)
        axi.spines['bottom'].set_linewidth(0.75)
        axi.spines['left'].set_linewidth(0.75)

    fig.tight_layout(h_pad=3)  #w_pad=0.6,h_pad=0.6    
    
    return fig

def sims_synaptic(N_units,k_low_com=0,k_high_com=0,k_high_nav=0,k_low_nav=0):
    # run 2 simulations varying contexts
    # outputs:
    # average of synaptic depression per synapse
    # average of synaptic weight per synapse 
    trials=N_units*20
    dur_ech=int(1.3311*10000)+1500
    dur_com=int(1.9266*10000)+1500

    masker=np.repeat(range(1,N_units+1),20)
    synapses=np.repeat(range(1,N_units+1),20)
    
    # SIMULATION 1: navigation context
    context=1
    exp=1
    exec(open("./input.py").read())  
    exec(open("./model.py").read()) 
    
    col=firing_context(ac,onset,context)
    colLF=firing_context(sp_low,onset,context)
    colHF=firing_context(sp_high,onset,context)
    masker=np.column_stack((masker, col))
    masker=np.column_stack((masker, colLF))
    masker=np.column_stack((masker, colHF))
    
    LF_xs=np.mean(mon.x_S[:trials,:dur_ech],1)
    HF_xs=np.mean(mon.x_S[trials:,:dur_ech],1)
    synapses=np.column_stack((synapses,LF_xs,HF_xs))
    LF_ge=np.mean(syn_w.g_e[:trials,:dur_ech]/psiemens,1)
    HF_ge=np.mean(syn_w.g_e[trials:,:dur_ech]/psiemens,1)
    synapses=np.column_stack((synapses,LF_ge,HF_ge))


    # SIMULATION 2: communication context
    context=2
    exp=1
    exec(open("./input.py").read())  
    exec(open("./model.py").read()) 
    
    col=firing_context(ac,onset,context)
    colLF=firing_context(sp_low,onset,context)
    colHF=firing_context(sp_high,onset,context)
    masker=np.column_stack((masker, col))
    masker=np.column_stack((masker, colLF))
    masker=np.column_stack((masker, colHF))
    
    LF_xs=np.mean(mon.x_S[:trials,:dur_com],1)
    HF_xs=np.mean(mon.x_S[trials:,:dur_com],1)
    synapses=np.column_stack((synapses,LF_xs,HF_xs))
    LF_ge=np.mean(syn_w.g_e[:trials,:dur_com]/psiemens,1)
    HF_ge=np.mean(syn_w.g_e[trials:,:dur_com]/psiemens,1)
    synapses=np.column_stack((synapses,LF_ge,HF_ge))

    
    return masker,synapses

def syn_variables(sims,con=0):
    # plot spike counts during context in: cortex, LF input and HF input
    # plot synaptic variables recorded: x_S and g_e
    # con=0 : results in response to echolocation context
    # con=1 : results in response to communication context
    parameters1=sims[0]
    n=len(parameters1)
    sims = np.delete(sims,[0], 0)
    if con==0: # in response to echo
        shift_sp=int(0)
        shift_syn=int(0)
    elif con==1: # in response to com
        shift_sp=int(3)
        shift_syn=int(4)


    # simulations changing selectivity of "con" context 
    sp_c=[]
    sp_LF=[]
    sp_HF=[]
    xs_av=[]
    xs_av_lf=[]
    xs_av_hf=[]
    ge_av=[]

    for s in range(n):
        s_i=sims[s]
        mean_i=np.mean(s_i,0) 
        sp_c_i=mean_i[1+shift_sp]    # cortical 
        sp_LF_i=mean_i[2+shift_sp]   # LF input
        sp_HF_i=mean_i[3+shift_sp]   # HF input
        xs_lf=(mean_i[8+shift_syn])
        xs_hf=(mean_i[9+shift_syn])
        xs_i=(mean_i[8+shift_syn]+mean_i[9+shift_syn])/2 # average between LF and HF
        ge_i=(mean_i[10+shift_syn]+mean_i[11+shift_syn])/2 # average between LF and HF
        sp_c.append(sp_c_i)
        sp_LF.append(sp_LF_i)
        sp_HF.append(sp_HF_i)
        xs_av.append(xs_i)
        xs_av_lf.append(xs_lf)
        xs_av_hf.append(xs_hf)
        ge_av.append(ge_i)

    # plotting
    fig, ax=plt.subplots(1,figsize=(3,2))    

    sum_inputs=[x + y for x, y in zip(sp_LF, sp_HF)]
#    ax.plot(parameters1,sp_c,'k')
    ax.plot(parameters1,sp_LF,'r')
    ax.plot(parameters1,sp_HF,'b')
    ax.plot(parameters1,sum_inputs,'k')

    ax2=ax.twinx()
#    ax2.plot(parameters1,ge_av/max(ge_av),'--g')
    ax2.plot(parameters1,[j*8 for j in xs_av_lf],'--r')
    ax2.plot(parameters1,[j*8 for j in xs_av_hf],'--b')
#    sum_xs=[x + y for x, y in zip(xs_av_lf, xs_av_hf)]
#    ax2.plot(parameters1,sum_xs,'--k')
    
    
    ax.set_xlabel('bandwidth of sound')    
    ax.set_ylabel('input spike count')
    ax2.set_ylabel('synaptic strength (nS)',color='k')
    ax2.set_ylim([5.8, 8])  #ax2.get_ylim()[0]


#    for tl in ax2.get_yticklabels():
#        tl.set_color('g')
    
    
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.spines['bottom'].set_linewidth(0.75)
    ax.spines['left'].set_linewidth(0.75)
    ax2.spines['top'].set_visible(False)
    ax2.spines['left'].set_visible(False)
    ax2.spines['bottom'].set_linewidth(0.75)
    ax2.spines['right'].set_linewidth(0.75)

    fig.tight_layout() 
    plt.locator_params(nbins=4)

    return fig

def compare_synweight(sims1,sims2):
    
    # plot spike counts during context in: cortex, LF input and HF input
    parameters1=sims1[0]
    parameters2=sims1[1]
    n_com=len(parameters1)
    n_nav=len(parameters2)
    sims1 = np.delete(sims1,[0,1], 0)
    sims2 = np.delete(sims2,[0,1], 0)
    
    # simulations changing selectivity of communication context 
    post_c_com1=[]
    post_c_com2=[]
    for s in range(n_com):
        s_i1=sims1[s]
        s_i2=sims2[s]
        mean_i1=np.mean(s_i1,0)
        mean_i2=np.mean(s_i2,0)
        post_com1=mean_i1[8]
        post_com2=mean_i2[8]
        post_c_com1.append(post_com1)
        post_c_com2.append(post_com2)

    # simulations changing selectivity of navigation context 
    post_n_nav1=[]
    post_n_nav2=[]
    for j in range(n_nav):
        s_i1=sims1[s+1+j]
        s_i2=sims2[s+1+j]
        mean_i1=np.mean(s_i1,0)
        mean_i2=np.mean(s_i2,0)
        post_nav1=mean_i1[6]
        post_nav2=mean_i2[6]
        post_n_nav1.append(post_nav1)
        post_n_nav2.append(post_nav2)
    
     # plotting the postsynaptic weight 
    fig, ax = plt.subplots(1,2,sharex=True,figsize=(4,2))
    parameters1=np.linspace(0,1,n_com)
    parameters2=np.linspace(0,1,n_nav)
    
    # communication change selectivity
    ax[0].plot(parameters1,post_c_com1,color='k',lw=0.75)
    ax[0].plot(parameters1,post_c_com2,color='r',lw=0.75)
    ax[0].set_xlabel('input selectivity to com') 
    ax[0].set_ylabel('postsynaptic conductance')
      

    # navigation change selectivity
    ax[1].plot(parameters2,post_n_nav1,color='k',lw=0.75)
    ax[1].plot(parameters2,post_n_nav2,color='r',lw=0.75)
    ax[1].set_xlabel('input selectivity to nav')   
    ax = ax.flatten() 
    for axi in ax:
        axi.spines['top'].set_visible(False)
        axi.spines['right'].set_visible(False)
        axi.spines['bottom'].set_linewidth(0.75)
        axi.spines['left'].set_linewidth(0.75)
        
    fig.tight_layout(h_pad=3)  #w_pad=0.6,h_pad=0.6 
    
    return fig 


def ce_synapses(sims): 
    # plot c.e changing synaptic weight
    parameters1=sims[0]
    n_com=len(parameters1)
    sims = np.delete(sims,[0], 0)
    
    # simulations changing selectivity of communication context 
    ce_nav_exp=[]
    ce_nav_unexp=[]
    ce_com_exp=[]
    ce_com_unexp=[]
    for s in range(n_com):
        s_i=sims[s]
        a=s_i[:,0]
        b=s_i[:,1]
        c=s_i[:,2]
        d=s_i[:,3]
        ce_nav_exp.append(a)
        ce_nav_unexp.append(b)
        ce_com_exp.append(c)
        ce_com_unexp.append(d)
 
    
    # plotting
    fig, ax = plt.subplots(1,2,sharex=True,sharey=True,figsize=(4,2))

    ax[0].axhline(y=0,color='k',lw=0.5,ls='--')
#    ax[0].set_ylim(-0.6, 0)  
    ax[0].set_ylabel('c.e.')
    ax[0].set_xlabel('synapse decrement (nav)')
    
    ax[0].plot(parameters1,np.mean(ce_nav_unexp,1),color='k',lw=0.75)
    ax[0].fill_between(parameters1, np.mean(ce_nav_unexp,1)-sem(ce_nav_unexp,1), np.mean(ce_nav_unexp,1)+sem(ce_nav_unexp,1), alpha=0.3, facecolor='k')
    ax[0].plot(parameters1,np.mean(ce_nav_exp,1),color='b',lw=0.75)
    ax[0].fill_between(parameters1, np.mean(ce_nav_exp,1)-sem(ce_nav_exp,1), np.mean(ce_nav_exp,1)+sem(ce_nav_exp,1), alpha=0.3, facecolor='b')

    ax[1].axhline(y=0,color='k',lw=0.5,ls='--')
    ax[1].set_xlabel('synapse decrement (com)')
    ax[1].plot(parameters1,np.mean(ce_com_unexp,1),color='k',lw=0.75)
    ax[1].fill_between(parameters1, np.mean(ce_com_unexp,1)-sem(ce_com_unexp,1), np.mean(ce_com_unexp,1)+sem(ce_com_unexp,1), alpha=0.3, facecolor='k')
    ax[1].plot(parameters1,np.mean(ce_com_exp,1),color='b',lw=0.75)
    ax[1].fill_between(parameters1, np.mean(ce_com_exp,1)-sem(ce_com_exp,1), np.mean(ce_com_exp,1)+sem(ce_nav_exp,1), alpha=0.3, facecolor='b')

    ax = ax.flatten() 
    for axi in ax:
        axi.spines['top'].set_visible(False)
        axi.spines['right'].set_visible(False)
        axi.spines['bottom'].set_linewidth(0.75)
        axi.spines['left'].set_linewidth(0.75)
    
    fig.tight_layout(h_pad=3)  #w_pad=0.6,h_pad=0.6    
    
    return fig

def plot_we_d(sims):
    parameters1=sims[0]
    parameters2=sims[1]
    sims = np.delete(sims,[0,1], 0)
    ntrials=len(sims[0])
    nx=len(parameters1)
    ny=len(parameters2)
    ntot=nx*ny
    N_units=int(ntrials/20)
    
    # calculate cliff delta from simulations per unit 
    sims_nav=[]
    sims_com=[]
    for s in range(ntot):
        s_i=np.mean(sims[s],0)
        # spikes in response to echolocation context
        sp_n_c=s_i[1]       # cortical neurons
        sp_n_LF=s_i[2]      # LF-tuned input
        sp_n_HF=s_i[3]      # HF-tuned input
        # spikes in response to communication context
        sp_c_c=s_i[4]
        sp_c_LF=s_i[5]
        sp_c_HF=s_i[6]
        
        sims_nav.append([sp_n_c,sp_n_LF,sp_n_HF])
        sims_com.append([sp_c_c,sp_c_LF,sp_c_HF])
    sims_nav=np.array(sims_nav)
    sims_com=np.array(sims_com)
    
    # convert into 2D matrix results and stats 
    # what to plot? cortical, LF input or HF input
    ind=1  # 1: cortical; 2: LF; 3: HF
    a=sims_nav[:,ind-1]
    b=sims_com[:,ind-1]
    av_nav=np.reshape(a,(nx,ny)).T
    av_com=np.reshape(b,(nx,ny)).T

    
    # plotting
    fig, ax = plt.subplots(1,2,sharex=True,sharey=True,figsize=(5.5,1.9))
    li=95      # max of colorbar

    nav=ax[0].imshow(av_nav,cmap='jet',origin='lower',vmin=0,vmax=li) 
    com=ax[1].imshow(av_com,cmap='jet',origin='lower',vmin=0,vmax=li)
    
    fig.colorbar(nav, ax=ax[0])
    fig.colorbar(com, ax=ax[1])
    fig.tight_layout() 
    
    fig2, ax = plt.subplots(1,3,sharex=True,sharey=True,figsize=(5.5,1.9))
    li=95      # max of colorbar

    fix_d=8
    ax[0].set_title('increasing synaptic weight')
    ax[0].plot(av_nav[:,fix_d],'b')
    ax[0].plot(av_com[:,fix_d]-20,'r')
    ax[2].set_title('increasing both (we and d)')
    ax[2].plot(np.fliplr(av_nav).diagonal(),'b')
    ax[2].plot(np.fliplr(av_com).diagonal()-20,'r')
    ax[1].set_title('increasing synaptic depression')
    ax[1].plot(av_nav[fix_d,:],'b')
    ax[1].plot(av_com[fix_d,:]-20,'r')
    
    fig2.tight_layout() 
    
    
    return fig, fig2 

def plot_v_alpha(sims,sims2):
    ratios=sims[0]
    vmax=sims[1]
    sims_com = np.delete(sims,[0,1], 0)
    sims_nav = np.delete(sims2,[0,1], 0)
    ntrials=len(sims_com[0])
    n_ratios=len(ratios)
    n_v=len(vmax)
    
    j=0
    ratios_com_c=[] # changing ratios ch1/ch2 in response to communication
    ratios_nav_c=[] # changing ratios ch1/ch2 in response to navigation
    ratios_com_LF=[] # changing ratios ch1/ch2 in response to communication
    ratios_nav_LF=[] # changing ratios ch1/ch2 in response to navigation
    ratios_com_HF=[] # changing ratios ch1/ch2 in response to communication
    ratios_nav_HF=[] # changing ratios ch1/ch2 in response to navigation
    for x in range(n_ratios):  # across simulations of ch1/ch2
        vs_per_ratio_com_c=[]
        vs_per_ratio_nav_c=[]
        vs_per_ratio_com_LF=[]
        vs_per_ratio_nav_LF=[]
        vs_per_ratio_com_HF=[]
        vs_per_ratio_nav_HF=[]
        for y in range(n_v):  # across simulations of v 
            # average across 1000 trials (50units*20trials)
            s_n=np.mean(sims_nav[j],0)
            s_c=np.mean(sims_com[j],0)
            # save average spikes in vectors per ratio
            vs_per_ratio_com_c.append(s_c[11]) # cortical neurons in response to com
            vs_per_ratio_nav_c.append(s_n[8])  # cortical neurons in response to nav
            vs_per_ratio_com_LF.append(s_c[12]) # LF input in response to com
            vs_per_ratio_nav_LF.append(s_n[9])  # LF input in response to nav
            vs_per_ratio_com_HF.append(s_c[13]) # HF input in response to com
            vs_per_ratio_nav_HF.append(s_n[10]) # HF input in response to nav
            j=j+1
        ratios_com_c.append(vs_per_ratio_com_c)
        ratios_nav_c.append(vs_per_ratio_nav_c)
        ratios_com_LF.append(vs_per_ratio_com_LF)
        ratios_nav_LF.append(vs_per_ratio_nav_LF)
        ratios_com_HF.append(vs_per_ratio_com_HF)
        ratios_nav_HF.append(vs_per_ratio_nav_HF)
                
    ## plotting spike counts during context across different ratios and v ##  
    dur_ech=1.3311
    dur_com=1.9266
    
    ratios_com_c=np.array(ratios_com_c)/dur_com
    ratios_nav_c=np.array(ratios_nav_c)/dur_ech
    
    
    fig1, ax=plt.subplots(2,3,sharex=True,figsize=(6.5,3))    

    colors = cm.jet(np.linspace(0,1,n_ratios))

    for i in range(n_ratios):
        ax[0,0].plot(vmax,ratios_com_c[i],color=colors[i],lw=1)
    ax[0,0].set_title('response to com')
    ax[0,0].set_ylabel('cortical')
    
    for i in range(n_ratios):
        ax[0,1].plot(vmax,ratios_com_LF[i],color=colors[i],lw=1)
    ax[0,1].set_title('response to com')
    ax[0,1].set_ylabel('LF inputs')


    for i in range(n_ratios):
        ax[0,2].plot(vmax,ratios_com_HF[i],color=colors[i],lw=1)
    ax[0,2].set_title('response to com')
    ax[0,2].set_ylabel('HF inputs')

    for i in range(n_ratios):
        ax[1,0].plot(vmax,ratios_nav_c[i],color=colors[i],lw=1)
    ax[1,0].set_title('response to nav')
    ax[1,0].set_ylabel('cortical')
    
    for i in range(n_ratios):
        ax[1,1].plot(vmax,ratios_nav_LF[i],color=colors[i],lw=1)
    ax[1,1].set_title('response to nav')
    ax[1,1].set_ylabel('LF inputs')


    for i in range(n_ratios):
        ax[1,2].plot(vmax,ratios_nav_HF[i],color=colors[i],lw=1)
    ax[1,2].set_title('response to nav')
    ax[1,2].set_ylabel('HF inputs')

     
    # add experimental data
#    ax[0,0].axhline(y=59,color='b',lw=0.5,ls='--')
    ax[0,0].axhline(y=81/dur_com,color='r',lw=0.5,ls='--')
#    ax[0,0].axhline(y=44,color='b',lw=0.5,ls='--')
    ax[0,0].axhline(y=72/dur_com,color='r',lw=0.5,ls='--')
    
    ax[1,0].axhline(y=59/dur_ech,color='b',lw=0.5,ls='--')
#    ax[1,0].axhline(y=81,color='r',lw=0.5,ls='--')
    ax[1,0].axhline(y=44/dur_ech,color='b',lw=0.5,ls='--')
#    ax[1,0].axhline(y=72,color='r',lw=0.5,ls='--')
    
    ax = ax.flatten() 
    for axi in ax:
        axi.spines['top'].set_visible(False)
        axi.spines['right'].set_visible(False)
        axi.spines['bottom'].set_linewidth(0.75)
        axi.spines['left'].set_linewidth(0.75)
        
    fig1.tight_layout()  
    
    ## plotting changes across different ratios for a fixed v ##
    # in response to communication
    fig2, ax=plt.subplots(1,figsize=(3,2))    
    v_inx=-1
    alpha_com_c=[vi[v_inx] for vi in ratios_com_c]
    alpha_com_LH=[vi[v_inx] for vi in ratios_com_LF]
    alpha_com_HF=[vi[v_inx] for vi in ratios_com_HF]
    sum_inputs_com=[x + y for x, y in zip(alpha_com_LH, alpha_com_HF)]
    
    ax.plot(ratios,alpha_com_LH,'r')
    ax.plot(ratios,alpha_com_HF,'b')
    ax.plot(ratios,sum_inputs_com,'k')
#    ax.plot(alpha_com_c,'k')
    
    ax.set_xlabel('bandwidth of communication')    
    ax.set_ylabel('input spike count')
    
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.spines['bottom'].set_linewidth(0.75)
    ax.spines['left'].set_linewidth(0.75)

    fig2.tight_layout()     
    
    # in response to navigation 
    fig3, ax=plt.subplots(1,figsize=(3,2))  
    alpha_nav_c=[vi[v_inx] for vi in ratios_nav_c]
    alpha_nav_LH=[vi[v_inx] for vi in ratios_nav_LF]
    alpha_nav_HF=[vi[v_inx] for vi in ratios_nav_HF]
    sum_inputs_nav=[x + y for x, y in zip(alpha_nav_LH, alpha_nav_HF)]
    
    ax.plot(ratios,alpha_nav_LH,'r')
    ax.plot(ratios,alpha_nav_HF,'b')
    ax.plot(ratios,sum_inputs_nav,'k')
    
    ax.set_ylabel('input spike count')
    ax.set_xlabel('bandwidth of sound')

    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.spines['bottom'].set_linewidth(0.75)
    ax.spines['left'].set_linewidth(0.75)
    
    fig3.tight_layout()   
    
#    fig4, ax=plt.subplots(1,1,figsize=(5,2.5))  
#    ax.plot(vmax[0:-1],np.diff(ratios_com_c[9])/np.diff(vmax),'r')
#    yy=[f for f in ratios_nav_c[0]]
#    ax.plot(vmax[0:-1],np.diff(yy)/np.diff(vmax),'b')
    
    fig5, ax=plt.subplots(1,figsize=(3,2))  
    ax.plot(ratios,alpha_nav_c,'k')
    
    ax.set_ylabel('output spike count')
    ax.set_xlabel('bandwidth of echolocation')

    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.spines['bottom'].set_linewidth(0.75)
    ax.spines['left'].set_linewidth(0.75)
    
    fig5.tight_layout()   
    
    return fig3
    
def first_approx(sims):
    # for the first plots on the paper first approximation to anesthesia effects 
    parameters1=sims[0]
    parameters2=sims[1]
    sims = np.delete(sims,[0,1], 0)
    nx=len(parameters1)
    ny=len(parameters2)
    dur_ech=1#1.3311
    dur_com=1#1.9266
    
    # for plotting 11 parameters average+- std for one fixed x 
    x_fixed=10 
    dd_nav_fixedx=[]
    dd_com_fixedx=[]
    di_nav_fixedx=[]
    di_iso_fixedx=[]
    di_com_fixedx=[]
    ce_n_e_fixedx=[]
    ce_n_u_fixedx=[]
    ce_c_e_fixedx=[]
    ce_c_u_fixedx=[]
    sp_nav_fixedx=[]
    sp_com_fixedx=[]
    # ...and for one fixed y
    y_fixed=0
    dd_nav_fixedy=[]
    dd_com_fixedy=[]
    di_nav_fixedy=[]
    di_iso_fixedy=[]
    di_com_fixedy=[]
    ce_n_e_fixedy=[]
    ce_n_u_fixedy=[]
    ce_c_e_fixedy=[]
    ce_c_u_fixedy=[]
    sp_nav_fixedy=[]
    sp_com_fixedy=[]    
    
    # for plotting 2D only the average across units 
    dd_nav=[]
    dd_com=[]
    di_nav=[]
    di_iso=[]
    di_com=[]
    ce_n_e=[]
    ce_n_u=[]
    ce_c_e=[]
    ce_c_u=[]
    rate_nav=[]
    rate_com=[]  

   
    j=0
    for x in range(nx):  # we or v 
        y_var_dd_nav=[]
        y_var_dd_com=[]
        y_var_di_nav=[]
        y_var_di_iso=[]
        y_var_di_com=[]
        y_var_ce_n_e=[]
        y_var_ce_n_u=[]
        y_var_ce_c_e=[]
        y_var_ce_c_u=[]
        y_var_sp_nav=[]
        y_var_sp_com=[]   
        for y in range(ny):   # d
            # j is a particular combinations of parameters
            # 1D needs all the values 
            if x==x_fixed:
                # context effect
                ce_fixed=context_effect(sims[j])
                ce_n_e_fixedx.append(ce_fixed[:,0])
                ce_n_u_fixedx.append(ce_fixed[:,1])
                ce_c_e_fixedx.append(ce_fixed[:,2])
                ce_c_u_fixedx.append(ce_fixed[:,3])
                # deviance detection
                dd_nav_fixed=(ce_fixed[:,1]-ce_fixed[:,0])/2
                dd_com_fixed=(ce_fixed[:,3]-ce_fixed[:,2])/2
                dd_nav_fixedx.append(dd_nav_fixed)
                dd_com_fixedx.append(dd_com_fixed)
                # discriminability index 
                di_fixed=discriminability_index(sims[j])
                di_nav_fixedx.append(di_fixed[:,0])
                di_iso_fixedx.append(di_fixed[:,1])
                di_com_fixedx.append(di_fixed[:,2])
                # firing rate context 
                sp_fixed=sims[j]
                sp_nav=np.reshape(sp_fixed[:,8],(50,20))  
                sp_nav_fixedx.append(np.mean(sp_nav,1)/dur_ech)
                sp_com=np.reshape(sp_fixed[:,11],(50,20))    
                sp_com_fixedx.append(np.mean(sp_com,1)/dur_com)
                
            if y==y_fixed:   
                # context effect
                ce_fixed=context_effect(sims[j])
                ce_n_e_fixedy.append(ce_fixed[:,0])
                ce_n_u_fixedy.append(ce_fixed[:,1])
                ce_c_e_fixedy.append(ce_fixed[:,2])
                ce_c_u_fixedy.append(ce_fixed[:,3])
                # deviance detection
                dd_nav_fixed=(ce_fixed[:,1]-ce_fixed[:,0])/2
                dd_com_fixed=(ce_fixed[:,3]-ce_fixed[:,2])/2
                dd_nav_fixedy.append(dd_nav_fixed)
                dd_com_fixedy.append(dd_com_fixed) 
                # discriminability index 
                di_fixed=discriminability_index(sims[j])
                di_nav_fixedy.append(di_fixed[:,0])
                di_iso_fixedy.append(di_fixed[:,1])
                di_com_fixedy.append(di_fixed[:,2])
                # firing rate context 
                sp_fixed=sims[j]
                sp_nav=np.reshape(sp_fixed[:,8],(50,20))  
                sp_nav_fixedy.append(np.mean(sp_nav,1)/dur_ech)
                sp_com=np.reshape(sp_fixed[:,11],(50,20))    
                sp_com_fixedy.append(np.mean(sp_com,1)/dur_com)
            
            # 2D only needs the average
            ce_i=np.mean(context_effect(sims[j]),0) 
            dd_nav_i=(ce_i[1]-ce_i[0])/2
            dd_com_i=(ce_i[3]-ce_i[2])/2
            y_var_dd_nav.append(dd_nav_i)
            y_var_dd_com.append(dd_com_i)
            # context effect 
            y_var_ce_n_e.append(ce_i[0])
            y_var_ce_n_u.append(ce_i[1])
            y_var_ce_c_e.append(ce_i[2])
            y_var_ce_c_u.append(ce_i[3])
            # discriminability index 
            di_i=np.mean(discriminability_index(sims[j]),0) 
            y_var_di_nav.append(di_i[0])
            y_var_di_iso.append(di_i[1])
            y_var_di_com.append(di_i[2])
            # firing rate 
            sp_contexts=sims[j]
            y_var_sp_nav.append(np.mean(sp_contexts[:,8])/dur_ech)
            y_var_sp_com.append(np.mean(sp_contexts[:,11])/dur_com)
            
            j=j+1
            
        dd_nav.append(y_var_dd_nav)
        dd_com.append(y_var_dd_com)
        di_nav.append(y_var_di_nav)
        di_iso.append(y_var_di_iso)
        di_com.append(y_var_di_com)
        ce_n_e.append(y_var_ce_n_e)
        ce_n_u.append(y_var_ce_n_u)
        ce_c_e.append(y_var_ce_c_e)
        ce_c_u.append(y_var_ce_c_u)
        rate_nav.append(y_var_sp_nav)
        rate_com.append(y_var_sp_com)

        
    dd_nav=np.array(dd_nav)
    dd_com=np.array(dd_com)
    di_nav=np.array(di_nav)
    di_iso=np.array(di_iso)
    di_com=np.array(di_com)
    ce_n_e=np.array(ce_n_e)
    ce_n_u=np.array(ce_n_u)
    ce_c_e=np.array(ce_c_e)
    ce_c_u=np.array(ce_c_u)
    rate_nav=np.array(rate_nav)
    rate_com=np.array(rate_com)
    
#    parameters1=np.arange(0.1,1.1,0.1)
    orr=0
    # plotting 1D with one parameter fixed
    fig1, ax = plt.subplots(1,figsize=(1.75,1.1))
    
    A = dd_nav_fixedx
    B = dd_com_fixedx
    X = parameters2
    
    ax.plot(X,np.mean(A,1),color='b',lw=0.75,marker="o",markersize=2)
    ax.fill_between(X, np.mean(A,1)-sem(A,1), np.mean(A,1)+sem(A,1), alpha=0.3, facecolor='b')
    ax.plot(X,np.mean(B,1),color='r',lw=0.75,marker="o",markersize=2)
    ax.fill_between(X, np.mean(B,1)-sem(B,1), np.mean(B,1)+sem(B,1), alpha=0.3, facecolor='r')

    # statistics
#    ax2 = ax.twinx()
#    stats=[comp_exp(par,'dd_nav','awake') for par in A]
#    ax2.plot(X,stats,'b--')
#    stats=[comp_exp(par,'dd_nav','anesthesia') for par in A]
#    ax2.plot(X,stats,'r--')
    
    # add errorbar for experimental data
    mat_contents = sio.loadmat(r'E:\Users\User\Documents\model-context\ei_data-an.mat')
    a=mat_contents['ei']
    exp_n_e = a[:,0]
    exp_n_u = a[:,1]
    exp_c_e = a[:,2]
    exp_c_u = a[:,3]
    exp_dd_navAN=(exp_n_u-exp_n_e)/2
    exp_dd_comAN=(exp_c_u-exp_c_e)/2
    mat_contents = sio.loadmat(r'E:\Users\User\Documents\model-context\ei_data.mat')
    a=mat_contents['ei']
    exp_n_e = a[:,0]
    exp_n_u = a[:,1]
    exp_c_e = a[:,2]
    exp_c_u = a[:,3]
    exp_dd_navAW=(exp_n_u-exp_n_e)/2
    exp_dd_comAW=(exp_c_u-exp_c_e)/2
#    ax.errorbar([2000,1200],[np.mean(exp_dd_navAW),np.mean(exp_dd_navAN)], yerr=[sem(exp_dd_navAW),sem(exp_dd_navAN)], fmt='o',color='b', mfc='w',ms=4)     
#    ax.errorbar([2000,1200],[np.mean(exp_dd_comAW),np.mean(exp_dd_comAN)], yerr=[sem(exp_dd_comAW),sem(exp_dd_comAN)], fmt='o',color='r', mfc='w',ms=4)


    
#    ax.set_ylabel('s.s.s')
#    ax.set_xlabel(r'inputs firing rate')
#    ax.set_xlabel(r'synaptic weight (nS)')
#    ax.set_xlabel(r'synaptic depression ($\Delta_{i,j}$)')
#    ax.set_xlabel(r'synaptic depression ($\Omega_{i,j}$)')
#    ax.set_xlabel(r'neuronal adaptation ($D_{th}$)')
#    ax.set_xlabel('bandwidth of sound')
#    ax.invert_xaxis()
#    ax.set_xticklabels([])
#    ax.axhline(0,linewidth=1, color='k',ls='--')
    
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.spines['bottom'].set_linewidth(0.75)
    ax.spines['left'].set_linewidth(0.75)

    fig1.tight_layout() 
    
    fig2=assymetries_plot(X,A,B,orr,'b','r',30)
    
    # another parameter
    fig3, ax = plt.subplots(1,figsize=(1.75,1.1))
    
    A = di_nav_fixedx
    B = di_iso_fixedx
    C = di_com_fixedx
    X = parameters2
    
    ax.plot(X,np.mean(A,1),color='b',lw=0.75,marker="o",markersize=2)
    ax.fill_between(X, np.mean(A,1)-sem(A,1), np.mean(A,1)+sem(A,1), alpha=0.3, facecolor='b')
    ax.plot(X,np.mean(B,1),color='K',lw=0.75,marker="o",markersize=2)
    ax.fill_between(X, np.mean(B,1)-sem(B,1), np.mean(B,1)+sem(B,1), alpha=0.3, facecolor='K')
    ax.plot(X,np.mean(C,1),color='r',lw=0.75,marker="o",markersize=2)
    ax.fill_between(X, np.mean(C,1)-sem(C,1), np.mean(C,1)+sem(C,1), alpha=0.3, facecolor='r')
    
    # statistics
#    ax2 = ax.twinx()
#    stats=[comp_exp(par,'di_iso','awake') for par in B]
#    ax2.plot(X,stats,'b--')
#    stats=[comp_exp(par,'di_iso','anesthesia') for par in B]
#    ax2.plot(X,stats,'r--')
    
#    ax.set_ylabel('d.i.')
#    ax.set_xlabel(r'inputs firing rate')
#    ax.set_xlabel(r'synaptic weight ($w_{e}$)')
#    ax.set_xlabel(r'synaptic depression ($\Delta_{i,j}$)')
#    ax.set_xlabel(r'synaptic depression ($\Omega_{i,j}$)')
#    ax.set_xlabel(r'neuronal adaptation ($D_{th}$)')
#    ax.invert_xaxis()
    
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.spines['bottom'].set_linewidth(0.75)
    ax.spines['left'].set_linewidth(0.75)

    ax.set_xticklabels([])

    fig3.tight_layout() 
    
    fig4=assymetries_plot(X,A,C,orr,'b','r',30)
    
    # another parameter    
    fig5, ax = plt.subplots(1,figsize=(1.75,1.1))
    
    A = sp_nav_fixedx
    B = sp_com_fixedx
    X = parameters2
    
    ax.plot(X,np.mean(A,1),color='b',lw=0.75,marker="o",markersize=2)
    ax.fill_between(X, np.mean(A,1)-sem(A,1), np.mean(A,1)+sem(A,1), alpha=0.3, facecolor='b')
    ax.plot(X,np.mean(B,1),color='r',lw=0.75,marker="o",markersize=2)
    ax.fill_between(X, np.mean(B,1)-sem(B,1), np.mean(B,1)+sem(B,1), alpha=0.3, facecolor='r')
    
#    ax2 = ax.twinx()
#    stats=[comp_exp(par,'rate_nav','awake') for par in A]
#    ax2.plot(X,stats,'b--')
#    stats=[comp_exp(par,'rate_nav','anesthesia') for par in A]
#    ax2.plot(X,stats,'r--')
    
#    ax.set_ylabel('firing context')
#    ax.set_xlabel(r'inputs firing rate')
#    ax.set_xlabel(r'synaptic weight ($w_{e}$)')
#    ax.set_xlabel(r'synaptic adaptation')  #($\Delta_{i,j}$)
#    ax.set_xlabel(r'synaptic depression ($\Omega_{i,j}$)')
#    ax.set_xlabel(r'neuronal adaptation')  #($D_{th}$)
#    ax.invert_xaxis()
    
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.spines['bottom'].set_linewidth(0.75)
    ax.spines['left'].set_linewidth(0.75)

    fig5.tight_layout() 
    
    fig6=assymetries_plot(X,A,B,orr,'b','r',45)
    
    # another parameter
    fig7, ax = plt.subplots(1,1,figsize=(1.75,1.1))
    
#    ax[0].set_title('echolocation')
    A = ce_n_e_fixedx
    B = ce_n_u_fixedx
    X = parameters2
    
    ax.plot(X,np.mean(A,1),color='cornflowerblue',lw=0.75,marker="o",markersize=2)
    ax.fill_between(X, np.mean(A,1)-sem(A,1), np.mean(A,1)+sem(A,1), alpha=0.3, facecolor='cornflowerblue')
    ax.plot(X,np.mean(B,1),color='darkblue',lw=0.75,marker="o",markersize=2)
    ax.fill_between(X, np.mean(B,1)-sem(B,1), np.mean(B,1)+sem(B,1), alpha=0.3, facecolor='darkblue')

#    ax2 = ax[0].twinx()
#    stats=[comp_exp(par,'ce_n_e','awake') for par in A]
#    ax2.plot(X,stats,'b--')
#    stats=[comp_exp(par,'ce_n_e','anesthesia') for par in A]
#    ax2.plot(X,stats,'r--')
    
#    ax.set_ylabel('echolocation effect')
#    ax.set_xlabel(r'inputs firing rate')
#    ax[0].set_xlabel(r'synaptic weight ($w_{e}$)')
#    ax[0].set_xlabel(r'synaptic depression ($\Delta_{i,j}$)')
#    ax[0].set_xlabel(r'synaptic depression ($\Omega_{i,j}$)')
#    ax[0].set_xlabel(r'neuronal adaptation ($D_{th}$)')
#    ax.invert_xaxis()
#    ax.set_xticklabels([])
#    ax[0].axhline(0,linewidth=1, color='k',ls='--')
    
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.spines['bottom'].set_linewidth(0.75)
    ax.spines['left'].set_linewidth(0.75)
    
#    ax[1].set_title('communication')
    fig7.tight_layout() 
    
    fig8=assymetries_plot(X,A,B,orr,'cornflowerblue','darkblue',30)
    
    fig9, ax = plt.subplots(1,1,figsize=(1.75,1.1))
    A = ce_c_e_fixedx
    B = ce_c_u_fixedx
    X = parameters2
    
    ax.plot(X,np.mean(A,1),color='indianred',lw=0.75,marker="o",markersize=2)
    ax.fill_between(X, np.mean(A,1)-sem(A,1), np.mean(A,1)+sem(A,1), alpha=0.3, facecolor='indianred')
    ax.plot(X,np.mean(B,1),color='darkred',lw=0.75,marker="o",markersize=2)
    ax.fill_between(X, np.mean(B,1)-sem(B,1), np.mean(B,1)+sem(B,1), alpha=0.3, facecolor='darkred')

#    ax2 = ax[1].twinx()
#    stats=[comp_exp(par,'ce_c_e','awake') for par in A]
#    ax2.plot(X,stats,'b--')
#    stats=[comp_exp(par,'ce_c_e','anesthesia') for par in A]
#    ax2.plot(X,stats,'r--')
    
#    ax.set_ylabel('communication effect')
#    ax.set_xlabel(r'inputs firing rate')
#    ax[1].set_xlabel(r'synaptic weight ($w_{e}$)')
#    ax[1].set_xlabel(r'synaptic depression ($\Delta_{i,j}$)')
#    ax[1].set_xlabel(r'synaptic depression ($\Omega_{i,j}$)')
#    ax[1].set_xlabel(r'neuronal adaptation ($D_{th}$)')
#    ax.invert_xaxis()
#    ax[1].axhline(0,linewidth=1, color='k',ls='--')
#    ax.set_xticklabels([])
    
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.spines['bottom'].set_linewidth(0.75)
    ax.spines['left'].set_linewidth(0.75)
    

    fig9.tight_layout() 
    
    fig10=assymetries_plot(X,A,B,orr,'indianred','darkred',30)
    
    
    # plottign 2D
    '''
    fig2, ax = plt.subplots(1,2,figsize=(5,2))
    A=rate_nav
    B=rate_com
    jump=2
    mm=min(np.min(A),np.min(B))
    MM=max(np.max(A),np.max(B))
    
    ax[0].imshow(A,cmap='jet',origin='lower',vmin=mm,vmax=MM)
    ax[0].set_xticks(np.arange(0,ny,jump))
    ax[0].set_xticklabels(['{:.1f}'.format(x) for x in parameters2[::jump]])
#    ax[0].invert_xaxis()
#    ax[0].set_yticks(np.arange(0,nx,jump))
    ax[0].set_yticklabels(['{:.1f}'.format(x) for x in parameters1[::jump]])
    ax[0].set_title('echolocation')
    
    im_com=ax[1].imshow(B,cmap='jet',origin='lower',vmin=mm,vmax=MM)
    ax[1].set_xticks(np.arange(0,ny,jump))
    ax[1].set_xticklabels(['{:.1f}'.format(x) for x in parameters2[::jump]])
#    ax[1].invert_xaxis()
#    ax[1].set_yticks(np.arange(0,nx,jump))
    ax[1].set_yticklabels(['{:.1f}'.format(x) for x in parameters1[::jump]])
    ax[1].set_title('communication')
    
#    ax[0].set_xlabel(r'synaptic weight ($w_{e}$)')
#    ax[1].set_xlabel(r'synaptic weight ($w_{e}$)')
    ax[1].set_xlabel(r'synaptic depression ($\Delta_{i,j}$)')
    ax[0].set_xlabel(r'synaptic depression ($\Delta_{i,j}$)')
#    ax[0].set_ylabel(r'input LF firing rate ($v$)')
    ax[0].set_ylabel(r'postsynaptic adaptation ($V_{th}$)')
    
    fig2.tight_layout() 
    
    fig2.subplots_adjust(right=0.85)
    cbar_ax = fig2.add_axes([0.88, 0.15, 0.04, 0.7])
    bar=fig2.colorbar(im_com, cax=cbar_ax)
    bar.ax.set_title('rate')
    '''
    
    return fig1,fig2,fig5,fig6,fig7,fig8,fig9,fig10


def comp_exp(distribution,parameter,state):
    # compare experimental data with modeling 
    # distribution: vector with simulated parameters 
    # parameter (11): 'dd_nav','di_nav', 'ce_n_e' or 'rate_nav' and others
    # state: 'awake' or 'anesthesia'
    # outputs: a (whether reject or no same distributions), p-values, effect size
    
    if (parameter[0:2]=='dd' or parameter[0:2]=='ce') and state=='awake':
        # import experimental awake data
        mat_contents = sio.loadmat(r'E:\Users\User\Documents\model-context\ei_data.mat')
        a=mat_contents['ei']
        exp_n_e = a[:,0]
        exp_n_u = a[:,1]
        exp_c_e = a[:,2]
        exp_c_u = a[:,3]
        exp_dd_nav=(exp_n_u-exp_n_e)/2
        exp_dd_com=(exp_c_u-exp_c_e)/2
        
        if parameter=='dd_nav':
            H,p=ranksums(distribution, exp_dd_nav)
            a=1 if p>0.05 else 0
            c,e = cliffsDelta(distribution, exp_dd_nav) 
        elif parameter=='dd_com':
            H,p=ranksums(distribution, exp_dd_com)
            a=1 if p>0.05 else 0 
            c,e = cliffsDelta(distribution, exp_dd_com)
        elif parameter=='ce_n_e':
            H,p=ranksums(distribution, exp_n_e)
            a=1 if p>0.05 else 0 
            c,e = cliffsDelta(distribution, exp_n_e)
        elif parameter=='ce_n_u':
            H,p=ranksums(distribution, exp_n_u)
            a=1 if p>0.05 else 0 
            c,e = cliffsDelta(distribution, exp_n_u)
        elif parameter=='ce_c_e':
            H,p=ranksums(distribution, exp_c_e)
            a=1 if p>0.05 else 0  
            c,e = cliffsDelta(distribution, exp_c_e)
        elif parameter=='ce_c_u':
            H,p=ranksums(distribution, exp_c_u)
            a=1 if p>0.05 else 0  
            c,e = cliffsDelta(distribution, exp_c_u)
        
    elif (parameter[0:2]=='dd' or parameter[0:2]=='ce') and state=='anesthesia':
        # import experimental anesthetized data 
        mat_contents = sio.loadmat(r'E:\Users\User\Documents\model-context\ei_data-an.mat')
        a=mat_contents['ei']
        exp_n_e = a[:,0]
        exp_n_u = a[:,1]
        exp_c_e = a[:,2]
        exp_c_u = a[:,3]
        exp_dd_nav=(exp_n_u-exp_n_e)/2
        exp_dd_com=(exp_c_u-exp_c_e)/2
        
        if parameter=='dd_nav':
            H,p=ranksums(distribution, exp_dd_nav)
            a=1 if p>0.05 else 0  
            c,e = cliffsDelta(distribution, exp_dd_nav)
        elif parameter=='dd_com':
            H,p=ranksums(distribution, exp_dd_com)
            a=1 if p>0.05 else 0  
            c,e = cliffsDelta(distribution, exp_dd_com)
        elif parameter=='ce_n_e':
            H,p=ranksums(distribution, exp_n_e)
            a=1 if p>0.05 else 0  
            c,e = cliffsDelta(distribution, exp_n_e)
        elif parameter=='ce_n_u':
            H,p=ranksums(distribution, exp_n_u)
            a=1 if p>0.05 else 0 
            c,e = cliffsDelta(distribution, exp_n_u)
        elif parameter=='ce_c_e':
            H,p=ranksums(distribution, exp_c_e)
            a=1 if p>0.05 else 0 
            c,e = cliffsDelta(distribution, exp_c_e)
        elif parameter=='ce_c_u':
            H,p=ranksums(distribution, exp_c_u)
            a=1 if p>0.05 else 0  
            c,e = cliffsDelta(distribution, exp_c_u)
        
    elif parameter[0:2]=='di' and state=='awake':
        # import experimental awake data - cliff delta values for 45 units 
        mat_contents = sio.loadmat(r'E:\Users\User\Documents\model-context\cliff_data.mat')
        a=mat_contents['cliff_data']
        exp_di_iso = a[:,0]
        exp_di_nav = a[:,1]
        exp_di_com = a[:,2]
        
        if parameter=='di_nav':
            H,p=ranksums(distribution, exp_di_nav)
            a=1 if p>0.05 else 0  
            c,e = cliffsDelta(distribution, exp_di_nav)
        elif parameter=='di_iso':
            H,p=ranksums(distribution, exp_di_iso)
            a=1 if p>0.05 else 0  
            c,e = cliffsDelta(distribution, exp_di_iso)
        elif parameter=='di_com':
            H,p=ranksums(distribution, exp_di_com)
            a=1 if p>0.05 else 0  
            c,e = cliffsDelta(distribution, exp_di_com)
            
    elif parameter[0:2]=='di' and state=='anesthesia':
        # import experimental anesthetized data - cliff delta values for 46 units 
        mat_contents = sio.loadmat(r'E:\Users\User\Documents\model-context\cliff_data-an.mat')
        a=mat_contents['cliff_data']
        exp_di_iso = a[:,0]
        exp_di_nav = a[:,1]
        exp_di_com = a[:,2]
        
        if parameter=='di_nav':
            H,p=ranksums(distribution, exp_di_nav)
            a=1 if p>0.05 else 0  
            c,e = cliffsDelta(distribution, exp_di_nav)
        elif parameter=='di_iso':
            H,p=ranksums(distribution, exp_di_iso)
            a=1 if p>0.05 else 0 
            c,e = cliffsDelta(distribution, exp_di_iso)
        elif parameter=='di_com':
            H,p=ranksums(distribution, exp_di_com)
            a=1 if p>0.05 else 0 
            c,e = cliffsDelta(distribution, exp_di_com)
            
    elif parameter[0:4]=='rate' and state=='awake':
        # import experimental awake data 
        mat_contents = sio.loadmat(r'E:\Users\User\Documents\model-context\rate_context.mat')
        a=mat_contents['rate']
        exp_rate_nav = a[:,0]
        exp_rate_com = a[:,1]
        
        if parameter=='rate_nav':
            H,p=ranksums(distribution, exp_rate_nav)
            a=1 if p>0.05 else 0  
            c,e = cliffsDelta(distribution, exp_rate_nav)
        elif parameter=='rate_com':
            H,p=ranksums(distribution, exp_rate_com)
            a=1 if p>0.05 else 0  
            c,e = cliffsDelta(distribution, exp_rate_com)
            
    elif parameter[0:4]=='rate' and state=='anesthesia':
        # import experimental awake data 
        mat_contents = sio.loadmat(r'E:\Users\User\Documents\model-context\rate_context-an.mat')
        a=mat_contents['rate']
        exp_rate_nav = a[:,0]
        exp_rate_com = a[:,1]
        
        if parameter=='rate_nav':
            H,p=ranksums(distribution, exp_rate_nav)
            a=1 if p>0.05 else 0  
            c,e = cliffsDelta(distribution, exp_rate_nav)
        elif parameter=='rate_com':
            H,p=ranksums(distribution, exp_rate_com)
            a=1 if p>0.05 else 0  
            c,e = cliffsDelta(distribution, exp_rate_com)
            
    return a #-np.log(p)            


def second_approx(sims,slope_rows=0,slope_cols=0,n_combis=0):
    # for the second plots on the paper caviets of parameters and solutions, anesthesia effects 
    parameters1=sims[0]
    parameters2=sims[1]
    sims = np.delete(sims,[0,1], 0)    
    L_rows=len(parameters1)
    L_cols=len(parameters2)
    dur_ech=1.3311
    dur_com=1.9266
    # create a matrix with the combination of parameters used in the same order 
    sims_par=[]
    for x in parameters1:   
        for y in parameters2:
            sims_par.append((x,y))        
    sims_par=np.array(sims_par)
    # picking certain parameters in a diagonal of the 2D parameters arrange
    # to change the slope of the diagonal 
#    slope_rows=4
#    slope_cols=1
#    n_combis=3
    
    ind=np.size(sims_par,0)-1 # position of the original combination -> awake
    diago=[ind]    # add this combination as the first of the list 
    for par in range(n_combis-1):
        ind_i=ind + slope_cols*L_cols- slope_rows
        diago.append(ind_i)
        ind=ind_i
    # to plot the location of the parameters chosen
    h=[0]*(L_rows*L_cols)
    for i in diago: 
        h[i]=1
    h=np.array(h)
    plt.figure()
    plt.imshow(h.reshape((L_rows,L_cols)),origin='lower')
    plt.xlabel('parameters 2')
    plt.ylabel('parameters 1')
    # to print the value of the parameters chosen 
    for c in diago:
        print('parameter1='+str(sims_par[c,0])+' parameter2='+str(sims_par[c,1]))
    # vector with parameters to  plot
    dd_nav=[]
    dd_com=[]
    di_nav=[]
    di_iso=[]
    di_com=[]
    ce_n_e=[]
    ce_n_u=[]
    ce_c_e=[]
    ce_c_u=[]
    rate_nav=[]
    rate_com=[]  
    
    for j in diago:  
        # context effect
        ce_fixed=context_effect(sims[j])
        ce_n_e.append(ce_fixed[:,0])
        ce_n_u.append(ce_fixed[:,1])
        ce_c_e.append(ce_fixed[:,2])
        ce_c_u.append(ce_fixed[:,3])
        # deviance detection
        dd_nav_fixed=(ce_fixed[:,1]-ce_fixed[:,0])/2
        dd_com_fixed=(ce_fixed[:,3]-ce_fixed[:,2])/2
        dd_nav.append(dd_nav_fixed)
        dd_com.append(dd_com_fixed)
        # discriminability index 
        di_fixed=discriminability_index(sims[j])
        di_nav.append(di_fixed[:,0])
        di_iso.append(di_fixed[:,1])
        di_com.append(di_fixed[:,2])
        # firing rate context 
        sp_fixed=sims[j]
        sp_nav=np.reshape(sp_fixed[:,8],(50,20))  
        rate_nav.append(np.mean(sp_nav,1))  #/dur_ech for firing rate 
        sp_com=np.reshape(sp_fixed[:,11],(50,20))    
        rate_com.append(np.mean(sp_com,1))  #/dur_com
        
        
    ## plotting a deviance detection (d.d) variation    
    fig, ax = plt.subplots(1,figsize=(3,2))
    
    A = dd_nav
    B = dd_com
    X = np.linspace(0,1,n_combis)
    
    ax.plot(X,np.mean(A,1),color='b',lw=0.75)
    ax.fill_between(X, np.mean(A,1)-sem(A,1), np.mean(A,1)+sem(A,1), alpha=0.3, facecolor='b')
    ax.plot(X,np.mean(B,1),color='r',lw=0.75)
    ax.fill_between(X, np.mean(B,1)-sem(B,1), np.mean(B,1)+sem(B,1), alpha=0.3, facecolor='r')
    
    ax.axhline(0,linewidth=1, color='k',ls='--')
    ax.set_ylabel('d.d.')
    ax.set_xlabel('anesthesia effect')
    
    # add errorbar for experimental data
    mat_contents = sio.loadmat(r'E:\Users\User\Documents\model-context\ei_data-an.mat')
    a=mat_contents['ei']
    exp_n_e = a[:,0]
    exp_n_u = a[:,1]
    exp_c_e = a[:,2]
    exp_c_u = a[:,3]
    exp_dd_navAN=(exp_n_u-exp_n_e)/2
    exp_dd_comAN=(exp_c_u-exp_c_e)/2
    mat_contents = sio.loadmat(r'E:\Users\User\Documents\model-context\ei_data.mat')
    a=mat_contents['ei']
    exp_n_e = a[:,0]
    exp_n_u = a[:,1]
    exp_c_e = a[:,2]
    exp_c_u = a[:,3]
    exp_dd_navAW=(exp_n_u-exp_n_e)/2
    exp_dd_comAW=(exp_c_u-exp_c_e)/2
#    ax[0].errorbar([0,1],[np.mean(exp_dd_navAW),np.mean(exp_dd_navAN)], yerr=[sem(exp_dd_navAW),sem(exp_dd_navAN)], fmt='o',color='k', mfc='w',ms=4)     
#    ax[1].errorbar([0,1],[np.mean(exp_dd_comAW),np.mean(exp_dd_comAN)], yerr=[sem(exp_dd_comAW),sem(exp_dd_comAN)], fmt='o',color='k', mfc='w',ms=4)

    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.spines['bottom'].set_linewidth(0.75)
    ax.spines['left'].set_linewidth(0.75)
        
    fig.tight_layout() 
    
    fig6=assymetries_plot(X,A,B,0,'b','r',20)
    
    ## plotting a discriminability index (d.i) variation 
    fig2, ax = plt.subplots(1,figsize=(3,2))
    
    A = di_nav
    B = di_com
    X = np.linspace(0,1,n_combis)
    
    ax.plot(X,np.mean(A,1),color='b',lw=0.75)
    ax.fill_between(X, np.mean(A,1)-sem(A,1), np.mean(A,1)+sem(A,1), alpha=0.3, facecolor='b')
    ax.plot(X,np.mean(B,1),color='r',lw=0.75)
    ax.fill_between(X, np.mean(B,1)-sem(B,1), np.mean(B,1)+sem(B,1), alpha=0.3, facecolor='r')
    
    ax.set_ylabel('d.i.')
    ax.set_xlabel('anesthesia effect')
  
    # add errorbar for experimental data
    mat_contents = sio.loadmat(r'E:\Users\User\Documents\model-context\cliff_data-an.mat')
    a=mat_contents['cliff_data']
#    exp_di_isoAN = a[:,0]
    exp_di_navAN = a[:,1]
    exp_di_comAN = a[:,2]
    mat_contents = sio.loadmat(r'E:\Users\User\Documents\model-context\cliff_data.mat')
    a=mat_contents['cliff_data']
#    exp_di_isoAW = a[:,0]
    exp_di_navAW = a[:,1]
    exp_di_comAW = a[:,2]
#    ax[0].errorbar([0,1],[np.mean(exp_di_navAW),np.mean(exp_di_navAN)], yerr=[np.std(exp_di_navAW),np.std(exp_di_navAN)], fmt='o',color='k', mfc='w',ms=4)       
#    ax[1].errorbar([0,1],[np.mean(exp_di_comAW),np.mean(exp_di_comAN)], yerr=[np.std(exp_di_comAW),np.std(exp_di_comAN)], fmt='o',color='k', mfc='w',ms=4)       

    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.spines['bottom'].set_linewidth(0.75)
    ax.spines['left'].set_linewidth(0.75)
        
    fig2.tight_layout() 
    
    fig7=assymetries_plot(X,A,B,0,'b','r',20)
    
    ## plotting a context effect (c.e) variation    
    fig3, ax = plt.subplots(1,2,figsize=(6,2))
    
    A = ce_n_e
    B = ce_n_u
    X = np.linspace(0,1,n_combis)
    
    ax[0].plot(X,np.mean(A,1),color='gray',lw=0.75)
    ax[0].plot(X,np.mean(A,1),'.',markersize=3,color='gray')
    ax[0].fill_between(X, np.mean(A,1)-sem(A,1), np.mean(A,1)+sem(A,1), alpha=0.3, facecolor='gray')
    ax[0].plot(X,np.mean(B,1),color='k',lw=0.75)
    ax[0].plot(X,np.mean(B,1),'.',markersize=3,color='k')
    ax[0].fill_between(X, np.mean(B,1)-sem(B,1), np.mean(B,1)+sem(B,1), alpha=0.3, facecolor='k')

    fig7=assymetries_plot(X,A,B,0,'gray','k',20)
    
    A = ce_c_e
    B = ce_c_u
    X = np.linspace(0,1,n_combis)
    
    ax[1].plot(X,np.mean(A,1),color='gray',lw=0.75)
    ax[1].plot(X,np.mean(A,1),'.',markersize=3,color='gray')
    ax[1].fill_between(X, np.mean(A,1)-sem(A,1), np.mean(A,1)+sem(A,1), alpha=0.3, facecolor='gray')
    ax[1].plot(X,np.mean(B,1),color='k',lw=0.75)
    ax[1].plot(X,np.mean(B,1),'.',markersize=3,color='k')
    ax[1].fill_between(X, np.mean(B,1)-sem(B,1), np.mean(B,1)+sem(B,1), alpha=0.3, facecolor='k')
    
    ax[0].set_ylabel('c.e.')
    ax[0].set_xlabel('anesthesia effect')
    ax[1].set_xlabel('anesthesia effect')
#    ax[0].axhline(0,linewidth=1, color='k',ls='--')
#    ax[1].axhline(0,linewidth=1, color='k',ls='--')
#    ax[1,1].set_ylim([-0.85, 0])
    
    # add errorbar for experimental data
    mat_contents = sio.loadmat(r'E:\Users\User\Documents\model-context\ei_data-an.mat')
    a=mat_contents['ei']
    exp_n_eAN = a[:,0]
    exp_n_uAN = a[:,1]
    exp_c_eAN = a[:,2]
    exp_c_uAN = a[:,3]
    mat_contents = sio.loadmat(r'E:\Users\User\Documents\model-context\ei_data.mat')
    a=mat_contents['ei']
    exp_n_eAW = a[:,0]
    exp_n_uAW = a[:,1]
    exp_c_eAW = a[:,2]
    exp_c_uAW = a[:,3]
#    ax[0,0].errorbar([0,1],[np.mean(exp_n_eAW),np.mean(exp_n_eAN)], yerr=[np.std(exp_n_eAW),np.std(exp_n_eAN)], fmt='o',color='k', mfc='w',ms=4)       
#    ax[1,0].errorbar([0,1],[np.mean(exp_n_uAW),np.mean(exp_n_uAN)], yerr=[np.std(exp_n_uAW),np.std(exp_n_uAN)], fmt='o',color='k', mfc='w',ms=4)       
#    ax[0,1].errorbar([0,1],[np.mean(exp_c_eAW),np.mean(exp_c_eAN)], yerr=[np.std(exp_c_eAW),np.std(exp_c_eAN)], fmt='o',color='k', mfc='w',ms=4)       
#    ax[1,1].errorbar([0,1],[np.mean(exp_c_uAW),np.mean(exp_c_uAN)], yerr=[np.std(exp_c_uAW),np.std(exp_c_uAN)], fmt='o',color='k', mfc='w',ms=4)       

    ax = ax.flatten() 
    for axi in ax:
        axi.spines['top'].set_visible(False)
        axi.spines['right'].set_visible(False)
        axi.spines['bottom'].set_linewidth(0.75)
        axi.spines['left'].set_linewidth(0.75)   
        
    fig3.tight_layout() 
    
    fig8=assymetries_plot(X,A,B,0,'gray','k',20)
    
    ## plotting a firing rate context(rate) variation    
    fig4, ax = plt.subplots(1,figsize=(3,2))
    
    A = rate_nav
    B = rate_com
    X = np.linspace(0,1,n_combis)
    
    ax.plot(X,np.mean(A,1),color='b',lw=0.75)
    ax.plot(X,np.mean(A,1),'.',markersize=3,color='b')
    ax.fill_between(X, np.mean(A,1)-np.std(A,1), np.mean(A,1)+np.std(A,1), alpha=0.3, facecolor='b')
    ax.plot(X,np.mean(B,1),color='r',lw=0.75)
    ax.plot(X,np.mean(B,1),'.',markersize=3,color='r')
    ax.fill_between(X, np.mean(B,1)-np.std(B,1), np.mean(B,1)+np.std(B,1), alpha=0.3, facecolor='r')
    
    ax.set_ylabel('spike count during context')
    ax.set_xlabel('anesthesia effect')
    
    # add errorbar for experimental data
    mat_contents = sio.loadmat(r'E:\Users\User\Documents\model-context\rate_context-an.mat')
    a=mat_contents['rate']
    exp_rate_navAN = a[:,0]
    exp_rate_comAN = a[:,1]
    mat_contents = sio.loadmat(r'E:\Users\User\Documents\model-context\rate_context.mat')
    a=mat_contents['rate']
    exp_rate_navAW = a[:,0]
    exp_rate_comAW = a[:,1]
    
#    ax[0].errorbar([0,1],[np.mean(exp_rate_navAW),np.mean(exp_rate_navAN)], yerr=[np.std(exp_rate_navAW),np.std(exp_rate_navAN)], fmt='o',color='k', mfc='w',ms=4)       
#    ax[1].errorbar([0,1],[np.mean(exp_rate_comAW),np.mean(exp_rate_comAN)], yerr=[np.std(exp_rate_comAW),np.std(exp_rate_comAN)], fmt='o',color='k', mfc='w',ms=4)       

    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.spines['bottom'].set_linewidth(0.75)
    ax.spines['left'].set_linewidth(0.75)
        
    fig4.tight_layout()   
    
    fig9=assymetries_plot(X,A,B,0,'b','r',45)
    
        
    return fig3, fig4  

def dif_2d(sims): 

    # for 2D showing differences between com and nav 
    parameters1=sims[0]
    parameters2=sims[1]
    sims = np.delete(sims,[0,1], 0)
    nx=len(parameters1)
    ny=len(parameters2)
    or_par1=1.0000000000000004
    or_par2=1#7.9999999999999964
    ind_par1=np.where(parameters1==or_par1)[0]
    ind_par2=np.where(parameters2==or_par2)[0]
#    dur_ech=1.3311
#    dur_com=1.9266
    
    # for plotting 2D difference between context averages 
    dd_nav=[]
    dd_com=[]
    di_nav=[]
    di_iso=[]
    di_com=[]
    ce_n_e=[]
    ce_n_u=[]
    ce_c_e=[]
    ce_c_u=[]
    rate_nav=[]
    rate_com=[]  

    for j in range(nx*ny):  
        # j is a particular combinations of parameters
        # 2D only needs the average
        ce_i=np.mean(context_effect(sims[j]),0) 
        dd_nav_i=(ce_i[1]-ce_i[0])/2
        dd_com_i=(ce_i[3]-ce_i[2])/2
        dd_nav.append(dd_nav_i)
        dd_com.append(dd_com_i)
        # context effect 
        ce_n_e.append(ce_i[0])
        ce_n_u.append(ce_i[1])
        ce_c_e.append(ce_i[2])
        ce_c_u.append(ce_i[3])
        # discriminability index 
        di_i=np.mean(discriminability_index(sims[j]),0) 
        di_nav.append(di_i[0])
        di_iso.append(di_i[1])
        di_com.append(di_i[2])
        # firing rate 
        sp_contexts=sims[j]
        rate_nav.append(np.mean(sp_contexts[:,8]))  #/dur_ech
        rate_com.append(np.mean(sp_contexts[:,11]))  #/dur_com
 
    dd_nav=np.array(dd_nav).reshape(nx,ny)
    dd_com=np.array(dd_com).reshape(nx,ny)
    di_nav=np.array(di_nav).reshape(nx,ny)
    di_iso=np.array(di_iso).reshape(nx,ny)
    di_com=np.array(di_com).reshape(nx,ny)
    ce_n_e=np.array(ce_n_e).reshape(nx,ny)
    ce_n_u=np.array(ce_n_u).reshape(nx,ny)
    ce_c_e=np.array(ce_c_e).reshape(nx,ny)
    ce_c_u=np.array(ce_c_u).reshape(nx,ny)
    rate_nav=np.array(rate_nav).reshape(nx,ny)
    rate_com=np.array(rate_com).reshape(nx,ny)
    

    ## echolocation

    # difference with awake state
    dd=dd_nav-dd_nav[ind_par1,ind_par2]
    ce_e=ce_n_e-ce_n_e[ind_par1,ind_par2]
    ce_u=ce_n_u-ce_n_u[ind_par1,ind_par2]
    rate=rate_nav-rate_nav[ind_par1,ind_par2]
    
    # smoothing data
    x=parameters2
    y=parameters1
    f_dd = interp2d(x,y,dd, kind='cubic')
    f_ce_e = interp2d(x,y,ce_e, kind='cubic')
    f_ce_u = interp2d(x,y,ce_u, kind='cubic')
    f_rate = interp2d(x,y,rate, kind='cubic')

    # increasing resolution axes
    alpha=5   # smoothing factor
    
    x_in=parameters2[0]
    x_end=parameters2[-1]
    x_n=len(parameters2)
    y_in=parameters1[0]
    y_end=parameters1[-1]
    y_n=len(parameters1)
    x_nn=x_n + (x_n-1)*(alpha-1)
    y_nn=y_n + (y_n-1)*(alpha-1)
    xnew = np.linspace(x_in,x_end,x_nn)
    ynew = np.linspace(y_in,y_end,y_nn)
    
    Xn, Yn = np.meshgrid(xnew, ynew)
    # generate smooth data for the new axes
    dd_smooth = f_dd(xnew,ynew)
    ce_e_smooth = f_ce_e(xnew,ynew)
    ce_u_smooth = f_ce_u(xnew,ynew)
    rate_smooth = f_rate(xnew,ynew)


    fig, ax = plt.subplots(1,2,sharex=True,sharey=True,figsize=(2.1,1))
    mm=0.2
    # c.e. match
    q1=ax[0].pcolor(Xn,Yn,ce_e_smooth,cmap='RdBu',vmin=-mm,vmax=mm)
#    ax[0].set_title('c.e. match')
    # c.e. mismatch
    q2=ax[1].pcolor(Xn,Yn,ce_u_smooth,cmap='RdBu',vmin=-mm,vmax=mm)
#    ax[1].set_title('c.e. mismatch')
  
    # adding colorbars
    divider = make_axes_locatable(ax[0])
    cax = divider.append_axes("right", size="5%", pad=0.075)
    fig.colorbar(q1,cax=cax)
    cax.remove() 
    divider = make_axes_locatable(ax[1])
    cax = divider.append_axes("right", size="5%", pad=0.075)
    fig.colorbar(q2,cax=cax)
    
    fig.tight_layout() 
    
    fig2, ax2 = plt.subplots(1,figsize=(1.15,1))   
    #  s.s.s    
    q1=ax2.pcolormesh(Xn,Yn,dd_smooth,cmap='RdBu',vmin=-mm,vmax=mm)
#    ax2[0].set_title('s.s.s.')
    divider = make_axes_locatable(ax2)
    cax = divider.append_axes("right", size="5%", pad=0.075)
    fig2.colorbar(q1,cax=cax)
    
    fig2.tight_layout() 
    
    fig3, ax3 = plt.subplots(1,figsize=(1.15,1))  
    # firing rate 
    mm=30 
    q1=ax3.pcolor(Xn,Yn,rate_smooth,cmap='RdBu',vmin=-mm,vmax=mm)
    ax3.contour(Xn,Yn,rate_smooth,levels = [-23],origin='lower',colors=['k'],linewidths=0.75)
#    ax2[3].set_title('firing rate')
    divider = make_axes_locatable(ax3)
    cax = divider.append_axes("right", size="5%", pad=0.075)
    fig3.colorbar(q1,cax=cax)
    
    fig3.tight_layout() 
    

    ## communication

    # difference with awake state
    dd=dd_com-dd_com[ind_par1,ind_par2]
    ce_e=ce_c_e-ce_c_e[ind_par1,ind_par2]
    ce_u=ce_c_u-ce_c_u[ind_par1,ind_par2]
    rate=rate_com-rate_com[ind_par1,ind_par2]
    
    #smoothing data
    f_dd = interp2d(x,y,dd, kind='cubic')
    f_ce_e = interp2d(x,y,ce_e, kind='cubic')
    f_ce_u = interp2d(x,y,ce_u, kind='cubic')
    f_rate = interp2d(x,y,rate, kind='cubic')
    
    dd_smooth = f_dd(xnew,ynew)
    ce_e_smooth = f_ce_e(xnew,ynew)
    ce_u_smooth = f_ce_u(xnew,ynew)
    rate_smooth = f_rate(xnew,ynew)

    fig4, ax4 = plt.subplots(1,2,sharex=True,sharey=True,figsize=(2.1,1))
    mm=0.2
    # c.e. match 
    q1=ax4[0].pcolor(Xn,Yn,ce_e_smooth,cmap='RdBu',vmin=-mm,vmax=mm)
#    ax4[0].set_title('c.e. match')
    # c.e. mismatch
    q2=ax4[1].pcolor(Xn,Yn,ce_u_smooth,cmap='RdBu',vmin=-mm,vmax=mm)
#    ax4[1].set_title('c.e. mismatch')
    divider = make_axes_locatable(ax4[0])
    cax = divider.append_axes("right", size="5%", pad=0.075)
    fig4.colorbar(q1,cax=cax)
    cax.remove() 
    divider = make_axes_locatable(ax4[1])
    cax = divider.append_axes("right", size="5%", pad=0.075)
    fig4.colorbar(q2,cax=cax)
    
    fig4.tight_layout() 
    
    fig5, ax5 = plt.subplots(1,figsize=(1.15,1))  
    # s.s.s    
    q1=ax5.pcolor(Xn,Yn,dd_smooth,cmap='RdBu',vmin=-mm,vmax=mm)
#    ax5.set_title('s.s.s.')
    divider = make_axes_locatable(ax5)
    cax = divider.append_axes("right", size="5%", pad=0.075)
    fig5.colorbar(q1,cax=cax)  
    
    fig5.tight_layout() 
    
    fig6, ax6 = plt.subplots(1,figsize=(1.15,1))  
    # firing rate
    mm=30 
    q1=ax6.pcolor(Xn,Yn,rate_smooth,cmap='RdBu',vmin=-mm,vmax=mm)
    ax6.contour(Xn,Yn,rate_smooth,levels = [-16],origin='lower',colors=['k'],linewidths=0.75)
#    ax6.set_title('firing rate')
    divider = make_axes_locatable(ax6)
    cax = divider.append_axes("right", size="5%", pad=0.075)
    fig6.colorbar(q1,cax=cax)

    fig6.tight_layout() 

    # plotting for difference between communication and echolocation
    jump=3
    fig7, ax7 = plt.subplots(2,2,figsize=(6,4))
    fig7.canvas.manager.set_window_title('Difference between contexts')
    # difference with awake state
    # communication
    dd2=dd_com-dd_com[ind_par1,ind_par2]
    ce_e2=ce_c_e-ce_c_e[ind_par1,ind_par2]
    ce_u2=ce_c_u-ce_c_u[ind_par1,ind_par2]
    rate2=rate_com-rate_com[ind_par1,ind_par2]
    # echolocation
    dd1=dd_nav-dd_nav[ind_par1,ind_par2]
    ce_e1=ce_n_e-ce_n_e[ind_par1,ind_par2]
    ce_u1=ce_n_u-ce_n_u[ind_par1,ind_par2]
    rate1=rate_nav-rate_nav[ind_par1,ind_par2]
    # difference between communication and echolocation
    dd=dd2-dd1
    ce_e=ce_e2-ce_e1
    ce_u=ce_u2-ce_u1
    rate=rate2-rate1
    # first plot: s.s.s  
    mm=0.2 
    q1=ax7[0,0].imshow(dd,cmap='RdBu',origin='lower',vmin=-mm,vmax=mm,interpolation='hamming')
    ax7[0,0].set_xticks(np.arange(0,ny,jump))
    ax7[0,0].set_xticklabels(['{:.1f}'.format(x) for x in parameters2[::jump]])
    ax7[0,0].set_yticks(np.arange(0,nx,jump))
    ax7[0,0].set_yticklabels(['{:.1f}'.format(x) for x in parameters1[::jump]])
    ax7[0,0].set_title('s.s.s.')
    # third plot: c.e. mismatch
    q3=ax7[1,0].imshow(ce_e,cmap='RdBu',origin='lower',vmin=-mm,vmax=mm,interpolation='hamming')
    ax7[1,0].set_xticks(np.arange(0,ny,jump))
    ax7[1,0].set_xticklabels(['{:.1f}'.format(x) for x in parameters2[::jump]])
    ax7[1,0].set_yticks(np.arange(0,nx,jump))
    ax7[1,0].set_yticklabels(['{:.1f}'.format(x) for x in parameters1[::jump]])
    ax7[1,0].set_title('c.e. match')
    # forth plot: firing rate
    q4=ax7[1,1].imshow(ce_u,cmap='RdBu',origin='lower',vmin=-mm,vmax=mm,interpolation='hamming')
    ax7[1,1].set_xticks(np.arange(0,ny,jump))
    ax7[1,1].set_xticklabels(['{:.1f}'.format(x) for x in parameters2[::jump]])
    ax7[1,1].set_yticks(np.arange(0,nx,jump))
    ax7[1,1].set_yticklabels(['{:.1f}'.format(x) for x in parameters1[::jump]])
    ax7[1,1].set_title('c.e. mismatch')
    # second plot: c.e. match 
    mm=30
    q2=ax7[0,1].imshow(rate,cmap='RdBu',origin='lower',vmin=-mm,vmax=mm,interpolation='hamming')
    ax7[0,1].set_xticks(np.arange(0,ny,jump))
    ax7[0,1].set_xticklabels(['{:.1f}'.format(x) for x in parameters2[::jump]])
    ax7[0,1].set_yticks(np.arange(0,nx,jump))
    ax7[0,1].set_yticklabels(['{:.1f}'.format(x) for x in parameters1[::jump]])
    ax7[0,1].set_title('firing context')
    
    # add colorbars
    fig7.colorbar(q1,ax=ax7[0,0])
    fig7.colorbar(q2,ax=ax7[0,1])
    fig7.colorbar(q3,ax=ax7[1,0])
    fig7.colorbar(q4,ax=ax7[1,1])
    
    fig7.tight_layout()     
    
    
    return fig,ax,fig2,ax2,fig3,ax3,fig4,ax4,fig5,ax5,fig6,ax6

def adapt_supp(N_units):
    delay=1000
    trials=20*N_units
    dt_=0.1
 
    # SIMULATION NAVIGATION CONTEXT AWAKE
    tau_th=550
#    d_LF=0.045 
    d_HF=0.040
    coefC=1
    
    context=1
    exp=1
    exec(open("./input.py").read())  
    exec(open("./model.py").read()) 
    nav_postsyn=pos.V_th*1000   # to mV 
    nav_presyn=mon.x_S
    nav_high=nav_presyn[trials:]
    nav_low=nav_presyn[:trials]
    nav_tot=np.shape(nav_presyn)[1]*dt_
    nav_time=np.arange(0,nav_tot,dt_)
    # average acros trials
    nav_post=np.mean(nav_postsyn,0)
    nav_pre_high=np.mean(nav_high,0)
    nav_pre_low=np.mean(nav_low,0)
    
    # SIMULATION COMMUNICATION CONTEXT AWAKE
    context=2
    exec(open("./input.py").read())  
    exec(open("./model.py").read()) 
    com_postsyn=pos.V_th*1000
    com_presyn=mon.x_S
    com_high=com_presyn[trials:]
    com_low=com_presyn[:trials]
    com_tot=np.shape(com_presyn)[1]*dt_
    com_time=np.arange(0,com_tot,dt_)
    # average across trials
    com_post=np.mean(com_postsyn,0)
    com_pre_high=np.mean(com_high,0)
    com_pre_low=np.mean(com_low,0)
    
    
    # SIMULATION NAVIGATION CONTEXT ANEST
    tau_th=2.0*550
#    d_LF=0.045 
    d_HF=1.8*0.040
    coefC=0.6
    
    context=1
    exp=1
    exec(open("./input.py").read())  
    exec(open("./model.py").read()) 
    nav_postsyn=pos.V_th*1000   # to mV 
    nav_presyn=mon.x_S
    nav_high=nav_presyn[trials:]
    nav_low=nav_presyn[:trials]
    nav_tot=np.shape(nav_presyn)[1]*dt_
    nav_time=np.arange(0,nav_tot,dt_)
    # average acros trials
    nav_post2=np.mean(nav_postsyn,0)
    nav_pre_high2=np.mean(nav_high,0)
    nav_pre_low2=np.mean(nav_low,0)
    
    # SIMULATION COMMUNICATION CONTEXT ANEST
    context=2
    exec(open("./input.py").read())  
    exec(open("./model.py").read()) 
    com_postsyn=pos.V_th*1000
    com_presyn=mon.x_S
    com_high=com_presyn[trials:]
    com_low=com_presyn[:trials]
    com_tot=np.shape(com_presyn)[1]*dt_
    com_time=np.arange(0,com_tot,dt_)
    # average across trials
    com_post2=np.mean(com_postsyn,0)
    com_pre_high2=np.mean(com_high,0)
    com_pre_low2=np.mean(com_low,0)
    
    

    ## plotting
    vm_i=-50
    vm_s=-20
    xs_i=-1.1
    xs_s=-0.0 

    fig, ax = plt.subplots(2,2,sharex=True,sharey='row',figsize=(3.54,2))
    
    # navigation postsyn
    ax[0,0].plot(nav_time,nav_post,c='k',lw=0.75)
    ax[0,0].plot(nav_time,nav_post2,c='k',lw=0.75,ls='--')
#    ax[0,0].vlines(onset[1],vm_i,vm_s,lw=0.75)
    ax[0,0].vlines(onset[1]-delay+60,vm_i,vm_s,lw=0.75)
    ax[0,0].set_ylim(vm_i,vm_s)
    ax[0,0].set_xlim(0,onset[1]+50)
    ax[0,0].set_ylabel('spiking\n threshold')
    
    # communication postsyn
    ax[0,1].plot(com_time,com_post,c='k',lw=0.75)
    ax[0,1].plot(com_time,com_post2,c='k',lw=0.75,ls='--')
#    ax[0,1].vlines(onset[2],vm_i,vm_s,lw=0.75)
    ax[0,1].vlines(onset[2]-delay+60,vm_i,vm_s,lw=0.75)
    ax[0,1].set_ylim(vm_i,vm_s)
    ax[0,1].set_xlim(0,onset[2]+50)
    
    # navigation presyn
    ax[1,0].plot(nav_time,-1*nav_pre_high,c='b',lw=0.75)
    ax[1,0].plot(nav_time,-1*nav_pre_high2,c='b',lw=0.75,ls='--')
    ax[1,0].plot(nav_time,-1*nav_pre_low,c='r',lw=0.75)
    ax[1,0].plot(nav_time,-1*nav_pre_low2,c='r',lw=0.75,ls='--')
#    ax[1,0].vlines(onset[1],xs_i,xs_s,lw=0.75)
    ax[1,0].vlines(onset[1]-delay+60,xs_i,xs_s,lw=0.75)
    ax[1,0].set_ylim(xs_i,xs_s)
    ax[1,0].set_xlim(0,onset[1]+50)
    ax[1,0].set_ylabel('strength\n synapse')
    
    # communication presyn
    ax[1,1].plot(com_time,-1*com_pre_high,c='b',lw=0.75)
    ax[1,1].plot(com_time,-1*com_pre_high2,c='b',lw=0.75,ls='--')
    ax[1,1].plot(com_time,-1*com_pre_low,c='r',lw=0.75)
    ax[1,1].plot(com_time,-1*com_pre_low2,c='r',lw=0.75,ls='--')
#    ax[1,1].vlines(onset[2],xs_i,xs_s,lw=0.75)
    ax[1,1].vlines(onset[2]-delay+60,xs_i,xs_s,lw=0.75)
    ax[1,1].set_ylim(xs_i,xs_s)
    ax[1,1].set_xlim(0,onset[2]+50)
    
    
    ax = ax.flatten() 
    for axi in ax:
        axi.spines['top'].set_visible(False)
        axi.spines['right'].set_visible(False)
        axi.spines['bottom'].set_linewidth(0.75)
        axi.spines['left'].set_linewidth(0.75)
   
    plt.figtext(0.5, 0.01, 'Time (ms)', rotation=0, verticalalignment='bottom')
    
    fig.tight_layout() 
    
    return fig 

def assymetries_plot(x,nav_data,com_data,original_pos,c1,c2,maxY):
    # comparison between distributions across several parameters and original awake distribution
    # original awake distribution
    original_nav=nav_data[original_pos]
    original_com=com_data[original_pos]
    # context: echolocation
    p_navs=[]
    for i in range(len(nav_data)): # across several parameters
#        h,p=ranksums(original_nav,nav_data[i])
        p,size = cliffsDelta(original_nav,nav_data[i]) 
        p_navs.append(abs(p))
    # context: communication
    p_coms=[]
    for i in range(len(com_data)):
#        h,p=ranksums(original_com,com_data[i])
        p,size = cliffsDelta(original_com,com_data[i]) 
        p_coms.append(abs(p))     
        
    # plotting     
    fig, ax = plt.subplots(1,1,figsize=(0.75,1))   # for insets: (1.3,0.75)

    # plot pvalues
#    ax.plot(x,-np.log(p_navs),linestyle="--",lw=0.5,marker="o",markersize=2,color=c1)
    ax.plot(x,p_navs,linestyle="-",lw=0.75,marker="o",markersize=2,color=c1)
#    ax.plot(x,-np.log(p_coms),linestyle="--",lw=0.5,marker="o",markersize=2,color=c2)
    ax.plot(x,p_coms,linestyle="-",lw=0.75,marker="o",markersize=2,color=c2)
    
    # plot significance level
#    minY=-0.5
#    maxY=40
#    ax.set_ylim(ax.get_ylim()[0],maxY)
    ax.set_ylim(-0.1,1)
#    ax.axhspan(ax.get_ylim()[0],-np.log(0.05),facecolor='0.95')
#    ax.axhspan(-np.log(0.05),ax.get_ylim()[1],facecolor='0.85')
    # plot cliff delta interpretation as background color
    neg=0.147
    small=0.33
    med=0.47
    larg=1
    upper_color='black'
    down_color='orange'
#    ax.axhline(0,ls='--',c='k',lw=0.75)
#    ax.axhspan(-neg,0,alpha=0.05,color=down_color)
    ax.axhspan(0,neg,alpha=0.05,color=upper_color)
#    ax.axhspan(-small,-neg,alpha=0.1,color=down_color)
    ax.axhspan(neg,small,alpha=0.1,color=upper_color)
#    ax.axhspan(-med,-small,alpha=0.2,color=down_color)
    ax.axhspan(small,med,alpha=0.2,color=upper_color)
#    ax.axhspan(-larg,-med,alpha=0.3,color=down_color)
    ax.axhspan(med,larg,alpha=0.3,color=upper_color)


#    ax.set_ylabel('effect size') #('-log(p-value)')
#    ax.set_xlabel(r'$\tau$ adaptation')
#    ax.set_xticklabels([])
#    ax.set_yticks([])
#    ax.invert_xaxis()
        
    ax.axis('off')
#    ax.set_xlim(x[0],x[-1])
#    ax.set_ylim(minY,maxY)

    fig.tight_layout()
    
    return fig
        
def ns_region(sims,fig1,ax1,fig2,ax2,fig3,ax3,fig4,ax4,fig5,ax5,fig6,ax6,sims_original=None):
    
    # for 2D non significant regions of anesthesia effect  
    parameters1=sims[0]
    parameters2=sims[1]
    sims = np.delete(sims,[0,1], 0)
    nx=len(parameters1)
    ny=len(parameters2)
    # original parameters (no anesthesia)
    or_par1=1.0000000000000004
    or_par2=1
    ind_par1=np.where(parameters1==or_par1)[0]
    ind_par2=np.where(parameters2==or_par2)[0]
    original_pos=int(ind_par1)*(ny)+int(ind_par2)
    # original distribution of output
    if sims_original is None:
        o_sims=sims[original_pos]
    else:
        sims_original = np.delete(sims_original,[0,1], 0)
        o_sims=sims_original[original_pos]
    # original context effect
    o_ce=context_effect(o_sims)
    o_ce_n_e=o_ce[:,0]
    o_ce_n_u=o_ce[:,1]
    o_ce_c_e=o_ce[:,2]
    o_ce_c_u=o_ce[:,3]
    # original deviance detection
    o_dd_nav=(o_ce_n_u-o_ce_n_e)/2
    o_dd_com=(o_ce_c_u-o_ce_c_e)/2
    # original spiking rate during context
    o_rate_n=o_sims[:,8]
    o_rate_c=o_sims[:,11]
    
    # duration nof context
#    dur_ech=1.3311
#    dur_com=1.9266
    
    # for plotting 2D non-significant regions between anesthesia and non-anesthesia
    ce_n_e=[]
    ce_n_u=[]
    ce_c_e=[]
    ce_c_u=[]
    dd_nav=[]
    dd_com=[]
    rate_nav=[]
    rate_com=[]  
                
    for j in range(nx*ny):  
        # j is a particular combinations of parameters
        j_sims=sims[j]
        # original context effect
        j_ce=context_effect(j_sims)
        j_ce_n_e=j_ce[:,0]
        j_ce_n_u=j_ce[:,1]
        j_ce_c_e=j_ce[:,2]
        j_ce_c_u=j_ce[:,3]
        # original deviance detection
        j_dd_nav=(j_ce_n_u-j_ce_n_e)/2
        j_dd_com=(j_ce_c_u-j_ce_c_e)/2
        # original spiking rate during context
        j_rate_n=j_sims[:,8]
        j_rate_c=j_sims[:,11]
        
        # comparison between anesthesia and no-anesthesia - p-value
#        h,p=ranksums(j_ce_n_e,o_ce_n_e)
#        ce_n_e.append(1)if p<0.05 else ce_n_e.append(0) 
#        h,p=ranksums(j_ce_n_u,o_ce_n_u)
#        ce_n_u.append(1)if p<0.05 else ce_n_u.append(0) 
#        h,p=ranksums(j_ce_c_e,o_ce_c_e)
#        ce_c_e.append(1)if p<0.05 else ce_c_e.append(0) 
#        h,p=ranksums(j_ce_c_u,o_ce_c_u)
#        ce_c_u.append(1) if p<0.05 else ce_c_u.append(0) 
#        
#        h,p=ranksums(j_dd_nav,o_dd_nav)
#        dd_nav.append(1) if p<0.05 else dd_nav.append(0)           
#        h,p=ranksums(j_dd_com,o_dd_com)
#        dd_com.append(1) if p<0.05 else dd_com.append(0) 
#        
#        h,p=ranksums(j_rate_n,o_rate_n)
#        rate_nav.append(1) if p<0.05 else rate_nav.append(0)
#        h,p=ranksums(j_rate_c,o_rate_c)
#        rate_com.append(1) if p<0.05 else rate_com.append(0) 
        
        # comparison between anesthesia and no-anesthesia - effect size 
        h,p=cliffsDelta(j_ce_n_e,o_ce_n_e)
        ce_n_e.append(h)
        h,p=cliffsDelta(j_ce_n_u,o_ce_n_u)
        ce_n_u.append(h)
        h,p=cliffsDelta(j_ce_c_e,o_ce_c_e)
        ce_c_e.append(h)
        h,p=cliffsDelta(j_ce_c_u,o_ce_c_u)
        ce_c_u.append(h) 
        
        h,p=cliffsDelta(j_dd_nav,o_dd_nav)
        dd_nav.append(h)
        h,p=cliffsDelta(j_dd_com,o_dd_com)
        dd_com.append(h) 
        
        h,p=cliffsDelta(j_rate_n,o_rate_n)
        rate_nav.append(h)
        h,p=cliffsDelta(j_rate_c,o_rate_c)
        rate_com.append(h) 
        
    ce_n_e=np.array(ce_n_e).reshape(nx,ny)
    ce_n_u=np.array(ce_n_u).reshape(nx,ny)
    ce_c_e=np.array(ce_c_e).reshape(nx,ny)
    ce_c_u=np.array(ce_c_u).reshape(nx,ny)
    dd_nav=np.array(dd_nav).reshape(nx,ny)
    dd_com=np.array(dd_com).reshape(nx,ny)
    rate_nav=np.array(rate_nav).reshape(nx,ny)
    rate_com=np.array(rate_com).reshape(nx,ny)
        
    ## plotting
#    fig, ax = plt.subplots(1,4,sharex=True,sharey=True,figsize=(5,1.5))
    Xn, Yn = np.meshgrid(parameters2,parameters1)

    c_='gray'
    ls_=0.75
#    cliff_interpret=[-1,-0.47,-0.33,-0.147,0.147,0.33,0.47,1]
    cliff_interpret=[-0.147,0.147]

    plt.figure(1,figsize=fig1.get_size_inches())
    ax1[0].contour(Xn,Yn,ce_n_e,levels=cliff_interpret,origin='lower',colors=[c_],linewidths=ls_)
    ax1[1].contour(Xn,Yn,ce_n_u,levels=cliff_interpret,origin='lower',colors=[c_],linewidths=ls_)

    plt.figure(2,figsize=fig2.get_size_inches())
    #    ax.imshow(dd_nav,origin='lower')
    ax2.contour(Xn,Yn,dd_nav,levels =cliff_interpret,origin='lower',colors=[c_],linewidths=ls_)
    
#    plt.figure(3,figsize=fig3.get_size_inches())
#    ax3.contour(Xn,Yn,rate_nav,levels=cliff_interpret,origin='lower',colors=[c_],linewidths=ls_)

    plt.figure(4,figsize=fig4.get_size_inches())
    ax4[0].contour(Xn,Yn,ce_c_e,levels=cliff_interpret,origin='lower',colors=[c_],linewidths=ls_)
    ax4[1].contour(Xn,Yn,ce_c_u,levels=cliff_interpret,origin='lower',colors=[c_],linewidths=ls_)

    plt.figure(5,figsize=fig5.get_size_inches())
    ax5.contour(Xn,Yn,dd_com,levels=cliff_interpret,origin='lower',colors=[c_],linewidths=ls_)
    
#    plt.figure(6,figsize=fig6.get_size_inches())
#    ax6.contour(Xn,Yn,rate_com,levels =cliff_interpret,origin='lower',colors=[c_],linewidths=ls_)
    return fig1,fig2,fig3,fig4,fig5,fig6

def effect_tauth(sims, sims_tauth): 

    # distribution of outputs for model awake using sims
    parameters1=sims[0]
    parameters2=sims[1]
    sims = np.delete(sims,[0,1], 0)
    nx=len(parameters1)
    ny=len(parameters2)
    or_par1=1.0000000000000004
    or_par2=1#7.9999999999999964
    ind_par1=np.where(parameters1==or_par1)[0]
    ind_par2=np.where(parameters2==or_par2)[0]
#    dur_ech=1.3311
#    dur_com=1.9266
    
    dd_nav=[]
    dd_com=[]
    di_nav=[]
    di_iso=[]
    di_com=[]
    ce_n_e=[]
    ce_n_u=[]
    ce_c_e=[]
    ce_c_u=[]
    rate_nav=[]
    rate_com=[]  

    for j in range(nx*ny):  
        # j is a particular combinations of parameters
        # 2D only needs the average
        ce_i=np.mean(context_effect(sims[j]),0) 
        dd_nav_i=(ce_i[1]-ce_i[0])/2
        dd_com_i=(ce_i[3]-ce_i[2])/2
        dd_nav.append(dd_nav_i)
        dd_com.append(dd_com_i)
        # context effect 
        ce_n_e.append(ce_i[0])
        ce_n_u.append(ce_i[1])
        ce_c_e.append(ce_i[2])
        ce_c_u.append(ce_i[3])
        # discriminability index 
        di_i=np.mean(discriminability_index(sims[j]),0) 
        di_nav.append(di_i[0])
        di_iso.append(di_i[1])
        di_com.append(di_i[2])
        # firing rate 
        sp_contexts=sims[j]
        rate_nav.append(np.mean(sp_contexts[:,8]))  #/dur_ech
        rate_com.append(np.mean(sp_contexts[:,11]))  #/dur_com
 
    dd_nav=np.array(dd_nav).reshape(nx,ny)
    dd_com=np.array(dd_com).reshape(nx,ny)
    di_nav=np.array(di_nav).reshape(nx,ny)
    di_iso=np.array(di_iso).reshape(nx,ny)
    di_com=np.array(di_com).reshape(nx,ny)
    ce_n_e=np.array(ce_n_e).reshape(nx,ny)
    ce_n_u=np.array(ce_n_u).reshape(nx,ny)
    ce_c_e=np.array(ce_c_e).reshape(nx,ny)
    ce_c_u=np.array(ce_c_u).reshape(nx,ny)
    rate_nav=np.array(rate_nav).reshape(nx,ny)
    rate_com=np.array(rate_com).reshape(nx,ny)
    
    ## original awake distributions 
    dd_nav_o=dd_nav[ind_par1,ind_par2]
    ce_n_e_o=ce_n_e[ind_par1,ind_par2]
    ce_n_u_o=ce_n_u[ind_par1,ind_par2]
    rate_nav_o=rate_nav[ind_par1,ind_par2]
    
    dd_com_o=dd_com[ind_par1,ind_par2]
    ce_c_e_o=ce_c_e[ind_par1,ind_par2]
    ce_c_u_o=ce_c_u[ind_par1,ind_par2]
    rate_com_o=rate_com[ind_par1,ind_par2]
    

    
    # distribution of outputs for model anesthesia using sims_tauth
    sims_tauth = np.delete(sims_tauth,[0,1], 0)
    
    dd_nav=[]
    dd_com=[]
    di_nav=[]
    di_iso=[]
    di_com=[]
    ce_n_e=[]
    ce_n_u=[]
    ce_c_e=[]
    ce_c_u=[]
    rate_nav=[]
    rate_com=[]  

    for j in range(nx*ny):  
        # j is a particular combinations of parameters
        # 2D only needs the average
        ce_i=np.mean(context_effect(sims_tauth[j]),0) 
        dd_nav_i=(ce_i[1]-ce_i[0])/2
        dd_com_i=(ce_i[3]-ce_i[2])/2
        dd_nav.append(dd_nav_i)
        dd_com.append(dd_com_i)
        # context effect 
        ce_n_e.append(ce_i[0])
        ce_n_u.append(ce_i[1])
        ce_c_e.append(ce_i[2])
        ce_c_u.append(ce_i[3])
        # discriminability index 
        di_i=np.mean(discriminability_index(sims_tauth[j]),0) 
        di_nav.append(di_i[0])
        di_iso.append(di_i[1])
        di_com.append(di_i[2])
        # firing rate 
        sp_contexts=sims_tauth[j]
        rate_nav.append(np.mean(sp_contexts[:,8]))  #/dur_ech
        rate_com.append(np.mean(sp_contexts[:,11]))  #/dur_com
 
    dd_nav=np.array(dd_nav).reshape(nx,ny)
    dd_com=np.array(dd_com).reshape(nx,ny)
    di_nav=np.array(di_nav).reshape(nx,ny)
    di_iso=np.array(di_iso).reshape(nx,ny)
    di_com=np.array(di_com).reshape(nx,ny)
    ce_n_e=np.array(ce_n_e).reshape(nx,ny)
    ce_n_u=np.array(ce_n_u).reshape(nx,ny)
    ce_c_e=np.array(ce_c_e).reshape(nx,ny)
    ce_c_u=np.array(ce_c_u).reshape(nx,ny)
    rate_nav=np.array(rate_nav).reshape(nx,ny)
    rate_com=np.array(rate_com).reshape(nx,ny)
    
    # difference with awake state for echolocation
    dd=dd_nav-dd_nav_o
    ce_e=ce_n_e-ce_n_e_o
    ce_u=ce_n_u-ce_n_u_o
    rate=rate_nav-rate_nav_o
    
    # smoothing data
    x=parameters2
    y=parameters1
    f_dd = interp2d(x,y,dd, kind='cubic')
    f_ce_e = interp2d(x,y,ce_e, kind='cubic')
    f_ce_u = interp2d(x,y,ce_u, kind='cubic')
    f_rate = interp2d(x,y,rate, kind='cubic')

    # increasing resolution axes
    alpha=5   # smoothing factor
    
    x_in=parameters2[0]
    x_end=parameters2[-1]
    x_n=len(parameters2)
    y_in=parameters1[0]
    y_end=parameters1[-1]
    y_n=len(parameters1)
    x_nn=x_n + (x_n-1)*(alpha-1)
    y_nn=y_n + (y_n-1)*(alpha-1)
    xnew = np.linspace(x_in,x_end,x_nn)
    ynew = np.linspace(y_in,y_end,y_nn)
    
    Xn, Yn = np.meshgrid(xnew, ynew)
    # generate smooth data for the new axes
    dd_smooth = f_dd(xnew,ynew)
    ce_e_smooth = f_ce_e(xnew,ynew)
    ce_u_smooth = f_ce_u(xnew,ynew)
    rate_smooth = f_rate(xnew,ynew)


    
    # c.e. matching and mismatching
    fig1, ax1 = plt.subplots(1,2,sharex=True,sharey=True,figsize=(2.1,1))
    mm=0.2 
    q1=ax1[0].pcolor(Xn,Yn,ce_e_smooth,cmap='RdBu',vmin=-mm,vmax=mm)
    q2=ax1[1].pcolor(Xn,Yn,ce_u_smooth,cmap='RdBu',vmin=-mm,vmax=mm)

    # adding colorbars
    divider = make_axes_locatable(ax1[0])
    cax = divider.append_axes("right", size="5%", pad=0.075)
    fig1.colorbar(q1,cax=cax)
    cax.remove() 
    divider = make_axes_locatable(ax1[1])
    cax = divider.append_axes("right", size="5%", pad=0.075)
    fig1.colorbar(q2,cax=cax)
    
    fig1.tight_layout() 
    
    # s.s.s. 
    fig2, ax2 = plt.subplots(1,figsize=(1.15,1))
    q1=ax2.pcolormesh(Xn,Yn,dd_smooth,cmap='RdBu',vmin=-mm,vmax=mm)

    divider = make_axes_locatable(ax2)
    cax = divider.append_axes("right", size="5%", pad=0.075)
    fig2.colorbar(q1,cax=cax)
    
    fig2.tight_layout() 
    
    # firing rate of inputs 
    fig3, ax3 = plt.subplots(1,figsize=(1.15,1))
    mm=30 
    q1=ax3.pcolor(Xn,Yn,rate_smooth,cmap='RdBu',vmin=-mm,vmax=mm)
    ax3.contour(Xn,Yn,rate_smooth,levels = [-23],origin='lower',colors=['k'],linewidths=0.75)

    divider = make_axes_locatable(ax3)
    cax = divider.append_axes("right", size="5%", pad=0.075)
    fig3.colorbar(q1,cax=cax)
    
    fig3.tight_layout() 
    

    # difference with awake state for communication
    dd=dd_com-dd_com_o
    ce_e=ce_c_e-ce_c_e_o
    ce_u=ce_c_u-ce_c_u_o
    rate=rate_com-rate_com_o
    
    #smoothing data
    f_dd = interp2d(x,y,dd, kind='cubic')
    f_ce_e = interp2d(x,y,ce_e, kind='cubic')
    f_ce_u = interp2d(x,y,ce_u, kind='cubic')
    f_rate = interp2d(x,y,rate, kind='cubic')
    
    dd_smooth = f_dd(xnew,ynew)
    ce_e_smooth = f_ce_e(xnew,ynew)
    ce_u_smooth = f_ce_u(xnew,ynew)
    rate_smooth = f_rate(xnew,ynew)


    # c.e matching and mismatching 
    fig4, ax4 = plt.subplots(1,2,sharex=True,sharey=True,figsize=(2.1,1))
    mm=0.2
    # c.e. match 
    q1=ax4[0].pcolor(Xn,Yn,ce_e_smooth,cmap='RdBu',vmin=-mm,vmax=mm)
#    ax4[0].set_title('c.e. match')
    # c.e. mismatch
    q2=ax4[1].pcolor(Xn,Yn,ce_u_smooth,cmap='RdBu',vmin=-mm,vmax=mm)
#    ax4[1].set_title('c.e. mismatch')
    divider = make_axes_locatable(ax4[0])
    cax = divider.append_axes("right", size="5%", pad=0.075)
    fig4.colorbar(q1,cax=cax)
    cax.remove() 
    divider = make_axes_locatable(ax4[1])
    cax = divider.append_axes("right", size="5%", pad=0.075)
    fig4.colorbar(q2,cax=cax)
    
    fig4.tight_layout() 
    
    fig5, ax5 = plt.subplots(1,figsize=(1.15,1))  
    # s.s.s    
    q1=ax5.pcolor(Xn,Yn,dd_smooth,cmap='RdBu',vmin=-mm,vmax=mm)
#    ax5.set_title('s.s.s.')
    divider = make_axes_locatable(ax5)
    cax = divider.append_axes("right", size="5%", pad=0.075)
    fig5.colorbar(q1,cax=cax)  
    
    fig5.tight_layout() 
    
    fig6, ax6 = plt.subplots(1,figsize=(1.15,1))  
    # firing rate
    mm=30 
    q1=ax6.pcolor(Xn,Yn,rate_smooth,cmap='RdBu',vmin=-mm,vmax=mm)
    ax6.contour(Xn,Yn,rate_smooth,levels = [-16],origin='lower',colors=['k'],linewidths=0.75)
#    ax6.set_title('firing rate')
    divider = make_axes_locatable(ax6)
    cax = divider.append_axes("right", size="5%", pad=0.075)
    fig6.colorbar(q1,cax=cax)

    fig6.tight_layout() 





    # plotting for difference between communication and echolocation
    jump=3
    
    fig7, ax7 = plt.subplots(2,2,figsize=(6,4))
    fig7.canvas.manager.set_window_title('Difference between contexts')
    # difference with awake state
    # communication
    dd2=dd_com-dd_com_o
    ce_e2=ce_c_e-ce_c_e_o
    ce_u2=ce_c_u-ce_c_u_o
    rate2=rate_com-rate_com_o
    # echolocation
    dd1=dd_nav-dd_nav_o
    ce_e1=ce_n_e-ce_n_e_o
    ce_u1=ce_n_u-ce_n_u_o
    rate1=rate_nav-rate_nav_o
    # difference between communication and echolocation
    dd=dd2-dd1
    ce_e=ce_e2-ce_e1
    ce_u=ce_u2-ce_u1
    rate=rate2-rate1
    # first plot: s.s.s  
    mm=0.2 
    q1=ax7[0,0].imshow(dd,cmap='RdBu',origin='lower',vmin=-mm,vmax=mm,interpolation='hamming')
    ax7[0,0].set_xticks(np.arange(0,ny,jump))
    ax7[0,0].set_xticklabels(['{:.1f}'.format(x) for x in parameters2[::jump]])
    ax7[0,0].set_yticks(np.arange(0,nx,jump))
    ax7[0,0].set_yticklabels(['{:.1f}'.format(x) for x in parameters1[::jump]])
    ax7[0,0].set_title('s.s.s.')
    # third plot: c.e. mismatch
    q3=ax7[1,0].imshow(ce_e,cmap='RdBu',origin='lower',vmin=-mm,vmax=mm,interpolation='hamming')
    ax7[1,0].set_xticks(np.arange(0,ny,jump))
    ax7[1,0].set_xticklabels(['{:.1f}'.format(x) for x in parameters2[::jump]])
    ax7[1,0].set_yticks(np.arange(0,nx,jump))
    ax7[1,0].set_yticklabels(['{:.1f}'.format(x) for x in parameters1[::jump]])
    ax7[1,0].set_title('c.e. match')
    # forth plot: firing rate
    q4=ax7[1,1].imshow(ce_u,cmap='RdBu',origin='lower',vmin=-mm,vmax=mm,interpolation='hamming')
    ax7[1,1].set_xticks(np.arange(0,ny,jump))
    ax7[1,1].set_xticklabels(['{:.1f}'.format(x) for x in parameters2[::jump]])
    ax7[1,1].set_yticks(np.arange(0,nx,jump))
    ax7[1,1].set_yticklabels(['{:.1f}'.format(x) for x in parameters1[::jump]])
    ax7[1,1].set_title('c.e. mismatch')
    # second plot: c.e. match 
    mm=30
    q2=ax7[0,1].imshow(rate,cmap='RdBu',origin='lower',vmin=-mm,vmax=mm,interpolation='hamming')
    ax7[0,1].set_xticks(np.arange(0,ny,jump))
    ax7[0,1].set_xticklabels(['{:.1f}'.format(x) for x in parameters2[::jump]])
    ax7[0,1].set_yticks(np.arange(0,nx,jump))
    ax7[0,1].set_yticklabels(['{:.1f}'.format(x) for x in parameters1[::jump]])
    ax7[0,1].set_title('firing context')
    
    # add colorbars
    fig7.colorbar(q1,ax=ax7[0,0])
    fig7.colorbar(q2,ax=ax7[0,1])
    fig7.colorbar(q3,ax=ax7[1,0])
    fig7.colorbar(q4,ax=ax7[1,1])
    
    fig7.tight_layout()     
    
    
    return fig1,ax1,fig2,ax2,fig3,ax3,fig4,ax4,fig5,ax5,fig6,ax6

def one_combi(sims,x,y,sims_original=None):
    # get output from a particular combination of parameters (x and y) from sims
    parameters1=sims[0]
    parameters2=sims[1]
    sims = np.delete(sims,[0,1], 0)
    ny=len(parameters2)
    or_par1=x
    or_par2=y  
    ind_par1=np.where(np.round(parameters1,2)==or_par1)[0]
    ind_par2=np.where(np.round(parameters2,2)==or_par2)[0]  
    # from 2D to 1D
    pos_2d=int(ind_par1)*(ny)+int(ind_par2)
    # particular combination in x, y 
    comb_sim=sims[pos_2d]
    # combination for awake model 
    ind_par1_awake=np.where(np.round(parameters1,2)==1.0)[0]
    ind_par2_awake=np.where(np.round(parameters2,2)==1.0)[0]
    pos_2d_awake=int(ind_par1_awake)*(ny)+int(ind_par2_awake) 
    if sims_original is None:
        # take the awake values from the same matrix sims
        comb_sim_awake=sims[pos_2d_awake]
        print('using same matrix to find awake state')
    else:
        print('using second matrix provided in the function')
        # take the awake values from the given matrix 'sims_original'
        sims_original = np.delete(sims_original,[0,1], 0)
        comb_sim_awake=sims_original[pos_2d_awake]
    
    if np.array_equal(comb_sim,comb_sim_awake):
        print('equals')

    # calculation of outputs for combination x,y
    ce_comb=context_effect(comb_sim)
    ce_n_e_comb=ce_comb[:,0]
    ce_n_u_comb=ce_comb[:,1]
    ce_c_e_comb=ce_comb[:,2]
    ce_c_u_comb=ce_comb[:,3]
    
    p_anest_nav=ranksums(ce_n_e_comb, ce_n_u_comb)[1]
    p_anest_com=ranksums(ce_c_e_comb, ce_c_u_comb)[1]

    dd_nav_comb=(ce_n_u_comb-ce_n_e_comb)/2
    dd_com_comb=(ce_c_u_comb-ce_c_e_comb)/2
    
    rate_nav_comb=np.mean(comb_sim[:,8].reshape(50,20),1)
    rate_com_comb=np.mean(comb_sim[:,11].reshape(50,20),1)
    
    # calculation of outputs for combination awake
    ce_awake=context_effect(comb_sim_awake)
    ce_n_e_awake=ce_awake[:,0]
    ce_n_u_awake=ce_awake[:,1]
    ce_c_e_awake=ce_awake[:,2]
    ce_c_u_awake=ce_awake[:,3]
    
    p_awake_nav=ranksums(ce_n_e_awake, ce_n_u_awake)[1]
    p_awake_com=ranksums(ce_c_e_awake, ce_c_u_awake)[1]

    dd_nav_awake=(ce_n_u_awake-ce_n_e_awake)/2
    dd_com_awake=(ce_c_u_awake-ce_c_e_awake)/2
    
    rate_nav_awake=np.mean(comb_sim_awake[:,8].reshape(50,20),1)
    rate_com_awake=np.mean(comb_sim_awake[:,11].reshape(50,20),1)
    
    # stats 
#    print('anesthedsia')
#    print(wilcoxon(dd_com_comb))
#    print('awake')
#    print(wilcoxon(dd_com_awake))
    if p_awake_nav<0.05:
        stats1='*'
    else:
        stats1='ns'
    if p_anest_nav<0.05:
        stats2='*'
    else:
        stats2='ns'
    if p_awake_com<0.05:
        stats3='*'
    else:
        stats3='ns'
    if p_anest_com<0.05:
        stats4='*'
        k=0
    else:
        stats4='ns'
        k=1
  
    ## plotting 
    red_ =(0.8500,0.3250, 0.0980)
    blue_=(0,0.4470,0.7410)
    aa=0.5
    red_light=(0.8500,0.3250*aa, 0.0980*aa)
    blue_light=(0,0.4470*aa,0.7410*aa)
    colors=[blue_,blue_light,red_,red_light]
    v=0     # capsize 
    t=1     # thick of the cap
    ss=2    # size of the marker 'o'
    fH=1.4  # figure size weight
    fW=1.7   # figure size height
    rot=30  # rotation of xlabels
    ew=1    # thickness of errorbars 
    
    # context effect echolocation
    fig1, ax = plt.subplots(1,figsize=(fH,fW))
    # add boxplot 
    colorsb1=[nc.to_rgba(i, alpha=0.5) for i in colors]
    b1=ax.boxplot([ce_n_e_awake,ce_n_u_awake,ce_n_e_comb,ce_n_u_comb],positions=[1,2,3,4], patch_artist=True,widths=0.5,notch=True)
    for item in ['whiskers', 'fliers', 'medians', 'caps','boxes']:
        for c in b1[item]:
            c.set(color='k', linewidth=0.5)     
    for patch, color in zip(b1['boxes'], colorsb1):
        patch.set_facecolor(color)
    for i in b1['fliers']:
        i.set_markersize(2)
        i.set_markeredgewidth(0.5)

    ax.errorbar(1, np.median(ce_n_e_awake), yerr = sem(ce_n_e_awake),c=colors[0], fmt='o',ms=ss, capsize=v, capthick=t,elinewidth=ew)
    ax.errorbar(2, np.median(ce_n_u_awake), yerr = sem(ce_n_u_awake),c=colors[1], fmt='o',ms=ss, capsize=v, capthick=t,elinewidth=ew)
    ax.errorbar(3, np.median(ce_n_e_comb), yerr = sem(ce_n_e_comb),c=colors[2], fmt='o',ms=ss, capsize=v, capthick=t,elinewidth=ew)
    ax.errorbar(4, np.median(ce_n_u_comb), yerr = sem(ce_n_u_comb),c=colors[3], fmt='o',ms=ss, capsize=v, capthick=t,elinewidth=ew)
    
    # add stats
    y=0.1 #-0.2
    h=0.01
    ax.plot([1, 1, 2, 2], [y, y+h, y+h, y], lw=0.75,c='k')
    ax.plot([3, 3, 4, 4], [y, y+h, y+h, y], lw=0.75,c='k')
    ax.text((1+2)*.5, y+h, stats1, ha='center', va='center',c='k',fontsize=12)
    ax.text((3+4)*.5, y+h, stats2, ha='center', va='center',c='k',fontsize=12)
    
    ax.set_ylabel('echolocation effect')
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.spines['bottom'].set_linewidth(1)
    ax.spines['left'].set_linewidth(1)
    ax.set_xlim([0.5, 4.5])
#    ax.set_ylim([-0.9, 0.175])
    ax.set_xticks(np.arange(1,5,1))
    ax.set_xticklabels(['match','mismatch','match','mismatch'],ha='right')
    ax.tick_params(axis='x', rotation=rot)
    fig1.tight_layout()
    
    # context effect communication
    fig2, ax = plt.subplots(1,figsize=(fH,fW))
    
    # add boxplot 
    colorsb1=[nc.to_rgba(i, alpha=0.5) for i in colors]
    b1=ax.boxplot([ce_c_e_awake,ce_c_u_awake,ce_c_e_comb,ce_c_u_comb],positions=[1,2,3,4], patch_artist=True,widths=0.5,notch=True)
    for item in ['whiskers', 'fliers', 'medians', 'caps','boxes']:
        for c in b1[item]:
            c.set(color='k', linewidth=0.5)     
    for patch, color in zip(b1['boxes'], colorsb1):
        patch.set_facecolor(color)
    for i in b1['fliers']:
        i.set_markersize(2)
        i.set_markeredgewidth(0.5)
          
    ax.errorbar(1, np.median(ce_c_e_awake), yerr = sem(ce_c_e_awake),c=colors[0], fmt='o',ms=ss, capsize=v, capthick=t,elinewidth=ew)
    ax.errorbar(2, np.median(ce_c_u_awake), yerr = sem(ce_c_u_awake),c=colors[1], fmt='o',ms=ss, capsize=v, capthick=t,elinewidth=ew)
    ax.errorbar(3, np.median(ce_c_e_comb), yerr = sem(ce_c_e_comb),c=colors[2], fmt='o',ms=ss, capsize=v, capthick=t,elinewidth=ew)
    ax.errorbar(4, np.median(ce_c_u_comb), yerr = sem(ce_c_u_comb),c=colors[3], fmt='o',ms=ss, capsize=v, capthick=t,elinewidth=ew)


    # add stats
    y=0.1 #-0.25 #-0.3
    h=0.01
    ax.plot([1, 1, 2, 2], [y, y+h, y+h, y], lw=0.75,c='k')
    ax.plot([3, 3, 4, 4], [y, y+h, y+h, y], lw=0.75,c='k')
    ax.text((1+2)*.5, y+h, stats3, ha='center', va='center',c='k',fontsize=12)
    ax.text((3+4)*.5, y+h, stats4, ha='center', va='center',c='k',fontsize=12)

    ax.set_ylabel('communication effect')    
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.spines['bottom'].set_linewidth(1)
    ax.spines['left'].set_linewidth(1)
    ax.set_xlim([0.5, 4.5])
#    ax.set_ylim([-0.5, -0.275])
    ax.set_xticks(np.arange(1,5,1))
    ax.set_xticklabels(['match','mismatch','match','mismatch'],ha='right')
    ax.tick_params(axis='x', rotation=rot)
    fig2.tight_layout()
    
    # deviance detection both context
    fig3, ax = plt.subplots(1,figsize=(fH,fW))
    ax.axhline(y=0,color='gray',lw=1,ls='--') 
    # add boxplot 
    colorsb1=[nc.to_rgba(i, alpha=0.5) for i in colors]
    b1=ax.boxplot([dd_nav_awake,dd_com_awake,dd_nav_comb,dd_com_comb],positions=[1,2,3,4], patch_artist=True,widths=0.5,notch=True)
    for item in ['whiskers', 'fliers', 'medians', 'caps','boxes']:
        for c in b1[item]:
            c.set(color='k', linewidth=0.5)     
    for patch, color in zip(b1['boxes'], colorsb1):
        patch.set_facecolor(color)
    for i in b1['fliers']:
        i.set_markersize(2)
        i.set_markeredgewidth(0.5)
        
    
    ax.errorbar(1, np.median(dd_nav_awake), yerr = sem(dd_nav_awake),c=colors[0], fmt='o',ms=ss, capsize=v, capthick=t,elinewidth=ew)
    ax.errorbar(2, np.median(dd_com_awake), yerr = sem(dd_nav_comb),c=colors[1], fmt='o',ms=ss, capsize=v, capthick=t,elinewidth=ew)
    ax.errorbar(3, np.median(dd_nav_comb), yerr = sem(dd_com_awake),c=colors[2], fmt='o',ms=ss, capsize=v, capthick=t,elinewidth=ew)
    ax.errorbar(4, np.median(dd_com_comb), yerr = sem(dd_com_comb),c=colors[3], fmt='o',ms=ss, capsize=v, capthick=t,elinewidth=ew)

    ax.set_ylabel('s.s.s')
#    ax.set_ylim([0,0.2])     
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.spines['bottom'].set_linewidth(1)
    ax.spines['left'].set_linewidth(1)
    ax.set_xlim([0.5, 4.5])
    ax.set_ylim([-0.2, 0.4])
    ax.set_xticks(np.arange(1,5,1))
    ax.set_xticklabels(['ech','com','ech','com'],ha='right')
    ax.tick_params(axis='x', rotation=rot)
    fig3.tight_layout()
    
    # rate firing both context 
    fig4, ax = plt.subplots(1,figsize=(fH,fW))
#    ax.errorbar(1, np.mean(rate_nav_awake), yerr = sem(rate_nav_awake),c=colors[0], fmt='o',ms=ss, capsize=v, capthick=t,elinewidth=ew)
#    ax.errorbar(2, np.mean(rate_nav_comb), yerr = sem(rate_nav_comb),c=colors[2], fmt='o',ms=ss, capsize=v, capthick=t,elinewidth=ew)
#    ax.errorbar(3, np.mean(rate_com_awake), yerr = sem(rate_com_awake),c=colors[1], fmt='o',ms=ss, capsize=v, capthick=t,elinewidth=ew)
#    ax.errorbar(4, np.mean(rate_com_comb), yerr = sem(rate_com_comb),c=colors[3], fmt='o',ms=ss, capsize=v, capthick=t,elinewidth=ew)

    ax.bar(1, np.mean(rate_nav_awake),color=colors[0],alpha=0.5,edgecolor=(0,0,0,1))
    ax.bar(2, np.mean(rate_nav_comb),color=colors[2],alpha=0.5,edgecolor=(0,0,0,1))
    ax.bar(4, np.mean(rate_com_awake),color=colors[1],alpha=0.5,edgecolor=(0,0,0,1))
    ax.bar(5, np.mean(rate_com_comb),color=colors[3],alpha=0.5,edgecolor=(0,0,0,1))
    
    print(p_awake_nav)               # (wilcoxon(dd_nav_awake))
    print(p_anest_nav)               #(wilcoxon(dd_nav_comb))
    print(p_awake_com)               #(wilcoxon(dd_com_awake))
    print(p_anest_com)               #(wilcoxon(dd_com_comb))
    
#    print(cliffsDelta(rate_nav_awake,rate_nav_comb)[0])
#    print(cliffsDelta(rate_com_awake,rate_com_comb)[0])
    

    ax.set_ylabel('spike count')
#    ax.set_ylim([-1,0]) 
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.spines['bottom'].set_linewidth(1)
    ax.spines['left'].set_linewidth(1)
    ax.set_xlim([0, 6])
    ax.set_xticks([1,2,4,5])
    ax.set_xticklabels(['ech-aw','ech-an','com-aw','com-an'],ha='right')
    ax.tick_params(axis='x', rotation=rot)
    fig4.tight_layout()
      
    plt.locator_params(axis='y',nbins=5)
    
    
    # plotting with boxplots
#    fig1, ax1 = plt.subplots(1,1,figsize=(1.5,1.5))
#    colors=[nc.to_rgba(i, alpha=0.5) for i in colors]
#    b1=ax1.boxplot([ce_c_e_awake,ce_c_u_awake,ce_c_e_comb,ce_c_u_comb],positions=[1,2,3,4], patch_artist=True,widths=0.5,notch=True)
#    for item in ['whiskers', 'fliers', 'medians', 'caps','boxes']:
#        for c in b1[item]:
#            c.set(color='k', linewidth=0.5)     
#    for patch, color in zip(b1['boxes'], colors):
#        patch.set_facecolor(color)
#    for i in b1['fliers']:
#        i.set_markersize(2)
#        i.set_markeredgewidth(0.5)


    return fig1, fig2, fig3, fig4, k