diff --git a/src/dsl/native.rs b/src/dsl/native.rs index f2fc1955..3ea2b8ce 100644 --- a/src/dsl/native.rs +++ b/src/dsl/native.rs @@ -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}; @@ -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(()) } } @@ -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 @@ -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; } diff --git a/src/error/mod.rs b/src/error/mod.rs index 80d1201e..f4de2e46 100644 --- a/src/error/mod.rs +++ b/src/error/mod.rs @@ -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 { diff --git a/src/simulator/equation/ode/closure.rs b/src/simulator/equation/ode/closure.rs index 14ae22e5..bbf4f888 100644 --- a/src/simulator/equation/ode/closure.rs +++ b/src/simulator/equation/ode/closure.rs @@ -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); @@ -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`. /// diff --git a/src/simulator/equation/ode/mod.rs b/src/simulator/equation/ode/mod.rs index c90128ab..29c67aa0 100644 --- a/src/simulator/equation/ode/mod.rs +++ b/src/simulator/equation/ode/mod.rs @@ -26,7 +26,7 @@ use crate::simulator::equation::Predictions; use closure::PMProblem; use diffsol::{ error::OdeSolverError, ode_solver::method::OdeSolverMethod, NalgebraContext, OdeBuilder, - OdeSolverStopReason, Vector, VectorHost, + OdeSolverConfig, OdeSolverStopReason, Vector, VectorHost, }; use nalgebra::DVector; use pharmsol_dsl::ModelKind; @@ -313,7 +313,7 @@ fn _simulate_subject_dense( Some(error_models) => Some(ode.bind_error_models(error_models)?), None => None, }; - let bound_error_models = bound_error_models.as_ref().map(|models| &**models); + let bound_error_models = bound_error_models.as_deref(); let mut output = SubjectPredictions::new(ode.nparticles()); @@ -560,7 +560,9 @@ impl EquationPriv for ODE { /// (right-continuous) RHS so a first-order restart predicts with the new /// dynamics instead of the pre-boundary ones; /// - `state_mut` marks the state as modified so the next step restarts the -/// multi-step method at first order. +/// multi-step method at first order; +/// - the step size is raised to at least [`restart_step_size_floor`] so a +/// step crushed by a nearby previous stop cannot doom the next segment. /// /// Shared with the DSL/JIT ODE path ([`crate::dsl::native::NativeOdeModel`]), /// which must apply the same discontinuity semantics as the closure-based @@ -569,6 +571,26 @@ pub(crate) fn reinitialize_at_boundary<'a, F, S>(solver: &mut S, dy_scratch: &mu where F: Fn(&V, &V, f64, &mut V, &V, &V, &Covariates) + 'a, S: OdeSolverMethod<'a, PMProblem<'a, F>>, +{ + restart_solver_at_current_state(solver, dy_scratch, None); +} + +/// Shared restart core for boundary reinitializations and rescue restarts. +/// +/// `step_size` replaces the current step size when given (rescue restarts); +/// otherwise the current step size is kept. Either way the result is floored: +/// diffsol shrinks `h` to land exactly on a close stop and *ignores* the +/// `StepSizeTooSmall` this may raise, and once `h` sits below half the +/// solver's minimum step every later step-size update — including growth, +/// which is capped at 2x — fails. Restarting is the safe place to undo that: +/// the state is already marked modified, so raising `h` costs nothing extra. +fn restart_solver_at_current_state<'a, F, S>( + solver: &mut S, + dy_scratch: &mut V, + step_size: Option, +) where + F: Fn(&V, &V, f64, &mut V, &V, &V, &Covariates) + 'a, + S: OdeSolverMethod<'a, PMProblem<'a, F>>, { let state = solver.state_clone(); solver.set_state(state); @@ -581,26 +603,355 @@ where .eqn .refresh_state_derivative(t, y, dy_scratch); } + let floor = restart_step_size_floor(t); let state = solver.state_mut(); state.dy.copy_from(dy_scratch); + *state.h = step_size.unwrap_or(*state.h).max(floor); +} + +/// Fresh-start step size for rescue restarts; matches the `h0` given to +/// `OdeBuilder`. +const RESTART_STEP_SIZE: f64 = 1e-3; +/// Cap on rescue restarts within one integration segment. +const MAX_RESCUES_PER_SEGMENT: usize = 16; +/// Cap on consecutive rescue restarts that fail to advance the state time. +const MAX_STALLED_RESCUES: usize = 8; + +/// Smallest step size a restarted segment may begin with. +/// +/// This is a starting value, not a hard limit: the error controller is free to +/// shrink below it afterwards (down to diffsol's minimum step). It only needs +/// to sit far enough above the minimum step (1e-13) that step-size growth, +/// capped at 2x per update, can never fail, and far enough above the ULP +/// spacing of `t` that steps still advance time. +fn restart_step_size_floor(t: f64) -> f64 { + (f64::EPSILON * t.abs() * 256.0).max(1e-9) +} + +/// Advance below this threshold counts a rescue restart as stalled. +fn progress_epsilon(t: f64) -> f64 { + f64::EPSILON * t.abs().max(1.0) * 100.0 +} + +/// Step failures that an in-place restart (fresh Jacobian, first order, small +/// step) can plausibly clear. +/// +/// These are the failure modes of a stiff transient or a mid-segment RHS +/// discontinuity (e.g. a covariate change): the failure budgets are exhausted +/// or the step collapses while the state itself is still the last accepted, +/// finite solution. Setup, interpolation, and ordering errors are not +/// recoverable by restarting. +fn recoverable_step_failure(error: &diffsol::error::DiffsolError) -> bool { + use diffsol::error::DiffsolError; + match error { + DiffsolError::OdeSolverError(ode_error) => matches!( + ode_error, + OdeSolverError::StepSizeTooSmall { .. } + | OdeSolverError::TooManyErrorTestFailures { .. } + | OdeSolverError::TooManyNonlinearSolverFailures { .. } + ), + DiffsolError::NonLinearSolverError(_) | DiffsolError::LaError(_) => true, + _ => false, + } } -/// Whether a requested solver stop is effectively at the current state time. +/// Restart step size for the next rescue attempt. /// -/// diffsol reports `StopTimeAtCurrentTime` not only for a stop exactly at the -/// current time, but also when its internal state time has landed a few ULPs -/// past the requested stop (adaptive steps may end slightly beyond a stop). -/// Dense output grids built with floating-point arithmetic routinely place -/// requested times a few ULPs away from event times (e.g. a `t += dt` -/// accumulation puts a point ~16 ULPs after a bolus at `t = 12`), so accept a -/// stop within a small relative tolerance of the current time instead of -/// erroring. The tolerance stays far below any meaningful time difference: -/// ~64-128 ULPs of the current time, i.e. at most ~1e-13 at `t = 12`. +/// Starts at the fresh-start size (bounded by the remaining segment) and +/// shrinks tenfold per stalled attempt, never below the restart floor. +fn rescue_step_size(state_time: f64, stop_time: f64, stalled_rescues: usize) -> f64 { + let floor = restart_step_size_floor(state_time); + let remaining = (stop_time - state_time).max(0.0); + let base = RESTART_STEP_SIZE.min(remaining * 0.5).max(floor); + (base * 10f64.powi(-(stalled_rescues.min(12) as i32))).max(floor) +} + +fn solver_state_is_finite<'a, F, S>(solver: &S) -> bool +where + F: Fn(&V, &V, f64, &mut V, &V, &V, &Covariates) + 'a, + S: OdeSolverMethod<'a, PMProblem<'a, F>>, +{ + let state = solver.state(); + state.t.is_finite() && state.y.as_slice().iter().all(|value| value.is_finite()) +} + +/// Gap below which the event loop coalesces a forward stop itself instead of +/// asking the solver to integrate it. +/// +/// diffsol cannot integrate a segment shorter than its minimum step: its +/// pre-step stop handling clamps `h` to the gap and either fails outright +/// (the Runge-Kutta family re-clamps on every retry, so restarts cannot +/// recover) or leaves behind a crushed step size that fails on the next +/// update. A gap of a few minimum steps is indistinguishable from the same +/// time for simulation purposes. +fn min_integrable_gap<'a, F, S>(solver: &S, state_time: f64) -> f64 +where + F: Fn(&V, &V, f64, &mut V, &V, &V, &Covariates) + 'a, + S: OdeSolverMethod<'a, PMProblem<'a, F>>, +{ + let minimum_timestep = *solver.config().as_base_ref().minimum_timestep; + (minimum_timestep * 4.0).max(f64::EPSILON * state_time.abs() * 4.0) +} + +/// Snap the solver's logical time forward to `stop_time` without integrating, +/// after verifying that doing so cannot omit material infusion input. /// -/// Shared with the DSL/JIT ODE path ([`crate::dsl::native::NativeOdeModel`]). -pub(crate) fn stop_time_reached(stop_time: f64, state_t: f64) -> bool { - let tolerance = f64::EPSILON * state_t.abs().max(1.0) * 64.0; - (stop_time - state_t).abs() <= tolerance +/// The state itself is left untouched; only stops closer than the solver can +/// resolve are coalesced, so the skipped interval is below the solution +/// tolerances by construction (and this is enforced against the scheduled +/// infusion rates). The caller must restart the solver before the next +/// segment because `state_mut` marks the state as modified. +fn coalesce_stop_to<'a, F, S>( + solver: &mut S, + stop_time: f64, + rtol: f64, + atol: f64, +) -> Result<(), PharmsolError> +where + F: Fn(&V, &V, f64, &mut V, &V, &V, &Covariates) + 'a, + S: OdeSolverMethod<'a, PMProblem<'a, F>>, +{ + let state_time = solver.state().t; + let skipped_infusion = solver + .problem() + .eqn + .infusion_amount_between(state_time, stop_time); + let state_scale = solver + .state() + .y + .as_slice() + .iter() + .fold(0.0_f64, |scale, value| scale.max(value.abs())); + let material_tolerance = atol.abs() + rtol.abs() * state_scale; + if skipped_infusion > material_tolerance { + return Err(PharmsolError::OtherError(format!( + "coalescing stop times from t = {state_time:.16e} to \ + t = {stop_time:.16e} would skip infusion amount \ + {skipped_infusion:.6e}, above tolerance \ + {material_tolerance:.6e}" + ))); + } + *solver.state_mut().t = stop_time; + Ok(()) +} + +/// Advance the solver to `next_event_time`, stopping at every infusion +/// boundary in between. +/// +/// This is the single implementation of the event-to-event integration loop +/// shared by the closure-based [`ODE`] equation and the DSL runtime ODE path +/// ([`crate::dsl::native::NativeOdeModel`]): stop selection, left-continuity +/// handling at infusion boundaries, coalescing of stops that diffsol reports +/// as already reached, and bounded rescue restarts after recoverable step +/// failures. +/// +/// `after_step` runs after every step; the DSL path uses it to surface +/// model-function errors raised inside the RHS callback. +#[allow(clippy::too_many_arguments)] +pub(crate) fn advance_solver_to_event<'a, F, S, H>( + solver: &mut S, + next_event_time: f64, + infusion_boundary_cursor: &mut usize, + pending_reinit: &mut bool, + dy_scratch: &mut V, + rtol: f64, + atol: f64, + after_step: &mut H, +) -> Result<(), PharmsolError> +where + F: Fn(&V, &V, f64, &mut V, &V, &V, &Covariates) + 'a, + S: OdeSolverMethod<'a, PMProblem<'a, F>>, + H: FnMut() -> Result<(), PharmsolError>, +{ + while next_event_time > solver.state().t { + let infusion_boundary_times = solver.problem().eqn.infusion_boundary_times(); + 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) = + match infusion_boundary_times.get(*infusion_boundary_cursor) { + Some(&boundary_time) if boundary_time <= next_event_time => { + *infusion_boundary_cursor += 1; + (boundary_time, true) + } + _ => (next_event_time, false), + }; + + // Gaps the solver cannot integrate are coalesced up front; letting + // diffsol clamp its step to such a gap either fails immediately or + // leaves a crushed step size behind. + let state_time = solver.state().t; + if stop_time > state_time && stop_time - state_time < min_integrable_gap(solver, state_time) + { + coalesce_stop_to(solver, stop_time, rtol, atol)?; + *pending_reinit = true; + continue; + } + + 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 { + reinitialize_at_boundary(solver, dy_scratch); + *pending_reinit = false; + } + integrate_to_stop( + solver, + stop_time, + is_infusion_boundary, + pending_reinit, + dy_scratch, + after_step, + )?; + } + Err(diffsol::error::DiffsolError::OdeSolverError( + OdeSolverError::StopTimeAtCurrentTime, + )) => { + solver.problem().eqn.set_left_continuity_time(None); + // The requested stop is within diffsol's round-off of the + // current time. Coalesce forward stops; an actually earlier + // stop is a genuine ordering error. Snapping the logical time + // matters: a state left a few ULPs before an infusion + // boundary would otherwise be restarted with the old rate. + let state_time = solver.state().t; + if state_time > stop_time { + return Err(PharmsolError::from_solver_error( + diffsol::error::DiffsolError::OdeSolverError( + OdeSolverError::StopTimeBeforeCurrentTime { + stop_time, + state_time, + }, + ), + stop_time, + )); + } + coalesce_stop_to(solver, stop_time, rtol, atol)?; + *pending_reinit = true; + } + Err(err) => { + solver.problem().eqn.set_left_continuity_time(None); + return Err(PharmsolError::from_solver_error(err, stop_time)); + } + } + } + Ok(()) +} + +/// Step the solver until the stop set by `set_stop_time`, rescuing +/// recoverable failures with bounded in-place restarts. +/// +/// A recoverable failure leaves the state parked at the last accepted step, +/// typically right before a stiff transient or a mid-segment RHS +/// discontinuity that is not an event (e.g. a covariate change). Restarting +/// there with a fresh Jacobian at first order and a small step resets +/// diffsol's failure budgets and lets the solver walk across the transient. +/// Rescues that do not advance the state shrink the restart step tenfold; the +/// attempt and stall caps bound the retry work, and exhaustion returns the +/// original error annotated with the restart count. +fn integrate_to_stop<'a, F, S, H>( + solver: &mut S, + stop_time: f64, + is_infusion_boundary: bool, + pending_reinit: &mut bool, + dy_scratch: &mut V, + after_step: &mut H, +) -> Result<(), PharmsolError> +where + F: Fn(&V, &V, f64, &mut V, &V, &V, &Covariates) + 'a, + S: OdeSolverMethod<'a, PMProblem<'a, F>>, + H: FnMut() -> Result<(), PharmsolError>, +{ + let mut rescues = 0usize; + let mut stalled_rescues = 0usize; + let mut last_rescue_time: Option = None; + + loop { + match solver.step() { + Ok(OdeSolverStopReason::InternalTimestep) => { + after_step()?; + } + Ok(OdeSolverStopReason::TstopReached) => { + after_step()?; + solver.problem().eqn.set_left_continuity_time(None); + if solver.state().t != stop_time { + *solver.state_mut().t = stop_time; + *pending_reinit = true; + } + if is_infusion_boundary { + *pending_reinit = true; + } + return Ok(()); + } + Ok(OdeSolverStopReason::RootFound(root_time, _)) => { + return Err(PharmsolError::OtherError(format!( + "solver stopped at an unexpected root at t = {:.4} \ + (root finding is not configured)", + root_time + ))); + } + Err(diffsol::error::DiffsolError::OdeSolverError( + OdeSolverError::StopTimeAtCurrentTime, + )) => { + // A restart re-arms the pending stop inside diffsol, and the + // state can sit within diffsol's round-off of it; the stop is + // then effectively reached. + after_step()?; + solver.problem().eqn.set_left_continuity_time(None); + if solver.state().t != stop_time { + *solver.state_mut().t = stop_time; + *pending_reinit = true; + } + if is_infusion_boundary { + *pending_reinit = true; + } + return Ok(()); + } + Err(err) => { + // A model-function error raised inside the RHS is the root + // cause when present; surface it over the solver error. + after_step()?; + let state_time = solver.state().t; + if !recoverable_step_failure(&err) + || rescues >= MAX_RESCUES_PER_SEGMENT + || !solver_state_is_finite(solver) + { + solver.problem().eqn.set_left_continuity_time(None); + return Err(PharmsolError::from_solver_error(err, stop_time) + .with_rescue_context(rescues)); + } + match last_rescue_time { + Some(previous) if state_time <= previous + progress_epsilon(previous) => { + stalled_rescues += 1; + if stalled_rescues > MAX_STALLED_RESCUES { + solver.problem().eqn.set_left_continuity_time(None); + return Err(PharmsolError::from_solver_error(err, stop_time) + .with_rescue_context(rescues)); + } + } + _ => stalled_rescues = 0, + } + last_rescue_time = Some(state_time); + rescues += 1; + let restart_step = rescue_step_size(state_time, stop_time, stalled_rescues); + tracing::debug!( + "rescuing ODE solve at t = {state_time:.6e} toward t = {stop_time:.6e} \ + (attempt {rescues}, restart step {restart_step:.3e}): {err}" + ); + restart_solver_at_current_state(solver, dy_scratch, Some(restart_step)); + } + } + } } impl ODE { @@ -627,7 +978,6 @@ impl ODE { F: Fn(&V, &V, f64, &mut V, &V, &V, &Covariates) + 'a, S: OdeSolverMethod<'a, PMProblem<'a, F>>, { - let infusion_boundary_times = solver.problem().eqn.infusion_boundary_times(); let mut infusion_boundary_cursor = 0usize; let mut index = 0usize; // Set when the previous event changed the state or the previous stop @@ -717,105 +1067,16 @@ impl ODE { // 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 { - reinitialize_at_boundary(solver, dy_scratch); - pending_reinit = false; - } - loop { - match solver.step() { - 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 = 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. Breaking early would skip - // the remaining interval and leave the solver - // state at the boundary when the observation is - // evaluated. - 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)); - } - } - } + advance_solver_to_event( + solver, + next_event.time(), + &mut infusion_boundary_cursor, + &mut pending_reinit, + dy_scratch, + self.rtol, + self.atol, + &mut || Ok(()), + )?; } index += 1; } @@ -1349,8 +1610,8 @@ mod tests { // The infusion ends exactly one ULP after the observation at t = 10. // After landing on the observation stop the solver is already within // diffsol's round-off of the end boundary, so `set_stop_time` reports - // `StopTimeAtCurrentTime` and the loop accepts it through the - // same-ULP check. The reached boundary is *before* the observation at + // `StopTimeAtCurrentTime`, which confirms the boundary is reached. + // The reached boundary is *before* the observation at // t = 20, so the event loop must keep integrating toward it; breaking // early would evaluate the observation with the state frozen at the // boundary and miss the exponential decay. diff --git a/tests/bolus_reinit_stop_time.rs b/tests/bolus_reinit_stop_time.rs new file mode 100644 index 00000000..95f4c292 --- /dev/null +++ b/tests/bolus_reinit_stop_time.rs @@ -0,0 +1,383 @@ +//! Regression coverage for accepted stop times around state and RHS discontinuities. +//! +//! A solver restart can land a few ULPs from a requested event or infusion +//! boundary while diffsol still correctly reports that stop as reached. The +//! event loop must accept diffsol's `StopTimeAtCurrentTime`, align the logical +//! state time with the accepted stop, and restart with the post-boundary RHS. + +use pharmsol::prelude::*; + +#[cfg(feature = "dsl-jit")] +use pharmsol::dsl::{ + compile_module_source_to_runtime, CompiledRuntimeModel, RuntimeCompilationTarget, +}; + +const OBSERVATION_TIMES: [f64; 15] = [ + 0.0, + 0.35, + 0.516666666666667, + 0.983333333333333, + 1.48333333333333, + 2.0, + 2.5, + 3.0, + 4.0, + 4.98333333333333, + 6.98333333333333, + 7.98333333333333, + 10.0, + 11.0, + 12.0, +]; + +fn solver_cases() -> [(&'static str, OdeSolver); 4] { + [ + ("BDF", OdeSolver::Bdf), + ("TSIT45", OdeSolver::ExplicitRk(ExplicitRkTableau::Tsit45)), + ("TRBDF2", OdeSolver::Sdirk(SdirkTableau::TrBdf2)), + ("ESDIRK34", OdeSolver::Sdirk(SdirkTableau::Esdirk34)), + ] +} + +const PARAMETERS: [(&str, f64); 5] = [ + ("ka", 3.6156922578811646), + ("cl0", 1.0289061069488525), + ("vc0", 187.13204860687256), + ("q0", 2.4602913856506348), + ("vp0", 58.32162380218506), +]; + +fn subject_with_bolus_history() -> Subject { + let mut builder = Subject::builder("g34"); + for dose_index in -16..=0 { + builder = builder.bolus(f64::from(dose_index) * 12.0, 750.0, "input_1"); + } + for time in OBSERVATION_TIMES { + builder = builder.missing_observation(time, "outeq_1"); + } + builder.build() +} + +fn closure_model(solver: OdeSolver) -> equation::ODE { + equation::ODE::new( + |x, p, _t, dx, b, rateiv, _cov| { + fetch_params!(p, ka, cl0, vc0, q0, vp0); + let ke = cl0 / vc0; + let k23 = q0 / vc0; + let k32 = q0 / vp0; + + dx[0] = b[0] - x[0] * ka; + dx[1] = rateiv[0] + x[0] * ka + x[2] * k32 - x[1] * (ke + k23); + dx[2] = x[1] * k23 - x[2] * k32; + }, + |_p, _t, _cov| lag! {}, + |_p, _t, _cov| fa! {}, + |_p, _t, _cov, _x| {}, + |x, p, _t, _cov, y| { + fetch_params!(p, _ka, _cl0, vc0, _q0, _vp0); + y[0] = x[1] / vc0; + }, + ) + .with_nstates(3) + .with_ndrugs(1) + .with_nout(1) + .with_solver(solver) + .with_metadata( + equation::metadata::new("bolus_reinit_stop_time") + .parameters(["ka", "cl0", "vc0", "q0", "vp0"]) + .states(["x1", "x2", "x3"]) + .outputs(["outeq_1"]) + .routes([equation::Route::bolus("input_1") + .to_state("x1") + .expect_explicit_input()]), + ) + .expect("regression model metadata should validate") +} + +#[test] +fn closure_solvers_accept_reached_stop_after_bolus_restarts( +) -> Result<(), Box> { + for (label, solver) in solver_cases() { + let model = closure_model(solver); + let parameters = Parameters::with_model(&model, PARAMETERS)?; + let predictions = model + .estimate_predictions_dense(&subject_with_bolus_history(), parameters.as_slice()) + .unwrap_or_else(|error| panic!("{label}: bolus-history simulation failed: {error}")); + + assert_eq!(predictions.predictions().len(), OBSERVATION_TIMES.len()); + assert!(predictions + .predictions() + .iter() + .all(|prediction| prediction.prediction().is_finite())); + } + Ok(()) +} + +#[cfg(feature = "dsl-jit")] +const DSL_MODEL: &str = r#" +name = bolus_reinit_stop_time +kind = ode +params = ka, cl0, vc0, q0, vp0 +states = x1, x2, x3 +outputs = outeq_1 + +bolus(input_1) -> x1 +infusion(input_1) -> x2 + +cl = cl0 +vc = vc0 +q = q0 +vp = vp0 +ke = cl / vc +k23 = q / vc +k32 = q / vp + +dx(x1) = -(x1 * ka) +dx(x2) = x1 * ka + (x3 * k32) - (x2 * (ke + k23)) +dx(x3) = x2 * k23 - (x3 * k32) + +out(outeq_1) = x2 / vc +"#; + +fn infusion_boundary_subject() -> Subject { + Subject::builder("accepted_infusion_boundary") + .infusion(5.0, 100.0, "input_1", 10.0_f64.next_up() - 5.0) + .missing_observation(10.0, "cp") + .missing_observation(20.0, "cp") + .build() +} + +fn stepped_infusion_boundary_subject() -> Subject { + let mut builder = + Subject::builder("stepped_infusion_boundary").infusion(0.0, 100.0, "input_1", 12.0); + for time in OBSERVATION_TIMES.into_iter().filter(|time| *time < 12.0) { + builder = builder.missing_observation(time, "cp"); + } + builder.missing_observation(20.0, "cp").build() +} + +fn closure_infusion_model(solver: OdeSolver) -> equation::ODE { + equation::ODE::new( + |x, _p, _t, dx, _b, rateiv, _cov| { + dx[0] = rateiv[0] - 0.5 * x[0]; + }, + |_p, _t, _cov| lag! {}, + |_p, _t, _cov| fa! {}, + |_p, _t, _cov, _x| {}, + |x, _p, _t, _cov, y| y[0] = x[0], + ) + .with_nstates(1) + .with_ndrugs(1) + .with_nout(1) + .with_solver(solver) + .with_metadata( + equation::metadata::new("accepted_infusion_boundary") + .states(["central"]) + .outputs(["cp"]) + .routes([equation::Route::infusion("input_1").to_state("central")]), + ) + .expect("infusion metadata should validate") +} + +fn accepted_boundary_expected() -> f64 { + let delivered = 100.0 * (1.0 - (-2.5_f64).exp()) / 2.5; + delivered * (-5.0_f64).exp() +} + +fn stepped_boundary_expected() -> f64 { + (100.0 / 12.0) / 0.5 * (1.0 - (-6.0_f64).exp()) * (-4.0_f64).exp() +} + +fn assert_post_infusion_decay( + label: &str, + predictions: &SubjectPredictions, + expected: f64, + maximum_relative_error: f64, +) { + assert!(!predictions.predictions().is_empty(), "{label}"); + let actual = predictions.predictions().last().unwrap().prediction(); + let relative_error = (actual - expected).abs() / expected; + assert!( + relative_error < maximum_relative_error, + "{label}: post-infusion prediction {actual:.16e}, expected {expected:.16e}, relative error {relative_error:.3e}" + ); +} + +#[test] +fn closure_solvers_restart_with_post_infusion_rhs() { + for (label, solver) in solver_cases() { + let model = closure_infusion_model(solver); + for (scenario, subject, expected) in [ + ( + "accepted", + infusion_boundary_subject(), + accepted_boundary_expected(), + ), + ( + "stepped", + stepped_infusion_boundary_subject(), + stepped_boundary_expected(), + ), + ] { + let predictions = model + .estimate_predictions_dense(&subject, &[]) + .unwrap_or_else(|error| { + panic!("closure {label} {scenario} infusion failed: {error}") + }); + let maximum_relative_error = if label == "TSIT45" { 1.0e-3 } else { 1.0e-2 }; + assert_post_infusion_decay( + &format!("closure {label} {scenario}"), + &predictions, + expected, + maximum_relative_error, + ); + } + } +} + +#[cfg(feature = "dsl-jit")] +const INFUSION_DSL_MODEL: &str = r#" +name = accepted_infusion_boundary +kind = ode +states = central +outputs = cp +infusion(input_1) -> central +dx(central) = -(0.5 * central) +out(cp) = central +"#; + +#[test] +#[cfg(feature = "dsl-jit")] +fn jit_solvers_restart_with_post_infusion_rhs() -> Result<(), Box> { + for (label, solver) in solver_cases() { + let compiled = compile_module_source_to_runtime( + INFUSION_DSL_MODEL, + Some("accepted_infusion_boundary"), + RuntimeCompilationTarget::Jit, + |_, _| {}, + )?; + let model = match compiled { + CompiledRuntimeModel::Ode(model) => { + CompiledRuntimeModel::Ode(model.with_solver(solver)) + } + _ => return Err("expected an ODE model".into()), + }; + for (scenario, subject, expected) in [ + ( + "accepted", + infusion_boundary_subject(), + accepted_boundary_expected(), + ), + ( + "stepped", + stepped_infusion_boundary_subject(), + stepped_boundary_expected(), + ), + ] { + let predictions = match &model { + CompiledRuntimeModel::Ode(model) => model + .estimate_predictions_dense(&subject, &[]) + .unwrap_or_else(|error| { + panic!("JIT {label} {scenario} infusion failed: {error}") + }), + _ => unreachable!(), + }; + let maximum_relative_error = if label == "TSIT45" { 1.0e-3 } else { 1.0e-2 }; + assert_post_infusion_decay( + &format!("JIT {label} {scenario}"), + &predictions, + expected, + maximum_relative_error, + ); + } + } + Ok(()) +} + +fn material_short_infusion_subject() -> Subject { + Subject::builder("material_short_infusion") + .infusion(1.0, 100.0, "input_1", 1.0_f64.next_up() - 1.0) + .missing_observation(2.0, "cp") + .build() +} + +fn assert_material_infusion_error(label: &str, error: PharmsolError) { + let message = error.to_string(); + assert!( + message.contains("would skip infusion amount"), + "{label}: unexpected error: {message}" + ); +} + +#[test] +fn closure_solvers_reject_material_coalesced_infusions() { + for (label, solver) in solver_cases() { + let model = closure_infusion_model(solver); + let error = model + .estimate_predictions_dense(&material_short_infusion_subject(), &[]) + .unwrap_err(); + assert_material_infusion_error(&format!("closure {label}"), error); + } +} + +#[test] +#[cfg(feature = "dsl-jit")] +fn jit_solvers_reject_material_coalesced_infusions() -> Result<(), Box> { + for (label, solver) in solver_cases() { + let compiled = compile_module_source_to_runtime( + INFUSION_DSL_MODEL, + Some("accepted_infusion_boundary"), + RuntimeCompilationTarget::Jit, + |_, _| {}, + )?; + let model = match compiled { + CompiledRuntimeModel::Ode(model) => { + CompiledRuntimeModel::Ode(model.with_solver(solver)) + } + _ => return Err("expected an ODE model".into()), + }; + let error = match &model { + CompiledRuntimeModel::Ode(model) => model + .estimate_predictions_dense(&material_short_infusion_subject(), &[]) + .unwrap_err(), + _ => unreachable!(), + }; + assert_material_infusion_error(&format!("JIT {label}"), error); + } + Ok(()) +} + +#[test] +#[cfg(feature = "dsl-jit")] +fn jit_solvers_accept_reached_stop_after_bolus_restarts() -> Result<(), Box> +{ + for (label, solver) in solver_cases() { + let compiled = compile_module_source_to_runtime( + DSL_MODEL, + Some("bolus_reinit_stop_time"), + RuntimeCompilationTarget::Jit, + |_, _| {}, + )?; + let model = match compiled { + CompiledRuntimeModel::Ode(model) => { + CompiledRuntimeModel::Ode(model.with_solver(solver)) + } + _ => return Err("expected an ODE model".into()), + }; + let parameters = Parameters::with_model(&model, PARAMETERS)?; + + let predictions = match &model { + CompiledRuntimeModel::Ode(model) => model + .estimate_predictions_dense(&subject_with_bolus_history(), parameters.as_slice()) + .unwrap_or_else(|error| panic!("{label}: JIT bolus-history failed: {error}")), + _ => unreachable!(), + }; + + assert_eq!(predictions.predictions().len(), OBSERVATION_TIMES.len()); + assert!(predictions + .predictions() + .iter() + .all(|prediction| prediction.prediction().is_finite())); + } + Ok(()) +} diff --git a/tests/extreme_stiffness.rs b/tests/extreme_stiffness.rs new file mode 100644 index 00000000..7794d352 --- /dev/null +++ b/tests/extreme_stiffness.rs @@ -0,0 +1,490 @@ +//! Extreme-stiffness regression tests for the ODE event loop. +//! +//! These target the failure modes seen in population fits (PMcore): sudden +//! RHS changes at infusion boundaries, lagged boluses, and steep covariate +//! changes drive the adaptive step size toward zero. diffsol also shrinks the +//! step to land exactly on stops that sit a few dozen ULPs apart and swallows +//! the resulting `StepSizeTooSmall`, after which every later step-size update +//! fails. The event loop must survive all of this — for every solver — by +//! restarting in place, and must return a descriptive error (not hang) when a +//! problem is genuinely beyond the selected solver. + +use pharmsol::prelude::*; + +/// Shared relative/absolute tolerance for analytical comparisons: well above +/// the solver tolerances (1e-4) but far below any qualitative difference. +const PRED_TOLERANCE: f64 = 1e-3; + +fn implicit_solvers() -> [(&'static str, OdeSolver); 3] { + [ + ("BDF", OdeSolver::Bdf), + ("TRBDF2", OdeSolver::Sdirk(SdirkTableau::TrBdf2)), + ("ESDIRK34", OdeSolver::Sdirk(SdirkTableau::Esdirk34)), + ] +} + +fn all_solvers() -> [(&'static str, OdeSolver); 4] { + [ + ("BDF", OdeSolver::Bdf), + ("TRBDF2", OdeSolver::Sdirk(SdirkTableau::TrBdf2)), + ("ESDIRK34", OdeSolver::Sdirk(SdirkTableau::Esdirk34)), + ("TSIT45", OdeSolver::ExplicitRk(ExplicitRkTableau::Tsit45)), + ] +} + +fn assert_close(label: &str, actual: f64, expected: f64) { + assert!( + actual.is_finite(), + "{label}: prediction is not finite (expected {expected:.6e})" + ); + let scale = expected.abs().max(1.0); + assert!( + (actual - expected).abs() <= PRED_TOLERANCE * scale, + "{label}: prediction {actual:.6e} differs from analytical {expected:.6e} \ + by more than {PRED_TOLERANCE:.0e} (scale {scale:.3e})" + ); +} + +fn predictions_for( + label: &str, + model: &equation::ODE, + subject: &Subject, + parameters: &[(&str, f64)], +) -> Vec { + let parameters = Parameters::with_model(model, parameters.iter().copied()) + .unwrap_or_else(|error| panic!("{label}: parameters should validate: {error}")); + let predictions = model + .estimate_predictions_dense(subject, parameters.as_slice()) + .unwrap_or_else(|error| panic!("{label}: simulation failed: {error}")); + predictions + .predictions() + .iter() + .map(|prediction| prediction.prediction()) + .collect() +} + +// --------------------------------------------------------------------------- +// Stiff infusion boundary vs the analytical one-compartment solution. +// --------------------------------------------------------------------------- + +fn infusion_model(solver: OdeSolver) -> equation::ODE { + equation::ODE::new( + |x, p, _t, dx, _b, rateiv, _cov| { + fetch_params!(p, ke); + dx[0] = rateiv[0] - ke * x[0]; + }, + |_p, _t, _cov| lag! {}, + |_p, _t, _cov| fa! {}, + |_p, _t, _cov, _x| {}, + |x, _p, _t, _cov, y| y[0] = x[0], + ) + .with_nstates(1) + .with_ndrugs(1) + .with_nout(1) + .with_solver(solver) + .with_metadata( + equation::metadata::new("stiff_infusion") + .parameters(["ke"]) + .states(["central"]) + .outputs(["cp"]) + .routes([equation::Route::infusion("iv").to_state("central")]), + ) + .expect("stiff infusion metadata should validate") +} + +/// Unit-height infusion: rate `ke` over [0, 1] makes the steady state 1.0. +fn stiff_infusion_subject(observation_times: &[f64], ke: f64) -> Subject { + let mut builder = Subject::builder("stiff_infusion").infusion(0.0, ke, "iv", 1.0); + for &time in observation_times { + builder = builder.missing_observation(time, "cp"); + } + builder.build() +} + +fn exact_unit_infusion(ke: f64, t: f64) -> f64 { + if t <= 1.0 { + 1.0 - (-ke * t).exp() + } else { + (1.0 - (-ke).exp()) * (-ke * (t - 1.0)).exp() + } +} + +/// Observations straddling the infusion end at t = 1, with post-boundary +/// offsets scaled to the elimination timescale so the fast transient itself +/// is sampled. +fn boundary_observation_times(ke: f64) -> Vec { + vec![ + 0.5, + 1.0, + 1.0 + 0.5 / ke, + 1.0 + 1.0 / ke, + 1.0 + 2.0 / ke, + 1.0 + 5.0 / ke, + 2.0, + ] +} + +fn assert_stiff_infusion_boundary(label: &str, solver: OdeSolver, ke: f64) { + let times = boundary_observation_times(ke); + let subject = stiff_infusion_subject(×, ke); + let model = infusion_model(solver); + let actual = predictions_for(label, &model, &subject, &[("ke", ke)]); + assert_eq!(actual.len(), times.len(), "{label}: prediction count"); + for (time, actual) in times.iter().zip(actual) { + let expected = exact_unit_infusion(ke, *time); + assert_close(&format!("{label} at t = {time}"), actual, expected); + } +} + +#[test] +fn stiff_infusion_boundary_matches_analytical_implicit() { + for (name, solver) in implicit_solvers() { + for ke in [1e2, 1e4, 1e6, 1e8] { + assert_stiff_infusion_boundary(&format!("{name} ke={ke:.0e}"), solver.clone(), ke); + } + } +} + +#[test] +fn stiff_infusion_boundary_matches_analytical_explicit() { + // Explicit RK is stability-limited; keep the stiffness within what it can + // integrate in reasonable time. + for ke in [1e2, 1e4] { + assert_stiff_infusion_boundary( + &format!("TSIT45 ke={ke:.0e}"), + OdeSolver::ExplicitRk(ExplicitRkTableau::Tsit45), + ke, + ); + } +} + +// --------------------------------------------------------------------------- +// Step size crushed by near-coincident stop times. +// +// diffsol shrinks the step to land exactly on a close stop and ignores the +// `StepSizeTooSmall` this raises; once the step sits below half the minimum, +// every later step-size update (even growth, capped at 2x) fails. These used +// to abort the simulation a few steps after the close pair. +// --------------------------------------------------------------------------- + +#[test] +fn ulp_scale_gap_between_observations_recovers() { + // Two observations 3e-14 apart at t = 0.1: too far apart for diffsol to + // report the second stop as already reached, close enough to crush the + // step size below recovery. + let close = 0.1_f64 + 3e-14; + assert!(close > 0.1, "gap must be representable"); + let ke = 0.5; + let times = [0.1, close, 2.0]; + + for (name, solver) in all_solvers() { + let label = format!("{name} ulp-gap observations"); + let subject = stiff_infusion_subject(×, ke); + let model = infusion_model(solver); + let actual = predictions_for(&label, &model, &subject, &[("ke", ke)]); + assert_eq!(actual.len(), times.len(), "{label}: prediction count"); + for (time, actual) in times.iter().zip(actual) { + let expected = exact_unit_infusion(ke, *time); + assert_close(&format!("{label} at t = {time}"), actual, expected); + } + } +} + +#[test] +fn ulp_scale_gap_between_observation_and_infusion_end_recovers() { + // The infusion ends 3e-14 after an observation, so the crushed step is + // carried into an infusion-boundary restart instead of a plain segment. + let duration = 0.1_f64 + 3e-14; + let ke = 0.5; + let times = [0.1, 2.0]; + + for (name, solver) in all_solvers() { + let label = format!("{name} ulp-gap infusion end"); + let mut builder = + Subject::builder("ulp_gap_infusion").infusion(0.0, ke * duration, "iv", duration); + for &time in × { + builder = builder.missing_observation(time, "cp"); + } + let subject = builder.build(); + let model = infusion_model(solver); + let actual = predictions_for(&label, &model, &subject, &[("ke", ke)]); + + // Same unit-height model, infusion just ends at `duration` instead of 1. + let exact = |t: f64| { + if t <= duration { + 1.0 - (-ke * t).exp() + } else { + (1.0 - (-ke * duration).exp()) * (-ke * (t - duration)).exp() + } + }; + assert_eq!(actual.len(), times.len(), "{label}: prediction count"); + for (time, actual) in times.iter().zip(actual) { + assert_close(&format!("{label} at t = {time}"), actual, exact(*time)); + } + } +} + +// --------------------------------------------------------------------------- +// Long free segment after a stiff infusion cut-off. +// +// With no observation near the boundary the restart begins with the large +// step size inherited from the smooth infusion phase, and the controller must +// reject its way down to the fast elimination timescale — the classic way to +// exhaust diffsol's error-test budget on very stiff parameter draws. +// --------------------------------------------------------------------------- + +#[test] +fn stiff_cutoff_with_distant_observation_recovers() { + for (name, solver) in implicit_solvers() { + for ke in [1e6, 1e9] { + let label = format!("{name} ke={ke:.0e} distant observation"); + // Rate ke over [0, 24] -> steady state 1.0; nothing stops the + // solver between the boundary at 24 and the observation at 48. + let subject = Subject::builder("stiff_cutoff") + .infusion(0.0, ke * 24.0, "iv", 24.0) + .missing_observation(48.0, "cp") + .build(); + let model = infusion_model(solver.clone()); + let actual = predictions_for(&label, &model, &subject, &[("ke", ke)]); + assert_eq!(actual.len(), 1, "{label}: prediction count"); + // Fully eliminated 24 time units after cut-off. + assert_close(&label, actual[0], 0.0); + } + } +} + +// --------------------------------------------------------------------------- +// Steep covariate ramp in the middle of a segment (no event, no stop time). +// --------------------------------------------------------------------------- + +fn covariate_clearance_model(solver: OdeSolver) -> equation::ODE { + equation::ODE::new( + |x, _p, t, dx, b, _rateiv, cov| { + fetch_cov!(cov, t, cl); + dx[0] = b[0] - cl * x[0]; + }, + |_p, _t, _cov| lag! {}, + |_p, _t, _cov| fa! {}, + |_p, _t, _cov, _x| {}, + |x, _p, _t, _cov, y| y[0] = x[0], + ) + .with_nstates(1) + .with_ndrugs(1) + .with_nout(1) + .with_solver(solver) + .with_metadata( + equation::metadata::new("covariate_clearance") + .parameters(["dummy"]) + .states(["central"]) + .outputs(["cp"]) + .routes([equation::Route::bolus("dose") + .to_state("central") + .expect_explicit_input()]), + ) + .expect("covariate clearance metadata should validate") +} + +fn assert_covariate_ramp(label: &str, solver: OdeSolver, lambda: f64) { + // Clearance ramps linearly from 1 to `lambda` over a window of width + // 2/lambda starting at t = 5 — a near-discontinuous covariate change in + // the middle of the [4, ...] segment. The integral of the ramp is + // (1 + lambda)/lambda, so a bolus of exp(5 + (1 + lambda)/lambda) makes + // the state exactly 1.0 at the end of the ramp. + let ramp_end = 5.0 + 2.0 / lambda; + let bolus = (5.0 + (1.0 + lambda) / lambda).exp(); + let times = [4.0, ramp_end + 1.0 / lambda, ramp_end + 3.0 / lambda]; + let expected = [bolus * (-4.0_f64).exp(), (-1.0_f64).exp(), (-3.0_f64).exp()]; + + let subject = { + let mut builder = Subject::builder("covariate_ramp") + .bolus(0.0, bolus, "dose") + .covariate("cl", 0.0, 1.0) + .covariate("cl", 5.0, 1.0) + .covariate("cl", ramp_end, lambda); + for &time in × { + builder = builder.missing_observation(time, "cp"); + } + builder.build() + }; + + let model = covariate_clearance_model(solver); + let actual = predictions_for(label, &model, &subject, &[("dummy", 1.0)]); + assert_eq!(actual.len(), times.len(), "{label}: prediction count"); + for ((time, actual), expected) in times.iter().zip(actual).zip(expected) { + assert_close(&format!("{label} at t = {time}"), actual, expected); + } +} + +#[test] +fn steep_covariate_ramp_mid_segment_recovers() { + for (name, solver) in all_solvers() { + assert_covariate_ramp(&format!("{name} lambda=1e4"), solver, 1e4); + } + for (name, solver) in implicit_solvers() { + assert_covariate_ramp(&format!("{name} lambda=1e6"), solver, 1e6); + } +} + +// --------------------------------------------------------------------------- +// Lagged bolus into a very stiff elimination. +// --------------------------------------------------------------------------- + +fn lagged_bolus_model(solver: OdeSolver) -> equation::ODE { + equation::ODE::new( + |x, p, _t, dx, b, _rateiv, _cov| { + fetch_params!(p, ke, _tlag); + dx[0] = b[0] - ke * x[0]; + }, + |p, _t, _cov| { + fetch_params!(p, _ke, tlag); + lag! {0 => tlag} + }, + |_p, _t, _cov| fa! {}, + |_p, _t, _cov, _x| {}, + |x, _p, _t, _cov, y| y[0] = x[0], + ) + .with_nstates(1) + .with_ndrugs(1) + .with_nout(1) + .with_solver(solver) + .with_metadata( + equation::metadata::new("lagged_bolus") + .parameters(["ke", "tlag"]) + .states(["central"]) + .outputs(["cp"]) + .routes([equation::Route::bolus("dose") + .to_state("central") + .expect_explicit_input()]), + ) + .expect("lagged bolus metadata should validate") +} + +#[test] +fn lagged_bolus_into_stiff_elimination_matches_analytical() { + let tlag = 0.5; + let dose_time = 1.0; + let arrival = dose_time + tlag; + let amount = 100.0; + + for (name, solver) in implicit_solvers() { + for ke in [1e4, 1e8] { + let label = format!("{name} ke={ke:.0e} lagged bolus"); + let times = [ + 1.0, + arrival + 0.5 / ke, + arrival + 1.0 / ke, + arrival + 3.0 / ke, + arrival + 10.0, + ]; + let mut builder = Subject::builder("lagged_bolus").bolus(dose_time, amount, "dose"); + for &time in × { + builder = builder.missing_observation(time, "cp"); + } + let subject = builder.build(); + let model = lagged_bolus_model(solver.clone()); + let actual = predictions_for(&label, &model, &subject, &[("ke", ke), ("tlag", tlag)]); + assert_eq!(actual.len(), times.len(), "{label}: prediction count"); + for (time, actual) in times.iter().zip(actual) { + let expected = if *time < arrival { + 0.0 + } else { + amount * (-ke * (time - arrival)).exp() + }; + assert_close(&format!("{label} at t = {time}"), actual, expected); + } + } + } +} + +// --------------------------------------------------------------------------- +// Michaelis-Menten depletion: stiffness switches on as the state crosses Km, +// far from any event. All solvers must agree with each other. +// --------------------------------------------------------------------------- + +fn michaelis_menten_model(solver: OdeSolver) -> equation::ODE { + equation::ODE::new( + |x, p, _t, dx, _b, rateiv, _cov| { + fetch_params!(p, vmax, km); + dx[0] = rateiv[0] - vmax * x[0] / (km + x[0]); + }, + |_p, _t, _cov| lag! {}, + |_p, _t, _cov| fa! {}, + |_p, _t, _cov, _x| {}, + |x, _p, _t, _cov, y| y[0] = x[0], + ) + .with_nstates(1) + .with_ndrugs(1) + .with_nout(1) + .with_solver(solver) + .with_metadata( + equation::metadata::new("michaelis_menten") + .parameters(["vmax", "km"]) + .states(["central"]) + .outputs(["cp"]) + .routes([equation::Route::infusion("iv").to_state("central")]), + ) + .expect("michaelis-menten metadata should validate") +} + +#[test] +fn michaelis_menten_depletion_agrees_across_solvers() { + // Saturated elimination consumes the 10-unit infusion at ~1/unit time, so + // the state crosses Km (1e-2) around t = 10 and the local stiffness jumps + // from ~0 to vmax/km = 1e2 mid-segment. Solver tolerances are tightened so + // cross-solver agreement is meaningful through the depletion corner. + let parameters = [("vmax", 1.0), ("km", 1e-2)]; + let times = [0.5, 1.0, 5.0, 9.5, 9.9, 10.0, 10.05, 10.5, 11.0]; + let subject = { + let mut builder = Subject::builder("michaelis_menten").infusion(0.0, 10.0, "iv", 1.0); + for &time in × { + builder = builder.missing_observation(time, "cp"); + } + builder.build() + }; + let model_for = |solver: OdeSolver| michaelis_menten_model(solver).with_tolerances(1e-6, 1e-6); + + let reference = predictions_for( + "BDF michaelis-menten", + &model_for(OdeSolver::Bdf), + &subject, + ¶meters, + ); + assert_eq!(reference.len(), times.len()); + + for (name, solver) in all_solvers() { + let label = format!("{name} michaelis-menten"); + let actual = predictions_for(&label, &model_for(solver), &subject, ¶meters); + assert_eq!(actual.len(), reference.len(), "{label}: prediction count"); + for ((time, actual), expected) in times.iter().zip(&actual).zip(&reference) { + assert_close(&format!("{label} at t = {time}"), *actual, *expected); + } + } +} + +// --------------------------------------------------------------------------- +// Beyond rescue: a problem stiffer than the solver's hard minimum step must +// fail with a descriptive error, not hang or panic. +// --------------------------------------------------------------------------- + +#[test] +fn impossibly_stiff_problem_returns_descriptive_error() { + // ke = 1e12 needs accuracy steps below diffsol's minimum step size; no + // number of restarts can integrate the transient. + let ke = 1e12; + let subject = Subject::builder("impossibly_stiff") + .infusion(0.0, ke * 24.0, "iv", 24.0) + .missing_observation(48.0, "cp") + .build(); + let model = infusion_model(OdeSolver::Bdf); + let parameters = + Parameters::with_model(&model, [("ke", ke)]).expect("parameters should validate"); + + let error = model + .estimate_predictions_dense(&subject, parameters.as_slice()) + .expect_err("a problem beyond the minimum step size should fail"); + let message = error.to_string(); + assert!( + message.contains("did not recover"), + "error should mention the exhausted restarts: {message}" + ); +}