def setup_sections(self):
        start = h.startsw()
        
        ###################################################
        # set up sections
        self.sections = []
        # old style, but it is need for section_name in hoc
        h(self.section_def_template % (self.name, len(self.tree)))
        for sec in h.allsec():
            self.sections.append(sec)


        ###################################################
        # connect sections
        for i,sec in enumerate(self.sections):
            parent = self.tree[i]
            #print "%d to %d" % (i, tree[i])
            if(parent != 0):
                sec.connect(self.sections[parent-1], 1, 0)

        self.num_compartment = 0
        for sec in h.allsec():
            self.num_compartment += 1

        self.setup_time += h.startsw() - start

def get_pos_data_short():
    """
    Get positions of all segments currently loaded in Neuron in a simple matrix.
    Section position information is not available.

    :returns: 
        Matrix (3 x nSegments) With x,y,z positions. 
    :rtype: :class:`~numpy.ndarray`

    Example:
        .. code-block:: python

            data = get_pos_data_short()
    """
    n = 0
    for sec in h.allsec():
        n += int(h.n3d())
    data = np.zeros([4, n])
    cnt = 0
    for sec in h.allsec():
        for i in xrange(int(h.n3d())):
            data[0, cnt] = h.x3d(i)
            data[1, cnt] = h.y3d(i)
            data[2, cnt] = h.z3d(i)
            data[3, cnt] = h.diam3d(i)
            cnt += 1
    return data

  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 get_pos_data():
    """
    Get positions x, y, z for all segments and their diameter. 

    :returns: 
        4 lists: x,y,z,d. One element per section where each element is
        a :class:`~numpy.ndarray`.

    Example:
        .. code-block:: python

            x,y,z,d = get_pos_data()
            for sec in xrange(len(x)):
                for seg in xrange(len(x[sec]):
                    print x[sec][seg], y[sec][seg], z[sec][seg]
    """
    x = []
    y = []
    z = []
    d = []
    for sec in h.allsec():
        n3d = int(h.n3d())
        x_i, y_i, z_i = np.zeros(n3d), np.zeros(n3d), np.zeros(n3d),
        d_i = np.zeros(n3d)
        for i in xrange(n3d):
            x_i[i] = h.x3d(i)
            y_i[i] = h.y3d(i)
            z_i[i] = h.z3d(i)
            d_i[i] = h.diam3d(i)
        x.append(x_i)
        y.append(y_i)
        z.append(z_i)
        d.append(d_i)
    return x, y, z, d

 def launch_visio(self):
     msg = "Plotting..."
     self.ui.statusbar.showMessage(msg, 3500)
     if self.visio == None:
         
         # Checking there are sections in the model.
         i = 0
         for sec in h.allsec():
             i += 1
         
         if i > 0:
             self.visio = Visio(self.ui.sec_info_label, self.manager)
             self.visio.draw_model()
             self.ui.selected_section.setEnabled(True)
         else:
             msg = """No model found, no section created. You need 
             to have at least one."""
             logger.warning(msg)
     else:
         #Raise the visio window
         self.visio.container.show()
     # Enabling the animation
     try:
         self.animation()
     except KeyError:
         # No simulation run an nothing loaded.
         # just pass
         pass

    def record_vectors(self, nrnManager):
        """Add a vecRef to record the vectors"""
        t_i_r = self.param['neuron_time_recording_interval']
        
        for spine_id in self.param['stimulated_spines']:
            spine = nrnManager.spines[spine_id]
            for syn in spine.synapses:
                pp = syn.chan     
                self.manager.create_time_record(time_interval_recording=t_i_r,
                                                point_process=pp)
        
        for var in self.param['var_to_plot']:
            for sec_rec in self.param['sec_to_rec']:
                if sec_rec == 'all':
                    self.manager.add_all_vecRef(var, 
                                                t_i_r)
                    break
                else:
                    for sec in h.allsec():
                        if sec.name() in self.param['sec_to_rec']:
                            self.manager.add_vecRef(var, 
                                                    sec, 
                                                    t_i_r)

        # Recording the synapses
        
        for spine_id in self.param['stimulated_spines']:
            spine = nrnManager.spines[spine_id]
            for syn in spine.synapses:
                self.manager.add_synVecRef(syn)    

 def compute_total_area(self):
     self.total_area = 0
     for sec in h.allsec():
         for seg in sec:
             self.total_area += h.area(seg.x, sec)
     if DEBUG:
         print('Total area: %.0f um^2.' % self.total_area)

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 get_var_data(self, var, time_point=0):
     """Retrieve the value of the `var` for the `time_point`.
     Prameters:
     var - variable to retrieve
     time_point - point in the simulation"""
     
     var_scalar = []
     for sec in h.allsec():
         var_value = 0
         #if self.manager.refs.has_key('VecRef'):
         #    for vecRef in self.manager.refs['VecRef']:
         #        if vecRef.sec.name() == sec.name():
         #            if vecRef.vecs.has_key(var):
         #                vec = vecRef.vecs[var]
         #                try:
         #                    var_value = vec[time_point]
         #                except IndexError:
         #                    pass # vector exist, but not initialized.
         sec_scalar = self.build_sec_scalar(sec, var_value)
         var_scalar.extend(sec_scalar)
             
                 
                     
     
     if len(var_scalar) == 0:
         logger.debug( "Var scalar 0 length. Var: %s point_time: %s" %(var, 
                                                               time_point))
     return np.array(var_scalar)

  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 set_kir_gkbar(self, gkbar):
     """Set the conductance of kir"""
     for sec in h.allsec():
         for seg in sec:
             for mech in seg:
                 if mech.name() == 'kir':
                     mech.gkbar = gkbar

    def draw_mayavi(self, x, y, z, d, edges):
        "Draw the surface the first time"

        # rendering disabled
        self.mayavi.visualization.scene.disable_render = True

        points = mlab.pipeline.scalar_scatter(x, y, z, d / 2.0)
        dataset = points.mlab_source.dataset
        dataset.point_data.get_array(0).name = "diameter"
        dataset.lines = np.vstack(edges)
        dataset.point_data.update()
        self.dataset = dataset

        # The tube
        src = mlab.pipeline.set_active_attribute(points, point_scalars="diameter")
        stripper = mlab.pipeline.stripper(src)
        tube = mlab.pipeline.tube(stripper, tube_sides=6, tube_radius=1)
        tube.filter.capping = True
        #        tube.filter.use_default_normal = False
        tube.filter.vary_radius = "vary_radius_by_absolute_scalar"
        self.tube = tube

        # Setting the voltage
        # Making room for the voltage
        v = []
        for sec in h.allsec():
            sec.push()
            v.extend(np.repeat(0.0, h.n3d()))
            h.pop_section()

        v = np.array(v)
        self.draw_surface(v, "v")

        # ReEnable the rendering
        self.mayavi.visualization.scene.disable_render = False

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 __init__(self):
     """Create a MorphologyDB with the current NEURON morphology"""
     self._children = {sec:[] for sec in h.allsec()}
     self._parents = {}
     self._connection_pts = {}
     for sec in h.allsec():
         parent_sec = parent(sec)
         if parent_sec is not None:
             self._children[parent_sec].append(sec)
             pt = (parent_sec, parent_loc(sec, parent))
             local_pt = (sec, h.section_orientation(sec=sec))
             if pt in self._connection_pts:
                 self._connection_pts[pt].append(local_pt)
             else:
                 self._connection_pts[pt] = [pt, local_pt]
         self._parents[sec] = parent_sec

