def simulate(self, simtime=simtime, dt=dt, plotif=False, **kwargs):

        self.dt = dt
        self.simtime = simtime
        self.T = np.ceil(simtime / dt)
        self.trange = np.arange(0, self.simtime + dt, dt)

        self._init_network(**kwargs)
        if plotif:
            self._init_plots()

        # moose simulation
        moose.useClock(0, "/network/syns", "process")
        moose.useClock(1, "/network", "process")
        moose.useClock(2, "/plotSpikes", "process")
        moose.useClock(3, "/plotVms", "process")
        moose.useClock(3, "/plotWeights", "process")
        moose.setClock(0, dt)
        moose.setClock(1, dt)
        moose.setClock(2, dt)
        moose.setClock(3, dt)
        moose.setClock(9, dt)
        t1 = time.time()
        print "reinit MOOSE"
        moose.reinit()
        print "reinit time t = ", time.time() - t1
        t1 = time.time()
        print "starting"
        moose.start(self.simtime)
        print "runtime, t = ", time.time() - t1

        if plotif:
            self._plot()

    def simulate(self,simtime=simtime,dt=dt,plotif=False,**kwargs):
        
        self.dt = dt
        self.simtime = simtime
        self.T = np.ceil(simtime/dt)
        self.trange = np.arange(0,self.simtime+dt,dt)   
        
        self._init_network(**kwargs)
        if plotif:
            self._init_plots()
        
        # moose simulation
        moose.useClock( 0, '/network/syns', 'process' )
        moose.useClock( 1, '/network', 'process' )
        moose.useClock( 2, '/plotSpikes', 'process' )
        moose.useClock( 3, '/plotVms', 'process' )
        moose.useClock( 3, '/plotWeights', 'process' )
        moose.setClock( 0, dt )
        moose.setClock( 1, dt )
        moose.setClock( 2, dt )
        moose.setClock( 3, dt )
        moose.setClock( 9, dt )
        t1 = time.time()
        print 'reinit MOOSE -- takes a while ~20s.'
        moose.reinit()
        print 'reinit time t = ', time.time() - t1
        t1 = time.time()
        print 'starting'
        moose.start(self.simtime)
        print 'runtime, t = ', time.time() - t1

        if plotif:
            self._plot()

def setup_hdf5_output(model, neuron, filename=None, compartments=DEFAULT_HDF5_COMPARTMENTS):
    # Make sure /hdf5 exists
    if not moose.exists(HDF5WRITER_NAME):
        print('creating', HDF5WRITER_NAME)
        writer = moose.HDF5DataWriter(HDF5WRITER_NAME)
        #writer = moose.NSDFWriter(HDF5WRITER_NAME)
        writer.mode = 2 # Truncate existing file
        if filename is not None:
            writer.filename = filename
        moose.useClock(8, HDF5WRITER_NAME, 'process')
    else:
        print('using', HDF5WRITER_NAME)
        writer = moose.element(HDF5WRITER_NAME)

    for typenum,neur_type in enumerate(neuron.keys()):
        for ii,compname in enumerate(compartments):  #neur_comps):
            comp=moose.element(neur_type+'/'+compname)
            moose.connect(writer, 'requestOut', comp, 'getVm')

    if model.calYN:
        for typenum,neur_type in enumerate(neuron.keys()):
            for ii,compname in enumerate(compartments):  #neur_comps):
                comp=moose.element(neur_type+'/'+compname)
                for child in comp.children:
                    if child.className in {"CaConc", "ZombieCaConc"}:
                        cal = moose.element(comp.path+'/'+child.name)
                        moose.connect(writer, 'requestOut', cal, 'getCa')
                    elif  child.className == 'DifShell':
                        cal = moose.element(comp.path+'/'+child.name)
                        moose.connect(writer, 'requestOut', cal, 'getC')
    return writer

