Add determinant methods in solve module

Author: jturner314

src/lapack_traits/solve.rs

@@ -10,6 +10,14 @@ use super::{Pivot, Transpose, into_result};
 
 /// Wraps `*getrf`, `*getri`, and `*getrs`
 pub trait Solve_: Sized {
+    /// Computes the LU factorization of a general `m x n` matrix `a` using
+    /// partial pivoting with row interchanges.
+    ///
+    /// If the result matches `Err(LinalgError::Lapack(LapackError {
+    /// return_code )) if return_code > 0`, then `U[(return_code,
+    /// return_code)]` is exactly zero. The factorization has been completed,
+    /// but the factor `U` is exactly singular, and division by zero will occur
+    /// if it is used to solve a system of equations.
     unsafe fn lu(MatrixLayout, a: &mut [Self]) -> Result<Pivot>;
     unsafe fn inv(MatrixLayout, a: &mut [Self], &Pivot) -> Result<()>;
     unsafe fn solve(MatrixLayout, Transpose, a: &[Self], &Pivot, b: &mut [Self]) -> Result<()>;

src/solve.rs

@@ -328,3 +328,100 @@ where
         f.inv_into()
     }
 }
+
+/// An interface for calculating determinants of matrix refs.
+pub trait Determinant<A: Scalar> {
+    /// Computes the determinant of the matrix.
+    fn det(&self) -> Result<A>;
+}
+
+/// An interface for calculating determinants of matrices.
+pub trait DeterminantInto<A: Scalar> {
+    /// Computes the determinant of the matrix.
+    fn det_into(self) -> Result<A>;
+}
+
+fn lu_det<'a, A, P, U>(ipiv_iter: P, u_diag_iter: U) -> A
+where
+    A: Scalar,
+    P: Iterator<Item = i32>,
+    U: Iterator<Item = &'a A>,
+{
+    let pivot_sign = if ipiv_iter
+        .enumerate()
+        .filter(|&(i, pivot)| pivot != i as i32 + 1)
+        .count() % 2 == 0
+    {
+        A::one()
+    } else {
+        -A::one()
+    };
+    let (upper_sign, ln_det) = u_diag_iter.fold((A::one(), A::zero()), |(upper_sign, ln_det), &elem| {
+        let abs_elem = elem.abs();
+        (
+            upper_sign * elem.div_real(abs_elem),
+            ln_det.add_real(abs_elem.ln()),
+        )
+    });
+    pivot_sign * upper_sign * ln_det.exp()
+}
+
+fn check_square<S: Data>(a: &ArrayBase<S, Ix2>) -> Result<()> {
+    if a.is_square() {
+        Ok(())
+    } else {
+        Err(NotSquareError::new(a.rows() as i32, a.cols() as i32).into())
+    }
+}
+
+impl<A, S> Determinant<A> for LUFactorized<S>
+where
+    A: Scalar,
+    S: Data<Elem = A>,
+{
+    fn det(&self) -> Result<A> {
+        check_square(&self.a)?;
+        Ok(lu_det(self.ipiv.iter().cloned(), self.a.diag().iter()))
+    }
+}
+
+impl<A, S> DeterminantInto<A> for LUFactorized<S>
+where
+    A: Scalar,
+    S: Data<Elem = A>,
+{
+    fn det_into(self) -> Result<A> {
+        check_square(&self.a)?;
+        Ok(lu_det(self.ipiv.into_iter(), self.a.into_diag().iter()))
+    }
+}
+
+impl<A, S> Determinant<A> for ArrayBase<S, Ix2>
+where
+    A: Scalar,
+    S: Data<Elem = A>,
+{
+    fn det(&self) -> Result<A> {
+        check_square(&self)?;
+        match self.factorize() {
+            Ok(fac) => fac.det(),
+            Err(LinalgError::Lapack(LapackError { return_code })) if return_code > 0 => Ok(A::zero()),
+            Err(err) => Err(err),
+        }
+    }
+}
+
+impl<A, S> DeterminantInto<A> for ArrayBase<S, Ix2>
+where
+    A: Scalar,
+    S: DataMut<Elem = A>,
+{
+    fn det_into(self) -> Result<A> {
+        check_square(&self)?;
+        match self.factorize_into() {
+            Ok(fac) => fac.det_into(),
+            Err(LinalgError::Lapack(LapackError { return_code })) if return_code > 0 => Ok(A::zero()),
+            Err(err) => Err(err),
+        }
+    }
+}

src/types.rs

@@ -22,6 +22,7 @@ pub use num_complex::Complex64 as c64;
 /// - [abs_sqr](trait.Absolute.html#tymethod.abs_sqr)
 /// - [sqrt](trait.SquareRoot.html#tymethod.sqrt)
 /// - [exp](trait.Exponential.html#tymethod.exp)
