diff --git a/examples/bestdose.rs b/examples/bestdose.rs index 04392f02c..35e9deb06 100644 --- a/examples/bestdose.rs +++ b/examples/bestdose.rs @@ -114,10 +114,25 @@ fn main() -> Result<()> { // Print results for (bias_weight, optimal) in &results { + let opt_doses = optimal + .optimal_subject + .iter() + .flat_map(|occ| { + occ.events() + .iter() + .filter_map(|event| match event { + Event::Bolus(bolus) => Some(bolus.amount()), + Event::Infusion(infusion) => Some(infusion.amount()), + _ => None, + }) + .collect::>() + }) + .collect::>(); + println!( "Bias weight: {:.2}\t\t Optimal dose: {:?}\t\tCost: {:.6}\t\tln Cost: {:.4}\t\tMethod: {}", bias_weight, - optimal.dose, + opt_doses, optimal.objf, optimal.objf.ln(), optimal.optimization_method diff --git a/examples/bestdose_auc.rs b/examples/bestdose_auc.rs index e7f2b15b3..c92ba55bc 100644 --- a/examples/bestdose_auc.rs +++ b/examples/bestdose_auc.rs @@ -83,8 +83,23 @@ fn main() -> Result<()> { println!("Optimizing dose...\n"); let optimal = problem.optimize()?; + let opt_doses = optimal + .optimal_subject + .iter() + .flat_map(|occ| { + occ.events() + .iter() + .filter_map(|event| match event { + Event::Bolus(bolus) => Some(bolus.amount()), + Event::Infusion(infusion) => Some(infusion.amount()), + _ => None, + }) + .collect::>() + }) + .collect::>(); + println!("=== RESULTS ==="); - println!("Optimal dose: {:.1} mg", optimal.dose[0]); + println!("Optimal dose: {:.1} mg", opt_doses[0]); println!("Cost function: {:.6}", optimal.objf); if let Some(auc_preds) = &optimal.auc_predictions { diff --git a/src/bestdose/cost.rs b/src/bestdose/cost.rs index cbaa12950..7915a7015 100644 --- a/src/bestdose/cost.rs +++ b/src/bestdose/cost.rs @@ -76,7 +76,7 @@ use pharmsol::Equation; /// /// # Dose Masking /// -/// When `problem.current_time` is set (past/future separation), only doses where +/// When `problem.time_offset` is set (past/future separation), only doses where /// `dose_optimization_mask[i] == true` are updated with values from `candidate_doses`. /// Past doses (mask == false) remain at their historical values. /// @@ -239,7 +239,7 @@ pub fn calculate_cost(problem: &BestDoseProblem, candidate_doses: &[f64]) -> Res .matrix() .row_iter() .zip(problem.posterior.iter()) // Posterior from NPAGFULL11 (patient-specific) - .zip(problem.prior_weights.iter()) + .zip(problem.population_weights.iter()) // Prior (population) { let spp = row.iter().copied().collect::>(); diff --git a/src/bestdose/mod.rs b/src/bestdose/mod.rs index 588d34310..c568adbf6 100644 --- a/src/bestdose/mod.rs +++ b/src/bestdose/mod.rs @@ -11,8 +11,8 @@ //! ```rust,no_run,ignore //! use pmcore::bestdose::{BestDoseProblem, Target, DoseRange}; //! -//! # fn example(prior_theta: pmcore::structs::theta::Theta, -//! # prior_weights: pmcore::structs::weights::Weights, +//! # 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, @@ -21,8 +21,8 @@ //! # -> anyhow::Result<()> { //! // Create optimization problem //! let problem = BestDoseProblem::new( -//! &prior_theta, // Population support points from NPAG -//! &prior_weights, // Population probabilities +//! &population_theta, // Population support points from NPAG +//! &population_weights, // Population probabilities //! Some(past_data), // Patient history (None = use prior) //! target, // Future template with targets //! eq, // PK/PD model @@ -147,8 +147,8 @@ //! use pmcore::bestdose::{BestDoseProblem, Target, DoseRange}; //! use pharmsol::prelude::Subject; //! -//! # fn example(prior_theta: pmcore::structs::theta::Theta, -//! # prior_weights: pmcore::structs::weights::Weights, +//! # 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, @@ -161,7 +161,7 @@ //! .build(); //! //! let problem = BestDoseProblem::new( -//! &prior_theta, &prior_weights, Some(past), target, eq, error_models, +//! &population_theta, &population_weights, Some(past), target, eq, error_models, //! DoseRange::new(10.0, 500.0), // 10-500 mg allowed //! 0.3, // Slight population emphasis //! settings, 500, Target::Concentration, @@ -179,8 +179,8 @@ //! use pmcore::bestdose::{BestDoseProblem, Target, DoseRange}; //! use pharmsol::prelude::Subject; //! -//! # fn example(prior_theta: pmcore::structs::theta::Theta, -//! # prior_weights: pmcore::structs::weights::Weights, +//! # 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, @@ -194,7 +194,7 @@ //! .build(); //! //! let problem = BestDoseProblem::new( -//! &prior_theta, &prior_weights, Some(past), target, None, eq, error_models, +//! &population_theta, &population_weights, Some(past), target, eq, error_models, //! DoseRange::new(50.0, 300.0), //! 0.0, // Full personalization //! settings, Target::AUCFromZero, // Cumulative AUC target! @@ -214,8 +214,8 @@ //! //! ```rust,no_run,ignore //! # use pmcore::bestdose::{BestDoseProblem, Target, DoseRange}; -//! # fn example(prior_theta: pmcore::structs::theta::Theta, -//! # prior_weights: pmcore::structs::weights::Weights, +//! # 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, @@ -223,7 +223,7 @@ //! # -> anyhow::Result<()> { //! // No patient history - use population prior directly //! let problem = BestDoseProblem::new( -//! &prior_theta, &prior_weights, +//! &population_theta, &population_weights, //! None, // No past data //! target, eq, error_models, //! DoseRange::new(0.0, 1000.0), @@ -264,8 +264,8 @@ //! For faster optimization: //! ```rust,no_run,ignore //! # use pmcore::bestdose::{BestDoseProblem, Target, DoseRange}; -//! # fn example(prior_theta: pmcore::structs::theta::Theta, -//! # prior_weights: pmcore::structs::weights::Weights, +//! # 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, @@ -273,7 +273,7 @@ //! # -> anyhow::Result<()> { //! // Reduce refinement cycles //! let problem = BestDoseProblem::new( -//! &prior_theta, &prior_weights, None, target, eq, error_models, +//! &population_theta, &population_weights, None, target, eq, error_models, //! DoseRange::new(0.0, 1000.0), 0.5, //! settings.clone(), //! 100, // Faster: 100 instead of 500 @@ -306,20 +306,20 @@ pub use types::{BestDoseProblem, BestDoseResult, DoseRange, Target}; /// /// 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 `current_time` +/// 2. Offsets all future subject event times by `time_offset` /// 3. Combines into single continuous subject /// /// # Arguments /// /// * `past` - Subject with past history (only doses will be used) /// * `future` - Subject template for future (all events: doses + observations) -/// * `current_time` - Time offset to apply to all future events +/// * `time_offset` - Time offset to apply to all future events /// /// # Returns /// /// Combined subject with: -/// - Past doses at original times [0, current_time) -/// - Future doses + observations at offset times [current_time, ∞) +/// - Past doses at original times [0, time_offset) +/// - Future doses + observations at offset times [time_offset, ∞) /// /// # Example /// @@ -336,7 +336,7 @@ pub use types::{BestDoseProblem, BestDoseResult, DoseRange, Target}; /// .observation(24.0, 10.0, 0) // Target at t=30 absolute /// .build(); /// -/// // Concatenate with current_time = 6.0 +/// // Concatenate with time_offset = 6.0 /// let combined = concatenate_past_and_future(&past, &future, 6.0); /// // Result: dose at t=0 (fixed, 500mg), dose at t=6 (optimizable, 100mg initial), /// // observation target at t=30 (10 mg/L) @@ -344,7 +344,7 @@ pub use types::{BestDoseProblem, BestDoseResult, DoseRange, Target}; fn concatenate_past_and_future( past: &pharmsol::prelude::Subject, future: &pharmsol::prelude::Subject, - current_time: f64, + time_offset: f64, ) -> pharmsol::prelude::Subject { use pharmsol::prelude::*; @@ -374,11 +374,11 @@ fn concatenate_past_and_future( match event { Event::Bolus(bolus) => { builder = - builder.bolus(bolus.time() + current_time, bolus.amount(), bolus.input()); + builder.bolus(bolus.time() + time_offset, bolus.amount(), bolus.input()); } Event::Infusion(inf) => { builder = builder.infusion( - inf.time() + current_time, + inf.time() + time_offset, inf.amount(), inf.input(), inf.duration(), @@ -387,7 +387,7 @@ fn concatenate_past_and_future( Event::Observation(obs) => { builder = match obs.value() { Some(val) => { - builder.observation(obs.time() + current_time, val, obs.outeq()) + builder.observation(obs.time() + time_offset, val, obs.outeq()) } None => builder, }; @@ -466,8 +466,8 @@ use crate::structs::weights::Weights; // Helper Functions for STAGE 1: Posterior Density Calculation // ═════════════════════════════════════════════════════════════════════════════ -/// Validate current_time parameter for past/future separation mode -fn validate_current_time(current_time: f64, past_data: &Option) -> Result<()> { +/// 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 .occasions() @@ -480,11 +480,11 @@ fn validate_current_time(current_time: f64, past_data: &Option) -> Resu }) .fold(0.0_f64, |max, time| max.max(time)); - if current_time < max_past_time { + if time_offset < max_past_time { return Err(anyhow::anyhow!( - "Invalid current_time: {} is before the last past_data event at time {}. \ - current_time must be >= the maximum time in past_data to avoid time travel!", - current_time, + "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 )); } @@ -517,10 +517,10 @@ fn validate_current_time(current_time: f64, past_data: &Option) -> Resu /// /// # Returns /// -/// Tuple: (posterior_theta, posterior_weights, filtered_prior_weights, past_subject) +/// Tuple: (posterior_theta, posterior_weights, filtered_population_weights, past_subject) fn calculate_posterior_density( - prior_theta: &Theta, - prior_weights: &Weights, + population_theta: &Theta, + population_weights: &Weights, past_data: Option<&Subject>, eq: &ODE, error_models: &ErrorModels, @@ -530,9 +530,9 @@ fn calculate_posterior_density( None => { tracing::info!(" No past data → using prior directly"); Ok(( - prior_theta.clone(), - prior_weights.clone(), - prior_weights.clone(), + population_theta.clone(), + population_weights.clone(), + population_weights.clone(), Subject::builder("Empty").build(), )) } @@ -548,9 +548,9 @@ fn calculate_posterior_density( if !has_observations { tracing::info!(" Past data has no observations → using prior directly"); Ok(( - prior_theta.clone(), - prior_weights.clone(), - prior_weights.clone(), + population_theta.clone(), + population_weights.clone(), + population_weights.clone(), past_subject.clone(), )) } else { @@ -561,10 +561,10 @@ fn calculate_posterior_density( let past_data_obj = Data::new(vec![past_subject.clone()]); - let (posterior_theta, posterior_weights, filtered_prior_weights) = + let (posterior_theta, posterior_weights, filtered_population_weights) = posterior::calculate_two_step_posterior( - prior_theta, - prior_weights, + population_theta, + population_weights, &past_data_obj, eq, error_models, @@ -574,7 +574,7 @@ fn calculate_posterior_density( Ok(( posterior_theta, posterior_weights, - filtered_prior_weights, + filtered_population_weights, past_subject.clone(), )) } @@ -590,9 +590,9 @@ fn calculate_posterior_density( fn prepare_target_subject( past_subject: Subject, target: Subject, - current_time: Option, + time_offset: Option, ) -> Result<(Subject, Subject)> { - match current_time { + match time_offset { None => { tracing::info!(" Mode: Standard (single subject)"); Ok((target, past_subject)) @@ -647,11 +647,11 @@ impl BestDoseProblem { /// /// # Parameters /// - /// * `prior_theta` - Population support points from NPAG - /// * `prior_weights` - Population probabilities + /// * `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 - /// * `current_time` - Optional time offset for concatenation (None = standard mode, Some(t) = Fortran mode) + /// * `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 @@ -665,11 +665,11 @@ impl BestDoseProblem { /// BestDoseProblem ready for `optimize()` #[allow(clippy::too_many_arguments)] pub fn new( - prior_theta: &Theta, - prior_weights: &Weights, + population_theta: &Theta, + population_weights: &Weights, past_data: Option, target: Subject, - current_time: Option, + time_offset: Option, eq: ODE, error_models: ErrorModels, doserange: DoseRange, @@ -683,17 +683,17 @@ impl BestDoseProblem { tracing::info!("╚══════════════════════════════════════════════════════════╝"); // Validate input if using past/future separation mode - if let Some(t) = current_time { - validate_current_time(t, &past_data)?; + if let Some(t) = time_offset { + validate_time_offset(t, &past_data)?; } // ═════════════════════════════════════════════════════════════ // STAGE 1: Calculate Posterior Density // ═════════════════════════════════════════════════════════════ - let (posterior_theta, posterior_weights, filtered_prior_weights, past_subject) = + let (posterior_theta, posterior_weights, filtered_population_weights, past_subject) = calculate_posterior_density( - prior_theta, - prior_weights, + population_theta, + population_weights, past_data.as_ref(), &eq, &error_models, @@ -702,7 +702,7 @@ impl BestDoseProblem { // Handle past/future concatenation if needed let (final_target, final_past_data) = - prepare_target_subject(past_subject, target, current_time)?; + prepare_target_subject(past_subject, target, time_offset)?; tracing::info!("╔══════════════════════════════════════════════════════════╗"); tracing::info!("║ Stage 1 Complete - Ready for Optimization ║"); @@ -715,8 +715,8 @@ impl BestDoseProblem { past_data: final_past_data, target: final_target, target_type, - prior_theta: prior_theta.clone(), - prior_weights: filtered_prior_weights, + population_theta: population_theta.clone(), + population_weights: filtered_population_weights, theta: posterior_theta, posterior: posterior_weights, eq, @@ -724,7 +724,7 @@ impl BestDoseProblem { settings, doserange, bias_weight, - current_time, + time_offset, }) } diff --git a/src/bestdose/optimization.rs b/src/bestdose/optimization.rs index 2732177df..1e1e02d1d 100644 --- a/src/bestdose/optimization.rs +++ b/src/bestdose/optimization.rs @@ -92,7 +92,7 @@ impl CostFunction for BestDoseProblem { /// /// This is a helper for the dual optimization approach. /// -/// When `problem.current_time` is set (past/future separation mode): +/// When `problem.time_offset` is set (past/future separation mode): /// - Only optimizes doses where `dose_optimization_mask[i] == true` /// - Creates a reduced-dimension simplex for future doses only /// - Maps optimized doses back to full vector (past doses unchanged) @@ -190,14 +190,14 @@ fn run_single_optimization( /// │ │ /// │ OPTIMIZATION 1: Posterior Weights │ /// │ Use NPAGFULL11 posterior probabilities │ -/// │ → (doses₁, cost₁) │ +/// │ → (doses₁, cost₁) │ /// │ │ /// │ OPTIMIZATION 2: Uniform Weights │ /// │ Use equal weights (1/M) for all points │ -/// │ → (doses₂, cost₂) │ +/// │ → (doses₂, cost₂) │ /// │ │ -/// │ SELECTION: Choose min(cost₁, cost₂) │ -/// │ → (optimal_doses, optimal_cost, method) │ +/// │ SELECTION: Choose min(cost₁, cost₂) │ +/// │ → (optimal_doses, optimal_cost, method) │ /// └────────────┬────────────────────────────────────┘ /// ↓ /// ┌─────────────────────────────────────────────────┐ @@ -261,34 +261,34 @@ pub fn dual_optimization(problem: &BestDoseProblem) -> Result { method ); + // Generate target subject with optimal doses + let mut optimal_subject = problem.target.clone(); + let mut dose_number = 0; + + for occasion in optimal_subject.iter_mut() { + for event in occasion.iter_mut() { + match event { + Event::Bolus(bolus) => { + bolus.set_amount(final_doses[dose_number]); + dose_number += 1; + } + Event::Infusion(infusion) => { + infusion.set_amount(final_doses[dose_number]); + dose_number += 1; + } + Event::Observation(_) => {} + } + } + } + let (preds, auc_predictions) = calculate_final_predictions(problem, &final_doses, &final_weights)?; - // Extract only the optimized doses (exclude fixed past doses) - let original_doses: Vec = problem - .target - .iter() - .flat_map(|occ| { - occ.iter().filter_map(|event| match event { - Event::Bolus(bolus) => Some(bolus.amount()), - Event::Infusion(infusion) => Some(infusion.amount()), - Event::Observation(_) => None, - }) - }) - .collect(); - - let optimized_doses: Vec = final_doses - .iter() - .zip(original_doses.iter()) - .filter(|(_, &orig)| orig == 0.0) // Only doses that were placeholders - .map(|(&opt, _)| opt) - .collect(); - tracing::info!(" ✓ Predictions complete"); tracing::info!("─────────────────────────────────────────────────────────────"); Ok(BestDoseResult { - dose: optimized_doses, + optimal_subject, objf: final_cost, status: "Converged".to_string(), preds, diff --git a/src/bestdose/posterior.rs b/src/bestdose/posterior.rs index a89774c13..9109f65c1 100644 --- a/src/bestdose/posterior.rs +++ b/src/bestdose/posterior.rs @@ -88,10 +88,10 @@ const KEEP_UNREFINED_POINTS: bool = true; /// /// Note: This uses only lambda filtering, NO QR decomposition or second burke call. /// -/// Returns: (filtered_theta, filtered_posterior_weights, filtered_prior_weights) +/// Returns: (filtered_theta, filtered_posterior_weights, filtered_population_weights) pub fn npagfull11_filter( - prior_theta: &Theta, - prior_weights: &Weights, + population_theta: &Theta, + population_weights: &Weights, past_data: &Data, eq: &ODE, error_models: &ErrorModels, @@ -99,7 +99,7 @@ pub fn npagfull11_filter( tracing::info!("Stage 1.1: NPAGFULL11 Bayesian filtering"); // Calculate psi matrix P(data|theta_i) for all support points - let psi = calculate_psi(eq, past_data, prior_theta, error_models, false, true)?; + let psi = calculate_psi(eq, past_data, population_theta, error_models, false, true)?; // First burke call to get initial posterior probabilities let (initial_weights, _) = burke(&psi)?; @@ -120,7 +120,7 @@ pub fn npagfull11_filter( .collect(); // Filter theta to keep only points above threshold - let mut filtered_theta = prior_theta.clone(); + let mut filtered_theta = population_theta.clone(); filtered_theta.filter_indices(&keep_lambda); // Filter and renormalize posterior weights @@ -130,10 +130,11 @@ pub fn npagfull11_filter( Weights::from_vec(filtered_weights.iter().map(|w| w / sum).collect()); // Also filter the prior weights to match the filtered theta - let filtered_prior_weights: Vec = keep_lambda.iter().map(|&i| prior_weights[i]).collect(); - let prior_sum: f64 = filtered_prior_weights.iter().sum(); - let final_prior_weights = Weights::from_vec( - filtered_prior_weights + let filtered_population_weights: Vec = + keep_lambda.iter().map(|&i| population_weights[i]).collect(); + let prior_sum: f64 = filtered_population_weights.iter().sum(); + let final_population_weights = Weights::from_vec( + filtered_population_weights .iter() .map(|w| w / prior_sum) .collect(), @@ -141,12 +142,16 @@ pub fn npagfull11_filter( tracing::info!( " {} → {} support points (lambda filter, threshold={:.0e})", - prior_theta.matrix().nrows(), + population_theta.matrix().nrows(), filtered_theta.matrix().nrows(), threshold * max_weight ); - Ok((filtered_theta, final_posterior_weights, final_prior_weights)) + Ok(( + filtered_theta, + final_posterior_weights, + final_population_weights, + )) } /// Step 1.2: NPAGFULL - Refine each filtered point with full NPAG optimization @@ -300,11 +305,11 @@ pub fn npagfull_refinement( /// /// This is the complete Stage 1 of the BestDose algorithm. /// -/// Returns (posterior_theta, posterior_weights, filtered_prior_weights) suitable for dose optimization. -/// The filtered_prior_weights are the original prior weights filtered to match the posterior support points. +/// Returns (posterior_theta, posterior_weights, filtered_population_weights) suitable for dose optimization. +/// The filtered_population_weights are the original prior weights filtered to match the posterior support points. pub fn calculate_two_step_posterior( - prior_theta: &Theta, - prior_weights: &Weights, + population_theta: &Theta, + population_weights: &Weights, past_data: &Data, eq: &ODE, error_models: &ErrorModels, @@ -313,8 +318,14 @@ pub fn calculate_two_step_posterior( tracing::info!("=== STAGE 1: Posterior Density Calculation ==="); // Step 1.1: NPAGFULL11 filtering (returns filtered posterior AND filtered prior) - let (filtered_theta, filtered_posterior_weights, filtered_prior_weights) = - npagfull11_filter(prior_theta, prior_weights, past_data, eq, error_models)?; + let (filtered_theta, filtered_posterior_weights, filtered_population_weights) = + npagfull11_filter( + population_theta, + population_weights, + past_data, + eq, + error_models, + )?; // Step 1.2: NPAGFULL refinement let (refined_theta, refined_weights) = npagfull_refinement( @@ -330,5 +341,5 @@ pub fn calculate_two_step_posterior( refined_theta.matrix().nrows() ); - Ok((refined_theta, refined_weights, filtered_prior_weights)) + Ok((refined_theta, refined_weights, filtered_population_weights)) } diff --git a/src/bestdose/types.rs b/src/bestdose/types.rs index 45ccc1957..7702b2bd8 100644 --- a/src/bestdose/types.rs +++ b/src/bestdose/types.rs @@ -208,8 +208,8 @@ impl Default for DoseRange { /// - `target_type`: [`Target::Concentration`] or [`Target::AUC`] /// /// ## Population Prior -/// - `prior_theta`: Support points from NPAG population model -/// - `prior_weights`: Probability weights for each support point +/// - `population_theta`: Support points from NPAG population model +/// - `population_weights`: Probability weights for each support point /// /// ## Patient-Specific Posterior /// - `theta`: Refined posterior support points (from NPAGFULL11 + NPAGFULL) @@ -229,8 +229,8 @@ impl Default for DoseRange { /// ```rust,no_run,ignore /// use pmcore::bestdose::{BestDoseProblem, Target, DoseRange}; /// -/// # fn example(prior_theta: pmcore::structs::theta::Theta, -/// # prior_weights: pmcore::structs::weights::Weights, +/// # 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, @@ -238,8 +238,8 @@ impl Default for DoseRange { /// # settings: pmcore::routines::settings::Settings) /// # -> anyhow::Result<()> { /// let problem = BestDoseProblem::new( -/// &prior_theta, -/// &prior_weights, +/// &population_theta, +/// &population_weights, /// Some(past), // Patient history /// target, // Dosing template with targets /// eq, @@ -258,26 +258,45 @@ impl Default for DoseRange { #[derive(Debug, Clone)] pub struct BestDoseProblem { // Input data - pub past_data: Subject, - pub target: Subject, - pub target_type: Target, + /// Past patient data for posterior calculation + /// + /// These observations are used to refine the population prior into a + /// patient-specific posterior, and will be used to inform dose optimization. + pub(crate) past_data: Subject, + /// 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(crate) target: Subject, + /// Target type for optimization + /// + /// Specifies whether to optimize for concentrations or AUC values. + pub(crate) target_type: Target, // Population prior - pub prior_theta: Theta, - pub prior_weights: Weights, + /// The population prior support points ([Theta]), representing your previous knowledge of the population parameter distribution. + pub(crate) population_theta: Theta, + /// 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 theta: Theta, - pub posterior: Weights, + pub(crate) theta: Theta, + pub(crate) posterior: Weights, // Model and settings - pub eq: ODE, - pub error_models: ErrorModels, - pub settings: Settings, + pub(crate) eq: ODE, + pub(crate) error_models: ErrorModels, + pub(crate) settings: Settings, // Optimization parameters - pub doserange: DoseRange, - pub bias_weight: f64, // λ: 0=personalized, 1=population + pub(crate) doserange: DoseRange, + pub(crate) bias_weight: f64, // λ: 0=personalized, 1=population /// Time offset between past and future data (used for concatenation) /// When Some(t): future events were offset by this time to create continuous simulation @@ -285,7 +304,21 @@ pub struct BestDoseProblem { /// /// This is used to track the boundary between past and future for reporting/debugging. /// The actual optimization mask is derived from dose amounts (0 = optimize, >0 = fixed). - pub current_time: Option, + pub(crate) time_offset: Option, +} + +impl BestDoseProblem { + /// Validate input + pub(crate) fn validate(&self) -> anyhow::Result<()> { + if self.bias_weight <= 0.0 || self.bias_weight >= 1.0 { + return Err(anyhow::anyhow!( + "Bias weight must be between 0.0 and 1.0, got {}", + self.bias_weight + )); + } + + Ok(()) + } } /// Result from BestDose optimization @@ -359,13 +392,13 @@ pub struct BestDoseProblem { /// # Ok(()) /// # } /// ``` -#[derive(Debug)] +#[derive(Debug, Clone)] pub struct BestDoseResult { /// Optimal dose amount(s) /// /// Vector contains one element per dose in the target subject. /// Order matches the dose events in the target subject. - pub dose: Vec, + pub optimal_subject: Subject, /// Final cost function value /// diff --git a/tests/bestdose_tests.rs b/tests/bestdose_tests.rs index 8ed3ac284..ddf368bb3 100644 --- a/tests/bestdose_tests.rs +++ b/tests/bestdose_tests.rs @@ -105,16 +105,42 @@ fn test_infusion_mask_inclusion() -> Result<()> { // We should get back 1 optimized dose (the infusion placeholder) assert_eq!( - result.dose.len(), + result + .optimal_subject + .iter() + .flat_map(|occ| { + occ.events() + .iter() + .filter_map(|event| match event { + Event::Infusion(inf) => Some(inf.amount()), + _ => None, + }) + .collect::>() + }) + .count(), 1, "Should have 1 optimized dose (the infusion)" ); + let optinf = result + .optimal_subject + .iter() + .flat_map(|occ| { + occ.events() + .iter() + .filter_map(|event| match event { + Event::Infusion(inf) => Some(inf.amount()), + _ => None, + }) + .collect::>() + }) + .collect::>(); + // The optimized dose should be reasonable (not NaN, not infinite) assert!( - result.dose[0].is_finite(), + optinf[0].is_finite(), "Optimized dose should be finite, got {}", - result.dose[0] + optinf[0] ); Ok(()) @@ -196,7 +222,22 @@ fn test_fixed_infusion_preservation() -> Result<()> { let result = problem.optimize()?; // Should only optimize the future bolus, not the past infusion - assert_eq!(result.dose.len(), 1, "Should have 1 optimized dose"); + let doses = result + .optimal_subject + .iter() + .flat_map(|occ| { + occ.events() + .iter() + .filter_map(|event| match event { + Event::Infusion(inf) if inf.amount() != 200.0 => Some(inf.amount()), + Event::Bolus(bol) => Some(bol.amount()), + _ => None, + }) + .collect::>() + }) + .collect::>(); + eprintln!("Optimized doses: {:?}", doses); + assert_eq!(doses.len(), 1, "Should have 1 optimized dose"); Ok(()) } @@ -424,7 +465,21 @@ fn test_basic_auc_mode() -> Result<()> { ); let result = result?; - assert_eq!(result.dose.len(), 1); + let doses = result + .optimal_subject + .iter() + .flat_map(|occ| { + occ.events() + .iter() + .filter_map(|event| match event { + Event::Infusion(inf) => Some(inf.amount()), + Event::Bolus(bol) => Some(bol.amount()), + _ => None, + }) + .collect::>() + }) + .collect::>(); + assert_eq!(doses.len(), 1); assert!(result.auc_predictions.is_some()); @@ -519,11 +574,25 @@ fn test_infusion_auc_mode() -> Result<()> { ); let result = result?; - - eprintln!("Optimized dose: {}", result.dose[0]); + let doses = result + .optimal_subject + .iter() + .flat_map(|occ| { + occ.events() + .iter() + .filter_map(|event| match event { + Event::Infusion(inf) => Some(inf.amount()), + Event::Bolus(bol) => Some(bol.amount()), + _ => None, + }) + .collect::>() + }) + .collect::>(); + + eprintln!("Optimized dose: {:?}", doses); // Should have 1 optimized dose (the infusion) - assert_eq!(result.dose.len(), 1, "Should have 1 optimized dose"); + assert_eq!(doses.len(), 1, "Should have 1 optimized dose"); // Should have AUC predictions assert!( @@ -697,8 +766,24 @@ fn test_multi_outeq_auc_optimization() -> Result<()> { ); let best_dose_result = result?; - assert_eq!(best_dose_result.dose.len(), 1); - assert!(best_dose_result.dose[0] > 0.0); + + let doses = best_dose_result + .optimal_subject + .iter() + .flat_map(|occ| { + occ.events() + .iter() + .filter_map(|event| match event { + Event::Infusion(inf) => Some(inf.amount()), + Event::Bolus(bol) => Some(bol.amount()), + _ => None, + }) + .collect::>() + }) + .collect::>(); + + assert_eq!(doses.len(), 1); + assert!(doses[0] > 0.0); assert!(best_dose_result.objf.is_finite()); assert!(best_dose_result.auc_predictions.is_some());