diff --git a/src/dsl/backend.rs b/src/dsl/backend.rs index 94abe0a2..2942670f 100644 --- a/src/dsl/backend.rs +++ b/src/dsl/backend.rs @@ -22,6 +22,7 @@ use super::model_info::{RuntimeModelInfo, RuntimeRouteInfo, RuntimeStateInfo}; use crate::{ data::error_model::AssayErrorModels, data::{Covariates, Infusion, InputLabel, OutputLabel}, + error::NonFiniteRhs, simulator::{ cache::{ BoundErrorModelCache, PredictionCache, SdeLikelihoodCache, @@ -1194,9 +1195,11 @@ impl RuntimeOdeModel { ))); } let function_error = RefCell::new(None::); + let non_finite_rhs = RefCell::new(None::); let diffeq_session = &session; let diffeq_error = &function_error; + let diffeq_non_finite = &non_finite_rhs; let diffeq = move |x: &V, p: &V, t: f64, @@ -1243,6 +1246,22 @@ impl RuntimeOdeModel { dx.as_mut_slice().fill(0.0); } else { shared.apply_route_inputs_to_rates(dx.as_mut_slice(), rateiv.as_slice()); + // Track whether the most recent evaluation was non-finite so a + // later solver failure can report it. Non-finite trial + // evaluations are normally recovered by the solver, so this + // never turns into an error on its own. + let mut slot = diffeq_non_finite.borrow_mut(); + if dx.as_slice().iter().all(|value| value.is_finite()) { + if slot.is_some() { + *slot = None; + } + } else { + *slot = Some(NonFiniteRhs { + time: t, + state: x.as_slice().to_vec(), + derivative: dx.as_slice().to_vec(), + }); + } } }; @@ -1292,6 +1311,7 @@ impl RuntimeOdeModel { &mut output, &session, &function_error, + &non_finite_rhs, )?; }}; } @@ -1327,6 +1347,7 @@ impl RuntimeOdeModel { output: &mut SubjectPredictions, session: &RefCell>, function_error: &RefCell>, + non_finite_rhs: &RefCell>, ) -> Result<(), PharmsolError> where F: Fn(&V, &V, f64, &mut V, &V, &V, &Covariates) + 'a, @@ -1393,7 +1414,25 @@ impl RuntimeOdeModel { Some(error) => Err(error), None => Ok(()), }, - )?; + ) + .map_err(|error| { + let states = self.shared.metadata().states(); + let state_name = |i: usize| { + states + .get(i) + .map(|state| state.name().to_string()) + .unwrap_or_else(|| format!("x[{i}]")) + }; + let rhs = non_finite_rhs.borrow(); + let eqn = &solver.problem().eqn; + error.with_solver_diagnostics( + eqn.absolute_time(solver.state().t), + solver.state().y.as_slice(), + self.atol, + &state_name, + rhs.as_ref().map(|rhs| (rhs, eqn.absolute_time(rhs.time))), + ) + })?; } index += 1; } diff --git a/src/error/mod.rs b/src/error/mod.rs index 80d1201e..976bec1c 100644 --- a/src/error/mod.rs +++ b/src/error/mod.rs @@ -115,6 +115,125 @@ impl PharmsolError { } } +/// Snapshot of the most recent right-hand-side evaluation that produced +/// non-finite derivatives. +/// +/// Recorded by the compiled-model dynamics callback and cleared as soon as a +/// later evaluation is finite, so it is only still present when the solver gave +/// up straight after a non-finite evaluation, which is the usual cause of +/// "step size collapsed" and "Newton solver failed" errors. +#[cfg_attr(not(feature = "dsl"), allow(dead_code))] +#[derive(Debug, Clone)] +pub(crate) struct NonFiniteRhs { + /// Solver-local time of the evaluation. + pub time: f64, + pub state: Vec, + pub derivative: Vec, +} + +/// Ratio of largest to smallest non-zero state magnitude above which a solver +/// failure is annotated with a rescaling hint. +#[cfg_attr(not(feature = "dsl"), allow(dead_code))] +const STATE_RANGE_HINT_RATIO: f64 = 1e8; +/// Maximum number of entries listed in a diagnostic note. +#[cfg_attr(not(feature = "dsl"), allow(dead_code))] +const MAX_LISTED_STATES: usize = 6; + +#[cfg_attr(not(feature = "dsl"), allow(dead_code))] +impl PharmsolError { + /// Append actionable context to a solver failure. + /// + /// Adds (when applicable): the non-finite derivatives of the last RHS + /// evaluation together with the state that produced them, non-finite entries + /// in the last accepted state, and a rescaling hint when the accepted state + /// spans many orders of magnitude. Only [`PharmsolError::DiffsolError`] is + /// augmented; all other variants pass through unchanged. + /// + /// `state_name` maps a state index to a display label. `rhs` pairs the + /// recorded evaluation with its absolute time. + pub(crate) fn with_solver_diagnostics( + self, + accepted_time: f64, + accepted_state: &[f64], + atol: f64, + state_name: &dyn Fn(usize) -> String, + rhs: Option<(&NonFiniteRhs, f64)>, + ) -> Self { + let PharmsolError::DiffsolError(message) = self else { + return self; + }; + let mut notes = Vec::new(); + + if let Some((rhs, absolute_time)) = rhs { + let bad = rhs + .derivative + .iter() + .enumerate() + .filter(|(_, value)| !value.is_finite()) + .map(|(i, value)| format!("d{}/dt = {value}", state_name(i))) + .take(MAX_LISTED_STATES) + .collect::>() + .join(", "); + let state = list_states(&rhs.state, state_name); + notes.push(format!( + "the last right-hand-side evaluation (t = {absolute_time:.4}) returned non-finite derivatives ({bad}) for state [{state}]; check for division by zero, negative bases of fractional powers, or overflow. A state that overshoots slightly below zero inside a stiff step can turn terms such as (c/x)^h into NaN" + )); + } + + let non_finite = accepted_state + .iter() + .enumerate() + .filter(|(_, value)| !value.is_finite()) + .map(|(i, value)| format!("{} = {value}", state_name(i))) + .take(MAX_LISTED_STATES) + .collect::>(); + if !non_finite.is_empty() { + notes.push(format!( + "the last accepted state (t = {accepted_time:.4}) contains non-finite values ({})", + non_finite.join(", ") + )); + } else { + let magnitudes = accepted_state + .iter() + .enumerate() + .map(|(i, value)| (i, value.abs())) + .filter(|(_, magnitude)| *magnitude > 0.0); + let largest = magnitudes.clone().max_by(|a, b| a.1.total_cmp(&b.1)); + let smallest = magnitudes.min_by(|a, b| a.1.total_cmp(&b.1)); + if let (Some((hi, max)), Some((lo, min))) = (largest, smallest) { + if max / min >= STATE_RANGE_HINT_RATIO { + notes.push(format!( + "the last accepted state (t = {accepted_time:.4}) spans {:.0} orders of magnitude ({} = {max:.3e}, {} = {min:.3e}) with atol = {atol:e}; consider rescaling states (e.g. counts in units of 1e6) or bounding unchecked exponential growth", + (max / min).log10().floor(), + state_name(hi), + state_name(lo), + )); + } + } + } + + if notes.is_empty() { + return PharmsolError::DiffsolError(message); + } + PharmsolError::DiffsolError(format!("{message} Diagnostics: {}.", notes.join("; "))) + } +} + +#[cfg_attr(not(feature = "dsl"), allow(dead_code))] +fn list_states(state: &[f64], state_name: &dyn Fn(usize) -> String) -> String { + let mut listed = state + .iter() + .enumerate() + .take(MAX_LISTED_STATES) + .map(|(i, value)| format!("{} = {value:.4e}", state_name(i))) + .collect::>() + .join(", "); + if state.len() > MAX_LISTED_STATES { + listed.push_str(", ..."); + } + listed +} + /// Render a ` (available: a, b, c)` suffix, or an empty string when the list is /// empty (e.g. index-only models without named labels). fn format_available(labels: &[&str]) -> String { @@ -195,3 +314,72 @@ fn describe_diffsol_error( }, } } + +#[cfg(test)] +mod solver_diagnostics_tests { + use super::*; + + fn solver_error() -> PharmsolError { + PharmsolError::DiffsolError("Newton solver failed".to_string()) + } + + fn name(i: usize) -> String { + ["drug", "bacteria"] + .get(i) + .map_or_else(|| format!("x[{i}]"), |n| n.to_string()) + } + + #[test] + fn wide_state_range_adds_rescaling_hint() { + let message = solver_error() + .with_solver_diagnostics(3.0, &[1.5, 1.6e10], 1e-4, &name, None) + .to_string(); + assert!(message.contains("10 orders of magnitude"), "{message}"); + assert!(message.contains("bacteria = 1.600e10"), "{message}"); + assert!(message.contains("rescaling"), "{message}"); + } + + #[test] + fn moderate_state_range_is_not_annotated() { + let original = solver_error().to_string(); + let message = solver_error() + .with_solver_diagnostics(3.0, &[1.5, 1.0e5], 1e-4, &name, None) + .to_string(); + assert_eq!(message, original); + } + + #[test] + fn non_finite_accepted_state_is_named() { + let message = solver_error() + .with_solver_diagnostics(3.0, &[1.0, f64::NAN], 1e-4, &name, None) + .to_string(); + assert!(message.contains("bacteria = NaN"), "{message}"); + } + + #[test] + fn non_finite_rhs_lists_derivatives_and_state() { + let rhs = NonFiniteRhs { + time: 2.0, + state: vec![1.0, -1e-9], + derivative: vec![0.5, f64::NAN], + }; + let message = solver_error() + .with_solver_diagnostics(2.0, &[1.0, 1.0], 1e-4, &name, Some((&rhs, 12.0))) + .to_string(); + assert!(message.contains("dbacteria/dt = NaN"), "{message}"); + assert!(message.contains("t = 12.0000"), "{message}"); + assert!(message.contains("bacteria = -1.0000e-9"), "{message}"); + } + + #[test] + fn non_solver_errors_pass_through() { + let error = PharmsolError::OtherError("boom".to_string()).with_solver_diagnostics( + 0.0, + &[f64::NAN], + 1e-4, + &name, + None, + ); + assert_eq!(error.to_string(), "Other error: boom"); + } +} diff --git a/tests/dsl_solver_diagnostics.rs b/tests/dsl_solver_diagnostics.rs new file mode 100644 index 00000000..a6db42f8 --- /dev/null +++ b/tests/dsl_solver_diagnostics.rs @@ -0,0 +1,41 @@ +//! Solver failures on compiled DSL models should explain themselves: the +//! message names the offending states and, when it applies, a rescaling hint. + +#![cfg(feature = "dsl")] + +use pharmsol::dsl::compile_module_source_to_runtime; +use pharmsol::{prelude::*, Parameters}; + +fn simulate(source: &str, name: &str) -> Result<(), String> { + let model = compile_module_source_to_runtime(source, Some(name), |_, _| {}).unwrap(); + let support_point = Parameters::with_model(&model, [("k", 1.0)]).unwrap(); + let subject = Subject::builder("s1") + .missing_observation(0.0, "y") + .missing_observation(5.0, "y") + .build(); + model + .estimate_predictions(&subject, &support_point) + .map(|_| ()) + .map_err(|error| error.to_string()) +} + +/// `sqrt` of a negative number once `amount` undershoots 0.5 yields NaN. +#[test] +fn non_finite_rhs_is_reported_with_state_names() { + let source = r#" +name = blowup +kind = ode + +params = k +states = amount +outputs = y + +dx(amount) = -k * sqrt(amount - 0.5) +init(amount) = 1 +out(y) = amount +"#; + let message = simulate(source, "blowup").expect_err("model must fail"); + assert!(message.contains("Diagnostics:"), "{message}"); + assert!(message.contains("damount/dt"), "{message}"); + assert!(message.contains("amount ="), "{message}"); +}