diff --git a/examples/bestdose.rs b/examples/bestdose.rs index 34ff95364..6968514a6 100644 --- a/examples/bestdose.rs +++ b/examples/bestdose.rs @@ -1,6 +1,5 @@ use anyhow::Result; -use pmcore::bestdose; // bestdose new - // use pmcore::bestdose::bestdose_old as bestdose; // bestdose old +use pmcore::bestdose::{BestDosePosterior, DoseRange, Target}; use pmcore::prelude::*; use pmcore::routines::initialization::parse_prior; @@ -78,22 +77,15 @@ fn main() -> Result<()> { ) .unwrap(); - // Example usage - using new() constructor which calculates NPAGFULL11 posterior - // max_cycles controls NPAGFULL refinement: - // 0 = NPAGFULL11 only (fast but less accurate) - // 100 = moderate refinement - // 500 = full refinement (Fortran default, slow but most accurate) - let problem = bestdose::BestDoseProblem::new( + // Example usage - two-stage API: + // Stage 1: Compute posterior (expensive, done once) + // Stage 2: Optimize doses (can be called multiple times with different params) + let posterior = BestDosePosterior::compute( &theta, &prior.unwrap(), Some(past_data.clone()), // Optional: past data for Bayesian updating - target_data.clone(), - None, eq.clone(), - bestdose::DoseRange::new(0.0, 300.0), - 0.0, settings.clone(), - bestdose::Target::Concentration, // Target concentrations (not AUCs) )?; println!("Optimizing dose..."); @@ -103,7 +95,13 @@ fn main() -> Result<()> { for bias_weight in &bias_weights { println!("Running optimization with bias weight: {}", bias_weight); - let optimal = problem.clone().with_bias_weight(*bias_weight).optimize()?; + let optimal = posterior.optimize( + target_data.clone(), + None, + DoseRange::new(0.0, 300.0), + *bias_weight, + Target::Concentration, + )?; results.push((bias_weight, optimal)); } @@ -124,7 +122,7 @@ fn main() -> Result<()> { // Print concentration-time predictions for the optimal dose let optimal = &results.last().unwrap().1; println!("\nConcentration-time predictions for optimal dose:"); - for pred in optimal.predictions().predictions().into_iter() { + for pred in optimal.predictions().predictions().iter() { println!( "Time: {:.2} h, Observed: {:.2}, (Pop Mean: {:.4}, Pop Median: {:.4}, Post Mean: {:.4}, Post Median: {:.4})", pred.time(), pred.obs().unwrap_or(0.0), pred.pop_mean(), pred.pop_median(), pred.post_mean(), pred.post_median() diff --git a/examples/bestdose_auc.rs b/examples/bestdose_auc.rs index 6f9b08c27..70aecb429 100644 --- a/examples/bestdose_auc.rs +++ b/examples/bestdose_auc.rs @@ -1,5 +1,5 @@ use anyhow::Result; -use pmcore::bestdose::{BestDoseProblem, DoseRange, Target}; +use pmcore::bestdose::{BestDosePosterior, DoseRange, Target}; use pmcore::prelude::*; use pmcore::routines::initialization::parse_prior; @@ -61,22 +61,23 @@ fn main() -> Result<()> { .observation(12.0, 80.0, 0) // Target AUC at 12h .build(); - println!("Creating BestDose problem with AUC targets..."); - let problem = BestDoseProblem::new( + println!("Creating BestDose posterior (no past data - use prior directly)..."); + let posterior = BestDosePosterior::compute( &theta, weights, None, // No past data - use prior directly - target_data.clone(), - None, eq.clone(), - DoseRange::new(100.0, 2000.0), // Wider range for AUC targets - 0.8, // for AUC targets higher bias_weight usually works best settings.clone(), - Target::AUCFromZero, // Cumulative AUC from time 0 )?; println!("Optimizing dose...\n"); - let optimal = problem.optimize()?; + let optimal = posterior.optimize( + target_data.clone(), + None, + DoseRange::new(100.0, 2000.0), // Wider range for AUC targets + 0.8, // for AUC targets higher bias_weight usually works best + Target::AUCFromZero, // Cumulative AUC from time 0 + )?; let opt_doses = optimal.doses(); @@ -137,22 +138,14 @@ fn main() -> Result<()> { .build(); println!("Creating BestDose problem with interval AUC target..."); - let problem_interval = BestDoseProblem::new( - &theta, - weights, - None, + let optimal_interval = posterior.optimize( target_interval.clone(), None, - eq.clone(), DoseRange::new(50.0, 500.0), 0.8, - settings.clone(), Target::AUCFromLastDose, // Interval AUC from last dose! )?; - println!("Optimizing maintenance dose...\n"); - let optimal_interval = problem_interval.optimize()?; - let doses: Vec = optimal_interval.doses(); println!("=== INTERVAL AUC RESULTS ==="); diff --git a/examples/bestdose_bounds.rs b/examples/bestdose_bounds.rs index d91709b39..4b2eab3bb 100644 --- a/examples/bestdose_bounds.rs +++ b/examples/bestdose_bounds.rs @@ -1,5 +1,5 @@ use anyhow::Result; -use pmcore::bestdose::{BestDoseProblem, DoseRange, Target}; +use pmcore::bestdose::{BestDosePosterior, DoseRange, Target}; use pmcore::prelude::*; use pmcore::routines::initialization::parse_prior; @@ -67,40 +67,20 @@ fn main() -> Result<()> { println!("{:<30} | {:>12} | {:>10}", "Range", "Optimal Dose", "Cost"); println!("{}", "-".repeat(60)); + // Compute posterior once, reuse for all dose ranges + let posterior = + BestDosePosterior::compute(&theta, weights, None, eq.clone(), settings.clone())?; + for (min, max, description) in dose_ranges { - let problem = BestDoseProblem::new( - &theta, - weights, - None, + let result = posterior.optimize( target_data.clone(), None, - eq.clone(), DoseRange::new(min, max), 0.5, - settings.clone(), Target::Concentration, )?; - let result = problem.optimize()?; - - let doses: Vec = result - .optimal_subject() - .iter() - .map(|occ| { - occ.iter() - .filter(|event| match event { - Event::Bolus(_) => true, - Event::Infusion(_) => true, - _ => false, - }) - .map(|event| match event { - Event::Bolus(bolus) => bolus.amount(), - Event::Infusion(infusion) => infusion.amount(), - _ => 0.0, - }) - }) - .flatten() - .collect(); + let doses: Vec = result.doses(); // Check if dose hit the bound let at_bound = if (doses[0] - max).abs() < 1.0 { diff --git a/src/bestdose/cost.rs b/src/bestdose/cost.rs index c3ba15aa0..b9dede94c 100644 --- a/src/bestdose/cost.rs +++ b/src/bestdose/cost.rs @@ -129,7 +129,7 @@ use pharmsol::Equation; /// - Model simulation fails /// - Prediction length doesn't match observation count /// - AUC calculation fails (for AUC targets) -pub fn calculate_cost(problem: &BestDoseProblem, candidate_doses: &[f64]) -> Result { +pub(crate) fn calculate_cost(problem: &BestDoseProblem, candidate_doses: &[f64]) -> Result { // Validate candidate_doses length matches expected optimizable dose count let expected_optimizable = problem .target @@ -233,7 +233,7 @@ pub fn calculate_cost(problem: &BestDoseProblem, candidate_doses: &[f64]) -> Res // Calculate variance (using posterior weights) and population mean (using prior weights) - for ((row, post_prob), prior_prob) in problem + for ((row, post_prob), _prior_prob) in problem .theta .matrix() .row_iter() @@ -509,8 +509,8 @@ pub fn calculate_cost(problem: &BestDoseProblem, candidate_doses: &[f64]) -> Res let pj = preds_i[j]; let se = (obs_val - pj).powi(2); sumsq_i += se; - // Calculate population mean using PRIOR probabilities - y_bar[j] += prior_prob * pj; + // Calculate population mean using POSTERIOR probabilities + y_bar[j] += post_prob * pj; } variance += post_prob * sumsq_i; // Weighted by posterior diff --git a/src/bestdose/mod.rs b/src/bestdose/mod.rs index f9588dc33..4502b817e 100644 --- a/src/bestdose/mod.rs +++ b/src/bestdose/mod.rs @@ -9,38 +9,37 @@ //! # Quick Start //! //! ```rust,no_run,ignore -//! use pmcore::bestdose::{BestDoseProblem, Target, DoseRange}; +//! use pmcore::bestdose::{BestDosePosterior, Target, DoseRange}; //! //! # fn example(population_theta: pmcore::structs::theta::Theta, //! # population_weights: pmcore::structs::weights::Weights, //! # past_data: pharmsol::prelude::Subject, //! # target: pharmsol::prelude::Subject, //! # eq: pharmsol::prelude::ODE, -//! # error_models: pharmsol::prelude::ErrorModels, //! # settings: pmcore::routines::settings::Settings) //! # -> anyhow::Result<()> { -//! // Create optimization problem -//! let problem = BestDoseProblem::new( -//! &population_theta, // Population support points from NPAG -//! &population_weights, // Population probabilities +//! // Stage 1: Compute posterior from patient history +//! let posterior = BestDosePosterior::compute( +//! &population_theta, // Population support points from NPAG +//! &population_weights, // Population probabilities //! Some(past_data), // Patient history (None = use prior) +//! eq, // PK/PD model +//! settings, // NPAG settings +//! )?; +//! +//! // Stage 2 & 3: Optimize doses and get predictions +//! let result = posterior.optimize( //! target, // Future template with targets //! None, // time_offset (None = standard mode) -//! eq, // PK/PD model -//! error_models, // Error specifications //! DoseRange::new(0.0, 1000.0), // Dose constraints (0-1000 mg) //! 0.5, // bias_weight: 0=personalized, 1=population -//! settings, // NPAG settings //! Target::Concentration, // Target type //! )?; //! -//! // Run optimization -//! let result = problem.optimize()?; -//! //! // Extract results -//! println!("Optimal dose: {:?} mg", result.dose); -//! println!("Final cost: {}", result.objf); -//! println!("Method: {}", result.optimization_method); // "posterior" or "uniform" +//! println!("Optimal dose: {:?} mg", result.doses()); +//! println!("Final cost: {}", result.objf()); +//! println!("Method: {}", result.optimization_method()); //! # Ok(()) //! # } //! ``` @@ -144,32 +143,33 @@ //! ## Single Dose Optimization //! //! ```rust,no_run,ignore -//! use pmcore::bestdose::{BestDoseProblem, Target, DoseRange}; +//! use pmcore::bestdose::{BestDosePosterior, Target, DoseRange}; //! use pharmsol::prelude::Subject; //! //! # fn example(population_theta: pmcore::structs::theta::Theta, //! # population_weights: pmcore::structs::weights::Weights, //! # past: pharmsol::prelude::Subject, //! # eq: pharmsol::prelude::ODE, -//! # error_models: pharmsol::prelude::ErrorModels, //! # settings: pmcore::routines::settings::Settings) //! # -> anyhow::Result<()> { //! // Define target: 5 mg/L at 24 hours //! let target = Subject::builder("patient_001") -//! .bolus(0.0, 100.0, 0) // Initial dose (will be optimized) +//! .bolus(0.0, 0.0, 0) // Dose placeholder (will be optimized) //! .observation(24.0, 5.0, 0) // Target: 5 mg/L at 24h //! .build(); //! -//! let problem = BestDoseProblem::new( -//! &population_theta, &population_weights, Some(past), target, None, -//! eq, error_models, +//! let posterior = BestDosePosterior::compute( +//! &population_theta, &population_weights, Some(past), eq, settings, +//! )?; +//! +//! let result = posterior.optimize( +//! target, None, //! DoseRange::new(10.0, 500.0), // 10-500 mg allowed //! 0.3, // Slight population emphasis -//! settings, Target::Concentration, +//! Target::Concentration, //! )?; //! -//! let result = problem.optimize()?; -//! println!("Optimal dose: {} mg", result.dose[0]); +//! println!("Optimal dose: {} mg", result.doses()[0]); //! # Ok(()) //! # } //! ``` @@ -177,35 +177,36 @@ //! ## Multiple Doses with AUC Target //! //! ```rust,no_run,ignore -//! use pmcore::bestdose::{BestDoseProblem, Target, DoseRange}; +//! use pmcore::bestdose::{BestDosePosterior, Target, DoseRange}; //! use pharmsol::prelude::Subject; //! //! # fn example(population_theta: pmcore::structs::theta::Theta, //! # population_weights: pmcore::structs::weights::Weights, //! # past: pharmsol::prelude::Subject, //! # eq: pharmsol::prelude::ODE, -//! # error_models: pharmsol::prelude::ErrorModels, //! # settings: pmcore::routines::settings::Settings) //! # -> anyhow::Result<()> { //! // Target: Achieve AUC₂₄ = 400 mg·h/L //! let target = Subject::builder("patient_002") -//! .bolus(0.0, 100.0, 0) // Dose 1 (optimized) -//! .bolus(12.0, 100.0, 0) // Dose 2 (optimized) +//! .bolus(0.0, 0.0, 0) // Dose 1 placeholder (optimized) +//! .bolus(12.0, 0.0, 0) // Dose 2 placeholder (optimized) //! .observation(24.0, 400.0, 0) // Target: AUC₂₄ = 400 //! .build(); //! -//! let problem = BestDoseProblem::new( -//! &population_theta, &population_weights, Some(past), target, None, -//! eq, error_models, +//! let posterior = BestDosePosterior::compute( +//! &population_theta, &population_weights, Some(past), eq, settings, +//! )?; +//! +//! let result = posterior.optimize( +//! target, None, //! DoseRange::new(50.0, 300.0), //! 0.0, // Full personalization -//! settings, Target::AUCFromZero, // Cumulative AUC target! +//! Target::AUCFromZero, // Cumulative AUC target! //! )?; //! -//! let result = problem.optimize()?; -//! println!("Dose 1: {} mg at t=0", result.dose[0]); -//! println!("Dose 2: {} mg at t=12", result.dose[1]); -//! if let Some(auc) = result.auc_predictions { +//! println!("Dose 1: {} mg at t=0", result.doses()[0]); +//! println!("Dose 2: {} mg at t=12", result.doses()[1]); +//! if let Some(auc) = result.auc_predictions() { //! println!("Predicted AUC₂₄: {} mg·h/L", auc[0].1); //! } //! # Ok(()) @@ -215,27 +216,26 @@ //! ## Population-Only Optimization //! //! ```rust,no_run,ignore -//! # use pmcore::bestdose::{BestDoseProblem, Target, DoseRange}; +//! # use pmcore::bestdose::{BestDosePosterior, Target, DoseRange}; //! # fn example(population_theta: pmcore::structs::theta::Theta, //! # population_weights: pmcore::structs::weights::Weights, //! # target: pharmsol::prelude::Subject, //! # eq: pharmsol::prelude::ODE, -//! # error_models: pharmsol::prelude::ErrorModels, //! # settings: pmcore::routines::settings::Settings) //! # -> anyhow::Result<()> { //! // No patient history - use population prior directly -//! let problem = BestDoseProblem::new( +//! let posterior = BestDosePosterior::compute( //! &population_theta, &population_weights, -//! None, // No past data -//! target, None, // time_offset -//! eq, error_models, +//! None, // No past data → use prior +//! eq, settings, +//! )?; +//! +//! let result = posterior.optimize( +//! target, None, //! DoseRange::new(0.0, 1000.0), //! 1.0, // Full population weighting -//! settings, //! Target::Concentration, //! )?; -//! -//! let result = problem.optimize()?; //! // Returns population-typical dose //! # Ok(()) //! # } @@ -261,67 +261,46 @@ //! - `Target::AUCFromZero`: Cumulative AUC from time 0 //! - `Target::AUCFromLastDose`: Interval AUC from last dose //! -//! ## Performance Tuning -//! -//! For faster optimization: -//! ```rust,no_run,ignore -//! # use pmcore::bestdose::{BestDoseProblem, Target, DoseRange}; -//! # fn example(population_theta: pmcore::structs::theta::Theta, -//! # population_weights: pmcore::structs::weights::Weights, -//! # target: pharmsol::prelude::Subject, -//! # eq: pharmsol::ODE, -//! # error_models: pharmsol::prelude::ErrorModels, -//! # mut settings: pmcore::routines::settings::Settings) -//! # -> anyhow::Result<()> { -//! // Reduce refinement cycles -//! let problem = BestDoseProblem::new( -//! &population_theta, &population_weights, None, target, None, -//! eq, error_models, -//! DoseRange::new(0.0, 1000.0), 0.5, -//! settings.clone(), -//! Target::Concentration, -//! )?; -//! -//! // For AUC: use coarser time grid -//! settings.predictions().idelta = 30.0; // 30-minute intervals -//! # Ok(()) -//! # } -//! ``` -//! //! # See Also //! -//! - [`BestDoseProblem`]: Main entry point for optimization +//! - [`BestDosePosterior`]: Two-stage API entry point (compute posterior, then optimize) //! - [`BestDoseResult`]: Output structure with optimal doses //! - [`Target`]: Enum for concentration vs AUC targets //! - [`DoseRange`]: Dose constraint specification -pub mod cost; +pub(crate) mod cost; mod optimization; mod posterior; -pub mod predictions; +pub(crate) mod predictions; mod types; // Re-export public API -pub use types::{BestDoseProblem, BestDoseResult, DoseRange, Target}; +pub use types::{ + BestDosePosterior, BestDoseResult, BestDoseStatus, DoseRange, OptimalMethod, Target, +}; /// Helper function to concatenate past and future subjects (Option 3: Fortran MAKETMP approach) /// /// This mimics Fortran's MAKETMP subroutine logic: /// 1. Takes doses (only doses, not observations) from past subject -/// 2. Offsets all future subject event times by `time_offset` +/// 2. Offsets all future subject event times by `effective_offset` (absolute) /// 3. Combines into single continuous subject /// +/// Note: This function receives the **effective** (absolute) offset, computed +/// by `optimize()` as `max_past_time + time_offset` where `time_offset` is the +/// user-facing gap parameter. +/// /// # Arguments /// /// * `past` - Subject with past history (only doses will be used) /// * `future` - Subject template for future (all events: doses + observations) -/// * `time_offset` - Time offset to apply to all future events +/// * `effective_offset` - Absolute time offset to apply to all future events /// /// # Returns /// /// Combined subject with: -/// - Past doses at original times [0, time_offset) -/// - Future doses + observations at offset times [time_offset, ∞) +/// - Past doses at original times [0, effective_offset) +/// - Future doses + observations at offset times [effective_offset, ∞) /// /// # Example /// @@ -329,24 +308,24 @@ pub use types::{BestDoseProblem, BestDoseResult, DoseRange, Target}; /// // Past: dose at t=0, observation at t=6 (patient has been on therapy 6 hours) /// let past = Subject::builder("patient") /// .bolus(0.0, 500.0, 0) -/// .observation(6.0, 15.0, 0) // 15 mg/L at 6 hours +/// .observation(6.0, 15.0, 0) // 15 mg/L at 6 hours (max_past_time = 6) /// .build(); /// /// // Future: dose at t=0 (relative), target at t=24 (relative) /// let future = Subject::builder("patient") -/// .bolus(0.0, 100.0, 0) // Dose to optimize, will be at t=6 absolute -/// .observation(24.0, 10.0, 0) // Target at t=30 absolute +/// .bolus(0.0, 100.0, 0) // At absolute t=6 (with gap=0) +/// .observation(24.0, 10.0, 0) // At absolute t=30 (with gap=0) /// .build(); /// -/// // Concatenate with time_offset = 6.0 +/// // effective_offset = max_past_time(6) + gap(0) = 6 /// let combined = concatenate_past_and_future(&past, &future, 6.0); -/// // Result: dose at t=0 (fixed, 500mg), dose at t=6 (optimizable, 100mg initial), +/// // Result: dose at t=0 (fixed, 500mg), dose at t=6 (optimizable), /// // observation target at t=30 (10 mg/L) /// ``` fn concatenate_past_and_future( past: &pharmsol::prelude::Subject, future: &pharmsol::prelude::Subject, - time_offset: f64, + effective_offset: f64, ) -> pharmsol::prelude::Subject { use pharmsol::prelude::*; @@ -370,17 +349,20 @@ fn concatenate_past_and_future( } } - // Add future events with time offset + // Add future events with effective offset for occasion in future.occasions() { for event in occasion.events() { match event { Event::Bolus(bolus) => { - builder = - builder.bolus(bolus.time() + time_offset, bolus.amount(), bolus.input()); + builder = builder.bolus( + bolus.time() + effective_offset, + bolus.amount(), + bolus.input(), + ); } Event::Infusion(inf) => { builder = builder.infusion( - inf.time() + time_offset, + inf.time() + effective_offset, inf.amount(), inf.input(), inf.duration(), @@ -389,7 +371,7 @@ fn concatenate_past_and_future( Event::Observation(obs) => { builder = match obs.value() { Some(val) => { - builder.observation(obs.time() + time_offset, val, obs.outeq()) + builder.observation(obs.time() + effective_offset, val, obs.outeq()) } None => builder, }; @@ -401,61 +383,6 @@ fn concatenate_past_and_future( builder.build() } -/// Calculate which doses are optimizable based on dose amounts -/// -/// Returns a boolean mask where: -/// - `true` = dose amount is 0 (placeholder, optimizable) -/// - `false` = dose amount > 0 (fixed past dose) -/// -/// This allows users to specify a combined subject with: -/// - Non-zero doses for past doses (e.g., 500 mg at t=0) - these are fixed -/// - Zero doses as placeholders for future doses (e.g., 0 mg at t=6) - these are optimized -/// -/// # Arguments -/// -/// * `subject` - The subject with both fixed and placeholder doses -/// -/// # Returns -/// -/// Vector of booleans, one per dose in the subject -/// -/// # Example -/// -/// ```rust,ignore -/// let subject = Subject::builder("patient") -/// .bolus(0.0, 500.0, 0) // Past dose (fixed) - mask[0] = false -/// .bolus(6.0, 0.0, 0) // Future dose (optimize) - mask[1] = true -/// .observation(30.0, 10.0, 0) -/// .build(); -/// let mask = calculate_dose_optimization_mask(&subject); -/// assert_eq!(mask, vec![false, true]); -/// ``` -fn calculate_dose_optimization_mask(subject: &pharmsol::prelude::Subject) -> Vec { - use pharmsol::prelude::*; - - let mut mask = Vec::new(); - - for occasion in subject.occasions() { - for event in occasion.events() { - match event { - Event::Bolus(bolus) => { - // Dose is optimizable if amount is 0 (placeholder) - mask.push(bolus.amount() == 0.0); - } - Event::Infusion(infusion) => { - // Infusion is optimizable if amount is 0 (placeholder) - mask.push(infusion.amount() == 0.0); - } - Event::Observation(_) => { - // Observations don't go in the mask - } - } - } - } - - mask -} - use anyhow::Result; use pharmsol::prelude::*; use pharmsol::ODE; @@ -464,36 +391,230 @@ use crate::routines::settings::Settings; use crate::structs::theta::Theta; use crate::structs::weights::Weights; +use types::BestDoseProblem; + // ═════════════════════════════════════════════════════════════════════════════ -// Helper Functions for STAGE 1: Posterior Density Calculation +// BestDosePosterior: Public two-stage API // ═════════════════════════════════════════════════════════════════════════════ -/// Validate time_offset parameter for past/future separation mode -fn validate_time_offset(time_offset: f64, past_data: &Option) -> Result<()> { - if let Some(past_subject) = past_data { - let max_past_time = past_subject +impl BestDosePosterior { + /// **Stage 1**: Compute the Bayesian posterior density from population prior and patient data + /// + /// This performs the expensive posterior calculation (NPAGFULL11 filtering + NPAGFULL refinement) + /// and returns a reusable `BestDosePosterior` that can be optimized multiple times. + /// + /// # Algorithm + /// + /// ```text + /// Prior (N support points) + /// ↓ + /// NPAGFULL11: Bayesian filtering + /// P(θᵢ|data) ∝ P(data|θᵢ) × P(θᵢ) + /// ↓ + /// Filtered posterior (M points) + /// ↓ + /// NPAGFULL: Local refinement (max_cycles iterations) + /// ↓ + /// Refined posterior (M points with updated weights) + /// ``` + /// + /// # Arguments + /// + /// * `population_theta` - Population support points from NPAG + /// * `population_weights` - Population probabilities + /// * `past_data` - Patient history (`None` = use prior directly) + /// * `eq` - Pharmacokinetic/pharmacodynamic model + /// * `settings` - NPAG settings (includes error models and posterior refinement config) + /// + /// # Example + /// + /// ```rust,no_run,ignore + /// let posterior = BestDosePosterior::compute( + /// &theta, &weights, + /// Some(past_subject), + /// eq, settings, + /// )?; + /// println!("Posterior has {} support points", posterior.n_support_points()); + /// ``` + pub fn compute( + population_theta: &Theta, + population_weights: &Weights, + past_data: Option, + eq: ODE, + settings: Settings, + ) -> Result { + tracing::info!("╔══════════════════════════════════════════════════════════╗"); + tracing::info!("║ BestDose Algorithm: STAGE 1 ║"); + tracing::info!("║ Posterior Density Calculation ║"); + tracing::info!("╚══════════════════════════════════════════════════════════╝"); + + let (posterior_theta, posterior_weights, filtered_population_weights, _past_subject) = + calculate_posterior_density( + population_theta, + population_weights, + past_data.as_ref(), + &eq, + &settings.errormodels, + &settings, + )?; + + tracing::info!("╔══════════════════════════════════════════════════════════╗"); + tracing::info!("║ Stage 1 Complete - Posterior Ready ║"); + tracing::info!("╚══════════════════════════════════════════════════════════╝"); + tracing::info!(" Support points: {}", posterior_theta.matrix().nrows()); + + Ok(BestDosePosterior { + theta: posterior_theta, + posterior: posterior_weights, + population_weights: filtered_population_weights, + past_data, + eq, + settings, + }) + } + + /// **Stage 2**: Optimize doses for target outcomes using the computed posterior + /// + /// This runs the dual optimization (posterior weights vs uniform weights) and + /// returns the best dosing regimen. Can be called multiple times on the same + /// posterior with different parameters. + /// + /// # Arguments + /// + /// * `target` - Future dosing template with target observations + /// * `time_offset` - Optional gap (in hours) between the last past event and the start of + /// the future target. 0 means the future starts immediately after the last past event. + /// The effective absolute offset is `max_past_time + time_offset`. + /// * `dose_range` - Allowable dose constraints + /// * `bias_weight` - λ in \[0,1\]: 0=personalized, 1=population + /// * `target_type` - Concentration or AUC targets + /// + /// # Example + /// + /// ```rust,no_run,ignore + /// // Try different bias weights + /// for &bw in &[0.0, 0.25, 0.5, 0.75, 1.0] { + /// let result = posterior.optimize( + /// target.clone(), + /// None, + /// DoseRange::new(0.0, 300.0), + /// bw, + /// Target::Concentration, + /// )?; + /// println!("λ={}: dose={:.1}", bw, result.doses()[0]); + /// } + /// ``` + pub fn optimize( + &self, + target: Subject, + time_offset: Option, + dose_range: DoseRange, + bias_weight: f64, + target_type: Target, + ) -> Result { + tracing::info!("╔══════════════════════════════════════════════════════════╗"); + tracing::info!("║ BestDose Algorithm: STAGE 2 & 3 ║"); + tracing::info!("║ Dual Optimization + Final Predictions ║"); + tracing::info!("╚══════════════════════════════════════════════════════════╝"); + tracing::info!(" Target type: {:?}", target_type); + tracing::info!(" Bias weight (λ): {}", bias_weight); + + // Validate and compute effective time_offset + // time_offset is a gap relative to the last past event: + // effective_offset = max_past_time + time_offset + // So time_offset=0 means "future starts right after last past event" + if let Some(t) = time_offset { + if t < 0.0 { + return Err(anyhow::anyhow!( + "Invalid time_offset: {} is negative. \ + time_offset must be >= 0 (it represents the gap after the last past event).", + t + )); + } + } + + // Compute the absolute offset for concatenation + let effective_offset = time_offset.map(|t| { + let max_past_time = self + .past_data + .as_ref() + .map(|past| { + past.occasions() + .iter() + .flat_map(|occ| occ.events()) + .map(|event| match event { + Event::Bolus(b) => b.time(), + Event::Infusion(i) => i.time(), + Event::Observation(o) => o.time(), + }) + .fold(0.0_f64, |max, time| max.max(time)) + }) + .unwrap_or(0.0); + max_past_time + t + }); + + // Handle past/future concatenation if needed + // When time_offset is provided, offset all target event times by the + // effective offset (max_past_time + gap) and prepend past doses. + let final_target = match effective_offset { + None => target, + Some(eff) => { + tracing::info!( + " Time offset gap: {:?} hours (effective absolute offset: {} hours)", + time_offset, + eff + ); + match &self.past_data { + Some(past) => { + tracing::info!(" Concatenating past doses with offset target events"); + concatenate_past_and_future(past, &target, eff) + } + None => { + tracing::info!(" No past data stored — offsetting target events only"); + concatenate_past_and_future( + &Subject::builder("empty").build(), + &target, + eff, + ) + } + } + } + }; + + // Validate that the target has observations + let has_observations = final_target .occasions() .iter() .flat_map(|occ| occ.events()) - .map(|event| match event { - Event::Bolus(b) => b.time(), - Event::Infusion(i) => i.time(), - Event::Observation(o) => o.time(), - }) - .fold(0.0_f64, |max, time| max.max(time)); - - if time_offset < max_past_time { + .any(|event| matches!(event, Event::Observation(_))); + if !has_observations { return Err(anyhow::anyhow!( - "Invalid time_offset: {} is before the last past_data event at time {}. \ - time_offset must be >= the maximum time in past_data to avoid time travel!", - time_offset, - max_past_time + "Target subject has no observations. At least one observation is required for dose optimization." )); } + + // Build the internal optimization problem + let problem = BestDoseProblem { + target: final_target, + target_type, + population_weights: self.population_weights.clone(), + theta: self.theta.clone(), + posterior: self.posterior.clone(), + eq: self.eq.clone(), + settings: self.settings.clone(), + doserange: dose_range, + bias_weight, + }; + + // Run dual optimization + final predictions + optimization::dual_optimization(&problem) } - Ok(()) } +// ═════════════════════════════════════════════════════════════════════════════ +// Helper Functions for STAGE 1: Posterior Density Calculation +// ═════════════════════════════════════════════════════════════════════════════ + /// Calculate posterior density (STAGE 1: Two-step process) /// /// # Algorithm Flow (Matches Diagram) @@ -583,236 +704,3 @@ fn calculate_posterior_density( } } } - -/// Prepare target subject by handling past/future concatenation if needed -/// -/// # Returns -/// -/// Tuple: (final_target, final_past_data) -fn prepare_target_subject( - past_subject: Subject, - target: Subject, - time_offset: Option, -) -> Result<(Subject, Subject)> { - match time_offset { - None => { - tracing::info!(" Mode: Standard (single subject)"); - Ok((target, past_subject)) - } - Some(t) => { - tracing::info!(" Mode: Past/Future separation (Fortran MAKETMP approach)"); - tracing::info!(" Current time boundary: {} hours", t); - tracing::info!(" Concatenating past and future subjects..."); - - let combined = concatenate_past_and_future(&past_subject, &target, t); - - // Log dose structure - let mask = calculate_dose_optimization_mask(&combined); - let num_fixed = mask.iter().filter(|&&x| !x).count(); - let num_optimizable = mask.iter().filter(|&&x| x).count(); - tracing::info!(" Fixed doses (from past): {}", num_fixed); - tracing::info!(" Optimizable doses (from future): {}", num_optimizable); - - Ok((combined, past_subject)) - } - } -} - -// ═════════════════════════════════════════════════════════════════════════════ - -impl BestDoseProblem { - /// Create a new BestDose problem with automatic posterior calculation - /// - /// This is the main entry point for the BestDose algorithm. - /// - /// # Algorithm Structure (Matches Flowchart) - /// - /// ```text - /// ┌─────────────────────────────────────────┐ - /// │ STAGE 1: Posterior Density Calculation │ - /// │ │ - /// │ Prior Density (N points) │ - /// │ ↓ │ - /// │ Has past data with observations? │ - /// │ ↓ Yes ↓ No │ - /// │ Step 1.1: Use prior │ - /// │ NPAGFULL11 directly │ - /// │ (Filter) │ - /// │ ↓ │ - /// │ Step 1.2: │ - /// │ NPAGFULL │ - /// │ (Refine) │ - /// │ ↓ │ - /// │ Posterior Density │ - /// └─────────────────────────────────────────┘ - /// ``` - /// - /// # Parameters - /// - /// * `population_theta` - Population support points from NPAG - /// * `population_weights` - Population probabilities - /// * `past_data` - Patient history (None = use prior directly) - /// * `target` - Future dosing template with targets - /// * `time_offset` - Optional time offset for concatenation (None = standard mode, Some(t) = Fortran mode) - /// * `eq` - Pharmacokinetic/pharmacodynamic model - /// * `error_models` - Error model specifications - /// * `doserange` - Allowable dose constraints - /// * `bias_weight` - λ ∈ [0,1]: 0=personalized, 1=population - /// * `settings` - NPAG settings for posterior refinement - /// * `max_cycles` - NPAGFULL cycles (0=skip refinement, 500=default) - /// * `target_type` - Concentration or AUC targets - /// - /// # Returns - /// - /// BestDoseProblem ready for `optimize()` - #[allow(clippy::too_many_arguments)] - pub fn new( - population_theta: &Theta, - population_weights: &Weights, - past_data: Option, - target: Subject, - time_offset: Option, - eq: ODE, - doserange: DoseRange, - bias_weight: f64, - settings: Settings, - target_type: Target, - ) -> Result { - tracing::info!("╔══════════════════════════════════════════════════════════╗"); - tracing::info!("║ BestDose Algorithm: STAGE 1 ║"); - tracing::info!("║ Posterior Density Calculation ║"); - tracing::info!("╚══════════════════════════════════════════════════════════╝"); - - // Validate input if using past/future separation mode - if let Some(t) = time_offset { - validate_time_offset(t, &past_data)?; - } - - // ═════════════════════════════════════════════════════════════ - // STAGE 1: Calculate Posterior Density - // ═════════════════════════════════════════════════════════════ - let (posterior_theta, posterior_weights, filtered_population_weights, past_subject) = - calculate_posterior_density( - population_theta, - population_weights, - past_data.as_ref(), - &eq, - &settings.errormodels, - &settings, - )?; - - // Handle past/future concatenation if needed - let (final_target, _) = prepare_target_subject(past_subject, target, time_offset)?; - - tracing::info!("╔══════════════════════════════════════════════════════════╗"); - tracing::info!("║ Stage 1 Complete - Ready for Optimization ║"); - tracing::info!("╚══════════════════════════════════════════════════════════╝"); - tracing::info!(" Support points: {}", posterior_theta.matrix().nrows()); - tracing::info!(" Target type: {:?}", target_type); - tracing::info!(" Bias weight (λ): {}", bias_weight); - - Ok(BestDoseProblem { - target: final_target, - target_type, - population_weights: filtered_population_weights, - theta: posterior_theta, - posterior: posterior_weights, - eq, - settings, - doserange, - bias_weight, - }) - } - - /// Run the complete BestDose optimization algorithm - /// - /// # Algorithm Flow (Matches Diagram!) - /// - /// ```text - /// ┌─────────────────────────────────────────┐ - /// │ STAGE 1: Posterior Calculation │ - /// │ [COMPLETED in new()] │ - /// └────────────┬────────────────────────────┘ - /// ↓ - /// ┌─────────────────────────────────────────┐ - /// │ STAGE 2: Dual Optimization │ - /// │ │ - /// │ Optimization 1: Posterior Weights │ - /// │ (Patient-specific) │ - /// │ ↓ │ - /// │ Result 1: (doses₁, cost₁) │ - /// │ │ - /// │ Optimization 2: Uniform Weights │ - /// │ (Population-based) │ - /// │ ↓ │ - /// │ Result 2: (doses₂, cost₂) │ - /// │ │ - /// │ Select: min(cost₁, cost₂) │ - /// └────────────┬────────────────────────────┘ - /// ↓ - /// ┌─────────────────────────────────────────┐ - /// │ STAGE 3: Final Predictions │ - /// │ │ - /// │ Calculate predictions with │ - /// │ optimal doses and winning weights │ - /// └─────────────────────────────────────────┘ - /// ``` - /// - /// # Returns - /// - /// `BestDoseResult` containing: - /// - `dose`: Optimal dose amount(s) - /// - `objf`: Final cost function value - /// - `preds`: Concentration-time predictions - /// - `auc_predictions`: AUC values (if target_type is AUC) - /// - `optimization_method`: "posterior" or "uniform" - pub fn optimize(self) -> Result { - tracing::info!("╔══════════════════════════════════════════════════════════╗"); - tracing::info!("║ BestDose Algorithm: STAGE 2 & 3 ║"); - tracing::info!("║ Dual Optimization + Final Predictions ║"); - tracing::info!("╚══════════════════════════════════════════════════════════╝"); - - // STAGE 2 & 3: Dual optimization + predictions - optimization::dual_optimization(&self) - } - - /// Set the bias weight (lambda parameter) - /// - /// - λ = 0.0 (default): Full personalization (minimize patient-specific variance) - /// - λ = 0.5: Balanced between individual and population - /// - λ = 1.0: Population-based (minimize deviation from population mean) - pub fn with_bias_weight(mut self, weight: f64) -> Self { - self.bias_weight = weight; - self - } - - /// Get a reference to the refined posterior support points (Θ) - pub fn posterior_theta(&self) -> &Theta { - &self.theta - } - - /// Get the posterior probability weights - pub fn posterior_weights(&self) -> &Weights { - &self.posterior - } - - /// Get the filtered population weights used for the bias term - pub fn population_weights(&self) -> &Weights { - &self.population_weights - } - - /// Get the prepared target subject - pub fn target_subject(&self) -> &Subject { - &self.target - } - - /// Get the currently configured bias weight (λ) - pub fn bias_weight(&self) -> f64 { - self.bias_weight - } - - /// Get the selected optimization target type - pub fn target_type(&self) -> Target { - self.target_type - } -} diff --git a/src/bestdose/predictions.rs b/src/bestdose/predictions.rs index 740a06690..369c1065a 100644 --- a/src/bestdose/predictions.rs +++ b/src/bestdose/predictions.rs @@ -267,7 +267,8 @@ pub fn calculate_interval_auc_per_observation( /// /// This generates the final NPPredictions structure with the optimal doses /// and appropriate weights (posterior or uniform depending on which optimization won). -pub fn calculate_final_predictions( +#[allow(clippy::type_complexity)] +pub(crate) fn calculate_final_predictions( problem: &BestDoseProblem, optimal_doses: &[f64], weights: &Weights, diff --git a/src/bestdose/types.rs b/src/bestdose/types.rs index e422cd0be..1e2a943c5 100644 --- a/src/bestdose/types.rs +++ b/src/bestdose/types.rs @@ -1,7 +1,7 @@ //! Core data types for the BestDose algorithm //! //! This module defines the main structures used throughout the BestDose optimization: -//! - [`BestDoseProblem`]: The complete optimization problem specification +//! - [`BestDosePosterior`]: Two-stage API entry point — compute posterior, then optimize //! - [`BestDoseResult`]: Output structure containing optimal doses and predictions //! - [`Target`]: Enum specifying concentration or AUC targets //! - [`DoseRange`]: Dose constraint specification @@ -183,135 +183,131 @@ impl Default for DoseRange { } } -/// The BestDose optimization problem +/// The computed Bayesian posterior for a patient /// -/// Contains all data needed for the three-stage BestDose algorithm. -/// Create via [`BestDoseProblem::new()`], then call [`.optimize()`](BestDoseProblem::optimize) -/// to run the full algorithm. +/// This is the main public entry point for the two-stage BestDose API: /// -/// # Three-Stage Algorithm -/// -/// 1. **Posterior Density Calculation** (automatic in `new()`) +/// 1. **Stage 1: Posterior computation** ([`BestDosePosterior::compute()`]) /// - NPAGFULL11: Bayesian filtering of prior support points /// - NPAGFULL: Local refinement of each filtered point /// -/// 2. **Dual Optimization** (automatic in `optimize()`) -/// - Optimization with posterior weights (patient-specific) -/// - Optimization with uniform weights (population-based) -/// - Selection of better result -/// -/// 3. **Final Predictions** (automatic in `optimize()`) -/// - Concentration or AUC predictions with optimal doses -/// -/// # Fields -/// -/// ## Input Data -/// - `target`: Future dosing template with target observations -/// - `target_type`: [`Target::Concentration`] or [`Target::AUC`] -/// -/// ## Population Prior -/// - `population_weights`: Filtered population probability weights (used for bias term) +/// 2. **Stage 2: Dose optimization** ([`BestDosePosterior::optimize()`]) +/// - Dual optimization (posterior vs uniform weights) +/// - Final predictions with optimal doses /// -/// ## Patient-Specific Posterior -/// - `theta`: Refined posterior support points (from NPAGFULL11 + NPAGFULL) -/// - `posterior`: Posterior probability weights -/// -/// ## Model Components -/// - `eq`: Pharmacokinetic/pharmacodynamic ODE model -/// - `settings`: NPAG configuration settings (used for prediction grid) -/// -/// ## Optimization Parameters -/// - `doserange`: Min/max dose constraints -/// - `bias_weight` (λ): Personalization parameter (0=personalized, 1=population) +/// The posterior can be reused across multiple `optimize()` calls with +/// different targets, dose ranges, or bias weights. /// /// # Example /// /// ```rust,no_run,ignore -/// use pmcore::bestdose::{BestDoseProblem, Target, DoseRange}; +/// use pmcore::bestdose::{BestDosePosterior, Target, DoseRange}; /// /// # fn example(population_theta: pmcore::structs::theta::Theta, /// # population_weights: pmcore::structs::weights::Weights, /// # past: pharmsol::prelude::Subject, /// # target: pharmsol::prelude::Subject, /// # eq: pharmsol::prelude::ODE, -/// # error_models: pharmsol::prelude::ErrorModels, /// # settings: pmcore::routines::settings::Settings) /// # -> anyhow::Result<()> { -/// let problem = BestDoseProblem::new( +/// // Stage 1: Compute posterior (expensive, done once) +/// let posterior = BestDosePosterior::compute( /// &population_theta, /// &population_weights, -/// Some(past), // Patient history -/// target, // Dosing template with targets +/// Some(past), /// eq, -/// error_models, -/// DoseRange::new(0.0, 1000.0), -/// 0.5, // Balanced personalization /// settings, -/// 500, // NPAGFULL cycles -/// Target::Concentration, /// )?; /// -/// let result = problem.optimize()?; +/// // Stage 2: Optimize doses (can be called multiple times) +/// let result = posterior.optimize( +/// target, +/// None, // No time offset +/// DoseRange::new(0.0, 1000.0), +/// 0.5, // bias_weight +/// Target::Concentration, +/// )?; /// # Ok(()) /// # } /// ``` #[derive(Debug, Clone)] -pub struct BestDoseProblem { - /// Target subject with dosing template and target observations - /// - /// This [Subject] defines the targets for optimization, including - /// dose events (with amounts to be optimized) and observation events - /// (with desired target values). - /// - /// For a `Target::Concentration`, observation values are target concentrations. - /// For a `Target::AUC`, observation values are target cumulative AUC. - /// - /// Only doses with a value of `0.0` will be optimized; non-zero doses remain fixed. +pub struct BestDosePosterior { + /// Refined posterior support points (from NPAGFULL11 + NPAGFULL) + pub(crate) theta: Theta, + /// Posterior probability weights + pub(crate) posterior: Weights, + /// Filtered population weights (used for bias term in cost function) + pub(crate) population_weights: Weights, + /// Past patient data (stored for use in optimize() with time_offset) + pub(crate) past_data: Option, + /// PK/PD model + pub(crate) eq: ODE, + /// Settings (used for prediction grid, error models, etc.) + pub(crate) settings: Settings, +} + +impl BestDosePosterior { + /// Get the refined posterior support points (Θ) + pub fn theta(&self) -> &Theta { + &self.theta + } + + /// Get the posterior probability weights + pub fn posterior_weights(&self) -> &Weights { + &self.posterior + } + + /// Get the filtered population weights used for the bias term + pub fn population_weights(&self) -> &Weights { + &self.population_weights + } + + /// Get the number of support points in the posterior + pub fn n_support_points(&self) -> usize { + self.theta.matrix().nrows() + } +} + +/// Internal optimization problem (not exposed in public API) +/// +/// Contains all data needed for dose optimization. +/// Created internally by [`BestDosePosterior::optimize()`]. +#[derive(Debug, Clone)] +pub(crate) struct BestDoseProblem { pub(crate) target: Subject, - /// Target type for optimization - /// - /// Specifies whether to optimize for concentrations or AUC values. pub(crate) target_type: Target, - - /// The population prior weights ([Weights]), representing the probability of each support point in the population. pub(crate) population_weights: Weights, - - // Patient-specific posterior (from NPAGFULL11 + NPAGFULL) pub(crate) theta: Theta, pub(crate) posterior: Weights, - - // Model and settings pub(crate) eq: ODE, pub(crate) settings: Settings, - - // Optimization parameters pub(crate) doserange: DoseRange, - pub(crate) bias_weight: f64, // λ: 0=personalized, 1=population + pub(crate) bias_weight: f64, } /// Result from BestDose optimization /// /// Contains the optimal doses and associated predictions from running -/// [`BestDoseProblem::optimize()`]. +/// [`BestDosePosterior::optimize()`]. /// /// # Fields /// -/// - `dose`: Optimal dose amount(s) in the same order as doses in target subject +/// - `doses`: Optimal dose amount(s) in the same order as doses in target subject /// - `objf`: Final cost function value at optimal doses -/// - `status`: Optimization status message (e.g., "converged", "max iterations") -/// - `preds`: Concentration-time predictions using optimal doses -/// - `auc_predictions`: AUC values at observation times (only for [`Target::AUC`]) -/// - `optimization_method`: Which method won: `"posterior"` or `"uniform"` +/// - `status`: Optimization status (converged or max iterations) +/// - `predictions`: Concentration-time predictions using optimal doses +/// - `auc_predictions`: AUC values at observation times (only for AUC targets) +/// - `optimization_method`: Which method won: `Posterior` or `Uniform` /// /// # Interpretation /// /// ## Optimization Method /// -/// - **"posterior"**: Patient-specific optimization won (uses posterior weights) +/// - **`Posterior`**: Patient-specific optimization won (uses posterior weights) /// - Indicates patient differs from population or has sufficient history /// - Doses are highly personalized /// -/// - **"uniform"**: Population-based optimization won (uses uniform weights) +/// - **`Uniform`**: Population-based optimization won (uses uniform weights) /// - Indicates patient is population-typical or has limited history /// - Doses are more conservative/robust /// @@ -327,32 +323,26 @@ pub struct BestDoseProblem { /// ## Extracting Results /// /// ```rust,no_run,ignore -/// # use pmcore::bestdose::BestDoseProblem; -/// # fn example(problem: BestDoseProblem) -> anyhow::Result<()> { -/// let result = problem.optimize()?; +/// # use pmcore::bestdose::{BestDosePosterior, Target, DoseRange, BestDoseResult}; +/// # fn example(posterior: BestDosePosterior, +/// # target: pharmsol::prelude::Subject) -> anyhow::Result<()> { +/// let result = posterior.optimize( +/// target, None, DoseRange::new(0.0, 1000.0), 0.5, Target::Concentration, +/// )?; /// /// // Single dose -/// println!("Optimal dose: {} mg", result.dose[0]); +/// println!("Optimal dose: {} mg", result.doses()[0]); /// /// // Multiple doses -/// for (i, &dose) in result.dose.iter().enumerate() { +/// for (i, dose) in result.doses().iter().enumerate() { /// println!("Dose {}: {} mg", i + 1, dose); /// } /// /// // Check which method was used -/// match result.optimization_method.as_str() { -/// "posterior" => println!("Patient-specific optimization"), -/// "uniform" => println!("Population-based optimization"), -/// _ => {} -/// } -/// -/// // Access predictions -/// for pred in result.preds.iter() { -/// println!("t={:.1}h: {:.2} mg/L", pred.time(), pred.prediction()); -/// } +/// println!("Method: {}", result.optimization_method()); /// /// // For AUC targets -/// if let Some(auc_values) = result.auc_predictions { +/// if let Some(auc_values) = result.auc_predictions() { /// for (time, auc) in auc_values { /// println!("AUC at t={:.1}h: {:.1} mg·h/L", time, auc); /// } diff --git a/tests/bestdose_tests.rs b/tests/bestdose_tests.rs index aa1202824..145bf11ae 100644 --- a/tests/bestdose_tests.rs +++ b/tests/bestdose_tests.rs @@ -1,5 +1,5 @@ use anyhow::Result; -use pmcore::bestdose::{BestDoseProblem, DoseRange, Target}; +use pmcore::bestdose::{BestDosePosterior, DoseRange, Target}; use pmcore::prelude::*; use pmcore::structs::theta::Theta; use pmcore::structs::weights::Weights; @@ -58,18 +58,13 @@ fn test_infusion_mask_inclusion() -> Result<()> { }; let prior_weights = Weights::uniform(1); - // Create BestDose problem - let problem = BestDoseProblem::new( + // Create BestDose posterior + let posterior = BestDosePosterior::compute( &prior_theta, &prior_weights, None, - target.clone(), - None, eq.clone(), - DoseRange::new(10.0, 300.0), - 0.5, settings.clone(), - Target::Concentration, )?; // Count optimizable doses in the target @@ -90,7 +85,13 @@ fn test_infusion_mask_inclusion() -> Result<()> { ); // Run optimization - it should not panic and should handle infusion - let result = problem.optimize(); + let result = posterior.optimize( + target.clone(), + None, + DoseRange::new(10.0, 300.0), + 0.5, + Target::Concentration, + ); // The optimization should succeed assert!( @@ -178,22 +179,24 @@ fn test_fixed_infusion_preservation() -> Result<()> { let prior_weights = Weights::uniform(1); // Use current_time to separate past and future - let problem = BestDoseProblem::new( + let posterior = BestDosePosterior::compute( &prior_theta, &prior_weights, Some(past), - target, - Some(2.0), // Current time = 2.0 hours eq.clone(), + settings.clone(), + )?; + + let result = posterior.optimize( + target, + Some(0.0), // No gap after past (past ends at t=2.0) DoseRange::new(0.0, 500.0), 0.5, - settings.clone(), Target::Concentration, )?; - let result = problem.optimize()?; - - // Should only optimize the future bolus, not the past infusion + // With time_offset, past doses are concatenated with future target. + // Result should have 2 doses: fixed past infusion + optimized future bolus. let doses = result.doses(); eprintln!("Optimized doses: {:?}", doses); assert_eq!( @@ -201,7 +204,11 @@ fn test_fixed_infusion_preservation() -> Result<()> { 2, "Should have 2 doses (past infusion + future bolus)" ); - assert_eq!(doses[0], 200.0, "Past infusion dose should be preserved"); + assert!( + (doses[0] - 200.0).abs() < 1e-6, + "Past infusion should remain fixed at 200.0, got {}", + doses[0] + ); assert!(doses[1] > 0.0, "Future bolus dose should be optimized"); Ok(()) @@ -210,8 +217,6 @@ fn test_fixed_infusion_preservation() -> Result<()> { /// Test that dose count validation works #[test] fn test_dose_count_validation() -> Result<()> { - use pmcore::bestdose::cost::calculate_cost; - let eq = equation::ODE::new( |x, p, _t, dx, b, _rateiv, _cov| { fetch_params!(p, ke, _v); @@ -258,30 +263,23 @@ fn test_dose_count_validation() -> Result<()> { }; let prior_weights = Weights::uniform(1); - let problem = BestDoseProblem::new( - &prior_theta, - &prior_weights, - None, + let posterior = BestDosePosterior::compute(&prior_theta, &prior_weights, None, eq, settings)?; + + // Optimize with the correct target (2 optimizable doses, 2 observations) - should succeed + let result = posterior.optimize( target, None, - eq, DoseRange::new(10.0, 300.0), 0.5, - settings, Target::Concentration, - )?; - - // Try with wrong number of doses - should fail - let result_wrong = calculate_cost(&problem, &[100.0]); // Only 1 dose, need 2 - assert!(result_wrong.is_err(), "Should fail with wrong dose count"); - assert!(result_wrong.unwrap_err().to_string().contains("mismatch")); - - // Try with correct number of doses - should succeed - let result_correct = calculate_cost(&problem, &[100.0, 150.0]); + ); assert!( - result_correct.is_ok(), - "Should succeed with correct dose count" + result.is_ok(), + "Should succeed with correct target: {:?}", + result.err() ); + let result = result?; + assert_eq!(result.doses().len(), 2, "Should have 2 optimized doses"); Ok(()) } @@ -289,8 +287,6 @@ fn test_dose_count_validation() -> Result<()> { /// Test that empty observations are caught #[test] fn test_empty_observations_validation() -> Result<()> { - use pmcore::bestdose::cost::calculate_cost; - let eq = equation::ODE::new( |x, p, _t, dx, b, _rateiv, _cov| { fetch_params!(p, ke, _v); @@ -332,21 +328,16 @@ fn test_empty_observations_validation() -> Result<()> { }; let prior_weights = Weights::uniform(1); - let problem = BestDoseProblem::new( - &prior_theta, - &prior_weights, - None, + let posterior = BestDosePosterior::compute(&prior_theta, &prior_weights, None, eq, settings)?; + + // Try to optimize - should fail with no observations + let result = posterior.optimize( target, None, - eq, DoseRange::new(10.0, 300.0), 0.5, - settings, Target::Concentration, - )?; - - // Try to calculate cost - should fail with no observations - let result = calculate_cost(&problem, &[100.0]); + ); assert!(result.is_err(), "Should fail with no observations"); assert!(result.unwrap_err().to_string().contains("no observations")); @@ -402,20 +393,15 @@ fn test_basic_auc_mode() -> Result<()> { }; let prior_weights = Weights::uniform(1); - let problem = BestDoseProblem::new( - &prior_theta, - &prior_weights, - None, + let posterior = BestDosePosterior::compute(&prior_theta, &prior_weights, None, eq, settings)?; + + let result = posterior.optimize( target, None, - eq, DoseRange::new(100.0, 2000.0), 0.8, - settings, Target::AUCFromZero, - )?; - - let result = problem.optimize(); + ); assert!( result.is_ok(), @@ -494,22 +480,17 @@ fn test_infusion_auc_mode() -> Result<()> { }; let prior_weights = Weights::uniform(1); - // Create BestDose problem in AUC mode - let problem = BestDoseProblem::new( - &prior_theta, - &prior_weights, - None, + // Create BestDose posterior and optimize in AUC mode + let posterior = BestDosePosterior::compute(&prior_theta, &prior_weights, None, eq, settings)?; + + // Run optimization + let result = posterior.optimize( target, None, - eq, DoseRange::new(100.0, 2000.0), - 0.8, // Higher bias weight typically works better for AUC targets - settings, + 0.8, // Higher bias weight typically works better for AUC targets Target::AUCFromZero, // AUC mode! - )?; - - // Run optimization - let result = problem.optimize(); + ); assert!( result.is_ok(), @@ -584,6 +565,7 @@ fn test_multi_outeq_auc_mode() -> Result<()> { settings.disable_output(); settings.set_cycles(0); + settings.set_idelta(30.0); // 30-minute intervals for AUC calculation // Subject with fixed dose and target observations at multiple outeqs let target = Subject::builder("test") @@ -603,22 +585,18 @@ fn test_multi_outeq_auc_mode() -> Result<()> { }; let prior_weights = Weights::uniform(1); - let _problem = BestDoseProblem::new( - &prior_theta, - &prior_weights, - None, + let posterior = BestDosePosterior::compute(&prior_theta, &prior_weights, None, eq, settings)?; + + let _result = posterior.optimize( target, None, - eq, DoseRange::new(0.0, 2000.0), 0.5, - settings, Target::AUCFromZero, )?; - // Just verify that problem was created successfully + // Just verify that posterior compute and optimize succeed // This tests that cost calculation works with multi-outeq - // (cost is calculated during problem validation) Ok(()) } @@ -672,20 +650,15 @@ fn test_multi_outeq_auc_optimization() -> Result<()> { }; let prior_weights = Weights::uniform(1); - let problem = BestDoseProblem::new( - &prior_theta, - &prior_weights, - None, + let posterior = BestDosePosterior::compute(&prior_theta, &prior_weights, None, eq, settings)?; + + let result = posterior.optimize( target, None, - eq, DoseRange::new(0.0, 2000.0), 0.5, - settings, Target::AUCFromZero, - )?; - - let result = problem.optimize(); + ); assert!( result.is_ok(), "Multi-outeq AUC optimization failed: {:?}", @@ -765,21 +738,16 @@ fn test_auc_from_zero_single_dose() -> Result<()> { }; let prior_weights = Weights::uniform(1); - let problem = BestDoseProblem::new( - &prior_theta, - &prior_weights, - None, + let posterior = BestDosePosterior::compute(&prior_theta, &prior_weights, None, eq, settings)?; + + let result = posterior.optimize( target, None, - eq, DoseRange::new(100.0, 1000.0), 0.8, - settings, Target::AUCFromZero, // Cumulative AUC from time 0 )?; - let result = problem.optimize()?; - let doses: Vec = result.doses(); // Verify we got a result @@ -856,20 +824,15 @@ fn test_auc_from_last_dose_maintenance() -> Result<()> { }; let prior_weights = Weights::uniform(1); - let problem = BestDoseProblem::new( - &prior_theta, - &prior_weights, - None, + let posterior = BestDosePosterior::compute(&prior_theta, &prior_weights, None, eq, settings)?; + + let result = posterior.optimize( target, None, - eq, DoseRange::new(50.0, 500.0), 0.8, - settings, Target::AUCFromLastDose, // Interval AUC from last dose )?; - - let result = problem.optimize()?; let doses = result.doses(); // Verify we got a result @@ -951,20 +914,21 @@ fn test_auc_modes_comparison() -> Result<()> { .observation(24.0, 100.0, 0) // Target: AUC₀₋₂₄ = 100 .build(); - let problem_zero = BestDoseProblem::new( + let posterior_zero = BestDosePosterior::compute( &prior_theta, &prior_weights, None, + eq.clone(), + settings.clone(), + )?; + + let result_zero = posterior_zero.optimize( target_zero, None, - eq.clone(), DoseRange::new(10.0, 2000.0), 0.8, - settings.clone(), Target::AUCFromZero, )?; - - let result_zero = problem_zero.optimize()?; // Extract only the second dose (the optimized one at t=12) let dose_zero = result_zero.doses()[1]; @@ -975,20 +939,16 @@ fn test_auc_modes_comparison() -> Result<()> { .observation(24.0, 100.0, 0) // Target: AUC₁₂₋₂₄ = 100 .build(); - let problem_last = BestDoseProblem::new( - &prior_theta, - &prior_weights, - None, + let posterior_last = + BestDosePosterior::compute(&prior_theta, &prior_weights, None, eq, settings)?; + + let result_last = posterior_last.optimize( target_last, None, - eq, DoseRange::new(10.0, 2000.0), 0.8, - settings, Target::AUCFromLastDose, )?; - - let result_last = problem_last.optimize()?; // Extract only the second dose (the optimized one at t=12) let dose_last = result_last.doses()[1]; @@ -1081,20 +1041,15 @@ fn test_auc_from_last_dose_multiple_observations() -> Result<()> { }; let prior_weights = Weights::uniform(1); - let problem = BestDoseProblem::new( - &prior_theta, - &prior_weights, - None, + let posterior = BestDosePosterior::compute(&prior_theta, &prior_weights, None, eq, settings)?; + + let result = posterior.optimize( target, None, - eq, DoseRange::new(50.0, 500.0), 0.8, - settings, Target::AUCFromLastDose, )?; - - let result = problem.optimize()?; let doses: Vec = result.doses(); // Should optimize 2 doses @@ -1178,20 +1133,15 @@ fn test_auc_from_last_dose_no_prior_dose() -> Result<()> { }; let prior_weights = Weights::uniform(1); - let problem = BestDoseProblem::new( - &prior_theta, - &prior_weights, - None, + let posterior = BestDosePosterior::compute(&prior_theta, &prior_weights, None, eq, settings)?; + + let result = posterior.optimize( target, None, - eq, DoseRange::new(50.0, 500.0), 0.8, - settings, Target::AUCFromLastDose, )?; - - let result = problem.optimize()?; let doses: Vec = result.doses(); assert_eq!(doses.len(), 1); @@ -1274,20 +1224,16 @@ fn test_dose_range_bounds_respected() -> Result<()> { // Set a narrow dose range: 50-200 mg let dose_range = DoseRange::new(50.0, 200.0); - let problem = BestDoseProblem::new( + let posterior = BestDosePosterior::compute( &prior_theta, &prior_weights, None, - target.clone(), - None, eq.clone(), - dose_range, - 0.0, settings.clone(), - Target::Concentration, )?; - let result = problem.optimize()?; + let result = + posterior.optimize(target.clone(), None, dose_range, 0.0, Target::Concentration)?; let doses: Vec = result.doses(); println!("Optimal dose: {:.1} mg", doses[0]); @@ -1315,3 +1261,597 @@ fn test_dose_range_bounds_respected() -> Result<()> { Ok(()) } + +// ═════════════════════════════════════════════════════════════════════════════ +// Tests for time_offset behavior +// ═════════════════════════════════════════════════════════════════════════════ + +/// Helper to build a simple one-compartment model used by multiple tests +fn one_compartment_model() -> pharmsol::ODE { + equation::ODE::new( + |x, p, _t, dx, b, _rateiv, _cov| { + fetch_params!(p, ke, _v); + dx[0] = -ke * x[0] + b[0]; + }, + |_p, _, _| lag! {}, + |_p, _, _| fa! {}, + |_p, _t, _cov, _x| {}, + |x, p, _t, _cov, y| { + fetch_params!(p, _ke, v); + y[0] = x[0] / v; + }, + ) +} + +/// Helper to build minimal settings for tests (no posterior refinement) +fn minimal_settings() -> Settings { + let params = Parameters::new() + .add("ke", 0.001, 3.0) + .add("v", 25.0, 250.0); + let ems = AssayErrorModels::new() + .add( + 0, + AssayErrorModel::additive(ErrorPoly::new(0.0, 0.20, 0.0, 0.0), 0.0), + ) + .unwrap(); + let mut settings = Settings::builder() + .set_algorithm(Algorithm::NPAG) + .set_parameters(params) + .set_error_models(ems) + .build(); + settings.disable_output(); + settings.set_cycles(0); + settings +} + +/// Helper to build a simple prior (single support point) +fn simple_prior(settings: &Settings) -> (Theta, Weights) { + let mat = faer::Mat::from_fn(1, 2, |_r, c| match c { + 0 => 0.3, // ke + 1 => 50.0, // v + _ => 0.0, + }); + let theta = Theta::from_parts(mat, settings.parameters().clone()).unwrap(); + let weights = Weights::uniform(1); + (theta, weights) +} + +/// Test that gap=0 and gap=12 produce different results +/// +/// When time_offset is applied as a gap after the last past event, +/// different gaps change when the future dose is given relative to +/// the past, affecting the PK simulation outcome. +#[test] +fn test_time_offset_zero_vs_nonzero_differ() -> Result<()> { + let eq = one_compartment_model(); + let settings = minimal_settings(); + let (theta, weights) = simple_prior(&settings); + + // Past data: dose at t=0, observation at t=6 + let past = Subject::builder("patient") + .bolus(0.0, 500.0, 0) + .observation(6.0, 5.0, 0) + .build(); + + let posterior = + BestDosePosterior::compute(&theta, &weights, Some(past), eq.clone(), settings.clone())?; + + // Target: optimizable dose at t=0 (relative), target conc at t=1 (relative) + // Short observation window so residual from past dose matters + let target = Subject::builder("patient") + .bolus(0.0, 0.0, 0) + .observation(1.0, 5.0, 0) // target: 5 mg/L at 1h after the future dose + .build(); + + // gap=0: target dose at t=6 absolute (right after past), obs at t=7 + // Past dose (500mg at t=0): C(7) = 500/50 * e^(-0.3*7) ≈ 1.22 mg/L residual + let result_gap0 = posterior.optimize( + target.clone(), + Some(0.0), + DoseRange::new(10.0, 1000.0), + 0.5, + Target::Concentration, + )?; + + // gap=12: target dose at t=18 absolute, obs at t=19 + // Past dose (500mg at t=0): C(19) = 500/50 * e^(-0.3*19) ≈ 0.003 mg/L (negligible) + let result_gap12 = posterior.optimize( + target, + Some(12.0), + DoseRange::new(10.0, 1000.0), + 0.5, + Target::Concentration, + )?; + + let doses_gap0 = result_gap0.doses(); + let doses_gap12 = result_gap12.doses(); + + eprintln!("Gap=0 doses: {:?}", doses_gap0); + eprintln!("Gap=12 doses: {:?}", doses_gap12); + + // With gap=0, there's still significant residual from the past dose (~1.2 mg/L), + // so the optimizer needs less future dose. With gap=12, the past dose is negligible, + // so it needs more future dose. The optimizable doses should differ. + assert!( + (doses_gap0.last().unwrap() - doses_gap12.last().unwrap()).abs() > 1e-3, + "gap=0 and gap=12 must produce different optimizable doses, \ + but got {:.4} vs {:.4}", + doses_gap0.last().unwrap(), + doses_gap12.last().unwrap() + ); + + Ok(()) +} + +/// Test that the first target event lands at last_past_time + gap +/// and subsequent target times are shifted correctly. +#[test] +fn test_time_offset_event_placement() -> Result<()> { + let eq = one_compartment_model(); + let settings = minimal_settings(); + let (theta, weights) = simple_prior(&settings); + + // Past: dose at t=0, observation at t=6 (last event at t=6) + let past = Subject::builder("patient") + .bolus(0.0, 500.0, 0) + .observation(6.0, 5.0, 0) + .build(); + + let posterior = + BestDosePosterior::compute(&theta, &weights, Some(past), eq.clone(), settings.clone())?; + + // Target: dose at t=0, dose at t=12, obs at t=24 (all relative) + let target = Subject::builder("patient") + .bolus(0.0, 0.0, 0) + .bolus(12.0, 0.0, 0) + .observation(24.0, 5.0, 0) + .build(); + + // gap=0: future starts immediately after last past event (t=6) + // effective_offset = 6 + 0 = 6 + let gap = 0.0; + let result = posterior.optimize( + target, + Some(gap), + DoseRange::new(10.0, 500.0), + 0.5, + Target::Concentration, + )?; + + // After concatenation we should have: + // past dose at t=0 (fixed 500mg) + // target dose at t=0+6=6 (optimizable) + // target dose at t=12+6=18 (optimizable) + // target obs at t=24+6=30 + + let optimal_subject = result.optimal_subject(); + let mut dose_times = Vec::new(); + let mut obs_times = Vec::new(); + + for occ in optimal_subject.occasions() { + for event in occ.events() { + match event { + Event::Bolus(b) => dose_times.push(b.time()), + Event::Infusion(i) => dose_times.push(i.time()), + Event::Observation(o) => obs_times.push(o.time()), + } + } + } + + eprintln!("Dose times: {:?}", dose_times); + eprintln!("Obs times: {:?}", obs_times); + + // Past dose at t=0 + assert!( + (dose_times[0] - 0.0).abs() < 1e-10, + "First dose (past) should be at t=0, got {}", + dose_times[0] + ); + // First target dose at t = 0 + 6 = 6 + assert!( + (dose_times[1] - 6.0).abs() < 1e-10, + "Second dose should be at t=0+effective_offset=6, got {}", + dose_times[1] + ); + // Second target dose at t = 12 + 6 = 18 + assert!( + (dose_times[2] - 18.0).abs() < 1e-10, + "Third dose should be at t=12+effective_offset=18, got {}", + dose_times[2] + ); + // Observation at t = 24 + 6 = 30 + assert!( + (obs_times[0] - 30.0).abs() < 1e-10, + "Observation should be at t=24+effective_offset=30, got {}", + obs_times[0] + ); + + // Past dose should remain fixed at 500 + let doses = result.doses(); + assert!( + (doses[0] - 500.0).abs() < 1e-6, + "Past dose should be fixed at 500, got {}", + doses[0] + ); + + Ok(()) +} + +/// Test that time_offset=None leaves target events unchanged +#[test] +fn test_time_offset_none_no_shift() -> Result<()> { + let eq = one_compartment_model(); + let settings = minimal_settings(); + let (theta, weights) = simple_prior(&settings); + + let posterior = + BestDosePosterior::compute(&theta, &weights, None, eq.clone(), settings.clone())?; + + let target = Subject::builder("patient") + .bolus(0.0, 0.0, 0) + .bolus(12.0, 0.0, 0) + .observation(24.0, 5.0, 0) + .build(); + + let result = posterior.optimize( + target, + None, // No offset + DoseRange::new(10.0, 500.0), + 0.5, + Target::Concentration, + )?; + + let optimal_subject = result.optimal_subject(); + let mut dose_times = Vec::new(); + let mut obs_times = Vec::new(); + + for occ in optimal_subject.occasions() { + for event in occ.events() { + match event { + Event::Bolus(b) => dose_times.push(b.time()), + Event::Infusion(i) => dose_times.push(i.time()), + Event::Observation(o) => obs_times.push(o.time()), + } + } + } + + // Without offset, times should be exactly as specified in target + assert!((dose_times[0] - 0.0).abs() < 1e-10); + assert!((dose_times[1] - 12.0).abs() < 1e-10); + assert!((obs_times[0] - 24.0).abs() < 1e-10); + + Ok(()) +} + +// ═════════════════════════════════════════════════════════════════════════════ +// Tests for multi-target / multi-dose optimization +// ═════════════════════════════════════════════════════════════════════════════ + +/// Test that multiple optimizable doses all get meaningful values +#[test] +fn test_multi_dose_all_optimized() -> Result<()> { + let eq = one_compartment_model(); + let settings = minimal_settings(); + let (theta, weights) = simple_prior(&settings); + + let posterior = + BestDosePosterior::compute(&theta, &weights, None, eq.clone(), settings.clone())?; + + // Two optimizable doses, two target concentrations + let target = Subject::builder("patient") + .bolus(0.0, 0.0, 0) + .bolus(12.0, 0.0, 0) + .observation(6.0, 5.0, 0) // Target 5 mg/L at t=6 + .observation(18.0, 5.0, 0) // Target 5 mg/L at t=18 + .build(); + + let result = posterior.optimize( + target, + None, + DoseRange::new(10.0, 500.0), + 0.5, + Target::Concentration, + )?; + + let doses = result.doses(); + eprintln!("Multi-dose optimization: {:?}", doses); + + assert_eq!(doses.len(), 2, "Should optimize 2 doses"); + + // Both doses should be meaningful (not collapsed to minimum) + assert!( + doses[0] > 10.0 + 1.0, + "Dose 1 should be above minimum bound, got {}", + doses[0] + ); + assert!( + doses[1] > 10.0 + 1.0, + "Dose 2 should be above minimum bound, got {}", + doses[1] + ); + + Ok(()) +} + +/// Test that changing target for dose 2 changes dose 2's result +#[test] +fn test_multi_target_second_dose_responds_to_target_change() -> Result<()> { + let eq = one_compartment_model(); + let settings = minimal_settings(); + let (theta, weights) = simple_prior(&settings); + + let posterior = + BestDosePosterior::compute(&theta, &weights, None, eq.clone(), settings.clone())?; + + // Scenario A: second target is LOW (2 mg/L) + let target_low = Subject::builder("patient") + .bolus(0.0, 0.0, 0) + .bolus(12.0, 0.0, 0) + .observation(6.0, 5.0, 0) + .observation(18.0, 2.0, 0) // Low second target + .build(); + + // Scenario B: second target is HIGH (15 mg/L) + let target_high = Subject::builder("patient") + .bolus(0.0, 0.0, 0) + .bolus(12.0, 0.0, 0) + .observation(6.0, 5.0, 0) + .observation(18.0, 15.0, 0) // High second target + .build(); + + let result_low = posterior.optimize( + target_low, + None, + DoseRange::new(10.0, 1000.0), + 0.5, + Target::Concentration, + )?; + + let result_high = posterior.optimize( + target_high, + None, + DoseRange::new(10.0, 1000.0), + 0.5, + Target::Concentration, + )?; + + let doses_low = result_low.doses(); + let doses_high = result_high.doses(); + + eprintln!("Low second target: doses = {:?}", doses_low); + eprintln!("High second target: doses = {:?}", doses_high); + + // The second dose should be higher when the second target is higher + assert!( + doses_high[1] > doses_low[1], + "Higher second target ({}) should produce higher second dose, \ + but got low={:.2} vs high={:.2}", + 15.0, + doses_low[1], + doses_high[1] + ); + + Ok(()) +} + +// ═════════════════════════════════════════════════════════════════════════════ +// Tests for BestDosePosterior and BestDoseResult API surface +// ═════════════════════════════════════════════════════════════════════════════ + +/// Test BestDosePosterior accessor methods +#[test] +fn test_posterior_accessors() -> Result<()> { + let eq = one_compartment_model(); + let settings = minimal_settings(); + let (theta, weights) = simple_prior(&settings); + + let posterior = + BestDosePosterior::compute(&theta, &weights, None, eq.clone(), settings.clone())?; + + // n_support_points should match the prior (no filtering with 0 cycles and no data) + assert!( + posterior.n_support_points() > 0, + "Posterior should have at least 1 support point" + ); + + // theta() should return a valid Theta with the correct number of rows + assert_eq!( + posterior.theta().matrix().nrows(), + posterior.n_support_points() + ); + + // posterior_weights() should sum to ~1 + let weight_sum: f64 = posterior.posterior_weights().iter().sum(); + assert!( + (weight_sum - 1.0).abs() < 1e-6, + "Posterior weights should sum to 1.0, got {}", + weight_sum + ); + + // population_weights() should also sum to ~1 + let pop_weight_sum: f64 = posterior.population_weights().iter().sum(); + assert!( + (pop_weight_sum - 1.0).abs() < 1e-6, + "Population weights should sum to 1.0, got {}", + pop_weight_sum + ); + + Ok(()) +} + +/// Test BestDoseResult accessor methods +#[test] +fn test_result_accessors() -> Result<()> { + let eq = one_compartment_model(); + let settings = minimal_settings(); + let (theta, weights) = simple_prior(&settings); + + let posterior = + BestDosePosterior::compute(&theta, &weights, None, eq.clone(), settings.clone())?; + + let target = Subject::builder("patient") + .bolus(0.0, 0.0, 0) + .observation(6.0, 5.0, 0) + .build(); + + let result = posterior.optimize( + target, + None, + DoseRange::new(10.0, 500.0), + 0.5, + Target::Concentration, + )?; + + // doses() should return 1 dose + assert_eq!(result.doses().len(), 1); + assert!(result.doses()[0].is_finite()); + + // objf() should be finite and non-negative + assert!(result.objf().is_finite()); + assert!(result.objf() >= 0.0, "Cost should be non-negative"); + + // status() should be Converged (1000 iterations is usually enough for 1D) + assert_eq!( + *result.status(), + pmcore::bestdose::BestDoseStatus::Converged + ); + + // predictions() should have predictions + assert!( + !result.predictions().predictions().is_empty(), + "Predictions should not be empty" + ); + + // optimization_method() should be Posterior or Uniform + let method = result.optimization_method(); + assert!( + method == pmcore::bestdose::OptimalMethod::Posterior + || method == pmcore::bestdose::OptimalMethod::Uniform + ); + + // auc_predictions() should be None for concentration targets + assert!( + result.auc_predictions().is_none(), + "AUC predictions should be None for concentration targets" + ); + + // optimal_subject() should have the optimized dose + let optimal_subj = result.optimal_subject(); + let mut found_dose = false; + for occ in optimal_subj.occasions() { + for event in occ.events() { + if let Event::Bolus(b) = event { + assert!(b.amount() > 0.0, "Optimized dose should be > 0"); + found_dose = true; + } + } + } + assert!( + found_dose, + "Should find at least one dose in optimal subject" + ); + + Ok(()) +} + +/// Test that negative time_offset is rejected +#[test] +fn test_negative_time_offset_rejected() -> Result<()> { + let eq = one_compartment_model(); + let settings = minimal_settings(); + let (theta, weights) = simple_prior(&settings); + + let posterior = + BestDosePosterior::compute(&theta, &weights, None, eq.clone(), settings.clone())?; + + let target = Subject::builder("patient") + .bolus(0.0, 0.0, 0) + .observation(6.0, 5.0, 0) + .build(); + + let result = posterior.optimize( + target, + Some(-1.0), // Negative offset should be rejected + DoseRange::new(10.0, 500.0), + 0.5, + Target::Concentration, + ); + + assert!(result.is_err(), "Negative time_offset should be rejected"); + assert!( + result.unwrap_err().to_string().contains("negative"), + "Error message should mention negative time_offset" + ); + + Ok(()) +} + +/// Test that posterior can be reused for multiple optimizations +/// This is the key new feature of the two-stage API +#[test] +fn test_posterior_reuse() -> Result<()> { + let eq = one_compartment_model(); + let settings = minimal_settings(); + let (theta, weights) = simple_prior(&settings); + + // Compute posterior once + let posterior = + BestDosePosterior::compute(&theta, &weights, None, eq.clone(), settings.clone())?; + + // Optimize with different dose ranges + let target = Subject::builder("patient") + .bolus(0.0, 0.0, 0) + .observation(6.0, 5.0, 0) + .build(); + + let result_narrow = posterior.optimize( + target.clone(), + None, + DoseRange::new(10.0, 100.0), + 0.5, + Target::Concentration, + )?; + + let result_wide = posterior.optimize( + target.clone(), + None, + DoseRange::new(10.0, 1000.0), + 0.5, + Target::Concentration, + )?; + + // Both should succeed + assert!(result_narrow.doses()[0].is_finite()); + assert!(result_wide.doses()[0].is_finite()); + + // Wide range should allow a potentially better (lower cost) result + assert!( + result_wide.objf() <= result_narrow.objf() + 1e-6, + "Wider dose range should give equal or better cost: wide={:.6} vs narrow={:.6}", + result_wide.objf(), + result_narrow.objf() + ); + + // Optimize with different bias weights + let result_personal = posterior.optimize( + target.clone(), + None, + DoseRange::new(10.0, 500.0), + 0.0, // Full personalization + Target::Concentration, + )?; + + let result_population = posterior.optimize( + target, + None, + DoseRange::new(10.0, 500.0), + 1.0, // Full population weighting + Target::Concentration, + )?; + + // Both should succeed + assert!(result_personal.doses()[0].is_finite()); + assert!(result_population.doses()[0].is_finite()); + + Ok(()) +}