begintemplate PropagationSearch
public simulation_database,find_cutoffs,n_simulations,init,prop_test,parameter_choice,pc_gNa,pc_ISgNa,kinetics
external test_axon,axon2,load_axon,stim1,stim2,run,tstop,calc_spikes,n_extstim,n_stim,n_prop,st_extstim,st_stim,st_prop,n_cells
objref simulation_index,simulation_gNa,simulation_database

proc init(){
  print "calling init"
  simulation_index = new Vector(10,-1)
  simulation_gNa = new Vector(10,-1)
  simulation_database = new Matrix(10,10)
  col_nexts = 1 // column with n_extstim
  col_nst = 2 // column with n_stim
  col_np = 3 // column wit n_prop
  n_simulations = 0
  
  kinetics = 4
  pc_gNa = 0
  pc_ISgNa = 1
}

proc set_to_default_state(){
    if(kinetics==1){
        test_axon.resetwhole_Jonas(100*test_axon.axon[0].g_pas)
        if(n_cells>1){
            axon2.resetwhole_Jonas(100*test_axon.axon[0].g_pas)
        }
    }else{
        test_axon.resetwhole_modTraub(100*test_axon.axon[0].g_pas)
        if(n_cells>1){
            axon2.resetwhole_modTraub(100*test_axon.axon[0].g_pas)
        }
    }
}

func get_axon_state(){
    if(parameter_choice == pc_gNa){
        if(kinetics==1){
            return test_axon.axon[test_axon.extstim_site].gbar_nafJonas
        }else{
            return test_axon.axon[test_axon.extstim_site].gbar_nafTraub
        }
    }else{
        if(kinetics==1){
            return test_axon.axon[test_axon.IS].gbar_nafJonas
        }else{
            return test_axon.axon[test_axon.IS].gbar_nafTraub
        }
    }
}

proc set_axon_state(){local gNa
    gNa = $1
    print "in set_axon_state(): parameter_choice = ", parameter_choice
    if(parameter_choice == pc_gNa){
        if(kinetics==1){
            test_axon.resetwhole_Jonas(gNa)
            if(n_cells>1){
                axon2.resetwhole_Jonas(gNa)
            }
        }else{
            test_axon.resetwhole_modTraub(gNa)
            if(n_cells>1){
                axon2.resetwhole_modTraub(gNa)
            }
        }
    }else{
        if(kinetics==1){
            test_axon.resetwhole_Jonas(0.03,gNa)
            if(n_cells>1){
                axon2.resetwhole_Jonas(0.03,gNa)
            }
        }else{
            test_axon.resetwhole_modTraub(200*test_axon.axon[0].g_pas,gNa)
            if(n_cells>1){
                axon2.resetwhole_modTraub(200*test_axon.axon[0].g_pas,gNa)
            }
        }
    }
}

proc run_and_record(){local i,total_spikes,start_ind
  // run simulation and then record results in simulation_database
  run()
  calc_spikes()
  
  n_simulations = n_simulations + 1
  if(simulation_database.nrow<n_simulations){
    simulation_database.resize(2*simulation_database.nrow,simulation_database.ncol)
  }
  
  simulation_database.x[n_simulations-1][0] = get_axon_state()
  simulation_database.x[n_simulations-1][1] = n_extstim
  simulation_database.x[n_simulations-1][2] = n_stim
  simulation_database.x[n_simulations-1][3] = n_prop
  
  total_spikes = st_extstim.size()+st_stim.size()+st_prop.size()
  if(simulation_database.ncol<4+total_spikes){
      simulation_database.resize(simulation_database.nrow,4+total_spikes)
  }
  for(i=0;i<st_extstim.size();i=i+1){
      simulation_database.x[n_simulations-1][4+i] = st_extstim.x[i]
  }
  start_ind = 4+st_extstim.size()
  for(i=0;i<st_stim.size();i=i+1){
      simulation_database.x[n_simulations-1][start_ind+i] = st_stim.x[i]
  }
  start_ind = start_ind + st_stim.size()
  for(i=0;i<st_prop.size();i=i+1){
      simulation_database.x[n_simulations-1][start_ind+i] = st_prop.x[i]
  }

}

proc reset_database_index(){
  simulation_gNa = simulation_database.getcol(0)
  simulation_gNa.resize(n_simulations)
  simulation_index = simulation_gNa.sortindex()
}

proc set_run_data(){
    // set the simulation_database to the desired matrix
    init()
    simulation_database = $o1
    n_simulations = simulation_database.nrow
    reset_database_index()
}

