Refactor sympification (#261)

Centralize sympification of PEtab and SBML math expressions.

Simplifies incorporating #260

Also fixes likely collisions of SBML model entities with sympy entities.

Author: dweindl

petab/calculate.py

@@ -7,11 +7,11 @@
 import numpy as np
 import pandas as pd
 import sympy
-from sympy.abc import _clash
 
 import petab
 
 from .C import *
+from .math import sympify_petab
 
 __all__ = [
     "calculate_residuals",
@@ -157,7 +157,7 @@ def get_symbolic_noise_formulas(observable_df) -> Dict[str, sympy.Expr]:
         if NOISE_FORMULA not in observable_df.columns:
             noise_formula = None
         else:
-            noise_formula = sympy.sympify(row.noiseFormula, locals=_clash)
+            noise_formula = sympify_petab(row.noiseFormula)
         noise_formulas[observable_id] = noise_formula
     return noise_formulas
 

petab/lint.py

@@ -10,12 +10,12 @@
 import numpy as np
 import pandas as pd
 import sympy as sp
-from sympy.abc import _clash
 
 import petab
 
 from . import core, measurements, parameters
 from .C import *  # noqa: F403
+from .math import sympify_petab
 from .models import Model
 
 logger = logging.getLogger(__name__)
@@ -317,7 +317,7 @@ def check_observable_df(observable_df: pd.DataFrame) -> None:
     for row in observable_df.itertuples():
         obs = getattr(row, OBSERVABLE_FORMULA)
         try:
-            sp.sympify(obs, locals=_clash)
+            sympify_petab(obs)
         except sp.SympifyError as e:
             raise AssertionError(
                 f"Cannot parse expression '{obs}' "
@@ -326,7 +326,7 @@ def check_observable_df(observable_df: pd.DataFrame) -> None:
 
         noise = getattr(row, NOISE_FORMULA)
         try:
-            sympified_noise = sp.sympify(noise, locals=_clash)
+            sympified_noise = sympify_petab(noise)
             if sympified_noise is None or (
                 sympified_noise.is_Number and not sympified_noise.is_finite
             ):

petab/math/__init__.py

@@ -0,0 +1,2 @@
+"""Functions for parsing and evaluating mathematical expressions."""
+from .sympify import sympify_petab  # noqa: F401

petab/math/sympify.py

@@ -0,0 +1,20 @@
+"""PEtab math to sympy conversion."""
+
+import sympy as sp
+from sympy.abc import _clash
+
+
+def sympify_petab(expr: str) -> sp.Expr:
+    """
+    Convert a PEtab math expression to a sympy expression.
+
+    Parameters
+    ----------
+    expr:
+        The PEtab math expression.
+
+    Returns
+    -------
+    The sympy expression corresponding to ``expr``.
+    """
+    return sp.sympify(expr, locals=_clash)

petab/models/sbml_model.py

@@ -6,6 +6,7 @@
 
 import libsbml
 import sympy as sp
+from sympy.abc import _clash
 
 from ..sbml import (
     get_sbml_model,
@@ -104,34 +105,17 @@ def get_free_parameter_ids_with_values(
             ar.getVariable() for ar in self.sbml_model.getListOfRules()
         }
 
-        parser_settings = libsbml.L3ParserSettings(
-            self.sbml_model,
-            libsbml.L3P_PARSE_LOG_AS_LOG10,
-            libsbml.L3P_EXPAND_UNARY_MINUS,
-            libsbml.L3P_NO_UNITS,
-            libsbml.L3P_AVOGADRO_IS_CSYMBOL,
-            libsbml.L3P_COMPARE_BUILTINS_CASE_INSENSITIVE,
-            None,
-            libsbml.L3P_MODULO_IS_PIECEWISE,
-        )
-
         def get_initial(p):
             # return the initial assignment value if there is one, and it is a
             # number; `None`, if there is a non-numeric initial assignment;
             # otherwise, the parameter value
             if ia := self.sbml_model.getInitialAssignmentBySymbol(p.getId()):
-                formula_str = libsbml.formulaToL3StringWithSettings(
-                    ia.getMath(), parser_settings
+                sym_expr = sympify_sbml(ia.getMath())
+                return (
+                    float(sym_expr.evalf())
+                    if sym_expr.evalf().is_Number
+                    else None
                 )
-                try:
-                    return float(formula_str)
-                except ValueError:
-                    sym_expr = sp.sympify(formula_str)
-                    return (
-                        float(sym_expr.evalf())
-                        if sym_expr.evalf().is_Number
-                        else None
-                    )
             return p.getValue()
 
         return (
@@ -200,3 +184,39 @@ def is_state_variable(self, id_: str) -> bool:
             or self.sbml_model.getCompartment(id_) is not None
             or self.sbml_model.getRuleByVariable(id_) is not None
         )
+
+
+def sympify_sbml(sbml_obj: libsbml.ASTNode | libsbml.SBase) -> sp.Expr:
+    """Convert SBML math expression to sympy expression.
+
+    Parameters
+    ----------
+    sbml_obj:
+        SBML math element or an SBML object with a math element.
+
+    Returns
+    -------
+    The sympy expression corresponding to ``sbml_obj``.
+    """
+    ast_node = (
+        sbml_obj
+        if isinstance(sbml_obj, libsbml.ASTNode)
+        else sbml_obj.getMath()
+    )
+
+    parser_settings = libsbml.L3ParserSettings(
+        ast_node.getParentSBMLObject().getModel(),
+        libsbml.L3P_PARSE_LOG_AS_LOG10,
+        libsbml.L3P_EXPAND_UNARY_MINUS,
+        libsbml.L3P_NO_UNITS,
+        libsbml.L3P_AVOGADRO_IS_CSYMBOL,
+        libsbml.L3P_COMPARE_BUILTINS_CASE_INSENSITIVE,
+        None,
+        libsbml.L3P_MODULO_IS_PIECEWISE,
+    )
+
+    formula_str = libsbml.formulaToL3StringWithSettings(
+        ast_node, parser_settings
+    )
+
+    return sp.sympify(formula_str, locals=_clash)

petab/observables.py

@@ -6,11 +6,10 @@
 from typing import List, Literal, Union
 
 import pandas as pd
-import sympy as sp
-from sympy.abc import _clash
 
 from . import core, lint
 from .C import *  # noqa: F403
+from .math import sympify_petab
 from .models import Model
 
 __all__ = [
@@ -102,15 +101,15 @@ def get_output_parameters(
 
     for formula in formulas:
         free_syms = sorted(
-            sp.sympify(formula, locals=_clash).free_symbols,
+            sympify_petab(formula).free_symbols,
             key=lambda symbol: symbol.name,
         )
         for free_sym in free_syms:
             sym = str(free_sym)
             if model.symbol_allowed_in_observable_formula(sym):
                 continue
 
-            # does it mapping to a model entity?
+            # does it map to a model entity?
             if (
                 mapping_df is not None
                 and sym in mapping_df.index