bayes_estimate 0.20.0

//! Numerical approximation of matrix reciprocal condition numbers 'rcond'.
//!
//! Matrices are well-conditioned if the reciprocal condition number is near 1 and ill-conditioned if it is near zero.
//!
//! For positive (semi-)definite matrices a simple definition is the minimum diagonal element / maximum diagonal element.
//! This is best applied to Cholesky factorisation by using the square of the factors diagonal or for a UdU' by using 'd'
//! directly.
//!
//! A rcond = 0 implies the matrix is semi definite, < 0 implies the matrix is negative.
//!
//! Required for all linear algebra in models and filters.

use nalgebra::{allocator::Allocator, DefaultAllocator, Dim, OMatrix, RealField};

/// Determine the reciprocal condition number of a symmetric matrix.
///
/// The 'rcond' is simply the minimum diagonal element / maximum diagonal element. A NaN result a negative result.
pub fn rcond_symmetric<N: Copy + RealField, R: Dim, C: Dim>(sm: &OMatrix<N, R, C>) -> N
where
    DefaultAllocator: Allocator<R, C>,
{
    // Special case an empty matrix
    let n = sm.nrows();
    if n == 0 {
        N::zero()
    } else {
        let mut mind = sm[(0, 0)];
        let mut maxd = mind;

        for i in 0..n {
            let d = sm[(i, i)];
            if d != d {
                // NaN
                mind = N::one().neg();
                break;
            }
            if d < mind {
                mind = d;
            }
            if d > maxd {
                maxd = d;
            }
        }

        rcond_min_max(mind, maxd)
    }
}

/// Determine the 'rcond' from minimum and maximum diagonal elements.
fn rcond_min_max<N: RealField>(mind: N, maxd: N) -> N {
    if mind < N::zero() {
        // matrix is negative
        mind // mind < 0 but does not represent a rcond
    } else {
        assert!(mind <= maxd); // check sanity

        // NOTE mind -0 will be handled here and propagate into rcond
        let rcond = mind / maxd; // rcond from min/max norm
        if !rcond.is_finite() {
            // NaN, singular due to (mind == maxd) == (zero or infinity)
            N::zero()
        } else {
            assert!(rcond <= N::one());
            rcond
        }
    }
}

#[cfg(test)]
mod test {
    use super::*;

    #[test]
    fn test_simple_rconds() {
        assert_eq!(rcond_min_max(1., 1.), 1.);
        assert_eq!(rcond_min_max(1., 2.), 0.5);
    }

    #[test]
    fn test_special_rconds() {
        assert_eq!(rcond_min_max(-1., 55.), -1.);
        assert_eq!(rcond_min_max(-2.,55.), -2.);
        assert_eq!(rcond_min_max(-0.,55.), -0.);
        assert_eq!(rcond_min_max(-0., 0.), 0.);
        assert_eq!(rcond_min_max(-0., -0.), 0.);
    }

    #[should_panic(expected = "mind <= maxd")]
    #[test]
    fn test_greater_one_panics() {
        rcond_min_max(2., 1.);
    }

}

/// Checks the reciprocal condition number is > 0 .
///
/// IEC 559 NaN values are never true
pub fn check_positive<N: RealField>(rcond: N, message: &str) -> Result<(), &str> {
    if rcond > N::zero() {
        Ok(())
    } else {
        Err(message)
    }
}

/// Checks the reciprocal condition number is >= 0 .
///
/// IEC 559 NaN values are never true
pub fn check_non_negative<N: RealField>(rcond: N, message: &str) -> Result<(), &str> {
    if rcond >= N::zero() {
        Ok(())
    } else {
        Err(message)
    }
}