func prop_test(){local beh
  // Takes in beh, the cutoff behavior you are looking for
  // and returns 1 if we have reached it
  beh = $1

  if(beh==0){ // looking for stim_site to fire 
    return n_stim>=1
  }else if(beh==1){ // looking for prop_site to fire 
    return n_prop>=1
  }else if(beh==2){ // looking for stim_site to fire at least twice
    return n_stim>=2
  }else if(beh==3){ // looking for prop_site to fire at least twice
    return n_prop>=2
  }else{
    print "unknown behavior to test"
    return 0
  }
}

obfunc refine_window(){local i,beh,dgNa,gNa localobj window, refwindow
    // inputs are: vector with two doubles for window, interval to jump
    // stim by, propagation cutoff we're looking for back propagation
    window = $o1
    beh = $2
    dgNa = $3
    refwindow = window.c()
    print "in refine window: dgNa=", dgNa, ", window=[", window.x[0][0], window.x[1][0], "]"
    i = 1
    while(window.x[0][0]+dgNa*i<window.x[1][0]){
	      gNa = window.x[0][0]+dgNa*i
        set_axon_state(gNa)
	      print "Looking for behavior: ", beh, get_axon_state()
	      run_and_record()
	      if(prop_test(beh)){
	          print "successfully found behavior ", beh
	          refwindow.x[1][0] = gNa
	          refwindow.x[1][1] = n_extstim
	          refwindow.x[1][2] = n_stim
	          refwindow.x[1][3] = n_prop
	          break
	      }else{
	          print "did not find behavior ", beh
	          refwindow.x[0][0] = gNa
	          refwindow.x[0][1] = n_extstim
	          refwindow.x[0][2] = n_stim
	          refwindow.x[0][3] = n_prop
	      }
	      i = i + 1
	      print "after iteration:"
	      print "i=",i," dgNa=",dgNa," gNa=",gNa,", window=[", window.x[0][0], window.x[1][0]
    }
    
    return refwindow
}

func intfir(){local i
    // check if axon is intrinsically firing. 
    
    stim1.del = 0 // make cell fire no matter what at beginning, see if it fires again afterward. Some cells fire initially just because of initial conditions, but aren't intrinsically firing.
    stim2.del = tstop
    
    run()
    calc_spikes()
    
    if(st_prop.indwhere(">",5)!=-1 && n_prop>1){ //axon is intrinsically firing
	      return 1
    }
    return 0
}

obfunc find_prop(){local depth localobj prop_info,window
    // Takes in ddel=$1, the granularity that you want to find T_c to
    // Outputs a vector with T_c, the value of t_f at stim_site at T_c
    depth = $1
    
    window = new Matrix(2,4)
    // columns stand for gNa, n_extstim, n_stim, n_prop
    // first row represents case with highest gNa and no propagation
    // second row represents case wit lowest gNa and propagation
    for(i=0;i<n_simulations;i=i+1){
        ind = simulation_index.x[i]
        if(simulation_database.x[ind][col_np]>=1){
            prev_ind = simulation_index.x[i-1]
            window.x[0][0] = simulation_database.x[prev_ind][0]
            window.x[0][1] = simulation_database.x[prev_ind][col_nexts]
            window.x[0][2] = simulation_database.x[prev_ind][col_nst]
            window.x[0][3] = simulation_database.x[prev_ind][col_np]
            window.x[1][0] = simulation_database.x[ind][0]
            window.x[1][1] = simulation_database.x[ind][col_nexts]
            window.x[1][2] = simulation_database.x[ind][col_nst]
            window.x[1][3] = simulation_database.x[ind][col_np]
            break
        }
    }
    
    if(window.x[0][0]==0){
        window.x[0][0] = 0
        window.x[0][1] = 0
        window.x[0][2] = 0
        window.x[1][0] = max_gNa
        window.x[1][1] = 0
        window.x[1][2] = 0
    }  
    
    window = refine_window(window,1,0.2)
    window = refine_window(window,1,0.1)
    if(simulation_database.getcol(col_np).indwhere(">=",1)==-1){
        print "need to try larger than gNa=",max_gNa," for propagation"
        prop_info = window.getrow(1)
        return prop_info
    }
    
    del_exp = 2
    while(del_exp<=depth){
        print "looking for propagation: ", del_exp, window.x[0][0], window.x[1][0]
        window = refine_window(window,1,2*10^(-del_exp))
        window = refine_window(window,1,1*10^(-del_exp))
        del_exp = del_exp + 1
    }
    
    prop_info = window.getrow(1)
    
    return prop_info
}

