Fix Heaviside variables of roots crossed back by event assignments

changelog.d/+heaviside-crossed-back.fix.md

@@ -0,0 +1,3 @@
+Fixed wrong simulation results or failures if an event assignment reverts a
+trigger or `piecewise` condition that changed at the same time, e.g.,
+`at x < thr: x = 1.5 * x` (#3425).

include/amici/forwardproblem.h

@@ -242,8 +242,8 @@ struct FwdSimWorkspace {
     /** Events that are waiting to be handled at the current timepoint. */
     EventQueue pending_events;
 
-    /** temporary rootval storage to check crossing in secondary event
-     * (dimension: ne) */
+    /** root function values before the latest event assignment, to check
+     * which root functions it changed (dimension: ne) */
     std::vector<realtype> rval_tmp;
 
     /** array of number of found roots for a certain event type
@@ -421,12 +421,25 @@ class EventHandlingSimulator {
      * @brief Check for secondary events.
      *
      * Updates `roots_found` for the roots crossed due to the preceding event
-     * assignment.
+     * assignment, in either direction, including roots that were crossed
+     * earlier in the same discontinuity.
      *
      * @return the number of secondary events found
      */
     int detect_secondary_events();
 
+    /**
+     * @brief Estimate the change of the root functions within the integration
+     * tolerances at the current state.
+     *
+     * Computes `|g(x + d) - g(x)|` with `d_i = atol + rtol * |x_i|` (likewise
+     * for `dx`). Afterwards, re-evaluates the root functions at the current
+     * state, which also restores the model's derived state.
+     *
+     * @return the estimated change for each root function
+     */
+    std::vector<realtype> compute_root_noise();
+
     /**
      * @brief Add the events that triggered in the latest discontinuity to the
      * queue of pending events.

include/amici/misc.h

@@ -263,20 +263,6 @@ template <class T> bool is_equal(T const& a, T const& b) {
     return true;
 }
 
-/**
- * @brief Check whether the value of a Heaviside function differs for the two
- * arguments.
- *
- * Assumes H(x) = (x >= 1 ? 1 : 0).
- *
- * @param a First value
- * @param b Second value
- * @return H(a) != H(b).
- */
-template <typename T> bool heaviside_differs(T a, T b) {
-    return (a < T{0} && b >= T{0}) || (a >= T{0} && b < T{0});
-}
-
 /**
  * @brief Tracks elapsed CPU time of the calling thread.
  *

python/tests/test_events.py

@@ -3443,3 +3443,133 @@ def test_dae_explicit_and_solver_tracked_events(tempdir):
         assert_allclose(
             rdatas[dae].by_id("y"), 2 * rdatas[dae].by_id("x"), rtol=1e-8
         )
+
+
+@skip_on_valgrind
+def test_event_crosses_own_trigger_back(tempdir):
+    """The Heaviside variables of roots that an event assignment crosses back
+    in the same discontinuity are updated (gh-3425)."""
+    model_name = "test_event_crosses_own_trigger_back"
+    antimony2amici(
+        r"""
+        k1 = 0.7; k2 = 0.3; a = 1.5; thr = 0.5
+        x = 1; y = 0.5
+        x' = piecewise(-k1 * x, x > thr, -2 * k1 * x)
+        y' = piecewise(-k2 * y, x > thr, -k2 * x * y)
+        E1: at x < thr, fromTrigger=false: x = x * a
+        """,
+        model_name=model_name,
+        output_dir=tempdir,
+        observation_model=[
+            MC("obs_x", formula="x"),
+            MC("obs_y", formula="y"),
+        ],
+    )
+    model = import_model_module(model_name, tempdir).get_model()
+    k1, k2, a, thr = 0.7, 0.3, 1.5, 0.5
+    t = np.array([0.0, 1.0, 2.0, 3.0, 4.0])
+    model.set_timepoints(t)
+
+    # `x` reaches `thr` at `t0`, `E1` moves it back above `thr`, from where it
+    # decays to `thr` again within `period`. So, `x > thr` holds throughout,
+    # and `x = a^m * exp(-k1 * t)` after `m` event occurrences.
+    t0 = np.log(1 / thr) / k1
+    period = np.log(a) / k1
+    m = np.where(t < t0, 0, np.floor((t - t0) / period) + 1)
+    assert np.sum(m >= 3) >= 2
+    x = a**m * np.exp(-k1 * t)
+    y = 0.5 * np.exp(-k2 * t)
+    # The event times depend on `thr`, but the state jumps compensate the
+    # shifted right-hand side discontinuities, so `dx/dthr = 0`.
+    sx = {"k1": -t * x, "k2": 0 * t, "a": m * x / a, "thr": 0 * t}
+    sy = {"k1": 0 * t, "k2": -t * y, "a": 0 * t, "thr": 0 * t}
+    par_ids = list(model.get_free_parameter_ids())
+    ix = model.get_state_ids().index("x")
+    iy = model.get_state_ids().index("y")
+
+    edata = ExpData(model)
+    measurements = np.vstack([x, y]).T + 0.1
+    edata.set_measurements(measurements.flatten())
+    edata.set_noise_scales(np.ones_like(measurements).flatten())
+    for sens_method in (SensitivityMethod.forward, SensitivityMethod.adjoint):
+        solver = model.create_solver()
+        solver.set_sensitivity_order(SensitivityOrder.first)
+        solver.set_sensitivity_method(sens_method)
+        solver.set_relative_tolerance(1e-12)
+        solver.set_absolute_tolerance(1e-12)
+        rdata = run_simulation(model, solver, edata=edata)
+        assert rdata.status == AMICI_SUCCESS
+        assert_allclose(rdata.by_id("x"), x, rtol=1e-8)
+        assert_allclose(rdata.by_id("y"), y, rtol=1e-8)
+        # the right-hand side at the final time point uses the `x > thr`
+        # branch
+        assert_allclose(rdata.xdot[ix], -k1 * x[-1], rtol=1e-8)
+
+        expected_sllh = [
+            np.sum(0.1 * (sx[par_id] + sy[par_id])) for par_id in par_ids
+        ]
+        assert_allclose(
+            rdata.sllh,
+            expected_sllh,
+            rtol=1e-6,
+            atol=1e-8,
+            err_msg=str(sens_method),
+        )
+        if sens_method == SensitivityMethod.forward:
+            for ip, par_id in enumerate(par_ids):
+                assert_allclose(
+                    rdata.sx[:, ip, ix], sx[par_id], rtol=1e-6, atol=1e-8
+                )
+                assert_allclose(
+                    rdata.sx[:, ip, iy], sy[par_id], rtol=1e-6, atol=1e-8
+                )
+
+
+@skip_on_valgrind
+def test_event_moves_own_trigger_time(tempdir):
+    """An event that moves its own explicit trigger time fires again
+    (gh-3425)."""
+    model_name = "test_event_moves_own_trigger_time"
+    antimony2amici(
+        r"""
+        s = 1; n = 0
+        x = 1; x' = -x
+        E: at time >= s: s = s + 1, n = n + 1
+        """,
+        model_name=model_name,
+        output_dir=tempdir,
+    )
+    model = import_model_module(model_name, tempdir).get_model()
+    assert model.ne_solver == 0
+    model.set_timepoints([0.0, 1.5, 2.5, 3.5])
+    rdata = run_simulation(model, model.create_solver())
+    assert rdata.status == AMICI_SUCCESS
+    assert_allclose(rdata.by_id("n"), [0.0, 1.0, 2.0, 3.0], rtol=0, atol=0)
+
+
+@skip_on_valgrind
+@pytest.mark.parametrize("x_new", ["xmax", "0.3"])
+def test_event_assignment_onto_threshold(tempdir, x_new):
+    """An event assignment that sets a just-crossed root function to zero,
+    exactly or up to rounding errors, does not cross it back.
+
+    `xmax` is computed at runtime, `0.1 + 0.2 = 0.30000000000000004`, so
+    `x = 0.3` is one rounding error below it.
+    """
+    model_name = f"test_event_assignment_onto_threshold_{x_new != 'xmax'}"
+    antimony2amici(
+        rf"""
+        r = 1; a = 0.1; b = 0.2
+        xmax := a + b
+        x = 0; x' = piecewise(r, x < xmax, 0)
+        E: at x > xmax: x = {x_new}
+        """,
+        model_name=model_name,
+        output_dir=tempdir,
+    )
+    model = import_model_module(model_name, tempdir).get_model()
+    model.set_timepoints([0.0, 1.0, 2.0])
+    rdata = run_simulation(model, model.create_solver())
+    assert rdata.status == AMICI_SUCCESS
+    x_expected = 0.1 + 0.2 if x_new == "xmax" else 0.3
+    assert_allclose(rdata.by_id("x"), [0.0, x_expected, x_expected], rtol=1e-8)

src/forwardproblem.cpp

@@ -838,24 +838,40 @@ int EventHandlingSimulator::detect_secondary_events() {
 
     // check whether we need to fire a secondary event
     model_->froot(ws_->sol.t, ws_->sol.x, ws_->sol.dx, ws_->rootvals);
+    auto const& h = model_->get_model_state().h;
+
+    // estimated change of the root functions within the integration
+    // tolerances, only computed if needed
+    std::vector<realtype> root_noise;
+
     for (int ie = 0; ie < model_->ne; ie++) {
-        // the same event should not trigger itself
-        if (ws_->roots_found.at(ie) == 0) {
-            // check whether the value of the Heaviside function changed
-            if (heaviside_differs(ws_->rval_tmp.at(ie), ws_->rootvals.at(ie))) {
-                ws_->roots_found.at(ie)
-                    = ws_->rval_tmp.at(ie) < ws_->rootvals.at(ie) ? 1 : -1;
-                if (ws_->roots_found.at(ie) == 1)
-                    triggered_event = true;
-                secondevent++;
-            } else {
-                ws_->roots_found.at(ie) = 0;
-            }
-        } else {
-            // don't fire the same event again
-            ws_->roots_found.at(ie) = 0;
+        ws_->roots_found.at(ie) = 0;
+        auto const rootval = ws_->rootvals.at(ie);
+        // Root functions not changed by the preceding event assignment keep
+        // their Heaviside value.
+        if (rootval == ws_->rval_tmp.at(ie) || std::isnan(rootval))
+            continue;
+        // The roots crossed at this time are close to zero, so their sign
+        // may change due to roundoff or, for DAEs, the re-computation of the
+        // algebraic states. Only consider them crossed (back) if they are
+        // clearly on the other side.
+        if (ws_->roots_for_ignore.at(ie) != 0) {
+            if (root_noise.empty())
+                root_noise = compute_root_noise();
+            if (std::abs(rootval) <= root_noise[ie])
+                continue;
+        }
+        // Compare to the Heaviside variable, not to the previous root value,
+        // to also detect roots that are crossed back
+        realtype const h_new = rootval < 0.0 ? 0.0 : 1.0;
+        if (h_new != h.at(ie)) {
+            ws_->roots_found.at(ie) = h_new > h.at(ie) ? 1 : -1;
+            triggered_event |= ws_->roots_found.at(ie) == 1;
+            secondevent++;
         }
     }
+    // the next event assignment is compared to the current values
+    ws_->rval_tmp = ws_->rootvals;
 
     // fire the secondary event?
     if (secondevent > 0) {
@@ -883,6 +899,24 @@ int EventHandlingSimulator::detect_secondary_events() {
     return secondevent;
 }
 
+std::vector<realtype> EventHandlingSimulator::compute_root_noise() {
+    auto const rtol = solver_->get_relative_tolerance();
+    auto const atol = solver_->get_absolute_tolerance();
+    AmiVector x(ws_->sol.x);
+    AmiVector dx(ws_->sol.dx);
+    for (int ix = 0; ix < x.size(); ++ix)
+        x[ix] += atol + rtol * std::abs(x[ix]);
+    for (int ix = 0; ix < dx.size(); ++ix)
+        dx[ix] += atol + rtol * std::abs(dx[ix]);
+    std::vector<realtype> root_noise(model_->ne);
+    model_->froot(ws_->sol.t, x, dx, root_noise);
+    // restore the model's derived state at the current state
+    model_->froot(ws_->sol.t, ws_->sol.x, ws_->sol.dx, ws_->rootvals);
+    for (int ie = 0; ie < model_->ne; ++ie)
+        root_noise[ie] = std::abs(root_noise[ie] - ws_->rootvals.at(ie));
+    return root_noise;
+}
+
 void EventHandlingSimulator::collect_roots(realtype const t_trigger) {
     // roots of the root functions tracked by the solver
     solver_->get_root_info(ws_->roots_found.data());