1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126

#![doc="Matrix rank computation methods
"]


// std imports

// local imports
use matrix::matrix::{MatrixF64};
use matrix::traits::{Shape, Search};
use matrix::eo::eo_traits::ECO;


/// Computes the rank of a matrix using elementary column operations
pub fn rank_eco(a : & MatrixF64) -> usize {
    let mut a = a.clone();
    let mut rank  = 0;
    let m = a.num_rows();
    let n = a.num_cols();
    // forward elimination row wise
    for k in 0..m{
        println!("a: {}", a);
        let (_, cc) = a.max_abs_scalar_in_row(k, k, n);
        if cc > k {
            // TODO : we can switch only part of column
            a.eco_switch(k, cc);
        }
        let mut v = a.view(k, k, m - k, n - k);
        // Pick the pivot
        let pivot  = v.get(0, 0).unwrap();
        if pivot == 0. {
            // Nothing to be done.
            continue;
        }
        // We have a non-zero pivot
        rank += 1;
        for c in 1..v.num_cols(){
            let first = v.get(0, c).unwrap();
            let factor = first/pivot;
            v.eco_scale_add(c, 0, -factor);
        }
    }
    rank
}


pub fn rank(a : & MatrixF64) -> usize{
    rank_eco(a)
}


/******************************************************
 *
 *   Unit tests follow.
 *
 *******************************************************/

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

    #[test]
    fn test_rank_eco_0(){
        let a = matrix_rw_f64(2, 2, &[
            1., 0.,
            1., 1.]);
        let r = rank(&a);
        assert_eq!(r, 2);
    }

    #[test]
    fn test_rank_eco_1(){
        let a = matrix_rw_f64(4, 4, &[
        4.0, 3.0, 4.0, 4.0,
        4.0, 1.0, 4.0, 2.0,
        1.0, 2.0, 1.0, 4.0,
        4.0, 3.0, 4.0, 1.0
        ]);
        let r = rank(&a);
        assert_eq!(r, 3);
    }

    #[test]
    fn test_rank_eco_2(){
        let a = matrix_rw_f64(3, 5, &[
        2.0, 1.0, 3.0, 2.0, 1.0,
        4.0, 2.0, 3.0, 3.0, 1.0,
        4.0, 2.0, 4.0, 3.0, 2.0
        ]);
        let r = rank(&a);
        assert_eq!(r, 3);
    }

    #[test]
    fn test_rank_eco_3(){
        let a = matrix_rw_f64(4, 5, &[
        2.0, 1.0, 3.0, 2.0, 1.0,
        4.0, 2.0, 3.0, 3.0, 1.0,
        4.0, 2.0, 4.0, 3.0, 2.0,
        10.0, 5.0, 9.0, 8.0, 3.0
        ]);
        let r = rank(&a);
        assert_eq!(r, 3);
    }

    #[test]
    fn test_rank_eco_4(){
        let a = matrix_rw_f64(2, 3, &[
        1.0000000000000, 1.0000000000000, 1.0000000000000,
        2.0000000000000, 2.0000000000000, 2.0000000000002
        ]);
        let r = rank(&a);
        assert_eq!(r, 2);
    }

    #[test]
    fn test_rank_eco_hilbert(){
        for i in 4..50{
            let m = hilbert(i);
            let r = rank(&m);
            assert_eq!(r, i);
        }
    }
}