Fix check for sign change, fix Heaviside update

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 {

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]

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