grassmann/functions/
rank.rs

1use crate::{Number, core::matrix::Matrix};
2
3
4
5pub fn rank<T: Number>(A: &Matrix<T>) -> u32 {
6
7    if A.columns > A.rows {
8
9        let At = A.transpose();
10
11        let lu = At.lu();
12
13        let rank = At.columns - lu.d.len();
14
15        rank as u32
16
17    } else {
18
19        let lu = A.lu();
20
21        let rank = A.columns - lu.d.len();
22
23        rank as u32
24    }
25}
26
27
28
29mod tests {
30
31    use crate::{ core::{matrix::{ Matrix }}, matrix, vector };
32
33    #[test]
34    fn rank() {
35        
36        let mut A1 = matrix![f64,
37            5.024026017784438, 2.858902178366669, 296.2138835869165;
38            7.929129970221636, 5.7210492203315795, 523.7802005055211;
39            8.85257084623291, 8.95057121546899, 704.1069012250204;
40        ];
41
42        let A2 = matrix![f64,
43            1.1908477166794595, 8.793086722414468, 194.77132992778556;
44            3.6478484951000456, 4.858421485429982, 187.58571816777294;
45            9.423462238282756, 8.321761784861303, 406.23378670237366;
46        ];
47
48        let lu = A1.lu();
49
50        let R: Matrix<f64> = &lu.L * &lu.U;
51
52        println!("\n d1 is {:?} \n A1 is {} \n R is {} \n L is {} \n U is {} \n diff is {} \n d is {:?} \n P is {} \n \n", lu.d, A1, R, lu.L, lu.U, &A1 - &R, lu.d, lu.P);
53
54        assert_eq!(A1.rank(), 2, "A1 - rank is 2");
55        
56        let lu = A2.lu();
57
58        let R: Matrix<f64> = &lu.L * &lu.U;
59
60        println!("\n d2 is {:?} \n A2 is {} \n R is {} \n L is {} \n U is {} \n diff is {} \n d is {:?} \n P is {} \n \n", lu.d, A2, R, lu.L, lu.U, &A2 - &R, lu.d, lu.P);
61        
62        assert_eq!(A2.rank(), 2, "A2 - rank is 2");
63        
64        let test = 50;
65
66        for i in 1..test {
67            let max = 1.;
68
69            let A: Matrix<f64> = Matrix::rand_shape(i, max);
70            let At: Matrix<f64> = A.transpose();
71            let AtA: Matrix<f64> = &At * &A;
72            let AAt: Matrix<f64> = &A * &At;
73
74            let rank_A = A.rank();
75            let rank_AAt = AAt.rank();
76            let rank_AtA = AtA.rank();
77
78            let eq1 = rank_A == rank_AAt;
79            let eq2 = rank_A == rank_AtA;
80
81            if !eq1 {
82                println!("\n rank A {}, rank AAt {} \n", rank_A, rank_AAt);
83                println!("\n A ({}, {}) is {} \n AAt ({}, {}) {} \n", A.rows, A.columns, A, AAt.rows, AAt.columns, AAt);
84
85                let luA = A.lu();
86                let luAAt = AAt.lu();
87
88                println!("\n U A is {} \n L A is {} \n d A is {:?} \n", luA.U, luA.L, luA.d);
89                println!("\n U AAt is {} \n L AAt is {} \n d AAt is {:?} \n", luAAt.U, luAAt.L, luAAt.d);
90            }
91
92            if !eq2 {
93                println!("\n rank A {}, rank AtA {} \n", rank_A, rank_AtA);
94                println!("\n A ({}, {}) is {} \n AtA ({}, {}) {} \n", A.rows, A.columns, A, AtA.rows, AtA.columns, AtA);
95
96                let luA = A.lu();
97                let luAtA = AtA.lu();
98
99                println!("\n U A is {} \n L A is {} \n d A is {:?} \n", luA.U, luA.L, luA.d);
100                println!("\n U AtA is {} \n L AtA is {} \n d AtA is {:?} \n", luAtA.U, luAtA.L, luAtA.d);
101            }
102            //TODO rank AtA < rank A ??? fix
103            assert!(eq1, "rank A and rank AAt should be equivalent");
104            assert!(eq2, "rank A and rank AtA should be equivalent");
105        }
106    }
107}
grassmann/functions/rank.rs

grassmann/functions/
rank.rs
