From a6dc5fc1640659df7997c83e739b07554671ac47 Mon Sep 17 00:00:00 2001 From: Markus Date: Sat, 25 Oct 2025 15:09:28 +0200 Subject: [PATCH 1/5] Documentation --- src/bestdose/types.rs | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/src/bestdose/types.rs b/src/bestdose/types.rs index b25180432..85b46f4fd 100644 --- a/src/bestdose/types.rs +++ b/src/bestdose/types.rs @@ -214,12 +214,31 @@ impl Default for DoseRange { #[derive(Debug, Clone)] pub struct BestDoseProblem { // Input data + /// 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 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 target: Subject, + /// Target type for optimization + /// + /// Specifies whether to optimize for concentrations or AUC values. pub target_type: Target, // Population prior + /// The population prior support points ([Theta]), representing your previous knowledge of the population parameter distribution. pub prior_theta: Theta, + /// The population prior weights ([Weights]), representing the probability of each support point in the population. pub prior_weights: Weights, // Patient-specific posterior (from NPAGFULL11 + NPAGFULL) From 9b36e886f468908bed12a6388e3078b7f6462283 Mon Sep 17 00:00:00 2001 From: Markus Date: Sat, 25 Oct 2025 15:50:05 +0200 Subject: [PATCH 2/5] Make fields of BestDoseProblem private and rename field --- src/bestdose/cost.rs | 2 +- src/bestdose/mod.rs | 46 ++++++++++++++++++------------------ src/bestdose/optimization.rs | 2 +- src/bestdose/types.rs | 26 ++++++++++---------- 4 files changed, 38 insertions(+), 38 deletions(-) diff --git a/src/bestdose/cost.rs b/src/bestdose/cost.rs index c3439bdcc..acd93a1f9 100644 --- a/src/bestdose/cost.rs +++ b/src/bestdose/cost.rs @@ -70,7 +70,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. /// diff --git a/src/bestdose/mod.rs b/src/bestdose/mod.rs index c239ea9d5..943260d76 100644 --- a/src/bestdose/mod.rs +++ b/src/bestdose/mod.rs @@ -303,20 +303,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 /// @@ -333,7 +333,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) @@ -341,7 +341,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::*; @@ -371,11 +371,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(), @@ -384,7 +384,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, }; @@ -463,8 +463,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() @@ -477,11 +477,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 )); } @@ -589,9 +589,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)) @@ -650,7 +650,7 @@ impl BestDoseProblem { /// * `prior_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 @@ -668,7 +668,7 @@ impl BestDoseProblem { prior_weights: &Weights, past_data: Option, target: Subject, - current_time: Option, + time_offset: Option, eq: ODE, error_models: ErrorModels, doserange: DoseRange, @@ -683,8 +683,8 @@ 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)?; } // ═════════════════════════════════════════════════════════════ @@ -703,7 +703,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 ║"); @@ -725,7 +725,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..a1ce81434 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) diff --git a/src/bestdose/types.rs b/src/bestdose/types.rs index 85b46f4fd..4c564432a 100644 --- a/src/bestdose/types.rs +++ b/src/bestdose/types.rs @@ -218,7 +218,7 @@ pub struct BestDoseProblem { /// /// These observations are used to refine the population prior into a /// patient-specific posterior, and will be used to inform dose optimization. - pub past_data: Subject, + pub(crate) past_data: Subject, /// Target subject with dosing template and target observations /// /// This [Subject] defines the targets for optimization, including @@ -229,30 +229,30 @@ pub struct BestDoseProblem { /// 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 target: Subject, + pub(crate) target: Subject, /// Target type for optimization /// /// Specifies whether to optimize for concentrations or AUC values. - pub target_type: Target, + pub(crate) target_type: Target, // Population prior /// The population prior support points ([Theta]), representing your previous knowledge of the population parameter distribution. - pub prior_theta: Theta, + pub(crate) prior_theta: Theta, /// The population prior weights ([Weights]), representing the probability of each support point in the population. - pub prior_weights: Weights, + pub(crate) prior_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 @@ -260,7 +260,7 @@ 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, } /// Result from BestDose optimization From ea57ea6ece01be9cc1f689d752a80fd566d36edc Mon Sep 17 00:00:00 2001 From: Markus Date: Sat, 25 Oct 2025 15:56:13 +0200 Subject: [PATCH 3/5] Rename fields --- src/bestdose/cost.rs | 2 +- src/bestdose/mod.rs | 76 +++++++++++++++++++-------------------- src/bestdose/posterior.rs | 47 ++++++++++++++---------- src/bestdose/types.rs | 30 +++++++++++----- 4 files changed, 90 insertions(+), 65 deletions(-) diff --git a/src/bestdose/cost.rs b/src/bestdose/cost.rs index acd93a1f9..f027f754b 100644 --- a/src/bestdose/cost.rs +++ b/src/bestdose/cost.rs @@ -187,7 +187,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 943260d76..ef681d6c0 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 @@ -145,8 +145,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, @@ -159,7 +159,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, @@ -177,8 +177,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, @@ -192,7 +192,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(50.0, 300.0), //! 0.0, // Full personalization //! settings, 500, Target::AUC, // AUC target! @@ -212,8 +212,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, @@ -221,7 +221,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), @@ -261,8 +261,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, @@ -270,7 +270,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 @@ -514,10 +514,10 @@ fn validate_time_offset(time_offset: f64, past_data: &Option) -> Result /// /// # 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, @@ -528,9 +528,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(), )) } @@ -546,9 +546,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 { @@ -559,10 +559,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, @@ -573,7 +573,7 @@ fn calculate_posterior_density( Ok(( posterior_theta, posterior_weights, - filtered_prior_weights, + filtered_population_weights, past_subject.clone(), )) } @@ -646,8 +646,8 @@ 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 /// * `time_offset` - Optional time offset for concatenation (None = standard mode, Some(t) = Fortran mode) @@ -664,8 +664,8 @@ 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, time_offset: Option, @@ -690,10 +690,10 @@ impl BestDoseProblem { // ═════════════════════════════════════════════════════════════ // 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, @@ -716,8 +716,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, diff --git a/src/bestdose/posterior.rs b/src/bestdose/posterior.rs index 9594df285..22324d9e7 100644 --- a/src/bestdose/posterior.rs +++ b/src/bestdose/posterior.rs @@ -87,10 +87,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, @@ -98,7 +98,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)?; @@ -119,7 +119,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 @@ -129,10 +129,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(), @@ -140,12 +141,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 @@ -292,11 +297,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, @@ -306,8 +311,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( @@ -324,5 +335,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 4c564432a..12d18eef8 100644 --- a/src/bestdose/types.rs +++ b/src/bestdose/types.rs @@ -164,8 +164,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) @@ -185,8 +185,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, @@ -194,8 +194,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, @@ -237,9 +237,9 @@ pub struct BestDoseProblem { // Population prior /// The population prior support points ([Theta]), representing your previous knowledge of the population parameter distribution. - pub(crate) prior_theta: Theta, + pub(crate) population_theta: Theta, /// The population prior weights ([Weights]), representing the probability of each support point in the population. - pub(crate) prior_weights: Weights, + pub(crate) population_weights: Weights, // Patient-specific posterior (from NPAGFULL11 + NPAGFULL) pub(crate) theta: Theta, @@ -263,6 +263,20 @@ pub struct BestDoseProblem { 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 /// /// Contains the optimal doses and associated predictions from running From 431338f024a7d64320768f3a392d47675b5229a6 Mon Sep 17 00:00:00 2001 From: Markus Date: Sun, 26 Oct 2025 12:31:54 +0100 Subject: [PATCH 4/5] Change output of result --- examples/bestdose.rs | 17 +++++++++++- examples/bestdose_auc.rs | 17 +++++++++++- src/bestdose/optimization.rs | 50 ++++++++++++++++++------------------ src/bestdose/types.rs | 2 +- 4 files changed, 58 insertions(+), 28 deletions(-) diff --git a/examples/bestdose.rs b/examples/bestdose.rs index 569502499..bf9ffe9ac 100644 --- a/examples/bestdose.rs +++ b/examples/bestdose.rs @@ -115,10 +115,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 0ad68707f..eb6231246 100644 --- a/examples/bestdose_auc.rs +++ b/examples/bestdose_auc.rs @@ -84,8 +84,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/optimization.rs b/src/bestdose/optimization.rs index a1ce81434..1e1e02d1d 100644 --- a/src/bestdose/optimization.rs +++ b/src/bestdose/optimization.rs @@ -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/types.rs b/src/bestdose/types.rs index 12d18eef8..b1b6abdf0 100644 --- a/src/bestdose/types.rs +++ b/src/bestdose/types.rs @@ -354,7 +354,7 @@ pub struct BestDoseResult { /// /// 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 /// From a267facf5fd3ff5b8b9a1a70e6052384761f387a Mon Sep 17 00:00:00 2001 From: Markus Date: Tue, 4 Nov 2025 18:01:55 +0100 Subject: [PATCH 5/5] Update tests --- src/bestdose/types.rs | 2 +- tests/bestdose_tests.rs | 105 ++++++++++++++++++++++++++++++++++++---- 2 files changed, 96 insertions(+), 11 deletions(-) diff --git a/src/bestdose/types.rs b/src/bestdose/types.rs index b1b6abdf0..258539053 100644 --- a/src/bestdose/types.rs +++ b/src/bestdose/types.rs @@ -348,7 +348,7 @@ impl BestDoseProblem { /// # Ok(()) /// # } /// ``` -#[derive(Debug)] +#[derive(Debug, Clone)] pub struct BestDoseResult { /// Optimal dose amount(s) /// diff --git a/tests/bestdose_tests.rs b/tests/bestdose_tests.rs index f7ca33e89..8992ed01f 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());