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