Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
23 changes: 5 additions & 18 deletions examples/bestdose.rs
Original file line number Diff line number Diff line change
Expand Up @@ -114,35 +114,22 @@ 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::<Vec<_>>()
})
.collect::<Vec<f64>>();
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
optimal.objf(),
optimal.objf().ln(),
optimal.optimization_method()
);
}

// 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.preds.predictions().into_iter() {
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()
Expand Down
48 changes: 9 additions & 39 deletions examples/bestdose_auc.rs
Original file line number Diff line number Diff line change
Expand Up @@ -83,33 +83,20 @@ 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::<Vec<_>>()
})
.collect::<Vec<f64>>();
let opt_doses = optimal.doses();

println!("=== RESULTS ===");
println!("Optimal dose: {:.1} mg", opt_doses[0]);
println!("Cost function: {:.6}", optimal.objf);
println!("Cost function: {:.6}", optimal.objf());

if let Some(auc_preds) = &optimal.auc_predictions {
if let Some(auc_preds) = &optimal.auc_predictions() {
println!("\nAUC Predictions:");
let mut total_error = 0.0;
for (time, auc) in auc_preds {
// Find the target AUC for this time
let target = if (*time - 6.0).abs() < 0.1 {
let target = if (time - 6.0).abs() < 0.1 {
50.0
} else if (*time - 12.0).abs() < 0.1 {
} else if (time - 12.0).abs() < 0.1 {
Comment thread
mhovd marked this conversation as resolved.
80.0
} else {
0.0
Expand All @@ -127,7 +114,7 @@ fn main() -> Result<()> {
);
} else {
println!("\nConcentration Predictions:");
for pred in optimal.preds.predictions() {
for pred in optimal.predictions().predictions() {
println!(
" Time: {:5.1}h | Target: {:6.1} | Predicted: {:6.2}",
pred.time(),
Expand Down Expand Up @@ -172,30 +159,13 @@ fn main() -> Result<()> {
println!("Optimizing maintenance dose...\n");
let optimal_interval = problem_interval.optimize()?;

let doses: Vec<f64> = optimal_interval
.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> = optimal_interval.doses();

println!("=== INTERVAL AUC RESULTS ===");
println!("Optimal maintenance dose (at t=12h): {:.1} mg", doses[0]);
println!("Cost function: {:.6}", optimal_interval.objf);
println!("Cost function: {:.6}", optimal_interval.objf());

if let Some(auc_preds) = &optimal_interval.auc_predictions {
if let Some(auc_preds) = &optimal_interval.auc_predictions() {
println!("\nInterval AUC Predictions:");
for (time, auc) in auc_preds {
let target = 60.0;
Expand Down
7 changes: 5 additions & 2 deletions examples/bestdose_bounds.rs
Original file line number Diff line number Diff line change
Expand Up @@ -89,7 +89,7 @@ fn main() -> Result<()> {
let result = problem.optimize()?;

let doses: Vec<f64> = result
.optimal_subject
.optimal_subject()
.iter()
.map(|occ| {
occ.iter()
Expand Down Expand Up @@ -118,7 +118,10 @@ fn main() -> Result<()> {

println!(
"{:<30} | {:>10.1} mg | {:>10.6}{}",
description, doses[0], result.objf, at_bound
description,
doses[0],
result.objf(),
at_bound
);
}

Expand Down
15 changes: 10 additions & 5 deletions src/bestdose/optimization.rs
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ use argmin::solver::neldermead::NelderMead;

use crate::bestdose::cost::calculate_cost;
use crate::bestdose::predictions::calculate_final_predictions;
use crate::bestdose::types::{BestDoseProblem, BestDoseResult};
use crate::bestdose::types::{BestDoseProblem, BestDoseResult, BestDoseStatus, OptimalMethod};
use crate::structs::weights::Weights;
use pharmsol::prelude::*;

Expand Down Expand Up @@ -244,10 +244,15 @@ pub fn dual_optimization(problem: &BestDoseProblem) -> Result<BestDoseResult> {

let (final_doses, final_cost, method, final_weights) = if cost1 <= cost2 {
tracing::info!(" → Winner: Posterior (lower cost) ✓");
(doses1, cost1, "posterior", problem.posterior.clone())
(
doses1,
cost1,
OptimalMethod::Posterior,
problem.posterior.clone(),
)
} else {
tracing::info!(" → Winner: Uniform (lower cost) ✓");
(doses2, cost2, "uniform", uniform_weights)
(doses2, cost2, OptimalMethod::Uniform, uniform_weights)
};

// ═════════════════════════════════════════════════════════════
Expand Down Expand Up @@ -290,9 +295,9 @@ pub fn dual_optimization(problem: &BestDoseProblem) -> Result<BestDoseResult> {
Ok(BestDoseResult {
optimal_subject,
objf: final_cost,
status: "Converged".to_string(),
status: BestDoseStatus::Converged,
preds,
auc_predictions,
optimization_method: method.to_string(),
optimization_method: method,
})
}
108 changes: 95 additions & 13 deletions src/bestdose/types.rs
Original file line number Diff line number Diff line change
Expand Up @@ -6,12 +6,15 @@
//! - [`Target`]: Enum specifying concentration or AUC targets
//! - [`DoseRange`]: Dose constraint specification

use std::fmt::Display;

use crate::prelude::*;
use crate::routines::output::predictions::NPPredictions;
use crate::routines::settings::Settings;
use crate::structs::theta::Theta;
use crate::structs::weights::Weights;
use pharmsol::prelude::*;
use serde::{Deserialize, Serialize};

/// Target type for dose optimization
///
Expand Down Expand Up @@ -49,7 +52,7 @@ use pharmsol::prelude::*;
/// - Automatically finds the most recent bolus/infusion before each observation
///
/// Both methods use trapezoidal rule on a dense time grid controlled by `settings.predictions().idelta`.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]
pub enum Target {
/// Target concentrations at observation times
///
Expand Down Expand Up @@ -133,10 +136,10 @@ pub enum Target {
/// ```rust,ignore
/// use pmcore::bestdose::DoseRange;
///
/// // Standard range: 0-1000 mg
/// // Large range: 0-1000 mg
/// let range = DoseRange::new(0.0, 1000.0);
///
/// // Narrow therapeutic window
/// // Narrow range: 50-150 mg
/// let range = DoseRange::new(50.0, 150.0);
///
/// // Access bounds
Expand Down Expand Up @@ -357,45 +360,124 @@ pub struct BestDoseProblem {
/// # Ok(())
/// # }
/// ```
#[derive(Debug, Clone)]
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct BestDoseResult {
/// Optimal dose amount(s)
/// Subject with optimal doses
///
/// Vector contains one element per dose in the target subject.
/// Order matches the dose events in the target subject.
pub optimal_subject: Subject,
/// The [Subject] contains the same events as the target subject,
/// but with the dose amounts updated to the optimal values.
pub(crate) optimal_subject: Subject,

/// Final cost function value
///
/// Lower is better. Represents the weighted combination of variance
/// (patient-specific error) and bias (deviation from population).
pub objf: f64,
pub(crate) objf: f64,

/// Optimization status message
///
/// Examples: "converged", "maximum iterations reached", etc.
pub status: String,
pub(crate) status: BestDoseStatus,

/// Concentration-time predictions for optimal doses
///
/// Contains predicted concentrations at observation times using the
/// optimal doses. Predictions use the weights from the winning optimization
/// method (posterior or uniform).
pub preds: NPPredictions,
pub(crate) preds: NPPredictions,

/// AUC values at observation times
///
/// Only populated when `target_type` is [`Target::AUC`].
/// Each tuple contains `(time, cumulative_auc)`.
///
/// For [`Target::Concentration`], this field is `None`.
pub auc_predictions: Option<Vec<(f64, f64)>>,
pub(crate) auc_predictions: Option<Vec<(f64, f64)>>,

/// Which optimization method produced the best result
///
/// - `"posterior"`: Patient-specific optimization (uses posterior weights)
/// - `"uniform"`: Population-based optimization (uses uniform weights)
///
/// The algorithm runs both optimizations and selects the one with lower cost.
pub optimization_method: String,
pub(crate) optimization_method: OptimalMethod,
}

impl BestDoseResult {
/// Get the optimized subject
pub fn optimal_subject(&self) -> &Subject {
&self.optimal_subject
}

/// Get the dose amounts of the optimized subject
///
/// This includes all doses (bolus and infusion) in the order they appear
/// in the optimal subject, and returns their amounts as a vector of f64.
pub fn doses(&self) -> Vec<f64> {
self.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::<Vec<_>>()
})
.collect::<Vec<f64>>()
}

/// Get the objective cost function value
pub fn objf(&self) -> f64 {
self.objf
}

/// Get the optimization status
pub fn status(&self) -> &BestDoseStatus {
&self.status
}

/// Get the concentration-time predictions
pub fn predictions(&self) -> &NPPredictions {
&self.preds
}

/// Get the AUC predictions, if available
pub fn auc_predictions(&self) -> Option<Vec<(f64, f64)>> {
self.auc_predictions.clone()
}

/// Get the optimization method used
pub fn optimization_method(&self) -> OptimalMethod {
self.optimization_method
}
}

/// Optimization method used in BestDose
///
/// This returns the type of optimization method that produced the best result:
/// - `Posterior`: Patient-specific optimization using posterior weights
/// - `Uniform`: Population-based optimization using uniform weights
#[derive(Debug, Clone, Serialize, Deserialize, Copy, PartialEq, Eq)]
pub enum OptimalMethod {
Posterior,
Uniform,
}

impl Display for OptimalMethod {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
match self {
OptimalMethod::Posterior => write!(f, "Posterior"),
OptimalMethod::Uniform => write!(f, "Uniform"),
}
}
}

/// Status of the BestDose optimization
#[derive(Debug, Clone, Serialize, Deserialize, Copy, PartialEq, Eq)]
pub enum BestDoseStatus {
Converged,
MaxIterations,
}
Loading
Loading