From 010d1da2835ffc6fba81c45479861a64a410f45c Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 16 Oct 2025 10:36:42 +0200 Subject: [PATCH] Fix SBML import for parameters with piecewise initial assignments Parameters that are targets of parameter-dependent initial assignments are currently implemented as `w` expressions. Those symbols are currently not supported in event triggers. Previously, `w` symbols were only substituted in a subset of root functions, resulting in compilation failures due to undefined symbols for such models. Now `w` symbols are substituted in all trigger functions. Fixes PEtab v2 import for `Smith_BMCSystBiol2013`. Also, add `NULL` to the list of reserved symbols, fixing macro-redefinition compiler warnings when importing `Smith_BMCSystBiol2013`. --- python/sdist/amici/de_model.py | 68 +++++++++++++++++------------- python/sdist/amici/import_utils.py | 12 +++++- python/tests/test_bngl.py | 9 +++- 3 files changed, 58 insertions(+), 31 deletions(-) diff --git a/python/sdist/amici/de_model.py b/python/sdist/amici/de_model.py index e423e10ec0..e62f14312b 100644 --- a/python/sdist/amici/de_model.py +++ b/python/sdist/amici/de_model.py @@ -1305,18 +1305,36 @@ def parse_events(self) -> None: and replaces the formulae of the found roots by identifiers of AMICI's Heaviside function implementation in the right-hand side """ + # toposorted w_sym -> w_expr for substitution of 'w' in trigger function + # do only once. `w` is not modified during this function. + w_toposorted = toposort_symbols( + dict( + zip( + [expr.get_id() for expr in self._expressions], + [expr.get_val() for expr in self._expressions], + strict=True, + ) + ) + ) + # Track all roots functions in the right-hand side roots = copy.deepcopy(self._events) for state in self._differential_states: - state.set_dt(self._process_heavisides(state.get_dt(), roots)) + state.set_dt( + self._process_heavisides(state.get_dt(), roots, w_toposorted) + ) for expr in self._expressions: - expr.set_val(self._process_heavisides(expr.get_val(), roots)) + expr.set_val( + self._process_heavisides(expr.get_val(), roots, w_toposorted) + ) # remove all possible Heavisides from roots, which may arise from # the substitution of `'w'` in `_collect_heaviside_roots` for root in roots: - root.set_val(self._process_heavisides(root.get_val(), roots)) + root.set_val( + self._process_heavisides(root.get_val(), roots, w_toposorted) + ) # Now add the found roots to the model components for root in roots: @@ -1326,6 +1344,11 @@ def parse_events(self) -> None: # add roots of heaviside functions self.add_component(root) + # Substitute 'w' expressions into root expressions, to avoid rewriting + # 'root.cpp' and 'stau.cpp' headers to include 'w.h'. + for event in self.events(): + event.set_val(event.get_val().subs(w_toposorted)) + # re-order events - first those that require root tracking, then the others constant_syms = set(self.sym("k")) | set(self.sym("p")) self._events = list( @@ -2391,7 +2414,7 @@ def _expr_is_time_dependent(self, expr: sp.Expr) -> bool: expr_syms = {str(sym) for sym in expr.free_symbols} # Check if the time variable is in the expression. - if "t" in expr_syms: + if amici_time_symbol.name in expr_syms: return True # Check if any time-dependent states are in the expression. @@ -2464,33 +2487,11 @@ def _collect_heaviside_roots( return root_funs - def _substitute_w_in_roots( - self, - root_funs: list[tuple[sp.Expr, sp.Expr]], - ) -> list[tuple[sp.Expr, sp.Expr]]: - """ - Substitute 'w' expressions into root expressions, to avoid rewriting - 'root.cpp' and 'stau.cpp' headers to include 'w.h'. - """ - w_sorted = toposort_symbols( - dict( - zip( - [expr.get_id() for expr in self._expressions], - [expr.get_val() for expr in self._expressions], - strict=True, - ) - ) - ) - root_funs = [ - (r[0].subs(w_sorted), r[1].subs(w_sorted)) for r in root_funs - ] - - return root_funs - def _process_heavisides( self, dxdt: sp.Expr, roots: list[Event], + w_toposorted: dict[sp.Symbol, sp.Expr], ) -> sp.Expr: """ Parses the RHS of a state variable, checks for Heaviside functions, @@ -2502,7 +2503,8 @@ def _process_heavisides( right-hand side of state variable :param roots: list of known root functions with identifier - + :param w_toposorted: + `w` symbols->expressions sorted in topological order :returns: dxdt with Heaviside functions replaced by amici helper variables """ @@ -2511,7 +2513,15 @@ def _process_heavisides( heavisides = [] # run through the expression tree and get the roots tmp_roots_old = self._collect_heaviside_roots((dxdt,)) - tmp_roots_old = self._substitute_w_in_roots(tmp_roots_old) + # substitute 'w' symbols in the root expression by their equations, + # because currently, + # 1) root functions must not depend on 'w' + # 2) the check for time-dependence currently assumes only state + # variables are implicitly time-dependent + tmp_roots_old = [ + (a.subs(w_toposorted), b.subs(w_toposorted)) + for a, b in tmp_roots_old + ] for tmp_root_old, tmp_x0_old in unique_preserve_order(tmp_roots_old): # we want unique identifiers for the roots tmp_root_new = self._get_unique_root(tmp_root_old, roots) diff --git a/python/sdist/amici/import_utils.py b/python/sdist/amici/import_utils.py index 3a8e2ae61b..dc516362ae 100644 --- a/python/sdist/amici/import_utils.py +++ b/python/sdist/amici/import_utils.py @@ -17,7 +17,17 @@ from sympy.logic.boolalg import BooleanAtom from toposort import toposort -RESERVED_SYMBOLS = ["x", "k", "p", "y", "w", "h", "t", "AMICI_EMPTY_BOLUS"] +RESERVED_SYMBOLS = [ + "x", + "k", + "p", + "y", + "w", + "h", + "t", + "AMICI_EMPTY_BOLUS", + "NULL", +] try: import pysb diff --git a/python/tests/test_bngl.py b/python/tests/test_bngl.py index b92ea566a5..5baf20d055 100644 --- a/python/tests/test_bngl.py +++ b/python/tests/test_bngl.py @@ -1,6 +1,5 @@ import os -import amici import numpy as np import pytest @@ -10,6 +9,7 @@ from amici.testing import TemporaryDirectoryWinSafe, skip_on_valgrind from pysb.importers.bngl import model_from_bngl from pysb.simulator import ScipyOdeSimulator +from contextlib import suppress tests = [ "CaOscillate_Func", @@ -39,6 +39,13 @@ @skip_on_valgrind @pytest.mark.parametrize("example", tests) def test_compare_to_pysb_simulation(example): + import amici.import_utils + + # allow "NULL" as model symbol + # (used in CaOscillate_Func and Repressilator examples) + with suppress(ValueError): + amici.import_utils.RESERVED_SYMBOLS.remove("NULL") + atol = 1e-8 rtol = 1e-8