tests: add property based tests

Author: schroedk

crates/approx-chol/tests/common/laplacian_prop.rs

@@ -1,4 +1,5 @@
 use proptest::prelude::*;
+use std::collections::VecDeque;
 
 pub type LaplacianCsr = (Vec<u32>, Vec<u32>, Vec<f64>, u32);
 
@@ -54,3 +55,95 @@ pub fn rhs_for_dimension(n: usize) -> Vec<f64> {
     }
     rhs
 }
+
+#[allow(dead_code)]
+pub fn csr_matvec(row_ptrs: &[u32], col_indices: &[u32], values: &[f64], x: &[f64]) -> Vec<f64> {
+    let n = row_ptrs.len() - 1;
+    let mut y = vec![0.0; n];
+    for i in 0..n {
+        let start = row_ptrs[i] as usize;
+        let end = row_ptrs[i + 1] as usize;
+        for k in start..end {
+            y[i] += values[k] * x[col_indices[k] as usize];
+        }
+    }
+    y
+}
+
+#[allow(dead_code)]
+pub fn is_connected(row_ptrs: &[u32], col_indices: &[u32], n: u32) -> bool {
+    if n <= 1 {
+        return true;
+    }
+    let n = n as usize;
+    let mut visited = vec![false; n];
+    let mut queue = VecDeque::new();
+    visited[0] = true;
+    queue.push_back(0usize);
+    while let Some(v) = queue.pop_front() {
+        let start = row_ptrs[v] as usize;
+        let end = row_ptrs[v + 1] as usize;
+        for &col in &col_indices[start..end] {
+            let u = col as usize;
+            if !visited[u] {
+                visited[u] = true;
+                queue.push_back(u);
+            }
+        }
+    }
+    visited.iter().all(|&v| v)
+}
+
+#[allow(dead_code)]
+pub fn norm2(v: &[f64]) -> f64 {
+    v.iter().map(|x| x * x).sum::<f64>().sqrt()
+}
+
+#[allow(dead_code)]
+pub fn random_zero_sum_rhs_strategy(n: usize) -> BoxedStrategy<Vec<f64>> {
+    if n <= 1 {
+        Just(vec![0.0; n]).boxed()
+    } else {
+        prop::collection::vec(-10.0f64..10.0, n)
+            .prop_map(|mut v| {
+                let mean = v.iter().sum::<f64>() / v.len() as f64;
+                for x in &mut v {
+                    *x -= mean;
+                }
+                v
+            })
+            .boxed()
+    }
+}
+
+#[allow(dead_code)]
+pub fn laplacian_with_rhs_strategy() -> impl Strategy<Value = (LaplacianCsr, Vec<f64>)> {
+    laplacian_csr_strategy().prop_flat_map(|(rp, ci, vals, n)| {
+        random_zero_sum_rhs_strategy(n as usize)
+            .prop_map(move |rhs| ((rp.clone(), ci.clone(), vals.clone(), n), rhs))
+    })
+}
+
+#[allow(dead_code)]
+pub fn sddm_csr_strategy() -> impl Strategy<Value = LaplacianCsr> {
+    (1usize..=8).prop_flat_map(|n| {
+        let pair_count = n * (n - 1) / 2;
+        (
+            prop::collection::vec(0u8..=4, pair_count),
+            prop::collection::vec(1u8..=5, n),
+        )
+            .prop_map(move |(edge_weights, surpluses)| {
+                let (rp, ci, mut vals, n_u32) = build_laplacian_csr(n, &edge_weights);
+                for i in 0..n {
+                    let start = rp[i] as usize;
+                    let end = rp[i + 1] as usize;
+                    for k in start..end {
+                        if ci[k] as usize == i {
+                            vals[k] += surpluses[i] as f64;
+                        }
+                    }
+                }
+                (rp, ci, vals, n_u32)
+            })
+    })
+}

crates/approx-chol/tests/property_factorization.rs

@@ -5,11 +5,18 @@ use panic_ok::OrPanic;
 mod laplacian_prop;
 
 use approx_chol::{factorize, factorize_with, Config, CsrRef};