def main():
    makeModel()
    gsolve = moose.Gsolve("/model/compartment/gsolve")
    stoich = moose.Stoich("/model/compartment/stoich")
    stoich.compartment = moose.element("/model/compartment")
    stoich.ksolve = gsolve
    stoich.path = "/model/compartment/##"
    # solver.method = "rk5"
    # mesh = moose.element( "/model/compartment/mesh" )
    # moose.connect( mesh, "remesh", solver, "remesh" )
    moose.setClock(5, 1.0)  # clock for the solver
    moose.useClock(5, "/model/compartment/gsolve", "process")

    moose.reinit()
    moose.start(100.0)  # Run the model for 100 seconds.

    a = moose.element("/model/compartment/a")
    b = moose.element("/model/compartment/b")

    # move most molecules over to bgsolve
    b.conc = b.conc + a.conc * 0.9
    a.conc = a.conc * 0.1
    moose.start(100.0)  # Run the model for 100 seconds.

    # move most molecules back to a
    a.conc = a.conc + b.conc * 0.99
    b.conc = b.conc * 0.01
    moose.start(100.0)  # Run the model for 100 seconds.

    # Iterate through all plots, dump their contents to data.plot.
    displayPlots()

    quit()

def main():
    """
    Example of Interpol object.
    """
    simtime = 1.0
    simdt = 0.001
    model = moose.Neutral('/model')
    data = moose.Neutral('/data')
    interpol = moose.Interpol('/model/sin')
    vec = np.sin(np.linspace(-3.14, 3.14, 100))
    interpol.vector = vec
    interpol.xmax = 3.14
    interpol.xmin = -3.14
    recorded = moose.Table('/data/output')
    moose.connect(recorded, 'requestOut', interpol, 'getY')

    stimtab = moose.StimulusTable('/model/x')
    stimtab.stepSize = 0.0
    # stimtab.startTime = 0.0
    # stimtab.stopTime = simtime
    stimtab.vector = np.linspace(-4, 4, 1000)
    print((stimtab.vector))
    print((interpol.vector))
    moose.connect(stimtab, 'output', interpol, 'input')


    moose.setClock(0, simdt)
    moose.useClock(0, '/data/##,/model/##', 'process')
    moose.reinit()
    moose.start(simtime)
    plt.plot(np.linspace(0, simtime, len(recorded.vector)), recorded.vector, 'b-+', label='interpolated')
    plt.plot(np.linspace(0, simtime, len(vec)), vec, 'r-+', label='original')
    plt.show()

def test_crossing_single():
    """This function creates an ematrix of two PulseGen elements and
    another ematrix of two Table elements.

    The two pulsegen elements have same amplitude but opposite phase.

    Table[0] is connected to PulseGen[1] and Table[1] to Pulsegen[0].

    In the plot you should see two square pulses of opposite phase.

    """
    size = 2
    pg = moose.PulseGen('pulsegen', size)
    for ix, ii in enumerate(pg.vec):
        pulse = moose.element(ii)
        pulse.delay[0] = 1.0
        pulse.width[0] = 2.0
        pulse.level[0] = (-1)**ix
    tab = moose.Table('table', size)
    moose.connect(tab.vec[0], 'requestOut', pg.vec[1], 'getOutputValue', 'Single')
    moose.connect(tab.vec[1], 'requestOut', pg.vec[0], 'getOutputValue', 'Single')
    print 'Neighbors:'
    for t in tab.vec:
        print t.path
        for n in moose.element(t).neighbors['requestOut']:
            print 'requestOut <-', n.path
    moose.setClock(0, 0.1)
    moose.useClock(0, '/##', 'process')
    moose.start(5)
    for ii in tab.vec:
        t = moose.Table(ii).vector
        print len(t)
        pylab.plot(t)
    pylab.show()

def example():
    """In this example we create a square-pulse generator object and
    record the output using a table.

    The steps are:

    1. Create a PulseGen element `pulse`.

    2. Set `delay[0]=1.0`, `width[0]=0.2`, `level[0]=0.5`, so it
       generates 0.2 s wide square pulses with 0.5 amplitude every 1 s.

    3. Create a Table element `tab`.

    4. Connect the `outputValue` field of `pulse` to `tab`.

    5. We set tick-interval of ticks 0 and 1 to 0.01 and schedule
       `pulse` on tick 0 and `tab` on tick 1.

    5. Run the simulation for 5 s and save data to the ascii file
       `output_tabledemo.csv`.

    """
    pg = moose.PulseGen('pulse')
    pg.delay[0] = 1.0
    pg.width[0] = 0.2
    pg.level[0] = 0.5
    tab = moose.Table('tab')
    moose.connect(tab, 'requestOut', pg, 'getOutputValue')
    moose.setClock(0, 0.01)
    moose.setClock(1, 0.01)
    moose.useClock(0, pg.path, 'process')
    moose.useClock(1, tab.path, 'process')
    moose.reinit()
    moose.start(5.0)
    tab.plainPlot('output_tabledemo.csv')

 def setupSolver(self, path = '/hsolve'):
     """Setting up HSolver """
     hsolve = moose.HSolve( path )
     hsolve.dt = self.simDt
     moose.setClock(1, self.simDt)
     moose.useClock(1, hsolve.path, 'process')
     hsolve.target = self.cablePath

