#cp fitterb8.py  fitterc.py # Use IPSC PPR in addition to EPSC PPR
#cp fitterb5.py fitterb8.py # larger ranges
#cp fitter.py fitter2.py # allow more parameters to be fitted.
#cp seriesfitbasal2.py seriesfitbasal3.py: use total cAMP level 0.1uM (Alberts et al. Molecular biology of the cell?)
# cp seriesfitbasal.py seriesfitbasal2.py: allow changing the target cAMP concentration
import matplotlib
matplotlib.use('Agg')
from pylab import *
import scipy.io
import os
import time
from os.path import exists
import re
import emoo
import random
from scipy.optimize import curve_fit
import mytools
import pickle

fixedBlock = 'Leak,AC5'
fixedBlockCoeff = '1.33,0.0'
fixedAltered = '1'
fixedAlteredCoeff = '1'

VARIABLES = [['k[0]',0.0,1000.0],                   #Gaba explusion rate
             ['k[1,9]',0.0001,1000.0],              #Gaba binding to GABABR. Since the Kd is known (110 nM), k[2,10] follow
             ['k[3,4,5,6,7,8]',0, 1000.0],          #Overall Gibg production rate (Gi binding+unbinding (to both Gaba-bound and -nonbound), Gia+Gibg dissolution)
             ['k[15,17,19,21,23]',0, 100.0],       #Gibg binding rate to GIRK and VGCC
             ['k[16,18,20,22,24]',0, 100.0],       #Gibg unbinding rate from GIRK and VGCC
             ['OnlyExp0_RGS',0.2,5.0],             #RGS concentration postsynaptically
             ['OnlyExp1_RGS',0.2,5.0],             #RGS concentration presynaptically, glutamatergic synapse that is activated at the same time as a nearby gabaergic synapse
             ['OnlyExp2_RGS',0.2,5.0],             #RGS concentration presynaptically, gabaergic synapse
             ['gaba_flux',0.,3000.0],              #Number of gaba particles/ms
            ]

mesh_input_file = open('mesh_general.out','r')
mesh_firstline = mesh_input_file.readline()
mesh_secondline = mesh_input_file.readline()
mesh_values = mesh_secondline.split()
my_volume = float(mesh_values[-2])*1e-15 #litres
mesh_input_file.close()

Nsamp = 2000

#These are fixed across experiments:
Duration         = 1e6+10000
tolerance        = 1e-7
gaba_input_onset   = 1e6
gaba_input_N       = 1
gaba_input_freq    = 1
gaba_input_dur     = 50.0 #Long duration due to assumed volume transmission
gaba_input_flux    = 5000.0 # Onset time constant (post-syn) should depend on concentration. This maybe corresponds to a saturating concentration
Ntrains          = 1
trainT           = 1
record_T         = 5

OBJECTIVES_DESCRIPTIONS = ['Onset and offset, post-synaptic', 'Onset and offset, pre-synaptic', 'Amplitudes']
OBJECTIVES = ['f'+str(i) for i in range(0,len(OBJECTIVES_DESCRIPTIONS))]

#Data for PPR
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_ydata_scaled = [data_exppr_scaled[i][1] for i in range(0,len(data_exppr_scaled))]
data_exppr_xdata = [data_exppr_scaled[i][0] for i in range(0,len(data_exppr_scaled))]
data_antiexppr_ydata_scaled = [1-data_exppr_ydata_scaled[i] for i in range(0,len(data_exppr_ydata_scaled))]                                                              #1-PPR(control)/PPR(GABAB_blocked) should closely reflect the temporal dynamics of GABAB activation

data_inppr_xdata = [350, 550, 750, 1000, 2000, 3000, 4000]
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_inppr_ydata_scaled = [data_inppr_ydata_nonscaled[i]/data_inppr_ydata_GABABblocked[i] for i in range(0,len(data_inppr_ydata_nonscaled))]                                 #Rationale: Certain fraction of PPR is caused by GABAB and the rest by vesicle depletion etc.
data_antiinppr_ydata_scaled = [1-data_inppr_ydata_scaled[i] for i in range(0,len(data_inppr_ydata_scaled))]                                                                  #1-PPR(control)/PPR(GABAB_blocked) should closely reflect the temporal dynamics of GABAB activation
#TODO: add these data to run_model and func_to_optimize

