def test_syn():
    precell = BallStick()
    postcell = BallStick()
    nc = precell.connect2target(postcell.synlist[0])
    nc.weight[0] = 0.01
    nc.delay = 0
    stim = h.IClamp(0.5, sec=precell.soma)
    stim.amp = 0.700
    stim.delay = 700
    stim.dur = 1000
    vec = {}
    for var in 't', 'pre', 'post':
        vec[var] = h.Vector()
    vec['t'].record(h._ref_t)
    vec['pre'].record(precell.soma(0.5)._ref_v)
    vec['post'].record(postcell.soma(0.5)._ref_v)
    cvode = h.CVode()
    cvode.active(1)
    h.finitialize(-65)
    tstop = 2000
    while h.t < tstop:
        h.fadvance()
    with open("vm.out", "w") as out:
        for time, vsoma in zip(vec['t'], vec['post']):
            out.write("%g %g\n" % (time, vsoma))
    try:
        import matplotlib.pyplot as plt
        plt.plot(vec['t'], vec['pre'], vec['t'], vec['post'])
        plt.show()
    except ImportError:
        pass

  def initialize():
      global Epas
      h.celsius = celsius
      for sec in h.soma:
          h.distance()    
 
      for sec in h.allsec():
          sec.v = Epas
          sec.e_pas = Epas
  
          sec.insert("pas")
          sec.e_pas = Epas
          sec.g_pas = 1/Rm
          sec.Ra = rall
          sec.cm = cap
          sec.gnabar_hh2 = 0 
          sec.gkbar_hh2 = 0
          dist = h.distance(0.5)
          # sec.gcabar_it2 = gcat_func(dist)
          sec.gcabar_it2 = gcat
      for sec in h.soma:
          sec.gnabar_hh2 = gna 
          sec.gkbar_hh2 = gkdr
          # sec.gcabar_it2 = 0.1054*gcat
          sec.gcabar_it2 = gcat
                                                                           
      h.finitialize()       
      h.fcurrent()     
      cvode.re_init()

def test_cell():
    cell = BallStick()
    stim = h.IClamp(0.5, sec=cell.soma)
    stim.amp = 0.620
    stim.delay = 700
    stim.dur = 1000
    tm = h.Vector()
    vm = h.Vector()
    ca = h.Vector()
    tm.record(h._ref_t)
    vm.record(cell.soma(0.5)._ref_v)
    ca.record(cell.soma(0.5)._ref_cai)
    cvode = h.CVode()
    cvode.active(1)
    h.finitialize(-65)
    tstop = 2000
    while h.t < tstop:
        h.fadvance()
    with open("vm.out", "w") as out:
        for time, vsoma in zip(tm, vm):
            out.write("%g %g\n" % (time, vsoma))
    with open("ca.out", "w") as out:
        for time, conc in zip(tm, ca):
            out.write("%g %g\n" % (time, conc))
    try:
        import matplotlib.pyplot as plt
        plt.plot(tm, vm)
        plt.show()
    except ImportError:
        pass

