[dweindl]

Fix event handling for event assignments that move a root function back across or onto its threshold, using post-hoc root filtering: the solver reports every root crossing, and those that don't change the Heaviside variables are discarded. Alternative to #3490 (root masking: the solver only reports crossings in the direction that would change the Heaviside variables), 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.

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.

The solver reports all root crossings. Those that don't change h are dropped (debug message ROOT_IGNORED) and the integration continues without handling a discontinuity. This covers a just-crossed root that is re-detected after reinitialization, and a root function whose h was set by its direction and that crosses zero a few steps later, which previously changed h to -1 or 2 and, e.g., triggered events twice. A dropped root that is exactly zero is deactivated until it becomes nonzero, as SUNDIALS does at a restart, since otherwise CVODES/IDAS fail with "Root found at and very near t" on the next step (SBML cases 00354, 00754–00756). This requires access to SUNDIALS internals (cv_gactive/ida_gactive).

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). That filter can't be made correct: from the root index, its direction, and the fact that it was found in the first step after reinitialization, a genuine fast re-crossing (e.g., the next bounce of a ball within that step) can't be distinguished from a re-detection, so it either misses events or leaves h inconsistent. With roots that don't change h dropped, it's redundant: all tests and the SBML cases below pass without it.

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, identical to #3490):

• 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 root masking (#3490):

• The forward-problem code is the same except for comments, and both pass the same tests.
• In Fix events crossed back or moved onto thresholds by event assignments, using root masking #3490, the solver doesn't stop at crossings that don't change h (CVodeSetRootDirection/IDASetRootDirection), so no SUNDIALS internals are needed for exact zeros, but h has to be up to date whenever the solver is called. Here, the solver stops at these crossings; this resumes the integration without reinitialization and didn't measurably affect performance.
• C++ logic lines (without comments and blank lines) vs. main: +210/−144 here, +232/−144 in Fix events crossed back or moved onto thresholds by event assignments, using root masking #3490.
• 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 99.43503% with 1 line in your changes missing coverage. Please review.
✅ Project coverage is 80.32%. Comparing base (b15fed8) to head (ece48c5).

Files with missing lines	Patch %	Lines
src/solver_idas.cpp	80.00%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #3491      +/-   ##
==========================================
+ Coverage   80.23%   80.32%   +0.08%     
==========================================
  Files         331      331              
  Lines       23420    23481      +61     
  Branches     1631     1645      +14     
==========================================
+ Hits        18792    18860      +68     
+ Misses       4622     4615       -7     
  Partials        6        6              

Flag	Coverage Δ
cpp	74.92% <99.43%> (+0.09%)	⬆️
cpp_python	34.89% <62.14%> (+0.23%)	⬆️
petab	46.87% <61.58%> (+0.22%)	⬆️
petab_sciml	16.55% <0.00%> (-0.06%)	⬇️
petab_sciml_benchmarks	15.15% <0.00%> (-0.05%)	⬇️
python	72.65% <90.39%> (+0.09%)	⬆️

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/forwardproblem.cpp	93.84% <100.00%> (+0.88%)	⬆️
src/solver_cvodes.cpp	75.00% <100.00%> (-0.96%)	⬇️
src/solver_idas.cpp	70.96% <80.00%> (-0.04%)	⬇️

🚀 New features to boost your workflow:

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