def main():
	"""
	This snippet shows the use of several objects.
	This snippet sets up a StimulusTable to control a RandSpike which
	sends its outputs to two places: to a SimpleSynHandler on an IntFire,
	which is used to monitor spike arrival, and to various Stats objects.
    Each of these are recorded and plotted.
	The StimulusTable has a sine-wave waveform.
	"""
        make_model()

        moose.reinit()
        moose.start( runtime )
        plots = moose.element( '/plots' )
        plot1 = moose.element( '/plot1' )
        plot2 = moose.element( '/plot2' )
        plotf = moose.element( '/plotf' )
        t = [i * dt for i in range( plot1.vector.size )]
        pylab.plot( t, plots.vector, label='stimulus' )
        pylab.plot( t, plot1.vector, label='spike rate mean' )
        pylab.plot( t, plot2.vector, label='Vm mean' )
        pylab.plot( t, plotf.vector, label='Vm' )
        pylab.legend()
        pylab.show()

	'''

def main():
    global synSpineList 
    global synDendList 
    numpy.random.seed( 1234 )
    rdes = buildRdesigneur( )
    for i in elecFileNames:
        print(i)
        rdes.cellProtoList = [ ['./cells/' + i, 'elec'] ]
        rdes.buildModel( )
        assert( moose.exists( '/model' ) )
        synSpineList = moose.wildcardFind( "/model/elec/#head#/glu,/model/elec/#head#/NMDA" )
        temp = set( moose.wildcardFind( "/model/elec/#/glu,/model/elec/#/NMDA" ) )

        synDendList = list( temp - set( synSpineList ) )
        moose.reinit()
        buildPlots( rdes )
        # Run for baseline, tetanus, and post-tetanic settling time 
        t1 = time.time()
        probeStimulus( baselineTime )
        tetanicStimulus( tetTime )
        probeStimulus( postTetTime )
        print(('real time = ', time.time() - t1))

        printPsd( i + ".fig5" )
        saveAndClearPlots( i + ".fig5" )
        moose.delete( '/model' )
        rdes.elecid = moose.element( '/' )

def run(nogui):
    
    reader = NML2Reader(verbose=True)

    filename = 'test_files/NML2_SingleCompHHCell.nml'
    print('Loading: %s'%filename)
    reader.read(filename, symmetric=True)
    
    
    msoma = reader.getComp(reader.doc.networks[0].populations[0].id,0,0)
    print(msoma)
    
    
    data = moose.Neutral('/data')
    
    pg = reader.getInput('pulseGen1')
    
    inj = moose.Table('%s/pulse' % (data.path))
    moose.connect(inj, 'requestOut', pg, 'getOutputValue')
    
    
    vm = moose.Table('%s/Vm' % (data.path))
    moose.connect(vm, 'requestOut', msoma, 'getVm')
    
    simdt = 1e-6
    plotdt = 1e-4
    simtime = 300e-3
    #moose.showmsg( '/clock' )
    for i in range(8):
        moose.setClock( i, simdt )
    moose.setClock( 8, plotdt )
    moose.reinit()
    moose.start(simtime)
    
    print("Finished simulation!")
    
    t = np.linspace(0, simtime, len(vm.vector))
    
    if not nogui:
        import matplotlib.pyplot as plt

        vfile = open('moose_v_hh.dat','w')

        for i in range(len(t)):
            vfile.write('%s\t%s\n'%(t[i],vm.vector[i]))
        vfile.close()
        plt.subplot(211)
        plt.plot(t, vm.vector * 1e3, label='Vm (mV)')
        plt.legend()
        plt.title('Vm')
        plt.subplot(212)
        plt.title('Input')
        plt.plot(t, inj.vector * 1e9, label='injected (nA)')
        #plt.plot(t, gK.vector * 1e6, label='K')
        #plt.plot(t, gNa.vector * 1e6, label='Na')
        plt.legend()
        plt.figure()
        test_channel_gates()
        plt.show()
        plt.close()

    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 auto-schedules
        #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' )
        for i in range(10):
            moose.setClock( i, 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 singleCompt( name, params ):
    mod = moose.copy( '/library/' + name + '/' + name, '/model' )
    A = moose.element( mod.path + '/A' )
    Z = moose.element( mod.path + '/Z' )
    Z.nInit = 1
    Ca = moose.element( mod.path + '/Ca' )
    CaStim = moose.element( Ca.path + '/CaStim' )
    runtime = params['preStimTime'] + params['stimWidth'] + params['postStimTime'] 
    steptime = 100

    CaStim.expr += ' + x2 * (t > ' + str( runtime ) + ' ) * ( t < ' + str( runtime + steptime ) +  ' )'
    print(CaStim.expr)
    tab = moose.Table2( '/model/' + name + '/Atab' )
    #for i in range( 10, 19 ):
        #moose.setClock( i, 0.01 )
    ampl = moose.element( mod.path + '/ampl' )
    phase = moose.element( mod.path + '/phase' )
    moose.connect( tab, 'requestOut', A, 'getN' )
    ampl.nInit = params['stimAmplitude'] * 1
    phase.nInit = params['preStimTime']

    ksolve = moose.Ksolve( mod.path + '/ksolve' )
    stoich = moose.Stoich( mod.path + '/stoich' )
    stoich.compartment = mod
    stoich.ksolve = ksolve
    stoich.path = mod.path + '/##'
    runtime += 2 * steptime

    moose.reinit()
    moose.start( runtime )
    t = np.arange( 0, runtime + 1e-9, tab.dt )
    return name, t, tab.vector

def main():
        """
        This example illustrates loading, and running a kinetic model 
        for a bistable positive feedback system, defined in kkit format. 
        This is based on Bhalla, Ram and Iyengar, Science 2002.

        The core of this model is a positive feedback loop comprising of
        the MAPK cascade, PLA2, and PKC. It receives PDGF and Ca2+ as 
        inputs.

        This model is quite a large one and due to some stiffness in its
        equations, it runs somewhat slowly. 

        The simulation illustrated here shows how the model starts out in
        a state of low activity. It is induced to 'turn on' when a 
        a PDGF stimulus is given for 400 seconds. 
        After it has settled to the new 'on' state, model is made to 
        'turn off'
        by setting the system calcium levels to zero for a while. This
        is a somewhat unphysiological manipulation!
        """
        solver = "gsl"  # Pick any of gsl, gssa, ee..
        #solver = "gssa"  # Pick any of gsl, gssa, ee..
	mfile = '../../genesis/acc35.g'
	runtime = 2000.0
	if ( len( sys.argv ) == 2 ):
                solver = sys.argv[1]
	modelId = moose.loadModel( mfile, 'model', solver )
        # Increase volume so that the stochastic solver gssa 
        # gives an interesting output
        compt = moose.element( '/model/kinetics' )
        compt.volume = 5e-19 

	moose.reinit()
	moose.start( 500 ) 
        moose.element( '/model/kinetics/PDGFR/PDGF' ).concInit = 0.0001
	moose.start( 400 ) 
        moose.element( '/model/kinetics/PDGFR/PDGF' ).concInit = 0.0
	moose.start( 2000 ) 
        moose.element( '/model/kinetics/Ca' ).concInit = 0.0
	moose.start( 500 ) 
        moose.element( '/model/kinetics/Ca' ).concInit = 0.00008
	moose.start( 2000 ) 

	# Display all plots.
        img = mpimg.imread( 'mapkFB.png' )
        fig = plt.figure( figsize=(12, 10 ) )
        png = fig.add_subplot( 211 )
        imgplot = plt.imshow( img )
        ax = fig.add_subplot( 212 )
	x = moose.wildcardFind( '/model/#graphs/conc#/#' )
        t = numpy.arange( 0, x[0].vector.size, 1 ) * x[0].dt
        ax.plot( t, x[0].vector, 'b-', label=x[0].name )
        ax.plot( t, x[1].vector, 'c-', label=x[1].name )
        ax.plot( t, x[2].vector, 'r-', label=x[2].name )
        ax.plot( t, x[3].vector, 'm-', label=x[3].name )
        plt.ylabel( 'Conc (mM)' )
        plt.xlabel( 'Time (seconds)' )
        pylab.legend()
        pylab.show()

def deliverStim(currTime):
	global injectionCurrent	
	global spineVm
	global somaVm
	if numpy.fabs( currTime - baselineTime ) < frameRunTime/2.0 :
		#start
		eList = moose.wildcardFind( '/model/elec/#soma#' )
		assert( len(eList) > 0 )
		eList[0].inject = injectionCurrent
		#print "1. injected current = ", injectionCurrent
		injectionCurrent += deltaCurrent
		#print "del stim first ", moose.element('/clock').currentTime
	if numpy.fabs( currTime - baselineTime - currPulseTime) < frameRunTime/2.0 :
		#end
		eList = moose.wildcardFind( '/model/elec/#soma#' )
		assert( len(eList) > 0 )
		eList[0].inject = 0.0
		#print "2. injected current = ", injectionCurrent
		#print "del stim second ", moose.element('/clock').currentTime
	if runtime - currTime < frameRunTime * 2.0 :
		#print "3. reinit-ing"
		somaVm.append( moose.element( '/graphs/VmTab' ).vector )
		spineVm.append( moose.element( '/graphs/eSpineVmTab' ).vector )
		iList.append(injectionCurrent)
		if injectionCurrent < maxCurrent :
			moose.reinit()	

 def resetSimulation(self, runTime, updateInterval, simulationInterval):
     self.runTime            = runTime
     self.updateInterval     = updateInterval
     self.simulationInterval = simulationInterval
     self.pause              = False
     moose.reinit()
     self.simulationReset.emit()

def test():
    global finish_all_
    os.environ['MOOSE_STREAMER_ADDRESS'] = 'http://127.0.0.1:%d'%port_
    done = mp.Value('d', 0)
    q = mp.Queue()
    client = mp.Process(target=socket_client, args=(q, done))
    client.start()

    time.sleep(0.1)

    print( '[INFO] Socket client is running now' )
    ts = models.simple_model_a()
    moose.reinit()
    time.sleep(0.1)
    # If TCP socket is created, some delay is often neccessary before start. Don't
    # know why. probably some latency in a fresh TCP socket. A TCP guru can
    # tell.
    moose.start(50)
    print( 'MOOSE is done' )

    time.sleep(0.5)
    done.value = 1
    res = q.get()
    client.join()

    if not res:
        raise RuntimeWarning('Nothing was streamed')
    for k in res:
        a = res[k][1::2]
        b = moose.element(k).vector
        print(k, len(a), len(b))
        assert( (a==b).all())
    print( 'Test 1 passed' )

 def resetAndStartSimulation(self):
     """TODO this should provide a clean scheduling through all kinds
     of simulation or default scheduling should be implemented in MOOSE
     itself. We need to define a policy for handling scheduling. It can
     be pushed to the plugin-developers who should have knowledge of
     the scheduling criteria for their domain."""
     settings = config.MooseSetting()
     try:
         simdt_kinetics = float(settings[config.KEY_KINETICS_SIMDT])
     except ValueError:
         simdt_kinetics = 0.1
     try:
         simdt_electrical = float(settings[config.KEY_ELECTRICAL_SIMDT])
     except ValueError:
         simdt_electrical = 0.25e-4
     try:
         plotdt_kinetics = float(settings[config.KEY_KINETICS_PLOTDT])
     except ValueError:
         plotdt_kinetics = 0.1
     try:
         plotdt_electrical = float(settings[config.KEY_ELECTRICAL_PLOTDT])
     except ValueError:
         plotdt_electrical = 0.25e-3
     try:
         simtime = float(settings[config.KEY_SIMTIME])
     except ValueError:
         simtime = 1.0
     moose.reinit()
     view = self.plugin.getRunView()
     moose.start(simtime)
     if view.getCentralWidget().plotAll:
         view.getCentralWidget().plotAllData()
     self.setCurrentView('run')        

def run_single_channel(channelname, Gbar, simtime, simdt=testutils.SIMDT, plotdt=testutils.PLOTDT):
    testId = uuid.uuid4().int
    container = moose.Neutral('test%d' % (testId))
    model_container = moose.Neutral('%s/model' % (container.path))
    data_container = moose.Neutral('%s/data' % (container.path))
    params = testutils.setup_single_compartment(
        model_container, data_container,
        channelbase.prototypes[channelname],
        Gbar)
    vm_data = params['Vm']
    gk_data = params['Gk']
    ik_data = params['Ik']
    testutils.setup_clocks(simdt, plotdt)
    testutils.assign_clocks(model_container, data_container)
    moose.reinit()
    print 'Starting simulation', testId, 'for', simtime, 's'
    moose.start(simtime)
    print 'Finished simulation'
    vm_file = 'data/%s_Vm.dat' % (channelname)
    gk_file = 'data/%s_Gk.dat' % (channelname)
    ik_file = 'data/%s_Ik.dat' % (channelname)
    tseries = np.array(range(len(vm_data.vec))) * simdt
    print 'Vm:', len(vm_data.vec), 'Gk', len(gk_data.vec), 'Ik', len(ik_data.vec)
    data = np.c_[tseries, vm_data.vec]
    np.savetxt(vm_file, data)
    print 'Saved Vm in', vm_file
    print len(gk_data.vec), len(vm_data.vec)
    data = np.c_[tseries, gk_data.vec]
    np.savetxt(gk_file, data)
    print 'Saved Gk in', gk_file
    data = np.c_[tseries, ik_data.vec]
    np.savetxt(ik_file, data)
    print 'Saved Gk in', ik_file
    return params

def main():
    """
    A toy compartmental neuronal + chemical model that causes bad things
    to happen to the hsolver, as of 28 May 2013. Hopefully this will
    become irrelevant soon.
    """
    fineDt = 1e-5
    coarseDt = 5e-5
    make_spiny_compt()
    make_plots()
    for i in range( 8 ):
        moose.setClock( i, fineDt )
    moose.setClock( 8, coarseDt )
    moose.reinit()
    moose.start( 0.1 )
    display_plots( 'instab.plot' )
    # make Hsolver and rerun
    hsolve = moose.HSolve( '/n/hsolve' )
    for i in range( 8 ):
        moose.setClock( i, coarseDt )
    hsolve.dt = coarseDt
    hsolve.target = '/n/compt'
    moose.reinit()
    moose.start( 0.1 )
    display_plots( 'h_instab.plot' )
    pylab.show()

def main():
		makeModel()
                '''
                '''
		ksolve = moose.Ksolve( '/model/compartment/ksolve' )
		stoich = moose.Stoich( '/model/compartment/stoich' )
		stoich.compartment = moose.element( '/model/compartment' )
		stoich.ksolve = ksolve
		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/ksolve', 'process' )
                '''

		moose.reinit()
		moose.start( 100.0 ) # Run the model for 100 seconds.
                func = moose.element( '/model/compartment/d/func' )
                if useY:
                    func.expr = "-y0 + 10*y1"
                else:
                    func.expr = "-x0 + 10*x1"
		moose.start( 100.0 ) # Run the model for 100 seconds.
                #moose.showfields( '/model/compartment/d' )
                #moose.showfields( '/model/compartment/d/func' )
                print func.x.value
                print moose.element( '/model/compartment/b' ).n

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

		quit()

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 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 main():
    """
    This shows the use of SynChan with Izhikevich neuron. This can be
    used for creating a network of Izhikevich neurons.
    """
    
    simtime = 200.0
    stepsize = 10.0
    model_dict = make_model()
    vm, inject, gk, spike = setup_data_recording(model_dict['neuron'],
                                          model_dict['pulse'],
                                          model_dict['synapse'],
                                          model_dict['spike_in'])
    mutils.setDefaultDt(elecdt=0.01, plotdt2=0.25)
    mutils.assignDefaultTicks(solver='ee')
    moose.reinit()
    mutils.stepRun(simtime, stepsize)
    pylab.subplot(411)
    pylab.plot(pylab.linspace(0, simtime, len(vm.vector)), vm.vector, label='Vm (mV)')
    pylab.legend()
    pylab.subplot(412)
    pylab.plot(pylab.linspace(0, simtime, len(inject.vector)), inject.vector, label='Inject (uA)')
    pylab.legend()
    pylab.subplot(413)
    pylab.plot(spike.vector, pylab.ones(len(spike.vector)), '|', label='input spike times')
    pylab.legend()
    pylab.subplot(414)
    pylab.plot(pylab.linspace(0, simtime, len(gk.vector)), gk.vector, label='Gk (mS)')
    pylab.legend()
    pylab.show()

