/*
Created by BluePyOpt(1.5.22) at 2017-08-25 10:55:58.013607
*/
{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", -65, v_init)
}

begintemplate CA1_PC
  public init, morphology, geom_nseg_fixed, geom_nsec, biophys,load_morphology
  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()

  //For compatibility with BBP CCells
  CellRef = this

  forall delete_section()

  if(numarg() >= 2) {
    load_morphology($s1, $s2)
  } else {
    load_morphology("morphologies", "dend-050921AM2_axon-051208AM2.asc")
  }
    dist()
    geom_nseg2()

    flag=1
    //flag=section_exists("a.*")
   if(flag==0){ 
        place_axon(1)
        geom_nseg2()
    }else{ 
      replace_axon()
   }
  insertChannel()
  biophys()
  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
}

external lambda_f

proc geom_nseg2() {
  forsec all { nseg = int((L/(0.1*lambda_f(100))+.9)/2)*2 + 1  }
}


proc dist(){
    this.soma[0] distance(0,0.5)
    }


proc insertChannel() {


  forsec this.all {
    insert pas

  }

  forsec this.axonal {
    insert kmdb
    insert nax
	insert kdr
	insert kap

  }

  forsec this.somatic {
    
	insert cad
    insert kmdb     
    insert hd
    insert caldb
    insert cat
    insert kca
    insert mykca     
    insert somacar
	insert can
    insert na3
    insert kdr 
	insert kap
	
	}
 

  forsec this.apical {
      insert car
      insert na3d
	  insert kdr   
	  insert kad
	  insert hd
	  insert caldb
	  insert cad
	  insert cat
	  insert can
	  insert kca
	  insert mykca
	  
    }
	
  forsec this.basal {
     insert car
     insert na3d
	 insert kdr 
     insert kap
     insert hd	  
	 insert caldb
	 insert cad
	 insert cat
	 insert can
	 insert kca
	 insert mykca
    }
}

proc biophys() {
  Rm=20000
   v_rest=-80
  
  forsec CellRef.all {
    cm = 1
 
  }
  
  forsec CellRef.apical {
  
   gcabar_car=0.00001
    gcanbar_can =0.0001
    gbar_mykca = 8.5e-5
    gbar_kca = 5e-5
    gcatbar_cat = 0.0001
    gcalbar_caldb = 0.0001
    gkdrbar_kdr =0.06
    gbar_na3d = 0.04
    Ra =600
    g_pas = 1/33000
	e_pas=v_rest-5
   
  }
  
  forsec CellRef.axonal {
    gbar_kmdb = 0.003
    gkabar_kap =0.02
    gbar_nax = 0.06
    gkdrbar_kdr = 0.3
    Ra = 150
    g_pas = 3e-5  
    e_pas =v_rest
   
  }


  forsec CellRef.basal {
    gcabar_car=0.00001
    gcanbar_can =0.0001
    gbar_mykca = 8.5e-5
    gbar_kca = 5e-5
    gcatbar_cat = 0.0001
    gcalbar_caldb = 0.0002
    ghdbar_hd = 0.00001
    gkabar_kap =0.015
    gkdrbar_kdr =0.06
    gbar_na3d = 0.05
    Ra = 400
    g_pas = 1/33000
    e_pas =v_rest
   }
  
  forsec CellRef.somatic {
	gkabar_kap =0.02
    gbar_kmdb = 0.001
	ghdbar_hd = 0.00001
	gcalbar_caldb = 0.0005
	gcanbar_can = 0.0001
    gcatbar_cat = 0.0001
    gbar_kca = 5e-5
    gbar_mykca = 8.5e-5
	gcabar_somacar =0.00001
	Ra = 250
    g_pas =3e-5
    e_pas =v_rest
    gkdrbar_kdr =0.1
    gbar_na3 = 0.05
  

  }
  
  forsec CellRef.myelinated {
  }
 
 distribute_distance(CellRef.apical, "ghdbar_hd", " (1+ 6/(1.0 + exp((280-%.17g)/50)))*1e-06")
 distribute_distance(CellRef.apical, "gkabar_kad", "(18./(1. + exp((160-%.17g)/40)))*0.015")

 
}

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
  }
}


  proc place_axon(){local L_chunk, L_target, chunkSize
   
     L_target = 60  // length of axon
     nseg0 = 5  // number of segments for each of the two axon sections

     nseg_total = nseg0 * 2
     chunkSize = L_target/nseg_total

            
     execute1("create axon[2]", CellRef)

         for i=0,1{
             access axon[i]
             L =  L_target/2
             nseg = nseg_total/2

             for (x) {
                 if (x > 0 && x < 1) {
                     if (i==0){ diam(x) = $1*2
                         }else { diam(x) = $1
                               }
                        }
             }

             all.append()
             axonal.append()
         }

         nSecAxonal = 2
         soma[0] connect axon[0](0), 1
         axon[0] connect axon[1](0), 1

         

 }



/*
 * Replace the axon built from the original morphology file with a stub axon
 */

      
    proc replace_axon(){local nSec, L_chunk, dist, i1, i2, count, L_target, chunkSize, L_real localobj diams, lens

     L_target = 60  // length of stub axon
     nseg0 = 5  // number of segments for each of the two axon sections

     nseg_total = nseg0 * 2
     chunkSize = L_target/nseg_total

     nSec = -1
     forsec axonal{nSec = nSec + 1}

     // Try to grab info from original axon
     if(nSec < 1){ //At least two axon sections have to be present!
            
        //execerror("Less than two axon sections are present! Add an axon to the morphology and try again!")
           access axon[0]
          ddiam=diam 
          forsec axonal{delete_section()}
          place_axon(ddiam)  

     } else {

         diams = new Vector()
         lens = new Vector()

         access axon[0]
         i1 = v(0.0001) // used when serializing sections prior to sim start

         access axon[1]
         i2 = v(0.0001) // used when serializing sections prior to sim start

         count = 0
         forsec axonal{ // loop through all axon sections

	     nseg = 1 + int(L/chunkSize/2.)*2  //nseg to get diameter

         for (x) {
             if (x > 0 && x < 1) {
                 count = count + 1
                 diams.resize(count)
                 diams.x[count-1] = diam(x)
                 lens.resize(count)
                 lens.x[count-1] = L/nseg
                 if( count == nseg_total ){
                     break
                 }
             }
         }
         if( count == nseg_total ){
             break
	 }
     }

         // get rid of the old axon
         forsec axonal{delete_section()}
         execute1("create axon[2]", CellRef)

         L_real = 0
         count = 0

         // new axon dependant on old diameters
         for i=0,1{
             access axon[i]
             L =  L_target/2
             nseg = nseg_total/2

             for (x) {
                 if (x > 0 && x < 1) {
                     diam(x) = diams.x[count]
                     L_real = L_real+lens.x[count]
                     count = count + 1
                 }
             }

             all.append()
             axonal.append()

             if (i == 0) {
                 v(0.0001) = i1
             } else {
                 v(0.0001) = i2
             }
         }

         nSecAxonal = 2
         soma[0] connect axon[0](0), 1
         axon[0] connect axon[1](0), 1

         print "Target stub axon length:", L_target, "um, equivalent length: ", L_real "um"
     }

 }

    



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()

  
}




endtemplate CA1_PC