/*
Created by BluePyOpt(1.6.42) at 2019-03-05 09:07:19.588315
Author: Elisabetta Iavarone @ Blue Brain Project
*/
{load_file("stdrun.hoc")}
{load_file("import3d.hoc")}
/*
 * Check that global parameters are the same as with the optimization
 */
proc check_parameter(/* name, expected_value, value */){
  strdef error
  if($2 != $3){
    sprint(error, "Parameter %s has different value %f != %f", $s1, $2, $3)
    execerror(error)
  }
}
proc check_simulator() {
  check_parameter("celsius", 34, celsius)
  check_parameter("v_init", -79, v_init)
}

begintemplate cAD_ltb
  public init, morphology, geom_nseg_fixed, geom_nsec, gid
  public channel_seed, channel_seed_set
  public soma, dend, apic, axon, myelin
  create soma[1], dend[1], apic[1], axon[1], myelin[1]

  objref this, CellRef, segCounts

  public all, somatic, apical, axonal, basal, myelinated, APC
  objref all, somatic, apical, axonal, basal, myelinated, APC

proc init(/* args: morphology_dir, morphology_name */) {
  all = new SectionList()
  apical = new SectionList()
  axonal = new SectionList()
  basal = new SectionList()
  somatic = new SectionList()
  myelinated = new SectionList()

  //gid in this case is only used for rng seeding
  gid = 0

  //For compatibility with BBP CCells
  CellRef = this

  forall delete_section()

  if(numarg() >= 2) {
    load_morphology($s1, $s2)
  } else {
    load_morphology($s1, "jy160728_A_idA.asc")
  }

  geom_nseg()
    replace_axon()
  insertChannel()
  biophys()

  // Initialize channel_seed_set to avoid accidents
  channel_seed_set = 0
  // Initialize random number generators
  re_init_rng()
}

proc load_morphology(/* morphology_dir, morphology_name */) {localobj morph, import, sf, extension
  strdef morph_path
  sprint(morph_path, "%s/%s", $s1, $s2)

  sf = new StringFunctions()
  extension = new String()

  sscanf(morph_path, "%s", extension.s)
  sf.right(extension.s, sf.len(extension.s)-4)

  if( strcmp(extension.s, ".asc") == 0 ) {
    morph = new Import3d_Neurolucida3()
  } else if( strcmp(extension.s, ".swc" ) == 0) {
    morph = new Import3d_SWC_read()
  } else {
    printf("Unsupported file format: Morphology file has to end with .asc or .swc" )
    quit()
  }

  morph.quiet = 1
  morph.input(morph_path)

  import = new Import3d_GUI(morph, 0)
  import.instantiate(this)
}

/*
 * Assignment of mechanism values based on distance from the soma
 * Matches the BluePyOpt method
 */
proc distribute_distance(){local x localobj sl
  strdef stmp, distfunc, mech

  sl = $o1
  mech = $s2
  distfunc = $s3
  this.soma[0] distance(0, 0.5)
  sprint(distfunc, "%%s %s(%%f) = %s", mech, distfunc)
  forsec sl for(x, 0) {
    sprint(stmp, distfunc, secname(), x, distance(x))
    execute(stmp)
  }
}

proc geom_nseg() {
  this.geom_nsec() //To count all sections
  //TODO: geom_nseg_fixed depends on segCounts which is calculated by
  //  geom_nsec.  Can this be collapsed?
  this.geom_nseg_fixed(40)
  this.geom_nsec() //To count all sections
}

proc insertChannel() {
  forsec this.all {
    insert pas
    insert TC_cad
  }
  forsec this.apical {
    insert TC_ih_Bud97
    insert TC_Nap_Et2
    insert TC_iA
    insert TC_iL
    insert SK_E2
    insert TC_HH
    insert TC_iT_Des98
  }
  forsec this.axonal {
    insert TC_HH
  }
  forsec this.basal {
    insert TC_ih_Bud97
    insert TC_Nap_Et2
    insert TC_iA
    insert TC_iL
    insert SK_E2
    insert TC_HH
    insert TC_iT_Des98
  }
  forsec this.somatic {
    insert TC_ih_Bud97
    insert TC_Nap_Et2
    insert TC_iA
    insert TC_iL
    insert SK_E2
    insert TC_HH
    insert TC_iT_Des98
  }
  forsec this.myelinated {
  }
}

