diff --git a/examples/bestdose.rs b/examples/bestdose.rs index 34ff95364..7d8ebf96c 100644 --- a/examples/bestdose.rs +++ b/examples/bestdose.rs @@ -73,7 +73,7 @@ fn main() -> Result<()> { .build(); let (theta, prior) = parse_prior( - &"examples/bimodal_ke/output/theta.csv".to_string(), + &"examples/bimodal_ke/prior.csv".to_string(), &settings, ) .unwrap(); @@ -131,5 +131,38 @@ fn main() -> Result<()> { ); } + // Print the posterior support points with their filtered population and posterior weights. + let posterior_theta = problem.posterior_theta(); + let posterior_weights = problem.posterior_weights(); + let population_weights = problem.population_weights(); + let param_names = posterior_theta.param_names(); + + println!("\n=== Support Points Summary ==="); + println!("Number of support points: {}", posterior_theta.nspp()); + + print!("\n{:<8} {:<15} {:<15}", "Point", "Prior Weight", "Posterior Weight"); + for name in ¶m_names { + print!(" {:<15}", name); + } + println!(); + println!("{}", "-".repeat(40 + 16 * param_names.len())); + + for point_idx in 0..posterior_theta.nspp() { + let row = posterior_theta.matrix().row(point_idx); + + print!( + "{:<8} {:<15.6e} {:<15.6e}", + point_idx, + population_weights[point_idx], + posterior_weights[point_idx] + ); + + for value in row.iter() { + print!(" {:<15.6}", value); + } + + println!(); + } + Ok(()) } diff --git a/examples/bestdose_posterior.rs b/examples/bestdose_posterior.rs new file mode 100644 index 000000000..972e031a5 --- /dev/null +++ b/examples/bestdose_posterior.rs @@ -0,0 +1,161 @@ +use anyhow::Result; +use pmcore::bestdose; // bestdose new + // use pmcore::bestdose::bestdose_old as bestdose; // bestdose old + +use pmcore::prelude::*; +use pmcore::routines::initialization::parse_prior; + +fn main() -> Result<()> { + // Example model + let eq = ode! { + diffeq: |x, p, _t, dx, b, _rateiv, _cov| { + // fetch_cov!(cov, t, wt); + fetch_params!(p, ke, _v); + dx[0] = -ke * x[0] + b[0]; + }, + out: |x, p, _t, _cov, y| { + fetch_params!(p, _ke, v); + y[0] = x[0] / v; + }, + }; + + 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), + )?; + + // Make settings + let mut settings = Settings::builder() + .set_algorithm(Algorithm::NPAG) + .set_parameters(params) + .set_error_models(ems.clone()) + .build(); + + settings.disable_output(); + + // Generate a patient with known parameters + // Ke = 0.5, V = 50 + // C(t) = Dose * exp(-ke * t) / V + + fn conc(t: f64, dose: f64) -> f64 { + let ke = 0.3406021231412888; // Elimination rate constant + let v = 99.99475717544556; // Volume of distribution + (dose * (-ke * t).exp()) / v + } + + // Some observed data + let subject = Subject::builder("Nikola Tesla") + .bolus(0.0, 150.0, 0) + .observation(2.0, conc(2.0, 150.0), 0) + .observation(4.0, conc(4.0, 150.0), 0) + .observation(6.0, conc(6.0, 150.0), 0) + .bolus(12.0, 75.0, 0) + .observation(14.0, conc(2.0, 75.0) + conc(14.0, 150.0), 0) + .observation(16.0, conc(4.0, 75.0) + conc(16.0, 150.0), 0) + .observation(18.0, conc(6.0, 75.0) + conc(18.0, 150.0), 0) + .build(); + + let past_data = subject.clone(); + + let target_data = Subject::builder("Thomas Edison") + .bolus(0.0, 0.0, 0) + .observation(2.0, conc(2.0, 150.0), 0) + .observation(4.0, conc(4.0, 150.0), 0) + .observation(6.0, conc(6.0, 150.0), 0) + .bolus(12.0, 0.0, 0) + .observation(14.0, conc(2.0, 75.0) + conc(14.0, 150.0), 0) + .observation(16.0, conc(4.0, 75.0) + conc(16.0, 150.0), 0) + .observation(18.0, conc(6.0, 75.0) + conc(18.0, 150.0), 0) + .build(); + + let (theta, prior) = parse_prior( + &"examples/bimodal_ke/prior.csv".to_string(), + &settings, + ) + .unwrap(); + + let problem = bestdose::BestDoseProblem::new( + &theta, + &prior.unwrap(), + Some(past_data.clone()), + target_data.clone(), + None, + eq.clone(), + bestdose::DoseRange::new(0.0, 300.0), + 0.0, + settings.clone(), + bestdose::Target::Concentration, + )? + .with_optimization_strategy(bestdose::OptimizationStrategy::PosteriorOnly); + + println!("Optimizing dose with posterior-only strategy..."); + + let bias_weights = vec![0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]; + let mut results = Vec::new(); + + for bias_weight in &bias_weights { + println!("Running optimization with bias weight: {}", bias_weight); + let optimal = problem.clone().with_bias_weight(*bias_weight).optimize()?; + results.push((bias_weight, optimal)); + } + + for (bias_weight, optimal) in &results { + let opt_doses = optimal.doses(); + + println!( + "Bias weight: {:.2}\t\t Optimal dose: {:?}\t\tCost: {:.6}\t\tln Cost: {:.4}\t\tMethod: {}", + bias_weight, + opt_doses, + optimal.objf(), + optimal.objf().ln(), + optimal.optimization_method() + ); + } + + let optimal = &results.last().unwrap().1; + println!("\nConcentration-time predictions for optimal dose:"); + for pred in optimal.predictions().predictions().into_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() + ); + } + + let posterior_theta = problem.posterior_theta(); + let posterior_weights = problem.posterior_weights(); + let population_weights = problem.population_weights(); + let param_names = posterior_theta.param_names(); + + println!("\n=== Support Points Summary ==="); + println!("Number of support points: {}", posterior_theta.nspp()); + + print!("\n{:<8} {:<15} {:<15}", "Point", "Prior Weight", "Posterior Weight"); + for name in ¶m_names { + print!(" {:<15}", name); + } + println!(); + println!("{}", "-".repeat(40 + 16 * param_names.len())); + + for point_idx in 0..posterior_theta.nspp() { + let row = posterior_theta.matrix().row(point_idx); + + print!( + "{:<8} {:<15.6e} {:<15.6e}", + point_idx, + population_weights[point_idx], + posterior_weights[point_idx] + ); + + for value in row.iter() { + print!(" {:<15.6}", value); + } + + println!(); + } + + Ok(()) +} \ No newline at end of file diff --git a/src/bestdose/cost.rs b/src/bestdose/cost.rs index c3ba15aa0..9c78cea08 100644 --- a/src/bestdose/cost.rs +++ b/src/bestdose/cost.rs @@ -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() @@ -510,7 +510,7 @@ pub fn calculate_cost(problem: &BestDoseProblem, candidate_doses: &[f64]) -> Res let se = (obs_val - pj).powi(2); sumsq_i += se; // Calculate population mean using PRIOR probabilities - y_bar[j] += prior_prob * pj; + 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..9500450d4 100644 --- a/src/bestdose/mod.rs +++ b/src/bestdose/mod.rs @@ -302,7 +302,7 @@ pub mod predictions; mod types; // Re-export public API -pub use types::{BestDoseProblem, BestDoseResult, DoseRange, Target}; +pub use types::{BestDoseProblem, BestDoseResult, DoseRange, OptimizationStrategy, Target}; /// Helper function to concatenate past and future subjects (Option 3: Fortran MAKETMP approach) /// @@ -710,6 +710,7 @@ impl BestDoseProblem { tracing::info!(" Support points: {}", posterior_theta.matrix().nrows()); tracing::info!(" Target type: {:?}", target_type); tracing::info!(" Bias weight (λ): {}", bias_weight); + tracing::info!(" Optimization strategy: {}", OptimizationStrategy::Dual); Ok(BestDoseProblem { target: final_target, @@ -721,6 +722,7 @@ impl BestDoseProblem { settings, doserange, bias_weight, + optimization_strategy: OptimizationStrategy::Dual, }) } @@ -767,13 +769,24 @@ impl BestDoseProblem { /// - `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!("╚══════════════════════════════════════════════════════════╝"); + match self.optimization_strategy { + OptimizationStrategy::Dual => { + tracing::info!("╔══════════════════════════════════════════════════════════╗"); + tracing::info!("║ BestDose Algorithm: STAGE 2 & 3 ║"); + tracing::info!("║ Dual Optimization + Final Predictions ║"); + tracing::info!("╚══════════════════════════════════════════════════════════╝"); + + optimization::dual_optimization(&self) + } + OptimizationStrategy::PosteriorOnly => { + tracing::info!("╔══════════════════════════════════════════════════════════╗"); + tracing::info!("║ BestDose Algorithm: STAGE 2 & 3 ║"); + tracing::info!("║ Posterior Optimization Only + Final Predictions ║"); + tracing::info!("╚══════════════════════════════════════════════════════════╝"); - // STAGE 2 & 3: Dual optimization + predictions - optimization::dual_optimization(&self) + optimization::posterior_optimization(&self) + } + } } /// Set the bias weight (lambda parameter) @@ -786,6 +799,12 @@ impl BestDoseProblem { self } + /// Set the Stage 2 optimization strategy. + pub fn with_optimization_strategy(mut self, strategy: OptimizationStrategy) -> Self { + self.optimization_strategy = strategy; + self + } + /// Get a reference to the refined posterior support points (Θ) pub fn posterior_theta(&self) -> &Theta { &self.theta @@ -811,6 +830,11 @@ impl BestDoseProblem { self.bias_weight } + /// Get the selected Stage 2 optimization strategy. + pub fn optimization_strategy(&self) -> OptimizationStrategy { + self.optimization_strategy + } + /// Get the selected optimization target type pub fn target_type(&self) -> Target { self.target_type diff --git a/src/bestdose/optimization.rs b/src/bestdose/optimization.rs index bd4056ca2..2caefd302 100644 --- a/src/bestdose/optimization.rs +++ b/src/bestdose/optimization.rs @@ -180,6 +180,78 @@ fn run_single_optimization( Ok((full_doses, final_cost)) } +fn finalize_optimization( + problem: &BestDoseProblem, + final_doses: Vec, + final_cost: f64, + method: OptimalMethod, + final_weights: Weights, +) -> Result { + // ═════════════════════════════════════════════════════════════ + // STAGE 3: Final Predictions + // ═════════════════════════════════════════════════════════════ + tracing::info!("─────────────────────────────────────────────────────────────"); + tracing::info!("STAGE 3: Final Predictions"); + tracing::info!("─────────────────────────────────────────────────────────────"); + tracing::info!( + " Calculating predictions with optimal doses and {} weights", + 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)?; + + tracing::info!(" ✓ Predictions complete"); + tracing::info!("─────────────────────────────────────────────────────────────"); + + Ok(BestDoseResult { + optimal_subject, + objf: final_cost, + status: BestDoseStatus::Converged, + preds, + auc_predictions, + optimization_method: method, + }) +} + +pub fn posterior_optimization(problem: &BestDoseProblem) -> Result { + tracing::info!("─────────────────────────────────────────────────────────────"); + tracing::info!("STAGE 2: Posterior Optimization Only"); + tracing::info!("─────────────────────────────────────────────────────────────"); + tracing::info!("│"); + tracing::info!("└─ Optimization: Posterior Weights (Patient-Specific)"); + + let (doses, cost) = run_single_optimization(problem, &problem.posterior, "Posterior")?; + + finalize_optimization( + problem, + doses, + cost, + OptimalMethod::Posterior, + problem.posterior.clone(), + ) +} + /// Stage 2 & 3: Dual optimization + Final predictions /// /// # Algorithm Flow (Matches Diagram) @@ -255,49 +327,5 @@ pub fn dual_optimization(problem: &BestDoseProblem) -> Result { (doses2, cost2, OptimalMethod::Uniform, uniform_weights) }; - // ═════════════════════════════════════════════════════════════ - // STAGE 3: Final Predictions - // ═════════════════════════════════════════════════════════════ - tracing::info!("─────────────────────────────────────────────────────────────"); - tracing::info!("STAGE 3: Final Predictions"); - tracing::info!("─────────────────────────────────────────────────────────────"); - tracing::info!( - " Calculating predictions with optimal doses and {} weights", - 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)?; - - tracing::info!(" ✓ Predictions complete"); - tracing::info!("─────────────────────────────────────────────────────────────"); - - Ok(BestDoseResult { - optimal_subject, - objf: final_cost, - status: BestDoseStatus::Converged, - preds, - auc_predictions, - optimization_method: method, - }) + finalize_optimization(problem, final_doses, final_cost, method, final_weights) } diff --git a/src/bestdose/posterior.rs b/src/bestdose/posterior.rs index 9506d390e..cc9440bf8 100644 --- a/src/bestdose/posterior.rs +++ b/src/bestdose/posterior.rs @@ -53,7 +53,6 @@ use anyhow::Result; use faer::Mat; -use crate::algorithms::npag::burke; use crate::algorithms::npag::NPAG; use crate::algorithms::Algorithms; use crate::algorithms::Status; @@ -86,7 +85,8 @@ const KEEP_UNREFINED_POINTS: bool = true; /// 2. Apply Bayes' rule to get P(θᵢ|data) /// 3. Filter: Keep points where P(θᵢ|data) > 1e-100 × max_weight /// -/// Note: This uses only lambda filtering, NO QR decomposition or second burke call. +/// Note: This uses only direct Bayesian filtering, with no QR decomposition or +/// NPAG weight optimization. /// /// Returns: (filtered_theta, filtered_posterior_weights, filtered_population_weights) pub fn npagfull11_filter( @@ -101,21 +101,81 @@ pub fn npagfull11_filter( // Calculate psi matrix P(data|theta_i) for all support points let psi = calculate_psi(eq, past_data, population_theta, error_models, false)?; - // First burke call to get initial posterior probabilities - let (initial_weights, _) = burke(&psi)?; + // NPAGFULL11 does a direct Bayes update on the incoming prior density: + // posterior_j ∝ prior_j * product_i P(data_i | theta_j). + // It does not run the NPAG/emint optimization path because the Fortran + // routine hard-codes MAXCYC = 0. + let n_points = population_theta.matrix().nrows(); + let mut log_joint_weights = vec![f64::NEG_INFINITY; n_points]; + + for point_idx in 0..n_points { + let prior_weight = population_weights[point_idx]; + if prior_weight <= 0.0 { + continue; + } + + let mut log_weight = prior_weight.ln(); + let mut is_zero = false; + + for subject_idx in 0..psi.matrix().nrows() { + let likelihood = psi.matrix()[(subject_idx, point_idx)]; + if likelihood <= 0.0 { + is_zero = true; + break; + } + log_weight += likelihood.ln(); + } + + if !is_zero { + log_joint_weights[point_idx] = log_weight; + } + } + + let max_log_weight = log_joint_weights + .iter() + .copied() + .fold(f64::NEG_INFINITY, f64::max); + + if !max_log_weight.is_finite() { + return Err(anyhow::anyhow!( + "NPAGFULL11 filtering produced zero joint probability for every support point" + )); + } + + let mut initial_weights: Vec = log_joint_weights + .iter() + .map(|&log_weight| { + if log_weight.is_finite() { + (log_weight - max_log_weight).exp() + } else { + 0.0 + } + }) + .collect(); + + let total_weight: f64 = initial_weights.iter().sum(); + if total_weight <= 0.0 { + return Err(anyhow::anyhow!( + "NPAGFULL11 filtering produced non-positive posterior mass" + )); + } + + for weight in &mut initial_weights { + *weight /= total_weight; + } // NPAGFULL11 filtering: Keep all points within 1e-100 of the maximum weight // This is different from NPAG's condensation - NO QR decomposition here! let max_weight = initial_weights .iter() - .fold(f64::NEG_INFINITY, |a, b| a.max(b)); + .fold(f64::NEG_INFINITY, |a, b| a.max(*b)); let threshold = 1e-100; // NPAGFULL11-specific threshold let keep_lambda: Vec = initial_weights .iter() .enumerate() - .filter(|(_, lam)| *lam > threshold * max_weight) + .filter(|(_, lam)| **lam > threshold * max_weight) .map(|(i, _)| i) .collect(); diff --git a/src/bestdose/types.rs b/src/bestdose/types.rs index e422cd0be..832a642e8 100644 --- a/src/bestdose/types.rs +++ b/src/bestdose/types.rs @@ -183,6 +183,33 @@ impl Default for DoseRange { } } +/// Strategy used for Stage 2 dose optimization. +/// +/// - [`OptimizationStrategy::Dual`]: Run both posterior and uniform optimizations, +/// then keep the lower-cost result. +/// - [`OptimizationStrategy::PosteriorOnly`]: Run only the posterior-weighted +/// optimization path. +#[derive(Debug, Clone, Copy, Serialize, Deserialize, PartialEq, Eq)] +pub enum OptimizationStrategy { + Dual, + PosteriorOnly, +} + +impl Default for OptimizationStrategy { + fn default() -> Self { + Self::Dual + } +} + +impl Display for OptimizationStrategy { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + OptimizationStrategy::Dual => write!(f, "Dual"), + OptimizationStrategy::PosteriorOnly => write!(f, "PosteriorOnly"), + } + } +} + /// The BestDose optimization problem /// /// Contains all data needed for the three-stage BestDose algorithm. @@ -287,6 +314,7 @@ pub struct BestDoseProblem { // Optimization parameters pub(crate) doserange: DoseRange, pub(crate) bias_weight: f64, // λ: 0=personalized, 1=population + pub(crate) optimization_strategy: OptimizationStrategy, } /// Result from BestDose optimization