No event trigger reinitialization after pre-equilibration and pre-simulation

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

include/amici/model.h

@@ -298,6 +298,22 @@ class Model : public AbstractModel, public ModelDimensions {
         std::vector<int>& 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<realtype> const& h_old, std::vector<int>& roots_found
+    );
+
     /**
      * @brief Re-compute the explicit roots.
      *

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]

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])
 

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,

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
             );
         }
     }

src/model.cpp

@@ -421,16 +421,30 @@ void Model::initializeStateSensitivities(
 void Model::initEvents(
     realtype t, AmiVector const& x, AmiVector const& dx,
     std::vector<int>& roots_found
+) {
+    // construct "old" Heaviside vector from initial values
+    std::vector<realtype> 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<realtype> const& h_old, std::vector<int>& roots_found
 ) {
     std::vector<realtype> 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;
             }

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;