Merge branch 'rcarray'

Author: termoshtt

Cargo.toml

@@ -14,7 +14,7 @@ lapack = "0.11.1"
 num-traits = "0.1.36"
 
 [dependencies.ndarray]
-version = "0.7"
+version = "0.7.3"
 features = ["blas"]
 
 [dev-dependencies]

src/hermite.rs

@@ -1,6 +1,6 @@
 //! Define trait for Hermite matrices
 
-use ndarray::{Ix2, Array};
+use ndarray::{Ix2, Array, RcArray};
 use lapack::c::Layout;
 
 use matrix::{Matrix, MFloat};
@@ -71,3 +71,21 @@ impl<A: HMFloat> HermiteMatrix for Array<A, Ix2> {
         Ok(rt * rt)
     }
 }
+
+impl<A: HMFloat> HermiteMatrix for RcArray<A, Ix2> {
+    fn eigh(self) -> Result<(Self::Vector, Self), LinalgError> {
+        let (e, v) = self.into_owned().eigh()?;
+        Ok((e.into_shared(), v.into_shared()))
+    }
+    fn ssqrt(self) -> Result<Self, LinalgError> {
+        let s = self.into_owned().ssqrt()?;
+        Ok(s.into_shared())
+    }
+    fn cholesky(self) -> Result<Self, LinalgError> {
+        let s = self.into_owned().cholesky()?;
+        Ok(s.into_shared())
+    }
+    fn deth(self) -> Result<Self::Scalar, LinalgError> {
+        self.into_owned().deth()
+    }
+}

src/matrix.rs

@@ -2,6 +2,7 @@
 
 use std::cmp::min;
 use ndarray::prelude::*;
+use ndarray::DataMut;
 use lapack::c::Layout;
 
 use error::{LinalgError, StrideError};
@@ -43,6 +44,35 @@ pub trait Matrix: Sized {
     }
 }
 
