def run(tstop=1000, dt=0, V=-65):
    h.load_file('stdrun.hoc')
    #h.finitialize(V)
    if dt > 0:
        h.dt = dt
    h.tstop = tstop
    h.run()

def passive_soma(quad, show=False):
  """
  Creates the model with basic pyramidal passive properties.
  """
  # Load the hoc into neuron
  h('xopen(%s)' %quad.hocfile)
  h.load_file('stdrun.hoc')
  seclist = list(h.allsec())
  for sec in seclist:
    sec.insert('pas')
    sec.Ra = 200.
  
  # Current injection into soma or tip
  soma_sec = return_soma_seg(quad, h)
  stim_loc = 0.
  stim = h.IClamp(stim_loc, sec=soma_sec)
  stim.delay = 1 # ms
  stim.dur = 1 # ms
  stim.amp = 20 # nA
  
  # Run sim and record data
  (v, labels) = ez_record(h) # PyNeuron to record all compartments
  t, I = h.Vector(), h.Vector()
  t.record(h._ref_t)
  I.record(stim._ref_i)
  h.init()
  h.tstop = 10 # s?
  h.run()
  v = ez_convert(v) # Convert v to numpy 2D array
  
  # If show, plot, else just return v
  if show:

def simulate(tstop=25):
    """Initialize and run a simulation.

    :param tstop: Duration of the simulation.
    """
    h.tstop = tstop
    h.run()

def test_simple_synaptic_input():
    soma = h.Section()
    soma.insert("hh")
    soma.insert("pas")
    tl = TargetLocation(soma, 0.5)
    proc = h.ExpSyn 

    ac = ArtificialCell()
    fire_times = np.linspace(0,50,10)
    firetimes = h.Vector(fire_times)
    ac.setFireTimes(firetimes)
    
    ic = InputConnection(ac, tl, proc)

    v = h.Vector()
    v.record(soma(0.5)._ref_v)
   
    h.tstop = 100
    h.run()

    test_result = np.array(v)
    correct = np.load("test_simple_synaptic_input_data.npy")

    print sum(abs(test_result - correct))

    #pyl.plot(test_result)
    #pyl.plot(correct)
    #pyl.show()

    assert np.all(test_result == correct)

def simulation(tstop, with_time = False):
    """
    runs the simulation and returns the current and
    time vectors as Numpy arrays
    """
    h.load_file('stdrun.hoc')
    h.v_init = -70
    
    h.tstop = tstop
    VC_patch.dur1 = tstop

    # define vectors
    current = h.Vector()
    current.record(VC_patch._ref_i)
   
    if with_time is True:
        time = h.Vector()
        time.record(h._ref_t)

    h.run()

    if with_time is True:
        return (time, np.array(current)*1000.)
    else:
        return np.array(current)*1000. 

def main():
    soma = h.Section()
    soma.insert('pas')
    soma.L = 100
    soma.diam = 100
    weight_min = 0.005
    weight_max = 0.05
    mu = (np.log(weight_min)+np.log(weight_max))/2
    sigma = (np.log(weight_max)-mu)/3
    weights = np.sort(np.exp(np.random.normal(mu,sigma,size=200)))
    synapses = [AMPASynapse(soma, 0.5, 0, w) for w in weights]
    for i,syn in enumerate(synapses):
        syn.set_presynaptic_spike_times([10+i*50])
    rec = {}
    for lbl in 't','v','g':
        rec[lbl] = h.Vector()
    rec['t'].record(h._ref_t)
    rec['v'].record(soma(0.5)._ref_v)
    rec['g'].record(syn.syn._ref_g)
    h.load_file('stdrun.hoc')
    h.v_init = -70
    h.celsius = 37
    h.tstop = len(weights)*50 + 100
    h.run()
    import pylab as p
    p.subplot(2,1,1)
    p.plot(rec['t'],rec['v'],'k')
    p.ylabel('Voltage (mV)')
    p.subplot(2,1,2)
    p.plot(rec['t'],rec['g'],'r')
    p.xlabel('Time (ms)')
    p.ylabel('Conductance (uS)')
    p.show()

    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 find_vrest(h, section_name):
    h.load_file("stdrun.hoc")
    tstop = 100
    h.dt = dt = 0.1
    soma, sec = fetch_soma_sec(section_name)
    h.init()
    h.cvode.re_init()
    t_vec, soma_vm, sec_vm = record(soma, sec)
    h.execute('tstop = 100')
    h.run()
    vrest = np.array(sec_vm)[-1]
    return vrest

def get_vRest(myEPas):
    """ 
    this functions returns the voltage in the given section
    as a function ePas. 
    
    Example:
    >>> get_vRest(soma(0.5), myEPas=-68)
    """
    mycell.soma.e_pas = myEPas
    
    h.run()
    return mycell.soma.v