def set_solve_type(domain=None, dimension=None, dx=None, nsubseg=None, method=None):
    """Specify the numerical discretization and solver options.
    
    domain -- a section or Python iterable of sections"""
    setting_default = False
    if domain is None:
        domain = h.allsec()
        setting_default = True
    elif isinstance(domain, nrn.Section):
        domain = [domain]
    
    # NOTE: These attributes are set on a per-nrn.Section basis; they cannot 
    #       assume Section1D objects exist because they might be specified before
    #       those objects are created
    
    # domain is now always an iterable (or invalid)
    if method is not None:
        raise RxDException('using set_solve_type to specify method is not yet implemented')
    if dimension is not None:
        if dimension not in (1, 3):
            raise RxDException('invalid option to set_solve_type: dimension must be 1 or 3')
        factory = lambda: dimension
        if setting_default:
            _dimensions.default_factory = factory
        for sec in domain:
            _dimensions[sec] = dimension 
    if dx is not None:
        raise RxDException('using set_solve_type to specify dx is not yet implemented')
    if nsubseg is not None:
        raise RxDException('using set_solve_type to specify nsubseg is not yet implemented')

 def _check_clean(self):
     """Check that all objects have been cleared from NEURON kernel.
     """
     # Release objects held by an internal buffer
     # See https://www.neuron.yale.edu/phpBB/viewtopic.php?f=2&t=3221
     h.Vector().size()    
     
     # Make sure nothing is hanging around in reference cycles 
     gc.collect()
     
     remaining = []
     
     # No sections left
     n = len(list(h.allsec()))
     if n > 0:
         remaining.append((n, 'Section'))
         
     # NetCon (and other object types?)
     for objtyp in ['NetCon']:
         n = len(h.List('NetCon'))
         if n > 0:
             remaining.append((n, 'NetCon'))
     
     # No point processes or artificial cells left
     for name, typ in Mechanism.all_mechanism_types().items():
         if typ['artificial_cell'] or typ['point_process']:
             n = len(h.List(name))
             if n > 0:
                 remaining.append((n, name))
         
     if len(remaining) > 0:
         msg = ("Cannot create new context--old objects have not been "
             "cleared: %s" % ', '.join(['%d %s' % rem for rem in remaining]))
         raise RuntimeError(msg)

 def build_subsets(self):
     '''
     NEURON staff
     adds sections in NEURON SectionList
     '''
     self.all = h.SectionList()
     for sec in h.allsec():
       self.all.append(sec=sec)  

 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 _create_sectionlists(self):
            '''Create section lists for different kinds of sections'''
            #list with all sections
            self.allsecnames = []
            self.allseclist = h.SectionList()
            for sec in h.allsec():
                if sec.name().find('Ib_') >= 0:
                    self.allsecnames.append(sec.name())
                    self.allseclist.append(sec=sec)

            #list of soma sections, assuming it is named on the format "soma*"
            self.nsomasec = 0
            self.somalist = h.SectionList()
            for sec in h.allsec():
                if sec.name().find('Ib_soma') >= 0:
                    self.somalist.append(sec=sec)
                    self.nsomasec += 1

    def insert_active_mech(self):
        ##### fast sodium and delayed rectifier potassium in all sections
        for sec in h.allsec():
            sec.insert('hh2')
            sec.ek = -90
            sec.ena = 50
            sec.vtraub_hh2 = -55
            sec.gnabar_hh2 = 0.05
            sec.gkbar_hh2 = 0.005

	##### SOMA
        for sec in self.soma:
            sec.insert('im')
            h.taumax_im = 1000
            sec.gkbar_im = 1e-3
            sec.insert('cad')
            sec.depth_cad = 1
            sec.taur_cad = 5
            sec.cainf_cad = 2.4e-4
            sec.kt_cad = 0
            sec.insert('it')
            sec.cai = 2.4e-4
            sec.cao = 2
            sec.eca = 120
            sec.gcabar_it = 0.002
            sec.insert('ical')
            sec.gcabar_ical = 1e-4
            sec.insert('KahpM95')
            sec.gbar_KahpM95 = 0.03
            sec.insert('kd')
            sec.gkdbar_kd = 1e-8
            sec.insert('napinst')
            sec.gbar_napinst = 1e-8

	##### DENDRITES
	for section in it.chain(self.basal,self.apical):
            section.insert('cad')
            section.depth_cad = 1
            section.taur_cad = 5
            section.cainf_cad = 2.4e-4
            section.kt_cad = 0
            section.cai = 2.4e-4
            section.cao = 2
            section.eca = 120
            section.insert('it')
            section.gcabar_it = 0.002
            section.insert('ical')
            section.gcabar_ical = 1e-4
            section.insert('KahpM95')
            section.gbar_KahpM95 = 0.02
            section.insert('kd')
            section.gkdbar_kd = 1e-8

        ##### AXON
        for section,distance in zip(self.axon,self.axon_length):
            if distance < 10: # AIS
                section.gnabar_hh2 = 0.25