feat: Implement parametric algorithms

.gitignore

@@ -2,13 +2,21 @@
 /Cargo.lock
 *.log
 theta.csv
-cycles.csv
-pred.csv
 obs.csv
 time.csv
 n_psi.csv
 psi.csv
 r.csv
+correlation.csv
+/docs
+diagnostics.json
+predictions.csv
+summary.csv
+summary.json
+iterations.csv
+population.csv
+shrinkage.csv
+statistics.csv
 posterior.csv
 simulation_output.csv
 /examples/rosuva/*
@@ -17,6 +25,8 @@ simulation_output.csv
 /examples/data/iohexol*
 /examples/data/rosuva*
 /examples/data/vori*
+/examples/paper_benchmarks
+/examples/*/output
 /.idea
 stop
 .vscode
@@ -28,7 +38,13 @@ settings.json
 log.txt
 op.csv
 *results.txt
-covs.csv
+covariates.csv
+individual_effects.csv
+individual_parameters.csv
+residual_error.csv
 error_theta.csv
 lcov.info
-Fortran/
+Fortran/
+paper/
+docs/
+examples/**/outputs/

Cargo.toml

@@ -32,15 +32,22 @@ pharmsol = "=0.26.1"
 anyhow = "1.0.100"
 rayon = "1.10.0"
 rand = "0.10.1"
+cobyla = "0.8.0"
+rand_distr = "0.6.0"
+statrs = "0.18"
 
 [features]
 default = []
+json = ["pharmsol/json"]
 exa = ["pharmsol/exa"]
 
 [profile.release]
 codegen-units = 1
 opt-level = 3
 
+[profile.test]
+inherits = "release"
+
 [dev-dependencies]
 approx = "0.5.1"
 criterion = { version = "0.8" }

benches/bimodal_ke.rs

@@ -20,10 +20,10 @@ fn create_equation() -> equation::ODE {
     )
 }
 
