From 5f8824fac95ded346a6035d0b31432a287f3f727 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 28 Aug 2025 12:08:03 +0200 Subject: [PATCH] Fix handling of initial assignments containing splines Currently, `x0` must not depend on splines. Therefore, substitute the spline symbol by it piecewise expression when compiling x0 expressions. Reported by @m-philipps. --- python/sdist/amici/sbml_import.py | 16 ++++++++++++++-- python/tests/test_sbml_import.py | 25 ++++++++++++++++++++++--- 2 files changed, 36 insertions(+), 5 deletions(-) diff --git a/python/sdist/amici/sbml_import.py b/python/sdist/amici/sbml_import.py index ec827a5ed4..1cdba3ae8c 100644 --- a/python/sdist/amici/sbml_import.py +++ b/python/sdist/amici/sbml_import.py @@ -1526,14 +1526,14 @@ def _process_rules(self) -> None: self.symbols[SymbolId.EXPRESSION], "value" ) - # expressions must not occur in the definition of x0 + # substitute symbols that must not occur in the definition of x0 + # allowed symbols: amici model parameters and time allowed_syms = ( set(self.symbols[SymbolId.PARAMETER]) | set(self.symbols[SymbolId.FIXED_PARAMETER]) | {sbml_time_symbol} ) for species in self.symbols[SymbolId.SPECIES].values(): - # only parameters are allowed as free symbols while True: species["init"] = species["init"].subs(self.compartments) sym_math, rateof_to_dummy = _rateof_to_dummy(species["init"]) @@ -2403,6 +2403,18 @@ def _make_initial( # no need to recurse here, as value is numeric init = sp.Float(element.getValue()) sym_math = sym_math.subs(var, init) + elif spline := [spl for spl in self.splines if spl.sbml_id == var]: + # x0 must not depend on splines -- substitute + assert len(spline) == 1 + spline = spline[0] + if spline.evaluate_at != sbml_time_symbol: + raise NotImplementedError( + "AMICI at the moment does not support splines " + "whose evaluation point is not the model time." + ) + sym_math = sym_math.subs( + var, spline.evaluate(sbml_time_symbol) + ) sym_math = _dummy_to_rateof(sym_math, rateof_to_dummy) diff --git a/python/tests/test_sbml_import.py b/python/tests/test_sbml_import.py index 3f124f8e11..f72773fc86 100644 --- a/python/tests/test_sbml_import.py +++ b/python/tests/test_sbml_import.py @@ -1150,6 +1150,7 @@ 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 + from amici.import_utils import amici_time_symbol ant_model = """ x1' = 1 @@ -1158,10 +1159,23 @@ def test_time_dependent_initial_assignment(compute_conservation_laws: bool): p1 = x1 x2 := x1 p2 = x2 + spline1 = 0 # replaced by actual spline below + x3' = 0 + x3 = spline1 """ + sbml_str = antimony2sbml(ant_model) + sbml_reader = libsbml.SBMLReader() + sbml_document = sbml_reader.readSBMLFromString(sbml_str) + sbml_model = sbml_document.getModel() + + spline = amici.splines.CubicHermiteSpline( + sbml_id="spline1", + nodes=[0, 1, 2], + values_at_nodes=[3, 4, 5], + ) + spline.add_to_sbml_model(sbml_model) + sbml_model.getElementBySId("spline1").setConstant(False) - 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"}}, @@ -1179,5 +1193,10 @@ def test_time_dependent_initial_assignment(compute_conservation_laws: bool): assert list(de_model.sym("x_rdata")) == [ symbol_with_assumptions("p2"), symbol_with_assumptions("x1"), + symbol_with_assumptions("x3"), + ] + assert list(de_model.eq("x0")) == [ + symbol_with_assumptions("p0"), + symbol_with_assumptions("p0"), + amici_time_symbol * 1.0 + 3.0, ] - assert list(de_model.eq("x0")) == [symbol_with_assumptions("p0")] * 2