Skip to content

Commit

Permalink
null space implemented for RlstScalar
Browse files Browse the repository at this point in the history
  • Loading branch information
@ignacia-fp
ignacia-fp committed Dec 20, 2024
1 parent ae4dfe7 commit 9d8a965
Showing 1 changed file with 76 additions and 82 deletions.
158 changes: 76 additions & 82 deletions src/dense/linalg/null_space.rs
View file
Original file line number Diff line number Diff line change
@@ -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.
///
Expand Down Expand Up @@ -35,28 +37,20 @@ pub trait MatrixNull: RlstScalar {
) -> RlstResult<NullSpace<Self>>;
}

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<Item = Self>,
>(
arr: Array<Self, ArrayImpl, 2>,
tol: <Self as RlstScalar>::Real
) -> RlstResult<NullSpace<Self>> {
NullSpace::<$scalar>::new(arr, tol)
}
}
};
}

implement_into_null!(f32);
implement_into_null!(f64);
implement_into_null!(c32);
implement_into_null!(c64);
impl <T: RlstScalar + MatrixSvd> MatrixNull for T {
fn into_null_alloc<
ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self>
+ Stride<2>
+ Shape<2>
+ RawAccessMut<Item = Self>,
>(
arr: Array<Self, ArrayImpl, 2>,
tol: <Self as RlstScalar>::Real
) -> RlstResult<NullSpace<Self>> {
NullSpace::<Self>::new(arr, tol)
}
}

impl<
Item: RlstScalar + MatrixNull,
Expand All @@ -74,73 +68,73 @@ impl<


type RealScalar<T> = <T as RlstScalar>::Real;
/// QR decomposition
/// Null Space
pub struct NullSpace<
Item: RlstScalar
> {
///Computed null space
pub null_space_arr: Array<Item, BaseArray<Item, VectorContainer<Item>, 2>, 2>,
}

macro_rules! implement_null_space {
($scalar:ty) => {
impl NullSpace<$scalar>
{

fn new<
ArrayImpl: UnsafeRandomAccessByValue<2, Item = $scalar>
+ Stride<2>
+ Shape<2>
+ RawAccessMut<Item = $scalar>,
>(arr: Array<$scalar, ArrayImpl, 2>, tol: RealScalar<$scalar>) -> RlstResult<Self> {

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<RealScalar<$scalar>, BaseArray<RealScalar<$scalar>, VectorContainer<RealScalar<$scalar>>, 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<RealScalar<$scalar>, BaseArray<RealScalar<$scalar>, VectorContainer<RealScalar<$scalar>>, 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::<Vec<_>>();
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<ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item>
+ Stride<2> + RawAccessMut<Item = Self::Item> + Shape<2>>
(arr: Array<Self::Item, ArrayImpl, 2>, tol: RealScalar<Self::Item>) -> RlstResult<Self> where Self: Sized;

///This function helps us to find the matrix rank
fn find_matrix_rank(singular_values: &mut DynamicArray<RealScalar<Self::Item>, 1>, dim: usize, tol:RealScalar<Self::Item>)->usize;
}


impl <T: RlstScalar + MatrixSvd> NullSpaceComputation for NullSpace<T> {
type Item = T;

fn new<ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item>
+ Stride<2>
+ RawAccessMut<Item = Self::Item>
+ Shape<2>>(arr: Array<Self::Item, ArrayImpl, 2>, tol: RealScalar<Self::Item>) -> RlstResult<Self> {

let shape: [usize; 2] = arr.shape();
let dim: usize = min(shape).unwrap();
let mut singular_values: DynamicArray<RealScalar<Self::Item>, 1> = rlst_dynamic_array1!(RealScalar<Self::Item>, [dim]);
let mode: SvdMode = SvdMode::Full;
let mut u: DynamicArray<Self::Item, 2> = rlst_dynamic_array2!(Self::Item, [shape[0], shape[0]]);
let mut vt: DynamicArray<Self::Item, 2> = 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<Self::Item, 2> = 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<RealScalar<Self::Item>, 1>, dim: usize, tol:<Self::Item as RlstScalar>::Real)->usize{
//We compute the rank of the matrix by expecting the values of the elements in the diagonal of R.
let max: RealScalar<Self::Item> = 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() > <RealScalar<Self::Item> as Zero>::zero(){
let alpha: RealScalar<Self::Item> = <RealScalar<Self::Item> as One>::one()/max;
singular_values.scale_inplace(alpha);
let aux_vec: Vec<RealScalar<Self::Item>> = singular_values.iter().filter(|el| el.abs() > tol ).collect::<Vec<RealScalar<Self::Item>>>();
rank = aux_vec.len();
}
};

rank

}
}

implement_null_space!(f64);
implement_null_space!(f32);
implement_null_space!(c64);
implement_null_space!(c32);

0 comments on commit 9d8a965

Please sign in to comment.