diff --git a/Cargo.toml b/Cargo.toml index d8a34c945..07d22858c 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -28,7 +28,7 @@ tracing-subscriber = { version = "0.3.19", features = [ "time", ] } faer = "0.24.0" -pharmsol = { git = "https://github.com/LAPKB/pharmsol.git", branch = "claude/sweet-davinci-o6m0ej" } +pharmsol = { git = "https://github.com/LAPKB/pharmsol.git", branch = "fix/ode-newton-budget" } anyhow = "1.0.100" rayon = "1.10.0" rand = "0.10.1" diff --git a/benches/bimodal_ke.rs b/benches/bimodal_ke.rs index 9c7cf1477..e3dd8d470 100644 --- a/benches/bimodal_ke.rs +++ b/benches/bimodal_ke.rs @@ -76,30 +76,15 @@ where } fn benchmark_bimodal_ke_npag(c: &mut Criterion) { - benchmark_algorithm( - c, - "bimodal_ke_npag", - setup_npag, - NonParametricAlgorithm::npag(), - ); + benchmark_algorithm(c, "bimodal_ke_npag", setup_npag, NpagConfig::new()); } fn benchmark_bimodal_ke_npod(c: &mut Criterion) { - benchmark_algorithm( - c, - "bimodal_ke_npod", - setup_npod, - NonParametricAlgorithm::npod(), - ); + benchmark_algorithm(c, "bimodal_ke_npod", setup_npod, NpodConfig::new()); } fn benchmark_bimodal_ke_postprob(c: &mut Criterion) { - benchmark_algorithm( - c, - "bimodal_ke_postprob", - setup_postprob, - NonParametricAlgorithm::npmap(), - ); + benchmark_algorithm(c, "bimodal_ke_postprob", setup_postprob, NpmapConfig::new()); } criterion_group! { diff --git a/examples/bestdose/main.rs b/examples/bestdose/main.rs index 504ff0d5e..227bef2b5 100644 --- a/examples/bestdose/main.rs +++ b/examples/bestdose/main.rs @@ -116,7 +116,7 @@ fn concentration_target( )? .fit_with(NcnpagConfig::default())?; - let post_theta = posterior.get_theta().clone(); + let post_theta = posterior.theta().clone(); let post_weights = posterior.weights().clone(); println!( "NCNPAG posterior: {} support points\n", diff --git a/examples/bimodal_ke/main.rs b/examples/bimodal_ke/main.rs index 6e8a035b5..7ee38ec86 100644 --- a/examples/bimodal_ke/main.rs +++ b/examples/bimodal_ke/main.rs @@ -33,9 +33,9 @@ fn main() -> Result<()> { let problem = EstimationProblem::nonparametric(eq, data, prior, error_models)?; - let _result = problem.fit_with(NonParametricAlgorithm::npag())?; + let _result = problem.fit_with(NpagConfig::new())?; - _result.write_outputs("outputs/bimodal_ke", 0.0, 0.0)?; + _result.write_outputs("outputs/bimodal_ke", PredictionGrid::new())?; Ok(()) } diff --git a/examples/bimodal_ke_backend_compare.rs b/examples/bimodal_ke_backend_compare.rs index eab7dab1b..726c4d586 100644 --- a/examples/bimodal_ke_backend_compare.rs +++ b/examples/bimodal_ke_backend_compare.rs @@ -115,7 +115,7 @@ fn run_case Result<()> { println!( "Chained NPAG→NPOD: OBJF = {:.2}, {} support points, {} total cycles", result.objf(), - result.get_theta().nspp(), + result.theta().nspp(), result.cycles(), ); diff --git a/examples/drusano/main.rs b/examples/drusano/main.rs index 2cf71f81a..fa100a5df 100644 --- a/examples/drusano/main.rs +++ b/examples/drusano/main.rs @@ -100,7 +100,6 @@ fn main() -> Result<()> { "4", AssayErrorModel::proportional(ErrorPoly::new(0.1, 0.1, 0.0, 0.0), 1.0), )?; - EstimationProblem::nonparametric(eq, data, prior, error_models)? - .fit_with(NonParametricAlgorithm::npag())?; + EstimationProblem::nonparametric(eq, data, prior, error_models)?.fit_with(NpagConfig::new())?; Ok(()) } diff --git a/examples/iov/main.rs b/examples/iov/main.rs index ecf4181b7..2feaac809 100644 --- a/examples/iov/main.rs +++ b/examples/iov/main.rs @@ -34,7 +34,7 @@ fn main() -> Result<()> { AssayErrorModel::additive(ErrorPoly::new(0.0, 0.0, 0.0, 0.0), 0.0000757575757576), )?; EstimationProblem::nonparametric(sde, data, prior, error_models)? - .fit_with(NonParametricAlgorithm::npag()) + .fit_with(NpagConfig::new()) .unwrap(); Ok(()) diff --git a/examples/meta/main.rs b/examples/meta/main.rs index 70a5cf104..4143274d1 100644 --- a/examples/meta/main.rs +++ b/examples/meta/main.rs @@ -47,8 +47,7 @@ fn main() -> Result<()> { "2", AssayErrorModel::proportional(ErrorPoly::new(1.0, 0.1, 0.0, 0.0), 5.0), )?; - EstimationProblem::nonparametric(eq, data, prior, error_models)? - .fit_with(NonParametricAlgorithm::npod())?; + EstimationProblem::nonparametric(eq, data, prior, error_models)?.fit_with(NpodConfig::new())?; Ok(()) } diff --git a/examples/neely/main.rs b/examples/neely/main.rs index d059e01cc..2d52a3c7d 100644 --- a/examples/neely/main.rs +++ b/examples/neely/main.rs @@ -69,8 +69,7 @@ fn main() -> Result<()> { "3", AssayErrorModel::proportional(ErrorPoly::new(1.0, 0.1, 0.0, 0.0), 5.0), )?; - EstimationProblem::nonparametric(eq, data, prior, error_models)? - .fit_with(NonParametricAlgorithm::npag())?; + EstimationProblem::nonparametric(eq, data, prior, error_models)?.fit_with(NpagConfig::new())?; Ok(()) } diff --git a/examples/new_iov/main.rs b/examples/new_iov/main.rs index 74d42c68c..4fb57810a 100644 --- a/examples/new_iov/main.rs +++ b/examples/new_iov/main.rs @@ -34,8 +34,7 @@ fn main() -> Result<()> { "1", AssayErrorModel::additive(ErrorPoly::new(-0.00119, 0.44379, -0.45864, 0.16537), 0.0), )?; - EstimationProblem::nonparametric(sde, data, prior, error_models)? - .fit_with(NonParametricAlgorithm::npag())?; + EstimationProblem::nonparametric(sde, data, prior, error_models)?.fit_with(NpagConfig::new())?; Ok(()) } diff --git a/examples/theophylline/main.rs b/examples/theophylline/main.rs index c4d9ad84e..ba1c536e5 100644 --- a/examples/theophylline/main.rs +++ b/examples/theophylline/main.rs @@ -26,7 +26,7 @@ fn main() -> Result<()> { AssayErrorModel::proportional(ErrorPoly::new(0.1, 0.1, 0.0, 0.0), 2.0), )?; EstimationProblem::nonparametric(analytical, data, prior, error_models)? - .fit_with(NonParametricAlgorithm::npag())?; + .fit_with(NpagConfig::new())?; Ok(()) } diff --git a/examples/two_eq_lag/main.rs b/examples/two_eq_lag/main.rs index 3b6f25235..0498e2807 100644 --- a/examples/two_eq_lag/main.rs +++ b/examples/two_eq_lag/main.rs @@ -29,8 +29,7 @@ fn main() -> Result<()> { "0", AssayErrorModel::additive(ErrorPoly::new(-0.00119, 0.44379, -0.45864, 0.16537), 0.0), )?; - EstimationProblem::nonparametric(eq, data, prior, error_models)? - .fit_with(NonParametricAlgorithm::npag())?; + EstimationProblem::nonparametric(eq, data, prior, error_models)?.fit_with(NpagConfig::new())?; Ok(()) } diff --git a/examples/vanco_sde/main.rs b/examples/vanco_sde/main.rs index d72db2d28..8e2675618 100644 --- a/examples/vanco_sde/main.rs +++ b/examples/vanco_sde/main.rs @@ -60,8 +60,7 @@ fn main() -> Result<()> { "1", AssayErrorModel::additive(ErrorPoly::new(0.00119, 0.20, 0.0, 0.0), 0.0), )?; - EstimationProblem::nonparametric(sde, data, prior, error_models)? - .fit_with(NonParametricAlgorithm::npag())?; + EstimationProblem::nonparametric(sde, data, prior, error_models)?.fit_with(NpagConfig::new())?; Ok(()) } diff --git a/src/algorithms/mod.rs b/src/algorithms/mod.rs index a6455ba31..7fa38bbdc 100644 --- a/src/algorithms/mod.rs +++ b/src/algorithms/mod.rs @@ -1,17 +1,8 @@ -use std::fs; -use std::path::Path; - -use crate::estimation::nonparametric::{NonParametricResult, Psi, Theta}; 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::{Predictions, Subject}; +use pharmsol::prelude::simulator::Equation; use serde::{Deserialize, Serialize}; /// Defines an algorithm that can fit an [`EstimationProblem`] to produce a result. @@ -19,7 +10,7 @@ use serde::{Deserialize, Serialize}; /// 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 { +pub trait Algorithm { /// The specific result struct (e.g. `NonParametricResult`). type Output: FitResult; @@ -43,259 +34,6 @@ impl EstimationProblem { } } -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, &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; - /// 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 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 [Status] of the algorithm - fn status(&self) -> &Status; - /// Set the current [Status] of the algorithm - fn set_status(&mut self, status: Status); - /// Evaluate convergence criteria and update status - 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(()) - } - fn estimation(&mut self) -> Result<()>; - /// 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<()>; - - /// 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<()>; - - /// 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<()>; - - /// 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. - fn next_cycle(&mut self) -> Result { - let cycle = self.increment_cycle(); - - if cycle > 1 { - self.expansion()?; - } - - let span = tracing::info_span!("", "{}", format!("Cycle {}", self.cycle())); - let _enter = span.enter(); - self.estimation()?; - self.condensation()?; - self.optimizations()?; - self.evaluation() - } - - /// Fit the model until convergence or stopping criteria are met - /// - /// 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() - } - - #[allow(clippy::wrong_self_convention)] - fn into_result(&self) -> Result>; -} - /// Where a fit stands: still running, or stopped (and why). #[derive(Debug, Clone, PartialEq, Eq, Serialize, Deserialize)] pub enum Status { @@ -346,8 +84,9 @@ pub enum StopReason { MaxCycles, /// A `stop` file was found on disk. StopFile, - /// Stopped from code — [`request_stop`](crate::algorithms::nonparametric::FitController::request_stop) - /// or an observer returning [`CycleFlow::Stop`](crate::algorithms::nonparametric::CycleFlow::Stop). + /// Stopped from code, with + /// [`FitController::request_stop`](crate::algorithms::nonparametric::FitController::request_stop) + /// or [`FitController::run_until`](crate::algorithms::nonparametric::FitController::run_until). Aborted, } diff --git a/src/algorithms/nonparametric/controller.rs b/src/algorithms/nonparametric/controller.rs index 22a714290..d679eeb28 100644 --- a/src/algorithms/nonparametric/controller.rs +++ b/src/algorithms/nonparametric/controller.rs @@ -1,144 +1,149 @@ //! Step through a fit cycle by cycle, rather than running it all at once. //! //! [`fit_with`](crate::estimation::EstimationProblem::fit_with) runs to the end and gives back -//! only the final result. For debugging or third-party monitoring, you can instead drive the fit yourself and inspect the state between cycles. -//! -//! - [`FitController`] — hold the fit and advance it one [`step`](FitController::step) at a -//! time, inspecting state in between. -//! - [`fit_with_observer`](crate::estimation::EstimationProblem::fit_with_observer) — let it -//! run, but get a callback after every cycle. -//! -//! Both drive the same underlying runner and finish at the same result. +//! only the final result. For debugging or monitoring, start a [`FitController`] with +//! [`fit_controller`](crate::estimation::EstimationProblem::fit_controller) instead, and either +//! advance it one [`step`](FitController::step) at a time or let it +//! [`run_until`](FitController::run_until) a condition of yours. Every route runs the same fit +//! and ends in the same result. use anyhow::Result; use pharmsol::prelude::{data::AssayErrorModels, simulator::Equation}; -use crate::algorithms::{NonParametricRunner, Status, StopReason}; -use crate::estimation::nonparametric::{NonParametricResult, Psi, Theta}; +use crate::algorithms::Status; +use crate::estimation::nonparametric::{NonParametricResult, Psi, Theta, Weights}; use crate::estimation::{EstimationProblem, NonParametric}; +use super::fit::Fit; use super::NonParametricAlgorithm; /// A running fit that you advance one cycle at a time. /// /// It's ready to go as soon as it's built, so the first [`step`](Self::step) runs cycle 1. -/// The accessors read the current state between steps, and [`finish`](Self::finish) or +/// The accessors read the current state between steps. Once the fit has stopped, whether on +/// its own or through [`request_stop`](Self::request_stop), it stays stopped: further steps do +/// nothing. [`finish`](Self::finish), [`run_until`](Self::run_until) and /// [`into_result`](Self::into_result) turn it into a result. -pub struct FitController { - runner: Box>, +pub struct FitController { + fit: Fit, } -impl FitController { +impl FitController { /// Build a controller and initialize it, ready for the first [`step`](Self::step). pub(crate) fn new( algorithm: NonParametricAlgorithm, problem: EstimationProblem, ) -> Result { - let mut runner = algorithm.into_runner(problem)?; - runner.initialize()?; - Ok(Self { runner }) + Ok(Self { + fit: algorithm.start(problem)?, + }) } /// Run one cycle and return the new [`Status`]. /// - /// [`Continue`](Status::Continue) means there's more to do; [`Stop`](Status::Stop) means the - /// algorithm is finished and further steps won't change the outcome. + /// [`Continue`](Status::Continue) means there's more to do. [`Stop`](Status::Stop) means the + /// fit is finished: from then on, `step` does nothing and returns the same status. pub fn step(&mut self) -> Result { - self.runner.next_cycle() + self.fit.step() } - /// Flag the fit as [`Aborted`](StopReason::Aborted) so the result shows an early stop. + /// Stop the fit now, as [`Aborted`](crate::algorithms::StopReason::Aborted). /// - /// 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). + /// Further steps do nothing. A fit that has already stopped keeps its original reason. pub fn request_stop(&mut self) { - self.runner.set_status(Status::Stop(StopReason::Aborted)); + self.fit.request_stop(); } /// The current cycle number (0 before the first [`step`](Self::step)). pub fn cycle(&self) -> usize { - self.runner.cycle() + self.fit.state().cycle() } /// The current [`Status`] of the fit. pub fn status(&self) -> &Status { - self.runner.status() + self.fit.state().status() } - /// The current support points and weights. + /// The current support points. pub fn theta(&self) -> &Theta { - self.runner.theta() + self.fit.state().solution().theta() + } + + /// The current weight of each support point in [`theta`](Self::theta). + pub fn weights(&self) -> &Weights { + self.fit.state().solution().weights() } /// The current likelihood matrix. pub fn psi(&self) -> &Psi { - self.runner.psi() + self.fit.state().solution().psi() } /// The current log-likelihood. - pub fn likelihood(&self) -> f64 { - self.runner.likelihood() + pub fn log_likelihood(&self) -> f64 { + self.fit.state().solution().log_likelihood() } - /// The current negative two log-likelihood (objective function). - pub fn n2ll(&self) -> f64 { - self.runner.n2ll() + /// The current objective function, -2 × log-likelihood, as reported by + /// [`NonParametricResult::objf`]. + pub fn objf(&self) -> f64 { + self.fit.state().objf() } /// The current error models (updated when the algorithm optimizes them). pub fn error_models(&self) -> &AssayErrorModels { - self.runner.error_models() + self.fit.state().error_models() } /// Run the fit to completion and return the result. /// /// Safe to call anytime: if the fit has already stopped, it just returns the result. pub fn finish(mut self) -> Result> { + self.fit.run()?; + self.fit.into_result() + } + + /// Run until the fit stops on its own or `stop` returns `true`, then return the result. + /// + /// `stop` is called after every cycle, the last one included, so it can also report + /// progress. If it returns `true` while the fit is still running, the fit stops as + /// [`Aborted`](crate::algorithms::StopReason::Aborted). + /// + /// ```no_run + /// use pmcore::prelude::*; + /// # fn run(problem: EstimationProblem) -> Result<()> { + /// let result = problem + /// .fit_controller(NpagConfig::new())? + /// .run_until(|fit| { + /// println!("cycle {} | -2LL {:.4}", fit.cycle(), fit.objf()); + /// fit.cycle() >= 50 + /// })?; + /// # let _ = result; Ok(()) + /// # } + /// ``` + pub fn run_until( + mut self, + mut stop: impl FnMut(&Self) -> bool, + ) -> Result> { while self.status().is_continue() { self.step()?; + if stop(&self) { + self.request_stop(); + } } - self.runner.into_result() + self.into_result() } /// Turn the current state into a result as-is, without running more cycles. pub fn into_result(self) -> Result> { - self.runner.into_result() - } -} - -/// An observer's verdict after a cycle: keep going, or stop. -#[derive(Debug, Clone, Copy, PartialEq, Eq)] -pub enum CycleFlow { - /// Keep going — the algorithm may still stop on its own. - Continue, - /// Stop now. - Stop, -} - -/// A callback run after every cycle of a fit. -/// -/// It receives the live [`FitController`] to read the current state, and returns a -/// [`CycleFlow`] to keep going or bail out. Any `FnMut(&FitController) -> CycleFlow` works, -/// so you can hand a plain closure to -/// [`fit_with_observer`](crate::estimation::EstimationProblem::fit_with_observer). -pub trait FitObserver { - fn on_cycle(&mut self, controller: &FitController) -> CycleFlow; -} - -impl FitObserver for F -where - E: Equation + Send + 'static, - F: FnMut(&FitController) -> CycleFlow, -{ - fn on_cycle(&mut self, controller: &FitController) -> CycleFlow { - self(controller) + self.fit.into_result() } } -impl EstimationProblem { +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 + /// Takes any non-parametric config (`NpagConfig`, `NpodConfig`, `NpmapConfig`, `NcnpagConfig`) or a /// [`NonParametricAlgorithm`]. /// /// ```no_run @@ -146,7 +151,7 @@ impl EstimationProblem { /// # fn run(problem: EstimationProblem) -> Result<()> { /// let mut controller = problem.fit_controller(NpagConfig::new())?; /// while controller.step()?.is_continue() { - /// println!("cycle {} | -2LL {:.4}", controller.cycle(), controller.n2ll()); + /// println!("cycle {} | -2LL {:.4}", controller.cycle(), controller.objf()); /// } /// let result = controller.into_result()?; /// # let _ = result; Ok(()) @@ -158,41 +163,4 @@ impl EstimationProblem { ) -> Result> { FitController::new(algorithm.into(), self) } - - /// Run a fit to completion, calling `observer` after every cycle. - /// - /// The observer sees every cycle, the last one included, and can return - /// [`CycleFlow::Stop`] to bail out early. Returns the same result as - /// [`fit_with`](Self::fit_with). - /// - /// ```no_run - /// use pmcore::prelude::*; - /// # fn run(problem: EstimationProblem) -> Result<()> { - /// let result = problem.fit_with_observer(NpagConfig::new(), |c: &FitController<_>| { - /// println!("cycle {} | -2LL {:.4}", c.cycle(), c.n2ll()); - /// CycleFlow::Continue - /// })?; - /// # let _ = result; Ok(()) - /// # } - /// ``` - pub fn fit_with_observer>( - self, - algorithm: impl Into, - mut observer: O, - ) -> Result> { - let mut controller = self.fit_controller(algorithm)?; - loop { - let status = controller.step()?; - let flow = observer.on_cycle(&controller); - // Stop on the algorithm's own criteria, or when the observer asks us to. - if status.is_stop() { - break; - } - if flow == CycleFlow::Stop { - controller.request_stop(); - break; - } - } - controller.into_result() - } } diff --git a/src/algorithms/nonparametric/error_optim.rs b/src/algorithms/nonparametric/error_optim.rs index 1fab17535..76123e8dd 100644 --- a/src/algorithms/nonparametric/error_optim.rs +++ b/src/algorithms/nonparametric/error_optim.rs @@ -1,13 +1,13 @@ //! Error-model factor optimization used by non-parametric algorithms. use anyhow::Result; -use pharmsol::prelude::{ - data::{AssayErrorModels, Data}, - simulator::Equation, -}; +use pharmsol::prelude::{data::AssayErrorModels, simulator::Equation}; +use pharmsol::OutputLabel; use serde::{Deserialize, Serialize}; -use crate::estimation::nonparametric::{calculate_psi, ipm::burke, Psi, Theta, Weights}; +use super::fit::FitState; +use super::solution::Solution; +use crate::estimation::nonparametric::calculate_psi; /// Configuration for the error-model factor (gamma/lambda) optimization. #[derive(Debug, Clone, PartialEq, Serialize, Deserialize)] @@ -59,110 +59,112 @@ impl ErrorOptimConfig { } } -/// Perform one pass of the error-model factor (gamma/lambda) optimization. +/// Optimizes the error-model factors (gamma/lambda) of a fit, one pass per cycle. /// -/// 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 -/// 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 -/// (that direction is skipped); if both directions fail the factor is left -/// unchanged for this pass. -#[allow(clippy::too_many_arguments)] -pub(crate) fn optimize_error_models( - equation: &E, - data: &Data, - theta: &Theta, - error_models: &mut AssayErrorModels, - gamma_delta: &mut [f64], - objf: &mut f64, - lambda: &mut Weights, - psi: &mut Psi, - config: &ErrorOptimConfig, -) -> Result<()> { - error_models - .clone() - .iter_mut() - .enumerate() - .filter_map(|(slot, (outeq, em))| { - if em.optimize() { - Some((slot, outeq, em)) - } else { - None - } - }) - .try_for_each(|(slot, outeq, em)| -> Result<()> { - // Optimize the best value of gamma/lambda - - let gamma_up = em.factor()? * (1.0 + gamma_delta[slot]); - let gamma_down = em.factor()? / (1.0 + gamma_delta[slot]); - - let mut error_model_up = error_models.clone(); - error_model_up.set_factor(outeq, gamma_up)?; - - let mut error_model_down = error_models.clone(); - error_model_down.set_factor(outeq, gamma_down)?; - - let psi_up = calculate_psi(equation, data, theta, &error_model_up, false)?; - let psi_down = calculate_psi(equation, data, theta, &error_model_down, false)?; - - // Treat errors in IPM, such as subjects with zero probability or non-finite - // 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)), - Err(err) => { - tracing::warn!( - "Error in IPM during optim (up) for outeq {}: {:?}", - outeq, - err - ); - None - } - }; - let down = match burke(&psi_down) { - Ok((lambda, objf)) => Some((lambda, objf)), - Err(err) => { - tracing::warn!( - "Error in IPM during optim (down) for outeq {}: {:?}", - outeq, - err - ); - None - } - }; - - // Select the best improving candidate (if any) over the current - // objective. Among the two directions, the one with the higher - // objective function 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_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((new_objf, new_lambda, new_psi, gamma)) = best { - error_models.set_factor(outeq, gamma)?; - *objf = new_objf; - gamma_delta[slot] *= config.growth; - *lambda = new_lambda; - *psi = new_psi; +/// Holds the per-output adaptive step that carries over between cycles. +#[derive(Debug)] +pub(crate) struct ErrorModelOptimizer { + config: ErrorOptimConfig, + /// Current relative perturbation step for each error model, in iteration order. + steps: Vec, +} + +impl ErrorModelOptimizer { + pub(crate) fn new(config: ErrorOptimConfig, error_models: &AssayErrorModels) -> Self { + let steps = vec![config.step; error_models.len()]; + Self { config, steps } + } + + /// Run one optimization pass over every optimizable error model in `state`. + /// + /// Each factor is perturbed up and down, and the current support points are re-solved + /// under both. The direction with the highest log-likelihood is adopted if it improves on + /// the current one. The step grows on improvement and shrinks every pass, resetting to + /// `config.step` once it falls to `config.min_step`. + /// + /// An IPM failure in one direction is logged and that direction is skipped; if both fail + /// the factor is left unchanged for this pass. + pub(crate) fn optimize(&mut self, state: &mut FitState) -> Result<()> { + let targets: Vec<(usize, OutputLabel, f64)> = state + .error_models() + .iter() + .enumerate() + .filter(|(_, (_, em))| em.optimize()) + .map(|(slot, (outeq, em))| Ok((slot, outeq.clone(), em.factor()?))) + .collect::>()?; + + for (slot, outeq, factor) in targets { + let step = &mut self.steps[slot]; + let gamma_up = factor * (1.0 + *step); + let gamma_down = factor / (1.0 + *step); + + let up = evaluate_factor(state, &outeq, gamma_up, "up")?; + let down = evaluate_factor(state, &outeq, gamma_down, "down")?; + + // Of the directions that improve on the current log-likelihood, keep the best. + let current = state.solution().log_likelihood(); + let best = [up, down] + .into_iter() + .flatten() + .filter(|candidate| candidate.solution.log_likelihood() > current) + .reduce(|best, candidate| { + if candidate.solution.log_likelihood() > best.solution.log_likelihood() { + candidate + } else { + best + } + }); + + if let Some(best) = best { + let weights = state.solution().weights().clone(); + state.adopt_error_models( + best.error_models, + best.solution.with_stale_weights(weights), + ); + *step *= self.config.growth; } - gamma_delta[slot] *= config.shrink; - if gamma_delta[slot] <= config.min_step { - gamma_delta[slot] = config.step; + *step *= self.config.shrink; + if *step <= self.config.min_step { + *step = self.config.step; } - Ok(()) - })?; + } + Ok(()) + } +} + +/// The error models with one factor set to a trial value, and the solution they give. +struct Candidate { + error_models: AssayErrorModels, + solution: Solution, +} - Ok(()) +/// Re-solve the current support points with the factor of `outeq` set to `factor`. +/// +/// Returns `None` (after logging a warning) if the IPM fails for this factor. +fn evaluate_factor( + state: &FitState, + outeq: &OutputLabel, + factor: f64, + direction: &str, +) -> Result> { + let mut error_models = state.error_models().clone(); + error_models.set_factor(outeq, factor)?; + let theta = state.solution().theta(); + let psi = calculate_psi(state.equation(), state.data(), theta, &error_models, false)?; + + match Solution::optimal(theta.clone(), psi) { + Ok(solution) => Ok(Some(Candidate { + error_models, + solution, + })), + Err(err) => { + tracing::warn!( + "Error in IPM during optim ({}) for outeq {}: {:?}", + direction, + outeq, + err + ); + Ok(None) + } + } } diff --git a/src/algorithms/nonparametric/fit.rs b/src/algorithms/nonparametric/fit.rs new file mode 100644 index 000000000..dbde7c603 --- /dev/null +++ b/src/algorithms/nonparametric/fit.rs @@ -0,0 +1,521 @@ +//! The machinery every non-parametric algorithm runs on. +//! +//! A fit is split in two: +//! +//! - [`FitState`] describes the fit at a point in time: the problem (equation, data, error +//! models, prior), the current [`Solution`], and the cycle bookkeeping. It is the same for +//! every algorithm. +//! - A [`Runner`] holds only what is specific to one algorithm (its configuration and any +//! adaptive state, such as NPAG's grid resolution) and knows how to run a single cycle. +//! +//! [`Fit`] pairs the two and owns everything around the cycle: the counter, the cycle cap, the +//! stop file, logging, and the cycle log. Both +//! [`fit_with`](crate::estimation::EstimationProblem::fit_with) and the stepping +//! [`FitController`](super::FitController) are thin wrappers around it. + +use std::fs; +use std::path::Path; + +use anyhow::{Context, Result}; +use pharmsol::prelude::{ + data::{AssayErrorModel, AssayErrorModels, Data}, + simulator::Equation, +}; +use pharmsol::{Predictions, Subject}; +use rayon::prelude::*; + +use super::solution::Solution; +use crate::algorithms::{Status, StopReason}; +use crate::estimation::nonparametric::{ + calculate_psi, qr, CycleLog, NPCycle, NonParametricResult, Psi, Theta, ZeroProbabilityError, +}; +use crate::estimation::{EstimationProblem, NonParametric}; + +/// A file with this name in the working directory stops a watched fit after the current cycle. +const STOP_FILE: &str = "stop"; + +/// What a runner concluded at the end of a cycle. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub(crate) enum CycleOutcome { + /// The algorithm wants another cycle. + Continue, + /// The algorithm's own convergence criterion is met. + Converged, +} + +/// One non-parametric algorithm: what a single cycle does to a [`FitState`]. +pub(crate) trait Runner: Send + Sync { + /// Run cycle [`state.cycle()`](FitState::cycle), replacing the solution (and, if it + /// optimizes them, the error models) as it goes. + /// + /// Everything else (the cycle cap, the stop file, logging, and the cycle log) is handled + /// by [`Fit`]. + fn step(&mut self, state: &mut FitState) -> Result; +} + +/// The state of a non-parametric fit, shared by every algorithm. +/// +/// Runners can read everything but only change the solution and the error models, and only +/// through [`set_solution`](Self::set_solution) and +/// [`adopt_error_models`](Self::adopt_error_models). The cycle bookkeeping belongs to [`Fit`]. +#[derive(Debug)] +pub(crate) struct FitState { + equation: E, + data: Data, + error_models: AssayErrorModels, + /// The support points the fit started from. + prior: Theta, + solution: Solution, + /// Log-likelihood at the end of the previous cycle, or `None` during the first one. + previous_log_likelihood: Option, + /// The current cycle, starting at 1 (0 before the first cycle runs). + cycle: usize, + status: Status, + cycle_log: CycleLog, +} + +impl FitState { + pub(crate) fn new(problem: EstimationProblem) -> Self { + Self::from_parts( + problem.model.equation, + problem.data, + problem.error_models, + problem.prior, + ) + } + + pub(crate) fn from_parts( + equation: E, + data: Data, + error_models: AssayErrorModels, + prior: Theta, + ) -> Self { + Self { + equation, + data, + error_models, + solution: Solution::unsolved(prior.clone()), + prior, + previous_log_likelihood: None, + cycle: 0, + status: Status::Continue, + cycle_log: CycleLog::new(), + } + } + + pub(crate) fn equation(&self) -> &E { + &self.equation + } + + pub(crate) fn data(&self) -> &Data { + &self.data + } + + pub(crate) fn error_models(&self) -> &AssayErrorModels { + &self.error_models + } + + pub(crate) fn solution(&self) -> &Solution { + &self.solution + } + + pub(crate) fn set_solution(&mut self, solution: Solution) { + self.solution = solution; + } + + /// Switch to new error models, together with the solution computed under them. + pub(crate) fn adopt_error_models( + &mut self, + error_models: AssayErrorModels, + solution: Solution, + ) { + self.error_models = error_models; + self.solution = solution; + } + + pub(crate) fn cycle(&self) -> usize { + self.cycle + } + + pub(crate) fn is_first_cycle(&self) -> bool { + self.cycle == 1 + } + + pub(crate) fn status(&self) -> &Status { + &self.status + } + + /// The objective function, -2 × log-likelihood. + pub(crate) fn objf(&self) -> f64 { + -2.0 * self.solution.log_likelihood() + } + + /// Absolute change in log-likelihood since the previous cycle (`None` in the first cycle). + pub(crate) fn objective_change(&self) -> Option { + self.previous_log_likelihood + .map(|previous| (previous - self.solution.log_likelihood()).abs()) + } + + /// Compute [`Psi`] for `theta` under the current error models. + /// + /// If some subjects have zero probability at every support point, likelihood diagnostics + /// for them are logged before the [`ZeroProbabilityError`] is returned. + pub(crate) fn compute_psi(&self, theta: &Theta, progress: bool) -> Result { + calculate_psi( + &self.equation, + &self.data, + theta, + &self.error_models, + progress, + ) + .inspect_err(|err| { + if let Some(zero) = err.downcast_ref::() { + for subject in self.data.subjects() { + if zero.subjects().contains(subject.id()) { + self.log_zero_probability_subject(subject, theta); + } + } + } + }) + } + + /// Evaluate `theta` and solve for its optimal weights. + pub(crate) fn estimate(&self, theta: Theta, progress: bool) -> Result { + let psi = self.compute_psi(&theta, progress)?; + self.check_zero_probability_subjects(&theta, &psi)?; + Solution::optimal(theta, psi).context("IPM failed during estimation") + } + + /// Drop support points that carry negligible weight or are linearly dependent, then solve + /// for the weights of the survivors. + /// + /// A point is dropped when its weight in `estimate` is at most `prune_threshold` × the + /// largest weight, or when its column of [`Psi`] is (numerically) a combination of the + /// others, as judged by a rank-revealing QR decomposition with tolerance `qr_tolerance`. + pub(crate) fn condense( + &self, + estimate: Solution, + prune_threshold: f64, + qr_tolerance: f64, + ) -> Result { + let (mut theta, mut psi, lambda, _) = estimate.into_parts(); + + let max_lambda = lambda.iter().fold(f64::NEG_INFINITY, f64::max); + let keep: Vec = lambda + .iter() + .enumerate() + .filter(|&(_, lam)| lam > max_lambda * prune_threshold) + .map(|(i, _)| i) + .collect(); + retain_support_points(&mut theta, &mut psi, &keep, "Lambda (max/1000)"); + + let (r, perm) = qr::qrd(&psi)?; + let rank_bound = psi.matrix().ncols().min(psi.matrix().nrows()); + let keep: Vec = (0..rank_bound) + .filter(|&i| (r.get(i, i) / r.col(i).norm_l2()).abs() >= qr_tolerance) + .map(|i| perm[i]) + .collect(); + retain_support_points(&mut theta, &mut psi, &keep, "QR decomposition"); + + self.check_zero_probability_subjects(&theta, &psi)?; + Solution::optimal(theta, psi).context("IPM failed during condensation") + } + + /// Log the objective function and error models at the end of a cycle, and warn if the + /// objective got worse. + fn log_cycle_summary(&self) { + tracing::info!("Objective function = {:.4}", self.objf()); + tracing::debug!("Support points: {}", self.solution.theta().nspp()); + + for (outeq, em) in self.error_models.iter() { + if *em != AssayErrorModel::None { + tracing::debug!( + "Error model for outeq {}: {:.4}", + outeq, + em.factor().unwrap_or_default() + ); + } + } + + // A falling log-likelihood signals instability or model misspecification. + if let Some(previous) = self.previous_log_likelihood { + if previous > self.solution.log_likelihood() + 1e-4 { + tracing::warn!( + "Objective function decreased from {:.4} to {:.4} (delta = {})", + -2.0 * previous, + self.objf(), + -2.0 * previous - self.objf() + ); + } + } + } + + /// Append the current cycle to the cycle log. + fn record_cycle(&mut self) { + let delta = self.objective_change().unwrap_or(f64::INFINITY); + let theta = self.solution.theta(); + self.cycle_log.push(NPCycle::new( + self.cycle, + self.objf(), + self.error_models.clone(), + theta.clone(), + self.solution.weights().clone(), + theta.nspp(), + delta, + self.status.clone(), + )); + self.previous_log_likelihood = Some(self.solution.log_likelihood()); + } + + fn into_result(self) -> Result> { + let objf = self.objf(); + let (theta, psi, weights, _) = self.solution.into_parts(); + NonParametricResult::new( + self.equation, + self.data, + self.error_models, + self.prior, + theta, + psi, + weights, + objf, + self.cycle, + self.status, + self.cycle_log, + ) + } + + /// Fail if any subject has zero (or non-finite) probability at the support points `theta`. + /// + /// Each row of `psi` holds a subject's likelihood across every support point, so a + /// subject's probability is the sum of its row. When any subject is flagged, per-subject + /// diagnostics are logged and an error naming the subjects is returned. + fn check_zero_probability_subjects(&self, theta: &Theta, psi: &Psi) -> Result<()> { + let psi = psi.matrix(); + + // Non-finite entries 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() + ); + } + + let subjects = self.data.subjects(); + let flagged: Vec = (0..psi.nrows()) + .filter(|&i| { + let probability: f64 = psi.row(i).iter().sum(); + !probability.is_finite() || probability == 0.0 + }) + .collect(); + + if flagged.is_empty() { + return Ok(()); + } + + for &i in &flagged { + self.log_zero_probability_subject(subjects[i], theta); + } + let ids = flagged.iter().map(|&i| subjects[i].id().clone()).collect(); + Err(ZeroProbabilityError::new(ids, psi.nrows()).into()) + } + + /// Log likelihood diagnostics for a subject with zero probability at the support points + /// `theta`. + fn log_zero_probability_subject(&self, subject: &Subject, theta: &Theta) { + tracing::debug!("Subject with zero probability: {}", subject.id()); + + // Simulate every support point for this subject. Points that fail to simulate count + // as having no likelihood. + let mut results: Vec<_> = (0..theta.nspp()) + .into_par_iter() + .map(|i| { + let support_point: Vec = theta.matrix().row(i).iter().copied().collect(); + match self + .equation + .simulate(subject, &support_point, Some(&self.error_models)) + { + Ok((pred, ll)) => (i, support_point, pred.get_predictions(), ll), + Err(_) => (i, support_point, Vec::new(), None), + } + }) + .collect(); + + let (mut nan, mut pos_inf, mut neg_inf, mut zero, mut valid) = (0, 0, 0, 0, 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!("====================="); + } +} + +/// Keep only the support points at `keep`, in `theta` and `psi` alike. +fn retain_support_points(theta: &mut Theta, psi: &mut Psi, keep: &[usize], step: &str) { + let dropped = theta.nspp() - keep.len(); + if dropped > 0 { + tracing::debug!("{} dropped {} support point(s)", step, dropped); + } + theta.filter_indices(keep); + psi.filter_column_indices(keep); +} + +/// A non-parametric fit in progress: the shared [`FitState`] plus the [`Runner`] that advances it. +pub(crate) struct Fit { + state: FitState, + runner: Box>, + max_cycles: usize, + watch_stop_file: bool, +} + +impl Fit { + /// Prepare a fit that stops on convergence or after `max_cycles` cycles. + pub(crate) fn new(state: FitState, runner: Box>, max_cycles: usize) -> Self { + Self { + state, + runner, + max_cycles, + watch_stop_file: false, + } + } + + /// Also stop after any cycle at which a [`STOP_FILE`] exists in the working directory. + /// + /// A stop file left over from an earlier run is removed now, so it can't end this one + /// before it starts. + pub(crate) fn watch_stop_file(mut self) -> Result { + if Path::new(STOP_FILE).exists() { + tracing::info!("Removing existing stop file prior to run"); + fs::remove_file(STOP_FILE).context("Unable to remove previous stop file")?; + } + self.watch_stop_file = true; + Ok(self) + } + + pub(crate) fn state(&self) -> &FitState { + &self.state + } + + /// Run the next cycle, record it in the cycle log, and return the new [`Status`]. + /// + /// Once the fit has stopped, this does nothing and returns the final status. + pub(crate) fn step(&mut self) -> Result { + if self.state.status.is_stop() { + return Ok(self.state.status.clone()); + } + self.state.cycle += 1; + let span = tracing::info_span!("", "{}", format!("Cycle {}", self.state.cycle)); + let _enter = span.enter(); + + let outcome = self.runner.step(&mut self.state)?; + self.state.log_cycle_summary(); + + self.state.status = match outcome { + CycleOutcome::Converged => { + tracing::info!("The model converged after {} cycles", self.state.cycle); + Status::Stop(StopReason::Converged) + } + CycleOutcome::Continue => self.limit_reached().map_or(Status::Continue, Status::Stop), + }; + self.state.record_cycle(); + Ok(self.state.status.clone()) + } + + /// Why the fit must stop although it has not converged, if it must. + fn limit_reached(&self) -> Option { + if self.state.cycle >= self.max_cycles { + tracing::warn!("Maximum number of cycles reached"); + return Some(StopReason::MaxCycles); + } + if self.watch_stop_file && Path::new(STOP_FILE).exists() { + tracing::warn!("Stopfile detected - breaking"); + return Some(StopReason::StopFile); + } + None + } + + /// Run cycles until the fit stops. Does nothing if it already has. + pub(crate) fn run(&mut self) -> Result<()> { + while self.state.status.is_continue() { + self.step()?; + } + Ok(()) + } + + /// Stop a running fit as [`Aborted`](StopReason::Aborted). A fit that has already stopped + /// keeps its reason. + pub(crate) fn request_stop(&mut self) { + if self.state.status.is_continue() { + self.state.status = Status::Stop(StopReason::Aborted); + } + } + + pub(crate) fn into_result(self) -> Result> { + self.state.into_result() + } +} diff --git a/src/algorithms/nonparametric/mod.rs b/src/algorithms/nonparametric/mod.rs index ef82222f9..ab65b2151 100644 --- a/src/algorithms/nonparametric/mod.rs +++ b/src/algorithms/nonparametric/mod.rs @@ -6,21 +6,37 @@ //! //! # Available Algorithms //! -//! - [`NPAG`](npag): Non-Parametric Adaptive Grid -//! - [`NPOD`](npod): Non-Parametric Optimal Design -//! - [`NPMAP`](npmap): Maximum a posteriori reweighting +//! - [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 //! //! Use the [`NonParametricAlgorithm`] enum to select and configure an algorithm. Each -//! variant wraps its algorithm-specific configuration struct (e.g. [`NpagConfig`]). The -//! internal execution state used while fitting implements the [`NonParametricRunner`] -//! trait, which defines the common interface for initialization, estimation, condensation, -//! expansion, and convergence evaluation. +//! variant wraps its algorithm-specific configuration struct (e.g. [`NpagConfig`]). +//! +//! # Structure +//! +//! Every algorithm runs on the same engine (see the `fit` module): the state shared by all +//! algorithms (support points, likelihoods, weights, cycle log) lives in one place, and each +//! algorithm only supplies what a single cycle does to it. +//! +//! # Stopping a running fit +//! +//! A fit started through [`fit_with`](crate::estimation::EstimationProblem::fit_with) or +//! [`fit_controller`](crate::estimation::EstimationProblem::fit_controller) stops after +//! the current cycle if a file named `stop` appears in the working directory, with +//! [`StopReason::StopFile`](crate::algorithms::StopReason::StopFile). A `stop` file left over +//! from an earlier run is removed when the fit starts. // Shared error-model factor optimization pub mod error_optim; +// The engine every algorithm runs on. +mod fit; +mod solution; + // Algorithm implementations pub mod ncnpag; pub mod npag; @@ -30,12 +46,6 @@ pub mod npod; // Incremental, observable fitting (stepping handle + per-cycle observers). pub mod controller; -// Re-export algorithm structs -pub use ncnpag::NCNPAG; -pub use npag::NPAG; -pub use npmap::NPMAP; -pub use npod::NPOD; - // Re-export per-algorithm configuration structs pub use error_optim::ErrorOptimConfig; pub use ncnpag::NcnpagConfig; @@ -44,40 +54,50 @@ pub use npmap::NpmapConfig; pub use npod::NpodConfig; // Re-export the incremental fitting API -pub use controller::{CycleFlow, FitController, FitObserver}; +pub use controller::FitController; -use crate::algorithms::{Algorithm, NonParametricRunner}; +use crate::algorithms::Algorithm; use crate::estimation::nonparametric::NonParametricResult; use crate::estimation::{EstimationProblem, NonParametric}; use anyhow::Result; +use fit::{Fit, FitState, Runner}; +use ncnpag::Ncnpag; +use npag::Npag; +use npmap::Npmap; +use npod::Npod; use pharmsol::prelude::simulator::Equation; /// The non-parametric algorithms supported by PMcore. /// -/// Use the constructors to select an algorithm with its default configuration: +/// To run an algorithm, pass its configuration struct straight to +/// [`fit_with`](crate::estimation::EstimationProblem::fit_with) or +/// [`fit_controller`](crate::estimation::EstimationProblem::fit_controller). Each one exposes +/// only the settings valid for its algorithm: /// /// ```no_run /// use pmcore::prelude::*; -/// -/// // Default NPAG configuration. -/// let algorithm = NonParametricAlgorithm::npag(); +/// # fn run(problem: EstimationProblem) -> Result<()> { +/// let result = problem.fit_with(NpagConfig::new().eps(0.1).max_cycles(500))?; +/// # let _ = result; Ok(()) +/// # } /// ``` /// -/// To customize an algorithm, build its configuration struct (which exposes only the -/// setters valid for that algorithm) and pass it directly to -/// [`fit_with`](crate::estimation::EstimationProblem::fit_with): +/// Use this enum when the algorithm is only known at runtime, for example when it comes from +/// a settings file. Every configuration converts into it, and it can be passed to the same +/// functions: /// /// ```no_run /// use pmcore::prelude::*; -/// -/// // NPAG with a tighter convergence criterion and a cycle cap. -/// let config = NpagConfig::new().eps(0.1).max_cycles(500); -/// // `problem.fit_with(config)` accepts the config directly. +/// # fn run(problem: EstimationProblem, use_npod: bool) -> Result<()> { +/// let algorithm: NonParametricAlgorithm = if use_npod { +/// NpodConfig::new().into() +/// } else { +/// NpagConfig::new().into() +/// }; +/// let result = problem.fit_with(algorithm)?; +/// # let _ = result; Ok(()) +/// # } /// ``` -/// -/// 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. #[derive(Debug, Clone)] pub enum NonParametricAlgorithm { /// Non-Parametric Adaptive Grid. @@ -142,91 +162,48 @@ impl NonParametricAlgorithm { Self::Ncnpag(NcnpagConfig::default()) } - /// Build the internal, mutable execution state (runner) for this algorithm. + /// Set up a fit of `problem` with this algorithm, ready to run its first cycle. /// /// Both [`fit`](Algorithm::fit) and the stepping /// [`FitController`](crate::algorithms::nonparametric::controller::FitController) - /// build on this primitive. - pub(crate) fn into_runner( + /// build on this. + pub(crate) fn start( self, problem: EstimationProblem, - ) -> 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 { - Self::Npag(config) => Box::new(NPAG::from_parts( - problem.model.equation, - problem.data, - problem.error_models, - problem.prior, - config, - )?), - Self::Npod(config) => Box::new(NPOD::from_parts( - problem.model.equation, - problem.data, - problem.error_models, - problem.prior, - config, - )?), - Self::Npmap(config) => Box::new(NPMAP::from_parts( - problem.model.equation, - problem.data, - problem.error_models, - problem.prior, - config, - )?), - Self::Ncnpag(config) => Box::new(NCNPAG::from_parts( - problem.model.equation, - problem.data, - problem.error_models, - problem.prior, - config, - )?), + ) -> Result> { + let state = FitState::new(problem); + let (runner, max_cycles): (Box>, usize) = match self { + Self::Npag(config) => { + let max_cycles = config.max_cycles; + let runner = Npag::new(config, state.solution().theta(), state.error_models()); + (Box::new(runner), max_cycles) + } + Self::Npod(config) => { + let max_cycles = config.max_cycles; + ( + Box::new(Npod::new(config, state.error_models())), + max_cycles, + ) + } + // Single-pass algorithms. + Self::Npmap(_) => (Box::new(Npmap), 1), + Self::Ncnpag(config) => (Box::new(Ncnpag::new(config)), 1), }; - Ok(runner) - } -} - -impl Algorithm for NonParametricAlgorithm { - type Output = NonParametricResult; - - fn fit(self, problem: EstimationProblem) -> Result { - let mut runner = self.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) + Fit::new(state, runner, max_cycles).watch_stop_file() } } -impl Algorithm for NcnpagConfig { +/// Every non-parametric algorithm (or its configuration) can be fitted directly. +impl Algorithm for A +where + E: Equation, + A: Into, +{ type Output = NonParametricResult; fn fit(self, problem: EstimationProblem) -> Result { - NonParametricAlgorithm::from(self).fit(problem) + let mut fit = self.into().start(problem)?; + fit.run()?; + fit.into_result() } } diff --git a/src/algorithms/nonparametric/ncnpag.rs b/src/algorithms/nonparametric/ncnpag.rs index 860968ce7..4360e7a25 100644 --- a/src/algorithms/nonparametric/ncnpag.rs +++ b/src/algorithms/nonparametric/ncnpag.rs @@ -1,22 +1,13 @@ -use crate::{ - algorithms::{ - nonparametric::{npag::NPAG, NpagConfig}, - NonParametricRunner, Status, StopReason, - }, - estimation::nonparametric::{ - calculate_psi, CycleLog, NPCycle, NonParametricResult, Psi, Theta, Weights, - }, -}; - use anyhow::Result; use faer::Mat; -use pharmsol::prelude::{ - data::{AssayErrorModels, Data}, - simulator::Equation, -}; - +use pharmsol::prelude::simulator::Equation; use serde::{Deserialize, Serialize}; +use super::fit::{CycleOutcome, Fit, FitState, Runner}; +use super::npag::{Npag, NpagConfig}; +use super::solution::Solution; +use crate::estimation::nonparametric::{Theta, Weights}; + /// Configuration options for the non-collapsing NPAG (NCNPAG) algorithm. #[derive(Debug, Clone, PartialEq, Serialize, Deserialize)] pub struct NcnpagConfig { @@ -73,55 +64,22 @@ impl NcnpagConfig { /// /// The result is returned as a standard [`NonParametricResult`], so the /// `(theta, weights)` can be consumed exactly like any other fit. -pub struct NCNPAG { - equation: E, - psi: Psi, - theta: Theta, - w: Weights, - objf: f64, - cycle: usize, - status: Status, - data: Data, - cyclelog: CycleLog, - error_models: AssayErrorModels, - prior: Theta, - cycles: usize, - progress: bool, +#[derive(Debug)] +pub(crate) struct Ncnpag { + config: NcnpagConfig, } -impl NCNPAG { - pub(crate) fn from_parts( - equation: E, - data: Data, - error_models: AssayErrorModels, - theta: Theta, - config: NcnpagConfig, - ) -> Result { - Ok(Self { - equation, - psi: Psi::new(), - theta: theta.clone(), - w: Weights::default(), - objf: f64::INFINITY, - cycle: 0, - status: Status::Continue, - data, - cyclelog: CycleLog::new(), - error_models, - prior: theta, - cycles: config.cycles, - progress: config.progress, - }) +impl Ncnpag { + pub(crate) fn new(config: NcnpagConfig) -> Self { + Self { config } } } /// Refine each support point with a full NPAG seeded from that single point, /// preserving the point's filter weight (NPAGFULL). Points whose refinement /// fails or produces no output are kept at their original location. -fn refine_points( - equation: &E, - data: &Data, - error_models: &AssayErrorModels, +fn refine_points( + state: &FitState, theta: &Theta, weights: &Weights, cycles: usize, @@ -142,31 +100,23 @@ fn refine_points( progress, ..Default::default() }; - let mut npag = NPAG::from_parts( - equation.clone(), - data.clone(), - error_models.clone(), + let npag = Npag::new(npag_config, &single_theta, state.error_models()); + let nested = FitState::from_parts( + state.equation().clone(), + state.data().clone(), + state.error_models().clone(), single_theta, - npag_config, - )?; - - #[allow(clippy::while_let_loop)] - let run = npag.initialize().and_then(|_| { - loop { - match npag.next_cycle()? { - Status::Continue => continue, - Status::Stop(_) => break, - } - } - Ok(()) - }); + ); + let mut fit = Fit::new(nested, Box::new(npag), cycles); + let run = fit.run().map(|()| fit); match run { - Ok(()) if npag.theta().matrix().nrows() > 0 => { - let refined: Vec = npag.theta().matrix().row(0).iter().copied().collect(); + Ok(fit) if fit.state().solution().theta().nspp() > 0 => { + let theta = fit.state().solution().theta(); + let refined: Vec = theta.matrix().row(0).iter().copied().collect(); refined_points.push(refined); } - Ok(()) => { + Ok(_) => { tracing::warn!( "NCNPAG: refinement produced no points for support point {} — keeping original", i + 1 @@ -199,96 +149,14 @@ fn refine_points( Ok((refined_theta, refined_weights)) } -/// Marginal log-likelihood of the data under a discrete `(psi, weights)` model. -fn marginal_loglik(psi: &Psi, w: &Weights) -> f64 { - let m = psi.matrix(); - (0..m.nrows()) - .map(|s| { - let acc: f64 = (0..m.ncols()).map(|j| *m.get(s, j) * w[j]).sum(); - acc.max(f64::MIN_POSITIVE).ln() - }) - .sum::() - + 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 - } - - fn equation(&self) -> &E { - &self.equation - } - - fn data(&self) -> &Data { - &self.data - } - - fn likelihood(&self) -> f64 { - self.objf - } - - fn increment_cycle(&mut self) -> usize { - 0 - } - - fn cycle(&self) -> usize { - 0 - } - - fn set_theta(&mut self, theta: Theta) { - self.theta = theta; - } - - fn theta(&self) -> &Theta { - &self.theta - } - - fn psi(&self) -> &Psi { - &self.psi - } - - fn set_status(&mut self, status: Status) { - self.status = status; - } - - fn status(&self) -> &Status { - &self.status - } - - fn evaluation(&mut self) -> Result { - self.status = Status::Stop(StopReason::Converged); - Ok(self.status.clone()) - } - - fn estimation(&mut self) -> Result<()> { +impl Runner for Ncnpag { + fn step(&mut self, state: &mut FitState) -> Result { // Likelihood of each fixed support point for the data. - let psi = calculate_psi( - &self.equation, - &self.data, - &self.theta, - &self.error_models, - false, - )?; + let mut theta = state.solution().theta().clone(); + let psi = state.compute_psi(&theta, false)?; // Flat (uniform) prior: postⱼ ∝ ∏ᵢ P(dataᵢ | θⱼ). Accumulate in log space. - let n_points = self.theta.matrix().nrows(); + let n_points = theta.nspp(); let mut log_weights = vec![f64::NEG_INFINITY; n_points]; for (j, slot) in log_weights.iter_mut().enumerate() { let mut log_weight = 0.0; // ln(uniform prior) is constant, drops out on normalization @@ -343,92 +211,25 @@ impl NonParametricRunner for NCNPAG { .collect(); // Filter theta and renormalize the surviving weights (NPAGFULL11). - self.theta.filter_indices(&keep); + theta.filter_indices(&keep); let kept: Vec = keep.iter().map(|&i| weights[i]).collect(); let sum: f64 = kept.iter().sum(); - self.w = Weights::from_vec(kept.iter().map(|w| w / sum).collect()); + let mut weights = 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 { - let (refined_theta, refined_weights) = refine_points( - &self.equation, - &self.data, - &self.error_models, - &self.theta, - &self.w, - self.cycles, - self.progress, + if self.config.cycles > 0 { + (theta, weights) = refine_points( + state, + &theta, + &weights, + self.config.cycles, + self.config.progress, )?; - self.theta = refined_theta; - self.w = refined_weights; } // Recompute psi over the final support points so psi/weights stay aligned. - self.psi = calculate_psi( - &self.equation, - &self.data, - &self.theta, - &self.error_models, - false, - )?; - - self.objf = 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)] -mod tests { - use super::*; - use approx::assert_relative_eq; - use ndarray::Array2; - - #[test] - fn marginal_loglik_restores_subject_scaling() -> anyhow::Result<()> { - let log_likelihoods = Array2::from_shape_vec((2, 2), vec![-1000.0, -1001.0, -2.0, -4.0])?; - let psi = Psi::from_log_likelihoods(log_likelihoods)?; - let weights = Weights::from_vec(vec![0.6, 0.4]); - - let expected = -1000.0 + (0.6 + 0.4 * (-1.0_f64).exp()).ln() - 2.0 - + (0.6 + 0.4 * (-2.0_f64).exp()).ln(); - - assert_relative_eq!(marginal_loglik(&psi, &weights), expected, epsilon = 1e-12); - Ok(()) + let psi = state.compute_psi(&theta, false)?; + state.set_solution(Solution::with_weights(theta, psi, weights)); + Ok(CycleOutcome::Converged) } } diff --git a/src/algorithms/nonparametric/npag.rs b/src/algorithms/nonparametric/npag.rs index 097b9e2d9..1c4603045 100644 --- a/src/algorithms/nonparametric/npag.rs +++ b/src/algorithms/nonparametric/npag.rs @@ -1,26 +1,11 @@ -use crate::algorithms::{NonParametricRunner, Status, StopReason}; -use crate::estimation::nonparametric::{ - calculate_psi, 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 pharmsol::prelude::{data::AssayErrorModels, simulator::Equation}; use serde::{Deserialize, Serialize}; +use super::error_optim::{ErrorModelOptimizer, ErrorOptimConfig}; +use super::fit::{CycleOutcome, FitState, Runner}; +use crate::estimation::nonparametric::{adaptative_grid, Theta}; + /// Configuration options for the Non-Parametric Adaptive Grid (NPAG) algorithm. #[derive(Debug, Clone, PartialEq, Serialize, Deserialize)] pub struct NpagConfig { @@ -109,332 +94,96 @@ impl NpagConfig { } } +/// NPAG's own state between cycles: its configuration, the adaptive grid resolution, and the +/// bookkeeping for its convergence test. #[derive(Debug)] -pub struct NPAG { - equation: E, +pub(crate) struct Npag { + config: NpagConfig, + /// Bounds of each parameter; the grid is refined within them. ranges: Vec<(f64, f64)>, - psi: Psi, - prior: Theta, - theta: Theta, - lambda: Weights, - w: Weights, + /// Grid resolution, as a fraction of each parameter's range. eps: f64, - last_objf: f64, - objf: f64, - f0: f64, - f1: f64, - cycle: usize, - gamma_delta: Vec, - error_models: AssayErrorModels, - status: Status, - cycle_log: CycleLog, - data: Data, - config: NpagConfig, + /// `Σ ln(Ψ·w)` the last time `eps` was refined down to `min_eps`. NPAG has converged when + /// two such values in a row agree to within `pyl_tolerance`. + last_fine_objective: f64, + error_optimizer: ErrorModelOptimizer, } -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`]. - pub(crate) fn from_parts( - equation: E, - data: Data, - error_models: AssayErrorModels, - theta: Theta, - config: NpagConfig, - ) -> Result { - let ranges = theta.parameters().finite_ranges(); - let gamma_delta = vec![config.error_optim.step; error_models.len()]; - let eps = config.eps; - - Ok(Self { - equation, - ranges, - psi: Psi::new(), - prior: theta.clone(), - theta, - lambda: Weights::default(), - w: Weights::default(), - eps, - last_objf: -1e30, - objf: f64::NEG_INFINITY, - f0: -1e30, - f1: f64::default(), - cycle: 0, - gamma_delta, - error_models, - status: Status::Continue, - cycle_log: CycleLog::new(), - data, +impl Npag { + pub(crate) fn new(config: NpagConfig, prior: &Theta, error_models: &AssayErrorModels) -> Self { + Self { + ranges: prior.parameters().finite_ranges(), + eps: config.eps, + last_fine_objective: -1e30, + error_optimizer: ErrorModelOptimizer::new(config.error_optim.clone(), error_models), config, - }) - } -} - -impl NonParametricRunner 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 error_models(&self) -> &AssayErrorModels { - &self.error_models - } - - fn data(&self) -> &Data { - &self.data - } - - fn likelihood(&self) -> f64 { - self.objf - } - - fn increment_cycle(&mut self) -> usize { - self.cycle += 1; - self.cycle - } - - fn cycle(&self) -> usize { - self.cycle - } - - fn set_theta(&mut self, theta: Theta) { - self.theta = theta; - } - - fn theta(&self) -> &Theta { - &self.theta - } - - fn psi(&self) -> &Psi { - &self.psi + } } - 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 {}: {:.2}", - outeq, - em.factor().unwrap_or_default() - ); - }); - + /// Whether the fit has converged. + /// + /// Every time the objective stalls, the grid is refined by halving `eps`. Once `eps` + /// reaches `min_eps`, the fit has converged if `Σ ln(Ψ·w)` matches its value at the + /// previous such point; otherwise `eps` is reset to start another round of refinement. + fn converged(&mut self, state: &FitState) -> bool { 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 - && self.eps > self.config.min_eps - { - self.eps /= 2.; - if self.eps <= self.config.min_eps { - 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()); - } else { - self.f0 = self.f1; - self.eps = self.config.eps; - } - } + let stalled = state + .objective_change() + .is_some_and(|change| change <= self.config.objective_tolerance); + if !stalled || self.eps <= self.config.min_eps { + return false; } - // 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()); + self.eps /= 2.0; + if self.eps > self.config.min_eps { + return false; } - // Stop if stopfile exists - if std::path::Path::new("stop").exists() { - tracing::warn!("Stopfile detected - breaking"); - self.set_status(Status::Stop(StopReason::StopFile)); - self.log_cycle_state(); - return Ok(self.status().clone()); + let solution = state.solution(); + let pyl = solution.psi().matrix() * solution.weights().weights(); + let objective: f64 = pyl.iter().map(|x| x.ln()).sum(); + if (objective - self.last_fine_objective).abs() <= self.config.pyl_tolerance { + return true; } - - // Continue with normal operation - self.set_status(Status::Continue); - self.log_cycle_state(); - Ok(self.status().clone()) + self.last_fine_objective = objective; + self.eps = self.config.eps; + false } +} - fn estimation(&mut self) -> Result<()> { - self.psi = calculate_psi( - &self.equation, - &self.data, - &self.theta, - &self.error_models, - self.cycle == 1 && self.config.progress, - )?; - - if let Err(err) = self.check_zero_probability_subjects() { - bail!(err); - } - - (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(), - ); +impl Runner for Npag { + fn step(&mut self, state: &mut FitState) -> Result { + let mut theta = state.solution().theta().clone(); + if !state.is_first_cycle() { + adaptative_grid( + &mut theta, + self.eps, + &self.ranges, + self.config.grid_tolerance, + )?; } - // 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.w = self.lambda.clone(); - Ok(()) - } - - fn optimizations(&mut self) -> Result<()> { - optimize_error_models( - &self.equation, - &self.data, - &self.theta, - &mut self.error_models, - &mut self.gamma_delta, - &mut self.objf, - &mut self.lambda, - &mut self.psi, - &self.config.error_optim, - ) - } - - fn expansion(&mut self) -> Result<()> { - adaptative_grid( - &mut self.theta, - self.eps, - &self.ranges, - self.config.grid_tolerance, + let estimate = state.estimate(theta, self.config.progress && state.is_first_cycle())?; + let solution = state.condense( + estimate, + self.config.prune_threshold, + self.config.qr_tolerance, )?; - Ok(()) - } - - fn set_status(&mut self, status: Status) { - self.status = status; - } + state.set_solution(solution); + self.error_optimizer.optimize(state)?; - 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; + Ok(if self.converged(state) { + CycleOutcome::Converged + } else { + CycleOutcome::Continue + }) } } #[cfg(test)] mod tests { + use crate::algorithms::{Status, StopReason}; use crate::prelude::*; use pharmsol::{fa, fetch_params, lag, Subject, SubjectBuilderExt}; @@ -492,7 +241,7 @@ mod tests { EstimationProblem::nonparametric(simple_equation(), simple_data(), prior, error_models) .expect("Failed to build problem"); - let result = problem.fit_with(NonParametricAlgorithm::npag()); + let result = problem.fit_with(NpagConfig::new()); assert!( result.is_ok(), @@ -500,4 +249,109 @@ mod tests { result.err() ); } + + #[test] + fn controller_exposes_a_consistent_solution_between_steps() -> anyhow::Result<()> { + let parameters = ParameterSpace::bounded() + .add("ke", 0.001, 3.0) + .add("v", 25.0, 250.0); + let error_models = AssayErrorModels::new().add( + "0", + AssayErrorModel::additive(ErrorPoly::new(0.0, 0.5, 0.0, 0.0), 0.0), + )?; + let problem = EstimationProblem::nonparametric( + simple_equation(), + simple_data(), + Theta::sobol_default(¶meters)?, + error_models, + )?; + + let mut controller = problem.fit_controller(NpagConfig::new().progress(false))?; + for _ in 0..3 { + controller.step()?; + let nspp = controller.theta().nspp(); + assert_eq!(controller.weights().len(), nspp); + assert_eq!(controller.psi().matrix().ncols(), nspp); + assert_eq!(controller.objf(), -2.0 * controller.log_likelihood()); + } + assert_eq!(controller.cycle(), 3); + Ok(()) + } + + fn problem() -> anyhow::Result> { + let parameters = ParameterSpace::bounded() + .add("ke", 0.001, 3.0) + .add("v", 25.0, 250.0); + let error_models = AssayErrorModels::new().add( + "0", + AssayErrorModel::additive(ErrorPoly::new(0.0, 0.5, 0.0, 0.0), 0.0), + )?; + EstimationProblem::nonparametric( + simple_equation(), + simple_data(), + Theta::sobol_default(¶meters)?, + error_models, + ) + } + + #[test] + fn a_stopped_fit_stays_stopped() -> anyhow::Result<()> { + let mut controller = problem()?.fit_controller(NpagConfig::new().progress(false))?; + controller.step()?; + controller.request_stop(); + let aborted = Status::Stop(StopReason::Aborted); + + assert_eq!(controller.step()?, aborted); + assert_eq!(controller.cycle(), 1); + + let result = controller.finish()?; + assert_eq!(result.cycles(), 1); + assert_eq!(result.cycle_log().cycles().len(), 1); + assert_eq!(result.stop_reason(), Some(&StopReason::Aborted)); + Ok(()) + } + + #[test] + fn request_stop_keeps_the_reason_of_a_finished_fit() -> anyhow::Result<()> { + let mut controller = + problem()?.fit_controller(NpagConfig::new().progress(false).max_cycles(2))?; + while controller.step()?.is_continue() {} + controller.request_stop(); + + let result = controller.into_result()?; + assert_eq!(result.status(), &Status::Stop(StopReason::MaxCycles)); + assert_eq!(result.cycles(), 2); + Ok(()) + } + + #[test] + fn run_until_sees_every_cycle_and_stops_when_asked() -> anyhow::Result<()> { + let mut seen = Vec::new(); + let result = problem()? + .fit_controller(NpagConfig::new().progress(false))? + .run_until(|fit| { + seen.push(fit.cycle()); + fit.cycle() >= 5 + })?; + + assert_eq!(seen, vec![1, 2, 3, 4, 5]); + assert_eq!(result.cycles(), 5); + assert_eq!(result.stop_reason(), Some(&StopReason::Aborted)); + Ok(()) + } + + #[test] + fn run_until_ends_with_the_algorithms_own_reason() -> anyhow::Result<()> { + let mut calls = 0; + let result = problem()? + .fit_controller(NpagConfig::new().progress(false).max_cycles(3))? + .run_until(|_| { + calls += 1; + false + })?; + + assert_eq!(calls, 3); + assert_eq!(result.stop_reason(), Some(&StopReason::MaxCycles)); + Ok(()) + } } diff --git a/src/algorithms/nonparametric/npmap.rs b/src/algorithms/nonparametric/npmap.rs index c5765a6a1..f24d2325a 100644 --- a/src/algorithms/nonparametric/npmap.rs +++ b/src/algorithms/nonparametric/npmap.rs @@ -1,19 +1,9 @@ -use crate::{ - algorithms::{NonParametricRunner, Status, StopReason}, - estimation::nonparametric::{ - calculate_psi, CycleLog, NPCycle, NonParametricResult, Psi, Theta, Weights, - }, -}; - -use anyhow::{Context, Result}; -use pharmsol::prelude::{ - data::{AssayErrorModels, Data}, - simulator::Equation, -}; - -use crate::estimation::nonparametric::ipm::burke; +use anyhow::Result; +use pharmsol::prelude::simulator::Equation; use serde::{Deserialize, Serialize}; +use super::fit::{CycleOutcome, FitState, Runner}; + /// Configuration options for the non-parametric maximum a posteriori (NPMAP) algorithm #[derive(Debug, Clone, PartialEq, Default, Serialize, Deserialize)] pub struct NpmapConfig {} @@ -24,161 +14,81 @@ impl NpmapConfig { } } -/// Non-parametric maximum a posteriori (NPMAP) algorithm +/// Non-parametric maximum a posteriori (NPMAP) algorithm. /// -/// This algorithm is a wrapper around the IPM algorithm that calculates the posterior probabilities of the support points -/// given a prior distribution and the likelihood of the data. +/// Reweights the prior support points by their likelihood given the data, in a single pass: +/// the support points themselves are left unchanged. #[derive(Debug)] -pub struct NPMAP { - equation: E, - psi: Psi, - theta: Theta, - w: Weights, - objf: f64, - cycle: usize, - status: Status, - data: Data, - cyclelog: CycleLog, - error_models: AssayErrorModels, - prior: Theta, -} +pub(crate) struct Npmap; -impl NPMAP { - pub(crate) fn from_parts( - equation: E, - data: Data, - error_models: AssayErrorModels, - theta: Theta, - _config: NpmapConfig, - ) -> Result { - Ok(Self { - equation, - psi: Psi::new(), - theta: theta.clone(), - w: Weights::default(), - objf: f64::INFINITY, - cycle: 0, - status: Status::Continue, - data, - cyclelog: CycleLog::new(), - error_models, - prior: theta, - }) +impl Runner for Npmap { + fn step(&mut self, state: &mut FitState) -> Result { + let solution = state.estimate(state.solution().theta().clone(), false)?; + state.set_solution(solution); + Ok(CycleOutcome::Converged) } } -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(), +#[cfg(test)] +mod tests { + use crate::estimation::nonparametric::burke; + use crate::prelude::*; + + use pharmsol::{fa, fetch_params, lag, Subject, SubjectBuilderExt}; + + fn equation() -> pharmsol::equation::ODE { + pharmsol::equation::ODE::new( + |x, p, _t, dx, b, _rateiv, _cov| { + fetch_params!(p, ke); + dx[0] = -ke * x[0] + b[0]; + }, + |_p, _t, _cov| lag! {}, + |_p, _t, _cov| fa! {}, + |_p, _t, _cov, _x| {}, + |x, p, _t, _cov, y| { + fetch_params!(p, v); + y[0] = x[0] / v; + }, ) - } - - fn error_models(&self) -> &AssayErrorModels { - &self.error_models - } - - fn equation(&self) -> &E { - &self.equation - } - - fn data(&self) -> &Data { - &self.data - } - - fn likelihood(&self) -> f64 { - self.objf - } - - fn increment_cycle(&mut self) -> usize { - 0 - } - - fn cycle(&self) -> usize { - 0 - } - - fn set_theta(&mut self, theta: Theta) { - self.theta = theta; - } - - fn theta(&self) -> &Theta { - &self.theta - } - - fn psi(&self) -> &Psi { - &self.psi - } - - fn set_status(&mut self, status: Status) { - self.status = status; - } - - fn status(&self) -> &Status { - &self.status - } - - fn evaluation(&mut self) -> Result { - self.status = Status::Stop(StopReason::Converged); - Ok(self.status.clone()) - } - - fn estimation(&mut self) -> Result<()> { - self.psi = calculate_psi( - &self.equation, - &self.data, - &self.theta, - &self.error_models, - false, + .with_nstates(1) + .with_ndrugs(1) + .with_nout(1) + .with_metadata( + pharmsol::equation::metadata::new("npmap_test") + .parameters(["ke", "v"]) + .states(["central"]) + .outputs(["0"]) + .route(pharmsol::equation::Route::bolus("0").to_state("central")), + ) + .expect("metadata attachment should validate") + } + + #[test] + fn npmap_reports_neg2ll_like_every_other_algorithm() -> anyhow::Result<()> { + let data = Data::new(vec![Subject::builder("1") + .bolus(0.0, 100.0, 0) + .observation(1.0, 10.0, 0) + .observation(2.0, 8.0, 0) + .build()]); + let parameters = ParameterSpace::bounded() + .add("ke", 0.001, 3.0) + .add("v", 25.0, 250.0); + let error_models = AssayErrorModels::new().add( + "0", + AssayErrorModel::additive(ErrorPoly::new(0.0, 0.5, 0.0, 0.0), 0.0), + )?; + let problem = EstimationProblem::nonparametric( + equation(), + data, + Theta::sobol(¶meters, 64)?, + error_models, )?; - (self.w, self.objf) = burke(&self.psi).context("Error in IPM")?; - Ok(()) - } - - fn condensation(&mut self) -> Result<()> { - Ok(()) - } - fn optimizations(&mut self) -> Result<()> { - Ok(()) - } + let result = problem.fit_with(NpmapConfig::new())?; - fn expansion(&mut self) -> Result<()> { + let (_, log_likelihood) = burke(result.psi())?; + approx::assert_relative_eq!(result.objf(), -2.0 * log_likelihood, epsilon = 1e-10); + assert_eq!(result.cycle_log().cycles()[0].objf(), result.objf()); 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..20ba99dd5 100644 --- a/src/algorithms/nonparametric/npod.rs +++ b/src/algorithms/nonparametric/npod.rs @@ -1,30 +1,28 @@ -use crate::{ - algorithms::{NonParametricRunner, Status, StopReason}, - estimation::nonparametric::{ - calculate_psi, ipm::burke, qr, 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 ndarray::Array1; -use rayon::prelude::{IntoParallelRefMutIterator, ParallelIterator}; +use pharmsol::prelude::{data::AssayErrorModels, simulator::Equation}; +use pharmsol::ParameterOptimizer; +use rayon::prelude::*; use serde::{Deserialize, Serialize}; -use super::error_optim::{optimize_error_models, ErrorOptimConfig}; - -const THETA_F: f64 = 1e-2; -const THETA_D: f64 = 1e-4; +use super::error_optim::{ErrorModelOptimizer, ErrorOptimConfig}; +use super::fit::{CycleOutcome, FitState, Runner}; +use crate::estimation::nonparametric::Theta; /// Configuration options for the Non-Parametric Optimal Design (NPOD) algorithm. #[derive(Debug, Clone, Serialize, Deserialize)] +#[serde(default)] pub struct NpodConfig { /// Maximum number of cycles to run the algorithm for. pub max_cycles: usize, + /// Converged once the log-likelihood changes by at most this much in a cycle. + pub objective_tolerance: f64, + /// Minimum distance between a new candidate point and the existing support points. + pub min_distance: f64, + /// Support points with at most this fraction of the largest weight are dropped. + pub prune_threshold: f64, + /// Tolerance of the rank-revealing QR decomposition that drops linearly dependent points. + pub qr_tolerance: f64, /// Configuration for the error-model factor (gamma/lambda) optimization. pub error_optim: ErrorOptimConfig, /// Whether to print progress information during the first cycle. @@ -41,6 +39,26 @@ impl NpodConfig { self } + pub fn objective_tolerance(mut self, tolerance: f64) -> Self { + self.objective_tolerance = tolerance; + self + } + + pub fn min_distance(mut self, distance: f64) -> Self { + self.min_distance = distance; + self + } + + pub fn prune_threshold(mut self, threshold: f64) -> Self { + self.prune_threshold = threshold; + self + } + + pub fn qr_tolerance(mut self, tolerance: f64) -> Self { + self.qr_tolerance = tolerance; + self + } + pub fn error_optim(mut self, config: ErrorOptimConfig) -> Self { self.error_optim = config; self @@ -56,311 +74,108 @@ impl Default for NpodConfig { fn default() -> Self { Self { max_cycles: 100, + objective_tolerance: 1e-2, + min_distance: 1e-4, + prune_threshold: 1e-3, + qr_tolerance: 1e-8, error_optim: ErrorOptimConfig::default(), progress: true, } } } +/// NPOD's own state between cycles: its configuration and the error-model optimizer. #[derive(Debug)] -pub struct NPOD { - equation: E, - psi: Psi, - prior: Theta, - theta: Theta, - lambda: Weights, - w: Weights, - last_objf: f64, - objf: f64, - cycle: usize, - gamma_delta: Vec, - error_models: AssayErrorModels, - converged: bool, - status: Status, - cycle_log: CycleLog, - data: Data, +pub(crate) struct Npod { config: NpodConfig, + error_optimizer: ErrorModelOptimizer, } -impl NPOD { - pub(crate) fn from_parts( - equation: E, - data: Data, - error_models: AssayErrorModels, - theta: Theta, - config: NpodConfig, - ) -> Result { - let gamma_delta = vec![config.error_optim.step; error_models.len()]; - - Ok(Self { - equation, - psi: Psi::new(), - prior: theta.clone(), - theta, - lambda: Weights::default(), - w: Weights::default(), - last_objf: -1e30, - objf: f64::NEG_INFINITY, - cycle: 0, - gamma_delta, - error_models, - converged: false, - status: Status::Continue, - cycle_log: CycleLog::new(), - data, +impl Npod { + pub(crate) fn new(config: NpodConfig, error_models: &AssayErrorModels) -> Self { + Self { + error_optimizer: ErrorModelOptimizer::new(config.error_optim.clone(), error_models), config, - }) - } -} - -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(), - ) - } - - fn equation(&self) -> &E { - &self.equation - } - - fn error_models(&self) -> &AssayErrorModels { - &self.error_models - } - - fn data(&self) -> &Data { - &self.data - } - - fn increment_cycle(&mut self) -> usize { - self.cycle += 1; - self.cycle - } - - fn cycle(&self) -> usize { - self.cycle - } - - fn set_theta(&mut self, theta: Theta) { - self.theta = theta; - } - - fn theta(&self) -> &Theta { - &self.theta - } - - fn psi(&self) -> &Psi { - &self.psi - } - - fn likelihood(&self) -> f64 { - self.objf - } - - 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; - } - - 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 - ); - } - - if (self.last_objf - self.objf).abs() <= THETA_F { - 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()); - } - - 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()); - } - - if std::path::Path::new("stop").exists() { - tracing::warn!("Stopfile detected - breaking"); - self.converged = true; - self.set_status(Status::Stop(StopReason::StopFile)); - self.log_cycle_state(); - return Ok(self.status.clone()); } - - self.status = Status::Continue; - self.log_cycle_state(); - Ok(self.status.clone()) } +} - fn estimation(&mut self) -> Result<()> { - let error_model: AssayErrorModels = self.error_models.clone(); +impl Runner for Npod { + fn step(&mut self, state: &mut FitState) -> Result { + let theta = if state.is_first_cycle() { + state.solution().theta().clone() + } else { + expand(state, self.config.min_distance)? + }; - self.psi = calculate_psi( - &self.equation, - &self.data, - &self.theta, - &error_model, - self.cycle == 1 && self.config.progress, + let estimate = state.estimate(theta, self.config.progress && state.is_first_cycle())?; + let solution = state.condense( + estimate, + self.config.prune_threshold, + self.config.qr_tolerance, )?; - - if let Err(err) = self.check_zero_probability_subjects() { - bail!(err); - } - - (self.lambda, _) = match burke(&self.psi) { - Ok((lambda, objf)) => (lambda, objf), - Err(err) => { - bail!(err); - } - }; - Ok(()) + state.set_solution(solution); + self.error_optimizer.optimize(state)?; + + let converged = state + .objective_change() + .is_some_and(|change| change <= self.config.objective_tolerance); + Ok(if converged { + CycleOutcome::Converged + } else { + CycleOutcome::Continue + }) } +} - 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()); - - (self.lambda, self.objf) = match burke(&self.psi) { - Ok((lambda, objf)) => (lambda, objf), - Err(err) => { - return Err(anyhow::anyhow!("Error in IPM: {:?}", err)); - } - }; - self.w = self.lambda.clone(); - Ok(()) - } +/// Move each support point towards a local maximum of the directional derivative of the +/// likelihood, and return the support points with the results added as new candidates. +fn expand(state: &FitState, min_distance: f64) -> Result { + let solution = state.solution(); + // The population likelihood of each subject, Ψ·w. `ParameterOptimizer` works with raw + // likelihoods, so undo the per-subject scaling of Ψ. + // TODO: Move NPOD and its parameter optimizer fully into log space. + let pyl: Array1 = (solution.psi().matrix() * solution.weights().weights()) + .iter() + .zip(solution.psi().row_log_scales()) + .map(|(pyl, row_log_scale)| pyl * row_log_scale.exp()) + .collect(); + + let candidates: Vec> = solution + .theta() + .matrix() + .row_iter() + .map(|row| row.iter().copied().collect()) + .collect::>() + .into_par_iter() + .map(|point| { + ParameterOptimizer::new(state.equation(), state.data(), state.error_models(), &pyl) + .optimize_point(point) + }) + .collect::>()?; - fn optimizations(&mut self) -> Result<()> { - optimize_error_models( - &self.equation, - &self.data, - &self.theta, - &mut self.error_models, - &mut self.gamma_delta, - &mut self.objf, - &mut self.lambda, - &mut self.psi, - &self.config.error_optim, - ) + let mut theta = solution.theta().clone(); + for candidate in candidates { + theta.suggest_point(&candidate.to_vec(), min_distance)?; } + Ok(theta) +} - fn expansion(&mut self) -> Result<()> { - let pyl_col = self.psi().matrix().as_ref() * self.w.weights().as_ref(); - let mut pyl: Array1 = pyl_col.iter().copied().collect(); - - // ParameterOptimizer currently evaluates raw likelihoods, so restore - // pyl to the same scale before calculating directional derivatives. - // TODO: Move NPOD and its parameter optimizer fully into log space. - for (pyl_i, row_log_scale) in pyl.iter_mut().zip(self.psi().row_log_scales()) { - *pyl_i *= row_log_scale.exp(); - } - - let error_model: AssayErrorModels = self.error_models.clone(); - - let mut candididate_points: Vec> = Vec::default(); - for spp in self.theta.matrix().row_iter() { - let candidate: Vec = spp.iter().cloned().collect(); - let spp = Array1::from(candidate); - candididate_points.push(spp.to_owned()); - } - candididate_points.par_iter_mut().for_each(|spp| { - let optimizer = ParameterOptimizer::new(&self.equation, &self.data, &error_model, &pyl); - let candidate_point = optimizer.optimize_point(spp.to_owned()).unwrap(); - *spp = candidate_point; - }); - for cp in candididate_points { - self.theta.suggest_point(cp.to_vec().as_slice(), THETA_D)?; - } +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn config_without_the_newer_fields_still_deserializes() -> anyhow::Result<()> { + let config: NpodConfig = serde_json::from_str(r#"{"max_cycles": 7, "progress": false}"#)?; + + let defaults = NpodConfig::default(); + assert_eq!(config.max_cycles, 7); + assert!(!config.progress); + assert_eq!(config.objective_tolerance, defaults.objective_tolerance); + assert_eq!(config.min_distance, defaults.min_distance); + assert_eq!(config.prune_threshold, defaults.prune_threshold); + assert_eq!(config.qr_tolerance, defaults.qr_tolerance); Ok(()) } } diff --git a/src/algorithms/nonparametric/solution.rs b/src/algorithms/nonparametric/solution.rs new file mode 100644 index 000000000..81d81861c --- /dev/null +++ b/src/algorithms/nonparametric/solution.rs @@ -0,0 +1,133 @@ +//! A set of support points together with everything derived from them. + +use anyhow::Result; + +use crate::estimation::nonparametric::{ipm::burke, Psi, Theta, Weights}; + +/// Support points with their likelihoods, weights, and the log-likelihood those give. +/// +/// The four parts always describe the same fit. A `Solution` is built either by solving for the +/// optimal weights ([`optimal`](Self::optimal)) or by computing the log-likelihood of given +/// weights ([`with_weights`](Self::with_weights)), and the fit only ever replaces it as a whole. +/// The one deliberate exception is [`with_stale_weights`](Self::with_stale_weights). +#[derive(Debug, Clone)] +pub(crate) struct Solution { + theta: Theta, + /// Likelihood of each subject (row) at each support point (column) of `theta`. + psi: Psi, + /// Probability of each support point in `theta`. + weights: Weights, + /// Log-likelihood of the data given `theta` and `weights`. + log_likelihood: f64, +} + +impl Solution { + /// The starting point of a fit: the prior support points, not yet evaluated. + pub(crate) fn unsolved(prior: Theta) -> Self { + Self { + theta: prior, + psi: Psi::new(), + weights: Weights::default(), + log_likelihood: f64::NEG_INFINITY, + } + } + + /// Solve for the weights that maximize the likelihood, using Burke's interior point method. + pub(crate) fn optimal(theta: Theta, psi: Psi) -> Result { + let (weights, log_likelihood) = burke(&psi)?; + Ok(Self { + theta, + psi, + weights, + log_likelihood, + }) + } + + /// Use the given weights and compute the log-likelihood they give. + pub(crate) fn with_weights(theta: Theta, psi: Psi, weights: Weights) -> Self { + let log_likelihood = marginal_log_likelihood(&psi, &weights); + Self { + theta, + psi, + weights, + log_likelihood, + } + } + + /// Replace the weights, keeping this solution's Ψ and log-likelihood. + /// + /// This breaks the invariant on purpose: the error-model optimization has always adopted a + /// better Ψ while keeping the weights from before it. Whether the weights should follow is + /// an open question (LAPKB/PMcore#334); dropping this call is all it takes to change that. + pub(crate) fn with_stale_weights(self, weights: Weights) -> Self { + Self { weights, ..self } + } + + pub(crate) fn theta(&self) -> &Theta { + &self.theta + } + + pub(crate) fn psi(&self) -> &Psi { + &self.psi + } + + pub(crate) fn weights(&self) -> &Weights { + &self.weights + } + + pub(crate) fn log_likelihood(&self) -> f64 { + self.log_likelihood + } + + pub(crate) fn into_parts(self) -> (Theta, Psi, Weights, f64) { + (self.theta, self.psi, self.weights, self.log_likelihood) + } +} + +/// Log-likelihood of the data under a discrete `(psi, weights)` model. +fn marginal_log_likelihood(psi: &Psi, weights: &Weights) -> f64 { + let m = psi.matrix(); + (0..m.nrows()) + .map(|s| { + let acc: f64 = (0..m.ncols()).map(|j| *m.get(s, j) * weights[j]).sum(); + acc.max(f64::MIN_POSITIVE).ln() + }) + .sum::() + + psi.log_scale() +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_relative_eq; + use ndarray::Array2; + + #[test] + fn marginal_log_likelihood_restores_subject_scaling() -> anyhow::Result<()> { + let log_likelihoods = Array2::from_shape_vec((2, 2), vec![-1000.0, -1001.0, -2.0, -4.0])?; + let psi = Psi::from_log_likelihoods(log_likelihoods)?; + let weights = Weights::from_vec(vec![0.6, 0.4]); + + let expected = -1000.0 + (0.6 + 0.4 * (-1.0_f64).exp()).ln() - 2.0 + + (0.6 + 0.4 * (-2.0_f64).exp()).ln(); + + assert_relative_eq!( + marginal_log_likelihood(&psi, &weights), + expected, + epsilon = 1e-12 + ); + Ok(()) + } + + #[test] + fn optimal_weights_give_the_reported_log_likelihood() -> anyhow::Result<()> { + let log_likelihoods = Array2::from_shape_vec((2, 2), vec![-1.0, -3.0, -4.0, -2.0])?; + let psi = Psi::from_log_likelihoods(log_likelihoods)?; + + let solution = Solution::optimal(Theta::new(), psi)?; + let recomputed = marginal_log_likelihood(solution.psi(), solution.weights()); + + assert_relative_eq!(solution.log_likelihood(), recomputed, epsilon = 1e-8); + Ok(()) + } +} diff --git a/src/algorithms/parametric/mod.rs b/src/algorithms/parametric/mod.rs index 9951345f1..cb1bd756d 100644 --- a/src/algorithms/parametric/mod.rs +++ b/src/algorithms/parametric/mod.rs @@ -49,8 +49,8 @@ use pharmsol::prelude::simulator::Equation; /// // `problem.fit_with(config)` accepts the config directly. /// ``` /// -/// [`SaemConfig`] implements [`Algorithm`] by delegating to the matching enum variant, so -/// configs can be passed to `fit_with` without converting them first. +/// Anything that converts into a `ParametricAlgorithm`, including [`SaemConfig`], implements +/// [`Algorithm`], so configs can be passed to `fit_with` without converting them first. #[derive(Debug, Clone)] pub enum ParametricAlgorithm { /// Stochastic Approximation Expectation-Maximization. @@ -77,24 +77,19 @@ impl ParametricAlgorithm { } } -impl Algorithm for ParametricAlgorithm { +/// Every parametric algorithm (or its configuration) can be fitted directly. +impl Algorithm for A +where + E: Equation, + 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/bestdose/mod.rs b/src/bestdose/mod.rs index b32d4104e..cc668b26a 100644 --- a/src/bestdose/mod.rs +++ b/src/bestdose/mod.rs @@ -28,11 +28,11 @@ //! let (theta, weights) = match past_data { //! Some(past) => { //! let post = EstimationProblem::nonparametric( -//! eq.clone(), data::Data::new(vec![past]), fit.get_theta().clone(), ems.clone())? +//! eq.clone(), data::Data::new(vec![past]), fit.theta().clone(), ems.clone())? //! .fit_with(NcnpagConfig::default())?; // or NpmapConfig::default() -//! (post.get_theta().clone(), post.weights().clone()) +//! (post.theta().clone(), post.weights().clone()) //! } -//! None => (fit.get_theta().clone(), fit.weights().clone()), +//! None => (fit.theta().clone(), fit.weights().clone()), //! }; //! //! // 3. Optimize doses. diff --git a/src/estimation/nonparametric/cycles.rs b/src/estimation/nonparametric/cycles.rs index a5333033c..bc2b51e9e 100644 --- a/src/estimation/nonparametric/cycles.rs +++ b/src/estimation/nonparametric/cycles.rs @@ -48,6 +48,7 @@ impl NPCycle { pub fn cycle(&self) -> usize { self.cycle } + /// The objective function at the end of this cycle, -2 × log-likelihood. pub fn objf(&self) -> f64 { self.objf } diff --git a/src/estimation/nonparametric/mod.rs b/src/estimation/nonparametric/mod.rs index 57c980202..ad599603b 100644 --- a/src/estimation/nonparametric/mod.rs +++ b/src/estimation/nonparametric/mod.rs @@ -18,9 +18,9 @@ pub use cycles::{CycleLog, NPCycle}; pub(crate) use expansion::adaptative_grid; pub use ipm::burke; pub use posterior::{posterior, Posterior}; -pub use predictions::{NPPredictionRow, NPPredictions}; +pub use predictions::{NPPredictionRow, NPPredictions, PredictionGrid}; pub(crate) use psi::calculate_psi; -pub use psi::Psi; +pub use psi::{Psi, ZeroProbabilityError}; pub use result::NonParametricResult; pub use statistics::{median, population_mean_median, posterior_mean_median, weighted_median}; pub use summaries::{fit_summary, individual_summaries, population_summary}; diff --git a/src/estimation/nonparametric/predictions.rs b/src/estimation/nonparametric/predictions.rs index 6cd431dfa..05189366d 100644 --- a/src/estimation/nonparametric/predictions.rs +++ b/src/estimation/nonparametric/predictions.rs @@ -60,6 +60,43 @@ impl NPPredictionRow { } } +/// The times at which to predict, besides the observation times. +/// +/// The default predicts at the observation times only. With an `idelta`, predictions are also +/// made every `idelta` time units, up to `tad` time units after the last dose of each occasion. +/// +/// ``` +/// use pmcore::prelude::*; +/// +/// let grid = PredictionGrid::new().idelta(0.5).tad(12.0); +/// ``` +#[derive(Debug, Clone, Copy, PartialEq, Default, Serialize, Deserialize)] +pub struct PredictionGrid { + /// Interval between the added prediction times; `0` adds none. + pub idelta: f64, + /// How far past the last dose of each occasion to keep predicting. + pub tad: f64, +} + +impl PredictionGrid { + /// Predictions at the observation times only. + pub fn new() -> Self { + Self::default() + } + + /// Also predict every `idelta` time units. + pub fn idelta(mut self, idelta: f64) -> Self { + self.idelta = idelta; + self + } + + /// Keep predicting `tad` time units past the last dose of each occasion. + pub fn tad(mut self, tad: f64) -> Self { + self.tad = tad; + self + } +} + #[derive(Debug, Clone, Serialize, Deserialize)] pub struct NPPredictions { predictions: Vec, @@ -118,15 +155,14 @@ impl NPPredictions { } pub fn calculate( - equation: &(impl pharmsol::prelude::simulator::Equation + Sync), + equation: &impl pharmsol::prelude::simulator::Equation, data: &Data, theta: &Theta, w: &Weights, posterior: &Posterior, - idelta: f64, - tad: f64, + grid: PredictionGrid, ) -> Result { - let data = data.clone().expand(idelta, tad); + let data = data.expand(grid.idelta, grid.tad); let subjects = data.subjects(); if subjects.len() != posterior.matrix().nrows() { diff --git a/src/estimation/nonparametric/psi.rs b/src/estimation/nonparametric/psi.rs index c6a744854..388d5934a 100644 --- a/src/estimation/nonparametric/psi.rs +++ b/src/estimation/nonparametric/psi.rs @@ -255,6 +255,42 @@ impl<'de> Deserialize<'de> for Psi { } } +/// Some subjects have zero probability given the model: their likelihood is zero (or not a +/// number) at every support point, so no weighting of the support points can explain their data. +/// +/// Fits return this inside an [`anyhow::Error`]; recover it with +/// [`downcast_ref`](anyhow::Error::downcast_ref) to see which subjects are affected. +#[derive(Debug, Clone, PartialEq, Eq)] +pub struct ZeroProbabilityError { + subjects: Vec, + total: usize, +} + +impl ZeroProbabilityError { + pub(crate) fn new(subjects: Vec, total: usize) -> Self { + Self { subjects, total } + } + + /// The IDs of the affected subjects. + pub fn subjects(&self) -> &[String] { + &self.subjects + } +} + +impl std::fmt::Display for ZeroProbabilityError { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + write!( + f, + "The probability of {}/{} subjects is zero given the model. Affected subjects: {:?}", + self.subjects.len(), + self.total, + self.subjects + ) + } +} + +impl std::error::Error for ZeroProbabilityError {} + pub(crate) fn calculate_psi( equation: &impl Equation, subjects: &Data, @@ -267,14 +303,65 @@ pub(crate) fn calculate_psi( let log_psi = log_likelihood_matrix(equation, subjects, &theta_ndarray, error_models, progress)?; + let zero = zero_probability_rows(&log_psi); + if !zero.is_empty() { + let all = subjects.subjects(); + let ids = zero.iter().map(|&i| all[i].id().clone()).collect(); + return Err(ZeroProbabilityError::new(ids, log_psi.nrows()).into()); + } + Psi::from_log_likelihoods(log_psi) } +/// Rows (subjects) without a single log-likelihood above `-inf`, ignoring NaN. +fn zero_probability_rows(log_likelihoods: &Array2) -> Vec { + log_likelihoods + .axis_iter(Axis(0)) + .enumerate() + .filter(|(_, row)| { + row.iter().copied().fold(f64::NEG_INFINITY, f64::max) == f64::NEG_INFINITY + }) + .map(|(i, _)| i) + .collect() +} + #[cfg(test)] mod tests { use super::*; use ndarray::Array2; + #[test] + fn subjects_without_a_finite_log_likelihood_have_zero_probability() { + let log_likelihoods = Array2::from_shape_vec( + (3, 2), + vec![ + -1.0, + -2.0, + f64::NEG_INFINITY, + f64::NAN, + -3.0, + f64::NEG_INFINITY, + ], + ) + .unwrap(); + + assert_eq!(zero_probability_rows(&log_likelihoods), vec![1]); + } + + #[test] + fn zero_probability_error_survives_anyhow() { + let err: anyhow::Error = ZeroProbabilityError::new(vec!["7".into()], 10).into(); + + let zero = err + .downcast_ref::() + .expect("should downcast"); + assert_eq!(zero.subjects(), ["7".to_string()]); + assert_eq!( + err.to_string(), + "The probability of 1/10 subjects is zero given the model. Affected subjects: [\"7\"]" + ); + } + #[test] fn log_likelihood_rows_are_scaled_before_exponentiation() -> Result<()> { let log_likelihoods = Array2::from_shape_vec((2, 2), vec![-1000.0, -1001.0, -2.0, -4.0])?; diff --git a/src/estimation/nonparametric/result.rs b/src/estimation/nonparametric/result.rs index 7af6ca340..1985de3c6 100644 --- a/src/estimation/nonparametric/result.rs +++ b/src/estimation/nonparametric/result.rs @@ -3,8 +3,10 @@ use std::path::Path; use pharmsol::Equation; use serde::Serialize; -use crate::algorithms::Status; -use crate::estimation::nonparametric::{CycleLog, NPPredictions, Posterior, Psi, Theta, Weights}; +use crate::algorithms::{Status, StopReason}; +use crate::estimation::nonparametric::{ + CycleLog, NPPredictions, Posterior, PredictionGrid, Psi, Theta, Weights, +}; use pharmsol::{AssayErrorModels, Data}; @@ -58,6 +60,7 @@ impl NonParametricResult { self.cycles } + /// The objective function, -2 × log-likelihood (`neg2ll` in `cycles.csv`). pub fn objf(&self) -> f64 { self.objf } @@ -66,14 +69,25 @@ impl NonParametricResult { self.status.converged() } - pub fn get_theta(&self) -> &Theta { + /// How the fit ended. + pub fn status(&self) -> &Status { + &self.status + } + + /// Why the fit stopped. + pub fn stop_reason(&self) -> Option<&StopReason> { + self.status.stop_reason() + } + + /// The optimized support points. + pub fn theta(&self) -> &Theta { &self.theta } /// The prior distribution ([`Theta`]) that seeded the algorithm. /// /// This is the initial set of support points, as opposed to the optimized - /// solution returned by [`get_theta`](Self::get_theta). + /// solution returned by [`theta`](Self::theta). pub fn prior(&self) -> &Theta { &self.prior } @@ -142,9 +156,9 @@ impl NonParametricResult { /// Compute predictions on demand. Nothing is cached on the result; callers /// that need the predictions repeatedly should hold on to the returned /// value themselves. - pub fn predictions(&self, idelta: f64, tad: f64) -> anyhow::Result { + pub fn predictions(&self, grid: PredictionGrid) -> anyhow::Result { let posterior = self.posterior()?; - self.predictions_with(&posterior, idelta, tad) + self.predictions_with(&posterior, grid) } /// Like [`predictions`](Self::predictions), but reuses an already-computed @@ -152,8 +166,7 @@ impl NonParametricResult { fn predictions_with( &self, posterior: &Posterior, - idelta: f64, - tad: f64, + grid: PredictionGrid, ) -> anyhow::Result { NPPredictions::calculate( &self.equation, @@ -161,8 +174,7 @@ impl NonParametricResult { &self.theta, &self.weights, posterior, - idelta, - tad, + grid, ) } @@ -316,10 +328,9 @@ impl NonParametricResult { /// Compute and write the population and posterior predictions to a CSV file /// readable by Pmetrics. /// - /// `idelta` is the interval used to densify the prediction grid and `tad` is - /// the additional time after the last event to simulate. - pub fn write_predictions(&self, path: &Path, idelta: f64, tad: f64) -> anyhow::Result<()> { - let predictions = self.predictions(idelta, tad)?; + /// `grid` sets the prediction times; see [`PredictionGrid`]. + pub fn write_predictions(&self, path: &Path, grid: PredictionGrid) -> anyhow::Result<()> { + let predictions = self.predictions(grid)?; predictions.write(path) } @@ -327,11 +338,11 @@ impl NonParametricResult { /// /// This includes the data, error models, prior, optimized support points, /// likelihoods, weights, objective function, status, cycle log, posterior - /// probabilities, and predictions. `idelta` and `tad` control the density of - /// the embedded predictions (see [`predictions`](Self::predictions)). - pub fn write_json(&self, path: &Path, idelta: f64, tad: f64) -> anyhow::Result<()> { + /// probabilities, and predictions. `grid` sets the prediction times; see + /// [`PredictionGrid`]. + pub fn write_json(&self, path: &Path, grid: PredictionGrid) -> anyhow::Result<()> { let posterior = self.posterior()?; - let predictions = self.predictions_with(&posterior, idelta, tad)?; + let predictions = self.predictions_with(&posterior, grid)?; self.write_json_with(path, &posterior, &predictions) } @@ -376,20 +387,18 @@ impl NonParametricResult { /// - `cycles.csv` — cycle log /// /// Use [`write_json`](Self::write_json) to explicitly write the complete result. - /// `idelta` and `tad` control the density of the predictions (see - /// [`predictions`](Self::predictions)). + /// `grid` sets the prediction times; see [`PredictionGrid`]. pub fn write_outputs( &self, directory: impl AsRef, - idelta: f64, - tad: f64, + grid: PredictionGrid, ) -> anyhow::Result<()> { let dir = directory.as_ref(); std::fs::create_dir_all(dir)?; // Compute the expensive artifacts a single time and reuse them. let posterior = self.posterior()?; - let predictions = self.predictions_with(&posterior, idelta, tad)?; + let predictions = self.predictions_with(&posterior, grid)?; self.write_theta(&dir.join("theta.csv"))?; self.write_posterior_with(&dir.join("posterior.csv"), &posterior)?; @@ -487,13 +496,13 @@ mod tests { .fit_with(NpagConfig::new().max_cycles(2)) .unwrap(); - let n_spp = r1.get_theta().nspp(); + let n_spp = r1.theta().nspp(); assert!(n_spp > 0, "NPAG should produce support points"); // Chain into NPOD — should complete without error let r2 = r1.chain(NpodConfig::new().max_cycles(1)).unwrap(); - assert!(r2.get_theta().nspp() > 0); + assert!(r2.theta().nspp() > 0); assert_eq!(r2.data().subjects().len(), 1); assert_eq!(r2.cycles(), 1); } diff --git a/src/estimation/nonparametric/summaries.rs b/src/estimation/nonparametric/summaries.rs index 6d979aaef..bded343f9 100644 --- a/src/estimation/nonparametric/summaries.rs +++ b/src/estimation/nonparametric/summaries.rs @@ -12,18 +12,18 @@ pub fn fit_summary(result: &NonParametricResult) -> FitSummary { iterations: result.cycles(), subject_count: result.data().subjects().len(), observation_count: count_observations(result.data()), - parameter_count: result.get_theta().parameters().len(), + parameter_count: result.theta().parameters().len(), } } pub fn population_summary(result: &NonParametricResult) -> PopulationSummary { - let theta_matrix = to_ndarray_matrix(result.get_theta().matrix()); + let theta_matrix = to_ndarray_matrix(result.theta().matrix()); let weights = Array1::from_iter(result.weights().iter()); let (mean, median) = population_mean_median(&theta_matrix, &weights) .expect("population summary should be derivable from theta and weights"); let parameters = result - .get_theta() + .theta() .parameters() .names() .into_iter() @@ -54,12 +54,12 @@ pub fn population_summary(result: &NonParametricResult) -> Popul pub fn individual_summaries( result: &NonParametricResult, ) -> Vec { - let theta_matrix = to_ndarray_matrix(result.get_theta().matrix()); + let theta_matrix = to_ndarray_matrix(result.theta().matrix()); let psi_matrix = to_ndarray_matrix(result.psi().matrix()); let weights = Array1::from_iter(result.weights().iter()); let (means, _) = posterior_mean_median(&theta_matrix, &psi_matrix, &weights) .expect("individual summaries should be derivable from theta, psi, and weights"); - let parameter_names = result.get_theta().parameters().names(); + let parameter_names = result.theta().parameters().names(); result .data() diff --git a/src/iov/mod.rs b/src/iov/mod.rs index 2ba67803c..83466d951 100644 --- a/src/iov/mod.rs +++ b/src/iov/mod.rs @@ -20,8 +20,8 @@ //! use pmcore::iov::DiffusionOptimize; //! use pmcore::iov::DiffusionConfig; //! -//! let r_ode = problem.fit_with(NPAG::default())?; -//! let mut joint = r_ode.get_theta() +//! let r_ode = problem.fit_with(NpagConfig::new())?; +//! let mut joint = r_ode.theta() //! .with_added_parameter("ske", 1e-6, 1.0, 0.01)?; //! //! let diff = sde.optimize_diffusion( diff --git a/src/lib.rs b/src/lib.rs index 58a6b8da7..b5a657401 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -72,9 +72,10 @@ pub mod prelude { BoundedParameter, Parameter, ParameterScale, ParameterSpace, UnboundedParameter, }; - pub use crate::algorithms::nonparametric::{CycleFlow, FitController, FitObserver}; + pub use crate::algorithms::nonparametric::FitController; pub use crate::estimation::nonparametric::{ - CycleLog, NPCycle, NPPredictions, NonParametricResult, Posterior, Psi, Theta, Weights, + CycleLog, NPCycle, NPPredictions, NonParametricResult, Posterior, PredictionGrid, Psi, + Theta, Weights, ZeroProbabilityError, }; pub use crate::iov::{DiffusionConfig, DiffusionOptimize, DiffusionResult}; pub use crate::model::{EquationMetadataSource, ModelMetadata}; diff --git a/tests/bestdose_tests.rs b/tests/bestdose_tests.rs index 35579df62..c04cff81c 100644 --- a/tests/bestdose_tests.rs +++ b/tests/bestdose_tests.rs @@ -357,7 +357,7 @@ fn ncnpag_individualizes_toward_matching_point() -> Result<()> { assert!((weight_sum - 1.0).abs() < 1e-6); // The dominant support point matches the data-generating parameters. - let point = best_point(posterior.get_theta(), posterior.weights()); + let point = best_point(posterior.theta(), posterior.weights()); assert!( (point[0] - matching[0]).abs() < 0.05, "expected ke ≈ {}, got {}", @@ -381,7 +381,7 @@ fn npmap_reweights_toward_matching_point() -> Result<()> { )? .fit_with(NpmapConfig::default())?; - let point = best_point(posterior.get_theta(), posterior.weights()); + let point = best_point(posterior.theta(), posterior.weights()); assert!( (point[0] - matching[0]).abs() < 0.05, "expected ke ≈ {}, got {}", @@ -409,7 +409,7 @@ fn ncnpag_posterior_feeds_dose_optimization() -> Result<()> { let problem = BestDoseProblem::new( bolus_model(), - posterior.get_theta().clone(), + posterior.theta().clone(), posterior.weights().clone(), )?; diff --git a/tests/onecomp.rs b/tests/onecomp.rs index c63509b52..65180c3cf 100644 --- a/tests/onecomp.rs +++ b/tests/onecomp.rs @@ -62,7 +62,7 @@ fn test_one_compartment_npag() -> Result<()> { AssayErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0), )?; let result = EstimationProblem::nonparametric(eq, data, prior, error_models)? - .fit_with(NonParametricAlgorithm::npag())?; + .fit_with(NpagConfig::new())?; // Check the results assert_eq!(result.cycles(), 31); @@ -71,7 +71,7 @@ fn test_one_compartment_npag() -> Result<()> { // The prior is preserved on the result and is distinct from the optimized // solution (which is condensed to far fewer support points). assert_eq!(result.prior().nspp(), 100); - assert!(result.get_theta().nspp() < result.prior().nspp()); + assert!(result.theta().nspp() < result.prior().nspp()); Ok(()) } @@ -128,7 +128,7 @@ fn test_one_compartment_npod() -> Result<()> { AssayErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0), )?; let result = EstimationProblem::nonparametric(eq, data, prior, error_models)? - .fit_with(NonParametricAlgorithm::npod())?; + .fit_with(NpodConfig::new())?; // Check the results assert_eq!(result.cycles(), 11); @@ -192,13 +192,14 @@ fn test_one_compartment_postprob() -> Result<()> { AssayErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0), )?; let result = EstimationProblem::nonparametric(eq, data, theta.clone(), error_models)? - .fit_with(NonParametricAlgorithm::npmap())?; + .fit_with(NpmapConfig::new())?; - // Check the results - assert_eq!(result.cycles(), 0); + // NPMAP is a single reweighting pass + 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()); + assert_eq!(result.theta().nspp(), theta.nspp()); Ok(()) } diff --git a/tests/results_summary_tests.rs b/tests/results_summary_tests.rs index 88363e4b9..1e7632928 100644 --- a/tests/results_summary_tests.rs +++ b/tests/results_summary_tests.rs @@ -49,7 +49,7 @@ fn test_nonparametric_fit_result_summary_surface() -> Result<()> { let error_models = AssayErrorModels::new().add("0", assay_error)?; let result = EstimationProblem::nonparametric(simple_equation(), simple_data(), prior, error_models)? - .fit_with(NonParametricAlgorithm::npag())?; + .fit_with(NpagConfig::new())?; let summary = result.summary();