obfunc find_ref(){local depth,ind localobj ref_info,window
    // Takes in ddel=$1, the granularity that you want to find T_c to
    // Outputs a vector with T_c, the value of t_f at stim_site at T_c
    depth = $1
    
    print "in find_ref()"
    window = new Matrix(2,4)
    // columns stand for gNa, n_extstim, n_stim, n_prop
    // first row represents case with highest gNa and no propagation
    // second row represents case with lowest gNa and propagation
    for(i=0;i<n_simulations;i=i+1){
        ind = simulation_index.x[i]
        print "i=", i, " ind=", ind
        if(simulation_database.x[ind][col_np]>=2){
            prev_ind = simulation_index.x[i-1]
            print "prev_ind=", prev_ind
            window.x[0][0] = simulation_database.x[prev_ind][0]
            window.x[0][1] = simulation_database.x[prev_ind][col_nexts]
            window.x[0][2] = simulation_database.x[prev_ind][col_nst]
            window.x[0][3] = simulation_database.x[prev_ind][col_np]
            window.x[1][0] = simulation_database.x[ind][0]
            window.x[1][1] = simulation_database.x[ind][col_nexts]
            window.x[1][2] = simulation_database.x[ind][col_nst]
            window.x[1][3] = simulation_database.x[ind][col_np]
            break
        }
    }
    
    if(window.x[0][0]==0){
        ind = simulation_index.x[n_simulations-1]
        window.x[0][0] = simulation_database.x[ind][0]
        window.x[0][1] = simulation_database.x[ind][col_nexts]
        window.x[0][2] = simulation_database.x[ind][col_nst]
        window.x[0][3] = simulation_database.x[ind][col_np]
        window.x[1][0] = max_gNa
        window.x[1][1] = 0
        window.x[1][2] = 0
        window.x[1][3] = 0
    }  
    
    window = refine_window(window,3,0.2)
    window = refine_window(window,3,0.1)
    if(simulation_database.getcol(col_np).indwhere(">=",2)==-1){
        print "need to try larger than gNa=",max_gNa," for t_r<5 ms"
        ref_info = window.getrow(1)
        return ref_info
    }
    
    print "one round of find_ref, depth=",depth
    del_exp = 2
    while(del_exp<=depth){
        print "looking for refractoriness: ", del_exp, window.x[0][0], window.x[1][0]
        window = refine_window(window,3,2*10^(-del_exp))
        window = refine_window(window,3,1*10^(-del_exp))
        del_exp = del_exp + 1
    }
    
    ref_info = window.getrow(1)
    
    return ref_info
}

obfunc find_cutoffs(){local depth,int_firing localobj cutoffs,prop_info,ref_info
    // Find the cutoff points for propagation from extstim_site to prop_site.
    // Vary parameter_choice (default g_Na within the axon),and return 
    // for the given depth number of decimal places:
    // (1) the lowest parameter_choice with propagation past prop_site
    // (2) the lowest parameter_choice with propagation past prop_site within 5 ms
    // Options for parameter_choice are: pc_gNa, pc_ISgNa
    
    parameter_choice = $1
    print "in find_cutoffs(): parameter_choice = ", parameter_choice
    if(parameter_choice == pc_gNa){
        max_gNa = 1
    }else{
        max_gNa = 2
    }
    depth = $2
    print "ps.kinetics = ", kinetics
    
    cutoffs = new Matrix(2,4)
    // contains matrix with external stimulation and firing times for
    // prop_info and ref_info,
    // columns are gNa, 
    
    if(numarg()>2){
        set_run_data($o3)
    }else{
        // reset simulation database
        n_simulations = 0
        simulation_index = new Vector(10,-1)
        simulation_gNa = new Vector(10,-1)
        simulation_database = new Matrix(10,10)
        
        // check to see if axon is intrinsically firing
        print "checking if the axon is intrinsically firing"
        set_to_default_state()
        int_firing = intfir()
        print "intrinsically firing?: ", int_firing
        if(int_firing){
	          return cutoffs
        }
    }
    
    stim1.del = 5
    stim2.del = 10
    print "tstop=",tstop
    
    // find cutoffs for propagation
    print "finding propagation cutoff"
    prop_info = find_prop(depth)
    print "prop_info: ", prop_info.printf()
    cutoffs.setrow(0,prop_info)
    if(prop_info.x[0]==1){
        return cutoffs
    }
    
  print "resetting database"
  // set database index
  reset_database_index()

  print "finding refractoriness cutoff"
  // find refractory cutoff
  ref_info = find_ref(depth)
  print "ref_info: ", ref_info.printf()
  cutoffs.setrow(1,ref_info)

  return cutoffs
}

endtemplate PropagationSearch