chore: Improve test coverage

Cargo.toml

@@ -40,6 +40,7 @@ codegen-units = 1
 opt-level = 3
 
 [dev-dependencies]
+approx = "0.5.1"
 criterion = { version = "0.5" }
 
 

src/algorithms/mod.rs

@@ -17,8 +17,6 @@ use pharmsol::{ErrorModel, Predictions, Subject};
 use postprob::POSTPROB;
 use serde::{Deserialize, Serialize};
 
-// use self::{data::Subject, simulator::Equation};
-
 pub mod npag;
 pub mod npod;
 pub mod postprob;

src/algorithms/npag.rs

@@ -227,7 +227,7 @@ impl<E: Equation> Algorithms<E> for NPAG<E> {
         self.psi.filter_column_indices(keep.as_slice());
 
         //Rank-Revealing Factorization
-        let (r, perm) = qr::calculate_r(&self.psi);
+        let (r, perm) = qr::qrd(&self.psi)?;
 
         let mut keep = Vec::<usize>::new();
 

src/algorithms/npod.rs

@@ -217,7 +217,7 @@ impl<E: Equation> Algorithms<E> for NPOD<E> {
         self.psi.filter_column_indices(keep.as_slice());
 
         //Rank-Revealing Factorization
-        let (r, perm) = qr::calculate_r(&self.psi);
+        let (r, perm) = qr::qrd(&self.psi)?;
 
         let mut keep = Vec::<usize>::new();
 

src/lib.rs

@@ -4,48 +4,12 @@
 //!
 //! # Configuration
 //!
-//! PMcore is configured using a TOML file, which specifies the settings for the algorithm. The settings file is divided into sections, each of which specifies a different aspect of the algorithm. They are further described in the [routines::settings] module.
-//!
-//! # Model definition
-//!
-//! As PMcore is provided as a library, the user must define the model to be used. Some algebraic models are provided, and more will be added, but the user is free to define their own model. The model must implement the [routines::simulation::predict] trait, which specifies the methods that the model must implement. For more information on how to define models with ordinary differential equations, please look at at the examples.
+//! PMcore is configured using [routines::settings::Settings], which specifies the settings for the program.
 //!
 //! # Data format
 //!
-//! Data is provided in a CSV file, and the format is described in the table below. For each subject, there must be at least one dose event.
-//!
-//! | Column | Description                                                         | Conditions                       |
-//! |--------|---------------------------------------------------------------------|----------------------------------|
-//! | `ID`   | Unique subject ID                                                   |                                  |
-//! | `EVID` | Event type; 0 = observation, 1 = dose, 4 = reset                    |                                  |
-//! | `TIME` | Time of event                                                       |                                  |
-//! | `DUR`  | Duration of an infusion                                             | Must be provided if EVID = 1     |
-//! | `DOSE` | The dose amount                                                     | Must be provided if EVID = 1     |
-//! | `ADDL` | The number of additional doses to be given at the interval `II`     |                                  |
-//! | `II`   | Interval for additional doses                                       |                                  |
-//! | `INPUT`| The compartment the dose should be delivered to                     |                                  |
-//! | `OUT`  | The observed value                                                  | Must be provided if EVID = 0     |
-//! | `OUTEQ`| Denotes the output equation for which `OUT` is provided             |                                  |
-//! | `C0`   | Optional override of the error polynomial coefficient               |                                  |
-//! | `C1`   | Optional override of the error polynomial coefficient               |                                  |
-//! | `C2`   | Optional override of the error polynomial coefficient               |                                  |
-//! | `C3`   | Optional override of the error polynomial coefficient               |                                  |
-//! | `COV...`| Any additional columns are assumed to be covariates, one per column | Must be present for the first event for each subject |
-//!
-//! # Examples
-//!
-//! A couple of examples are provided in the `examples` directory. The `settings.toml` file contains the settings for the algorithm, and the `data.csv` file contains the data.
-//!
-//! They can be run using the following command
-//!
-//! ```sh
-//! cargo run --release --example `example_name`
-//! ```
-//!
-//! where `example_name` is the name of the example to run. Currently, the following examples are available:
+//! PMcore is heavily linked to [pharmsol], which provides the data structures and routines for handling pharmacokinetic data. The data is stored in a [pharmsol::Data] structure, and can either be read from a CSV file, using [pharmsol::data::parse_pmetrics::read_pmetrics], or created dynamically using the [pharmsol::data::builder::SubjectBuilder].
 //!
-//! - `bimodal_ke`: A simple, one-compartmental model following an intravenous infusion. The example is named by the bimodal distribution of one of two parameters, `Ke`, the elimination rate constant. The example is designed to demonstrate the ability of the algorithm to handle multimodal distributions, and detect outliers.
-//! - `simple_covariates`: A simple example with a single subject and a single dose event, with covariates.
 
 /// Provides the various algorithms used within the framework
 // pub mod algorithms;

src/routines/evaluation/ipm.rs

@@ -30,23 +30,20 @@ use rayon::prelude::*;
 /// 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)> {
-    // Get the underlying matrix. (Assume psi.matrix() returns an ndarray-compatible matrix.)
     let mut psi = psi.matrix().to_owned();
 
     // Ensure all entries are finite and make them non-negative.
-    psi.row_iter_mut()
-        .try_for_each(|row| {
-            row.iter_mut().try_for_each(|x| {
-                if !x.is_finite() {
-                    bail!("Input matrix must have finite entries")
-                } else {
-                    // Coerce negatives to non-negative (could alternatively return an error)
-                    *x = x.abs();
-                    Ok(())
-                }
-            })
+    psi.row_iter_mut().try_for_each(|row| {
+        row.iter_mut().try_for_each(|x| {
+            if !x.is_finite() {
+                bail!("Input matrix must have finite entries")
+            } else {
+                // Coerce negatives to non-negative (could alternatively return an error)
+                *x = x.abs();
+                Ok(())
+            }
         })
-        .unwrap();
+    })?;
 
     // Let psi be of shape (n_sub, n_point)
     let (n_sub, n_point) = psi.shape();
@@ -280,9 +277,154 @@ pub fn burke(psi: &Psi) -> anyhow::Result<(Col<f64>, f64)> {
     Ok((lam, obj))
 }
 
-// fn pprint(x: &Mat<f64>, name: &str) {
-//     println!("Matrix: {}", name);
-//     x.row_iter().for_each(|row| {
-//         println!("{:.unwrap()}", row);
-//     });
-// }
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use approx::assert_relative_eq;
+    use faer::Mat;
+
+    #[test]
+    fn test_burke_identity() {
+        // Test with a small identity matrix
+        // For an identity matrix, each support point should have equal weight
+        let n = 100;
+        let mat = Mat::identity(n, n);
+        let psi = Psi::from(mat);
+
+        let (lam, _) = burke(&psi).unwrap();
+
+        // For identity matrix, all lambda values should be equal
+        let expected = 1.0 / n as f64;
+        for i in 0..n {
+            assert_relative_eq!(lam[i], expected, epsilon = 1e-10);
+        }
+
+        // Check that lambda sums to 1
+        assert_relative_eq!(lam.iter().sum::<f64>(), 1.0, epsilon = 1e-10);
+    }
+
+    #[test]
+    fn test_burke_uniform_square() {
+        // Test with a matrix of all ones
+        // This should also result in uniform weights
+        let n_sub = 10;
+        let n_point = 10;
+        let mat = Mat::from_fn(n_sub, n_point, |_, _| 1.0);
+        let psi = Psi::from(mat);
+
+        let (lam, _) = burke(&psi).unwrap();
+
+        // Check that lambda sums to 1
+        assert_relative_eq!(lam.iter().sum::<f64>(), 1.0, epsilon = 1e-10);
+
+        // For uniform matrix, all lambda values should be equal
+        let expected = 1.0 / n_point as f64;
+        for i in 0..n_point {
+            assert_relative_eq!(lam[i], expected, epsilon = 1e-10);
+        }
+    }
+
+    #[test]
+    fn test_burke_uniform_wide() {
+        // Test with a matrix of all ones
+        // This should also result in uniform weights
+        let n_sub = 10;
+        let n_point = 100;
+        let mat = Mat::from_fn(n_sub, n_point, |_, _| 1.0);
+        let psi = Psi::from(mat);
+
+        let (lam, _) = burke(&psi).unwrap();
+
+        // Check that lambda sums to 1
+        assert_relative_eq!(lam.iter().sum::<f64>(), 1.0, epsilon = 1e-10);
+
+        // For uniform matrix, all lambda values should be equal
+        let expected = 1.0 / n_point as f64;
+        for i in 0..n_point {
+            assert_relative_eq!(lam[i], expected, epsilon = 1e-10);
+        }
+    }
+
+    #[test]
+    fn test_burke_uniform_long() {
+        // Test with a matrix of all ones
+        // This should also result in uniform weights
+        let n_sub = 100;
+        let n_point = 10;
+        let mat = Mat::from_fn(n_sub, n_point, |_, _| 1.0);
+        let psi = Psi::from(mat);
+
+        let (lam, _) = burke(&psi).unwrap();
+
+        // Check that lambda sums to 1
+        assert_relative_eq!(lam.iter().sum::<f64>(), 1.0, epsilon = 1e-10);
+
+        // For uniform matrix, all lambda values should be equal
+        let expected = 1.0 / n_point as f64;
+        for i in 0..n_point {
+            assert_relative_eq!(lam[i], expected, epsilon = 1e-10);
+        }
+    }
+
+    #[test]
+    fn test_burke_with_non_uniform_matrix() {
+        // Test with a non-uniform matrix
+        // Create a matrix where one column is clearly better
+        let n_sub = 3;
+        let n_point = 4;
+        let mat = Mat::from_fn(n_sub, n_point, |_, j| if j == 0 { 10.0 } else { 1.0 });
+        let psi = Psi::from(mat);
+
+        let (lam, _) = burke(&psi).unwrap();
+
+        // Check that lambda sums to 1
+        assert_relative_eq!(lam.iter().sum::<f64>(), 1.0, epsilon = 1e-10);
+
+        // First support point should have highest weight
+        assert!(lam[0] > lam[1]);
+        assert!(lam[0] > lam[2]);
+        assert!(lam[0] > lam[3]);
+    }
+
+    #[test]
+    fn test_burke_with_negative_values() {
+        // The algorithm should handle negative values by taking their absolute value
+        let n_sub = 2;
+        let n_point = 3;
+        let mat = Mat::from_fn(
+            n_sub,
+            n_point,
+            |i, j| if i == 0 && j == 0 { -5.0 } else { 1.0 },
+        );
+        let psi = Psi::from(mat);
+
+        let result = burke(&psi);
+        assert!(result.is_ok());
+
+        let (lam, _) = result.unwrap();
+        // Check that lambda sums to 1
+        assert_relative_eq!(lam.iter().sum::<f64>(), 1.0, epsilon = 1e-10);
+
+        // First support point should have highest weight due to the high absolute value
+        assert!(lam[0] > lam[1]);
+        assert!(lam[0] > lam[2]);
+    }
+
+    #[test]
+    fn test_burke_with_non_finite_values() {
+        // The algorithm should return an error for non-finite values
+        let n_sub = 10;
+        let n_point = 10;
+        let mat = Mat::from_fn(n_sub, n_point, |i, j| {
+            if i == 0 && j == 0 {
+                f64::NAN
+            } else {
+                1.0
+            }
+        });
+        let psi = Psi::from(mat);
+
+        let result = burke(&psi);
+        assert!(result.is_err());
+    }
+}

src/routines/evaluation/qr.rs

@@ -1,16 +1,32 @@
 use crate::structs::psi::Psi;
+use anyhow::{bail, Result};
 use faer::linalg::solvers::ColPivQr;
 use faer::Mat;
 
-pub fn calculate_r(psi: &Psi) -> (Mat<f64>, Vec<usize>) {
-    // Clone the matrix, as we will modify it
-    let mut mat = psi.matrix().clone();
+/// Perform a QR decomposition on the Psi matrix
+///
+/// Normalizes each row of the matrix to sum to 1 before decomposition.
+/// Returns the R matrix from QR decomposition and the column permutation vector.
+///
+/// # Arguments
+/// * `psi` - The Psi matrix to decompose
+///
+/// # Returns
+/// * Tuple containing the R matrix (as [faer::Mat<f64>]) and permutation vector (as [Vec<usize>])
+/// * Error if any row in the matrix sums to zero
+pub fn qrd(psi: &Psi) -> Result<(Mat<f64>, Vec<usize>)> {
+    let mut mat = psi.matrix().to_owned();
 
     // Normalize the rows to sum to 1
-    mat.row_iter_mut().for_each(|row| {
+    for (index, row) in mat.row_iter_mut().enumerate() {
         let row_sum: f64 = row.as_ref().iter().sum();
+
+        // Check if the row sum is zero
+        if row_sum.abs() == 0.0 {
+            bail!("In psi, the row with index {} sums to zero", index);
+        }
         row.iter_mut().for_each(|x| *x /= row_sum);
-    });
+    }
 
     // Perform column pivoted QR decomposition
     let qr: ColPivQr<f64> = mat.col_piv_qr();
@@ -20,5 +36,63 @@ pub fn calculate_r(psi: &Psi) -> (Mat<f64>, Vec<usize>) {
 
     // Get the permutation information
     let perm = qr.P().arrays().0.to_vec();
-    (r_mat, perm)
+    Ok((r_mat, perm))
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn test_identity() {
+        // Create a 2x2 identity matrix
+        let mat: Mat<f64> = Mat::identity(10, 10);
+        let psi = Psi::from(mat);
+
+        // Perform the QR decomposition
+        let (r_mat, perm) = qrd(&psi).unwrap();
+
+        // Check that R is an identity matrix
+        let expected_r_mat: Mat<f64> = Mat::identity(10, 10);
+        assert_eq!(r_mat, expected_r_mat);
+
+        // Check that the permutation is the identity
+        assert_eq!(perm, vec![0, 1, 2, 3, 4, 5, 6, 7, 8, 9]);
+    }
+
+    #[test]
+    fn test_with_zero_row_sum() {
+        // Create a test matrix with a row that sums to zero
+        let mat = Mat::from_fn(2, 2, |i, j| {
+            match (i, j) {
+                (0, 0) => 1.0,
+                (0, 1) => 2.0,
+                (1, 0) => 0.0, // Row that sums to zero
+                (1, 1) => 0.0,
+                _ => 0.0,
+            }
+        });
+        let psi = Psi::from(mat);
+
+        // Perform the QR decomposition
+        let result = qrd(&psi);
+
+        // Confirm that the function returns an error
+        assert!(result.is_err(), "Expected an error due to zero row sum");
+    }
+
+    #[test]
+    fn test_empty_matrix() {
+        // Create an empty Psi
+        let mat = Mat::<f64>::new();
+        let psi = Psi::from(mat);
+
+        // Should not panic
+        let (r_mat, perm) = qrd(&psi).unwrap();
+
+        // Empty matrix should produce empty results
+        assert_eq!(r_mat.nrows(), 0);
+        assert_eq!(r_mat.ncols(), 0);
+        assert_eq!(perm.len(), 0);
+    }
 }

tests/onecomp.rs

@@ -55,9 +55,6 @@ fn test_one_compartment() -> Result<()> {
     });
 
     let data = data::Data::new(subjects);
-    data.get_subjects().iter().for_each(|subject| {
-        println!("{}", subject);
-    });
 
     // Run the algorithm
     let mut algorithm = dispatch_algorithm(settings, eq, data)?;