objref recv_hold

proc refit_ihold() {
    // save original values to tstop and injclamp in case these are <1s
    injclamp_del = injclamp.del
    injclamp.del = 1000
    tstop_orig = tstop

    // Fit holding current
    tstop=injclamp.del // only run for initial period
    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
          printf("*** FATAL ERROR: holding current fitting algorithm found inappropriate model. ABORTING...\n")
          return
        } 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
              leeway = leeway * 2/3
              amp_increment = amp_increment / 2
              printf("*** GOING TO FINER resolution of amp_increment: %g and leeway: %g\n", amp_increment, leeway)
            }
          } 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)
          }
        }
      }
    }
    tstop = tstop_orig
    injclamp.del = injclamp_del
}