[dweindl]

Fix event handling for event assignments that move a root function back across or onto its threshold, using root masking as suggested in #3439: the solver only reports root crossings in the direction that would change the Heaviside variables. Alternative to #3491 (post-hoc root filtering: the solver reports every crossing and the irrelevant ones are discarded), which has the same functionality and passes the same tests; one of them supersedes #3439. Fixes #3425.

If an event assignment moved a just-crossed root function back across zero, e.g., E1: at x < thr: x = 1.5 * x, the Heaviside variable kept the post-crossing value, and later crossings incremented it further (h = 2, 3, …), so piecewise expressions took the wrong branch. Likewise, an event moving its own explicit trigger time, E: at time >= s: s = s + 1, fired only once, and an assignment exactly onto a threshold, e.g., a bouncing ball (E: at x <= 0: x = 0, v = -e * v) or onto another event's threshold, left stale Heaviside variables, as did an event that changes the direction of a root function sitting at its threshold.

The solver now only reports root crossings that change the Heaviside variables: before each Solver::run/Solver::step, CVodeSetRootDirection/IDASetRootDirection are set from h (h = 1: only decreasing, h = 0: only increasing). This replaces the ignore-after-reinit machinery from #3266 (Solver::ignore_roots_after_reinit, the root filtering in reinit_post_process_f, the ROOT_AFTER_REINIT warning), since a just-crossed root can no longer be re-detected. SUNDIALS still reports a root that becomes exactly zero right after a step regardless of its direction; such roots are dropped (debug message ROOT_IGNORED) and the integration continues.

With masking, h must be consistent with the root functions, otherwise crossings are missed. After each event assignment, the root functions are compared to h. Root functions within their uncertainty of zero, estimated from a state perturbation within the integration tolerances and from the time tolerance of root-finding, are resolved as follows: an assignment that moves a root function onto zero triggers the event (the trigger is true at zero); afterwards, the direction in which the root function moves decides (central finite difference in time). The same holds for root functions that crossed zero at the current time, for unchanged root functions at zero whose direction may have changed, and, at the start of a simulation or steady-state computation, for root functions at zero. A switch that flips back and forth at the same time, e.g., x' = piecewise(-1, x > 0, 1) at x = 0, now fails with an error instead of looping or giving wrong results.

With adjoint sensitivities, the event-time terms are now passed through secondary discontinuities without event assignments, and forward and adjoint sensitivities use the last rising root, or else the last crossing root, for the event time, so that secondary events triggered at a discontinuity with only falling roots get event-time sensitivities (not covered by a test).

Tests (test_events.py):

• from Fix Heaviside variables of roots crossed back by event assignments #3439: test_event_crosses_own_trigger_back, test_event_moves_own_trigger_time, test_event_assignment_onto_threshold
• test_bouncing_ball: event counts, trajectory, and forward/adjoint sensitivities against the analytical solution
• test_sawtooth_many_events: 1000 periods of events that reverse their own root function; event counts and forward/adjoint sensitivities against the analytical solution (fails on main with adjoint sensitivities)
• test_root_at_threshold_changes_direction (fails on main), test_assignment_onto_other_threshold, test_trigger_leaves_threshold (after an assignment / at t0), test_event_moves_own_trigger_back_within_tolerance, test_chattering_switch_fails
• test_root_after_reinit_direction_is_checked is renamed to test_event_assignment_just_below_other_trigger, and the ROOT_AFTER_REINIT checks are removed.

The SBML semantic test suite cases with events or piecewise give the same results as on main (627 passed, including forward and adjoint sensitivity checks).

Comparison with post-hoc root filtering (#3491):

• The forward-problem code is the same except for comments, and both pass the same tests. Both remove the ignore-after-reinit filter, which can't distinguish a fast genuine re-crossing in the first step after reinitialization from a re-detection, and is redundant once roots that don't change h are dropped.
• Here, the solver doesn't stop at crossings that don't change h. In Fix events crossed back or moved onto thresholds by event assignments, using post-hoc root filtering #3491, it does, and a dropped root that is exactly zero has to be deactivated via SUNDIALS internals (cv_gactive/ida_gactive), as the solver doesn't continue otherwise ("Root found at and very near t", SBML cases 00354, 00754–00756). This requires h to be up to date whenever the solver is called, which Solver::apply_root_directions documents.
• C++ logic lines (without comments and blank lines) vs. main: +232/−144 here, +210/−144 in Fix events crossed back or moved onto thresholds by event assignments, using post-hoc root filtering #3491.
• Performance: no difference beyond noise and identical step counts on the benchmark collection models with events (Smith_BMCSystBiol2013, Liu_IFACPapersOnLine2025, Alkan_SciSignal2018, Isensee_JCB2018, Giordano_Nature2020, Brannmark_JBC2010) and on toy models with many events. For Smith_BMCSystBiol2013, the likelihood with forward sensitivities on main differed from the one without sensitivities (−896682 vs. −888198); both branches give −888198 in all modes.

Known limitations (both): The direction of a root function within its uncertainty of zero is estimated by finite differences, which is unreliable for very small slopes, and for DAEs uses the algebraic states' derivatives only approximately. An event that moves its own root function back by less than that uncertainty, e.g., at x >= 1e6: x = x - 0.005 with default tolerances, does not trigger again, and a root function that dips below zero by less than it within a single step is not detected (as on main). Event assignments are applied as state changes, so x = 1 may result in x one rounding error off (as on main).

🤖 Generated with Claude Code

[codecov]

Codecov Report

❌ Patch coverage is 94.24084% with 11 lines in your changes missing coverage. Please review.
✅ Project coverage is 80.29%. Comparing base (b15fed8) to head (c4c82af).

Files with missing lines	Patch %	Lines
src/forwardproblem.cpp	95.23%	7 Missing ⚠️
src/solver_cvodes.cpp	71.42%	2 Missing ⚠️
src/solver_idas.cpp	71.42%	2 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #3490      +/-   ##
==========================================
+ Coverage   80.23%   80.29%   +0.05%     
==========================================
  Files         331      331              
  Lines       23420    23495      +75     
  Branches     1631     1647      +16     
==========================================
+ Hits        18792    18865      +73     
- Misses       4622     4624       +2     
  Partials        6        6              

Flag	Coverage Δ
cpp	74.90% <94.24%> (+0.07%)	⬆️
cpp_python	34.90% <61.25%> (+0.24%)	⬆️
petab	46.90% <63.87%> (+0.26%)	⬆️
petab_sciml	16.54% <0.00%> (-0.07%)	⬇️
petab_sciml_benchmarks	15.13% <0.00%> (-0.07%)	⬇️
python	72.60% <83.76%> (+0.03%)	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

Files with missing lines	Coverage Δ
include/amici/forwardproblem.h	100.00% <100.00%> (ø)
include/amici/misc.h	91.83% <100.00%> (+1.83%)	⬆️
include/amici/solver.h	100.00% <ø> (ø)
include/amici/solver_cvodes.h	100.00% <ø> (ø)
include/amici/solver_idas.h	100.00% <ø> (ø)
src/backwardproblem.cpp	94.72% <100.00%> (-0.04%)	⬇️
src/solver.cpp	79.88% <100.00%> (+0.26%)	⬆️
src/solver_cvodes.cpp	74.28% <71.42%> (-1.68%)	⬇️
src/solver_idas.cpp	70.61% <71.42%> (-0.40%)	⬇️
src/forwardproblem.cpp	93.00% <95.23%> (+0.04%)	⬆️

... and 2 files with indirect coverage changes

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.