//
// Entry point for the model. There are many options and switches here, some defined in 
// included hoc files. See startup_*.sh scripts in the same directory for examples of
// how to run this file. Most parameters are specified on the command-line.
//

VERBOSE=0
DEBUG=0
if (name_declared("verbose")==5) {
    if (verbose >= 1) {
      printf("*** verbose enabled!\n")
      VERBOSE = 1
    }
    if (verbose >= 2) {
      printf("*** debug enabled!\n")
      DEBUG = 1
    }
}

/*
 * Before anything, figure out if we're running Richy (Cell 1) or
 * Starfish (Cell 2).
 */
if (name_declared("cell") != 5) {
  printf("*** FATAL ERROR: no cell type! ABORTING...\n") 
  quit()
}
if (cell != 1 && cell != 2) {
  printf("*** FATAL ERROR: cell type must be 1 or 2! ABORTING...\n")
  quit()
}

load_file("nrngui.hoc")
if (cell == 1) {
  load_file("initactivesdRichy.hoc")
} else {
  load_file("initactivesdStarfish.hoc")
}
load_file("rn.hoc")
load_file("iprotocol.hoc")
load_file("params.hoc")
load_file("nsegfuncs.hoc")
load_file("param_active_cond.hoc")
load_file("active_mods.hoc")    // for changing gh on the fly across model
load_file("ihold_funcs.hoc")    // for refitting ihold after changing gh (above)
load_file("dendIh_funcs.hoc")

v_init = -74
cvode_active(1)
cvode.atol(0.005)

print "DC input resistance of cell is ", rn(), "megohms"

// For some reason this needs to be after calling rn()
geom_nseg() // need this as Cm changed from loading fitted ses file

// Define maximum conductance densities as per command-line args
err = 0
if (set_gbars()) { err = 1 }
if (err) {
  printf("*** FATAL ERROR. ABORTING... \n")
  quit()
}

// Fit holding current, only if ihold not specified as commandline argument
tstop=1000 // only run for initial period
objref recv_hold
if (name_declared("ihold") == 5) {
  holdclamp.amp = ihold
} else {
  injclamp_del = injclamp.del // probably 180, see in iprotocol.hoc
  injclamp.del=1000 // so that we can do the fitting properly, then set it back afterwards

  recv_hold = new Vector()
  recv_hold.record(&soma[0].v(0.5)) // first need to find spike, so use variable dt from cvode
  target_Vm = -60 - 13.8 // == -73.8mV; NB: -13.8mV for JP correction
  //leeway = 2 // leeway in fitting holding current; start off at 2mV
  leeway = 5 // leeway in fitting holding current; start off at 5mV
  //target_error = 0.05 // mV
  //target_error = 0.1 // mV - changed it as 0.05 wouldn't always converge
  target_error = 0.15 // mV - changed it as 0.1 wouldn't always converge


  // start from -15pA, we will then either go more positive (likely), or more negative, as needed
  start_amp = -0.015
  holdclamp.amp = start_amp // holdclamp defined in iprotocol.hoc
  stopspike_amp = start_amp // amount of amp to stop spike with
  amp_increment = 0.001 // amount to increment amp, at first
  found_spike = 0
  found_hold = 0
  while (!found_hold) {
    printf("*** FITTING holding clamp amp: trying %g nA\n", holdclamp.amp)
    run()

    /*
      Two phases: 1) Find amp to stop spiking at, using variable dt; 2) Find amp
      to get avg_Vm to target_Vm, using fixed dt. We need both phases because of
      Vector.record() functioning. If we just have one phase, and set variable dt,
      then we won't get accurate avg_Vm in the last fifth of initial period as there
      will be too few data points, so won't converge properly. If we have fixed dt,
      then likely won't catch the spike in the case where the spike threshold is
      lower than target_Vm. In this case, the fitting procedure will run indefinitely
      because the spiking will throw off avg_Vm.  
        So our solution are these two phases. The first is to definitively find
      the spike threshold (rather, amount of holding current needed to stop spiking).
      Then we start the second phase with fixed dt in recv_hold.record() vector. This
      is where we look for the holding current to keep avg_Vm close to target_Vm.
      *But* if in the course of this second phase we notice that our holding clamp
      amp is getting to be smaller (i.e., higher nA value) than the current needed to
      stop the spike, then we know that this is one of the model cells where spiking
      occurs even before we've reached target_Vm from below.  In this case we stop
      with the holding clamp amp set to stopspike_amp.
    */
    if (!found_spike) {
      if (recv_hold.max() >= -10) {
        stopspike_amp = holdclamp.amp - amp_increment
        printf("*** GENERATED SPIKE!!! FOUND holding clamp amp to stop spike at %g nA\n", stopspike_amp)
        found_spike = 1
      } else if (holdclamp.amp >= 0.015) {
        printf("*** CELL NOT INTRINSICALLY SPIKING; reached upper bound on check. Continuing with holding clamp fitting...\n")
        stopspike_amp = holdclamp.amp
        found_spike = 1
      } else {
        holdclamp.amp = holdclamp.amp + amp_increment
      }
   
      if (found_spike) {
        // now re-init vector and use fixed dt in order to get good avg_Vm calculation
        recv_hold = new Vector()
        recv_hold.record(&soma[0].v(0.5),20) // this will give us 50 timepoints in recv_hold
        holdclamp.amp = start_amp // reset holdclamp.amp too
      }
    } else {

      if (holdclamp.amp > stopspike_amp) {
        holdclamp.amp = holdclamp.amp - amp_increment // just to output what we actually tested for
        if (want_output) {
          write_err_and_quit()
        } else {
          printf("*** FATAL ERROR: holding current fitting algorithm found inappropriate model. ABORTING...\n")
          quit()
        }
      } else {

        // Get average of somatic Vm in last 200ms of initial period, so last 5th of 1s-long trace
        end = recv_hold.size() - 1
        start = (end) * 4/5 + 1
        avg_Vm = recv_hold.sum(start, end) / (end - start + 1)
        printf("*** AVG Vm %g mV\n", avg_Vm)

        if (avg_Vm >= target_Vm - leeway && avg_Vm <= target_Vm + leeway) {
          if (abs(abs(avg_Vm) - abs(target_Vm)) <= target_error) {
            // we found our holding current amp!
            printf("*** FOUND holding clamp amp at %g nA\n", holdclamp.amp)
            found_hold = 1
          } else {
            leeway = leeway * 2/3
            amp_increment = amp_increment / 2
            printf("*** GOING TO FINER resolution of amp_increment: %g and leeway: %g avg_Vm: %g target_Vm: %g\n", amp_increment, leeway, avg_Vm, target_Vm)
          }
        } else if (avg_Vm < target_Vm - leeway) {
          // still hyperpolarized wrt target Vm, so inject *less* hyperpolarizing current
          holdclamp.amp = holdclamp.amp + amp_increment
        } else if (avg_Vm > target_Vm - leeway) {
          // we are depolarized wrt target Vm, so inject *more* hyperpolarizing current
          holdclamp.amp = holdclamp.amp - amp_increment
        } else {
          // umm, something broke
          printf("*** FATAL ERROR: impossible control flow in holding current code; avg_Vm %g holdclamp.amp %g", avg_Vm, holdclamp.amp)
          quit()
        }
      }
    }
  }

  injclamp.del = injclamp_del
}
  
tstop=1300
run()