From bbda6a6c945231087d742af1c5a73e747ca57100 Mon Sep 17 00:00:00 2001 From: Markus <66058642+mhovd@users.noreply.github.com> Date: Fri, 11 Sep 2026 14:31:51 +0200 Subject: [PATCH] wip: Update algorithm trait --- benches/bimodal_ke.rs | 2 +- examples/bimodal_ke_backend_compare.rs | 2 +- examples/chain_algorithms.rs | 4 +- src/algorithms/diagnostics.rs | 186 +++++++++++ src/algorithms/mod.rs | 338 +++++++------------- src/algorithms/nonparametric/controller.rs | 20 +- src/algorithms/nonparametric/error_optim.rs | 34 +- src/algorithms/nonparametric/mod.rs | 61 ++-- src/algorithms/nonparametric/ncnpag.rs | 139 ++++---- src/algorithms/nonparametric/npag.rs | 271 ++++++---------- src/algorithms/nonparametric/npmap.rs | 126 +++----- src/algorithms/nonparametric/npod.rs | 281 +++++++--------- src/algorithms/nonparametric/stages.rs | 150 +++++++++ src/algorithms/parametric/mod.rs | 23 +- src/algorithms/stop.rs | 28 ++ src/estimation/nonparametric/cycles.rs | 40 ++- src/estimation/nonparametric/result.rs | 45 ++- src/estimation/nonparametric/summaries.rs | 2 +- src/lib.rs | 2 +- src/results/fit_result.rs | 13 +- tests/onecomp.rs | 9 +- 21 files changed, 941 insertions(+), 835 deletions(-) create mode 100644 src/algorithms/diagnostics.rs create mode 100644 src/algorithms/nonparametric/stages.rs create mode 100644 src/algorithms/stop.rs diff --git a/benches/bimodal_ke.rs b/benches/bimodal_ke.rs index 97d805c91..5787f1992 100644 --- a/benches/bimodal_ke.rs +++ b/benches/bimodal_ke.rs @@ -68,7 +68,7 @@ fn setup_postprob() -> Result> { fn benchmark_algorithm(c: &mut Criterion, bench_name: &str, setup_fn: F, config: A) where F: Fn() -> Result>, - A: Algorithm + Clone, + A: Algorithm> + Clone, { c.bench_function(bench_name, |b| { b.iter_with_setup( diff --git a/examples/bimodal_ke_backend_compare.rs b/examples/bimodal_ke_backend_compare.rs index 0a5488aae..298d4ef3b 100644 --- a/examples/bimodal_ke_backend_compare.rs +++ b/examples/bimodal_ke_backend_compare.rs @@ -128,7 +128,7 @@ fn run_case Result<()> { .chain(NpodConfig::new().max_cycles(5))?; println!( - "Chained NPAG→NPOD: OBJF = {:.2}, {} support points, {} total cycles", - result.objf(), + "Chained NPAG→NPOD: -2LL = {:.2}, {} support points, {} total cycles", + result.n2ll(), result.get_theta().nspp(), result.cycles(), ); diff --git a/src/algorithms/diagnostics.rs b/src/algorithms/diagnostics.rs new file mode 100644 index 000000000..b064d8c0f --- /dev/null +++ b/src/algorithms/diagnostics.rs @@ -0,0 +1,186 @@ +//! Diagnostics shared by the non-parametric algorithms. +//! +//! These are free functions rather than methods on the algorithm state, so they +//! can be applied to any combination of equation, data and support points, and +//! so that logging helpers stay out of the algorithm trait. + +use anyhow::Result; +use ndarray::parallel::prelude::{IntoParallelIterator, ParallelIterator}; + +use pharmsol::prelude::data::{AssayErrorModels, Data}; +use pharmsol::prelude::simulator::Equation; +use pharmsol::{Predictions, Subject}; + +use crate::estimation::nonparametric::{Psi, Theta}; + +/// Identify subjects whose total probability given the model is zero or +/// non-finite. +/// +/// Each row of [`Psi`] holds the likelihood of a subject across every +/// support point, so a subject's probability is the sum across its row. A +/// subject is flagged when that sum is zero or not finite, meaning the model +/// cannot explain the subject's data. When any subject is flagged, detailed +/// per-subject diagnostics are logged and an error is returned. +pub fn check_zero_probability_subjects( + equation: &E, + data: &Data, + error_models: &AssayErrorModels, + theta: &Theta, + psi: &Psi, +) -> Result<()> { + let matrix = psi.matrix(); + + // Report non-finite entries; these propagate into the row sums below. + let nonfinite = matrix + .row_iter() + .flat_map(|row| row.iter().copied()) + .filter(|v| !v.is_finite()) + .count(); + if nonfinite > 0 { + tracing::warn!( + "Psi matrix contains {} non-finite value(s) of {} total", + nonfinite, + matrix.nrows() * matrix.ncols() + ); + } + + // A subject's probability is the sum across its row. + let subjects = data.subjects(); + let flagged: Vec = (0..matrix.nrows()) + .filter(|&i| { + let probability: f64 = (0..matrix.ncols()).map(|j| matrix[(i, j)]).sum(); + !probability.is_finite() || probability == 0.0 + }) + .collect(); + + if flagged.is_empty() { + return Ok(()); + } + + tracing::error!( + "{}/{} subjects have zero probability given the model", + flagged.len(), + matrix.nrows() + ); + + for &i in &flagged { + log_zero_probability_subject(equation, error_models, theta, subjects[i]); + } + + let ids: Vec<&String> = flagged.iter().map(|&i| subjects[i].id()).collect(); + Err(anyhow::anyhow!( + "The probability of {}/{} subjects is zero given the model. Affected subjects: {:?}", + flagged.len(), + matrix.nrows(), + ids + )) +} + +/// Log detailed likelihood diagnostics for a single subject whose +/// probability given the model is zero or non-finite. +/// +/// This is best-effort: support points that fail to simulate are reported at +/// debug level and skipped, since the caller is already handling an error. +fn log_zero_probability_subject( + equation: &E, + error_models: &AssayErrorModels, + theta: &Theta, + subject: &Subject, +) { + tracing::debug!("Subject with zero probability: {}", subject.id()); + + // Simulate every support point for this subject in parallel. + let mut results: Vec<_> = theta + .matrix() + .row_iter() + .enumerate() + .collect::>() + .into_par_iter() + .filter_map(|(i, spp)| { + let support_point: Vec = spp.iter().copied().collect(); + match equation.simulate_subject_dense(subject, &support_point, Some(error_models)) { + Ok((pred, ll)) => Some((i, support_point, pred.get_predictions(), ll)), + Err(err) => { + tracing::debug!("Support point #{} could not be simulated: {}", i, err); + None + } + } + }) + .collect(); + + if results.is_empty() { + tracing::debug!("\tNo support point could be simulated for this subject"); + return; + } + + // Summarise the distribution of likelihood values. + let mut nan = 0; + let mut pos_inf = 0; + let mut neg_inf = 0; + let mut zero = 0; + let mut valid = 0; + for (_, _, _, ll) in &results { + match ll { + Some(v) if v.is_nan() => nan += 1, + Some(v) if v.is_infinite() && v.is_sign_positive() => pos_inf += 1, + Some(v) if v.is_infinite() => neg_inf += 1, + Some(v) if *v == 0.0 => zero += 1, + Some(_) => valid += 1, + None => nan += 1, + } + } + + let total = results.len(); + let pct = |n: usize| 100.0 * n as f64 / total as f64; + tracing::debug!( + "\tLikelihood analysis for subject {} ({} support points):", + subject.id(), + total + ); + tracing::debug!("\tNaN likelihoods: {} ({:.1}%)", nan, pct(nan)); + tracing::debug!("\t+Inf likelihoods: {} ({:.1}%)", pos_inf, pct(pos_inf)); + tracing::debug!("\t-Inf likelihoods: {} ({:.1}%)", neg_inf, pct(neg_inf)); + tracing::debug!("\tZero likelihoods: {} ({:.1}%)", zero, pct(zero)); + tracing::debug!("\tValid likelihoods: {} ({:.1}%)", valid, pct(valid)); + + // Show the most likely support points to aid debugging. + results.sort_by(|a, b| { + b.3.unwrap_or(f64::NEG_INFINITY) + .partial_cmp(&a.3.unwrap_or(f64::NEG_INFINITY)) + .unwrap_or(std::cmp::Ordering::Equal) + }); + + const TAKE: usize = 3; + tracing::debug!("Top {} most likely support points:", TAKE); + for (i, support_point, preds, ll) in results.iter().take(TAKE) { + tracing::debug!("\tSupport point #{}: {:?}", i, support_point); + tracing::debug!("\t\tLog-likelihood: {:?}", ll); + tracing::debug!( + "\t\tTimes: {:?}", + preds.iter().map(|x| x.time()).collect::>() + ); + tracing::debug!( + "\t\tObservations: {:?}", + preds + .iter() + .map(|x| x.observation()) + .collect::>>() + ); + tracing::debug!( + "\t\tPredictions: {:?}", + preds.iter().map(|x| x.prediction()).collect::>() + ); + tracing::debug!( + "\t\tOuteqs: {:?}", + preds.iter().map(|x| x.outeq()).collect::>() + ); + tracing::debug!( + "\t\tStates: {:?}", + preds + .iter() + .map(|x| x.state().to_vec()) + .collect::>>() + ); + } + tracing::debug!("====================="); +} diff --git a/src/algorithms/mod.rs b/src/algorithms/mod.rs index 30a6c8d3d..31166e8e2 100644 --- a/src/algorithms/mod.rs +++ b/src/algorithms/mod.rs @@ -1,272 +1,177 @@ -use std::fs; -use std::path::Path; - -use crate::estimation::nonparametric::{NonParametricResult, Psi, Theta}; +use crate::estimation::nonparametric::{NPCycle, Psi, Theta, Weights}; use crate::estimation::{EstimationProblem, Framework}; use crate::results::FitResult; -use anyhow::Context; use anyhow::Result; -use ndarray::parallel::prelude::{IntoParallelIterator, ParallelIterator}; -use pharmsol::prelude::{data::Data, simulator::Equation}; +use pharmsol::prelude::data::{AssayErrorModels, Data}; +use pharmsol::prelude::simulator::Equation; -use pharmsol::{Predictions, Subject}; use serde::{Deserialize, Serialize}; +pub mod diagnostics; +pub mod nonparametric; +pub mod parametric; +mod stop; + /// Defines an algorithm that can fit an [`EstimationProblem`] to produce a result. /// /// Implementors are the lightweight, user-facing configuration structs (e.g. -/// `NpagConfig`). The heavy, mutable execution state used while fitting is an -/// internal implementation detail. -pub trait Algorithm { +/// [`NpagConfig`](nonparametric::NpagConfig)). The heavy, mutable execution +/// state used while fitting is an internal implementation detail. +pub trait Algorithm

{ /// The specific result struct (e.g. `NonParametricResult`). type Output: FitResult; /// Consumes the configuration and the problem, runs the optimization to /// completion, and returns the strictly-typed result. - fn fit(self, problem: EstimationProblem) -> Result; + fn fit(self, problem: P) -> Result; } -// Module organization for algorithm types -pub mod nonparametric; -pub mod parametric; - impl EstimationProblem { /// Consumes the problem and an algorithm configuration, runs the fit to /// completion, and returns the result. pub fn fit_with(self, algorithm: A) -> Result where - A: Algorithm, + A: Algorithm, { algorithm.fit(self) } } -pub trait NonParametricRunner: Sync + Send + 'static { - /// Identify subjects whose total probability given the model is zero or - /// non-finite. - /// - /// Each row of [`Psi`] holds the likelihood of a subject across every - /// support point, so a subject's probability is the sum across its row. A - /// subject is flagged when that sum is zero or not finite, meaning the model - /// cannot explain the subject's data. When any subject is flagged, detailed - /// per-subject diagnostics are logged and an error is returned. - fn check_zero_probability_subjects(&self) -> Result<()> { - let psi = self.psi().matrix(); - - // Report non-finite entries; these propagate into the row sums below. - let nonfinite = psi - .row_iter() - .flat_map(|row| row.iter().copied()) - .filter(|v| !v.is_finite()) - .count(); - if nonfinite > 0 { - tracing::warn!( - "Psi matrix contains {} non-finite value(s) of {} total", - nonfinite, - psi.nrows() * psi.ncols() - ); - } - - // A subject's probability is the sum across its row. - let subjects = self.data().subjects(); - let flagged: Vec = (0..psi.nrows()) - .filter(|&i| { - let probability: f64 = (0..psi.ncols()).map(|j| psi[(i, j)]).sum(); - !probability.is_finite() || probability == 0.0 - }) - .collect(); - - if flagged.is_empty() { - return Ok(()); - } - - tracing::error!( - "{}/{} subjects have zero probability given the model", - flagged.len(), - psi.nrows() - ); - - for &i in &flagged { - self.log_zero_probability_subject(subjects[i]); - } - - let ids: Vec<&String> = flagged.iter().map(|&i| subjects[i].id()).collect(); - Err(anyhow::anyhow!( - "The probability of {}/{} subjects is zero given the model. Affected subjects: {:?}", - flagged.len(), - psi.nrows(), - ids - )) - } - - /// Log detailed likelihood diagnostics for a single subject whose - /// probability given the model is zero or non-finite. - fn log_zero_probability_subject(&self, subject: &Subject) { - tracing::debug!("Subject with zero probability: {}", subject.id()); - - let error_model = self.error_models().clone(); - - // Simulate every support point for this subject in parallel. - let mut results: Vec<_> = self - .theta() - .matrix() - .row_iter() - .enumerate() - .collect::>() - .into_par_iter() - .map(|(i, spp)| { - let support_point: Vec = spp.iter().copied().collect(); - let (pred, ll) = self - .equation() - .simulate_subject_dense(subject, &support_point, Some(&error_model)) - .unwrap(); //TODO: Handle error - (i, support_point, pred.get_predictions(), ll) - }) - .collect(); - - // Summarise the distribution of likelihood values. - let mut nan = 0; - let mut pos_inf = 0; - let mut neg_inf = 0; - let mut zero = 0; - let mut valid = 0; - for (_, _, _, ll) in &results { - match ll { - Some(v) if v.is_nan() => nan += 1, - Some(v) if v.is_infinite() && v.is_sign_positive() => pos_inf += 1, - Some(v) if v.is_infinite() => neg_inf += 1, - Some(v) if *v == 0.0 => zero += 1, - Some(_) => valid += 1, - None => nan += 1, - } - } - - let total = results.len(); - let pct = |n: usize| 100.0 * n as f64 / total as f64; - tracing::debug!( - "\tLikelihood analysis for subject {} ({} support points):", - subject.id(), - total - ); - tracing::debug!("\tNaN likelihoods: {} ({:.1}%)", nan, pct(nan)); - tracing::debug!("\t+Inf likelihoods: {} ({:.1}%)", pos_inf, pct(pos_inf)); - tracing::debug!("\t-Inf likelihoods: {} ({:.1}%)", neg_inf, pct(neg_inf)); - tracing::debug!("\tZero likelihoods: {} ({:.1}%)", zero, pct(zero)); - tracing::debug!("\tValid likelihoods: {} ({:.1}%)", valid, pct(valid)); - - // Show the most likely support points to aid debugging. - results.sort_by(|a, b| { - b.3.unwrap_or(f64::NEG_INFINITY) - .partial_cmp(&a.3.unwrap_or(f64::NEG_INFINITY)) - .unwrap_or(std::cmp::Ordering::Equal) - }); - - const TAKE: usize = 3; - tracing::debug!("Top {} most likely support points:", TAKE); - for (i, support_point, preds, ll) in results.iter().take(TAKE) { - tracing::debug!("\tSupport point #{}: {:?}", i, support_point); - tracing::debug!("\t\tLog-likelihood: {:?}", ll); - tracing::debug!( - "\t\tTimes: {:?}", - preds.iter().map(|x| x.time()).collect::>() - ); - tracing::debug!( - "\t\tObservations: {:?}", - preds - .iter() - .map(|x| x.observation()) - .collect::>>() - ); - tracing::debug!( - "\t\tPredictions: {:?}", - preds.iter().map(|x| x.prediction()).collect::>() - ); - tracing::debug!( - "\t\tOuteqs: {:?}", - preds.iter().map(|x| x.outeq()).collect::>() - ); - tracing::debug!( - "\t\tStates: {:?}", - preds - .iter() - .map(|x| x.state().to_vec()) - .collect::>>() - ); - } - tracing::debug!("====================="); - } - - fn error_models(&self) -> &pharmsol::prelude::data::AssayErrorModels; +/// Read-only view of the state of a running fit. +/// +/// A runner exposes this state between cycles. Log-likelihoods are reported on +/// their natural scale: [`log_likelihood`](Self::log_likelihood) returns the +/// log-likelihood of the data under the current model (higher is better), and +/// [`n2ll`](Self::n2ll) the objective function actually minimized, +/// `-2 × log_likelihood` (lower is better). Runners that converge on a +/// different statistic (NCNPAG's marginal log-likelihood) report that statistic +/// instead, so `n2ll` is always the quantity being optimized. +pub trait FitState { /// Get the equation used in the algorithm fn equation(&self) -> &E; /// Get the data used in the algorithm fn data(&self) -> &Data; - - /// Increment the cycle counter and return the new value - fn increment_cycle(&mut self) -> usize; - /// Get the current cycle number - fn cycle(&self) -> usize; - /// Set the current [Theta] - fn set_theta(&mut self, theta: Theta); + /// Get the error models used in the algorithm + fn error_models(&self) -> &AssayErrorModels; /// Get the current [Theta] fn theta(&self) -> Θ /// Get the current [Psi] fn psi(&self) -> Ψ - /// Get the current likelihood - fn likelihood(&self) -> f64; - /// Get the current negative two log-likelihood - fn n2ll(&self) -> f64 { - -2.0 * self.likelihood() - } + /// Get the current support point weights + fn weights(&self) -> &Weights; + /// Get the current cycle number + fn cycle(&self) -> usize; /// Get the current [Status] of the algorithm fn status(&self) -> &Status; + /// Get the current log-likelihood (higher is better) + fn log_likelihood(&self) -> f64; + /// Get the current objective function, `-2 × log_likelihood` (lower is better) + fn n2ll(&self) -> f64 { + -2.0 * self.log_likelihood() + } +} + +/// The mutable state of a non-parametric fit, advanced one cycle at a time. +/// +/// Most algorithms run the same cycle: estimate the likelihood of every support +/// point, condense the support points, optimize the error models, then evaluate +/// the convergence criteria. Stages that don't apply default to no-ops, so an +/// algorithm only implements what it needs. Single-pass algorithms (NPMAP, +/// NCNPAG) implement `estimation` and `evaluation`, and report +/// [`StopReason::Converged`] once their single pass is done. +/// +/// Only [`next_cycle`](Self::next_cycle) writes the [`Status`], so the status an +/// observer sees always matches the state of the fit. +pub trait NonParametricRunner: FitState + Send + 'static { + /// The result produced when the fit finishes. + type Output: FitResult; + /// Set the current [Status] of the algorithm fn set_status(&mut self, status: Status); - /// Evaluate convergence criteria and update status + /// Increment the cycle counter and return the new value + fn increment_cycle(&mut self) -> usize; + /// Evaluate the convergence criteria and return the resulting [Status] + /// + /// Implementations must not store the status they return: + /// [`next_cycle`](Self::next_cycle) is the single writer. fn evaluation(&mut self) -> Result; - - /// Create and log a cycle state with the current algorithm state - fn log_cycle_state(&mut self); - - /// Initialize the algorithm, setting up initial [Theta] and [Status] - fn initialize(&mut self) -> Result<()> { - // If a stop file exists in the current directory, remove it - if Path::new("stop").exists() { - tracing::info!("Removing existing stop file prior to run"); - fs::remove_file("stop").context("Unable to remove previous stop file")?; - } - self.set_status(Status::Continue); - - Ok(()) - } + /// Estimate the likelihood of every support point for every subject fn estimation(&mut self) -> Result<()>; + /// Consume the runner and build the final result + #[allow(clippy::wrong_self_convention)] + fn into_result(self: Box) -> Result; + /// Append a cycle to the cycle log + fn push_cycle(&mut self, cycle: NPCycle); + /// Performs condensation of [Theta] and updates [Psi] /// /// This step reduces the number of support points in [Theta] based on the current weights, /// and updates the [Psi] matrix accordingly to reflect the new set of support points. /// It is typically performed after the estimation step in each cycle of the algorithm. - fn condensation(&mut self) -> Result<()>; - + fn condensation(&mut self) -> Result<()> { + Ok(()) + } /// Performs optimizations on the current `AssayErrorModels` and updates [Psi] accordingly /// /// This step refines the error model parameters to better fit the data, /// and subsequently updates the [Psi] matrix to reflect these changes. - fn optimizations(&mut self) -> Result<()>; - + fn optimizations(&mut self) -> Result<()> { + Ok(()) + } /// Performs expansion of [Theta] /// /// This step increases the number of support points in [Theta] based on the current distribution, /// allowing for exploration of the parameter space. - fn expansion(&mut self) -> Result<()>; + fn expansion(&mut self) -> Result<()> { + Ok(()) + } + /// Change in the objective function since the last logged cycle + /// + /// Called once per logged cycle, before the cycle is pushed to the log. + /// Implementations may use this to update their own bookkeeping. + fn objective_delta(&mut self) -> f64 { + 0.0 + } + + /// Initialize the algorithm, setting up initial [Theta] and [Status] + fn initialize(&mut self) -> Result<()> { + // If a stop file exists in the current directory, remove it + crate::algorithms::stop::remove_stop_file()?; + self.set_status(Status::Continue); + + Ok(()) + } + + /// Create and log a cycle state with the current algorithm state + fn log_cycle_state(&mut self) { + let delta = self.objective_delta(); + let state = NPCycle::new( + self.cycle(), + self.n2ll(), + self.error_models().clone(), + self.theta().clone(), + self.weights().clone(), + delta, + self.status().clone(), + ); + self.push_cycle(state); + } /// Proceed to the next cycle of the algorithm /// /// This method increments the cycle counter, performs expansion if necessary, /// and then runs the estimation, condensation, optimization, logging, and evaluation steps - /// in sequence. It returns the current [Status] of the algorithm after completing these steps. + /// in sequence. It returns the [Status] of the algorithm after completing these steps. + /// + /// A stop is final: once the runner has stopped, further cycles return the + /// same [Status] without touching the fit state. fn next_cycle(&mut self) -> Result { + if self.status().is_stop() { + return Ok(self.status().clone()); + } + let cycle = self.increment_cycle(); if cycle > 1 { @@ -278,7 +183,10 @@ pub trait NonParametricRunner: Sync + Send + 'stat self.estimation()?; self.condensation()?; self.optimizations()?; - self.evaluation() + let status = self.evaluation()?; + self.set_status(status.clone()); + self.log_cycle_state(); + Ok(status) } /// Fit the model until convergence or stopping criteria are met @@ -286,14 +194,12 @@ pub trait NonParametricRunner: Sync + Send + 'stat /// This method runs the full fitting process, starting with initialization, /// followed by iterative cycles of estimation, condensation, optimization, and evaluation /// until the algorithm converges or meets a stopping criteria. - fn fit(&mut self) -> Result> { - self.initialize()?; - while let Status::Continue = self.next_cycle()? {} - self.into_result() + fn fit(self: Box) -> Result { + let mut runner = self; + runner.initialize()?; + while runner.next_cycle()?.is_continue() {} + runner.into_result() } - - #[allow(clippy::wrong_self_convention)] - fn into_result(&self) -> Result>; } /// Where a fit stands: still running, or stopped (and why). diff --git a/src/algorithms/nonparametric/controller.rs b/src/algorithms/nonparametric/controller.rs index 22a714290..24e8342ab 100644 --- a/src/algorithms/nonparametric/controller.rs +++ b/src/algorithms/nonparametric/controller.rs @@ -25,7 +25,7 @@ use super::NonParametricAlgorithm; /// The accessors read the current state between steps, and [`finish`](Self::finish) or /// [`into_result`](Self::into_result) turn it into a result. pub struct FitController { - runner: Box>, + runner: Box>>, } impl FitController { @@ -49,8 +49,9 @@ impl FitController { /// Flag the fit as [`Aborted`](StopReason::Aborted) so the result shows an early stop. /// - /// This only sets the status; it won't halt a later [`step`](Self::step), which runs a full - /// cycle and overwrites it. Call it just before [`into_result`](Self::into_result). + /// A stop is final: the status is remembered, so the fit stays stopped even if you keep + /// calling [`step`](Self::step) — later steps return [`Stop`](Status::Stop) without + /// advancing the fit. pub fn request_stop(&mut self) { self.runner.set_status(Status::Stop(StopReason::Aborted)); } @@ -75,12 +76,13 @@ impl FitController { self.runner.psi() } - /// The current log-likelihood. - pub fn likelihood(&self) -> f64 { - self.runner.likelihood() + /// The current log-likelihood (higher is better). + pub fn log_likelihood(&self) -> f64 { + self.runner.log_likelihood() } - /// The current negative two log-likelihood (objective function). + /// The current negative two log-likelihood, the objective function that is minimized + /// (lower is better). pub fn n2ll(&self) -> f64 { self.runner.n2ll() } @@ -138,8 +140,8 @@ where impl EstimationProblem { /// Start a fit you drive yourself, one [`step`](FitController::step) at a time. /// - /// Takes any non-parametric config (`NpagConfig`, `NpodConfig`, `NpmapConfig`) or a - /// [`NonParametricAlgorithm`]. + /// Takes any non-parametric config (`NpagConfig`, `NpodConfig`, `NpmapConfig`, + /// `NcnpagConfig`) or a [`NonParametricAlgorithm`]. /// /// ```no_run /// use pmcore::prelude::*; diff --git a/src/algorithms/nonparametric/error_optim.rs b/src/algorithms/nonparametric/error_optim.rs index 33e8a5b26..495c43b51 100644 --- a/src/algorithms/nonparametric/error_optim.rs +++ b/src/algorithms/nonparametric/error_optim.rs @@ -63,9 +63,9 @@ impl ErrorOptimConfig { /// /// For each optimizable error model, the factor is perturbed up and down, /// `psi` is recomputed for each perturbation, and the IPM ([`burke`]) is run. -/// The direction with the highest objective function that also improves on the -/// current `objf` is adopted, updating `error_models`, `objf`, `lambda`, and -/// `psi` in place. The per-output adaptive step in `gamma_delta` is grown on +/// The direction with the highest log-likelihood that also improves on the +/// current one is adopted, updating `error_models`, `log_likelihood`, `lambda`, +/// and `psi` in place. The per-output adaptive step in `gamma_delta` is grown on /// improvement and shrunk every pass, resetting once it falls to `config.min_step`. /// /// A failure of the IPM for one perturbation direction is treated as a warning @@ -78,7 +78,7 @@ pub(crate) fn optimize_error_models( theta: &Theta, error_models: &mut AssayErrorModels, gamma_delta: &mut [f64], - objf: &mut f64, + log_likelihood: &mut f64, lambda: &mut Weights, psi: &mut Psi, config: &ErrorOptimConfig, @@ -112,7 +112,7 @@ pub(crate) fn optimize_error_models( // likelihoods as warnings, and continue with the other direction. If both // directions fail, we will not update the error model factor. let up = match burke(&psi_up) { - Ok((lambda, objf)) => Some((lambda, objf)), + Ok((lambda, ll)) => Some((lambda, ll)), Err(err) => { tracing::warn!( "Error in IPM during optim (up) for outeq {}: {:?}", @@ -123,7 +123,7 @@ pub(crate) fn optimize_error_models( } }; let down = match burke(&psi_down) { - Ok((lambda, objf)) => Some((lambda, objf)), + Ok((lambda, ll)) => Some((lambda, ll)), Err(err) => { tracing::warn!( "Error in IPM during optim (down) for outeq {}: {:?}", @@ -135,23 +135,23 @@ pub(crate) fn optimize_error_models( }; // Select the best improving candidate (if any) over the current - // objective. Among the two directions, the one with the higher - // objective function wins. + // log-likelihood. Among the two directions, the one with the higher + // log-likelihood wins. let mut best: Option<(f64, Weights, Psi, f64)> = None; - if let Some((lambda_up, objf_up)) = up { - if objf_up > *objf { - best = Some((objf_up, lambda_up, psi_up, gamma_up)); + if let Some((lambda_up, ll_up)) = up { + if ll_up > *log_likelihood { + best = Some((ll_up, lambda_up, psi_up, gamma_up)); } } - if let Some((lambda_down, objf_down)) = down { - let threshold = best.as_ref().map_or(*objf, |(o, ..)| *o); - if objf_down > threshold { - best = Some((objf_down, lambda_down, psi_down, gamma_down)); + if let Some((lambda_down, ll_down)) = down { + let threshold = best.as_ref().map_or(*log_likelihood, |(o, ..)| *o); + if ll_down > threshold { + best = Some((ll_down, lambda_down, psi_down, gamma_down)); } } - if let Some((new_objf, new_lambda, new_psi, gamma)) = best { + if let Some((new_log_likelihood, new_lambda, new_psi, gamma)) = best { error_models.set_factor(outeq, gamma)?; - *objf = new_objf; + *log_likelihood = new_log_likelihood; gamma_delta[outeq] *= config.growth; *lambda = new_lambda; *psi = new_psi; diff --git a/src/algorithms/nonparametric/mod.rs b/src/algorithms/nonparametric/mod.rs index ef82222f9..1945ce191 100644 --- a/src/algorithms/nonparametric/mod.rs +++ b/src/algorithms/nonparametric/mod.rs @@ -9,6 +9,7 @@ //! - [`NPAG`](npag): Non-Parametric Adaptive Grid //! - [`NPOD`](npod): Non-Parametric Optimal Design //! - [`NPMAP`](npmap): Maximum a posteriori reweighting +//! - [`NCNPAG`](ncnpag): Non-collapsing NPAG //! //! # Algorithm Selection //! @@ -27,6 +28,9 @@ pub mod npag; pub mod npmap; pub mod npod; +// Cycle stages shared by the non-parametric algorithms +mod stages; + // Incremental, observable fitting (stepping handle + per-cycle observers). pub mod controller; @@ -75,9 +79,9 @@ use pharmsol::prelude::simulator::Equation; /// // `problem.fit_with(config)` accepts the config directly. /// ``` /// -/// Each configuration type ([`NpagConfig`], [`NpodConfig`], [`NpmapConfig`]) implements -/// [`Algorithm`] by delegating to the matching enum variant, so configs can be passed to -/// `fit_with` without converting them first. +/// Each configuration type ([`NpagConfig`], [`NpodConfig`], [`NpmapConfig`], +/// [`NcnpagConfig`]) converts into this enum and therefore implements [`Algorithm`], so +/// configs can be passed to `fit_with` without converting them first. #[derive(Debug, Clone)] pub enum NonParametricAlgorithm { /// Non-Parametric Adaptive Grid. @@ -150,10 +154,10 @@ impl NonParametricAlgorithm { pub(crate) fn into_runner( self, problem: EstimationProblem, - ) -> Result>> { + ) -> Result>>> { // `problem.prior` is the prior `Theta` (which also carries the parameter // space) and `problem.error_models` is strictly `AssayErrorModels`. - let runner: Box> = match self { + let runner: Box>> = match self { Self::Npag(config) => Box::new(NPAG::from_parts( problem.model.equation, problem.data, @@ -187,46 +191,19 @@ impl NonParametricAlgorithm { } } -impl Algorithm for NonParametricAlgorithm { +/// Every configuration that converts into a [`NonParametricAlgorithm`] is itself an +/// [`Algorithm`], so the algorithm-specific config structs can be passed straight to +/// [`fit_with`](crate::estimation::EstimationProblem::fit_with) without converting them +/// first. +impl Algorithm> for A +where + E: Equation + Send + 'static, + A: Into, +{ type Output = NonParametricResult; fn fit(self, problem: EstimationProblem) -> Result { - let mut runner = self.into_runner(problem)?; + let runner = self.into().into_runner(problem)?; runner.fit() } } - -// Each configuration struct delegates to its matching `NonParametricAlgorithm` variant so it -// can be passed directly to `fit_with`. This keeps the variant-specific setters on the config -// types (compile-time checked) while the enum remains the single source of fitting logic. -impl Algorithm for NpagConfig { - type Output = NonParametricResult; - - fn fit(self, problem: EstimationProblem) -> Result { - NonParametricAlgorithm::from(self).fit(problem) - } -} - -impl Algorithm for NpodConfig { - type Output = NonParametricResult; - - fn fit(self, problem: EstimationProblem) -> Result { - NonParametricAlgorithm::from(self).fit(problem) - } -} - -impl Algorithm for NpmapConfig { - type Output = NonParametricResult; - - fn fit(self, problem: EstimationProblem) -> Result { - NonParametricAlgorithm::from(self).fit(problem) - } -} - -impl Algorithm for NcnpagConfig { - type Output = NonParametricResult; - - fn fit(self, problem: EstimationProblem) -> Result { - NonParametricAlgorithm::from(self).fit(problem) - } -} diff --git a/src/algorithms/nonparametric/ncnpag.rs b/src/algorithms/nonparametric/ncnpag.rs index 860968ce7..ee069d6bd 100644 --- a/src/algorithms/nonparametric/ncnpag.rs +++ b/src/algorithms/nonparametric/ncnpag.rs @@ -1,7 +1,7 @@ use crate::{ algorithms::{ nonparametric::{npag::NPAG, NpagConfig}, - NonParametricRunner, Status, StopReason, + FitState, NonParametricRunner, Status, StopReason, }, estimation::nonparametric::{ calculate_psi, CycleLog, NPCycle, NonParametricResult, Psi, Theta, Weights, @@ -78,14 +78,14 @@ pub struct NCNPAG { psi: Psi, theta: Theta, w: Weights, - objf: f64, + log_likelihood: f64, cycle: usize, status: Status, data: Data, cyclelog: CycleLog, error_models: AssayErrorModels, prior: Theta, - cycles: usize, + refinement_cycles: usize, progress: bool, } @@ -102,14 +102,14 @@ impl NCNPAG { psi: Psi::new(), theta: theta.clone(), w: Weights::default(), - objf: f64::INFINITY, + log_likelihood: f64::INFINITY, cycle: 0, status: Status::Continue, data, cyclelog: CycleLog::new(), error_models, prior: theta, - cycles: config.cycles, + refinement_cycles: config.cycles, progress: config.progress, }) } @@ -211,27 +211,7 @@ fn marginal_loglik(psi: &Psi, w: &Weights) -> f64 { + psi.log_scale() } -impl NonParametricRunner for NCNPAG { - fn into_result(&self) -> Result> { - NonParametricResult::new( - self.equation.clone(), - self.data.clone(), - self.error_models.clone(), - self.prior.clone(), - self.theta.clone(), - self.psi.clone(), - self.w.clone(), - self.objf, - self.cycle, - self.status.clone(), - self.cyclelog.clone(), - ) - } - - fn error_models(&self) -> &AssayErrorModels { - &self.error_models - } - +impl FitState for NCNPAG { fn equation(&self) -> &E { &self.equation } @@ -240,41 +220,74 @@ impl NonParametricRunner for NCNPAG { &self.data } - fn likelihood(&self) -> f64 { - self.objf + fn error_models(&self) -> &AssayErrorModels { + &self.error_models } - fn increment_cycle(&mut self) -> usize { - 0 + fn theta(&self) -> &Theta { + &self.theta } - fn cycle(&self) -> usize { - 0 + fn psi(&self) -> &Psi { + &self.psi } - fn set_theta(&mut self, theta: Theta) { - self.theta = theta; + fn weights(&self) -> &Weights { + &self.w } - fn theta(&self) -> &Theta { - &self.theta + fn cycle(&self) -> usize { + self.cycle } - fn psi(&self) -> &Psi { - &self.psi + fn status(&self) -> &Status { + &self.status + } + + fn log_likelihood(&self) -> f64 { + self.log_likelihood } +} + +impl NonParametricRunner for NCNPAG { + type Output = NonParametricResult; fn set_status(&mut self, status: Status) { self.status = status; } - fn status(&self) -> &Status { - &self.status + fn increment_cycle(&mut self) -> usize { + self.cycle += 1; + self.cycle + } + + fn into_result(self: Box) -> Result { + let this = *self; + let n2ll = this.n2ll(); + + NonParametricResult::new( + this.equation, + this.data, + this.error_models, + this.prior, + this.theta, + this.psi, + this.w, + n2ll, + this.cycle, + this.status, + this.cyclelog, + ) } + fn push_cycle(&mut self, cycle: NPCycle) { + self.cyclelog.push(cycle); + } + + /// NCNPAG is a single-pass reweighting of the prior support points, so the + /// only cycle stops as soon as it has been evaluated. fn evaluation(&mut self) -> Result { - self.status = Status::Stop(StopReason::Converged); - Ok(self.status.clone()) + Ok(Status::Stop(StopReason::Converged)) } fn estimation(&mut self) -> Result<()> { @@ -349,14 +362,14 @@ impl NonParametricRunner for NCNPAG { self.w = Weights::from_vec(kept.iter().map(|w| w / sum).collect()); // NPAGFULL: refine each surviving point with a full NPAG seeded from it. - if self.cycles > 0 { + if self.refinement_cycles > 0 { let (refined_theta, refined_weights) = refine_points( &self.equation, &self.data, &self.error_models, &self.theta, &self.w, - self.cycles, + self.refinement_cycles, self.progress, )?; self.theta = refined_theta; @@ -372,45 +385,9 @@ impl NonParametricRunner for NCNPAG { false, )?; - self.objf = marginal_loglik(&self.psi, &self.w); + self.log_likelihood = marginal_loglik(&self.psi, &self.w); Ok(()) } - - fn condensation(&mut self) -> Result<()> { - Ok(()) - } - - fn optimizations(&mut self) -> Result<()> { - Ok(()) - } - - fn expansion(&mut self) -> Result<()> { - Ok(()) - } - - fn log_cycle_state(&mut self) { - let state = NPCycle::new( - self.cycle, - self.objf, - self.error_models.clone(), - self.theta.clone(), - self.w.clone(), - self.theta.nspp(), - 0.0, - self.status.clone(), - ); - self.cyclelog.push(state); - } - - /// NCNPAG is a single-pass reweighting: it evaluates the likelihood of the - /// fixed prior support points once, rather than iterating cycles. - fn fit(&mut self) -> Result> { - self.estimation()?; - self.evaluation()?; - self.log_cycle_state(); - - self.into_result() - } } #[cfg(test)] diff --git a/src/algorithms/nonparametric/npag.rs b/src/algorithms/nonparametric/npag.rs index 097b9e2d9..9e554e77c 100644 --- a/src/algorithms/nonparametric/npag.rs +++ b/src/algorithms/nonparametric/npag.rs @@ -1,23 +1,18 @@ -use crate::algorithms::{NonParametricRunner, Status, StopReason}; +use crate::algorithms::{FitState, NonParametricRunner, Status, StopReason}; use crate::estimation::nonparametric::{ - calculate_psi, CycleLog, NPCycle, NonParametricResult, Psi, Theta, Weights, + CycleLog, NPCycle, NonParametricResult, Psi, Theta, Weights, }; -pub(crate) use crate::estimation::nonparametric::ipm::burke; -pub(crate) use crate::estimation::nonparametric::qr; - -use anyhow::bail; use anyhow::Result; use pharmsol::prelude::{ data::{AssayErrorModels, Data}, simulator::Equation, }; -use pharmsol::prelude::AssayErrorModel; - use crate::estimation::nonparametric::adaptative_grid; use super::error_optim::{optimize_error_models, ErrorOptimConfig}; +use super::stages::{self, CondensationOptions}; use serde::{Deserialize, Serialize}; @@ -119,8 +114,8 @@ pub struct NPAG { lambda: Weights, w: Weights, eps: f64, - last_objf: f64, - objf: f64, + last_log_likelihood: f64, + log_likelihood: f64, f0: f64, f1: f64, cycle: usize, @@ -136,8 +131,8 @@ impl NPAG { /// Construct an `NPAG` instance from explicit parts. /// /// The `parameter_space` is used solely to derive the finite bounds for the - /// adaptive grid. Initial support points can be supplied separately via - /// [`NonParametricRunner::set_theta`]. + /// adaptive grid. The initial support points come from `theta`, which is + /// also kept as the prior. pub(crate) fn from_parts( equation: E, data: Data, @@ -158,8 +153,8 @@ impl NPAG { lambda: Weights::default(), w: Weights::default(), eps, - last_objf: -1e30, - objf: f64::NEG_INFINITY, + last_log_likelihood: -1e30, + log_likelihood: f64::NEG_INFINITY, f0: -1e30, f1: f64::default(), cycle: 0, @@ -173,89 +168,97 @@ impl NPAG { } } -impl NonParametricRunner for NPAG { +impl FitState for NPAG { fn equation(&self) -> &E { &self.equation } - fn into_result(&self) -> Result> { - NonParametricResult::new( - self.equation.clone(), - self.data.clone(), - self.error_models.clone(), - self.prior.clone(), - self.theta.clone(), - self.psi.clone(), - self.w.clone(), - -2. * self.objf, - self.cycle, - self.status.clone(), - self.cycle_log.clone(), - ) + fn data(&self) -> &Data { + &self.data } fn error_models(&self) -> &AssayErrorModels { &self.error_models } - fn data(&self) -> &Data { - &self.data + fn theta(&self) -> &Theta { + &self.theta } - fn likelihood(&self) -> f64 { - self.objf + fn psi(&self) -> &Psi { + &self.psi } - fn increment_cycle(&mut self) -> usize { - self.cycle += 1; - self.cycle + fn weights(&self) -> &Weights { + &self.w } fn cycle(&self) -> usize { self.cycle } - fn set_theta(&mut self, theta: Theta) { - self.theta = theta; + fn status(&self) -> &Status { + &self.status } - fn theta(&self) -> &Theta { - &self.theta + fn log_likelihood(&self) -> f64 { + self.log_likelihood + } +} + +impl NonParametricRunner for NPAG { + type Output = NonParametricResult; + + fn set_status(&mut self, status: Status) { + self.status = status; } - fn psi(&self) -> &Psi { - &self.psi + fn increment_cycle(&mut self) -> usize { + self.cycle += 1; + self.cycle } - fn evaluation(&mut self) -> Result { - tracing::info!("Objective function = {:.4}", -2.0 * self.objf); - tracing::debug!("Support points: {}", self.theta.nspp()); + fn into_result(self: Box) -> Result { + let this = *self; + let n2ll = this.n2ll(); - self.error_models.iter().for_each(|(outeq, em)| { - if AssayErrorModel::None == *em { - return; - } - tracing::debug!( - "Error model for outeq {}: {:.2}", - outeq, - em.factor().unwrap_or_default() - ); - }); + NonParametricResult::new( + this.equation, + this.data, + this.error_models, + this.prior, + this.theta, + this.psi, + this.w, + n2ll, + this.cycle, + this.status, + this.cycle_log, + ) + } + fn push_cycle(&mut self, cycle: NPCycle) { + self.cycle_log.push(cycle); + } + + fn objective_delta(&mut self) -> f64 { + let delta = (self.last_log_likelihood - self.log_likelihood).abs(); + self.last_log_likelihood = self.log_likelihood; + delta + } + + fn evaluation(&mut self) -> Result { + stages::log_objective( + self.log_likelihood, + self.last_log_likelihood, + &self.theta, + &self.error_models, + ); tracing::debug!("EPS = {:.4}", self.eps); - // Increasing objf signals instability or model misspecification. - if self.last_objf > self.objf + 1e-4 { - tracing::warn!( - "Objective function decreased from {:.4} to {:.4} (delta = {})", - -2.0 * self.last_objf, - -2.0 * self.objf, - -2.0 * self.last_objf - -2.0 * self.objf - ); - } let psi = self.psi.matrix(); let w = &self.w; - if (self.last_objf - self.objf).abs() <= self.config.objective_tolerance + if (self.last_log_likelihood - self.log_likelihood).abs() <= self.config.objective_tolerance && self.eps > self.config.min_eps { self.eps /= 2.; @@ -263,10 +266,8 @@ impl NonParametricRunner for NPAG { let pyl = psi * w.weights(); self.f1 = pyl.iter().map(|x| x.ln()).sum(); if (self.f1 - self.f0).abs() <= self.config.pyl_tolerance { - tracing::info!("The model converged after {} cycles", self.cycle,); - self.set_status(Status::Stop(StopReason::Converged)); - self.log_cycle_state(); - return Ok(self.status().clone()); + tracing::info!("The model converged after {} cycles", self.cycle); + return Ok(Status::Stop(StopReason::Converged)); } else { self.f0 = self.f1; self.eps = self.config.eps; @@ -277,111 +278,54 @@ impl NonParametricRunner for NPAG { // Stop if we have reached maximum number of cycles if self.cycle >= self.config.max_cycles { tracing::warn!("Maximum number of cycles reached"); - self.set_status(Status::Stop(StopReason::MaxCycles)); - self.log_cycle_state(); - return Ok(self.status().clone()); + return Ok(Status::Stop(StopReason::MaxCycles)); } // Stop if stopfile exists - if std::path::Path::new("stop").exists() { + if crate::algorithms::stop::stop_file_present() { tracing::warn!("Stopfile detected - breaking"); - self.set_status(Status::Stop(StopReason::StopFile)); - self.log_cycle_state(); - return Ok(self.status().clone()); + return Ok(Status::Stop(StopReason::StopFile)); } // Continue with normal operation - self.set_status(Status::Continue); - self.log_cycle_state(); - Ok(self.status().clone()) + Ok(Status::Continue) } fn estimation(&mut self) -> Result<()> { - self.psi = calculate_psi( + let (psi, lambda, log_likelihood) = stages::estimate( &self.equation, &self.data, - &self.theta, &self.error_models, + &self.theta, self.cycle == 1 && self.config.progress, )?; - if let Err(err) = self.check_zero_probability_subjects() { - bail!(err); - } + self.psi = psi; + self.lambda = lambda; + self.log_likelihood = log_likelihood; - (self.lambda, _) = match burke(&self.psi) { - Ok((lambda, objf)) => (lambda, objf), - Err(err) => { - bail!("Error in IPM during estimation: {:?}", err); - } - }; Ok(()) } fn condensation(&mut self) -> Result<()> { - // Filter out the support points with lambda < max(lambda)/1000 - - let max_lambda = self - .lambda - .iter() - .fold(f64::NEG_INFINITY, |acc, x| x.max(acc)); - - let mut keep = Vec::::new(); - for (index, lam) in self.lambda.iter().enumerate() { - if lam > max_lambda * self.config.prune_threshold { - keep.push(index); - } - } - if self.psi.matrix().ncols() != keep.len() { - tracing::debug!( - "Lambda (max/1000) dropped {} support point(s)", - self.psi.matrix().ncols() - keep.len(), - ); - } - - self.theta.filter_indices(keep.as_slice()); - self.psi.filter_column_indices(keep.as_slice()); - - //Rank-Revealing Factorization - let (r, perm) = qr::qrd(&self.psi)?; - - let mut keep = Vec::::new(); - - // The minimum between the number of subjects and the actual number of support points - let keep_n = self.psi.matrix().ncols().min(self.psi.matrix().nrows()); - for i in 0..keep_n { - let test = r.col(i).norm_l2(); - let r_diag_val = r.get(i, i); - let ratio = r_diag_val / test; - if ratio.abs() >= self.config.qr_tolerance { - keep.push(*perm.get(i).unwrap()); - } - } - - // If a support point is dropped, log it as a debug message - if self.psi.matrix().ncols() != keep.len() { - tracing::debug!( - "QR decomposition dropped {} support point(s)", - self.psi.matrix().ncols() - keep.len(), - ); - } + let (lambda, log_likelihood) = stages::condense( + &self.equation, + &self.data, + &self.error_models, + &mut self.theta, + &mut self.psi, + &self.lambda, + &CondensationOptions { + prune_threshold: self.config.prune_threshold, + qr_tolerance: self.config.qr_tolerance, + check_zero_probability: true, + }, + )?; - // Filter to keep only the support points (rows) that are in the `keep` vector - self.theta.filter_indices(keep.as_slice()); - // Filter to keep only the support points (columns) that are in the `keep` vector - self.psi.filter_column_indices(keep.as_slice()); - - self.check_zero_probability_subjects()?; - (self.lambda, self.objf) = match burke(&self.psi) { - Ok((lambda, objf)) => (lambda, objf), - Err(err) => { - return Err(anyhow::anyhow!( - "Error in IPM during condensation: {:?}", - err - )); - } - }; + self.lambda = lambda; + self.log_likelihood = log_likelihood; self.w = self.lambda.clone(); + Ok(()) } @@ -392,7 +336,7 @@ impl NonParametricRunner for NPAG { &self.theta, &mut self.error_models, &mut self.gamma_delta, - &mut self.objf, + &mut self.log_likelihood, &mut self.lambda, &mut self.psi, &self.config.error_optim, @@ -408,29 +352,6 @@ impl NonParametricRunner for NPAG { )?; Ok(()) } - - fn set_status(&mut self, status: Status) { - self.status = status; - } - - fn status(&self) -> &Status { - &self.status - } - - fn log_cycle_state(&mut self) { - let state = NPCycle::new( - self.cycle, - -2. * self.objf, - self.error_models.clone(), - self.theta.clone(), - self.w.clone(), - self.theta.nspp(), - (self.last_objf - self.objf).abs(), - self.status.clone(), - ); - self.cycle_log.push(state); - self.last_objf = self.objf; - } } #[cfg(test)] diff --git a/src/algorithms/nonparametric/npmap.rs b/src/algorithms/nonparametric/npmap.rs index c5765a6a1..2c28a7bbd 100644 --- a/src/algorithms/nonparametric/npmap.rs +++ b/src/algorithms/nonparametric/npmap.rs @@ -1,5 +1,5 @@ use crate::{ - algorithms::{NonParametricRunner, Status, StopReason}, + algorithms::{FitState, NonParametricRunner, Status, StopReason}, estimation::nonparametric::{ calculate_psi, CycleLog, NPCycle, NonParametricResult, Psi, Theta, Weights, }, @@ -67,27 +67,7 @@ impl NPMAP { } } -impl NonParametricRunner for NPMAP { - fn into_result(&self) -> Result> { - NonParametricResult::new( - self.equation.clone(), - self.data.clone(), - self.error_models.clone(), - self.prior.clone(), - self.theta.clone(), - self.psi.clone(), - self.w.clone(), - self.objf, - self.cycle, - self.status.clone(), - self.cyclelog.clone(), - ) - } - - fn error_models(&self) -> &AssayErrorModels { - &self.error_models - } - +impl FitState for NPMAP { fn equation(&self) -> &E { &self.equation } @@ -96,41 +76,74 @@ impl NonParametricRunner for NPMAP { &self.data } - fn likelihood(&self) -> f64 { - self.objf + fn error_models(&self) -> &AssayErrorModels { + &self.error_models } - fn increment_cycle(&mut self) -> usize { - 0 + fn theta(&self) -> &Theta { + &self.theta } - fn cycle(&self) -> usize { - 0 + fn psi(&self) -> &Psi { + &self.psi } - fn set_theta(&mut self, theta: Theta) { - self.theta = theta; + fn weights(&self) -> &Weights { + &self.w } - fn theta(&self) -> &Theta { - &self.theta + fn cycle(&self) -> usize { + self.cycle } - fn psi(&self) -> &Psi { - &self.psi + fn status(&self) -> &Status { + &self.status + } + + fn log_likelihood(&self) -> f64 { + self.objf } +} + +impl NonParametricRunner for NPMAP { + type Output = NonParametricResult; fn set_status(&mut self, status: Status) { self.status = status; } - fn status(&self) -> &Status { - &self.status + fn increment_cycle(&mut self) -> usize { + self.cycle += 1; + self.cycle + } + + fn into_result(self: Box) -> Result { + let this = *self; + let objf = this.n2ll(); + + NonParametricResult::new( + this.equation, + this.data, + this.error_models, + this.prior, + this.theta, + this.psi, + this.w, + objf, + this.cycle, + this.status, + this.cyclelog, + ) } + fn push_cycle(&mut self, cycle: NPCycle) { + self.cyclelog.push(cycle); + } + + /// NPMAP is a single-pass reweighting of the prior support points, so the + /// only cycle stops as soon as it has been evaluated. fn evaluation(&mut self) -> Result { - self.status = Status::Stop(StopReason::Converged); - Ok(self.status.clone()) + Ok(Status::Stop(StopReason::Converged)) } fn estimation(&mut self) -> Result<()> { @@ -144,41 +157,4 @@ impl NonParametricRunner for NPMAP { (self.w, self.objf) = burke(&self.psi).context("Error in IPM")?; Ok(()) } - - fn condensation(&mut self) -> Result<()> { - Ok(()) - } - - fn optimizations(&mut self) -> Result<()> { - Ok(()) - } - - fn expansion(&mut self) -> Result<()> { - Ok(()) - } - - fn log_cycle_state(&mut self) { - // NPMAP doesn't track last_objf, so we use 0.0 as the delta - let state = NPCycle::new( - self.cycle, - self.objf, - self.error_models.clone(), - self.theta.clone(), - self.w.clone(), - self.theta.nspp(), - 0.0, - self.status.clone(), - ); - self.cyclelog.push(state); - } - - /// NPMAP is a single-pass reweighting: it evaluates the likelihood of the - /// fixed prior support points once, rather than iterating cycles. - fn fit(&mut self) -> Result> { - self.estimation()?; - self.evaluation()?; - self.log_cycle_state(); - - self.into_result() - } } diff --git a/src/algorithms/nonparametric/npod.rs b/src/algorithms/nonparametric/npod.rs index 32aaa3372..6b6e5d8e7 100644 --- a/src/algorithms/nonparametric/npod.rs +++ b/src/algorithms/nonparametric/npod.rs @@ -1,23 +1,20 @@ use crate::{ - algorithms::{NonParametricRunner, Status, StopReason}, - estimation::nonparametric::{ - calculate_psi, ipm::burke, qr, CycleLog, NPCycle, NonParametricResult, Psi, Theta, Weights, - }, + algorithms::{FitState, NonParametricRunner, Status, StopReason}, + estimation::nonparametric::{CycleLog, NPCycle, NonParametricResult, Psi, Theta, Weights}, }; use pharmsol::ParameterOptimizer; -use anyhow::bail; use anyhow::Result; use pharmsol::prelude::{data::Data, simulator::Equation}; -use pharmsol::{prelude::AssayErrorModel, AssayErrorModels}; +use pharmsol::AssayErrorModels; use ndarray::Array1; use rayon::prelude::{IntoParallelRefMutIterator, ParallelIterator}; use serde::{Deserialize, Serialize}; use super::error_optim::{optimize_error_models, ErrorOptimConfig}; +use super::stages::{self, CondensationOptions}; -const THETA_F: f64 = 1e-2; const THETA_D: f64 = 1e-4; /// Configuration options for the Non-Parametric Optimal Design (NPOD) algorithm. @@ -29,6 +26,12 @@ pub struct NpodConfig { pub error_optim: ErrorOptimConfig, /// Whether to print progress information during the first cycle. pub progress: bool, + /// Convergence tolerance on the change in the objective function. + pub objective_tolerance: f64, + /// Support points whose weight is below `max_weight × prune_threshold` are dropped. + pub prune_threshold: f64, + /// Minimum `|r_ii| / ‖r_i‖` ratio for a support point to survive the QR decomposition. + pub qr_tolerance: f64, } impl NpodConfig { @@ -50,6 +53,26 @@ impl NpodConfig { self.progress = progress; self } + + /// Set the convergence tolerance on the change in the objective function. + pub fn objective_tolerance(mut self, tolerance: f64) -> Self { + self.objective_tolerance = tolerance; + self + } + + /// Set the weight below which support points are dropped, relative to the + /// largest weight. + pub fn prune_threshold(mut self, threshold: f64) -> Self { + self.prune_threshold = threshold; + self + } + + /// Set the minimum `|r_ii| / ‖r_i‖` ratio for a support point to survive the + /// QR decomposition. + pub fn qr_tolerance(mut self, tolerance: f64) -> Self { + self.qr_tolerance = tolerance; + self + } } impl Default for NpodConfig { @@ -58,6 +81,9 @@ impl Default for NpodConfig { max_cycles: 100, error_optim: ErrorOptimConfig::default(), progress: true, + objective_tolerance: 1e-2, + prune_threshold: 1e-3, + qr_tolerance: 1e-8, } } } @@ -70,12 +96,11 @@ pub struct NPOD { theta: Theta, lambda: Weights, w: Weights, - last_objf: f64, - objf: f64, + last_log_likelihood: f64, + log_likelihood: f64, cycle: usize, gamma_delta: Vec, error_models: AssayErrorModels, - converged: bool, status: Status, cycle_log: CycleLog, data: Data, @@ -99,12 +124,11 @@ impl NPOD { theta, lambda: Weights::default(), w: Weights::default(), - last_objf: -1e30, - objf: f64::NEG_INFINITY, + last_log_likelihood: -1e30, + log_likelihood: f64::NEG_INFINITY, cycle: 0, gamma_delta, error_models, - converged: false, status: Status::Continue, cycle_log: CycleLog::new(), data, @@ -113,210 +137,147 @@ impl NPOD { } } -impl NonParametricRunner for NPOD { - fn into_result(&self) -> Result> { - NonParametricResult::new( - self.equation.clone(), - self.data.clone(), - self.error_models.clone(), - self.prior.clone(), - self.theta.clone(), - self.psi.clone(), - self.w.clone(), - -2. * self.objf, - self.cycle, - self.status.clone(), - self.cycle_log.clone(), - ) - } - +impl FitState for NPOD { fn equation(&self) -> &E { &self.equation } + fn data(&self) -> &Data { + &self.data + } + fn error_models(&self) -> &AssayErrorModels { &self.error_models } - fn data(&self) -> &Data { - &self.data + fn theta(&self) -> &Theta { + &self.theta } - fn increment_cycle(&mut self) -> usize { - self.cycle += 1; - self.cycle + fn psi(&self) -> &Psi { + &self.psi + } + + fn weights(&self) -> &Weights { + &self.w } fn cycle(&self) -> usize { self.cycle } - fn set_theta(&mut self, theta: Theta) { - self.theta = theta; + fn status(&self) -> &Status { + &self.status } - fn theta(&self) -> &Theta { - &self.theta + fn log_likelihood(&self) -> f64 { + self.log_likelihood } +} - fn psi(&self) -> &Psi { - &self.psi +impl NonParametricRunner for NPOD { + type Output = NonParametricResult; + + fn set_status(&mut self, status: Status) { + self.status = status; } - fn likelihood(&self) -> f64 { - self.objf + fn increment_cycle(&mut self) -> usize { + self.cycle += 1; + self.cycle } - fn set_status(&mut self, status: Status) { - self.status = status; + fn into_result(self: Box) -> Result { + let this = *self; + let n2ll = this.n2ll(); + + NonParametricResult::new( + this.equation, + this.data, + this.error_models, + this.prior, + this.theta, + this.psi, + this.w, + n2ll, + this.cycle, + this.status, + this.cycle_log, + ) } - fn status(&self) -> &Status { - &self.status + fn push_cycle(&mut self, cycle: NPCycle) { + self.cycle_log.push(cycle); } - fn log_cycle_state(&mut self) { - let state = NPCycle::new( - self.cycle, - -2. * self.objf, - self.error_models.clone(), - self.theta.clone(), - self.w.clone(), - self.theta.nspp(), - (self.last_objf - self.objf).abs(), - self.status.clone(), - ); - self.cycle_log.push(state); - self.last_objf = self.objf; + fn objective_delta(&mut self) -> f64 { + let delta = (self.last_log_likelihood - self.log_likelihood).abs(); + self.last_log_likelihood = self.log_likelihood; + delta } fn evaluation(&mut self) -> Result { - tracing::info!("Objective function = {:.4}", -2.0 * self.objf); - tracing::debug!("Support points: {}", self.theta.nspp()); - self.error_models.iter().for_each(|(outeq, em)| { - if AssayErrorModel::None == *em { - return; - } - tracing::debug!( - "Error model for outeq {}: {:.16}", - outeq, - em.factor().unwrap_or_default() - ); - }); - if self.last_objf > self.objf + 1e-4 { - tracing::warn!( - "Objective function decreased from {:.4} to {:.4} (delta = {})", - -2.0 * self.last_objf, - -2.0 * self.objf, - -2.0 * self.last_objf - -2.0 * self.objf - ); - } + stages::log_objective( + self.log_likelihood, + self.last_log_likelihood, + &self.theta, + &self.error_models, + ); - if (self.last_objf - self.objf).abs() <= THETA_F { + if (self.last_log_likelihood - self.log_likelihood).abs() <= self.config.objective_tolerance + { tracing::info!("Objective function convergence reached"); - self.converged = true; - self.set_status(Status::Stop(StopReason::Converged)); - self.log_cycle_state(); - return Ok(self.status.clone()); + return Ok(Status::Stop(StopReason::Converged)); } if self.cycle >= self.config.max_cycles { tracing::warn!("Maximum number of cycles reached"); - self.converged = true; - self.set_status(Status::Stop(StopReason::MaxCycles)); - self.log_cycle_state(); - return Ok(self.status.clone()); + return Ok(Status::Stop(StopReason::MaxCycles)); } - if std::path::Path::new("stop").exists() { + if crate::algorithms::stop::stop_file_present() { tracing::warn!("Stopfile detected - breaking"); - self.converged = true; - self.set_status(Status::Stop(StopReason::StopFile)); - self.log_cycle_state(); - return Ok(self.status.clone()); + return Ok(Status::Stop(StopReason::StopFile)); } - self.status = Status::Continue; - self.log_cycle_state(); - Ok(self.status.clone()) + Ok(Status::Continue) } fn estimation(&mut self) -> Result<()> { - let error_model: AssayErrorModels = self.error_models.clone(); - - self.psi = calculate_psi( + let (psi, lambda, log_likelihood) = stages::estimate( &self.equation, &self.data, + &self.error_models, &self.theta, - &error_model, self.cycle == 1 && self.config.progress, )?; - if let Err(err) = self.check_zero_probability_subjects() { - bail!(err); - } + self.psi = psi; + self.lambda = lambda; + self.log_likelihood = log_likelihood; - (self.lambda, _) = match burke(&self.psi) { - Ok((lambda, objf)) => (lambda, objf), - Err(err) => { - bail!(err); - } - }; Ok(()) } fn condensation(&mut self) -> Result<()> { - let max_lambda = self - .lambda - .iter() - .fold(f64::NEG_INFINITY, |acc, x| x.max(acc)); - - let mut keep = Vec::::new(); - for (index, lam) in self.lambda.iter().enumerate() { - if lam > max_lambda / 1000_f64 { - keep.push(index); - } - } - if self.psi.matrix().ncols() != keep.len() { - tracing::debug!( - "Lambda (max/1000) dropped {} support point(s)", - self.psi.matrix().ncols() - keep.len(), - ); - } - - self.theta.filter_indices(keep.as_slice()); - self.psi.filter_column_indices(keep.as_slice()); - - let (r, perm) = qr::qrd(&self.psi)?; - - let mut keep = Vec::::new(); - let keep_n = self.psi.matrix().ncols().min(self.psi.matrix().nrows()); - for i in 0..keep_n { - let test = r.col(i).norm_l2(); - let r_diag_val = r.get(i, i); - let ratio = r_diag_val / test; - if ratio.abs() >= 1e-8 { - keep.push(*perm.get(i).unwrap()); - } - } - - if self.psi.matrix().ncols() != keep.len() { - tracing::debug!( - "QR decomposition dropped {} support point(s)", - self.psi.matrix().ncols() - keep.len(), - ); - } - - self.theta.filter_indices(keep.as_slice()); - self.psi.filter_column_indices(keep.as_slice()); + let (lambda, log_likelihood) = stages::condense( + &self.equation, + &self.data, + &self.error_models, + &mut self.theta, + &mut self.psi, + &self.lambda, + &CondensationOptions { + prune_threshold: self.config.prune_threshold, + qr_tolerance: self.config.qr_tolerance, + check_zero_probability: true, + }, + )?; - (self.lambda, self.objf) = match burke(&self.psi) { - Ok((lambda, objf)) => (lambda, objf), - Err(err) => { - return Err(anyhow::anyhow!("Error in IPM: {:?}", err)); - } - }; + self.lambda = lambda; + self.log_likelihood = log_likelihood; self.w = self.lambda.clone(); + Ok(()) } @@ -327,7 +288,7 @@ impl NonParametricRunner for NPOD { &self.theta, &mut self.error_models, &mut self.gamma_delta, - &mut self.objf, + &mut self.log_likelihood, &mut self.lambda, &mut self.psi, &self.config.error_optim, diff --git a/src/algorithms/nonparametric/stages.rs b/src/algorithms/nonparametric/stages.rs new file mode 100644 index 000000000..6cc0809ba --- /dev/null +++ b/src/algorithms/nonparametric/stages.rs @@ -0,0 +1,150 @@ +//! Cycle stages shared by the non-parametric algorithms. +//! +//! NPAG and NPOD run the same pipeline: estimate the likelihood of every +//! support point for every subject, condense the support points, then solve for +//! the weights. Only the tolerances and a few extra checks differ, so the +//! implementation lives here and each runner passes its configuration in. + +use anyhow::Result; + +use pharmsol::prelude::data::{AssayErrorModels, Data}; +use pharmsol::prelude::simulator::Equation; +use pharmsol::prelude::AssayErrorModel; + +use crate::algorithms::diagnostics; +use crate::estimation::nonparametric::{calculate_psi, ipm::burke, qr, Psi, Theta, Weights}; + +/// Estimate the likelihood of every support point for every subject, then solve +/// for the weights with the interior point method. +/// +/// Returns the likelihood matrix, the support point weights, and the +/// log-likelihood of the data under the weighted support points (higher is +/// better) — the quantity that [`FitState::n2ll`](crate::algorithms::FitState::n2ll) +/// negates and doubles. +pub(crate) fn estimate( + equation: &E, + data: &Data, + error_models: &AssayErrorModels, + theta: &Theta, + show_progress: bool, +) -> Result<(Psi, Weights, f64)> { + let psi = calculate_psi(equation, data, theta, error_models, show_progress)?; + diagnostics::check_zero_probability_subjects(equation, data, error_models, theta, &psi)?; + let (lambda, log_likelihood) = + burke(&psi).map_err(|err| anyhow::anyhow!("Error in IPM during estimation: {err:?}"))?; + + Ok((psi, lambda, log_likelihood)) +} + +/// How aggressively to drop support points during condensation. +pub(crate) struct CondensationOptions { + /// Support points whose weight is below `max_weight × prune_threshold` are dropped. + pub prune_threshold: f64, + /// Minimum `|r_ii| / ‖r_i‖` ratio for a support point to survive the QR decomposition. + pub qr_tolerance: f64, + /// Whether to re-check for zero-probability subjects after condensing. + pub check_zero_probability: bool, +} + +/// Drop support points that carry negligible weight or that are numerically +/// redundant, then solve for the weights of the survivors. +/// +/// `theta` and `psi` are filtered in place and stay aligned; the returned +/// weights and log-likelihood describe the surviving support points. +pub(crate) fn condense( + equation: &E, + data: &Data, + error_models: &AssayErrorModels, + theta: &mut Theta, + psi: &mut Psi, + lambda: &Weights, + options: &CondensationOptions, +) -> Result<(Weights, f64)> { + // Drop the support points with lambda < max(lambda) * prune_threshold + let max_lambda = lambda.iter().fold(f64::NEG_INFINITY, |acc, x| x.max(acc)); + + let mut keep = Vec::::new(); + for (index, lam) in lambda.iter().enumerate() { + if lam > max_lambda * options.prune_threshold { + keep.push(index); + } + } + if psi.matrix().ncols() != keep.len() { + tracing::debug!( + "Weight pruning (max × {}) dropped {} support point(s)", + options.prune_threshold, + psi.matrix().ncols() - keep.len(), + ); + } + + theta.filter_indices(keep.as_slice()); + psi.filter_column_indices(keep.as_slice()); + + // Rank-revealing factorization + let (r, perm) = qr::qrd(psi)?; + + let mut keep = Vec::::new(); + + // The minimum between the number of subjects and the actual number of support points + let keep_n = psi.matrix().ncols().min(psi.matrix().nrows()); + for i in 0..keep_n { + let test = r.col(i).norm_l2(); + let ratio = r.get(i, i) / test; + if ratio.abs() >= options.qr_tolerance { + keep.push(*perm.get(i).unwrap()); + } + } + + // If a support point is dropped, log it as a debug message + if psi.matrix().ncols() != keep.len() { + tracing::debug!( + "QR decomposition dropped {} support point(s)", + psi.matrix().ncols() - keep.len(), + ); + } + + theta.filter_indices(keep.as_slice()); + psi.filter_column_indices(keep.as_slice()); + + if options.check_zero_probability { + diagnostics::check_zero_probability_subjects(equation, data, error_models, theta, psi)?; + } + + let (lambda, log_likelihood) = + burke(psi).map_err(|err| anyhow::anyhow!("Error in IPM during condensation: {err:?}"))?; + + Ok((lambda, log_likelihood)) +} + +/// Log the objective function, the number of support points and the current +/// error models, warning when the log-likelihood decreased since the last cycle. +pub(crate) fn log_objective( + log_likelihood: f64, + last_log_likelihood: f64, + theta: &Theta, + error_models: &AssayErrorModels, +) { + tracing::info!("Objective function = {:.4}", -2.0 * log_likelihood); + tracing::debug!("Support points: {}", theta.nspp()); + error_models.iter().for_each(|(outeq, em)| { + if AssayErrorModel::None == *em { + return; + } + tracing::debug!( + "Error model for outeq {}: {:.4}", + outeq, + em.factor().unwrap_or_default() + ); + }); + + // A decreasing log-likelihood signals instability or model misspecification. + if last_log_likelihood > log_likelihood + 1e-4 { + tracing::warn!( + "Log-likelihood decreased from {:.4} to {:.4} (objective function increased from {:.4} to {:.4})", + last_log_likelihood, + log_likelihood, + -2.0 * last_log_likelihood, + -2.0 * log_likelihood + ); + } +} diff --git a/src/algorithms/parametric/mod.rs b/src/algorithms/parametric/mod.rs index 9951345f1..9fc8f3162 100644 --- a/src/algorithms/parametric/mod.rs +++ b/src/algorithms/parametric/mod.rs @@ -77,24 +77,21 @@ impl ParametricAlgorithm { } } -impl Algorithm for ParametricAlgorithm { +/// Every configuration that converts into a [`ParametricAlgorithm`] is itself an +/// [`Algorithm`], so [`SaemConfig`] can be passed straight to +/// [`fit_with`](crate::estimation::EstimationProblem::fit_with). +impl Algorithm> for A +where + E: Equation + Send + 'static, + A: Into, +{ type Output = ParametricResult; fn fit(self, _problem: EstimationProblem) -> Result { - match self { - Self::Saem(_config) => { + match self.into() { + ParametricAlgorithm::Saem(_config) => { unimplemented!("SAEM fitting is not yet implemented") } } } } - -// `SaemConfig` delegates to the matching `ParametricAlgorithm` variant so it can be passed -// directly to `fit_with`, keeping its setters compile-time checked. -impl Algorithm for SaemConfig { - type Output = ParametricResult; - - fn fit(self, problem: EstimationProblem) -> Result { - ParametricAlgorithm::from(self).fit(problem) - } -} diff --git a/src/algorithms/stop.rs b/src/algorithms/stop.rs new file mode 100644 index 000000000..028dfd558 --- /dev/null +++ b/src/algorithms/stop.rs @@ -0,0 +1,28 @@ +//! Out-of-band stopping of a running fit. +//! +//! PMcore follows the classic NONMEM/Pmetrics convention of a `stop` file in the +//! working directory: creating one asks the running fit to stop after the +//! current cycle. [`FitController::request_stop`](crate::algorithms::nonparametric::FitController::request_stop) +//! is the in-process equivalent. + +use std::path::Path; + +use anyhow::{Context, Result}; + +/// Name of the file whose presence asks a running fit to stop. +const STOP_FILE: &str = "stop"; + +/// Whether a stop has been requested through the stop file. +pub(crate) fn stop_file_present() -> bool { + Path::new(STOP_FILE).exists() +} + +/// Remove a stale stop file left over from a previous run. +pub(crate) fn remove_stop_file() -> Result<()> { + if stop_file_present() { + tracing::info!("Removing existing stop file prior to run"); + std::fs::remove_file(STOP_FILE).context("Unable to remove previous stop file")?; + } + + Ok(()) +} diff --git a/src/estimation/nonparametric/cycles.rs b/src/estimation/nonparametric/cycles.rs index 9c20d3539..884fc40c6 100644 --- a/src/estimation/nonparametric/cycles.rs +++ b/src/estimation/nonparametric/cycles.rs @@ -10,37 +10,43 @@ use crate::{ estimation::nonparametric::{median, theta::Theta, weights::Weights}, }; +/// A snapshot of the fit state after one cycle. +/// +/// The objective function is stored on the scale that is actually minimized, +/// `-2 × log-likelihood`, so that cycles logged by different algorithms are +/// directly comparable. #[derive(Debug, Clone, Serialize)] pub struct NPCycle { cycle: usize, - objf: f64, + n2ll: f64, error_models: AssayErrorModels, theta: Theta, weights: Weights, nspp: usize, - delta_objf: f64, + delta_log_likelihood: f64, status: Status, } impl NPCycle { pub fn new( cycle: usize, - objf: f64, + n2ll: f64, error_models: AssayErrorModels, theta: Theta, weights: Weights, - nspp: usize, - delta_objf: f64, + delta_log_likelihood: f64, status: Status, ) -> Self { + let nspp = theta.nspp(); + Self { cycle, - objf, + n2ll, error_models, theta, weights, nspp, - delta_objf, + delta_log_likelihood, status, } } @@ -48,8 +54,9 @@ impl NPCycle { pub fn cycle(&self) -> usize { self.cycle } - pub fn objf(&self) -> f64 { - self.objf + /// The objective function minimized by the algorithm, `-2 × log-likelihood`. + pub fn n2ll(&self) -> f64 { + self.n2ll } pub fn error_models(&self) -> &AssayErrorModels { &self.error_models @@ -63,8 +70,9 @@ impl NPCycle { pub fn nspp(&self) -> usize { self.nspp } - pub fn delta_objf(&self) -> f64 { - self.delta_objf + /// Change in log-likelihood since the previous cycle. + pub fn delta_log_likelihood(&self) -> f64 { + self.delta_log_likelihood } pub fn status(&self) -> &Status { &self.status @@ -73,12 +81,12 @@ impl NPCycle { pub fn placeholder() -> Self { Self { cycle: 0, - objf: 0.0, + n2ll: 0.0, error_models: AssayErrorModels::default(), theta: Theta::new(), weights: Weights::default(), nspp: 0, - delta_objf: 0.0, + delta_log_likelihood: 0.0, status: Status::Continue, } } @@ -149,10 +157,8 @@ impl CycleLog { writer.write_field(format!("{}", cycle.cycle))?; writer.write_field(format!("{}", cycle.status.converged()))?; writer.write_field(format!("{}", cycle.status))?; - writer.write_field(format!("{}", cycle.objf))?; - writer - .write_field(format!("{}", cycle.theta.nspp())) - .unwrap(); + writer.write_field(format!("{}", cycle.n2ll))?; + writer.write_field(format!("{}", cycle.nspp))?; cycle.error_models.iter().try_for_each( |(_, errmod): (usize, &AssayErrorModel)| -> Result<()> { diff --git a/src/estimation/nonparametric/result.rs b/src/estimation/nonparametric/result.rs index 7af6ca340..c526ab58d 100644 --- a/src/estimation/nonparametric/result.rs +++ b/src/estimation/nonparametric/result.rs @@ -18,7 +18,7 @@ pub struct NonParametricResult { theta: Theta, psi: Psi, weights: Weights, - objf: f64, + n2ll: f64, cycles: usize, status: Status, cyclelog: CycleLog, @@ -34,7 +34,7 @@ impl NonParametricResult { theta: Theta, psi: Psi, weights: Weights, - objf: f64, + n2ll: f64, cycles: usize, status: Status, cyclelog: CycleLog, @@ -47,21 +47,37 @@ impl NonParametricResult { theta, psi, weights, - objf, + n2ll, cycles, status, cyclelog, }) } + /// The number of cycles the algorithm ran + /// + /// This is the number of entries in the [cycle log](Self::cycle_log): every + /// cycle that was evaluated is logged, including the final one. Single-pass + /// algorithms such as NPMAP and NCNPAG therefore report `1`. pub fn cycles(&self) -> usize { self.cycles } - pub fn objf(&self) -> f64 { - self.objf + /// The objective function the algorithm minimized, `-2 × log-likelihood` + /// (lower is better). + /// + /// Reported on the same scale for every algorithm, so results from different + /// algorithms can be compared directly. + pub fn n2ll(&self) -> f64 { + self.n2ll + } + + /// The log-likelihood of the data under the fitted model (higher is better). + pub fn log_likelihood(&self) -> f64 { + -0.5 * self.n2ll } + /// Whether the fit stopped because it converged, rather than being cut short. pub fn converged(&self) -> bool { self.status.converged() } @@ -105,8 +121,7 @@ impl NonParametricResult { pub fn chain(self, algorithm: A) -> anyhow::Result> where A: crate::algorithms::Algorithm< - E, - crate::estimation::NonParametric, + crate::estimation::EstimationProblem, Output = NonParametricResult, >, E: crate::model::EquationMetadataSource, @@ -352,7 +367,7 @@ impl NonParametricResult { theta: &self.theta, psi: &self.psi, weights: &self.weights, - objf: self.objf, + n2ll: self.n2ll, cycles: self.cycles, status: &self.status, cyclelog: &self.cyclelog, @@ -413,7 +428,7 @@ struct NonParametricResultJson<'a> { theta: &'a Theta, psi: &'a Psi, weights: &'a Weights, - objf: f64, + n2ll: f64, cycles: usize, status: &'a Status, cyclelog: &'a CycleLog, @@ -499,7 +514,7 @@ mod tests { } #[test] - fn chain_npag_to_npag_maintains_or_improves_objf() { + fn chain_npag_to_npag_maintains_or_improves_n2ll() { let ode = minimal_ode(); let data = minimal_data(); let params = ParameterSpace::bounded() @@ -518,17 +533,17 @@ mod tests { .fit_with(NpagConfig::new().max_cycles(5)) .unwrap(); - let objf1 = r1.objf(); + let n2ll1 = r1.n2ll(); // Chain into another NPAG run let r2 = r1.chain(NpagConfig::new().max_cycles(3)).unwrap(); // Second run should not regress significantly assert!( - r2.objf() <= objf1 + 0.5, - "OBJF regressed: {} -> {}", - objf1, - r2.objf() + r2.n2ll() <= n2ll1 + 0.5, + "-2LL regressed: {} -> {}", + n2ll1, + r2.n2ll() ); } diff --git a/src/estimation/nonparametric/summaries.rs b/src/estimation/nonparametric/summaries.rs index 6d979aaef..a1de6616d 100644 --- a/src/estimation/nonparametric/summaries.rs +++ b/src/estimation/nonparametric/summaries.rs @@ -7,7 +7,7 @@ use crate::results::{FitSummary, IndividualSummary, ParameterSummary, Population pub fn fit_summary(result: &NonParametricResult) -> FitSummary { FitSummary { - objective_function: result.objf(), + objective_function: result.n2ll(), converged: result.converged(), iterations: result.cycles(), subject_count: result.data().subjects().len(), diff --git a/src/lib.rs b/src/lib.rs index 824aae3f9..ce1676b5b 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -58,7 +58,7 @@ pub mod prelude { pub use super::HashMap; pub use super::Result; pub use crate::algorithms; - pub use crate::algorithms::Algorithm; + pub use crate::algorithms::{Algorithm, FitState, NonParametricRunner, Status, StopReason}; pub use crate::estimation::NonParametric; pub use crate::estimation::Parametric; diff --git a/src/results/fit_result.rs b/src/results/fit_result.rs index eadcfeaaf..b9e22564d 100644 --- a/src/results/fit_result.rs +++ b/src/results/fit_result.rs @@ -5,7 +5,10 @@ use crate::results::{FitSummary, IndividualSummary, PopulationSummary}; /// A shared trait for the output of any estimation algorithm. pub trait FitResult { - fn objf(&self) -> f64; + /// The objective function minimized by the algorithm, `-2 × log-likelihood` + /// for likelihood-based algorithms (lower is better). + fn n2ll(&self) -> f64; + /// Whether the fit stopped because it converged, rather than being cut short. fn converged(&self) -> bool; fn summary(&self) -> FitSummary; fn population_summary(&self) -> PopulationSummary; @@ -20,7 +23,7 @@ pub struct ParametricResult { } impl FitResult for ParametricResult { - fn objf(&self) -> f64 { + fn n2ll(&self) -> f64 { unimplemented!("Parametric result not yet implemented") } fn converged(&self) -> bool { @@ -40,12 +43,12 @@ impl FitResult for ParametricResult { use crate::estimation::nonparametric; impl FitResult for NonParametricResult { - fn objf(&self) -> f64 { - self.objf() // Assuming the struct has this native method + fn n2ll(&self) -> f64 { + NonParametricResult::n2ll(self) } fn converged(&self) -> bool { - self.converged() // Assuming the struct has this native method + NonParametricResult::converged(self) } fn summary(&self) -> FitSummary { diff --git a/tests/onecomp.rs b/tests/onecomp.rs index c63509b52..44560dc1a 100644 --- a/tests/onecomp.rs +++ b/tests/onecomp.rs @@ -66,7 +66,7 @@ fn test_one_compartment_npag() -> Result<()> { // Check the results assert_eq!(result.cycles(), 31); - assert!(result.objf() - 565.7749 < 0.01); + assert!(result.n2ll() - 565.7749 < 0.01); // The prior is preserved on the result and is distinct from the optimized // solution (which is condensed to far fewer support points). @@ -132,7 +132,7 @@ fn test_one_compartment_npod() -> Result<()> { // Check the results assert_eq!(result.cycles(), 11); - assert!(result.objf() - 565.7749 < 0.01); + assert!(result.n2ll() - 565.7749 < 0.01); Ok(()) } @@ -194,8 +194,9 @@ fn test_one_compartment_postprob() -> Result<()> { let result = EstimationProblem::nonparametric(eq, data, theta.clone(), error_models)? .fit_with(NonParametricAlgorithm::npmap())?; - // Check the results - assert_eq!(result.cycles(), 0); + // NPMAP is single-pass: exactly one cycle is evaluated and logged + assert_eq!(result.cycles(), 1); + assert_eq!(result.cycle_log().cycles().len(), 1); // Should be 100 points in theta (no change in points) assert_eq!(result.get_theta().nspp(), theta.nspp());