Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add fixed species support #364

Merged
merged 9 commits into from Feb 5, 2019
@@ -8,7 +8,6 @@
# exported from PySB model 'robertson'

import numpy
import weave
import scipy.integrate
import collections
import itertools
@@ -18,8 +17,11 @@
_use_inline = False
# try to inline a C statement to see if inline is functional
try:
import weave
weave.inline('int i;', force=1)
_use_inline = True
except ImportError:
pass
except distutils.errors.CompileError:
pass

@@ -41,7 +43,7 @@ def __init__(self):
self.sim_param_values = numpy.empty(6)
self.parameters = [None] * 6
self.observables = [None] * 3
self.initial_conditions = [None] * 3
self.initials = [None] * 3

self.parameters[0] = Parameter('k1', 0.040000000000000001)
self.parameters[1] = Parameter('k2', 30000000)
@@ -54,28 +56,28 @@ def __init__(self):
self.observables[1] = Observable('B_total', [1], [1])
self.observables[2] = Observable('C_total', [2], [1])

self.initial_conditions[0] = Initial(3, 0)
self.initial_conditions[1] = Initial(4, 1)
self.initial_conditions[2] = Initial(5, 2)
self.initials[0] = Initial(3, 0)
self.initials[1] = Initial(4, 1)
self.initials[2] = Initial(5, 2)

if _use_inline:

def ode_rhs(self, t, y, p):
ydot = self.ydot
weave.inline(r'''
ydot[0] = -p[0]*y[0] + p[2]*y[1]*y[2];
ydot[1] = p[0]*y[0] - p[1]*pow(y[1], 2) - p[2]*y[1]*y[2];
ydot[2] = p[1]*pow(y[1], 2);
ydot[0] = y[1]*y[2]*p[2] + (y[0]*p[0])*(-1);
ydot[1] = y[0]*p[0] + (pow(y[1], 2)*p[1])*(-1) + (y[1]*y[2]*p[2])*(-1);
ydot[2] = pow(y[1], 2)*p[1];
''', ['ydot', 't', 'y', 'p'])
return ydot

else:

def ode_rhs(self, t, y, p):
ydot = self.ydot
ydot[0] = -p[0]*y[0] + p[2]*y[1]*y[2]
ydot[1] = p[0]*y[0] - p[1]*pow(y[1], 2) - p[2]*y[1]*y[2]
ydot[2] = p[1]*pow(y[1], 2)
ydot[0] = y[1]*y[2]*p[2] + (y[0]*p[0])*(-1)
ydot[1] = y[0]*p[0] + (pow(y[1], 2)*p[1])*(-1) + (y[1]*y[2]*p[2])*(-1)
ydot[2] = pow(y[1], 2)*p[1]
return ydot


@@ -90,14 +92,14 @@ def simulate(self, tspan, param_values=None, view=False):
# create parameter vector from the values in the model
self.sim_param_values[:] = [p.value for p in self.parameters]
self.y0.fill(0)
for ic in self.initial_conditions:
for ic in self.initials:
self.y0[ic.species_index] = self.sim_param_values[ic.param_index]
if self.y is None or len(tspan) != len(self.y):
self.y = numpy.empty((len(tspan), len(self.y0)))
if len(self.observables):
self.yobs = numpy.ndarray(len(tspan),
zip((obs.name for obs in self.observables),
itertools.repeat(float)))
list(zip((obs.name for obs in self.observables),
itertools.repeat(float))))
else:
self.yobs = numpy.ndarray((len(tspan), 0))
self.yobs_view = self.yobs.view(float).reshape(len(self.yobs),
@@ -452,8 +452,8 @@ one shown (output shown below the ``'>>>'``` prompts)::
{'obsC8': <pysb.core.Observable object at 0x104b2c4d0>,
'obsBid': <pysb.core.Observable object at 0x104b2c5d0>,
'obstBid': <pysb.core.Observable object at 0x104b2c6d0>}
>>> m.model.initial_conditions
[(C8(b=None), Parameter(name='C8_0', value=1000)), (Bid(b=None, S=u), Parameter(name='Bid_0', value=10000))]
>>> m.model.initials
[Initial(C8(b=None), C8_0), Initial(Bid(b=None, S='u'), Bid_0)]
>>> m.model.rules
{'C8_Bid_bind': Rule(name='C8_Bid_bind', reactants=C8(b=None) +
Bid(b=None, S=None), products=C8(b=1) % Bid(b=1, S=None),
@@ -96,8 +96,10 @@ def _check_model(self):
"""
if not self.model.rules:
raise NoRulesError()
if not self.model.initial_conditions and not any(r.is_synth() for r
in self.model.rules):
if (
not self.model.initials
and not any(r.is_synth() for r in self.model.rules)
):
raise NoInitialConditionsError()

@classmethod
@@ -778,7 +780,7 @@ def _parse_species(model, line):
monomer_patterns = []
for ms in monomer_strings:
monomer_name, site_strings, monomer_compartment_name = \
re.match(r'(\w+)\(([^)]*)\)(?:@(\w+))?', ms).groups()
re.match(r'\$?(\w+)\(([^)]*)\)(?:@(\w+))?', ms).groups()
site_conditions = {}
if len(site_strings):
for ss in site_strings.split(','):
@@ -172,9 +172,11 @@ def expression(self, *args, **kwargs):
self.model.add_component(e)
return e

def initial(self, *args):
def initial(self, *args, **kwargs):
"""Adds an initial condition to the Builder's model instance."""
self.model.initial(*args)
i = Initial(*args, _export=False, **kwargs)
self.model.add_initial(i)
return i

def __getitem__(self, index):
"""Returns the component with the given string index
ProTip! Use n and p to navigate between commits in a pull request.
You can’t perform that action at this time.