diff --git a/benches/bimodal_ke.rs b/benches/bimodal_ke.rs
index 25ba91ea8..b6f12d869 100644
--- a/benches/bimodal_ke.rs
+++ b/benches/bimodal_ke.rs
@@ -1,6 +1,6 @@
 use anyhow::Result;
 use criterion::{criterion_group, criterion_main, Criterion};
-use pmcore::prelude::*;
+use pmcore::{algorithms::Algorithm, prelude::*};
 
 use std::hint::black_box;
 
@@ -22,74 +22,87 @@ fn create_equation() -> equation::ODE {
     }
 }
 
-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_npag() -> Result<EstimationProblem<equation::ODE>> {
+fn setup_npag() -> Result<EstimationProblem<equation::ODE, NonParametric>> {
     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(Npag::new())
-        .error("outeq_1", create_error_model())?
-        .progress(false)
-        .prior(Prior::sobol(2048, 22))
-        .build()
+    let parameters = ParameterSpace::bounded()
+        .add("ke", 0.001, 3.0)
+        .add("v", 25.0, 250.0);
+    let prior = Theta::sobol_default(&parameters)?;
+    let error_models = AssayErrorModels::new().add(
+        "outeq_1",
+        AssayErrorModel::additive(ErrorPoly::new(0.0, 0.5, 0.0, 0.0), 0.0),
+    )?;
+    EstimationProblem::nonparametric(create_equation(), data, prior, error_models)
 }
 
-fn setup_npod() -> Result<EstimationProblem<equation::ODE>> {
+fn setup_npod() -> Result<EstimationProblem<equation::ODE, NonParametric>> {
     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())?
-        .progress(false)
-        .prior(Prior::sobol(2048, 22))
-        .build()
+    let parameters = ParameterSpace::bounded()
+        .add("ke", 0.001, 3.0)
+        .add("v", 25.0, 250.0);
+    let prior = Theta::sobol_default(&parameters)?;
+    let error_models = AssayErrorModels::new().add(
+        "outeq_1",
+        AssayErrorModel::additive(ErrorPoly::new(0.0, 0.5, 0.0, 0.0), 0.0),
+    )?;
+    EstimationProblem::nonparametric(create_equation(), data, prior, error_models)
 }
 
-fn setup_postprob() -> Result<EstimationProblem<equation::ODE>> {
+fn setup_postprob() -> Result<EstimationProblem<equation::ODE, NonParametric>> {
     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())?
-        .progress(false)
-        .prior(Prior::sobol(2048, 22))
-        .build()
+    let parameters = ParameterSpace::bounded()
+        .add("ke", 0.001, 3.0)
+        .add("v", 25.0, 250.0);
+    let prior = Theta::sobol_default(&parameters)?;
+    let error_models = AssayErrorModels::new().add(
+        "outeq_1",
+        AssayErrorModel::additive(ErrorPoly::new(0.0, 0.5, 0.0, 0.0), 0.0),
+    )?;
+    EstimationProblem::nonparametric(create_equation(), data, prior, error_models)
 }
 
-fn benchmark_algorithm<F>(c: &mut Criterion, bench_name: &str, setup_fn: F)
+fn benchmark_algorithm<F, A>(c: &mut Criterion, bench_name: &str, setup_fn: F, config: A)
 where
-    F: Fn() -> Result<EstimationProblem<equation::ODE>>,
+    F: Fn() -> Result<EstimationProblem<equation::ODE, NonParametric>>,
+    A: Algorithm<equation::ODE, NonParametric> + Clone,
 {
-    let problem = setup_fn().unwrap();
-
     c.bench_function(bench_name, |b| {
         b.iter_with_setup(
-            || problem.clone(),
-            |problem| black_box(problem.fit().unwrap()),
+            || setup_fn().unwrap(),
+            |problem| black_box(problem.fit_with(config.clone()).unwrap()),
         )
     });
 }
 
 fn benchmark_bimodal_ke_npag(c: &mut Criterion) {
-    benchmark_algorithm(c, "bimodal_ke_npag", setup_npag);
+    benchmark_algorithm(
+        c,
+        "bimodal_ke_npag",
+        setup_npag,
+        NonParametricAlgorithm::npag(),
+    );
 }
 
 fn benchmark_bimodal_ke_npod(c: &mut Criterion) {
-    benchmark_algorithm(c, "bimodal_ke_npod", setup_npod);
+    benchmark_algorithm(
+        c,
+        "bimodal_ke_npod",
+        setup_npod,
+        NonParametricAlgorithm::npod(),
+    );
 }
 
 fn benchmark_bimodal_ke_postprob(c: &mut Criterion) {
-    benchmark_algorithm(c, "bimodal_ke_postprob", setup_postprob);
+    benchmark_algorithm(
+        c,
+        "bimodal_ke_postprob",
+        setup_postprob,
+        NonParametricAlgorithm::npmap(),
+    );
 }
 
 criterion_group! {
@@ -98,3 +111,4 @@ criterion_group! {
     targets = benchmark_bimodal_ke_npag, benchmark_bimodal_ke_npod, benchmark_bimodal_ke_postprob
 }
 criterion_main!(benches);
+
diff --git a/examples/bestdose.rs b/examples/bestdose.rs
index ebbd341ff..9fd9aed2e 100644
--- a/examples/bestdose.rs
+++ b/examples/bestdose.rs
@@ -20,9 +20,9 @@ fn main() -> Result<()> {
         },
     };
 
-    let parameter_space = ParameterSpace::new()
-        .add(Parameter::bounded("ke", 0.001, 3.0))
-        .add(Parameter::bounded("v", 25.0, 250.0));
+    let parameter_space = ParameterSpace::<BoundedParameter>::new()
+        .add("ke", 0.001, 3.0)
+        .add("v", 25.0, 250.0);
 
     let ems = AssayErrorModels::new().add(
         0,
@@ -65,7 +65,8 @@ fn main() -> Result<()> {
         .observation(18.0, conc(6.0, 75.0) + conc(18.0, 150.0), 0)
         .build();
 
-    let (theta, prior) = read_prior("examples/bimodal_ke/output/theta.csv", &parameter_space)?;
+    let (theta, prior) =
+        Theta::from_file("examples/bimodal_ke/output/theta.csv", &parameter_space)?;
 
     let posterior = BestDosePosterior::compute(
         &theta,
diff --git a/examples/bestdose_auc.rs b/examples/bestdose_auc.rs
index df179004f..0bcacc0e4 100644
--- a/examples/bestdose_auc.rs
+++ b/examples/bestdose_auc.rs
@@ -28,9 +28,9 @@ fn main() -> Result<()> {
     };
 
     // Minimal parameter ranges
-    let parameter_space = ParameterSpace::new()
-        .add(Parameter::bounded("ke", 0.001, 3.0))
-        .add(Parameter::bounded("v", 25.0, 250.0));
+    let parameter_space = ParameterSpace::<BoundedParameter>::new()
+        .add("ke", 0.001, 3.0)
+        .add("v", 25.0, 250.0);
 
     let ems = AssayErrorModels::new().add(
         0,
@@ -43,7 +43,10 @@ fn main() -> Result<()> {
 
     // Load realistic prior from previous NPAG run (47 support points)
     println!("Loading prior from bimodal_ke example...");
-    let (theta, prior) = read_prior("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());
@@ -118,9 +121,8 @@ fn main() -> Result<()> {
         }
     }
 
-    // =========================================================================
     // EXAMPLE 2: Interval AUC (AUCFromLastDose)
-    // =========================================================================
+
     println!("\n\n");
     println!("════════════════════════════════════════════════════════");
     println!("  EXAMPLE 2: Interval AUC (AUCFromLastDose)");
diff --git a/examples/bestdose_bounds.rs b/examples/bestdose_bounds.rs
index 91b668b86..286251618 100644
--- a/examples/bestdose_bounds.rs
+++ b/examples/bestdose_bounds.rs
@@ -27,9 +27,9 @@ fn main() -> Result<()> {
         },
     };
 
-    let parameter_space = ParameterSpace::new()
-        .add(Parameter::bounded("ke", 0.001, 3.0))
-        .add(Parameter::bounded("v", 25.0, 250.0));
+    let parameter_space = ParameterSpace::<BoundedParameter>::new()
+        .add("ke", 0.001, 3.0)
+        .add("v", 25.0, 250.0);
 
     let ems = AssayErrorModels::new().add(
         0,
@@ -40,7 +40,10 @@ fn main() -> Result<()> {
 
     // Load realistic prior from previous NPAG run
     println!("Loading prior from bimodal_ke example...");
-    let (theta, prior) = read_prior("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());
diff --git a/examples/bestdose_cov.rs b/examples/bestdose_cov.rs
index 07065d818..1cb6ac3e2 100644
--- a/examples/bestdose_cov.rs
+++ b/examples/bestdose_cov.rs
@@ -1,7 +1,5 @@
 use anyhow::Result;
-use pmcore::bestdose; // bestdose new
-                      // use pmcore::bestdose::bestdose_old as bestdose; // bestdose old
-
+use pmcore::bestdose::{BestDoseConfig, BestDosePosterior, DoseRange, Target};
 use pmcore::prelude::*;
 
 fn main() -> Result<()> {
@@ -23,17 +21,16 @@ fn main() -> Result<()> {
         },
     };
 
-    let parameter_space = ParameterSpace::new()
-        .add(Parameter::bounded("ke", 0.001, 3.0))
-        .add(Parameter::bounded("v", 25.0 / 70.0, 250.0 / 70.0));
+    let parameter_space = ParameterSpace::<BoundedParameter>::new()
+        .add("ke", 0.001, 3.0)
+        .add("v", 25.0 / 70.0, 250.0 / 70.0);
 
     let ems = AssayErrorModels::new().add(
         0,
         AssayErrorModel::additive(ErrorPoly::new(0.0, 0.20, 0.0, 0.0), 0.0),
     )?;
 
-    let config =
-        bestdose::BestDoseConfig::new(parameter_space.clone(), ems.clone()).with_progress(false);
+    let config = BestDoseConfig::new(parameter_space.clone(), ems.clone()).with_progress(false);
 
     // Generate a patient with known parameters
     // Ke = 0.5, V = 50
@@ -93,29 +90,22 @@ fn main() -> Result<()> {
         .observation(18.0, conc(6.0, 75.0) + conc(18.0, 150.0), 0)
         .build();
 
-    let (mut theta, prior) = read_prior("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() {
         m_t[(i, 1)] /= 70.0;
     }
 
-    // Example usage - using new() constructor which calculates NPAGFULL11 posterior
-    // max_cycles controls NPAGFULL refinement:
-    //   0 = NPAGFULL11 only (fast but less accurate)
-    //   100 = moderate refinement
-    //   500 = full refinement (Fortran default, slow but most accurate)
-    let problem = bestdose::BestDoseProblem::new(
+    let posterior = BestDosePosterior::compute(
         &theta,
         &prior.unwrap(),
         Some(past_data.clone()), // Optional: past data for Bayesian updating
-        target_data.clone(),
-        None,
         eq.clone(),
-        bestdose::DoseRange::new(0.0, 300.0),
-        0.0,
         config.clone(),
-        bestdose::Target::Concentration, // Target concentrations (not AUCs)
     )?;
 
     println!("Optimizing dose...");
@@ -125,7 +115,13 @@ fn main() -> Result<()> {
 
     for bias_weight in &bias_weights {
         println!("Running optimization with bias weight: {}", bias_weight);
-        let optimal = problem.clone().with_bias_weight(*bias_weight).optimize()?;
+        let optimal = posterior.optimize(
+            target_data.clone(),
+            None,
+            DoseRange::new(0.0, 300.0),
+            *bias_weight,
+            Target::Concentration,
+        )?;
         results.push((bias_weight, optimal));
     }
 
diff --git a/examples/bimodal_ke/main.rs b/examples/bimodal_ke/main.rs
index 447e8cc17..bc6fe7604 100644
--- a/examples/bimodal_ke/main.rs
+++ b/examples/bimodal_ke/main.rs
@@ -1,8 +1,8 @@
 use anyhow::Result;
-use pmcore::{output::logging::Logger, prelude::*};
+use pmcore::prelude::*;
 
 fn main() -> Result<()> {
-    Logger::new().init()?;
+    Logger::new().stdout(true).init()?;
 
     let eq = ode! {
         name: "bimodal_ke",
@@ -23,15 +23,20 @@ fn main() -> Result<()> {
 
     let data = data::read_pmetrics("examples/bimodal_ke/bimodal_ke.csv")?;
 
-    let _result = EstimationProblem::builder(eq, data)
-        .method(Npag::default())
-        .parameter(Parameter::bounded("ke", 0.001, 3.0))?
-        .parameter(Parameter::bounded("v", 25.0, 250.0))?
-        .error(
-            "outeq_1",
-            AssayErrorModel::additive(ErrorPoly::new(0.0, 0.5, 0.0, 0.0), 0.0),
-        )?
-        .fit()?;
+    let parameters = ParameterSpace::bounded()
+        .add("ke", 0.001, 3.0)
+        .add("v", 25.0, 250.0);
+
+    let prior = Theta::sobol_default(&parameters)?;
+
+    let error_models = AssayErrorModels::new().add(
+        "outeq_1",
+        AssayErrorModel::additive(ErrorPoly::new(0.0, 0.5, 0.0, 0.0), 0.0),
+    )?;
+
+    let problem = EstimationProblem::nonparametric(eq, data, prior, error_models)?;
+
+    let _result = problem.fit_with(NonParametricAlgorithm::npag())?;
 
     Ok(())
 }
diff --git a/examples/bimodal_ke_backend_compare.rs b/examples/bimodal_ke_backend_compare.rs
index 52e99ad13..5aa08b8bc 100644
--- a/examples/bimodal_ke_backend_compare.rs
+++ b/examples/bimodal_ke_backend_compare.rs
@@ -202,17 +202,16 @@ fn run_case<E: pharmsol::Equation + Clone + Send + 'static + EquationMetadataSou
 ) -> Result<ComparisonResult> {
     let data = data::read_pmetrics(DATA_PATH)?;
     let fit_started = Instant::now();
-    let result = EstimationProblem::builder(equation, data)
-        .parameter(Parameter::bounded("ke", 0.001, 3.0))?
-        .parameter(Parameter::bounded("v", 25.0, 250.0))?
-        .method(Npag::new())
-        .error(
-            "outeq_1",
-            AssayErrorModel::additive(ErrorPoly::new(0.0, 0.5, 0.0, 0.0), 0.0),
-        )?
-        .cache(true)
-        .progress(false)
-        .fit()?;
+    let parameters = ParameterSpace::bounded()
+        .add("ke", 0.001, 3.0)
+        .add("v", 25.0, 250.0);
+    let prior = Theta::sobol_default(&parameters)?;
+    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(equation, data, prior, error_models)?
+        .fit_with(NonParametricAlgorithm::npag())?;
     let fit_time = fit_started.elapsed();
 
     summarize_result(label, compile_time, fit_time, &result)
@@ -222,12 +221,8 @@ fn summarize_result<E: pharmsol::Equation>(
     label: &'static str,
     compile_time: Duration,
     fit_time: Duration,
-    result: &FitResult<E>,
+    result: &NonParametricResult<E>,
 ) -> Result<ComparisonResult> {
-    result
-        .as_nonparametric()
-        .ok_or_else(|| anyhow!("expected nonparametric result for {label}"))?;
-
     Ok(ComparisonResult {
         label,
         compile_time,
@@ -318,3 +313,4 @@ fn print_summary(results: &[ComparisonResult]) -> Result<()> {
 
     Ok(())
 }
+
diff --git a/examples/drusano/main.rs b/examples/drusano/main.rs
index 56cc0c590..240b69656 100644
--- a/examples/drusano/main.rs
+++ b/examples/drusano/main.rs
@@ -64,49 +64,51 @@ fn main() -> Result<()> {
     };
 
     let data = data::read_pmetrics("examples/drusano/data.csv")?;
-    EstimationProblem::builder(eq, data)
-        .parameter(Parameter::bounded("v1", 5.0, 160.0))?
-        .parameter(Parameter::bounded("cl1", 4.0, 9.0))?
-        .parameter(Parameter::bounded("v2", 100.0, 200.0))?
-        .parameter(Parameter::bounded("cl2", 25.0, 35.0))?
-        .parameter(Parameter::bounded("popmax", 100000000.0, 100000000000.0))?
-        .parameter(Parameter::bounded("kgs", 0.01, 0.25))?
-        .parameter(Parameter::bounded("kks", 0.01, 0.5))?
-        .parameter(Parameter::bounded("e50_1s", 0.1, 2.5))?
-        .parameter(Parameter::bounded("e50_2s", 0.1, 10.0))?
-        .parameter(Parameter::bounded("alpha_s", -8.0, 5.0))?
-        .parameter(Parameter::bounded("kgr1", 0.004, 0.1))?
-        .parameter(Parameter::bounded("kkr1", 0.08, 0.4))?
-        .parameter(Parameter::bounded("e50_1r1", 8.0, 17.0))?
-        .parameter(Parameter::bounded("alpha_r1", -8.0, 5.0))?
-        .parameter(Parameter::bounded("kgr2", 0.004, 0.3))?
-        .parameter(Parameter::bounded("kkr2", 0.1, 0.5))?
-        .parameter(Parameter::bounded("e50_2r2", 5.0, 8.0))?
-        .parameter(Parameter::bounded("alpha_r2", -5.0, 5.0))?
-        .parameter(Parameter::bounded("init_4", -1.0, 4.0))?
-        .parameter(Parameter::bounded("init_5", -1.0, 3.0))?
-        .parameter(Parameter::bounded("h1s", 0.5, 8.0))?
-        .parameter(Parameter::bounded("h2s", 0.1, 4.0))?
-        .parameter(Parameter::bounded("h1r1", 5.0, 25.0))?
-        .parameter(Parameter::bounded("h2r2", 10.0, 22.0))?
-        .method(Npag::new())
-        .error(
+    let parameters = ParameterSpace::bounded()
+        .add("v1", 5.0, 160.0)
+        .add("cl1", 4.0, 9.0)
+        .add("v2", 100.0, 200.0)
+        .add("cl2", 25.0, 35.0)
+        .add("popmax", 100000000.0, 100000000000.0)
+        .add("kgs", 0.01, 0.25)
+        .add("kks", 0.01, 0.5)
+        .add("e50_1s", 0.1, 2.5)
+        .add("e50_2s", 0.1, 10.0)
+        .add("alpha_s", -8.0, 5.0)
+        .add("kgr1", 0.004, 0.1)
+        .add("kkr1", 0.08, 0.4)
+        .add("e50_1r1", 8.0, 17.0)
+        .add("alpha_r1", -8.0, 5.0)
+        .add("kgr2", 0.004, 0.3)
+        .add("kkr2", 0.1, 0.5)
+        .add("e50_2r2", 5.0, 8.0)
+        .add("alpha_r2", -5.0, 5.0)
+        .add("init_4", -1.0, 4.0)
+        .add("init_5", -1.0, 3.0)
+        .add("h1s", 0.5, 8.0)
+        .add("h2s", 0.1, 4.0)
+        .add("h1r1", 5.0, 25.0)
+        .add("h2r2", 10.0, 22.0);
+    let prior = Theta::sobol_default(&parameters)?;
+    let error_models = AssayErrorModels::new()
+        .add(
             "outeq_1",
             AssayErrorModel::proportional(ErrorPoly::new(0.1, 0.1, 0.0, 0.0), 1.0),
         )?
-        .error(
+        .add(
             "outeq_2",
             AssayErrorModel::proportional(ErrorPoly::new(0.1, 0.1, 0.0, 0.0), 1.0),
         )?
-        .error(
+        .add(
             "outeq_3",
             AssayErrorModel::proportional(ErrorPoly::new(0.1, 0.1, 0.0, 0.0), 1.0),
         )?
-        .error(
+        .add(
             "outeq_4",
             AssayErrorModel::proportional(ErrorPoly::new(0.1, 0.1, 0.0, 0.0), 1.0),
-        )?
-        .prior(Prior::sobol(212900, 347))
-        .fit()?;
+        )?;
+    EstimationProblem::nonparametric(eq, data, prior, error_models)?
+        .fit_with(NonParametricAlgorithm::npag())?;
     Ok(())
 }
+
diff --git a/examples/iov/main.rs b/examples/iov/main.rs
index 5d95ba56a..7e0f5ba4c 100644
--- a/examples/iov/main.rs
+++ b/examples/iov/main.rs
@@ -27,16 +27,16 @@ fn main() -> Result<()> {
     };
 
     let data = data::read_pmetrics("examples/iov/test.csv").unwrap();
-    EstimationProblem::builder(sde, data)
-        .parameter(Parameter::bounded("ke0", 0.001, 2.0))?
-        .method(Npag::new())
-        .error(
-            "outeq_1",
-            AssayErrorModel::additive(ErrorPoly::new(0.0, 0.0, 0.0, 0.0), 0.0000757575757576),
-        )?
-        .prior(Prior::sobol(100, 347))
-        .fit()
+    let parameters = ParameterSpace::bounded().add("ke0", 0.001, 2.0);
+    let prior = Theta::sobol_default(&parameters)?;
+    let error_models = AssayErrorModels::new().add(
+        "outeq_1",
+        AssayErrorModel::additive(ErrorPoly::new(0.0, 0.0, 0.0, 0.0), 0.0000757575757576),
+    )?;
+    EstimationProblem::nonparametric(sde, data, prior, error_models)?
+        .fit_with(NonParametricAlgorithm::npag())
         .unwrap();
 
     Ok(())
 }
+
diff --git a/examples/meta/main.rs b/examples/meta/main.rs
index 644cba6a6..ea89bdfd4 100644
--- a/examples/meta/main.rs
+++ b/examples/meta/main.rs
@@ -4,7 +4,7 @@
 
 use pmcore::prelude::*;
 
-fn main() {
+fn main() -> Result<()> {
     let eq = ode! {
         name: "meta",
         params: [cls, fm, k20, relv, theta1, theta2, vs],
@@ -31,33 +31,28 @@ fn main() {
         },
     };
 
-    let data = data::read_pmetrics("examples/meta/meta.csv").unwrap();
-    EstimationProblem::builder(eq, data)
-        .parameter(Parameter::bounded("cls", 0.1, 10.0))
-        .unwrap()
-        .parameter(Parameter::bounded("fm", 0.0, 1.0))
-        .unwrap()
-        .parameter(Parameter::bounded("k20", 0.01, 1.0))
-        .unwrap()
-        .parameter(Parameter::bounded("relv", 0.1, 1.0))
-        .unwrap()
-        .parameter(Parameter::bounded("theta1", 0.1, 10.0))
-        .unwrap()
-        .parameter(Parameter::bounded("theta2", 0.1, 10.0))
-        .unwrap()
-        .parameter(Parameter::bounded("vs", 1.0, 10.0))
-        .unwrap()
-        .method(Npod::new())
-        .error(
+    let data = data::read_pmetrics("examples/meta/meta.csv")?;
+    let parameters = ParameterSpace::bounded()
+        .add("cls", 0.1, 10.0)
+        .add("fm", 0.0, 1.0)
+        .add("k20", 0.01, 1.0)
+        .add("relv", 0.1, 1.0)
+        .add("theta1", 0.1, 10.0)
+        .add("theta2", 0.1, 10.0)
+        .add("vs", 1.0, 10.0);
+    let prior = Theta::sobol_default(&parameters)?;
+    let error_models = AssayErrorModels::new()
+        .add(
             "outeq_1",
             AssayErrorModel::proportional(ErrorPoly::new(1.0, 0.1, 0.0, 0.0), 5.0),
-        )
-        .unwrap()
-        .error(
+        )?
+        .add(
             "outeq_2",
             AssayErrorModel::proportional(ErrorPoly::new(1.0, 0.1, 0.0, 0.0), 5.0),
-        )
-        .unwrap()
-        .fit()
-        .unwrap();
+        )?;
+    EstimationProblem::nonparametric(eq, data, prior, error_models)?
+        .fit_with(NonParametricAlgorithm::npod())?;
+
+    Ok(())
 }
+
diff --git a/examples/neely/main.rs b/examples/neely/main.rs
index 5f5e7d7fa..922547721 100644
--- a/examples/neely/main.rs
+++ b/examples/neely/main.rs
@@ -1,6 +1,6 @@
 use pmcore::prelude::*;
 
-fn main() {
+fn main() -> Result<()> {
     let eq = ode! {
         name: "neely",
         params: [cls, k30, k40, qs, vps, vs, fm1, fm2, theta1, theta2],
@@ -46,45 +46,35 @@ fn main() {
         },
     };
 
-    let data = data::read_pmetrics("examples/neely/data.csv").unwrap();
-    EstimationProblem::builder(eq, data)
-        .parameter(Parameter::bounded("cls", 0.0, 0.4))
-        .unwrap()
-        .parameter(Parameter::bounded("k30", 0.0, 0.5))
-        .unwrap()
-        .parameter(Parameter::bounded("k40", 0.3, 1.5))
-        .unwrap()
-        .parameter(Parameter::bounded("qs", 0.0, 0.5))
-        .unwrap()
-        .parameter(Parameter::bounded("vps", 0.0, 5.0))
-        .unwrap()
-        .parameter(Parameter::bounded("vs", 0.0, 2.0))
-        .unwrap()
-        .parameter(Parameter::bounded("fm1", 0.0, 0.2))
-        .unwrap()
-        .parameter(Parameter::bounded("fm2", 0.0, 0.1))
-        .unwrap()
-        .parameter(Parameter::bounded("theta1", -4.0, 2.0))
-        .unwrap()
-        .parameter(Parameter::bounded("theta2", -2.0, 0.5))
-        .unwrap()
-        .method(Npag::new())
-        .error(
+    let data = data::read_pmetrics("examples/neely/data.csv")?;
+    let parameters = ParameterSpace::bounded()
+        .add("cls", 0.0, 0.4)
+        .add("k30", 0.0, 0.5)
+        .add("k40", 0.3, 1.5)
+        .add("qs", 0.0, 0.5)
+        .add("vps", 0.0, 5.0)
+        .add("vs", 0.0, 2.0)
+        .add("fm1", 0.0, 0.2)
+        .add("fm2", 0.0, 0.1)
+        .add("theta1", -4.0, 2.0)
+        .add("theta2", -2.0, 0.5);
+    let prior = Theta::sobol_default(&parameters)?;
+    let error_models = AssayErrorModels::new()
+        .add(
             "outeq_1",
             AssayErrorModel::proportional(ErrorPoly::new(1.0, 0.1, 0.0, 0.0), 5.0),
-        )
-        .unwrap()
-        .error(
+        )?
+        .add(
             "outeq_2",
             AssayErrorModel::proportional(ErrorPoly::new(1.0, 0.1, 0.0, 0.0), 5.0),
-        )
-        .unwrap()
-        .error(
+        )?
+        .add(
             "outeq_3",
             AssayErrorModel::proportional(ErrorPoly::new(1.0, 0.1, 0.0, 0.0), 5.0),
-        )
-        .unwrap()
-        .prior(Prior::sobol(2028, 22))
-        .fit()
-        .unwrap();
+        )?;
+    EstimationProblem::nonparametric(eq, data, prior, error_models)?
+        .fit_with(NonParametricAlgorithm::npag())?;
+
+    Ok(())
 }
+
diff --git a/examples/new_iov/main.rs b/examples/new_iov/main.rs
index 550899020..9d1ea6eb6 100644
--- a/examples/new_iov/main.rs
+++ b/examples/new_iov/main.rs
@@ -1,6 +1,6 @@
 use pmcore::prelude::*;
 
-fn main() {
+fn main() -> Result<()> {
     let sde = sde! {
         name: "new_iov",
         params: [ke0, ske],
@@ -25,19 +25,18 @@ fn main() {
         },
     };
 
-    let data = data::read_pmetrics("examples/new_iov/data.csv").unwrap();
-    EstimationProblem::builder(sde, data)
-        .parameter(Parameter::bounded("ke0", 0.0001, 2.4))
-        .unwrap()
-        .parameter(Parameter::bounded("ske", 0.0001, 0.2))
-        .unwrap()
-        .method(Npag::new())
-        .error(
-            "outeq_1",
-            AssayErrorModel::additive(ErrorPoly::new(-0.00119, 0.44379, -0.45864, 0.16537), 0.0),
-        )
-        .unwrap()
-        .prior(Prior::sobol(100, 347))
-        .fit()
-        .unwrap();
+    let data = data::read_pmetrics("examples/new_iov/data.csv")?;
+    let parameters = ParameterSpace::bounded()
+        .add("ke0", 0.0001, 2.4)
+        .add("ske", 0.0001, 0.2);
+    let prior = Theta::sobol_default(&parameters)?;
+    let error_models = AssayErrorModels::new().add(
+        "outeq_1",
+        AssayErrorModel::additive(ErrorPoly::new(-0.00119, 0.44379, -0.45864, 0.16537), 0.0),
+    )?;
+    EstimationProblem::nonparametric(sde, data, prior, error_models)?
+        .fit_with(NonParametricAlgorithm::npag())?;
+
+    Ok(())
 }
+
diff --git a/examples/theophylline/main.rs b/examples/theophylline/main.rs
index 73f994bc2..1d48aaa97 100644
--- a/examples/theophylline/main.rs
+++ b/examples/theophylline/main.rs
@@ -1,6 +1,6 @@
 use pmcore::prelude::*;
 
-fn main() {
+fn main() -> Result<()> {
     let analytical = analytical! {
         name: "theophylline",
         params: [ka, ke, v],
@@ -15,20 +15,19 @@ fn main() {
         },
     };
 
-    let data = data::read_pmetrics("examples/theophylline/theophylline.csv").unwrap();
-    EstimationProblem::builder(analytical, data)
-        .parameter(Parameter::bounded("ka", 0.001, 3.0))
-        .unwrap()
-        .parameter(Parameter::bounded("ke", 0.001, 3.0))
-        .unwrap()
-        .parameter(Parameter::bounded("v", 0.001, 50.0))
-        .unwrap()
-        .method(Npag::new())
-        .error(
-            "outeq_0",
-            AssayErrorModel::proportional(ErrorPoly::new(0.1, 0.1, 0.0, 0.0), 2.0),
-        )
-        .unwrap()
-        .fit()
-        .unwrap();
+    let data = data::read_pmetrics("examples/theophylline/theophylline.csv")?;
+    let parameters = ParameterSpace::bounded()
+        .add("ka", 0.001, 3.0)
+        .add("ke", 0.001, 3.0)
+        .add("v", 0.001, 50.0);
+    let prior = Theta::sobol_default(&parameters)?;
+    let error_models = AssayErrorModels::new().add(
+        "outeq_0",
+        AssayErrorModel::proportional(ErrorPoly::new(0.1, 0.1, 0.0, 0.0), 2.0),
+    )?;
+    EstimationProblem::nonparametric(analytical, data, prior, error_models)?
+        .fit_with(NonParametricAlgorithm::npag())?;
+
+    Ok(())
 }
+
diff --git a/examples/two_eq_lag/main.rs b/examples/two_eq_lag/main.rs
index e137d6a70..c6185517d 100644
--- a/examples/two_eq_lag/main.rs
+++ b/examples/two_eq_lag/main.rs
@@ -1,6 +1,6 @@
 use pmcore::prelude::*;
 
-fn main() {
+fn main() -> Result<()> {
     let eq = ode! {
         name: "two_eq_lag",
         params: [ka, ke, tlag, v],
@@ -21,22 +21,20 @@ fn main() {
         },
     };
 
-    let data = data::read_pmetrics("examples/two_eq_lag/two_eq_lag.csv").unwrap();
-    EstimationProblem::builder(eq, data)
-        .parameter(Parameter::bounded("ka", 0.1, 0.9))
-        .unwrap()
-        .parameter(Parameter::bounded("ke", 0.001, 0.1))
-        .unwrap()
-        .parameter(Parameter::bounded("tlag", 0.0, 4.0))
-        .unwrap()
-        .parameter(Parameter::bounded("v", 30.0, 120.0))
-        .unwrap()
-        .method(Npag::new())
-        .error(
-            "outeq_0",
-            AssayErrorModel::additive(ErrorPoly::new(-0.00119, 0.44379, -0.45864, 0.16537), 0.0),
-        )
-        .unwrap()
-        .fit()
-        .unwrap();
+    let data = data::read_pmetrics("examples/two_eq_lag/two_eq_lag.csv")?;
+    let parameters = ParameterSpace::bounded()
+        .add("ka", 0.1, 0.9)
+        .add("ke", 0.001, 0.1)
+        .add("tlag", 0.0, 4.0)
+        .add("v", 30.0, 120.0);
+    let prior = Theta::sobol_default(&parameters)?;
+    let error_models = AssayErrorModels::new().add(
+        "outeq_0",
+        AssayErrorModel::additive(ErrorPoly::new(-0.00119, 0.44379, -0.45864, 0.16537), 0.0),
+    )?;
+    EstimationProblem::nonparametric(eq, data, prior, error_models)?
+        .fit_with(NonParametricAlgorithm::npag())?;
+
+    Ok(())
 }
+
diff --git a/examples/vanco.rs b/examples/vanco.rs
index be35b5eaa..469d20bca 100644
--- a/examples/vanco.rs
+++ b/examples/vanco.rs
@@ -1,5 +1,5 @@
 use pmcore::prelude::*;
-fn main() {
+fn main() -> Result<()> {
     let eq = ode! {
         name: "vanco_two_compartment",
         params: [ke, kcp, kpc],
@@ -41,8 +41,7 @@ fn main() {
         .build();
 
     let op = eq
-        .simulate_subject_dense(&subject, &[0.3, 0.2, 0.5], None)
-        .unwrap()
+        .simulate_subject_dense(&subject, &[0.3, 0.2, 0.5], None)?
         .0;
 
     let times = op.flat_times();
@@ -51,4 +50,6 @@ fn main() {
     for (t, p) in times.iter().zip(pred.iter()) {
         println!("{}, {}", t, p);
     }
+
+    Ok(())
 }
diff --git a/examples/vanco_sde/main.rs b/examples/vanco_sde/main.rs
index 5b9e22ec4..36e77449a 100644
--- a/examples/vanco_sde/main.rs
+++ b/examples/vanco_sde/main.rs
@@ -1,6 +1,6 @@
 use pmcore::prelude::*;
 
-fn main() {
+fn main() -> Result<()> {
     let sde = sde! {
         name: "vanco_sde",
         params: [ka, ke0, kcp, kpc, vol, ske],
@@ -47,27 +47,22 @@ fn main() {
     //     (3, 1),
     // );
 
-    let data = data::read_pmetrics("examples/vanco_sde/vanco_clean.csv").unwrap();
-    EstimationProblem::builder(sde, data)
-        .parameter(Parameter::bounded("ka", 0.0001, 2.4))
-        .unwrap()
-        .parameter(Parameter::bounded("ke0", 0.0001, 2.7))
-        .unwrap()
-        .parameter(Parameter::bounded("kcp", 0.0001, 2.4))
-        .unwrap()
-        .parameter(Parameter::bounded("kpc", 0.0001, 2.4))
-        .unwrap()
-        .parameter(Parameter::bounded("vol", 0.2, 12.0))
-        .unwrap()
-        .parameter(Parameter::bounded("ske", 0.0001, 0.2))
-        .unwrap()
-        .method(Npag::new())
-        .error(
-            "outeq_1",
-            AssayErrorModel::additive(ErrorPoly::new(0.00119, 0.20, 0.0, 0.0), 0.0),
-        )
-        .unwrap()
-        .prior(Prior::sobol(100, 347))
-        .fit()
-        .unwrap();
+    let data = data::read_pmetrics("examples/vanco_sde/vanco_clean.csv")?;
+    let parameters = ParameterSpace::bounded()
+        .add("ka", 0.0001, 2.4)
+        .add("ke0", 0.0001, 2.7)
+        .add("kcp", 0.0001, 2.4)
+        .add("kpc", 0.0001, 2.4)
+        .add("vol", 0.2, 12.0)
+        .add("ske", 0.0001, 0.2);
+    let prior = Theta::sobol_default(&parameters)?;
+    let error_models = AssayErrorModels::new().add(
+        "outeq_1",
+        AssayErrorModel::additive(ErrorPoly::new(0.00119, 0.20, 0.0, 0.0), 0.0),
+    )?;
+    EstimationProblem::nonparametric(sde, data, prior, error_models)?
+        .fit_with(NonParametricAlgorithm::npag())?;
+
+    Ok(())
 }
+
diff --git a/examples/vanco_sde/vanco_sde.zip b/examples/vanco_sde/vanco_sde.zip
deleted file mode 100644
index f8a7bdde4..000000000
Binary files a/examples/vanco_sde/vanco_sde.zip and /dev/null differ
diff --git a/src/algorithms/mod.rs b/src/algorithms/mod.rs
index d0b11adb9..2c7c6c9e0 100644
--- a/src/algorithms/mod.rs
+++ b/src/algorithms/mod.rs
@@ -1,320 +1,201 @@
 use std::fs;
 use std::path::Path;
-use std::time::Instant;
 
-use crate::api::estimation_problem::NonparametricMethod;
-use crate::api::{OutputPlan, RuntimeOptions};
-use crate::estimation::nonparametric::{NonparametricWorkspace, Prior, Psi, Theta};
-use crate::model::{ModelDefinition, ParameterSpace};
-use crate::output::shared::RunConfiguration;
+use crate::estimation::nonparametric::{NonParametricResult, Psi, Theta};
+use crate::estimation::{EstimationProblem, Framework};
+use crate::results::FitResult;
+
 use anyhow::Context;
 use anyhow::Result;
 use ndarray::parallel::prelude::{IntoParallelIterator, ParallelIterator};
-use ndarray::{Array, ArrayBase, Dim, OwnedRepr};
-use nonparametric::npag::*;
-use nonparametric::npod::NPOD;
-use nonparametric::postprob::POSTPROB;
+
 use pharmsol::prelude::{data::Data, simulator::Equation};
-use pharmsol::AssayErrorModels;
+
 use pharmsol::{Predictions, Subject};
 use serde::{Deserialize, Serialize};
 
-// Module organization for algorithm types
-pub mod nonparametric;
-
-#[derive(Debug, Clone)]
-pub(crate) struct NonparametricAlgorithmInput<E: Equation> {
-    pub method: NonparametricMethod,
-    pub equation: E,
-    pub data: Data,
-    pub parameter_space: ParameterSpace,
-    pub error_models: AssayErrorModels,
-    pub output: OutputPlan,
-    pub runtime: RuntimeOptions,
+/// Defines an algorithm that can fit an [`EstimationProblem`] to produce a result.
+///
+/// Implementors are the lightweight, user-facing configuration structs (e.g.
+/// `NpagConfig`). The heavy, mutable execution state used while fitting is an
+/// internal implementation detail.
+pub trait Algorithm<E: Equation, F: crate::estimation::Framework> {
+    /// The specific result struct (e.g. `NonParametricResult<E>`).
+    type Output: FitResult;
+
+    /// Consumes the configuration and the problem, runs the optimization to
+    /// completion, and returns the strictly-typed result.
+    fn fit(self, problem: EstimationProblem<E, F>) -> Result<Self::Output>;
 }
 
-#[derive(Debug, Clone)]
-pub(crate) struct NativeNonparametricConfig {
-    pub parameter_space: ParameterSpace,
-    pub ranges: Vec<(f64, f64)>,
-    pub prior: Prior,
-    pub max_cycles: usize,
-    pub progress: bool,
-    pub run_configuration: RunConfiguration,
+// Module organization for algorithm types
+pub mod nonparametric;
+pub mod parametric;
+
+impl<E: Equation, F: Framework> EstimationProblem<E, F> {
+    /// Consumes the problem and an algorithm configuration, runs the fit to
+    /// completion, and returns the result.
+    pub fn fit_with<A>(self, algorithm: A) -> Result<A::Output>
+    where
+        A: Algorithm<E, F>,
+    {
+        algorithm.fit(self)
+    }
 }
 
-impl<E: Equation> NonparametricAlgorithmInput<E> {
-    pub(crate) fn new(
-        method: NonparametricMethod,
-        model: ModelDefinition<E>,
-        data: Data,
-        error_models: AssayErrorModels,
-        output: OutputPlan,
-        runtime: RuntimeOptions,
-    ) -> Self {
-        let ModelDefinition {
-            equation,
-            parameters,
-            ..
-        } = model;
-
-        Self {
-            method,
-            equation,
-            data,
-            parameter_space: parameters,
-            error_models,
-            output,
-            runtime,
+pub trait NonParametricRunner<E: Equation + Send + 'static>: Sync + Send + 'static {
+    /// Identify subjects whose total probability given the model is zero or
+    /// non-finite.
+    ///
+    /// Each row of [`Psi`] holds the likelihood of a subject across every
+    /// support point, so a subject's probability is the sum across its row. A
+    /// subject is flagged when that sum is zero or not finite, meaning the model
+    /// cannot explain the subject's data. When any subject is flagged, detailed
+    /// per-subject diagnostics are logged and an error is returned.
+    fn check_zero_probability_subjects(&self) -> Result<()> {
+        let psi = self.psi().matrix();
+
+        // Report non-finite entries; these propagate into the row sums below.
+        let nonfinite = psi
+            .row_iter()
+            .flat_map(|row| row.iter().copied())
+            .filter(|v| !v.is_finite())
+            .count();
+        if nonfinite > 0 {
+            tracing::warn!(
+                "Psi matrix contains {} non-finite value(s) of {} total",
+                nonfinite,
+                psi.nrows() * psi.ncols()
+            );
         }
-    }
-
-    pub(crate) fn algorithm(&self) -> Algorithm {
-        self.method.algorithm()
-    }
 
-    pub(crate) fn error_models(&self) -> &pharmsol::prelude::data::AssayErrorModels {
-        &self.error_models
-    }
-
-    pub(crate) fn max_cycles(&self) -> usize {
-        self.runtime.cycles
-    }
-
-    pub(crate) fn progress_enabled(&self) -> bool {
-        self.runtime.progress
-    }
+        // A subject's probability is the sum across its row.
+        let subjects = self.data().subjects();
+        let flagged: Vec<usize> = (0..psi.nrows())
+            .filter(|&i| {
+                let probability: f64 = (0..psi.ncols()).map(|j| psi[(i, j)]).sum();
+                !probability.is_finite() || probability == 0.0
+            })
+            .collect();
+
+        if flagged.is_empty() {
+            return Ok(());
+        }
 
-    pub(crate) fn prior(&self) -> Prior {
-        self.runtime.prior.clone().unwrap_or_default()
-    }
+        tracing::error!(
+            "{}/{} subjects have zero probability given the model",
+            flagged.len(),
+            psi.nrows()
+        );
 
-    pub(crate) fn run_configuration(&self) -> RunConfiguration {
-        RunConfiguration::new(
-            self.algorithm(),
-            &self.output,
-            &self.runtime,
-            self.parameter_space
-                .iter()
-                .map(|parameter| parameter.name.clone())
-                .collect(),
-        )
-    }
+        for &i in &flagged {
+            self.log_zero_probability_subject(subjects[i]);
+        }
 
-    pub(crate) fn native_config(&self) -> Result<NativeNonparametricConfig> {
-        Ok(NativeNonparametricConfig {
-            ranges: self.parameter_space.finite_ranges()?,
-            parameter_space: self.parameter_space.clone(),
-            prior: self.prior(),
-            max_cycles: self.max_cycles(),
-            progress: self.progress_enabled(),
-            run_configuration: self.run_configuration(),
-        })
+        let ids: Vec<&String> = flagged.iter().map(|&i| subjects[i].id()).collect();
+        Err(anyhow::anyhow!(
+            "The probability of {}/{} subjects is zero given the model. Affected subjects: {:?}",
+            flagged.len(),
+            psi.nrows(),
+            ids
+        ))
     }
-}
 
-/// Algorithm type enumeration
-///
-/// This enum represents the baseline algorithms available on the structure branch.
-#[derive(Debug, Clone, Copy, Serialize, Deserialize, PartialEq)]
-pub enum Algorithm {
-    /// Non-Parametric Adaptive Grid
-    NPAG,
-    /// Non-Parametric Optimal Design
-    NPOD,
-    /// Posterior Probability calculation
-    POSTPROB,
-}
+    /// Log detailed likelihood diagnostics for a single subject whose
+    /// probability given the model is zero or non-finite.
+    fn log_zero_probability_subject(&self, subject: &Subject) {
+        tracing::debug!("Subject with zero probability: {}", subject.id());
 
-impl Algorithm {
-    /// Check if this is a non-parametric algorithm
-    pub fn is_nonparametric(&self) -> bool {
-        matches!(
-            self,
-            Algorithm::NPAG | Algorithm::NPOD | Algorithm::POSTPROB
-        )
-    }
-}
+        let error_model = self.error_models().clone();
 
-pub trait Algorithms<E: Equation + Send + 'static>: Sync + Send + 'static {
-    fn validate_psi(&mut self) -> Result<()> {
-        // Count problematic values in psi
-        let mut nan_count = 0;
-        let mut inf_count = 0;
-
-        let psi = self.psi().to_ndarray();
-        // First coerce all NaN and infinite in psi to 0.0
-        for i in 0..psi.nrows() {
-            for j in 0..self.psi().matrix().ncols() {
-                let val = psi.get((i, j)).unwrap();
-                if val.is_nan() {
-                    nan_count += 1;
-                    // *val = 0.0;
-                } else if val.is_infinite() {
-                    inf_count += 1;
-                    // *val = 0.0;
-                }
+        // Simulate every support point for this subject in parallel.
+        let mut results: Vec<_> = self
+            .theta()
+            .matrix()
+            .row_iter()
+            .enumerate()
+            .collect::<Vec<_>>()
+            .into_par_iter()
+            .map(|(i, spp)| {
+                let support_point: Vec<f64> = spp.iter().copied().collect();
+                let (pred, ll) = self
+                    .equation()
+                    .simulate_subject_dense(subject, &support_point, Some(&error_model))
+                    .unwrap(); //TODO: Handle error
+                (i, support_point, pred.get_predictions(), ll)
+            })
+            .collect();
+
+        // Summarise the distribution of likelihood values.
+        let mut nan = 0;
+        let mut pos_inf = 0;
+        let mut neg_inf = 0;
+        let mut zero = 0;
+        let mut valid = 0;
+        for (_, _, _, ll) in &results {
+            match ll {
+                Some(v) if v.is_nan() => nan += 1,
+                Some(v) if v.is_infinite() && v.is_sign_positive() => pos_inf += 1,
+                Some(v) if v.is_infinite() => neg_inf += 1,
+                Some(v) if *v == 0.0 => zero += 1,
+                Some(_) => valid += 1,
+                None => nan += 1,
             }
         }
 
-        if nan_count + inf_count > 0 {
-            tracing::warn!(
-                "Psi matrix contains {} NaN, {} Infinite values of {} total values",
-                nan_count,
-                inf_count,
-                psi.ncols() * psi.nrows()
+        let total = results.len();
+        let pct = |n: usize| 100.0 * n as f64 / total as f64;
+        tracing::debug!(
+            "\tLikelihood analysis for subject {} ({} support points):",
+            subject.id(),
+            total
+        );
+        tracing::debug!("\tNaN likelihoods: {} ({:.1}%)", nan, pct(nan));
+        tracing::debug!("\t+Inf likelihoods: {} ({:.1}%)", pos_inf, pct(pos_inf));
+        tracing::debug!("\t-Inf likelihoods: {} ({:.1}%)", neg_inf, pct(neg_inf));
+        tracing::debug!("\tZero likelihoods: {} ({:.1}%)", zero, pct(zero));
+        tracing::debug!("\tValid likelihoods: {} ({:.1}%)", valid, pct(valid));
+
+        // Show the most likely support points to aid debugging.
+        results.sort_by(|a, b| {
+            b.3.unwrap_or(f64::NEG_INFINITY)
+                .partial_cmp(&a.3.unwrap_or(f64::NEG_INFINITY))
+                .unwrap_or(std::cmp::Ordering::Equal)
+        });
+
+        const TAKE: usize = 3;
+        tracing::debug!("Top {} most likely support points:", TAKE);
+        for (i, support_point, preds, ll) in results.iter().take(TAKE) {
+            tracing::debug!("\tSupport point #{}: {:?}", i, support_point);
+            tracing::debug!("\t\tLog-likelihood: {:?}", ll);
+            tracing::debug!(
+                "\t\tTimes: {:?}",
+                preds.iter().map(|x| x.time()).collect::<Vec<f64>>()
             );
-        }
-
-        let (_, col) = psi.dim();
-        let ecol: ArrayBase<OwnedRepr<f64>, Dim<[usize; 1]>> = Array::ones(col);
-        let plam = psi.dot(&ecol);
-        let w = 1. / &plam;
-
-        // Get the index of each element in `w` that is NaN or infinite
-        let indices: Vec<usize> = w
-            .iter()
-            .enumerate()
-            .filter(|(_, x)| x.is_nan() || x.is_infinite())
-            .map(|(i, _)| i)
-            .collect::<Vec<_>>();
-
-        if !indices.is_empty() {
-            let subject: Vec<&Subject> = self.data().subjects();
-            let zero_probability_subjects: Vec<&String> =
-                indices.iter().map(|&i| subject[i].id()).collect();
-
-            tracing::error!(
-                "{}/{} subjects have zero probability given the model",
-                indices.len(),
-                psi.nrows()
+            tracing::debug!(
+                "\t\tObservations: {:?}",
+                preds
+                    .iter()
+                    .map(|x| x.observation())
+                    .collect::<Vec<Option<f64>>>()
+            );
+            tracing::debug!(
+                "\t\tPredictions: {:?}",
+                preds.iter().map(|x| x.prediction()).collect::<Vec<f64>>()
+            );
+            tracing::debug!(
+                "\t\tOuteqs: {:?}",
+                preds.iter().map(|x| x.outeq()).collect::<Vec<usize>>()
+            );
+            tracing::debug!(
+                "\t\tStates: {:?}",
+                preds
+                    .iter()
+                    .map(|x| x.state().to_vec())
+                    .collect::<Vec<Vec<f64>>>()
             );
-
-            // For each problematic subject
-            for index in &indices {
-                tracing::debug!("Subject with zero probability: {}", subject[*index].id());
-
-                let error_model = self.error_models().clone();
-
-                // Simulate all support points in parallel
-                let spp_results: Vec<_> = self
-                    .theta()
-                    .matrix()
-                    .row_iter()
-                    .enumerate()
-                    .collect::<Vec<_>>()
-                    .into_par_iter()
-                    .map(|(i, spp)| {
-                        let support_point: Vec<f64> = spp.iter().copied().collect();
-                        let (pred, ll) = self
-                            .equation()
-                            .simulate_subject_dense(
-                                subject[*index],
-                                &support_point,
-                                Some(&error_model),
-                            )
-                            .unwrap(); //TODO: Handle error
-                        (i, support_point, pred.get_predictions(), ll)
-                    })
-                    .collect();
-
-                // Count problematic likelihoods for this subject
-                let mut nan_ll = 0;
-                let mut inf_pos_ll = 0;
-                let mut inf_neg_ll = 0;
-                let mut zero_ll = 0;
-                let mut valid_ll = 0;
-
-                for (_, _, _, ll) in &spp_results {
-                    match ll {
-                        Some(ll_val) if ll_val.is_nan() => nan_ll += 1,
-                        Some(ll_val) if ll_val.is_infinite() && ll_val.is_sign_positive() => {
-                            inf_pos_ll += 1
-                        }
-                        Some(ll_val) if ll_val.is_infinite() && ll_val.is_sign_negative() => {
-                            inf_neg_ll += 1
-                        }
-                        Some(ll_val) if *ll_val == 0.0 => zero_ll += 1,
-                        Some(_) => valid_ll += 1,
-                        None => nan_ll += 1,
-                    }
-                }
-
-                tracing::debug!(
-                    "\tLikelihood analysis for subject {} ({} support points):",
-                    subject[*index].id(),
-                    spp_results.len()
-                );
-                tracing::debug!(
-                    "\tNaN likelihoods: {} ({:.1}%)",
-                    nan_ll,
-                    100.0 * nan_ll as f64 / spp_results.len() as f64
-                );
-                tracing::debug!(
-                    "\t+Inf likelihoods: {} ({:.1}%)",
-                    inf_pos_ll,
-                    100.0 * inf_pos_ll as f64 / spp_results.len() as f64
-                );
-                tracing::debug!(
-                    "\t-Inf likelihoods: {} ({:.1}%)",
-                    inf_neg_ll,
-                    100.0 * inf_neg_ll as f64 / spp_results.len() as f64
-                );
-                tracing::debug!(
-                    "\tZero likelihoods: {} ({:.1}%)",
-                    zero_ll,
-                    100.0 * zero_ll as f64 / spp_results.len() as f64
-                );
-                tracing::debug!(
-                    "\tValid likelihoods: {} ({:.1}%)",
-                    valid_ll,
-                    100.0 * valid_ll as f64 / spp_results.len() as f64
-                );
-
-                // Sort and show top 10 most likely support points
-                let mut sorted_results = spp_results;
-                sorted_results.sort_by(|a, b| {
-                    b.3.unwrap_or(f64::NEG_INFINITY)
-                        .partial_cmp(&a.3.unwrap_or(f64::NEG_INFINITY))
-                        .unwrap_or(std::cmp::Ordering::Equal)
-                });
-                let take = 3;
-
-                tracing::debug!("Top {} most likely support points:", take);
-                for (i, support_point, preds, ll) in sorted_results.iter().take(take) {
-                    tracing::debug!("\tSupport point #{}: {:?}", i, support_point);
-                    tracing::debug!("\t\tLog-likelihood: {:?}", ll);
-
-                    let times = preds.iter().map(|x| x.time()).collect::<Vec<f64>>();
-                    let observations = preds
-                        .iter()
-                        .map(|x| x.observation())
-                        .collect::<Vec<Option<f64>>>();
-                    let predictions = preds.iter().map(|x| x.prediction()).collect::<Vec<f64>>();
-                    let outeqs = preds.iter().map(|x| x.outeq()).collect::<Vec<usize>>();
-                    let states = preds
-                        .iter()
-                        .map(|x| x.state().to_vec())
-                        .collect::<Vec<Vec<f64>>>();
-
-                    tracing::debug!("\t\tTimes: {:?}", times);
-                    tracing::debug!("\t\tObservations: {:?}", observations);
-                    tracing::debug!("\t\tPredictions: {:?}", predictions);
-                    tracing::debug!("\t\tOuteqs: {:?}", outeqs);
-                    tracing::debug!("\t\tStates: {:?}", states);
-                }
-                tracing::debug!("=====================");
-            }
-
-            return Err(anyhow::anyhow!(
-                    "The probability of {}/{} subjects is zero given the model. Affected subjects: {:?}",
-                    indices.len(),
-                    self.psi().matrix().nrows(),
-                    zero_probability_subjects
-                ));
         }
-
-        Ok(())
+        tracing::debug!("=====================");
     }
 
     fn error_models(&self) -> &pharmsol::prelude::data::AssayErrorModels;
@@ -322,7 +203,7 @@ pub trait Algorithms<E: Equation + Send + 'static>: Sync + Send + 'static {
     fn equation(&self) -> &E;
     /// Get the data used in the algorithm
     fn data(&self) -> &Data;
-    fn get_prior(&self) -> Theta;
+
     /// Increment the cycle counter and return the new value
     fn increment_cycle(&mut self) -> usize;
     /// Get the current cycle number
@@ -357,7 +238,7 @@ pub trait Algorithms<E: Equation + Send + 'static>: Sync + Send + 'static {
             fs::remove_file("stop").context("Unable to remove previous stop file")?;
         }
         self.set_status(Status::Continue);
-        self.set_theta(self.get_prior());
+
         Ok(())
     }
     fn estimation(&mut self) -> Result<()>;
@@ -405,76 +286,14 @@ pub trait Algorithms<E: Equation + Send + 'static>: Sync + Send + 'static {
     /// This method runs the full fitting process, starting with initialization,
     /// followed by iterative cycles of estimation, condensation, optimization, and evaluation
     /// until the algorithm converges or meets a stopping criteria.
-    fn fit(&mut self) -> Result<NonparametricWorkspace<E>> {
-        self.initialize().unwrap();
+    fn fit(&mut self) -> Result<NonParametricResult<E>> {
+        self.initialize()?;
         while let Status::Continue = self.next_cycle()? {}
-        self.into_workspace()
+        self.into_result()
     }
 
     #[allow(clippy::wrong_self_convention)]
-    fn into_workspace(&self) -> Result<NonparametricWorkspace<E>>;
-}
-
-pub(crate) fn dispatch_nonparametric_algorithm<E: Equation + Send + 'static>(
-    input: NonparametricAlgorithmInput<E>,
-) -> Result<Box<dyn Algorithms<E>>> {
-    match input.method {
-        NonparametricMethod::Npag(_) => {
-            let algorithm: Box<dyn Algorithms<E>> = NPAG::from_input(input)?;
-            Ok(algorithm)
-        }
-        NonparametricMethod::Npod(_) => {
-            let algorithm: Box<dyn Algorithms<E>> = NPOD::from_input(input)?;
-            Ok(algorithm)
-        }
-        NonparametricMethod::PostProb(_) => {
-            let algorithm: Box<dyn Algorithms<E>> = POSTPROB::from_input(input)?;
-            Ok(algorithm)
-        }
-    }
-}
-
-pub(crate) fn run_nonparametric_algorithm_with_progress<E, F>(
-    input: NonparametricAlgorithmInput<E>,
-    mut on_progress: F,
-) -> Result<NonparametricWorkspace<E>>
-where
-    E: Equation + Send + 'static,
-    F: FnMut(usize, f64, Option<f64>, u64, Status),
-{
-    let mut algorithm = dispatch_nonparametric_algorithm(input)?;
-    algorithm.initialize()?;
-
-    let mut previous_objective: Option<f64> = None;
-
-    loop {
-        let cycle_started = Instant::now();
-        let status = algorithm.next_cycle()?;
-        let objective = algorithm.n2ll();
-        let objective_delta = previous_objective.map(|prev| objective - prev);
-        previous_objective = Some(objective);
-
-        on_progress(
-            algorithm.cycle(),
-            objective,
-            objective_delta,
-            cycle_started.elapsed().as_millis() as u64,
-            status.clone(),
-        );
-
-        if matches!(status, Status::Stop(_)) {
-            break;
-        }
-    }
-
-    algorithm.into_workspace()
-}
-
-pub(crate) fn run_nonparametric_algorithm<E: Equation + Send + 'static>(
-    input: NonparametricAlgorithmInput<E>,
-) -> Result<NonparametricWorkspace<E>> {
-    let mut algorithm = dispatch_nonparametric_algorithm(input)?;
-    algorithm.fit()
+    fn into_result(&self) -> Result<NonParametricResult<E>>;
 }
 
 /// Represents the status/result of the algorithm
diff --git a/src/algorithms/nonparametric/mod.rs b/src/algorithms/nonparametric/mod.rs
index db8459a86..8e8eb809c 100644
--- a/src/algorithms/nonparametric/mod.rs
+++ b/src/algorithms/nonparametric/mod.rs
@@ -10,21 +10,177 @@
 //! - [`NPOD`](npod): Non-Parametric Optimal Design
 //! - [`POSTPROB`](postprob): Posterior probability reweighting
 //!
-//! # Algorithm Trait
+//! # Algorithm Selection
 //!
-//! All non-parametric algorithms implement the [`NPAlgorithm`] trait (aliased from `Algorithms`)
-//! which defines the common interface for initialization, estimation, condensation, expansion,
-//! and convergence evaluation.
+//! Use the [`NonParametricAlgorithm`] enum to select and configure an algorithm. Each
+//! variant wraps its algorithm-specific configuration struct (e.g. [`NpagConfig`]). The
+//! internal execution state used while fitting implements the [`NonParametricRunner`]
+//! trait, which defines the common interface for initialization, estimation, condensation,
+//! expansion, and convergence evaluation.
 
 // Algorithm implementations
 pub mod npag;
+pub mod npmap;
 pub mod npod;
-pub mod postprob;
 
 // Re-export algorithm structs
 pub use npag::NPAG;
+pub use npmap::NPMAP;
 pub use npod::NPOD;
-pub use postprob::POSTPROB;
 
-// Re-export the NP algorithm trait from parent
-pub use super::Algorithms as NPAlgorithm;
+// Re-export per-algorithm configuration structs
+pub use npag::NpagConfig;
+pub use npmap::NpmapConfig;
+pub use npod::NpodConfig;
+
+use crate::algorithms::{Algorithm, NonParametricRunner};
+use crate::estimation::nonparametric::NonParametricResult;
+use crate::estimation::{EstimationProblem, NonParametric};
+use anyhow::Result;
+use pharmsol::prelude::simulator::Equation;
+
+/// The non-parametric algorithms supported by PMcore.
+///
+/// Use the constructors to select an algorithm with its default configuration:
+///
+/// ```no_run
+/// use pmcore::prelude::*;
+///
+/// // Default NPAG configuration.
+/// let algorithm = NonParametricAlgorithm::npag();
+/// ```
+///
+/// To customize an algorithm, build its configuration struct (which exposes only the
+/// setters valid for that algorithm) and pass it directly to
+/// [`fit_with`](crate::estimation::EstimationProblem::fit_with):
+///
+/// ```no_run
+/// use pmcore::prelude::*;
+///
+/// // NPAG with a tighter convergence criterion and a cycle cap.
+/// let config = NpagConfig::new().eps(0.1).max_cycles(500);
+/// // `problem.fit_with(config)` accepts the config directly.
+/// ```
+///
+/// Each configuration type ([`NpagConfig`], [`NpodConfig`], [`NpmapConfig`]) implements
+/// [`Algorithm`] by delegating to the matching enum variant, so configs can be passed to
+/// `fit_with` without converting them first.
+#[derive(Debug, Clone)]
+pub enum NonParametricAlgorithm {
+    /// Non-Parametric Adaptive Grid.
+    Npag(NpagConfig),
+    /// Non-Parametric Optimal Design.
+    Npod(NpodConfig),
+    /// Non-parametric maximum a posteriori (posterior probability reweighting).
+    Npmap(NpmapConfig),
+}
+
+impl Default for NonParametricAlgorithm {
+    fn default() -> Self {
+        Self::npag()
+    }
+}
+
+impl From<NpagConfig> for NonParametricAlgorithm {
+    fn from(config: NpagConfig) -> Self {
+        Self::Npag(config)
+    }
+}
+
+impl From<NpodConfig> for NonParametricAlgorithm {
+    fn from(config: NpodConfig) -> Self {
+        Self::Npod(config)
+    }
+}
+
+impl From<NpmapConfig> for NonParametricAlgorithm {
+    fn from(config: NpmapConfig) -> Self {
+        Self::Npmap(config)
+    }
+}
+
+impl NonParametricAlgorithm {
+    /// The Non-Parametric Adaptive Grid (NPAG) algorithm with its default configuration.
+    pub fn npag() -> Self {
+        Self::Npag(NpagConfig::default())
+    }
+
+    /// The Non-Parametric Optimal Design (NPOD) algorithm with its default configuration.
+    pub fn npod() -> Self {
+        Self::Npod(NpodConfig::default())
+    }
+
+    /// The non-parametric maximum a posteriori (NPMAP) algorithm with its default
+    /// configuration.
+    pub fn npmap() -> Self {
+        Self::Npmap(NpmapConfig::default())
+    }
+}
+
+impl<E: Equation + Send + 'static> Algorithm<E, NonParametric> for NonParametricAlgorithm {
+    type Output = NonParametricResult<E>;
+
+    fn fit(self, problem: EstimationProblem<E, NonParametric>) -> Result<Self::Output> {
+        match self {
+            Self::Npag(config) => {
+                // `problem.prior` is the prior `Theta` (which also carries the parameter
+                // space) and `problem.error_models` is strictly `AssayErrorModels`.
+                let mut runner = NPAG::from_parts(
+                    problem.model.equation,
+                    problem.data,
+                    problem.error_models,
+                    problem.prior,
+                    config,
+                )?;
+                NonParametricRunner::fit(&mut runner)
+            }
+            Self::Npod(config) => {
+                let mut runner = NPOD::from_parts(
+                    problem.model.equation,
+                    problem.data,
+                    problem.error_models,
+                    problem.prior,
+                    config,
+                )?;
+                NonParametricRunner::fit(&mut runner)
+            }
+            Self::Npmap(config) => {
+                let mut runner = NPMAP::from_parts(
+                    problem.model.equation,
+                    problem.data,
+                    problem.error_models,
+                    problem.prior,
+                    config,
+                )?;
+                NonParametricRunner::fit(&mut runner)
+            }
+        }
+    }
+}
+
+// Each configuration struct delegates to its matching `NonParametricAlgorithm` variant so it
+// can be passed directly to `fit_with`. This keeps the variant-specific setters on the config
+// types (compile-time checked) while the enum remains the single source of fitting logic.
+impl<E: Equation + Send + 'static> Algorithm<E, NonParametric> for NpagConfig {
+    type Output = NonParametricResult<E>;
+
+    fn fit(self, problem: EstimationProblem<E, NonParametric>) -> Result<Self::Output> {
+        NonParametricAlgorithm::from(self).fit(problem)
+    }
+}
+
+impl<E: Equation + Send + 'static> Algorithm<E, NonParametric> for NpodConfig {
+    type Output = NonParametricResult<E>;
+
+    fn fit(self, problem: EstimationProblem<E, NonParametric>) -> Result<Self::Output> {
+        NonParametricAlgorithm::from(self).fit(problem)
+    }
+}
+
+impl<E: Equation + Send + 'static> Algorithm<E, NonParametric> for NpmapConfig {
+    type Output = NonParametricResult<E>;
+
+    fn fit(self, problem: EstimationProblem<E, NonParametric>) -> Result<Self::Output> {
+        NonParametricAlgorithm::from(self).fit(problem)
+    }
+}
diff --git a/src/algorithms/nonparametric/npag.rs b/src/algorithms/nonparametric/npag.rs
index e29c004b8..efe0caf0e 100644
--- a/src/algorithms/nonparametric/npag.rs
+++ b/src/algorithms/nonparametric/npag.rs
@@ -1,12 +1,7 @@
-use crate::algorithms::{
-    NativeNonparametricConfig, NonparametricAlgorithmInput, Status, StopReason,
-};
-use crate::api::estimation_problem::NonparametricMethod;
-use crate::api::Npag;
+use crate::algorithms::{NonParametricRunner, Status, StopReason};
 use crate::estimation::nonparametric::{
-    calculate_psi, CycleLog, NPCycle, NonparametricWorkspace, Psi, Theta, Weights,
+    calculate_psi, CycleLog, NPCycle, NonParametricResult, Psi, Theta, Weights,
 };
-use crate::prelude::algorithms::Algorithms;
 
 pub(crate) use crate::estimation::nonparametric::ipm::burke;
 pub(crate) use crate::estimation::nonparametric::qr;
@@ -20,15 +15,125 @@ use pharmsol::prelude::{
 
 use pharmsol::prelude::AssayErrorModel;
 
-use crate::estimation::nonparametric::sample_space_for_parameters;
-
 use crate::estimation::nonparametric::adaptative_grid;
 
+use serde::{Deserialize, Serialize};
+
+/// Configuration options for the Non-Parametric Adaptive Grid (NPAG) algorithm.
+#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
+pub struct NpagConfig {
+    pub eps: f64,
+    pub min_eps: f64,
+    pub objective_tolerance: f64,
+    pub pyl_tolerance: f64,
+    pub prune_threshold: f64,
+    pub qr_tolerance: f64,
+    pub grid_tolerance: f64,
+    pub error_step: f64,
+    pub min_error_step: f64,
+    pub error_step_growth: f64,
+    pub error_step_shrink: f64,
+    pub max_cycles: usize,
+    pub progress: bool,
+}
+
+impl Default for NpagConfig {
+    fn default() -> Self {
+        Self {
+            eps: 0.2,
+            min_eps: 1e-4,
+            objective_tolerance: 1e-4,
+            pyl_tolerance: 1e-2,
+            prune_threshold: 1e-3,
+            qr_tolerance: 1e-8,
+            grid_tolerance: 1e-4,
+            error_step: 0.1,
+            min_error_step: 0.01,
+            error_step_growth: 4.0,
+            error_step_shrink: 0.5,
+            max_cycles: 1000,
+            progress: true,
+        }
+    }
+}
+
+impl NpagConfig {
+    pub fn new() -> Self {
+        Self::default()
+    }
+
+    pub fn eps(mut self, eps: f64) -> Self {
+        self.eps = eps;
+        self
+    }
+
+    pub fn min_eps(mut self, min_eps: f64) -> Self {
+        self.min_eps = min_eps;
+        self
+    }
+
+    pub fn objective_tolerance(mut self, tolerance: f64) -> Self {
+        self.objective_tolerance = tolerance;
+        self
+    }
+
+    pub fn pyl_tolerance(mut self, tolerance: f64) -> Self {
+        self.pyl_tolerance = tolerance;
+        self
+    }
+
+    pub fn prune_threshold(mut self, threshold: f64) -> Self {
+        self.prune_threshold = threshold;
+        self
+    }
+
+    pub fn qr_tolerance(mut self, tolerance: f64) -> Self {
+        self.qr_tolerance = tolerance;
+        self
+    }
+
+    pub fn grid_tolerance(mut self, tolerance: f64) -> Self {
+        self.grid_tolerance = tolerance;
+        self
+    }
+
+    pub fn error_step(mut self, step: f64) -> Self {
+        self.error_step = step;
+        self
+    }
+
+    pub fn min_error_step(mut self, step: f64) -> Self {
+        self.min_error_step = step;
+        self
+    }
+
+    pub fn error_step_growth(mut self, factor: f64) -> Self {
+        self.error_step_growth = factor;
+        self
+    }
+
+    pub fn error_step_shrink(mut self, factor: f64) -> Self {
+        self.error_step_shrink = factor;
+        self
+    }
+
+    pub fn max_cycles(mut self, cycles: usize) -> Self {
+        self.max_cycles = cycles;
+        self
+    }
+
+    pub fn progress(mut self, progress: bool) -> Self {
+        self.progress = progress;
+        self
+    }
+}
+
 #[derive(Debug)]
 pub struct NPAG<E: Equation + Send + 'static> {
     equation: E,
     ranges: Vec<(f64, f64)>,
     psi: Psi,
+    prior: Theta,
     theta: Theta,
     lambda: Weights,
     w: Weights,
@@ -43,29 +148,35 @@ pub struct NPAG<E: Equation + Send + 'static> {
     status: Status,
     cycle_log: CycleLog,
     data: Data,
-    config: NativeNonparametricConfig,
-    settings: Npag,
+    config: NpagConfig,
 }
 
 impl<E: Equation + Send + 'static> NPAG<E> {
-    pub(crate) fn from_config(
+    /// Construct an `NPAG` instance from explicit parts.
+    ///
+    /// The `parameter_space` is used solely to derive the finite bounds for the
+    /// adaptive grid. Initial support points can be supplied separately via
+    /// [`NonParametricRunner::set_theta`].
+    pub(crate) fn from_parts(
         equation: E,
         data: Data,
         error_models: AssayErrorModels,
-        config: NativeNonparametricConfig,
-        settings: Npag,
-    ) -> Box<Self> {
-        let ranges = config.ranges.clone();
-        let gamma_delta = vec![settings.error_step; error_models.len()];
-
-        Box::new(Self {
+        theta: Theta,
+        config: NpagConfig,
+    ) -> Result<Self> {
+        let ranges = theta.parameters().finite_ranges();
+        let gamma_delta = vec![config.error_step; error_models.len()];
+        let eps = config.eps;
+
+        Ok(Self {
             equation,
             ranges,
             psi: Psi::new(),
-            theta: Theta::new(),
+            prior: theta.clone(),
+            theta,
             lambda: Weights::default(),
             w: Weights::default(),
-            eps: settings.eps,
+            eps,
             last_objf: -1e30,
             objf: f64::NEG_INFINITY,
             f0: -1e30,
@@ -77,45 +188,27 @@ impl<E: Equation + Send + 'static> NPAG<E> {
             cycle_log: CycleLog::new(),
             data,
             config,
-            settings,
         })
     }
-
-    pub(crate) fn from_input(input: NonparametricAlgorithmInput<E>) -> Result<Box<Self>> {
-        let method = match input.method {
-            NonparametricMethod::Npag(method) => method,
-            _ => unreachable!("NPAG::from_input requires an NPAG method"),
-        };
-        let config = input.native_config()?;
-        let error_models = input.error_models().clone();
-        let equation = input.equation;
-        let data = input.data;
-
-        Ok(Self::from_config(
-            equation,
-            data,
-            error_models,
-            config,
-            method,
-        ))
-    }
 }
 
-impl<E: Equation + Send + 'static> Algorithms<E> for NPAG<E> {
+impl<E: Equation + Send + 'static> NonParametricRunner<E> for NPAG<E> {
     fn equation(&self) -> &E {
         &self.equation
     }
-    fn into_workspace(&self) -> Result<NonparametricWorkspace<E>> {
-        NonparametricWorkspace::new(
+
+    fn into_result(&self) -> Result<NonParametricResult<E>> {
+        NonParametricResult::new(
             self.equation.clone(),
             self.data.clone(),
+            self.error_models.clone(),
+            self.prior.clone(),
             self.theta.clone(),
             self.psi.clone(),
             self.w.clone(),
             -2. * self.objf,
             self.cycle,
             self.status.clone(),
-            self.config.run_configuration.clone(),
             self.cycle_log.clone(),
         )
     }
@@ -128,10 +221,6 @@ impl<E: Equation + Send + 'static> Algorithms<E> for NPAG<E> {
         &self.data
     }
 
-    fn get_prior(&self) -> Theta {
-        sample_space_for_parameters(&self.config.parameter_space, &self.config.prior).unwrap()
-    }
-
     fn likelihood(&self) -> f64 {
         self.objf
     }
@@ -185,21 +274,21 @@ impl<E: Equation + Send + 'static> Algorithms<E> for NPAG<E> {
 
         let psi = self.psi.matrix();
         let w = &self.w;
-        if (self.last_objf - self.objf).abs() <= self.settings.objective_tolerance
-            && self.eps > self.settings.min_eps
+        if (self.last_objf - self.objf).abs() <= self.config.objective_tolerance
+            && self.eps > self.config.min_eps
         {
             self.eps /= 2.;
-            if self.eps <= self.settings.min_eps {
+            if self.eps <= self.config.min_eps {
                 let pyl = psi * w.weights();
                 self.f1 = pyl.iter().map(|x| x.ln()).sum();
-                if (self.f1 - self.f0).abs() <= self.settings.pyl_tolerance {
+                if (self.f1 - self.f0).abs() <= self.config.pyl_tolerance {
                     tracing::info!("The model converged after {} cycles", self.cycle,);
                     self.set_status(Status::Stop(StopReason::Converged));
                     self.log_cycle_state();
                     return Ok(self.status().clone());
                 } else {
                     self.f0 = self.f1;
-                    self.eps = self.settings.eps;
+                    self.eps = self.config.eps;
                 }
             }
         }
@@ -235,7 +324,7 @@ impl<E: Equation + Send + 'static> Algorithms<E> for NPAG<E> {
             self.cycle == 1 && self.config.progress,
         )?;
 
-        if let Err(err) = self.validate_psi() {
+        if let Err(err) = self.check_zero_probability_subjects() {
             bail!(err);
         }
 
@@ -258,7 +347,7 @@ impl<E: Equation + Send + 'static> Algorithms<E> for NPAG<E> {
 
         let mut keep = Vec::<usize>::new();
         for (index, lam) in self.lambda.iter().enumerate() {
-            if lam > max_lambda * self.settings.prune_threshold {
+            if lam > max_lambda * self.config.prune_threshold {
                 keep.push(index);
             }
         }
@@ -283,7 +372,7 @@ impl<E: Equation + Send + 'static> Algorithms<E> for NPAG<E> {
             let test = r.col(i).norm_l2();
             let r_diag_val = r.get(i, i);
             let ratio = r_diag_val / test;
-            if ratio.abs() >= self.settings.qr_tolerance {
+            if ratio.abs() >= self.config.qr_tolerance {
                 keep.push(*perm.get(i).unwrap());
             }
         }
@@ -301,7 +390,7 @@ impl<E: Equation + Send + 'static> Algorithms<E> for NPAG<E> {
         // Filter to keep only the support points (columns) that are in the `keep` vector
         self.psi.filter_column_indices(keep.as_slice());
 
-        self.validate_psi()?;
+        self.check_zero_probability_subjects()?;
         (self.lambda, self.objf) = match burke(&self.psi) {
             Ok((lambda, objf)) => (lambda, objf),
             Err(err) => {
@@ -368,20 +457,20 @@ impl<E: Equation + Send + 'static> Algorithms<E> for NPAG<E> {
                 if objf_up > self.objf {
                     self.error_models.set_factor(outeq, gamma_up)?;
                     self.objf = objf_up;
-                    self.gamma_delta[outeq] *= self.settings.error_step_growth;
+                    self.gamma_delta[outeq] *= self.config.error_step_growth;
                     self.lambda = lambda_up;
                     self.psi = psi_up;
                 }
                 if objf_down > self.objf {
                     self.error_models.set_factor(outeq, gamma_down)?;
                     self.objf = objf_down;
-                    self.gamma_delta[outeq] *= self.settings.error_step_growth;
+                    self.gamma_delta[outeq] *= self.config.error_step_growth;
                     self.lambda = lambda_down;
                     self.psi = psi_down;
                 }
-                self.gamma_delta[outeq] *= self.settings.error_step_shrink;
-                if self.gamma_delta[outeq] <= self.settings.min_error_step {
-                    self.gamma_delta[outeq] = self.settings.error_step;
+                self.gamma_delta[outeq] *= self.config.error_step_shrink;
+                if self.gamma_delta[outeq] <= self.config.min_error_step {
+                    self.gamma_delta[outeq] = self.config.error_step;
                 }
                 Ok(())
             })?;
@@ -394,7 +483,7 @@ impl<E: Equation + Send + 'static> Algorithms<E> for NPAG<E> {
             &mut self.theta,
             self.eps,
             &self.ranges,
-            self.settings.grid_tolerance,
+            self.config.grid_tolerance,
         )?;
         Ok(())
     }
@@ -424,10 +513,9 @@ impl<E: Equation + Send + 'static> Algorithms<E> for NPAG<E> {
 
 #[cfg(test)]
 mod tests {
-    use super::*;
-    use crate::api::{OutputPlan, RuntimeOptions};
-    use crate::model::{ModelDefinition, Parameter};
-    use pharmsol::{fa, fetch_params, lag, AssayErrorModel, ErrorPoly, Subject, SubjectBuilderExt};
+    use crate::prelude::*;
+
+    use pharmsol::{fa, fetch_params, lag, Subject, SubjectBuilderExt};
 
     fn simple_equation() -> pharmsol::equation::ODE {
         pharmsol::equation::ODE::new(
@@ -467,44 +555,27 @@ mod tests {
     }
 
     #[test]
-    fn from_input_uses_npag_method_settings() -> Result<()> {
-        let method = Npag::new()
-            .eps(0.125)
-            .min_eps(0.0025)
-            .objective_tolerance(3e-5)
-            .pyl_tolerance(4e-3)
-            .prune_threshold(2e-3)
-            .qr_tolerance(9e-7)
-            .grid_tolerance(8e-4)
-            .error_step(0.2)
-            .min_error_step(0.03)
-            .error_step_growth(3.0)
-            .error_step_shrink(0.25);
-
-        let model = ModelDefinition::builder(simple_equation())
-            .parameter(Parameter::bounded("ke", 0.1, 1.0))?
-            .parameter(Parameter::bounded("v", 1.0, 20.0))?
-            .build()?;
-
-        let error_models = AssayErrorModels::new().add(
-            0,
-            AssayErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0),
-        )?;
-
-        let input = NonparametricAlgorithmInput::new(
-            NonparametricMethod::Npag(method),
-            model,
-            simple_data(),
-            error_models,
-            OutputPlan::disabled(),
-            RuntimeOptions::default(),
+    fn npag_runs_without_error() {
+        let parameters = ParameterSpace::bounded()
+            .add("ke", 0.001, 3.0)
+            .add("v", 25.0, 250.0);
+        let prior = Theta::sobol_default(&parameters).expect("Failed to build prior");
+        let error_models = AssayErrorModels::new()
+            .add(
+                "0",
+                AssayErrorModel::additive(ErrorPoly::new(0.0, 0.5, 0.0, 0.0), 0.0),
+            )
+            .expect("Failed to build error models");
+        let problem =
+            EstimationProblem::nonparametric(simple_equation(), simple_data(), prior, error_models)
+                .expect("Failed to build problem");
+
+        let result = problem.fit_with(NonParametricAlgorithm::npag());
+
+        assert!(
+            result.is_ok(),
+            "NPAG algorithm should run without error, but got: {:?}",
+            result.err()
         );
-
-        let algorithm = NPAG::from_input(input)?;
-
-        assert_eq!(algorithm.settings, method);
-        assert_eq!(algorithm.eps, method.eps);
-        assert_eq!(algorithm.gamma_delta, vec![method.error_step]);
-        Ok(())
     }
 }
diff --git a/src/algorithms/nonparametric/postprob.rs b/src/algorithms/nonparametric/npmap.rs
similarity index 62%
rename from src/algorithms/nonparametric/postprob.rs
rename to src/algorithms/nonparametric/npmap.rs
index 3e778b009..92b122f64 100644
--- a/src/algorithms/nonparametric/postprob.rs
+++ b/src/algorithms/nonparametric/npmap.rs
@@ -1,23 +1,34 @@
 use crate::{
-    algorithms::{NativeNonparametricConfig, NonparametricAlgorithmInput, Status, StopReason},
+    algorithms::{NonParametricRunner, Status, StopReason},
     estimation::nonparametric::{
-        calculate_psi, CycleLog, NPCycle, NonparametricWorkspace, Psi, Theta, Weights,
+        calculate_psi, CycleLog, NPCycle, NonParametricResult, Psi, Theta, Weights,
     },
-    prelude::algorithms::Algorithms,
 };
-use anyhow::{Context, Result};
 
+use anyhow::{Context, Result};
 use pharmsol::prelude::{
     data::{AssayErrorModels, Data},
     simulator::Equation,
 };
 
 use crate::estimation::nonparametric::ipm::burke;
-use crate::estimation::nonparametric::sample_space_for_parameters;
+use serde::{Deserialize, Serialize};
+
+/// Configuration options for the non-parametric maximum a posteriori (NPMAP) algorithm
+#[derive(Debug, Clone, PartialEq, Default, Serialize, Deserialize)]
+pub struct NpmapConfig {}
 
-/// Posterior probability algorithm
-/// Reweights the prior probabilities to the observed data and error model
-pub struct POSTPROB<E: Equation + Send + 'static> {
+impl NpmapConfig {
+    pub fn new() -> Self {
+        Self::default()
+    }
+}
+
+/// Non-parametric maximum a posteriori (NPMAP) algorithm
+///
+/// This algorithm is a wrapper around the IPM algorithm that calculates the posterior probabilities of the support points
+/// given a prior distribution and the likelihood of the data.
+pub struct NPMAP<E: Equation + Send + 'static> {
     equation: E,
     psi: Psi,
     theta: Theta,
@@ -26,26 +37,52 @@ pub struct POSTPROB<E: Equation + Send + 'static> {
     cycle: usize,
     status: Status,
     data: Data,
-    config: NativeNonparametricConfig,
     cyclelog: CycleLog,
     error_models: AssayErrorModels,
+    prior: Theta,
 }
 
-impl<E: Equation + Send + 'static> Algorithms<E> for POSTPROB<E> {
-    fn into_workspace(&self) -> Result<NonparametricWorkspace<E>> {
-        NonparametricWorkspace::new(
+impl<E: Equation + Send + 'static> NPMAP<E> {
+    pub(crate) fn from_parts(
+        equation: E,
+        data: Data,
+        error_models: AssayErrorModels,
+        theta: Theta,
+        _config: NpmapConfig,
+    ) -> Result<Self> {
+        Ok(Self {
+            equation,
+            psi: Psi::new(),
+            theta: theta.clone(),
+            w: Weights::default(),
+            objf: f64::INFINITY,
+            cycle: 0,
+            status: Status::Continue,
+            data,
+            cyclelog: CycleLog::new(),
+            error_models,
+            prior: theta,
+        })
+    }
+}
+
+impl<E: Equation + Send + 'static> NonParametricRunner<E> for NPMAP<E> {
+    fn into_result(&self) -> Result<NonParametricResult<E>> {
+        NonParametricResult::new(
             self.equation.clone(),
             self.data.clone(),
+            self.error_models.clone(),
+            self.prior.clone(),
             self.theta.clone(),
             self.psi.clone(),
             self.w.clone(),
             self.objf,
             self.cycle,
             self.status.clone(),
-            self.config.run_configuration.clone(),
             self.cyclelog.clone(),
         )
     }
+
     fn error_models(&self) -> &AssayErrorModels {
         &self.error_models
     }
@@ -58,10 +95,6 @@ impl<E: Equation + Send + 'static> Algorithms<E> for POSTPROB<E> {
         &self.data
     }
 
-    fn get_prior(&self) -> Theta {
-        sample_space_for_parameters(&self.config.parameter_space, &self.config.prior).unwrap()
-    }
-
     fn likelihood(&self) -> f64 {
         self.objf
     }
@@ -114,6 +147,7 @@ impl<E: Equation + Send + 'static> Algorithms<E> for POSTPROB<E> {
     fn condensation(&mut self) -> Result<()> {
         Ok(())
     }
+
     fn optimizations(&mut self) -> Result<()> {
         Ok(())
     }
@@ -135,27 +169,14 @@ impl<E: Equation + Send + 'static> Algorithms<E> for POSTPROB<E> {
         );
         self.cyclelog.push(state);
     }
-}
 
-impl<E: Equation + Send + 'static> POSTPROB<E> {
-    pub(crate) fn from_input(input: NonparametricAlgorithmInput<E>) -> Result<Box<Self>> {
-        let config = input.native_config()?;
-        let error_models = input.error_models().clone();
-        let equation = input.equation;
-        let data = input.data;
+    /// POSTPROB is a single-pass reweighting: it evaluates the likelihood of the
+    /// fixed prior support points once, rather than iterating cycles.
+    fn fit(&mut self) -> Result<NonParametricResult<E>> {
+        self.estimation()?;
+        self.evaluation()?;
+        self.log_cycle_state();
 
-        Ok(Box::new(Self {
-            equation,
-            psi: Psi::new(),
-            theta: Theta::new(),
-            w: Weights::default(),
-            objf: f64::INFINITY,
-            cycle: 0,
-            status: Status::Continue,
-            data,
-            config,
-            cyclelog: CycleLog::new(),
-            error_models,
-        }))
+        self.into_result()
     }
 }
diff --git a/src/algorithms/nonparametric/npod.rs b/src/algorithms/nonparametric/npod.rs
index 83476b89e..634c673f8 100644
--- a/src/algorithms/nonparametric/npod.rs
+++ b/src/algorithms/nonparametric/npod.rs
@@ -1,29 +1,61 @@
-use crate::algorithms::{NativeNonparametricConfig, NonparametricAlgorithmInput, StopReason};
-use crate::estimation::nonparametric::ipm::burke;
-use crate::estimation::nonparametric::qr;
-use crate::estimation::nonparametric::{
-    calculate_psi, CycleLog, NPCycle, NonparametricWorkspace, Psi, Theta, Weights,
+use crate::{
+    algorithms::{NonParametricRunner, Status, StopReason},
+    estimation::nonparametric::{
+        calculate_psi, ipm::burke, qr, CycleLog, NPCycle, NonParametricResult, Psi, Theta, Weights,
+    },
 };
-use crate::{algorithms::Status, prelude::algorithms::Algorithms};
 use pharmsol::ParameterOptimizer;
 
 use anyhow::bail;
 use anyhow::Result;
+use pharmsol::prelude::{data::Data, simulator::Equation};
 use pharmsol::{prelude::AssayErrorModel, AssayErrorModels};
-use pharmsol::{
-    prelude::{data::Data, simulator::Equation},
-    Subject,
-};
 
 use ndarray::Array1;
 use rayon::prelude::{IntoParallelRefMutIterator, ParallelIterator};
+use serde::{Deserialize, Serialize};
 
 const THETA_F: f64 = 1e-2;
 const THETA_D: f64 = 1e-4;
 
+/// Configuration options for the Non-Parametric Optimal Design (NPOD) algorithm.
+#[derive(Debug, Clone, Serialize, Deserialize)]
+pub struct NpodConfig {
+    /// Maximum number of cycles to run the algorithm for.
+    pub max_cycles: usize,
+    /// Whether to print progress information during the first cycle.
+    pub progress: bool,
+}
+
+impl NpodConfig {
+    pub fn new() -> Self {
+        Self::default()
+    }
+
+    pub fn max_cycles(mut self, cycles: usize) -> Self {
+        self.max_cycles = cycles;
+        self
+    }
+
+    pub fn progress(mut self, progress: bool) -> Self {
+        self.progress = progress;
+        self
+    }
+}
+
+impl Default for NpodConfig {
+    fn default() -> Self {
+        Self {
+            max_cycles: 100,
+            progress: true,
+        }
+    }
+}
+
 pub struct NPOD<E: Equation + Send + 'static> {
     equation: E,
     psi: Psi,
+    prior: Theta,
     theta: Theta,
     lambda: Weights,
     w: Weights,
@@ -36,21 +68,53 @@ pub struct NPOD<E: Equation + Send + 'static> {
     status: Status,
     cycle_log: CycleLog,
     data: Data,
-    config: NativeNonparametricConfig,
+    config: NpodConfig,
+}
+
+impl<E: Equation + Send + 'static> NPOD<E> {
+    pub(crate) fn from_parts(
+        equation: E,
+        data: Data,
+        error_models: AssayErrorModels,
+        theta: Theta,
+        config: NpodConfig,
+    ) -> Result<Self> {
+        let gamma_delta = vec![0.1; error_models.len()];
+
+        Ok(Self {
+            equation,
+            psi: Psi::new(),
+            prior: theta.clone(),
+            theta: theta,
+            lambda: Weights::default(),
+            w: Weights::default(),
+            last_objf: -1e30,
+            objf: f64::NEG_INFINITY,
+            cycle: 0,
+            gamma_delta,
+            error_models,
+            converged: false,
+            status: Status::Continue,
+            cycle_log: CycleLog::new(),
+            data,
+            config,
+        })
+    }
 }
 
-impl<E: Equation + Send + 'static> Algorithms<E> for NPOD<E> {
-    fn into_workspace(&self) -> Result<NonparametricWorkspace<E>> {
-        NonparametricWorkspace::new(
+impl<E: Equation + Send + 'static> NonParametricRunner<E> for NPOD<E> {
+    fn into_result(&self) -> Result<NonParametricResult<E>> {
+        NonParametricResult::new(
             self.equation.clone(),
             self.data.clone(),
+            self.error_models.clone(),
+            self.prior.clone(),
             self.theta.clone(),
             self.psi.clone(),
             self.w.clone(),
             -2. * self.objf,
             self.cycle,
             self.status.clone(),
-            self.config.run_configuration.clone(),
             self.cycle_log.clone(),
         )
     }
@@ -67,14 +131,6 @@ impl<E: Equation + Send + 'static> Algorithms<E> for NPOD<E> {
         &self.data
     }
 
-    fn get_prior(&self) -> Theta {
-        crate::estimation::nonparametric::sample_space_for_parameters(
-            &self.config.parameter_space,
-            &self.config.prior,
-        )
-        .unwrap()
-    }
-
     fn increment_cycle(&mut self) -> usize {
         self.cycle += 1;
         self.cycle
@@ -184,7 +240,7 @@ impl<E: Equation + Send + 'static> Algorithms<E> for NPOD<E> {
             self.cycle == 1 && self.config.progress,
         )?;
 
-        if let Err(err) = self.validate_psi() {
+        if let Err(err) = self.check_zero_probability_subjects() {
             bail!(err);
         }
 
@@ -347,73 +403,3 @@ impl<E: Equation + Send + 'static> Algorithms<E> for NPOD<E> {
         Ok(())
     }
 }
-
-impl<E: Equation + Send + 'static> NPOD<E> {
-    pub(crate) fn from_input(input: NonparametricAlgorithmInput<E>) -> Result<Box<Self>> {
-        let config = input.native_config()?;
-        let error_models = input.error_models().clone();
-        let gamma_delta = vec![0.1; error_models.len()];
-        let equation = input.equation;
-        let data = input.data;
-
-        Ok(Box::new(Self {
-            equation,
-            psi: Psi::new(),
-            theta: Theta::new(),
-            lambda: Weights::default(),
-            w: Weights::default(),
-            last_objf: -1e30,
-            objf: f64::NEG_INFINITY,
-            cycle: 0,
-            gamma_delta,
-            error_models,
-            converged: false,
-            status: Status::Continue,
-            cycle_log: CycleLog::new(),
-            data,
-            config,
-        }))
-    }
-
-    fn validate_psi(&mut self) -> Result<()> {
-        let mut psi = self.psi().matrix().to_owned();
-        // First coerce all NaN and infinite in psi to 0.0
-        let mut has_bad_values = false;
-        for i in 0..psi.nrows() {
-            for j in 0..psi.ncols() {
-                let val = psi[(i, j)];
-                if val.is_nan() || val.is_infinite() {
-                    has_bad_values = true;
-                    psi[(i, j)] = 0.0;
-                }
-            }
-        }
-        if has_bad_values {
-            tracing::warn!("Psi contains NaN or Inf values, coercing to 0.0");
-        }
-
-        // Calculate row sums and check for zero-probability subjects
-        let nrows = psi.nrows();
-        let ncols = psi.ncols();
-        let indices: Vec<usize> = (0..nrows)
-            .filter(|&i| {
-                let row_sum: f64 = (0..ncols).map(|j| psi[(i, j)]).sum();
-                let w: f64 = 1.0 / row_sum;
-                w.is_nan() || w.is_infinite()
-            })
-            .collect();
-
-        // If any elements in `w` are NaN or infinite, return the subject IDs for each index
-        if !indices.is_empty() {
-            let subject: Vec<&Subject> = self.data.subjects();
-            let zero_probability_subjects: Vec<&String> =
-                indices.iter().map(|&i| subject[i].id()).collect();
-
-            return Err(anyhow::anyhow!(
-                "The probability of one or more subjects, given the model, is zero. The following subjects have zero probability: {:?}", zero_probability_subjects
-            ));
-        }
-
-        Ok(())
-    }
-}
diff --git a/src/algorithms/parametric/mod.rs b/src/algorithms/parametric/mod.rs
new file mode 100644
index 000000000..9951345f1
--- /dev/null
+++ b/src/algorithms/parametric/mod.rs
@@ -0,0 +1,100 @@
+//! Parametric algorithm implementations.
+//!
+//! Parametric algorithms estimate the population distribution as a parametric model (for
+//! example a multivariate normal) whose parameters are fitted to the data.
+//!
+//! # Algorithm Selection
+//!
+//! Use the [`ParametricAlgorithm`] enum to select and configure an algorithm. Each variant
+//! wraps its algorithm-specific configuration struct (e.g. [`SaemConfig`]).
+//!
+//! Note: the parametric fitting machinery is not yet implemented. Constructing a problem and
+//! calling [`fit_with`](crate::estimation::EstimationProblem::fit_with) with a
+//! [`ParametricAlgorithm`] type-checks today, but running it will panic until the SAEM solver
+//! is implemented.
+
+pub mod saem_config;
+
+pub use saem_config::SaemConfig;
+
+use crate::algorithms::Algorithm;
+use crate::estimation::{EstimationProblem, Parametric};
+use crate::results::ParametricResult;
+use anyhow::Result;
+use pharmsol::prelude::simulator::Equation;
+
+/// The parametric algorithms supported by PMcore.
+///
+/// Use the constructors to select an algorithm with its default configuration:
+///
+/// ```no_run
+/// use pmcore::prelude::*;
+///
+/// // Default SAEM configuration.
+/// let algorithm = ParametricAlgorithm::saem();
+/// ```
+///
+/// To customize an algorithm, build its configuration struct (which exposes only the
+/// setters valid for that algorithm) and pass it directly to
+/// [`fit_with`](crate::estimation::EstimationProblem::fit_with):
+///
+/// ```no_run
+/// use pmcore::prelude::*;
+///
+/// // SAEM with a custom iteration schedule and seed.
+/// let config = SaemConfig::new()
+///     .k1_iterations(500)
+///     .k2_iterations(200)
+///     .seed(42);
+/// // `problem.fit_with(config)` accepts the config directly.
+/// ```
+///
+/// [`SaemConfig`] implements [`Algorithm`] by delegating to the matching enum variant, so
+/// configs can be passed to `fit_with` without converting them first.
+#[derive(Debug, Clone)]
+pub enum ParametricAlgorithm {
+    /// Stochastic Approximation Expectation-Maximization.
+    Saem(SaemConfig),
+}
+
+impl Default for ParametricAlgorithm {
+    fn default() -> Self {
+        Self::saem()
+    }
+}
+
+impl From<SaemConfig> for ParametricAlgorithm {
+    fn from(config: SaemConfig) -> Self {
+        Self::Saem(config)
+    }
+}
+
+impl ParametricAlgorithm {
+    /// The Stochastic Approximation Expectation-Maximization (SAEM) algorithm with its
+    /// default configuration.
+    pub fn saem() -> Self {
+        Self::Saem(SaemConfig::default())
+    }
+}
+
+impl<E: Equation + Send + 'static> Algorithm<E, Parametric> for ParametricAlgorithm {
+    type Output = ParametricResult<E>;
+
+    fn fit(self, _problem: EstimationProblem<E, Parametric>) -> Result<Self::Output> {
+        match self {
+            Self::Saem(_config) => {
+                unimplemented!("SAEM fitting is not yet implemented")
+            }
+        }
+    }
+}
+
+// `SaemConfig` delegates to the matching `ParametricAlgorithm` variant so it can be passed
+// directly to `fit_with`, keeping its setters compile-time checked.
+impl<E: Equation + Send + 'static> Algorithm<E, Parametric> for SaemConfig {
+    type Output = ParametricResult<E>;
+
+    fn fit(self, problem: EstimationProblem<E, Parametric>) -> Result<Self::Output> {
+        ParametricAlgorithm::from(self).fit(problem)
+    }
+}
diff --git a/src/api/saem_config.rs b/src/algorithms/parametric/saem_config.rs
similarity index 77%
rename from src/api/saem_config.rs
rename to src/algorithms/parametric/saem_config.rs
index 7f73f373c..289c41750 100644
--- a/src/api/saem_config.rs
+++ b/src/algorithms/parametric/saem_config.rs
@@ -61,6 +61,47 @@ impl Default for SaemConfig {
 }
 
 impl SaemConfig {
+    /// Creates a new `SaemConfig` with default values.
+    pub fn new() -> Self {
+        Self::default()
+    }
+
+    /// Number of exploration-phase (K1) iterations.
+    pub fn k1_iterations(mut self, iterations: usize) -> Self {
+        self.k1_iterations = iterations;
+        self
+    }
+
+    /// Number of smoothing-phase (K2) iterations.
+    pub fn k2_iterations(mut self, iterations: usize) -> Self {
+        self.k2_iterations = iterations;
+        self
+    }
+
+    /// Number of burn-in iterations.
+    pub fn burn_in(mut self, burn_in: usize) -> Self {
+        self.burn_in = burn_in;
+        self
+    }
+
+    /// Number of MCMC chains.
+    pub fn n_chains(mut self, n_chains: usize) -> Self {
+        self.n_chains = n_chains;
+        self
+    }
+
+    /// MCMC step size.
+    pub fn mcmc_step_size(mut self, step_size: f64) -> Self {
+        self.mcmc_step_size = step_size;
+        self
+    }
+
+    /// Random-number-generator seed.
+    pub fn seed(mut self, seed: u64) -> Self {
+        self.seed = seed;
+        self
+    }
+
     pub fn total_iterations(&self) -> usize {
         self.k1_iterations + self.k2_iterations
     }
diff --git a/src/api/estimation_problem.rs b/src/api/estimation_problem.rs
deleted file mode 100644
index b82af45ec..000000000
--- a/src/api/estimation_problem.rs
+++ /dev/null
@@ -1,566 +0,0 @@
-use anyhow::Result;
-use pharmsol::{AssayErrorModel, AssayErrorModels, Data, Equation};
-use serde::{Deserialize, Serialize};
-
-use crate::algorithms::Algorithm;
-use crate::api::error_models::ErrorModels;
-use crate::api::SaemConfig;
-use crate::estimation::nonparametric::Prior;
-use crate::model::{
-    CovariateSpec, EquationMetadataSource, ModelDefinition, ModelDefinitionBuilder, ModelMetadata,
-    Parameter, VariabilityModel,
-};
-
-// =============================================================================
-// Method selection
-// =============================================================================
-
-pub trait MethodSpec {
-    type Builder<E: Equation>;
-
-    fn attach<E: Equation>(self, builder: EstimationProblemBuilder<E>) -> Self::Builder<E>;
-}
-
-#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
-pub struct Npag {
-    pub eps: f64,
-    pub min_eps: f64,
-    pub objective_tolerance: f64,
-    pub pyl_tolerance: f64,
-    pub prune_threshold: f64,
-    pub qr_tolerance: f64,
-    pub grid_tolerance: f64,
-    pub error_step: f64,
-    pub min_error_step: f64,
-    pub error_step_growth: f64,
-    pub error_step_shrink: f64,
-}
-
-impl Default for Npag {
-    fn default() -> Self {
-        Self {
-            eps: 0.2,
-            min_eps: 1e-4,
-            objective_tolerance: 1e-4,
-            pyl_tolerance: 1e-2,
-            prune_threshold: 1e-3,
-            qr_tolerance: 1e-8,
-            grid_tolerance: 1e-4,
-            error_step: 0.1,
-            min_error_step: 0.01,
-            error_step_growth: 4.0,
-            error_step_shrink: 0.5,
-        }
-    }
-}
-
-impl Npag {
-    pub fn new() -> Self {
-        Self::default()
-    }
-
-    pub fn eps(mut self, eps: f64) -> Self {
-        self.eps = eps;
-        self
-    }
-
-    pub fn min_eps(mut self, min_eps: f64) -> Self {
-        self.min_eps = min_eps;
-        self
-    }
-
-    pub fn objective_tolerance(mut self, tolerance: f64) -> Self {
-        self.objective_tolerance = tolerance;
-        self
-    }
-
-    pub fn pyl_tolerance(mut self, tolerance: f64) -> Self {
-        self.pyl_tolerance = tolerance;
-        self
-    }
-
-    pub fn prune_threshold(mut self, threshold: f64) -> Self {
-        self.prune_threshold = threshold;
-        self
-    }
-
-    pub fn qr_tolerance(mut self, tolerance: f64) -> Self {
-        self.qr_tolerance = tolerance;
-        self
-    }
-
-    pub fn grid_tolerance(mut self, tolerance: f64) -> Self {
-        self.grid_tolerance = tolerance;
-        self
-    }
-
-    pub fn error_step(mut self, step: f64) -> Self {
-        self.error_step = step;
-        self
-    }
-
-    pub fn min_error_step(mut self, step: f64) -> Self {
-        self.min_error_step = step;
-        self
-    }
-
-    pub fn error_step_growth(mut self, factor: f64) -> Self {
-        self.error_step_growth = factor;
-        self
-    }
-
-    pub fn error_step_shrink(mut self, factor: f64) -> Self {
-        self.error_step_shrink = factor;
-        self
-    }
-}
-
-#[derive(Debug, Clone, Copy, Default, PartialEq, Eq, Serialize, Deserialize)]
-pub struct Npod;
-
-impl Npod {
-    pub fn new() -> Self {
-        Self
-    }
-}
-
-#[derive(Debug, Clone, Copy, Default, PartialEq, Eq, Serialize, Deserialize)]
-pub struct PostProb;
-
-impl PostProb {
-    pub fn new() -> Self {
-        Self
-    }
-}
-
-#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
-pub(crate) enum NonparametricMethod {
-    Npag(Npag),
-    Npod(Npod),
-    PostProb(PostProb),
-}
-
-impl NonparametricMethod {
-    pub fn algorithm(self) -> Algorithm {
-        match self {
-            Self::Npag(_) => Algorithm::NPAG,
-            Self::Npod(_) => Algorithm::NPOD,
-            Self::PostProb(_) => Algorithm::POSTPROB,
-        }
-    }
-}
-
-impl MethodSpec for Npag {
-    type Builder<E: Equation> = NonparametricEstimationProblemBuilder<E>;
-
-    fn attach<E: Equation>(self, builder: EstimationProblemBuilder<E>) -> Self::Builder<E> {
-        NonparametricEstimationProblemBuilder::new(builder, NonparametricMethod::Npag(self))
-    }
-}
-
-impl MethodSpec for Npod {
-    type Builder<E: Equation> = NonparametricEstimationProblemBuilder<E>;
-
-    fn attach<E: Equation>(self, builder: EstimationProblemBuilder<E>) -> Self::Builder<E> {
-        NonparametricEstimationProblemBuilder::new(builder, NonparametricMethod::Npod(self))
-    }
-}
-
-impl MethodSpec for PostProb {
-    type Builder<E: Equation> = NonparametricEstimationProblemBuilder<E>;
-
-    fn attach<E: Equation>(self, builder: EstimationProblemBuilder<E>) -> Self::Builder<E> {
-        NonparametricEstimationProblemBuilder::new(builder, NonparametricMethod::PostProb(self))
-    }
-}
-
-#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
-#[serde(default)]
-pub struct OutputPlan {
-    pub write: bool,
-    pub path: Option<String>,
-}
-
-impl OutputPlan {
-    pub fn disabled() -> Self {
-        Self {
-            write: false,
-            path: None,
-        }
-    }
-}
-
-impl Default for OutputPlan {
-    fn default() -> Self {
-        Self::disabled()
-    }
-}
-
-#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
-#[serde(default)]
-pub struct ConvergenceOptions {
-    pub likelihood: f64,
-    pub pyl: f64,
-    pub eps: f64,
-}
-
-impl Default for ConvergenceOptions {
-    fn default() -> Self {
-        Self {
-            likelihood: 1e-4,
-            pyl: 1e-2,
-            eps: 1e-2,
-        }
-    }
-}
-
-#[derive(Debug, Clone, Serialize, Deserialize)]
-#[serde(default)]
-pub struct AlgorithmTuning {
-    pub min_distance: f64,
-    pub nm_steps: usize,
-    pub tolerance: f64,
-    pub saem: SaemConfig,
-}
-
-impl Default for AlgorithmTuning {
-    fn default() -> Self {
-        Self {
-            min_distance: 1e-4,
-            nm_steps: 100,
-            tolerance: 1e-6,
-            saem: SaemConfig::default(),
-        }
-    }
-}
-
-#[derive(Debug, Clone, Serialize, Deserialize)]
-#[serde(default)]
-pub struct RuntimeOptions {
-    pub cycles: usize,
-    pub cache: bool,
-    pub progress: bool,
-    pub idelta: f64,
-    pub tad: f64,
-    pub prior: Option<Prior>,
-    pub convergence: ConvergenceOptions,
-    pub tuning: AlgorithmTuning,
-}
-
-impl Default for RuntimeOptions {
-    fn default() -> Self {
-        Self {
-            cycles: 100,
-            cache: true,
-            progress: true,
-            idelta: 0.12,
-            tad: 0.0,
-            prior: None,
-            convergence: ConvergenceOptions::default(),
-            tuning: AlgorithmTuning::default(),
-        }
-    }
-}
-
-#[derive(Debug, Clone)]
-pub struct EstimationProblem<E: Equation> {
-    pub(crate) model: ModelDefinition<E>,
-    pub(crate) data: Data,
-    pub(crate) error_models: ErrorModels,
-    pub(crate) method: NonparametricMethod,
-    pub(crate) output: OutputPlan,
-    pub(crate) runtime: RuntimeOptions,
-}
-
-impl<E: Equation> EstimationProblem<E> {
-    pub fn builder(equation: E, data: Data) -> EstimationProblemBuilder<E> {
-        EstimationProblemBuilder {
-            model: ModelDefinition::builder(equation),
-            data,
-            output: Some(OutputPlan::default()),
-            runtime: Some(RuntimeOptions::default()),
-        }
-    }
-}
-
-impl<E: Equation + Clone + Send + 'static + EquationMetadataSource> EstimationProblem<E> {
-    pub fn fit(self) -> Result<crate::results::FitResult<E>> {
-        crate::api::fit(self)
-    }
-
-    pub fn fit_with_progress<F>(self, on_progress: F) -> Result<crate::results::FitResult<E>>
-    where
-        F: FnMut(crate::api::FitProgress),
-    {
-        crate::api::fit_with_progress(self, on_progress)
-    }
-}
-
-pub struct EstimationProblemBuilder<E: Equation> {
-    model: ModelDefinitionBuilder<E>,
-    data: Data,
-    output: Option<OutputPlan>,
-    runtime: Option<RuntimeOptions>,
-}
-
-impl<E: Equation> EstimationProblemBuilder<E> {
-    pub fn method<M: MethodSpec>(self, method: M) -> M::Builder<E> {
-        method.attach(self)
-    }
-
-    pub fn output_dir(self, path: impl Into<String>) -> Self {
-        self.with_output_plan(|output| {
-            output.write = true;
-            output.path = Some(path.into());
-        })
-    }
-
-    pub fn no_output(self) -> Self {
-        self.with_output_plan(|output| {
-            output.write = false;
-            output.path = None;
-        })
-    }
-
-    pub fn cycles(self, cycles: usize) -> Self {
-        self.with_runtime_options(|runtime| runtime.cycles = cycles)
-    }
-
-    pub fn cache(self, enabled: bool) -> Self {
-        self.with_runtime_options(|runtime| runtime.cache = enabled)
-    }
-
-    pub fn progress(self, enabled: bool) -> Self {
-        self.with_runtime_options(|runtime| runtime.progress = enabled)
-    }
-
-    pub fn idelta(self, value: f64) -> Self {
-        self.with_runtime_options(|runtime| runtime.idelta = value)
-    }
-
-    pub fn tad(self, value: f64) -> Self {
-        self.with_runtime_options(|runtime| runtime.tad = value)
-    }
-
-    pub fn prior(self, prior: Prior) -> Self {
-        self.with_runtime_options(|runtime| runtime.prior = Some(prior))
-    }
-
-    pub fn convergence(self, convergence: ConvergenceOptions) -> Self {
-        self.with_runtime_options(|runtime| runtime.convergence = convergence)
-    }
-
-    pub fn tuning(self, tuning: AlgorithmTuning) -> Self {
-        self.with_runtime_options(|runtime| runtime.tuning = tuning)
-    }
-}
-
-impl<E: Equation + EquationMetadataSource> EstimationProblemBuilder<E> {
-    pub fn parameter(self, parameter: Parameter) -> Result<Self> {
-        self.map_model_builder(|model| model.parameter(parameter))
-    }
-
-    pub fn variability(self, variability: VariabilityModel) -> Self {
-        self.with_model_builder(|model| model.variability(variability))
-    }
-
-    pub fn covariates(self, covariates: CovariateSpec) -> Self {
-        self.with_model_builder(|model| model.covariates(covariates))
-    }
-
-    pub fn metadata(self, metadata: ModelMetadata) -> Self {
-        self.with_model_builder(|model| model.metadata(metadata))
-    }
-
-    fn map_model_builder(
-        self,
-        map: impl FnOnce(ModelDefinitionBuilder<E>) -> Result<ModelDefinitionBuilder<E>>,
-    ) -> Result<Self> {
-        let EstimationProblemBuilder {
-            model,
-            data,
-            output,
-            runtime,
-        } = self;
-
-        Ok(Self {
-            model: map(model)?,
-            data,
-            output,
-            runtime,
-        })
-    }
-
-    fn with_model_builder(
-        self,
-        map: impl FnOnce(ModelDefinitionBuilder<E>) -> ModelDefinitionBuilder<E>,
-    ) -> Self {
-        let EstimationProblemBuilder {
-            model,
-            data,
-            output,
-            runtime,
-        } = self;
-
-        Self {
-            model: map(model),
-            data,
-            output,
-            runtime,
-        }
-    }
-}
-
-impl<E: Equation> EstimationProblemBuilder<E> {
-    fn with_output_plan(mut self, map: impl FnOnce(&mut OutputPlan)) -> Self {
-        map(self.output.get_or_insert_with(OutputPlan::default));
-        self
-    }
-
-    fn with_runtime_options(mut self, map: impl FnOnce(&mut RuntimeOptions)) -> Self {
-        map(self.runtime.get_or_insert_with(RuntimeOptions::default));
-        self
-    }
-}
-
-pub struct NonparametricEstimationProblemBuilder<E: Equation> {
-    builder: EstimationProblemBuilder<E>,
-    method: NonparametricMethod,
-    error_models: Option<AssayErrorModels>,
-}
-
-impl<E: Equation> NonparametricEstimationProblemBuilder<E> {
-    fn new(builder: EstimationProblemBuilder<E>, method: NonparametricMethod) -> Self {
-        Self {
-            builder,
-            method,
-            error_models: None,
-        }
-    }
-
-    pub fn cache(self, enabled: bool) -> Self {
-        self.with_builder(|builder| builder.cache(enabled))
-    }
-
-    pub fn progress(self, enabled: bool) -> Self {
-        self.with_builder(|builder| builder.progress(enabled))
-    }
-
-    pub fn idelta(self, value: f64) -> Self {
-        self.with_builder(|builder| builder.idelta(value))
-    }
-
-    pub fn tad(self, value: f64) -> Self {
-        self.with_builder(|builder| builder.tad(value))
-    }
-
-    pub fn prior(self, prior: Prior) -> Self {
-        self.with_builder(|builder| builder.prior(prior))
-    }
-
-    pub fn convergence(self, convergence: ConvergenceOptions) -> Self {
-        self.with_builder(|builder| builder.convergence(convergence))
-    }
-
-    pub fn tuning(self, tuning: AlgorithmTuning) -> Self {
-        self.with_builder(|builder| builder.tuning(tuning))
-    }
-
-    fn with_builder(
-        self,
-        map: impl FnOnce(EstimationProblemBuilder<E>) -> EstimationProblemBuilder<E>,
-    ) -> Self {
-        let NonparametricEstimationProblemBuilder {
-            builder,
-            method,
-            error_models,
-        } = self;
-
-        Self {
-            builder: map(builder),
-            method,
-            error_models,
-        }
-    }
-}
-
-impl<E: Equation + EquationMetadataSource> NonparametricEstimationProblemBuilder<E> {
-    pub fn parameter(self, parameter: Parameter) -> Result<Self> {
-        self.map_builder(|builder| builder.parameter(parameter))
-    }
-
-    pub fn variability(self, variability: VariabilityModel) -> Self {
-        self.with_builder(|builder| builder.variability(variability))
-    }
-
-    pub fn covariates(self, covariates: CovariateSpec) -> Self {
-        self.with_builder(|builder| builder.covariates(covariates))
-    }
-
-    pub fn metadata(self, metadata: ModelMetadata) -> Self {
-        self.with_builder(|builder| builder.metadata(metadata))
-    }
-
-    pub fn error(mut self, name: &str, model: AssayErrorModel) -> Result<Self> {
-        let outeq = self
-            .builder
-            .model
-            .output_index(name)
-            .ok_or_else(|| anyhow::anyhow!("unknown equation output label: {name}"))?;
-
-        self.error_models = Some(match self.error_models.take() {
-            None => AssayErrorModels::new().add(outeq, model)?,
-            Some(models) => models.add(outeq, model)?,
-        });
-
-        Ok(self)
-    }
-
-    pub fn build(self) -> Result<EstimationProblem<E>> {
-        let model = self.builder.model.build()?;
-        let error_models = self
-            .error_models
-            .ok_or_else(|| anyhow::anyhow!("error models are required"))?;
-
-        Ok(EstimationProblem {
-            model,
-            data: self.builder.data,
-            error_models: ErrorModels::Nonparametric(error_models),
-            method: self.method,
-            output: self.builder.output.unwrap_or_default(),
-            runtime: self.builder.runtime.unwrap_or_default(),
-        })
-    }
-
-    fn map_builder(
-        self,
-        map: impl FnOnce(EstimationProblemBuilder<E>) -> Result<EstimationProblemBuilder<E>>,
-    ) -> Result<Self> {
-        let NonparametricEstimationProblemBuilder {
-            builder,
-            method,
-            error_models,
-        } = self;
-
-        Ok(Self {
-            builder: map(builder)?,
-            method,
-            error_models,
-        })
-    }
-}
-
-impl<E: Equation + Clone + Send + 'static + EquationMetadataSource>
-    NonparametricEstimationProblemBuilder<E>
-{
-    pub fn fit(self) -> Result<crate::results::FitResult<E>> {
-        self.build()?.fit()
-    }
-
-    pub fn fit_with_progress<F>(self, on_progress: F) -> Result<crate::results::FitResult<E>>
-    where
-        F: FnMut(crate::api::FitProgress),
-    {
-        self.build()?.fit_with_progress(on_progress)
-    }
-}
diff --git a/src/api/fit.rs b/src/api/fit.rs
deleted file mode 100644
index e9f0a553b..000000000
--- a/src/api/fit.rs
+++ /dev/null
@@ -1,109 +0,0 @@
-use anyhow::Result;
-use pharmsol::equation::Equation;
-
-use crate::api::estimation_problem::EstimationProblem;
-use crate::api::progress::FitProgress;
-use crate::estimation::nonparametric;
-use crate::model::EquationMetadataSource;
-use crate::results::FitResult;
-
-pub fn fit<E: Equation + Clone + Send + 'static + EquationMetadataSource>(
-    problem: EstimationProblem<E>,
-) -> Result<FitResult<E>> {
-    let compiled = problem.compile()?;
-    let mut result = nonparametric::fit(compiled)?;
-    result.write_outputs()?;
-    Ok(result)
-}
-
-pub fn fit_with_progress<E, F>(
-    problem: EstimationProblem<E>,
-    mut on_progress: F,
-) -> Result<FitResult<E>>
-where
-    E: Equation + Clone + Send + 'static + EquationMetadataSource,
-    F: FnMut(FitProgress),
-{
-    let compiled = problem.compile()?;
-    let mut result = nonparametric::fit_with_progress(compiled, |event| {
-        on_progress(FitProgress::NonparametricCycle(event));
-    })?;
-    result.write_outputs()?;
-    Ok(result)
-}
-
-#[cfg(test)]
-mod tests {
-    use anyhow::Result;
-    use pharmsol::{AssayErrorModel, ErrorPoly, Subject};
-
-    use super::{fit_with_progress, FitProgress};
-    use crate::algorithms::Status;
-    use crate::api::{EstimationProblem, Npag};
-    use crate::prelude::*;
-
-    fn equation() -> equation::ODE {
-        equation::ODE::new(
-            |x, p, _t, dx, b, _rateiv, _cov| {
-                fetch_params!(p, ke);
-                dx[0] = -ke * x[0] + b[0];
-            },
-            |_p, _t, _cov| lag! {},
-            |_p, _t, _cov| fa! {},
-            |_p, _t, _cov, _x| {},
-            |x, p, _t, _cov, y| {
-                fetch_params!(p, v);
-                y[0] = x[0] / v;
-            },
-        )
-        .with_nstates(1)
-        .with_ndrugs(1)
-        .with_nout(1)
-        .with_metadata(
-            equation::metadata::new("fit_progress")
-                .parameters(["ke", "v"])
-                .states(["central"])
-                .outputs(["0"])
-                .route(equation::Route::bolus("0").to_state("central")),
-        )
-        .expect("metadata attachment should validate")
-    }
-
-    #[test]
-    fn fit_with_progress_reports_nonparametric_cycles() -> Result<()> {
-        let data = pharmsol::Data::new(vec![Subject::builder("1")
-            .bolus(0.0, 100.0, 0)
-            .observation(1.0, 10.0, 0)
-            .observation(2.0, 7.0, 0)
-            .build()]);
-
-        let assay_error = AssayErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0);
-
-        let problem = EstimationProblem::builder(equation(), data)
-            .parameter(Parameter::bounded("ke", 0.05, 1.0))?
-            .parameter(Parameter::bounded("v", 5.0, 50.0))?
-            .method(Npag::new())
-            .error("0", assay_error)?
-            .progress(false)
-            .prior(Prior::sobol(8, 7))
-            .build()?;
-
-        let mut progress_events = Vec::new();
-        let result = fit_with_progress(problem, |event| progress_events.push(event))?;
-
-        let workspace = result
-            .as_nonparametric()
-            .expect("expected nonparametric fit result");
-
-        assert!(!progress_events.is_empty());
-        assert_eq!(progress_events.len(), workspace.cycle_log().cycles().len());
-        assert_eq!(progress_events.len(), workspace.cycles());
-
-        let last_event = progress_events.last().cloned().expect("missing last event");
-        let FitProgress::NonparametricCycle(last_cycle) = last_event;
-        assert!(matches!(last_cycle.status, Status::Stop(_)));
-        assert_eq!(last_cycle.cycle, workspace.cycles());
-
-        Ok(())
-    }
-}
diff --git a/src/api/mod.rs b/src/api/mod.rs
deleted file mode 100644
index 811abc0c5..000000000
--- a/src/api/mod.rs
+++ /dev/null
@@ -1,16 +0,0 @@
-pub mod error_models;
-pub mod estimation_problem;
-pub mod fit;
-pub mod model_definition;
-pub mod progress;
-pub mod saem_config;
-
-pub use error_models::ErrorModels;
-pub use estimation_problem::{
-    AlgorithmTuning, ConvergenceOptions, EstimationProblem, EstimationProblemBuilder, MethodSpec,
-    NonparametricEstimationProblemBuilder, Npag, Npod, OutputPlan, PostProb, RuntimeOptions,
-};
-pub use fit::{fit, fit_with_progress};
-pub use model_definition::{ModelDefinition, ModelDefinitionBuilder};
-pub use progress::{FitProgress, NonparametricCycleProgress};
-pub use saem_config::SaemConfig;
diff --git a/src/api/model_definition.rs b/src/api/model_definition.rs
deleted file mode 100644
index f596d46c3..000000000
--- a/src/api/model_definition.rs
+++ /dev/null
@@ -1 +0,0 @@
-pub use crate::model::{ModelDefinition, ModelDefinitionBuilder};
diff --git a/src/bestdose/posterior.rs b/src/bestdose/posterior.rs
index 5e06e7dfe..3c7d30a7f 100644
--- a/src/bestdose/posterior.rs
+++ b/src/bestdose/posterior.rs
@@ -55,11 +55,11 @@ use faer::Mat;
 
 use crate::algorithms::nonparametric::npag::burke;
 use crate::algorithms::nonparametric::npag::NPAG;
-use crate::algorithms::Algorithms;
-use crate::algorithms::NativeNonparametricConfig;
+use crate::algorithms::NonParametricRunner;
+
 use crate::algorithms::Status;
 use crate::bestdose::types::BestDoseConfig;
-use crate::estimation::nonparametric::{calculate_psi, Prior, Theta, Weights};
+use crate::estimation::nonparametric::{calculate_psi, Theta, Weights};
 use crate::prelude::*;
 use pharmsol::prelude::*;
 
@@ -188,17 +188,6 @@ pub fn npagfull_refinement(
     let mut kept_weights: Vec<f64> = Vec::new();
     let num_points = filtered_theta.matrix().nrows();
     let parameter_space = config.parameter_space().clone();
-    let runtime = RuntimeOptions {
-        cycles: config.refinement_cycles(),
-        cache: true,
-        progress: config.progress(),
-        idelta: config.prediction_interval(),
-        tad: 0.0,
-        prior: None,
-
-        convergence: ConvergenceOptions::default(),
-        tuning: AlgorithmTuning::default(),
-    };
 
     for i in 0..num_points {
         tracing::debug!("  Refining point {}/{}", i + 1, num_points);
@@ -212,26 +201,19 @@ pub fn npagfull_refinement(
         let single_point_theta = Theta::from_parts(single_point_matrix, parameter_space.clone())?;
 
         // Create and run NPAG
-        let mut npag = NPAG::from_config(
+        let npag_config = NpagConfig {
+            max_cycles: config.refinement_cycles(),
+            progress: config.progress(),
+            ..Default::default()
+        };
+
+        let mut npag = NPAG::from_parts(
             eq.clone(),
             past_data.clone(),
             config.error_models().clone(),
-            NativeNonparametricConfig {
-                ranges: parameter_space.finite_ranges()?,
-                parameter_space: parameter_space.clone(),
-                prior: Prior::Theta(single_point_theta.clone()),
-                max_cycles: config.refinement_cycles(),
-                progress: config.progress(),
-                run_configuration: crate::output::shared::RunConfiguration::new(
-                    Algorithm::NPAG,
-                    &OutputPlan::disabled(),
-                    &runtime,
-                    config.parameter_names(),
-                ),
-            },
-            Npag::default(),
-        );
-        npag.set_theta(single_point_theta);
+            single_point_theta,
+            npag_config,
+        )?;
 
         // Run NPAG optimization
         let refinement_result = npag.initialize().and_then(|_| {
diff --git a/src/bestdose/types.rs b/src/bestdose/types.rs
index d759d9d7d..5c08de9c5 100644
--- a/src/bestdose/types.rs
+++ b/src/bestdose/types.rs
@@ -9,7 +9,8 @@
 
 use std::fmt::Display;
 
-use crate::estimation::nonparametric::{NPPredictions, Prior, Theta, Weights};
+use crate::estimation::nonparametric::{NPPredictions, Theta, Weights};
+use crate::model::{BoundedParameter, ParameterSpace};
 use crate::prelude::*;
 use pharmsol::prelude::*;
 use serde::{Deserialize, Serialize};
@@ -185,31 +186,24 @@ impl Default for DoseRange {
 
 #[derive(Debug, Clone)]
 pub struct BestDoseConfig {
-    pub(crate) parameter_space: ParameterSpace,
+    pub(crate) parameter_space: ParameterSpace<BoundedParameter>,
     pub(crate) error_models: AssayErrorModels,
-    pub(crate) prior: Prior,
     pub(crate) refinement_cycles: usize,
     pub(crate) progress: bool,
     pub(crate) prediction_interval: f64,
 }
 
 impl BestDoseConfig {
-    pub fn new(parameter_space: ParameterSpace, error_models: AssayErrorModels) -> Self {
+    pub fn new(parameter_space: ParameterSpace<BoundedParameter>, error_models: AssayErrorModels) -> Self {
         Self {
             parameter_space,
             error_models,
-            prior: Prior::default(),
             refinement_cycles: 500,
             progress: true,
             prediction_interval: 0.12,
         }
     }
 
-    pub fn with_prior(mut self, prior: Prior) -> Self {
-        self.prior = prior;
-        self
-    }
-
     pub fn with_refinement_cycles(mut self, refinement_cycles: usize) -> Self {
         self.refinement_cycles = refinement_cycles;
         self
@@ -225,7 +219,7 @@ impl BestDoseConfig {
         self
     }
 
-    pub fn parameter_space(&self) -> &ParameterSpace {
+    pub fn parameter_space(&self) -> &ParameterSpace<BoundedParameter> {
         &self.parameter_space
     }
 
@@ -233,10 +227,6 @@ impl BestDoseConfig {
         &self.error_models
     }
 
-    pub fn prior(&self) -> &Prior {
-        &self.prior
-    }
-
     pub fn refinement_cycles(&self) -> usize {
         self.refinement_cycles
     }
@@ -249,7 +239,7 @@ impl BestDoseConfig {
         self.prediction_interval
     }
 
-    pub(crate) fn parameter_names(&self) -> Vec<String> {
+    pub fn parameter_names(&self) -> Vec<String> {
         self.parameter_space
             .iter()
             .map(|parameter| parameter.name.clone())
diff --git a/src/compile/caches.rs b/src/compile/caches.rs
deleted file mode 100644
index 1ceeb2961..000000000
--- a/src/compile/caches.rs
+++ /dev/null
@@ -1,6 +0,0 @@
-use serde::{Deserialize, Serialize};
-
-#[derive(Debug, Clone, Default, Serialize, Deserialize, PartialEq, Eq)]
-pub struct ExecutionCaches {
-    pub prediction_cache_enabled: bool,
-}
diff --git a/src/compile/compiled_problem.rs b/src/compile/compiled_problem.rs
deleted file mode 100644
index ec010b449..000000000
--- a/src/compile/compiled_problem.rs
+++ /dev/null
@@ -1,72 +0,0 @@
-use pharmsol::equation::Equation;
-
-use crate::algorithms::Algorithm;
-use crate::api::estimation_problem::NonparametricMethod;
-use crate::api::{ErrorModels, OutputPlan, RuntimeOptions};
-use crate::compile::{DesignContext, ExecutionCaches, ObservationIndex};
-use crate::model::ModelDefinition;
-use pharmsol::Data;
-
-#[derive(Debug, Clone)]
-pub struct CompiledProblem<E: Equation> {
-    pub model: ModelDefinition<E>,
-    pub data: Data,
-    error_models: ErrorModels,
-    method: NonparametricMethod,
-    output: OutputPlan,
-    runtime: RuntimeOptions,
-    pub design: DesignContext,
-    pub observation_index: ObservationIndex,
-    pub caches: ExecutionCaches,
-}
-
-impl<E: Equation> CompiledProblem<E> {
-    #[allow(clippy::too_many_arguments)]
-    pub(crate) fn new(
-        model: ModelDefinition<E>,
-        data: Data,
-        error_models: ErrorModels,
-        method: NonparametricMethod,
-        output: OutputPlan,
-        runtime: RuntimeOptions,
-        design: DesignContext,
-        observation_index: ObservationIndex,
-        caches: ExecutionCaches,
-    ) -> Self {
-        Self {
-            model,
-            data,
-            error_models,
-            method,
-            output,
-            runtime,
-            design,
-            observation_index,
-            caches,
-        }
-    }
-
-    pub(crate) fn method(&self) -> NonparametricMethod {
-        self.method
-    }
-
-    pub fn algorithm(&self) -> Algorithm {
-        self.method.algorithm()
-    }
-
-    pub fn error_models(&self) -> &ErrorModels {
-        &self.error_models
-    }
-
-    pub fn output_plan(&self) -> &OutputPlan {
-        &self.output
-    }
-
-    pub fn runtime_options(&self) -> &RuntimeOptions {
-        &self.runtime
-    }
-
-    pub fn into_parts(self) -> (ModelDefinition<E>, Data) {
-        (self.model, self.data)
-    }
-}
diff --git a/src/compile/design_context.rs b/src/compile/design_context.rs
deleted file mode 100644
index 9927e1b7f..000000000
--- a/src/compile/design_context.rs
+++ /dev/null
@@ -1,59 +0,0 @@
-use serde::{Deserialize, Serialize};
-
-#[derive(Debug, Clone, Serialize, Deserialize, PartialEq)]
-pub struct DesignContext {
-    pub parameter_names: Vec<String>,
-    pub subjects: Vec<SubjectDesign>,
-    pub occasions: Vec<OccasionDesign>,
-    pub structured_covariates: StructuredCovariateDesign,
-}
-
-impl DesignContext {
-    pub fn subject_count(&self) -> usize {
-        self.subjects.len()
-    }
-
-    pub fn occasion_count(&self) -> usize {
-        self.occasions.len()
-    }
-}
-
-#[derive(Debug, Clone, Serialize, Deserialize, PartialEq, Eq)]
-pub struct SubjectDesign {
-    pub subject_index: usize,
-    pub id: String,
-    pub occasion_count: usize,
-    pub observation_count: usize,
-}
-
-#[derive(Debug, Clone, Serialize, Deserialize, PartialEq, Eq)]
-pub struct OccasionDesign {
-    pub subject_index: usize,
-    pub occasion_index: usize,
-    pub event_count: usize,
-    pub observation_count: usize,
-}
-
-#[derive(Debug, Clone, Default, Serialize, Deserialize, PartialEq)]
-pub struct StructuredCovariateDesign {
-    pub subject_columns: Vec<String>,
-    pub subject_rows: Vec<SubjectCovariateRow>,
-    pub occasion_columns: Vec<String>,
-    pub occasion_rows: Vec<OccasionCovariateRow>,
-}
-
-#[derive(Debug, Clone, Serialize, Deserialize, PartialEq)]
-pub struct SubjectCovariateRow {
-    pub subject_index: usize,
-    pub id: String,
-    pub anchor_time: f64,
-    pub values: Vec<Option<f64>>,
-}
-
-#[derive(Debug, Clone, Serialize, Deserialize, PartialEq)]
-pub struct OccasionCovariateRow {
-    pub subject_index: usize,
-    pub occasion_index: usize,
-    pub anchor_time: f64,
-    pub values: Vec<Option<f64>>,
-}
diff --git a/src/compile/mod.rs b/src/compile/mod.rs
deleted file mode 100644
index 9484b559f..000000000
--- a/src/compile/mod.rs
+++ /dev/null
@@ -1,251 +0,0 @@
-use anyhow::Result;
-use pharmsol::{Data, Equation, Event};
-
-use crate::api::EstimationProblem;
-use crate::model::{CovariateSpec, EquationMetadataSource, ParameterSpace};
-
-mod caches;
-mod compiled_problem;
-mod design_context;
-mod observation_index;
-mod validation;
-
-pub use caches::ExecutionCaches;
-pub use compiled_problem::CompiledProblem;
-pub use design_context::{
-    DesignContext, OccasionCovariateRow, OccasionDesign, StructuredCovariateDesign,
-    SubjectCovariateRow, SubjectDesign,
-};
-pub use observation_index::{ObservationIndex, ObservationRecord};
-pub use validation::validate_problem;
-
-impl<E: Equation + Clone + EquationMetadataSource> EstimationProblem<E> {
-    pub fn compile(self) -> Result<CompiledProblem<E>> {
-        compile_problem(self)
-    }
-}
-
-pub fn compile_problem<E: Equation + Clone + EquationMetadataSource>(
-    problem: EstimationProblem<E>,
-) -> Result<CompiledProblem<E>> {
-    validate_problem(&problem)?;
-
-    let design = build_design_context(
-        &problem.model.parameters,
-        &problem.model.covariates,
-        &problem.data,
-    );
-    let observation_index = build_observation_index(&problem.data)?;
-    let caches = ExecutionCaches {
-        prediction_cache_enabled: problem.runtime.cache,
-    };
-
-    Ok(CompiledProblem::new(
-        problem.model,
-        problem.data,
-        problem.error_models,
-        problem.method,
-        problem.output,
-        problem.runtime,
-        design,
-        observation_index,
-        caches,
-    ))
-}
-
-fn build_design_context(
-    parameter_space: &ParameterSpace,
-    covariates: &CovariateSpec,
-    data: &Data,
-) -> DesignContext {
-    let subjects = data.subjects();
-
-    let subject_design = subjects
-        .iter()
-        .enumerate()
-        .map(|(subject_index, subject)| {
-            let occasions = subject.occasions();
-            let observation_count = occasions
-                .iter()
-                .map(|occasion| {
-                    occasion
-                        .events()
-                        .iter()
-                        .filter(|event| matches!(event, Event::Observation(_)))
-                        .count()
-                })
-                .sum();
-
-            SubjectDesign {
-                subject_index,
-                id: subject.id().clone(),
-                occasion_count: occasions.len(),
-                observation_count,
-            }
-        })
-        .collect::<Vec<_>>();
-
-    let occasion_design = subjects
-        .iter()
-        .enumerate()
-        .flat_map(|(subject_index, subject)| {
-            subject.occasions().iter().map(move |occasion| {
-                let events = occasion.events();
-                let observation_count = events
-                    .iter()
-                    .filter(|event| matches!(event, Event::Observation(_)))
-                    .count();
-
-                OccasionDesign {
-                    subject_index,
-                    occasion_index: occasion.index(),
-                    event_count: events.len(),
-                    observation_count,
-                }
-            })
-        })
-        .collect::<Vec<_>>();
-
-    let structured_covariates = match covariates {
-        CovariateSpec::InEquation => StructuredCovariateDesign::default(),
-        CovariateSpec::Structured(spec) => build_structured_covariate_design(
-            &spec.subject_columns(),
-            &spec.occasion_columns(),
-            data,
-        ),
-    };
-
-    DesignContext {
-        parameter_names: parameter_space
-            .iter()
-            .map(|item| item.name.clone())
-            .collect(),
-        subjects: subject_design,
-        occasions: occasion_design,
-        structured_covariates,
-    }
-}
-
-fn build_structured_covariate_design(
-    subject_columns: &[String],
-    occasion_columns: &[String],
-    data: &Data,
-) -> StructuredCovariateDesign {
-    let subject_rows = data
-        .subjects()
-        .iter()
-        .enumerate()
-        .map(|(subject_index, subject)| {
-            let anchor_time = subject_anchor_time(subject);
-            let values = subject_columns
-                .iter()
-                .map(|name| subject_covariate_value(subject, name))
-                .collect();
-
-            SubjectCovariateRow {
-                subject_index,
-                id: subject.id().clone(),
-                anchor_time,
-                values,
-            }
-        })
-        .collect();
-
-    let occasion_rows = data
-        .subjects()
-        .iter()
-        .enumerate()
-        .flat_map(|(subject_index, subject)| {
-            subject.occasions().iter().map(move |occasion| {
-                let anchor_time = occasion_anchor_time(occasion);
-                let values = occasion_columns
-                    .iter()
-                    .map(|name| {
-                        occasion
-                            .covariates()
-                            .get_covariate(name)
-                            .and_then(|covariate| covariate.interpolate(anchor_time).ok())
-                    })
-                    .collect();
-
-                OccasionCovariateRow {
-                    subject_index,
-                    occasion_index: occasion.index(),
-                    anchor_time,
-                    values,
-                }
-            })
-        })
-        .collect();
-
-    StructuredCovariateDesign {
-        subject_columns: subject_columns.to_vec(),
-        subject_rows,
-        occasion_columns: occasion_columns.to_vec(),
-        occasion_rows,
-    }
-}
-
-fn subject_anchor_time(subject: &pharmsol::Subject) -> f64 {
-    subject
-        .occasions()
-        .iter()
-        .find_map(|occasion| occasion.events().first().map(|event| event.time()))
-        .unwrap_or(0.0)
-}
-
-fn subject_covariate_value(subject: &pharmsol::Subject, name: &str) -> Option<f64> {
-    subject.occasions().iter().find_map(|occasion| {
-        let anchor_time = occasion_anchor_time(occasion);
-        occasion
-            .covariates()
-            .get_covariate(name)
-            .and_then(|covariate| covariate.interpolate(anchor_time).ok())
-    })
-}
-
-fn occasion_anchor_time(occasion: &pharmsol::Occasion) -> f64 {
-    occasion
-        .events()
-        .first()
-        .map(|event| event.time())
-        .unwrap_or(0.0)
-}
-
-fn build_observation_index(data: &Data) -> Result<ObservationIndex> {
-    let records =
-        data.subjects()
-            .iter()
-            .enumerate()
-            .flat_map(|(subject_index, subject)| {
-                subject.occasions().iter().flat_map(move |occasion| {
-                    occasion
-                        .events()
-                        .iter()
-                        .enumerate()
-                        .filter_map(move |(event_index, event)| match event {
-                            Event::Observation(observation) => Some(
-                                observation
-                                    .outeq_index()
-                                    .map(|outeq| ObservationRecord {
-                                        subject_index,
-                                        occasion_index: occasion.index(),
-                                        event_index,
-                                        outeq,
-                                        time: observation.time(),
-                                    })
-                                    .ok_or_else(|| {
-                                        anyhow::anyhow!(
-                                            "Compilation requires numeric observation output labels; got `{}`",
-                                            observation.outeq()
-                                        )
-                                    }),
-                            ),
-                            _ => None,
-                        })
-                })
-            })
-            .collect::<Result<Vec<_>>>()?;
-
-    Ok(ObservationIndex { records })
-}
diff --git a/src/compile/observation_index.rs b/src/compile/observation_index.rs
deleted file mode 100644
index 7cd989de9..000000000
--- a/src/compile/observation_index.rs
+++ /dev/null
@@ -1,25 +0,0 @@
-use serde::{Deserialize, Serialize};
-
-#[derive(Debug, Clone, Default, Serialize, Deserialize, PartialEq)]
-pub struct ObservationIndex {
-    pub records: Vec<ObservationRecord>,
-}
-
-impl ObservationIndex {
-    pub fn len(&self) -> usize {
-        self.records.len()
-    }
-
-    pub fn is_empty(&self) -> bool {
-        self.records.is_empty()
-    }
-}
-
-#[derive(Debug, Clone, Serialize, Deserialize, PartialEq)]
-pub struct ObservationRecord {
-    pub subject_index: usize,
-    pub occasion_index: usize,
-    pub event_index: usize,
-    pub outeq: usize,
-    pub time: f64,
-}
diff --git a/src/compile/validation.rs b/src/compile/validation.rs
deleted file mode 100644
index 8d145d238..000000000
--- a/src/compile/validation.rs
+++ /dev/null
@@ -1,30 +0,0 @@
-use anyhow::{bail, Result};
-use pharmsol::equation::Equation;
-
-use crate::api::EstimationProblem;
-use crate::model::EquationMetadataSource;
-
-pub fn validate_problem<E: Equation + EquationMetadataSource>(
-    problem: &EstimationProblem<E>,
-) -> Result<()> {
-    if problem.model.parameters.is_empty() {
-        bail!("estimation problem requires at least one parameter");
-    }
-
-    if problem.runtime.cycles == 0 {
-        bail!("runtime cycles must be greater than zero");
-    }
-
-    if problem.model.output_count() == 0 {
-        bail!("at least one equation output is required");
-    }
-
-    let error_models = problem.error_models.models();
-    if error_models.iter().next().is_none() {
-        bail!("at least one nonparametric error model is required");
-    }
-
-    problem.model.parameters.finite_ranges()?;
-
-    Ok(())
-}
diff --git a/src/api/error_models.rs b/src/estimation/error_models.rs
similarity index 80%
rename from src/api/error_models.rs
rename to src/estimation/error_models.rs
index 20e9767c8..4dcc27b5b 100644
--- a/src/api/error_models.rs
+++ b/src/estimation/error_models.rs
@@ -18,3 +18,9 @@ impl ErrorModels {
         self.models().iter().next().is_none()
     }
 }
+
+impl Into<AssayErrorModels> for ErrorModels {
+    fn into(self) -> AssayErrorModels {
+        self.models().clone()
+    }
+}
diff --git a/src/estimation/mod.rs b/src/estimation/mod.rs
index d97ede1fc..2bf9f0104 100644
--- a/src/estimation/mod.rs
+++ b/src/estimation/mod.rs
@@ -1 +1,14 @@
+pub mod error_models;
 pub mod nonparametric;
+pub mod problem;
+pub mod progress;
+
+pub use crate::algorithms::nonparametric::{
+    NonParametricAlgorithm, NpagConfig, NpmapConfig, NpodConfig,
+};
+pub use crate::algorithms::parametric::{ParametricAlgorithm, SaemConfig};
+pub use error_models::ErrorModels;
+pub use problem::{
+    EstimationProblem, Framework, NonParametric, Parametric, ParametricBuilder,
+};
+pub use progress::{FitProgress, NonparametricCycleProgress};
diff --git a/src/estimation/nonparametric/cycles.rs b/src/estimation/nonparametric/cycles.rs
index 6a60e4fd6..8115a9bff 100644
--- a/src/estimation/nonparametric/cycles.rs
+++ b/src/estimation/nonparametric/cycles.rs
@@ -1,3 +1,5 @@
+use std::{fs::File, path::Path};
+
 use anyhow::Result;
 use csv::WriterBuilder;
 use pharmsol::{AssayErrorModel, AssayErrorModels};
@@ -7,7 +9,6 @@ use crate::{
     algorithms::{Status, StopReason},
     estimation::nonparametric::median,
     estimation::nonparametric::theta::Theta,
-    output::OutputFile,
 };
 
 #[derive(Debug, Clone, Serialize)]
@@ -95,12 +96,12 @@ impl CycleLog {
         self.cycles.push(cycle);
     }
 
-    pub fn write(&self, folder: &str, parameter_names: &[String]) -> Result<()> {
+    pub fn write(&self, path: &Path) -> Result<()> {
         tracing::debug!("Writing cycles...");
-        let outputfile = OutputFile::new(folder, "iterations.csv")?;
-        let mut writer = WriterBuilder::new()
-            .has_headers(false)
-            .from_writer(outputfile.file());
+
+        std::fs::create_dir_all(path)?;
+        let file = File::create(path)?;
+        let mut writer = WriterBuilder::new().has_headers(true).from_writer(file);
 
         writer.write_field("cycle")?;
         writer.write_field("converged")?;
@@ -124,7 +125,13 @@ impl CycleLog {
             )?;
         }
 
-        for param_name in parameter_names {
+        let names = self
+            .cycles
+            .first()
+            .map(|cycle| cycle.theta.param_names())
+            .expect("No cycles");
+
+        for param_name in names {
             writer.write_field(format!("{}.mean", param_name))?;
             writer.write_field(format!("{}.median", param_name))?;
             writer.write_field(format!("{}.sd", param_name))?;
@@ -174,7 +181,8 @@ impl CycleLog {
             writer.write_record(None::<&[u8]>)?;
         }
         writer.flush()?;
-        tracing::debug!("Cycles written to {:?}", &outputfile.relative_path());
+
+        tracing::debug!("Cycles written to {:?}", path);
         Ok(())
     }
 }
diff --git a/src/estimation/nonparametric/engine.rs b/src/estimation/nonparametric/engine.rs
deleted file mode 100644
index 3f07e8607..000000000
--- a/src/estimation/nonparametric/engine.rs
+++ /dev/null
@@ -1,80 +0,0 @@
-use anyhow::Result;
-use pharmsol::Equation;
-
-use crate::algorithms::{
-    run_nonparametric_algorithm, run_nonparametric_algorithm_with_progress,
-    NonparametricAlgorithmInput,
-};
-use crate::api::NonparametricCycleProgress;
-use crate::compile::CompiledProblem;
-use crate::estimation::nonparametric::workspace::NonparametricWorkspace;
-use crate::results::FitResult;
-
-#[derive(Debug, Default, Clone, Copy)]
-pub struct NonparametricEngine;
-
-impl NonparametricEngine {
-    pub fn fit<E: Equation + Clone + Send + 'static>(
-        problem: CompiledProblem<E>,
-    ) -> Result<NonparametricWorkspace<E>> {
-        let input = input_from_compiled_problem(problem)?;
-        run_nonparametric_algorithm(input)
-    }
-
-    pub fn fit_with_progress<E, F>(
-        problem: CompiledProblem<E>,
-        mut on_progress: F,
-    ) -> Result<NonparametricWorkspace<E>>
-    where
-        E: Equation + Clone + Send + 'static,
-        F: FnMut(NonparametricCycleProgress),
-    {
-        let input = input_from_compiled_problem(problem)?;
-        run_nonparametric_algorithm_with_progress(
-            input,
-            |cycle, objective, objective_delta, elapsed_ms, status| {
-                on_progress(NonparametricCycleProgress {
-                    cycle,
-                    objective,
-                    objective_delta,
-                    elapsed_ms,
-                    status,
-                });
-            },
-        )
-    }
-}
-
-fn input_from_compiled_problem<E: Equation + Clone + Send + 'static>(
-    problem: CompiledProblem<E>,
-) -> Result<NonparametricAlgorithmInput<E>> {
-    let method = problem.method();
-    let error_models = problem.error_models().models().clone();
-    let output = problem.output_plan().clone();
-    let runtime = problem.runtime_options().clone();
-    let (model, data) = problem.into_parts();
-    Ok(NonparametricAlgorithmInput::new(
-        method,
-        model,
-        data,
-        error_models,
-        output,
-        runtime,
-    ))
-}
-
-pub fn fit<E: Equation + Clone + Send + 'static>(
-    problem: CompiledProblem<E>,
-) -> Result<FitResult<E>> {
-    let workspace = NonparametricEngine::fit(problem)?;
-    Ok(workspace.into_fit_result())
-}
-
-pub fn fit_with_progress<E, F>(problem: CompiledProblem<E>, on_progress: F) -> Result<FitResult<E>>
-where
-    E: Equation + Clone + Send + 'static,
-    F: FnMut(NonparametricCycleProgress),
-{
-    let workspace = NonparametricEngine::fit_with_progress(problem, on_progress)?;
-    Ok(workspace.into_fit_result())
-}
diff --git a/src/estimation/nonparametric/expansion.rs b/src/estimation/nonparametric/expansion.rs
index 63d058a1f..7ca892ae0 100644
--- a/src/estimation/nonparametric/expansion.rs
+++ b/src/estimation/nonparametric/expansion.rs
@@ -45,14 +45,14 @@ pub fn adaptative_grid(
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::model::{Parameter, ParameterSpace};
+    use crate::model::{BoundedParameter, ParameterSpace};
     use faer::mat;
 
     #[test]
     fn adaptive_grid_expands_points_within_bounds() {
-        let parameters = ParameterSpace::new()
-            .add(Parameter::bounded("x", 0.0, 1.0))
-            .add(Parameter::bounded("y", 0.0, 1.0));
+        let parameters = ParameterSpace::<BoundedParameter>::new()
+            .add("x", 0.0, 1.0)
+            .add("y", 0.0, 1.0);
         let mut theta = Theta::from_parts(mat![[0.5, 0.5]], parameters).unwrap();
         let ranges = [(0.0, 1.0), (0.0, 1.0)];
 
diff --git a/src/estimation/nonparametric/mod.rs b/src/estimation/nonparametric/mod.rs
index b74278d98..c3cdb91cb 100644
--- a/src/estimation/nonparametric/mod.rs
+++ b/src/estimation/nonparametric/mod.rs
@@ -1,30 +1,28 @@
 mod cycles;
-mod engine;
+
 mod expansion;
 pub(crate) mod ipm;
 mod posterior;
 mod predictions;
-mod prior;
+
+pub mod sampling;
 mod psi;
 pub(crate) mod qr;
+mod result;
 mod statistics;
 mod summaries;
 mod theta;
 mod weights;
-mod workspace;
 
 pub use cycles::{CycleLog, NPCycle};
-pub use engine::{fit, fit_with_progress, NonparametricEngine};
 pub(crate) use expansion::adaptative_grid;
 pub use ipm::burke;
 pub use posterior::{posterior, Posterior};
 pub use predictions::{NPPredictionRow, NPPredictions};
-pub(crate) use prior::sample_space_for_parameters;
-pub use prior::{read_prior, Prior};
 pub(crate) use psi::calculate_psi;
 pub use psi::Psi;
+pub use result::NonParametricResult;
 pub use statistics::{median, population_mean_median, posterior_mean_median, weighted_median};
 pub use summaries::{fit_summary, individual_summaries, population_summary};
 pub use theta::Theta;
 pub use weights::Weights;
-pub use workspace::NonparametricWorkspace;
diff --git a/src/estimation/nonparametric/prior.rs b/src/estimation/nonparametric/prior.rs
deleted file mode 100644
index 8e0488f2e..000000000
--- a/src/estimation/nonparametric/prior.rs
+++ /dev/null
@@ -1,325 +0,0 @@
-use std::fs::File;
-
-use crate::estimation::nonparametric::{Theta, Weights};
-use crate::model::{ParameterDomain, ParameterSpace};
-use anyhow::{bail, Context, Result};
-use faer::Mat;
-use serde::{Deserialize, Serialize};
-
-pub mod latin;
-pub mod sobol;
-
-#[derive(Debug, Clone)]
-pub enum Prior {
-    Sobol(usize, usize),
-    Latin(usize, usize),
-    File(String),
-    Theta(Theta),
-}
-
-impl Serialize for Prior {
-    fn serialize<S>(&self, serializer: S) -> Result<S::Ok, S::Error>
-    where
-        S: serde::Serializer,
-    {
-        #[derive(Serialize)]
-        #[serde(tag = "kind", rename_all = "snake_case")]
-        enum Ref<'a> {
-            Sobol { points: usize, seed: usize },
-            Latin { points: usize, seed: usize },
-            File { path: &'a str },
-        }
-        let wire = match self {
-            Prior::Sobol(points, seed) => Ref::Sobol {
-                points: *points,
-                seed: *seed,
-            },
-            Prior::Latin(points, seed) => Ref::Latin {
-                points: *points,
-                seed: *seed,
-            },
-            Prior::File(path) => Ref::File { path },
-            Prior::Theta(_) => {
-                return Err(serde::ser::Error::custom(
-                    "Prior::Theta cannot be serialized",
-                ))
-            }
-        };
-        wire.serialize(serializer)
-    }
-}
-
-impl<'de> Deserialize<'de> for Prior {
-    fn deserialize<D>(deserializer: D) -> Result<Self, D::Error>
-    where
-        D: serde::Deserializer<'de>,
-    {
-        #[derive(Deserialize)]
-        #[serde(tag = "kind", rename_all = "snake_case")]
-        enum Wire {
-            Sobol { points: usize, seed: usize },
-            Latin { points: usize, seed: usize },
-            File { path: String },
-        }
-        match Wire::deserialize(deserializer)? {
-            Wire::Sobol { points, seed } => Ok(Prior::Sobol(points, seed)),
-            Wire::Latin { points, seed } => Ok(Prior::Latin(points, seed)),
-            Wire::File { path } => Ok(Prior::File(path)),
-        }
-    }
-}
-
-impl Prior {
-    pub fn sobol(points: usize, seed: usize) -> Prior {
-        Prior::Sobol(points, seed)
-    }
-
-    pub fn points(&self) -> Option<usize> {
-        match self {
-            Prior::Sobol(points, _) => Some(*points),
-            Prior::Latin(points, _) => Some(*points),
-            Prior::File(_) => None,
-            Prior::Theta(theta) => Some(theta.nspp()),
-        }
-    }
-
-    pub fn seed(&self) -> Option<usize> {
-        match self {
-            Prior::Sobol(_, seed) => Some(*seed),
-            Prior::Latin(_, seed) => Some(*seed),
-            Prior::File(_) => None,
-            Prior::Theta(_) => None,
-        }
-    }
-}
-
-impl Default for Prior {
-    fn default() -> Self {
-        Prior::Sobol(2028, 22)
-    }
-}
-
-pub fn read_prior(
-    path: impl AsRef<str>,
-    parameters: impl Into<ParameterSpace>,
-) -> Result<(Theta, Option<Weights>)> {
-    let path = path.as_ref().to_string();
-    parse_prior_for_parameters(&path, parameters)
-}
-
-pub(crate) fn sample_space_for_parameters(
-    parameters: impl Into<ParameterSpace>,
-    prior: &Prior,
-) -> Result<Theta> {
-    let parameter_space = parameters.into();
-
-    for parameter in parameter_space.iter() {
-        let (lower, upper) = match parameter.domain {
-            ParameterDomain::Bounded { lower, upper } => (lower, upper),
-            ParameterDomain::Positive {
-                lower: Some(lower),
-                upper: Some(upper),
-            }
-            | ParameterDomain::Unbounded {
-                lower: Some(lower),
-                upper: Some(upper),
-            } => (lower, upper),
-            _ => bail!(
-                "Parameter '{}' is missing finite bounds required for nonparametric initialization",
-                parameter.name
-            ),
-        };
-
-        if lower.is_infinite() || upper.is_infinite() {
-            bail!(
-                "Parameter '{}' has infinite bounds: [{}, {}]",
-                parameter.name,
-                lower,
-                upper
-            );
-        }
-
-        if lower >= upper {
-            bail!(
-                "Parameter '{}' has invalid bounds: [{}, {}]. Lower bound must be less than upper bound.",
-                parameter.name,
-                lower,
-                upper
-            );
-        }
-    }
-
-    let prior = match prior {
-        Prior::Sobol(points, seed) => sobol::generate(&parameter_space, *points, *seed)?,
-        Prior::Latin(points, seed) => latin::generate(&parameter_space, *points, *seed)?,
-        Prior::File(path) => parse_prior_for_parameters(path, &parameter_space)?.0,
-        Prior::Theta(theta) => return Ok(theta.clone()),
-    };
-    Ok(prior)
-}
-
-pub(crate) fn parse_prior_for_parameters(
-    path: &String,
-    parameters: impl Into<ParameterSpace>,
-) -> Result<(Theta, Option<Weights>)> {
-    let parameters = parameters.into();
-    tracing::info!("Reading prior from {}", path);
-    let file = File::open(path).context(format!("Unable to open the prior file '{}'", path))?;
-    let mut reader = csv::ReaderBuilder::new()
-        .has_headers(true)
-        .from_reader(file);
-
-    let mut parameter_names: Vec<String> = reader
-        .headers()?
-        .clone()
-        .into_iter()
-        .map(|s| s.trim().to_owned())
-        .collect();
-
-    let prob_index = parameter_names.iter().position(|name| name == "prob");
-    if let Some(index) = prob_index {
-        parameter_names.remove(index);
-    }
-
-    let random_names: Vec<String> = parameters.names();
-
-    let mut reordered_indices: Vec<usize> = Vec::new();
-    for random_name in &random_names {
-        match parameter_names.iter().position(|name| name == random_name) {
-            Some(index) => {
-                let adjusted_index = if let Some(prob_idx) = prob_index {
-                    if index >= prob_idx {
-                        index + 1
-                    } else {
-                        index
-                    }
-                } else {
-                    index
-                };
-                reordered_indices.push(adjusted_index);
-            }
-            None => bail!("Parameter {} is not present in the CSV file.", random_name),
-        }
-    }
-
-    if parameter_names.len() > random_names.len() {
-        let extra_parameters: Vec<&String> = parameter_names.iter().collect();
-        bail!(
-            "Found parameters in the prior not present in configuration: {:?}",
-            extra_parameters
-        );
-    }
-
-    let mut theta_values = Vec::new();
-    let mut prob_values = Vec::new();
-
-    for result in reader.records() {
-        let record = result.unwrap();
-        let values: Vec<f64> = reordered_indices
-            .iter()
-            .map(|&i| record[i].parse::<f64>().unwrap())
-            .collect();
-        theta_values.push(values);
-
-        if let Some(prob_idx) = prob_index {
-            let prob_value: f64 = record[prob_idx].parse::<f64>().unwrap();
-            prob_values.push(prob_value);
-        }
-    }
-
-    let n_points = theta_values.len();
-    let n_params = random_names.len();
-    let theta_values: Vec<f64> = theta_values.into_iter().flatten().collect();
-    let theta_matrix: Mat<f64> =
-        Mat::from_fn(n_points, n_params, |i, j| theta_values[i * n_params + j]);
-
-    let theta = Theta::from_parts(theta_matrix, parameters.clone())?;
-    let weights = if !prob_values.is_empty() {
-        Some(Weights::from_vec(prob_values))
-    } else {
-        None
-    };
-
-    Ok((theta, weights))
-}
-
-#[cfg(test)]
-mod tests {
-    use super::*;
-    use crate::model::{Parameter, ParameterSpace};
-    use std::fs;
-
-    fn parameter_space() -> ParameterSpace {
-        ParameterSpace::new()
-            .add(Parameter::bounded("ke", 0.1, 1.0))
-            .add(Parameter::bounded("v", 5.0, 50.0))
-    }
-
-    fn temp_csv_path() -> String {
-        format!("test_temp_prior_{}.csv", rand::random::<u32>())
-    }
-
-    #[test]
-    fn prior_metadata_accessors() {
-        let sobol = Prior::sobol(100, 42);
-        assert_eq!(sobol.points(), Some(100));
-        assert_eq!(sobol.seed(), Some(42));
-
-        let latin = Prior::Latin(50, 7);
-        assert_eq!(latin.points(), Some(50));
-        assert_eq!(latin.seed(), Some(7));
-
-        let file = Prior::File("prior.csv".to_string());
-        assert_eq!(file.points(), None);
-        assert_eq!(file.seed(), None);
-    }
-
-    #[test]
-    fn sample_space_generates_expected_shape() {
-        let theta = sample_space_for_parameters(parameter_space(), &Prior::sobol(10, 42)).unwrap();
-        assert_eq!(theta.nspp(), 10);
-        assert_eq!(theta.matrix().ncols(), 2);
-    }
-
-    #[test]
-    fn sample_space_returns_custom_theta_verbatim() {
-        let parameters = parameter_space();
-        let matrix = Mat::from_fn(3, 2, |i, j| (i + j) as f64);
-        let custom = Theta::from_parts(matrix, parameters).unwrap();
-
-        let theta =
-            sample_space_for_parameters(parameter_space(), &Prior::Theta(custom.clone())).unwrap();
-        assert_eq!(theta.matrix(), custom.matrix());
-    }
-
-    #[test]
-    fn read_prior_parses_weights_and_reorders_columns() {
-        let path = temp_csv_path();
-        fs::write(&path, "v,ke,prob\n10.0,0.5,0.3\n15.0,0.7,0.7\n").unwrap();
-
-        let (theta, weights) = read_prior(&path, parameter_space()).unwrap();
-        let _ = fs::remove_file(&path);
-
-        assert_eq!(theta.nspp(), 2);
-        assert_eq!(theta.matrix()[(0, 0)], 0.5);
-        assert_eq!(theta.matrix()[(0, 1)], 10.0);
-
-        let weights = weights.expect("weights should be parsed from prob column");
-        assert_eq!(weights.len(), 2);
-        assert_eq!(weights[0], 0.3);
-        assert_eq!(weights[1], 0.7);
-    }
-
-    #[test]
-    fn read_prior_rejects_extra_parameters() {
-        let path = temp_csv_path();
-        fs::write(&path, "ke,v,extra\n0.5,10.0,1.0\n").unwrap();
-
-        let err = read_prior(&path, parameter_space()).unwrap_err();
-        let _ = fs::remove_file(&path);
-
-        assert!(err
-            .to_string()
-            .contains("Found parameters in the prior not present in configuration"));
-    }
-}
diff --git a/src/estimation/nonparametric/workspace.rs b/src/estimation/nonparametric/result.rs
similarity index 67%
rename from src/estimation/nonparametric/workspace.rs
rename to src/estimation/nonparametric/result.rs
index 16463fa67..96cc75bd3 100644
--- a/src/estimation/nonparametric/workspace.rs
+++ b/src/estimation/nonparametric/result.rs
@@ -1,58 +1,55 @@
+use std::path::Path;
+
 use pharmsol::Equation;
 
 use crate::algorithms::{Status, StopReason};
-use crate::estimation::nonparametric::{
-    posterior, CycleLog, NPPredictions, Posterior, Psi, Theta, Weights,
-};
-use crate::output::shared::RunConfiguration;
-use crate::results::FitResult;
-use pharmsol::Data;
+use crate::estimation::nonparametric::{CycleLog, NPPredictions, Posterior, Psi, Theta, Weights};
+
+use pharmsol::{AssayErrorModels, Data};
 
+/// Contains the results of a nonparametric estimation, including the final parameter
 #[derive(Debug)]
-pub struct NonparametricWorkspace<E: Equation> {
+pub struct NonParametricResult<E: Equation> {
     equation: E,
     data: Data,
+    error_models: AssayErrorModels,
+    prior: Theta,
     theta: Theta,
     psi: Psi,
     weights: Weights,
     objf: f64,
     cycles: usize,
     status: Status,
-    run_configuration: RunConfiguration,
     cyclelog: CycleLog,
-    predictions: Option<NPPredictions>,
-    posterior: Posterior,
 }
 
-impl<E: Equation> NonparametricWorkspace<E> {
+impl<E: Equation> NonParametricResult<E> {
     #[allow(clippy::too_many_arguments)]
     pub(crate) fn new(
         equation: E,
         data: Data,
+        error_models: AssayErrorModels,
+        prior: Theta,
         theta: Theta,
         psi: Psi,
         weights: Weights,
         objf: f64,
         cycles: usize,
         status: Status,
-        run_configuration: RunConfiguration,
         cyclelog: CycleLog,
     ) -> anyhow::Result<Self> {
-        let posterior = posterior::posterior(&psi, &weights)?;
-
         Ok(Self {
             equation,
             data,
+            error_models,
+            prior,
             theta,
             psi,
             weights,
             objf,
             cycles,
             status,
-            run_configuration,
             cyclelog,
-            predictions: None,
-            posterior,
         })
     }
 
@@ -72,6 +69,14 @@ impl<E: Equation> NonparametricWorkspace<E> {
         &self.theta
     }
 
+    /// The prior distribution ([`Theta`]) that seeded the algorithm.
+    ///
+    /// This is the initial set of support points, as opposed to the optimized
+    /// solution returned by [`get_theta`](Self::get_theta).
+    pub fn prior(&self) -> &Theta {
+        &self.prior
+    }
+
     pub fn data(&self) -> &Data {
         &self.data
     }
@@ -84,33 +89,6 @@ impl<E: Equation> NonparametricWorkspace<E> {
         &self.cyclelog
     }
 
-    pub(crate) fn run_configuration(&self) -> &RunConfiguration {
-        &self.run_configuration
-    }
-
-    pub(crate) fn algorithm(&self) -> crate::algorithms::Algorithm {
-        self.run_configuration.algorithm
-    }
-
-    pub(crate) fn output_folder(&self) -> &str {
-        self.run_configuration.output_path()
-    }
-
-    pub(crate) fn should_write_outputs(&self) -> bool {
-        self.run_configuration.should_write_outputs()
-    }
-
-    pub(crate) fn prediction_interval(&self) -> (f64, f64) {
-        (
-            self.run_configuration.runtime.idelta,
-            self.run_configuration.runtime.tad,
-        )
-    }
-
-    pub fn predictions(&self) -> Option<&NPPredictions> {
-        self.predictions.as_ref()
-    }
-
     pub fn psi(&self) -> &Psi {
         &self.psi
     }
@@ -119,26 +97,35 @@ impl<E: Equation> NonparametricWorkspace<E> {
         &self.weights
     }
 
-    pub fn posterior(&self) -> &Posterior {
-        &self.posterior
+    pub fn error_models(&self) -> &AssayErrorModels {
+        &self.error_models
+    }
+
+    /// Compute the posterior probabilities on demand from [`Psi`] and the
+    /// [`Weights`]. This is a cheap matrix operation and is intentionally not
+    /// cached on the result.
+    pub fn posterior(&self) -> anyhow::Result<Posterior> {
+        Posterior::calculate(&self.psi, &self.weights)
     }
 
-    pub fn calculate_predictions(&mut self, idelta: f64, tad: f64) -> anyhow::Result<()> {
-        let predictions = NPPredictions::calculate(
+    /// Compute predictions on demand. Nothing is cached on the result; callers
+    /// that need the predictions repeatedly should hold on to the returned
+    /// value themselves.
+    pub fn predictions(&self, idelta: f64, tad: f64) -> anyhow::Result<NPPredictions> {
+        let posterior = self.posterior()?;
+        NPPredictions::calculate(
             &self.equation,
             &self.data,
             &self.theta,
             &self.weights,
-            &self.posterior,
+            &posterior,
             idelta,
             tad,
-        )?;
-        self.predictions = Some(predictions);
-        Ok(())
+        )
     }
 
-    pub fn write_theta(&self) -> anyhow::Result<()> {
-        use anyhow::{bail, Context};
+    pub fn write_theta(&self, path: &Path) -> anyhow::Result<()> {
+        use anyhow::bail;
         use csv::WriterBuilder;
 
         tracing::debug!("Writing population parameter distribution...");
@@ -152,14 +139,11 @@ impl<E: Equation> NonparametricWorkspace<E> {
             );
         }
 
-        let outputfile = crate::output::OutputFile::new(self.output_folder(), "theta.csv")
-            .context("Failed to create output file for theta")?;
-
-        let mut writer = WriterBuilder::new()
-            .has_headers(true)
-            .from_writer(outputfile.file());
+        std::fs::create_dir_all(path)?;
+        let file = std::fs::File::create(path)?;
+        let mut writer = WriterBuilder::new().has_headers(true).from_writer(file);
 
-        let mut theta_header = self.run_configuration.parameter_names.clone();
+        let mut theta_header = self.theta.param_names().clone();
         theta_header.push("prob".to_string());
         writer.write_record(&theta_header)?;
 
@@ -172,16 +156,14 @@ impl<E: Equation> NonparametricWorkspace<E> {
         Ok(())
     }
 
-    pub fn write_posterior(&self) -> anyhow::Result<()> {
+    pub fn write_posterior(&self, path: &str) -> anyhow::Result<()> {
         use csv::WriterBuilder;
 
         tracing::debug!("Writing posterior parameter probabilities...");
 
-        let outputfile = crate::output::OutputFile::new(self.output_folder(), "posterior.csv")?;
-
-        let mut writer = WriterBuilder::new()
-            .has_headers(true)
-            .from_writer(outputfile.file());
+        std::fs::create_dir_all(path)?;
+        let file = std::fs::File::create(path)?;
+        let mut writer = WriterBuilder::new().has_headers(true).from_writer(file);
 
         writer.write_field("id")?;
         writer.write_field("point")?;
@@ -191,8 +173,9 @@ impl<E: Equation> NonparametricWorkspace<E> {
         writer.write_field("prob")?;
         writer.write_record(None::<&[u8]>)?;
 
+        let posterior = self.posterior()?;
         let subjects = self.data.subjects();
-        self.posterior
+        posterior
             .matrix()
             .row_iter()
             .enumerate()
@@ -217,15 +200,14 @@ impl<E: Equation> NonparametricWorkspace<E> {
         Ok(())
     }
 
-    pub fn write_covariates(&self) -> anyhow::Result<()> {
+    pub fn write_covariates(&self, path: &Path) -> anyhow::Result<()> {
         use csv::WriterBuilder;
         use pharmsol::Event;
 
         tracing::debug!("Writing covariates...");
-        let outputfile = crate::output::OutputFile::new(self.output_folder(), "covariates.csv")?;
-        let mut writer = WriterBuilder::new()
-            .has_headers(true)
-            .from_writer(outputfile.file());
+        std::fs::create_dir_all(path)?;
+        let file = std::fs::File::create(path)?;
+        let mut writer = WriterBuilder::new().has_headers(true).from_writer(file);
 
         let mut covariate_names = std::collections::HashSet::new();
         for subject in self.data.subjects() {
@@ -279,8 +261,4 @@ impl<E: Equation> NonparametricWorkspace<E> {
         writer.flush()?;
         Ok(())
     }
-
-    pub fn into_fit_result(self) -> FitResult<E> {
-        FitResult::Nonparametric(self)
-    }
 }
diff --git a/src/estimation/nonparametric/prior/latin.rs b/src/estimation/nonparametric/sampling/latin.rs
similarity index 68%
rename from src/estimation/nonparametric/prior/latin.rs
rename to src/estimation/nonparametric/sampling/latin.rs
index a58c0e94c..b6f7f4f32 100644
--- a/src/estimation/nonparametric/prior/latin.rs
+++ b/src/estimation/nonparametric/sampling/latin.rs
@@ -4,15 +4,10 @@ use rand::prelude::*;
 use rand::rngs::StdRng;
 
 use crate::estimation::nonparametric::Theta;
-use crate::model::ParameterSpace;
-
-pub fn generate(
-    parameters: impl Into<ParameterSpace>,
-    points: usize,
-    seed: usize,
-) -> Result<Theta> {
-    let parameters = parameters.into();
-    let ranges = parameters.finite_ranges()?;
+use crate::model::{BoundedParameter, ParameterSpace};
+
+pub fn generate(parameters: &ParameterSpace<BoundedParameter>, points: usize, seed: usize) -> Result<Theta> {
+    let ranges = parameters.finite_ranges();
     let mut rng = StdRng::seed_from_u64(seed as u64);
 
     let mut intervals = Vec::new();
@@ -36,14 +31,14 @@ pub fn generate(
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::model::{Parameter, ParameterSpace};
+    use crate::model::{BoundedParameter, ParameterSpace};
 
     #[test]
     fn latin_generate_produces_requested_shape() {
-        let params = ParameterSpace::new()
-            .add(Parameter::bounded("a", 0.0, 1.0))
-            .add(Parameter::bounded("b", 0.0, 1.0))
-            .add(Parameter::bounded("c", 0.0, 1.0));
+        let params = ParameterSpace::<BoundedParameter>::new()
+            .add("a", 0.0, 1.0)
+            .add("b", 0.0, 1.0)
+            .add("c", 0.0, 1.0);
 
         let theta = generate(&params, 10, 22).unwrap();
         assert_eq!(theta.nspp(), 10);
diff --git a/src/estimation/nonparametric/sampling/mod.rs b/src/estimation/nonparametric/sampling/mod.rs
new file mode 100644
index 000000000..899edc1cf
--- /dev/null
+++ b/src/estimation/nonparametric/sampling/mod.rs
@@ -0,0 +1,10 @@
+//! Methods for generating the initial set of support points ([`Theta`]).
+
+pub mod latin;
+pub mod sobol;
+
+/// Default seed used for quasi-random sampling when none is given.
+pub const DEFAULT_SEED: usize = 22;
+
+/// Default number of support points sampled for a quasi-random starting grid.
+pub const DEFAULT_POINTS: usize = 2028;
diff --git a/src/estimation/nonparametric/prior/sobol.rs b/src/estimation/nonparametric/sampling/sobol.rs
similarity index 65%
rename from src/estimation/nonparametric/prior/sobol.rs
rename to src/estimation/nonparametric/sampling/sobol.rs
index 61a162fb3..1dc8f2af0 100644
--- a/src/estimation/nonparametric/prior/sobol.rs
+++ b/src/estimation/nonparametric/sampling/sobol.rs
@@ -3,16 +3,15 @@ use faer::Mat;
 use sobol_burley::sample;
 
 use crate::estimation::nonparametric::Theta;
-use crate::model::ParameterSpace;
+use crate::model::{BoundedParameter, ParameterSpace};
 
 pub fn generate(
-    parameters: impl Into<ParameterSpace>,
+    parameters: &ParameterSpace<BoundedParameter>,
     points: usize,
     seed: usize,
 ) -> Result<Theta> {
-    let parameters = parameters.into();
     let seed = seed as u32;
-    let ranges = parameters.finite_ranges()?;
+    let ranges = parameters.finite_ranges();
 
     let rand_matrix = Mat::from_fn(points, ranges.len(), |i, j| {
         let unscaled = sample((i).try_into().unwrap(), j.try_into().unwrap(), seed) as f64;
@@ -26,14 +25,13 @@ pub fn generate(
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::model::{Parameter, ParameterSpace};
 
     #[test]
     fn sobol_generate_produces_requested_shape() {
-        let params = ParameterSpace::new()
-            .add(Parameter::bounded("a", 0.0, 1.0))
-            .add(Parameter::bounded("b", 0.0, 1.0))
-            .add(Parameter::bounded("c", 0.0, 1.0));
+        let params = ParameterSpace::<BoundedParameter>::new()
+            .add("a", 0.0, 1.0)
+            .add("b", 0.0, 1.0)
+            .add("c", 0.0, 1.0);
 
         let theta = generate(&params, 10, 22).unwrap();
         assert_eq!(theta.nspp(), 10);
diff --git a/src/estimation/nonparametric/summaries.rs b/src/estimation/nonparametric/summaries.rs
index e0d145af9..6d979aaef 100644
--- a/src/estimation/nonparametric/summaries.rs
+++ b/src/estimation/nonparametric/summaries.rs
@@ -1,11 +1,11 @@
 use ndarray::{Array1, Array2};
 use pharmsol::{Data, Equation, Event};
 
-use crate::estimation::nonparametric::NonparametricWorkspace;
+use crate::estimation::nonparametric::NonParametricResult;
 use crate::estimation::nonparametric::{population_mean_median, posterior_mean_median};
 use crate::results::{FitSummary, IndividualSummary, ParameterSummary, PopulationSummary};
 
-pub fn fit_summary<E: Equation>(result: &NonparametricWorkspace<E>) -> FitSummary {
+pub fn fit_summary<E: Equation>(result: &NonParametricResult<E>) -> FitSummary {
     FitSummary {
         objective_function: result.objf(),
         converged: result.converged(),
@@ -13,11 +13,10 @@ pub fn fit_summary<E: Equation>(result: &NonparametricWorkspace<E>) -> FitSummar
         subject_count: result.data().subjects().len(),
         observation_count: count_observations(result.data()),
         parameter_count: result.get_theta().parameters().len(),
-        algorithm: format!("{:?}", result.algorithm()),
     }
 }
 
-pub fn population_summary<E: Equation>(result: &NonparametricWorkspace<E>) -> PopulationSummary {
+pub fn population_summary<E: Equation>(result: &NonParametricResult<E>) -> PopulationSummary {
     let theta_matrix = to_ndarray_matrix(result.get_theta().matrix());
     let weights = Array1::from_iter(result.weights().iter());
     let (mean, median) = population_mean_median(&theta_matrix, &weights)
@@ -53,7 +52,7 @@ pub fn population_summary<E: Equation>(result: &NonparametricWorkspace<E>) -> Po
 }
 
 pub fn individual_summaries<E: Equation>(
-    result: &NonparametricWorkspace<E>,
+    result: &NonParametricResult<E>,
 ) -> Vec<IndividualSummary> {
     let theta_matrix = to_ndarray_matrix(result.get_theta().matrix());
     let psi_matrix = to_ndarray_matrix(result.psi().matrix());
diff --git a/src/estimation/nonparametric/theta.rs b/src/estimation/nonparametric/theta.rs
index 2d5052bd3..abe0c0a70 100644
--- a/src/estimation/nonparametric/theta.rs
+++ b/src/estimation/nonparametric/theta.rs
@@ -1,11 +1,12 @@
-use std::fmt::Debug;
+use std::{fmt::Debug, fs::File, path::Path};
 
-use anyhow::{bail, Result};
+use anyhow::{bail, Context, Result};
 use faer::Mat;
 use serde::{Deserialize, Serialize};
 
+use super::sampling::{self, latin, sobol};
 use super::weights::Weights;
-use crate::model::ParameterSpace;
+use crate::model::{BoundedParameter, ParameterSpace};
 
 /// [Theta] is a structure that holds the support points
 /// These represent the joint population parameter distribution
@@ -14,14 +15,14 @@ use crate::model::ParameterSpace;
 #[derive(Clone, PartialEq)]
 pub struct Theta {
     matrix: Mat<f64>,
-    parameters: ParameterSpace,
+    parameters: ParameterSpace<BoundedParameter>,
 }
 
 impl Default for Theta {
     fn default() -> Self {
         Theta {
             matrix: Mat::new(),
-            parameters: ParameterSpace::new(),
+            parameters: ParameterSpace::<BoundedParameter>::new(),
         }
     }
 }
@@ -37,8 +38,10 @@ impl Theta {
     /// in the [ParameterSpace]
     ///
     /// The order of parameters in the [ParameterSpace] should match the order of columns in the matrix
-    pub fn from_parts(matrix: Mat<f64>, parameters: impl Into<ParameterSpace>) -> Result<Self> {
-        let parameters = parameters.into();
+    pub fn from_parts(
+        matrix: Mat<f64>,
+        parameters: ParameterSpace<BoundedParameter>,
+    ) -> Result<Self> {
         if matrix.ncols() != parameters.len() {
             bail!(
                 "Number of columns in matrix ({}) does not match number of parameters ({})",
@@ -63,12 +66,12 @@ impl Theta {
     }
 
     /// Get the [ParameterSpace] associated with this [Theta]
-    pub fn parameters(&self) -> &ParameterSpace {
+    pub fn parameters(&self) -> &ParameterSpace<BoundedParameter> {
         &self.parameters
     }
 
     /// Get a mutable reference to the [ParameterSpace]
-    pub fn parameters_mut(&mut self) -> &mut ParameterSpace {
+    pub fn parameters_mut(&mut self) -> &mut ParameterSpace<BoundedParameter> {
         &mut self.parameters
     }
 
@@ -124,10 +127,7 @@ impl Theta {
             return true;
         }
 
-        let limits = self
-            .parameters
-            .finite_ranges()
-            .expect("theta requires finite parameter bounds");
+        let limits = self.parameters.finite_ranges();
 
         for row_idx in 0..self.matrix.nrows() {
             let mut squared_dist = 0.0;
@@ -226,10 +226,143 @@ impl Theta {
         }
 
         let mat = Mat::from_fn(nrows, ncols, |i, j| rows[i][j]);
-        let parameters = ParameterSpace::new();
+        let parameters = ParameterSpace::<BoundedParameter>::new();
 
         Theta::from_parts(mat, parameters)
     }
+
+    /// Generate a starting grid of `points` support points over `parameters`
+    /// using a Sobol sequence and the default seed ([`sampling::DEFAULT_SEED`]).
+    ///
+    /// The returned [Theta] carries `parameters`, so the chosen grid is explicit
+    /// and self-describing.
+    pub fn sobol(parameters: &ParameterSpace<BoundedParameter>, points: usize) -> Result<Self> {
+        Self::sobol_with_seed(parameters, points, sampling::DEFAULT_SEED)
+    }
+
+    /// Generate a starting grid over `parameters` using a Sobol sequence with the
+    /// default number of support points ([`sampling::DEFAULT_POINTS`]) and the
+    /// default seed ([`sampling::DEFAULT_SEED`]).
+    pub fn sobol_default(parameters: &ParameterSpace<BoundedParameter>) -> Result<Self> {
+        Self::sobol(parameters, sampling::DEFAULT_POINTS)
+    }
+
+    /// Like [`Theta::sobol`], with an explicit seed for the quasi-random sequence.
+    pub fn sobol_with_seed(
+        parameters: &ParameterSpace<BoundedParameter>,
+        points: usize,
+        seed: usize,
+    ) -> Result<Self> {
+        validate_bounds(parameters)?;
+        sobol::generate(parameters, points, seed)
+    }
+
+    /// Generate a starting grid of `points` support points over `parameters`
+    /// using Latin Hypercube Sampling and the default seed ([`sampling::DEFAULT_SEED`]).
+    ///
+    /// The returned [Theta] carries `parameters`, so the chosen grid is explicit
+    /// and self-describing.
+    pub fn latin(parameters: &ParameterSpace<BoundedParameter>, points: usize) -> Result<Self> {
+        Self::latin_with_seed(parameters, points, sampling::DEFAULT_SEED)
+    }
+
+    /// Like [`Theta::latin`], with an explicit seed for the quasi-random sequence.
+    pub fn latin_with_seed(
+        parameters: &ParameterSpace<BoundedParameter>,
+        points: usize,
+        seed: usize,
+    ) -> Result<Self> {
+        validate_bounds(parameters)?;
+        latin::generate(parameters, points, seed)
+    }
+
+    pub fn from_file(
+        path: impl AsRef<Path>,
+        parameters: &ParameterSpace<BoundedParameter>,
+    ) -> Result<(Theta, Option<Weights>)> {
+        let path = path.as_ref();
+        tracing::info!("Reading prior from {}", path.display());
+        let file = File::open(path).context(format!(
+            "Unable to open the prior file '{}'",
+            path.display()
+        ))?;
+        let mut reader = csv::ReaderBuilder::new()
+            .has_headers(true)
+            .from_reader(file);
+
+        let mut parameter_names: Vec<String> = reader
+            .headers()?
+            .clone()
+            .into_iter()
+            .map(|s| s.trim().to_owned())
+            .collect();
+
+        let prob_index = parameter_names.iter().position(|name| name == "prob");
+        if let Some(index) = prob_index {
+            parameter_names.remove(index);
+        }
+
+        let random_names: Vec<String> = parameters.names();
+
+        let mut reordered_indices: Vec<usize> = Vec::new();
+        for random_name in &random_names {
+            match parameter_names.iter().position(|name| name == random_name) {
+                Some(index) => {
+                    let adjusted_index = if let Some(prob_idx) = prob_index {
+                        if index >= prob_idx {
+                            index + 1
+                        } else {
+                            index
+                        }
+                    } else {
+                        index
+                    };
+                    reordered_indices.push(adjusted_index);
+                }
+                None => bail!("Parameter {} is not present in the CSV file.", random_name),
+            }
+        }
+
+        if parameter_names.len() > random_names.len() {
+            let extra_parameters: Vec<&String> = parameter_names.iter().collect();
+            bail!(
+                "Found parameters in the prior not present in configuration: {:?}",
+                extra_parameters
+            );
+        }
+
+        let mut theta_values = Vec::new();
+        let mut prob_values = Vec::new();
+
+        for result in reader.records() {
+            let record = result.unwrap();
+            let values: Vec<f64> = reordered_indices
+                .iter()
+                .map(|&i| record[i].parse::<f64>().unwrap())
+                .collect();
+            theta_values.push(values);
+
+            if let Some(prob_idx) = prob_index {
+                let prob_value: f64 = record[prob_idx].parse::<f64>().unwrap();
+                prob_values.push(prob_value);
+            }
+        }
+
+        let n_points = theta_values.len();
+        let n_params = random_names.len();
+        let theta_values: Vec<f64> = theta_values.into_iter().flatten().collect();
+        let theta_matrix: Mat<f64> =
+            Mat::from_fn(n_points, n_params, |i, j| theta_values[i * n_params + j]);
+
+        let theta = Theta::from_parts(theta_matrix, parameters.clone())?;
+        let weights = if !prob_values.is_empty() {
+            Some(Weights::from_vec(prob_values))
+        } else {
+            None
+        };
+
+        Ok((theta, weights))
+    }
 }
 
 impl Debug for Theta {
@@ -281,7 +414,7 @@ impl<'de> Deserialize<'de> for Theta {
         #[derive(Deserialize)]
         struct ThetaSerde {
             matrix: Vec<Vec<f64>>,
-            parameters: ParameterSpace,
+            parameters: ParameterSpace<BoundedParameter>,
         }
 
         let decoded = ThetaSerde::deserialize(deserializer)?;
@@ -310,3 +443,86 @@ impl<'de> Deserialize<'de> for Theta {
         Self::from_parts(matrix, decoded.parameters).map_err(serde::de::Error::custom)
     }
 }
+
+/// Validates that every parameter has a strictly-ordered, finite bound interval.
+fn validate_bounds(parameters: &ParameterSpace<BoundedParameter>) -> Result<()> {
+    for parameter in parameters.iter() {
+        if parameter.lower >= parameter.upper {
+            bail!(
+                "Parameter '{}' has invalid bounds: [{}, {}]. Lower bound must be less than upper bound.",
+                parameter.name,
+                parameter.lower,
+                parameter.upper
+            );
+        }
+    }
+    Ok(())
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use std::fs;
+
+    fn parameters() -> ParameterSpace<BoundedParameter> {
+        ParameterSpace::<BoundedParameter>::new()
+            .add("ke", 0.1, 1.0)
+            .add("v", 5.0, 50.0)
+    }
+
+    fn temp_csv_path() -> String {
+        format!("test_temp_theta_{}.csv", rand::random::<u32>())
+    }
+
+    #[test]
+    fn sobol_generates_expected_shape() {
+        let theta = Theta::sobol_with_seed(&parameters(), 10, 42).unwrap();
+        assert_eq!(theta.nspp(), 10);
+        assert_eq!(theta.matrix().ncols(), 2);
+    }
+
+    #[test]
+    fn latin_generates_expected_shape() {
+        let theta = Theta::latin(&parameters(), 10).unwrap();
+        assert_eq!(theta.nspp(), 10);
+        assert_eq!(theta.matrix().ncols(), 2);
+    }
+
+    #[test]
+    fn sampling_rejects_invalid_bounds() {
+        let bad = ParameterSpace::<BoundedParameter>::new().add("ke", 1.0, 1.0);
+        let err = Theta::sobol(&bad, 10).unwrap_err();
+        assert!(err.to_string().contains("invalid bounds"));
+    }
+
+    #[test]
+    fn from_file_parses_weights_and_reorders_columns() {
+        let path = temp_csv_path();
+        fs::write(&path, "v,ke,prob\n10.0,0.5,0.3\n15.0,0.7,0.7\n").unwrap();
+
+        let (theta, weights) = Theta::from_file(&path, &parameters()).unwrap();
+        let _ = fs::remove_file(&path);
+
+        assert_eq!(theta.nspp(), 2);
+        assert_eq!(theta.matrix()[(0, 0)], 0.5);
+        assert_eq!(theta.matrix()[(0, 1)], 10.0);
+
+        let weights = weights.expect("weights should be parsed from prob column");
+        assert_eq!(weights.len(), 2);
+        assert_eq!(weights[0], 0.3);
+        assert_eq!(weights[1], 0.7);
+    }
+
+    #[test]
+    fn from_file_rejects_extra_parameters() {
+        let path = temp_csv_path();
+        fs::write(&path, "ke,v,extra\n0.5,10.0,1.0\n").unwrap();
+
+        let err = Theta::from_file(&path, &parameters()).unwrap_err();
+        let _ = fs::remove_file(&path);
+
+        assert!(err
+            .to_string()
+            .contains("Found parameters in the prior not present in configuration"));
+    }
+}
diff --git a/src/estimation/problem.rs b/src/estimation/problem.rs
new file mode 100644
index 000000000..ef66805c2
--- /dev/null
+++ b/src/estimation/problem.rs
@@ -0,0 +1,378 @@
+use anyhow::{anyhow, Result};
+use pharmsol::{
+    AssayErrorModel, AssayErrorModels, Data, Equation, Event, ResidualErrorModel,
+    ResidualErrorModels,
+};
+use std::collections::{BTreeSet, HashSet};
+
+use crate::estimation::nonparametric::Theta;
+use crate::model::parameter_space::{BoundedParameter, ParameterSpace, UnboundedParameter};
+use crate::model::{EquationMetadataSource, Model, ModelBuilder};
+
+pub trait Framework {
+    type ErrorModels;
+    /// The prior that seeds the algorithm.
+    ///
+    /// For the non-parametric framework this is a [`Theta`] (a discrete prior
+    /// distribution that also carries the parameter space). For the parametric
+    /// framework it is the [`ParameterSpace`] of unbounded parameters.
+    type Prior;
+}
+
+pub struct Parametric;
+
+impl Framework for Parametric {
+    type ErrorModels = ResidualErrorModels;
+    type Prior = ParameterSpace<UnboundedParameter>;
+}
+
+pub struct NonParametric;
+
+impl Framework for NonParametric {
+    type ErrorModels = AssayErrorModels;
+    type Prior = Theta;
+}
+
+#[derive(Debug, Clone)]
+pub struct EstimationProblem<E: Equation, F: Framework> {
+    pub(crate) model: Model<E>,
+    pub(crate) data: Data,
+    pub(crate) error_models: F::ErrorModels,
+    /// The prior that seeds the algorithm.
+    ///
+    /// For the non-parametric framework this is the prior [`Theta`], which also
+    /// carries the parameter space. The parameter space is therefore not stored
+    /// separately.
+    pub(crate) prior: F::Prior,
+}
+
+impl<E: Equation + EquationMetadataSource> EstimationProblem<E, NonParametric> {
+    /// Creates a non-parametric estimation problem.
+    ///
+    /// The `prior` is a [`Theta`] holding the prior distribution (the initial
+    /// set of support points) together with the [`ParameterSpace`] it was built
+    /// from. The parameter space is taken directly from the prior, so there is
+    /// no separate parameter-declaration step.
+    pub fn nonparametric(
+        equation: E,
+        data: Data,
+        prior: Theta,
+        error_models: AssayErrorModels,
+    ) -> Result<Self> {
+        let model_builder = Model::builder(equation);
+
+        validate_nonparametric_parameters(&model_builder, prior.parameters())?;
+
+        let model = model_builder.build()?;
+
+        validate_nonparametric_error_models(&model, &data, &error_models)?;
+
+        Ok(EstimationProblem {
+            model,
+            data,
+            error_models,
+            prior,
+        })
+    }
+}
+
+impl<E: Equation> EstimationProblem<E, Parametric> {
+    /// Begins building a parametric estimation problem.
+    pub fn parametric(equation: E, data: Data) -> ParametricBuilder<E> {
+        ParametricBuilder {
+            model: Model::builder(equation),
+            data,
+            parameters: ParameterSpace::<UnboundedParameter>::new(),
+            error_models: Vec::new(),
+        }
+    }
+
+    /// Returns the parameter space defined for this problem.
+    pub fn parameters(&self) -> &ParameterSpace<UnboundedParameter> {
+        &self.prior
+    }
+}
+
+impl<E: Equation> EstimationProblem<E, NonParametric> {
+    /// Returns the parameter space carried by the prior [`Theta`].
+    pub fn parameters(&self) -> &ParameterSpace<BoundedParameter> {
+        self.prior.parameters()
+    }
+}
+
+pub struct ParametricBuilder<E: Equation> {
+    model: ModelBuilder<E>,
+    data: Data,
+    parameters: ParameterSpace<UnboundedParameter>,
+    error_models: Vec<(String, ResidualErrorModel)>,
+}
+
+impl<E: Equation> ParametricBuilder<E> {
+    pub fn parameter(mut self, parameter: impl Into<UnboundedParameter>) -> Self {
+        self.parameters.push(parameter.into());
+        self
+    }
+
+    pub fn parameters<P, I>(mut self, parameters: I) -> Self
+    where
+        P: Into<UnboundedParameter>,
+        I: IntoIterator<Item = P>,
+    {
+        for param in parameters {
+            self.parameters.push(param.into());
+        }
+        self
+    }
+
+    pub fn error_model(mut self, name: impl Into<String>, model: ResidualErrorModel) -> Self {
+        self.error_models.push((name.into(), model));
+        self
+    }
+}
+
+impl<E: Equation + EquationMetadataSource> ParametricBuilder<E> {
+    pub fn build(self) -> Result<EstimationProblem<E, Parametric>> {
+        validate_parametric_parameters(&self.model, &self.parameters)?;
+        validate_parametric_error_models(&self.model, &self.error_models)?;
+
+        let mut all_errors = ResidualErrorModels::new();
+        for (name, error_model) in self.error_models {
+            let outeq = self
+                .model
+                .output_index(&name)
+                .ok_or_else(|| anyhow!("unknown equation output label: {name}"))?;
+
+            all_errors = all_errors.add(outeq, error_model);
+        }
+
+        Ok(EstimationProblem {
+            model: self.model.build()?,
+            data: self.data,
+            error_models: all_errors,
+            prior: self.parameters,
+        })
+    }
+}
+
+fn validate_nonparametric_parameters<E: Equation + EquationMetadataSource>(
+    model: &ModelBuilder<E>,
+    parameters: &ParameterSpace<BoundedParameter>,
+) -> Result<()> {
+    if parameters.is_empty() {
+        anyhow::bail!("at least one parameter is required for non-parametric models");
+    }
+
+    for parameter in parameters.iter() {
+        if !parameter.lower.is_finite() || !parameter.upper.is_finite() {
+            anyhow::bail!(
+                "invalid bounds for parameter '{}': bounds must be finite numbers",
+                parameter.name
+            );
+        }
+
+        if parameter.lower >= parameter.upper {
+            anyhow::bail!(
+                "invalid bounds for parameter '{}': lower bound ({}) must be strictly less than upper bound ({})",
+                parameter.name,
+                parameter.lower,
+                parameter.upper
+            );
+        }
+    }
+
+    let names: Vec<String> = parameters
+        .iter()
+        .map(|parameter| parameter.name.clone())
+        .collect();
+    validate_parameter_declarations(model, &names)
+}
+
+fn validate_parametric_parameters<E: Equation + EquationMetadataSource>(
+    model: &ModelBuilder<E>,
+    parameters: &ParameterSpace<UnboundedParameter>,
+) -> Result<()> {
+    if parameters.is_empty() {
+        anyhow::bail!("at least one parameter is required for parametric models");
+    }
+
+    let names: Vec<String> = parameters
+        .iter()
+        .map(|parameter| parameter.name.clone())
+        .collect();
+    validate_parameter_declarations(model, &names)
+}
+
+fn validate_parameter_declarations<E: Equation + EquationMetadataSource>(
+    model: &ModelBuilder<E>,
+    provided_names: &[String],
+) -> Result<()> {
+    let mut seen: HashSet<&str> = HashSet::new();
+    let mut duplicates: Vec<String> = Vec::new();
+    for name in provided_names {
+        if !seen.insert(name.as_str()) {
+            duplicates.push(name.clone());
+        }
+    }
+
+    if !duplicates.is_empty() {
+        duplicates.sort();
+        duplicates.dedup();
+        anyhow::bail!(
+            "duplicate parameter declarations found: {}",
+            duplicates.join(", ")
+        );
+    }
+
+    let declared = model.parameter_names();
+
+    let unknown: Vec<String> = provided_names
+        .iter()
+        .filter(|name| model.parameter_index(name).is_none())
+        .cloned()
+        .collect();
+    if !unknown.is_empty() {
+        anyhow::bail!(
+            "unknown parameter name(s): {}. Valid parameters are: {}",
+            unknown.join(", "),
+            declared.join(", ")
+        );
+    }
+
+    let provided: HashSet<&str> = provided_names.iter().map(|name| name.as_str()).collect();
+    let missing: Vec<String> = declared
+        .iter()
+        .filter(|name| !provided.contains(name.as_str()))
+        .cloned()
+        .collect();
+
+    if !missing.is_empty() {
+        anyhow::bail!("missing parameter declaration(s): {}", missing.join(", "));
+    }
+
+    Ok(())
+}
+
+fn validate_nonparametric_error_models<E: Equation + EquationMetadataSource>(
+    model: &Model<E>,
+    data: &Data,
+    error_models: &AssayErrorModels,
+) -> Result<()> {
+    // Bind the (label-first) error models to the equation. This resolves and
+    // validates that every declared output label maps to a valid model output.
+    let bound = model
+        .equation
+        .bind_error_models(error_models)
+        .map_err(|e| anyhow!("invalid assay error model output(s): {e}"))?;
+
+    // Collect the set of model output indices that are actually observed in the
+    // data, resolving each observation's output label the same way the simulator
+    // does (exact name, then the `outeq_<N>` numeric alias).
+    let mut observed_outputs: BTreeSet<usize> = BTreeSet::new();
+    let mut unresolved_labels: BTreeSet<String> = BTreeSet::new();
+    for subject in data.subjects() {
+        for occasion in subject.occasions() {
+            for event in occasion.events() {
+                if let Event::Observation(obs) = event {
+                    let label = obs.outeq().to_string();
+                    match resolve_output_index(model, &label) {
+                        Some(outeq) => {
+                            observed_outputs.insert(outeq);
+                        }
+                        None => {
+                            unresolved_labels.insert(label);
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    if !unresolved_labels.is_empty() {
+        let labels: Vec<String> = unresolved_labels.into_iter().collect();
+        anyhow::bail!(
+            "the data references output label(s) that are not defined by the model: {}",
+            labels.join(", ")
+        );
+    }
+
+    if observed_outputs.is_empty() {
+        anyhow::bail!("the data contains no observations to fit");
+    }
+
+    // Every observed output must have a (non-`None`) assay error model.
+    for &outeq in &observed_outputs {
+        let has_model = matches!(
+            bound.error_model(outeq),
+            Ok(error_model) if *error_model != AssayErrorModel::None
+        );
+
+        if !has_model {
+            let label = model
+                .output_name(outeq)
+                .map(|name| name.to_string())
+                .unwrap_or_else(|| outeq.to_string());
+            anyhow::bail!(
+                "no assay error model defined for output '{}' (index {}), which is observed in the data",
+                label,
+                outeq
+            );
+        }
+    }
+
+    Ok(())
+}
+
+/// Resolves an observation output `label` to a model output index, mirroring the
+/// simulator: first by exact output name, then via the `outeq_<N>` numeric alias.
+fn resolve_output_index<E: Equation + EquationMetadataSource>(
+    model: &Model<E>,
+    label: &str,
+) -> Option<usize> {
+    model.output_index(label).or_else(|| {
+        if !label.is_empty() && label.bytes().all(|b| b.is_ascii_digit()) {
+            model.output_index(&format!("outeq_{label}"))
+        } else {
+            None
+        }
+    })
+}
+
+fn validate_parametric_error_models<E: Equation + EquationMetadataSource>(
+    model: &ModelBuilder<E>,
+    error_models: &[(String, ResidualErrorModel)],
+) -> Result<()> {
+    if error_models.is_empty() {
+        anyhow::bail!("at least one residual error model is required");
+    }
+
+    validate_error_model_labels(model, error_models.iter().map(|(name, _)| name.as_str()))
+}
+
+fn validate_error_model_labels<'a, E, I>(model: &ModelBuilder<E>, labels: I) -> Result<()>
+where
+    E: Equation + EquationMetadataSource,
+    I: IntoIterator<Item = &'a str>,
+{
+    let valid_outputs = model.output_names();
+    let mut seen_output_indexes: HashSet<usize> = HashSet::new();
+
+    for name in labels {
+        let outeq = model.output_index(name).ok_or_else(|| {
+            anyhow!(
+                "unknown equation output label: {}. Valid outputs are: {}",
+                name,
+                valid_outputs.join(", ")
+            )
+        })?;
+
+        if !seen_output_indexes.insert(outeq) {
+            anyhow::bail!(
+                "duplicate error model declaration for output '{}' (index {})",
+                name,
+                outeq
+            );
+        }
+    }
+
+    Ok(())
+}
diff --git a/src/api/progress.rs b/src/estimation/progress.rs
similarity index 100%
rename from src/api/progress.rs
rename to src/estimation/progress.rs
diff --git a/src/lib.rs b/src/lib.rs
index 7018af0c0..1ef54e274 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -11,11 +11,11 @@
 //! - NPOD (Non-Parametric Optimal Design)
 //! - POSTPROB (Posterior probability reweighting)
 //!
-//! # Public API
+//! # Public Interface
 //!
-//! PMcore centers on the model/problem API in [api]. Models are defined with
-//! [api::ModelDefinition], configured with [api::EstimationProblem], and executed through
-//! [api::fit].
+//! PMcore centers on the estimation interface in [estimation]. Models are defined in
+//! [model], configured with [estimation::EstimationProblem], and then executed with the
+//! selected algorithm.
 //!
 //! # Data format
 //!
@@ -28,12 +28,6 @@
 /// Provides the various algorithms used within the framework
 pub mod algorithms;
 
-/// New public modeling and execution API.
-pub mod api;
-
-/// Shared preprocessing and compilation layer.
-pub mod compile;
-
 /// Estimation family boundaries for the new architecture.
 pub mod estimation;
 
@@ -43,8 +37,8 @@ pub mod model;
 /// Shared result and summary types for the new API.
 pub mod results;
 
-/// Shared output writers for the new API.
-pub mod output;
+/// Logs
+pub mod logs;
 
 // Re-export commonly used items
 pub use anyhow::Result;
@@ -54,40 +48,43 @@ pub use std::collections::HashMap;
 pub mod bestdose;
 
 /// A collection of commonly used items to simplify imports.
-#[allow(ambiguous_glob_reexports)]
 pub mod prelude {
+    pub use super::logs::Logger;
     pub use super::HashMap;
     pub use super::Result;
     pub use crate::algorithms;
     pub use crate::algorithms::Algorithm;
-    pub use crate::api::fit;
-    pub use crate::api::fit_with_progress;
-    pub use crate::api::{
-        AlgorithmTuning, ConvergenceOptions, ErrorModels, EstimationProblem, FitProgress,
-        MethodSpec, ModelDefinition, NonparametricCycleProgress, Npag, Npod, OutputPlan, PostProb,
-        RuntimeOptions,
+
+    pub use crate::estimation::NonParametric;
+    pub use crate::estimation::Parametric;
+    pub use crate::estimation::{
+        ErrorModels, EstimationProblem, FitProgress, NonParametricAlgorithm,
+        NonparametricCycleProgress, NpagConfig, NpmapConfig, NpodConfig, ParametricAlgorithm,
+        SaemConfig,
     };
-    pub use crate::compile::{CompiledProblem, DesignContext, ObservationIndex};
-    pub use crate::estimation::nonparametric::{
-        CycleLog, NPCycle, NPPredictions, NonparametricEngine, NonparametricWorkspace, Posterior,
-        Psi, Theta, Weights,
+
+    pub use crate::model::parameter_space::{
+        BoundedParameter, Parameter, ParameterScale, ParameterSpace, UnboundedParameter,
     };
-    pub use crate::model::{
-        CovariateEffectsSpec, CovariateModel, CovariateSpec, EquationMetadataSource, ModelMetadata,
-        Parameter, ParameterDomain, ParameterSpace, ParameterTransform as ModelParameterTransform,
-        ParameterVariability, RandomEffectsSpec, VariabilityModel,
+
+    pub use crate::estimation::nonparametric::{
+        CycleLog, NPCycle, NPPredictions, NonParametricResult, Posterior, Psi, Theta, Weights,
     };
+    pub use crate::model::{EquationMetadataSource, ModelMetadata};
     pub use crate::results::{
-        ArtifactIndex, DiagnosticsBundle, FitResult, FitSummary, IndividualSummary,
-        ParameterSummary, PopulationSummary, PredictionsBundle,
+        FitResult, FitSummary, IndividualSummary, ParameterSummary, PopulationSummary,
     };
-    pub use pharmsol::optimize::effect::get_e2;
 
+    // pharmsol: re-export the crate itself and its curated prelude.
     pub use pharmsol;
+    pub use pharmsol::prelude::*;
 
-    pub use crate::api::SaemConfig;
-    pub use crate::estimation::nonparametric::{read_prior, Prior};
+    // Items required by downstream code that are not part of `pharmsol::prelude`.
+    pub use pharmsol::equation::{EquationTypes, Predictions};
+    pub use pharmsol::optimize::effect::get_e2;
+    pub use pharmsol::{ODE, SDE};
 
+    // Organized submodules mirroring pharmsol's grouping.
     pub mod simulator {
         pub use pharmsol::prelude::simulator::*;
     }
@@ -98,20 +95,6 @@ pub mod prelude {
         pub use pharmsol::prelude::models::*;
     }
 
-    //traits
-    pub use pharmsol::data::*;
-    pub use pharmsol::equation::Equation;
-    pub use pharmsol::equation::EquationTypes;
-    pub use pharmsol::equation::Predictions;
-    pub use pharmsol::equation::*;
-    pub use pharmsol::prelude::*;
-    pub use pharmsol::simulator::*;
-    pub use pharmsol::ODE;
-    pub use pharmsol::SDE;
-
-    //macros
-    pub use pharmsol::fa;
-    pub use pharmsol::fetch_cov;
-    pub use pharmsol::fetch_params;
-    pub use pharmsol::lag;
+    // macros
+    pub use pharmsol::{fa, fetch_cov, fetch_params, lag};
 }
diff --git a/src/output/logging.rs b/src/logs.rs
similarity index 99%
rename from src/output/logging.rs
rename to src/logs.rs
index 9f9698ea0..69abca8f4 100644
--- a/src/output/logging.rs
+++ b/src/logs.rs
@@ -172,7 +172,7 @@ where
 ///
 /// ```no_run
 /// use std::time::Instant;
-/// use pmcore::output::Logger::format_elapsed;
+/// use pmcore::logs::format_elapsed;
 ///
 /// let start = Instant::now();
 /// // ... do work ...
diff --git a/src/model/covariate_model.rs b/src/model/covariate_model.rs
deleted file mode 100644
index 0de14c163..000000000
--- a/src/model/covariate_model.rs
+++ /dev/null
@@ -1,444 +0,0 @@
-//! Structured covariate model for population parameter regression.
-
-use anyhow::{bail, Result};
-use faer::{Col, Mat};
-use serde::Serialize;
-use std::collections::HashMap;
-
-/// Covariate model specification.
-#[derive(Debug, Clone)]
-pub struct CovariateModel {
-    param_names: Vec<String>,
-    covariate_names: Vec<String>,
-    covariate_mask: Vec<Vec<bool>>,
-    beta: Col<f64>,
-    estimate_beta: Vec<bool>,
-    reference_values: HashMap<String, f64>,
-}
-
-impl CovariateModel {
-    pub fn new(
-        param_names: Vec<impl Into<String>>,
-        covariate_names: Vec<impl Into<String>>,
-        covariate_mask: Vec<Vec<bool>>,
-    ) -> Result<Self> {
-        let param_names: Vec<String> = param_names.into_iter().map(|s| s.into()).collect();
-        let covariate_names: Vec<String> = covariate_names.into_iter().map(|s| s.into()).collect();
-
-        let n_params = param_names.len();
-        let n_covs = covariate_names.len();
-
-        if covariate_mask.len() != n_params {
-            bail!(
-                "Covariate mask rows ({}) must match number of parameters ({})",
-                covariate_mask.len(),
-                n_params
-            );
-        }
-
-        for (i, row) in covariate_mask.iter().enumerate() {
-            if row.len() != n_covs {
-                bail!(
-                    "Covariate mask row {} has {} columns, expected {}",
-                    i,
-                    row.len(),
-                    n_covs
-                );
-            }
-        }
-
-        let n_beta = Self::count_beta_coefficients(&covariate_mask, n_params);
-
-        Ok(Self {
-            param_names,
-            covariate_names,
-            covariate_mask,
-            beta: Col::zeros(n_beta),
-            estimate_beta: vec![true; n_beta],
-            reference_values: HashMap::new(),
-        })
-    }
-
-    pub fn intercept_only(param_names: Vec<impl Into<String>>) -> Result<Self> {
-        let param_names: Vec<String> = param_names.into_iter().map(|s| s.into()).collect();
-        let n_params = param_names.len();
-
-        Ok(Self {
-            param_names,
-            covariate_names: Vec::new(),
-            covariate_mask: vec![Vec::new(); n_params],
-            beta: Col::zeros(n_params),
-            estimate_beta: vec![true; n_params],
-            reference_values: HashMap::new(),
-        })
-    }
-
-    pub fn from_saemix_matrix(
-        param_names: Vec<impl Into<String>>,
-        covariate_names: Vec<impl Into<String>>,
-        matrix: &[f64],
-    ) -> Result<Self> {
-        let param_names: Vec<String> = param_names.into_iter().map(|s| s.into()).collect();
-        let covariate_names: Vec<String> = covariate_names.into_iter().map(|s| s.into()).collect();
-
-        let n_params = param_names.len();
-        let n_covs = covariate_names.len();
-
-        if matrix.len() != n_params * n_covs {
-            bail!(
-                "Matrix length ({}) doesn't match n_params × n_covs ({} × {} = {})",
-                matrix.len(),
-                n_params,
-                n_covs,
-                n_params * n_covs
-            );
-        }
-
-        let covariate_mask: Vec<Vec<bool>> = (0..n_params)
-            .map(|i| (0..n_covs).map(|j| matrix[i * n_covs + j] != 0.0).collect())
-            .collect();
-
-        Self::new(param_names, covariate_names, covariate_mask)
-    }
-
-    pub fn set_beta(&mut self, beta: Col<f64>) -> Result<()> {
-        let expected = self.n_beta();
-        if beta.nrows() != expected {
-            bail!(
-                "Beta length ({}) doesn't match expected ({})",
-                beta.nrows(),
-                expected
-            );
-        }
-        self.beta = beta;
-        Ok(())
-    }
-
-    pub fn set_intercepts(&mut self, intercepts: &[f64]) -> Result<()> {
-        if intercepts.len() != self.n_params() {
-            bail!(
-                "Intercepts length ({}) doesn't match n_params ({})",
-                intercepts.len(),
-                self.n_params()
-            );
-        }
-
-        let mut idx = 0;
-        for (i, &intercept) in intercepts.iter().enumerate() {
-            self.beta[idx] = intercept;
-            idx += 1;
-            idx += self.covariate_mask[i].iter().filter(|&&x| x).count();
-        }
-
-        Ok(())
-    }
-
-    pub fn set_estimate_beta(&mut self, estimate: Vec<bool>) -> Result<()> {
-        if estimate.len() != self.beta.nrows() {
-            bail!(
-                "estimate_beta length ({}) doesn't match n_beta ({})",
-                estimate.len(),
-                self.beta.nrows()
-            );
-        }
-        self.estimate_beta = estimate;
-        Ok(())
-    }
-
-    pub fn fix_intercept(&mut self, param_idx: usize) -> Result<()> {
-        let beta_idx = self.intercept_beta_index(param_idx)?;
-        self.estimate_beta[beta_idx] = false;
-        Ok(())
-    }
-
-    pub fn set_reference(&mut self, covariate: &str, value: f64) -> Result<()> {
-        if !self.covariate_names.contains(&covariate.to_string()) {
-            bail!("Unknown covariate: {}", covariate);
-        }
-        self.reference_values.insert(covariate.to_string(), value);
-        Ok(())
-    }
-
-    pub fn compute_mu(&self, covariates: &HashMap<String, f64>) -> Col<f64> {
-        let n_params = self.param_names.len();
-        let mut mu = Col::zeros(n_params);
-        let mut beta_idx = 0;
-
-        for i in 0..n_params {
-            mu[i] = self.beta[beta_idx];
-            beta_idx += 1;
-
-            for (j, cov_name) in self.covariate_names.iter().enumerate() {
-                if self.covariate_mask[i][j] {
-                    let cov_value = covariates.get(cov_name).copied().unwrap_or(0.0);
-                    let reference = self.reference_values.get(cov_name).copied().unwrap_or(0.0);
-                    mu[i] += self.beta[beta_idx] * (cov_value - reference);
-                    beta_idx += 1;
-                }
-            }
-        }
-
-        mu
-    }
-
-    pub fn build_design_row(&self, covariates: &HashMap<String, f64>) -> Col<f64> {
-        let n_beta = self.n_beta();
-        let mut x = Col::zeros(n_beta);
-        let mut beta_idx = 0;
-
-        for i in 0..self.n_params() {
-            x[beta_idx] = 1.0;
-            beta_idx += 1;
-
-            for (j, cov_name) in self.covariate_names.iter().enumerate() {
-                if self.covariate_mask[i][j] {
-                    let cov_value = covariates.get(cov_name).copied().unwrap_or(0.0);
-                    let reference = self.reference_values.get(cov_name).copied().unwrap_or(0.0);
-                    x[beta_idx] = cov_value - reference;
-                    beta_idx += 1;
-                }
-            }
-        }
-
-        x
-    }
-
-    pub fn build_design_matrix(&self, all_covariates: &[HashMap<String, f64>]) -> Mat<f64> {
-        let n_subjects = all_covariates.len();
-        let n_params = self.n_params();
-        let n_beta = self.n_beta();
-        let n_rows = n_subjects * n_params;
-        let mut x = Mat::zeros(n_rows, n_beta);
-
-        for (subject_idx, covs) in all_covariates.iter().enumerate() {
-            let mut beta_idx = 0;
-            for param_idx in 0..n_params {
-                let row_idx = subject_idx * n_params + param_idx;
-                x[(row_idx, beta_idx)] = 1.0;
-                beta_idx += 1;
-
-                for (j, cov_name) in self.covariate_names.iter().enumerate() {
-                    if self.covariate_mask[param_idx][j] {
-                        let cov_value = covs.get(cov_name).copied().unwrap_or(0.0);
-                        let reference = self.reference_values.get(cov_name).copied().unwrap_or(0.0);
-                        x[(row_idx, beta_idx)] = cov_value - reference;
-                        beta_idx += 1;
-                    }
-                }
-            }
-        }
-
-        x
-    }
-
-    pub fn n_params(&self) -> usize {
-        self.param_names.len()
-    }
-
-    pub fn n_covariates(&self) -> usize {
-        self.covariate_names.len()
-    }
-
-    pub fn n_beta(&self) -> usize {
-        Self::count_beta_coefficients(&self.covariate_mask, self.param_names.len())
-    }
-
-    pub fn n_beta_estimated(&self) -> usize {
-        self.estimate_beta.iter().filter(|&&x| x).count()
-    }
-
-    pub fn param_names(&self) -> &[String] {
-        &self.param_names
-    }
-
-    pub fn covariate_names(&self) -> &[String] {
-        &self.covariate_names
-    }
-
-    pub fn covariate_mask(&self) -> &[Vec<bool>] {
-        &self.covariate_mask
-    }
-
-    pub fn beta(&self) -> &Col<f64> {
-        &self.beta
-    }
-
-    pub fn beta_mut(&mut self) -> &mut Col<f64> {
-        &mut self.beta
-    }
-
-    pub fn estimate_beta(&self) -> &[bool] {
-        &self.estimate_beta
-    }
-
-    pub fn intercept(&self, param_idx: usize) -> Option<f64> {
-        let beta_idx = self.intercept_beta_index(param_idx).ok()?;
-        Some(self.beta[beta_idx])
-    }
-
-    pub fn has_covariates(&self, param_idx: usize) -> bool {
-        param_idx < self.covariate_mask.len() && self.covariate_mask[param_idx].iter().any(|&x| x)
-    }
-
-    fn count_beta_coefficients(mask: &[Vec<bool>], n_params: usize) -> usize {
-        let mut count = n_params;
-        for row in mask {
-            count += row.iter().filter(|&&x| x).count();
-        }
-        count
-    }
-
-    fn intercept_beta_index(&self, param_idx: usize) -> Result<usize> {
-        if param_idx >= self.n_params() {
-            bail!("Parameter index {} out of range", param_idx);
-        }
-
-        let mut idx = 0;
-        for i in 0..param_idx {
-            idx += 1;
-            idx += self.covariate_mask[i].iter().filter(|&&x| x).count();
-        }
-        Ok(idx)
-    }
-
-    pub fn estimated_beta_indices(&self) -> Vec<usize> {
-        self.estimate_beta
-            .iter()
-            .enumerate()
-            .filter_map(|(i, &est)| if est { Some(i) } else { None })
-            .collect()
-    }
-}
-
-impl Default for CovariateModel {
-    fn default() -> Self {
-        Self {
-            param_names: Vec::new(),
-            covariate_names: Vec::new(),
-            covariate_mask: Vec::new(),
-            beta: Col::zeros(0),
-            estimate_beta: Vec::new(),
-            reference_values: HashMap::new(),
-        }
-    }
-}
-
-impl Serialize for CovariateModel {
-    fn serialize<S>(&self, serializer: S) -> std::result::Result<S::Ok, S::Error>
-    where
-        S: serde::Serializer,
-    {
-        use serde::ser::SerializeStruct;
-
-        let mut state = serializer.serialize_struct("CovariateModel", 6)?;
-        state.serialize_field("param_names", &self.param_names)?;
-        state.serialize_field("covariate_names", &self.covariate_names)?;
-        state.serialize_field("covariate_mask", &self.covariate_mask)?;
-        let beta_vec: Vec<f64> = (0..self.beta.nrows()).map(|i| self.beta[i]).collect();
-        state.serialize_field("beta", &beta_vec)?;
-        state.serialize_field("estimate_beta", &self.estimate_beta)?;
-        state.serialize_field("reference_values", &self.reference_values)?;
-        state.end()
-    }
-}
-
-#[cfg(test)]
-mod tests {
-    use super::*;
-
-    #[test]
-    fn test_intercept_only() {
-        let model = CovariateModel::intercept_only(vec!["CL", "V"]).unwrap();
-
-        assert_eq!(model.n_params(), 2);
-        assert_eq!(model.n_covariates(), 0);
-        assert_eq!(model.n_beta(), 2);
-
-        let mut model = model;
-        model.set_intercepts(&[5.0, 50.0]).unwrap();
-
-        let mu = model.compute_mu(&HashMap::new());
-        assert_eq!(mu[0], 5.0);
-        assert_eq!(mu[1], 50.0);
-    }
-
-    #[test]
-    fn test_with_covariates() {
-        let model = CovariateModel::new(
-            vec!["CL", "V"],
-            vec!["WT", "SEX"],
-            vec![vec![true, false], vec![true, true]],
-        )
-        .unwrap();
-
-        assert_eq!(model.n_params(), 2);
-        assert_eq!(model.n_covariates(), 2);
-        assert_eq!(model.n_beta(), 5);
-    }
-
-    #[test]
-    fn test_compute_mu() {
-        let mut model =
-            CovariateModel::new(vec!["CL", "V"], vec!["WT"], vec![vec![true], vec![true]]).unwrap();
-
-        model
-            .set_beta(Col::from_fn(4, |i| match i {
-                0 => 5.0,
-                1 => 0.1,
-                2 => 50.0,
-                3 => 1.0,
-                _ => 0.0,
-            }))
-            .unwrap();
-
-        let mut covs = HashMap::new();
-        covs.insert("WT".to_string(), 70.0);
-        let mu = model.compute_mu(&covs);
-
-        assert!((mu[0] - 12.0).abs() < 1e-10);
-        assert!((mu[1] - 120.0).abs() < 1e-10);
-    }
-
-    #[test]
-    fn test_centering() {
-        let mut model = CovariateModel::new(vec!["CL"], vec!["WT"], vec![vec![true]]).unwrap();
-        model.set_reference("WT", 70.0).unwrap();
-        model
-            .set_beta(Col::from_fn(2, |i| if i == 0 { 5.0 } else { 0.1 }))
-            .unwrap();
-
-        let mut covs = HashMap::new();
-        covs.insert("WT".to_string(), 70.0);
-        let mu = model.compute_mu(&covs);
-        assert!((mu[0] - 5.0).abs() < 1e-10);
-
-        covs.insert("WT".to_string(), 80.0);
-        let mu = model.compute_mu(&covs);
-        assert!((mu[0] - 6.0).abs() < 1e-10);
-    }
-
-    #[test]
-    fn test_from_saemix_matrix() {
-        let model = CovariateModel::from_saemix_matrix(
-            vec!["CL", "V"],
-            vec!["WT", "SEX"],
-            &[1.0, 0.0, 1.0, 1.0],
-        )
-        .unwrap();
-
-        assert!(model.covariate_mask[0][0]);
-        assert!(!model.covariate_mask[0][1]);
-        assert!(model.covariate_mask[1][0]);
-        assert!(model.covariate_mask[1][1]);
-    }
-
-    #[test]
-    fn test_fix_intercept() {
-        let mut model = CovariateModel::intercept_only(vec!["CL", "V"]).unwrap();
-
-        assert!(model.estimate_beta[0]);
-        model.fix_intercept(0).unwrap();
-        assert!(!model.estimate_beta[0]);
-        assert!(model.estimate_beta[1]);
-    }
-}
diff --git a/src/model/covariates.rs b/src/model/covariates.rs
deleted file mode 100644
index fb514add8..000000000
--- a/src/model/covariates.rs
+++ /dev/null
@@ -1,33 +0,0 @@
-use serde::Serialize;
-
-use crate::model::CovariateModel;
-
-#[allow(clippy::large_enum_variant)]
-#[derive(Debug, Clone, Default, Serialize)]
-pub enum CovariateSpec {
-    #[default]
-    InEquation,
-    Structured(CovariateEffectsSpec),
-}
-
-#[derive(Debug, Clone, Default, Serialize)]
-pub struct CovariateEffectsSpec {
-    pub subject_effects: Option<CovariateModel>,
-    pub occasion_effects: Option<CovariateModel>,
-}
-
-impl CovariateEffectsSpec {
-    pub fn subject_columns(&self) -> Vec<String> {
-        self.subject_effects
-            .as_ref()
-            .map(|model| model.covariate_names().to_vec())
-            .unwrap_or_default()
-    }
-
-    pub fn occasion_columns(&self) -> Vec<String> {
-        self.occasion_effects
-            .as_ref()
-            .map(|model| model.covariate_names().to_vec())
-            .unwrap_or_default()
-    }
-}
diff --git a/src/model/mod.rs b/src/model/mod.rs
index 23271eaf3..26b22787e 100644
--- a/src/model/mod.rs
+++ b/src/model/mod.rs
@@ -1,43 +1,30 @@
-use anyhow::{anyhow, bail, Result};
+use anyhow::Result;
 use pharmsol::equation::Equation;
 use pharmsol::{Analytical, ValidatedModelMetadata, ODE, SDE};
 
-pub mod covariate_model;
-pub mod covariates;
 pub mod metadata;
 pub mod parameter_space;
-pub mod variability;
 
-pub use covariate_model::CovariateModel;
-pub use covariates::{CovariateEffectsSpec, CovariateSpec};
 pub use metadata::ModelMetadata;
+// Re-exporting the new typestate parameter elements for easy access downstream
 pub use parameter_space::{
-    Parameter, ParameterDomain, ParameterSpace, ParameterTransform, ParameterVariability,
+    BoundedParameter, Parameter, ParameterMeta, ParameterScale, ParameterSpace, UnboundedParameter,
 };
-pub use variability::{CovarianceStructure, RandomEffectsSpec, VariabilityModel};
 
 #[derive(Debug, Clone)]
-pub struct ModelDefinition<E: Equation> {
+pub struct Model<E: Equation> {
     pub equation: E,
-    pub parameters: ParameterSpace,
-    pub variability: VariabilityModel,
-    pub covariates: CovariateSpec,
-    pub metadata: ModelMetadata,
+    // Note: No 'parameters' field here! It is now managed by EstimationProblem<E, F>.
+    // This struct is ready to hold Covariates or Variability specs if you add them later.
 }
 
-impl<E: Equation> ModelDefinition<E> {
-    pub fn builder(equation: E) -> ModelDefinitionBuilder<E> {
-        ModelDefinitionBuilder {
-            equation,
-            parameters: ParameterSpace::new(),
-            variability: Some(VariabilityModel::default()),
-            covariates: Some(CovariateSpec::InEquation),
-            metadata: Some(ModelMetadata::default()),
-        }
+impl<E: Equation> Model<E> {
+    pub fn builder(equation: E) -> ModelBuilder<E> {
+        ModelBuilder { equation }
     }
 }
 
-impl<E: EquationMetadataSource> ModelDefinition<E> {
+impl<E: EquationMetadataSource> Model<E> {
     pub fn parameter_count(&self) -> usize {
         self.equation
             .equation_metadata()
@@ -80,66 +67,51 @@ impl<E: EquationMetadataSource> ModelDefinition<E> {
     }
 }
 
-pub struct ModelDefinitionBuilder<E: Equation> {
+pub struct ModelBuilder<E: Equation> {
     equation: E,
-    parameters: ParameterSpace,
-    variability: Option<VariabilityModel>,
-    covariates: Option<CovariateSpec>,
-    metadata: Option<ModelMetadata>,
 }
 
-impl<E: Equation> ModelDefinitionBuilder<E> {
-    pub fn parameter(mut self, parameter: Parameter) -> Result<Self>
-    where
-        E: EquationMetadataSource,
-    {
-        validate_parameter(&self.equation, &self.parameters, &parameter)?;
-        self.parameters.push(parameter);
-        Ok(self)
-    }
-
-    pub fn variability(mut self, variability: VariabilityModel) -> Self {
-        self.variability = Some(variability);
-        self
+impl<E: Equation> ModelBuilder<E> {
+    pub fn build(self) -> Result<Model<E>> {
+        Ok(Model {
+            equation: self.equation,
+        })
     }
+}
 
-    pub fn covariates(mut self, covariates: CovariateSpec) -> Self {
-        self.covariates = Some(covariates);
-        self
+impl<E: EquationMetadataSource> ModelBuilder<E> {
+    pub(crate) fn parameter_index(&self, name: &str) -> Option<usize> {
+        self.equation
+            .equation_metadata()
+            .and_then(|metadata| metadata.parameter_index(name))
     }
 
-    pub fn metadata(mut self, metadata: ModelMetadata) -> Self {
-        self.metadata = Some(metadata);
-        self
+    pub(crate) fn parameter_names(&self) -> Vec<String> {
+        self.equation
+            .equation_metadata()
+            .map(|metadata| {
+                metadata
+                    .parameters()
+                    .iter()
+                    .map(|parameter| parameter.name().to_string())
+                    .collect()
+            })
+            .unwrap_or_default()
     }
 
-    pub fn build(self) -> Result<ModelDefinition<E>>
-    where
-        E: EquationMetadataSource,
-    {
-        let ModelDefinitionBuilder {
-            equation,
-            parameters,
-            variability,
-            covariates,
-            metadata,
-        } = self;
-
-        if parameters.is_empty() {
-            bail!("model parameters cannot be empty");
-        }
-
-        Ok(ModelDefinition {
-            equation,
-            parameters,
-            variability: variability.unwrap_or_default(),
-            covariates: covariates.unwrap_or_default(),
-            metadata: metadata.unwrap_or_default(),
-        })
+    pub(crate) fn output_names(&self) -> Vec<String> {
+        self.equation
+            .equation_metadata()
+            .map(|metadata| {
+                metadata
+                    .outputs()
+                    .iter()
+                    .map(|output| output.name().to_string())
+                    .collect()
+            })
+            .unwrap_or_default()
     }
-}
 
-impl<E: EquationMetadataSource> ModelDefinitionBuilder<E> {
     pub(crate) fn output_index(&self, name: &str) -> Option<usize> {
         self.equation.equation_metadata().and_then(|metadata| {
             metadata
@@ -150,93 +122,41 @@ impl<E: EquationMetadataSource> ModelDefinitionBuilder<E> {
     }
 }
 
-fn validate_parameter<E: EquationMetadataSource>(
-    equation: &E,
-    existing: &ParameterSpace,
-    parameter: &Parameter,
-) -> Result<()> {
-    let metadata = equation
-        .equation_metadata()
-        .ok_or_else(|| anyhow!("equation metadata is required to define model parameters"))?;
-
-    if parameter.name.trim().is_empty() {
-        bail!("model parameter name cannot be empty");
-    }
-
-    if metadata.parameter_index(&parameter.name).is_none() {
-        bail!("unknown equation parameter: {}", parameter.name);
-    }
-
-    if existing
-        .iter()
-        .any(|existing_parameter| existing_parameter.name == parameter.name)
-    {
-        bail!("duplicate model parameter: {}", parameter.name);
-    }
-
-    Ok(())
-}
-
 pub trait EquationMetadataSource: Equation {
     fn equation_metadata(&self) -> Option<&ValidatedModelMetadata>;
 }
 
-impl EquationMetadataSource for ODE {
-    fn equation_metadata(&self) -> Option<&ValidatedModelMetadata> {
-        self.metadata()
-    }
-}
-
-impl EquationMetadataSource for Analytical {
-    fn equation_metadata(&self) -> Option<&ValidatedModelMetadata> {
-        self.metadata()
-    }
-}
-
-impl EquationMetadataSource for SDE {
-    fn equation_metadata(&self) -> Option<&ValidatedModelMetadata> {
-        self.metadata()
-    }
+// Macro for standard pharmsol equations
+macro_rules! impl_metadata_opt {
+    ($($t:ty),+) => {
+        $(impl EquationMetadataSource for $t {
+            fn equation_metadata(&self) -> Option<&ValidatedModelMetadata> {
+                self.metadata()
+            }
+        })+
+    };
 }
-
-#[cfg(any(
-    feature = "dsl-jit",
-    all(feature = "dsl-aot", feature = "dsl-aot-load"),
-    all(
-        feature = "dsl-wasm",
-        not(all(target_arch = "wasm32", target_os = "unknown"))
-    )
-))]
-impl EquationMetadataSource for pharmsol::dsl::RuntimeOdeModel {
-    fn equation_metadata(&self) -> Option<&ValidatedModelMetadata> {
-        Some(self.metadata())
-    }
-}
-
-#[cfg(any(
-    feature = "dsl-jit",
-    all(feature = "dsl-aot", feature = "dsl-aot-load"),
-    all(
-        feature = "dsl-wasm",
-        not(all(target_arch = "wasm32", target_os = "unknown"))
-    )
-))]
-impl EquationMetadataSource for pharmsol::dsl::RuntimeAnalyticalModel {
-    fn equation_metadata(&self) -> Option<&ValidatedModelMetadata> {
-        Some(self.metadata())
-    }
-}
-
-#[cfg(any(
-    feature = "dsl-jit",
-    all(feature = "dsl-aot", feature = "dsl-aot-load"),
-    all(
-        feature = "dsl-wasm",
-        not(all(target_arch = "wasm32", target_os = "unknown"))
-    )
-))]
-impl EquationMetadataSource for pharmsol::dsl::RuntimeSdeModel {
-    fn equation_metadata(&self) -> Option<&ValidatedModelMetadata> {
-        Some(self.metadata())
-    }
+impl_metadata_opt!(ODE, Analytical, SDE);
+
+// Macro for runtime/JIT models
+macro_rules! impl_metadata_some {
+    ($($t:ty),+) => {
+        $(
+            #[cfg(any(
+                feature = "dsl-jit",
+                all(feature = "dsl-aot", feature = "dsl-aot-load"),
+                all(feature = "dsl-wasm", not(all(target_arch = "wasm32", target_os = "unknown")))
+            ))]
+            impl EquationMetadataSource for $t {
+                fn equation_metadata(&self) -> Option<&ValidatedModelMetadata> {
+                    Some(self.metadata())
+                }
+            }
+        )+
+    };
 }
+impl_metadata_some!(
+    pharmsol::dsl::RuntimeOdeModel,
+    pharmsol::dsl::RuntimeAnalyticalModel,
+    pharmsol::dsl::RuntimeSdeModel
+);
diff --git a/src/model/parameter_space.rs b/src/model/parameter_space.rs
index 1efeb676b..b92785bea 100644
--- a/src/model/parameter_space.rs
+++ b/src/model/parameter_space.rs
@@ -1,24 +1,21 @@
-use anyhow::{bail, Result};
 use serde::{Deserialize, Serialize};
 
-#[derive(Debug, Clone, Serialize, Deserialize, PartialEq)]
-pub struct ParameterSpace {
-    pub items: Vec<Parameter>,
+/// Ordered collection of parameters.
+///
+/// Use `ParameterSpace<BoundedParameter>` for non-parametric problems and
+/// `ParameterSpace<UnboundedParameter>` for parametric problems.
+#[derive(Debug, Clone, Default, Serialize, Deserialize, PartialEq)]
+pub struct ParameterSpace<T> {
+    pub items: Vec<T>,
 }
 
-impl ParameterSpace {
+impl<T> ParameterSpace<T> {
     pub fn new() -> Self {
         Self { items: Vec::new() }
     }
 
-    pub fn push(&mut self, item: Parameter) {
-        self.items.push(item);
-    }
-
-    #[allow(clippy::should_implement_trait)]
-    pub fn add(mut self, item: Parameter) -> Self {
-        self.push(item);
-        self
+    pub fn push(&mut self, item: impl Into<T>) {
+        self.items.push(item.into());
     }
 
     pub fn len(&self) -> usize {
@@ -29,123 +26,233 @@ impl ParameterSpace {
         self.items.is_empty()
     }
 
-    pub fn iter(&self) -> std::slice::Iter<'_, Parameter> {
+    pub fn iter(&self) -> std::slice::Iter<'_, T> {
         self.items.iter()
     }
+}
 
+impl<T: ParameterMeta> ParameterSpace<T> {
     pub fn names(&self) -> Vec<String> {
-        self.items.iter().map(|item| item.name.clone()).collect()
-    }
-
-    pub fn finite_ranges(&self) -> Result<Vec<(f64, f64)>> {
-        self.items
-            .iter()
-            .map(|parameter| match parameter.domain {
-                ParameterDomain::Bounded { lower, upper } => Ok((lower, upper)),
-                ParameterDomain::Positive {
-                    lower: Some(lower),
-                    upper: Some(upper),
-                }
-                | ParameterDomain::Unbounded {
-                    lower: Some(lower),
-                    upper: Some(upper),
-                } => Ok((lower, upper)),
-                _ => bail!(
-                    "nonparametric initialization requires finite lower/upper bounds for parameter '{}'",
-                    parameter.name
-                ),
-            })
-            .collect()
-    }
-}
-
-impl Default for ParameterSpace {
-    fn default() -> Self {
+        self.items.iter().map(|p| p.name().to_string()).collect()
+    }
+}
+
+/// Helpers for bounded parameter spaces.
+impl ParameterSpace<BoundedParameter> {
+    /// Creates an empty bounded parameter space (for non-parametric problems).
+    ///
+    /// Prefer this over `ParameterSpace::<BoundedParameter>::new()` to avoid the
+    /// turbofish:
+    ///
+    /// ```ignore
+    /// let space = ParameterSpace::bounded()
+    ///     .add("ke", 0.1, 1.0)
+    ///     .add("v", 1.0, 20.0);
+    /// ```
+    pub fn bounded() -> Self {
         Self::new()
     }
+
+    /// Adds a bounded parameter with the given `name`, `lower`, and `upper` bounds.
+    #[allow(clippy::should_implement_trait)]
+    pub fn add(mut self, name: impl Into<String>, lower: f64, upper: f64) -> Self {
+        self.items.push(BoundedParameter::new(name, lower, upper));
+        self
+    }
+
+    /// Returns `(lower, upper)` for each parameter.
+    pub fn finite_ranges(&self) -> Vec<(f64, f64)> {
+        self.items.iter().map(|p| (p.lower, p.upper)).collect()
+    }
 }
 
-impl From<&ParameterSpace> for ParameterSpace {
-    fn from(parameter_space: &ParameterSpace) -> Self {
-        parameter_space.clone()
+/// Helpers for unbounded parameter spaces.
+impl ParameterSpace<UnboundedParameter> {
+    /// Creates an empty unbounded parameter space (for parametric problems).
+    ///
+    /// Prefer this over `ParameterSpace::<UnboundedParameter>::new()` to avoid
+    /// the turbofish.
+    pub fn unbounded() -> Self {
+        Self::new()
+    }
+
+    /// Adds an unbounded parameter to the space.
+    #[allow(clippy::should_implement_trait)]
+    pub fn add(mut self, parameter: impl Into<UnboundedParameter>) -> Self {
+        self.items.push(parameter.into());
+        self
     }
 }
 
-fn default_estimate() -> bool {
-    true
+impl<T> FromIterator<T> for ParameterSpace<T> {
+    fn from_iter<I: IntoIterator<Item = T>>(iter: I) -> Self {
+        Self {
+            items: iter.into_iter().collect(),
+        }
+    }
+}
+
+impl<T> IntoIterator for ParameterSpace<T> {
+    type Item = T;
+    type IntoIter = std::vec::IntoIter<T>;
+
+    fn into_iter(self) -> Self::IntoIter {
+        self.items.into_iter()
+    }
+}
+
+impl<'a, T> IntoIterator for &'a ParameterSpace<T> {
+    type Item = &'a T;
+    type IntoIter = std::slice::Iter<'a, T>;
+
+    fn into_iter(self) -> Self::IntoIter {
+        self.items.iter()
+    }
+}
+
+/// Common metadata exposed by parameter types.
+pub trait ParameterMeta {
+    fn name(&self) -> &str;
 }
 
 #[derive(Debug, Clone, Serialize, Deserialize, PartialEq)]
-pub struct Parameter {
+pub struct BoundedParameter {
     pub name: String,
-    pub domain: ParameterDomain,
-    #[serde(default)]
-    pub transform: ParameterTransform,
-    #[serde(default)]
-    pub initial: Option<f64>,
-    #[serde(default = "default_estimate")]
-    pub estimate: bool,
-    #[serde(default)]
-    pub variability: ParameterVariability,
+    pub lower: f64,
+    pub upper: f64,
 }
 
-impl Parameter {
-    pub fn bounded(name: impl Into<String>, lower: f64, upper: f64) -> Self {
+impl BoundedParameter {
+    pub fn new(name: impl Into<String>, lower: f64, upper: f64) -> Self {
         Self {
             name: name.into(),
-            domain: ParameterDomain::Bounded { lower, upper },
-            transform: ParameterTransform::Identity,
+            lower,
+            upper,
+        }
+    }
+}
+
+impl ParameterMeta for BoundedParameter {
+    fn name(&self) -> &str {
+        &self.name
+    }
+}
+
+/// Converts a bounded parameter into a parametric parameter with logit scaling.
+impl From<BoundedParameter> for UnboundedParameter {
+    fn from(p: BoundedParameter) -> Self {
+        UnboundedParameter {
+            name: p.name,
+            scale: ParameterScale::Logit {
+                lower: p.lower,
+                upper: p.upper,
+            },
             initial: None,
             estimate: true,
-            variability: ParameterVariability::Subject,
         }
     }
+}
 
-    pub fn positive(name: impl Into<String>) -> Self {
+/// Parametric parameter with an optional scale transform.
+#[derive(Debug, Clone, Serialize, Deserialize, PartialEq)]
+pub struct UnboundedParameter {
+    pub name: String,
+    pub scale: ParameterScale,
+    pub initial: Option<f64>,
+    pub estimate: bool,
+}
+
+impl UnboundedParameter {
+    /// Creates a parameter with an explicit scale.
+    pub fn new(name: impl Into<String>, scale: ParameterScale) -> Self {
         Self {
             name: name.into(),
-            domain: ParameterDomain::Positive {
-                lower: Some(0.0),
-                upper: None,
-            },
-            transform: ParameterTransform::LogNormal,
+            scale,
             initial: None,
             estimate: true,
-            variability: ParameterVariability::Subject,
         }
     }
+
+    /// Creates a parameter on identity scale.
+    pub fn real(name: impl Into<String>) -> Self {
+        Self::new(name, ParameterScale::Identity)
+    }
+
+    /// Sets an initial value.
+    pub fn with_initial(mut self, value: f64) -> Self {
+        self.initial = Some(value);
+        self
+    }
 }
 
-#[derive(Debug, Clone, Serialize, Deserialize, PartialEq)]
-pub enum ParameterDomain {
-    Positive {
-        lower: Option<f64>,
-        upper: Option<f64>,
-    },
-    Unbounded {
-        lower: Option<f64>,
-        upper: Option<f64>,
-    },
-    Bounded {
-        lower: f64,
-        upper: f64,
-    },
-}
-
-#[derive(Debug, Clone, Serialize, Deserialize, PartialEq, Eq, Default)]
-pub enum ParameterTransform {
-    #[default]
+impl ParameterMeta for UnboundedParameter {
+    fn name(&self) -> &str {
+        &self.name
+    }
+}
+
+/// Scale transform for parametric parameters.
+#[derive(Debug, Clone, Copy, Serialize, Deserialize, PartialEq)]
+pub enum ParameterScale {
+    /// Identity transform.
     Identity,
-    LogNormal,
-    Probit,
-    Logit,
-}
-
-#[derive(Debug, Clone, Serialize, Deserialize, PartialEq, Eq, Default)]
-pub enum ParameterVariability {
-    FixedOnly,
-    #[default]
-    Subject,
-    Occasion,
-    SubjectAndOccasion,
+    /// Log transform.
+    Log,
+    /// Logistic transform on `(lower, upper)`.
+    Logit { lower: f64, upper: f64 },
+    /// Probit transform on `(lower, upper)`.
+    Probit { lower: f64, upper: f64 },
+}
+
+impl Default for ParameterScale {
+    fn default() -> Self {
+        ParameterScale::Identity
+    }
+}
+
+/// Entry point for building parameter declarations.
+///
+/// ```ignore
+/// use pmcore::prelude::*;
+///
+/// // Non-parametric: only bounded parameters are accepted.
+/// builder.parameter(Parameter::bounded("ke", 0.001, 3.0));
+///
+/// // Parametric: pick the scale explicitly.
+/// builder.parameter(Parameter::log("ke"));
+/// builder.parameter(Parameter::logit("frac", 0.0, 1.0));
+/// builder.parameter(Parameter::bounded("v", 25.0, 250.0)); // mapped to Logit
+/// ```
+pub struct Parameter;
+
+impl Parameter {
+    /// Creates a bounded parameter.
+    pub fn bounded(name: impl Into<String>, lower: f64, upper: f64) -> BoundedParameter {
+        BoundedParameter::new(name, lower, upper)
+    }
+
+    /// Creates a parametric parameter on identity scale.
+    pub fn real(name: impl Into<String>) -> UnboundedParameter {
+        UnboundedParameter::real(name)
+    }
+
+    /// Creates a parametric parameter with an explicit scale.
+    pub fn scaled(name: impl Into<String>, scale: ParameterScale) -> UnboundedParameter {
+        UnboundedParameter::new(name, scale)
+    }
+
+    /// Creates a parametric parameter on log scale.
+    pub fn log(name: impl Into<String>) -> UnboundedParameter {
+        UnboundedParameter::new(name, ParameterScale::Log)
+    }
+
+    /// Creates a parametric parameter on logit scale.
+    pub fn logit(name: impl Into<String>, lower: f64, upper: f64) -> UnboundedParameter {
+        UnboundedParameter::new(name, ParameterScale::Logit { lower, upper })
+    }
+
+    /// Creates a parametric parameter on probit scale.
+    pub fn probit(name: impl Into<String>, lower: f64, upper: f64) -> UnboundedParameter {
+        UnboundedParameter::new(name, ParameterScale::Probit { lower, upper })
+    }
 }
diff --git a/src/model/variability.rs b/src/model/variability.rs
deleted file mode 100644
index 20cab7638..000000000
--- a/src/model/variability.rs
+++ /dev/null
@@ -1,28 +0,0 @@
-use serde::{Deserialize, Serialize};
-
-#[derive(Debug, Clone, Serialize, Deserialize, PartialEq, Eq, Default)]
-pub struct VariabilityModel {
-    pub subject: RandomEffectsSpec,
-    pub occasion: Option<RandomEffectsSpec>,
-}
-
-#[derive(Debug, Clone, Serialize, Deserialize, PartialEq, Eq)]
-pub struct RandomEffectsSpec {
-    pub enabled_for: Vec<bool>,
-    pub covariance: CovarianceStructure,
-}
-
-impl Default for RandomEffectsSpec {
-    fn default() -> Self {
-        Self {
-            enabled_for: Vec::new(),
-            covariance: CovarianceStructure::Diagonal,
-        }
-    }
-}
-
-#[derive(Debug, Clone, Serialize, Deserialize, PartialEq, Eq)]
-pub enum CovarianceStructure {
-    Diagonal,
-    Full,
-}
diff --git a/src/output/file.rs b/src/output/file.rs
deleted file mode 100644
index e48210f4a..000000000
--- a/src/output/file.rs
+++ /dev/null
@@ -1,45 +0,0 @@
-use anyhow::{Context, Result};
-use std::fs::{create_dir_all, File, OpenOptions};
-use std::path::{Path, PathBuf};
-
-/// Contains all the necessary information of an output file.
-#[derive(Debug)]
-pub struct OutputFile {
-    file: File,
-    relative_path: PathBuf,
-}
-
-impl OutputFile {
-    pub fn new(folder: &str, file_name: &str) -> Result<Self> {
-        let relative_path = Path::new(folder).join(file_name);
-
-        if let Some(parent) = relative_path.parent() {
-            create_dir_all(parent)
-                .with_context(|| format!("Failed to create directories for {:?}", parent))?;
-        }
-
-        let file = OpenOptions::new()
-            .write(true)
-            .create(true)
-            .truncate(true)
-            .open(&relative_path)
-            .with_context(|| format!("Failed to open file: {:?}", relative_path))?;
-
-        Ok(Self {
-            file,
-            relative_path,
-        })
-    }
-
-    pub fn file(&self) -> &File {
-        &self.file
-    }
-
-    pub fn file_owned(self) -> File {
-        self.file
-    }
-
-    pub fn relative_path(&self) -> &Path {
-        &self.relative_path
-    }
-}
diff --git a/src/output/mod.rs b/src/output/mod.rs
deleted file mode 100644
index ffa994aa2..000000000
--- a/src/output/mod.rs
+++ /dev/null
@@ -1,8 +0,0 @@
-mod file;
-pub mod logging;
-pub mod nonparametric;
-pub mod shared;
-pub mod writer;
-
-pub use file::OutputFile;
-pub use writer::write_result;
diff --git a/src/output/nonparametric.rs b/src/output/nonparametric.rs
deleted file mode 100644
index f0fa9d38b..000000000
--- a/src/output/nonparametric.rs
+++ /dev/null
@@ -1,44 +0,0 @@
-use anyhow::Result;
-
-use crate::estimation::nonparametric::NonparametricWorkspace;
-use crate::output::shared::shared_output_file_names;
-
-pub(crate) fn output_file_names<E: pharmsol::Equation>(
-    result: &NonparametricWorkspace<E>,
-) -> Vec<String> {
-    let mut files = shared_output_file_names();
-    files.extend(
-        ["iterations.csv", "theta.csv", "posterior.csv"]
-            .into_iter()
-            .map(str::to_string),
-    );
-
-    let has_covariates = result.data().subjects().iter().any(|subject| {
-        subject
-            .occasions()
-            .iter()
-            .any(|occasion| !occasion.covariates().covariates().is_empty())
-    });
-    if has_covariates {
-        files.push("covariates.csv".to_string());
-    }
-
-    files.sort();
-    files.dedup();
-    files
-}
-
-pub fn write_nonparametric_outputs<E: pharmsol::Equation>(
-    result: &mut NonparametricWorkspace<E>,
-) -> Result<()> {
-    let parameter_names = result.get_theta().parameters().names();
-    result
-        .cycle_log()
-        .write(result.output_folder(), &parameter_names)?;
-    result.write_theta()?;
-    result.write_covariates()?;
-    result.write_posterior()?;
-    let (idelta, tad) = result.prediction_interval();
-    result.calculate_predictions(idelta, tad)?;
-    Ok(())
-}
diff --git a/src/output/shared.rs b/src/output/shared.rs
deleted file mode 100644
index 4836d84fc..000000000
--- a/src/output/shared.rs
+++ /dev/null
@@ -1,113 +0,0 @@
-use anyhow::Result;
-use csv::WriterBuilder;
-use serde::Serialize;
-
-use crate::algorithms::Algorithm;
-use crate::api::{OutputPlan, RuntimeOptions};
-use crate::output::OutputFile;
-use crate::results::{DiagnosticsBundle, FitSummary};
-
-pub(crate) fn shared_output_file_names() -> Vec<String> {
-    vec![
-        "settings.json",
-        "summary.json",
-        "summary.csv",
-        "diagnostics.json",
-        "predictions.csv",
-    ]
-    .into_iter()
-    .map(str::to_string)
-    .collect()
-}
-
-#[derive(Debug, Clone, Serialize)]
-pub(crate) struct RunConfiguration {
-    pub algorithm: Algorithm,
-    pub output: OutputPlan,
-    pub runtime: RuntimeOptions,
-    pub parameter_names: Vec<String>,
-}
-
-impl RunConfiguration {
-    pub(crate) fn new(
-        algorithm: Algorithm,
-        output: &OutputPlan,
-        runtime: &RuntimeOptions,
-        parameter_names: Vec<String>,
-    ) -> Self {
-        Self {
-            algorithm,
-            output: output.clone(),
-            runtime: runtime.clone(),
-            parameter_names,
-        }
-    }
-
-    pub(crate) fn output_path(&self) -> &str {
-        self.output.path.as_deref().unwrap_or("outputs/")
-    }
-
-    pub(crate) fn should_write_outputs(&self) -> bool {
-        self.output.write
-    }
-}
-
-pub(crate) fn write_settings(folder: &str, configuration: &RunConfiguration) -> Result<()> {
-    let outputfile = OutputFile::new(folder, "settings.json")?;
-    let mut file = outputfile.file_owned();
-    let serialized = serde_json::to_string_pretty(configuration)?;
-    std::io::Write::write_all(&mut file, serialized.as_bytes())?;
-    Ok(())
-}
-
-pub fn write_summary(folder: &str, summary: &FitSummary) -> Result<()> {
-    let outputfile = OutputFile::new(folder, "summary.json")?;
-    let mut file = outputfile.file_owned();
-    let serialized = serde_json::to_string_pretty(summary)?;
-    std::io::Write::write_all(&mut file, serialized.as_bytes())?;
-
-    let outputfile = OutputFile::new(folder, "summary.csv")?;
-    let mut writer = WriterBuilder::new()
-        .has_headers(true)
-        .from_writer(outputfile.file_owned());
-    writer.write_record(["metric", "value"])?;
-    writer.write_record([
-        "objective_function",
-        &summary.objective_function.to_string(),
-    ])?;
-    writer.write_record(["converged", &summary.converged.to_string()])?;
-    writer.write_record(["iterations", &summary.iterations.to_string()])?;
-    writer.write_record(["subject_count", &summary.subject_count.to_string()])?;
-    writer.write_record(["observation_count", &summary.observation_count.to_string()])?;
-    writer.write_record(["parameter_count", &summary.parameter_count.to_string()])?;
-    writer.write_record(["algorithm", &summary.algorithm])?;
-    writer.flush()?;
-
-    Ok(())
-}
-
-pub fn write_diagnostics(folder: &str, diagnostics: &DiagnosticsBundle) -> Result<()> {
-    let outputfile = OutputFile::new(folder, "diagnostics.json")?;
-    let mut file = outputfile.file_owned();
-    let serialized = serde_json::to_string_pretty(diagnostics)?;
-    std::io::Write::write_all(&mut file, serialized.as_bytes())?;
-    Ok(())
-}
-
-pub fn write_csv_rows<T: Serialize>(
-    folder: &str,
-    file_name: &str,
-    rows: impl IntoIterator<Item = T>,
-) -> Result<()> {
-    let outputfile = OutputFile::new(folder, file_name)?;
-    let mut writer = WriterBuilder::new()
-        .has_headers(true)
-        .from_writer(outputfile.file_owned());
-
-    for row in rows {
-        writer.serialize(row)?;
-    }
-
-    writer.flush()?;
-    Ok(())
-}
diff --git a/src/output/writer.rs b/src/output/writer.rs
deleted file mode 100644
index 60ef89a9f..000000000
--- a/src/output/writer.rs
+++ /dev/null
@@ -1,72 +0,0 @@
-use anyhow::Result;
-use pharmsol::{Censor, Equation};
-use serde::Serialize;
-
-use crate::estimation::nonparametric as np_estimation;
-use crate::estimation::nonparametric::NonparametricWorkspace;
-use crate::output::{nonparametric as np_output, shared};
-use crate::results::FitResult;
-use crate::results::{nonparametric_diagnostics, FitSummary};
-
-#[derive(Debug, Clone, Serialize)]
-struct SharedPredictionRow {
-    id: String,
-    time: f64,
-    outeq: usize,
-    block: usize,
-    obs: Option<f64>,
-    cens: Censor,
-    pred_population: f64,
-    pred_individual: f64,
-    residual_population: Option<f64>,
-    residual_individual: Option<f64>,
-    source_method: String,
-}
-
-pub fn write_result<E: Equation>(result: &mut FitResult<E>) -> Result<()> {
-    match result {
-        FitResult::Nonparametric(inner) => write_nonparametric_result(inner)?,
-    }
-
-    Ok(())
-}
-
-pub fn write_nonparametric_result<E: Equation>(
-    result: &mut NonparametricWorkspace<E>,
-) -> Result<()> {
-    if !result.should_write_outputs() {
-        return Ok(());
-    }
-
-    let folder = result.output_folder().to_string();
-    shared::write_settings(&folder, result.run_configuration())?;
-    shared::write_summary(&folder, &nonparametric_summary(result))?;
-    shared::write_diagnostics(&folder, &nonparametric_diagnostics(result))?;
-    np_output::write_nonparametric_outputs(result)?;
-
-    if let Some(predictions) = result.predictions() {
-        let rows = predictions
-            .predictions()
-            .iter()
-            .map(|row| SharedPredictionRow {
-                id: row.id().to_string(),
-                time: row.time(),
-                outeq: row.outeq(),
-                block: row.block(),
-                obs: row.obs(),
-                cens: row.censoring(),
-                pred_population: row.pop_mean(),
-                pred_individual: row.post_mean(),
-                residual_population: row.obs().map(|obs| obs - row.pop_mean()),
-                residual_individual: row.obs().map(|obs| obs - row.post_mean()),
-                source_method: "nonparametric".to_string(),
-            });
-        shared::write_csv_rows(&folder, "predictions.csv", rows)?;
-    }
-
-    Ok(())
-}
-
-fn nonparametric_summary<E: Equation>(result: &NonparametricWorkspace<E>) -> FitSummary {
-    np_estimation::fit_summary(result)
-}
diff --git a/src/results/artifacts.rs b/src/results/artifacts.rs
deleted file mode 100644
index abd817c3e..000000000
--- a/src/results/artifacts.rs
+++ /dev/null
@@ -1,88 +0,0 @@
-use std::fs;
-use std::path::Path;
-
-use pharmsol::Equation;
-use serde::{Deserialize, Serialize};
-
-use crate::estimation::nonparametric::NonparametricWorkspace;
-use crate::output::shared::shared_output_file_names;
-
-#[derive(Debug, Clone, Default, Serialize, Deserialize, PartialEq, Eq)]
-pub struct ArtifactIndex {
-    pub files: Vec<String>,
-    pub expected_files: Vec<String>,
-    pub missing_files: Vec<String>,
-    pub shared_expected_files: Vec<String>,
-    pub method_specific_expected_files: Vec<String>,
-}
-
-pub(crate) fn nonparametric_artifacts<E: Equation>(
-    result: &NonparametricWorkspace<E>,
-) -> ArtifactIndex {
-    artifact_index(
-        result.output_folder(),
-        result.should_write_outputs(),
-        crate::output::nonparametric::output_file_names(result),
-    )
-}
-
-fn artifact_index(
-    folder: &str,
-    should_write_outputs: bool,
-    mut expected_files: Vec<String>,
-) -> ArtifactIndex {
-    if !should_write_outputs {
-        return ArtifactIndex::default();
-    }
-
-    expected_files.sort();
-    expected_files.dedup();
-
-    let shared_output_files = shared_output_file_names();
-    let shared_expected_files = expected_files
-        .iter()
-        .filter(|file| shared_output_files.contains(*file))
-        .cloned()
-        .collect::<Vec<_>>();
-    let method_specific_expected_files = expected_files
-        .iter()
-        .filter(|file| !shared_output_files.contains(*file))
-        .cloned()
-        .collect::<Vec<_>>();
-
-    let path = Path::new(folder);
-    if !path.exists() {
-        return ArtifactIndex {
-            files: Vec::new(),
-            missing_files: expected_files.clone(),
-            expected_files,
-            shared_expected_files,
-            method_specific_expected_files,
-        };
-    }
-
-    let mut files = expected_files
-        .iter()
-        .filter(|file| {
-            fs::metadata(path.join(file))
-                .map(|meta| meta.is_file())
-                .unwrap_or(false)
-        })
-        .cloned()
-        .collect::<Vec<_>>();
-    files.sort();
-
-    let missing_files = expected_files
-        .iter()
-        .filter(|file| !files.contains(file))
-        .cloned()
-        .collect();
-
-    ArtifactIndex {
-        files,
-        expected_files,
-        missing_files,
-        shared_expected_files,
-        method_specific_expected_files,
-    }
-}
diff --git a/src/results/diagnostics.rs b/src/results/diagnostics.rs
deleted file mode 100644
index cfeaf2b5a..000000000
--- a/src/results/diagnostics.rs
+++ /dev/null
@@ -1,59 +0,0 @@
-use std::collections::BTreeMap;
-
-use pharmsol::Equation;
-use serde::{Deserialize, Serialize};
-
-use crate::estimation::nonparametric::NonparametricWorkspace;
-
-#[derive(Debug, Clone, Default, Serialize, Deserialize, PartialEq, Eq)]
-pub struct DiagnosticsBundle {
-    pub warnings: Vec<String>,
-    pub deferred_features: Vec<String>,
-    pub convergence_notes: Vec<String>,
-    pub estimator_metadata: BTreeMap<String, String>,
-}
-
-pub(crate) fn nonparametric_diagnostics<E: Equation>(
-    result: &NonparametricWorkspace<E>,
-) -> DiagnosticsBundle {
-    let mut convergence_notes = Vec::new();
-    if result.converged() {
-        convergence_notes.push("Estimator reported convergence.".to_string());
-    } else {
-        convergence_notes.push("Estimator stopped without convergence.".to_string());
-    }
-
-    let status = result
-        .cycle_log()
-        .cycles()
-        .last()
-        .map(|cycle| format!("{:?}", cycle.status()))
-        .unwrap_or_else(|| "Continue".to_string());
-
-    let mut estimator_metadata = BTreeMap::new();
-    estimator_metadata.insert("algorithm".to_string(), format!("{:?}", result.algorithm()));
-    estimator_metadata.insert("status".to_string(), status);
-    estimator_metadata.insert(
-        "outputs_requested".to_string(),
-        result.should_write_outputs().to_string(),
-    );
-    estimator_metadata.insert(
-        "support_point_count".to_string(),
-        result.get_theta().nspp().to_string(),
-    );
-    estimator_metadata.insert(
-        "prediction_cache".to_string(),
-        if result.predictions().is_some() {
-            "available".to_string()
-        } else {
-            "not_materialized".to_string()
-        },
-    );
-
-    DiagnosticsBundle {
-        warnings: Vec::new(),
-        deferred_features: Vec::new(),
-        convergence_notes,
-        estimator_metadata,
-    }
-}
diff --git a/src/results/fit_result.rs b/src/results/fit_result.rs
index c38532429..3e8470b52 100644
--- a/src/results/fit_result.rs
+++ b/src/results/fit_result.rs
@@ -1,74 +1,63 @@
-use anyhow::Result;
 use pharmsol::Equation;
 
-use crate::estimation::nonparametric;
-use crate::estimation::nonparametric::NonparametricWorkspace;
-use crate::results::{
-    nonparametric_artifacts, nonparametric_diagnostics, nonparametric_predictions, ArtifactIndex,
-    DiagnosticsBundle, FitSummary, IndividualSummary, PopulationSummary, PredictionsBundle,
-};
+use crate::estimation::nonparametric::NonParametricResult;
+use crate::results::{FitSummary, IndividualSummary, PopulationSummary};
+
+/// A shared trait for the output of any estimation algorithm.
+pub trait FitResult {
+    fn objf(&self) -> f64;
+    fn converged(&self) -> bool;
+    fn summary(&self) -> FitSummary;
+    fn population_summary(&self) -> PopulationSummary;
+    fn individual_summaries(&self) -> Vec<IndividualSummary>;
+}
 
+// Placeholder for your future Parametric implementation
 #[derive(Debug)]
-pub enum FitResult<E: Equation> {
-    Nonparametric(NonparametricWorkspace<E>),
+#[allow(unused)]
+//TODO: Implement ParametricResult
+pub struct ParametricResult<E: Equation> {
+    _phantom: std::marker::PhantomData<E>,
 }
 
-impl<E: Equation> FitResult<E> {
-    pub fn objf(&self) -> f64 {
-        match self {
-            Self::Nonparametric(result) => result.objf(),
-        }
+impl<E: Equation> FitResult for ParametricResult<E> {
+    fn objf(&self) -> f64 {
+        unimplemented!("Parametric result not yet implemented")
     }
-
-    pub fn converged(&self) -> bool {
-        match self {
-            Self::Nonparametric(result) => result.converged(),
-        }
+    fn converged(&self) -> bool {
+        unimplemented!()
     }
-
-    pub(crate) fn write_outputs(&mut self) -> Result<()> {
-        crate::output::write_result(self)
+    fn summary(&self) -> FitSummary {
+        unimplemented!()
     }
-
-    pub fn summary(&self) -> FitSummary {
-        match self {
-            Self::Nonparametric(result) => nonparametric::fit_summary(result),
-        }
+    fn population_summary(&self) -> PopulationSummary {
+        unimplemented!()
     }
-
-    pub fn population_summary(&self) -> PopulationSummary {
-        match self {
-            Self::Nonparametric(result) => nonparametric::population_summary(result),
-        }
+    fn individual_summaries(&self) -> Vec<IndividualSummary> {
+        unimplemented!()
     }
+}
+
+use crate::estimation::nonparametric;
 
-    pub fn individual_summaries(&self) -> Vec<IndividualSummary> {
-        match self {
-            Self::Nonparametric(result) => nonparametric::individual_summaries(result),
-        }
+impl<E: Equation> FitResult for NonParametricResult<E> {
+    fn objf(&self) -> f64 {
+        self.objf() // Assuming the struct has this native method
     }
 
-    pub fn diagnostics(&self) -> DiagnosticsBundle {
-        match self {
-            Self::Nonparametric(result) => nonparametric_diagnostics(result),
-        }
+    fn converged(&self) -> bool {
+        self.converged() // Assuming the struct has this native method
     }
 
-    pub fn predictions(&self) -> PredictionsBundle {
-        match self {
-            Self::Nonparametric(result) => nonparametric_predictions(result),
-        }
+    fn summary(&self) -> FitSummary {
+        nonparametric::fit_summary(self)
     }
 
-    pub fn artifacts(&self) -> ArtifactIndex {
-        match self {
-            Self::Nonparametric(result) => nonparametric_artifacts(result),
-        }
+    fn population_summary(&self) -> PopulationSummary {
+        nonparametric::population_summary(self)
     }
 
-    pub fn as_nonparametric(&self) -> Option<&NonparametricWorkspace<E>> {
-        match self {
-            Self::Nonparametric(result) => Some(result),
-        }
+    fn individual_summaries(&self) -> Vec<IndividualSummary> {
+        nonparametric::individual_summaries(self)
     }
 }
diff --git a/src/results/mod.rs b/src/results/mod.rs
index 3ff684906..5a0bed71e 100644
--- a/src/results/mod.rs
+++ b/src/results/mod.rs
@@ -1,15 +1,7 @@
-mod artifacts;
-mod diagnostics;
 mod fit_result;
-mod predictions;
+
 mod summary;
 
-pub use artifacts::ArtifactIndex;
-pub use diagnostics::DiagnosticsBundle;
-pub use fit_result::FitResult;
-pub use predictions::PredictionsBundle;
-pub use summary::{FitSummary, IndividualSummary, ParameterSummary, PopulationSummary};
+pub use fit_result::{FitResult, ParametricResult};
 
-pub(crate) use artifacts::nonparametric_artifacts;
-pub(crate) use diagnostics::nonparametric_diagnostics;
-pub(crate) use predictions::nonparametric_predictions;
+pub use summary::{FitSummary, IndividualSummary, ParameterSummary, PopulationSummary};
diff --git a/src/results/predictions.rs b/src/results/predictions.rs
deleted file mode 100644
index 8f3118999..000000000
--- a/src/results/predictions.rs
+++ /dev/null
@@ -1,38 +0,0 @@
-use pharmsol::Equation;
-use serde::{Deserialize, Serialize};
-
-use crate::estimation::nonparametric::NonparametricWorkspace;
-use crate::results::nonparametric_artifacts;
-
-#[derive(Debug, Clone, Default, Serialize, Deserialize, PartialEq, Eq)]
-pub struct PredictionsBundle {
-    pub available: bool,
-    pub row_count: Option<usize>,
-    pub source: Option<String>,
-    pub artifact: Option<String>,
-}
-
-pub(crate) fn nonparametric_predictions<E: Equation>(
-    result: &NonparametricWorkspace<E>,
-) -> PredictionsBundle {
-    let artifact = nonparametric_artifacts(result)
-        .files
-        .into_iter()
-        .find(|file| file == "predictions.csv");
-
-    if let Some(predictions) = result.predictions() {
-        return PredictionsBundle {
-            available: true,
-            row_count: Some(predictions.predictions().len()),
-            source: Some("in_memory".to_string()),
-            artifact,
-        };
-    }
-
-    PredictionsBundle {
-        available: artifact.is_some(),
-        row_count: None,
-        source: artifact.as_ref().map(|_| "artifact".to_string()),
-        artifact,
-    }
-}
diff --git a/src/results/summary.rs b/src/results/summary.rs
index 1ad5de295..c6769f38e 100644
--- a/src/results/summary.rs
+++ b/src/results/summary.rs
@@ -8,7 +8,6 @@ pub struct FitSummary {
     pub subject_count: usize,
     pub observation_count: usize,
     pub parameter_count: usize,
-    pub algorithm: String,
 }
 
 #[derive(Debug, Clone, Serialize, Deserialize, PartialEq)]
diff --git a/tests/acceptance_baseline_tests.rs b/tests/acceptance_baseline_tests.rs
deleted file mode 100644
index 933c3198e..000000000
--- a/tests/acceptance_baseline_tests.rs
+++ /dev/null
@@ -1,76 +0,0 @@
-use anyhow::Result;
-use pmcore::prelude::*;
-
-fn assert_close(actual: f64, expected: f64, tolerance: f64, label: &str) {
-    let difference = (actual - expected).abs();
-    assert!(
-        difference <= tolerance,
-        "{label}: expected {expected}, got {actual}, abs diff {difference} > {tolerance}"
-    );
-}
-
-fn bimodal_ode_equation() -> equation::ODE {
-    ode! {
-        name: "acceptance_baseline_bimodal_ke",
-        params: [ke, v],
-        states: [central],
-        outputs: [1],
-        routes: [
-            infusion(1) -> central,
-        ],
-        diffeq: |x, _t, dx| {
-            dx[central] = -ke * x[central];
-        },
-        out: |x, _t, y| {
-            y[1] = x[central] / v;
-        },
-    }
-    .with_solver(OdeSolver::ExplicitRk(ExplicitRkTableau::Tsit45))
-}
-
-fn bimodal_data() -> Result<Data> {
-    Ok(data::read_pmetrics("examples/bimodal_ke/bimodal_ke.csv")?)
-}
-
-#[test]
-fn test_acceptance_baseline_npag_bimodal_ke() -> Result<()> {
-    let result = EstimationProblem::builder(bimodal_ode_equation(), bimodal_data()?)
-        .parameter(Parameter::bounded("ke", 0.001, 3.0))?
-        .parameter(Parameter::bounded("v", 25.0, 250.0))?
-        .method(Npag::new())
-        .error(
-            "1",
-            AssayErrorModel::additive(ErrorPoly::new(0.0, 0.5, 0.0, 0.0), 0.0),
-        )?
-        .progress(false)
-        .fit()?;
-    let summary = result.summary();
-    let population = result.population_summary();
-    let result = result
-        .as_nonparametric()
-        .expect("NPAG acceptance baseline should yield a nonparametric result");
-
-    // This is the canonical rewrite-blocking nonparametric baseline for the bimodal_ke example.
-    assert_close(
-        summary.objective_function,
-        -425.60904902364695,
-        1e-6,
-        "npag.objf",
-    );
-    assert!(summary.converged);
-    assert_eq!(summary.iterations, 288);
-    assert_eq!(result.get_theta().nspp(), 46);
-    assert_close(
-        population.parameters[0].mean,
-        0.187047284678325,
-        1e-6,
-        "npag.ke.mean",
-    );
-    assert_close(
-        population.parameters[1].mean,
-        107.94241284196241,
-        1e-6,
-        "npag.v.mean",
-    );
-    Ok(())
-}
diff --git a/tests/api_smoke_tests.rs b/tests/api_smoke_tests.rs
deleted file mode 100644
index db82f8655..000000000
--- a/tests/api_smoke_tests.rs
+++ /dev/null
@@ -1,169 +0,0 @@
-use anyhow::Result;
-use pharmsol::{AssayErrorModel, ErrorPoly};
-use pmcore::prelude::*;
-
-fn simple_equation() -> equation::ODE {
-    equation::ODE::new(
-        |x, p, _t, dx, b, _rateiv, _cov| {
-            fetch_params!(p, ke);
-            dx[0] = -ke * x[0] + b[0];
-        },
-        |_p, _t, _cov| lag! {},
-        |_p, _t, _cov| fa! {},
-        |_p, _t, _cov, _x| {},
-        |x, p, _t, _cov, y| {
-            fetch_params!(p, v);
-            y[0] = x[0] / v;
-        },
-    )
-    .with_nstates(1)
-    .with_ndrugs(1)
-    .with_nout(1)
-    .with_metadata(
-        equation::metadata::new("simple_equation")
-            .parameters(["ke", "v"])
-            .states(["central"])
-            .outputs(["0"])
-            .route(equation::Route::bolus("0").to_state("central")),
-    )
-    .expect("metadata attachment should validate")
-}
-
-fn simple_data() -> Data {
-    let subject = Subject::builder("1")
-        .bolus(0.0, 100.0, 0)
-        .observation(1.0, 10.0, 0)
-        .observation(2.0, 8.0, 0)
-        .build();
-
-    Data::new(vec![subject])
-}
-
-fn metadata_equation() -> equation::ODE {
-    ode! {
-        name: "metadata_derived_outputs",
-        params: [ke, v],
-        states: [central],
-        outputs: [1],
-        routes: [
-            infusion(1) -> central,
-        ],
-        diffeq: |x, _t, dx| {
-            dx[central] = -ke * x[central];
-        },
-        out: |x, _t, y| {
-            y[1] = x[central] / v;
-        },
-    }
-}
-
-#[test]
-fn test_model_definition_builder() -> Result<()> {
-    let model = ModelDefinition::builder(simple_equation())
-        .parameter(Parameter::bounded("ke", 0.1, 1.0))?
-        .parameter(Parameter::bounded("v", 1.0, 20.0))?
-        .build()?;
-
-    assert_eq!(model.parameters.len(), 2);
-    assert_eq!(model.parameter_count(), 2);
-    assert_eq!(model.parameter_name(0), Some("ke"));
-    assert_eq!(model.output_count(), 1);
-    assert_eq!(model.output_index("0"), Some(0));
-    Ok(())
-}
-
-#[test]
-fn test_model_definition_rejects_unknown_parameter_name() {
-    let err = ModelDefinition::builder(simple_equation())
-        .parameter(Parameter::bounded("clearance", 0.1, 1.0))
-        .err()
-        .expect("parameter name should be validated against equation metadata");
-
-    assert!(err
-        .to_string()
-        .contains("unknown equation parameter: clearance"));
-}
-
-#[test]
-fn test_model_definition_derives_observations_from_equation_metadata() -> Result<()> {
-    let model = ModelDefinition::builder(metadata_equation())
-        .parameter(Parameter::bounded("ke", 0.1, 1.0))?
-        .parameter(Parameter::bounded("v", 1.0, 20.0))?
-        .build()?;
-
-    assert_eq!(model.output_count(), 1);
-    assert_eq!(model.output_name(0), Some("1"));
-    assert_eq!(model.output_index("1"), Some(0));
-
-    Ok(())
-}
-
-#[test]
-fn test_unified_fit_nonparametric_smoke() -> Result<()> {
-    let assay_error = AssayErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0);
-    let result = EstimationProblem::builder(simple_equation(), simple_data())
-        .parameter(Parameter::bounded("ke", 0.1, 1.0))?
-        .parameter(Parameter::bounded("v", 1.0, 20.0))?
-        .method(Npag::new())
-        .error("0", assay_error)?
-        .progress(false)
-        .fit()?;
-
-    assert!(result.objf().is_finite());
-    assert_eq!(result.summary().parameter_count, 2);
-    assert_eq!(result.population_summary().parameters.len(), 2);
-    assert_eq!(result.individual_summaries().len(), 1);
-    Ok(())
-}
-
-#[test]
-fn test_problem_compile_preserves_runtime_configuration() -> Result<()> {
-    let assay_error = AssayErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0);
-    let convergence = ConvergenceOptions {
-        likelihood: 1e-5,
-        pyl: 5e-3,
-        eps: 2e-3,
-    };
-    let tuning = AlgorithmTuning {
-        min_distance: 2e-4,
-        nm_steps: 222,
-        tolerance: 3e-6,
-        saem: SaemConfig {
-            k1_iterations: 111,
-            k2_iterations: 22,
-            ..SaemConfig::default()
-        },
-    };
-
-    let compiled = EstimationProblem::builder(simple_equation(), simple_data())
-        .parameter(Parameter::bounded("ke", 0.1, 1.0))?
-        .parameter(Parameter::bounded("v", 1.0, 20.0))?
-        .method(Npag::new())
-        .error("0", assay_error)?
-        .cache(false)
-        .progress(false)
-        .idelta(0.5)
-        .tad(24.0)
-        .convergence(convergence)
-        .tuning(tuning)
-        .build()?
-        .compile()?;
-
-    assert_eq!(compiled.algorithm(), Algorithm::NPAG);
-    assert!(!compiled.output_plan().write);
-    assert_eq!(compiled.runtime_options().cycles, 7);
-    assert!(!compiled.runtime_options().cache);
-    assert!(!compiled.runtime_options().progress);
-    assert_eq!(compiled.runtime_options().idelta, 0.5);
-    assert_eq!(compiled.runtime_options().tad, 24.0);
-
-    assert_eq!(compiled.runtime_options().convergence.likelihood, 1e-5);
-    assert_eq!(compiled.runtime_options().convergence.pyl, 5e-3);
-    assert_eq!(compiled.runtime_options().convergence.eps, 2e-3);
-    assert_eq!(compiled.runtime_options().tuning.min_distance, 2e-4);
-    assert_eq!(compiled.runtime_options().tuning.nm_steps, 222);
-    assert_eq!(compiled.runtime_options().tuning.tolerance, 3e-6);
-    assert_eq!(compiled.runtime_options().tuning.saem.k1_iterations, 111);
-    assert_eq!(compiled.runtime_options().tuning.saem.k2_iterations, 22);
-    Ok(())
-}
diff --git a/tests/bestdose_tests.rs b/tests/bestdose_tests.rs
index 270d09fbe..ee4a268b2 100644
--- a/tests/bestdose_tests.rs
+++ b/tests/bestdose_tests.rs
@@ -1,16 +1,22 @@
 use anyhow::Result;
 use pmcore::bestdose::{BestDoseConfig, BestDosePosterior, BestDoseProblem, DoseRange, Target};
 use pmcore::estimation::nonparametric::{Theta, Weights};
+use pmcore::model::{BoundedParameter, ParameterSpace};
 use pmcore::prelude::*;
 
-fn pk_parameter_space(ke_lower: f64, ke_upper: f64, v_lower: f64, v_upper: f64) -> ParameterSpace {
-    ParameterSpace::new()
-        .add(Parameter::bounded("ke", ke_lower, ke_upper))
-        .add(Parameter::bounded("v", v_lower, v_upper))
+fn pk_parameter_space(
+    ke_lower: f64,
+    ke_upper: f64,
+    v_lower: f64,
+    v_upper: f64,
+) -> ParameterSpace<BoundedParameter> {
+    ParameterSpace::<BoundedParameter>::new()
+        .add("ke", ke_lower, ke_upper)
+        .add("v", v_lower, v_upper)
 }
 
 fn bestdose_config(
-    params: &ParameterSpace,
+    params: &ParameterSpace<BoundedParameter>,
     error_models: AssayErrorModels,
     refinement_cycles: usize,
     prediction_interval: f64,
diff --git a/tests/compile_layer_tests.rs b/tests/compile_layer_tests.rs
deleted file mode 100644
index 46f6900bb..000000000
--- a/tests/compile_layer_tests.rs
+++ /dev/null
@@ -1,141 +0,0 @@
-use anyhow::Result;
-use pharmsol::{AssayErrorModel, ErrorPoly};
-use pmcore::prelude::*;
-
-fn simple_equation() -> equation::ODE {
-    equation::ODE::new(
-        |x, p, _t, dx, b, _rateiv, _cov| {
-            fetch_params!(p, ke);
-            dx[0] = -ke * x[0] + b[0];
-        },
-        |_p, _t, _cov| lag! {},
-        |_p, _t, _cov| fa! {},
-        |_p, _t, _cov, _x| {},
-        |x, p, _t, _cov, y| {
-            fetch_params!(p, v);
-            y[0] = x[0] / v;
-        },
-    )
-    .with_nstates(1)
-    .with_ndrugs(1)
-    .with_nout(1)
-    .with_metadata(
-        equation::metadata::new("compile_layer")
-            .parameters(["ke", "v"])
-            .states(["central"])
-            .outputs(["0"])
-            .route(equation::Route::bolus("0").to_state("central")),
-    )
-    .expect("metadata attachment should validate")
-}
-
-fn multi_subject_data() -> Data {
-    let first = Subject::builder("1")
-        .bolus(0.0, 100.0, 0)
-        .observation(1.0, 10.0, 0)
-        .observation(2.0, 8.0, 0)
-        .build();
-
-    let second = Subject::builder("2")
-        .bolus(0.0, 100.0, 0)
-        .observation(1.0, 9.0, 0)
-        .observation(3.0, 7.0, 0)
-        .build();
-
-    Data::new(vec![first, second])
-}
-
-fn structured_covariate_data() -> Data {
-    let subject = Subject::builder("1")
-        .covariate("wt", 0.0, 70.0)
-        .covariate("study_day", 0.0, 1.0)
-        .bolus(0.0, 100.0, 0)
-        .observation(1.0, 10.0, 0)
-        .reset()
-        .covariate("wt", 0.0, 72.0)
-        .covariate("study_day", 0.0, 2.0)
-        .bolus(0.0, 100.0, 0)
-        .observation(1.5, 8.0, 0)
-        .build();
-
-    Data::new(vec![subject])
-}
-
-fn simple_problem() -> Result<EstimationProblem<equation::ODE>> {
-    let assay_error = AssayErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0);
-    EstimationProblem::builder(simple_equation(), multi_subject_data())
-        .parameter(Parameter::bounded("ke", 0.1, 1.0))?
-        .parameter(Parameter::bounded("v", 1.0, 20.0))?
-        .method(Npag::new())
-        .error("0", assay_error)?
-        .build()
-}
-
-#[test]
-fn test_compile_problem_builds_indexes() -> Result<()> {
-    let compiled = simple_problem()?.compile()?;
-
-    assert_eq!(compiled.design.subject_count(), 2);
-    assert_eq!(compiled.design.occasion_count(), 2);
-    assert_eq!(compiled.observation_index.len(), 4);
-    assert_eq!(compiled.design.parameter_names, vec!["ke", "v"]);
-    Ok(())
-}
-
-#[test]
-fn test_compile_problem_builds_algorithm_settings() -> Result<()> {
-    let compiled = simple_problem()?.compile()?;
-
-    assert_eq!(compiled.algorithm(), Algorithm::NPAG);
-    assert_eq!(compiled.design.parameter_names.len(), 2);
-    assert!(!compiled.output_plan().write);
-    Ok(())
-}
-
-#[test]
-fn test_compile_problem_extracts_structured_covariate_values() -> Result<()> {
-    let assay_error = AssayErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0);
-    let compiled = EstimationProblem::builder(simple_equation(), structured_covariate_data())
-        .parameter(Parameter::bounded("ke", 0.1, 1.0))?
-        .parameter(Parameter::bounded("v", 1.0, 20.0))?
-        .covariates(CovariateSpec::Structured(CovariateEffectsSpec {
-            subject_effects: Some(CovariateModel::new(
-                vec!["ke", "v"],
-                vec!["wt"],
-                vec![vec![true], vec![false]],
-            )?),
-            occasion_effects: Some(CovariateModel::new(
-                vec!["ke", "v"],
-                vec!["study_day"],
-                vec![vec![true], vec![false]],
-            )?),
-        }))
-        .method(Npag::new())
-        .error("0", assay_error)?
-        .build()?
-        .compile()?;
-
-    assert_eq!(
-        compiled.design.structured_covariates.subject_columns,
-        vec!["wt"]
-    );
-    assert_eq!(
-        compiled.design.structured_covariates.occasion_columns,
-        vec!["study_day"]
-    );
-    assert_eq!(compiled.design.structured_covariates.subject_rows.len(), 1);
-    assert_eq!(compiled.design.structured_covariates.occasion_rows.len(), 2);
-    assert_eq!(
-        compiled.design.structured_covariates.subject_rows[0].values,
-        vec![Some(70.0)]
-    );
-    assert_eq!(
-        compiled.design.structured_covariates.occasion_rows[0].values,
-        vec![Some(1.0)]
-    );
-    assert_eq!(
-        compiled.design.structured_covariates.occasion_rows[1].values,
-        vec![Some(2.0)]
-    );
-    Ok(())
-}
diff --git a/tests/nonparametric_engine_tests.rs b/tests/nonparametric_engine_tests.rs
deleted file mode 100644
index ba2b12ef5..000000000
--- a/tests/nonparametric_engine_tests.rs
+++ /dev/null
@@ -1,62 +0,0 @@
-use anyhow::Result;
-use pharmsol::{AssayErrorModel, ErrorPoly};
-use pmcore::prelude::*;
-
-fn simple_equation() -> equation::ODE {
-    equation::ODE::new(
-        |x, p, _t, dx, b, _rateiv, _cov| {
-            fetch_params!(p, ke);
-            dx[0] = -ke * x[0] + b[0];
-        },
-        |_p, _t, _cov| lag! {},
-        |_p, _t, _cov| fa! {},
-        |_p, _t, _cov, _x| {},
-        |x, p, _t, _cov, y| {
-            fetch_params!(p, v);
-            y[0] = x[0] / v;
-        },
-    )
-    .with_nstates(1)
-    .with_ndrugs(1)
-    .with_nout(1)
-    .with_metadata(
-        equation::metadata::new("nonparametric_engine")
-            .parameters(["ke", "v"])
-            .states(["central"])
-            .outputs(["0"])
-            .route(equation::Route::bolus("0").to_state("central")),
-    )
-    .expect("metadata attachment should validate")
-}
-
-fn simple_data() -> Data {
-    let subject = Subject::builder("1")
-        .bolus(0.0, 100.0, 0)
-        .observation(1.0, 10.0, 0)
-        .observation(2.0, 8.0, 0)
-        .build();
-
-    Data::new(vec![subject])
-}
-
-#[test]
-fn test_nonparametric_engine_returns_workspace() -> Result<()> {
-    let assay_error = AssayErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0);
-    let compiled = EstimationProblem::builder(simple_equation(), simple_data())
-        .parameter(Parameter::bounded("ke", 0.1, 1.0))?
-        .parameter(Parameter::bounded("v", 1.0, 20.0))?
-        .method(Npag::new())
-        .error("0", assay_error)?
-        .progress(false)
-        .build()?
-        .compile()?;
-
-    let workspace = NonparametricEngine::fit(compiled)?;
-    assert!(workspace.objf().is_finite());
-    assert_eq!(workspace.get_theta().parameters().len(), 2);
-
-    let fit_result = workspace.into_fit_result();
-    assert_eq!(fit_result.population_summary().parameters.len(), 2);
-    assert_eq!(fit_result.individual_summaries().len(), 1);
-    Ok(())
-}
diff --git a/tests/onecomp.rs b/tests/onecomp.rs
index 7064c8d53..85c6a27be 100644
--- a/tests/onecomp.rs
+++ b/tests/onecomp.rs
@@ -52,24 +52,27 @@ fn test_one_compartment_npag() -> Result<()> {
 
     let data = data::Data::new(subjects);
 
-    let result = EstimationProblem::builder(eq, data)
-        .parameter(Parameter::bounded("ke", 0.1, 1.0))?
-        .parameter(Parameter::bounded("v", 1.0, 20.0))?
-        .method(Npag::new())
-        .error(
-            "0",
-            AssayErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0),
-        )?
-        .prior(Prior::sobol(64, 22))
-        .fit()?;
-    let result = result
-        .as_nonparametric()
-        .expect("NPAG should yield a nonparametric result");
+    let parameters = ParameterSpace::<BoundedParameter>::new()
+        .add("ke", 0.1, 1.0)
+        .add("v", 1.0, 20.0);
+
+    let prior = Theta::sobol(&parameters, 100)?;
+    let error_models = AssayErrorModels::new().add(
+        "0",
+        AssayErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0),
+    )?;
+    let result = EstimationProblem::nonparametric(eq, data, prior, error_models)?
+        .fit_with(NonParametricAlgorithm::npag())?;
 
     // Check the results
-    assert_eq!(result.cycles(), 32);
+    assert_eq!(result.cycles(), 31);
     assert!(result.objf() - 565.7749 < 0.01);
 
+    // The prior is preserved on the result and is distinct from the optimized
+    // solution (which is condensed to far fewer support points).
+    assert_eq!(result.prior().nspp(), 100);
+    assert!(result.get_theta().nspp() < result.prior().nspp());
+
     Ok(())
 }
 
@@ -116,19 +119,16 @@ fn test_one_compartment_npod() -> Result<()> {
 
     let data = data::Data::new(subjects);
 
-    let result = EstimationProblem::builder(eq, data)
-        .parameter(Parameter::bounded("ke", 0.1, 1.0))?
-        .parameter(Parameter::bounded("v", 1.0, 20.0))?
-        .method(Npod::new())
-        .error(
-            "0",
-            AssayErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0),
-        )?
-        .prior(Prior::sobol(64, 22))
-        .fit()?;
-    let result = result
-        .as_nonparametric()
-        .expect("NPOD should yield a nonparametric result");
+    let parameters = ParameterSpace::<BoundedParameter>::new()
+        .add("ke", 0.1, 1.0)
+        .add("v", 1.0, 20.0);
+    let prior = Theta::sobol_default(&parameters)?;
+    let error_models = AssayErrorModels::new().add(
+        "0",
+        AssayErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0),
+    )?;
+    let result = EstimationProblem::nonparametric(eq, data, prior, error_models)?
+        .fit_with(NonParametricAlgorithm::npod())?;
 
     // Check the results
     assert_eq!(result.cycles(), 11);
@@ -180,25 +180,26 @@ fn test_one_compartment_postprob() -> Result<()> {
 
     let data = data::Data::new(subjects);
 
-    let result = EstimationProblem::builder(eq, data)
-        .parameter(Parameter::bounded("ke", 0.1, 1.0))?
-        .parameter(Parameter::bounded("v", 1.0, 20.0))?
-        .method(PostProb::new())
-        .error(
-            "0",
-            AssayErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0),
-        )?
-        .prior(Prior::sobol(64, 22))
-        .fit()?;
-    let result = result
-        .as_nonparametric()
-        .expect("POSTPROB should yield a nonparametric result");
+    // Generate a prior distribution to test against
+    let parameters = ParameterSpace::<BoundedParameter>::new()
+        .add("ke", 0.1, 1.0)
+        .add("v", 1.0, 20.0);
+
+    let theta = Theta::sobol(&parameters, 100)?;
+
+    let error_models = AssayErrorModels::new().add(
+        "0",
+        AssayErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0),
+    )?;
+    let result = EstimationProblem::nonparametric(eq, data, theta.clone(), error_models)?
+        .fit_with(NonParametricAlgorithm::npmap())?;
 
     // Check the results
     assert_eq!(result.cycles(), 0);
 
-    // Should be 64 points in theta (no change in points)
-    assert_eq!(result.get_theta().nspp(), 64);
+    // Should be 100 points in theta (no change in points)
+    assert_eq!(result.get_theta().nspp(), theta.nspp());
 
     Ok(())
 }
+
diff --git a/tests/results_summary_tests.rs b/tests/results_summary_tests.rs
index 44dd13fe8..6edcc8a89 100644
--- a/tests/results_summary_tests.rs
+++ b/tests/results_summary_tests.rs
@@ -42,13 +42,11 @@ fn simple_data() -> Data {
 #[test]
 fn test_nonparametric_fit_result_summary_surface() -> Result<()> {
     let assay_error = AssayErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0);
-    let result = EstimationProblem::builder(simple_equation(), simple_data())
-        .parameter(Parameter::bounded("ke", 0.1, 1.0))?
-        .parameter(Parameter::bounded("v", 1.0, 20.0))?
-        .method(Npag::new())
-        .error("0", assay_error)?
-        .progress(false)
-        .fit()?;
+    let parameters = ParameterSpace::bounded().add("ke", 0.1, 1.0).add("v", 1.0, 20.0);
+    let prior = Theta::sobol_default(&parameters)?;
+    let error_models = AssayErrorModels::new().add("0", assay_error)?;
+    let result = EstimationProblem::nonparametric(simple_equation(), simple_data(), prior, error_models)?
+        .fit_with(NonParametricAlgorithm::npag())?;
 
     let summary = result.summary();
 
@@ -57,20 +55,7 @@ fn test_nonparametric_fit_result_summary_surface() -> Result<()> {
     assert_eq!(summary.observation_count, 2);
     assert_eq!(result.population_summary().parameters.len(), 2);
     assert_eq!(result.individual_summaries().len(), 1);
-    let diagnostics = result.diagnostics();
-    assert_eq!(
-        diagnostics.estimator_metadata.get("algorithm"),
-        Some(&"NPAG".to_string())
-    );
-    assert_eq!(
-        diagnostics.estimator_metadata.get("outputs_requested"),
-        Some(&"false".to_string())
-    );
-    assert!(!diagnostics.convergence_notes.is_empty());
-    assert!(diagnostics.deferred_features.is_empty());
-    let predictions = result.predictions();
-    assert!(!predictions.available);
-    assert!(result.artifacts().files.is_empty());
-    assert!(result.artifacts().expected_files.is_empty());
+
     Ok(())
 }
+