from pylab import *
import mytools
import scipy.io

def boxoff(ax,whichxoff='top'):
    ax.spines[whichxoff].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.get_xaxis().tick_bottom()
    ax.get_yaxis().tick_left()

FACILITATION_INDEX_DATA_ISAACSON = array([
  [-0.09293747952584752, 1.9327628880989067],
  [1.5117506460732102, 1.9408752017083024],
  [3.036200321115576, 2.0736007813542665],
  [4.578687470426228, 1.8721139030101552],
  [6.066091570513986, 1.927796010725422],
  [7.497765537100172, 1.9936681792261681],
  [9.124050100500277, 1.900952123851931],
  [10.558150633131522, 1.955495262129797],
  [12.067151171020328, 1.910349000861431],
  [13.503192956487624, 1.9558289149609929],
  [14.954036794829797, 1.9322017447009863],
  [16.380372316116844, 2.022997779692069],
  [17.95165471582889, 1.937073076036447],
  [19.442698664984245, 1.9757616384173935],
  [21.005002770329572, 1.9317543465864282],
  [22.460214427552845, 1.8877349219252375],
  [23.976252006972327, 1.8097344730105192],
  [25.593558276005705, 1.7589358294609385],
  [27.052895095505587, 1.6956570534208515],
  [28.48917953757739, 1.7400040644981254],
  [29.994540226398602, 1.7118513485640796],
  [31.59340459348952, 1.7471533347083876],
  [33.09415480682512, 1.7405257761978137],
  [34.64869390082623, 1.7327713810800647],
  [36.023181794283815, 1.8156310285000181],
  [37.58742715246518, 1.7625605124907489],
  [39.07046343367183, 1.8386348746071997]])



stim_T = 50.0
stim_T_GABAB = 200.0 # GABAB stimulation at 1000.0, el. stimulation at 1000.0+stim_T_GABAB and 1000.0+stim_T_GABAB+stim_T
Nstims = 3
decay = 1.0
axonPasLen = 100
axonActLen = 20

Vesmax = 0.001    #Concentration of vesicles (mM)
Fmax = 0.001 #Concentration of fusion factor (mM)
initfile = 'None'
tolstochange = 'Gi,GiaGDP,Gibg'.split(',')
tolschange = [float(x) for x in '2,2,2'.split(',')]
gaba_input_dur = 3.0
isyntype = 1 #1 for E->E connection (OnlyExp1_RGS used), 2 for I->E connection (OnlyExp2_RGS used)

spineLen = 0.8
spineDiam = 0.4
stimAmp = 0.3
kcoeff_f = 5e-06 
kcoeff_b = 1.0
kcoeff_2f = 100.0
kcoeff_2b = 1.0
kcoeff_3 = 1.0
kcoeff_4 = 1.0

gnabar = 3.9  #L5PCbiophysics4.hoc
gkbar = 0.077 #L5PCbiophysics4.hoc
gl = 0.00003  #L5PCbiophysics4.hoc
el = -90.0    #L5PCbiophysics4.hoc
gCaNtypebar_single = 2.7e-6 #uS (https://www.nature.com/articles/nn.2657) #13e-6 #uS (https://www.ncbi.nlm.nih.gov/books/NBK6181/)
minCa = 50e-6
params_all = ['107.034,514.651,167.633,11.026,66.931,3.009,1.904,2.241,1507.740','91.805,846.811,67.856,9.858,71.217,5.000,1.817,2.170,3000.000','133.092,528.522,185.581,8.574,25.868,3.186,2.125,2.729,2815.426','108.557,729.183,292.710,2.473,16.923,4.406,2.127,2.973,696.729','111.531,708.687,165.616,10.797,93.218,2.903,1.728,2.235,2797.567','109.869,469.365,114.137,59.053,73.263,3.573,2.385,3.028,2852.993','84.387,57.393,126.013,4.253,57.544,5.000,1.722,2.122,929.204','126.993,135.047,348.945,3.733,18.257,3.800,2.647,3.237,2804.842','105.438,642.913,142.107,19.602,40.482,3.599,2.391,2.974,1386.488','91.403,459.570,97.658,3.443,65.833,5.000,1.540,1.942,2447.114','71.410,380.140,143.521,2.991,45.796,10.000,2.044,2.836,2126.204','146.055,702.846,310.725,3.041,13.076,3.334,2.351,3.166,2112.405','90.417,993.718,98.403,6.879,95.287,3.944,1.572,2.204,1413.654','119.681,591.014,91.960,89.687,88.948,3.453,2.300,2.755,2950.763','123.170,616.560,117.913,74.842,46.625,3.952,2.609,3.283,3000.000']


params_str = 'p0'

if len(sys.argv) > 1:
  gnabar = float(sys.argv[1])
if len(sys.argv) > 2:
  gkbar = float(sys.argv[2])
if len(sys.argv) > 3:
  gl = float(sys.argv[3])
if len(sys.argv) > 4:
  el = float(sys.argv[4])
if len(sys.argv) > 5:
  decay = float(sys.argv[5])
if len(sys.argv) > 6:
  minCa = float(sys.argv[6])