def assign_clocks(model_container_list, simdt, plotdt):
    """
    Assign clocks to elements under the listed paths.

    This should be called only after all model components have been
    created. Anything created after this will not be scheduled.

    """
    global inited
    # `inited` is for avoiding double scheduling of the same object
    if not inited:
        print(('SimDt=%g, PlotDt=%g' % (simdt, plotdt)))
        moose.setClock(0, simdt)
        moose.setClock(1, simdt)
        moose.setClock(2, simdt)
        moose.setClock(3, simdt)
        moose.setClock(4, plotdt)
        for path in model_container_list:
            print(('Scheduling elements under:', path))
            moose.useClock(0, '%s/##[TYPE=Compartment]' % (path), 'init')
            moose.useClock(1, '%s/##[TYPE=Compartment]' % (path), 'process')
            moose.useClock(2, '%s/##[TYPE=SynChan],%s/##[TYPE=HHChannel]' % (path, path), 'process')
            moose.useClock(3, '%s/##[TYPE=SpikeGen],%s/##[TYPE=PulseGen]' % (path, path), 'process')
        moose.useClock(4, '/data/##[TYPE=Table]', 'process')
        inited = True
    moose.reinit()

def testElecAlone():
    makeSpinyCompt()
    makeElecPlots()
    head2 = moose.element( '/n/head2' )
    kchan = moose.element( '/n/compt/K' )
    moose.setClock( 0, 2e-6 )
    moose.setClock( 1, 2e-6 )
    moose.setClock( 2, 2e-6 )
    moose.setClock( 8, 0.1e-3 )
    #moose.useClock( 0, '/n/#', 'init' )
    #moose.useClock( 1, '/n/#', 'process' )
    #moose.useClock( 2, '/n/#/#', 'process' )
    #print moose.wildcardFind( '/n/##[ISA=SpikeGen]' )
    moose.useClock( 0, '/n/##[ISA=Compartment]', 'init' )
    moose.useClock( 1, '/n/##[ISA=Compartment]', 'process' )
    moose.useClock( 2, '/n/##[ISA=ChanBase],/n/##[ISA=SynBase],/n/##[ISA=CaConc],/n/##[ISA=SpikeGen]','process')
    moose.useClock( 8, '/graphs/elec/#', 'process' )
    moose.reinit()
    moose.start( 0.1 )
    dumpPlots( 'instab.plot' )
    # make Hsolver and rerun
    hsolve = moose.HSolve( '/n/hsolve' )
    moose.useClock( 1, '/n/hsolve', 'process' )
    moose.setClock( 0, 2e-5 )
    moose.setClock( 1, 2e-5 )
    moose.setClock( 2, 2e-5 )
    hsolve.dt = 2e-5
    hsolve.target = '/n/compt'
    moose.reinit()
    #print kchan, ', Gbar = ', kchan.Gbar
    #kchan.Gbar = 0.1e-3
    #print 'Gbar = ', kchan.Gbar
    moose.start( 0.11 )
    dumpPlots( 'h_instab.plot' )

def main():
    makeModel()
    solver = moose.GslStoich("/model/compartment/solver")
    solver.path = "/model/compartment/##"
    solver.method = "rk5"
    mesh = moose.element("/model/compartment/mesh")
    moose.connect(mesh, "remesh", solver, "remesh")
    moose.setClock(5, 1.0)  # clock for the solver
    moose.useClock(5, "/model/compartment/solver", "process")

    moose.reinit()
    moose.start(100.0)  # Run the model for 100 seconds.

    a = moose.element("/model/compartment/a")
    b = moose.element("/model/compartment/b")

    # move most molecules over to b
    b.conc = b.conc + a.conc * 0.9
    a.conc = a.conc * 0.1
    moose.start(100.0)  # Run the model for 100 seconds.

    # move most molecules back to a
    a.conc = a.conc + b.conc * 0.99
    b.conc = b.conc * 0.01
    moose.start(100.0)  # Run the model for 100 seconds.

    # Iterate through all plots, dump their contents to data.plot.
    for x in moose.wildcardFind("/model/graphs/conc#"):
        moose.element(x[0]).xplot("scriptKineticSolver.plot", x[0].name)

    quit()

