#cp simsteadystate_someonly.py simsteadystate_li2020.py
#simsteadystate.py: Simulate the Ca binding and protein activation time courses for steady Ca inputs with varying Ca input fluxes
#Input: the amplitude of the fluxes of other ligands: beta-adrenergic, glutamate and ACh
#Tuomo Maki-Marttunen, 2019-2020
#CC BY 4.0
import matplotlib
matplotlib.use('Agg')
import numpy as np
import emoo
from pylab import *
import pickle
import mytools
import time
import random
import sys
import scipy.io
from os.path import exists
import os
import calcconds

globalcounter = 0
             
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()

conds_hom1 = [12.4, 18.9]
conds_hom2 = 2.2
conds_het = 2.5

for iargv in range(0,len(sys.argv)):
  print "argv["+str(iargv)+"] = "+sys.argv[iargv]

Duration = 4940000
tolerance = 1e-6
Ca_input_onset   = 40000.0
Ca_input_N       = 1
Ca_input_freq    = 1.0
Ca_input_dur     = 300000.0

Ca_input_fluxes = [0.000025,0.00005,0.0001,0.00025,0.0005,0.001,0.0015,0.002]+[0.025*i for i in range(1,10)]+[0.25,0.3,0.35,0.4,0.45]+[0.25*i for i in range(2,10)] +[2.5,3.0,3.5,4.0,4.5]+ [2.5*i for i in range(2,10)]
Ntrains          = 1
cols = ['#0000FF','#FF0000','#00BB00']
trainT = 100000.0
nrnfactor = 6.022e23*my_volume*1e-9 #how many molecules is 1nM in the given volume (0.5 fl)

MAXERR = 1e8 # Maximum error for missing data

fb_coeff = 1.0
CaM_tot = 5.0
calbin_coeff = 1.0

if len(sys.argv) > 1:
  fb_coeff = float(sys.argv[1])
if len(sys.argv) > 2:
  CaM_tot = float(sys.argv[2])
if len(sys.argv) > 3:
  calbin_coeff = float(sys.argv[3])


