rmath/src/hludecomposition.rs

use num_traits::One;
use num_traits::Zero;
use std::ops::Add;
use std::ops::Div;
use std::ops::DivAssign;
use std::ops::Mul;
use std::ops::MulAssign;
use std::ops::Neg;
use std::ops::SubAssign;

use crate::err::RmathError;

use super::hpivot::HPivot;
use super::matrix::LUDecompositionImpl;
use super::matrix::NumericalMatrixImpl;
use super::HMatrix;
use super::LUDecomposition;
use super::Matrix;
use super::NumericalMatrix;
use super::Pivot;
use super::SingularMatrixError;

use super::misc::EnhancedOption;

pub struct HLUDecomposition<TYPE> {
    matrix: HMatrix<TYPE>,
    pivot: HPivot,
}

impl<TYPE> HLUDecomposition<TYPE> 
where TYPE : Zero + SubAssign<TYPE> + Clone,
    for<'a> &'a TYPE : Mul<&'a TYPE, Output = TYPE>,
    for<'a> TYPE : DivAssign<&'a TYPE> {
    pub (crate) fn new(matrix: HMatrix<TYPE>, pivot: Option<HPivot>) -> Self {
        let size = matrix.rows();
        HLUDecomposition {
            matrix,
            pivot: pivot.otherwise(|| HPivot::new(size)),
        }
    }

    pub fn solve(&self, b : &HMatrix<TYPE>) -> Result<HMatrix<TYPE>, Box<dyn RmathError>> {
        let mut x = HMatrix::<TYPE>::new(b.rows(), b.columns(), |_| TYPE::zero());
        for n in 0..b.columns() {
            for i in 0..self.matrix.rows() {
                x[(i, n)] = b[(self.pivot[i], n)].clone();
                for k in 0..i {
                    let sub = &self.matrix[(i, k)] * &x[(k, n)];
                    x[(i, n)] -= sub;
                }
            }
            for i in (0..self.matrix.rows()).rev() {
                for k in (i + 1)..self.matrix.rows() {
                    let sub = &self.matrix[(i, k)] * &x[(k, n)];
                    x[(i, n)] -= sub;
                }
                if self.matrix[(i, i)].is_zero() {
                    return Err(Box::<SingularMatrixError>::default());
                } else {
                    x[(i, n)] /= &self.matrix[(i, i)];
                }
            }
        }
        Ok(x)
    }
}

impl<TYPE> LUDecompositionImpl<TYPE, HMatrix<TYPE>, HPivot> for HLUDecomposition<TYPE>
where
    for<'a> TYPE: 'a + One + MulAssign<&'a TYPE> + Neg<Output = TYPE>,
{
    fn get_matrix(&self) -> &HMatrix<TYPE> {
        &self.matrix
    }

    fn get_pivot(&self) -> &HPivot {
        &self.pivot
    }
}

impl<TYPE> LUDecomposition<TYPE, HMatrix<TYPE>, HPivot> for HLUDecomposition<TYPE>
where
    for<'b> TYPE: 'b
        + Clone
        + Zero
        + One
        + PartialEq
        + PartialOrd
        + Add<&'b TYPE, Output = TYPE>
        + Neg<Output = TYPE>
        + Div<Output = TYPE>
        + SubAssign<TYPE>
        + SubAssign<&'b TYPE>
        + DivAssign<TYPE>
        + DivAssign<&'b TYPE>
        + MulAssign<&'b TYPE>,
    for<'c> &'c TYPE:
        Mul<Output = TYPE> + Add<Output = TYPE> + Div<Output = TYPE> + Neg<Output = TYPE>,
{
    fn l(&self) -> HMatrix<TYPE> {
        self.matrix.clone().tril_replace(TYPE::one())
    }

    fn u(&self) -> HMatrix<TYPE> {
        self.matrix.clone().triu()
    }

    fn invert(&self) -> HMatrix<TYPE> {
        let mut result =
            HMatrix::<TYPE>::new(self.matrix.rows(), self.matrix.columns(), |_| TYPE::zero());
        NumericalMatrixImpl::<TYPE, HPivot, HLUDecomposition<TYPE>>::lu_invert(
            &self.matrix,
            &mut result,
            &mut self.pivot.clone(),
        );
        result
    }

    fn pivot(&self, m: HMatrix<TYPE>) -> HMatrix<TYPE> {
        &self.pivot * m
    }

    // fn p<'a>(&'a self) -> &'a HPivot {
    //     &self.pivot
    // }
}


#[cfg(test)]
mod test {
    use rand::rngs::StdRng;
    use rand::RngCore;
    use rand::SeedableRng;

    use crate::NumericalMatrix;
    use crate::Rational;
    use crate::HMatrix;
    use num_traits::Zero;

    use rand;

    #[test]
    fn solve_linear_system() {
        let mut rand = {
            let seed = [
                1,0,1,3, 
                2,5,0,0, 
                200,1,0,0, 
                210,30,0,0,
                78,134,31,0, 
                253,11,7,0, 
                120,169,89,48, 
                200,0,202,0
            ];
            StdRng::from_seed(seed)
        };
        let mtx = HMatrix::<Rational<i64>>::new(5, 5, |_| {
            Rational::new(i64::try_from(rand.next_u32() % 40).unwrap() - 20, 1)
        });
        let b = HMatrix::<Rational<i64>>::new(5, 5, |_| {
            Rational::new(i64::try_from(rand.next_u32() % 20).unwrap() - 10, 1)
        });
        let lu = mtx.clone().lu().unwrap();
        let x = lu.solve(&b).unwrap();
        assert!((mtx * x - b).is_zero());
    }

}


#[cfg(test)]
mod test_big_int {
    use num_traits::Zero;
    use rand::rngs::StdRng;
    use rand::RngCore;
    use rand::SeedableRng;

    use crate::NumericalMatrix;
    use crate::Rational;
    use crate::HMatrix;

    use rand;

    #[test]
    fn solve_linear_system() {
        use num_bigint::BigInt;
        use num_traits::One;

        let mut rand = {
            let seed = [
                1,0,1,3, 
                2,5,0,0, 
                200,1,0,0, 
                210,30,0,0,
                78,134,31,0, 
                253,11,7,0, 
                120,169,89,48, 
                200,0,202,0
            ];
            StdRng::from_seed(seed)
        };
        let mtx = HMatrix::<Rational<BigInt>>::new(10, 10, |_| {
            Rational::new(BigInt::from(rand.next_u32() % 200) - 100, BigInt::one())
        });
        let b = HMatrix::<Rational<BigInt>>::new(10, 10, |_| {
            Rational::new(BigInt::from(rand.next_u32() % 100) - 50, BigInt::one())
        });
        let lu = mtx.clone().lu().unwrap();
        let x = lu.solve(&b).unwrap();

        assert!((mtx * x - b).is_zero());

    }

}