+fn check_layout(strides: &[Ixs]) -> Result<Layout, StrideError> {
+    if min(strides[0], strides[1]) != 1 {
+        return Err(StrideError {
+            s0: strides[0],
+            s1: strides[1],
+        });;
+    }
+    if strides[0] < strides[1] {
+        Ok(Layout::ColumnMajor)
+    } else {
+        Ok(Layout::RowMajor)
+    }
+}
+
+fn permutate<A: NdFloat, S>(mut a: &mut ArrayBase<S, Ix2>, ipiv: &Vec<i32>)
+    where S: DataMut<Elem = A>
+{
+    let m = a.cols();
+    for (i, j_) in ipiv.iter().enumerate().rev() {
+        let j = (j_ - 1) as usize;
+        if i == j {
+            continue;
+        }
+        for k in 0..m {
+            a.swap((i, k), (j, k));
+        }
+    }
+}
+
 impl<A: MFloat> Matrix for Array<A, Ix2> {
     type Scalar = A;
     type Vector = Array<A, Ix1>;
@@ -52,18 +82,7 @@ impl<A: MFloat> Matrix for Array<A, Ix2> {
         (self.rows(), self.cols())
     }
     fn layout(&self) -> Result<Layout, StrideError> {
-        let strides = self.strides();
-        if min(strides[0], strides[1]) != 1 {
-            return Err(StrideError {
-                s0: strides[0],
-                s1: strides[1],
-            });;
-        }
-        if strides[0] < strides[1] {
-            Ok(Layout::ColumnMajor)
-        } else {
-            Ok(Layout::RowMajor)
-        }
+        check_layout(self.strides())
     }
     fn opnorm_1(&self) -> Self::Scalar {
         let (m, n) = self.size();
@@ -159,15 +178,45 @@ impl<A: MFloat> Matrix for Array<A, Ix2> {
         Ok((p, lm, am))
     }
     fn permutate(&mut self, ipiv: &Self::Permutator) {
-        let (_, m) = self.size();
-        for (i, j_) in ipiv.iter().enumerate().rev() {
-            let j = (j_ - 1) as usize;
-            if i == j {
-                continue;
-            }
-            for k in 0..m {
-                self.swap((i, k), (j, k));
-            }
-        }
+        permutate(self, ipiv);
+    }
+}
+
+impl<A: MFloat> Matrix for RcArray<A, Ix2> {
+    type Scalar = A;
+    type Vector = RcArray<A, Ix1>;
+    type Permutator = Vec<i32>;
+    fn size(&self) -> (usize, usize) {
+        (self.rows(), self.cols())
+    }
+    fn layout(&self) -> Result<Layout, StrideError> {
+        check_layout(self.strides())
+    }
+    fn opnorm_1(&self) -> Self::Scalar {
+        // XXX unnecessary clone
+        self.to_owned().opnorm_1()
+    }
+    fn opnorm_i(&self) -> Self::Scalar {
+        // XXX unnecessary clone
+        self.to_owned().opnorm_i()
+    }
+    fn opnorm_f(&self) -> Self::Scalar {
+        // XXX unnecessary clone
+        self.to_owned().opnorm_f()
+    }
+    fn svd(self) -> Result<(Self, Self::Vector, Self), LinalgError> {
+        let (u, s, v) = self.into_owned().svd()?;
+        Ok((u.into_shared(), s.into_shared(), v.into_shared()))
+    }
+    fn qr(self) -> Result<(Self, Self), LinalgError> {
+        let (q, r) = self.into_owned().qr()?;
+        Ok((q.into_shared(), r.into_shared()))
+    }
+    fn lu(self) -> Result<(Self::Permutator, Self, Self), LinalgError> {
+        let (p, l, u) = self.into_owned().lu()?;
+        Ok((p, l.into_shared(), u.into_shared()))
+    }
+    fn permutate(&mut self, ipiv: &Self::Permutator) {
+        permutate(self, ipiv);
     }
 }

src/prelude.rs

@@ -4,3 +4,4 @@ pub use matrix::Matrix;
 pub use square::SquareMatrix;
 pub use hermite::HermiteMatrix;
 pub use triangular::{TriangularMatrix, drop_lower, drop_upper};
+pub use util::{all_close_l1, all_close_l2, all_close_max};

src/square.rs

@@ -1,6 +1,6 @@
 //! Define trait for Hermite matrices
 
-use ndarray::{Ix2, Array};
+use ndarray::{Ix2, Array, RcArray, ArrayBase, Data};
 use lapack::c::Layout;
 
 use matrix::{Matrix, MFloat};
@@ -32,6 +32,13 @@ pub trait SquareMatrix: Matrix {
     }
 }
 
+fn trace<A: MFloat, S>(a: &ArrayBase<S, Ix2>) -> A
+    where S: Data<Elem = A>
+{
+    let n = a.rows();
+    (0..n).fold(A::zero(), |sum, i| sum + a[(i, i)])
+}
+
 impl<A: MFloat> SquareMatrix for Array<A, Ix2> {
     fn inv(self) -> Result<Self, LinalgError> {
         self.check_square()?;
@@ -47,7 +54,18 @@ impl<A: MFloat> SquareMatrix for Array<A, Ix2> {
     }
     fn trace(&self) -> Result<Self::Scalar, LinalgError> {
         self.check_square()?;
-        let (n, _) = self.size();
-        Ok((0..n).fold(A::zero(), |sum, i| sum + self[(i, i)]))
+        Ok(trace(self))
+    }
+}
+
+impl<A: MFloat> SquareMatrix for RcArray<A, Ix2> {
+    fn inv(self) -> Result<Self, LinalgError> {
+        // XXX unnecessary clone (should use into_owned())
+        let i = self.to_owned().inv()?;
+        Ok(i.into_shared())
+    }
+    fn trace(&self) -> Result<Self::Scalar, LinalgError> {
+        self.check_square()?;
+        Ok(trace(self))
     }
 }

src/triangular.rs

@@ -1,5 +1,5 @@
 
-use ndarray::{Ix2, Array};
+use ndarray::{Ix2, Array, RcArray, NdFloat, ArrayBase, DataMut};
 
 use matrix::{Matrix, MFloat};
 use square::SquareMatrix;
@@ -32,19 +32,36 @@ impl<A: MFloat> TriangularMatrix for Array<A, Ix2> {
     }
 }
 
