Pythonize symbol parsing (#1296)

* refactor * format * simplify * fixup * fix code smells * refactor llh funs * fixup * refactor code + fix typehints * refactor species initial assignment * denest * fix initialization * fixup identifiers * fix doc format * make doc nicer * refactor * fixups * fixup * fixup compartment rate rules * refactor SBML testing * fixup dtype * fix assert * fix whitespaces * fixup parameter assignment rules with hasOnlySubstanceUnits * fixup compartment assignment rules * add all compartments to observables * Update sbml_import.py * Update sbml_import.py * fixup measurement symbols, species references, initial assignments * fixup * fixup, add support for parameter initial assignments * fix observable names, conversion factors and expose initial assignments * fixups reaction ids, initial assignments * fix float initial assignments * fixup species references * fixup species references * denest * add documentation * Update python_interface.rst * Apply suggestions from code review Co-authored-by: Daniel Weindl <[email protected]> * address review comments * fixup * pythonize species parsing * pythonization of remaining symbols * fix symbol resetting * implement enum for fields, fix symbol resetting * fix doc * don't creat duplicate initial assignments in petab import * cleanup bad merge * fixup make initials * fixup * sanitize petab output ... * fixup sbml testsuite * fix observables * fixup * fixup * refactor field names * pythonize compartments and get rid of rate rule dicts * fixup * add support for parameter rate rules, fixes #1278 * fixup * fixup * format + efficiency * fix volume conversion and reserved symbol replacement * change handling of 'amount' species, do not simulate concentrations * add more compute time tracking * fixup initials and comp time derivative * iterable constant parameters for backwards compatibility (#1303) * Apply suggestions from code review Co-authored-by: dilpath <[email protected]> * fix compartments * Apply suggestions from code review Co-authored-by: dilpath <[email protected]> * refactor & adress review comments * fixups * fixup * fixup volume/concentration, add usage hint with pytest * fix code smell * update testsuite stats * better check for species in stoichiometric matrix * Apply suggestions from code review Co-authored-by: dilpath <[email protected]> * fix local symbols * update SBML test statistic * update doc * Update python/amici/sbml_import.py Co-authored-by: Daniel Weindl <[email protected]> * move code comment to doctstring Co-authored-by: Daniel Weindl <[email protected]> Co-authored-by: dilpath <[email protected]>
AMICI-dev · Oct 20, 2020 · 9da5cfd · 9da5cfd
1 parent 49134ae
commit 9da5cfd
Show file tree

Hide file tree

Showing 8 changed files with 527 additions and 641 deletions.
diff --git a/.gitignore b/.gitignore
@@ -188,3 +188,4 @@ documentation/generated/*
 tests/performance/CS_Signalling_ERBB_RAS_AKT/
 tests/performance/CS_Signalling_ERBB_RAS_AKT_petab/*
 coverage_SBMLSuite.xml
+Benchmark-Models-PEtab/
diff --git a/documentation/python_interface.rst b/documentation/python_interface.rst
@@ -25,7 +25,7 @@ AMICI can import :term:`SBML` models via the
 Status of SBML support in Python-AMICI
 ++++++++++++++++++++++++++++++++++++++
 
-Python-AMICI currently **passes 725 out of the 1780 (~41%) test cases** from
+Python-AMICI currently **passes 750 out of the 1780 (~42%) test cases** from
 the semantic
 `SBML Test Suite <https://github.com/sbmlteam/sbml-test-suite/>`_
 (`current status <https://github.com/AMICI-dev/AMICI/actions>`_).
@@ -44,7 +44,6 @@ However, the following features are unlikely to be supported:
 - SBML extensions
 - `factorial()`, `ceil()`, `floor()`, due to incompatibility with
   symbolic sensitivity computations
-- initial assignments for parameters
 - `delay()` due to missing SUNDIALS solver support
 
 

diff --git a/python/amici/constants.py b/python/amici/constants.py
@@ -0,0 +1,28 @@
+"""
+Constants
+-----------
+This module provides a central place to define native python enums and
+constants that are used in multiple other modules
+"""
+
+import enum
+
+
+class SymbolId(str, enum.Enum):
+    """
+    Defines the different fields in the symbol dict to which sbml entities
+    get parsed to.
+
+    .. note:: This class inherits from str enabling direct comparison to
+        strings, which means that the species symbols can be accessed as
+        symbols['species'], which is convenient for debugging and symbols[
+        SymbolId.SPECIES], which is how the field should be accessed
+        programmatically.
+    """
+    SPECIES = 'species'
+    PARAMETER = 'parameter'
+    FIXED_PARAMETER = 'fixed_parameter'
+    OBSERVABLE = 'observable'
+    EXPRESSION = 'expression'
+    SIGMAY = 'sigmay'
+    LLHY = 'llhy'
diff --git a/python/amici/ode_export.py b/python/amici/ode_export.py
@@ -4,9 +4,9 @@
 This module provides all necessary functionality specify an ODE model and
 generate executable C++ simulation code. The user generally won't have to
 directly call any function from this module as this will be done by
-:func:`amici.pysb_import.pysb2amici`,
-:meth:`amici.sbml_import.SbmlImporter.sbml2amici` and
-:func:`amici.petab_import.import_model`
+:py:func:`amici.pysb_import.pysb2amici`,
+:py:func:`amici.sbml_import.SbmlImporter.sbml2amici` and
+:py:func:`amici.petab_import.import_model`
 """
 import sympy as sp
 import numpy as np
@@ -40,6 +40,7 @@
     sbml_import
 )
 from .logging import get_logger, log_execution_time, set_log_level
+from .constants import SymbolId
 
 # Template for model simulation main.cpp file
 CXX_MAIN_TEMPLATE_FILE = os.path.join(amiciSrcPath, 'main.template.cpp')
@@ -275,7 +276,7 @@ class ModelQuantity:
     def __init__(self,
                  identifier: sp.Symbol,
                  name: str,
-                 value: Union[SupportsFloat, numbers.Number, sp.Basic]):
+                 value: Union[SupportsFloat, numbers.Number, sp.Expr]):
         """
         Create a new ModelQuantity instance.
 
@@ -301,8 +302,8 @@ def __init__(self,
         if isinstance(value, sp.RealNumber) \
                 or isinstance(value, numbers.Number):
             value = float(value)
-        if not isinstance(value, sp.Basic) and not isinstance(value, float):
-            raise TypeError(f'value must be sympy.Symbol or float, was '
+        if not isinstance(value, sp.Expr) and not isinstance(value, float):
+            raise TypeError(f'value must be sympy.Expr or float, was '
                             f'{type(value)}')
         self._value = value
 
@@ -363,8 +364,8 @@ class State(ModelQuantity):
     def __init__(self,
                  identifier: sp.Symbol,
                  name: str,
-                 value: sp.Basic,
-                 dt: sp.Basic):
+                 init: sp.Expr,
+                 dt: sp.Expr):
         """
         Create a new State instance. Extends :meth:`ModelQuantity.__init__`
         by dt
@@ -375,13 +376,13 @@ def __init__(self,
         :param name:
             individual name of the state (does not need to be unique)
 
-        :param value:
+        :param init:
             initial value
 
         :param dt:
             time derivative
         """
-        super(State, self).__init__(identifier, name, value)
+        super(State, self).__init__(identifier, name, init)
         if not isinstance(dt, sp.Expr):
             raise TypeError(f'dt must have type sympy.Expr, was '
                             f'{type(dt)}')
@@ -641,13 +642,13 @@ def __init__(self,
 
 # defines the type of some attributes in ODEModel
 symbol_to_type = {
-    'species': State,
-    'parameter': Parameter,
-    'fixed_parameter': Constant,
-    'observable': Observable,
-    'sigmay': SigmaY,
-    'llhy': LogLikelihood,
-    'expression': Expression,
+    SymbolId.SPECIES: State,
+    SymbolId.PARAMETER: Parameter,
+    SymbolId.FIXED_PARAMETER: Constant,
+    SymbolId.OBSERVABLE: Observable,
+    SymbolId.SIGMAY: SigmaY,
+    SymbolId.LLHY: LogLikelihood,
+    SymbolId.EXPRESSION: Expression,
 }
 
 
@@ -921,32 +922,21 @@ def import_from_sbml_importer(self,
 
         dxdotdw_updates = []
 
-        def dx_dt(x_index, dxdt):
+        def transform_dxdt_to_concentration(specie_id, dxdt):
             """
             Produces the appropriate expression for the first derivative of a
             species with respect to time, for species that reside in
             compartments with a constant volume, or a volume that is defined by
             an assignment or rate rule.
 
-            :param x_index:
-                The index (not identifier) of the species in the variables
-                (generated in "sbml_import.py") that describe the model.
+            :param specie_id:
+                The identifier of the species (generated in "sbml_import.py").
 
             :param dxdt:
                 The element-wise product of the row in the stoichiometric
                 matrix that corresponds to the species (row x_index) and the
                 flux (kinetic laws) vector. Ignored in the case of rate rules.
             """
-            x_id = symbols['species']['identifier'][x_index]
-
-            # Rate rules specify dx_dt.
-            # Note that the rate rule of species may describe amount, not
-            # concentration.
-            if x_id in si.compartment_rate_rules:
-                return si.compartment_rate_rules[x_id]
-            elif x_id in si.species_rate_rules:
-                return si.species_rate_rules[x_id]
-
             # The derivation of the below return expressions can be found in
             # the documentation. They are found by rearranging
             # $\frac{d}{dt} (vx) = Sw$ for $\frac{dx}{dt}$, where $v$ is the
@@ -956,49 +946,87 @@ def dx_dt(x_index, dxdt):
             # species in (i) compartments with a rate rule, (ii) compartments
             # with an assignment rule, and (iii) compartments with a constant
             # volume, respectively.
-            v_name = si.species_compartment[x_index]
-            if v_name in si.compartment_rate_rules:
-                dv_dt = si.compartment_rate_rules[v_name]
-                xdot = (dxdt - dv_dt * x_id) / v_name
-                for w_index, flux in enumerate(fluxes):
-                    dxdotdw_updates.append((x_index, w_index, xdot.diff(flux)))
+            specie = si.symbols[SymbolId.SPECIES][specie_id]
+
+            comp = specie['compartment']
+            x_index = specie['index']
+            if comp in si.symbols[SymbolId.SPECIES]:
+                dv_dt = si.symbols[SymbolId.SPECIES][comp]['dt']
+                xdot = (dxdt - dv_dt * specie_id) / comp
+                dxdotdw_updates.extend(
+                    (x_index, w_index, xdot.diff(r_flux))
+                    for w_index, r_flux in enumerate(fluxes)
+                )
                 return xdot
-            elif v_name in si.compartment_assignment_rules:
-                v = si.compartment_assignment_rules[v_name]
+            elif comp in si.compartment_assignment_rules:
+                v = si.compartment_assignment_rules[comp]
                 dv_dt = v.diff(si.amici_time_symbol)
-                dv_dx = v.diff(x_id)
-                xdot = (dxdt - dv_dt * x_id) / (dv_dx * x_id + v)
-                for w_index, flux in enumerate(fluxes):
-                    dxdotdw_updates.append((x_index, w_index, xdot.diff(flux)))
+                # we may end up with a time derivative of the compartment
+                # volume due to parameter rate rules
+                comp_rate_vars = [p for p in v.free_symbols
+                                  if p in si.symbols[SymbolId.SPECIES]]
+                for var in comp_rate_vars:
+                    dv_dt += \
+                        v.diff(var) * si.symbols[SymbolId.SPECIES][var]['dt']
+                dv_dx = v.diff(specie_id)
+                xdot = (dxdt - dv_dt * specie_id) / (dv_dx * specie_id + v)
+                dxdotdw_updates.extend(
+                    (x_index, w_index, xdot.diff(r_flux))
+                    for w_index, r_flux in enumerate(fluxes)
+                )
                 return xdot
             else:
-                v = si.compartment_volume[list(si.compartment_symbols).index(
-                    si.species_compartment[x_index])]
+                v = si.compartments[comp]
 
                 if v == 1.0:
                     return dxdt
 
-                for w_index, flux in enumerate(fluxes):
-                    if si.stoichiometric_matrix[x_index, w_index] != 0:
-                        dxdotdw_updates.append((x_index, w_index,
-                            si.stoichiometric_matrix[x_index, w_index] / v))
+                dxdotdw_updates.extend(
+                    (x_index, w_index,
+                     si.stoichiometric_matrix[x_index, w_index] / v)
+                    for w_index in range(si.stoichiometric_matrix.shape[1])
+                    if si.stoichiometric_matrix[x_index, w_index] != 0
+                )
+
                 return dxdt / v
 
         # create dynamics without respecting conservation laws first
         dxdt = smart_multiply(si.stoichiometric_matrix,
                               MutableDenseMatrix(fluxes))
-        symbols['species']['dt'] = sp.Matrix([
-            dx_dt(x_index, dxdt[x_index])
-            for x_index in range(dxdt.rows)
-        ])
+        for ix, ((specie_id, specie), formula) in enumerate(zip(
+                symbols[SymbolId.SPECIES].items(),
+                dxdt
+        )):
+            assert ix == specie['index']  # check that no reordering occurred
+            # rate rules and amount species don't need to be updated
+            if 'dt' in specie:
+                continue
+            if specie['amount']:
+                specie['dt'] = formula
+            else:
+                specie['dt'] = transform_dxdt_to_concentration(specie_id,
+                                                               formula)
 
         # create all basic components of the ODE model and add them.
-        for symbol in [s for s in symbols if s != 'my']:
+        for symbol_name in symbols:
             # transform dict of lists into a list of dicts
-            protos = [dict(zip(symbols[symbol], t))
-                      for t in zip(*symbols[symbol].values())]
+            args = ['name', 'identifier']
+
+            if symbol_name == SymbolId.SPECIES:
+                args += ['dt', 'init']
+            else:
+                args += ['value']
+
+            protos = [
+                {
+                    'identifier': var_id,
+                    **{k: v for k, v in var.items() if k in args}
+                }
+                for var_id, var in symbols[symbol_name].items()
+            ]
+
             for proto in protos:
-                self.add_component(symbol_to_type[symbol](**proto))
+                self.add_component(symbol_to_type[symbol_name](**proto))
 
         # process conservation laws
         if compute_cls:

diff --git a/python/amici/petab_import.py b/python/amici/petab_import.py
@@ -10,6 +10,7 @@
 import logging
 import math
 import os
+import re
 import shutil
 import tempfile
 from _collections import OrderedDict
@@ -526,8 +527,10 @@ def import_model_sbml(
             init_par = sbml_model.createParameter()
             init_par.setId(init_par_id)
             init_par.setName(init_par_id)
-        assignment = sbml_model.createInitialAssignment()
-        assignment.setSymbol(assignee_id)
+        assignment = sbml_model.getInitialAssignment(assignee_id)
+        if assignment is None:
+            assignment = sbml_model.createInitialAssignment()
+            assignment.setSymbol(assignee_id)
         formula = f'{PREEQ_INDICATOR_ID} * {init_par_id_preeq} ' \
                   f'+ (1 - {PREEQ_INDICATOR_ID}) * {init_par_id_sim}'
         math_ast = libsbml.parseL3Formula(formula)
@@ -584,9 +587,11 @@ def get_observation_model(observable_df: pd.DataFrame
 
     for _, observable in observable_df.iterrows():
         oid = observable.name
-        name = observable.get(OBSERVABLE_NAME, "")
-        formula_obs = observable[OBSERVABLE_FORMULA]
-        formula_noise = observable[NOISE_FORMULA]
+        # need to sanitize due to https://github.com/PEtab-dev/PEtab/issues/447
+        pat = r'^[nN]a[nN]$'
+        name = re.sub(pat, '', str(observable.get(OBSERVABLE_NAME, '')))
+        formula_obs = re.sub(pat, '', str(observable[OBSERVABLE_FORMULA]))
+        formula_noise = re.sub(pat, '', str(observable[NOISE_FORMULA]))
         observables[oid] = {'name': name, 'formula': formula_obs}
         sigmas[oid] = formula_noise