ndarray_linalg/
generate.rs

1//! Generator functions for matrices
2
3use ndarray::*;
4use rand::prelude::*;
5
6use super::convert::*;
7use super::error::*;
8use super::qr::*;
9use super::types::*;
10
11/// Hermite conjugate matrix
12pub fn conjugate<A, Si, So>(a: &ArrayBase<Si, Ix2>) -> ArrayBase<So, Ix2>
13where
14    A: Scalar,
15    Si: Data<Elem = A>,
16    So: DataOwned<Elem = A> + DataMut,
17{
18    let mut a: ArrayBase<So, Ix2> = replicate(&a.t());
19    for val in a.iter_mut() {
20        *val = val.conj();
21    }
22    a
23}
24
25/// Generate random array
26pub fn random<A, S, Sh, D>(sh: Sh) -> ArrayBase<S, D>
27where
28    A: Scalar,
29    S: DataOwned<Elem = A>,
30    D: Dimension,
31    Sh: ShapeBuilder<Dim = D>,
32{
33    let mut rng = thread_rng();
34    ArrayBase::from_shape_fn(sh, |_| A::rand(&mut rng))
35}
36
37/// Generate random unitary matrix using QR decomposition
38///
39/// Be sure that this it **NOT** a uniform distribution. Use it only for test purpose.
40pub fn random_unitary<A>(n: usize) -> Array2<A>
41where
42    A: Scalar + Lapack,
43{
44    let a: Array2<A> = random((n, n));
45    let (q, _r) = a.qr_into().unwrap();
46    q
47}
48
49/// Generate random regular matrix
50///
51/// Be sure that this it **NOT** a uniform distribution. Use it only for test purpose.
52pub fn random_regular<A>(n: usize) -> Array2<A>
53where
54    A: Scalar + Lapack,
55{
56    let a: Array2<A> = random((n, n));
57    let (q, mut r) = a.qr_into().unwrap();
58    for i in 0..n {
59        r[(i, i)] = A::one() + A::from_real(r[(i, i)].abs());
60    }
61    q.dot(&r)
62}
63
64/// Random Hermite matrix
65pub fn random_hermite<A, S>(n: usize) -> ArrayBase<S, Ix2>
66where
67    A: Scalar,
68    S: DataOwned<Elem = A> + DataMut,
69{
70    let mut a: ArrayBase<S, Ix2> = random((n, n));
71    for i in 0..n {
72        a[(i, i)] = a[(i, i)] + a[(i, i)].conj();
73        for j in (i + 1)..n {
74            a[(i, j)] = a[(j, i)].conj();
75        }
76    }
77    a
78}
79
80/// Random Hermite Positive-definite matrix
81///
82/// - Eigenvalue of matrix must be larger than 1 (thus non-singular)
83///
84pub fn random_hpd<A, S>(n: usize) -> ArrayBase<S, Ix2>
85where
86    A: Scalar,
87    S: DataOwned<Elem = A> + DataMut,
88{
89    let a: Array2<A> = random((n, n));
90    let ah: Array2<A> = conjugate(&a);
91    ArrayBase::eye(n) + &ah.dot(&a)
92}
93
94/// construct matrix from diag
95pub fn from_diag<A>(d: &[A]) -> Array2<A>
96where
97    A: Scalar,
98{
99    let n = d.len();
100    let mut e = Array::zeros((n, n));
101    for i in 0..n {
102        e[(i, i)] = d[i];
103    }
104    e
105}
106
107/// stack vectors into matrix horizontally
108pub fn hstack<A, S>(xs: &[ArrayBase<S, Ix1>]) -> Result<Array<A, Ix2>>
109where
110    A: Scalar,
111    S: Data<Elem = A>,
112{
113    let views: Vec<_> = xs.iter().map(|x| x.view()).collect();
114    stack(Axis(1), &views).map_err(Into::into)
115}
116
117/// stack vectors into matrix vertically
118pub fn vstack<A, S>(xs: &[ArrayBase<S, Ix1>]) -> Result<Array<A, Ix2>>
119where
120    A: Scalar,
121    S: Data<Elem = A>,
122{
123    let views: Vec<_> = xs.iter().map(|x| x.view()).collect();
124    stack(Axis(0), &views).map_err(Into::into)
125}