def simulate(pool, tstop=1000, vinit=-55):
    ''' simulation control 
    Parameters
    ----------
    cell: NEURON cell
        cell for simulation
    tstop: int (ms)
        simulation time
    vinit: int (mV)
        initialized voltage
    '''
    h.finitialize(vinit)
    for i in pool:
        cell = pc.gid2cell(i)
        balance(cell)
    if h.cvode.active():
        h.cvode.active()
    else:
        h.fcurrent()
    h.frecord_init()
    h.tstop = tstop
    h.v_init = vinit
    pc.set_maxstep(0.5)
    h.stdinit()
    pc.psolve(tstop)

  def initialize(Tdist):
      global Epas
      h.celsius = celsius
      for sec in h.soma:
          h.distance()    
 
      for sec in h.allsec():
          sec.v = Epas
          sec.e_pas = Epas
          sec.insert("pas")
          sec.e_pas = Epas
          sec.g_pas = 1/Rm
          sec.Ra = rall
          sec.cm = cap
          sec.gnabar_hh2 = 0 
          sec.gkbar_hh2 = 0
          for seg in sec:
              if Tdist == 1:
                  seg.gcabar_it2 = gcat
              if Tdist == 2:
                  seg.gcabar_it2 = gcat * (1 + 0.04 * (h.distance(0) + sec.L * seg.x)) * 0.10539397661220173
              
      for sec in h.soma:
          sec.gnabar_hh2 = gna 
          sec.gkbar_hh2 = gkdr
          if Tdist == 1:
              seg.gcabar_it2 = gcat
          if Tdist == 2:
              seg.gcabar_it2 = gcat * 0.10539397661220173
                                                                           
      h.finitialize()       
      h.fcurrent()     
      cvode.re_init()

 def test_initializer_initialize(self):
     init = simulator.initializer
     orig_initialize = init._initialize
     init._initialize = Mock()
     h.finitialize(-65)
     self.assertTrue(init._initialize.called)
     init._initialize = orig_initialize

 def test_initializer_initialize(self):
     init = simulator.initializer
     orig_initialize = init._initialize
     init._initialize = Mock()
     h.finitialize(-65)
     init._initialize.assert_called()
     init._initialize = orig_initialize

 def initialise(self, vrest=-65):
     """
     Initialise the model, to launch before each simulations
     """
     for sec in h.allsec():
         h.finitialize(vrest, sec)
         h.fcurrent(sec)
     h.frecord_init()

 def reset(self):
     """Reset the state of the current network to time t = 0."""
     self.running = False
     self.t = 0
     self.tstop = 0
     self.t_start = 0
     self.segment_counter += 1
     h.finitialize()

    def go(self):
        self.set_recording()
        h.dt = self.dt
        h.tstop = self.sim_time
        h.finitialize(-60)#self.v_init)
        h.init()
        h.run()

        self.rec_i = self.rec_ina.to_python()

def ivcurve(mechanism_name, i_type, vmin=-100, vmax=100, deltav=1, transient_time=50, test_time=50, rs=1, vinit=-665):
    """
    Returns the (peak) current-voltage relationship for an ion channel.

    Args:
        mechanism_name = name of the mechanism (e.g. hh)
        i_type = which current to monitor (e.g. ik, ina)
        vmin = minimum voltage step to test
        vmax = maximum voltage step to test
        deltav = increment of voltage
        transient_time = how long to ignore for initial conditions to stabilize (ms)
        test_time = duration of the voltage clamp tests (ms)
        rs = resistance of voltage clamp in MOhm
        vinit = initialization voltage

    Returns:
        i = iterable of peak currents (in mA/cm^2)
        v = iterable of corresponding test voltages

    Note:
        The initialization potential (vinit) may affect the result. For example, consider
        the Hodgkin-Huxley sodium channel; a large fraction are inactivated at rest. Using a
        strongly hyperpolarizing vinit will uninactivate many channels, leading to more
        current.
    """
    from neuron import h
    import numpy
    h.load_file('stdrun.hoc')
    sec = h.Section()
    sec.insert(mechanism_name)
    sec.L = 1
    sec.diam = 1
    seclamp = h.SEClamp(sec(0.5))
    seclamp.amp1 = vinit
    seclamp.dur1 = transient_time
    seclamp.dur2 = test_time 
    seclamp.rs = rs
    i_record = h.Vector()
    i_record.record(sec(0.5).__getattribute__('_ref_' + i_type))
    result_i = []
    result_v = numpy.arange(vmin, vmax, deltav)
    for test_v in result_v:
        seclamp.amp2 = test_v
        h.finitialize(vinit)
        h.continuerun(transient_time)
        num_transient_points = len(i_record)
        h.continuerun(test_time + transient_time)
        i_record2 = i_record.as_numpy()[num_transient_points:]
        baseline_i = i_record2[0]
        i_record_shift = i_record2 - baseline_i
        max_i = max(i_record_shift)
        min_i = min(i_record_shift)
        peak_i = max_i if abs(max_i) > abs(min_i) else min_i
        peak_i += baseline_i
        result_i.append(peak_i)
    return result_i, result_v