def main():
    make_spiny_compt()
    make_plots()
    moose.setClock( 0, 1e-5 )
    moose.setClock( 1, 1e-5 )
    moose.setClock( 2, 1e-5 )
    moose.setClock( 8, 0.1e-3 )
    moose.useClock( 0, '/n/#', 'init' )
    moose.useClock( 1, '/n/#', 'process' )
    moose.useClock( 2, '/n/#/#', 'process' )
    moose.useClock( 8, '/graphs/#', 'process' )
    moose.reinit()
    moose.start( 0.1 )
    dump_plots( 'instab.plot' )
    # make Hsolver and rerun
    hsolve = moose.HSolve( '/n/hsolve' )
    moose.useClock( 1, '/n/hsolve', 'process' )
    moose.setClock( 0, 2e-5 )
    moose.setClock( 1, 2e-5 )
    moose.setClock( 2, 2e-5 )
    hsolve.dt = 2e-5
    hsolve.target = '/n/compt'
    moose.reinit()
    moose.start( 0.1 )
    dump_plots( 'h_instab.plot' )

def test_elec_alone():
    eeDt = 2e-6
    hSolveDt = 2e-5
    runTime = 0.02

    make_spiny_compt()
    make_elec_plots()
    head2 = moose.element( '/n/head2' )
    moose.setClock( 0, 2e-6 )
    moose.setClock( 1, 2e-6 )
    moose.setClock( 2, 2e-6 )
    moose.setClock( 8, 0.1e-3 )
    moose.useClock( 0, '/n/##[ISA=Compartment]', 'init' )
    moose.useClock( 1, '/n/##[ISA=Compartment]', 'process' )
    moose.useClock( 2, '/n/##[ISA=ChanBase],/n/##[ISA=SynBase],/n/##[ISA=CaConc],/n/##[ISA=SpikeGen]','process')
    moose.useClock( 8, '/graphs/elec/#', 'process' )
    moose.reinit()
    moose.start( runTime )
    dump_plots( 'instab.plot' )
    print "||||", len(moose.wildcardFind('/##[ISA=HHChannel]'))
    # make Hsolver and rerun
    hsolve = moose.HSolve( '/n/hsolve' )
    moose.useClock( 1, '/n/hsolve', 'process' )
    hsolve.dt = 20e-6
    hsolve.target = '/n/compt'
    moose.le( '/n' )
    for dt in ( 20e-6, 50e-6, 100e-6 ):
        print 'running at dt =', dt
        moose.setClock( 0, dt )
        moose.setClock( 1, dt )
        moose.setClock( 2, dt )
        hsolve.dt = dt
        moose.reinit()
        moose.start( runTime )
        dump_plots( 'h_instab' + str( dt ) + '.plot' )

def testNeuroMeshMultiscale():
        useHsolve = 0
        runtime = 0.5
	elecDt = 10e-6
	chemDt = 0.005
	ePlotDt = 0.5e-3
	cPlotDt = 0.005
	plotName = 'nm.plot'

	makeNeuroMeshModel()
	print "after model is completely done"
	for i in moose.wildcardFind( '/model/chem/#/#/#/transloc#' ):
		print i[0].name, i[0].Kf, i[0].Kb, i[0].kf, i[0].kb

	makeChemPlots()
	makeElecPlots()
	makeCaPlots()
        for i in range( 10 ):
            moose.setClock( i, elecDt )
        for i in range( 10, 20 ):
            moose.setClock( i, chemDt )
	moose.setClock( 8, ePlotDt )
	moose.setClock( 18, cPlotDt )
        if useHsolve:
	    hsolve = moose.HSolve( '/model/elec/hsolve' )
	    #moose.useClock( 1, '/model/elec/hsolve', 'process' )
	    hsolve.dt = elecDt
	    hsolve.target = '/model/elec/compt'
	    moose.reinit()
        #soma = moose.element( '/model/elec/soma' )
        '''
        else:
	    moose.useClock( 0, '/model/elec/##[ISA=Compartment]', 'init' )
	    moose.useClock( 1, '/model/elec/##[ISA=Compartment]', 'process' )
	    moose.useClock( 2, '/model/elec/##[ISA=ChanBase],/model/##[ISA=SynBase],/model/##[ISA=CaConc]','process')

	moose.useClock( 1, '/model/elec/##[ISA=SpikeGen]', 'process' )
	moose.useClock( 2, '/model/##[ISA=SynBase],/model/##[ISA=CaConc]','process')
	#moose.useClock( 5, '/model/chem/##[ISA=PoolBase],/model/##[ISA=ReacBase],/model/##[ISA=EnzBase]', 'process' )
	#moose.useClock( 4, '/model/chem/##[ISA=Adaptor]', 'process' )
	moose.useClock( 4, '/model/chem/#/dsolve', 'process' )
	moose.useClock( 5, '/model/chem/#/ksolve', 'process' )
	moose.useClock( 6, '/model/chem/spine/adaptCa', 'process' )
	moose.useClock( 6, '/model/chem/dend/DEND/adaptCa', 'process' )
        '''

	moose.useClock( 18, '/graphs/chem/#', 'process' )
	moose.useClock( 8, '/graphs/elec/#,/graphs/ca/#', 'process' )
        moose.element( '/model/elec/soma' ).inject = 2e-10
        moose.element( '/model/chem/psd/Ca' ).concInit = 0.001
        moose.element( '/model/chem/spine/Ca' ).concInit = 0.002
        moose.element( '/model/chem/dend/DEND/Ca' ).concInit = 0.003
	moose.reinit()

	moose.start( runtime )
