From e64ef5d4cd4e26f13a2fb02fddcb3bdc1fbbdc51 Mon Sep 17 00:00:00 2001
From: Daniel Weindl <daniel.weindl@uni-bonn.de>
Date: Mon, 28 Jul 2025 14:32:18 +0200
Subject: [PATCH] Fix check for sign change, fix Heaviside update

The old check for sign changes did not work if one values was zero.

Heaviside state has to be updated independent of whether the event occurs at the initial timepoint or not.

Fixes #2926.
---
 include/amici/misc.h        | 14 ++++++++++++++
 python/tests/test_events.py | 33 +++++++++++++++++++++++++++++++++
 src/forwardproblem.cpp      |  8 +++-----
 3 files changed, 50 insertions(+), 5 deletions(-)

diff --git a/include/amici/misc.h b/include/amici/misc.h
index a9821c1cfc..3f9caf8b0c 100644
--- a/include/amici/misc.h
+++ b/include/amici/misc.h
@@ -252,6 +252,20 @@ 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});
+}
+
 #ifdef BOOST_CHRONO_HAS_THREAD_CLOCK
 /** Tracks elapsed CPU time using boost::chrono::thread_clock. */
 class CpuTimer {
diff --git a/python/tests/test_events.py b/python/tests/test_events.py
index efc22f81c1..8a16b7cba4 100644
--- a/python/tests/test_events.py
+++ b/python/tests/test_events.py
@@ -1260,3 +1260,36 @@ def test_preeq_presim_preserve_heaviside_state(tempdir):
     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]
+
+
+@skip_on_valgrind
+def test_gh2926(tempdir):
+    """Two simultaneous events. Event `E1` changes the root function
+    for the piecewise-switch from 0 to <0."""
+    from amici.antimony_import import antimony2amici
+
+    model_name = "test_gh2926"
+    antimony2amici(
+        r"""
+        some_time = 0
+        some_time' = 1
+
+        threshold = 0
+
+        x1 := piecewise(2, (some_time >= threshold), 1)
+        E1: at some_time >= 0, t0 = false: threshold = 10
+        """,
+        model_name=model_name,
+        output_dir=tempdir,
+    )
+
+    model_module = import_model_module(model_name, tempdir)
+    model = model_module.get_model()
+    # check output just after t=10, otherwise the solver stops at `10 - eps`
+    #  and the event triggers only after the output was recorded
+    model.setTimepoints([0, 5, 10.1, 11])
+    solver = model.getSolver()
+
+    rdata = amici.runAmiciSimulation(model, solver)
+    assert rdata.status == amici.AMICI_SUCCESS
+    assert rdata.by_id("x1").tolist() == [1.0, 1.0, 2.0, 2.0]
diff --git a/src/forwardproblem.cpp b/src/forwardproblem.cpp
index e2cd16e85e..934117c349 100644
--- a/src/forwardproblem.cpp
+++ b/src/forwardproblem.cpp
@@ -495,9 +495,7 @@ void EventHandlingSimulator::handle_events(
 
             store_pre_event_info(true);
 
-            if (!initial_event) {
-                model_->updateHeaviside(ws_->roots_found);
-            }
+            model_->updateHeaviside(ws_->roots_found);
         }
     }
     store_post_event_info();
@@ -616,8 +614,8 @@ int EventHandlingSimulator::detect_secondary_events() {
     for (int ie = 0; ie < model_->ne; ie++) {
         // the same event should not trigger itself
         if (ws_->roots_found.at(ie) == 0) {
-            // check whether there was a zero-crossing
-            if (0 > ws_->rval_tmp.at(ie) * ws_->rootvals.at(ie)) {
+            // check whether the value of the Heaviside function changed
+            if (heaviside_differs(ws_->rval_tmp.at(ie), ws_->rootvals.at(ie))) {
                 if (ws_->rval_tmp.at(ie) < ws_->rootvals.at(ie)) {
                     ws_->roots_found.at(ie) = 1;
                     auto const& event = model_->get_event(ie);