Norm for ArrayBase<S, D>

src/lib.rs

@@ -48,6 +48,8 @@ pub mod triangular;
 pub mod qr;
 pub mod svd;
 pub mod eigh;
-pub mod norm;
+pub mod opnorm;
 pub mod solve;
 pub mod cholesky;
+
+pub mod util;

src/matrix.rs

@@ -7,11 +7,11 @@ use lapack::c::Layout;
 use error::{LinalgError, StrideError};
 use qr::ImplQR;
 use svd::ImplSVD;
-use norm::ImplNorm;
+use opnorm::ImplOpNorm;
 use solve::ImplSolve;
 
-pub trait MFloat: ImplQR + ImplSVD + ImplNorm + ImplSolve + NdFloat {}
-impl<A: ImplQR + ImplSVD + ImplNorm + ImplSolve + NdFloat> MFloat for A {}
+pub trait MFloat: ImplQR + ImplSVD + ImplOpNorm + ImplSolve + NdFloat {}
+impl<A: ImplQR + ImplSVD + ImplOpNorm + ImplSolve + NdFloat> MFloat for A {}
 
 /// Methods for general matrices
 pub trait Matrix: Sized {
@@ -23,11 +23,11 @@ pub trait Matrix: Sized {
     /// Layout (C/Fortran) of matrix
     fn layout(&self) -> Result<Layout, StrideError>;
     /// Operator norm for L-1 norm
-    fn norm_1(&self) -> Self::Scalar;
+    fn opnorm_1(&self) -> Self::Scalar;
     /// Operator norm for L-inf norm
-    fn norm_i(&self) -> Self::Scalar;
+    fn opnorm_i(&self) -> Self::Scalar;
     /// Frobenius norm
-    fn norm_f(&self) -> Self::Scalar;
+    fn opnorm_f(&self) -> Self::Scalar;
     /// singular-value decomposition (SVD)
     fn svd(self) -> Result<(Self, Self::Vector, Self), LinalgError>;
     /// QR decomposition
@@ -65,27 +65,27 @@ impl<A: MFloat> Matrix for Array<A, Ix2> {
             Ok(Layout::RowMajor)
         }
     }
-    fn norm_1(&self) -> Self::Scalar {
+    fn opnorm_1(&self) -> Self::Scalar {
         let (m, n) = self.size();
         let strides = self.strides();
         if strides[0] > strides[1] {
-            ImplNorm::norm_i(n, m, self.clone().into_raw_vec())
+            ImplOpNorm::opnorm_i(n, m, self.clone().into_raw_vec())
         } else {
-            ImplNorm::norm_1(m, n, self.clone().into_raw_vec())
+            ImplOpNorm::opnorm_1(m, n, self.clone().into_raw_vec())
         }
     }
-    fn norm_i(&self) -> Self::Scalar {
+    fn opnorm_i(&self) -> Self::Scalar {
         let (m, n) = self.size();
         let strides = self.strides();
         if strides[0] > strides[1] {
-            ImplNorm::norm_1(n, m, self.clone().into_raw_vec())
+            ImplOpNorm::opnorm_1(n, m, self.clone().into_raw_vec())
         } else {
-            ImplNorm::norm_i(m, n, self.clone().into_raw_vec())
+            ImplOpNorm::opnorm_i(m, n, self.clone().into_raw_vec())
         }
     }
-    fn norm_f(&self) -> Self::Scalar {
+    fn opnorm_f(&self) -> Self::Scalar {
         let (m, n) = self.size();
-        ImplNorm::norm_f(m, n, self.clone().into_raw_vec())
+        ImplOpNorm::opnorm_f(m, n, self.clone().into_raw_vec())
     }
     fn svd(self) -> Result<(Self, Self::Vector, Self), LinalgError> {
         let (n, m) = self.size();

src/opnorm.rs

@@ -0,0 +1,27 @@
+//! Implement Norms for matrices
+
+use lapack::c::*;
+
+pub trait ImplOpNorm: Sized {
+    fn opnorm_1(m: usize, n: usize, a: Vec<Self>) -> Self;
+    fn opnorm_i(m: usize, n: usize, a: Vec<Self>) -> Self;
+    fn opnorm_f(m: usize, n: usize, a: Vec<Self>) -> Self;
+}
+
+macro_rules! impl_opnorm {
+    ($scalar:ty, $lange:path) => {
+impl ImplOpNorm for $scalar {
+    fn opnorm_1(m: usize, n: usize, mut a: Vec<Self>) -> Self {
+        $lange(Layout::ColumnMajor, b'o', m as i32, n as i32, &mut a, m as i32)
+    }
+    fn opnorm_i(m: usize, n: usize, mut a: Vec<Self>) -> Self {
+        $lange(Layout::ColumnMajor, b'i', m as i32, n as i32, &mut a, m as i32)
+    }
+    fn opnorm_f(m: usize, n: usize, mut a: Vec<Self>) -> Self {
+        $lange(Layout::ColumnMajor, b'f', m as i32, n as i32, &mut a, m as i32)
+    }
+}
+}} // end macro_rules
+
+impl_opnorm!(f64, dlange);
+impl_opnorm!(f32, slange);

src/util.rs

@@ -0,0 +1,43 @@
+//! module for topologcal vector space
+//!
+
+use std::iter::Sum;
+use ndarray::{ArrayBase, Data, Dimension, LinalgScalar};
+use num_traits::Float;
+use super::vector::*;
+
+pub fn all_close_max<A, Tol, S1, S2, D>(test: &ArrayBase<S1, D>,
+                                        truth: &ArrayBase<S2, D>,
+                                        atol: Tol)
+                                        -> Result<Tol, Tol>
+    where A: LinalgScalar + Squared<Output = Tol>,
+          Tol: Float + Sum,
+          S1: Data<Elem = A>,
+          S2: Data<Elem = A>,
+          D: Dimension
+{
+    let tol = (test - truth).norm_max();
+    if tol < atol { Ok(tol) } else { Err(tol) }
+}
+
+pub fn all_close_l1<A, Tol, S1, S2, D>(test: &ArrayBase<S1, D>, truth: &ArrayBase<S2, D>, rtol: Tol) -> Result<Tol, Tol>
+    where A: LinalgScalar + Squared<Output = Tol>,
+          Tol: Float + Sum,
+          S1: Data<Elem = A>,
+          S2: Data<Elem = A>,
+          D: Dimension
+{
+    let tol = (test - truth).norm_l1() / truth.norm_l1();
+    if tol < rtol { Ok(tol) } else { Err(tol) }
+}
+
+pub fn all_close_l2<A, Tol, S1, S2, D>(test: &ArrayBase<S1, D>, truth: &ArrayBase<S2, D>, rtol: Tol) -> Result<Tol, Tol>
+    where A: LinalgScalar + Squared<Output = Tol>,
+          Tol: Float + Sum,
+          S1: Data<Elem = A>,
+          S2: Data<Elem = A>,
+          D: Dimension
+{
+    let tol = (test - truth).norm_l2() / truth.norm_l2();
+    if tol < rtol { Ok(tol) } else { Err(tol) }
+}

src/vector.rs

@@ -1,18 +1,57 @@
 //! Define trait for vectors
 
-use ndarray::{LinalgScalar, Array, Ix1};
+use std::iter::Sum;
+use ndarray::{LinalgScalar, ArrayBase, Data, Dimension};
 use num_traits::float::Float;
 
 /// Methods for vectors
 pub trait Vector {
     type Scalar;
+    /// rename of norm_l2
+    fn norm(&self) -> Self::Scalar {
+        self.norm_l2()
+    }
+    /// L-1 norm
+    fn norm_l1(&self) -> Self::Scalar;
     /// L-2 norm
-    fn norm(&self) -> Self::Scalar;
+    fn norm_l2(&self) -> Self::Scalar;
+    /// maximum norm
+    fn norm_max(&self) -> Self::Scalar;
 }
 
-impl<A: Float + LinalgScalar> Vector for Array<A, Ix1> {
-    type Scalar = A;
-    fn norm(&self) -> Self::Scalar {
-        self.dot(&self).sqrt()
+impl<A, S, D, T> Vector for ArrayBase<S, D>
+    where A: LinalgScalar + Squared<Output = T>,
+          T: Float + Sum,
+          S: Data<Elem = A>,
+          D: Dimension
+{
+    type Scalar = T;
+    fn norm_l1(&self) -> Self::Scalar {
+        self.iter().map(|x| x.sq_abs()).sum()
+    }
+    fn norm_l2(&self) -> Self::Scalar {
+        self.iter().map(|x| x.squared()).sum::<T>().sqrt()
+    }
+    fn norm_max(&self) -> Self::Scalar {
+        self.iter().fold(T::zero(), |f, &val| {
+            let v = val.sq_abs();
+            if f > v { f } else { v }
+        })
+    }
+}
+
+pub trait Squared {
+    type Output;
+    fn squared(&self) -> Self::Output;
+    fn sq_abs(&self) -> Self::Output;
+}
+
+impl<A: Float> Squared for A {
+    type Output = A;
+    fn squared(&self) -> A {
+        *self * *self
+    }
+    fn sq_abs(&self) -> A {
+        self.abs()
     }
 }

tests/opnorm.rs

@@ -0,0 +1,55 @@
+include!("header.rs");
+
+#[test]
+fn matrix_opnorm_square() {
+    let a = Array::range(1., 10., 1.).into_shape((3, 3)).unwrap();
+    a.opnorm_1().assert_close(18.0, 1e-7);
+    a.opnorm_i().assert_close(24.0, 1e-7);
+    a.opnorm_f().assert_close(285.0.sqrt(), 1e-7);
+}
+
+#[test]
+fn matrix_opnorm_square_t() {
+    let a = Array::range(1., 10., 1.).into_shape((3, 3)).unwrap().reversed_axes();
+    a.opnorm_1().assert_close(24.0, 1e-7);
+    a.opnorm_i().assert_close(18.0, 1e-7);
+    a.opnorm_f().assert_close(285.0.sqrt(), 1e-7);
+}
+
+#[test]
+fn matrix_opnorm_3x4() {
+    let a = Array::range(1., 13., 1.).into_shape((3, 4)).unwrap();
+    a.opnorm_1().assert_close(24.0, 1e-7);
+    a.opnorm_i().assert_close(42.0, 1e-7);
+    a.opnorm_f().assert_close(650.0.sqrt(), 1e-7);
+}
+
+#[test]
+fn matrix_opnorm_3x4_t() {
+    let a = Array::range(1., 13., 1.)
+        .into_shape((3, 4))
+        .unwrap()
+        .reversed_axes();
+    a.opnorm_1().assert_close(42.0, 1e-7);
+    a.opnorm_i().assert_close(24.0, 1e-7);
+    a.opnorm_f().assert_close(650.0.sqrt(), 1e-7);
+}
+
+#[test]
+fn matrix_opnorm_4x3() {
+    let a = Array::range(1., 13., 1.).into_shape((4, 3)).unwrap();
+    a.opnorm_1().assert_close(30.0, 1e-7);
+    a.opnorm_i().assert_close(33.0, 1e-7);
+    a.opnorm_f().assert_close(650.0.sqrt(), 1e-7);
+}
+
+#[test]
+fn matrix_opnorm_4x3_t() {
+    let a = Array::range(1., 13., 1.)
+        .into_shape((4, 3))
+        .unwrap()
+        .reversed_axes();
+    a.opnorm_1().assert_close(33.0, 1e-7);
+    a.opnorm_i().assert_close(30.0, 1e-7);
+    a.opnorm_f().assert_close(650.0.sqrt(), 1e-7);
+}

tests/vector.rs

@@ -0,0 +1,23 @@
+include!("header.rs");
+
+#[test]
+fn vector_norm() {
+    let a = Array::range(1., 10., 1.);
+    a.norm().assert_close(285.0.sqrt(), 1e-7);
+}
+
+#[test]
+fn vector_norm_l1() {
+    let a = arr1(&[1.0, -1.0]);
+    a.norm_l1().assert_close(2.0, 1e-7);
+    let b = arr2(&[[0.0, -1.0], [1.0, 0.0]]);
+    b.norm_l1().assert_close(2.0, 1e-7);
+}
+
+#[test]
+fn vector_norm_max() {
+    let a = arr1(&[1.0, 1.0, -3.0]);
+    a.norm_max().assert_close(3.0, 1e-7);
+    let b = arr2(&[[1.0, 3.0], [1.0, -4.0]]);
+    b.norm_max().assert_close(4.0, 1e-7);
+}