#cp drawpresynGABAB_actpasaxon.py drawfig_fitting_composite.py
#cp presynAP_test_gCaNtype_axonspine_train_withrelease_testgamma_withreluctant_varygbars.py presynGABAB.py
from neuron import h,rxd
from pylab import *
import mytools
import scipy.io

Nstims = 2
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
blockedStr = 'None'
blockeds = []
block_factors = []
isyntype = 1 #1 for E->E connection (OnlyExp1_RGS used), 2 for I->E connection (OnlyExp2_RGS used)

spineLen = 0.8
spineDiam = 0.4
doInterpolate = 1
doDraw = 0
fileOut = 'None'
stimAmp = 0.3
kcoeff_f = 5e-6
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','117.497,596.380,166.400,10.972,56.312,2.905,1.798,2.486,2569.532','133.092,528.522,185.581,8.574,25.868,3.186,2.125,2.729,2815.426','91.979,819.933,169.073,5.239,22.576,4.547,2.153,2.962,525.594','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','127.868,614.980,360.316,2.912,15.822,3.778,2.492,3.119,1872.786','119.681,591.014,91.960,89.687,88.948,3.453,2.300,2.755,2950.763','170.128,194.903,421.922,2.138,22.729,2.763,1.791,2.151,2061.804']


data_inppr_ydata_nonscaled = [0.4978165938864628, 0.6462882096069869, 0.6965065502183405, 0.7467248908296943, 0.8275109170305677, 0.8668122270742359, 0.8973799126637554,  ] #https://link.springer.com/content/pdf/10.1007/BF00169135.pdf (Olpe et al. 1994)
data_inppr_ydata_GABABblocked = [0.8078602620087335, 0.7991266375545851, 0.7991266375545851, 0.8296943231441048, 0.8362445414847161, 0.868995633187773, 0.9301310043668121 ] #

data_exppr_scaled = [[147.51475085724883, 0.7605217633066633], #Local and diffuse synaptic actions of GABA in the hippocampus. J.S.Isaacson, J.M.Solis, R.A.Nicoll, 1993
                         [198.375088814062, 0.7229294739118345],
                         [299.67563559976526, 0.664576318309598],
                         [497.7047357202435, 0.7310354947329398],
                         [693.0153532482777, 0.8064038800160637],
                         [892.7743968366749, 0.903557196255908],
                         [1295.6040900806277, 0.9651926106700441],
                         [1591.7704117883293, 1.0000293472552593]]
data_exppr_nonscaled = [[153.11004784688993, 97.40181268882175],
                        [202.10526315789474, 100],
                        [303.157894736842, 100.12084592145015],
                        [500.6698564593301, 100.96676737160121],
                        [900.287081339713, 101.99395770392749],
                        [1601.531100478469, 98.18731117824774],
                        [150.0478468899521, 76.13293051359517],
                        [199.0430622009569, 72.44712990936556],
                        [297.0334928229666, 66.58610271903322],
                        [499.13875598086116, 73.05135951661632],
                        [898.7559808612441, 90.33232628398791],
                        [1600, 100]]
cols = ['#99AA00','#6600FF']
redcols = ['#FF6600','#CC0044']


f,axarr = subplots(1,2)
cols = mytools.colorsredtolila(5,1.0)+mytools.colorsredtolila(5,0.5)+mytools.colorsredtolila(5,0.9,0.7)
axnew = []
for iax in range(0,2):
  axarr[iax].set_position([0.07+0.48*iax,0.1,0.41,0.8])
  axnew.append(f.add_axes([0.07+0.48*iax+0.2,0.15,0.2,0.15]))

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

for iax in range(0,2):
  axarr[iax].tick_params(axis='both', which='major', labelsize=5)
  axnew[iax].tick_params(axis='both', which='major', labelsize=5)
  boxoff(axarr[iax])
  boxoff(axnew[iax])