data_postsyn_scaled = [[77.25316584201036, 0.12989552390326564],
                       [97.62754955662777, 0.27964900975021945],
                       [120.74912053813799, 0.5282774023543535],
                       [138.750612423254, 0.7288741564007608],
                       [156.63910945934614, 0.8842729708376185],
                       [196.90764878022938, 0.9916886991909393],
                       [234.23125984842804, 0.9211331264703673],
                       [271.3853786430907, 0.7827806443354708],
                       [311.0253841162788, 0.6387828336106707],
                       [355.6724738368372, 0.4976187218340474],
                       [392.88309005601195, 0.38186520950392605],
                       [432.6007794879039, 0.26894098226069135],
                       [469.88907966578245, 0.18426105341213533],
                       [577.1282535675034, 0.07993061410052119],
                       [679.5227733173257, 0.03773851402945802],
                       [731.8817614837636, -0.018666219395389358]]
data_postsyn_xdata = [data_postsyn_scaled[i][0] for i in range(0,len(data_postsyn_scaled))]
data_postsyn_ydata_scaled = [data_postsyn_scaled[i][1] for i in range(0,len(data_postsyn_scaled))]

nrnfactor = 6.022e23*my_volume*1e-9 #how many molecules is 1nM in the given volume (0.5 fl)
experiments_blockeds = [ [['VGCC', 0.0]], [['GIRK', 0.0]], [['GIRK', 0.0]] ] #Three experiment types, where in the first one the only concentration change is setting VGCC to zero, and the second and third one where the only change is setting GIRK to zero
experiments_altereds = [ ['', ''], ['', ''], ['', ''] ]               #All rates same in the two experiments

MAXERR = 1e8 # Maximum error for missing data

def run_model(parameters,deleteFiles=True,rankID=0):
  data = []
  gaba_flux = gaba_input_flux
  thisAltered = ''
  thisAlteredCoeff = ''
  thisBlock = ''
  thisBlockCoeff = ''
  paramkeys = list(parameters.keys())
  experimentwiseBlockeds = []
  for iparam in range(0,len(paramkeys)):
    if paramkeys[iparam].find('k[') > -1:
      iks = paramkeys[iparam][2:-1].split(',')
      thisAltered = thisAltered + ','.join(iks)+','
      thisAlteredCoeff = thisAlteredCoeff + ','.join([str(parameters[paramkeys[iparam]]) for i in iks])+','
      if paramkeys[iparam].find('k[1,9]') > -1: #since we know the affinity of gaba + GABABR <-> gabaGABABR, use it to calculate the backward rate
        thisAltered = thisAltered + '2,10,'
        thisAlteredCoeff = thisAlteredCoeff + str(5.555*parameters[paramkeys[iparam]]*0.00011/0.005)+','+str(5.555*parameters[paramkeys[iparam]]*0.00011/0.005)+','
    elif paramkeys[iparam].find('gaba_flux') > -1:
      gaba_flux = parameters[paramkeys[iparam]]
    elif 'OnlyExp' in paramkeys[iparam]:
      blockediexp = int(paramkeys[iparam][7])
      blockeds = paramkeys[iparam][9:].split(',')
      experimentwiseBlockeds.append([blockediexp,paramkeys[iparam][9:],str(parameters[paramkeys[iparam]])])
    else:
      blockeds = paramkeys[iparam].split(',')
      thisBlock = thisBlock +','.join(blockeds)+','
      thisBlockCoeff = thisBlockCoeff + ','.join([str(parameters[paramkeys[iparam]]) for i in blockeds])+','

  if len(thisAltered) > 0 and thisAltered[-1] == ',':
    thisAltered = thisAltered[0:-1]
    thisAlteredCoeff = thisAlteredCoeff[0:-1]
  #if len(thisBlock) > 0 and thisBlock[0] == ',':
  #  thisBlock = thisBlock[1:]
  #  thisBlockCoeff = thisBlockCoeff[1:]
  
  DataAll = []
  filenames = []
  timenow = time.time()
  timetmp = time.time()

  nplotted_experiment = 0
  for iexperiment_type in range(0,3):
    Blocked = thisBlock+','.join([experiments_blockeds[iexperiment_type][i][0] for i in range(0,len(experiments_blockeds[iexperiment_type]))])
    BlockedCoeff = thisBlockCoeff+','.join([str(experiments_blockeds[iexperiment_type][i][1]) for i in range(0,len(experiments_blockeds[iexperiment_type]))])
    Altered = thisAltered + experiments_altereds[iexperiment_type][0]
    AlteredCoeff = thisAlteredCoeff + experiments_altereds[iexperiment_type][1]
    for iexperimentwiseBlock in range(0,len(experimentwiseBlockeds)):
      if experimentwiseBlockeds[iexperimentwiseBlock][0] == iexperiment_type:
        blockeds = experimentwiseBlockeds[iexperimentwiseBlock][1].split(',')
        Blocked = Blocked + ','+experimentwiseBlockeds[iexperimentwiseBlock][1]
        BlockedCoeff = BlockedCoeff + ','+','.join([str(experimentwiseBlockeds[iexperimentwiseBlock][2]) for x in blockeds])
    if len(Altered) > 0 and Altered[0] == ',':
      Altered = Altered[1:]
      AlteredCoeff = AlteredCoeff[1:]

    randomString = ''.join(random.choice('ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789') for _ in range(0,10))+'_'+str(rankID)
    print('python3 model_nrn_extfilename.py '+str(Duration)+' '+str(tolerance)+' '+str(gaba_input_onset)+' '+str(gaba_input_N)+' '+str(gaba_input_freq)+' '+str(gaba_input_dur)+' '+
       str(gaba_flux)+' '+str(Ntrains)+' '+str(trainT)+' None '+Blocked+' '+BlockedCoeff+' '+Altered+' '+AlteredCoeff+' fit'+randomString+'.mat '+str(record_T))
    os.system('python3 model_nrn_extfilename.py '+str(Duration)+' '+str(tolerance)+' '+str(gaba_input_onset)+' '+str(gaba_input_N)+' '+str(gaba_input_freq)+' '+str(gaba_input_dur)+' '+
         str(gaba_flux)+' '+str(Ntrains)+' '+str(trainT)+' None '+Blocked+' '+BlockedCoeff+' '+Altered+' '+AlteredCoeff+' fit'+randomString+'.mat '+str(record_T))
    print('Exp.  '+str(iexperiment_type)+', ID='+str(rankID)+' done in '+str(time.time()-timenow)+' sec')

    timetmp = time.time()
    filename = 'fit'+randomString+'.mat'
    filenames.append(filename)
    if not exists(filename):
      print('Error: filename = '+filename+' does not exists, Exp. ='+str(iexperiment))
      DataAll.append([])
      continue

    A = scipy.io.loadmat(filename)
    DataAll.append(A['DATA'])
    times = A['DATA'][0]
    headers = A['headers']    
  if deleteFiles:
    for filename in filenames:
      os.system('rm '+filename)
  return [DataAll,headers]



