From 720c59490bc254a393f89d7bd9a3ce33b3a49521 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Mon, 18 Aug 2025 14:37:09 +0200 Subject: [PATCH 1/7] Fix incorrect initialization for non-constant initial assignments Fixes a bug where targets of initial assignments are not initialized correctly if the initial assignment expression is implicitly time-dependent. Fixes #2936. --- python/sdist/amici/sbml_import.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/python/sdist/amici/sbml_import.py b/python/sdist/amici/sbml_import.py index 8d4db42e91..b379d43e85 100644 --- a/python/sdist/amici/sbml_import.py +++ b/python/sdist/amici/sbml_import.py @@ -1392,7 +1392,9 @@ def _process_parameters( and not ia.is_Number and not self.is_rate_rule_target(par) ): - if not ia.has(sbml_time_symbol): + if not ia.has(sbml_time_symbol) and not ( + ia.free_symbols - set(self.symbols[SymbolId.PARAMETER]) + ): self.symbols[SymbolId.EXPRESSION][ _get_identifier_symbol(par) ] = { @@ -1407,6 +1409,10 @@ def _process_parameters( # We can't represent that as expression, since the # initial simulation time is only known at the time of the # simulation, so we can't substitute it. + # Also, any parameter with an initial assignment + # that expression that is implicitly time-dependent + # must be converted to a species to avoid re-evaluating + # the initial assignment at every time step. self.symbols[SymbolId.SPECIES][ _get_identifier_symbol(par) ] = { From 27e9d1bf68b5761d2fa7ee0a47a70cfac21a2463 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Mon, 18 Aug 2025 16:11:30 +0200 Subject: [PATCH 2/7] .. --- python/sdist/amici/sbml_import.py | 28 +++++++++++++++++++++------- 1 file changed, 21 insertions(+), 7 deletions(-) diff --git a/python/sdist/amici/sbml_import.py b/python/sdist/amici/sbml_import.py index b379d43e85..6a200934fe 100644 --- a/python/sdist/amici/sbml_import.py +++ b/python/sdist/amici/sbml_import.py @@ -1386,15 +1386,20 @@ def _process_parameters( # Parameters that need to be turned into expressions or species # so far, this concerns parameters with symbolic initial assignments # (those have been skipped above) that are not rate rule targets + + # Set of symbols in initial assignments that still allows handling them + # via amici expressions + syms_allowed_in_expr_ia = set(self.symbols[SymbolId.PARAMETER]) | set( + self.symbols[SymbolId.PARAMETER] + ) + for par in self.sbml.getListOfParameters(): if ( (ia := par_id_to_ia.get(par.getId())) is not None and not ia.is_Number and not self.is_rate_rule_target(par) ): - if not ia.has(sbml_time_symbol) and not ( - ia.free_symbols - set(self.symbols[SymbolId.PARAMETER]) - ): + if not (ia.free_symbols - syms_allowed_in_expr_ia): self.symbols[SymbolId.EXPRESSION][ _get_identifier_symbol(par) ] = { @@ -1522,12 +1527,21 @@ def _process_rules(self) -> None: ) # expressions must not occur in definition of x0 + allowed_syms = ( + set(self.symbols[SymbolId.PARAMETER]) + | set(self.symbols[SymbolId.FIXED_PARAMETER]) + | {amici_time_symbol} + ) for species in self.symbols[SymbolId.SPECIES].values(): - species["init"] = self._make_initial( - smart_subs_dict( - species["init"], self.symbols[SymbolId.EXPRESSION], "value" + # only parameters are allowed as free symbols + while species["init"].free_symbols - allowed_syms: + species["init"] = self._make_initial( + smart_subs_dict( + species["init"], + self.symbols[SymbolId.EXPRESSION], + "value", + ) ) - ) def _process_rule_algebraic(self, rule: libsbml.AlgebraicRule): formula = self._sympify(rule) From 08c8244b9b091029e82bab798143bb71aab95489 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Mon, 18 Aug 2025 22:02:12 +0200 Subject: [PATCH 3/7] .. --- python/sdist/amici/de_model.py | 9 +++++++++ python/sdist/amici/sbml_import.py | 12 ++++++++++-- 2 files changed, 19 insertions(+), 2 deletions(-) diff --git a/python/sdist/amici/de_model.py b/python/sdist/amici/de_model.py index 4aaead9770..3be456e0cf 100644 --- a/python/sdist/amici/de_model.py +++ b/python/sdist/amici/de_model.py @@ -455,12 +455,21 @@ def get_rate(symbol: sp.Symbol): self._eqs["xdot"] = smart_subs_dict(self.eq("xdot"), subs) # replace rateOf-instances in x0 by xdot equation + made_substitutions = False for i_state in range(len(self.eq("x0"))): new, replacement = self._eqs["x0"][i_state].replace( rate_of_func, get_rate, map=True ) if replacement: self._eqs["x0"][i_state] = new + made_substitutions = True + if made_substitutions: + # replace any newly introduced state variables + # by their x0 expressions + subs = toposort_symbols( + dict(zip(self.sym("x_rdata"), self.eq("x0"), strict=True)) + ) + self._eqs["x0"] = smart_subs_dict(self.eq("x0"), subs) # replace rateOf-instances in w by xdot equation # here we may need toposort, as xdot may depend on w diff --git a/python/sdist/amici/sbml_import.py b/python/sdist/amici/sbml_import.py index 6a200934fe..cf300b499f 100644 --- a/python/sdist/amici/sbml_import.py +++ b/python/sdist/amici/sbml_import.py @@ -1526,7 +1526,7 @@ def _process_rules(self) -> None: self.symbols[SymbolId.EXPRESSION], "value" ) - # expressions must not occur in definition of x0 + # expressions must not occur in the definition of x0 allowed_syms = ( set(self.symbols[SymbolId.PARAMETER]) | set(self.symbols[SymbolId.FIXED_PARAMETER]) @@ -1534,7 +1534,15 @@ def _process_rules(self) -> None: ) for species in self.symbols[SymbolId.SPECIES].values(): # only parameters are allowed as free symbols - while species["init"].free_symbols - allowed_syms: + while True: + sym_math, rateof_to_dummy = _rateof_to_dummy(species["init"]) + if ( + sym_math.free_symbols + - allowed_syms + - set(rateof_to_dummy.values()) + == set() + ): + break species["init"] = self._make_initial( smart_subs_dict( species["init"], From 76dc01517da513060870b470294529e6ac459c55 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 19 Aug 2025 10:57:35 +0200 Subject: [PATCH 4/7] .. --- python/sdist/amici/sbml_import.py | 1 + 1 file changed, 1 insertion(+) diff --git a/python/sdist/amici/sbml_import.py b/python/sdist/amici/sbml_import.py index cf300b499f..eb12065d08 100644 --- a/python/sdist/amici/sbml_import.py +++ b/python/sdist/amici/sbml_import.py @@ -1533,6 +1533,7 @@ def _process_rules(self) -> None: | {amici_time_symbol} ) for species in self.symbols[SymbolId.SPECIES].values(): + species["init"] = species["init"].subs(self.compartments) # only parameters are allowed as free symbols while True: sym_math, rateof_to_dummy = _rateof_to_dummy(species["init"]) From c6b934a91dc2a4bf3ab144c31cb3ee5a39094e10 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 19 Aug 2025 13:38:40 +0200 Subject: [PATCH 5/7] .. --- pytest.ini | 2 +- python/sdist/amici/sbml_import.py | 10 +++++-- python/tests/test_sbml_import.py | 45 ++++++++++++++++++++++++++++--- 3 files changed, 51 insertions(+), 6 deletions(-) diff --git a/pytest.ini b/pytest.ini index 481b60e04f..5a478bc1f9 100644 --- a/pytest.ini +++ b/pytest.ini @@ -30,4 +30,4 @@ filterwarnings = ignore:jax.* is deprecated:DeprecationWarning -norecursedirs = .git amici_models build doc documentation matlab models ThirdParty amici sdist examples +norecursedirs = .git amici_models build doc documentation matlab models ThirdParty amici sdist examples *build* diff --git a/python/sdist/amici/sbml_import.py b/python/sdist/amici/sbml_import.py index eb12065d08..a217eeb813 100644 --- a/python/sdist/amici/sbml_import.py +++ b/python/sdist/amici/sbml_import.py @@ -1390,7 +1390,7 @@ def _process_parameters( # Set of symbols in initial assignments that still allows handling them # via amici expressions syms_allowed_in_expr_ia = set(self.symbols[SymbolId.PARAMETER]) | set( - self.symbols[SymbolId.PARAMETER] + self.symbols[SymbolId.FIXED_PARAMETER] ) for par in self.sbml.getListOfParameters(): @@ -1530,13 +1530,15 @@ def _process_rules(self) -> None: allowed_syms = ( set(self.symbols[SymbolId.PARAMETER]) | set(self.symbols[SymbolId.FIXED_PARAMETER]) - | {amici_time_symbol} + | {sbml_time_symbol} ) for species in self.symbols[SymbolId.SPECIES].values(): species["init"] = species["init"].subs(self.compartments) # only parameters are allowed as free symbols + old_init = None while True: sym_math, rateof_to_dummy = _rateof_to_dummy(species["init"]) + old_init = species["init"] if ( sym_math.free_symbols - allowed_syms @@ -1551,6 +1553,10 @@ def _process_rules(self) -> None: "value", ) ) + if species["init"] == old_init: + raise AssertionError( + "Infinite loop detected in _process_rules." + ) def _process_rule_algebraic(self, rule: libsbml.AlgebraicRule): formula = self._sympify(rule) diff --git a/python/tests/test_sbml_import.py b/python/tests/test_sbml_import.py index eaa6896fab..3f124f8e11 100644 --- a/python/tests/test_sbml_import.py +++ b/python/tests/test_sbml_import.py @@ -5,16 +5,16 @@ import sys from numbers import Number from pathlib import Path - import amici import libsbml import numpy as np import pytest from amici.gradient_check import check_derivatives -from amici.sbml_import import SbmlImporter -from amici.testing import skip_on_valgrind +from amici.sbml_import import SbmlImporter, SymbolId +from amici.import_utils import symbol_with_assumptions from numpy.testing import assert_allclose, assert_array_equal from amici import import_model_module +from amici.testing import skip_on_valgrind from amici.testing import TemporaryDirectoryWinSafe as TemporaryDirectory from conftest import MODEL_STEADYSTATE_SCALED_XML import sympy as sp @@ -1142,3 +1142,42 @@ def test_contains_periodic_subexpression(): assert cps(sp.sin(t), t) is True assert cps(sp.cos(t), t) is True assert cps(t + sp.sin(t), t) is True + + +@skip_on_valgrind +@pytest.mark.parametrize("compute_conservation_laws", [True, False]) +def test_time_dependent_initial_assignment(compute_conservation_laws: bool): + """Check that dynamic expressions for initial assignments are only + evaluated at t=t0.""" + from amici.antimony_import import antimony2sbml + + ant_model = """ + x1' = 1 + x1 = p0 + p0 = 1 + p1 = x1 + x2 := x1 + p2 = x2 + """ + + sbml_model = antimony2sbml(ant_model) + print(sbml_model) + si = SbmlImporter(sbml_model, from_file=False) + de_model = si._build_ode_model( + observables={"obs_p1": {"formula": "p1"}, "obs_p2": {"formula": "p2"}}, + compute_conservation_laws=compute_conservation_laws, + ) + # "species", because the initial assignment expression is time-dependent + assert symbol_with_assumptions("p2") in si.symbols[SymbolId.SPECIES].keys() + # "species", because differential state + assert symbol_with_assumptions("x1") in si.symbols[SymbolId.SPECIES].keys() + + assert "p0" in [str(p.get_id()) for p in de_model.parameters()] + assert "p1" not in [str(p.get_id()) for p in de_model.parameters()] + assert "p2" not in [str(p.get_id()) for p in de_model.parameters()] + + assert list(de_model.sym("x_rdata")) == [ + symbol_with_assumptions("p2"), + symbol_with_assumptions("x1"), + ] + assert list(de_model.eq("x0")) == [symbol_with_assumptions("p0")] * 2 From 0868fa604e8239a9b8e9e4182486d81951dbd4de Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 19 Aug 2025 14:57:25 +0200 Subject: [PATCH 6/7] .. --- python/sdist/amici/de_model.py | 37 +++++++++++++++++-------------- python/sdist/amici/sbml_import.py | 9 +++++--- 2 files changed, 26 insertions(+), 20 deletions(-) diff --git a/python/sdist/amici/de_model.py b/python/sdist/amici/de_model.py index 3be456e0cf..82cc81922c 100644 --- a/python/sdist/amici/de_model.py +++ b/python/sdist/amici/de_model.py @@ -454,23 +454,6 @@ def get_rate(symbol: sp.Symbol): ) self._eqs["xdot"] = smart_subs_dict(self.eq("xdot"), subs) - # replace rateOf-instances in x0 by xdot equation - made_substitutions = False - for i_state in range(len(self.eq("x0"))): - new, replacement = self._eqs["x0"][i_state].replace( - rate_of_func, get_rate, map=True - ) - if replacement: - self._eqs["x0"][i_state] = new - made_substitutions = True - if made_substitutions: - # replace any newly introduced state variables - # by their x0 expressions - subs = toposort_symbols( - dict(zip(self.sym("x_rdata"), self.eq("x0"), strict=True)) - ) - self._eqs["x0"] = smart_subs_dict(self.eq("x0"), subs) - # replace rateOf-instances in w by xdot equation # here we may need toposort, as xdot may depend on w made_substitutions = False @@ -518,6 +501,26 @@ def get_rate(symbol: sp.Symbol): self._syms["w"] = sp.Matrix(topo_expr_syms) self._eqs["w"] = sp.Matrix(list(w_sorted.values())) + # replace rateOf-instances in x0 by xdot equation + made_substitutions = False + for i_state in range(len(self.eq("x0"))): + new, replacement = self._eqs["x0"][i_state].replace( + rate_of_func, get_rate, map=True + ) + if replacement: + self._eqs["x0"][i_state] = new + made_substitutions = True + if made_substitutions: + # Replace any newly introduced state variables + # by their x0 expressions. + # Also replace any newly introduced `w` symbols by their + # expressions (after `w` was toposorted above). + subs = toposort_symbols( + dict(zip(self.sym("x_rdata"), self.eq("x0"), strict=True)) + ) + subs = dict(zip(self._syms["w"], self.eq("w"), strict=True)) | subs + self._eqs["x0"] = smart_subs_dict(self.eq("x0"), subs) + for component in chain( self.observables(), self.events(), diff --git a/python/sdist/amici/sbml_import.py b/python/sdist/amici/sbml_import.py index a217eeb813..ec827a5ed4 100644 --- a/python/sdist/amici/sbml_import.py +++ b/python/sdist/amici/sbml_import.py @@ -1533,10 +1533,9 @@ def _process_rules(self) -> None: | {sbml_time_symbol} ) for species in self.symbols[SymbolId.SPECIES].values(): - species["init"] = species["init"].subs(self.compartments) # only parameters are allowed as free symbols - old_init = None while True: + species["init"] = species["init"].subs(self.compartments) sym_math, rateof_to_dummy = _rateof_to_dummy(species["init"]) old_init = species["init"] if ( @@ -1555,7 +1554,7 @@ def _process_rules(self) -> None: ) if species["init"] == old_init: raise AssertionError( - "Infinite loop detected in _process_rules." + f"Infinite loop detected in _process_rules {species}." ) def _process_rule_algebraic(self, rule: libsbml.AlgebraicRule): @@ -2394,6 +2393,10 @@ def _make_initial( sym_math = sym_math.subs( var, self.symbols[SymbolId.SPECIES][var]["init"] ) + elif var in self.symbols[SymbolId.ALGEBRAIC_STATE]: + sym_math = sym_math.subs( + var, self.symbols[SymbolId.ALGEBRAIC_STATE][var]["init"] + ) elif ( element := self.sbml.getElementBySId(element_id) ) and self.is_rate_rule_target(element): From 6996b7e70e811a1c8dbef978d52fe1bff40c1e30 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 26 Aug 2025 09:54:17 +0200 Subject: [PATCH 7/7] only substitute in modified expressions --- python/sdist/amici/de_model.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/python/sdist/amici/de_model.py b/python/sdist/amici/de_model.py index 82cc81922c..538bd33f69 100644 --- a/python/sdist/amici/de_model.py +++ b/python/sdist/amici/de_model.py @@ -502,15 +502,16 @@ def get_rate(symbol: sp.Symbol): self._eqs["w"] = sp.Matrix(list(w_sorted.values())) # replace rateOf-instances in x0 by xdot equation - made_substitutions = False + # indices of state variables whose x0 was modified + changed_indices = [] for i_state in range(len(self.eq("x0"))): new, replacement = self._eqs["x0"][i_state].replace( rate_of_func, get_rate, map=True ) if replacement: self._eqs["x0"][i_state] = new - made_substitutions = True - if made_substitutions: + changed_indices.append(i_state) + if changed_indices: # Replace any newly introduced state variables # by their x0 expressions. # Also replace any newly introduced `w` symbols by their @@ -519,7 +520,10 @@ def get_rate(symbol: sp.Symbol): dict(zip(self.sym("x_rdata"), self.eq("x0"), strict=True)) ) subs = dict(zip(self._syms["w"], self.eq("w"), strict=True)) | subs - self._eqs["x0"] = smart_subs_dict(self.eq("x0"), subs) + for i_state in changed_indices: + self._eqs["x0"][i_state] = smart_subs_dict( + self._eqs["x0"][i_state], subs + ) for component in chain( self.observables(),