for params_str in ['p'+str(i) for i in range(0,15)]:
  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])}

  gaba_input_flux = params['gaba_flux']

  styles = ['-','--']
  widths = [0.3,0.6]

  toplot = []
  meanerrors = []
  for isyntype in [0,1]:
   if isyntype == 1:
    Ts = [150.0, 200.0, 300.0, 500.0, 700.0, 900.0, 1300.0, 1600.0]
    if iparam == 0:
      axarr[isyntype].plot(Ts, array(data_exppr_scaled)[:,1], 'k.',zorder=10,label='Exp.') #. data, with GABAB\nrelative to GABAB-blocked')
    myexpdata = array(data_exppr_scaled)[:,1]
   else:
    Ts = [350.0, 550.0, 750.0, 1000.0, 2000.0, 3000.0, 4000.0]
    if iparam == 0:
      axarr[isyntype].plot(Ts, [data_inppr_ydata_nonscaled[i]/data_inppr_ydata_GABABblocked[i] for i in range(0,len(data_inppr_ydata_nonscaled))],'k.',zorder=10,label='Exp.') #'Exp. data, with GABAB\nrelative to GABAB-blocked')
    myexpdata = [data_inppr_ydata_nonscaled[i]/data_inppr_ydata_GABABblocked[i] for i in range(0,len(data_inppr_ydata_nonscaled))]

   NT_amps_all = []
   NT_firstamps_all = []
   times_all_all = []
   maxPoint = 0

   for iGABAB in [0,1]:
    NTs_all = []
    times_all = []
    NT_amps = []
    NT_firstamps = []
    for iT in range(0,len(Ts)):
      addition = ('_no1stglu_EE' if isyntype == 1 else '_IE') + ('' if iGABAB else '_GABABRx0.0')
      A = scipy.io.loadmat("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)+"_"+str(Ts[iT])+"ms_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('Loeaded '+"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)+"_"+str(Ts[iT])+"ms_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')
      times = A['times'][0]
      NT = A['NT'][0]
      NT_amps.append(max([NT[i] for i in range(0,len(times)) if times[i] >= 1000+Ts[iT]]))
      NT_firstamps.append(max([NT[i] for i in range(0,len(times)) if times[i] < 1000+Ts[iT]]))
      NTs_all.append(NT[:])

    NT_amps_all.append(NT_amps[:])
    NT_firstamps_all.append(NT_firstamps[:])
    maxPoint = max(maxPoint,max(NT_amps))
    #print("iGABAB="+str(iGABAB)+", min max = "+str(min([NT_amps_all[iGABAB][i] for i in range(0,len(NT_amps_all[0]))]))+" "+str(max([NT_amps_all[iGABAB][i] for i in range(0,len(NT_amps_all[0]))])))
   meanerror = mean([abs(NT_amps_all[1][i]/NT_amps_all[0][i] - myexpdata[i]) for i in range(0,len(NT_amps_all[0]))])
   print(params_str+', isyntype='+str(isyntype)+', mean PPR error = '+str(meanerror))
   toplot.append([Ts,[NT_amps_all[1][i]/NT_amps_all[0][i] for i in range(0,len(NT_amps_all[0]))]])
   axnew[isyntype].bar(iparam,meanerror,facecolor=cols[iparam])
   meanerrors.append(meanerror)
  for isyntype in [0,1]:
   axarr[isyntype].plot(toplot[isyntype][0],toplot[isyntype][1],'k-' if mean(meanerrors) < 0.05 else 'k--',lw=0.9+0.9*(iparam==0),color=cols[iparam],label='Pred., param '+str(iparam))
  axarr[isyntype].set_ylim([0.6,1.02])

  f.savefig("fig2_suppAB.pdf")
 
axarr[0].set_ylabel('GABABR-mediated PPR',fontsize=8)
axarr[0].set_xlabel('ISI (ms)',fontsize=8)
axarr[1].set_xlabel('ISI (ms)',fontsize=8)
axarr[0].set_title('I->E synapse')
axarr[1].set_title('E->E synapse')
for iax in range(0,2):
  axarr[iax].legend(fontsize=4,frameon=False,loc=5)
  axnew[iax].set_xticks([])
  axnew[iax].set_ylabel('Mean PPR error',fontsize=6)
  pos = axarr[iax].get_position()
  f.text(pos.x0 - 0.05, pos.y1 - 0.01+0.02, chr(ord('A')+iax), fontsize=12)

  
f.savefig("fig2_suppAB.pdf")