Skip to content

fix: Add recovery for stiff ODE problems #346

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
mhovd wants to merge 3 commits into from
Closed

fix: Add recovery for stiff ODE problems #346

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
97f52ad
fix: handle close ODE stop times safely
Siel Aug 20, 2026
733f986
fix: recover ODE solves from stiff transients and crushed step sizes
mhovd Aug 21, 2026
6fa3b8e
Address copilot review items
mhovd Aug 21, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
131 changes: 18 additions & 113 deletions src/dsl/native.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,8 +3,7 @@ use std::collections::HashMap;
use std::sync::Arc;

use diffsol::{
error::OdeSolverError, ode_solver::method::OdeSolverMethod, NalgebraContext, OdeBuilder,
OdeSolverStopReason, Vector, VectorHost,
ode_solver::method::OdeSolverMethod, NalgebraContext, OdeBuilder, Vector, VectorHost,
};
use nalgebra::DVector;
use ndarray::{concatenate, Array2, Axis};
Expand Down Expand Up @@ -234,7 +233,7 @@ impl FunctionSession for NativeFunctionSession<'_> {
))
})?;

function(time, states, params, covariates, routes, derived, out);
unsafe { function(time, states, params, covariates, routes, derived, out) };
Ok(())
}
}
Expand Down Expand Up @@ -1349,13 +1348,9 @@ impl NativeOdeModel {
F: Fn(&V, &V, f64, &mut V, &V, &V, &Covariates) + 'a,
S: OdeSolverMethod<'a, PMProblem<'a, F>>,
{
// Mirror the closure-based ODE event loop: stop at every infusion
// start and end boundary in addition to subject events, using the
// left-continuous rate while integrating toward a boundary and the
// right-continuous rate after reaching it. This keeps the JIT
// implementation numerically consistent with the reference [`ODE`]
// path (see `ode::run_events`).
let infusion_boundary_times = solver.problem().eqn.infusion_boundary_times();
// The event-to-event integration loop is shared with the closure-based
// [`ODE`] equation so compiled models stay numerically consistent with
// the reference path by construction.
let mut infusion_boundary_cursor = 0usize;
let mut index = 0usize;
// Set when the previous event changed the state or the previous stop
Expand Down Expand Up @@ -1402,109 +1397,19 @@ impl NativeOdeModel {

// Advance to the next event time if it exists.
if let Some(next_event) = next_event {
let next_event_time = next_event.time();
while next_event_time > solver.state().t {
while infusion_boundary_cursor < infusion_boundary_times.len()
&& infusion_boundary_times[infusion_boundary_cursor] <= solver.state().t
{
infusion_boundary_cursor += 1;
}

let (stop_time, is_infusion_boundary) = if let Some(stop_time) =
infusion_boundary_times.get(infusion_boundary_cursor)
{
if *stop_time <= next_event_time {
infusion_boundary_cursor += 1;
(*stop_time, true)
} else {
(next_event_time, false)
}
} else {
(next_event_time, false)
};

solver
.problem()
.eqn
.set_left_continuity_time(if is_infusion_boundary {
Some(stop_time)
} else {
None
});

match solver.set_stop_time(stop_time) {
Ok(_) => {
if pending_reinit {
crate::simulator::equation::ode::reinitialize_at_boundary(
solver, dy_scratch,
);
pending_reinit = false;
}
loop {
match solver.step() {
Ok(_) if function_error.borrow().is_some() => {
return Err(function_error.borrow_mut().take().unwrap());
}
Ok(OdeSolverStopReason::InternalTimestep) => continue,
Ok(OdeSolverStopReason::TstopReached) => {
solver.problem().eqn.set_left_continuity_time(None);
if is_infusion_boundary {
pending_reinit = true;
}
break;
}
Ok(OdeSolverStopReason::RootFound(_, _)) => {
return Err(PharmsolError::OtherError(format!(
"solver stopped at an unexpected root at t = {:.4} \
(root finding is not configured)",
stop_time
)));
}
Err(err) => {
return Err(PharmsolError::from_solver_error(
err, stop_time,
));
}
}
}
}
Err(diffsol::error::DiffsolError::OdeSolverError(
OdeSolverError::StopTimeAtCurrentTime,
)) => {
solver.problem().eqn.set_left_continuity_time(None);
let state_t = solver.state().t;
let stop_reached = crate::simulator::equation::ode::stop_time_reached(
stop_time, state_t,
);

if stop_reached {
if is_infusion_boundary {
pending_reinit = true;
}
// The requested stop is the current time within
// a small relative tolerance. If it is an
// infusion boundary before the next subject
// event, keep integrating toward the event;
// break only when the reached stop is the
// event time itself.
if stop_time < next_event_time {
continue;
}
break;
}
return Err(PharmsolError::from_solver_error(
diffsol::error::DiffsolError::OdeSolverError(
OdeSolverError::StopTimeAtCurrentTime,
),
stop_time,
));
}
Err(err) => {
solver.problem().eqn.set_left_continuity_time(None);
return Err(PharmsolError::from_solver_error(err, stop_time));
}
}
}
crate::simulator::equation::ode::advance_solver_to_event(
solver,
next_event.time(),
&mut infusion_boundary_cursor,
&mut pending_reinit,
dy_scratch,
self.rtol,
self.atol,
&mut || match function_error.borrow_mut().take() {
Some(error) => Err(error),
None => Ok(()),
},
)?;
}
index += 1;
}
Expand Down
17 changes: 17 additions & 0 deletions src/error/mod.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,23 @@ impl PharmsolError {
pub fn from_solver_error(error: diffsol::error::DiffsolError, target_time: f64) -> Self {
PharmsolError::DiffsolError(describe_diffsol_error(&error, Some(target_time)))
}

/// Annotate a solver error with the number of automatic in-place restarts
/// that were attempted before giving up. No-op for zero attempts.
pub(crate) fn with_rescue_context(self, rescues: usize) -> Self {
if rescues == 0 {
return self;
}
match self {
PharmsolError::DiffsolError(msg) => PharmsolError::DiffsolError(format!(
"{msg} ({rescues} automatic solver restart(s) did not recover; \
the problem is too stiff for the selected solver at this point — \
if using an explicit solver, switch to an implicit one such as BDF; \
otherwise check for extreme or implausible parameter values)"
)),
other => other,
}
}
}

impl PharmsolError {
Expand Down
18 changes: 18 additions & 0 deletions src/simulator/equation/ode/closure.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -182,6 +182,20 @@ impl InfusionSchedule {
&self.boundary_times
}

/// Absolute infusion input omitted when two adjacent stops are coalesced.
/// The event loop calls this only before the next infusion boundary, so the
/// active rate is constant over the interval.
fn infusion_amount_between(&self, from: f64, to: f64) -> f64 {
if to <= from {
return 0.0;
}
let duration = to - from;
self.tracks
.iter()
.map(|track| track.rate_at_left(to).abs() * duration)
.sum()
}

fn fill_rate_vector(&self, time: f64, rateiv: &mut V) {
let left_continuity_time = self.left_continuity_time.get();
rateiv.fill(0.0);
Expand Down Expand Up @@ -414,6 +428,10 @@ where
self.infusion_schedule.infusion_boundary_times()
}

pub(crate) fn infusion_amount_between(&self, from: f64, to: f64) -> f64 {
self.infusion_schedule.infusion_amount_between(from, to)
}

/// Evaluate the full RHS (including the currently scheduled infusion
/// rates) at time `t` into `dx`.
///
Expand Down
Loading
Loading