diff --git a/.gitignore b/.gitignore index b5a16c6..fdf592a 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,6 @@ /target /.cspellcache -/coverage/* +/coverage +/build* +/cobertura.xml +.DS_Store diff --git a/Cargo.lock b/Cargo.lock index ed2107b..ed0fe09 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -47,6 +47,27 @@ version = "1.5.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "c08606f8c3cbf4ce6ec8e28fb0014a2c086708fe954eaa885384a6165172e7e8" +[[package]] +name = "bit-set" +version = "0.8.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "08807e080ed7f9d5433fa9b275196cfc35414f66a0c79d864dc51a0d825231a3" +dependencies = [ + "bit-vec", +] + +[[package]] +name = "bit-vec" +version = "0.8.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5e764a1d40d510daf35e07be9eb06e75770908c27d411ee6c92109c9840eaaf7" + +[[package]] +name = "bitflags" +version = "2.10.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "812e12b5285cc515a9c72a5c1d3b6d46a19dac5acfef5265968c166106e31dd3" + [[package]] name = "bumpalo" version = "3.19.0" @@ -205,12 +226,46 @@ version = "1.15.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "48c757948c5ede0e46177b7add2e67155f70e33c07fea8284df6576da70b3719" +[[package]] +name = "errno" +version = "0.3.14" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "39cab71617ae0d63f51a36d69f866391735b51691dbda63cf6f96d042b63efeb" +dependencies = [ + "libc", + "windows-sys", +] + +[[package]] +name = "fastrand" +version = "2.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "37909eebbb50d72f9059c3b6d82c0463f2ff062c9e95845c43a6c9c0355411be" + [[package]] name = "find-msvc-tools" version = "0.1.5" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "3a3076410a55c90011c298b04d0cfa770b00fa04e1e3c97d3f6c9de105a03844" +[[package]] +name = "fnv" +version = "1.0.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3f9eec918d3f24069decb9af1554cad7c880e2da24a9afd88aca000531ab82c1" + +[[package]] +name = "getrandom" +version = "0.3.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "899def5c37c4fd7b2664648c28120ecec138e4d395b459e5ca34f9cce2dd77fd" +dependencies = [ + "cfg-if", + "libc", + "r-efi", + "wasip2", +] + [[package]] name = "glam" version = "0.14.0" @@ -350,6 +405,8 @@ dependencies = [ "approx", "criterion", "nalgebra", + "pastey", + "proptest", ] [[package]] @@ -358,6 +415,12 @@ version = "0.2.178" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "37c93d8daa9d8a012fd8ab92f088405fb202ea0b6ab73ee2482ae66af4f42091" +[[package]] +name = "linux-raw-sys" +version = "0.11.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "df1d3c3b53da64cf5760482273a98e575c651a67eec7f77df96b5b642de8f039" + [[package]] name = "matrixmultiply" version = "0.3.10" @@ -493,6 +556,12 @@ version = "1.0.15" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "57c0d7b74b563b49d38dae00a0c37d4d6de9b432382b2892f0574ddcae73fd0a" +[[package]] +name = "pastey" +version = "0.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "57d6c094ee800037dff99e02cab0eaf3142826586742a270ab3d7a62656bd27a" + [[package]] name = "plotters" version = "0.3.7" @@ -521,6 +590,15 @@ dependencies = [ "plotters-backend", ] +[[package]] +name = "ppv-lite86" +version = "0.2.21" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "85eae3c4ed2f50dcfe72643da4befc30deadb458a9b590d720cde2f2b1e97da9" +dependencies = [ + "zerocopy", +] + [[package]] name = "proc-macro2" version = "1.0.103" @@ -530,6 +608,31 @@ dependencies = [ "unicode-ident", ] +[[package]] +name = "proptest" +version = "1.9.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bee689443a2bd0a16ab0348b52ee43e3b2d1b1f931c8aa5c9f8de4c86fbe8c40" +dependencies = [ + "bit-set", + "bit-vec", + "bitflags", + "num-traits", + "rand", + "rand_chacha", + "rand_xorshift", + "regex-syntax", + "rusty-fork", + "tempfile", + "unarray", +] + +[[package]] +name = "quick-error" +version = "1.2.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a1d01941d82fa2ab50be1e79e6714289dd7cde78eba4c074bc5a4374f650dfe0" + [[package]] name = "quote" version = "1.0.42" @@ -539,6 +642,50 @@ dependencies = [ "proc-macro2", ] +[[package]] +name = "r-efi" +version = "5.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "69cdb34c158ceb288df11e18b4bd39de994f6657d83847bdffdbd7f346754b0f" + +[[package]] +name = "rand" +version = "0.9.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6db2770f06117d490610c7488547d543617b21bfa07796d7a12f6f1bd53850d1" +dependencies = [ + "rand_chacha", + "rand_core", +] + +[[package]] +name = "rand_chacha" +version = "0.9.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d3022b5f1df60f26e1ffddd6c66e8aa15de382ae63b3a0c1bfc0e4d3e3f325cb" +dependencies = [ + "ppv-lite86", + "rand_core", +] + +[[package]] +name = "rand_core" +version = "0.9.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "99d9a13982dcf210057a8a78572b2217b667c3beacbf3a0d8b454f6f82837d38" +dependencies = [ + "getrandom", +] + +[[package]] +name = "rand_xorshift" +version = "0.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "513962919efc330f829edb2535844d1b912b0fbe2ca165d613e4e8788bb05a5a" +dependencies = [ + "rand_core", +] + [[package]] name = "rawpointer" version = "0.2.1" @@ -594,12 +741,37 @@ version = "0.8.8" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "7a2d987857b319362043e95f5353c0535c1f58eec5336fdfcf626430af7def58" +[[package]] +name = "rustix" +version = "1.1.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cd15f8a2c5551a84d56efdc1cd049089e409ac19a3072d5037a17fd70719ff3e" +dependencies = [ + "bitflags", + "errno", + "libc", + "linux-raw-sys", + "windows-sys", +] + [[package]] name = "rustversion" version = "1.0.22" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "b39cdef0fa800fc44525c84ccb54a029961a8215f9619753635a9c0d2538d46d" +[[package]] +name = "rusty-fork" +version = "0.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cc6bf79ff24e648f6da1f8d1f011e9cac26491b619e6b9280f2b47f1774e6ee2" +dependencies = [ + "fnv", + "quick-error", + "tempfile", + "wait-timeout", +] + [[package]] name = "ryu" version = "1.0.20" @@ -697,6 +869,19 @@ dependencies = [ "unicode-ident", ] +[[package]] +name = "tempfile" +version = "3.23.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2d31c77bdf42a745371d260a26ca7163f1e0924b64afa0b688e61b5a9fa02f16" +dependencies = [ + "fastrand", + "getrandom", + "once_cell", + "rustix", + "windows-sys", +] + [[package]] name = "tinytemplate" version = "1.2.1" @@ -713,12 +898,27 @@ version = "1.19.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "562d481066bde0658276a35467c4af00bdc6ee726305698a55b86e61d7ad82bb" +[[package]] +name = "unarray" +version = "0.1.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "eaea85b334db583fe3274d12b4cd1880032beab409c0d774be044d4480ab9a94" + [[package]] name = "unicode-ident" version = "1.0.22" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "9312f7c4f6ff9069b165498234ce8be658059c6728633667c526e27dc2cf1df5" +[[package]] +name = "wait-timeout" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "09ac3b126d3914f9849036f826e054cbabdc8519970b8998ddaf3b5bd3c65f11" +dependencies = [ + "libc", +] + [[package]] name = "walkdir" version = "2.5.0" @@ -729,6 +929,15 @@ dependencies = [ "winapi-util", ] +[[package]] +name = "wasip2" +version = "1.0.1+wasi-0.2.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0562428422c63773dad2c345a1882263bbf4d65cf3f42e90921f787ef5ad58e7" +dependencies = [ + "wit-bindgen", +] + [[package]] name = "wasm-bindgen" version = "0.2.106" @@ -840,6 +1049,12 @@ dependencies = [ "windows-link", ] +[[package]] +name = "wit-bindgen" +version = "0.46.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f17a85883d4e6d00e8a97c586de764dabcc06133f7f1d55dce5cdc070ad7fe59" + [[package]] name = "zerocopy" version = "0.8.31" diff --git a/Cargo.toml b/Cargo.toml index de2a357..6daaf94 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -5,6 +5,7 @@ edition = "2024" rust-version = "1.92" license = "BSD-3-Clause" description = "Small, stack-allocated linear algebra for fixed dimensions" +readme = "README.md" repository = "https://github.com/acgetchell/la-stack" categories = ["mathematics", "science"] keywords = ["linear-algebra", "geometry", "const-generics"] @@ -13,9 +14,11 @@ keywords = ["linear-algebra", "geometry", "const-generics"] # Intentionally empty [dev-dependencies] -approx = "0.5" -criterion = { version = "0.8", features = ["html_reports"] } -nalgebra = "0.34" +approx = "0.5.1" +criterion = { version = "0.8.1", features = ["html_reports"] } +nalgebra = "0.34.1" +pastey = "0.2.0" +proptest = "1.9.0" [[bench]] name = "vs_nalgebra" diff --git a/README.md b/README.md index d0707f0..bbafc55 100644 --- a/README.md +++ b/README.md @@ -2,7 +2,7 @@ [![Crates.io](https://img.shields.io/crates/v/la-stack.svg)](https://crates.io/crates/la-stack) [![Downloads](https://img.shields.io/crates/d/la-stack.svg)](https://crates.io/crates/la-stack) -[![License](https://img.shields.io/crates/l/la-stack.svg)](LICENSE) +[![License](https://img.shields.io/crates/l/la-stack.svg)](./LICENSE) [![Docs.rs](https://docs.rs/la-stack/badge.svg)](https://docs.rs/la-stack) [![CI](https://github.com/acgetchell/la-stack/actions/workflows/ci.yml/badge.svg)](https://github.com/acgetchell/la-stack/actions/workflows/ci.yml) [![rust-clippy analyze](https://github.com/acgetchell/la-stack/actions/workflows/rust-clippy.yml/badge.svg)](https://github.com/acgetchell/la-stack/actions/workflows/rust-clippy.yml) @@ -11,7 +11,8 @@ Fast, stack-allocated linear algebra for fixed dimensions in Rust. -This crate exists primarily to support the [`delaunay`](https://crates.io/crates/delaunay) crate’s needs while keeping the API intentionally small and explicit. +This crate grew from the need to support [`delaunay`](https://crates.io/crates/delaunay) with fast, stack-allocated linear algebra primitives and algorithms +while keeping the API intentionally small and explicit. ## πŸ“ Introduction @@ -45,21 +46,30 @@ Add this to your `Cargo.toml`: la-stack = "0.1" ``` -Solve a 3Γ—3 system via LU: +Solve a 5Γ—5 system via LU: ```rust -use la_stack::{LaError, Matrix, Vector, DEFAULT_PIVOT_TOL}; - -fn main() -> Result<(), LaError> { - // Requires pivoting (a[0][0] = 0), so it's a good LU demo. - let a = Matrix::<3>::from_rows([[0.0, 1.0, 1.0], [1.0, 0.0, 1.0], [1.0, 1.0, 0.0]]); - let b = Vector::<3>::new([5.0, 4.0, 3.0]); - - let lu = a.lu(DEFAULT_PIVOT_TOL)?; - let x = lu.solve_vec(b)?.into_array(); - - println!("x = {x:?}"); - Ok(()) +use la_stack::prelude::*; + +// This system requires pivoting (a[0][0] = 0), so it's a good LU demo. +// A = J - I: zeros on diagonal, ones elsewhere. +let a = Matrix::<5>::from_rows([ + [0.0, 1.0, 1.0, 1.0, 1.0], + [1.0, 0.0, 1.0, 1.0, 1.0], + [1.0, 1.0, 0.0, 1.0, 1.0], + [1.0, 1.0, 1.0, 0.0, 1.0], + [1.0, 1.0, 1.0, 1.0, 0.0], +]); + +let b = Vector::<5>::new([14.0, 13.0, 12.0, 11.0, 10.0]); + +let lu = a.lu(DEFAULT_PIVOT_TOL).unwrap(); +let x = lu.solve_vec(b).unwrap().into_array(); + +// Floating-point rounding is expected; compare with a tolerance. +let expected = [1.0, 2.0, 3.0, 4.0, 5.0]; +for (x_i, e_i) in x.iter().zip(expected.iter()) { + assert!((*x_i - *e_i).abs() <= 1e-12); } ``` @@ -80,8 +90,8 @@ The `examples/` directory contains small, runnable programs: ```bash just examples # or: -cargo run --example solve_3x3 -cargo run --example det_3x3 +cargo run --example solve_5x5 +cargo run --example det_5x5 ``` ## 🀝 Contributing @@ -98,4 +108,4 @@ For the full set of developer commands, see `just --list` and `WARP.md`. ## πŸ“„ License -BSD 3-Clause License. See [LICENSE](LICENSE). +BSD 3-Clause License. See [LICENSE](./LICENSE). diff --git a/WARP.md b/WARP.md index 4f209ee..334430f 100644 --- a/WARP.md +++ b/WARP.md @@ -21,7 +21,7 @@ When making changes in this repo, prioritize (in order): - Lint (Clippy): `cargo clippy --all-targets --all-features -- -D warnings` (or `just clippy`) - Spell check: `just spell-check` (uses `cspell.json` at repo root; keep the `words` list sorted lexicographically) - Run benchmarks: `cargo bench` (or `just bench`) -- Run examples: `just examples` (or `cargo run --example det_3x3`) +- Run examples: `just examples` (or `cargo run --example det_5x5` / `cargo run --example solve_5x5`) - CI simulation (lint + tests + bench compile): `just ci` - Pre-commit validation: `just commit-check` diff --git a/cspell.json b/cspell.json index 3662968..fc25c02 100644 --- a/cspell.json +++ b/cspell.json @@ -6,6 +6,7 @@ "clippy", "codacy", "const", + "doctests", "elif", "f128", "f32", @@ -14,18 +15,24 @@ "Getchell", "Justfile", "justfile", + "laerror", "lu", "markdownlint", "MSRV", "msvc", "mult", "nalgebra", + "nonfinite", + "pastey", "pipefail", "pivoting", "println", + "proptest", + "proptests", "rug", "RUSTDOCFLAGS", "shellcheck", + "tridiagonal", "usize" ], "ignorePaths": [ diff --git a/examples/det_3x3.rs b/examples/det_3x3.rs deleted file mode 100644 index a7e7522..0000000 --- a/examples/det_3x3.rs +++ /dev/null @@ -1,14 +0,0 @@ -//! Compute the determinant of a 3Γ—3 matrix via explicit LU factorization. - -use la_stack::{DEFAULT_PIVOT_TOL, LaError, Matrix}; - -fn main() -> Result<(), LaError> { - let a = Matrix::<3>::from_rows([[1.0, 2.0, 3.0], [0.0, 4.0, 5.0], [1.0, 0.0, 6.0]]); - - // Compute via explicit LU factorization. - let lu = a.lu(DEFAULT_PIVOT_TOL)?; - let det = lu.det(); - - println!("det = {det}"); - Ok(()) -} diff --git a/examples/det_5x5.rs b/examples/det_5x5.rs new file mode 100644 index 0000000..899151d --- /dev/null +++ b/examples/det_5x5.rs @@ -0,0 +1,22 @@ +//! Compute the determinant of a 5Γ—5 matrix via explicit LU factorization. + +use la_stack::prelude::*; + +fn main() -> Result<(), LaError> { + // 5Γ—5 matrix with zeros on diagonal and ones elsewhere (J - I). + // det(J - I) = (D - 1) * (-1)^(D-1) = 4 for D=5. + let a = Matrix::<5>::from_rows([ + [0.0, 1.0, 1.0, 1.0, 1.0], + [1.0, 0.0, 1.0, 1.0, 1.0], + [1.0, 1.0, 0.0, 1.0, 1.0], + [1.0, 1.0, 1.0, 0.0, 1.0], + [1.0, 1.0, 1.0, 1.0, 0.0], + ]); + + // Compute via explicit LU factorization. + let lu = a.lu(DEFAULT_PIVOT_TOL)?; + let det = lu.det(); + + println!("det = {det}"); + Ok(()) +} diff --git a/examples/solve_3x3.rs b/examples/solve_3x3.rs deleted file mode 100644 index 4048a83..0000000 --- a/examples/solve_3x3.rs +++ /dev/null @@ -1,15 +0,0 @@ -//! Solve a 3Γ—3 linear system via LU factorization (with pivoting). - -use la_stack::{DEFAULT_PIVOT_TOL, LaError, Matrix, Vector}; - -fn main() -> Result<(), LaError> { - // This system requires pivoting (a[0][0] = 0), so it's a good LU demo. - let a = Matrix::<3>::from_rows([[0.0, 1.0, 1.0], [1.0, 0.0, 1.0], [1.0, 1.0, 0.0]]); - let b = Vector::<3>::new([5.0, 4.0, 3.0]); - - let lu = a.lu(DEFAULT_PIVOT_TOL)?; - let x = lu.solve_vec(b)?.into_array(); - - println!("x = [{:.6}, {:.6}, {:.6}]", x[0], x[1], x[2]); - Ok(()) -} diff --git a/examples/solve_5x5.rs b/examples/solve_5x5.rs new file mode 100644 index 0000000..8515146 --- /dev/null +++ b/examples/solve_5x5.rs @@ -0,0 +1,24 @@ +//! Solve a 5Γ—5 linear system via LU factorization (with pivoting). + +use la_stack::prelude::*; + +fn main() -> Result<(), LaError> { + // This system requires pivoting (a[0][0] = 0), so it's a good LU demo. + // A = J - I: zeros on diagonal, ones elsewhere. + let a = Matrix::<5>::from_rows([ + [0.0, 1.0, 1.0, 1.0, 1.0], + [1.0, 0.0, 1.0, 1.0, 1.0], + [1.0, 1.0, 0.0, 1.0, 1.0], + [1.0, 1.0, 1.0, 0.0, 1.0], + [1.0, 1.0, 1.0, 1.0, 0.0], + ]); + + // Choose x = [1, 2, 3, 4, 5]. Then b = A x = [14, 13, 12, 11, 10]. + let b = Vector::<5>::new([14.0, 13.0, 12.0, 11.0, 10.0]); + + let lu = a.lu(DEFAULT_PIVOT_TOL)?; + let x = lu.solve_vec(b)?.into_array(); + + println!("x = {x:?}"); + Ok(()) +} diff --git a/justfile b/justfile index bd77637..21010ba 100644 --- a/justfile +++ b/justfile @@ -42,7 +42,7 @@ check: cargo check # CI simulation (matches delaunay's `just ci` shape) -ci: lint test test-integration bench-compile +ci: lint test-all bench-compile @echo "🎯 CI simulation complete!" clean: @@ -54,7 +54,7 @@ clippy: # Pre-commit workflow: comprehensive validation before committing # Runs: linting + all Rust tests (lib + doc + integration) + examples -commit-check: lint test-all examples +commit-check: lint test-all bench-compile examples @echo "πŸš€ Ready to commit! All checks passed!" # Coverage (cargo-tarpaulin) @@ -101,8 +101,8 @@ doc-check: # Examples examples: - cargo run --quiet --example det_3x3 - cargo run --quiet --example solve_3x3 + cargo run --quiet --example det_5x5 + cargo run --quiet --example solve_5x5 fmt: cargo fmt --all diff --git a/src/lib.rs b/src/lib.rs index e9f5341..bd4b323 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,12 +1,35 @@ #![forbid(unsafe_code)] #![warn(missing_docs)] +#![doc = include_str!("../README.md")] -//! Small, stack-allocated linear algebra for fixed dimensions. -//! -//! This crate intentionally focuses on a minimal, explicit API surface: -//! - const-generic sizes (no dynamic dimensions) -//! - stack-backed `Copy` types -//! - explicit algorithms (LU, solve, det) +#[cfg(doc)] +mod readme_doctests { + //! Executable versions of README examples. + /// ```rust + /// use la_stack::prelude::*; + /// + /// // This system requires pivoting (a[0][0] = 0), so it's a good LU demo. + /// let a = Matrix::<5>::from_rows([ + /// [0.0, 1.0, 1.0, 1.0, 1.0], + /// [1.0, 0.0, 1.0, 1.0, 1.0], + /// [1.0, 1.0, 0.0, 1.0, 1.0], + /// [1.0, 1.0, 1.0, 0.0, 1.0], + /// [1.0, 1.0, 1.0, 1.0, 0.0], + /// ]); + /// + /// let b = Vector::<5>::new([14.0, 13.0, 12.0, 11.0, 10.0]); + /// + /// let lu = a.lu(DEFAULT_PIVOT_TOL).unwrap(); + /// let x = lu.solve_vec(b).unwrap().into_array(); + /// + /// // Floating-point rounding is expected; compare with a tolerance. + /// let expected = [1.0, 2.0, 3.0, 4.0, 5.0]; + /// for (x_i, e_i) in x.iter().zip(expected.iter()) { + /// assert!((*x_i - *e_i).abs() <= 1e-12); + /// } + /// ``` + fn solve_5x5_example() {} +} mod lu; mod matrix; @@ -63,3 +86,47 @@ pub use vector::Vector; pub mod prelude { pub use crate::{DEFAULT_PIVOT_TOL, LaError, Lu, Matrix, Vector}; } + +#[cfg(test)] +mod tests { + use super::*; + + use approx::assert_abs_diff_eq; + + #[test] + fn default_pivot_tol_is_expected() { + assert_abs_diff_eq!(DEFAULT_PIVOT_TOL, 1e-12, epsilon = 0.0); + } + + #[test] + fn laerror_display_formats_singular() { + let err = LaError::Singular { pivot_col: 3 }; + assert_eq!(err.to_string(), "singular matrix at pivot column 3"); + } + + #[test] + fn laerror_display_formats_nonfinite() { + let err = LaError::NonFinite { pivot_col: 2 }; + assert_eq!( + err.to_string(), + "non-finite value encountered at pivot column 2" + ); + } + + #[test] + fn laerror_is_std_error_with_no_source() { + let err = LaError::Singular { pivot_col: 0 }; + let e: &dyn std::error::Error = &err; + assert!(e.source().is_none()); + } + + #[test] + fn prelude_reexports_compile_and_work() { + use crate::prelude::*; + + // Use the items so we know they are in scope and usable. + let m = Matrix::<2>::identity(); + let v = Vector::<2>::new([1.0, 2.0]); + let _ = m.lu(DEFAULT_PIVOT_TOL).unwrap().solve_vec(v).unwrap(); + } +} diff --git a/src/lu.rs b/src/lu.rs index 813cbf6..7c04bfe 100644 --- a/src/lu.rs +++ b/src/lu.rs @@ -180,34 +180,271 @@ mod tests { use super::*; use crate::DEFAULT_PIVOT_TOL; - fn assert_approx(a: f64, b: f64, eps: f64) { - assert!((a - b).abs() <= eps, "{a} !~= {b} (eps={eps})"); + use core::hint::black_box; + + use approx::assert_abs_diff_eq; + use pastey::paste; + + macro_rules! gen_public_api_pivoting_solve_vec_and_det_tests { + ($d:literal) => { + paste! { + #[test] + fn []() { + // Public API path under test: + // Matrix::lu (pub) -> Lu::solve_vec (pub). + + // Permutation matrix that swaps the first two basis vectors. + // This forces pivoting in column 0 for any D >= 2. + let mut rows = [[0.0f64; $d]; $d]; + for i in 0..$d { + rows[i][i] = 1.0; + } + rows.swap(0, 1); + + let a = Matrix::<$d>::from_rows(black_box(rows)); + let lu_fn: fn(Matrix<$d>, f64) -> Result, LaError> = + black_box(Matrix::<$d>::lu); + let lu = lu_fn(a, DEFAULT_PIVOT_TOL).unwrap(); + + // Pick a simple RHS with unique entries, so the expected swap is obvious. + let b_arr = { + let mut arr = [0.0f64; $d]; + let mut val = 1.0f64; + for dst in arr.iter_mut() { + *dst = val; + val += 1.0; + } + arr + }; + let mut expected = b_arr; + expected.swap(0, 1); + let b = Vector::<$d>::new(black_box(b_arr)); + + let solve_fn: fn(&Lu<$d>, Vector<$d>) -> Result, LaError> = + black_box(Lu::<$d>::solve_vec); + let x = solve_fn(&lu, b).unwrap().into_array(); + + for i in 0..$d { + assert_abs_diff_eq!(x[i], expected[i], epsilon = 1e-12); + } + } + + #[test] + fn []() { + // Public API path under test: + // Matrix::lu (pub) -> Lu::det (pub). + + // Permutation matrix that swaps the first two basis vectors. + let mut rows = [[0.0f64; $d]; $d]; + for i in 0..$d { + rows[i][i] = 1.0; + } + rows.swap(0, 1); + + let a = Matrix::<$d>::from_rows(black_box(rows)); + let lu_fn: fn(Matrix<$d>, f64) -> Result, LaError> = + black_box(Matrix::<$d>::lu); + let lu = lu_fn(a, DEFAULT_PIVOT_TOL).unwrap(); + + // Row swap β‡’ determinant sign flip. + let det_fn: fn(&Lu<$d>) -> f64 = black_box(Lu::<$d>::det); + assert_abs_diff_eq!(det_fn(&lu), -1.0, epsilon = 1e-12); + } + } + }; + } + + gen_public_api_pivoting_solve_vec_and_det_tests!(2); + gen_public_api_pivoting_solve_vec_and_det_tests!(3); + gen_public_api_pivoting_solve_vec_and_det_tests!(4); + gen_public_api_pivoting_solve_vec_and_det_tests!(5); + + macro_rules! gen_public_api_tridiagonal_smoke_solve_vec_and_det_tests { + ($d:literal) => { + paste! { + #[test] + fn []() { + // Public API path under test: + // Matrix::lu (pub) -> Lu::solve_vec (pub). + + // Classic SPD tridiagonal: 2 on diagonal, -1 on sub/super-diagonals. + #[allow(clippy::large_stack_arrays)] + let mut rows = [[0.0f64; $d]; $d]; + for i in 0..$d { + rows[i][i] = 2.0; + if i > 0 { + rows[i][i - 1] = -1.0; + } + if i + 1 < $d { + rows[i][i + 1] = -1.0; + } + } + + let a = Matrix::<$d>::from_rows(black_box(rows)); + let lu_fn: fn(Matrix<$d>, f64) -> Result, LaError> = + black_box(Matrix::<$d>::lu); + let lu = lu_fn(a, DEFAULT_PIVOT_TOL).unwrap(); + + // Choose x = 1, so b = A x is simple: [1, 0, 0, ..., 0, 1]. + let mut b_arr = [0.0f64; $d]; + b_arr[0] = 1.0; + b_arr[$d - 1] = 1.0; + let b = Vector::<$d>::new(black_box(b_arr)); + + let solve_fn: fn(&Lu<$d>, Vector<$d>) -> Result, LaError> = + black_box(Lu::<$d>::solve_vec); + let x = solve_fn(&lu, b).unwrap().into_array(); + + for &x_i in &x { + assert_abs_diff_eq!(x_i, 1.0, epsilon = 1e-9); + } + } + + #[test] + fn []() { + // Public API path under test: + // Matrix::lu (pub) -> Lu::det (pub). + + // Classic SPD tridiagonal: 2 on diagonal, -1 on sub/super-diagonals. + // Determinant is known exactly: det = D + 1. + #[allow(clippy::large_stack_arrays)] + let mut rows = [[0.0f64; $d]; $d]; + for i in 0..$d { + rows[i][i] = 2.0; + if i > 0 { + rows[i][i - 1] = -1.0; + } + if i + 1 < $d { + rows[i][i + 1] = -1.0; + } + } + + let a = Matrix::<$d>::from_rows(black_box(rows)); + let lu_fn: fn(Matrix<$d>, f64) -> Result, LaError> = + black_box(Matrix::<$d>::lu); + let lu = lu_fn(a, DEFAULT_PIVOT_TOL).unwrap(); + + let det_fn: fn(&Lu<$d>) -> f64 = black_box(Lu::<$d>::det); + assert_abs_diff_eq!(det_fn(&lu), f64::from($d) + 1.0, epsilon = 1e-8); + } + } + }; + } + + gen_public_api_tridiagonal_smoke_solve_vec_and_det_tests!(16); + gen_public_api_tridiagonal_smoke_solve_vec_and_det_tests!(32); + gen_public_api_tridiagonal_smoke_solve_vec_and_det_tests!(64); + + #[test] + fn solve_1x1() { + let a = Matrix::<1>::from_rows(black_box([[2.0]])); + let lu = (black_box(Lu::<1>::factor))(a, DEFAULT_PIVOT_TOL).unwrap(); + + let b = Vector::<1>::new(black_box([6.0])); + let solve_fn: fn(&Lu<1>, Vector<1>) -> Result, LaError> = + black_box(Lu::<1>::solve_vec); + let x = solve_fn(&lu, b).unwrap().into_array(); + assert_abs_diff_eq!(x[0], 3.0, epsilon = 1e-12); + + let det_fn: fn(&Lu<1>) -> f64 = black_box(Lu::<1>::det); + assert_abs_diff_eq!(det_fn(&lu), 2.0, epsilon = 0.0); } #[test] fn solve_2x2_basic() { - let a = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]); - let lu = Lu::factor(a, DEFAULT_PIVOT_TOL).unwrap(); - let b = Vector::<2>::new([5.0, 11.0]); - let x = lu.solve_vec(b).unwrap().into_array(); - assert_approx(x[0], 1.0, 1e-12); - assert_approx(x[1], 2.0, 1e-12); + let a = Matrix::<2>::from_rows(black_box([[1.0, 2.0], [3.0, 4.0]])); + let lu = (black_box(Lu::<2>::factor))(a, DEFAULT_PIVOT_TOL).unwrap(); + let b = Vector::<2>::new(black_box([5.0, 11.0])); + + let solve_fn: fn(&Lu<2>, Vector<2>) -> Result, LaError> = + black_box(Lu::<2>::solve_vec); + let x = solve_fn(&lu, b).unwrap().into_array(); + + assert_abs_diff_eq!(x[0], 1.0, epsilon = 1e-12); + assert_abs_diff_eq!(x[1], 2.0, epsilon = 1e-12); + } + + #[test] + fn det_2x2_basic() { + let a = Matrix::<2>::from_rows(black_box([[1.0, 2.0], [3.0, 4.0]])); + let lu = a.lu(DEFAULT_PIVOT_TOL).unwrap(); + + let det_fn: fn(&Lu<2>) -> f64 = black_box(Lu::<2>::det); + assert_abs_diff_eq!(det_fn(&lu), -2.0, epsilon = 1e-12); + } + + #[test] + fn det_requires_pivot_sign() { + // Row swap β‡’ determinant sign flip. + let a = Matrix::<2>::from_rows(black_box([[0.0, 1.0], [1.0, 0.0]])); + let lu = a.lu(DEFAULT_PIVOT_TOL).unwrap(); + + let det_fn: fn(&Lu<2>) -> f64 = black_box(Lu::<2>::det); + assert_abs_diff_eq!(det_fn(&lu), -1.0, epsilon = 0.0); } #[test] fn solve_requires_pivoting() { - let a = Matrix::<2>::from_rows([[0.0, 1.0], [1.0, 0.0]]); + let a = Matrix::<2>::from_rows(black_box([[0.0, 1.0], [1.0, 0.0]])); let lu = a.lu(DEFAULT_PIVOT_TOL).unwrap(); - let b = Vector::<2>::new([1.0, 2.0]); - let x = lu.solve_vec(b).unwrap().into_array(); - assert_approx(x[0], 2.0, 1e-12); - assert_approx(x[1], 1.0, 1e-12); + let b = Vector::<2>::new(black_box([1.0, 2.0])); + + let solve_fn: fn(&Lu<2>, Vector<2>) -> Result, LaError> = + black_box(Lu::<2>::solve_vec); + let x = solve_fn(&lu, b).unwrap().into_array(); + + assert_abs_diff_eq!(x[0], 2.0, epsilon = 1e-12); + assert_abs_diff_eq!(x[1], 1.0, epsilon = 1e-12); } #[test] fn singular_detected() { - let a = Matrix::<2>::from_rows([[1.0, 2.0], [2.0, 4.0]]); + let a = Matrix::<2>::from_rows(black_box([[1.0, 2.0], [2.0, 4.0]])); let err = a.lu(DEFAULT_PIVOT_TOL).unwrap_err(); assert_eq!(err, LaError::Singular { pivot_col: 1 }); } + + #[test] + fn singular_due_to_tolerance_at_first_pivot() { + // Not exactly singular, but below DEFAULT_PIVOT_TOL. + let a = Matrix::<2>::from_rows(black_box([[1e-13, 0.0], [0.0, 1.0]])); + let err = a.lu(DEFAULT_PIVOT_TOL).unwrap_err(); + assert_eq!(err, LaError::Singular { pivot_col: 0 }); + } + + #[test] + fn nonfinite_detected_on_pivot_entry() { + let a = Matrix::<2>::from_rows([[f64::NAN, 0.0], [0.0, 1.0]]); + let err = a.lu(DEFAULT_PIVOT_TOL).unwrap_err(); + assert_eq!(err, LaError::NonFinite { pivot_col: 0 }); + } + + #[test] + fn nonfinite_detected_in_pivot_column_scan() { + let a = Matrix::<2>::from_rows([[1.0, 0.0], [f64::INFINITY, 1.0]]); + let err = a.lu(DEFAULT_PIVOT_TOL).unwrap_err(); + assert_eq!(err, LaError::NonFinite { pivot_col: 0 }); + } + + #[test] + fn solve_vec_nonfinite_forward_substitution_overflow() { + // L has a -1 multiplier, and a large RHS makes forward substitution overflow. + let a = Matrix::<3>::from_rows([[1.0, 0.0, 0.0], [-1.0, 1.0, 0.0], [0.0, 0.0, 1.0]]); + let lu = a.lu(DEFAULT_PIVOT_TOL).unwrap(); + + let b = Vector::<3>::new([1.0e308, 1.0e308, 0.0]); + let err = lu.solve_vec(b).unwrap_err(); + assert_eq!(err, LaError::NonFinite { pivot_col: 1 }); + } + + #[test] + fn solve_vec_nonfinite_back_substitution_overflow() { + // Make x[1] overflow during back substitution, then ensure it is detected on the next row. + let a = Matrix::<2>::from_rows([[1.0, 1.0], [0.0, 2.0e-12]]); + let lu = a.lu(DEFAULT_PIVOT_TOL).unwrap(); + + let b = Vector::<2>::new([0.0, 1.0e300]); + let err = lu.solve_vec(b).unwrap_err(); + assert_eq!(err, LaError::NonFinite { pivot_col: 0 }); + } } diff --git a/src/matrix.rs b/src/matrix.rs index bb3e11e..9acf314 100644 --- a/src/matrix.rs +++ b/src/matrix.rs @@ -200,47 +200,103 @@ mod tests { use super::*; use crate::DEFAULT_PIVOT_TOL; - fn assert_approx(a: f64, b: f64, eps: f64) { - assert!((a - b).abs() <= eps, "{a} !~= {b} (eps={eps})"); - } + use approx::assert_abs_diff_eq; + use pastey::paste; - #[test] - fn get_set_bounds_checked() { - let mut m = Matrix::<2>::zero(); - assert!(m.set(0, 0, 1.0)); - assert_eq!(m.get(0, 0), Some(1.0)); + macro_rules! gen_public_api_matrix_tests { + ($d:literal) => { + paste! { + #[test] + fn []() { + let mut rows = [[0.0f64; $d]; $d]; + rows[0][0] = 1.0; + rows[$d - 1][$d - 1] = -2.0; - assert!(!m.set(2, 0, 1.0)); - assert_eq!(m.get(2, 0), None); - } + let mut m = Matrix::<$d>::from_rows(rows); - #[test] - fn inf_norm_max_row_sum() { - let m = Matrix::<2>::from_rows([[1.0, -2.0], [3.0, 4.0]]); - assert_approx(m.inf_norm(), 7.0, 0.0); - } + assert_eq!(m.get(0, 0), Some(1.0)); + assert_eq!(m.get($d - 1, $d - 1), Some(-2.0)); - #[test] - fn det_identity_is_one() { - let det = Matrix::<3>::identity().det(DEFAULT_PIVOT_TOL).unwrap(); - assert_approx(det, 1.0, 1e-12); - } + // Out-of-bounds is None. + assert_eq!(m.get($d, 0), None); + + // Out-of-bounds set fails. + assert!(!m.set($d, 0, 3.0)); + + // In-bounds set works. + assert!(m.set(0, $d - 1, 3.0)); + assert_eq!(m.get(0, $d - 1), Some(3.0)); + } + + #[test] + fn []() { + let z = Matrix::<$d>::zero(); + assert_abs_diff_eq!(z.inf_norm(), 0.0, epsilon = 0.0); + + let d = Matrix::<$d>::default(); + assert_abs_diff_eq!(d.inf_norm(), 0.0, epsilon = 0.0); + } + + #[test] + fn []() { + let mut rows = [[0.0f64; $d]; $d]; + + // Row 0 has absolute row sum = D. + for c in 0..$d { + rows[0][c] = -1.0; + } - #[test] - fn identity_has_ones_on_diag_and_zeros_off_diag() { - let m = Matrix::<3>::identity(); + // Row 1 has smaller absolute row sum. + for c in 0..$d { + rows[1][c] = 0.5; + } - for r in 0..3 { - for c in 0..3 { - let expected = if r == c { 1.0 } else { 0.0 }; - assert_approx(m.get(r, c).unwrap(), expected, 0.0); + let m = Matrix::<$d>::from_rows(rows); + assert_abs_diff_eq!(m.inf_norm(), f64::from($d), epsilon = 0.0); + } + + #[test] + fn []() { + let m = Matrix::<$d>::identity(); + + // Identity has ones on diag and zeros off diag. + for r in 0..$d { + for c in 0..$d { + let expected = if r == c { 1.0 } else { 0.0 }; + assert_abs_diff_eq!(m.get(r, c).unwrap(), expected, epsilon = 0.0); + } + } + + // Determinant is 1. + let det = m.det(DEFAULT_PIVOT_TOL).unwrap(); + assert_abs_diff_eq!(det, 1.0, epsilon = 1e-12); + + // LU solve on identity returns the RHS. + let lu = m.lu(DEFAULT_PIVOT_TOL).unwrap(); + + let b_arr = { + let mut arr = [0.0f64; $d]; + let values = [1.0f64, 2.0, 3.0, 4.0, 5.0]; + for (dst, src) in arr.iter_mut().zip(values.iter()) { + *dst = *src; + } + arr + }; + + let b = crate::Vector::<$d>::new(b_arr); + let x = lu.solve_vec(b).unwrap().into_array(); + + for (x_i, b_i) in x.iter().zip(b_arr.iter()) { + assert_abs_diff_eq!(*x_i, *b_i, epsilon = 1e-12); + } + } } - } + }; } - #[test] - fn default_is_zero() { - let m = Matrix::<3>::default(); - assert_approx(m.inf_norm(), 0.0, 0.0); - } + // Mirror delaunay-style multi-dimension tests. + gen_public_api_matrix_tests!(2); + gen_public_api_matrix_tests!(3); + gen_public_api_matrix_tests!(4); + gen_public_api_matrix_tests!(5); } diff --git a/src/vector.rs b/src/vector.rs index e4171e3..e63c3c2 100644 --- a/src/vector.rs +++ b/src/vector.rs @@ -122,41 +122,112 @@ impl Default for Vector { mod tests { use super::*; - fn assert_approx(a: f64, b: f64, eps: f64) { - assert!((a - b).abs() <= eps, "{a} !~= {b} (eps={eps})"); + use core::hint::black_box; + + use approx::assert_abs_diff_eq; + use pastey::paste; + + macro_rules! gen_public_api_vector_tests { + ($d:literal) => { + paste! { + #[test] + fn []() { + let arr = { + let mut arr = [0.0f64; $d]; + let values = [1.0f64, 2.0, 3.0, 4.0, 5.0]; + for (dst, src) in arr.iter_mut().zip(values.iter()) { + *dst = *src; + } + arr + }; + + let v = Vector::<$d>::new(arr); + + for i in 0..$d { + assert_abs_diff_eq!(v.as_array()[i], arr[i], epsilon = 0.0); + } + + let out = v.into_array(); + for i in 0..$d { + assert_abs_diff_eq!(out[i], arr[i], epsilon = 0.0); + } + } + + #[test] + fn []() { + let z = Vector::<$d>::zero(); + for &x in z.as_array() { + assert_abs_diff_eq!(x, 0.0, epsilon = 0.0); + } + for x in z.into_array() { + assert_abs_diff_eq!(x, 0.0, epsilon = 0.0); + } + + let d = Vector::<$d>::default(); + for x in d.into_array() { + assert_abs_diff_eq!(x, 0.0, epsilon = 0.0); + } + } + + #[test] + fn []() { + // Use black_box to avoid constant-folding/inlining eliminating the actual dot loop, + // which can make coverage tools report the mul_add line as uncovered. + + let a_arr = { + let mut arr = [0.0f64; $d]; + let values = [1.0f64, 2.0, 3.0, 4.0, 5.0]; + for (dst, src) in arr.iter_mut().zip(values.iter()) { + *dst = black_box(*src); + } + arr + }; + let b_arr = { + let mut arr = [0.0f64; $d]; + let values = [-2.0f64, 0.5, 4.0, -1.0, 2.0]; + for (dst, src) in arr.iter_mut().zip(values.iter()) { + *dst = black_box(*src); + } + arr + }; + + let expected_dot = { + let mut acc = 0.0; + let mut i = 0; + while i < $d { + acc = a_arr[i].mul_add(b_arr[i], acc); + i += 1; + } + acc + }; + let expected_norm2_sq = { + let mut acc = 0.0; + let mut i = 0; + while i < $d { + acc = a_arr[i].mul_add(a_arr[i], acc); + i += 1; + } + acc + }; + + let a = Vector::<$d>::new(black_box(a_arr)); + let b = Vector::<$d>::new(black_box(b_arr)); + + // Call via (black_boxed) fn pointers to discourage inlining, improving line-level coverage + // attribution for the loop body. + let dot_fn: fn(Vector<$d>, Vector<$d>) -> f64 = black_box(Vector::<$d>::dot); + let norm2_sq_fn: fn(Vector<$d>) -> f64 = black_box(Vector::<$d>::norm2_sq); + + assert_abs_diff_eq!(dot_fn(black_box(a), black_box(b)), expected_dot, epsilon = 1e-14); + assert_abs_diff_eq!(norm2_sq_fn(black_box(a)), expected_norm2_sq, epsilon = 1e-14); + } + } + }; } - #[test] - fn dot_and_norm2_sq() { - // Use black_box to avoid constant-folding/inlining eliminating the actual dot loop, - // which can make coverage tools report the mul_add line as uncovered. - use core::hint::black_box; - - let a = Vector::<3>::new([black_box(1.0), black_box(2.0), black_box(3.0)]); - let b = Vector::<3>::new([black_box(-2.0), black_box(0.5), black_box(4.0)]); - - // Call via (black_boxed) fn pointers to discourage inlining, improving line-level coverage - // attribution for the loop body. - let dot_fn: fn(Vector<3>, Vector<3>) -> f64 = black_box(Vector::<3>::dot); - let norm2_sq_fn: fn(Vector<3>) -> f64 = black_box(Vector::<3>::norm2_sq); - - assert_approx(dot_fn(black_box(a), black_box(b)), 11.0, 0.0); - assert_approx(norm2_sq_fn(black_box(a)), 14.0, 0.0); - } - - #[test] - fn zero_as_array_into_array_and_default() { - let z = Vector::<3>::zero(); - for &x in z.as_array() { - assert_approx(x, 0.0, 0.0); - } - for x in z.into_array() { - assert_approx(x, 0.0, 0.0); - } - - let d = Vector::<3>::default(); - for x in d.into_array() { - assert_approx(x, 0.0, 0.0); - } - } + // Mirror delaunay-style multi-dimension tests. + gen_public_api_vector_tests!(2); + gen_public_api_vector_tests!(3); + gen_public_api_vector_tests!(4); + gen_public_api_vector_tests!(5); } diff --git a/tests/proptest_matrix.rs b/tests/proptest_matrix.rs new file mode 100644 index 0000000..c6795bf --- /dev/null +++ b/tests/proptest_matrix.rs @@ -0,0 +1,110 @@ +//! Property-based tests for the `Matrix` public API. + +use approx::assert_abs_diff_eq; +use pastey::paste; +use proptest::prelude::*; + +use la_stack::prelude::*; + +fn small_f64() -> impl Strategy { + (-1000i16..=1000i16).prop_map(|x| f64::from(x) / 10.0) +} + +fn small_nonzero_f64() -> impl Strategy { + prop_oneof![(-1000i16..=-1i16), (1i16..=1000i16)].prop_map(|x| f64::from(x) / 10.0) +} + +macro_rules! gen_public_api_matrix_proptests { + ($d:literal) => { + paste! { + proptest! { + #![proptest_config(ProptestConfig::with_cases(64))] + + #[test] + fn []( + rows in proptest::array::[]( + proptest::array::[](small_f64()), + ), + ) { + let m = Matrix::<$d>::from_rows(rows); + + for r in 0..$d { + for c in 0..$d { + assert_abs_diff_eq!(m.get(r, c).unwrap(), rows[r][c], epsilon = 0.0); + } + } + + // Out-of-bounds is None. + prop_assert_eq!(m.get($d, 0), None); + prop_assert_eq!(m.get(0, $d), None); + } + + #[test] + fn []( + r in 0usize..$d, + c in 0usize..$d, + v in small_f64(), + ) { + let mut m = Matrix::<$d>::zero(); + prop_assert!(m.set(r, c, v)); + assert_abs_diff_eq!(m.get(r, c).unwrap(), v, epsilon = 0.0); + } + + #[test] + fn []( + rows in proptest::array::[]( + proptest::array::[](small_f64()), + ), + ) { + let m = Matrix::<$d>::from_rows(rows); + + let expected = rows + .iter() + .map(|row| row.iter().map(|&x| x.abs()).sum::()) + .fold(0.0f64, f64::max); + + assert_abs_diff_eq!(m.inf_norm(), expected, epsilon = 0.0); + prop_assert!(m.inf_norm() >= 0.0); + } + + #[test] + fn []( + diag in proptest::array::[](small_nonzero_f64()), + b_arr in proptest::array::[](small_f64()), + ) { + // Diagonal matrix: det is product of diagonal, and solve is element-wise division. + let mut rows = [[0.0f64; $d]; $d]; + for i in 0..$d { + rows[i][i] = diag[i]; + } + let a = Matrix::<$d>::from_rows(rows); + + let det = a.det(DEFAULT_PIVOT_TOL).unwrap(); + let expected_det = { + let mut acc = 1.0; + for i in 0..$d { + acc *= diag[i]; + } + acc + }; + assert_abs_diff_eq!(det, expected_det, epsilon = 1e-12); + + let lu = a.lu(DEFAULT_PIVOT_TOL).unwrap(); + let b = Vector::<$d>::new(b_arr); + let x = lu.solve_vec(b).unwrap().into_array(); + + for i in 0..$d { + let expected_x = b_arr[i] / diag[i]; + assert_abs_diff_eq!(x[i], expected_x, epsilon = 1e-12); + } + } + } + } + }; +} + +// Mirror delaunay-style multi-dimension tests. +gen_public_api_matrix_proptests!(2); +gen_public_api_matrix_proptests!(3); +gen_public_api_matrix_proptests!(4); +gen_public_api_matrix_proptests!(5); diff --git a/tests/proptest_vector.rs b/tests/proptest_vector.rs new file mode 100644 index 0000000..11d754b --- /dev/null +++ b/tests/proptest_vector.rs @@ -0,0 +1,66 @@ +//! Property-based tests for the `Vector` public API. + +use approx::assert_abs_diff_eq; +use pastey::paste; +use proptest::prelude::*; + +use la_stack::prelude::*; + +fn small_f64() -> impl Strategy { + (-1000i16..=1000i16).prop_map(|x| f64::from(x) / 10.0) +} + +macro_rules! gen_public_api_vector_proptests { + ($d:literal) => { + paste! { + proptest! { + #![proptest_config(ProptestConfig::with_cases(64))] + + #[test] + fn []( + arr in proptest::array::[](small_f64()), + ) { + let v = Vector::<$d>::new(arr); + + for i in 0..$d { + assert_abs_diff_eq!(v.as_array()[i], arr[i], epsilon = 0.0); + } + + let out = v.into_array(); + for i in 0..$d { + assert_abs_diff_eq!(out[i], arr[i], epsilon = 0.0); + } + } + + #[test] + fn []( + a_arr in proptest::array::[](small_f64()), + b_arr in proptest::array::[](small_f64()), + ) { + let a = Vector::<$d>::new(a_arr); + let b = Vector::<$d>::new(b_arr); + + let dot_ab = a.dot(b); + let dot_ba = b.dot(a); + assert_abs_diff_eq!(dot_ab, dot_ba, epsilon = 1e-14); + + let dot_aa = a.dot(a); + assert_abs_diff_eq!(a.norm2_sq(), dot_aa, epsilon = 0.0); + + // Squared norm is always non-negative for finite inputs. + prop_assert!(a.norm2_sq() >= 0.0); + + // Dot with zero vector is zero. + let z = Vector::<$d>::zero(); + assert_abs_diff_eq!(a.dot(z), 0.0, epsilon = 1e-14); + } + } + } + }; +} + +// Mirror delaunay-style multi-dimension tests. +gen_public_api_vector_proptests!(2); +gen_public_api_vector_proptests!(3); +gen_public_api_vector_proptests!(4); +gen_public_api_vector_proptests!(5);