fix: Add recovery for stiff ODE problems - #346
Closed
mhovd wants to merge 3 commits into
Closed
Conversation
Stiff models with sudden rate changes (infusion start/end, lagged boluses, covariate-driven clearance jumps) could kill a whole simulation through two diffsol traps: 1. Crushed-step doom loop: when approaching a stop, diffsol clamps h to the gap and silently ignores the StepSizeTooSmall this can raise. Once h sits below half the minimum step, every later step-size update -- including growth, capped at 2x -- fails, even though the state is healthy. Near-coincident stops (infusion end ULPs from an observation) trigger exactly this. 2. Stiff-restart rejection cascades: discontinuity restarts keep the old step size. A large h from a smooth segment meeting a stiff transient must reject its way down at <= 2x shrink per rejection within a budget of 40 consecutive rejections; very stiff transients exhaust the budget (TooManyErrorTestFailures / StepSizeTooSmall). The fix is layered, works for all solver families (BDF, SDIRK, RK), and uses no unsafe code: - Coalesce un-integrable gaps up front: stops closer than a few minimum timesteps are snapped to without integrating, after verifying against the scheduled infusion rates that the skip cannot omit material input (bounded by atol + rtol*||y||). The same check accepts StopTimeAtCurrentTime for a stop a few ULPs ahead and snaps the clock so the next segment cannot restart with the pre-boundary infusion rate. - Restart hygiene: every restart (boundary reinit or rescue) now floors the step size at max(256*eps*|t|, 1e-9), repairing any h crushed by diffsol's stop handling. The floor is a starting value only; the error controller may shrink below it afterwards. - Bounded rescue restarts: recoverable step failures (step-size collapse, error-test/Newton budget exhaustion, linear-solve failures) restart in place at first order with a fresh Jacobian and a small step, resetting diffsol's failure budgets so the solver can walk across the transient. Restart steps shrink 10x per stalled attempt; budgets of 16 rescues per segment and 8 consecutive stalls bound the work. Exhaustion returns the original descriptive error annotated with the restart count and a pointer toward implicit solvers. Each rescue emits a tracing debug event. This also crosses mid-segment RHS discontinuities that are not events, e.g. carry-forward covariate changes. - Deduplicate the event loop: the advance-to-next-event logic (stop selection, left-continuity at infusion boundaries, coalescing, rescues) now lives in one shared advance_solver_to_event(); the closure-based ODE path and the DSL runtime path both call it. The DSL path injects its model-function error surfacing through an after_step hook, which runs before rescue classification so user-code errors are reported as the root cause instead of the solver error they provoked. This removes the duplicated loop in dsl/native.rs that let the two paths diverge in the first place. tests/extreme_stiffness.rs exercises the failure modes across all solver families against closed-form solutions: infusion-end boundaries at ke up to 1e6 observed inside the transient, ULP-scale gaps between stops, stiff cutoff followed by a distant observation, a steep covariate clearance ramp mid-segment, a lagged bolus into stiff elimination, Michaelis-Menten depletion with tiny Km, and an impossibly stiff explicit-solver case asserting the descriptive error. The new tests fail without the fixes.
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #346 +/- ##
==========================================
+ Coverage 83.82% 84.02% +0.19%
==========================================
Files 82 82
Lines 33983 34113 +130
==========================================
+ Hits 28486 28662 +176
+ Misses 5497 5451 -46
🚀 New features to boost your workflow:
|
Contributor
There was a problem hiding this comment.
Pull request overview
This PR hardens the ODE “event loop” against extreme stiffness and floating-point edge cases (near-coincident stop times, discontinuity restarts), aiming to prevent otherwise-healthy simulations from aborting or hanging, while producing clearer failure messages when the chosen solver truly can’t handle the problem.
Changes:
- Adds bounded “rescue restart” logic (restart-in-place with fresh Jacobian / small step) and step-size floor handling on restarts to recover from diffsol step-size collapse and rejection cascades.
- Coalesces un-integrable tiny stop gaps (with an infusion-amount safety check) and centralizes the event-advance loop into a shared
advance_solver_to_event()used by both closure and DSL runtime paths. - Introduces new regression tests covering extreme stiffness scenarios and stop-time edge cases across solver families.
Reviewed changes
Copilot reviewed 6 out of 6 changed files in this pull request and generated 3 comments.
Show a summary per file
| File | Description |
|---|---|
| tests/extreme_stiffness.rs | Adds comprehensive regression tests for stiff discontinuities, near-ULP stop gaps, and “beyond rescue” descriptive failures across solvers. |
| tests/bolus_reinit_stop_time.rs | Adds regression tests for accepting StopTimeAtCurrentTime and ensuring correct post-boundary RHS semantics in both closure and DSL/JIT paths. |
| src/simulator/equation/ode/mod.rs | Implements shared event-loop advancement, stop-gap coalescing, restart step-size flooring, and bounded rescue restarts. |
| src/simulator/equation/ode/closure.rs | Adds infusion-amount computation used to validate safe coalescing of tiny stop gaps. |
| src/error/mod.rs | Adds rescue-attempt context to solver error messages. |
| src/dsl/native.rs | Removes duplicated DSL-native event loop in favor of the shared advance_solver_to_event() path; adjusts function invocation. |
💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.
Comment on lines
27
to
30
| use diffsol::{ | ||
| error::OdeSolverError, ode_solver::method::OdeSolverMethod, NalgebraContext, OdeBuilder, | ||
| OdeSolverStopReason, Vector, VectorHost, | ||
| OdeSolverConfig, OdeSolverStopReason, Vector, VectorHost, | ||
| }; |
Comment on lines
+896
to
+902
| Ok(OdeSolverStopReason::RootFound(_, _)) => { | ||
| return Err(PharmsolError::OtherError(format!( | ||
| "solver stopped at an unexpected root at t = {:.4} \ | ||
| (root finding is not configured)", | ||
| stop_time | ||
| ))); | ||
| } |
Comment on lines
+71
to
+75
| PharmsolError::DiffsolError(msg) => PharmsolError::DiffsolError(format!( | ||
| "{msg} ({rescues} automatic solver restart(s) did not recover; \ | ||
| the system may be too stiff for the selected solver — stiff \ | ||
| problems need an implicit solver such as BDF)" | ||
| )), |
Contributor
|
Bencher Report| Project | pharmsol |
| Branch | fix/ode-stiff-recovery |
| Testbed | mhovd-pgx |
⚠️ WARNING: Truncated view!The full continuous benchmarking report exceeds the maximum length allowed on this platform.
🚨 2 Alerts
🐰 View full continuous benchmarking report in Bencher
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Note: this PR was drafted by Fable 5, with manual review.
The LLM notes are below
Stiff models with sudden rate changes (infusion start/end, lagged
boluses, covariate-driven clearance jumps) could kill a whole
simulation through two diffsol traps:
Crushed-step doom loop: when approaching a stop, diffsol clamps
h to the gap and silently ignores the StepSizeTooSmall this can
raise. Once h sits below half the minimum step, every later
step-size update -- including growth, capped at 2x -- fails, even
though the state is healthy. Near-coincident stops (infusion end
ULPs from an observation) trigger exactly this.
Stiff-restart rejection cascades: discontinuity restarts keep the
old step size. A large h from a smooth segment meeting a stiff
transient must reject its way down at <= 2x shrink per rejection
within a budget of 40 consecutive rejections; very stiff transients
exhaust the budget (TooManyErrorTestFailures / StepSizeTooSmall).
The fix is layered, works for all solver families (BDF, SDIRK, RK),
and uses no unsafe code:
Coalesce un-integrable gaps up front: stops closer than a few
minimum timesteps are snapped to without integrating, after
verifying against the scheduled infusion rates that the skip
cannot omit material input (bounded by atol + rtol*||y||). The
same check accepts StopTimeAtCurrentTime for a stop a few ULPs
ahead and snaps the clock so the next segment cannot restart
with the pre-boundary infusion rate.
Restart hygiene: every restart (boundary reinit or rescue) now
floors the step size at max(256eps|t|, 1e-9), repairing any h
crushed by diffsol's stop handling. The floor is a starting value
only; the error controller may shrink below it afterwards.
Bounded rescue restarts: recoverable step failures (step-size
collapse, error-test/Newton budget exhaustion, linear-solve
failures) restart in place at first order with a fresh Jacobian
and a small step, resetting diffsol's failure budgets so the
solver can walk across the transient. Restart steps shrink 10x
per stalled attempt; budgets of 16 rescues per segment and 8
consecutive stalls bound the work. Exhaustion returns the
original descriptive error annotated with the restart count and
a pointer toward implicit solvers. Each rescue emits a tracing
debug event. This also crosses mid-segment RHS discontinuities
that are not events, e.g. carry-forward covariate changes.
Deduplicate the event loop: the advance-to-next-event logic
(stop selection, left-continuity at infusion boundaries,
coalescing, rescues) now lives in one shared
advance_solver_to_event(); the closure-based ODE path and the
DSL runtime path both call it. The DSL path injects its
model-function error surfacing through an after_step hook, which
runs before rescue classification so user-code errors are
reported as the root cause instead of the solver error they
provoked. This removes the duplicated loop in dsl/native.rs that
let the two paths diverge in the first place.
tests/extreme_stiffness.rs exercises the failure modes across all
solver families against closed-form solutions: infusion-end
boundaries at ke up to 1e6 observed inside the transient, ULP-scale
gaps between stops, stiff cutoff followed by a distant observation,
a steep covariate clearance ramp mid-segment, a lagged bolus into
stiff elimination, Michaelis-Menten depletion with tiny Km, and an
impossibly stiff explicit-solver case asserting the descriptive
error. The new tests fail without the fixes.