Refactor ScipyOdeSimulator code generation (#531)

This moves all codegen out of the simulator class. A new abstract class
RhsBuilder is introduced to manage codegen, and Python and Cython concrete
implementations are provided. The codegen itself is now performed using sympy
expression manipulation and sympy's own codegen rather than hand-rolled string
manipulations and Cython.inline. Performance is improved since Cython.inline's
overhead on long code strings was really hurting us.

pysb/core.py

@@ -1964,6 +1964,10 @@ def all_components(self):
     def components(self):
         return self.all_components()
 
+    def parameters_all(self):
+        """Return a ComponentSet of all parameters and derived parameters."""
+        return self.parameters | self._derived_parameters
+
     def parameters_rules(self):
         """Return a ComponentSet of the parameters used in rules."""
         # rate_reverse is None for irreversible rules, so we'll need to filter those out
@@ -1999,16 +2003,21 @@ def parameters_unused(self):
         cset_used = (self.parameters_rules() | self.parameters_initial_conditions() |
                      self.parameters_compartments() | self.parameters_expressions())
         return self.parameters - cset_used
-    
-    def expressions_constant(self):
+
+    def expressions_constant(self, include_derived=False):
         """Return a ComponentSet of constant expressions."""
-        cset = ComponentSet(e for e in self.expressions
-                            if e.is_constant_expression())
+        expressions = self.expressions
+        if include_derived:
+            expressions = expressions | self._derived_expressions
+        cset = ComponentSet(e for e in expressions if e.is_constant_expression())
         return cset
 
-    def expressions_dynamic(self, include_local=True):
+    def expressions_dynamic(self, include_local=True, include_derived=False):
         """Return a ComponentSet of non-constant expressions."""
-        cset = self.expressions - self.expressions_constant()
+        expressions = self.expressions
+        if include_derived:
+            expressions = expressions | self._derived_expressions
+        cset = expressions - self.expressions_constant(include_derived)
         if not include_local:
             cset = ComponentSet(e for e in cset if not e.is_local)
         return cset

pysb/examples/run_earm_1_0.py

@@ -10,13 +10,15 @@
 
 
 # saturating level of ligand (corresponding to ~1000 ng/ml SuperKiller TRAIL)
-Lsat = 6E4; 
+Lsat = 6E4;
 
 # relationship of ligand concentration in the model (in # molecules/cell) to actual TRAIL concentration (in ng/ml)
 Lfactor = model.parameters['L_0'].value / 50;
 
 L_0_baseline = model.parameters['L_0'].value
 
+sim = ScipyOdeSimulator(model)
+
 
 def fig_4a():
     print("Simulating model for figure 4A...")
@@ -31,17 +33,15 @@ def fig_4a():
 
     CVenv = 0.2
     # num steps was originally 40, but 15 is plenty smooth enough for screen display
-    Ls = floor(logspace(1,5,15)) 
+    Ls = floor(logspace(1,5,15))
 
     fs = empty_like(Ls)
     Ts = empty_like(Ls)
     Td = empty_like(Ls)
     print("Scanning over %d values of L_0" % len(Ls))
-    for i in range(len(Ls)):
-        model.parameters['L_0'].value = Ls[i]
-
-        print("  L_0 = %g" % Ls[i])
-        x = ScipyOdeSimulator(model).run(tspan=t).all
+    for i, L_0 in enumerate(Ls):
+        print("  L_0 = %g" % L_0)
+        x = sim.run(tspan=t, param_values={"L_0": L_0}).all
 
         fs[i] = (x['PARP_unbound'][0] - x['PARP_unbound'][-1]) / x['PARP_unbound'][0]
         dP = 60 * (x['PARP_unbound'][:-1] - x['PARP_unbound'][1:]) / (dt * x['PARP_unbound'][0])  # in minutes
@@ -68,7 +68,7 @@ def fig_4b():
     print("Simulating model for figure 4B...")
 
     t = linspace(0, 6*3600, 6*60+1)  # 6 hours
-    x = ScipyOdeSimulator(model).run(tspan=t).all
+    x = sim.run(tspan=t).all
 
     x_norm = c_[x['Bid_unbound'], x['PARP_unbound'], x['mSmac_unbound']]
     x_norm = 1 - x_norm / x_norm[0, :]  # gets away without max() since first values are largest

pysb/logging.py

@@ -191,10 +191,10 @@ def get_logger(logger_name=BASE_LOGGER_NAME, model=None, log_level=None,
     if model is None:
         return logger
     else:
-        return PySBModelLoggerAdapter(logger, {'model': model})
+        return PySBModelLoggerAdapter(logger, {'model_name': model.name})
 
 
 class PySBModelLoggerAdapter(logging.LoggerAdapter):
     """ A logging adapter to prepend a model's name to log entries """
     def process(self, msg, kwargs):
-        return '[%s] %s' % (self.extra['model'].name, msg), kwargs
+        return '[%s] %s' % (self.extra['model_name'], msg), kwargs

pysb/simulator/base.py

@@ -190,7 +190,7 @@ def _update_initials_dict(self, initials_dict, initials_source, subs=None):
                         and all(isinstance(v, numbers.Number) for v in value_obj):
                     value = value_obj
                 elif isinstance(value_obj, Expression):
-                    value = [value_obj.expand_expr().evalf(subs=subs[sim]) for sim in range(len(subs))]
+                    value = [value_obj.expand_expr().xreplace(subs[sim]) for sim in range(len(subs))]
                 elif isinstance(value_obj, Parameter):
                     # Set parameter using param_values
                     pi = self._model.parameters.index(value_obj)