# Parameters:
# N: size of population
# C: size of capacity 
N_samples = Nsamp
C_samples = 2*Nsamp
N_generations = 25

def plotmybox(ax,ys,x=0,w=0.5,lw=0.5):
  ax.plot([x-w,x+w,x,x,x-w,x-w,x+w,x+w,x,nan,x-w,x+w,nan,x,x,x-w,x+w],[ys[0],ys[0],ys[0],ys[1],ys[1],ys[3],ys[3],ys[1],ys[1],nan,ys[2],ys[2],nan,ys[3],ys[4],ys[4],ys[4]],'k-',lw=lw)
def plotmybox2(ax,ys,x=0,w=0.5,w2=0.7,lw=0.5):
  ax.plot([x-w,x+w,x,x,x-w,x-w,x+w,x+w,x,nan,x-w2,x+w2,nan,x,x,x-w,x+w],[ys[0],ys[0],ys[0],ys[1],ys[1],ys[3],ys[3],ys[1],ys[1],nan,ys[2],ys[2],nan,ys[3],ys[4],ys[4],ys[4]],'k-',lw=lw)

#fvals = [[sum([y[j] for j in objective_groups_with_midpoints[i]]) for i in range(0,len(objective_groups_with_midpoints))] for y in params_f]

f,axarr = subplots(5,1)
for iax in range(0,2):
  axarr[iax].set_position([0.1+0.48*iax,0.75,0.4,0.19])
for iax in range(0,3):
  axarr[2+iax].set_position([0.14,0.5-0.21*iax,0.84,0.16])
for iax in range(0,len(axarr)):
  axarr[iax].tick_params(axis='both', which='major', labelsize=4)


fvals_all = []
params_all = []
medians_all = []
mins_all = []
maxs_all = []
prc25s_all = []
prc75s_all = []

