Generate specific matrices

Cargo.toml

@@ -15,6 +15,7 @@ openblas = ["blas/openblas", "lapack/openblas"]
 netlib   = ["blas/netlib", "lapack/netlib"]
 
 [dependencies]
+rand = "0.3"
 derive-new = "0.4"
 enum-error-derive = "0.1"
 num-traits  = "0.1"

src/assert.rs

@@ -4,6 +4,7 @@ use std::iter::Sum;
 use num_traits::Float;
 use ndarray::*;
 
+use super::types::*;
 use super::vector::*;
 
 pub fn rclose<A, Tol>(test: A, truth: A, rtol: Tol) -> Result<Tol, Tol>

src/generate.rs

@@ -0,0 +1,103 @@
+
+use ndarray::*;
+use std::ops::*;
+use rand::*;
+
+use super::layout::*;
+use super::types::*;
+use super::error::*;
+
+pub fn conjugate<A, Si, So>(a: &ArrayBase<Si, Ix2>) -> ArrayBase<So, Ix2>
+    where A: Conjugate,
+          Si: Data<Elem = A>,
+          So: DataOwned<Elem = A> + DataMut
+{
+    let mut a = replicate(&a.t());
+    for val in a.iter_mut() {
+        *val = Conjugate::conj(*val);
+    }
+    a
+}
+
+/// Random matrix
+pub fn random<A, S>(n: usize, m: usize) -> ArrayBase<S, Ix2>
+    where A: RandNormal,
+          S: DataOwned<Elem = A>
+{
+    let mut rng = thread_rng();
+    let v: Vec<A> = (0..n * m).map(|_| A::randn(&mut rng)).collect();
+    ArrayBase::from_shape_vec((n, m), v).unwrap()
+}
+
+/// Random square matrix
+pub fn random_square<A, S>(n: usize) -> ArrayBase<S, Ix2>
+    where A: RandNormal,
+          S: DataOwned<Elem = A>
+{
+    random(n, n)
+}
+
+/// Random Hermite matrix
+pub fn random_hermite<A, S>(n: usize) -> ArrayBase<S, Ix2>
+    where A: RandNormal + Conjugate + Add<Output = A>,
+          S: DataOwned<Elem = A> + DataMut
+{
+    let mut a = random_square(n);
+    for i in 0..n {
+        a[(i, i)] = a[(i, i)] + Conjugate::conj(a[(i, i)]);
+        for j in (i + 1)..n {
+            a[(i, j)] = Conjugate::conj(a[(j, i)])
+        }
+    }
+    a
+}
+
+/// Random Hermite Positive-definite matrix
+pub fn random_hpd<A, S>(n: usize) -> ArrayBase<S, Ix2>
+    where A: RandNormal + Conjugate + LinalgScalar,
+          S: DataOwned<Elem = A> + DataMut
+{
+    let a: Array2<A> = random_square(n);
+    let ah: Array2<A> = conjugate(&a);
+    replicate(&ah.dot(&a))
+}
+
+/// construct matrix from diag
+pub fn from_diag<A>(d: &[A]) -> Array2<A>
+    where A: LinalgScalar
+{
+    let n = d.len();
+    let mut e = Array::zeros((n, n));
+    for i in 0..n {
+        e[(i, i)] = d[i];
+    }
+    e
+}
+
+/// stack vectors into matrix horizontally
+pub fn hstack<A, S>(xs: &[ArrayBase<S, Ix1>]) -> Result<Array<A, Ix2>>
+    where A: LinalgScalar,
+          S: Data<Elem = A>
+{
+    let views: Vec<_> = xs.iter()
+        .map(|x| {
+            let n = x.len();
+            x.view().into_shape((n, 1)).unwrap()
+        })
+        .collect();
+    stack(Axis(1), &views).map_err(|e| e.into())
+}
+
+/// stack vectors into matrix vertically
+pub fn vstack<A, S>(xs: &[ArrayBase<S, Ix1>]) -> Result<Array<A, Ix2>>
+    where A: LinalgScalar,
+          S: Data<Elem = A>
+{
+    let views: Vec<_> = xs.iter()
+        .map(|x| {
+            let n = x.len();
+            x.view().into_shape((1, n)).unwrap()
+        })
+        .collect();
+    stack(Axis(0), &views).map_err(|e| e.into())
+}

