nalgebra-lapack: run tests with proptest instead of quickcheck.

This commit is contained in:
Crozet Sébastien 2021-02-28 18:39:18 +01:00
parent 6cfd2bca14
commit 74f4b0ba4d
11 changed files with 224 additions and 236 deletions

View File

@ -18,9 +18,11 @@ maintenance = { status = "actively-developed" }
[features]
serde-serialize = [ "serde", "serde_derive" ]
proptest-support = [ "nalgebra/proptest-support" ]
arbitrary = [ "nalgebra/arbitrary" ]
# For BLAS/LAPACK
default = ["openblas"]
default = ["netlib"]
openblas = ["lapack-src/openblas"]
netlib = ["lapack-src/netlib"]
accelerate = ["lapack-src/accelerate"]
@ -39,6 +41,7 @@ lapack-src = { version = "0.6", default-features = false }
[dev-dependencies]
nalgebra = { version = "0.24", features = [ "arbitrary" ], path = ".." }
proptest = { version = "1", default-features = false, features = ["std"] }
quickcheck = "1"
approx = "0.4"
rand = "0.8"

View File

@ -1,8 +1,14 @@
#[macro_use]
extern crate approx;
#[cfg(not(feature = "proptest-support"))]
compile_error!("Tests must be run with `proptest-support`");
extern crate nalgebra as na;
extern crate nalgebra_lapack as nl;
#[macro_use]
extern crate quickcheck;
extern crate lapack;
extern crate lapack_src;
mod linalg;
#[path = "../../tests/proptest/mod.rs"]
mod proptest;

View File

