nalgebra/tests/matrix.rs

720 lines
20 KiB
Rust
Raw Blame History

This file contains invisible Unicode characters

This file contains invisible Unicode characters that are indistinguishable to humans but may be processed differently by a computer. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

#[cfg(feature = "arbitrary")]
#[macro_use]
extern crate quickcheck;
#[macro_use]
extern crate approx;
extern crate num_traits as num;
extern crate alga;
extern crate nalgebra as na;
use num::{Zero, One};
use std::fmt::Display;
use alga::linear::FiniteDimInnerSpace;
use na::{U8, U15,
DVector, DMatrix,
Vector1, Vector2, Vector3, Vector4, Vector5, Vector6,
RowVector4, RowVector5,
Matrix1, Matrix2, Matrix3, Matrix4, Matrix5, Matrix6,
MatrixNM, Matrix2x3, Matrix3x2, Matrix3x4, Matrix4x3, Matrix2x4, Matrix4x5, Matrix4x6};
#[test]
fn is_column_major() {
let a = Matrix2x3::new(1.0, 2.0, 3.0,
4.0, 5.0, 6.0);
let expected = &[ 1.0, 4.0, 2.0, 5.0, 3.0, 6.0 ];
assert_eq!(a.as_slice(), expected);
let a = Matrix2x3::from_row_slice(&[1.0, 2.0, 3.0,
4.0, 5.0, 6.0]);
assert_eq!(a.as_slice(), expected);
let a = Matrix2x3::from_column_slice(&[1.0, 4.0,
2.0, 5.0,
3.0, 6.0]);
assert_eq!(a.as_slice(), expected);
}
#[test]
fn linear_index() {
let a = Matrix2x3::new(1, 2, 3,
4, 5, 6);
assert_eq!(a[0], 1);
assert_eq!(a[1], 4);
assert_eq!(a[2], 2);
assert_eq!(a[3], 5);
assert_eq!(a[4], 3);
assert_eq!(a[5], 6);
let b = Vector4::new(1, 2, 3, 4);
assert_eq!(b[0], 1);
assert_eq!(b[1], 2);
assert_eq!(b[2], 3);
assert_eq!(b[3], 4);
let c = RowVector4::new(1, 2, 3, 4);
assert_eq!(c[0], 1);
assert_eq!(c[1], 2);
assert_eq!(c[2], 3);
assert_eq!(c[3], 4);
}
#[test]
fn identity() {
let id1 = Matrix3::<f64>::identity();
let id2 = Matrix3x4::new(1.0, 0.0, 0.0, 0.0,
0.0, 1.0, 0.0, 0.0,
0.0, 0.0, 1.0, 0.0);
let id2bis = Matrix3x4::identity();
let id3 = Matrix4x3::new(1.0, 0.0, 0.0,
0.0, 1.0, 0.0,
0.0, 0.0, 1.0,
0.0, 0.0, 0.0);
let id3bis = Matrix4x3::identity();
let not_id1 = Matrix3::identity() * 2.0;
let not_id2 = Matrix3x4::new(1.0, 0.0, 0.0, 0.0,
0.0, 1.0, 0.0, 0.0,
0.0, 0.0, 1.0, 1.0);
let not_id3 = Matrix4x3::new(1.0, 0.0, 0.0,
0.0, 1.0, 0.0,
0.0, 0.0, 1.0,
0.0, 1.0, 0.0);
assert_eq!(id2, id2bis);
assert_eq!(id3, id3bis);
assert!(id1.is_identity(0.0));
assert!(id2.is_identity(0.0));
assert!(id3.is_identity(0.0));
assert!(!not_id1.is_identity(0.0));
assert!(!not_id2.is_identity(0.0));
assert!(!not_id3.is_identity(0.0));
}
#[test]
fn coordinates() {
let a = Matrix3x4::new(11, 12, 13, 14,
21, 22, 23, 24,
31, 32, 33, 34);
assert_eq!(a.m11, 11);
assert_eq!(a.m12, 12);
assert_eq!(a.m13, 13);
assert_eq!(a.m14, 14);
assert_eq!(a.m21, 21);
assert_eq!(a.m22, 22);
assert_eq!(a.m23, 23);
assert_eq!(a.m24, 24);
assert_eq!(a.m31, 31);
assert_eq!(a.m32, 32);
assert_eq!(a.m33, 33);
assert_eq!(a.m34, 34);
}
#[test]
fn from_diagonal() {
let diag = Vector3::new(1, 2, 3);
let expected = Matrix3::new(
1, 0, 0,
0, 2, 0,
0, 0, 3);
let a = Matrix3::from_diagonal(&diag);
assert_eq!(a, expected);
}
#[test]
fn from_rows() {
let rows = &[
RowVector4::new(11, 12, 13, 14),
RowVector4::new(21, 22, 23, 24),
RowVector4::new(31, 32, 33, 34)
];
let expected = Matrix3x4::new(
11, 12, 13, 14,
21, 22, 23, 24,
31, 32, 33, 34);
let a = Matrix3x4::from_rows(rows);
assert_eq!(a, expected);
}
#[test]
fn from_columns() {
let columns = &[
Vector3::new(11, 21, 31),
Vector3::new(12, 22, 32),
Vector3::new(13, 23, 33),
Vector3::new(14, 24, 34)
];
let expected = Matrix3x4::new(
11, 12, 13, 14,
21, 22, 23, 24,
31, 32, 33, 34);
let a = Matrix3x4::from_columns(columns);
assert_eq!(a, expected);
}
#[test]
fn from_columns_dynamic() {
let columns = &[
DVector::from_row_slice(3, &[11, 21, 31]),
DVector::from_row_slice(3, &[12, 22, 32]),
DVector::from_row_slice(3, &[13, 23, 33]),
DVector::from_row_slice(3, &[14, 24, 34])
];
let expected = DMatrix::from_row_slice(3, 4,
&[ 11, 12, 13, 14,
21, 22, 23, 24,
31, 32, 33, 34 ]);
let a = DMatrix::from_columns(columns);
assert_eq!(a, expected);
}
#[test]
#[should_panic]
fn from_too_many_rows() {
let rows = &[
RowVector4::new(11, 12, 13, 14),
RowVector4::new(21, 22, 23, 24),
RowVector4::new(31, 32, 33, 34),
RowVector4::new(31, 32, 33, 34)
];
let _ = Matrix3x4::from_rows(rows);
}
#[test]
#[should_panic]
fn from_not_enough_columns() {
let columns = &[
Vector3::new(11, 21, 31),
Vector3::new(14, 24, 34)
];
let _ = Matrix3x4::from_columns(columns);
}
#[test]
#[should_panic]
fn from_rows_with_different_dimensions() {
let columns = &[
DVector::from_row_slice(3, &[11, 21, 31]),
DVector::from_row_slice(3, &[12, 22, 32, 33])
];
let _ = DMatrix::from_columns(columns);
}
#[test]
fn to_homogeneous() {
let a = Vector3::new(1.0, 2.0, 3.0);
let expected_a = Vector4::new(1.0, 2.0, 3.0, 0.0);
let b = DVector::from_row_slice(3, &[1.0, 2.0, 3.0]);
let expected_b = DVector::from_row_slice(4, &[1.0, 2.0, 3.0, 0.0]);
assert_eq!(a.to_homogeneous(), expected_a);
assert_eq!(b.to_homogeneous(), expected_b);
}
#[test]
fn simple_add() {
let a = Matrix2x3::new(1.0, 2.0, 3.0,
4.0, 5.0, 6.0);
let b = Matrix2x3::new(10.0, 20.0, 30.0,
40.0, 50.0, 60.0);
let c = DMatrix::from_row_slice(2, 3, &[ 10.0, 20.0, 30.0,
40.0, 50.0, 60.0 ]);
let expected = Matrix2x3::new(11.0, 22.0, 33.0,
44.0, 55.0, 66.0);
assert_eq!(expected, &a + &b);
assert_eq!(expected, &a + b);
assert_eq!(expected, a + &b);
assert_eq!(expected, a + b);
// Sum of a static matrix with a dynamic one.
assert_eq!(expected, &a + &c);
assert_eq!(expected, a + &c);
assert_eq!(expected, &c + &a);
assert_eq!(expected, &c + a);
}
#[test]
fn simple_scalar_mul() {
let a = Matrix2x3::new(1.0, 2.0, 3.0,
4.0, 5.0, 6.0);
let expected = Matrix2x3::new(10.0, 20.0, 30.0,
40.0, 50.0, 60.0);
assert_eq!(expected, a * 10.0);
assert_eq!(expected, &a * 10.0);
assert_eq!(expected, 10.0 * a);
assert_eq!(expected, 10.0 * &a);
}
#[test]
fn simple_mul() {
let a = Matrix2x3::new(1.0, 2.0, 3.0,
4.0, 5.0, 6.0);
let b = Matrix3x4::new(10.0, 20.0, 30.0, 40.0,
50.0, 60.0, 70.0, 80.0,
90.0, 100.0, 110.0, 120.0);
let expected = Matrix2x4::new(380.0, 440.0, 500.0, 560.0,
830.0, 980.0, 1130.0, 1280.0);
assert_eq!(expected, &a * &b);
assert_eq!(expected, a * &b);
assert_eq!(expected, &a * b);
assert_eq!(expected, a * b);
}
#[test]
fn simple_scalar_conversion() {
let a = Matrix2x3::new(1.0, 2.0, 3.0,
4.0, 5.0, 6.0);
let expected = Matrix2x3::new(1, 2, 3,
4, 5, 6);
let a_u32: Matrix2x3<u32> = na::try_convert(a).unwrap(); // f32 -> u32
let a_f32: Matrix2x3<f32> = na::convert(a_u32); // u32 -> f32
assert_eq!(a, a_f32);
assert_eq!(expected, a_u32);
}
#[test]
#[should_panic]
fn trace_panic() {
let m = DMatrix::<f32>::new_random(2, 3);
let _ = m.trace();
}
#[test]
fn trace() {
let m = Matrix2::new(1.0, 20.0,
30.0, 4.0);
assert_eq!(m.trace(), 5.0);
}
#[test]
fn simple_transpose() {
let a = Matrix2x3::new(1.0, 2.0, 3.0,
4.0, 5.0, 6.0);
let expected = Matrix3x2::new(1.0, 4.0,
2.0, 5.0,
3.0, 6.0);
assert_eq!(a.transpose(), expected);
}
#[test]
fn simple_transpose_mut() {
let mut a = Matrix3::new(1.0, 2.0, 3.0,
4.0, 5.0, 6.0,
7.0, 8.0, 9.0);
let expected = Matrix3::new(1.0, 4.0, 7.0,
2.0, 5.0, 8.0,
3.0, 6.0, 9.0);
a.transpose_mut();
assert_eq!(a, expected);
}
#[test]
fn vector_index_mut() {
let mut v = Vector3::new(1, 2, 3);
assert_eq!(v[0], 1);
assert_eq!(v[1], 2);
assert_eq!(v[2], 3);
v[0] = 10;
v[1] = 20;
v[2] = 30;
assert_eq!(v, Vector3::new(10, 20, 30));
}
#[test]
fn components_mut() {
let mut m2 = Matrix2::from_element(1.0);
let mut m3 = Matrix3::from_element(1.0);
let mut m4 = Matrix4::from_element(1.0);
let mut m5 = Matrix5::from_element(1.0);
let mut m6 = Matrix6::from_element(1.0);
m2.m11 = 0.0; m2.m12 = 0.0;
m2.m21 = 0.0; m2.m22 = 0.0;
m3.m11 = 0.0; m3.m12 = 0.0; m3.m13 = 0.0;
m3.m21 = 0.0; m3.m22 = 0.0; m3.m23 = 0.0;
m3.m31 = 0.0; m3.m32 = 0.0; m3.m33 = 0.0;
m4.m11 = 0.0; m4.m12 = 0.0; m4.m13 = 0.0; m4.m14 = 0.0;
m4.m21 = 0.0; m4.m22 = 0.0; m4.m23 = 0.0; m4.m24 = 0.0;
m4.m31 = 0.0; m4.m32 = 0.0; m4.m33 = 0.0; m4.m34 = 0.0;
m4.m41 = 0.0; m4.m42 = 0.0; m4.m43 = 0.0; m4.m44 = 0.0;
m5.m11 = 0.0; m5.m12 = 0.0; m5.m13 = 0.0; m5.m14 = 0.0; m5.m15 = 0.0;
m5.m21 = 0.0; m5.m22 = 0.0; m5.m23 = 0.0; m5.m24 = 0.0; m5.m25 = 0.0;
m5.m31 = 0.0; m5.m32 = 0.0; m5.m33 = 0.0; m5.m34 = 0.0; m5.m35 = 0.0;
m5.m41 = 0.0; m5.m42 = 0.0; m5.m43 = 0.0; m5.m44 = 0.0; m5.m45 = 0.0;
m5.m51 = 0.0; m5.m52 = 0.0; m5.m53 = 0.0; m5.m54 = 0.0; m5.m55 = 0.0;
m6.m11 = 0.0; m6.m12 = 0.0; m6.m13 = 0.0; m6.m14 = 0.0; m6.m15 = 0.0; m6.m16 = 0.0;
m6.m21 = 0.0; m6.m22 = 0.0; m6.m23 = 0.0; m6.m24 = 0.0; m6.m25 = 0.0; m6.m26 = 0.0;
m6.m31 = 0.0; m6.m32 = 0.0; m6.m33 = 0.0; m6.m34 = 0.0; m6.m35 = 0.0; m6.m36 = 0.0;
m6.m41 = 0.0; m6.m42 = 0.0; m6.m43 = 0.0; m6.m44 = 0.0; m6.m45 = 0.0; m6.m46 = 0.0;
m6.m51 = 0.0; m6.m52 = 0.0; m6.m53 = 0.0; m6.m54 = 0.0; m6.m55 = 0.0; m6.m56 = 0.0;
m6.m61 = 0.0; m6.m62 = 0.0; m6.m63 = 0.0; m6.m64 = 0.0; m6.m65 = 0.0; m6.m66 = 0.0;
assert!(m2.is_zero());
assert!(m3.is_zero());
assert!(m4.is_zero());
assert!(m5.is_zero());
assert!(m6.is_zero());
let mut v1 = Vector1::from_element(1.0);
let mut v2 = Vector2::from_element(1.0);
let mut v3 = Vector3::from_element(1.0);
let mut v4 = Vector4::from_element(1.0);
let mut v5 = Vector5::from_element(1.0);
let mut v6 = Vector6::from_element(1.0);
v1.x = 0.0;
v2.x = 0.0; v2.y = 0.0;
v3.x = 0.0; v3.y = 0.0; v3.z = 0.0;
v4.x = 0.0; v4.y = 0.0; v4.z = 0.0; v4.w = 0.0;
v5.x = 0.0; v5.y = 0.0; v5.z = 0.0; v5.w = 0.0; v5.a = 0.0;
v6.x = 0.0; v6.y = 0.0; v6.z = 0.0; v6.w = 0.0; v6.a = 0.0; v6.b = 0.0;
assert!(v1.is_zero());
assert!(v2.is_zero());
assert!(v3.is_zero());
assert!(v4.is_zero());
assert!(v5.is_zero());
assert!(v6.is_zero());
// Check that the components order is correct.
m3.m11 = 11.0; m3.m12 = 12.0; m3.m13 = 13.0;
m3.m21 = 21.0; m3.m22 = 22.0; m3.m23 = 23.0;
m3.m31 = 31.0; m3.m32 = 32.0; m3.m33 = 33.0;
let expected_m3 = Matrix3::new(11.0, 12.0, 13.0,
21.0, 22.0, 23.0,
31.0, 32.0, 33.0);
assert_eq!(expected_m3, m3);
}
#[test]
fn kronecker() {
let a = Matrix2x3::new(
11, 12, 13,
21, 22, 23);
let b = Matrix4x5::new(
110, 120, 130, 140, 150,
210, 220, 230, 240, 250,
310, 320, 330, 340, 350,
410, 420, 430, 440, 450);
let expected = MatrixNM::<_, U8, U15>::from_row_slice(&[
1210, 1320, 1430, 1540, 1650, 1320, 1440, 1560, 1680, 1800, 1430, 1560, 1690, 1820, 1950,
2310, 2420, 2530, 2640, 2750, 2520, 2640, 2760, 2880, 3000, 2730, 2860, 2990, 3120, 3250,
3410, 3520, 3630, 3740, 3850, 3720, 3840, 3960, 4080, 4200, 4030, 4160, 4290, 4420, 4550,
4510, 4620, 4730, 4840, 4950, 4920, 5040, 5160, 5280, 5400, 5330, 5460, 5590, 5720, 5850,
2310, 2520, 2730, 2940, 3150, 2420, 2640, 2860, 3080, 3300, 2530, 2760, 2990, 3220, 3450,
4410, 4620, 4830, 5040, 5250, 4620, 4840, 5060, 5280, 5500, 4830, 5060, 5290, 5520, 5750,
6510, 6720, 6930, 7140, 7350, 6820, 7040, 7260, 7480, 7700, 7130, 7360, 7590, 7820, 8050,
8610, 8820, 9030, 9240, 9450, 9020, 9240, 9460, 9680, 9900, 9430, 9660, 9890, 10120, 10350 ]);
let computed = a.kronecker(&b);
assert_eq!(computed, expected);
let a = Vector2::new(1, 2);
let b = Vector3::new(10, 20, 30);
let expected = Vector6::new(10, 20, 30, 20, 40, 60);
assert_eq!(a.kronecker(&b), expected);
let a = Vector2::new(1, 2);
let b = RowVector4::new(10, 20, 30, 40);
let expected = Matrix2x4::new(
10, 20, 30, 40,
20, 40, 60, 80);
assert_eq!(a.kronecker(&b), expected);
}
#[test]
fn set_row_column() {
let a = Matrix4x5::new(
11, 12, 13, 14, 15,
21, 22, 23, 24, 25,
31, 32, 33, 34, 35,
41, 42, 43, 44, 45);
let expected1 = Matrix4x5::new(
11, 12, 13, 14, 15,
42, 43, 44, 45, 46,
31, 32, 33, 34, 35,
41, 42, 43, 44, 45);
let expected2 = Matrix4x5::new(
11, 12, 100, 14, 15,
42, 43, 101, 45, 46,
31, 32, 102, 34, 35,
41, 42, 103, 44, 45);
let row = RowVector5::new(42, 43, 44, 45, 46);
let col = Vector4::new(100, 101, 102, 103);
let mut computed = a;
computed.set_row(1, &row);
assert_eq!(expected1, computed);
computed.set_column(2, &col);
assert_eq!(expected2, computed);
}
#[cfg(feature = "arbitrary")]
quickcheck!{
/*
*
* Transposition.
*
*/
fn transpose_transpose_is_self(m: Matrix2x3<f64>) -> bool {
m.transpose().transpose() == m
}
fn transpose_mut_transpose_mut_is_self(m: Matrix3<f64>) -> bool {
let mut mm = m;
mm.transpose_mut();
mm.transpose_mut();
m == mm
}
fn transpose_transpose_is_id_dyn(m: DMatrix<f64>) -> bool {
m.transpose().transpose() == m
}
fn check_transpose_components_dyn(m: DMatrix<f64>) -> bool {
let tr = m.transpose();
let (nrows, ncols) = m.shape();
if nrows != tr.shape().1 || ncols != tr.shape().0 {
return false
}
for i in 0 .. nrows {
for j in 0 .. ncols {
if m[(i, j)] != tr[(j, i)] {
return false
}
}
}
true
}
fn tr_mul_is_transpose_then_mul(m: Matrix4x6<f64>, v: Vector4<f64>) -> bool {
m.transpose() * v == m.tr_mul(&v)
}
/*
*
*
* Inversion.
*
*
*/
fn self_mul_inv_is_id_dim1(m: Matrix1<f64>) -> bool {
if let Some(im) = m.try_inverse() {
let id = Matrix1::one();
relative_eq!(im * m, id, epsilon = 1.0e-7) &&
relative_eq!(m * im, id, epsilon = 1.0e-7)
}
else {
true
}
}
fn self_mul_inv_is_id_dim2(m: Matrix2<f64>) -> bool {
if let Some(im) = m.try_inverse() {
let id = Matrix2::one();
relative_eq!(im * m, id, epsilon = 1.0e-7) &&
relative_eq!(m * im, id, epsilon = 1.0e-7)
}
else {
true
}
}
fn self_mul_inv_is_id_dim3(m: Matrix3<f64>) -> bool {
if let Some(im) = m.try_inverse() {
let id = Matrix3::one();
relative_eq!(im * m, id, epsilon = 1.0e-7) &&
relative_eq!(m * im, id, epsilon = 1.0e-7)
}
else {
true
}
}
fn self_mul_inv_is_id_dim6(m: Matrix6<f64>) -> bool {
if let Some(im) = m.try_inverse() {
let id = Matrix6::one();
relative_eq!(im * m, id, epsilon = 1.0e-7) &&
relative_eq!(m * im, id, epsilon = 1.0e-7)
}
else {
true
}
}
/*
*
* Normalization.
*
*/
fn normalized_vec_norm_is_one(v: Vector3<f64>) -> bool {
if let Some(nv) = v.try_normalize(1.0e-10) {
relative_eq!(nv.norm(), 1.0, epsilon = 1.0e-7)
}
else {
true
}
}
fn normalized_vec_norm_is_one_dyn(v: DVector<f64>) -> bool {
if let Some(nv) = v.try_normalize(1.0e-10) {
relative_eq!(nv.norm(), 1.0)
}
else {
true
}
}
}
// FIXME: move this to alga ?
macro_rules! finite_dim_inner_space_test(
($($Vector: ident, $orthonormal_subspace: ident, $orthonormalization: ident);* $(;)*) => {$(
#[cfg(feature = "arbitrary")]
quickcheck!{
fn $orthonormal_subspace(vs: Vec<$Vector<f64>>) -> bool {
let mut given_basis = vs.clone();
let given_basis_dim = $Vector::orthonormalize(&mut given_basis[..]);
let mut ortho_basis = Vec::new();
$Vector::orthonormal_subspace_basis(
&given_basis[.. given_basis_dim],
|e| { ortho_basis.push(*e); true }
);
if !is_subspace_basis(&ortho_basis[..]) {
return false;
}
for v in vs {
for b in &ortho_basis {
if !relative_eq!(v.dot(b), 0.0, epsilon = 1.0e-7) {
println!("Found dot product: {} · {} = {}", v, b, v.dot(b));
return false;
}
}
}
true
}
fn $orthonormalization(vs: Vec<$Vector<f64>>) -> bool {
let mut basis = vs.clone();
let subdim = $Vector::orthonormalize(&mut basis[..]);
if !is_subspace_basis(&basis[.. subdim]) {
return false;
}
for mut e in vs {
for b in &basis[.. subdim] {
e -= e.dot(b) * b
}
// Any element of `e` must be a linear combination of the basis elements.
if !relative_eq!(e.norm(), 0.0, epsilon = 1.0e-7) {
println!("Orthonormalization; element decomposition failure: {}", e);
println!("... the non-zero norm is: {}", e.norm());
return false;
}
}
true
}
}
)*}
);
finite_dim_inner_space_test!(
Vector1, orthonormal_subspace_basis1, orthonormalize1;
Vector2, orthonormal_subspace_basis2, orthonormalize2;
Vector3, orthonormal_subspace_basis3, orthonormalize3;
Vector4, orthonormal_subspace_basis4, orthonormalize4;
Vector5, orthonormal_subspace_basis5, orthonormalize5;
Vector6, orthonormal_subspace_basis6, orthonormalize6;
);
/*
*
* Helper functions.
*
*/
fn is_subspace_basis<T: FiniteDimInnerSpace<Real = f64> + Display>(vs: &[T]) -> bool {
for i in 0 .. vs.len() {
// Basis elements must be normalized.
if !relative_eq!(vs[i].norm(), 1.0, epsilon = 1.0e-7) {
println!("Non-zero basis element norm: {}", vs[i].norm());
return false;
}
for j in 0 .. i {
// Basis elements must be orthogonal.
if !relative_eq!(vs[i].dot(&vs[j]), 0.0, epsilon = 1.0e-7) {
println!("Non-orthogonal basis elements: {} · {} = {}", vs[i], vs[j], vs[i].dot(&vs[j]));
return false
}
}
}
true
}