def run_single_simulation(config, interactive):
    axon = Axon(config)
    axon.insert_stim(config['stim_position'], config['stim_amplitude'],
            config['stim_start_time'], config['stim_duration'])

    # set up recording vectors for python plots and the csv file
    t = h.Vector()
    t.record(h._ref_t)

    num_v_traces = config['num_v_traces']
    v_traces = []
    for i in range(num_v_traces):
        v = h.Vector()
        v.record(axon.section_at_f(
            # record at num_v_traces points along the axon, equally spaced
            # from eachother and from the end points (since we don't care
            # about things like the impedance mismatch at the ends)
            (i + 1) * 1.0 / (num_v_traces + 1))
            (Axon.middle)._ref_v)
        v_traces.append(v)

    # set up NEURON plotting code (if we're in an interactive session)
    if interactive:
        g = h.Graph()
        g.size(0, config['integration_time'], -80, 55)
        for i in range(num_v_traces):
            g.addvar('v(0.5)',
                    sec=axon.section_at_f((i+1) * 1.0 / (num_v_traces + 1)))

    # initialize the simulation
    h.dt = config['max_time_step']
    tstop = config['integration_time']
    h.finitialize(config['initial_membrane_potential'])
    h.fcurrent()

    # run the simulation
    if interactive:
        g.begin()
        while h.t < tstop:
            h.fadvance()
            g.plot(h.t)
        g.flush()
    else:
        while h.t < tstop:
            h.fadvance()

    # save the data as a csv
    with open(config['csv_filename'], 'w') as csv_file:
        # start with a header of the form "t_ms, V0_mV, V1_mv, V2_mV,..."
        csv_file.write(", ".join(
            ["t_ms"] + ["V{0}_mV".format(i) for i in range(num_v_traces)]
            ) + "\n")

        # write the time and each of the recorded voltages at that time
        for row in zip(t, *v_traces):
            csv_file.write(", ".join([str(x) for x in row]) + "\n")

    def run(self, v_init=-60, tstop=20000., dt=0.1,
                                        cvode=True, ga_use_half=False):
        '''
        Simulates this cell and all desired vectors are recorded. Uses fixed
        or variable timestep depending on the `cvode=` boolean.

        Parameters:
        ----------
        v_init : int, float
            The starting voltage of the simulation.
        tstop : int, float
            The maximum time of integration
        dt : float
            The desired integration step.
        cvode : bool
            Selects variable time step integration. Default is False.
        ga_use_half : bool
            Will only use the 2nd have of recordings for GA
        '''
        
        h.load_file('stdrun.hoc')
        h.v_init = v_init
        h.tstop = tstop
        h.dt = dt
        #set the recording into the vecs dictionary
        #the _ref dictionary contain the hoc object attribute references
        for key in self._ref.keys():
            #This makes sure we overwrite any vectors if we already ran a sim
            if isinstance(self.vecs['time'], np.ndarray):
                self.vecs[key] = h.Vector()
                self.vecs[key].record(self._ref[key])
            else:
                self.vecs[key].record(self._ref[key])

        
        
        if cvode:
            solver = h.CVode()
            solver.active(1)
            h.finitialize(h.v_init)
            solver.solve(h.tstop)
        else:
            h.CVode().active(0)
            h.finitialize()
            for t in range(0, int(h.tstop/h.dt)):
                h.fadvance()

        for key, val in self.vecs.iteritems():
            self.vecs[key] = np.array(val)

        if ga_use_half:
            for key, val in self.vecs.iteritems():
                self.vecs[key] = val[(val.size)/2:]

        return self.vecs

def integrate(a, b, c, t):
    # g.begin()
    k = 0
    h.finitialize()
    while h.t < tstop:
        h.fadvance()
        a[k - 1] = cell.ek
        b[k - 1] = cell.ko
        c[k - 1] = cell.ik
        t[k - 1] = h.dt * (k - 1)
        k = k + 1

