feat: Chain fit methods using result (#284)

* Makes it possible to directly chain different algorithms, or the same algorithms with different configurations.

Author: Siel

.gitignore

@@ -47,4 +47,6 @@ lcov.info
 Fortran/
 paper/
 docs/
-examples/**/outputs/
+examples/**/outputs/
+examples/iov_synthetic/
+examples/iov_*.rs

examples/chain_algorithms.rs

@@ -0,0 +1,61 @@
+//! 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, _t, dx| {
+            dx[central] = -ke * x[central];
+        },
+        out: |x, _t, y| {
+            y[outeq_1] = x[central] / v;
+        },
+    };
+
+    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),
+    )?;
+
+    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(())
+}

src/algorithms/nonparametric/npod.rs

@@ -95,7 +95,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/estimation/error_models.rs

@@ -19,8 +19,8 @@ impl ErrorModels {
     }
 }
 
-impl Into<AssayErrorModels> for ErrorModels {
-    fn into(self) -> AssayErrorModels {
-        self.models().clone()
+impl From<ErrorModels> for AssayErrorModels {
+    fn from(val: ErrorModels) -> Self {
+        val.models().clone()
     }
 }

src/estimation/nonparametric/result.rs

@@ -89,6 +89,36 @@ impl<E: Equation> NonParametricResult<E> {
         &self.cyclelog
     }
 
+    /// Chain this result into a new fit with a different algorithm.
+    ///
+    /// Uses the optimized theta as the prior for the new run. Keeps the same
+    /// equation, data, and error models. Consumes `self`.
+    ///
+    /// # Example
+    ///
+    /// ```ignore
+    /// let result = problem
+    ///     .fit_with(NpagConfig::new().max_cycles(100))?
+    ///     .chain(NpodConfig::new().max_cycles(50))?;
+    /// ```
+    pub fn chain<A>(self, algorithm: A) -> anyhow::Result<NonParametricResult<E>>
+    where
+        A: crate::algorithms::Algorithm<
+            E,
+            crate::estimation::NonParametric,
+            Output = NonParametricResult<E>,
+        >,
+        E: crate::model::EquationMetadataSource,
+    {
+        crate::estimation::EstimationProblem::nonparametric(
+            self.equation,
+            self.data,
+            self.theta,
+            self.error_models,
+        )?
+        .fit_with(algorithm)
+    }
+
     pub fn psi(&self) -> &Psi {
         &self.psi
     }
@@ -262,3 +292,144 @@ impl<E: Equation> NonParametricResult<E> {
         Ok(())
     }
 }
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::algorithms::nonparametric::{NpagConfig, NpodConfig};
+    use crate::estimation::EstimationProblem;
+    use crate::model::ParameterSpace;
+    use pharmsol::equation::metadata;
+    use pharmsol::prelude::data::{AssayErrorModel, AssayErrorModels};
+    use pharmsol::{ErrorPoly, SubjectBuilderExt};
+
+    fn minimal_ode() -> pharmsol::ODE {
+        pharmsol::equation::ODE::new(
+            |x, p, _t, dx, b, _rateiv, _cov| {
+                let ke = p[0];
+                dx[0] = -ke * x[0] + b[0];
+            },
+            |_p, _t, _cov| pharmsol::lag! {},
+            |_p, _t, _cov| pharmsol::fa! {},
+            |_p, _t, _cov, _x| {},
+            |x, p, _t, _cov, y| {
+                let v = p[1];
+                y[0] = x[0] / v;
+            },
+        )
+        .with_nstates(1)
+        .with_ndrugs(1)
+        .with_nout(1)
+        .with_metadata(
+            metadata::new("chain_test")
+                .parameters(["ke", "v"])
+                .states(["central"])
+                .outputs(["0"])
+                .route(pharmsol::equation::Route::bolus("0").to_state("central")),
+        )
+        .expect("metadata should validate")
+    }
+
+    fn minimal_data() -> pharmsol::Data {
+        let subject = pharmsol::Subject::builder("1")
+            .bolus(0.0, 100.0, 0)
+            .observation(1.0, 10.0, 0)
+            .observation(2.0, 8.0, 0)
+            .build();
+        pharmsol::Data::new(vec![subject])
+    }
+
+    #[test]
+    fn chain_npag_to_npod_preserves_support_points() {
+        let ode = minimal_ode();
+        let data = minimal_data();
+        let params = ParameterSpace::bounded()
+            .add("ke", 0.001, 3.0)
+            .add("v", 25.0, 250.0);
+        let prior = Theta::sobol_default(&params).unwrap();
+        let err = AssayErrorModels::new()
+            .add(
+                "0",
+                AssayErrorModel::additive(ErrorPoly::new(0.0, 0.5, 0.0, 0.0), 0.0),
+            )
+            .unwrap();
+
+        let r1 = EstimationProblem::nonparametric(ode, data, prior, err)
+            .unwrap()
+            .fit_with(NpagConfig::new().max_cycles(2))
+            .unwrap();
+
+        let n_spp = r1.get_theta().nspp();
+        assert!(n_spp > 0, "NPAG should produce support points");
+
+        // Chain into NPOD — should complete without error
+        let r2 = r1.chain(NpodConfig::new().max_cycles(1)).unwrap();
+
+        assert!(r2.get_theta().nspp() > 0);
+        assert_eq!(r2.data().subjects().len(), 1);
+        assert_eq!(r2.cycles(), 1);
+    }
+
+    #[test]
+    fn chain_npag_to_npag_maintains_or_improves_objf() {
+        let ode = minimal_ode();
+        let data = minimal_data();
+        let params = ParameterSpace::bounded()
+            .add("ke", 0.001, 3.0)
+            .add("v", 25.0, 250.0);
+        let prior = Theta::sobol_with_seed(&params, 5, 42).unwrap();
+        let err = AssayErrorModels::new()
+            .add(
+                "0",
+                AssayErrorModel::additive(ErrorPoly::new(0.0, 0.5, 0.0, 0.0), 0.0),
+            )
+            .unwrap();
+
+        let r1 = EstimationProblem::nonparametric(ode, data, prior, err)
+            .unwrap()
+            .fit_with(NpagConfig::new().max_cycles(5))
+            .unwrap();
+
+        let objf1 = r1.objf();
+
+        // Chain into another NPAG run
+        let r2 = r1.chain(NpagConfig::new().max_cycles(3)).unwrap();
+
+        // Second run should not regress significantly
+        assert!(
+            r2.objf() <= objf1 + 0.5,
+            "OBJF regressed: {} -> {}",
+            objf1,
+            r2.objf()
+        );
+    }
+
+    #[test]
+    fn chain_npag_to_npod_with_unconverged_result() {
+        let ode = minimal_ode();
+        let data = minimal_data();
+        let params = ParameterSpace::bounded()
+            .add("ke", 0.001, 3.0)
+            .add("v", 25.0, 250.0);
+        let prior = Theta::sobol_with_seed(&params, 5, 42).unwrap();
+        let err = AssayErrorModels::new()
+            .add(
+                "0",
+                AssayErrorModel::additive(ErrorPoly::new(0.0, 0.5, 0.0, 0.0), 0.0),
+            )
+            .unwrap();
+
+        // Run only 1 cycle — will be unconverged
+        let r1 = EstimationProblem::nonparametric(ode, data, prior, err)
+            .unwrap()
+            .fit_with(NpagConfig::new().max_cycles(1))
+            .unwrap();
+
+        // Chaining from unconverged should still work
+        let r2 = r1.chain(NpodConfig::new().max_cycles(1));
+        assert!(
+            r2.is_ok(),
+            "Chaining from unconverged result should succeed"
+        );
+    }
+}

