From bd2e7ac23cb46e60c5b40142a489da2666d263ec Mon Sep 17 00:00:00 2001 From: Markus Date: Mon, 14 Apr 2025 12:18:30 +0200 Subject: [PATCH 1/6] tests: Add tests for IPM --- Cargo.toml | 1 + src/routines/evaluation/ipm.rs | 180 +++++++++++++++++++++++++++++---- 2 files changed, 162 insertions(+), 19 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 356ac27e..4240cd67 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -40,6 +40,7 @@ codegen-units = 1 opt-level = 3 [dev-dependencies] +approx = "0.5.1" criterion = { version = "0.5" } diff --git a/src/routines/evaluation/ipm.rs b/src/routines/evaluation/ipm.rs index 98a633df..2b2a1153 100644 --- a/src/routines/evaluation/ipm.rs +++ b/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)> { - // 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)> { Ok((lam, obj)) } -// fn pprint(x: &Mat, 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::(), 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::(), 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::(), 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::(), 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::(), 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::(), 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()); + } +} From 1bafb4810eaaf7797674b5dd6264c41e916e8f2d Mon Sep 17 00:00:00 2001 From: Markus Date: Mon, 14 Apr 2025 16:21:41 +0200 Subject: [PATCH 2/6] Tests and docs for QR decomposition Renames calculate_R to qrd --- src/algorithms/npag.rs | 2 +- src/algorithms/npod.rs | 2 +- src/routines/evaluation/qr.rs | 84 ++++++++++++++++++++++++++++++++--- 3 files changed, 80 insertions(+), 8 deletions(-) diff --git a/src/algorithms/npag.rs b/src/algorithms/npag.rs index f2490bb9..85b94ddc 100644 --- a/src/algorithms/npag.rs +++ b/src/algorithms/npag.rs @@ -227,7 +227,7 @@ impl Algorithms for NPAG { 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::::new(); diff --git a/src/algorithms/npod.rs b/src/algorithms/npod.rs index ec81a775..c63bfe67 100644 --- a/src/algorithms/npod.rs +++ b/src/algorithms/npod.rs @@ -217,7 +217,7 @@ impl Algorithms for NPOD { 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::::new(); diff --git a/src/routines/evaluation/qr.rs b/src/routines/evaluation/qr.rs index ebb6b8ff..1613fb13 100644 --- a/src/routines/evaluation/qr.rs +++ b/src/routines/evaluation/qr.rs @@ -1,16 +1,30 @@ use crate::structs::psi::Psi; +use anyhow::{bail, Result}; use faer::linalg::solvers::ColPivQr; use faer::Mat; -pub fn calculate_r(psi: &Psi) -> (Mat, Vec) { - // 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]) and permutation vector (as [Vec]) +/// * Error if any row in the matrix sums to zero +pub fn qrd(psi: &Psi) -> Result<(Mat, Vec)> { + 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(); + if row_sum == 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 = mat.col_piv_qr(); @@ -20,5 +34,63 @@ pub fn calculate_r(psi: &Psi) -> (Mat, Vec) { // 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 = 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 = 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::::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); + } } From a633f54dc6c2e5dac4eb4b1b9b7a527585c8f247 Mon Sep 17 00:00:00 2001 From: Markus Date: Mon, 14 Apr 2025 16:27:02 +0200 Subject: [PATCH 3/6] Update ancient lib.rs --- src/lib.rs | 40 ++-------------------------------------- 1 file changed, 2 insertions(+), 38 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index 86656daa..1416e5ec 100644 --- a/src/lib.rs +++ b/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; From 07249fbb1af2659838b4b58549080914864b705c Mon Sep 17 00:00:00 2001 From: Markus Date: Mon, 14 Apr 2025 16:27:52 +0200 Subject: [PATCH 4/6] Remove comment --- src/algorithms/mod.rs | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/algorithms/mod.rs b/src/algorithms/mod.rs index b8e21b52..3b9e6eae 100644 --- a/src/algorithms/mod.rs +++ b/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; From dd4692ba890f5d63c9f836955e6a0ee158cffac3 Mon Sep 17 00:00:00 2001 From: Markus Hovd Date: Mon, 14 Apr 2025 16:31:30 +0200 Subject: [PATCH 5/6] Update src/routines/evaluation/qr.rs Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- src/routines/evaluation/qr.rs | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/routines/evaluation/qr.rs b/src/routines/evaluation/qr.rs index 1613fb13..e74732b1 100644 --- a/src/routines/evaluation/qr.rs +++ b/src/routines/evaluation/qr.rs @@ -20,7 +20,8 @@ pub fn qrd(psi: &Psi) -> Result<(Mat, Vec)> { // Normalize the rows to sum to 1 for (index, row) in mat.row_iter_mut().enumerate() { let row_sum: f64 = row.as_ref().iter().sum(); - if row_sum == 0.0 { + const EPSILON: f64 = 1e-10; + if row_sum.abs() < EPSILON { bail!("In psi, the row with index {} sums to zero", index); } row.iter_mut().for_each(|x| *x /= row_sum); From e97148124aab219e27336d998ce97f82976cf6b9 Mon Sep 17 00:00:00 2001 From: Markus Date: Mon, 14 Apr 2025 17:44:38 +0200 Subject: [PATCH 6/6] Fix comparison of zero Because psi can contain very small values, we can check for exactly zero --- src/routines/evaluation/qr.rs | 5 +++-- tests/onecomp.rs | 3 --- 2 files changed, 3 insertions(+), 5 deletions(-) diff --git a/src/routines/evaluation/qr.rs b/src/routines/evaluation/qr.rs index e74732b1..317fabf6 100644 --- a/src/routines/evaluation/qr.rs +++ b/src/routines/evaluation/qr.rs @@ -20,8 +20,9 @@ pub fn qrd(psi: &Psi) -> Result<(Mat, Vec)> { // Normalize the rows to sum to 1 for (index, row) in mat.row_iter_mut().enumerate() { let row_sum: f64 = row.as_ref().iter().sum(); - const EPSILON: f64 = 1e-10; - if row_sum.abs() < EPSILON { + + // 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); diff --git a/tests/onecomp.rs b/tests/onecomp.rs index 5047e8a7..739cc8c6 100644 --- a/tests/onecomp.rs +++ b/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)?;