#        moose.element( '/model/elec/soma' ).inject = 0
#	moose.start( 0.25 )
        makeGraphics( cPlotDt, ePlotDt )

def main():
    """
    This demo shows how to start, stop, and continue a simulation. 
    This is commonly done when we want to run a model till settling, then
    change a parameter or deliver a stimulus, and then continue the
    simulation.

    Here, the model is just the output of a PulseGen object which 
    generates periodic pulses.
    The demo shows how to start the simulation. using the 
    *moose.reinit* command to reset the model to its initial state,
    and *moose.start* command to run the model for the specified duration. 
    We issue multiple *moose.start* commands and do different things to
    the model between them. First, we change the delay of the pulseGen.
    Then we show a number of ways to assign the timestep (dt) to the
    table object in the simulation.
    Note that throughout this simulation the pulsegen is going at a 
    uniform rate, it is just being sampled by the output table at
    different intervals.
    """

    dt = 0.1
    steps = 100
    simtime = dt * steps
    # Pulsegen is on tick 0, we can pre-emptively set its dt.
    moose.setClock(0, dt)
    table = setup_model()
    pulse = moose.element("/model/pulse")
    # The 'tick' field is on every object, we can use this to set its dt.
    moose.setClock(table.tick, dt)
    moose.reinit()
    clock = moose.element("/clock")
    print dt
    print "dt = ", dt, ", Total simulation time = ", simtime
    print "Running simulation for", simtime, "seconds"
    moose.start(simtime)
    print "Simulator time:", clock.currentTime
    # Here we change the pulse delay and then run again.
    pulse.delay[0] = 1.0
    moose.start(simtime)

    # We change the table tick and use a different dt for it:
    table.tick = 2
    moose.setClock(table.tick, dt * 2)
    moose.start(simtime)

    # Here is yet another way to change clocks used by the table
    moose.useClock(9, "/model/pulse/tab", "process")
    print table.tick
    moose.setClock(9, dt / 2.0)
    moose.start(simtime)

    # Finally, here we change the pulse delay to 1 second and run again.

    print "Simulator time at end of simulation", clock.currentTime
    pylab.plot(pylab.linspace(0, clock.currentTime, len(table.vector)), table.vector)
    pylab.show()