if len(sys.argv) > 7:
  params_str = sys.argv[7]
if len(sys.argv) > 9:
  stim_T = float(sys.argv[9])
if len(sys.argv) > 10:
  stim_T_GABAB = float(sys.argv[10])
if len(sys.argv) > 11:
  axonPasLen = int(float(sys.argv[11]))
if len(sys.argv) > 12:
  axonActLen = int(float(sys.argv[12]))
if len(sys.argv) > 14:
  stimAmp = float(sys.argv[14])
if len(sys.argv) > 15:
  kcoeff_f = float(sys.argv[15])
if len(sys.argv) > 16:
  kcoeff_b = float(sys.argv[16])
if len(sys.argv) > 17:
  kcoeff_2f = float(sys.argv[17])
if len(sys.argv) > 18:
  kcoeff_2b = float(sys.argv[18])
if len(sys.argv) > 19:
  kcoeff_3 = float(sys.argv[19])
if len(sys.argv) > 20:
  kcoeff_4 = float(sys.argv[20])
if len(sys.argv) > 21:
  Fmax = float(sys.argv[21])
if len(sys.argv) > 22:
  Vesmax = float(sys.argv[22])
if len(sys.argv) > 23:
  doInterpolate = int(float(sys.argv[23]))
if len(sys.argv) > 24:
  doDraw = int(float(sys.argv[24]))
if len(sys.argv) > 25:
  fileOut = sys.argv[25]

if params_str[0] == 'p':
  iparam = int(params_str[1:])
  params = {'k[0]': float(params_all[iparam].split(',')[0]), 'k[1,9]': float(params_all[iparam].split(',')[1]), 'k[3,4,5,6,7,8]': float(params_all[iparam].split(',')[2]), 'k[15,17,19,21,23]': float(params_all[iparam].split(',')[3]), 'k[16,18,20,22,24]': float(params_all[iparam].split(',')[4]), 'OnlyExp0_RGS': float(params_all[iparam].split(',')[5]), 'OnlyExp1_RGS': float(params_all[iparam].split(',')[6]), 'OnlyExp2_RGS': float(params_all[iparam].split(',')[7]), 'gaba_flux': float(params_all[iparam].split(',')[8])}
else:
  params = {'k[0]': float(params_str.split(',')[0]), 'k[1,9]': float(params_str.split(',')[1]), 'k[3,4,5,6,7,8]': float(params_str.split(',')[2]), 'k[15,17,19,21,23]': float(params_str.split(',')[3]), 'k[16,18,20,22,24]': float(params_str.split(',')[4]), 'OnlyExp0_RGS': float(params_str.split(',')[5]), 'OnlyExp1_RGS': float(params_str.split(',')[6]), 'OnlyExp2_RGS': float(params_str.split(',')[7]), 'gaba_flux': float(params_str.split(',')[8])}

gaba_input_flux = params['gaba_flux']

blockeds = ['GABABRx0.0','None']
f,axarr = subplots(1,7)
for iax in range(0,2):
  axarr[iax].set_position([0.07+0.4*iax,0.6,0.32,0.3])
  axarr[2+iax].set_position([0.07+0.4*iax+0.17,0.6+0.1,0.04,0.16])
  axarr[5+iax].set_position([0.07+0.4*iax,0.48,0.32,0.07])
axarr[4].set_position([0.91,0.6,0.04,0.27])
for iax in range(0,7):
  axarr[iax].tick_params(axis='both', which='major', labelsize=5)
  boxoff(axarr[iax])

from matplotlib.collections import PatchCollection
for i in [0,1]:
  for iburst in [0]:
    polygon = Polygon(array([[1000+20*iburst,1003+20*iburst,1003+20*iburst,1000+20*iburst],[0,0,5,5]]).T)
    p = PatchCollection([polygon], cmap=matplotlib.cm.jet)
    p.set_facecolor('#BBBBBB')
    p.set_edgecolor(None)
    axarr[i].add_collection(p)

