#[cfg(feature = "serde-serialize")]
use serde::{Deserialize, Serialize};
use num::Signed;
use std::cmp;
use na::allocator::Allocator;
use na::dimension::{Dim, DimMin, DimMinimum, U1};
use na::storage::Storage;
use na::{DefaultAllocator, Matrix, MatrixMN, MatrixN, Scalar, VectorN};
use lapack;
#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
#[cfg_attr(
feature = "serde-serialize",
serde(bound(
serialize = "DefaultAllocator: Allocator<N, DimMinimum<R, C>> +
Allocator<N, R, R> +
Allocator<N, C, C>,
MatrixN<N, R>: Serialize,
MatrixN<N, C>: Serialize,
VectorN<N, DimMinimum<R, C>>: Serialize"
))
)]
#[cfg_attr(
feature = "serde-serialize",
serde(bound(
serialize = "DefaultAllocator: Allocator<N, DimMinimum<R, C>> +
Allocator<N, R, R> +
Allocator<N, C, C>,
MatrixN<N, R>: Deserialize<'de>,
MatrixN<N, C>: Deserialize<'de>,
VectorN<N, DimMinimum<R, C>>: Deserialize<'de>"
))
)]
#[derive(Clone, Debug)]
pub struct SVD<N: Scalar, R: DimMin<C>, C: Dim>
where DefaultAllocator: Allocator<N, R, R> + Allocator<N, DimMinimum<R, C>> + Allocator<N, C, C>
{
pub u: MatrixN<N, R>,
pub vt: MatrixN<N, C>,
pub singular_values: VectorN<N, DimMinimum<R, C>>,
}
impl<N: Scalar, R: DimMin<C>, C: Dim> Copy for SVD<N, R, C>
where
DefaultAllocator: Allocator<N, C, C> + Allocator<N, R, R> + Allocator<N, DimMinimum<R, C>>,
MatrixMN<N, R, R>: Copy,
MatrixMN<N, C, C>: Copy,
VectorN<N, DimMinimum<R, C>>: Copy,
{}
pub trait SVDScalar<R: DimMin<C>, C: Dim>: Scalar
where DefaultAllocator: Allocator<Self, R, R>
+ Allocator<Self, R, C>
+ Allocator<Self, DimMinimum<R, C>>
+ Allocator<Self, C, C>
{
fn compute(m: MatrixMN<Self, R, C>) -> Option<SVD<Self, R, C>>;
}
impl<N: SVDScalar<R, C>, R: DimMin<C>, C: Dim> SVD<N, R, C>
where DefaultAllocator: Allocator<N, R, R>
+ Allocator<N, R, C>
+ Allocator<N, DimMinimum<R, C>>
+ Allocator<N, C, C>
{
pub fn new(m: MatrixMN<N, R, C>) -> Option<Self> {
N::compute(m)
}
}
macro_rules! svd_impl(
($t: ty, $lapack_func: path) => (
impl<R: Dim, C: Dim> SVDScalar<R, C> for $t
where R: DimMin<C>,
DefaultAllocator: Allocator<$t, R, C> +
Allocator<$t, R, R> +
Allocator<$t, C, C> +
Allocator<$t, DimMinimum<R, C>> {
fn compute(mut m: MatrixMN<$t, R, C>) -> Option<SVD<$t, R, C>> {
let (nrows, ncols) = m.data.shape();
if nrows.value() == 0 || ncols.value() == 0 {
return None;
}
let job = b'A';
let lda = nrows.value() as i32;
let mut u = unsafe { Matrix::new_uninitialized_generic(nrows, nrows) };
let mut s = unsafe { Matrix::new_uninitialized_generic(nrows.min(ncols), U1) };
let mut vt = unsafe { Matrix::new_uninitialized_generic(ncols, ncols) };
let ldu = nrows.value();
let ldvt = ncols.value();
let mut work = [ 0.0 ];
let mut lwork = -1 as i32;
let mut info = 0;
let mut iwork = unsafe { ::uninitialized_vec(8 * cmp::min(nrows.value(), ncols.value())) };
unsafe {
$lapack_func(job, nrows.value() as i32, ncols.value() as i32, m.as_mut_slice(),
lda, &mut s.as_mut_slice(), u.as_mut_slice(), ldu as i32, vt.as_mut_slice(),
ldvt as i32, &mut work, lwork, &mut iwork, &mut info);
}
lapack_check!(info);
lwork = work[0] as i32;
let mut work = unsafe { ::uninitialized_vec(lwork as usize) };
unsafe {
$lapack_func(job, nrows.value() as i32, ncols.value() as i32, m.as_mut_slice(),
lda, &mut s.as_mut_slice(), u.as_mut_slice(), ldu as i32, vt.as_mut_slice(),
ldvt as i32, &mut work, lwork, &mut iwork, &mut info);
}
lapack_check!(info);
Some(SVD { u: u, singular_values: s, vt: vt })
}
}
impl<R: DimMin<C>, C: Dim> SVD<$t, R, C>
where DefaultAllocator: Allocator<$t, R, C> +
Allocator<$t, C, R> +
Allocator<$t, U1, R> +
Allocator<$t, U1, C> +
Allocator<$t, R, R> +
Allocator<$t, DimMinimum<R, C>> +
Allocator<$t, DimMinimum<R, C>, R> +
Allocator<$t, DimMinimum<R, C>, C> +
Allocator<$t, R, DimMinimum<R, C>> +
Allocator<$t, C, C> {
#[inline]
pub fn recompose(self) -> MatrixMN<$t, R, C> {
let nrows = self.u.data.shape().0;
let ncols = self.vt.data.shape().1;
let min_nrows_ncols = nrows.min(ncols);
let mut res: MatrixMN<_, R, C> = Matrix::zeros_generic(nrows, ncols);
{
let mut sres = res.generic_slice_mut((0, 0), (min_nrows_ncols, ncols));
sres.copy_from(&self.vt.rows_generic(0, min_nrows_ncols));
for i in 0 .. min_nrows_ncols.value() {
let eigval = self.singular_values[i];
let mut row = sres.row_mut(i);
row *= eigval;
}
}
self.u * res
}
#[inline]
pub fn pseudo_inverse(&self, epsilon: $t) -> MatrixMN<$t, C, R> {
let nrows = self.u.data.shape().0;
let ncols = self.vt.data.shape().1;
let min_nrows_ncols = nrows.min(ncols);
let mut res: MatrixMN<_, C, R> = Matrix::zeros_generic(ncols, nrows);
{
let mut sres = res.generic_slice_mut((0, 0), (min_nrows_ncols, nrows));
self.u.columns_generic(0, min_nrows_ncols).transpose_to(&mut sres);
for i in 0 .. min_nrows_ncols.value() {
let eigval = self.singular_values[i];
let mut row = sres.row_mut(i);
if eigval.abs() > epsilon {
row /= eigval
}
else {
row.fill(0.0);
}
}
}
self.vt.tr_mul(&res)
}
#[inline]
pub fn rank(&self, epsilon: $t) -> usize {
let mut i = 0;
for e in self.singular_values.as_slice().iter() {
if e.abs() > epsilon {
i += 1;
}
}
i
}
}
);
);
svd_impl!(f32, lapack::sgesdd);
svd_impl!(f64, lapack::dgesdd);