-fn create_parameters() -> Parameters {
-    Parameters::new()
-        .add("ke", 0.001, 3.0)
-        .add("v", 25.0, 250.0)
+fn create_parameter_space() -> ParameterSpace {
+    ParameterSpace::new()
+        .add(ParameterSpec::bounded("ke", 0.001, 3.0))
+        .add(ParameterSpec::bounded("v", 25.0, 250.0))
 }
 
 fn create_error_models() -> Result<AssayErrorModels> {
@@ -37,51 +37,50 @@ fn load_data() -> Result<data::Data> {
     Ok(data::read_pmetrics("examples/bimodal_ke/bimodal_ke.csv")?)
 }
 
-fn setup_with_algorithm(algorithm: Algorithm) -> Result<(Settings, equation::ODE, data::Data)> {
-    let params = create_parameters();
+fn setup_with_algorithm(method: NonparametricMethod) -> Result<EstimationProblem<equation::ODE>> {
     let ems = create_error_models()?;
-
-    let mut settings = Settings::builder()
-        .set_algorithm(algorithm)
-        .set_parameters(params)
-        .set_error_models(ems)
-        .build();
-
-    settings.set_cycles(1000);
-    settings.set_prior(Prior::sobol(2048, 22));
-    settings.disable_output();
-    settings.set_progress(false);
-
+    let observations = ObservationSpec::new()
+        .add_channel(ObservationChannel::continuous(1, "cp"))
+        .with_assay_error_models(ems);
+    let model = ModelDefinition::builder(create_equation())
+        .parameters(create_parameter_space())
+        .observations(observations)
+        .build()?;
     let data = load_data()?;
-    Ok((settings, create_equation(), data))
+    EstimationProblem::builder(model, data)
+        .method(EstimationMethod::Nonparametric(method))
+        .output(OutputPlan::disabled())
+        .runtime(RuntimeOptions {
+            cycles: 1000,
+            progress: false,
+            prior: Some(Prior::sobol(2048, 22)),
+            ..RuntimeOptions::default()
+        })
+        .build()
 }
 
-fn setup_npag() -> Result<(Settings, equation::ODE, data::Data)> {
-    setup_with_algorithm(Algorithm::NPAG)
+fn setup_npag() -> Result<EstimationProblem<equation::ODE>> {
+    setup_with_algorithm(NonparametricMethod::Npag(NpagOptions))
 }
 
-fn setup_npod() -> Result<(Settings, equation::ODE, data::Data)> {
-    setup_with_algorithm(Algorithm::NPOD)
+fn setup_npod() -> Result<EstimationProblem<equation::ODE>> {
+    setup_with_algorithm(NonparametricMethod::Npod(NpodOptions))
 }
 
-fn setup_postprob() -> Result<(Settings, equation::ODE, data::Data)> {
-    setup_with_algorithm(Algorithm::POSTPROB)
+fn setup_postprob() -> Result<EstimationProblem<equation::ODE>> {
+    setup_with_algorithm(NonparametricMethod::Postprob(PostProbOptions))
 }
 
 fn benchmark_algorithm<F>(c: &mut Criterion, bench_name: &str, setup_fn: F)
 where
-    F: Fn() -> Result<(Settings, equation::ODE, data::Data)>,
+    F: Fn() -> Result<EstimationProblem<equation::ODE>>,
 {
-    let (settings, eq, data) = setup_fn().unwrap();
+    let problem = setup_fn().unwrap();
 
     c.bench_function(bench_name, |b| {
         b.iter_with_setup(
-            || (settings.clone(), eq.clone(), data.clone()),
-            |(s, e, d)| {
-                let mut algorithm = dispatch_algorithm(s, e, d).unwrap();
-                let result = algorithm.fit().unwrap();
-                black_box(result)
-            },
+            || problem.clone(),
+            |problem| black_box(problem.run().unwrap()),
         )
     });
 }

examples/bestdose.rs

@@ -1,8 +1,6 @@
 use anyhow::Result;
-use pmcore::bestdose::{BestDosePosterior, DoseRange, Target};
-
+use pmcore::bestdose::{BestDoseConfig, BestDosePosterior, DoseRange, Target};
 use pmcore::prelude::*;
-use pmcore::routines::initialization::parse_prior;
 
 fn main() -> Result<()> {
     // Example model
@@ -18,23 +16,15 @@ fn main() -> Result<()> {
         },
     };
 
-    let params = Parameters::new()
-        .add("ke", 0.001, 3.0)
-        .add("v", 25.0, 250.0);
+    let parameter_space = ParameterSpace::new()
+        .add(ParameterSpec::bounded("ke", 0.001, 3.0))
+        .add(ParameterSpec::bounded("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();
+    let config = BestDoseConfig::new(parameter_space.clone(), ems.clone()).with_progress(false);
 
     // Generate a patient with known parameters
     // Ke = 0.5, V = 50
@@ -71,21 +61,14 @@ fn main() -> Result<()> {
         .observation(18.0, conc(6.0, 75.0) + conc(18.0, 150.0), 0)
         .build();
 
-    let (theta, prior) = parse_prior(
-        &"examples/bimodal_ke/output/theta.csv".to_string(),
-        &settings,
-    )
-    .unwrap();
+    let (theta, prior) = read_prior("examples/bimodal_ke/output/theta.csv", &parameter_space)?;
 
-    // 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
         eq.clone(),
-        settings.clone(),
+        config.clone(),
     )?;
 
     println!("Optimizing dose...");

examples/bestdose_auc.rs

@@ -1,10 +1,11 @@
 use anyhow::Result;
-use pmcore::bestdose::{BestDosePosterior, DoseRange, Target};
+use pmcore::bestdose::{BestDoseConfig, BestDosePosterior, DoseRange, Target};
 use pmcore::prelude::*;
-use pmcore::routines::initialization::parse_prior;
 
 fn main() -> Result<()> {
-    tracing_subscriber::fmt::init();
+    tracing_subscriber::fmt()
+        .with_env_filter(tracing_subscriber::EnvFilter::new("info,diffsol=off"))
+        .init();
 
     println!("BestDose AUC Target - Minimal Example\n");
     println!("======================================\n");
@@ -22,30 +23,22 @@ fn main() -> Result<()> {
     };
 
     // Minimal parameter ranges
-    let params = Parameters::new()
-        .add("ke", 0.001, 3.0)
-        .add("v", 25.0, 250.0);
+    let parameter_space = ParameterSpace::new()
+        .add(ParameterSpec::bounded("ke", 0.001, 3.0))
+        .add(ParameterSpec::bounded("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),
     )?;
 
-    let mut settings = Settings::builder()
-        .set_algorithm(Algorithm::NPAG)
-        .set_parameters(params)
-        .set_error_models(ems.clone())
-        .build();
-
-    settings.disable_output();
-    settings.set_idelta(60.0); // 1 hour intervals for AUC calculation
+    let config = BestDoseConfig::new(parameter_space.clone(), ems.clone())
+        .with_progress(false)
+        .with_prediction_interval(60.0);
 
     // Load realistic prior from previous NPAG run (47 support points)
     println!("Loading prior from bimodal_ke example...");
-    let (theta, prior) = parse_prior(
-        &"examples/bimodal_ke/output/theta.csv".to_string(),
-        &settings,
-    )?;
+    let (theta, prior) = read_prior("examples/bimodal_ke/output/theta.csv", &parameter_space)?;
     let weights = prior.as_ref().unwrap();
 
     println!("Prior: {} support points\n", theta.matrix().nrows());
@@ -67,16 +60,16 @@ fn main() -> Result<()> {
         weights,
         None, // No past data - use prior directly
         eq.clone(),
-        settings.clone(),
+        config.clone(),
     )?;
 
     println!("Optimizing dose...\n");
     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
+        DoseRange::new(100.0, 2000.0),
+        0.8,
+        Target::AUCFromZero,
     )?;
 
     let opt_doses = optimal.doses();
@@ -143,7 +136,7 @@ fn main() -> Result<()> {
         None,
         DoseRange::new(50.0, 500.0),
         0.8,
-        Target::AUCFromLastDose, // Interval AUC from last dose!
+        Target::AUCFromLastDose,
     )?;
 
     let doses: Vec<f64> = optimal_interval.doses();