@ -1,37 +1,36 @@
use std::cmp;
use na::{DMatrix, DVector, Matrix3, Matrix4, Matrix4x3, Vector4};
use na::{DMatrix, DVector, Matrix4x3, Vector4};
use nl::Cholesky;
quickcheck! {
fn cholesky(m: DMatrix<f64>) -> bool {
if m.len() != 0 {
use crate::proptest::*;
use proptest::{prop_assert, proptest};
proptest! {
#[test]
fn cholesky(m in dmatrix()) {
let m = &m * m.transpose();
if let Some(chol) = Cholesky::new(m.clone()) {
let l = chol.unpack();
let reconstructed_m = &l * l.transpose();
return relative_eq!(reconstructed_m, m, epsilon = 1.0e-7)
prop_assert!(relative_eq!(reconstructed_m, m, epsilon = 1.0e-7));
}
}
return true
}
fn cholesky_static(m: Matrix3<f64>) -> bool {
#[test]
fn cholesky_static(m in matrix3()) {
let m = &m * m.transpose();
if let Some(chol) = Cholesky::new(m) {
let l = chol.unpack();
let reconstructed_m = &l * l.transpose();
relative_eq!(reconstructed_m, m, epsilon = 1.0e-7)
}
else {
false
prop_assert!(relative_eq!(reconstructed_m, m, epsilon = 1.0e-7))
}
}
fn cholesky_solve(n: usize, nb: usize) -> bool {
if n != 0 {
#[test]
fn cholesky_solve(n in PROPTEST_MATRIX_DIM, nb in PROPTEST_MATRIX_DIM) {
let n = cmp::min(n, 15); // To avoid slowing down the test too much.
let nb = cmp::min(nb, 15); // To avoid slowing down the test too much.
let m = DMatrix::<f64>::new_random(n, n);
@ -44,33 +43,28 @@ quickcheck! {
let sol1 = chol.solve(&b1).unwrap();
let sol2 = chol.solve(&b2).unwrap();
return relative_eq!(&m * sol1, b1, epsilon = 1.0e-6) &&
relative_eq!(&m * sol2, b2, epsilon = 1.0e-6)
prop_assert!(relative_eq!(&m * sol1, b1, epsilon = 1.0e-6));
prop_assert!(relative_eq!(&m * sol2, b2, epsilon = 1.0e-6));
}
}
return true;
}
fn cholesky_solve_static(m: Matrix4<f64>) -> bool {
#[test]
fn cholesky_solve_static(m in matrix4()) {
let m = &m * m.transpose();
match Cholesky::new(m) {
Some(chol) => {
if let Some(chol) = Cholesky::new(m) {
let b1 = Vector4::new_random();
let b2 = Matrix4x3::new_random();
let sol1 = chol.solve(&b1).unwrap();
let sol2 = chol.solve(&b2).unwrap();
relative_eq!(m * sol1, b1, epsilon = 1.0e-7) &&
relative_eq!(m * sol2, b2, epsilon = 1.0e-7)
},
None => true
prop_assert!(relative_eq!(m * sol1, b1, epsilon = 1.0e-7));
prop_assert!(relative_eq!(m * sol2, b2, epsilon = 1.0e-7));
}
}
fn cholesky_inverse(n: usize) -> bool {
if n != 0 {
#[test]
fn cholesky_inverse(n in PROPTEST_MATRIX_DIM) {
let n = cmp::min(n, 15); // To avoid slowing down the test too much.
let m = DMatrix::<f64>::new_random(n, n);
let m = &m * m.transpose();
@ -79,23 +73,18 @@ quickcheck! {
let id1 = &m * &m1;
let id2 = &m1 * &m;
return id1.is_identity(1.0e-6) && id2.is_identity(1.0e-6);
prop_assert!(id1.is_identity(1.0e-6) && id2.is_identity(1.0e-6));
}
}
return true;
}
fn cholesky_inverse_static(m: Matrix4<f64>) -> bool {
#[test]
fn cholesky_inverse_static(m in matrix4()) {
let m = m * m.transpose();
match Cholesky::new(m.clone()).unwrap().inverse() {
Some(m1) => {
if let Some(m1) = Cholesky::new(m.clone()).unwrap().inverse() {
let id1 = &m * &m1;
let id2 = &m1 * &m;
id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5)
},
None => true
prop_assert!(id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5))
}
}
}

View File

@ -1,38 +1,32 @@
use std::cmp;
use nl::Hessenberg;
use na::{DMatrix, Matrix4};
use nl::Hessenberg;
quickcheck!{
fn hessenberg(n: usize) -> bool {
if n != 0 {
use crate::proptest::*;
use proptest::{prop_assert, proptest};
proptest! {
#[test]
fn hessenberg(n in PROPTEST_MATRIX_DIM) {
let n = cmp::min(n, 25);
let m = DMatrix::<f64>::new_random(n, n);
match Hessenberg::new(m.clone()) {
Some(hess) => {
if let Some(hess) = Hessenberg::new(m.clone()) {
let h = hess.h();
let p = hess.p();
relative_eq!(m, &p * h * p.transpose(), epsilon = 1.0e-7)
},
None => true
}
}
else {
true
prop_assert!(relative_eq!(m, &p * h * p.transpose(), epsilon = 1.0e-7))
}
}
fn hessenberg_static(m: Matrix4<f64>) -> bool {
match Hessenberg::new(m) {
Some(hess) => {
#[test]
fn hessenberg_static(m in matrix4()) {
if let Some(hess) = Hessenberg::new(m) {
let h = hess.h();
let p = hess.p();
relative_eq!(m, p * h * p.transpose(), epsilon = 1.0e-7)
},
None => true
prop_assert!(relative_eq!(m, p * h * p.transpose(), epsilon = 1.0e-7))
}
}
}

View File

@ -1,11 +1,14 @@
use std::cmp;
use na::{DMatrix, DVector, Matrix3x4, Matrix4, Matrix4x3, Vector4};
use na::{DMatrix, DVector, Matrix4x3, Vector4};
use nl::LU;
quickcheck! {
fn lup(m: DMatrix<f64>) -> bool {
if m.len() != 0 {
use crate::proptest::*;
use proptest::{prop_assert, proptest};
proptest! {
#[test]
fn lup(m in dmatrix()) {
let lup = LU::new(m.clone());
let l = lup.l();
let u = lup.u();
@ -14,15 +17,12 @@ quickcheck! {
let computed2 = lup.p() * l * u;
relative_eq!(computed1, m, epsilon = 1.0e-7) &&
relative_eq!(computed2, m, epsilon = 1.0e-7)
}
else {
true
}
prop_assert!(relative_eq!(computed1, m, epsilon = 1.0e-7));
prop_assert!(relative_eq!(computed2, m, epsilon = 1.0e-7));
}
fn lu_static(m: Matrix3x4<f64>) -> bool {
#[test]
fn lu_static(m in matrix3x5()) {
let lup = LU::new(m);
let l = lup.l();
let u = lup.u();
@ -31,12 +31,12 @@ quickcheck! {
let computed2 = lup.p() * l * u;
relative_eq!(computed1, m, epsilon = 1.0e-7) &&
relative_eq!(computed2, m, epsilon = 1.0e-7)
prop_assert!(relative_eq!(computed1, m, epsilon = 1.0e-7));
prop_assert!(relative_eq!(computed2, m, epsilon = 1.0e-7));
}
fn lu_solve(n: usize, nb: usize) -> bool {
if n != 0 {
#[test]
fn lu_solve(n in PROPTEST_MATRIX_DIM, nb in PROPTEST_MATRIX_DIM) {
let n = cmp::min(n, 25); // To avoid slowing down the test too much.
let nb = cmp::min(nb, 25); // To avoid slowing down the test too much.
let m = DMatrix::<f64>::new_random(n, n);
@ -51,17 +51,14 @@ quickcheck! {
let tr_sol1 = lup.solve_transpose(&b1).unwrap();
let tr_sol2 = lup.solve_transpose(&b2).unwrap();
relative_eq!(&m * sol1, b1, epsilon = 1.0e-7) &&
relative_eq!(&m * sol2, b2, epsilon = 1.0e-7) &&
relative_eq!(m.transpose() * tr_sol1, b1, epsilon = 1.0e-7) &&
relative_eq!(m.transpose() * tr_sol2, b2, epsilon = 1.0e-7)
}
else {
true
}
prop_assert!(relative_eq!(&m * sol1, b1, epsilon = 1.0e-7));
prop_assert!(relative_eq!(&m * sol2, b2, epsilon = 1.0e-7));
prop_assert!(relative_eq!(m.transpose() * tr_sol1, b1, epsilon = 1.0e-7));
prop_assert!(relative_eq!(m.transpose() * tr_sol2, b2, epsilon = 1.0e-7));
}
fn lu_solve_static(m: Matrix4<f64>) -> bool {
#[test]
fn lu_solve_static(m in matrix4()) {
let lup = LU::new(m);
let b1 = Vector4::new_random();
let b2 = Matrix4x3::new_random();
@ -71,14 +68,14 @@ quickcheck! {
let tr_sol1 = lup.solve_transpose(&b1).unwrap();
let tr_sol2 = lup.solve_transpose(&b2).unwrap();
relative_eq!(m * sol1, b1, epsilon = 1.0e-7) &&
relative_eq!(m * sol2, b2, epsilon = 1.0e-7) &&
relative_eq!(m.transpose() * tr_sol1, b1, epsilon = 1.0e-7) &&
relative_eq!(m.transpose() * tr_sol2, b2, epsilon = 1.0e-7)
prop_assert!(relative_eq!(m * sol1, b1, epsilon = 1.0e-7));
prop_assert!(relative_eq!(m * sol2, b2, epsilon = 1.0e-7));
prop_assert!(relative_eq!(m.transpose() * tr_sol1, b1, epsilon = 1.0e-7));
prop_assert!(relative_eq!(m.transpose() * tr_sol2, b2, epsilon = 1.0e-7));
}
fn lu_inverse(n: usize) -> bool {
if n != 0 {
#[test]
fn lu_inverse(n in PROPTEST_MATRIX_DIM) {
let n = cmp::min(n, 25); // To avoid slowing down the test too much.
let m = DMatrix::<f64>::new_random(n, n);
@ -86,22 +83,17 @@ quickcheck! {
let id1 = &m * &m1;
let id2 = &m1 * &m;
return id1.is_identity(1.0e-7) && id2.is_identity(1.0e-7);
prop_assert!(id1.is_identity(1.0e-7) && id2.is_identity(1.0e-7));
}
}
return true;
}
fn lu_inverse_static(m: Matrix4<f64>) -> bool {
match LU::new(m.clone()).inverse() {
Some(m1) => {
#[test]
fn lu_inverse_static(m in matrix4()) {
if let Some(m1) = LU::new(m.clone()).inverse() {
let id1 = &m * &m1;
let id2 = &m1 * &m;
id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5)
},
None => true
prop_assert!(id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5))
}
}
}

View File

@ -1,20 +1,24 @@
use na::{DMatrix, Matrix4x3};
use nl::QR;
quickcheck! {
fn qr(m: DMatrix<f64>) -> bool {
use crate::proptest::*;
use proptest::{prop_assert, proptest};
proptest! {
#[test]
fn qr(m in dmatrix()) {
let qr = QR::new(m.clone());
let q = qr.q();
let r = qr.r();
relative_eq!(m, q * r, epsilon = 1.0e-7)
prop_assert!(relative_eq!(m, q * r, epsilon = 1.0e-7))
}
fn qr_static(m: Matrix4x3<f64>) -> bool {
#[test]
fn qr_static(m in matrix5x3()) {
let qr = QR::new(m);
let q = qr.q();
let r = qr.r();
relative_eq!(m, q * r, epsilon = 1.0e-7)
prop_assert!(relative_eq!(m, q * r, epsilon = 1.0e-7))
}
}

View File

@ -3,14 +3,16 @@ use std::cmp;
use na::{DMatrix, Matrix4};
use nl::Eigen;
quickcheck! {
fn eigensystem(n: usize) -> bool {
if n != 0 {
use crate::proptest::*;
use proptest::{prop_assert, proptest};
proptest! {
#[test]
fn eigensystem(n in PROPTEST_MATRIX_DIM) {
let n = cmp::min(n, 25);
let m = DMatrix::<f64>::new_random(n, n);
match Eigen::new(m.clone(), true, true) {
Some(eig) => {
if let Some(eig) = Eigen::new(m.clone(), true, true) {
let eigvals = DMatrix::from_diagonal(&eig.eigenvalues);
let transformed_eigvectors = &m * eig.eigenvectors.as_ref().unwrap();
let scaled_eigvectors = eig.eigenvectors.as_ref().unwrap() * &eigvals;
@ -18,20 +20,14 @@ quickcheck! {
let transformed_left_eigvectors = m.transpose() * eig.left_eigenvectors.as_ref().unwrap();
let scaled_left_eigvectors = eig.left_eigenvectors.as_ref().unwrap() * &eigvals;
relative_eq!(transformed_eigvectors, scaled_eigvectors, epsilon = 1.0e-7) &&
relative_eq!(transformed_left_eigvectors, scaled_left_eigvectors, epsilon = 1.0e-7)
},
None => true
}
}
else {
true
prop_assert!(relative_eq!(transformed_eigvectors, scaled_eigvectors, epsilon = 1.0e-7));
prop_assert!(relative_eq!(transformed_left_eigvectors, scaled_left_eigvectors, epsilon = 1.0e-7));
}
}
fn eigensystem_static(m: Matrix4<f64>) -> bool {
match Eigen::new(m, true, true) {
Some(eig) => {
#[test]
fn eigensystem_static(m in matrix4()) {
if let Some(eig) = Eigen::new(m, true, true) {
let eigvals = Matrix4::from_diagonal(&eig.eigenvalues);
let transformed_eigvectors = m * eig.eigenvectors.unwrap();
let scaled_eigvectors = eig.eigenvectors.unwrap() * eigvals;
@ -39,10 +35,8 @@ quickcheck! {
let transformed_left_eigvectors = m.transpose() * eig.left_eigenvectors.unwrap();
let scaled_left_eigvectors = eig.left_eigenvectors.unwrap() * eigvals;
relative_eq!(transformed_eigvectors, scaled_eigvectors, epsilon = 1.0e-7) &&
relative_eq!(transformed_left_eigvectors, scaled_left_eigvectors, epsilon = 1.0e-7)
},
None => true
prop_assert!(relative_eq!(transformed_eigvectors, scaled_eigvectors, epsilon = 1.0e-7));
prop_assert!(relative_eq!(transformed_left_eigvectors, scaled_left_eigvectors, epsilon = 1.0e-7));
}
}
}

View File

@ -1,20 +1,24 @@
use na::{DMatrix, Matrix4};
use na::DMatrix;
use nl::Schur;
use std::cmp;
quickcheck! {
fn schur(n: usize) -> bool {
use crate::proptest::*;
use proptest::{prop_assert, proptest};
proptest! {
#[test]
fn schur(n in PROPTEST_MATRIX_DIM) {
let n = cmp::max(1, cmp::min(n, 10));
let m = DMatrix::<f64>::new_random(n, n);
let (vecs, vals) = Schur::new(m.clone()).unpack();
relative_eq!(&vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7)
prop_assert!(relative_eq!(&vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7))
}
fn schur_static(m: Matrix4<f64>) -> bool {
#[test]
fn schur_static(m in matrix4()) {
let (vecs, vals) = Schur::new(m.clone()).unpack();
relative_eq!(vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7)
prop_assert!(relative_eq!(vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7))
}
}

View File

@ -1,57 +1,53 @@
use na::{DMatrix, Matrix3x4};
use na::{DMatrix, Matrix3x5};
use nl::SVD;
quickcheck! {
fn svd(m: DMatrix<f64>) -> bool {
if m.nrows() != 0 && m.ncols() != 0 {
use crate::proptest::*;
use proptest::{prop_assert, proptest};
proptest! {
#[test]
fn svd(m in dmatrix()) {
let svd = SVD::new(m.clone()).unwrap();
let sm = DMatrix::from_partial_diagonal(m.nrows(), m.ncols(), svd.singular_values.as_slice());
let reconstructed_m = &svd.u * sm * &svd.vt;
let reconstructed_m2 = svd.recompose();
relative_eq!(reconstructed_m, m, epsilon = 1.0e-7) &&
relative_eq!(reconstructed_m2, reconstructed_m, epsilon = 1.0e-7)
}
else {
true
}
prop_assert!(relative_eq!(reconstructed_m, m, epsilon = 1.0e-7));
prop_assert!(relative_eq!(reconstructed_m2, reconstructed_m, epsilon = 1.0e-7));
}
fn svd_static(m: Matrix3x4<f64>) -> bool {
#[test]
fn svd_static(m in matrix3x5()) {
let svd = SVD::new(m).unwrap();
let sm = Matrix3x4::from_partial_diagonal(svd.singular_values.as_slice());
let sm = Matrix3x5::from_partial_diagonal(svd.singular_values.as_slice());
let reconstructed_m = &svd.u * &sm * &svd.vt;
let reconstructed_m2 = svd.recompose();
relative_eq!(reconstructed_m, m, epsilon = 1.0e-7) &&
relative_eq!(reconstructed_m2, m, epsilon = 1.0e-7)
}
fn pseudo_inverse(m: DMatrix<f64>) -> bool {
if m.nrows() == 0 || m.ncols() == 0 {
return true;
prop_assert!(relative_eq!(reconstructed_m, m, epsilon = 1.0e-7));
prop_assert!(relative_eq!(reconstructed_m2, m, epsilon = 1.0e-7));
}
#[test]
fn pseudo_inverse(m in dmatrix()) {
let svd = SVD::new(m.clone()).unwrap();
let im = svd.pseudo_inverse(1.0e-7);
if m.nrows() <= m.ncols() {
return (&m * &im).is_identity(1.0e-7)
prop_assert!((&m * &im).is_identity(1.0e-7));
}
if m.nrows() >= m.ncols() {
return (im * m).is_identity(1.0e-7)
prop_assert!((im * m).is_identity(1.0e-7));
}
}
return true;
}
fn pseudo_inverse_static(m: Matrix3x4<f64>) -> bool {
#[test]
fn pseudo_inverse_static(m in matrix3x5()) {
let svd = SVD::new(m).unwrap();
let im = svd.pseudo_inverse(1.0e-7);
(m * im).is_identity(1.0e-7)
prop_assert!((m * im).is_identity(1.0e-7))
}
}

View File

@ -1,20 +1,25 @@
use std::cmp;
use na::{DMatrix, Matrix4};
use na::DMatrix;
use nl::SymmetricEigen;
quickcheck! {
fn symmetric_eigen(n: usize) -> bool {
use crate::proptest::*;
use proptest::{prop_assert, proptest};
proptest! {
#[test]
fn symmetric_eigen(n in PROPTEST_MATRIX_DIM) {
let n = cmp::max(1, cmp::min(n, 10));
let m = DMatrix::<f64>::new_random(n, n);
let eig = SymmetricEigen::new(m.clone());
let recomp = eig.recompose();
relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)
prop_assert!(relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5))
}
fn symmetric_eigen_static(m: Matrix4<f64>) -> bool {
#[test]
fn symmetric_eigen_static(m in matrix4()) {
let eig = SymmetricEigen::new(m);
let recomp = eig.recompose();
relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)
prop_assert!(relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5))
}
}

View File

@ -1,4 +1,7 @@
//! Tests for proptest-related functionality.
#![allow(dead_code)]
use nalgebra::allocator::Allocator;
use nalgebra::base::dimension::*;
use nalgebra::proptest::{DimRange, MatrixStrategy};
@ -118,12 +121,10 @@ where
macro_rules! define_strategies(
($($strategy_: ident $strategy: ident<$nrows: ident, $ncols: ident>),*) => {$(
#[allow(dead_code)]
pub fn $strategy() -> impl Strategy<Value = MatrixMN<f64, $nrows, $ncols>> {
matrix(PROPTEST_F64, $nrows, $ncols)
}
#[allow(dead_code)]
pub fn $strategy_<ScalarStrategy>(scalar_strategy: ScalarStrategy) -> impl Strategy<Value = MatrixMN<ScalarStrategy::Value, $nrows, $ncols>>
where
ScalarStrategy: Strategy + Clone + 'static,