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 fetch_soma_sec(section_name):
    cell_model = 'Hayton.hoc'
    h.load_file(cell_model)
    cell = h.L5PC
    soma = cell.soma[0]
    exec('sec = cell.' + section_name)
    return soma, sec

    def __init__(self, sim_vars, dt=0.1):

        h.load_file('stdrun.hoc')
        self.dt = dt
        self.sim_time = 300 # this will be rewritten in set_SEClamp
        h.celsius = 22

        #load known/default parameters
        params = yaml.load(open('params_example_start.yaml'))

        soma = h.Section(name='soma')
        soma.L = 15
        soma.diam = 15
        soma.cm = 1

        soma.insert('Narsg')
        #set known/default parameters
        for p in params['Channel']:
            cmd = 'soma(0.5).%s_Narsg = %s' % (p, params['Channel'][p])
            exec(cmd)
        #assign passed variables
        for sv in sim_vars:
            cmd = 'soma(0.5).%s_Narsg = %s' % (sv, sim_vars[sv])
            exec(cmd)
        self.recording_section = soma
        self.soma = soma

def importCell (fileName, cellName, cellArgs = None):
    h.initnrn()

    if cellArgs is None: cellArgs = [] # Define as empty list if not otherwise defined

    ''' Import cell from HOC template or python file into framework format (dict of sections, with geom, topol, mechs, syns)'''
    if fileName.endswith('.hoc'):
        h.load_file(fileName)
        if isinstance(cellArgs, dict):
            cell = getattr(h, cellName)(**cellArgs)  # create cell using template, passing dict with args
        else:
            cell = getattr(h, cellName)(*cellArgs) # create cell using template, passing list with args
    elif fileName.endswith('.py'):
        filePath,fileNameOnly = os.path.split(fileName)  # split path from filename
        if filePath not in sys.path:  # add to path if not there (need to import module)
            sys.path.insert(0, filePath)
        moduleName = fileNameOnly.split('.py')[0]  # remove .py to obtain module name
        exec('import ' + moduleName + ' as tempModule') in globals(), locals() # import module dynamically
        modulePointer = tempModule
        if isinstance(cellArgs, dict):
            cell = getattr(modulePointer, cellName)(**cellArgs) # create cell using template, passing dict with args
        else:
            cell = getattr(modulePointer, cellName)(*cellArgs)  # create cell using template, passing list with args
        sys.path.remove(filePath)
    else:
        print "File name should be either .hoc or .py file"
        return

    secDic, secListDic, synMechs = getCellParams(cell)
    return secDic, secListDic, synMechs

def coarse(hoc_filename,cube_length,save_filename):
	# INPUT: NEURON .hoc filename to import (str), voxel cube side length (fl/str for gcd), name of file to create for mesh output (str)
	# 		>> cube_length: 'gcd' (gcd of box dims) OR floating-point (must be common factor of all box dims)
	# This function reads in NEURON data and passes their info to
	# coarse_gen(), then associates the tets of the STEPS Tetmesh object returned
	# to the NEURON sections which exist inside of them. 
	# Returns a tet_hoc dictionary -- tet_hoc[tet_index] = [encapsulated hoc section references] -- as well as the Tetmesh object

	## GET HOC SECTION INFO ##
	h.load_file(hoc_filename)

	allp = [[],[],[]]
	for s in h.allsec():
		for j in range(int(h.n3d())):
			allp[0].append(h.x3d(j))
			allp[1].append(h.y3d(j))
			allp[2].append(h.z3d(j))

	maxl = [max(allp[0]),max(allp[1]),max(allp[2])]
	minl = [min(allp[0]),min(allp[1]),min(allp[2])]
	bdim = [ maxl[0] - minl[0], maxl[1] - minl[1], maxl[2] - minl[2] ]
	print "dims: ", bdim
	print "mins: ", minl

	## CREATE COARSE MESH ##
	if (cube_length == 'gcd'):
		gcd = fractions.gcd(fractions.gcd(bdim[0],bdim[1]),fractions.gcd(bdim[2],bdim[1]))
		print "GCD: ", gcd
		cube_length = gcd

	sm = coarse_gen(cube_length,bdim,minl,save_filename)

	## ASSOCIATE HOC SECTIONS WITH THEIR TETS ##
	tet_hoc = tet_associate(sm[0])
	return tet_hoc, sm[0]

 def enable_threads(self, n_threads, multisplit_on=True):
     """Enable threads in neuron Using the parall Neuron"""
     h.load_file('parcom.hoc')
     pc =h.ParallelComputeTool()
     pc.nthread(n_threads, 1)
     if multisplit_on:
         pc.multisplit(1)

    def set_up(self):

        h.load_file(1, 'NEURON_stuff/HayStuff/set_up_passive.hoc')

        # structuring the sections
        for sec in h.L5PC.basal:
            self.secs["basal"].append(sec)
            
        for sec in h.L5PC.somatic:
            self.secs["soma"].append(sec)

        for sec in h.L5PC.axonal:
            self.secs["axon"].append(sec)

        
        hoc_tuft = h.SectionList()
        hoc_tuft.subtree(sec=h.L5PC.apic[36])
        hoc_trunk = h.SectionList()
        for sec in h.L5PC.apical:
            hoc_trunk.append(sec=sec)
        
        for sec in hoc_tuft:
            if sec.name() != h.L5PC.apic[36].name():
                self.secs["tuft"].append(sec)

        for sec in self.secs["tuft"]:
            hoc_trunk.remove(sec=sec)

        for sec in hoc_trunk:
            self.secs["trunk"].append(sec)

        hoc_tuft = None     # making sure the object gets destroyed.
        self.bifurcation_info = (self.sections("trunk")[36].name(), 1)