src/lib.rs

@@ -35,6 +35,7 @@ extern crate blas;
 extern crate lapack;
 extern crate num_traits;
 extern crate num_complex;
+extern crate rand;
 #[macro_use(s)]
 extern crate ndarray;
 #[macro_use]
@@ -61,5 +62,6 @@ pub mod square;
 pub mod triangular;
 
 pub mod util;
+pub mod generate;
 pub mod assert;
 pub mod prelude;

src/prelude.rs

@@ -3,6 +3,8 @@ pub use matrix::Matrix;
 pub use square::SquareMatrix;
 pub use triangular::*;
 pub use util::*;
+pub use types::*;
+pub use generate::*;
 pub use assert::*;
 
 pub use qr::*;

src/triangular.rs

@@ -7,7 +7,7 @@ use super::impl2::UPLO;
 use super::matrix::{Matrix, MFloat};
 use super::square::SquareMatrix;
 use super::error::LinalgError;
-use super::util::hstack;
+use super::generate::hstack;
 use super::impls::solve::ImplSolve;
 
 pub trait SolveTriangular<Rhs>: Matrix + SquareMatrix {

src/types.rs

@@ -1,8 +1,11 @@
 
 pub use num_complex::Complex32 as c32;
 pub use num_complex::Complex64 as c64;
+use num_complex::Complex;
 use num_traits::Float;
 use std::ops::*;
+use rand::Rng;
+use rand::distributions::*;
 
 pub trait AssociatedReal: Sized {
     type Real: Float + Mul<Self, Output = Self>;
@@ -11,21 +14,91 @@ pub trait AssociatedComplex: Sized {
     type Complex;
 }
 
-macro_rules! impl_assoc {
+/// Field with norm
+pub trait Absolute {
+    type Output: Float;
+    fn squared(&self) -> Self::Output;
+    fn abs(&self) -> Self::Output {
+        self.squared().sqrt()
+    }
+}
+
+pub trait Conjugate: Copy {
+    fn conj(self) -> Self;
+}
+
+pub trait RandNormal {
+    fn randn<R: Rng>(&mut R) -> Self;
+}
+
+macro_rules! impl_traits {
     ($real:ty, $complex:ty) => {
+
 impl AssociatedReal for $real {
     type Real = $real;
 }
+
 impl AssociatedReal for $complex {
     type Real = $real;
 }
+
 impl AssociatedComplex for $real {
     type Complex = $complex;
 }
+
 impl AssociatedComplex for $complex {
     type Complex = $complex;
 }
-}} // impl_assoc!
 
-impl_assoc!(f64, c64);
-impl_assoc!(f32, c32);
+impl Absolute for $real {
+    type Output = Self;
+    fn squared(&self) -> Self::Output {
+        *self * *self
+    }
+    fn abs(&self) -> Self::Output {
+        Float::abs(*self)
+    }
+}
+
+impl Absolute for $complex {
+    type Output = $real;
+    fn squared(&self) -> Self::Output {
+        self.norm_sqr()
+    }
+    fn abs(&self) -> Self::Output {
+        self.norm()
+    }
+}
+
+impl Conjugate for $real {
+    fn conj(self) -> Self {
+        self
+    }
+}
+
+impl Conjugate for $complex {
+    fn conj(self) -> Self {
+        Complex::conj(&self)
+    }
+}
+
+impl RandNormal for $real {
+    fn randn<R: Rng>(rng: &mut R) -> Self {
+        let dist = Normal::new(0., 1.);
+        dist.ind_sample(rng) as $real
+    }
+}
+
+impl RandNormal for $complex {
+    fn randn<R: Rng>(rng: &mut R) -> Self {
+        let dist = Normal::new(0., 1.);
+        let re = dist.ind_sample(rng) as $real;
+        let im = dist.ind_sample(rng) as $real;
+        Self::new(re, im)
+    }
+}
+
+}} // impl_traits!
+
+impl_traits!(f64, c64);
+impl_traits!(f32, c32);

src/util.rs

@@ -3,48 +3,10 @@
 use std::iter::Sum;
 use ndarray::*;
 use num_traits::Float;
-use super::vector::*;
 use std::ops::Div;
 
-/// construct matrix from diag
-pub fn from_diag<A>(d: &[A]) -> Array2<A>
-    where A: LinalgScalar
-{
-    let n = d.len();
-    let mut e = Array::zeros((n, n));
-    for i in 0..n {
-        e[(i, i)] = d[i];
-    }
-    e
-}
-
-/// stack vectors into matrix horizontally
-pub fn hstack<A, S>(xs: &[ArrayBase<S, Ix1>]) -> Result<Array<A, Ix2>, ShapeError>
-    where A: LinalgScalar,
-          S: Data<Elem = A>
-{
-    let views: Vec<_> = xs.iter()
-        .map(|x| {
-            let n = x.len();
-            x.view().into_shape((n, 1)).unwrap()
-        })
-        .collect();
-    stack(Axis(1), &views)
-}
-
-/// stack vectors into matrix vertically
-pub fn vstack<A, S>(xs: &[ArrayBase<S, Ix1>]) -> Result<Array<A, Ix2>, ShapeError>
-    where A: LinalgScalar,
-          S: Data<Elem = A>
-{
-    let views: Vec<_> = xs.iter()
-        .map(|x| {
-            let n = x.len();
-            x.view().into_shape((1, n)).unwrap()
-        })
-        .collect();
-    stack(Axis(0), &views)
-}
+use super::types::*;
+use super::vector::*;
 
 pub enum NormalizeAxis {
     Row = 0,

src/vector.rs

@@ -43,40 +43,3 @@ impl<A, S, D, T> Norm for ArrayBase<S, D>
         })
     }
 }
