From 788c1e031900e2d5b20368612f422ec02c3f125e Mon Sep 17 00:00:00 2001
From: Markus <66058642+mhovd@users.noreply.github.com>
Date: Tue, 17 Mar 2026 08:38:44 +0100
Subject: [PATCH 1/9] WIP

---
 Cargo.toml                 |  4 +--
 src/algorithms/mod.rs      | 11 +++---
 src/algorithms/npag.rs     |  8 ++---
 src/algorithms/npod.rs     | 69 ++++++++++++++++++--------------------
 src/algorithms/postprob.rs |  4 +--
 src/bestdose/mod.rs        |  2 +-
 src/bestdose/posterior.rs  |  4 +--
 src/routines/output/mod.rs | 32 ++++--------------
 src/routines/settings.rs   | 16 ++++-----
 src/structs/psi.rs         |  4 +--
 10 files changed, 65 insertions(+), 89 deletions(-)

diff --git a/Cargo.toml b/Cargo.toml
index 65ef68a6b..5dfb2883d 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -16,7 +16,7 @@ exclude = [".github/*", ".vscode/*"]
 
 [dependencies]
 csv = "1.3.1"
-ndarray = { version = "0.16.1", features = ["rayon"] }
+ndarray = { version = "0.17.2", features = ["rayon"] }
 serde = "1.0.188"
 serde_json = "1.0.66"
 sobol_burley = "0.5.0"
@@ -29,7 +29,7 @@ tracing-subscriber = { version = "0.3.19", features = [
 ] }
 faer = "0.23.1"
 faer-ext = { version = "0.7.1", features = ["nalgebra", "ndarray"] }
-pharmsol = "=0.22.1"
+pharmsol = "=0.23.0"
 rand = "0.9.0"
 anyhow = "1.0.100"
 rayon = "1.10.0"
diff --git a/src/algorithms/mod.rs b/src/algorithms/mod.rs
index ab5c16912..3636184ef 100644
--- a/src/algorithms/mod.rs
+++ b/src/algorithms/mod.rs
@@ -9,7 +9,7 @@ use anyhow::Context;
 use anyhow::Result;
 use faer_ext::IntoNdarray;
 use ndarray::parallel::prelude::{IntoParallelIterator, ParallelIterator};
-use ndarray::{Array, ArrayBase, Dim, OwnedRepr};
+
 use npag::*;
 use npod::NPOD;
 use pharmsol::prelude::{data::Data, simulator::Equation};
@@ -61,10 +61,11 @@ pub trait Algorithms<E: Equation + Send + 'static>: Sync + Send + 'static {
             );
         }
 
-        let (_, col) = psi.dim();
-        let ecol: ArrayBase<OwnedRepr<f64>, Dim<[usize; 1]>> = Array::ones(col);
-        let plam = psi.dot(&ecol);
-        let w = 1. / &plam;
+        let (row, col) = psi.dim();
+        let plam: Vec<f64> = (0..row)
+            .map(|i| (0..col).map(|j| psi[(i, j)]).sum::<f64>())
+            .collect();
+        let w: Vec<f64> = plam.iter().map(|&x| 1.0 / x).collect();
 
         // Get the index of each element in `w` that is NaN or infinite
         let indices: Vec<usize> = w
diff --git a/src/algorithms/npag.rs b/src/algorithms/npag.rs
index 68ed04693..ac5a4b94e 100644
--- a/src/algorithms/npag.rs
+++ b/src/algorithms/npag.rs
@@ -13,11 +13,11 @@ use crate::structs::weights::Weights;
 use anyhow::bail;
 use anyhow::Result;
 use pharmsol::prelude::{
-    data::{Data, ErrorModels},
+    data::{AssayErrorModels, Data},
     simulator::Equation,
 };
 
-use pharmsol::prelude::ErrorModel;
+use pharmsol::prelude::AssayErrorModel;
 
 use crate::routines::initialization;
 
@@ -43,7 +43,7 @@ pub struct NPAG<E: Equation + Send + 'static> {
     f1: f64,
     cycle: usize,
     gamma_delta: Vec<f64>,
-    error_models: ErrorModels,
+    error_models: AssayErrorModels,
     status: Status,
     cycle_log: CycleLog,
     data: Data,