def multisplit():
    h.load_file("parcom.hoc")
    parcom = h.ParallelComputeTool()
    parcom.multisplit(1)
    if settings.rank == 0:
        lb = parcom.lb
        print ('multisplit rank 0: %d pieces  Load imbalance %.1f%%' % (lb.npiece, (lb.thread_cxbal_ -1)*100))

    def _load_geom(self, filename):
        """Load the geometry of the model"""
        
        h5f = tables.openFile(filename)
        node = "/%s/%s" %(self.geometry_root, self.geometry_node_name) 
        geom = h5f.getNode(node)
        
        # List check for legacy code
        xml_data = None
        geom_data = geom.read()
        if isinstance(geom_data, list):
            xml_data = geom_data[0]  # get the string.
        else:
            xml_data = geom_data # The node is directly a string
        logger.debug(type (xml_data))
        logger.debug("xml_data is a list: %s" %isinstance(xml_data, list))
        tmp_file = 'temp.xml'
        f = open(tmp_file, 'w')
        f.write(xml_data)
        f.close()
        
#        import rdxml # This has to go ASAP they fix NEURON install
        h.load_file('celbild.hoc')
        cb = h.CellBuild(0)
        cb.manage.neuroml(tmp_file)
        cb.cexport(1)
        
        os.remove(tmp_file)

    def __init__(self,f_name, El, Rm, Ra, Cm, min_distance = 0., convert_to_3pt_soma = True):
        """
        El is the reversal potential of the leak conductance.
        Rm is the input resistance of the cell, in MOhm.
        Cm is the membrane capacitance.
        Ra is the axial resistance of the cell.
        All three variables are dictionaries, with the following keys:
           dend - value to be used for the dendrites
           soma - value to be used for the soma
           axon - value to be used for the axon
        """
        h.load_file('stdlib.hoc') # contains the notorious lambda rule

        # the path of the SWC file
        if convert_to_3pt_soma:
            self.swc_filename = '.'.join(f_name.split('.')[:-1]) + '_converted.swc'
            convert_morphology(f_name, self.swc_filename)
        else:
            self.swc_filename = f_name

        # parameters
        self.El = El
        self.Rm = Rm
        self.Ra = Ra
        self.Cm = Cm
        self.min_distance = min_distance

        self.load_morphology()
        self.compute_measures()
        self.insert_passive_mech()
        self.insert_active_mech()

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 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 load_hoc_model(self, model_dir, hoc_file):
     """Load an hoc files. It compiles the mod 
     before loading the hoc."""
     try:
         os.path.isfile(os.path.join (model_dir, hoc_file))
     except IOError:
         logger.error("Not existing file: %s" %e.value)
         
     old_dir = os.path.abspath(os.getcwd())
     logger.info("Path changed to %s" %(os.path.abspath(model_dir)))
     if model_dir != '' :
         os.chdir(model_dir)
     try:
         # Add all mod files into current directory
         self.find_mod_files()
     
         # If windows
         if os.name == 'nt':                
             self.windows_compile_mod_files('.')
             from neuron import h
             h.nrn_load_dll('./nrnmech.dll')
         else: # Anything else.
             call(['nrnivmodl'])
             import neuron            
             neuron.load_mechanisms('.') # Auto loading. Not needed anymore.
         from neuron import gui # to not freeze neuron gui
         from neuron import h
         logger.info("Loading model in %s from %s"%(model_dir, hoc_file))
         h.load_file(hoc_file)
     except Exception as e:
         logger.warning("Error running model: " + e.message)
     logger.info("Path changed back to %s" %old_dir)
     os.chdir(old_dir)
     return True

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 instantiate_swc(filename):
    """load an swc file and instantiate it"""
    h.load_file('stdgui.hoc')
    h.load_file('import3d.hoc')
    cell = h.Import3d_SWC_read()
    cell.input(filename)
    i3d = h.Import3d_GUI(cell, 0)
    i3d.instantiate(None)
    return i3d

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 generate_solinas2010():
    # Load a shor version of the Start.hoc file
    # that does not run the sim
    # it saves the network structure to the Sim_data dir

    if not os.path.exists('SimData'):
        os.mkdir('SimData')

    from neuron import h
    h.load_file('Start_test.hoc')

    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 __init__(self):
        
        self.groups = {}
#        self.vecRefs = [] 
#        self.synVecRefs = []
        self.refs = {}
        self.results_root = 'results'
        self.geometry_root = 'geometry'
        self.geometry_node_name = 'geom'        
        # Load the std run for NEURON
        h.load_file("stdrun.hoc")

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