From 9d8a9659604d1f8a5a2fc9af52116b98b52cb188 Mon Sep 17 00:00:00 2001 From: ignacia-fp Date: Fri, 20 Dec 2024 18:47:54 +0000 Subject: [PATCH] null space implemented for RlstScalar --- src/dense/linalg/null_space.rs | 158 ++++++++++++++++----------------- 1 file changed, 76 insertions(+), 82 deletions(-) diff --git a/src/dense/linalg/null_space.rs b/src/dense/linalg/null_space.rs index 7d5d626..94fb546 100644 --- a/src/dense/linalg/null_space.rs +++ b/src/dense/linalg/null_space.rs @@ -1,8 +1,10 @@ use crate::dense::array::Array; use crate::dense::traits::{RawAccessMut, Shape, Stride, UnsafeRandomAccessByValue, DefaultIterator}; -use crate::dense::types::{c32, c64, RlstResult, RlstScalar}; -use crate::{BaseArray, VectorContainer, rlst_dynamic_array1, rlst_dynamic_array2, empty_array}; +use crate::dense::types::{RlstResult, RlstScalar}; +use crate::{empty_array, rlst_dynamic_array1, rlst_dynamic_array2, BaseArray, DynamicArray, MatrixSvd, SvdMode, VectorContainer}; use itertools::min; +use num::{One, Zero}; + /// Compute the matrix nullspace. /// @@ -35,28 +37,20 @@ pub trait MatrixNull: RlstScalar { ) -> RlstResult>; } -macro_rules! implement_into_null { - ($scalar:ty) => { - impl MatrixNull for $scalar { - fn into_null_alloc< - ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self> - + Stride<2> - + Shape<2> - + RawAccessMut, - >( - arr: Array, - tol: ::Real - ) -> RlstResult> { - NullSpace::<$scalar>::new(arr, tol) - } - } - }; -} -implement_into_null!(f32); -implement_into_null!(f64); -implement_into_null!(c32); -implement_into_null!(c64); +impl MatrixNull for T { + fn into_null_alloc< + ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self> + + Stride<2> + + Shape<2> + + RawAccessMut, + >( + arr: Array, + tol: ::Real + ) -> RlstResult> { + NullSpace::::new(arr, tol) + } +} impl< Item: RlstScalar + MatrixNull, @@ -74,7 +68,7 @@ impl< type RealScalar = ::Real; -/// QR decomposition +/// Null Space pub struct NullSpace< Item: RlstScalar > { @@ -82,65 +76,65 @@ pub struct NullSpace< pub null_space_arr: Array, 2>, 2>, } -macro_rules! implement_null_space { - ($scalar:ty) => { - impl NullSpace<$scalar> - { - - fn new< - ArrayImpl: UnsafeRandomAccessByValue<2, Item = $scalar> - + Stride<2> - + Shape<2> - + RawAccessMut, - >(arr: Array<$scalar, ArrayImpl, 2>, tol: RealScalar<$scalar>) -> RlstResult { - - let shape = arr.shape(); - let dim: usize = min(shape).unwrap(); - let mut singular_values = rlst_dynamic_array1!(RealScalar<$scalar>, [dim]); - let mode = crate::dense::linalg::svd::SvdMode::Full; - let mut u = rlst_dynamic_array2!($scalar, [shape[0], shape[0]]); - let mut vt = rlst_dynamic_array2!($scalar, [shape[1], shape[1]]); - - arr.into_svd_alloc(u.view_mut(), vt.view_mut(), singular_values.data_mut(), mode).unwrap(); - - //For a full rank rectangular matrix, then rank = dim. - //find_matrix_rank checks if the matrix is full rank and recomputes the rank. - let rank: usize = Self::find_matrix_rank(singular_values, dim, tol); - - //The null space is given by the last shape[1]-rank columns of V - let mut null_space_arr = empty_array(); - null_space_arr.fill_from_resize(vt.conj().transpose().into_subview([0, rank], [shape[1], shape[1]-rank])); - - Ok(Self {null_space_arr}) - - } - - fn find_matrix_rank(singular_values: Array, BaseArray, VectorContainer>, 1>, 1>, dim: usize, tol:<$scalar as RlstScalar>::Real)->usize{ - //We compute the rank of the matrix by expecting the values of the elements in the diagonal of R. - let mut singular_values_copy: Array, BaseArray, VectorContainer>, 1>, 1> = rlst_dynamic_array1!(RealScalar<$scalar>, [dim]); - singular_values_copy.fill_from(singular_values.view()); - - let max: RealScalar<$scalar> = singular_values.view().iter().max_by(|a, b| a.abs().total_cmp(&b.abs())).unwrap().abs().into(); - let rank: usize; - - if max.re() > 0.0{ - let alpha: RealScalar<$scalar> = (1.0/max); - singular_values_copy.scale_inplace(alpha); - let aux_vec = singular_values_copy.iter().filter(|el| el.abs() > tol ).collect::>(); - rank = aux_vec.len(); - } - else{ - rank = dim; - } - rank - - } +///NullSpaceComputation creates the null space decomposition and saves it in the NullSpace Struct +pub trait NullSpaceComputation{ + ///This trait is implemented for RlstScalar (ie. f32, f64, c32, c64) + type Item: RlstScalar; + + ///We create the null space decomposition + fn new + + Stride<2> + RawAccessMut + Shape<2>> + (arr: Array, tol: RealScalar) -> RlstResult where Self: Sized; + + ///This function helps us to find the matrix rank + fn find_matrix_rank(singular_values: &mut DynamicArray, 1>, dim: usize, tol:RealScalar)->usize; +} + + +impl NullSpaceComputation for NullSpace { + type Item = T; + + fn new + + Stride<2> + + RawAccessMut + + Shape<2>>(arr: Array, tol: RealScalar) -> RlstResult { + + let shape: [usize; 2] = arr.shape(); + let dim: usize = min(shape).unwrap(); + let mut singular_values: DynamicArray, 1> = rlst_dynamic_array1!(RealScalar, [dim]); + let mode: SvdMode = SvdMode::Full; + let mut u: DynamicArray = rlst_dynamic_array2!(Self::Item, [shape[0], shape[0]]); + let mut vt: DynamicArray = rlst_dynamic_array2!(Self::Item, [shape[1], shape[1]]); + + arr.into_svd_alloc(u.view_mut(), vt.view_mut(), singular_values.data_mut(), mode).unwrap(); + + //For a full rank rectangular matrix, then rank = dim. + //find_matrix_rank checks if the matrix is full rank and recomputes the rank. + let rank: usize = Self::find_matrix_rank(&mut singular_values, dim, tol); + + //The null space is given by the last shape[1]-rank columns of V + let mut null_space_arr: DynamicArray = empty_array(); + null_space_arr.fill_from_resize(vt.conj().transpose().into_subview([0, rank], [shape[1], shape[1]-rank])); + + Ok(Self{null_space_arr}) + + } + + fn find_matrix_rank(singular_values: &mut DynamicArray, 1>, dim: usize, tol:::Real)->usize{ + //We compute the rank of the matrix by expecting the values of the elements in the diagonal of R. + let max: RealScalar = singular_values.view().iter().max_by(|a, b| (a.abs().partial_cmp(&b.abs())).unwrap()).unwrap().abs(); + let mut rank: usize = dim; + + if max.re() > as Zero>::zero(){ + let alpha: RealScalar = as One>::one()/max; + singular_values.scale_inplace(alpha); + let aux_vec: Vec> = singular_values.iter().filter(|el| el.abs() > tol ).collect::>>(); + rank = aux_vec.len(); } - }; + + rank + + } } -implement_null_space!(f64); -implement_null_space!(f32); -implement_null_space!(c64); -implement_null_space!(c32);