From 2fba371aa31f9e333e1020eb9e102598cb816fe0 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Wed, 14 Jun 2017 16:13:00 +0900 Subject: [PATCH 1/6] Merge Absolute traint into types.rs --- src/types.rs | 42 ++++++++++++++++++++++++++++++++++++++---- src/vector.rs | 37 ------------------------------------- 2 files changed, 38 insertions(+), 41 deletions(-) diff --git a/src/types.rs b/src/types.rs index 2d0393d8..7cda5b65 100644 --- a/src/types.rs +++ b/src/types.rs @@ -11,21 +11,55 @@ 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() + } +} + +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_traits! + +impl_traits!(f64, c64); +impl_traits!(f32, c32); diff --git a/src/vector.rs b/src/vector.rs index 46b88f88..e770e7ed 100644 --- a/src/vector.rs +++ b/src/vector.rs @@ -43,40 +43,3 @@ impl Norm for ArrayBase }) } } - -/// 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); From 24555e1f17173550293ae78b4e0c3f5c9931bfc3 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Wed, 14 Jun 2017 16:13:18 +0900 Subject: [PATCH 2/6] Add Conjugate trait --- src/types.rs | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/src/types.rs b/src/types.rs index 7cda5b65..3bac3d0b 100644 --- a/src/types.rs +++ b/src/types.rs @@ -1,6 +1,7 @@ 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::*; @@ -20,6 +21,10 @@ pub trait Absolute { } } +pub trait Conjugate: Copy { + fn conj(self) -> Self; +} + macro_rules! impl_traits { ($real:ty, $complex:ty) => { @@ -59,6 +64,18 @@ impl Absolute for $complex { } } +impl Conjugate for $real { + fn conj(self) -> Self { + self + } +} + +impl Conjugate for $complex { + fn conj(self) -> Self { + Complex::conj(&self) + } +} + }} // impl_traits! impl_traits!(f64, c64); From ae27f6fdd7b2e36858eae4e99690d6bbe0e59bb1 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Wed, 14 Jun 2017 16:13:43 +0900 Subject: [PATCH 3/6] Create generate mod --- Cargo.toml | 1 + src/assert.rs | 1 + src/generate.rs | 70 +++++++++++++++++++++++++++++++++++++++++++++++ src/lib.rs | 2 ++ src/triangular.rs | 2 +- src/util.rs | 42 ++-------------------------- 6 files changed, 77 insertions(+), 41 deletions(-) create mode 100644 src/generate.rs diff --git a/Cargo.toml b/Cargo.toml index bc60de63..69535a92 100644 --- a/Cargo.toml +++ b/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" diff --git a/src/assert.rs b/src/assert.rs index 3a1f3fd8..880dba1b 100644 --- a/src/assert.rs +++ b/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(test: A, truth: A, rtol: Tol) -> Result diff --git a/src/generate.rs b/src/generate.rs new file mode 100644 index 00000000..dd99c863 --- /dev/null +++ b/src/generate.rs @@ -0,0 +1,70 @@ + +use ndarray::*; +use rand::*; +use std::ops::*; + +use super::types::*; +use super::error::*; + +pub fn random_square(n: usize) -> ArrayBase + where A: Rand, + S: DataOwned +{ + let mut rng = thread_rng(); + let v: Vec = (0..n * n).map(|_| rng.gen()).collect(); + ArrayBase::from_shape_vec((n, n), v).unwrap() +} + +pub fn random_hermite(n: usize) -> ArrayBase + where A: Rand + Conjugate + Add, + S: DataOwned + 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..n { + a[(i, j)] = a[(j, i)]; + } + } + a +} + +/// construct matrix from diag +pub fn from_diag(d: &[A]) -> Array2 + 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(xs: &[ArrayBase]) -> Result> + where A: LinalgScalar, + S: Data +{ + 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(xs: &[ArrayBase]) -> Result> + where A: LinalgScalar, + S: Data +{ + 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()) +} diff --git a/src/lib.rs b/src/lib.rs index df9e17be..b45c369b 100644 --- a/src/lib.rs +++ b/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; diff --git a/src/triangular.rs b/src/triangular.rs index 99bf7220..cae6b64a 100644 --- a/src/triangular.rs +++ b/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: Matrix + SquareMatrix { diff --git a/src/util.rs b/src/util.rs index aee89d65..03cb9482 100644 --- a/src/util.rs +++ b/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(d: &[A]) -> Array2 - 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(xs: &[ArrayBase]) -> Result, ShapeError> - where A: LinalgScalar, - S: Data -{ - 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(xs: &[ArrayBase]) -> Result, ShapeError> - where A: LinalgScalar, - S: Data -{ - 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, From 7917fd365bbdab1be7cdb5725ec5b600c3f4a83a Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Wed, 14 Jun 2017 17:00:36 +0900 Subject: [PATCH 4/6] Add random_hpd --- src/generate.rs | 29 +++++++++++++++++++++++++++-- src/prelude.rs | 2 ++ tests/cholesky.rs | 14 ++------------ 3 files changed, 31 insertions(+), 14 deletions(-) diff --git a/src/generate.rs b/src/generate.rs index dd99c863..8eff951f 100644 --- a/src/generate.rs +++ b/src/generate.rs @@ -3,9 +3,23 @@ use ndarray::*; use rand::*; use std::ops::*; +use super::layout::*; use super::types::*; use super::error::*; +pub fn conjugate(a: &ArrayBase) -> ArrayBase + where A: Conjugate, + Si: Data, + So: DataOwned + DataMut +{ + let mut a = replicate(&a.t()); + for val in a.iter_mut() { + *val = Conjugate::conj(*val); + } + a +} + +/// Random square matrix pub fn random_square(n: usize) -> ArrayBase where A: Rand, S: DataOwned @@ -15,6 +29,7 @@ pub fn random_square(n: usize) -> ArrayBase ArrayBase::from_shape_vec((n, n), v).unwrap() } +/// Random Hermite matrix pub fn random_hermite(n: usize) -> ArrayBase where A: Rand + Conjugate + Add, S: DataOwned + DataMut @@ -22,13 +37,23 @@ pub fn random_hermite(n: usize) -> ArrayBase 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..n { - a[(i, j)] = a[(j, 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(n: usize) -> ArrayBase + where A: Rand + Conjugate + LinalgScalar, + S: DataOwned + DataMut +{ + let a: Array2 = random_square(n); + let ah: Array2 = conjugate(&a); + replicate(&ah.dot(&a)) +} + /// construct matrix from diag pub fn from_diag(d: &[A]) -> Array2 where A: LinalgScalar diff --git a/src/prelude.rs b/src/prelude.rs index e9899e16..d3f0da73 100644 --- a/src/prelude.rs +++ b/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::*; diff --git a/tests/cholesky.rs b/tests/cholesky.rs index b74fec0b..526d4959 100644 --- a/tests/cholesky.rs +++ b/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 { - let r_dist = RealNormal::new(0., 1.); - let a = Array::::random((n, n), r_dist); - a.dot(&a.t()) -} - #[test] fn cholesky() { - let a = random_hermite(3); + let a: Array2 = 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 = random_hpd(3); println!("a = \n{:?}", a); let c: Array2<_> = (&a).cholesky(UPLO::Upper).unwrap(); println!("c = \n{:?}", c); From 0c28823bcc49a1c361bca291820425befbf7ef37 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Wed, 14 Jun 2017 17:22:14 +0900 Subject: [PATCH 5/6] Add randn(&mut rng) --- src/generate.rs | 10 +++++----- src/types.rs | 22 ++++++++++++++++++++++ 2 files changed, 27 insertions(+), 5 deletions(-) diff --git a/src/generate.rs b/src/generate.rs index 8eff951f..230d4d5d 100644 --- a/src/generate.rs +++ b/src/generate.rs @@ -1,7 +1,7 @@ use ndarray::*; -use rand::*; use std::ops::*; +use rand::*; use super::layout::*; use super::types::*; @@ -21,17 +21,17 @@ pub fn conjugate(a: &ArrayBase) -> ArrayBase /// Random square matrix pub fn random_square(n: usize) -> ArrayBase - where A: Rand, + where A: RandNormal, S: DataOwned { let mut rng = thread_rng(); - let v: Vec = (0..n * n).map(|_| rng.gen()).collect(); + let v: Vec = (0..n * n).map(|_| A::randn(&mut rng)).collect(); ArrayBase::from_shape_vec((n, n), v).unwrap() } /// Random Hermite matrix pub fn random_hermite(n: usize) -> ArrayBase - where A: Rand + Conjugate + Add, + where A: RandNormal + Conjugate + Add, S: DataOwned + DataMut { let mut a = random_square(n); @@ -46,7 +46,7 @@ pub fn random_hermite(n: usize) -> ArrayBase /// Random Hermite Positive-definite matrix pub fn random_hpd(n: usize) -> ArrayBase - where A: Rand + Conjugate + LinalgScalar, + where A: RandNormal + Conjugate + LinalgScalar, S: DataOwned + DataMut { let a: Array2 = random_square(n); diff --git a/src/types.rs b/src/types.rs index 3bac3d0b..d8118b69 100644 --- a/src/types.rs +++ b/src/types.rs @@ -4,6 +4,8 @@ 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; @@ -25,6 +27,10 @@ pub trait Conjugate: Copy { fn conj(self) -> Self; } +pub trait RandNormal { + fn randn(&mut R) -> Self; +} + macro_rules! impl_traits { ($real:ty, $complex:ty) => { @@ -76,6 +82,22 @@ impl Conjugate for $complex { } } +impl RandNormal for $real { + fn randn(rng: &mut R) -> Self { + let dist = Normal::new(0., 1.); + dist.ind_sample(rng) as $real + } +} + +impl RandNormal for $complex { + fn randn(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); From a56fbc4260ac83542572b657f86e01b11c622e8a Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Wed, 14 Jun 2017 17:26:18 +0900 Subject: [PATCH 6/6] Add generate::random(n, m) --- src/generate.rs | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/src/generate.rs b/src/generate.rs index 230d4d5d..54211c9b 100644 --- a/src/generate.rs +++ b/src/generate.rs @@ -19,14 +19,22 @@ pub fn conjugate(a: &ArrayBase) -> ArrayBase a } +/// Random matrix +pub fn random(n: usize, m: usize) -> ArrayBase + where A: RandNormal, + S: DataOwned +{ + let mut rng = thread_rng(); + let v: Vec = (0..n * m).map(|_| A::randn(&mut rng)).collect(); + ArrayBase::from_shape_vec((n, m), v).unwrap() +} + /// Random square matrix pub fn random_square(n: usize) -> ArrayBase where A: RandNormal, S: DataOwned { - let mut rng = thread_rng(); - let v: Vec = (0..n * n).map(|_| A::randn(&mut rng)).collect(); - ArrayBase::from_shape_vec((n, n), v).unwrap() + random(n, n) } /// Random Hermite matrix