feat: Update representation of Predictions

Cargo.toml

@@ -17,7 +17,6 @@ csv = "1.3.0"
 diffsol = "=0.7.0"
 libloading = { version = "0.8.6", optional = true, features = [] }
 nalgebra = "0.34.1"
-ndarray = { version = "0.17.1", features = ["rayon"] }
 rand = "0.9.0"
 rand_distr = "0.5.0"
 rayon = "1.10.0"
@@ -28,6 +27,7 @@ thiserror = "2.0.11"
 argmin = "0.11.0"
 argmin-math = "0.5.1"
 tracing = "0.1.41"
+faer = "0.23.2"
 
 [dev-dependencies]
 criterion = { version = "0.7.0", features = ["html_reports"] }

examples/gendata.rs

@@ -63,7 +63,7 @@ fn main() {
     let mut data = vec![];
     for (i, spp) in support_points.iter().enumerate() {
         let trajectories = sde.estimate_predictions(&subject, spp).unwrap();
-        let trajectory = trajectories.row(0);
+        let trajectory = trajectories.row(0).unwrap();
         // dbg!(&trajectory);
         let mut sb = data::Subject::builder(format!("id{}", i)).bolus(0.0, 20.0, 0);
         for (t, point) in trajectory.iter().enumerate() {

examples/likelihood.rs

@@ -1,12 +1,10 @@
-use ndarray::Array1;
 const FRAC_1_SQRT_2PI: f64 =
     std::f64::consts::FRAC_2_SQRT_PI * std::f64::consts::FRAC_1_SQRT_2 / 2.0;
-fn main() {
-    use ndarray::array;
 
-    let predictions = array![1.0, 2.0, 3.0];
-    let observations = array![1.1, 1.9, 3.2];
-    let sigma = array![0.1, 0.2, 0.3];
+fn main() {
+    let predictions = vec![1.0, 2.0, 3.0];
+    let observations = vec![1.1, 1.9, 3.2];
+    let sigma = vec![0.1, 0.2, 0.3];
 
     let result1 = normal_likelihood1(&predictions, &observations, &sigma);
     let result2 = normal_likelihood2(&predictions, &observations, &sigma);
@@ -15,32 +13,30 @@ fn main() {
     println!("Result 2: {}", result2);
 }
 
-fn normal_likelihood1(
-    predictions: &Array1<f64>,
-    observations: &Array1<f64>,
-    sigma: &Array1<f64>,
-) -> f64 {
+fn normal_likelihood1(predictions: &[f64], observations: &[f64], sigma: &[f64]) -> f64 {
     const FRAC_1_SQRT_2PI: f64 =
         std::f64::consts::FRAC_2_SQRT_PI * std::f64::consts::FRAC_1_SQRT_2 / 2.0;
-    let diff = (observations - predictions).mapv(|x| x.powi(2));
-    let two_sigma_sq = sigma.mapv(|x| 2.0 * x * x);
-    let exponent = (-&diff / two_sigma_sq).mapv(|x| x.exp());
-    let aux_vec = FRAC_1_SQRT_2PI * exponent / sigma;
-    aux_vec.product()
+
+    predictions
+        .iter()
+        .zip(observations.iter())
+        .zip(sigma.iter())
+        .map(|((pred, obs), sig)| {
+            let diff = (obs - pred).powi(2);
+            let two_sigma_sq = 2.0 * sig * sig;
+            let exponent = (-diff / two_sigma_sq).exp();
+            FRAC_1_SQRT_2PI * exponent / sig
+        })
+        .product()
 }
 
-fn normal_likelihood2(
-    predictions: &Array1<f64>,
-    observations: &Array1<f64>,
-    sigma: &Array1<f64>,
-) -> f64 {
-    let aux_vec = predictions
+fn normal_likelihood2(predictions: &[f64], observations: &[f64], sigma: &[f64]) -> f64 {
+    predictions
         .iter()
         .zip(observations.iter())
         .zip(sigma.iter())
         .map(|((pred, obs), sig)| normpdf(*obs, *pred, *sig))
-        .collect::<Array1<f64>>();
-    aux_vec.product()
+        .product()
 }
 
 /// Probability density function

src/error/mod.rs

@@ -5,16 +5,12 @@ use crate::data::parser::pmetrics::PmetricsError;
 
 use crate::CovariateError;
 
-use ndarray::ShapeError;
-
 #[derive(Error, Debug, Clone)]
 pub enum PharmsolError {
     #[error("Error in the error model: {0}")]
     ErrorModelError(#[from] ErrorModelError),
     #[error("Covariate error: {0}")]
     CovariateError(#[from] CovariateError),
-    #[error("Shape error: {0}")]
-    NdarrayShapeError(#[from] ShapeError),
     #[error("Error parsing Pmetrics datafile: {0}")]
     PmetricsError(#[from] PmetricsError),
     #[error("Diffsol error: {0}")]

src/lib.rs

@@ -30,7 +30,7 @@ pub mod prelude {
         pub use crate::simulator::{
             equation,
             equation::Equation,
-            likelihood::{psi, PopulationPredictions, Prediction, SubjectPredictions},
+            likelihood::{psi, Prediction, PredictionMatrix, SubjectPredictions},
         };
     }
     pub mod models {

src/optimize/spp.rs

@@ -3,22 +3,23 @@ use argmin::{
     solver::neldermead::NelderMead,
 };
 
-use ndarray::{Array1, Axis};
+use faer::Mat;
 
 use crate::{prelude::simulator::psi, Data, Equation, ErrorModels};
 
 pub struct SppOptimizer<'a, E: Equation> {
     equation: &'a E,
     data: &'a Data,
     sig: &'a ErrorModels,
-    pyl: &'a Array1<f64>,
+    pyl: &'a Vec<f64>,
 }
 
 impl<E: Equation> CostFunction for SppOptimizer<'_, E> {
     type Param = Vec<f64>;
     type Output = f64;
     fn cost(&self, spp: &Self::Param) -> Result<Self::Output, Error> {
-        let theta = Array1::from(spp.clone()).insert_axis(Axis(0));
+        // Create a 1xN matrix (1 support point with N parameters)
+        let theta = Mat::from_fn(1, spp.len(), |_, j| spp[j]);
 
         let psi = psi(self.equation, self.data, &theta, self.sig, false, false)?;
 
@@ -34,35 +35,32 @@ impl<E: Equation> CostFunction for SppOptimizer<'_, E> {
         }
         let nsub = psi.nrows() as f64;
         let mut sum = -nsub;
-        for (p_i, pyl_i) in psi.iter().zip(self.pyl.iter()) {
+        // Iterate through each row of psi (which has only 1 column)
+        for (i, pyl_i) in self.pyl.iter().enumerate() {
+            let p_i = psi[(i, 0)];
             sum += p_i / pyl_i;
         }
         Ok(-sum)
     }
 }
 
 impl<'a, E: Equation> SppOptimizer<'a, E> {
-    pub fn new(
-        equation: &'a E,
-        data: &'a Data,
-        sig: &'a ErrorModels,
-        pyl: &'a Array1<f64>,
-    ) -> Self {
+    pub fn new(equation: &'a E, data: &'a Data, sig: &'a ErrorModels, pyl: &'a Vec<f64>) -> Self {
         Self {
             equation,
             data,
             sig,
             pyl,
         }
     }
-    pub fn optimize_point(self, spp: Array1<f64>) -> Result<Array1<f64>, Error> {
-        let simplex = create_initial_simplex(&spp.to_vec());
+    pub fn optimize_point(self, spp: Vec<f64>) -> Result<Vec<f64>, Error> {
+        let simplex = create_initial_simplex(&spp);
         let solver: NelderMead<Vec<f64>, f64> = NelderMead::new(simplex).with_sd_tolerance(1e-2)?;
         let res = Executor::new(self, solver)
             .configure(|state| state.max_iters(5))
             // .add_observer(SlogLogger::term(), ObserverMode::Always)
             .run()?;
-        Ok(Array1::from(res.state.best_param.unwrap()))
+        Ok(res.state.best_param.unwrap())
     }
 }
 

src/simulator/equation/sde/mod.rs

@@ -2,7 +2,6 @@ mod em;
 
 use diffsol::{NalgebraContext, Vector};
 use nalgebra::DVector;
-use ndarray::{concatenate, Array2, Axis};
 use rand::{rng, Rng};
 use rayon::prelude::*;
 
@@ -13,7 +12,7 @@ use crate::{
     data::{Covariates, Infusion},
     error_model::ErrorModels,
     prelude::simulator::Prediction,
-    simulator::{Diffusion, Drift, Fa, Init, Lag, Neqs, Out, V},
+    simulator::{likelihood::PredictionMatrix, Diffusion, Drift, Fa, Init, Lag, Neqs, Out, V},
     Subject,
 };
 
@@ -151,9 +150,9 @@ impl State for Vec<DVector<f64>> {
 /// Predictions implementation for particle-based SDE simulation outputs.
 ///
 /// This implementation manages and processes predictions from multiple particles.
-impl Predictions for Array2<Prediction> {
+impl Predictions for PredictionMatrix {
     fn new(nparticles: usize) -> Self {
-        Array2::from_shape_fn((nparticles, 0), |_| Prediction::default())
+        PredictionMatrix::new(nparticles, 0)
     }
     fn squared_error(&self) -> f64 {
         unimplemented!();
@@ -171,11 +170,11 @@ impl Predictions for Array2<Prediction> {
 
             let mean_prediction: f64 = column
                 .iter()
-                .map(|pred: &Prediction| pred.prediction())
+                .map(|pred: &&Prediction| pred.prediction())
                 .sum::<f64>()
                 / self.nrows() as f64;
 
-            let mut prediction = column.first().unwrap().clone();
+            let mut prediction = (*column.first().unwrap()).clone();
             prediction.set_prediction(mean_prediction);
             result.push(prediction);
         }
@@ -186,7 +185,7 @@ impl Predictions for Array2<Prediction> {
 
 impl EquationTypes for SDE {
     type S = Vec<DVector<f64>>; // Vec -> particles, DVector -> state
-    type P = Array2<Prediction>; // Rows -> particles, Columns -> time
+    type P = PredictionMatrix; // Rows -> particles, Columns -> time
 }
 
 impl EquationPriv for SDE {
@@ -286,8 +285,10 @@ impl EquationPriv for SDE {
             );
             *p = observation.to_prediction(y[observation.outeq()], x[i].as_slice().to_vec());
         });
-        let out = Array2::from_shape_vec((self.nparticles, 1), pred.clone())?;
-        *output = concatenate(Axis(1), &[output.view(), out.view()]).unwrap();
+
+        // Append the new column of predictions
+        output.append_column(pred.clone())?;
+
         //e = y[t] .- x[:,1]
         // q = pdf.(Distributions.Normal(0, 0.5), e)
         if let Some(em) = error_models {