@@ -134,7 +134,7 @@ impl<E: Equation + Send + 'static> Algorithms<E> for NPAG<E> {
         tracing::debug!("Support points: {}", self.theta.nspp());
 
         self.error_models.iter().for_each(|(outeq, em)| {
-            if ErrorModel::None == *em {
+            if AssayErrorModel::None == *em {
                 return;
             }
             tracing::debug!(
diff --git a/src/algorithms/npod.rs b/src/algorithms/npod.rs
index ed962d971..9fde065df 100644
--- a/src/algorithms/npod.rs
+++ b/src/algorithms/npod.rs
@@ -20,17 +20,14 @@ use pharmsol::SppOptimizer;
 
 use anyhow::bail;
 use anyhow::Result;
-use faer_ext::IntoNdarray;
-use pharmsol::{prelude::ErrorModel, ErrorModels};
+use pharmsol::{prelude::AssayErrorModel, AssayErrorModels};
 use pharmsol::{
     prelude::{data::Data, simulator::Equation},
     Subject,
 };
 
-use ndarray::{
-    parallel::prelude::{IntoParallelRefMutIterator, ParallelIterator},
-    Array, Array1, ArrayBase, Dim, OwnedRepr,
-};
+use ndarray::Array1;
+use rayon::prelude::{IntoParallelRefMutIterator, ParallelIterator};
 
 const THETA_F: f64 = 1e-2;
 const THETA_D: f64 = 1e-4;
@@ -45,7 +42,7 @@ pub struct NPOD<E: Equation + Send + 'static> {
     objf: f64,
     cycle: usize,
     gamma_delta: Vec<f64>,
-    error_models: ErrorModels,
+    error_models: AssayErrorModels,
     converged: bool,
     status: Status,
     cycle_log: CycleLog,
@@ -155,7 +152,7 @@ impl<E: Equation + Send + 'static> Algorithms<E> for NPOD<E> {
         tracing::info!("Objective function = {:.4}", -2.0 * self.objf);
         tracing::debug!("Support points: {}", self.theta.nspp());
         self.error_models.iter().for_each(|(outeq, em)| {
-            if ErrorModel::None == *em {
+            if AssayErrorModel::None == *em {
                 return;
             }
             tracing::debug!(
@@ -207,7 +204,7 @@ impl<E: Equation + Send + 'static> Algorithms<E> for NPOD<E> {
     }
 
     fn estimation(&mut self) -> Result<()> {
-        let error_model: ErrorModels = self.error_models.clone();
+        let error_model: AssayErrorModels = self.error_models.clone();
 
         self.psi = calculate_psi(
             &self.equation,
@@ -295,7 +292,7 @@ impl<E: Equation + Send + 'static> Algorithms<E> for NPOD<E> {
             .clone()
             .iter_mut()
             .filter_map(|(outeq, em)| {
-                if *em == ErrorModel::None || em.is_factor_fixed().unwrap_or(true) {
+                if *em == AssayErrorModel::None || em.is_factor_fixed().unwrap_or(true) {
                     None
                 } else {
                     Some((outeq, em))
@@ -367,13 +364,12 @@ impl<E: Equation + Send + 'static> Algorithms<E> for NPOD<E> {
     }
 
     fn expansion(&mut self) -> Result<()> {
-        // If no stop signal, add new point to theta based on the optimization of the D function
-        let psi = self.psi().matrix().as_ref().into_ndarray().to_owned();
-        let w: Array1<f64> = self.w.clone().iter().collect();
-        let pyl = psi.dot(&w);
+        // Compute pyl = psi * w using faer native operations
+        let pyl_col = self.psi().matrix().as_ref() * self.w.weights().as_ref();
+        let pyl: Array1<f64> = pyl_col.iter().copied().collect();
 
         // Add new point to theta based on the optimization of the D function
-        let error_model: ErrorModels = self.error_models.clone();
+        let error_model: AssayErrorModels = self.error_models.clone();
 
         let mut candididate_points: Vec<Array1<f64>> = Vec::default();
         for spp in self.theta.matrix().row_iter() {
@@ -400,33 +396,32 @@ impl<E: Equation + Send + 'static> Algorithms<E> for NPOD<E> {
 
 impl<E: Equation + Send + 'static> NPOD<E> {
     fn validate_psi(&mut self) -> Result<()> {
-        let mut psi = self.psi().matrix().as_ref().into_ndarray().to_owned();
+        let mut psi = self.psi().matrix().to_owned();
         // First coerce all NaN and infinite in psi to 0.0
-        if psi.iter().any(|x| x.is_nan() || x.is_infinite()) {
-            tracing::warn!("Psi contains NaN or Inf values, coercing to 0.0");
-            for i in 0..psi.nrows() {
-                for j in 0..psi.ncols() {
-                    let val = psi.get_mut((i, j)).unwrap();
-                    if val.is_nan() || val.is_infinite() {
-                        *val = 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 the sum of each column in psi
-        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<_>>();
+        // 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() {
diff --git a/src/algorithms/postprob.rs b/src/algorithms/postprob.rs
index 496b36e28..de2fa0307 100644
--- a/src/algorithms/postprob.rs
+++ b/src/algorithms/postprob.rs
@@ -10,7 +10,7 @@ use crate::{
 use anyhow::{Context, Result};
 
 use pharmsol::prelude::{
-    data::{Data, ErrorModels},
+    data::{AssayErrorModels, Data},
     simulator::Equation,
 };
 
@@ -32,7 +32,7 @@ pub struct POSTPROB<E: Equation + Send + 'static> {
     data: Data,
     settings: Settings,
     cyclelog: CycleLog,
-    error_models: ErrorModels,
+    error_models: AssayErrorModels,
 }
 
 impl<E: Equation + Send + 'static> Algorithms<E> for POSTPROB<E> {
diff --git a/src/bestdose/mod.rs b/src/bestdose/mod.rs
index 56626aef0..f9588dc33 100644
--- a/src/bestdose/mod.rs
+++ b/src/bestdose/mod.rs
@@ -525,7 +525,7 @@ fn calculate_posterior_density(
     population_weights: &Weights,
     past_data: Option<&Subject>,
     eq: &ODE,
-    error_models: &ErrorModels,
+    error_models: &AssayErrorModels,
     settings: &Settings,
 ) -> Result<(Theta, Weights, Weights, Subject)> {
     match past_data {
diff --git a/src/bestdose/posterior.rs b/src/bestdose/posterior.rs
index 9109f65c1..e86a4ae07 100644
--- a/src/bestdose/posterior.rs
+++ b/src/bestdose/posterior.rs
@@ -94,7 +94,7 @@ pub fn npagfull11_filter(
     population_weights: &Weights,
     past_data: &Data,
     eq: &ODE,
-    error_models: &ErrorModels,
+    error_models: &AssayErrorModels,
 ) -> Result<(Theta, Weights, Weights)> {
     tracing::info!("Stage 1.1: NPAGFULL11 Bayesian filtering");
 
@@ -312,7 +312,7 @@ pub fn calculate_two_step_posterior(
     population_weights: &Weights,
     past_data: &Data,
     eq: &ODE,
-    error_models: &ErrorModels,
+    error_models: &AssayErrorModels,
     settings: &Settings,
 ) -> Result<(Theta, Weights, Weights)> {
     tracing::info!("=== STAGE 1: Posterior Density Calculation ===");
diff --git a/src/routines/output/mod.rs b/src/routines/output/mod.rs
index 2a9e43d69..3835904d2 100644
--- a/src/routines/output/mod.rs
+++ b/src/routines/output/mod.rs
@@ -9,8 +9,6 @@ use crate::structs::theta::Theta;
 use crate::structs::weights::Weights;
 use anyhow::{bail, Context, Result};
 use csv::WriterBuilder;
-use faer::linalg::zip::IntoView;
-use faer_ext::IntoNdarray;
 use ndarray::{Array, Array1, Array2, Axis};
 use pharmsol::prelude::data::*;
 use pharmsol::prelude::simulator::Equation;
@@ -171,22 +169,11 @@ impl<E: Equation> NPResult<E> {
             post_median: f64,
         }
 
-        let theta: Array2<f64> = self
-            .theta
-            .matrix()
-            .clone()
-            .as_mut()
-            .into_ndarray()
-            .to_owned();
-        let w: Array1<f64> = self
-            .w
-            .weights()
-            .clone()
-            .into_view()
-            .iter()
-            .cloned()
-            .collect();
-        let psi: Array2<f64> = self.psi.matrix().as_ref().into_ndarray().to_owned();
+        let tm = self.theta.matrix();
+        let theta = Array2::from_shape_fn((tm.nrows(), tm.ncols()), |(i, j)| tm[(i, j)]);
+        let w: Array1<f64> = self.w.iter().collect();
+        let pm = self.psi.matrix();
+        let psi = Array2::from_shape_fn((pm.nrows(), pm.ncols()), |(i, j)| pm[(i, j)]);
 
         let (post_mean, post_median) = posterior_mean_median(&theta, &psi, &w)
             .context("Failed to calculate posterior mean and median")?;
@@ -299,14 +286,7 @@ impl<E: Equation> NPResult<E> {
         tracing::debug!("Writing population parameter distribution...");
 
         let theta = &self.theta;
-        let w: Vec<f64> = self
-            .w
-            .weights()
-            .clone()
-            .into_view()
-            .iter()
-            .cloned()
-            .collect();
+        let w: Vec<f64> = self.w.to_vec();
 
         if w.len() != theta.matrix().nrows() {
             bail!(
diff --git a/src/routines/settings.rs b/src/routines/settings.rs
index 3471a683f..45383b287 100644
--- a/src/routines/settings.rs
+++ b/src/routines/settings.rs
@@ -2,7 +2,7 @@ use crate::algorithms::Algorithm;
 use crate::routines::initialization::Prior;
 use crate::routines::output::OutputFile;
 use anyhow::{bail, Result};
-use pharmsol::prelude::data::ErrorModels;
+use pharmsol::prelude::data::AssayErrorModels;
 
 use serde::{Deserialize, Serialize};
 use serde_json;
@@ -18,7 +18,7 @@ pub struct Settings {
     /// Parameters to be estimated
     pub(crate) parameters: Parameters,
     /// Defines the error models and polynomials to be used
-    pub(crate) errormodels: ErrorModels,
+    pub(crate) errormodels: AssayErrorModels,
     /// Configuration for predictions
     pub(crate) predictions: Predictions,
     /// Configuration for logging
@@ -48,7 +48,7 @@ impl Settings {
         &self.parameters
     }
 
-    pub fn errormodels(&self) -> &ErrorModels {
+    pub fn errormodels(&self) -> &AssayErrorModels {
         &self.errormodels
     }
 
@@ -441,7 +441,7 @@ impl Default for Output {
 pub struct SettingsBuilder<State> {
     config: Option<Config>,
     parameters: Option<Parameters>,
-    errormodels: Option<ErrorModels>,
+    errormodels: Option<AssayErrorModels>,
     predictions: Option<Predictions>,
     log: Option<Log>,
     prior: Option<Prior>,
@@ -524,7 +524,7 @@ impl SettingsBuilder<AlgorithmSet> {
 
 // Parameters are set, move to defining error model
 impl SettingsBuilder<ParametersSet> {
-    pub fn set_error_models(self, ems: ErrorModels) -> SettingsBuilder<ErrorSet> {
+    pub fn set_error_models(self, ems: AssayErrorModels) -> SettingsBuilder<ErrorSet> {
         SettingsBuilder {
             config: self.config,
             parameters: self.parameters,
@@ -575,7 +575,7 @@ fn parse_output_folder(path: String) -> String {
 #[cfg(test)]
 
 mod tests {
-    use pharmsol::{ErrorModel, ErrorPoly};
+    use pharmsol::{AssayErrorModel, AssayErrorModels, ErrorPoly};
 
     use super::*;
     use crate::algorithms::Algorithm;
@@ -584,10 +584,10 @@ mod tests {
     fn test_builder() {
         let parameters = Parameters::new().add("Ke", 0.0, 5.0).add("V", 10.0, 200.0);
 
-        let ems = ErrorModels::new()
+        let ems = AssayErrorModels::new()
             .add(
                 0,
-                ErrorModel::Proportional {
+                AssayErrorModel::Proportional {
                     gamma: pharmsol::Factor::Variable(5.0),
                     poly: ErrorPoly::new(0.0, 0.1, 0.0, 0.0),
                 },
diff --git a/src/structs/psi.rs b/src/structs/psi.rs
index c63756cae..0fda454e1 100644
--- a/src/structs/psi.rs
+++ b/src/structs/psi.rs
@@ -5,6 +5,7 @@ use faer_ext::IntoFaer;
 use faer_ext::IntoNdarray;
 use ndarray::{Array2, ArrayView2};
 use pharmsol::prelude::simulator::psi;
+use pharmsol::AssayErrorModels;
 use pharmsol::Data;
 use pharmsol::Equation;
 use pharmsol::ErrorModels;
@@ -220,7 +221,7 @@ pub(crate) fn calculate_psi(
     equation: &impl Equation,
     subjects: &Data,
     theta: &Theta,
-    error_models: &ErrorModels,
+    error_models: &AssayErrorModels,
     progress: bool,
     cache: bool,
 ) -> Result<Psi> {
@@ -230,7 +231,6 @@ pub(crate) fn calculate_psi(
         &theta.matrix().clone().as_ref().into_ndarray().to_owned(),
         error_models,
         progress,
-        cache,
     )?;
 
     Ok(psi_ndarray.view().into())

From 2dc947e72900000cdfd39094c5af1b6fe36e4c5c Mon Sep 17 00:00:00 2001
From: Markus <66058642+mhovd@users.noreply.github.com>
Date: Tue, 17 Mar 2026 09:00:44 +0100
Subject: [PATCH 2/9] Replace faer_ext operations with faer

---
 src/algorithms/mod.rs |  9 ++++----
 src/structs/psi.rs    | 52 +++++++++++++++++--------------------------
 2 files changed, 24 insertions(+), 37 deletions(-)

diff --git a/src/algorithms/mod.rs b/src/algorithms/mod.rs
index 3636184ef..1db36a506 100644
--- a/src/algorithms/mod.rs
+++ b/src/algorithms/mod.rs
@@ -7,7 +7,6 @@ use crate::structs::psi::Psi;
 use crate::structs::theta::Theta;
 use anyhow::Context;
 use anyhow::Result;
-use faer_ext::IntoNdarray;
 use ndarray::parallel::prelude::{IntoParallelIterator, ParallelIterator};
 
 use npag::*;
@@ -37,11 +36,11 @@ pub trait Algorithms<E: Equation + Send + 'static>: Sync + Send + 'static {
         let mut nan_count = 0;
         let mut inf_count = 0;
 
-        let psi = self.psi().matrix().as_ref().into_ndarray();
+        let psi = self.psi().matrix();
         // 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();
+            for j in 0..psi.ncols() {
+                let val = psi[(i, j)];
                 if val.is_nan() {
                     nan_count += 1;
                     // *val = 0.0;
@@ -61,7 +60,7 @@ pub trait Algorithms<E: Equation + Send + 'static>: Sync + Send + 'static {
             );
         }
 
-        let (row, col) = psi.dim();
+        let (row, col) = (psi.nrows(), psi.ncols());
         let plam: Vec<f64> = (0..row)
             .map(|i| (0..col).map(|j| psi[(i, j)]).sum::<f64>())
             .collect();
diff --git a/src/structs/psi.rs b/src/structs/psi.rs
index 0fda454e1..ea6fb0131 100644
--- a/src/structs/psi.rs
+++ b/src/structs/psi.rs
@@ -1,14 +1,11 @@
 use anyhow::bail;
 use anyhow::Result;
 use faer::Mat;
-use faer_ext::IntoFaer;
-use faer_ext::IntoNdarray;
-use ndarray::{Array2, ArrayView2};
-use pharmsol::prelude::simulator::psi;
+use ndarray::Array2;
+use pharmsol::prelude::simulator::log_likelihood_matrix;
 use pharmsol::AssayErrorModels;
 use pharmsol::Data;
 use pharmsol::Equation;
-use pharmsol::ErrorModels;
 use serde::{Deserialize, Serialize};
 
 use super::theta::Theta;
@@ -114,7 +111,7 @@ impl Default for Psi {
 
 impl From<Array2<f64>> for Psi {
     fn from(array: Array2<f64>) -> Self {
-        let matrix = array.view().into_faer().to_owned();
+        let matrix = Mat::from_fn(array.nrows(), array.ncols(), |i, j| array[(i, j)]);
         Psi { matrix }
     }
 }
@@ -125,16 +122,9 @@ impl From<Mat<f64>> for Psi {
     }
 }
 
-impl From<ArrayView2<'_, f64>> for Psi {
-    fn from(array_view: ArrayView2<'_, f64>) -> Self {
-        let matrix = array_view.into_faer().to_owned();
-        Psi { matrix }
-    }
-}
-
 impl From<&Array2<f64>> for Psi {
     fn from(array: &Array2<f64>) -> Self {
-        let matrix = array.view().into_faer().to_owned();
+        let matrix = Mat::from_fn(array.nrows(), array.ncols(), |i, j| array[(i, j)]);
         Psi { matrix }
     }
 }
@@ -225,15 +215,13 @@ pub(crate) fn calculate_psi(
     progress: bool,
     cache: bool,
 ) -> Result<Psi> {
-    let psi_ndarray = psi(
-        equation,
-        subjects,
-        &theta.matrix().clone().as_ref().into_ndarray().to_owned(),
-        error_models,
-        progress,
-    )?;
-
-    Ok(psi_ndarray.view().into())
+    let tm = theta.matrix();
+    let theta_ndarray = Array2::from_shape_fn((tm.nrows(), tm.ncols()), |(i, j)| tm[(i, j)]);
+    let log_psi =
+        log_likelihood_matrix(equation, subjects, &theta_ndarray, error_models, progress)?;
+    let psi_ndarray = log_psi.mapv(f64::exp);
+
+    Ok(Psi::from(psi_ndarray))
 }
 
 #[cfg(test)]
@@ -252,11 +240,11 @@ mod tests {
         assert_eq!(psi.nspp(), 2);
         assert_eq!(psi.nsub(), 3);
 
-        // Check values by converting back to ndarray and comparing
-        let result_array = psi.matrix().as_ref().into_ndarray();
+        // Check values using faer matrix directly
+        let m = psi.matrix();
         for i in 0..2 {
             for j in 0..3 {
-                assert_eq!(result_array[[i, j]], array[[i, j]]);
+                assert_eq!(m[(i, j)], array[[i, j]]);
             }
         }
     }
@@ -273,11 +261,11 @@ mod tests {
         assert_eq!(psi.nspp(), 3);
         assert_eq!(psi.nsub(), 2);
 
-        // Check values by converting back to ndarray and comparing
-        let result_array = psi.matrix().as_ref().into_ndarray();
+        // Check values using faer matrix directly
+        let m = psi.matrix();
         for i in 0..3 {
             for j in 0..2 {
-                assert_eq!(result_array[[i, j]], array[[i, j]]);
+                assert_eq!(m[(i, j)], array[[i, j]]);
             }
         }
     }
@@ -350,12 +338,12 @@ mod tests {
         assert_eq!(psi_from_owned.nsub(), psi_from_ref.nsub());
 
         // And the same values
-        let owned_array = psi_from_owned.matrix().as_ref().into_ndarray();
-        let ref_array = psi_from_ref.matrix().as_ref().into_ndarray();
+        let owned_m = psi_from_owned.matrix();
+        let ref_m = psi_from_ref.matrix();
 
         for i in 0..2 {
             for j in 0..3 {
-                assert_eq!(owned_array[[i, j]], ref_array[[i, j]]);
+                assert_eq!(owned_m[(i, j)], ref_m[(i, j)]);
             }
         }
     }

From 92e5553a333d6d7b5d9ef062bd95856efff6e378 Mon Sep 17 00:00:00 2001
From: Markus <markushh@uio.no>
Date: Sun, 22 Mar 2026 21:32:43 +0100
Subject: [PATCH 3/9] Update error model

---
 Cargo.toml                         |  1 -
 benches/bimodal_ke.rs              |  6 ++--
 examples/bestdose.rs               |  4 +--
 examples/bestdose_auc.rs           |  4 +--
 examples/bestdose_bounds.rs        |  4 +--
 examples/bimodal_ke/main.rs        |  6 ++--
 examples/drusano/main.rs           | 12 +++----
 examples/iov/main.rs               |  4 +--
 examples/meta/main.rs              |  6 ++--
 examples/neely/main.rs             |  8 ++---
 examples/new_iov/main.rs           |  4 +--
 examples/theophylline/main.rs      |  4 +--
 examples/two_eq_lag/main.rs        |  4 +--
 examples/vanco_sde/main.rs         |  4 +--
 src/algorithms/npag.rs             |  3 --
 src/algorithms/npod.rs             |  3 --
 src/algorithms/postprob.rs         |  1 -
 src/bestdose/posterior.rs          |  2 +-
 src/routines/initialization/mod.rs | 14 ++++----
 src/routines/output/cycles.rs      | 26 +++++++-------
 src/structs/psi.rs                 |  1 -
 tests/bestdose_tests.rs            | 56 +++++++++++++++---------------
 tests/onecomp.rs                   | 12 +++----
 tests/settings_tests.rs            | 34 +++++++++---------
 24 files changed, 107 insertions(+), 116 deletions(-)

diff --git a/Cargo.toml b/Cargo.toml
index 5dfb2883d..dad6db1aa 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -28,7 +28,6 @@ tracing-subscriber = { version = "0.3.19", features = [
     "time",
 ] }
 faer = "0.23.1"
-faer-ext = { version = "0.7.1", features = ["nalgebra", "ndarray"] }
 pharmsol = "=0.23.0"
 rand = "0.9.0"
 anyhow = "1.0.100"
diff --git a/benches/bimodal_ke.rs b/benches/bimodal_ke.rs
index c2518e636..9655cdc28 100644
--- a/benches/bimodal_ke.rs
+++ b/benches/bimodal_ke.rs
@@ -27,10 +27,10 @@ fn create_parameters() -> Parameters {
         .add("v", 25.0, 250.0)
 }
 
-fn create_error_models() -> Result<ErrorModels> {
-    Ok(ErrorModels::new().add(
+fn create_error_models() -> Result<AssayErrorModels> {
+    Ok(AssayErrorModels::new().add(
         0,
-        ErrorModel::additive(ErrorPoly::new(0.0, 0.5, 0.0, 0.0), 0.0),
+        AssayErrorModel::additive(ErrorPoly::new(0.0, 0.5, 0.0, 0.0), 0.0),
     )?)
 }
 
diff --git a/examples/bestdose.rs b/examples/bestdose.rs
index 7761d8433..8272f1a2c 100644
--- a/examples/bestdose.rs
+++ b/examples/bestdose.rs
@@ -27,9 +27,9 @@ fn main() -> Result<()> {
         .add("ke", 0.001, 3.0)
         .add("v", 25.0, 250.0);
 
-    let ems = ErrorModels::new().add(
+    let ems = AssayErrorModels::new().add(
         0,
-        ErrorModel::additive(ErrorPoly::new(0.0, 0.20, 0.0, 0.0), 0.0),
+        AssayErrorModel::additive(ErrorPoly::new(0.0, 0.20, 0.0, 0.0), 0.0),
     )?;
 
     // Make settings
diff --git a/examples/bestdose_auc.rs b/examples/bestdose_auc.rs
index c1f44d2e2..0ffa444e3 100644
--- a/examples/bestdose_auc.rs
+++ b/examples/bestdose_auc.rs
@@ -30,9 +30,9 @@ fn main() -> Result<()> {
         .add("ke", 0.001, 3.0)
         .add("v", 25.0, 250.0);
 
-    let ems = ErrorModels::new().add(
+    let ems = AssayErrorModels::new().add(
         0,
-        ErrorModel::additive(ErrorPoly::new(0.0, 0.20, 0.0, 0.0), 0.0),
+        AssayErrorModel::additive(ErrorPoly::new(0.0, 0.20, 0.0, 0.0), 0.0),
     )?;
 
     let mut settings = Settings::builder()
diff --git a/examples/bestdose_bounds.rs b/examples/bestdose_bounds.rs
index f442278db..e256ef62b 100644
--- a/examples/bestdose_bounds.rs
+++ b/examples/bestdose_bounds.rs
@@ -29,9 +29,9 @@ fn main() -> Result<()> {
         .add("ke", 0.001, 3.0)
         .add("v", 25.0, 250.0);
 
-    let ems = ErrorModels::new().add(
+    let ems = AssayErrorModels::new().add(
         0,
-        ErrorModel::additive(ErrorPoly::new(0.0, 0.20, 0.0, 0.0), 0.0),
+        AssayErrorModel::additive(ErrorPoly::new(0.0, 0.20, 0.0, 0.0), 0.0),
     )?;
 
     let mut settings = Settings::builder()
diff --git a/examples/bimodal_ke/main.rs b/examples/bimodal_ke/main.rs
index bbd92febe..2049e531c 100644
--- a/examples/bimodal_ke/main.rs
+++ b/examples/bimodal_ke/main.rs
@@ -22,13 +22,13 @@ fn main() -> Result<()> {
         .add("ke", 0.001, 3.0)
         .add("v", 25.0, 250.0);
 
-    let ems = ErrorModels::new()
+    let ems = AssayErrorModels::new()
         .add(
             0,
-            ErrorModel::additive(ErrorPoly::new(0.0, 0.5, 0.0, 0.0), 0.0),
+            AssayErrorModel::additive(ErrorPoly::new(0.0, 0.5, 0.0, 0.0), 0.0),
         )
         .unwrap()
-        .add(1, ErrorModel::None)
+        .add(1, AssayErrorModel::None)
         .unwrap();
 
     let mut settings = Settings::builder()
diff --git a/examples/drusano/main.rs b/examples/drusano/main.rs
index 2262b249a..a9ff9f982 100644
--- a/examples/drusano/main.rs
+++ b/examples/drusano/main.rs
@@ -100,26 +100,26 @@ fn main() -> Result<()> {
         .add("h1r1", 5.0, 25.0)
         .add("h2r2", 10.0, 22.0);
 
-    let ems = ErrorModels::new()
+    let ems = AssayErrorModels::new()
         .add(
             0,
-            ErrorModel::proportional(ErrorPoly::new(0.1, 0.1, 0.0, 0.0), 1.0),
+            AssayErrorModel::proportional(ErrorPoly::new(0.1, 0.1, 0.0, 0.0), 1.0),
         )?
         .add(
             1,
-            ErrorModel::proportional(ErrorPoly::new(0.1, 0.1, 0.0, 0.0), 1.0),
+            AssayErrorModel::proportional(ErrorPoly::new(0.1, 0.1, 0.0, 0.0), 1.0),
         )?
         .add(
             2,
-            ErrorModel::proportional(ErrorPoly::new(0.1, 0.1, 0.0, 0.0), 1.0),
+            AssayErrorModel::proportional(ErrorPoly::new(0.1, 0.1, 0.0, 0.0), 1.0),
         )?
         .add(
             3,
-            ErrorModel::proportional(ErrorPoly::new(0.1, 0.1, 0.0, 0.0), 1.0),
+            AssayErrorModel::proportional(ErrorPoly::new(0.1, 0.1, 0.0, 0.0), 1.0),
         )?
         .add(
             4,
-            ErrorModel::proportional(ErrorPoly::new(0.1, 0.1, 0.0, 0.0), 1.0),
+            AssayErrorModel::proportional(ErrorPoly::new(0.1, 0.1, 0.0, 0.0), 1.0),
         )?;
 
     let mut settings = SettingsBuilder::new()
diff --git a/examples/iov/main.rs b/examples/iov/main.rs
index ffa3f39d2..82d095751 100644
--- a/examples/iov/main.rs
+++ b/examples/iov/main.rs
@@ -31,9 +31,9 @@ fn main() -> Result<()> {
 
     let params = Parameters::new().add("ke0", 0.001, 2.0);
 
-    let ems = ErrorModels::new().add(
+    let ems = AssayErrorModels::new().add(
         0,
-        ErrorModel::additive(ErrorPoly::new(0.0, 0.0, 0.0, 0.0), 0.0000757575757576),
+        AssayErrorModel::additive(ErrorPoly::new(0.0, 0.0, 0.0, 0.0), 0.0000757575757576),
     )?;
 
     let mut settings = Settings::builder()
diff --git a/examples/meta/main.rs b/examples/meta/main.rs
index f2784a7c3..817ffad20 100644
--- a/examples/meta/main.rs
+++ b/examples/meta/main.rs
@@ -41,15 +41,15 @@ fn main() {
         .add("theta2", 0.1, 10.0)
         .add("vs", 1.0, 10.0);
 
-    let ems = ErrorModels::new()
+    let ems = AssayErrorModels::new()
         .add(
             0,
-            ErrorModel::proportional(ErrorPoly::new(1.0, 0.1, 0.0, 0.0), 5.0),
+            AssayErrorModel::proportional(ErrorPoly::new(1.0, 0.1, 0.0, 0.0), 5.0),
         )
         .unwrap()
         .add(
             1,
-            ErrorModel::proportional(ErrorPoly::new(1.0, 0.1, 0.0, 0.0), 5.0),
+            AssayErrorModel::proportional(ErrorPoly::new(1.0, 0.1, 0.0, 0.0), 5.0),
         )
         .unwrap();
 
diff --git a/examples/neely/main.rs b/examples/neely/main.rs
index 47a8894a6..3bcb332a4 100644
--- a/examples/neely/main.rs
+++ b/examples/neely/main.rs
@@ -66,20 +66,20 @@ fn main() {
         .add("theta1", -4.0, 2.0)
         .add("theta2", -2.0, 0.5);
 
-    let ems = ErrorModels::new()
+    let ems = AssayErrorModels::new()
         .add(
             0,
-            ErrorModel::proportional(ErrorPoly::new(1.0, 0.1, 0.0, 0.0), 5.0),
+            AssayErrorModel::proportional(ErrorPoly::new(1.0, 0.1, 0.0, 0.0), 5.0),
         )
         .unwrap()
         .add(
             1,
-            ErrorModel::proportional(ErrorPoly::new(1.0, 0.1, 0.0, 0.0), 5.0),
+            AssayErrorModel::proportional(ErrorPoly::new(1.0, 0.1, 0.0, 0.0), 5.0),
         )
         .unwrap()
         .add(
             2,
-            ErrorModel::proportional(ErrorPoly::new(1.0, 0.1, 0.0, 0.0), 5.0),
+            AssayErrorModel::proportional(ErrorPoly::new(1.0, 0.1, 0.0, 0.0), 5.0),
         )
         .unwrap();
     let mut settings = Settings::builder()
diff --git a/examples/new_iov/main.rs b/examples/new_iov/main.rs
index e23ff50a6..f60d0d683 100644
--- a/examples/new_iov/main.rs
+++ b/examples/new_iov/main.rs
@@ -33,10 +33,10 @@ fn main() {
         .add("ke0", 0.0001, 2.4)
         .add("ske", 0.0001, 0.2);
 
-    let ems = ErrorModels::new()
+    let ems = AssayErrorModels::new()
         .add(
             0,
-            ErrorModel::additive(ErrorPoly::new(-0.00119, 0.44379, -0.45864, 0.16537), 0.0),
+            AssayErrorModel::additive(ErrorPoly::new(-0.00119, 0.44379, -0.45864, 0.16537), 0.0),
         )
         .unwrap();
 
diff --git a/examples/theophylline/main.rs b/examples/theophylline/main.rs
index 67e2b5724..774df5a69 100644
--- a/examples/theophylline/main.rs
+++ b/examples/theophylline/main.rs
@@ -35,10 +35,10 @@ fn main() {
         .add("ke", 0.001, 3.0)
         .add("v", 0.001, 50.0);
 
-    let ems = ErrorModels::new()
+    let ems = AssayErrorModels::new()
         .add(
             0,
-            ErrorModel::proportional(ErrorPoly::new(0.1, 0.1, 0.0, 0.0), 2.0),
+            AssayErrorModel::proportional(ErrorPoly::new(0.1, 0.1, 0.0, 0.0), 2.0),
         )
         .unwrap();
 
diff --git a/examples/two_eq_lag/main.rs b/examples/two_eq_lag/main.rs
index b9226003d..eebeef1c1 100644
--- a/examples/two_eq_lag/main.rs
+++ b/examples/two_eq_lag/main.rs
@@ -73,10 +73,10 @@ fn main() {
         .add("tlag", 0.0, 4.0)
         .add("v", 30.0, 120.0);
 
-    let ems = ErrorModels::new()
+    let ems = AssayErrorModels::new()
         .add(
             0,
-            ErrorModel::additive(ErrorPoly::new(-0.00119, 0.44379, -0.45864, 0.16537), 0.0),
+            AssayErrorModel::additive(ErrorPoly::new(-0.00119, 0.44379, -0.45864, 0.16537), 0.0),
         )
         .unwrap();
 
diff --git a/examples/vanco_sde/main.rs b/examples/vanco_sde/main.rs
index 94459d93e..29533be4e 100644
--- a/examples/vanco_sde/main.rs
+++ b/examples/vanco_sde/main.rs
@@ -56,10 +56,10 @@ fn main() {
         .add("vol", 0.2, 12.0)
         .add("ske", 0.0001, 0.2);
 
-    let ems = ErrorModels::new()
+    let ems = AssayErrorModels::new()
         .add(
             0,
-            ErrorModel::additive(ErrorPoly::new(0.00119, 0.20, 0.0, 0.0), 0.0),
+            AssayErrorModel::additive(ErrorPoly::new(0.00119, 0.20, 0.0, 0.0), 0.0),
         )
         .unwrap();
 
diff --git a/src/algorithms/npag.rs b/src/algorithms/npag.rs
index ac5a4b94e..3f17989e8 100644
--- a/src/algorithms/npag.rs
+++ b/src/algorithms/npag.rs
@@ -203,7 +203,6 @@ impl<E: Equation + Send + 'static> Algorithms<E> for NPAG<E> {
             &self.theta,
             &self.error_models,
             self.cycle == 1 && self.settings.config().progress,
-            self.cycle != 1,
         )?;
 
         if let Err(err) = self.validate_psi() {
@@ -315,7 +314,6 @@ impl<E: Equation + Send + 'static> Algorithms<E> for NPAG<E> {
                     &self.theta,
                     &error_model_up,
                     false,
-                    true,
                 )?;
                 let psi_down = calculate_psi(
                     &self.equation,
@@ -323,7 +321,6 @@ impl<E: Equation + Send + 'static> Algorithms<E> for NPAG<E> {
                     &self.theta,
                     &error_model_down,
                     false,
-                    true,
                 )?;
 
                 let (lambda_up, objf_up) = match burke(&psi_up) {
diff --git a/src/algorithms/npod.rs b/src/algorithms/npod.rs
index 9fde065df..a5756ffe8 100644
--- a/src/algorithms/npod.rs
+++ b/src/algorithms/npod.rs
@@ -212,7 +212,6 @@ impl<E: Equation + Send + 'static> Algorithms<E> for NPOD<E> {
             &self.theta,
             &error_model,
             self.cycle == 1 && self.settings.config().progress,
-            self.cycle != 1,
         )?;
 
         if let Err(err) = self.validate_psi() {
@@ -316,7 +315,6 @@ impl<E: Equation + Send + 'static> Algorithms<E> for NPOD<E> {
                     &self.theta,
                     &error_model_up,
                     false,
-                    true,
                 )?;
                 let psi_down = calculate_psi(
                     &self.equation,
@@ -324,7 +322,6 @@ impl<E: Equation + Send + 'static> Algorithms<E> for NPOD<E> {
                     &self.theta,
                     &error_model_down,
                     false,
-                    true,
                 )?;
 
                 let (lambda_up, objf_up) = match burke(&psi_up) {
diff --git a/src/algorithms/postprob.rs b/src/algorithms/postprob.rs
index de2fa0307..19d030f66 100644
--- a/src/algorithms/postprob.rs
+++ b/src/algorithms/postprob.rs
@@ -125,7 +125,6 @@ impl<E: Equation + Send + 'static> Algorithms<E> for POSTPROB<E> {
             &self.theta,
             &self.error_models,
             false,
-            false,
         )?;
         (self.w, self.objf) = burke(&self.psi).context("Error in IPM")?;
         Ok(())
diff --git a/src/bestdose/posterior.rs b/src/bestdose/posterior.rs
index e86a4ae07..4f73ea6de 100644
--- a/src/bestdose/posterior.rs
+++ b/src/bestdose/posterior.rs
@@ -99,7 +99,7 @@ pub fn npagfull11_filter(
     tracing::info!("Stage 1.1: NPAGFULL11 Bayesian filtering");
 
     // Calculate psi matrix P(data|theta_i) for all support points
-    let psi = calculate_psi(eq, past_data, population_theta, error_models, false, true)?;
+    let psi = calculate_psi(eq, past_data, population_theta, error_models, false)?;
 
     // First burke call to get initial posterior probabilities
     let (initial_weights, _) = burke(&psi)?;
diff --git a/src/routines/initialization/mod.rs b/src/routines/initialization/mod.rs
index cff7e6046..fa03c4442 100644
--- a/src/routines/initialization/mod.rs
+++ b/src/routines/initialization/mod.rs
@@ -206,14 +206,14 @@ pub fn parse_prior(path: &String, settings: &Settings) -> Result<(Theta, Option<
 mod tests {
     use super::*;
     use crate::prelude::*;
-    use pharmsol::{ErrorModel, ErrorModels, ErrorPoly};
+    use pharmsol::{AssayErrorModel, AssayErrorModels, ErrorPoly};
     use std::fs;
 
     fn create_test_settings() -> Settings {
         let parameters = Parameters::new().add("ke", 0.1, 1.0).add("v", 5.0, 50.0);
 
-        let em = ErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0);
-        let ems = ErrorModels::new().add(0, em).unwrap();
+        let em = AssayErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0);
+        let ems = AssayErrorModels::new().add(0, em).unwrap();
 
         Settings::builder()
             .set_algorithm(Algorithm::NPAG)
@@ -319,8 +319,8 @@ mod tests {
             .add("ke", f64::NEG_INFINITY, 1.0) // Invalid: infinite lower bound
             .add("v", 5.0, 50.0);
 
-        let em = ErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0);
-        let ems = ErrorModels::new().add(0, em).unwrap();
+        let em = AssayErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0);
+        let ems = AssayErrorModels::new().add(0, em).unwrap();
 
         let mut settings = Settings::builder()
             .set_algorithm(Algorithm::NPAG)
@@ -341,8 +341,8 @@ mod tests {
             .add("ke", 1.0, 0.5) // Invalid: lower bound >= upper bound
             .add("v", 5.0, 50.0);
 
-        let em = ErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0);
-        let ems = ErrorModels::new().add(0, em).unwrap();
+        let em = AssayErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0);
+        let ems = AssayErrorModels::new().add(0, em).unwrap();
 
         let mut settings = Settings::builder()
             .set_algorithm(Algorithm::NPAG)
diff --git a/src/routines/output/cycles.rs b/src/routines/output/cycles.rs
index f65720aa9..24af587e0 100644
--- a/src/routines/output/cycles.rs
+++ b/src/routines/output/cycles.rs
@@ -1,6 +1,6 @@
 use anyhow::Result;
 use csv::WriterBuilder;
-use pharmsol::{ErrorModel, ErrorModels};
+use pharmsol::{AssayErrorModel, AssayErrorModels};
 use serde::Serialize;
 
 use crate::{
@@ -23,7 +23,7 @@ use crate::{
 pub struct NPCycle {
     cycle: usize,
     objf: f64,
-    error_models: ErrorModels,
+    error_models: AssayErrorModels,
     theta: Theta,
     nspp: usize,
     delta_objf: f64,
@@ -34,7 +34,7 @@ impl NPCycle {
     pub fn new(
         cycle: usize,
         objf: f64,
-        error_models: ErrorModels,
+        error_models: AssayErrorModels,
         theta: Theta,
         nspp: usize,
         delta_objf: f64,
@@ -57,7 +57,7 @@ impl NPCycle {
     pub fn objf(&self) -> f64 {
         self.objf
     }
-    pub fn error_models(&self) -> &ErrorModels {
+    pub fn error_models(&self) -> &AssayErrorModels {
         &self.error_models
     }
     pub fn theta(&self) -> &Theta {
@@ -77,7 +77,7 @@ impl NPCycle {
         Self {
             cycle: 0,
             objf: 0.0,
-            error_models: ErrorModels::default(),
+            error_models: AssayErrorModels::default(),
             theta: Theta::new(),
             nspp: 0,
             delta_objf: 0.0,
@@ -120,15 +120,15 @@ impl CycleLog {
         writer.write_field("nspp")?;
         if let Some(first_cycle) = self.cycles.first() {
             first_cycle.error_models.iter().try_for_each(
-                |(outeq, errmod): (usize, &ErrorModel)| -> Result<(), csv::Error> {
+                |(outeq, errmod): (usize, &AssayErrorModel)| -> Result<(), csv::Error> {
                     match errmod {
-                        ErrorModel::Additive { .. } => {
+                        AssayErrorModel::Additive { .. } => {
                             writer.write_field(format!("gamlam.{}", outeq))?;
                         }
-                        ErrorModel::Proportional { .. } => {
+                        AssayErrorModel::Proportional { .. } => {
                             writer.write_field(format!("gamlam.{}", outeq))?;
                         }
-                        ErrorModel::None => {}
+                        AssayErrorModel::None => {}
                     }
                     Ok(())
                 },
@@ -158,15 +158,15 @@ impl CycleLog {
 
             // Write the error models
             cycle.error_models.iter().try_for_each(
-                |(_, errmod): (usize, &ErrorModel)| -> Result<()> {
+                |(_, errmod): (usize, &AssayErrorModel)| -> Result<()> {
                     match errmod {
-                        ErrorModel::Additive { lambda: _, poly: _ } => {
+                        AssayErrorModel::Additive { lambda: _, poly: _ } => {
                             writer.write_field(format!("{:.5}", errmod.factor()?))?;
                         }
-                        ErrorModel::Proportional { gamma: _, poly: _ } => {
+                        AssayErrorModel::Proportional { gamma: _, poly: _ } => {
                             writer.write_field(format!("{:.5}", errmod.factor()?))?;
                         }
-                        ErrorModel::None => {}
+                        AssayErrorModel::None => {}
                     }
                     Ok(())
                 },
diff --git a/src/structs/psi.rs b/src/structs/psi.rs
index ea6fb0131..e50dfd0c7 100644
--- a/src/structs/psi.rs
+++ b/src/structs/psi.rs
@@ -213,7 +213,6 @@ pub(crate) fn calculate_psi(
     theta: &Theta,
     error_models: &AssayErrorModels,
     progress: bool,
-    cache: bool,
 ) -> Result<Psi> {
     let tm = theta.matrix();
     let theta_ndarray = Array2::from_shape_fn((tm.nrows(), tm.ncols()), |(i, j)| tm[(i, j)]);
diff --git a/tests/bestdose_tests.rs b/tests/bestdose_tests.rs
index 53ac2618e..6033a6d93 100644
--- a/tests/bestdose_tests.rs
+++ b/tests/bestdose_tests.rs
@@ -26,9 +26,9 @@ fn test_infusion_mask_inclusion() -> Result<()> {
 
     let params = Parameters::new().add("ke", 0.1, 0.5).add("v", 40.0, 60.0);
 
-    let ems = ErrorModels::new().add(
+    let ems = AssayErrorModels::new().add(
         0,
-        ErrorModel::additive(ErrorPoly::new(0.0, 0.20, 0.0, 0.0), 0.0),
+        AssayErrorModel::additive(ErrorPoly::new(0.0, 0.20, 0.0, 0.0), 0.0),
     )?;
 
     let mut settings = Settings::builder()
@@ -143,9 +143,9 @@ fn test_fixed_infusion_preservation() -> Result<()> {
         .add("ke", 0.001, 3.0)
         .add("v", 25.0, 250.0);
 
-    let ems = ErrorModels::new().add(
+    let ems = AssayErrorModels::new().add(
         0,
-        ErrorModel::additive(ErrorPoly::new(0.0, 0.20, 0.0, 0.0), 0.0),
+        AssayErrorModel::additive(ErrorPoly::new(0.0, 0.20, 0.0, 0.0), 0.0),
     )?;
 
     let mut settings = Settings::builder()
@@ -230,9 +230,9 @@ fn test_dose_count_validation() -> Result<()> {
     );
 
     let params = Parameters::new().add("ke", 0.1, 0.5).add("v", 40.0, 60.0);
-    let ems = ErrorModels::new().add(
+    let ems = AssayErrorModels::new().add(
         0,
-        ErrorModel::additive(ErrorPoly::new(0.0, 0.20, 0.0, 0.0), 0.0),
+        AssayErrorModel::additive(ErrorPoly::new(0.0, 0.20, 0.0, 0.0), 0.0),
     )?;
 
     let mut settings = Settings::builder()
@@ -310,9 +310,9 @@ fn test_empty_observations_validation() -> Result<()> {
     );
 
     let params = Parameters::new().add("ke", 0.1, 0.5).add("v", 40.0, 60.0);
-    let ems = ErrorModels::new().add(
+    let ems = AssayErrorModels::new().add(
         0,
-        ErrorModel::additive(ErrorPoly::new(0.0, 0.20, 0.0, 0.0), 0.0),
+        AssayErrorModel::additive(ErrorPoly::new(0.0, 0.20, 0.0, 0.0), 0.0),
     )?;
 
     let mut settings = Settings::builder()
@@ -377,9 +377,9 @@ fn test_basic_auc_mode() -> Result<()> {
 
     let params = Parameters::new().add("ke", 0.1, 0.5).add("v", 40.0, 60.0);
 
-    let ems = ErrorModels::new().add(
+    let ems = AssayErrorModels::new().add(
         0,
-        ErrorModel::additive(ErrorPoly::new(0.0, 0.20, 0.0, 0.0), 0.0),
+        AssayErrorModel::additive(ErrorPoly::new(0.0, 0.20, 0.0, 0.0), 0.0),
     )?;
 
     let mut settings = Settings::builder()
@@ -468,9 +468,9 @@ fn test_infusion_auc_mode() -> Result<()> {
 
     let params = Parameters::new().add("ke", 0.1, 0.5).add("v", 40.0, 60.0);
 
-    let ems = ErrorModels::new().add(
+    let ems = AssayErrorModels::new().add(
         0,
-        ErrorModel::additive(ErrorPoly::new(0.0, 0.20, 0.0, 0.0), 0.0),
+        AssayErrorModel::additive(ErrorPoly::new(0.0, 0.20, 0.0, 0.0), 0.0),
     )?;
 
     let mut settings = Settings::builder()
@@ -578,8 +578,8 @@ fn test_multi_outeq_auc_mode() -> Result<()> {
 
     let params = Parameters::new().add("ke", 0.1, 0.5).add("v", 40.0, 60.0);
 
-    let error_model = ErrorModel::additive(ErrorPoly::new(0.0, 5.0, 0.0, 0.0), 0.0);
-    let ems = ErrorModels::new()
+    let error_model = AssayErrorModel::additive(ErrorPoly::new(0.0, 5.0, 0.0, 0.0), 0.0);
+    let ems = AssayErrorModels::new()
         .add(0, error_model.clone())?
         .add(1, error_model)?;
 
@@ -651,8 +651,8 @@ fn test_multi_outeq_auc_optimization() -> Result<()> {
     );
 
     let params = Parameters::new().add("ke", 0.1, 0.5).add("v", 40.0, 60.0);
-    let error_model = ErrorModel::additive(ErrorPoly::new(0.0, 5.0, 0.0, 0.0), 0.0);
-    let ems = ErrorModels::new()
+    let error_model = AssayErrorModel::additive(ErrorPoly::new(0.0, 5.0, 0.0, 0.0), 0.0);
+    let ems = AssayErrorModels::new()
         .add(0, error_model.clone())?
         .add(1, error_model)?;
 
@@ -743,9 +743,9 @@ fn test_auc_from_zero_single_dose() -> Result<()> {
 
     let params = Parameters::new().add("ke", 0.2, 0.4).add("v", 40.0, 60.0);
 
-    let ems = ErrorModels::new().add(
+    let ems = AssayErrorModels::new().add(
         0,
-        ErrorModel::additive(ErrorPoly::new(0.0, 0.20, 0.0, 0.0), 0.0),
+        AssayErrorModel::additive(ErrorPoly::new(0.0, 0.20, 0.0, 0.0), 0.0),
     )?;
 
     let mut settings = Settings::builder()
@@ -833,9 +833,9 @@ fn test_auc_from_last_dose_maintenance() -> Result<()> {
 
     let params = Parameters::new().add("ke", 0.2, 0.4).add("v", 40.0, 60.0);
 
-    let ems = ErrorModels::new().add(
+    let ems = AssayErrorModels::new().add(
         0,
-        ErrorModel::additive(ErrorPoly::new(0.0, 0.20, 0.0, 0.0), 0.0),
+        AssayErrorModel::additive(ErrorPoly::new(0.0, 0.20, 0.0, 0.0), 0.0),
     )?;
 
     let mut settings = Settings::builder()
@@ -927,9 +927,9 @@ fn test_auc_modes_comparison() -> Result<()> {
 
     let params = Parameters::new().add("ke", 0.3, 0.3).add("v", 50.0, 50.0);
 
-    let ems = ErrorModels::new().add(
+    let ems = AssayErrorModels::new().add(
         0,
-        ErrorModel::additive(ErrorPoly::new(0.0, 0.20, 0.0, 0.0), 0.0),
+        AssayErrorModel::additive(ErrorPoly::new(0.0, 0.20, 0.0, 0.0), 0.0),
     )?;
 
     let mut settings = Settings::builder()
@@ -1060,9 +1060,9 @@ fn test_auc_from_last_dose_multiple_observations() -> Result<()> {
 
     let params = Parameters::new().add("ke", 0.2, 0.4).add("v", 40.0, 60.0);
 
-    let ems = ErrorModels::new().add(
+    let ems = AssayErrorModels::new().add(
         0,
-        ErrorModel::additive(ErrorPoly::new(0.0, 0.20, 0.0, 0.0), 0.0),
+        AssayErrorModel::additive(ErrorPoly::new(0.0, 0.20, 0.0, 0.0), 0.0),
     )?;
 
     let mut settings = Settings::builder()
@@ -1160,9 +1160,9 @@ fn test_auc_from_last_dose_no_prior_dose() -> Result<()> {
 
     let params = Parameters::new().add("ke", 0.2, 0.4).add("v", 40.0, 60.0);
 
-    let ems = ErrorModels::new().add(
+    let ems = AssayErrorModels::new().add(
         0,
-        ErrorModel::additive(ErrorPoly::new(0.0, 0.20, 0.0, 0.0), 0.0),
+        AssayErrorModel::additive(ErrorPoly::new(0.0, 0.20, 0.0, 0.0), 0.0),
     )?;
 
     let mut settings = Settings::builder()
@@ -1255,9 +1255,9 @@ fn test_dose_range_bounds_respected() -> Result<()> {
 
     let params = Parameters::new().add("ke", 0.1, 0.5).add("v", 40.0, 60.0);
 
-    let ems = ErrorModels::new().add(
+    let ems = AssayErrorModels::new().add(
         0,
-        ErrorModel::additive(ErrorPoly::new(0.0, 0.20, 0.0, 0.0), 0.0),
+        AssayErrorModel::additive(ErrorPoly::new(0.0, 0.20, 0.0, 0.0), 0.0),
     )?;
 
     let mut settings = Settings::builder()
diff --git a/tests/onecomp.rs b/tests/onecomp.rs
index c9cbd2020..7ab278dd6 100644
--- a/tests/onecomp.rs
+++ b/tests/onecomp.rs
@@ -22,8 +22,8 @@ fn test_one_compartment_npag() -> Result<()> {
     // Define parameters
     let params = Parameters::new().add("ke", 0.1, 1.0).add("v", 1.0, 20.0);
 
-    let em = ErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0);
-    let ems = ErrorModels::new().add(0, em).unwrap();
+    let em = AssayErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0);
+    let ems = AssayErrorModels::new().add(0, em).unwrap();
 
     // Create settings
     let mut settings = Settings::builder()
@@ -90,8 +90,8 @@ fn test_one_compartment_npod() -> Result<()> {
     // Define parameters
     let params = Parameters::new().add("ke", 0.1, 1.0).add("v", 1.0, 20.0);
 
-    let em = ErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0);
-    let ems = ErrorModels::new().add(0, em).unwrap();
+    let em = AssayErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0);
+    let ems = AssayErrorModels::new().add(0, em).unwrap();
 
     // Create settings
     let mut settings = Settings::builder()
@@ -158,8 +158,8 @@ fn test_one_compartment_postprob() -> Result<()> {
     // Define parameters
     let params = Parameters::new().add("ke", 0.1, 1.0).add("v", 1.0, 20.0);
 
-    let em = ErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0);
-    let ems = ErrorModels::new().add(0, em).unwrap();
+    let em = AssayErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0);
+    let ems = AssayErrorModels::new().add(0, em).unwrap();
 
     // Create settings
     let mut settings = Settings::builder()
diff --git a/tests/settings_tests.rs b/tests/settings_tests.rs
index 188ebd832..23992c34b 100644
--- a/tests/settings_tests.rs
+++ b/tests/settings_tests.rs
@@ -6,8 +6,8 @@ use pmcore::prelude::*;
 fn test_settings_builder_basic() -> Result<()> {
     let params = Parameters::new().add("ke", 0.1, 1.0).add("v", 1.0, 20.0);
 
-    let em = ErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0);
-    let ems = ErrorModels::new().add(0, em)?;
+    let em = AssayErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0);
+    let ems = AssayErrorModels::new().add(0, em)?;
 
     let settings = Settings::builder()
         .set_algorithm(Algorithm::NPAG)
@@ -27,8 +27,8 @@ fn test_settings_builder_basic() -> Result<()> {
 fn test_settings_serialization() -> Result<()> {
     let params = Parameters::new().add("ke", 0.1, 1.0).add("v", 5.0, 15.0);
 
-    let em = ErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0);
-    let ems = ErrorModels::new().add(0, em)?;
+    let em = AssayErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0);
+    let ems = AssayErrorModels::new().add(0, em)?;
 
     let settings = Settings::builder()
         .set_algorithm(Algorithm::NPOD)
@@ -54,8 +54,8 @@ fn test_settings_serialization() -> Result<()> {
 #[test]
 fn test_settings_algorithms() -> Result<()> {
     let params = Parameters::new().add("ke", 0.1, 1.0);
-    let em = ErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0);
-    let ems = ErrorModels::new().add(0, em)?;
+    let em = AssayErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0);
+    let ems = AssayErrorModels::new().add(0, em)?;
 
     // Test NPAG
     let settings_npag = Settings::builder()
@@ -80,8 +80,8 @@ fn test_settings_algorithms() -> Result<()> {
 #[test]
 fn test_settings_setters() -> Result<()> {
     let params = Parameters::new().add("ke", 0.1, 1.0);
-    let em = ErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0);
-    let ems = ErrorModels::new().add(0, em)?;
+    let em = AssayErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0);
+    let ems = AssayErrorModels::new().add(0, em)?;
 
     let mut settings = Settings::builder()
         .set_algorithm(Algorithm::NPAG)
@@ -116,8 +116,8 @@ fn test_settings_setters() -> Result<()> {
 #[test]
 fn test_settings_with_prior() -> Result<()> {
     let params = Parameters::new().add("ke", 0.1, 1.0);
-    let em = ErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0);
-    let ems = ErrorModels::new().add(0, em)?;
+    let em = AssayErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0);
+    let ems = AssayErrorModels::new().add(0, em)?;
 
     let mut settings = Settings::builder()
         .set_algorithm(Algorithm::NPAG)
@@ -139,8 +139,8 @@ fn test_settings_with_prior() -> Result<()> {
 #[test]
 fn test_settings_latin_prior() -> Result<()> {
     let params = Parameters::new().add("ke", 0.1, 1.0).add("v", 5.0, 15.0);
-    let em = ErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0);
-    let ems = ErrorModels::new().add(0, em)?;
+    let em = AssayErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0);
+    let ems = AssayErrorModels::new().add(0, em)?;
 
     let mut settings = Settings::builder()
         .set_algorithm(Algorithm::NPAG)
@@ -181,10 +181,10 @@ fn test_parameters() {
 /// Test ErrorModels construction
 #[test]
 fn test_error_models() -> Result<()> {
-    let em1 = ErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0);
-    let em2 = ErrorModel::proportional(ErrorPoly::new(0.0, 0.0, 0.15, 0.0), 2.0);
+    let em1 = AssayErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0);
+    let em2 = AssayErrorModel::proportional(ErrorPoly::new(0.0, 0.0, 0.15, 0.0), 2.0);
 
-    let mut ems = ErrorModels::new();
+    let mut ems = AssayErrorModels::new();
     ems = ems.add(0, em1)?;
     ems = ems.add(1, em2)?;
 
@@ -198,8 +198,8 @@ fn test_error_models() -> Result<()> {
 #[test]
 fn test_config_accessors() -> Result<()> {
     let params = Parameters::new().add("ke", 0.1, 1.0);
-    let em = ErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0);
-    let ems = ErrorModels::new().add(0, em)?;
+    let em = AssayErrorModel::additive(ErrorPoly::new(0.0, 0.10, 0.0, 0.0), 2.0);
+    let ems = AssayErrorModels::new().add(0, em)?;
 
     let settings = Settings::builder()
         .set_algorithm(Algorithm::NPAG)

From 2cf3d31baef1e0e62e6b68c90535fd209826c460 Mon Sep 17 00:00:00 2001
From: Markus <markushh@uio.no>
Date: Sun, 22 Mar 2026 21:40:35 +0100
Subject: [PATCH 4/9] Clippy

---
 src/algorithms/mod.rs     | 3 ++-
 src/bestdose/posterior.rs | 1 +
 src/structs/weights.rs    | 5 +++++
 3 files changed, 8 insertions(+), 1 deletion(-)

diff --git a/src/algorithms/mod.rs b/src/algorithms/mod.rs
index 1db36a506..949ded9c0 100644
--- a/src/algorithms/mod.rs
+++ b/src/algorithms/mod.rs
@@ -298,13 +298,14 @@ pub trait Algorithms<E: Equation + Send + 'static>: Sync + Send + 'static {
     /// until the algorithm converges or meets a stopping criteria.
     fn fit(&mut self) -> Result<NPResult<E>> {
         self.initialize().unwrap();
+        #[allow(clippy::while_let_loop)]
         loop {
             match self.next_cycle()? {
                 Status::Continue => continue,
                 Status::Stop(_) => break,
             }
         }
-        Ok(self.into_npresult()?)
+        self.into_npresult()
     }
 
     #[allow(clippy::wrong_self_convention)]
diff --git a/src/bestdose/posterior.rs b/src/bestdose/posterior.rs
index 4f73ea6de..9506d390e 100644
--- a/src/bestdose/posterior.rs
+++ b/src/bestdose/posterior.rs
@@ -213,6 +213,7 @@ pub fn npagfull_refinement(
 
         // Run NPAG optimization
         let refinement_result = npag.initialize().and_then(|_| {
+            #[allow(clippy::while_let_loop)]
             loop {
                 match npag.next_cycle()? {
                     Status::Continue => continue,
diff --git a/src/structs/weights.rs b/src/structs/weights.rs
index e84974443..483817d48 100644
--- a/src/structs/weights.rs
+++ b/src/structs/weights.rs
@@ -58,6 +58,11 @@ impl Weights {
         self.weights.nrows()
     }
 
+    // Check if there are no weights.
+    pub fn is_empty(&self) -> bool {
+        self.len() == 0
+    }
+
     /// Get a vector representation of the weights.
     pub fn to_vec(&self) -> Vec<f64> {
         self.weights.iter().cloned().collect()

From de58a5f9ab38a6b0f3a090620765b309e6f90941 Mon Sep 17 00:00:00 2001
From: Markus <66058642+mhovd@users.noreply.github.com>
Date: Mon, 23 Mar 2026 16:14:50 +0100
Subject: [PATCH 5/9] Silence diffsol logs

---
 src/routines/logger.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/routines/logger.rs b/src/routines/logger.rs
index 2b5d411cc..d9ba24e6a 100644
--- a/src/routines/logger.rs
+++ b/src/routines/logger.rs
@@ -27,7 +27,7 @@ pub(crate) fn setup_log(settings: &mut Settings) -> Result<()> {
 
     // Use the log level defined in configuration file
     let log_level = settings.log().level.clone();
-    let env_filter = EnvFilter::new(log_level);
+    let env_filter = EnvFilter::new(format!("{},diffsol=off", log_level));
 
     let timestamper = CompactTimestamp {
         start: Instant::now(),

From 22a7cbf056fdc07e7c86697049e42a9d3013271c Mon Sep 17 00:00:00 2001
From: Markus <66058642+mhovd@users.noreply.github.com>
Date: Mon, 23 Mar 2026 16:14:55 +0100
Subject: [PATCH 6/9] Update BKE example

---
 examples/bimodal_ke/main.rs | 10 ++++------
 1 file changed, 4 insertions(+), 6 deletions(-)

diff --git a/examples/bimodal_ke/main.rs b/examples/bimodal_ke/main.rs
index 2049e531c..fe72dfb96 100644
--- a/examples/bimodal_ke/main.rs
+++ b/examples/bimodal_ke/main.rs
@@ -6,16 +6,16 @@ fn main() -> Result<()> {
         |x, p, _t, dx, b, rateiv, _cov| {
             // fetch_cov!(cov, t, wt);
             fetch_params!(p, ke, _v);
-            dx[0] = -ke * x[0] + rateiv[0] + b[0];
+            dx[0] = -ke * x[0] + rateiv[1] + b[1];
         },
         |_p, _t, _cov| lag! {},
         |_p, _t, _cov| fa! {},
         |_p, _t, _cov, _x| {},
         |x, p, _t, _cov, y| {
             fetch_params!(p, _ke, v);
-            y[0] = x[0] / v;
+            y[1] = x[0] / v;
         },
-        (1, 1),
+        (2, 2),
     );
 
     let params = Parameters::new()
@@ -24,11 +24,9 @@ fn main() -> Result<()> {
 
     let ems = AssayErrorModels::new()
         .add(
-            0,
+            1,
             AssayErrorModel::additive(ErrorPoly::new(0.0, 0.5, 0.0, 0.0), 0.0),
         )
-        .unwrap()
-        .add(1, AssayErrorModel::None)
         .unwrap();
 
     let mut settings = Settings::builder()

From af85309e78274129813089284c713d9297f95e19 Mon Sep 17 00:00:00 2001
From: Markus <66058642+mhovd@users.noreply.github.com>
Date: Mon, 23 Mar 2026 16:15:36 +0100
Subject: [PATCH 7/9] Update benchmark

---
 benches/bimodal_ke.rs | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/benches/bimodal_ke.rs b/benches/bimodal_ke.rs
index 9655cdc28..f1037b61b 100644
--- a/benches/bimodal_ke.rs
+++ b/benches/bimodal_ke.rs
@@ -8,16 +8,16 @@ fn create_equation() -> equation::ODE {
     equation::ODE::new(
         |x, p, _t, dx, b, rateiv, _cov| {
             fetch_params!(p, ke, _v);
-            dx[0] = -ke * x[0] + rateiv[0] + b[0];
+            dx[0] = -ke * x[0] + rateiv[1] + b[1];
         },
         |_p, _t, _cov| lag! {},
         |_p, _t, _cov| fa! {},
         |_p, _t, _cov, _x| {},
         |x, p, _t, _cov, y| {
             fetch_params!(p, _ke, v);
-            y[0] = x[0] / v;
+            y[1] = x[0] / v;
         },
-        (1, 1),
+        (2, 2),
     )
 }
 
@@ -29,7 +29,7 @@ fn create_parameters() -> Parameters {
 
 fn create_error_models() -> Result<AssayErrorModels> {
     Ok(AssayErrorModels::new().add(
-        0,
+        1,
         AssayErrorModel::additive(ErrorPoly::new(0.0, 0.5, 0.0, 0.0), 0.0),
     )?)
 }

From 0ce22393935badfa4e2400325e2db75c1a9d3b19 Mon Sep 17 00:00:00 2001
From: Markus <66058642+mhovd@users.noreply.github.com>
Date: Tue, 24 Mar 2026 12:13:47 +0100
Subject: [PATCH 8/9] Update to pharmsol 0.24

---
 Cargo.toml                    |  2 +-
 benches/bimodal_ke.rs         |  1 -
 examples/bestdose.rs          |  1 -
 examples/bestdose_auc.rs      |  1 -
 examples/bestdose_bounds.rs   |  1 -
 examples/bimodal_ke/main.rs   |  1 -
 examples/drusano/main.rs      |  1 -
 examples/iov/main.rs          |  1 -
 examples/meta/main.rs         |  1 -
 examples/neely/main.rs        |  1 -
 examples/new_iov/main.rs      |  1 -
 examples/theophylline/main.rs | 21 ++-------------------
 examples/two_eq_lag/main.rs   |  1 -
 examples/vanco.rs             |  1 -
 examples/vanco_sde/main.rs    |  1 -
 tests/bestdose_tests.rs       | 14 --------------
 tests/onecomp.rs              |  3 ---
 17 files changed, 3 insertions(+), 50 deletions(-)

diff --git a/Cargo.toml b/Cargo.toml
index dad6db1aa..df35f74df 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -28,7 +28,7 @@ tracing-subscriber = { version = "0.3.19", features = [
     "time",
 ] }
 faer = "0.23.1"
-pharmsol = "=0.23.0"
+pharmsol = "=0.24.0"
 rand = "0.9.0"
 anyhow = "1.0.100"
 rayon = "1.10.0"
diff --git a/benches/bimodal_ke.rs b/benches/bimodal_ke.rs
index f1037b61b..b8ce6c55e 100644
--- a/benches/bimodal_ke.rs
+++ b/benches/bimodal_ke.rs
@@ -17,7 +17,6 @@ fn create_equation() -> equation::ODE {
             fetch_params!(p, _ke, v);
             y[1] = x[0] / v;
         },
-        (2, 2),
     )
 }
 
diff --git a/examples/bestdose.rs b/examples/bestdose.rs
index 8272f1a2c..3114ea3bc 100644
--- a/examples/bestdose.rs
+++ b/examples/bestdose.rs
@@ -20,7 +20,6 @@ fn main() -> Result<()> {
             fetch_params!(p, _ke, v);
             y[0] = x[0] / v;
         },
-        (1, 1),
     );
 
     let params = Parameters::new()
diff --git a/examples/bestdose_auc.rs b/examples/bestdose_auc.rs
index 0ffa444e3..75ac3ee12 100644
--- a/examples/bestdose_auc.rs
+++ b/examples/bestdose_auc.rs
@@ -22,7 +22,6 @@ fn main() -> Result<()> {
             fetch_params!(p, _ke, v);
             y[0] = x[0] / v;
         },
-        (1, 1),
     );
 
     // Minimal parameter ranges
diff --git a/examples/bestdose_bounds.rs b/examples/bestdose_bounds.rs
index e256ef62b..6058b1cf3 100644
--- a/examples/bestdose_bounds.rs
+++ b/examples/bestdose_bounds.rs
@@ -22,7 +22,6 @@ fn main() -> Result<()> {
             fetch_params!(p, _ke, v);
             y[0] = x[0] / v;
         },
-        (1, 1),
     );
 
     let params = Parameters::new()
diff --git a/examples/bimodal_ke/main.rs b/examples/bimodal_ke/main.rs
index fe72dfb96..dfeac2b22 100644
--- a/examples/bimodal_ke/main.rs
+++ b/examples/bimodal_ke/main.rs
@@ -15,7 +15,6 @@ fn main() -> Result<()> {
             fetch_params!(p, _ke, v);
             y[1] = x[0] / v;
         },
-        (2, 2),
     );
 
     let params = Parameters::new()
diff --git a/examples/drusano/main.rs b/examples/drusano/main.rs
index a9ff9f982..0e030302d 100644
--- a/examples/drusano/main.rs
+++ b/examples/drusano/main.rs
@@ -71,7 +71,6 @@ fn main() -> Result<()> {
             y[3] = x[3].log10();
             y[4] = x[4].log10();
         },
-        (5, 5),
     );
 
     let params = Parameters::new()
diff --git a/examples/iov/main.rs b/examples/iov/main.rs
index 82d095751..023762ab9 100644
--- a/examples/iov/main.rs
+++ b/examples/iov/main.rs
@@ -25,7 +25,6 @@ fn main() -> Result<()> {
             fetch_params!(p, _ke0);
             y[0] = x[0] / 50.0;
         },
-        (2, 1),
         10000,
     );
 
diff --git a/examples/meta/main.rs b/examples/meta/main.rs
index 817ffad20..20fbb51c9 100644
--- a/examples/meta/main.rs
+++ b/examples/meta/main.rs
@@ -29,7 +29,6 @@ fn main() {
             y[0] = x[0] / v;
             y[1] = x[1] / v2;
         },
-        (2, 2),
     );
 
     let params = Parameters::new()
diff --git a/examples/neely/main.rs b/examples/neely/main.rs
index 3bcb332a4..d74d4ede4 100644
--- a/examples/neely/main.rs
+++ b/examples/neely/main.rs
@@ -52,7 +52,6 @@ fn main() {
             y[1] = x[2] / vm1;
             y[2] = x[3] / vm2;
         },
-        (4, 3),
     );
     let params = Parameters::new()
         .add("cls", 0.0, 0.4)
diff --git a/examples/new_iov/main.rs b/examples/new_iov/main.rs
index f60d0d683..4ac298e83 100644
--- a/examples/new_iov/main.rs
+++ b/examples/new_iov/main.rs
@@ -25,7 +25,6 @@ fn main() {
             fetch_params!(p, _ke0, _ske);
             y[0] = x[0] / 50.0;
         },
-        (2, 1),
         11,
     );
 
diff --git a/examples/theophylline/main.rs b/examples/theophylline/main.rs
index 774df5a69..9addb5e2f 100644
--- a/examples/theophylline/main.rs
+++ b/examples/theophylline/main.rs
@@ -1,23 +1,7 @@
 use pmcore::prelude::*;
 
 fn main() {
-    // let eq = Equation::new_ode(
-    //     |x, p, _t, dx, rateiv, _cov| {
-    //         // fetch_cov!(cov, t, wt);
-    //         fetch_params!(p, ka, ke, _v);
-    //         dx[0] = -ka * x[0];
-    //         dx[1] = ka * x[0] - ke * x[1];
-    //     },
-    //     |_p, _t, _cov| lag! {},
-    //     |_p, _t, _cov| fa! {},
-    //     |_p, _t, _cov, _x| {},
-    //     |x, p, _t, _cov, y| {
-    //         fetch_params!(p, _ka, _ke, v);
-    //         y[0] = x[1] * 1000.0 / v;
-    //     },
-    //     (2, 1),
-    // );
-    let eq = equation::Analytical::new(
+    let analytical = equation::Analytical::new(
         one_compartment_with_absorption,
         |_p, _t, _cov| {},
         |_p, _t, _cov| lag! {},
@@ -27,7 +11,6 @@ fn main() {
             fetch_params!(p, _ka, _ke, v);
             y[0] = x[1] * 1000.0 / v;
         },
-        (2, 1),
     );
 
     let params = Parameters::new()
@@ -50,7 +33,7 @@ fn main() {
 
     settings.initialize_logs().unwrap();
     let data = data::read_pmetrics("examples/theophylline/theophylline.csv").unwrap();
-    let mut algorithm = dispatch_algorithm(settings, eq, data).unwrap();
+    let mut algorithm = dispatch_algorithm(settings, analytical, data).unwrap();
     // let result = algorithm.fit().unwrap();
     algorithm.initialize().unwrap();
     let mut result = algorithm.fit().unwrap();
diff --git a/examples/two_eq_lag/main.rs b/examples/two_eq_lag/main.rs
index eebeef1c1..feb4f4766 100644
--- a/examples/two_eq_lag/main.rs
+++ b/examples/two_eq_lag/main.rs
@@ -22,7 +22,6 @@ fn main() {
             fetch_params!(p, _ka, _ke, _tlag, v);
             y[0] = x[1] / v;
         },
-        (3, 1),
     );
     // let eq = Equation::new_analytical(
     //     one_compartment_with_absorption,
diff --git a/examples/vanco.rs b/examples/vanco.rs
index 3a8a0d1d7..f17721cfd 100644
--- a/examples/vanco.rs
+++ b/examples/vanco.rs
@@ -14,7 +14,6 @@ fn main() {
         |x, _p, _t, _cov, y| {
             y[0] = x[1];
         },
-        (2, 1),
     );
     // same eq but analytical
     // let eq = Equation::new_analytical(
diff --git a/examples/vanco_sde/main.rs b/examples/vanco_sde/main.rs
index 29533be4e..024659c4e 100644
--- a/examples/vanco_sde/main.rs
+++ b/examples/vanco_sde/main.rs
@@ -26,7 +26,6 @@ fn main() {
             fetch_cov!(cov, t, wt);
             y[0] = x[1] / (vol * wt);
         },
-        (4, 1),
         100,
     );
 
diff --git a/tests/bestdose_tests.rs b/tests/bestdose_tests.rs
index 6033a6d93..aa1202824 100644
--- a/tests/bestdose_tests.rs
+++ b/tests/bestdose_tests.rs
@@ -21,7 +21,6 @@ fn test_infusion_mask_inclusion() -> Result<()> {
             fetch_params!(p, _ke, v);
             y[0] = x[0] / v;
         },
-        (1, 1),
     );
 
     let params = Parameters::new().add("ke", 0.1, 0.5).add("v", 40.0, 60.0);
@@ -136,7 +135,6 @@ fn test_fixed_infusion_preservation() -> Result<()> {
             fetch_params!(p, _ke, v);
             y[0] = x[0] / v;
         },
-        (1, 1),
     );
 
     let params = Parameters::new()
@@ -226,7 +224,6 @@ fn test_dose_count_validation() -> Result<()> {
             fetch_params!(p, _ke, v);
             y[0] = x[0] / v;
         },
-        (1, 1),
     );
 
     let params = Parameters::new().add("ke", 0.1, 0.5).add("v", 40.0, 60.0);
@@ -306,7 +303,6 @@ fn test_empty_observations_validation() -> Result<()> {
             fetch_params!(p, _ke, v);
             y[0] = x[0] / v;
         },
-        (1, 1),
     );
 
     let params = Parameters::new().add("ke", 0.1, 0.5).add("v", 40.0, 60.0);
@@ -372,7 +368,6 @@ fn test_basic_auc_mode() -> Result<()> {
             fetch_params!(p, _ke, v);
             y[0] = x[0] / v;
         },
-        (1, 1),
     );
 
     let params = Parameters::new().add("ke", 0.1, 0.5).add("v", 40.0, 60.0);
@@ -463,7 +458,6 @@ fn test_infusion_auc_mode() -> Result<()> {
             fetch_params!(p, _ke, v);
             y[0] = x[0] / v;
         },
-        (1, 1),
     );
 
     let params = Parameters::new().add("ke", 0.1, 0.5).add("v", 40.0, 60.0);
@@ -573,7 +567,6 @@ fn test_multi_outeq_auc_mode() -> Result<()> {
             y[0] = x[0] / v; // outeq 0: concentration
             y[1] = x[0]; // outeq 1: amount
         },
-        (1, 2),
     );
 
     let params = Parameters::new().add("ke", 0.1, 0.5).add("v", 40.0, 60.0);
@@ -647,7 +640,6 @@ fn test_multi_outeq_auc_optimization() -> Result<()> {
             y[0] = x[0] / v;
             y[1] = x[0];
         },
-        (1, 2),
     );
 
     let params = Parameters::new().add("ke", 0.1, 0.5).add("v", 40.0, 60.0);
@@ -738,7 +730,6 @@ fn test_auc_from_zero_single_dose() -> Result<()> {
             fetch_params!(p, _ke, v);
             y[0] = x[0] / v;
         },
-        (1, 1),
     );
 
     let params = Parameters::new().add("ke", 0.2, 0.4).add("v", 40.0, 60.0);
@@ -828,7 +819,6 @@ fn test_auc_from_last_dose_maintenance() -> Result<()> {
             fetch_params!(p, _ke, v);
             y[0] = x[0] / v;
         },
-        (1, 1),
     );
 
     let params = Parameters::new().add("ke", 0.2, 0.4).add("v", 40.0, 60.0);
@@ -922,7 +912,6 @@ fn test_auc_modes_comparison() -> Result<()> {
             fetch_params!(p, _ke, v);
             y[0] = x[0] / v;
         },
-        (1, 1),
     );
 
     let params = Parameters::new().add("ke", 0.3, 0.3).add("v", 50.0, 50.0);
@@ -1055,7 +1044,6 @@ fn test_auc_from_last_dose_multiple_observations() -> Result<()> {
             fetch_params!(p, _ke, v);
             y[0] = x[0] / v;
         },
-        (1, 1),
     );
 
     let params = Parameters::new().add("ke", 0.2, 0.4).add("v", 40.0, 60.0);
@@ -1155,7 +1143,6 @@ fn test_auc_from_last_dose_no_prior_dose() -> Result<()> {
             fetch_params!(p, _ke, v);
             y[0] = x[0] / v;
         },
-        (1, 1),
     );
 
     let params = Parameters::new().add("ke", 0.2, 0.4).add("v", 40.0, 60.0);
@@ -1250,7 +1237,6 @@ fn test_dose_range_bounds_respected() -> Result<()> {
             fetch_params!(p, _ke, v);
             y[0] = x[0] / v;
         },
-        (1, 1),
     );
 
     let params = Parameters::new().add("ke", 0.1, 0.5).add("v", 40.0, 60.0);
diff --git a/tests/onecomp.rs b/tests/onecomp.rs
index 7ab278dd6..58230aca0 100644
--- a/tests/onecomp.rs
+++ b/tests/onecomp.rs
@@ -16,7 +16,6 @@ fn test_one_compartment_npag() -> Result<()> {
             fetch_params!(p, v);
             y[0] = x[0] / v;
         },
-        (1, 1),
     );
 
     // Define parameters
@@ -84,7 +83,6 @@ fn test_one_compartment_npod() -> Result<()> {
             fetch_params!(p, v);
             y[0] = x[0] / v;
         },
-        (1, 1),
     );
 
     // Define parameters
@@ -152,7 +150,6 @@ fn test_one_compartment_postprob() -> Result<()> {
             fetch_params!(p, v);
             y[0] = x[0] / v;
         },
-        (1, 1),
     );
 
     // Define parameters

From 3d717ac5e077330ab69c98fa472542af285ea0a7 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Juli=C3=A1n=20D=2E=20Ot=C3=A1lvaro?=
 <juliandavid347@gmail.com>
Date: Tue, 24 Mar 2026 21:07:35 +0000
Subject: [PATCH 9/9] All ODE examples are updated to use the new ode! macro

---
 Cargo.toml                  |  2 +-
 examples/bestdose.rs        | 11 ++++-------
 examples/bestdose_auc.rs    | 11 ++++-------
 examples/bestdose_bounds.rs | 11 ++++-------
 examples/bimodal_ke/main.rs | 12 +++++-------
 examples/drusano/main.rs    | 12 +++++-------
 examples/iov/main.rs        |  5 ++++-
 examples/meta/main.rs       | 11 ++++-------
 examples/neely/main.rs      | 15 ++++-----------
 examples/new_iov/main.rs    |  5 ++++-
 examples/two_eq_lag/main.rs | 12 +++++-------
 examples/vanco.rs           | 14 ++++++--------
 12 files changed, 50 insertions(+), 71 deletions(-)

diff --git a/Cargo.toml b/Cargo.toml
index df35f74df..94406d884 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -28,7 +28,7 @@ tracing-subscriber = { version = "0.3.19", features = [
     "time",
 ] }
 faer = "0.23.1"
-pharmsol = "=0.24.0"
+pharmsol = "=0.24.1"
 rand = "0.9.0"
 anyhow = "1.0.100"
 rayon = "1.10.0"
diff --git a/examples/bestdose.rs b/examples/bestdose.rs
index 3114ea3bc..34ff95364 100644
--- a/examples/bestdose.rs
+++ b/examples/bestdose.rs
@@ -7,20 +7,17 @@ use pmcore::routines::initialization::parse_prior;
 
 fn main() -> Result<()> {
     // Example model
-    let eq = equation::ODE::new(
-        |x, p, _t, dx, b, _rateiv, _cov| {
+    let eq = ode! {
+        diffeq: |x, p, _t, dx, b, _rateiv, _cov| {
             // fetch_cov!(cov, t, wt);
             fetch_params!(p, ke, _v);
             dx[0] = -ke * x[0] + b[0];
         },
-        |_p, _, _| lag! {},
-        |_p, _, _| fa! {},
-        |_p, _t, _cov, _x| {},
-        |x, p, _t, _cov, y| {
+        out: |x, p, _t, _cov, y| {
             fetch_params!(p, _ke, v);
             y[0] = x[0] / v;
         },
-    );
+    };
 
     let params = Parameters::new()
         .add("ke", 0.001, 3.0)
diff --git a/examples/bestdose_auc.rs b/examples/bestdose_auc.rs
index 75ac3ee12..6f9b08c27 100644
--- a/examples/bestdose_auc.rs
+++ b/examples/bestdose_auc.rs
@@ -10,19 +10,16 @@ fn main() -> Result<()> {
     println!("======================================\n");
 
     // Simple one-compartment PK model
-    let eq = equation::ODE::new(
-        |x, p, _t, dx, b, _rateiv, _cov| {
+    let eq = ode! {
+        diffeq: |x, p, _t, dx, b, _rateiv, _cov| {
             fetch_params!(p, ke, _v);
             dx[0] = -ke * x[0] + b[0];
         },
-        |_p, _, _| lag! {},
-        |_p, _, _| fa! {},
-        |_p, _t, _cov, _x| {},
-        |x, p, _t, _cov, y| {
+        out: |x, p, _t, _cov, y| {
             fetch_params!(p, _ke, v);
             y[0] = x[0] / v;
         },
-    );
+    };
 
     // Minimal parameter ranges
     let params = Parameters::new()
diff --git a/examples/bestdose_bounds.rs b/examples/bestdose_bounds.rs
index 6058b1cf3..d91709b39 100644
--- a/examples/bestdose_bounds.rs
+++ b/examples/bestdose_bounds.rs
@@ -10,19 +10,16 @@ fn main() -> Result<()> {
     println!("==========================================\n");
 
     // Simple one-compartment PK model
-    let eq = equation::ODE::new(
-        |x, p, _t, dx, b, _rateiv, _cov| {
+    let eq = ode! {
+        diffeq: |x, p, _t, dx, b, _rateiv, _cov| {
             fetch_params!(p, ke, _v);
             dx[0] = -ke * x[0] + b[0];
         },
-        |_p, _, _| lag! {},
-        |_p, _, _| fa! {},
-        |_p, _t, _cov, _x| {},
-        |x, p, _t, _cov, y| {
+        out: |x, p, _t, _cov, y| {
             fetch_params!(p, _ke, v);
             y[0] = x[0] / v;
         },
-    );
+    };
 
     let params = Parameters::new()
         .add("ke", 0.001, 3.0)
diff --git a/examples/bimodal_ke/main.rs b/examples/bimodal_ke/main.rs
index dfeac2b22..0ea9bb68c 100644
--- a/examples/bimodal_ke/main.rs
+++ b/examples/bimodal_ke/main.rs
@@ -2,20 +2,18 @@ use anyhow::Result;
 use pmcore::prelude::*;
 
 fn main() -> Result<()> {
-    let eq = equation::ODE::new(
-        |x, p, _t, dx, b, rateiv, _cov| {
+    let eq = ode! {
+        diffeq: |x, p, _t, dx, b, rateiv, _cov| {
             // fetch_cov!(cov, t, wt);
             fetch_params!(p, ke, _v);
             dx[0] = -ke * x[0] + rateiv[1] + b[1];
         },
-        |_p, _t, _cov| lag! {},
-        |_p, _t, _cov| fa! {},
-        |_p, _t, _cov, _x| {},
-        |x, p, _t, _cov, y| {
+        out: |x, p, _t, _cov, y| {
             fetch_params!(p, _ke, v);
             y[1] = x[0] / v;
         },
-    );
+    }
+    .with_solver(OdeSolver::ExplicitRk(ExplicitRkTableau::Tsit45));
 
     let params = Parameters::new()
         .add("ke", 0.001, 3.0)
diff --git a/examples/drusano/main.rs b/examples/drusano/main.rs
index 0e030302d..640fb20eb 100644
--- a/examples/drusano/main.rs
+++ b/examples/drusano/main.rs
@@ -2,8 +2,8 @@ use anyhow::Result;
 use pmcore::prelude::*;
 #[allow(unused_variables)]
 fn main() -> Result<()> {
-    let eq = equation::ODE::new(
-        |x, p, t, dx, b, rateiv, _cov| {
+    let eq = ode! {
+        diffeq:|x, p, t, dx, b, rateiv, _cov| {
             fetch_params!(
                 p, v1, cl1, v2, cl2, popmax, kgs, kks, e50_1s, e50_2s, alpha_s, kgr1, kkr1,
                 e50_1r1, alpha_r1, kgr2, kkr2, e50_2r2, alpha_r2, init_4, init_5, h1s, h2s, h1r1,
@@ -44,9 +44,7 @@ fn main() -> Result<()> {
             let xm0best = get_e2(u_r2, v_r2, w_r2, 1.0 / h1r2, 1.0 / h2r2, alpha_s);
             dx[4] = xnr2 * (kgr2 * e - kkr2 * xm0best);
         },
-        |_p, _t, _cov| lag! {},
-        |_p, _t, _cov| fa! {},
-        |p, t, cov, x| {
+        init: |p, t, cov, x| {
             fetch_params!(
                 p, v1, cl1, v2, cl2, popmax, kgs, kks, e50_1s, e50_2s, alpha_s, kgr1, kkr1,
                 e50_1r1, alpha_r1, kgr2, kkr2, e50_2r2, alpha_r2, init_4, init_5, h1s, h2s, h1r1,
@@ -59,7 +57,7 @@ fn main() -> Result<()> {
             x[3] = 10.0_f64.powf(init_4);
             x[4] = 10.0_f64.powf(init_5);
         },
-        |x, p, _t, _cov, y| {
+        out: |x, p, _t, _cov, y| {
             fetch_params!(
                 p, v1, cl1, v2, cl2, popmax, kgs, kks, e50_1s, e50_2s, alpha_s, kgr1, kkr1,
                 e50_1r1, alpha_r1, kgr2, kkr2, e50_2r2, alpha_r2, init_4, init_5, h1s, h2s, h1r1,
@@ -71,7 +69,7 @@ fn main() -> Result<()> {
             y[3] = x[3].log10();
             y[4] = x[4].log10();
         },
-    );
+    };
 
     let params = Parameters::new()
         .add("v1", 5.0, 160.0)
diff --git a/examples/iov/main.rs b/examples/iov/main.rs
index 023762ab9..222015dfc 100644
--- a/examples/iov/main.rs
+++ b/examples/iov/main.rs
@@ -26,7 +26,10 @@ fn main() -> Result<()> {
             y[0] = x[0] / 50.0;
         },
         10000,
-    );
+    )
+    .with_nstates(2)
+    .with_ndrugs(1)
+    .with_nout(1);
 
     let params = Parameters::new().add("ke0", 0.001, 2.0);
 
diff --git a/examples/meta/main.rs b/examples/meta/main.rs
index 20fbb51c9..4c3fa264d 100644
--- a/examples/meta/main.rs
+++ b/examples/meta/main.rs
@@ -5,8 +5,8 @@
 use pmcore::{prelude::*, routines::settings};
 
 fn main() {
-    let eq = equation::ODE::new(
-        |x, p, t, dx, b, rateiv, cov| {
+    let eq = ode! {
+        diffeq: |x, p, t, dx, b, rateiv, cov| {
             fetch_cov!(cov, t, wt, pkvisit);
             fetch_params!(p, cls, fm, k20, relv, theta1, theta2, vs);
             let cl = cls * ((pkvisit - 1.0) * theta1).exp() * (wt / 70.0).powf(0.75);
@@ -16,10 +16,7 @@ fn main() {
             dx[0] = rateiv[0] - ke * x[0] * (1.0 - fm) - fm * x[0] + b[0];
             dx[1] = fm * x[0] - k20 * x[1];
         },
-        |_p, _t, _cov| lag! {},
-        |_p, _t, _cov| fa! {},
-        |_p, _t, _cov, _x| {},
-        |x, p, t, cov, y| {
+        out: |x, p, t, cov, y| {
             fetch_cov!(cov, t, wt, pkvisit);
             fetch_params!(p, cls, fm, k20, relv, theta1, theta2, vs);
             let cl = cls * ((pkvisit - 1.0) * theta1).exp() * (wt / 70.0).powf(0.75);
@@ -29,7 +26,7 @@ fn main() {
             y[0] = x[0] / v;
             y[1] = x[1] / v2;
         },
-    );
+    };
 
     let params = Parameters::new()
         .add("cls", 0.1, 10.0)
diff --git a/examples/neely/main.rs b/examples/neely/main.rs
index d74d4ede4..98b47cbee 100644
--- a/examples/neely/main.rs
+++ b/examples/neely/main.rs
@@ -1,7 +1,7 @@
 use pmcore::prelude::*;
 fn main() {
-    let ode = equation::ODE::new(
-        |x, p, t, dx, b, rateiv, cov| {
+    let ode = ode! {
+        diffeq: |x, p, t, dx, b, rateiv, cov| {
             fetch_params!(p, cls, k30, k40, qs, vps, vs, fm1, fm2, theta1, theta2);
             fetch_cov!(cov, t, wt, pkvisit);
 
@@ -25,14 +25,7 @@ fn main() {
             dx[2] = fm1 * x[0] - k30 * x[2];
             dx[3] = fm2 * x[0] - k40 * x[3];
         },
-        |_p, _t, _cov| {
-            lag! {}
-        },
-        |_p, _t, _cov| {
-            fa! {}
-        },
-        |_p, _t, _cov, _x| {},
-        |x, p, t, cov, y| {
+        out: |x, p, t, cov, y| {
             fetch_params!(p, cls, _k30, _k40, qs, vps, vs, _fm1, _fm2, theta1, theta2);
             fetch_cov!(cov, t, wt, pkvisit);
 
@@ -52,7 +45,7 @@ fn main() {
             y[1] = x[2] / vm1;
             y[2] = x[3] / vm2;
         },
-    );
+    };
     let params = Parameters::new()
         .add("cls", 0.0, 0.4)
         .add("k30", 0.0, 0.5)
diff --git a/examples/new_iov/main.rs b/examples/new_iov/main.rs
index 4ac298e83..40490b7f2 100644
--- a/examples/new_iov/main.rs
+++ b/examples/new_iov/main.rs
@@ -26,7 +26,10 @@ fn main() {
             y[0] = x[0] / 50.0;
         },
         11,
-    );
+    )
+    .with_nstates(2)
+    .with_ndrugs(1)
+    .with_nout(1);
 
     let params = Parameters::new()
         .add("ke0", 0.0001, 2.4)
diff --git a/examples/two_eq_lag/main.rs b/examples/two_eq_lag/main.rs
index feb4f4766..ec1c06d65 100644
--- a/examples/two_eq_lag/main.rs
+++ b/examples/two_eq_lag/main.rs
@@ -5,24 +5,22 @@
 use pmcore::prelude::*;
 
 fn main() {
-    let eq = equation::ODE::new(
-        |x, p, _t, dx, b, rateiv, _cov| {
+    let eq = ode! {
+        diffeq: |x, p, _t, dx, b, rateiv, _cov| {
             fetch_cov!(cov, t,);
             fetch_params!(p, ka, ke);
             dx[0] = -ka * x[0] + b[0];
             dx[1] = ka * x[0] - ke * x[1];
         },
-        |p, _t, _cov| {
+        lag: |p, _t, _cov| {
             fetch_params!(p, _ka, _ke, tlag, _v);
             lag! {0=>tlag}
         },
-        |_p, _t, _cov| fa! {},
-        |_p, _t, _cov, _x| {},
-        |x, p, _t, _cov, y| {
+        out: |x, p, _t, _cov, y| {
             fetch_params!(p, _ka, _ke, _tlag, v);
             y[0] = x[1] / v;
         },
-    );
+    };
     // let eq = Equation::new_analytical(
     //     one_compartment_with_absorption,
     //     |_p, _cov| {},
diff --git a/examples/vanco.rs b/examples/vanco.rs
index f17721cfd..f76b60c67 100644
--- a/examples/vanco.rs
+++ b/examples/vanco.rs
@@ -1,20 +1,18 @@
-use pmcore::prelude::{simulator::Equation, *};
+use pmcore::prelude::*;
 fn main() {
-    let eq = equation::ODE::new(
-        |x, p, _t, dx, b, _rateiv, _cov| {
+    let eq = ode! {
+        diffeq: |x, p, _t, dx, b, _rateiv, _cov| {
             fetch_params!(p, ke, kcp, kpc);
             dx[0] = -ke * x[0] - kcp * x[0] + kpc * x[1] + b[0];
             dx[1] = -kpc * x[1] + kcp * x[0];
         },
-        |_p, _t, _cov| lag! {},
-        |_p, _t, _cov| fa! {},
-        |_p, _t, _cov, x| {
+        init: |_p, _t, _cov, x| {
             x[0] = 500.0;
         },
-        |x, _p, _t, _cov, y| {
+        out: |x, _p, _t, _cov, y| {
             y[0] = x[1];
         },
-    );
+    };
     // same eq but analytical
     // let eq = Equation::new_analytical(
     //     two_compartments,