def loadAndRun(solver=True):
    simtime = 500e-3
    model = moose.loadModel('../data/h10.CNG.swc', '/cell')
    comp = moose.element('/cell/apical_e_177_0')
    soma = moose.element('/cell/soma')
    for i in range(10):
        moose.setClock(i, dt)
    if solver:
        solver = moose.HSolve('/cell/solver')
        solver.target = soma.path
        solver.dt = dt
    stim = moose.PulseGen('/cell/stim')
    stim.delay[0] = 50e-3
    stim.delay[1] = 1e9
    stim.level[0] = 1e-9
    stim.width[0] = 2e-3    
    moose.connect(stim, 'output', comp, 'injectMsg')
    tab = moose.Table('/cell/Vm')
    moose.connect(tab, 'requestOut', comp, 'getVm')
    tab_soma = moose.Table('/cell/Vm_soma')
    moose.connect(tab_soma, 'requestOut', soma, 'getVm')
    moose.reinit()
    print('[INFO] Running for %s' % simtime)
    moose.start(simtime )
    vec = tab_soma.vector 
    moose.delete( '/cell' )
    return vec

def main():
    ######## Put your favourite cell model here ######
    ##This one is from PMID 22730554: Suo et al J. Mol Cell Biol 2012
    filename = 'cells/ko20x-07.CNG.swc'
    moose.Neutral( '/library' )
    rdes = rd.rdesigneur( \
            cellProto = [[ filename, 'elec' ] ],\
            passiveDistrib = passiveDistrib_,
            chanProto = chanProto_,
            chanDistrib = chanDistrib_
        )
    rdes.buildModel( '/model' )
    moose.reinit()

    ################## Now we store plots ########################
    somaVm = moose.Table( '/somaVm' )
    moose.connect( somaVm, 'requestOut', rdes.soma, 'getVm' )
    ################## Now we set up the display ########################
    compts = moose.wildcardFind( "/model/elec/#[ISA=CompartmentBase]" )
    compts[0].inject = inject

    print("Setting Up 3D Display")
    app = QtGui.QApplication(sys.argv)
    vm_viewer = create_vm_viewer(rdes, somaVm)
    vm_viewer.show()
    vm_viewer.start()
    return app.exec_()