+/// - [ln](trait.NaturalLogarithm.html#tymethod.ln)
 /// - [conj](trait.Conjugate.html#tymethod.conj)
 /// - [randn](trait.RandNormal.html#tymethod.randn)
 ///
@@ -33,6 +34,7 @@ pub trait Scalar
     + Absolute
     + SquareRoot
     + Exponential
+    + NaturalLogarithm
     + Conjugate
     + RandNormal
     + Neg<Output = Self>
@@ -118,6 +120,11 @@ pub trait Exponential {
     fn exp(&self) -> Self;
 }
 
+/// Define `ln()` more generally
+pub trait NaturalLogarithm {
+    fn ln(&self) -> Self;
+}
+
 /// Complex conjugate value
 pub trait Conjugate: Copy {
     fn conj(self) -> Self;
@@ -207,6 +214,18 @@ impl Exponential for $complex {
     }
 }
 
+impl NaturalLogarithm for $real {
+    fn ln(&self) -> Self {
+        Float::ln(*self)
+    }
+}
+
+impl NaturalLogarithm for $complex {
+    fn ln(&self) -> Self {
+        Complex::ln(self)
+    }
+}
+
 impl Conjugate for $real {
     fn conj(self) -> Self {
         self

tests/solve.rs

@@ -0,0 +1,89 @@
+extern crate ndarray;
+#[macro_use]
+extern crate ndarray_linalg;
+extern crate num_traits;
+
+use ndarray::*;
+use ndarray_linalg::*;
+use num_traits::Zero;
+
+fn det_3x3<A, S>(a: ArrayBase<S, Ix2>) -> A
+where
+    A: Scalar,
+    S: Data<Elem = A>,
+{
+    a[(0, 0)] * a[(1, 1)] * a[(2, 2)] + a[(0, 1)] * a[(1, 2)] * a[(2, 0)] + a[(0, 2)] * a[(1, 0)] * a[(2, 1)] -
+        a[(0, 2)] * a[(1, 1)] * a[(2, 0)] - a[(0, 1)] * a[(1, 0)] * a[(2, 2)] - a[(0, 0)] * a[(1, 2)] * a[(2, 1)]
+}
+
+#[test]
+fn det_zero() {
+    macro_rules! det_zero {
+        ($elem:ty) => {
+            let a: Array2<$elem> = array![[Zero::zero()]];
+            assert_eq!(a.det().unwrap(), Zero::zero());
+            assert_eq!(a.det_into().unwrap(), Zero::zero());
+        }
+    }
+    det_zero!(f64);
+    det_zero!(f32);
+    det_zero!(c64);
+    det_zero!(c32);
+}
+
+#[test]
+fn det_zero_nonsquare() {
+    macro_rules! det_zero_nonsquare {
+        ($elem:ty, $shape:expr) => {
+            let a: Array2<$elem> = Array2::zeros($shape);
+            assert!(a.det().is_err());
+            assert!(a.det_into().is_err());
+        }
+    }
+    for &shape in &[(1, 2).into_shape(), (1, 2).f()] {
+        det_zero_nonsquare!(f64, shape);
+        det_zero_nonsquare!(f32, shape);
+        det_zero_nonsquare!(c64, shape);
+        det_zero_nonsquare!(c32, shape);
+    }
+}
+
+#[test]
+fn det() {
+    macro_rules! det {
+        ($elem:ty, $shape:expr, $rtol:expr) => {
+            let a: Array2<$elem> = random($shape);
+            println!("a = \n{:?}", a);
+            let det = det_3x3(a.view());
+            assert_rclose!(a.factorize().unwrap().det().unwrap(), det, $rtol);
+            assert_rclose!(a.factorize().unwrap().det_into().unwrap(), det, $rtol);
+            assert_rclose!(a.det().unwrap(), det, $rtol);
+            assert_rclose!(a.det_into().unwrap(), det, $rtol);
+        }
+    }
+    for &shape in &[(3, 3).into_shape(), (3, 3).f()] {
+        det!(f64, shape, 1e-9);
+        det!(f32, shape, 1e-4);
+        det!(c64, shape, 1e-9);
+        det!(c32, shape, 1e-4);
+    }
+}
+
+#[test]
+fn det_nonsquare() {
+    macro_rules! det_nonsquare {
+        ($elem:ty, $shape:expr) => {
+            let a: Array2<$elem> = random($shape);
+            assert!(a.factorize().unwrap().det().is_err());
+            assert!(a.factorize().unwrap().det_into().is_err());
+            assert!(a.det().is_err());
+            assert!(a.det_into().is_err());
+        }
+    }
+    for &shape in &[(1, 2).into_shape(), (1, 2).f(), (2, 1).into_shape(), (2, 1).f()] {
+        det_nonsquare!(f64, shape);
+        det_nonsquare!(f32, shape);
+        det_nonsquare!(c64, shape);
+        det_nonsquare!(c32, shape);
+    }
+}