"""
Module containing a class for converting a PySB model to a set of ordinary
differential equations for integration or analysis in Mathematica.
For information on how to use the model exporters, see the documentation
for :py:mod:`pysb.export`.
Output for the Robertson example model
======================================
The Mathematica code produced will follow the form as given below for
``pysb.examples.robertson``::
(*
A simple three-species chemical kinetics system known as "Robertson's
example", as presented in:
H. H. Robertson, The solution of a set of reaction rate equations, in Numerical
Analysis: An Introduction, J. Walsh, ed., Academic Press, 1966, pp. 178-182.
Mathematica model definition file for model robertson.
Generated by pysb.export.mathematica.MathematicaExporter.
Run with (for example):
tmax = 10
soln = NDSolve[Join[odes, initconds], slist, {t, 0, tmax}]
Plot[s0[t] /. soln, {t, 0, tmax}, PlotRange -> All]
*)
(* Parameters *)
k1 = 0.040000000000000001;
k2 = 30000000;
k3 = 10000;
A0 = 1;
B0 = 0;
C0 = 0;
(* List of Species *)
(* s0[t] = A() *)
(* s1[t] = B() *)
(* s2[t] = C() *)
(* ODEs *)
odes = {
s0'[t] == -k1*s0[t] + k3*s1[t]*s2[t],
s1'[t] == k1*s0[t] - k2*s1[t]^2 - k3*s1[t]*s2[t],
s2'[t] == k2*s1[t]^2
}
(* Initial Conditions *)
initconds = {
s0[0] == A0,
s1[0] == B0,
s2[0] == C0
}
(* List of Variables (e.g., as an argument to NDSolve) *)
solvelist = {
s0[t],
s1[t],
s2[t]
}
(* Run the simulation -- example *)
tmax = 100
soln = NDSolve[Join[odes, initconds], solvelist, {t, 0, tmax}]
(* Observables *)
Atotal = (s0[t] * 1) /. soln
Btotal = (s1[t] * 1) /. soln
Ctotal = (s2[t] * 1) /. soln
The output consists of a block of commands that define the ODEs, parameters,
species and other variables for the model, along with a set of descriptive
comments. The sections are as follows:
* The header comments identify the model and show an example of how to
integrate the ODEs in Mathematica.
* The parameters block defines the numerical values of the named parameters.
* The list of species gives the mapping between the indexed species (``s0``,
``s1``, ``s2``) and their representation in PySB (``A()``, ``B()``, ``C()``).
* The ODEs block defines the set of ordinary differential equations and assigns
the set of equations to the variable ``odes``.
* The initial conditions block defines the initial values for each species and
assigns the set of conditions to the variable ``initconds``.
* The "list of variables" block enumerates all of the species in the model
(``s0[t]``, ``s1[t]``, ``s2[t]``) and assigns them to the variable
``solvelist``; this list can be passed to the Mathematica command ``NDSolve``
to indicate the variables to be solved for.
* This is followed by an example of how to call ``NDSolve`` to integrate the
equations.
* Finally, the observables block enumerates the observables in the model,
expressing each one as a linear combination of the appropriate species in
the model. The interpolating functions returned by ``NDSolve`` are substituted
in from the solution variable ``soln``, allowing the observables to be
plotted.
Note that Mathematica does not permit underscores in variable names, so
any underscores used in PySB variables will be removed (e.g., ``A_total`` will
be converted to ``Atotal``).
"""
import pysb
import pysb.bng
import sympy
import re
from io import StringIO
from pysb.export import Exporter, ExpressionsNotSupported, \
CompartmentsNotSupported
[docs]class MathematicaExporter(Exporter):
"""A class for returning the ODEs for a given PySB model for use in
Mathematica.
Inherits from :py:class:`pysb.export.Exporter`, which implements
basic functionality for all exporters.
"""
[docs] def export(self):
"""Generate the corresponding Mathematica ODEs for the PySB model
associated with the exporter.
Returns
-------
string
String containing the Mathematica code for the model's ODEs.
"""
if self.model.expressions:
raise ExpressionsNotSupported()
if self.model.compartments:
raise CompartmentsNotSupported()
output = StringIO()
pysb.bng.generate_equations(self.model)
# Add docstring if there is one
if self.docstring:
output.write('(*\n' + self.docstring + '\n')
else:
output.write("(*\n")
# Header comment
output.write("Mathematica model definition file for ")
output.write("model " + self.model.name + ".\n")
output.write("Generated by " \
"pysb.export.mathematica.MathematicaExporter.\n")
output.write("\n")
output.write("Run with (for example):\n")
output.write("tmax = 10\n")
output.write("soln = NDSolve[Join[odes, initconds], slist, " \
"{t, 0, tmax}]\n")
output.write("Plot[s0[t] /. soln, {t, 0, tmax}, PlotRange -> All]\n")
output.write("*)\n\n")
# PARAMETERS
# Note that in Mathematica, underscores are not allowed in variable
# names, so we simply strip them out here
params_str = ''
for i, p in enumerate(self.model.parameters):
# Remove underscores
pname = p.name.replace('_', '')
# Convert parameter values to scientific notation
# If the parameter is 0, don't take the log!
if p.value == 0:
params_str += '%s = %g;\n' % (pname, p.value)
# Otherwise, take the log (base 10) and format accordingly
else:
val_str = '%.17g' % p.value
if 'e' in val_str:
(mantissa, exponent) = val_str.split('e')
params_str += '%s = %s * 10^%s;\n' % \
(pname, mantissa, exponent)
else:
params_str += '%s = %s;\n' % (pname, val_str)
## ODEs ###
odes_str = 'odes = {\n'
# Concatenate the equations
odes_str += ',\n'.join(['s%d == %s' %
(i, sympy.ccode(self.model.odes[i]))
for i in range(len(self.model.odes))])
# Replace, e.g., s0 with s[0]
odes_str = re.sub(r's(\d+)', lambda m: 's%s[t]' % (int(m.group(1))),
odes_str)
# Add the derivative symbol ' to the left hand sides
odes_str = re.sub(r's(\d+)\[t\] ==', r"s\1'[t] ==", odes_str)
# Correct the exponentiation syntax
odes_str = re.sub(r'pow\(([^,]+), ([^)]+)\)', r'\1^\2', odes_str)
odes_str += '\n}'
#c_code = odes_str
# Eliminate underscores from parameter names in equations
for i, p in enumerate(self.model.parameters):
odes_str = re.sub(r'\b(%s)\b' % p.name, p.name.replace('_', ''),
odes_str)
## INITIAL CONDITIONS
ic_values = ['0'] * len(self.model.odes)
for i, ic in enumerate(self.model.initials):
idx = self.model.get_species_index(ic.pattern)
ic_values[idx] = ic.value.name.replace('_', '')
init_conds_str = 'initconds = {\n'
init_conds_str += ',\n'.join(['s%s[0] == %s' % (i, val)
for i, val in enumerate(ic_values)])
init_conds_str += '\n}'
## SOLVE LIST
solvelist_str = 'solvelist = {\n'
solvelist_str += ',\n'.join(['s%s[t]' % (i)
for i in range(0, len(self.model.odes))])
solvelist_str += '\n}'
## OBSERVABLES
observables_str = ''
for obs in self.model.observables:
# Remove underscores
observables_str += obs.name.replace('_', '') + ' = '
#groups = self.model.observable_groups[obs_name]
observables_str += ' + '.join(['(s%s[t] * %d)' % (s, c)
for s, c in zip(obs.species, obs.coefficients)])
observables_str += ' /. soln\n'
# Add comments identifying the species
species_str = '\n'.join(['(* s%d[t] = %s *)' % (i, s) for i, s in
enumerate(self.model.species)])
output.write('(* Parameters *)\n')
output.write(params_str + "\n")
output.write('(* List of Species *)\n')
output.write(species_str + "\n\n")
output.write('(* ODEs *)\n')
output.write(odes_str + "\n\n")
output.write('(* Initial Conditions *)\n')
output.write(init_conds_str + "\n\n")
output.write('(* List of Variables (e.g., as an argument to NDSolve) ' \
'*)\n')
output.write(solvelist_str + '\n\n')
output.write('(* Run the simulation -- example *)\n')
output.write('tmax = 100\n')
output.write('soln = NDSolve[Join[odes, initconds], ' \
'solvelist, {t, 0, tmax}]\n\n')
output.write('(* Observables *)\n')
output.write(observables_str + '\n')
return output.getvalue()