-
-/// Field with norm
-pub trait Absolute {
-    type Output: Float;
-    fn squared(&self) -> Self::Output;
-    fn abs(&self) -> Self::Output {
-        self.squared().sqrt()
-    }
-}
-
-macro_rules! impl_abs {
-    ($f:ty, $c:ty) => {
-
-impl Absolute for $f {
-    type Output = Self;
-    fn squared(&self) -> Self::Output {
-        *self * *self
-    }
-    fn abs(&self) -> Self::Output {
-        Float::abs(*self)
-    }
-}
-
-impl Absolute for $c {
-    type Output = $f;
-    fn squared(&self) -> Self::Output {
-        self.norm_sqr()
-    }
-    fn abs(&self) -> Self::Output {
-        self.norm()
-    }
-}
-
-}} // impl_abs!
-
-impl_abs!(f64, c64);
-impl_abs!(f32, c32);

tests/cholesky.rs

@@ -1,24 +1,14 @@
 
-extern crate rand_extra;
 extern crate ndarray;
-extern crate ndarray_rand;
 #[macro_use]
 extern crate ndarray_linalg;
 
-use rand_extra::*;
 use ndarray::*;
-use ndarray_rand::RandomExt;
 use ndarray_linalg::prelude::*;
 
-pub 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())
-}
-
 #[test]
 fn cholesky() {
-    let a = random_hermite(3);
+    let a: Array2<f64> = random_hpd(3);
     println!("a = \n{:?}", a);
     let c: Array2<_> = (&a).cholesky(UPLO::Upper).unwrap();
     println!("c = \n{:?}", c);
@@ -28,7 +18,7 @@ fn cholesky() {
 
 #[test]
 fn cholesky_t() {
-    let a = random_hermite(3);
+    let a: Array2<f64> = random_hpd(3);
     println!("a = \n{:?}", a);
     let c: Array2<_> = (&a).cholesky(UPLO::Upper).unwrap();
     println!("c = \n{:?}", c);