for igen in range(0,25):
 fvals_thisgen = []
 params_thisgen = []
 for myseed in range(1,16):
  filename = 'fitfiles/fitc_seed'+str(myseed)+'_N'+str(N_samples)
  if exists(filename+"_tmp"+str(igen)+".sav"):
    unpicklefile = open(filename+"_tmp"+str(igen)+".sav", 'rb')
  else:
    print(filename+"_tmp"+str(igen)+".sav not found")
    continue
  unpickledlist = pickle.load(unpicklefile,encoding='bytes')
  unpicklefile.close()
  fvals = [[unpickledlist[0][j][i] for i in range(len(VARIABLES),len(VARIABLES)+2)] for j in range(0,N_samples)]
  params = [[unpickledlist[0][j][i] for i in range(0,len(VARIABLES)+3)] for j in range(0,N_samples)]
  fvals_thisgen = fvals_thisgen + fvals
  params_thisgen = params_thisgen + params
 fvals_all.append(fvals_thisgen[:])
 params_all.append(params_thisgen[:])
 medians_all.append([median([fvals[i][ifval] for i in range(0,N_samples)]) for ifval in range(0,2)])
 mins_all.append([min([fvals[i][ifval] for i in range(0,N_samples)]) for ifval in range(0,2)])
 maxs_all.append([max([fvals[i][ifval] for i in range(0,N_samples)]) for ifval in range(0,2)])
 prc25s_all.append([percentile([fvals[i][ifval] for i in range(0,N_samples)],25.) for ifval in range(0,2)])
 prc75s_all.append([percentile([fvals[i][ifval] for i in range(0,N_samples)],75.) for ifval in range(0,2)])
 print(str(unpickledlist[1]))
 print(str(len(fvals)))
 sumvals = [fvals_all[-1][i][0]/medians_all[-1][0] + fvals_all[-1][i][1]/medians_all[-1][1] for i in range(0,len(fvals_all[-1]))]
 #sumvals = [fvals_all[-1][i][0]/medians_all[-1][0] + 0*fvals_all[-1][i][1]/medians_all[-1][1] for i in range(0,len(fvals_all[-1]))]
 ibest = argmin(sumvals)
 print('best params: '+str(params_all[-1][ibest])+", best_f0="+str(fvals_all[-1][ibest][0])+", min_f0="+str(min([fvals_all[-1][i][0] for i in range(0,len(fvals_all[-1]))])))
 #ichosen = [i for i in range(0,len(sumvals)) if abs(params_all[-1][i][len(VARIABLES)]-0.4840479) < 0.00001 and abs(params_all[-1][i][len(VARIABLES)+1]-0.8611440) < 0.00001]
 #ichosen = [i for i in range(0,len(sumvals)) if abs(params_all[-1][i][len(VARIABLES)]-0.428515837) < 0.00001 and abs(params_all[-1][i][len(VARIABLES)+1]-1.06957820) < 0.00001]
 #vars_chosen = [106.712,760.671,100.983,8.947,47.145,3.772,1.699,2.329,2693.683] #p1, see choosefinals.sh
 #vars_chosen = [91.805,846.811,67.856,9.858,71.217,5.000,1.817,2.170,3000.000] #p1, see choosefinals.sh
 vars_chosen = [107.034,514.651,167.633,11.026,66.931,3.009,1.904,2.241,1507.740] #p0, see choosefinals.sh
 ichosen = [i for i in range(0,len(sumvals)) if sum([abs(params_all[-1][i][j]-vars_chosen[j]) for j in range(0,len(vars_chosen))]) < 0.02]
 if len(ichosen) > 0:
  for ifval in range(0,2):
   axarr[ifval].plot(igen+1,params_all[-1][ichosen[0]][len(VARIABLES)+ifval],'k.',zorder=4)
  print("Found ichosen="+str(ichosen)+" for igen="+str(igen)+", params_all[-1][ichosen[0]][len(VARIABLES):len(VARIABLES)+2] = "+str(params_all[-1][ichosen[0]][len(VARIABLES):len(VARIABLES)+2]))
 else:
  print("No ichosen found for igen="+str(igen)+", min diff = "+str(min([sum([abs(params_all[-1][i][j]-vars_chosen[j]) for j in range(0,len(vars_chosen))]) for i in range(0,len(sumvals))])))