def run():
    """
    Runs sim and plots figures, saves to outfig
    """
    t_vec.record(h._ref_t)
    v_vec.record(dend(0.5)._ref_v)
    h.tstop = 20000       #stop sim at 20000 ms or 20 s

    h.run()  # run sim
    outfig = "v" + date + ".png"
    pp.plot(t_vec, v_vec)    # plot time versus voltage
    pp.savefig(outfig)       # save figure using outfig

def simulate_voltage(tstop, list_seg):
    """ runs a simulation and recall the voltage vectors
    defined in segments entered in **args
    return a list of HocVector objects
    """
    myvector_list = list()
    for seg in list_seg:
        myvector_list.append(HocVector(seg))

    h.tstop = tstop
    h.run()

    # return a list of HocVector objects
    return myvector_list

def custom_cable(length=526., tiprad=1.4, somarad=15.4, inj=1.2, 
                 tinj=1., injloc=1., Ra=35.4):
  """
  Simulate a simple passive current injection. Length=um, rads=um,
  inj=nA, tinj=ms, injloc=1 (tip). 
  Recorded from all 11 segments.
  """
  # Set up model
  h.load_file('stdrun.hoc')
  cell = h.Section()
  cell.nseg = 11   # It is a good idea to have nseg be an odd number
  cell.Ra = 35.4   # Ohm*cm
  cell.insert('pas')
  
  # Create structure
  h.pt3dadd(0,0,0,somarad,sec=cell)
  h.pt3dadd(length,0,0,tiprad,sec=cell)
  stim = h.IClamp(injloc, sec=cell)
  stim.delay = 5 # ms
  stim.dur = tinj # ms
  stim.amp = inj # nA
  print("Stim: %.2f nA, %.2f ms, at location %.2f" %(stim.amp, stim.dur,
                                                     injloc))
  
  # Segment positions, equall spaced from 0 to 1
  seg_positions = np.linspace(0,1,cell.nseg)
  # Use toolbox to record v
  # ez_record records v in all compartments by default
  (v,v_labels) = ez_record(h)
  # Manually record time and current
  t, I = h.Vector(), h.Vector()
  t.record(h._ref_t)
  I.record(stim._ref_i)
  
  # Run the simulation 
  h.init()
  h.tstop = 30
  h.run()
  # Use toolbox convert v into a numpy 2D array
  v = ez_convert(v) 
  
  # Plotting options
  fig = plt.figure()
  for i in range(cell.nseg):
    t = [v[u][i] for u in range(len(v))]
    ax = fig.add_subplot(cell.nseg, 1, i+1)
    ax.plot(t)
  plt.show()
  return h, v

def main(args):
    global nseg, nchan, simulator, _args
    _args = args
    loadModel(args.swc_file, args)
    print("Done loading")

    h.init()
    print("[INFO] Running NEURON for %s sec" % args.sim_time)
    t1 = time.time()
    h.tstop = 1e3 * float(args.sim_time)
    h.run()
    t = time.time() - t1
    print("Time taken by neuron: %s sec" % t)

    makePlots()

def passive_cable(stimloc=1., stimdur=5., show=False):
  """
  This simulates a pulse injected into a simple passive cable.
  """
  # Set up model
  h.load_file('stdrun.hoc')
  cell = h.Section()
  cell.nseg = 11   # It is a good idea to have nseg be an odd number
  cell.Ra = 35.4   # Ohm*cm
  cell.insert('pas')
  
  # create 3d structure
  h.pt3dadd(0,0,0,1.0,sec=cell)
  h.pt3dadd(1732,1732,1732,1.0,sec=cell)
  # Specify current injection
  stim = h.IClamp(stimloc,sec=cell) # Stim @ 1th end of segment
  stim.delay = 5   # ms
  stim.dur = stimdur   # ms
  stim.amp = 0.2   # nA
  print("Stim: %.2f nA, %.2f ms, at location %.2f" %(stim.amp, stim.dur,
                                                     stimloc))
  
  # Segment positions, equall spaced from 0 to 1
  seg_positions = np.linspace(0,1,cell.nseg)
  # Use toolbox to record v
  # ez_record records v in all compartments by default
  (v,v_labels) = ez_record(h)
  # Manually record time and current
  t, I = h.Vector(), h.Vector()
  t.record(h._ref_t)
  I.record(stim._ref_i)
  
  # Run the simulation 
  h.init()
  h.tstop = 30
  h.run()
  # Use toolbox convert v into a numpy 2D array
  v = ez_convert(v)
  
  # Plotting options
  if show:
    fig = plt.figure()
    for i in range(cell.nseg):
      t = [v[u][i] for u in range(len(v))]
      ax = fig.add_subplot(cell.nseg, 1, i+1)
      ax.plot(t)
    plt.show()
  return h, v