def benchmark_cell():
    cell = BallStick()
    stim = h.IClamp(0.5, sec=cell.soma)
    stim.amp = 0.700
    stim.delay = 0
    stim.dur = 1e6
    cvode = h.CVode()
    cvode.active(1)
    h.finitialize(-65)
    tstop = 1e6
    while h.t < tstop:
        h.fadvance()

 def _pre_run(self):
     if not self.running:
         self.running = True
         local_minimum_delay = self.parallel_context.set_maxstep(self.default_maxstep)
         if state.vargid_offsets:
             logger.info("Setting up transfer on MPI process {}".format(state.mpi_rank))
             state.parallel_context.setup_transfer()
         h.finitialize()
         self.tstop = 0
         logger.debug("default_maxstep on host #%d = %g" % (self.mpi_rank, self.default_maxstep ))
         logger.debug("local_minimum_delay on host #%d = %g" % (self.mpi_rank, local_minimum_delay))
         if self.num_processes > 1:
             assert local_minimum_delay >= self.min_delay, \
                    "There are connections with delays (%g) shorter than the minimum delay (%g)" % (local_minimum_delay, self.min_delay)

    def go(self, sim_time=None):
        """
        Start the simulation once it's been intialized
        """

        self.set_recording()
        h.dt = self.dt
        
        h.finitialize(self.v_init)
        neuron.init()
        if sim_time:
            neuron.run(sim_time)
        else:
            neuron.run(self.sim_time)
        self.go_already = True

 def run(self, simtime):
     """Advance the simulation for a certain time."""
     if not self.running:
         self.running = True
         local_minimum_delay = self.parallel_context.set_maxstep(self.default_maxstep)
         h.finitialize()
         self.tstop = 0
         logger.debug("default_maxstep on host #%d = %g" % (self.mpi_rank, self.default_maxstep ))
         logger.debug("local_minimum_delay on host #%d = %g" % (self.mpi_rank, local_minimum_delay))
         if self.num_processes > 1:
             assert local_minimum_delay >= self.min_delay, \
                    "There are connections with delays (%g) shorter than the minimum delay (%g)" % (local_minimum_delay, self.min_delay)
     self.tstop += simtime
     logger.info("Running the simulation for %g ms" % simtime)
     self.parallel_context.psolve(self.tstop)

def trivial_ecs(scale):
    from neuron import h, crxd as rxd
    import numpy
    import warnings
    warnings.simplefilter("ignore", UserWarning)
    h.load_file('stdrun.hoc')
    tstop = 10
    if scale:   #variable step case
        h.CVode().active(True)
        h.CVode().event(tstop)
    else:           #fixed step case
        h.dt = 0.1

    sec = h.Section() #NEURON requires at least 1 section

    # enable extracellular RxD
    rxd.options.enable.extracellular = True

    # simulation parameters
    dx = 1.0    # voxel size
    L = 9.0     # length of initial cube
    Lecs = 21.0 # lengths of ECS

    # define the extracellular region
    extracellular = rxd.Extracellular(-Lecs/2., -Lecs/2., -Lecs/2.,
                                      Lecs/2., Lecs/2., Lecs/2., dx=dx,
                                      volume_fraction=0.2, tortuosity=1.6)

    # define the extracellular species
    k_rxd = rxd.Species(extracellular, name='k', d=2.62, charge=1,
                        atolscale=scale, initial=lambda nd: 1.0 if 
                        abs(nd.x3d) <= L/2. and abs(nd.y3d) <= L/2. and 
                        abs(nd.z3d) <= L/2. else 0.0)

    # record the concentration at (0,0,0)
    ecs_vec = h.Vector()
    ecs_vec.record(k_rxd[extracellular].node_by_location(0, 0, 0)._ref_value)
    h.finitialize()
    h.continuerun(tstop) #run the simulation
    
    # compare with previous solution 
    ecs_vec.sub(h.Vector(trivial_ecs_data[scale]))
    ecs_vec.abs()
    if ecs_vec.sum() > 1e-9:
        return -1
    return 0

    def min_sim(self, TSTOP=100):
        """
        Launch a minimal simulation to test the model and determine its resting potential empirically
        """
        vrec = h.Vector()
        vrec.record(self.soma(0.5)._ref_v)

        for sec in h.allsec():
            h.finitialize(-65, sec)
            h.fcurrent(sec)
        h.frecord_init()

        while h.t < TSTOP: #Launch a simulation
            h.fadvance()

        vrest = np.array(vrec)[-1]

        return vrest