feat: Chain fit methods using result

benches/bimodal_ke.rs

@@ -111,4 +111,3 @@ criterion_group! {
     targets = benchmark_bimodal_ke_npag, benchmark_bimodal_ke_npod, benchmark_bimodal_ke_postprob
 }
 criterion_main!(benches);
-

examples/bestdose_auc.rs

@@ -43,10 +43,8 @@ fn main() -> Result<()> {
 
     // Load realistic prior from previous NPAG run (47 support points)
     println!("Loading prior from bimodal_ke example...");
-    let (theta, prior) = Theta::from_file(
-        "examples/bimodal_ke/output/theta.csv",
-        &parameter_space,
-    )?;
+    let (theta, prior) =
+        Theta::from_file("examples/bimodal_ke/output/theta.csv", &parameter_space)?;
     let weights = prior.as_ref().unwrap();
 
     println!("Prior: {} support points\n", theta.matrix().nrows());

examples/bestdose_bounds.rs

@@ -40,10 +40,8 @@ fn main() -> Result<()> {
 
     // Load realistic prior from previous NPAG run
     println!("Loading prior from bimodal_ke example...");
-    let (theta, prior) = Theta::from_file(
-        "examples/bimodal_ke/output/theta.csv",
-        &parameter_space,
-    )?;
+    let (theta, prior) =
+        Theta::from_file("examples/bimodal_ke/output/theta.csv", &parameter_space)?;
     let weights = prior.as_ref().unwrap();
 
     println!("Prior: {} support points\n", theta.matrix().nrows());

examples/bestdose_cov.rs

@@ -90,10 +90,8 @@ fn main() -> Result<()> {
         .observation(18.0, conc(6.0, 75.0) + conc(18.0, 150.0), 0)
         .build();
 
-    let (mut theta, prior) = Theta::from_file(
-        "examples/bimodal_ke/output/theta.csv",
-        &parameter_space,
-    )?;
+    let (mut theta, prior) =
+        Theta::from_file("examples/bimodal_ke/output/theta.csv", &parameter_space)?;
 
     let m_t = theta.matrix_mut();
     for i in 0..m_t.nrows() {

examples/bimodal_ke_backend_compare.rs

@@ -313,4 +313,3 @@ fn print_summary(results: &[ComparisonResult]) -> Result<()> {
 
     Ok(())
 }
-

examples/bootstrap_model.rs

@@ -0,0 +1,104 @@
+//! Model complexity bootstrapping: expand a 1-compartment posterior into a
+//! prior for a 2-compartment model using `bootstrap_theta()`.
+//!
+//! Stage 1 fits (ke, v) with NPAG. Stage 2 expands the posterior into a
+//! prior over (ke, v, kcp, kpc) and fits the 2-compartment model.
+//!
+//! Uses inline synthetic data to keep the example self-contained.
+
+use pharmsol::SubjectBuilderExt;
+use pmcore::prelude::*;
+
+fn main() -> Result<()> {
+    // ── Synthetic data ──
+    let subject = Subject::builder("1")
+        .bolus(0.0, 100.0, 1)
+        .observation(1.0, 8.0, 1)
+        .observation(2.0, 5.0, 1)
+        .observation(4.0, 2.0, 1)
+        .build();
+    let data = Data::new(vec![subject]);
+
+    // ── Stage 1: 1-compartment ODE ──
+
+    let ode_1cmt = ode! {
+        name: "one_cmt_bootstrap",
+        params: [ke, v],
+        states: [central],
+        outputs: [outeq_1],
+        routes: [bolus(input_1) -> central],
+        diffeq: |x, p, _t, dx, _cov| {
+            fetch_params!(p, ke);
+            dx[0] = -ke * x[0];
+        },
+        out: |x, p, _t, _cov, y| {
+            fetch_params!(p, v);
+            y[0] = x[0] / v;
+        },
+    };
+
+    let params_1cmt = ParameterSpace::bounded()
+        .add("ke", 0.001, 3.0)
+        .add("v", 5.0, 50.0);
+    let prior = Theta::sobol(&params_1cmt, 100)?;
+    let error_models = AssayErrorModels::new().add(
+        "outeq_1",
+        AssayErrorModel::additive(ErrorPoly::new(0.0, 0.5, 0.0, 0.0), 0.0),
+    )?;
+
+    println!("Stage 1: Fitting 1-cmt model...");
+    let r1 = EstimationProblem::nonparametric(ode_1cmt, data.clone(), prior, error_models.clone())?
+        .fit_with(NonParametricAlgorithm::npag())?;
+
+    println!(
+        "  {} support points, OBJF = {:.2}",
+        r1.get_theta().nspp(),
+        r1.objf(),
+    );
+
+    // ── Stage 2: bootstrap into 2-compartment prior ──
+
+    let ode_2cmt = ode! {
+        name: "two_cmt_bootstrap",
+        params: [ke, v, kcp, kpc],
+        states: [central, peripheral],
+        outputs: [outeq_1],
+        routes: [bolus(input_1) -> central],
+        diffeq: |x, p, _t, dx, _cov| {
+            fetch_params!(p, ke, _v, kcp, kpc);
+            dx[0] = -ke * x[0] - kcp * x[0] + kpc * x[1];
+            dx[1] = -kpc * x[1] + kcp * x[0];
+        },
+        out: |x, p, _t, _cov, y| {
+            fetch_params!(p, _ke, v);
+            y[0] = x[0] / v;
+        },
+    };
+
+    let params_2cmt = ParameterSpace::bounded()
+        .add("ke", 0.001, 3.0)
+        .add("v", 5.0, 50.0)
+        .add("kcp", 0.01, 5.0)
+        .add("kpc", 0.01, 5.0);
+
+    // Expand posterior: common params keep their values, new ones are Sobol-sampled
+    let prior_2cmt = bootstrap_theta(&r1, &ode_2cmt, &params_2cmt, 500)?;
+
+    println!(
+        "Prior expanded: {} support points over {:?}",
+        prior_2cmt.nspp(),
+        prior_2cmt.parameters().names(),
+    );
+
+    println!("Stage 2: Fitting 2-cmt model...");
+    let r2 = EstimationProblem::nonparametric(ode_2cmt, data, prior_2cmt, error_models)?
+        .fit_with(NonParametricAlgorithm::npag())?;
+
+    println!(
+        "  {} support points, OBJF = {:.2}",
+        r2.get_theta().nspp(),
+        r2.objf(),
+    );
+
+    Ok(())
+}

examples/chain_algorithms.rs

@@ -0,0 +1,65 @@
+//! Algorithm chaining: NPAG(10) → NPOD(5) on a simple 1-compartment IV bolus model.
+//!
+//! Demonstrates `NonParametricResult::chain()`.
+//! Uses inline synthetic data to keep the example self-contained.
+
+use pharmsol::SubjectBuilderExt;
+use pmcore::prelude::*;
+
+fn main() -> Result<()> {
+    let eq = ode! {
+        name: "chain_example",
+        params: [ke, v],
+        states: [central],
+        outputs: [outeq_1],
+        routes: [bolus(input_1) -> central],
+        diffeq: |x, p, _t, dx, _cov| {
+            fetch_params!(p, ke);
+            dx[0] = -ke * x[0];
+        },
+        out: |x, p, _t, _cov, y| {
+            fetch_params!(p, v);
+            y[0] = x[0] / v;
+        },
+    };
+
+    // Small synthetic dataset: 2 subjects, 3 obs each
+    let subject1 = Subject::builder("1")
+        .bolus(0.0, 100.0, 1)
+        .observation(1.0, 8.0, 1)
+        .observation(2.0, 5.0, 1)
+        .observation(4.0, 2.0, 1)
+        .build();
+
+    let subject2 = Subject::builder("2")
+        .bolus(0.0, 80.0, 1)
+        .observation(1.0, 6.0, 1)
+        .observation(2.0, 4.0, 1)
+        .observation(4.0, 1.5, 1)
+        .build();
+
+    let data = Data::new(vec![subject1, subject2]);
+
+    let parameters = ParameterSpace::bounded()
+        .add("ke", 0.001, 3.0)
+        .add("v", 5.0, 50.0);
+    let prior = Theta::sobol(&parameters, 100)?;
+    let error_models = AssayErrorModels::new().add(
+        "outeq_1",
+        AssayErrorModel::additive(ErrorPoly::new(0.0, 0.5, 0.0, 0.0), 0.0),
+    )?;
+
+    // Chain: NPAG runs first, then NPOD continues from its result
+    let result = EstimationProblem::nonparametric(eq, data, prior, error_models)?
+        .fit_with(NpagConfig::new().max_cycles(10))?
+        .chain(NpodConfig::new().max_cycles(5))?;
+
+    println!(
+        "Chained NPAG→NPOD: OBJF = {:.2}, {} support points, {} total cycles",
+        result.objf(),
+        result.get_theta().nspp(),
+        result.cycles(),
+    );
+
+    Ok(())
+}

examples/drusano/main.rs

@@ -111,4 +111,3 @@ fn main() -> Result<()> {
         .fit_with(NonParametricAlgorithm::npag())?;
     Ok(())
 }
-

examples/iov/main.rs

@@ -39,4 +39,3 @@ fn main() -> Result<()> {
 
     Ok(())
 }
-

examples/iov_analysis.rs

@@ -0,0 +1,116 @@
+//! SDE Inter-Occasion Variability (IOV) analysis.
+//!
+//! Stage 1: NPAG fit with ODE to find support points.
+//! Stage 2: Optimize SDE diffusion parameter (ske) per support point.
+//!
+//! Uses small synthetic data and low particle count (50) to keep the example fast.
+
+use pharmsol::SubjectBuilderExt;
+use pmcore::iov::DiffusionConfig;
+use pmcore::prelude::*;
+
+fn main() -> Result<()> {
+    // ── Synthetic data ──
+    let subject = Subject::builder("1")
+        .bolus(0.0, 100.0, 1)
+        .observation(1.0, 10.0, 1)
+        .observation(2.0, 6.0, 1)
+        .observation(4.0, 2.5, 1)
+        .build();
+    let data = Data::new(vec![subject]);
+
+    // ── Stage 1: ODE fit ──
+
+    let ode = ode! {
+        name: "iov_stage1_ode",
+        params: [ke, v],
+        states: [central],
+        outputs: [outeq_1],
+        routes: [bolus(input_1) -> central],
+        diffeq: |x, p, _t, dx, _cov| {
+            fetch_params!(p, ke);
+            dx[0] = -ke * x[0];
+        },
+        out: |x, p, _t, _cov, y| {
+            fetch_params!(p, v);
+            y[0] = x[0] / v;
+        },
+    };
+
+    let parameters = ParameterSpace::bounded()
+        .add("ke", 0.001, 2.0)
+        .add("v", 5.0, 50.0);
+    let prior = Theta::sobol(&parameters, 50)?;
+    let error_models = AssayErrorModels::new().add(
+        "outeq_1",
+        AssayErrorModel::additive(ErrorPoly::new(0.0, 0.5, 0.0, 0.0), 0.0),
+    )?;
+
+    println!("Stage 1: Fitting ODE with NPAG...");
+    let r_ode = EstimationProblem::nonparametric(ode, data.clone(), prior, error_models.clone())?
+        .fit_with(NonParametricAlgorithm::npag())?;
+
+    let n_spp = r_ode.get_theta().nspp();
+    println!("  {} support points, OBJF = {:.2}", n_spp, r_ode.objf());
+
+    // ── Stage 2: SDE IOV ──
+
+    let sde = sde! {
+        name: "iov_stage2_sde",
+        params: [ke, v, ske],
+        states: [central, ke0],
+        outputs: [outeq_1],
+        particles: 50,
+        routes: [bolus(input_1) -> central],
+        drift: |x, p, _t, dx, _cov| {
+            fetch_params!(p, ke);
+            dx[ke0] = -x[ke0] + ke;
+            dx[central] = -x[ke0] * x[central];
+        },
+        diffusion: |p, sigma| {
+            fetch_params!(p, _ke, _v, ske);
+            sigma[ke0] = ske;
+        },
+        init: |p, _t, _cov, x| {
+            fetch_params!(p, ke);
+            x[ke0] = ke;
+        },
+        out: |x, p, _t, _cov, y| {
+            fetch_params!(p, _ke, v);
+            y[0] = x[central] / v;
+        },
+    };
+
+    let mut joint = r_ode
+        .get_theta()
+        .with_added_parameter("ske", 1e-6, 1.0, 0.01)?;
+
+    println!(
+        "Stage 2: Optimizing sigma for {} support points...",
+        joint.nspp()
+    );
+
+    let diff = sde.optimize_diffusion(
+        r_ode.data(),
+        &mut joint,
+        &["ske".to_string()],
+        r_ode.error_models(),
+        DiffusionConfig {
+            max_iter: 10,
+            ..DiffusionConfig::default()
+        },
+    )?;
+
+    let n_converged = diff.per_point_converged.iter().filter(|&&c| c).count();
+    let mean_ll: f64 =
+        diff.per_point_likelihood.iter().sum::<f64>() / diff.per_point_likelihood.len() as f64;
+
+    println!(
+        "  Converged: {}/{} points, mean log-likelihood: {:.4}",
+        n_converged,
+        diff.per_point_converged.len(),
+        mean_ll,
+    );
+
+    Ok(())
+}

examples/meta/main.rs

@@ -55,4 +55,3 @@ fn main() -> Result<()> {
 
     Ok(())
 }
-

examples/neely/main.rs

@@ -77,4 +77,3 @@ fn main() -> Result<()> {
 
     Ok(())
 }
-

examples/new_iov/main.rs

@@ -39,4 +39,3 @@ fn main() -> Result<()> {
 
     Ok(())
 }
-

examples/theophylline/main.rs

@@ -30,4 +30,3 @@ fn main() -> Result<()> {
 
     Ok(())
 }
-

examples/two_eq_lag/main.rs

@@ -37,4 +37,3 @@ fn main() -> Result<()> {
 
     Ok(())
 }
-

examples/vanco_sde/main.rs

@@ -65,4 +65,3 @@ fn main() -> Result<()> {
 
     Ok(())
 }
-

src/algorithms/nonparametric/npod.rs

@@ -85,7 +85,7 @@ impl<E: Equation + Send + 'static> NPOD<E> {
             equation,
             psi: Psi::new(),
             prior: theta.clone(),
-            theta: theta,
+            theta,
             lambda: Weights::default(),
             w: Weights::default(),
             last_objf: -1e30,

src/bestdose/types.rs

@@ -194,7 +194,10 @@ pub struct BestDoseConfig {
 }
 
 impl BestDoseConfig {
-    pub fn new(parameter_space: ParameterSpace<BoundedParameter>, error_models: AssayErrorModels) -> Self {
+    pub fn new(
+        parameter_space: ParameterSpace<BoundedParameter>,
+        error_models: AssayErrorModels,
+    ) -> Self {
         Self {
             parameter_space,
             error_models,