From de72ad975aec9b98af76be1eea06db37f6d912c5 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Mon, 21 Jul 2025 10:17:40 +0200 Subject: [PATCH] No event trigger reinitialization after pre-equilibration and pre-simulation Meanwhile, events are handled during pre-equilibration and pre-simulation. So far, SBML's `event.trigger.initialValue` was applied at the beginning of each period. However, this should only be applied at the beginning of the very first simulation period. This is fixed and tested here. Closes #2918. --- doc/glossary.rst | 3 + include/amici/model.h | 16 ++++++ python/tests/test_events.py | 81 +++++++++++++++++++++++++++ python/tests/test_preequilibration.py | 8 +-- python/tests/test_sbml_import.py | 20 +------ src/forwardproblem.cpp | 12 +++- src/model.cpp | 18 +++++- swig/model.i | 1 + 8 files changed, 133 insertions(+), 26 deletions(-) diff --git a/doc/glossary.rst b/doc/glossary.rst index 39767080cc..942ab048f6 100644 --- a/doc/glossary.rst +++ b/doc/glossary.rst @@ -93,6 +93,9 @@ Glossary be specified for each period, except for post-equilibration, which always uses the fixed parameters of the main simulation period. + For SBML models, initial assignments and initial values of event + triggers are only applied once at the start of the first period. + steady state In AMICI, a model is considered to be at steady state if the time derivative of the state variables is zero, up to the specified diff --git a/include/amici/model.h b/include/amici/model.h index d60222442a..5a032d765f 100644 --- a/include/amici/model.h +++ b/include/amici/model.h @@ -298,6 +298,22 @@ class Model : public AbstractModel, public ModelDimensions { std::vector& roots_found ); + /** + * @brief Re-initialize the Heaviside variables `h` at the beginning + * of the second or subsequent simulation. + * + * @param t Timepoint + * @param x Reference to state variables + * @param dx Reference to time derivative of states (DAE only) + * @param h_old Previous values of Heaviside variables + * @param roots_found boolean indicators indicating whether roots were found + * at t0 by this fun + */ + void reinit_events( + realtype t, AmiVector const& x, AmiVector const& dx, + std::vector const& h_old, std::vector& roots_found + ); + /** * @brief Re-compute the explicit roots. * diff --git a/python/tests/test_events.py b/python/tests/test_events.py index 76509713a3..edd29413fe 100644 --- a/python/tests/test_events.py +++ b/python/tests/test_events.py @@ -1180,3 +1180,84 @@ def test_posteq_events_are_handled(tempdir): epsilon=1e-8, skip_fields=["res"], ) + + +@skip_on_valgrind +def test_preeq_presim_preserve_heaviside_state(tempdir): + """Test the state of event trigger functions is preserved after + pre-equilibration and presimulation. + + I.e., the trigger.initialValue is only applied in the very beginning. + """ + from amici.antimony_import import antimony2amici + + model_name = "test_preeq_presim_preserve_heaviside_state" + antimony2amici( + r""" + some_time = 0 + some_time' = piecewise(1, (time <= 10), 0) + + target1 = 0 + target2 = 0 + + # random parameter to be able to enable + # pre-equilibration and presimulation + k1 = 42 + + # this should only trigger at the beginning of the first period + E1: at true, t0=false: + target1 = target1 + 1; + # this will only trigger at the beginning of the second period + # if the trigger is not re-initialized to `true` + E2: at some_time >= 10 and time < 1, t0 = true: + target2 = target2 + 1; + """, + constant_parameters=["k1"], + model_name=model_name, + output_dir=tempdir, + ) + + model_module = import_model_module(model_name, tempdir) + model = model_module.get_model() + model.setTimepoints(np.linspace(0, 2, 3)) + model.setSteadyStateComputationMode( + amici.SteadyStateComputationMode.integrationOnly + ) + solver = model.getSolver() + + # Only main simulation. E1 triggers, E2 does not. + rdata = amici.runAmiciSimulation(model, solver) + assert rdata.status == amici.AMICI_SUCCESS + assert list(rdata.by_id("target1")) == [1.0, 1.0, 1.0] + assert list(rdata.by_id("target2")) == [0.0, 0.0, 0.0] + + # Pre-equilibration + main simulation. Both E1 and E2 trigger. + edata = amici.ExpData(rdata, 1, 0, 0) + edata.fixedParametersPreequilibration = [1] + rdata = amici.runAmiciSimulation(model, solver, edata=edata) + assert rdata.status == amici.AMICI_SUCCESS + assert list(rdata.by_id("target1")) == [1.0, 1.0, 1.0] + assert list(rdata.by_id("target2")) == [1.0, 1.0, 1.0] + + # Pre-simulation + main simulation. Both E1 and E2 trigger. + # FIXME: this is currently not supported + # (root-after-reinitialization when switching from pre-simulation# + # to main simulation) + # edata = amici.ExpData(rdata, 1, 0, 0) + # edata.fixedParametersPresimulation = [1] + # edata.t_presim = 10 + # rdata = amici.runAmiciSimulation(model, solver, edata=edata) + # assert rdata.status == amici.AMICI_SUCCESS + # assert list(rdata.by_id("target1")) == [1.0, 1.0, 1.0] + # assert list(rdata.by_id("target2")) == [1.0, 1.0, 1.0] + + # Pre-equilibration + pre-simulation + main simulation. + # Both E1 and E2 trigger. + edata = amici.ExpData(rdata, 1, 0, 0) + edata.fixedParametersPreequilibration = [1] + edata.fixedParametersPresimulation = [1] + edata.t_presim = 10 + rdata = amici.runAmiciSimulation(model, solver, edata=edata) + assert rdata.status == amici.AMICI_SUCCESS + assert list(rdata.by_id("target1")) == [1.0, 1.0, 1.0] + assert list(rdata.by_id("target2")) == [1.0, 1.0, 1.0] diff --git a/python/tests/test_preequilibration.py b/python/tests/test_preequilibration.py index afe3fcba47..49422fbe0a 100644 --- a/python/tests/test_preequilibration.py +++ b/python/tests/test_preequilibration.py @@ -761,8 +761,8 @@ def test_preequilibration_events(tempdir): bolus2 = 1 bolus3 = 1 bolus4 = 1 - # E1 & E2 will both trigger during pre-equilibration and main - # simulation (Heaviside is reset after pre-equilibration) + # E1 & E2 will only trigger during pre-equilibration + # (initialValue is only used once at the start of pre-equilibration) E1: at some_time >= 0, t0 = false: target1 = target1 + bolus1 E2: at time >= 0, t0 = false: target2 = target2 + bolus2 # requires early time point @@ -815,8 +815,8 @@ def test_preequilibration_events(tempdir): assert rdata.x_ss[target2_idx] == 1 assert rdata.x_ss[target3_idx] == 1 assert rdata.x_ss[target4_idx] == 1 - assert np.all(rdata.x[:, target1_idx] == [2, 2]) - assert np.all(rdata.x[:, target2_idx] == [2, 2]) + assert np.all(rdata.x[:, target1_idx] == [1, 1]) + assert np.all(rdata.x[:, target2_idx] == [1, 1]) assert np.all(rdata.x[:, target3_idx] == [1, 1]) assert np.all(rdata.x[:, target4_idx] == [1, 2]) diff --git a/python/tests/test_sbml_import.py b/python/tests/test_sbml_import.py index 09e6ad4cc1..eaa6896fab 100644 --- a/python/tests/test_sbml_import.py +++ b/python/tests/test_sbml_import.py @@ -378,9 +378,7 @@ def test_presimulation_events_and_sensitivities(tempdir): xx = 0 xx' = piecewise(k_pre, time < 0, k_main) - # this will trigger twice, once in presimulation - # and once in the main simulation - at time >= -1 , t0=false: some_time = some_time + bolus + at time < -1 , t0=false: some_time = some_time + bolus """, model_name=model_name, output_dir=tempdir, @@ -412,22 +410,10 @@ def test_presimulation_events_and_sensitivities(tempdir): assert rdata.status == amici.AMICI_SUCCESS assert_allclose( - rdata.by_id("some_time"), np.array([0, 1, 2]) + 2.1, atol=1e-14 + rdata.by_id("some_time"), np.array([0, 1, 2]) + 1.1, atol=1e-14 ) - if sensi_method == amici.SensitivityMethod.forward: - model.requireSensitivitiesForAllParameters() - else: - # FIXME ASA with events: - # https://github.com/AMICI-dev/AMICI/pull/1539 - model.setParameterList( - [ - i - for i, p in enumerate(model.getParameterIds()) - if p != "bolus" - ] - ) - + model.requireSensitivitiesForAllParameters() check_derivatives( model, solver, diff --git a/src/forwardproblem.cpp b/src/forwardproblem.cpp index 9f00707113..52f7a4786f 100644 --- a/src/forwardproblem.cpp +++ b/src/forwardproblem.cpp @@ -297,7 +297,11 @@ void ForwardProblem::handlePresimulation() { ws_.roots_found ); } else if (model->ne) { - model->initEvents(ws_.sol.t, ws_.sol.x, ws_.sol.dx, ws_.roots_found); + // copy, since model state will be updated in reinit_events + auto h_old = model->getModelState().h; + model->reinit_events( + ws_.sol.t, ws_.sol.x, ws_.sol.dx, h_old, ws_.roots_found + ); } solver->setup(ws_.sol.t, model, ws_.sol.x, ws_.sol.dx, ws_.sol.sx, ws_.sdx); solver->updateAndReinitStatesAndSensitivities(model); @@ -338,8 +342,10 @@ void ForwardProblem::handleMainSimulation() { // Reset the time and re-initialize events for the main simulation solver->updateAndReinitStatesAndSensitivities(model); if (model->ne) { - model->initEvents( - model->t0(), ws_.sol.x, ws_.sol.dx, ws_.roots_found + // copy, since model state will be updated in reinit_events + auto h_old = model->getModelState().h; + model->reinit_events( + ws_.sol.t, ws_.sol.x, ws_.sol.dx, h_old, ws_.roots_found ); } } diff --git a/src/model.cpp b/src/model.cpp index f193931ccb..fc29c52ec5 100644 --- a/src/model.cpp +++ b/src/model.cpp @@ -421,16 +421,30 @@ void Model::initializeStateSensitivities( void Model::initEvents( realtype t, AmiVector const& x, AmiVector const& dx, std::vector& roots_found +) { + // construct "old" Heaviside vector from initial values + std::vector h_old(ne, 1.0); + for (int ie = 0; ie < ne; ie++) { + if (events_.at(ie).get_initial_value() == false) + h_old.at(ie) = 0.0; + } + + reinit_events(t, x, dx, h_old, roots_found); +} + +void Model::reinit_events( + realtype t, AmiVector const& x, AmiVector const& dx, + std::vector const& h_old, std::vector& roots_found ) { std::vector rootvals(ne, 0.0); froot(t, x, dx, rootvals); std::ranges::fill(roots_found, 0); for (int ie = 0; ie < ne; ie++) { - if (rootvals.at(ie) < 0) { + if (rootvals.at(ie) < 0.0) { state_.h.at(ie) = 0.0; } else { state_.h.at(ie) = 1.0; - if (pythonGenerated && !events_.at(ie).get_initial_value()) { + if (h_old.at(ie) <= 0.0) { // only false->true triggers event roots_found.at(ie) = 1; } diff --git a/swig/model.i b/swig/model.i index eb0ad4d7a2..6ad594c315 100644 --- a/swig/model.i +++ b/swig/model.i @@ -34,6 +34,7 @@ using namespace amici; %ignore getObservableSensitivity; %ignore getExpression; %ignore initEvents; +%ignore reinit_events; %ignore initialize; %ignore initializeB; %ignore initializeStateSensitivities;