Merge pull request #153 from johnbachman/fix_expressions

Fix expressions

pysb/core.py

@@ -786,9 +786,9 @@ def func(self):
 
     def __repr__(self):
         return  '%s(%s, %s)' % (self.__class__.__name__, repr(self.name), repr(self.value))
-    
+
     def __str__(self):
-        return  '%s(%s, %s)' % (self.__class__.__name__, repr(self.name), repr(self.value))
+        return  repr(self)
 
 
 
@@ -844,7 +844,6 @@ def __repr__(self):
             (self.__class__.__name__, repr(self.name), repr(self.parent), repr(self.dimension), repr(self.size))
 
 
-
 class Rule(Component):
 
     """
@@ -1010,15 +1009,9 @@ def __repr__(self):
             ret += ', match=%s' % repr(self.match)
         ret += ')'
         return ret
-
-    def __str__(self):
-        ret = '%s(%s, %s' % (self.__class__.__name__, repr(self.name),
-                              repr(self.reaction_pattern))
-        if self.match != 'molecules':
-            ret += ', match=%s' % repr(self.match)
-        ret += ')'
-        return ret
 
+    def __str__(self):
+        return repr(self)
 
 
 class Expression(Component, sympy.Symbol):
@@ -1060,7 +1053,7 @@ def is_constant_expression(self):
                    (isinstance(a, Expression) and a.is_constant_expression()) or
                    isinstance(a, sympy.Number)
                    for a in self.expr.atoms())
-        
+
     def get_value(self):
         return self.expr.evalf()
 
@@ -1074,12 +1067,9 @@ def __repr__(self):
         ret = '%s(%s, %s)' % (self.__class__.__name__, repr(self.name),
                               repr(self.expr))
         return ret
-
-    def __str__(self):
-        ret = '%s(%s, %s)' % (self.__class__.__name__, repr(self.name),
-                              repr(self.expr))
-        return ret
 
+    def __str__(self):
+        return repr(self)
 
 
 class Model(object):

pysb/generator/bng.py

@@ -39,7 +39,7 @@ def generate_parameters(self):
                                (p.name, p.value))
         for e in exprs:
             self.__content += (("  %-" + str(max_length) + "s   %s\n") %
-                               (e.name, sympy_to_muparser(e.expr)))
+                               (e.name, expression_to_muparser(e)))
         self.__content += "end parameters\n\n"
 
     def generate_compartments(self):
@@ -121,7 +121,7 @@ def generate_functions(self):
         for i, e in enumerate(exprs):
             signature = e.name + '()'
             self.__content += ("  %-" + str(max_length) + "s   %s\n") % \
-                (signature, sympy_to_muparser(e.expr))
+                (signature, expression_to_muparser(e))
         self.__content += "end functions\n\n"
 
     def generate_species(self):
@@ -210,10 +210,14 @@ def warn_caller(message):
         module = inspect.getmodule(caller_frame)
     warnings.warn(message, stacklevel=stacklevel)
 
-def sympy_to_muparser(expr):
-#     code = sympy.fcode(expr)
-#     code = sympy.ccode(expr)
-    code = str(expr)
+def expression_to_muparser(expression):
+    """Render the Expression as a muparser-compatible string."""
+
+    # sympy.printing.sstr is the preferred way to render an Expression as a
+    # string (rather than, e.g., str(Expression.expr) or repr(Expression.expr).
+    # Note: "For large expressions where speed is a concern, use the setting
+    # order='none'"
+    code = sympy.printing.sstr(expression.expr, order='none')
     code = code.replace('\n     @', '')
     code = code.replace('**', '^')
     return code

pysb/integrate.py

@@ -335,8 +335,9 @@ def run(self, param_values=None, y0=None):
                     raise IndexError("param_values dictionary has unknown "
                                      "parameter name (%s)" % key)
                 param_values[pi] = param_values_dict[key]
-
-        subs = dict((p.name, param_values[i]) for i, p in
+
+        # The substitution dict must have Symbols as keys, not strings
+        subs = dict((p, param_values[i]) for i, p in
                     enumerate(self.model.parameters))
 
         if y0 is not None:
@@ -356,7 +357,7 @@ def run(self, param_values=None, y0=None):
                     pi = self.model.parameters.index(value_obj)
                     value = param_values[pi]
                 elif value_obj in self.model.expressions:
-                    value = value_obj.expand_expr(self.model).evalf(subs=subs)
+                    value = value_obj.expand_expr().evalf(subs=subs)
                 else:
                     raise ValueError("Unexpected initial condition value type")
                 si = self.model.get_species_index(cp)

pysb/tests/test_bng.py

@@ -41,4 +41,5 @@ def test_bidirectional_rules():
     ok_(len(model.reactions_bidirectional)==1)
     ok_(len(model.reactions_bidirectional[0]['rule'])==3)
     ok_(model.reactions_bidirectional[0]['reversible'])
-    #TODO Check that 'rate' has 4 terms
+    #TODO Check that 'rate' has 4 terms
+

pysb/tests/test_expressions.py

@@ -0,0 +1,59 @@
+from pysb import *
+from pysb.testing import *
+from pysb.bng import *
+from pysb.integrate import Solver
+import numpy as np
+
+@with_model
+def test_ic_expression_with_two_parameters():
+    Monomer('A')
+    Parameter('k1', 1)
+    Parameter('k2', 2)
+    Expression('e1', k1*k2)
+    Rule('A_deg', A() >> None, k1)
+    Initial(A(), e1)
+    generate_equations(model)
+    t = np.linspace(0, 1000, 100)
+    sol = Solver(model, t, use_analytic_jacobian=True)
+    sol.run()
+
+@with_model
+def test_ic_expression_with_one_parameter():
+    Monomer('A')
+    Parameter('k1', 1)
+    Expression('e1', k1)
+    Rule('A_deg', A() >> None, k1)
+    Initial(A(), e1)
+    generate_equations(model)
+    t = np.linspace(0, 1000, 100)
+    sol = Solver(model, t, use_analytic_jacobian=True)
+    sol.run()
+
+@with_model
+def test_expressions_with_one_observable():
+    Monomer('A')
+    Parameter('k1', 1)
+    Observable('o1', A())
+    Expression('e1', o1)
+    Rule('A_deg', A() >> None, k1)
+    Initial(A(), k1)
+    generate_equations(model)
+    t = np.linspace(0, 1000, 100)
+    sol = Solver(model, t, use_analytic_jacobian=True)
+    sol.run()
+
+@with_model
+def test_nested_expression():
+    Monomer('A')
+    Monomer('B')
+    Parameter('k1', 1)
+    Observable('o1', B())
+    Expression('e1', o1*10)
+    Expression('e2', e1+5)
+    Initial(A(), Parameter('A_0', 1000))
+    Initial(B(), Parameter('B_0', 1))
+    Rule('A_deg', A() >> None, e2)
+    generate_equations(model)
+    t = np.linspace(0, 1000, 100)
+    sol = Solver(model, t, use_analytic_jacobian=True)
+    sol.run()

pysb/tests/test_integrate.py

@@ -125,7 +125,7 @@ def test_integrate_with_expression():
     Observable('s16_obs', s16())
     Observable('s20_obs', s20())
 
-    Expression('keff', sympify("ka*ka20/(ka20+s9_obs)"))
+    Expression('keff', (ka*ka20)/(ka20+s9_obs))
 
     Rule('R1', None >> s16(), ka)
     Rule('R2', None >> s20(), ka)