def makeModel():
		# create container for model
		r0 = 1e-6	# m
		r1 = 0.5e-6	# m. Note taper.
		num = 200
		diffLength = 1e-6 # m
		comptLength = num * diffLength	# m
		diffConst = 20e-12 # m^2/sec
		concA = 1 # millimolar
		dt4 = 0.02  # for the diffusion
		dt5 = 0.2   # for the reaction
                mfile = '../../Genesis_files/M1719.g'

		model = moose.Neutral( 'model' )
		compartment = moose.CylMesh( '/model/kinetics' )

		# load in model
                modelId = moose.loadModel( mfile, '/model', 'ee' )
		a = moose.element( '/model/kinetics/a' )
		b = moose.element( '/model/kinetics/b' )
		c = moose.element( '/model/kinetics/c' )

                ac = a.concInit
                bc = b.concInit
                cc = c.concInit

		compartment.r0 = r0
		compartment.r1 = r1
		compartment.x0 = 0
		compartment.x1 = comptLength
		compartment.diffLength = diffLength
		assert( compartment.numDiffCompts == num )

		# Assign parameters
                for x in moose.wildcardFind( '/model/kinetics/##[ISA=PoolBase]' ):
                    #print 'pools: ', x, x.name
                    x.diffConst = diffConst

		# Make solvers
		ksolve = moose.Ksolve( '/model/kinetics/ksolve' )
		dsolve = moose.Dsolve( '/model/dsolve' )
                # Set up clocks. The dsolver to know before assigning stoich
		moose.setClock( 4, dt4 )
		moose.setClock( 5, dt5 )
		moose.useClock( 4, '/model/dsolve', 'process' )
                # Ksolve must be scheduled after dsolve.
		moose.useClock( 5, '/model/kinetics/ksolve', 'process' )

		stoich = moose.Stoich( '/model/kinetics/stoich' )
		stoich.compartment = compartment
		stoich.ksolve = ksolve
		stoich.dsolve = dsolve
		stoich.path = "/model/kinetics/##"
                print 'dsolve.numPools, num = ', dsolve.numPools, num
                b.vec[num-1].concInit *= 1.01 # Break symmetry.

def main():
		# Schedule the whole lot
		moose.setClock( 4, 0.1 ) # for the computational objects
		moose.setClock( 5, 0.2 ) # clock for the solver
		moose.setClock( 8, 1.0 ) # for the plots
		# The wildcard uses # for single level, and ## for recursive.
		compartment = makeModel()
		ksolve = moose.Ksolve( '/model/compartment/ksolve' )
		stoich = moose.Stoich( '/model/compartment/stoich' )
		stoich.compartment = compartment
		stoich.ksolve = ksolve
		ksolve.stoich = stoich
		stoich.path = "/model/compartment/##"
		state = moose.SteadyState( '/model/compartment/state' )
		#solver.method = "rk5"
		#mesh = moose.element( "/model/compartment/mesh" )
		#moose.connect( mesh, "remesh", solver, "remesh" )
		#moose.useClock( 4, '/model/compartment/##', 'process' )
		moose.useClock( 5, '/model/compartment/ksolve', 'process' )
		moose.useClock( 8, '/model/graphs/#', 'process' )

		moose.reinit()
		state.stoich = stoich
		state.showMatrices()
		state.convergenceCriterion = 1e-7

		for i in range( 0, 100 ):
			getState( ksolve, state )
			# state.randomInit()
			# moose.start( 50.0 ) # Run the model for 100 seconds.
			# state.resettle()
			# print ksolve.nVec[0], sum( ksolve.nVec[0] )
			# print state.nIter, state.status, state.stateType, state.nNegEigenvalues, state.nPosEigenvalues, state.solutionStatus
		
		moose.start( 100.0 ) # Run the model for 100 seconds.

		a = moose.element( '/model/compartment/a' )
		b = moose.element( '/model/compartment/b' )

		# move most molecules over to b
		b.conc = b.conc + a.conc * 0.9
		a.conc = a.conc * 0.1
		moose.start( 100.0 ) # Run the model for 100 seconds.

		# move most molecules back to a
		a.conc = a.conc + b.conc * 0.99
		b.conc = b.conc * 0.01
		moose.start( 100.0 ) # Run the model for 100 seconds.

		# Iterate through all plots, dump their contents to data.plot.
		displayPlots()

		quit()

def main():
    utils.parser
    nml.loadNeuroML_L123('./two_cells_nml_1.8/two_cells.nml')
    #mumbl.loadMumbl("./two_cells_nml_1.8/mumbl.xml")
    table1 = utils.recordTarget('/tableA', '/cells/purkinjeGroup_0/Dend_37_41', 'vm')
    table2 = utils.recordTarget('/tableB', '/cells/granuleGroup_0/Soma_0', 'vm')
    moose.setClock(0, 5e-6)
    moose.useClock(0, '/##', 'process')
    moose.useClock(0, '/##', 'init')
    moose.reinit()
    utils.run(0.1)
    graphviz.writeGraphviz(__file__+".dot", ignore='/library')

