Extract unified platform structure

.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

@@ -28,19 +28,31 @@ tracing-subscriber = { version = "0.3.19", features = [
     "time",
 ] }
 faer = "0.24.0"
-pharmsol = "=0.26.1"
+pharmsol = { version = "=0.27.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 = []
-exa = ["pharmsol/exa"]
+json = []
+exa = []
+dsl-core = ["pharmsol/dsl-core"]
+dsl-jit = ["pharmsol/dsl-jit"]
+dsl-aot = ["pharmsol/dsl-aot"]
+dsl-aot-load = ["pharmsol/dsl-aot-load"]
+dsl-wasm = ["pharmsol/dsl-wasm"]
 
 [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

@@ -5,83 +5,80 @@ use pmcore::prelude::*;
 use std::hint::black_box;
 
 fn create_equation() -> equation::ODE {
-    equation::ODE::new(
-        |x, p, _t, dx, b, rateiv, _cov| {
-            fetch_params!(p, ke, _v);
-            dx[0] = -ke * x[0] + rateiv[1] + b[1];
+    ode! {
+        name: "bimodal_ke",
+        params: [ke, v],
+        states: [central],
+        outputs: [outeq_1],
+        routes: [
+            infusion(input_1) -> central,
+        ],
+        diffeq: |x, _t, dx| {
+            dx[central] = -ke * x[central];
         },
-        |_p, _t, _cov| lag! {},
-        |_p, _t, _cov| fa! {},
-        |_p, _t, _cov, _x| {},
-        |x, p, _t, _cov, y| {
-            fetch_params!(p, _ke, v);
-            y[1] = x[0] / v;
+        out: |x, _t, y| {
+            y[outeq_1] = x[central] / v;
         },
-    )
+    }
 }
 
-fn create_parameters() -> Parameters {
-    Parameters::new()
-        .add("ke", 0.001, 3.0)
-        .add("v", 25.0, 250.0)
-}
-
-fn create_error_models() -> Result<AssayErrorModels> {
-    Ok(AssayErrorModels::new().add(
-        1,
-        AssayErrorModel::additive(ErrorPoly::new(0.0, 0.5, 0.0, 0.0), 0.0),
-    )?)
+fn create_error_model() -> AssayErrorModel {
+    AssayErrorModel::additive(ErrorPoly::new(0.0, 0.5, 0.0, 0.0), 0.0)
 }
 
 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();
-    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);
-
+fn setup_npag() -> Result<EstimationProblem<equation::ODE>> {
     let data = load_data()?;
-    Ok((settings, create_equation(), data))
+    EstimationProblem::builder(create_equation(), data)
+        .parameter(Parameter::bounded("ke", 0.001, 3.0))?
+        .parameter(Parameter::bounded("v", 25.0, 250.0))?
+        .method(Npag::new())
+        .error("outeq_1", create_error_model())?
+        .cycles(1000)
+        .progress(false)
+        .prior(Prior::sobol(2048, 22))
+        .build()
 }
 
-fn setup_npag() -> Result<(Settings, equation::ODE, data::Data)> {
-    setup_with_algorithm(Algorithm::NPAG)
-}
-
-fn setup_npod() -> Result<(Settings, equation::ODE, data::Data)> {
-    setup_with_algorithm(Algorithm::NPOD)
+fn setup_npod() -> Result<EstimationProblem<equation::ODE>> {
+    let data = load_data()?;
+    EstimationProblem::builder(create_equation(), data)
+        .parameter(Parameter::bounded("ke", 0.001, 3.0))?
+        .parameter(Parameter::bounded("v", 25.0, 250.0))?
+        .method(Npod::new())
+        .error("outeq_1", create_error_model())?
+        .cycles(1000)
+        .progress(false)
+        .prior(Prior::sobol(2048, 22))
+        .build()
 }
 
-fn setup_postprob() -> Result<(Settings, equation::ODE, data::Data)> {
-    setup_with_algorithm(Algorithm::POSTPROB)
+fn setup_postprob() -> Result<EstimationProblem<equation::ODE>> {
+    let data = load_data()?;
+    EstimationProblem::builder(create_equation(), data)
+        .parameter(Parameter::bounded("ke", 0.001, 3.0))?
+        .parameter(Parameter::bounded("v", 25.0, 250.0))?
+        .method(PostProb::new())
+        .error("outeq_1", create_error_model())?
+        .cycles(1000)
+        .progress(false)
+        .prior(Prior::sobol(2048, 22))
+        .build()
 }
 
 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.fit().unwrap()),
         )
     });
 }

examples/bestdose.rs

@@ -1,40 +1,34 @@
 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
     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];
+        name: "bestdose_one_compartment",
+        params: [ke, v],
+        states: [central],
+        outputs: [cp],
+        routes: [
+            bolus(dose) -> central,
+        ],
+        diffeq: |x, _t, dx| {
+            dx[central] = -ke * x[central];
         },
-        out: |x, p, _t, _cov, y| {
-            fetch_params!(p, _ke, v);
-            y[0] = x[0] / v;
+        out: |x, _t, y| {
+            y[cp] = x[central] / v;
         },
     };
 
-    let params = Parameters::new()
-        .add("ke", 0.001, 3.0)
-        .add("v", 25.0, 250.0);
+    let parameter_space = ParameterSpace::new()
+        .add(Parameter::bounded("ke", 0.001, 3.0))
+        .add(Parameter::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 +65,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...");