for ifval in range(0,2):
 for igen in range(0,len(fvals_all)):
   plotmybox2(axarr[ifval],[mins_all[igen][ifval],prc25s_all[igen][ifval],medians_all[igen][ifval],prc75s_all[igen][ifval],maxs_all[igen][ifval]],igen+1,0.28,0.35)
 axarr[ifval].set_xlim([0.5,len(fvals_all)+0.5])
 if ifval < 5:
   axarr[ifval].set_yscale("log", nonpositive='clip')
   axarr[ifval].set_xticks([5,10,15,20,25])

 if len(fvals_all) < 20:
   pos = axarr[ifval].get_position()
   axarr[ifval].set_position([pos.x0, pos.y0, pos.width*len(fvals_all)/20, pos.height])

 axarr[ifval].tick_params(axis='both', which='major', labelsize=7)

params_all_all = []
for iparams in range(0,len(params_all)):
  params_all_all = params_all_all + params_all[iparams]

minNow = inf
params_all_chosen = []
for myseed in range(1,16):
  filename = 'fitc_seed'+str(myseed)+'_N'+str(N_samples)
  gensdone = -1
  coeffs = rand(len(VARIABLES),)
  for gen in range(25,-1,-1):
    if exists('fitfiles/'+filename+'_tmp'+str(gen)+'.sav'):
      gensdone = gen
      unpicklefile = open('fitfiles/'+filename+'_tmp'+str(gen)+'.sav', 'rb')
      unpickledlist = pickle.load(unpicklefile,encoding='bytes')
      unpicklefile.close()
      params_all = unpickledlist[0]
      columns = unpickledlist[1]
      break
    else:
      print('file '+'fitfiles/'+filename+'_tmp'+str(gen)+'.sav not found')
      
  if gensdone == -1:
    continue

  param_fdims = range(len(VARIABLES), len(VARIABLES)+len(OBJECTIVES))
  fs_all = [params_all[i][param_fdims[0]:param_fdims[0]+3] for i in range(0,len(params_all))]
  params_all_0to1 = params_all+0
  for ifdim in range(0,len(OBJECTIVES)):
    params_all_0to1[:,param_fdims[ifdim]] = (params_all[:,param_fdims[ifdim]] - min(params_all[:,param_fdims[ifdim]]))/(max(params_all[:,param_fdims[ifdim]])-min(params_all[:,param_fdims[ifdim]]))
  medians = [median([y[i] for y in params_all_0to1]) for i in param_fdims]
  params_f =  [[1.0*y[param_fdims[j]]/medians[j] for j in range(0,len(param_fdims))] for y in params_all_0to1]

  fvals = [y[0]+y[1] for y in params_f]
  ibest = argmin(fvals)

  print('f0-1best = '+str(params_all[ibest,param_fdims[0:2]]))
  print('params_best = '+','.join(['{:.3f}'.format(params_all[ibest,j]) for j in range(0,len(VARIABLES))]))
  params_all_chosen.append(params_all[ibest,0:param_fdims[0]])

params_all = params_all_chosen[0]

paramdict = {}
for iparam in range(0,len(VARIABLES)):
  paramdict[VARIABLES[iparam][0]] = params_all[iparam]

print('paramdict = '+str(paramdict))
A = run_model(paramdict,True,0)

DataAll = A[0]
headers = A[1]
timesAll = [DataAll[i][0] for i in range(0,len(DataAll))]
VGCCboundsAll = [0*x for x in timesAll]
for i in range(0,len(headers)):
  if headers[i].find(' ') != -1:
    headers[i] = headers[i][0:headers[i].find(' ')]
  if 'VGCCGibg' in headers[i]:
    VGCCboundsAll = [DataAll[iexp][i] for iexp in range(0,len(timesAll))]
iGIRK = [i for i in range(0,len(headers)) if 'GIRKGibg' in headers[i]]
nGibgs = [1 if headers[i][-1] == 'g' else int(headers[i][-1]) for i in iGIRK]
coeffs = [0.0, 0.01, 0.06, 0.26, 1.0]
GIRKcondsAll = [0*x for x in timesAll]
GIRK1boundsAll = [0*x for x in timesAll]
for iiGIRK in range(0,len(iGIRK)):
  for iexp in range(0,len(timesAll)):
    GIRKcondsAll[iexp] = GIRKcondsAll[iexp] + coeffs[nGibgs[iiGIRK]]*DataAll[iexp][iGIRK[iiGIRK]]
    if iiGIRK == 0:
      GIRK1boundsAll[iexp] = GIRK1boundsAll[iexp] + DataAll[iexp][iGIRK[0]]

mydict = {}

amps = [max(GIRK1boundsAll[0]), max(VGCCboundsAll[1]), max(VGCCboundsAll[2])]