def main():
    """ This example illustrates loading, running of an SBML model defined in XML format.\n
	The model 00001-sbml-l3v1.xml is taken from l3v1 SBML testcase.\n
	Plots are setup.\n
	Model is run for 20sec.\n
	As a general rule we created model under '/path/model' and plots under '/path/graphs'.\n
    """

    mfile =  os.path.join( script_dir, 'chem_models/00001-sbml-l3v1.xml')
    runtime = 20.0
        
    # Loading the sbml file into MOOSE, models are loaded in path/model
    sbmlId = moose.readSBML(mfile,'sbml')
    

    s1 = moose.element('/sbml/model/compartment/S1')
    s2= moose.element('/sbml/model/compartment/S2')
                      
    # Creating MOOSE Table, Table2 is for the chemical model
    graphs = moose.Neutral( '/sbml/graphs' )
    outputs1 = moose.Table2 ( '/sbml/graphs/concS1')
    outputs2 = moose.Table2 ( '/sbml/graphs/concS2')

    # connect up the tables
    moose.connect( outputs1,'requestOut', s1, 'getConc' );
    moose.connect( outputs2,'requestOut', s2, 'getConc' );

        
    # Reset and Run
    moose.reinit()
    moose.start(runtime)

def main():
        """
        This example illustrates loading and running a reaction system that
        spans two volumes, that is, is in different compartments. It uses a
        kkit model file. You can tell if it is working if you see nice 
        relaxation oscillations.
        """
        # the kkit reader doesn't know how to do multicompt solver setup.
        solver = "ee"  
	mfile = '../Genesis_files/OSC_diff_vols.g'
	runtime = 3000.0
        simDt = 1.0
	modelId = moose.loadModel( mfile, 'model', solver )
        #moose.delete( '/model/kinetics/A/Stot' )
        compt0 = moose.element( '/model/kinetics' )
        compt1 = moose.element( '/model/compartment_1' )
        assert( deq( compt0.volume, 2e-20 ) )
        assert( deq( compt1.volume, 1e-20 ) )
        dy = compt0.dy
        compt1.y1 += dy
        compt1.y0 = dy
        assert( deq( compt1.volume, 1e-20 ) )
        # We now have two cubes adjacent to each other. Compt0 has 2x vol.
        # Compt1 touches it.
        stoich0 = moose.Stoich( '/model/kinetics/stoich' )
        stoich1 = moose.Stoich( '/model/compartment_1/stoich' )
        ksolve0 = moose.Ksolve( '/model/kinetics/ksolve' )
        ksolve1 = moose.Ksolve( '/model/compartment_1/ksolve' )
        stoich0.compartment = compt0
        stoich0.ksolve = ksolve0
        stoich0.path = '/model/kinetics/##'
        stoich1.compartment = compt1
        stoich1.ksolve = ksolve1
        stoich1.path = '/model/compartment_1/##'
        #stoich0.buildXreacs( stoich1 )
        print ksolve0.numLocalVoxels, ksolve0.numPools, stoich0.numAllPools
        assert( ksolve0.numLocalVoxels == 1 )
        assert( ksolve0.numPools == 7 )
        assert( stoich0.numAllPools == 6 )
        print len( stoich0.proxyPools[stoich1] ),
        print len( stoich1.proxyPools[stoich0] )
        assert( len( stoich0.proxyPools[stoich1] ) == 1 )
        assert( len( stoich1.proxyPools[stoich0] ) == 1 )
        print ksolve1.numLocalVoxels, ksolve1.numPools, stoich1.numAllPools
        assert( ksolve1.numLocalVoxels == 1 )
        assert( ksolve1.numPools == 6 )
        assert( stoich1.numAllPools == 5 )
        stoich0.buildXreacs( stoich1 )
        print moose.element( '/model/kinetics/endo' )
        print moose.element( '/model/compartment_1/exo' )
        moose.le( '/model/compartment_1' )
	moose.reinit()
	moose.start( runtime ) 

	# Display all plots.
	for x in moose.wildcardFind( '/model/#graphs/conc#/#' ):
            t = numpy.arange( 0, x.vector.size, 1 ) * simDt
            pylab.plot( t, x.vector, label=x.name )
        pylab.legend()
        pylab.show()