-use laplacian_prop::{laplacian_csr_strategy, rhs_for_dimension};
+use laplacian_prop::{
+    csr_matvec, is_connected, laplacian_csr_strategy, laplacian_with_rhs_strategy, norm2,
+    rhs_for_dimension, sddm_csr_strategy,
+};
 use proptest::prelude::*;
 use std::panic::{catch_unwind, AssertUnwindSafe};
 
 proptest! {
+    // -----------------------------------------------------------------------
+    // Existing: panic-freedom and finiteness
+    // -----------------------------------------------------------------------
+
     #[test]
     fn default_factorization_solve_is_panic_free_and_finite(
         (row_ptrs, col_indices, values, n) in laplacian_csr_strategy()
@@ -79,4 +86,210 @@ proptest! {
         let work = run.or_panic("checked above");
         prop_assert!(work.iter().all(|x| x.is_finite()));
     }
+
+    // -----------------------------------------------------------------------
+    // Solution quality: residual ||Ax - b|| / ||b|| is bounded
+    // -----------------------------------------------------------------------
+
+    #[test]
+    fn ac_residual_is_bounded(
+        (row_ptrs, col_indices, values, n) in laplacian_csr_strategy()
+    ) {
+        prop_assume!(n >= 2);
+        prop_assume!(is_connected(&row_ptrs, &col_indices, n));
+
+        let csr = CsrRef::new(&row_ptrs, &col_indices, &values, n)
+            .or_panic("valid CSR");
+        let factor = factorize(csr).or_panic("factorization");
+        let rhs = rhs_for_dimension(n as usize);
+        let x = factor.solve(&rhs).or_panic("solve");
+
+        let ax = csr_matvec(&row_ptrs, &col_indices, &values, &x);
+        let residual: Vec<f64> = ax.iter().zip(rhs.iter()).map(|(a, b)| a - b).collect();
+        let r_norm = norm2(&residual);
+        let b_norm = norm2(&rhs);
+
+        prop_assert!(
+            b_norm < 1e-15 || r_norm / b_norm < 100.0,
+            "AC relative residual too large: {:.4e} (r_norm={:.4e}, b_norm={:.4e})",
+            r_norm / b_norm, r_norm, b_norm
+        );
+    }
+
+    #[test]
+    fn ac2_residual_is_bounded(
+        (row_ptrs, col_indices, values, n) in laplacian_csr_strategy()
+    ) {
+        prop_assume!(n >= 2);
+        prop_assume!(is_connected(&row_ptrs, &col_indices, n));
+
+        let csr = CsrRef::new(&row_ptrs, &col_indices, &values, n)
+            .or_panic("valid CSR");
+        let factor = factorize_with(
+            csr,
+            Config {
+                seed: 7,
+                split_merge: Some(2),
+            },
+        )
+        .or_panic("AC2 factorization");
+        let rhs = rhs_for_dimension(n as usize);
+        let x = factor.solve(&rhs).or_panic("solve");
+
+        let ax = csr_matvec(&row_ptrs, &col_indices, &values, &x);
+        let residual: Vec<f64> = ax.iter().zip(rhs.iter()).map(|(a, b)| a - b).collect();
+        let r_norm = norm2(&residual);
+        let b_norm = norm2(&rhs);
+
+        prop_assert!(
+            b_norm < 1e-15 || r_norm / b_norm < 100.0,
+            "AC2 relative residual too large: {:.4e} (r_norm={:.4e}, b_norm={:.4e})",
+            r_norm / b_norm, r_norm, b_norm
+        );
+    }
+
+    #[test]
+    fn random_rhs_residual_is_bounded(
+        ((row_ptrs, col_indices, values, n), rhs) in laplacian_with_rhs_strategy()
+    ) {
+        prop_assume!(n >= 2);
+        prop_assume!(is_connected(&row_ptrs, &col_indices, n));
+        let b_norm = norm2(&rhs);
+        prop_assume!(b_norm > 1e-10);
+
+        let csr = CsrRef::new(&row_ptrs, &col_indices, &values, n)
+            .or_panic("valid CSR");
+        let factor = factorize(csr).or_panic("factorization");
+        let x = factor.solve(&rhs).or_panic("solve");
+
+        let ax = csr_matvec(&row_ptrs, &col_indices, &values, &x);
+        let residual: Vec<f64> = ax.iter().zip(rhs.iter()).map(|(a, b)| a - b).collect();
+        let r_norm = norm2(&residual);
+
+        prop_assert!(
+            r_norm / b_norm < 100.0,
+            "random-RHS relative residual too large: {:.4e}", r_norm / b_norm
+        );
+    }
+
+    // -----------------------------------------------------------------------
+    // SDDM matrices (Gremban augmentation path)
+    // -----------------------------------------------------------------------
+
+    #[test]
+    fn sddm_factorization_is_panic_free_and_finite(
+        (row_ptrs, col_indices, values, n) in sddm_csr_strategy()
+    ) {
+        let run = catch_unwind(AssertUnwindSafe(|| {
+            let csr = CsrRef::new(&row_ptrs, &col_indices, &values, n)
+                .or_panic("valid SDDM CSR");
+            let factor = factorize(csr).or_panic("factorization");
+            let mut rhs = vec![0.0_f64; n as usize];
+            if n >= 2 {
+                rhs[0] = 1.0;
+                rhs[(n as usize) - 1] = -1.0;
+            }
+            factor.solve(&rhs).or_panic("solve")
+        }));
+
+        prop_assert!(run.is_ok(), "SDDM factorization or solve panicked");
+        let x = run.or_panic("checked above");
+        prop_assert!(x.iter().all(|v| v.is_finite()), "SDDM solution has non-finite values");
+    }
+
+    // -----------------------------------------------------------------------
+    // Determinism: same seed + same input → identical output
+    // -----------------------------------------------------------------------
+
+    #[test]
+    fn deterministic_with_fixed_seed(
+        (row_ptrs, col_indices, values, n) in laplacian_csr_strategy()
+    ) {
+        let config = Config { seed: 42, ..Default::default() };
+        let rhs = rhs_for_dimension(n as usize);
+
+        let csr1 = CsrRef::new(&row_ptrs, &col_indices, &values, n)
+            .or_panic("valid CSR");
+        let x1 = factorize_with(csr1, config).or_panic("factorize 1")
+            .solve(&rhs).or_panic("solve 1");
+
+        let csr2 = CsrRef::new(&row_ptrs, &col_indices, &values, n)
+            .or_panic("valid CSR");
+        let x2 = factorize_with(csr2, config).or_panic("factorize 2")
+            .solve(&rhs).or_panic("solve 2");
+
+        prop_assert_eq!(x1.len(), x2.len());
+        for (a, b) in x1.iter().zip(x2.iter()) {
+            prop_assert!(
+                a.to_bits() == b.to_bits(),
+                "non-deterministic: {} vs {}", a, b
+            );
+        }
+    }
+
+    // -----------------------------------------------------------------------
+    // Factor dimensions are consistent
+    // -----------------------------------------------------------------------
+
+    #[test]
+    fn factor_dimensions_are_consistent(
+        (row_ptrs, col_indices, values, n) in laplacian_csr_strategy()
+    ) {
+        let csr = CsrRef::new(&row_ptrs, &col_indices, &values, n)
+            .or_panic("valid CSR");
+        let factor = factorize(csr).or_panic("factorization");
+
+        prop_assert_eq!(
+            factor.original_n(), n as usize,
+            "original_n must match input dimension"
+        );
+        prop_assert!(
+            factor.n() >= n as usize,
+            "factor.n() must be >= input dimension"
+        );
+        // Pure Laplacians should not be augmented
+        prop_assert_eq!(
+            factor.n(), n as usize,
+            "pure Laplacian should not trigger Gremban augmentation"
+        );
+    }
+
+    #[test]
+    fn sddm_factor_dimensions_are_consistent(
+        (row_ptrs, col_indices, values, n) in sddm_csr_strategy()
+    ) {
+        let csr = CsrRef::new(&row_ptrs, &col_indices, &values, n)
+            .or_panic("valid SDDM CSR");
+        let factor = factorize(csr).or_panic("factorization");
+
+        prop_assert_eq!(
+            factor.original_n(), n as usize,
+            "original_n must match input dimension"
+        );
+        prop_assert!(
+            factor.n() > n as usize,
+            "SDDM should trigger Gremban augmentation (factor.n() must be > n)"
+        );
+    }
+
+    // -----------------------------------------------------------------------
+    // f32 support
+    // -----------------------------------------------------------------------
+
+    #[test]
+    fn f32_factorization_is_finite(
+        (row_ptrs, col_indices, values_f64, n) in laplacian_csr_strategy()
+    ) {
+        let values_f32: Vec<f32> = values_f64.iter().map(|&v| v as f32).collect();
+        let csr = CsrRef::new(&row_ptrs, &col_indices, &values_f32, n)
+            .or_panic("valid f32 CSR");
+        let factor = factorize(csr).or_panic("f32 factorization");
+        let mut rhs = vec![0.0_f32; n as usize];
+        if n >= 2 {
+            rhs[0] = 1.0;
+            rhs[(n as usize) - 1] = -1.0;
+        }
+        let x = factor.solve(&rhs).or_panic("f32 solve");
+        prop_assert!(x.iter().all(|v| v.is_finite()), "f32 solution has non-finite values");
+    }
 }