feat: Refactor outputs

src/algorithms/npag.rs

@@ -5,9 +5,10 @@ pub use crate::routines::evaluation::ipm::burke;
 pub use crate::routines::evaluation::qr;
 use crate::routines::settings::Settings;
 
-use crate::routines::output::{CycleLog, NPCycle, NPResult};
+use crate::routines::output::{cycles::CycleLog, cycles::NPCycle, NPResult};
 use crate::structs::psi::{calculate_psi, Psi};
 use crate::structs::theta::Theta;
+use crate::structs::weights::Weights;
 
 use anyhow::bail;
 use anyhow::Result;
@@ -18,8 +19,6 @@ use pharmsol::prelude::{
 
 use pharmsol::prelude::ErrorModel;
 
-use faer::Col;
-
 use crate::routines::initialization;
 
 use crate::routines::expansion::adaptative_grid::adaptative_grid;
@@ -35,8 +34,8 @@ pub struct NPAG<E: Equation> {
     ranges: Vec<(f64, f64)>,
     psi: Psi,
     theta: Theta,
-    lambda: Col<f64>,
-    w: Col<f64>,
+    lambda: Weights,
+    w: Weights,
     eps: f64,
     last_objf: f64,
     objf: f64,
@@ -59,8 +58,8 @@ impl<E: Equation> Algorithms<E> for NPAG<E> {
             ranges: settings.parameters().ranges(),
             psi: Psi::new(),
             theta: Theta::new(),
-            lambda: Col::zeros(0),
-            w: Col::zeros(0),
+            lambda: Weights::default(),
+            w: Weights::default(),
             eps: 0.2,
             last_objf: -1e30,
             objf: f64::NEG_INFINITY,
@@ -138,7 +137,7 @@ impl<E: Equation> Algorithms<E> for NPAG<E> {
         if (self.last_objf - self.objf).abs() <= THETA_G && self.eps > THETA_E {
             self.eps /= 2.;
             if self.eps <= THETA_E {
-                let pyl = psi * w;
+                let pyl = psi * w.weights();
                 self.f1 = pyl.iter().map(|x| x.ln()).sum();
                 if (self.f1 - self.f0).abs() <= THETA_F {
                     tracing::info!("The model converged after {} cycles", self.cycle,);
@@ -165,15 +164,15 @@ impl<E: Equation> Algorithms<E> for NPAG<E> {
         }
 
         // Create state object
-        let state = NPCycle {
-            cycle: self.cycle,
-            objf: -2. * self.objf,
-            delta_objf: (self.last_objf - self.objf).abs(),
-            nspp: self.theta.nspp(),
-            theta: self.theta.clone(),
-            error_models: self.error_models.clone(),
-            status: self.status.clone(),
-        };
+        let state = NPCycle::new(
+            self.cycle,
+            -2. * self.objf,
+            self.error_models.clone(),
+            self.theta.clone(),
+            self.theta.nspp(),
+            (self.last_objf - self.objf).abs(),
+            self.status.clone(),
+        );
 
         // Write cycle log
         self.cycle_log.push(state);
@@ -199,7 +198,7 @@ impl<E: Equation> Algorithms<E> for NPAG<E> {
         }
 
         (self.lambda, _) = match burke(&self.psi) {
-            Ok((lambda, objf)) => (lambda, objf),
+            Ok((lambda, objf)) => (lambda.into(), objf),
             Err(err) => {
                 bail!("Error in IPM during evaluation: {:?}", err);
             }
@@ -213,11 +212,11 @@ impl<E: Equation> Algorithms<E> for NPAG<E> {
         let max_lambda = self
             .lambda
             .iter()
-            .fold(f64::NEG_INFINITY, |acc, &x| x.max(acc));
+            .fold(f64::NEG_INFINITY, |acc, x| x.max(acc));
 
         let mut keep = Vec::<usize>::new();
         for (index, lam) in self.lambda.iter().enumerate() {
-            if *lam > max_lambda / 1000_f64 {
+            if lam > max_lambda / 1000_f64 {
                 keep.push(index);
             }
         }
@@ -262,15 +261,15 @@ impl<E: Equation> Algorithms<E> for NPAG<E> {
 
         self.validate_psi()?;
         (self.lambda, self.objf) = match burke(&self.psi) {
-            Ok((lambda, objf)) => (lambda, objf),
+            Ok((lambda, objf)) => (lambda.into(), objf),
             Err(err) => {
                 return Err(anyhow::anyhow!(
                     "Error in IPM during condensation: {:?}",
                     err
                 ));
             }
         };
-        self.w = self.lambda.clone();
+        self.w = self.lambda.clone().into();
         Ok(())
     }
 

src/algorithms/npod.rs

@@ -1,10 +1,11 @@
+use crate::routines::output::{cycles::CycleLog, cycles::NPCycle, NPResult};
+use crate::structs::weights::Weights;
 use crate::{
     algorithms::Status,
     prelude::{
         algorithms::Algorithms,
         routines::{
             evaluation::{ipm::burke, qr},
-            output::{CycleLog, NPCycle, NPResult},
             settings::Settings,
         },
     },
@@ -13,9 +14,9 @@ use crate::{
         theta::Theta,
     },
 };
+
 use anyhow::bail;
 use anyhow::Result;
-use faer::Col;
 use faer_ext::IntoNdarray;
 use pharmsol::{prelude::ErrorModel, ErrorModels};
 use pharmsol::{
@@ -37,8 +38,8 @@ pub struct NPOD<E: Equation> {
     equation: E,
     psi: Psi,
     theta: Theta,
-    lambda: Col<f64>,
-    w: Col<f64>,
+    lambda: Weights,
+    w: Weights,
     last_objf: f64,
     objf: f64,
     cycle: usize,
@@ -57,8 +58,8 @@ impl<E: Equation> Algorithms<E> for NPOD<E> {
             equation,
             psi: Psi::new(),
             theta: Theta::new(),
-            lambda: Col::zeros(0),
-            w: Col::zeros(0),
+            lambda: Weights::default(),
+            w: Weights::default(),
             last_objf: -1e30,
             objf: f64::NEG_INFINITY,
             cycle: 0,
@@ -157,15 +158,15 @@ impl<E: Equation> Algorithms<E> for NPOD<E> {
         }
 
         // Create state object
-        let state = NPCycle {
-            cycle: self.cycle,
-            objf: -2. * self.objf,
-            delta_objf: (self.last_objf - self.objf).abs(),
-            nspp: self.theta.nspp(),
-            theta: self.theta.clone(),
-            error_models: self.error_models.clone(),
-            status: self.status.clone(),
-        };
+        let state = NPCycle::new(
+            self.cycle,
+            -2. * self.objf,
+            self.error_models.clone(),
+            self.theta.clone(),
+            self.theta.nspp(),
+            (self.last_objf - self.objf).abs(),
+            self.status.clone(),
+        );
 
         // Write cycle log
         self.cycle_log.push(state);
@@ -205,11 +206,11 @@ impl<E: Equation> Algorithms<E> for NPOD<E> {
         let max_lambda = self
             .lambda
             .iter()
-            .fold(f64::NEG_INFINITY, |acc, &x| x.max(acc));
+            .fold(f64::NEG_INFINITY, |acc, x| x.max(acc));
 
         let mut keep = Vec::<usize>::new();
         for (index, lam) in self.lambda.iter().enumerate() {
-            if *lam > max_lambda / 1000_f64 {
+            if lam > max_lambda / 1000_f64 {
                 keep.push(index);
             }
         }
@@ -367,7 +368,7 @@ impl<E: Equation> 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().cloned().collect();
+        let w: Array1<f64> = self.w.clone().iter().collect();
         let pyl = psi.dot(&w);
 
         // Add new point to theta based on the optimization of the D function

src/algorithms/postprob.rs

@@ -4,18 +4,19 @@ use crate::{
     structs::{
         psi::{calculate_psi, Psi},
         theta::Theta,
+        weights::Weights,
     },
 };
 use anyhow::{Context, Result};
-use faer::Col;
+
 use pharmsol::prelude::{
     data::{Data, ErrorModels},
     simulator::Equation,
 };
 
 use crate::routines::evaluation::ipm::burke;
 use crate::routines::initialization;
-use crate::routines::output::{CycleLog, NPResult};
+use crate::routines::output::{cycles::CycleLog, NPResult};
 use crate::routines::settings::Settings;
 
 /// Posterior probability algorithm
@@ -24,7 +25,7 @@ pub struct POSTPROB<E: Equation> {
     equation: E,
     psi: Psi,
     theta: Theta,
-    w: Col<f64>,
+    w: Weights,
     objf: f64,
     cycle: usize,
     status: Status,
@@ -40,7 +41,7 @@ impl<E: Equation> Algorithms<E> for POSTPROB<E> {
             equation,
             psi: Psi::new(),
             theta: Theta::new(),
-            w: Col::zeros(0),
+            w: Weights::default(),
             objf: f64::INFINITY,
             cycle: 0,
             status: Status::Starting,

src/routines/evaluation/ipm.rs

@@ -1,4 +1,5 @@
 use crate::structs::psi::Psi;
+use crate::structs::weights::Weights;
 use anyhow::bail;
 use faer::linalg::triangular_solve::solve_lower_triangular_in_place;
 use faer::linalg::triangular_solve::solve_upper_triangular_in_place;
@@ -21,15 +22,15 @@ use rayon::prelude::*;
 ///
 /// # Returns
 ///
-/// On success, returns a tuple `(lam, obj)` where:
-///   - `lam` is a faer::Col<f64> containing the computed probability vector,
+/// On success, returns a tuple `(weights, obj)` where:
+///   - [Weights] contains the optimized weights (probabilities) for each support point.
 ///   - `obj` is the value of the objective function at the solution.
 ///
 /// # Errors
 ///
 /// This function returns an error if any step in the optimization (e.g. Cholesky factorization)
 /// fails.
-pub fn burke(psi: &Psi) -> anyhow::Result<(Col<f64>, f64)> {
+pub fn burke(psi: &Psi) -> anyhow::Result<(Weights, f64)> {
     let mut psi = psi.matrix().to_owned();
 
     // Ensure all entries are finite and make them non-negative.
@@ -274,7 +275,7 @@ pub fn burke(psi: &Psi) -> anyhow::Result<(Col<f64>, f64)> {
     let lam_sum: f64 = lam.iter().sum();
     lam = &lam / lam_sum;
 
-    Ok((lam, obj))
+    Ok((lam.into(), obj))
 }
 
 #[cfg(test)]
@@ -465,7 +466,11 @@ mod tests {
         // distribution depends on the optimization algorithm's convergence
 
         // Just verify that no single weight dominates excessively (basic sanity check)
-        let max_weight = lam.iter().cloned().fold(f64::NEG_INFINITY, f64::max);
+        let max_weight = lam
+            .weights()
+            .iter()
+            .cloned()
+            .fold(f64::NEG_INFINITY, f64::max);
         assert!(
             max_weight < 0.1,
             "No single weight should dominate in uniform matrix (max weight: {})",

src/routines/logger.rs

@@ -48,7 +48,7 @@ pub(crate) fn setup_log(settings: &mut Settings) -> Result<()> {
     let file_layer = match settings.log().write {
         true => {
             let layer = fmt::layer()
-                .with_writer(outputfile.file)
+                .with_writer(outputfile.file_owned())
                 .with_ansi(false)
                 .with_timer(timestamper.clone());