proc biophys() {
  
  forsec CellRef.all {
    cm = 1
    Ra = 100
    ena = 50
    ek = -90
    e_pas = -80
    g_pas = 3.2505679140951987e-05
  }
  
  forsec CellRef.apical {
    pcabar_TC_iT_Des98 = 8.9487559513902622e-05
    gh_max_TC_ih_Bud97 = 8.4228903849141138e-06
    gNap_Et2bar_TC_Nap_Et2 = 7.4653712158998865e-05
    gk_max_TC_iA = 0.0041818662141791902
    gk_max_TC_HH = 0.0099243833236075541
    gna_max_TC_HH = 0.0052508466616232969
    pcabar_TC_iL = 1.1024317581168538e-05
    gSK_E2bar_SK_E2 = 0.00021799609115766169
    taur_TC_cad = 9.5512406128698544
    gamma_TC_cad = 0.00053950618656117702
  }
  
  forsec CellRef.axonal {
    gk_max_TC_HH = 0.11490723847692205
    gna_max_TC_HH = 0.21874510090978222
  }
  
  forsec CellRef.basal {
    pcabar_TC_iT_Des98 = 8.9487559513902622e-05
    gh_max_TC_ih_Bud97 = 8.4228903849141138e-06
    gNap_Et2bar_TC_Nap_Et2 = 7.4653712158998865e-05
    gk_max_TC_iA = 0.0041818662141791902
    gk_max_TC_HH = 0.0099243833236075541
    gna_max_TC_HH = 0.0052508466616232969
    pcabar_TC_iL = 1.1024317581168538e-05
    gSK_E2bar_SK_E2 = 0.00021799609115766169
    taur_TC_cad = 9.5512406128698544
    gamma_TC_cad = 0.00053950618656117702
  }
  
  forsec CellRef.somatic {
    pcabar_TC_iT_Des98 = 8.9487559513902622e-05
    gh_max_TC_ih_Bud97 = 4.7504774088700974e-05
    gNap_Et2bar_TC_Nap_Et2 = 1.3389527019193626e-05
    gk_max_TC_iA = 0.063745893099437026
    gk_max_TC_HH = 0.11616929624507591
    gna_max_TC_HH = 0.091136959874091761
    pcabar_TC_iL = 0.00049811286552489979
    gSK_E2bar_SK_E2 = 0.0012832790164235054
    taur_TC_cad = 11.005585762426819
    gamma_TC_cad = 0.0067401226271938839
  }
  
  forsec CellRef.myelinated {
  }
  
}

func sec_count(/* SectionList */) { local nSec
  nSec = 0
  forsec $o1 {
      nSec += 1
  }
  return nSec
}

/*
 * Iterate over the section and compute how many segments should be allocate to
 * each.
 */
proc geom_nseg_fixed(/* chunkSize */) { local secIndex, chunkSize
  chunkSize = $1
  soma area(.5) // make sure diam reflects 3d points
  secIndex = 0
  forsec all {
    nseg = 1 + 2*int(L/chunkSize)
    segCounts.x[secIndex] = nseg
    secIndex += 1
  }
}

/*
 * Count up the number of sections
 */
proc geom_nsec() { local nSec
  nSecAll = sec_count(all)
  nSecSoma = sec_count(somatic)
  nSecApical = sec_count(apical)
  nSecBasal = sec_count(basal)
  nSecMyelinated = sec_count(myelinated)
  nSecAxonalOrig = nSecAxonal = sec_count(axonal)

  segCounts = new Vector()
  segCounts.resize(nSecAll)
  nSec = 0
  forsec all {
    segCounts.x[nSec] = nseg
    nSec += 1
  }
}

/*
 * Replace the axon built from the original morphology file with a stub axon
 */
    
proc replace_axon(){ local nSec, D1, D2
  // preserve the number of original axonal sections
  nSec = sec_count(axonal)

  // Try to grab info from original axon
  if(nSec == 0) { //No axon section present
    D1 = D2 = 1
  } else if(nSec == 1) {
    access axon[0]
    D1 = D2 = diam
  } else {
    access axon[0]
    D1 = diam
    access soma distance() //to calculate distance from soma
    forsec axonal{
      //if section is longer than 60um then store diam and exit from loop
      if(distance(0.5) > 60){
        D2 = diam
        break
      }
    }
  }

  // get rid of the old axon
  forsec axonal{
    delete_section()
  }

  create axon[2]

  access axon[0] {
    L = 30
    diam = D1
    nseg = 1 + 2*int(L/40)
    all.append()
    axonal.append()
  }
  access axon[1] {
    L = 30
    diam = D2
    nseg = 1 + 2*int(L/40)
    all.append()
    axonal.append()
  }
  nSecAxonal = 2
  soma[0] connect axon[0](0), 1
  axon[0] connect axon[1](0), 1
}
        



func hash_str() {localobj sf strdef right
  sf = new StringFunctions()

  right = $s1

  n_of_c = sf.len(right)

  hash = 0
  char_int = 0
  for i = 0, n_of_c - 1 {
     sscanf(right, "%c", & char_int)
     hash = (hash * 31 + char_int) % (2 ^ 31 - 1)
     sf.right(right, 1)
  }

  return hash
}

proc re_init_rng() {localobj sf
    strdef full_str, name

    sf = new StringFunctions()

    if(numarg() == 1) {
        // We received a third seed
        channel_seed = $1
        channel_seed_set = 1
    } else {
        channel_seed_set = 0
    }


}


endtemplate cAD_ltb