def iclamp(cell, sec, i_inj, v_init, tstop, dt, celsius=35, pos_i=0.5, pos_v=0.5):
    """
    Runs a NEURON simulation of the cell for the given parameters.

    :param sec: List with 1st entry the name of the section and 2nd entry the index (or None in case of soma)
    :type sec: list[str, int]
    :param i_inj: Amplitude of the injected current for all times t.
    :type i_inj: array_like
    :param v_init: Initial membrane potential of the cell.
    :type v_init: float
    :param tstop: Duration of a whole run.
    :type tstop: float
    :param dt: Time step.
    :type dt: float
    :param celsius: Temperature during the simulation (affects ion channel kinetics).
    :type celsius: float
    :param pos_i: Position of the IClamp on the Section (number between 0 and 1).
    :type pos_i: float
    :param pos_v: Position of the recording electrode on the Section (number between 0 and 1).
    :type pos_v: float
    :return: Membrane potential of the cell and time recorded at each time step.
    :rtype: tuple of three ndarrays
    """

    section = cell.substitute_section(sec[0], sec[1])

    # time
    t = np.arange(0, tstop + dt, dt)

    # insert an IClamp with the current trace from the experiment
    stim, i_vec, t_vec = section.play_current(i_inj, t, pos_i)

    # record the membrane potential
    v = section.record('v', pos_v)
    t = h.Vector()
    t.record(h._ref_t)

    # run simulation
    h.celsius = celsius
    h.v_init = v_init
    h.tstop = tstop
    h.steps_per_ms = 1 / dt  # change steps_per_ms before dt, otherwise dt not changed properly
    h.dt = dt
    h.run()

    return np.array(v), np.array(t)

def main():
    fixedInput = {'probability': 0.5, 'spikeTimes': [100,215]}

    neurons = [KhaliqRaman(i) for i in range(2)]
    for n in neurons:
        n.addFixedInput(fixedInput['probability'], fixedInput['spikeTimes'])
        n.addPoissonInputs(clusterSize=5, stimulusProps={'frequency': 0.1, 'noise': 1.})
        n.addSomaticVoltageRecorder()

    time = h.Vector()
    time.record(h._ref_t)

    h.load_file('stdrun.hoc')
    h.tstop = 400
    h.run()
    for n in neurons:
        #print('Neuron [%d] emitted %d spikes.' % (n.ID, len(n.spikeTimes)))
        saveNeuron('neuron.h5', n)

def runsim (izhm):
  global Iin, ax1, recvecs
  h.dt = 0.025
  f1 = []
  ax1=None
  if fig1 is not None: plt.clf()
  for Iin in IinRange:
    izhm.Iin=Iin
    h.init()

    if burstMode:
            izhm.Iin=IinHyper
            h.init()
            h.continuerun(120)
            izhm.Iin=Iin
            h.run()
    else:
            h.run()
    plotall()

def simulate(tstop):
    """
    simulate the somatic EPSC and
    returns the time (in ms) and current (in pA)
    """
    
    # simulation details
    h.load_file('stdrun.hoc')
    h.v_init = -70
    time, current = h.Vector(), h.Vector()

    time.record(h._ref_t)
    current.record(VC_patch._ref_i)
    
    h.tstop = tstop
    VC_patch.dur1 = tstop
    
    h.run()
    # note that current is in pA
    return (time, np.array(current)*1000.)

def run_sim(h, section_name, v_peak, tau_raise, tau_fall, onset=100):
    tstop = 500
    h.dt = dt = 0.1
    h.load_file("stdrun.hoc")
    soma, sec = fetch_soma_sec(section_name)
    v_rest = -75.711 # find_vrest(h, section_name)
    h.init()
    h.cvode.re_init()
    s_v, a_v = fetch_soma_apic_pots()
    vv = voltage_clamp(tstop, dt, v_rest, v_peak, tau_raise, tau_fall, onset)
    vc = h.SEClamp(sec(0.5))
    vc.rs = 0.001
    vc.dur1 = tstop
    vamp = h.Vector(vv)
    vamp.play(vc._ref_amp1, h.dt)
    t_vec, soma_vm, sec_vm = record(soma, sec)
    h.execute('tstop = ' + str(tstop))
    h.run()
    diff_v = np.array(a_v) - np.array(s_v)
    return t_vec, soma_vm, sec_vm, diff_v, vv

def main():
    pkj = h.DSB94()
    h.dt = 0.02
    Vrest = -74
    h.celsius = 37
    rec={}
    rec['v'] = h.Vector()
    rec['t'] = h.Vector()
    rec['v'].record(pkj.soma(0.5)._ref_v)  
    rec['t'].record(h._ref_t)
    iclamp = h.IClamp(pkj.soma(0.5))
    iclamp.amp = 0.300
    iclamp.dur = 1000
    iclamp.delay = 0
    h.finitialize(Vrest)
    h.tstop = 700
    h.run()
    plt.plot(rec['t'],rec['v'])
    plt.show()
    from ipdb import set_trace
    set_trace()