GIRKcond_interpolated = mytools.interpolate(timesAll[0],GIRKcondsAll[0]/max(GIRKcondsAll[0]),[1e6+x for x in data_postsyn_xdata])
VGCCbound_ex_interpolated = mytools.interpolate(timesAll[1],VGCCboundsAll[1]/max(VGCCboundsAll[1]),[1e6+x for x in data_exppr_xdata])
VGCCbound_in_interpolated = mytools.interpolate(timesAll[2],VGCCboundsAll[2]/max(VGCCboundsAll[2]),[1e6+x for x in data_inppr_xdata])
mydict['f0'] = sum([abs(GIRKcond_interpolated[i]/max(GIRKcondsAll[0]) - data_postsyn_ydata_scaled[i]) for i in range(0,len(data_postsyn_ydata_scaled))])
mydict['f1'] = sum([abs(VGCCbound_ex_interpolated[i] - data_antiexppr_ydata_scaled[i]/max(data_antiexppr_ydata_scaled)) for i in range(0,len(data_antiexppr_ydata_scaled))]) # data_antiexppr_ydata_scaled is re-scaled because the maximum of the PPR need not be 1 at the maximum point
mydict['f1'] = mydict['f1'] + sum([abs(VGCCbound_in_interpolated[i] - data_antiinppr_ydata_scaled[i]/max(data_antiinppr_ydata_scaled)) for i in range(0,len(data_antiinppr_ydata_scaled))]) # data_antiinppr_ydata_scaled is re-scaled because the maximum of the PPR need not be 1 at the maximum point
mydict['f2'] = -amps[0]*amps[1]*amps[2]
for i in [0,1,2]:
  if isnan(mydict['f'+str(i)]):
    mydict['f'+str(i)] = MAXERR
axarr[2+0].plot(timesAll[0]-1e6,GIRKcondsAll[0]/max(GIRKcondsAll[0]),'k-')
#axarr[2+0].plot(data_postsyn_xdata,GIRKcond_interpolated,'b.')
axarr[2+0].plot(data_postsyn_xdata,data_postsyn_ydata_scaled,'k.')
axarr[2+1].plot(timesAll[1]-1e6,VGCCboundsAll[1]/max(VGCCboundsAll[1]),'k-')
#axarr[2+1].plot(data_exppr_xdata,VGCCbound_ex_interpolated,'k.')
axarr[2+1].plot(data_exppr_xdata,[x*max(VGCCboundsAll[1]/max(VGCCboundsAll[1]))/max(data_antiexppr_ydata_scaled) for x in data_antiexppr_ydata_scaled],'k.')
axarr[2+2].plot(timesAll[2]-1e6,VGCCboundsAll[2]/max(VGCCboundsAll[2]),'k-')
#axarr[2+2].plot(data_inppr_xdata,VGCCbound_in_interpolated,'k.')
axarr[2+2].plot(data_inppr_xdata,[x*max(VGCCboundsAll[2]/max(VGCCboundsAll[2]))/max(data_antiinppr_ydata_scaled) for x in data_antiinppr_ydata_scaled],'k.')
myparstr = str(paramdict)
myparstr2 = myparstr[:]
commapos = 0
ncommafound = 0
while myparstr2.find(', ') > -1:
  commapos = commapos + myparstr2.find(', ')+1
  ncommafound = ncommafound + 1
  myparstr2 = myparstr2[myparstr2.find(', ')+1:]
  if ncommafound == 5:
    break
if commapos > 0:
  myparstr = myparstr[0:commapos-1]+',\n'+myparstr[commapos:]
print('paramdict='+str(paramdict)+', mydict='+str(mydict))

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

axarr[0].set_ylabel('Error function 1',fontsize=7)
axarr[1].set_ylabel('Error function 2',fontsize=7)
axarr[0].set_xlabel("Generation",fontsize=7)
axarr[1].set_xlabel("Generation",fontsize=7)

axarr[2].set_ylabel('Normalized\n GIRK cond.', fontsize=7)
axarr[3].set_ylabel('Normalized\n VGCC binding\n(excitatory)', fontsize=7)
axarr[4].set_ylabel('Normalized\n VGCC binding\n(inhibitory)', fontsize=7)
axarr[4].set_xlabel('time (ms)',fontsize=7)

#axarr[0].set_xlim([
axarr[2].set_xlim([0,4500])
axarr[3].set_xlim([0,4500])
axarr[4].set_xlim([0,4500])

f.savefig('fig1.pdf')