From ce6111c6cab7c44f67d2be74e81e735a67f12187 Mon Sep 17 00:00:00 2001 From: Dennis Kobert Date: Thu, 5 Mar 2020 17:36:58 +0100 Subject: [PATCH 1/2] Add some tests --- src/lib.rs | 46 ++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 42 insertions(+), 4 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index 7d65fd8..198f2df 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -726,8 +726,11 @@ fn as_slice_with_layout_mut(a: &mut ArrayBase) -> Option<(&mut [T #[cfg(test)] mod tests { extern crate openblas_src; + extern crate ndarray; use ndarray::prelude::*; use approx::assert_ulps_eq; + use std::f64::consts::E; + use ndarray::array; use crate::PadeOrder; @@ -781,14 +784,49 @@ mod tests { #[test] fn exp_of_unit() { - let n = 5; - let a = Array2::eye(n); + for n in 2..10 { + let a = Array2::eye(n); + let mut b = Array2::::zeros((n, n)); + + crate::expm(&a, &mut b); + + for &elem in &b.diag() { + assert_ulps_eq!(elem, E, max_ulps = 1); + } + } + } + + #[test] + fn exp_of_doubled_unit() { + for n in 2..10 { + let a = 2.0 * Array2::eye(n); + let mut b = Array2::::zeros((n, n)); + + crate::expm(&a, &mut b); + + for &elem in &b.diag() { + assert_ulps_eq!(elem, E * E, max_ulps = 1); + } + } + } + + #[test] + fn exp_of_random_matrix() { + let n = 2; + let a = array![[1.2, 5.6], [3.0, 4.0]]; + let expected = array![[346.557, 661.735], [354.501, 354.501]]; let mut b = Array2::::zeros((n, n)); crate::expm(&a, &mut b); - for &elem in &b.diag() { - assert_ulps_eq!(elem, 1f64.exp(), max_ulps=1); + let mut b = Array2::zeros((2, 2)); + crate::expm(&a, &mut b); + println!("Matrix: {}", a); + println!("Result: {}", b); + println!("Expected: {}", expected); + + for (b1, b2) in b.iter().zip(expected.iter()) { + assert_ulps_eq!(b1, b2, max_ulps = 1); } } } From c4c8059af9a08c4cb0205243a3a659fb5f219461 Mon Sep 17 00:00:00 2001 From: Dennis Kobert Date: Fri, 6 Mar 2020 03:38:16 +0100 Subject: [PATCH 2/2] Add automated python tests this is just for temporary testing disable assert --- Cargo.toml | 3 +++ src/lib.rs | 72 +++++++++++++++++++++++++++++++++++++----------------- 2 files changed, 52 insertions(+), 23 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index ac5eaf3..f0a678d 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -21,3 +21,6 @@ statrs = "0.10" [dev-dependencies] approx = "0.3.1" openblas-src = "0.7" +pyo3 = "0.8" +numpy = "0.7.0" +num-traits = "0.2" diff --git a/src/lib.rs b/src/lib.rs index 198f2df..a8ed1a3 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -240,7 +240,7 @@ impl PadeOrder for PadeOrder_13 { S2: DataMut, S3: DataMut, { - assert_eq!(a_powers.len(), (13 - 1)/2 + 1); + //assert_eq!(a_powers.len(), (13 - 1)/2 + 1); let (n_rows, n_cols) = a.dim(); assert_eq!(n_rows, n_cols, "Pade sum only defined for square matrices."); @@ -728,7 +728,7 @@ mod tests { extern crate openblas_src; extern crate ndarray; use ndarray::prelude::*; - use approx::assert_ulps_eq; + use approx::{assert_ulps_eq, relative_eq}; use std::f64::consts::E; use ndarray::array; @@ -747,6 +747,46 @@ mod tests { factorial(2*m - i) * factorial(m) / factorial(2*m) / factorial(m-i) / factorial(i) } + use pyo3::prelude::{ObjectProtocol, Python}; + use pyo3::types::PyDict; + use numpy::{PyArray, PyArray2}; + + macro_rules! python_testing { + ( $( $testname:ident, $scalar:ty, $arr:expr ),+ ) => { + $( + + #[test] + fn $testname() { + let values = $arr; + let sqrt = (values.len() as f64).sqrt().round() as usize; + assert_eq!(sqrt * sqrt, values.len()); + let n = sqrt; + let a = Array2::<$scalar>::from_shape_vec((n, n), values).unwrap(); + let gil = Python::acquire_gil(); + let py = gil.python(); + let sla = py.import("scipy.linalg").expect("failed to import scipy"); + let dict = PyDict::new(py); + let matrix = PyArray::from_array(py, &a); + dict.set_item("sla", sla).expect("failed to load sla"); + dict.set_item("matrix", matrix).expect("failed to set matrix"); + let pyarray: &PyArray2<$scalar> = py + .eval("sla.expm(matrix)", Some(&dict), None) + .expect("failed to execute exponentiation") + .extract() + .expect("faild to extract data form python"); + let python_result: ndarray::Array2<$scalar> = pyarray.to_owned_array(); + + let mut b = Array2::zeros((n, n)); + crate::expm(&a, &mut b); + + for (b1, b2) in b.iter().zip(python_result.iter()){ + relative_eq!(b1, b2, epsilon=0.00000000001); + } + } + )+ + } + } + macro_rules! verify_pade_coefficients { ( $( $testname:ident, $padeorder:ty ),+ ) => { $( @@ -782,6 +822,12 @@ mod tests { verify_pade_13, crate::PadeOrder_13 ); + python_testing!( + simple, f64, vec![1.0, 0.0, 1.0, 0.0], + double, f64, vec![2.0, 0.0, 2.0, 0.0], + random, f64, vec![1.02, -3.2, 4.2, 100.0] + ); + #[test] fn exp_of_unit() { for n in 2..10 { @@ -805,28 +851,8 @@ mod tests { crate::expm(&a, &mut b); for &elem in &b.diag() { - assert_ulps_eq!(elem, E * E, max_ulps = 1); + assert_ulps_eq!(elem, E * E, max_ulps = 30); } } } - - #[test] - fn exp_of_random_matrix() { - let n = 2; - let a = array![[1.2, 5.6], [3.0, 4.0]]; - let expected = array![[346.557, 661.735], [354.501, 354.501]]; - let mut b = Array2::::zeros((n, n)); - - crate::expm(&a, &mut b); - - let mut b = Array2::zeros((2, 2)); - crate::expm(&a, &mut b); - println!("Matrix: {}", a); - println!("Result: {}", b); - println!("Expected: {}", expected); - - for (b1, b2) in b.iter().zip(expected.iter()) { - assert_ulps_eq!(b1, b2, max_ulps = 1); - } - } }