amps = []
cols = ['#99AA00','#6600FF']
for iblocked in [0,1]:
  addition = '_no1stglu_recVGCC_EE'
  blockedStr = blockeds[iblocked]
  if len(blockedStr) > 0 and blockedStr != 'None':
    addition = addition+'_'+blockedStr

  filename = "presynGABA_rxdrel6"+addition+"_gnabar"+str(gnabar)+"_gkbar"+str(gkbar)+"_gl"+str(gl)+"_el"+str(el)+"_decay"+str(decay)+'_minCa'+str(minCa)+"_"+params_str+"_Ves"+str(Vesmax)+"_Fmax"+str(Fmax)+"_200.0ms_stim_T_last50.0ms_N"+str(Nstims)+"_len"+str(axonPasLen)+","+str(axonActLen)+'_amp'+str(stimAmp)+'_kcoeff'+str(kcoeff_f)+"_"+str(kcoeff_b)+"_"+str(kcoeff_2f)+"_"+str(kcoeff_2b)+"_"+str(kcoeff_3)+"_"+str(kcoeff_4)+'.mat'

  print('Loading '+filename)
  A = scipy.io.loadmat(filename)
  
  axarr[iblocked].plot(A['times'][0],A['NT'][0],lw=0.5,color=cols[iblocked])
  amps.append([max([A['NT'][0][i] for i in range(0,len(A['times'][0])) if A['times'][0][i] < 1000+stim_T_GABAB+0.5*stim_T]),max([A['NT'][0][i] for i in range(0,len(A['times'][0])) if A['times'][0][i] >= 1000+stim_T_GABAB+0.5*stim_T])])
  axarr[5+iblocked].plot(A['times'][0],A['VGCCrec'][0]/A['VGCCrec'][0][0],'k-',lw=0.5,color=cols[iblocked])
  axarr[5+iblocked].plot(A['times'][0],A['VGCCGibgrec'][0]/A['VGCCrec'][0][0],'k--',lw=0.5,dashes=(1,1),color=cols[iblocked])

  axarr[iblocked].set_xlim([990,1280])
  axarr[5+iblocked].set_xlim([990,1280])
  axarr[iblocked].set_ylim([0,5])
  axarr[iblocked].plot([1150,1200],[4.8,4.8],'k-',lw=1)
  axarr[iblocked].text(1175,4.85,'50 ms',fontsize=5,ha='center',va='bottom')
  axarr[5+iblocked].plot([1150,1200],[0.88-0.29*(iblocked==0)]*2,'k-',lw=1)
  axarr[5+iblocked].text(1175,0.93-0.29*(iblocked==0),'50 ms',fontsize=5,ha='center',va='bottom')
  axarr[2+iblocked].bar(1,amps[-1][0],facecolor=cols[iblocked])
  axarr[2+iblocked].bar(2,amps[-1][1],facecolor=cols[iblocked])
  for ibar in [0,1]:
    axarr[2+iblocked].text(1+ibar,amps[-1][ibar]+0.1,'#'+str(1+ibar),ha='center',va='bottom',fontsize=6)
  axarr[2+iblocked].set_ylim([0,5])
print(str(amps))
axarr[4].bar(1,amps[0][1]/amps[0][0],facecolor=cols[0])
axarr[4].bar(2,amps[1][1]/amps[1][0],facecolor=cols[1])
DATA_STD_approximation = sqrt((std(FACILITATION_INDEX_DATA_ISAACSON[17:25,1])/mean(FACILITATION_INDEX_DATA_ISAACSON[0:16,1]))**2 + (mean(FACILITATION_INDEX_DATA_ISAACSON[17:25,1])*std(FACILITATION_INDEX_DATA_ISAACSON[0:16,1])/mean(FACILITATION_INDEX_DATA_ISAACSON[0:16,1])**2)**2)
print("Facilitation index in control "+str((amps[0][1]/amps[0][0])/(amps[1][1]/amps[1][0])*100)+"% from that in GABABR blocked, in the data: "+str(mean(FACILITATION_INDEX_DATA_ISAACSON[0:16,1]))+" vs "+str(mean(FACILITATION_INDEX_DATA_ISAACSON[17:25,1]))+" = "+str(mean(FACILITATION_INDEX_DATA_ISAACSON[17:25,1])/mean(FACILITATION_INDEX_DATA_ISAACSON[0:16,1])*100)+"% +- "+str(DATA_STD_approximation*100))


for i in range(0,2):
  axarr[i].set_ylabel('[NT] (mM)',fontsize=7)
  axarr[2+i].set_ylabel('peak [NT] (mM)',fontsize=7)
  myleg0 = mytools.mylegend(f,[0.07+0.4*i,0.48+0.07,0.15,0.02],['-','--'],['Willing','Reluctant'],nx=2,dx=2,yplus=0.5,yplustext=0.35,colors=[cols[i]]*2,dashes=[(20,0),(1,1)],linewidths=[0.5,0.5],myfontsize=5)
  for q in ['top','bottom','left','right']:
    myleg0.spines[q].set_linewidth(0.0)
  axarr[5+i].set_ylabel('VGCC\nfraction',fontsize=6)
  axarr[5+i].set_xticks([])
  
axarr[4].set_ylabel('facilitation index\n(second peak/first peak)',fontsize=7)
for i in range(0,5):
  axarr[i].set_xticks([])
axarr[0].set_title('GABABR blocked',fontsize=7)
axarr[1].set_title('Control',fontsize=7)

myleg = mytools.mylegend(f,[0.9,0.88,0.09,0.07],['-','-'],['GABABR\nblocked','Control'],nx=1,dx=2,yplus=0.5,yplustext=0.35,colors=cols,dashes=[],linewidths=[],myfontsize=6)
for q in ['top','bottom','left','right']:
  myleg.spines[q].set_linewidth(0.0)

for iax in range(0,3):
  pos = axarr[iax+2*(iax==2)].get_position()
  f.text(pos.x0 - 0.05, pos.y1 + 0.01+0.02, chr(ord('C')+iax), fontsize=12)

f.text(0.5,0.999, 'Effects of GABABR activation on paired-pulse facilitation in E->E synapse',ha='center',va='top')

f.savefig("fig2_suppCDE.pdf")