-pub fn drop_upper(mut a: Array<f64, Ix2>) -> Array<f64, Ix2> {
+impl<A: MFloat> TriangularMatrix for RcArray<A, Ix2> {
+    fn solve_upper(&self, b: Self::Vector) -> Result<Self::Vector, LinalgError> {
+        // XXX unnecessary clone
+        let x = self.to_owned().solve_upper(b.to_owned())?;
+        Ok(x.into_shared())
+    }
+    fn solve_lower(&self, b: Self::Vector) -> Result<Self::Vector, LinalgError> {
+        // XXX unnecessary clone
+        let x = self.to_owned().solve_lower(b.to_owned())?;
+        Ok(x.into_shared())
+    }
+}
+
+pub fn drop_upper<A: NdFloat, S>(mut a: ArrayBase<S, Ix2>) -> ArrayBase<S, Ix2>
+    where S: DataMut<Elem = A>
+{
     for ((i, j), val) in a.indexed_iter_mut() {
         if i < j {
-            *val = 0.0;
+            *val = A::zero();
         }
     }
     a
 }
 
-pub fn drop_lower(mut a: Array<f64, Ix2>) -> Array<f64, Ix2> {
+pub fn drop_lower<A: NdFloat, S>(mut a: ArrayBase<S, Ix2>) -> ArrayBase<S, Ix2>
+    where S: DataMut<Elem = A>
+{
     for ((i, j), val) in a.indexed_iter_mut() {
         if i > j {
-            *val = 0.0;
+            *val = A::zero();
         }
     }
     a

tests/cholesky.rs

@@ -1,25 +1,30 @@
 include!("header.rs");
 
-#[test]
-fn cholesky() {
-    let r_dist = RealNormal::new(0., 1.);
-    let mut a = Array::<f64, _>::random((3, 3), r_dist);
-    a = a.dot(&a.t());
-    println!("a = \n{:?}", a);
-    let c = a.clone().cholesky().unwrap();
-    println!("c = \n{:?}", c);
-    println!("cc = \n{:?}", c.t().dot(&c));
-    c.t().dot(&c).assert_allclose_l2(&a, 1e-7);
+macro_rules! impl_test {
+    ($modname:ident, $clone:ident) => {
+mod $modname {
+    use super::random_hermite;
+    use ndarray_linalg::prelude::*;
+    #[test]
+    fn cholesky() {
+        let a = random_hermite(3);
+        println!("a = \n{:?}", a);
+        let c = a.$clone().cholesky().unwrap();
+        println!("c = \n{:?}", c);
+        println!("cc = \n{:?}", c.t().dot(&c));
+        all_close_l2(&c.t().dot(&c), &a, 1e-7).unwrap();
+    }
+    #[test]
+    fn cholesky_t() {
+        let a = random_hermite(3);
+        println!("a = \n{:?}", a);
+        let c = a.$clone().cholesky().unwrap();
+        println!("c = \n{:?}", c);
+        println!("cc = \n{:?}", c.t().dot(&c));
+        all_close_l2(&c.t().dot(&c), &a, 1e-7).unwrap();
+    }
 }
+}} // impl_test
 
-#[test]
-fn cholesky_t() {
-    let r_dist = RealNormal::new(0., 1.);
-    let mut a = Array::<f64, _>::random((3, 3), r_dist);
-    a = a.dot(&a.t()).reversed_axes();
-    println!("a = \n{:?}", a);
-    let c = a.clone().cholesky().unwrap();
-    println!("c = \n{:?}", c);
-    println!("cc = \n{:?}", c.t().dot(&c));
-    c.t().dot(&c).assert_allclose_l2(&a, 1e-7);
-}
+impl_test!(owned, clone);
+impl_test!(shared, to_shared);

tests/det.rs

@@ -1,16 +1,21 @@
 include!("header.rs");
 
-fn random_hermite(n: usize) -> Array<f64, Ix2> {
-    let r_dist = RealNormal::new(0., 1.);
-    let a = Array::<f64, _>::random((n, n), r_dist);
-    a.dot(&a.t())
+macro_rules! impl_test{
+    ($modname:ident, $clone:ident) => {
+mod $modname {
+    use super::random_hermite;
+    use ndarray_linalg::prelude::*;
+    use ndarray_numtest::prelude::*;
+    #[test]
+    fn deth() {
+        let a = random_hermite(3);
+        let (e, _) = a.$clone().eigh().unwrap();
+        let deth = a.$clone().deth().unwrap();
+        let det_eig = e.iter().fold(1.0, |x, y| x * y);
+        deth.assert_close(det_eig, 1.0e-7);
+    }
 }
+}} // impl_test
 
-#[test]
-fn deth() {
-    let a = random_hermite(3);
-    let (e, _) = a.clone().eigh().unwrap();
-    let deth = a.clone().deth().unwrap();
-    let det_eig = e.iter().fold(1.0, |x, y| x * y);
-    deth.assert_close(det_eig, 1.0e-7);
-}
+impl_test!(owned, clone);
+impl_test!(shared, to_shared);

tests/eigh.rs

@@ -1,23 +1,33 @@
 include!("header.rs");
 
-#[test]
-fn eigen_vector_manual() {
-    let a = arr2(&[[3.0, 1.0, 1.0], [1.0, 3.0, 1.0], [1.0, 1.0, 3.0]]);
-    let (e, vecs) = a.clone().eigh().unwrap();
-    assert!(e.all_close(&arr1(&[2.0, 2.0, 5.0]), 1.0e-7));
-    for (i, v) in vecs.axis_iter(Axis(1)).enumerate() {
-        let av = a.dot(&v);
-        let ev = v.mapv(|x| e[i] * x);
-        assert!(av.all_close(&ev, 1.0e-7));
+macro_rules! impl_test {
+    ($modname:ident, $clone:ident) => {
+mod $modname {
+    use ndarray::prelude::*;
+    use ndarray_linalg::prelude::*;
+    use ndarray_numtest::prelude::*;
+    #[test]
+    fn eigen_vector_manual() {
+        let a = arr2(&[[3.0, 1.0, 1.0], [1.0, 3.0, 1.0], [1.0, 1.0, 3.0]]);
+        let (e, vecs) = a.$clone().eigh().unwrap();
+        all_close_l2(&e, &arr1(&[2.0, 2.0, 5.0]), 1.0e-7).unwrap();
+        for (i, v) in vecs.axis_iter(Axis(1)).enumerate() {
+            let av = a.dot(&v);
+            let ev = v.mapv(|x| e[i] * x);
+            all_close_l2(&av, &ev, 1.0e-7).unwrap();
+        }
     }
-}
-
-#[test]
-fn diagonalize() {
-    let a = arr2(&[[3.0, 1.0, 1.0], [1.0, 3.0, 1.0], [1.0, 1.0, 3.0]]);
-    let (e, vecs) = a.clone().eigh().unwrap();
-    let s = vecs.t().dot(&a).dot(&vecs);
-    for i in 0..3 {
-        e[i].assert_close(s[(i, i)], 1e-7);
+    #[test]
+    fn diagonalize() {
+        let a = arr2(&[[3.0, 1.0, 1.0], [1.0, 3.0, 1.0], [1.0, 1.0, 3.0]]);
+        let (e, vecs) = a.$clone().eigh().unwrap();
+        let s = vecs.t().dot(&a).dot(&vecs);
+        for i in 0..3 {
+            e[i].assert_close(s[(i, i)], 1e-7);
+        }
     }
 }
+}} // impl_test
+
+impl_test!(owned, clone);
+impl_test!(shared, to_shared);