Fix bestdose posterior

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 &param_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(())
 }

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 &param_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(())
+}

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

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<BestDoseResult> {
-        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