def run_model(parameters,saveFig="",deleteFiles=False,rankID=0):
  data = []
  
  paramkeys = parameters.keys()
  L_input_flux   = 0.0
  Glu_input_flux = 0.0
  ACh_input_flux = 0.0


  species_alts = [['CaOut', 1.0],
                  ['Leak', 0.0],
                  ['Calbin', calbin_coeff],
                  ['CalbinC', calbin_coeff],
                  ['LOut', 0.0],
                  ['Epac1', 0.0],
                  ['PMCA', 0.0],
                  ['NCX', 0.0],
                  ['L', 0.0],
                  ['R', 0.0],
                  ['Gs', 0.0],
                  ['Gi', 0.0],
                  ['AC1', 0.0],
                  ['ATP', 0.0],
                  ['AC8', 0.0],
                  ['PDE1', 0.0],
                  ['AMP', 0.0],
                  ['Ng', 2.5],            # originally 20 uM
                  ['CaM', CaM_tot/60.0],  # originally 60 uM
                  ['PP2B', 21.739],       # originally 2.3 uM
                  ['CK', 2.1739],         # originally 23 uM
                  ['PKA', 0.0],
                  ['I1', 0.0],
                  ['PP1', 0.0],
                  ['GluR1', 0.0],
                  ['GluR1_memb', 0.0],
                  ['PDE4', 0.0],
                  ['fixedbuffer', fb_coeff],
                  ['MGluR', 0.0],
                  ['GluOut', 0.0],
                  ['Gqabg', 0.0],
                  ['PLC', 0.0],
                  ['Pip2', 0.0],
                  ['PIkinase', 0.0],
                  ['Ip3degPIk', 0.0],
                  ['PKC', 0.0],
                  ['DAG', 0.0],
                  ['DAGK', 0.0],
                  ['DGL', 0.0],
                  ['CaDGL', 0.0],
                  ['DAGCaDGL', 0.0],
                  ['Ip3degrad', 0.0],
                  ['GluR2', 0.0],
                  ['GluR2_memb', 0.0],
                  ['PP2A', 0.0],
                  ['M1R', 0.0],
                  ['PLA2', 0.0]]

  addition_IC = ''
  addition_IC_values = ''
  for ispecalt in range(0,len(species_alts)):
    if species_alts[ispecalt][1] != 1.0:
      addition_IC = addition_IC + species_alts[ispecalt][0] + '-'
      addition_IC_values = addition_IC_values + str(species_alts[ispecalt][1]) + '-'
  addition_IC = addition_IC[0:-1]
  addition_IC_values = addition_IC_values[0:-1]

  timesAll = []
  timeCoursesAll = []
  maxValsAll = []
  if len(saveFig) > 0:
    close("all")
    f,axarr = subplots(len(Ca_input_fluxes),2)
  filenames = []
  timenow = time.time()
  DATA_all_all_all = []
  for iexperiment in range(0,len(Ca_input_fluxes)):
    myString = 'nrn_oldCaM_li2020_'+str(fb_coeff)+'_'+str(CaM_tot)+'_'+str(calbin_coeff)+'_caflux'+str(iexperiment)

    filename = myString
    print "thisfile = "+filename
    filenames.append(filename)
    #print 'loading filename = '+filename
    if not exists(filename+'.mat'):
      print('python model_nrn_oldCaM_altered_noU_extfilename_lowmem_recall.py '+str(Duration)+' '+str(tolerance)+' '+str(Ca_input_onset)+' '+str(Ca_input_N)+' '+str(Ca_input_freq)+' '+str(Ca_input_dur)+' '+
             str(Ca_input_fluxes[iexperiment])+' '+str(L_input_flux)+' '+str(Glu_input_flux)+' '+str(ACh_input_flux)+' '+
             str(Ntrains)+' '+str(trainT)+' None '+myString+' '+addition_IC+' '+addition_IC_values)
      os.system('python model_nrn_oldCaM_altered_noU_extfilename_lowmem_recall.py '+str(Duration)+' '+str(tolerance)+' '+str(Ca_input_onset)+' '+str(Ca_input_N)+' '+str(Ca_input_freq)+' '+str(Ca_input_dur)+' '+
             str(Ca_input_fluxes[iexperiment])+' '+str(L_input_flux)+' '+str(Glu_input_flux)+' '+str(ACh_input_flux)+' '+
             str(Ntrains)+' '+str(trainT)+' None '+myString+' '+addition_IC+' '+addition_IC_values)
      print('Exp. '+str(iexperiment)+', ID='+str(rankID)+' done in '+str(time.time()-timenow)+' sec')
      if not exists(filename+'.mat'):
        print 'Error: filename = '+filename+'.mat does not exists, Exp. ='+str(iexperiment)
        timesAll.append([])
        timeCoursesAll.append([])
        maxValsAll.append([])
        DATA_all_all_all.append([])
        continue
    
    conds, times = calcconds.calcconds_nrn(filename+'.mat')

    DATA_all_all = scipy.io.loadmat(filename+'.mat')
    if DATA_all_all['maxDATA'].shape[0] != 1:
      DATA_all_all['maxDATA'] = DATA_all_all['maxDATA'].T
    DATA_all = {}
    maxDATA_all = {}
    for ikey in range(0,len(DATA_all_all['headers'])):
      first_space = DATA_all_all['headers'][ikey].find(' ')
      mykey = DATA_all_all['headers'][ikey]
      if first_space > -1:
        mykey = DATA_all_all['headers'][ikey][0:first_space]
      DATA_all[mykey] = DATA_all_all['DATA'][ikey]
      maxDATA_all[mykey] = DATA_all_all['maxDATA'][0][ikey]

    itime = argmin(abs(times - 4339000))
    timesAll.append(times[itime])
    timeCoursesAll.append(conds[itime])
    maxValsAll.append(maxDATA_all['Ca'])
    DATA_all_all_all.append([DATA_all_all['DATA'][ikey][itime] for ikey in range(0,len(DATA_all_all['DATA']))])
  if deleteFiles:
    for filename in filenames:
      os.system('rm '+filename+'.mat')
      #print('rm '+filename+'.mat')
  return [timesAll, timeCoursesAll, maxValsAll,  DATA_all_all_all, DATA_all_all['headers']]

if True:
  A = run_model({}, "", False, 0)
  picklelist = [A, Ca_input_fluxes]
  file=open('steadystate_new_oldCaM_li2020_'+str(fb_coeff)+'_'+str(CaM_tot)+'_'+str(calbin_coeff)+'_processed.sav', 'w')
  pickle.dump(picklelist,file)
  file.close()