src/estimation/nonparametric/theta.rs

@@ -144,6 +144,61 @@ impl Theta {
         true
     }
 
+    /// Create a new Theta with an additional parameter column.
+    ///
+    /// All existing rows keep their values. The new column is filled with
+    /// `initial_value` for every row. Returns a new `Theta` — does not
+    /// mutate in place.
+    ///
+    /// # Errors
+    ///
+    /// - `name` already exists in the parameter space
+    /// - `lower` or `upper` are non-finite
+    /// - `lower >= upper`
+    pub fn with_added_parameter(
+        &self,
+        name: &str,
+        lower: f64,
+        upper: f64,
+        initial_value: f64,
+    ) -> Result<Theta> {
+        // Validate uniqueness
+        if self.parameters().iter().any(|p| p.name.as_str() == name) {
+            bail!("parameter '{}' already exists in theta", name);
+        }
+
+        // Validate bounds
+        if !lower.is_finite() || !upper.is_finite() {
+            bail!(
+                "bounds must be finite for parameter '{}': [{}, {}]",
+                name,
+                lower,
+                upper
+            );
+        }
+        if lower >= upper {
+            bail!(
+                "lower bound ({}) must be strictly less than upper bound ({}) for parameter '{}'",
+                lower,
+                upper,
+                name
+            );
+        }
+
+        let (nrows, ncols) = (self.matrix().nrows(), self.matrix().ncols());
+        let new_matrix = faer::Mat::from_fn(nrows, ncols + 1, |r, c| {
+            if c < ncols {
+                self.matrix()[(r, c)]
+            } else {
+                initial_value
+            }
+        });
+
+        let new_params = self.parameters().clone().add(name, lower, upper);
+
+        Theta::from_parts(new_matrix, new_params)
+    }
+
     /// Write the matrix to a CSV file
     pub fn write(&self, path: &str) {
         let mut writer = csv::Writer::from_path(path).unwrap();