Merge branch 'determinant'

Author: termoshtt

Cargo.toml

@@ -20,3 +20,4 @@ features = ["blas"]
 [dev-dependencies]
 rand = "0.3.14"
 ndarray-rand = "0.3"
+float-cmp = "0.2.3"

src/hermite.rs

@@ -3,6 +3,7 @@
 use ndarray::{Ix2, Array, LinalgScalar};
 use std::fmt::Debug;
 use num_traits::float::Float;
+use num_traits::One;
 use lapack::c::Layout;
 
 use matrix::Matrix;
@@ -23,10 +24,12 @@ pub trait HermiteMatrix: SquareMatrix + Matrix {
     fn ssqrt(self) -> Result<Self, LinalgError>;
     /// Cholesky factorization
     fn cholesky(self) -> Result<Self, LinalgError>;
+    /// calc determinant using Cholesky factorization
+    fn deth(self) -> Result<Self::Scalar, LinalgError>;
 }
 
 impl<A> HermiteMatrix for Array<A, Ix2>
-    where A: ImplQR + ImplSVD + ImplNorm + ImplSolve + ImplEigh + ImplCholesky + LinalgScalar + Float + Debug
+    where A: ImplQR + ImplSVD + ImplNorm + ImplSolve + ImplEigh + ImplCholesky + LinalgScalar + Float + Debug + One
 {
     fn eigh(self) -> Result<(Self::Vector, Self), LinalgError> {
         self.check_square()?;
@@ -67,4 +70,10 @@ impl<A> HermiteMatrix for Array<A, Ix2>
         }
         Ok(c)
     }
+    fn deth(self) -> Result<Self::Scalar, LinalgError> {
+        let (n, _) = self.size();
+        let c = self.cholesky()?;
+        let rt = (0..n).map(|i| c[(i, i)]).fold(A::one(), |det, c| det * c);
+        Ok(rt * rt)
+    }
 }

tests/det.rs

@@ -0,0 +1,36 @@
+
+extern crate ndarray;
+extern crate ndarray_linalg;
+extern crate ndarray_rand;
+extern crate rand;
+extern crate float_cmp;
+
+use ndarray::prelude::*;
+use ndarray_linalg::prelude::*;
+use rand::distributions::*;
+use ndarray_rand::RandomExt;
+use float_cmp::ApproxEqRatio;
+
+fn approx_eq(val: f64, truth: f64, ratio: f64) {
+    if !val.approx_eq_ratio(&truth, ratio) {
+        panic!("Not almost equal! val={:?}, truth={:?}, ratio={:?}",
+               val,
+               truth,
+               ratio);
+    }
+}
+
+fn random_hermite(n: usize) -> Array<f64, Ix2> {
+    let r_dist = Range::new(0., 1.);
+    let a = Array::<f64, _>::random((n, n), r_dist);
+    a.dot(&a.t())
+}
+
+#[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);
+    approx_eq(deth, det_eig, 1.0e-7);
+}