#########################################################################
## This program is part of 'MOOSE', the
## Messaging Object Oriented Simulation Environment.
## Copyright (C) 2014 Upinder S. Bhalla. and NCBS
## It is made available under the terms of the
## GNU Lesser General Public License version 2.1
## See the file COPYING.LIB for the full notice.
#########################################################################
import moose
import pylab
import numpy
import matplotlib.pyplot as plt
import sys
def runAndSavePlots( name ):
runtime = 20.0
moose.reinit()
moose.start( runtime )
pa = moose.Neutral( '/model/graphs/' + name )
for x in moose.wildcardFind( '/model/#graphs/conc#/#' ):
if ( x.tick != -1 ):
tabname = '/model/graphs/' + name + '/' + x.name + '.' + name
y = moose.Table( tabname )
y.vector = x.vector
y.tick = -1
# Takes args ee, gsl, or gssa
def switchSolvers( solver ):
if ( moose.exists( 'model/kinetics/stoich' ) ):
moose.delete( '/model/kinetics/stoich' )
moose.delete( '/model/kinetics/ksolve' )
compt = moose.element( '/model/kinetics' )
if ( solver == 'gsl' ):
ksolve = moose.Ksolve( '/model/kinetics/ksolve' )
if ( solver == 'gssa' ):
ksolve = moose.Gsolve( '/model/kinetics/ksolve' )
if ( solver != 'ee' ):
stoich = moose.Stoich( '/model/kinetics/stoich' )
stoich.compartment = compt
stoich.ksolve = ksolve
stoich.path = "/model/kinetics/##"
[docs]def main():
"""
At zero order, you can select the solver you want to use within the
function moose.loadModel( filename, modelpath, solver ).
Having loaded in the model, you can change the solver to use on it.
This example illustrates how to assign and change solvers for a
kinetic model. This process is necessary in two situations:
* If we want to change the numerical method employed, for example,
from deterministic to stochastic.
* If we are already using a solver, and we have changed the reaction
network by adding or removing molecules or reactions.
Note that we do not have to change the solvers if the volume or
reaction rates change.
In this example the model is loaded in with a gsl solver. The
sequence of solver calculations is:
#. gsl
#. ee
#. gsl
#. gssa
#. gsl
If you're removing the solvers, you just delete the stoichiometry
object and the associated ksolve/gsolve. Should there be diffusion
(a dsolve)then you should delete that too. If you're
building the solvers up again, then you must do the following
steps in order:
#. build up the ksolve/gsolve and stoich (any order)
#. Assign stoich.ksolve
#. Assign stoich.path.
See the Reaction-diffusion section should you want to do diffusion
as well.
"""
solver = "gsl" # Pick any of gsl, gssa, ee..
mfile = '../genesis/kkit_objects_example.g'
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 = 1e-19
runAndSavePlots( 'gsl' )
#########################################################
switchSolvers( 'ee' )
runAndSavePlots( 'ee' )
#########################################################
switchSolvers( 'gsl' )
runAndSavePlots( 'gsl2' )
#########################################################
switchSolvers( 'gssa' )
runAndSavePlots( 'gssa' )
#########################################################
switchSolvers( 'gsl' )
runAndSavePlots( 'gsl3' )
#########################################################
# Display all plots.
fig = plt.figure( figsize = (12, 10) )
orig = fig.add_subplot( 511 )
gsl = fig.add_subplot( 512 )
ee = fig.add_subplot( 513 )
gsl2 = fig.add_subplot( 514 )
gssa = fig.add_subplot( 515 )
plotdt = moose.element( '/clock' ).tickDt[18]
for x in moose.wildcardFind( '/model/#graphs/conc#/#' ):
t = numpy.arange( 0, x.vector.size, 1 ) * plotdt
orig.plot( t, x.vector, label=x.name )
for x in moose.wildcardFind( '/model/graphs/gsl/#' ):
t = numpy.arange( 0, x.vector.size, 1 ) * plotdt
gsl.plot( t, x.vector, label=x.name )
for x in moose.wildcardFind( '/model/graphs/ee/#' ):
t = numpy.arange( 0, x.vector.size, 1 ) * plotdt
ee.plot( t, x.vector, label=x.name )
for x in moose.wildcardFind( '/model/graphs/gsl2/#' ):
t = numpy.arange( 0, x.vector.size, 1 ) * plotdt
gsl2.plot( t, x.vector, label=x.name )
for x in moose.wildcardFind( '/model/graphs/gssa/#' ):
t = numpy.arange( 0, x.vector.size, 1 ) * plotdt
gssa.plot( t, x.vector, label=x.name )
plt.legend()
pylab.show()
quit()
# Run the 'main' if this script is executed standalone.
if __name__ == '__main__':
main()