def testNeuroMeshMultiscale():
        runtime = 0.5
	#elecDt = 0.2e-6
	elecDt = 10e-6
	chemDt = 0.0025
	ePlotDt = 0.5e-3
	cPlotDt = 0.0025
	plotName = 'nm.plot'

	makeNeuroMeshModel()
	print "after model is completely done"
	for i in moose.wildcardFind( '/model/chem/#/#/#/transloc#' ):
		print i[0].name, i[0].Kf, i[0].Kb, i[0].kf, i[0].kb

	makeChemPlots()
	makeElecPlots()
	makeCaPlots()
        for i in range (10):
	    moose.setClock( i, elecDt )
        for i in range ( 10, 20 ):
	    moose.setClock( i, chemDt )
	moose.setClock( 8, ePlotDt )
	moose.setClock( 18, cPlotDt )
	moose.useClock( 8, '/graphs/elec/#,/graphs/ca/#', 'process' )
	moose.useClock( 18, '/graphs/chem/#', 'process' )
	hsolve = moose.HSolve( '/model/elec/hsolve' )
	hsolve.dt = elecDt
	hsolve.target = '/model/elec/compt'
        plotlist = makeGraphics()
	moose.reinit()
        moose.element( '/model/elec/soma' ).inject = 2e-10
        moose.element( '/model/chem/psd/Ca' ).concInit = 0.001
        moose.element( '/model/chem/spine/Ca' ).concInit = 0.002
        moose.element( '/model/chem/dend/DEND/Ca' ).concInit = 0.003
	moose.reinit()

        numDivs = 200
        partialRuntime = runtime / numDivs
	max_voxel = find_max_voxel()
	voxel_val_dict = {'spine':numpy.zeros((max_voxel, numDivs)), 
			  'dend':numpy.zeros((max_voxel, numDivs)), 
			  'elec':numpy.zeros((max_voxel, numDivs)), 
			  'spineCaM':numpy.zeros((max_voxel, numDivs)), 
			  'psdCaM':numpy.zeros((max_voxel, numDivs))}
        for i in range( numDivs ):
	    moose.start( partialRuntime )
            voxel_val_dict = updateGraphics( plotlist, voxel_val_dict, i ) #Edited by Chaitanya
	    
#        moose.element( '/model/elec/soma' ).inject = 0
#	moose.start( 0.25 )
	save_NSDF(cPlotDt, ePlotDt, voxel_val_dict, [max_voxel, numDivs, partialRuntime]) #Edited by Chaitanya
	finalizeGraphics( plotlist, cPlotDt, ePlotDt ) 

def main():
		# Schedule the whole lot
		moose.setClock( 4, 0.1 ) # for the computational objects
		moose.setClock( 5, 0.1 ) # clock for the solver
		moose.setClock( 8, 1.0 ) # for the plots
		# The wildcard uses # for single level, and ## for recursive.
		#compartment = makeModel()
                moose.loadModel( '../Genesis_files/M1719.cspace', '/model', 'ee' )
                compartment = moose.element( 'model/kinetics' )
                compartment.name = 'compartment'
		ksolve = moose.Ksolve( '/model/compartment/ksolve' )
		stoich = moose.Stoich( '/model/compartment/stoich' )
		stoich.compartment = compartment
		stoich.ksolve = ksolve
		#ksolve.stoich = stoich
		stoich.path = "/model/compartment/##"
		state = moose.SteadyState( '/model/compartment/state' )
		moose.useClock( 5, '/model/compartment/ksolve', 'process' )
		moose.useClock( 8, '/model/graphs/#', 'process' )

		moose.reinit()
		state.stoich = stoich
		#state.showMatrices()
		state.convergenceCriterion = 1e-7

                moose.le( '/model/graphs' )
                a = moose.element( '/model/compartment/a' )
                b = moose.element( '/model/compartment/b' )
                c = moose.element( '/model/compartment/c' )

		for i in range( 0, 100 ):
			getState( ksolve, state )
		
		moose.start( 100.0 ) # Run the model for 100 seconds.

		b = moose.element( '/model/compartment/b' )
		c = moose.element( '/model/compartment/c' )

		# move most molecules over to b
		b.conc = b.conc + c.conc * 0.95
		c.conc = c.conc * 0.05
		moose.start( 100.0 ) # Run the model for 100 seconds.

		# move most molecules back to a
		c.conc = c.conc + b.conc * 0.95
		b.conc = b.conc * 0.05
		moose.start( 100.0 ) # Run the model for 100 seconds.

		# Iterate through all plots, dump their contents to data.plot.
		displayPlots()

		quit()