feat: new bestdose API

examples/bestdose.rs

@@ -1,6 +1,5 @@
 use anyhow::Result;
-use pmcore::bestdose; // bestdose new
-                      // use pmcore::bestdose::bestdose_old as bestdose; // bestdose old
+use pmcore::bestdose::{BestDosePosterior, DoseRange, Target};
 
 use pmcore::prelude::*;
 use pmcore::routines::initialization::parse_prior;
@@ -78,22 +77,15 @@ fn main() -> Result<()> {
     )
     .unwrap();
 
-    // Example usage - using new() constructor which calculates NPAGFULL11 posterior
-    // max_cycles controls NPAGFULL refinement:
-    //   0 = NPAGFULL11 only (fast but less accurate)
-    //   100 = moderate refinement
-    //   500 = full refinement (Fortran default, slow but most accurate)
-    let problem = bestdose::BestDoseProblem::new(
+    // Example usage - two-stage API:
+    // Stage 1: Compute posterior (expensive, done once)
+    // Stage 2: Optimize doses (can be called multiple times with different params)
+    let posterior = BestDosePosterior::compute(
         &theta,
         &prior.unwrap(),
         Some(past_data.clone()), // Optional: past data for Bayesian updating
-        target_data.clone(),
-        None,
         eq.clone(),
-        bestdose::DoseRange::new(0.0, 300.0),
-        0.0,
         settings.clone(),
-        bestdose::Target::Concentration, // Target concentrations (not AUCs)
     )?;
 
     println!("Optimizing dose...");
@@ -103,7 +95,13 @@ fn main() -> Result<()> {
 
     for bias_weight in &bias_weights {
         println!("Running optimization with bias weight: {}", bias_weight);
-        let optimal = problem.clone().with_bias_weight(*bias_weight).optimize()?;
+        let optimal = posterior.optimize(
+            target_data.clone(),
+            None,
+            DoseRange::new(0.0, 300.0),
+            *bias_weight,
+            Target::Concentration,
+        )?;
         results.push((bias_weight, optimal));
     }
 
@@ -124,7 +122,7 @@ fn main() -> Result<()> {
     // Print concentration-time predictions for the optimal dose
     let optimal = &results.last().unwrap().1;
     println!("\nConcentration-time predictions for optimal dose:");
-    for pred in optimal.predictions().predictions().into_iter() {
+    for pred in optimal.predictions().predictions().iter() {
         println!(
             "Time: {:.2} h, Observed: {:.2}, (Pop Mean: {:.4}, Pop Median: {:.4}, Post Mean: {:.4}, Post Median: {:.4})",
             pred.time(), pred.obs().unwrap_or(0.0), pred.pop_mean(), pred.pop_median(), pred.post_mean(), pred.post_median()

examples/bestdose_auc.rs

@@ -1,5 +1,5 @@
 use anyhow::Result;
-use pmcore::bestdose::{BestDoseProblem, DoseRange, Target};
+use pmcore::bestdose::{BestDosePosterior, DoseRange, Target};
 use pmcore::prelude::*;
 use pmcore::routines::initialization::parse_prior;
 
@@ -61,22 +61,23 @@ fn main() -> Result<()> {
         .observation(12.0, 80.0, 0) // Target AUC at 12h
         .build();
 
-    println!("Creating BestDose problem with AUC targets...");
-    let problem = BestDoseProblem::new(
+    println!("Creating BestDose posterior (no past data - use prior directly)...");
+    let posterior = BestDosePosterior::compute(
         &theta,
         weights,
         None, // No past data - use prior directly
-        target_data.clone(),
-        None,
         eq.clone(),
-        DoseRange::new(100.0, 2000.0), // Wider range for AUC targets
-        0.8,                           // for AUC targets higher bias_weight usually works best
         settings.clone(),
-        Target::AUCFromZero, // Cumulative AUC from time 0
     )?;
 
     println!("Optimizing dose...\n");
-    let optimal = problem.optimize()?;
+    let optimal = posterior.optimize(
+        target_data.clone(),
+        None,
+        DoseRange::new(100.0, 2000.0), // Wider range for AUC targets
+        0.8,                           // for AUC targets higher bias_weight usually works best
+        Target::AUCFromZero,           // Cumulative AUC from time 0
+    )?;
 
     let opt_doses = optimal.doses();
 
@@ -137,22 +138,14 @@ fn main() -> Result<()> {
         .build();
 
     println!("Creating BestDose problem with interval AUC target...");
-    let problem_interval = BestDoseProblem::new(
-        &theta,
-        weights,
-        None,
+    let optimal_interval = posterior.optimize(
         target_interval.clone(),
         None,
-        eq.clone(),
         DoseRange::new(50.0, 500.0),
         0.8,
-        settings.clone(),
         Target::AUCFromLastDose, // Interval AUC from last dose!
     )?;
 
-    println!("Optimizing maintenance dose...\n");
-    let optimal_interval = problem_interval.optimize()?;
-
     let doses: Vec<f64> = optimal_interval.doses();
 
     println!("=== INTERVAL AUC RESULTS ===");

examples/bestdose_bounds.rs

@@ -1,5 +1,5 @@
 use anyhow::Result;
-use pmcore::bestdose::{BestDoseProblem, DoseRange, Target};
+use pmcore::bestdose::{BestDosePosterior, DoseRange, Target};
 use pmcore::prelude::*;
 use pmcore::routines::initialization::parse_prior;
 
@@ -67,40 +67,20 @@ fn main() -> Result<()> {
     println!("{:<30} | {:>12} | {:>10}", "Range", "Optimal Dose", "Cost");
     println!("{}", "-".repeat(60));
 
+    // Compute posterior once, reuse for all dose ranges
+    let posterior =
+        BestDosePosterior::compute(&theta, weights, None, eq.clone(), settings.clone())?;
+
     for (min, max, description) in dose_ranges {
-        let problem = BestDoseProblem::new(
-            &theta,
-            weights,
-            None,
+        let result = posterior.optimize(
             target_data.clone(),
             None,
-            eq.clone(),
             DoseRange::new(min, max),
             0.5,
-            settings.clone(),
             Target::Concentration,
         )?;
 
-        let result = problem.optimize()?;
-
-        let doses: Vec<f64> = result
-            .optimal_subject()
-            .iter()
-            .map(|occ| {
-                occ.iter()
-                    .filter(|event| match event {
-                        Event::Bolus(_) => true,
-                        Event::Infusion(_) => true,
-                        _ => false,
-                    })
-                    .map(|event| match event {
-                        Event::Bolus(bolus) => bolus.amount(),
-                        Event::Infusion(infusion) => infusion.amount(),
-                        _ => 0.0,
-                    })
-            })
-            .flatten()
-            .collect();
+        let doses: Vec<f64> = result.doses();
 
         // Check if dose hit the bound
         let at_bound = if (doses[0] - max).abs() < 1.0 {

src/bestdose/cost.rs

@@ -129,7 +129,7 @@ use pharmsol::Equation;
 /// - Model simulation fails
 /// - Prediction length doesn't match observation count
 /// - AUC calculation fails (for AUC targets)
-pub fn calculate_cost(problem: &BestDoseProblem, candidate_doses: &[f64]) -> Result<f64> {
+pub(crate) fn calculate_cost(problem: &BestDoseProblem, candidate_doses: &[f64]) -> Result<f64> {
     // Validate candidate_doses length matches expected optimizable dose count
     let expected_optimizable = problem
         .target
@@ -233,7 +233,7 @@ pub fn calculate_cost(problem: &BestDoseProblem, candidate_doses: &[f64]) -> Res
 
     // Calculate variance (using posterior weights) and population mean (using prior weights)
 
-    for ((row, post_prob), prior_prob) in problem
+    for ((row, post_prob), _prior_prob) in problem
         .theta
         .matrix()
         .row_iter()
@@ -509,8 +509,8 @@ pub fn calculate_cost(problem: &BestDoseProblem, candidate_doses: &[f64]) -> Res
             let pj = preds_i[j];
             let se = (obs_val - pj).powi(2);
             sumsq_i += se;
-            // Calculate population mean using PRIOR probabilities
-            y_bar[j] += prior_prob * pj;
+            // Calculate population mean using POSTERIOR probabilities
+            y_bar[j] += post_prob * pj;
         }
 
         variance += post_prob * sumsq_i; // Weighted by posterior