test: gaussian and gauus jordan elimination

src/lib.rs

@@ -1,4 +1,5 @@
 use fractions::{Fraction, FractionError};
+use std::cmp::min;
 use std::fmt::{Debug, Display};
 use std::ops::Add;
 use std::ops::Mul;
@@ -232,10 +233,10 @@ impl Matrix {
         }
 
         let mut max_row = col;
-        let mut max_value = self.get(col, col).unwrap().abs();
+        let mut max_value = self.get(col, col)?.abs();
 
         for r in (col + 1)..self.rows {
-            let val = self.get(r, col).unwrap().abs();
+            let val = self.get(r, col)?.abs();
             if val > max_value {
                 max_value = val;
                 max_row = r;
@@ -260,24 +261,41 @@ impl Matrix {
             rows: self.rows,
             data: self.data.clone(),
         };
-        let mut sign = Fraction::new(1, 1).unwrap();
+        let mut sign = Fraction::new(1, 1)?;
 
-        for i in 0..self.columns {
+        for i in 0..min(self.rows, self.columns) {
             // We do partial pivoting for better efifiency and security
             let pivot_exists = trig_matrix.partial_pivoting(i, &mut sign)?;
 
             // If there ain't no other thing but 0 then we're
             // fucked, determinant is zero
             if !pivot_exists {
-                return Err(MatrixError::FailedGauss);
+                let mut all_zero = true;
+
+                for r in i..trig_matrix.rows {
+                    if !trig_matrix.get(r, i)?.is_zero() {
+                        all_zero = false;
+                        break;
+                    }
                 }
 
-            let pivot = *trig_matrix.get(i, i).unwrap();
+                if all_zero {
+                    return Err(MatrixError::FailedGauss); // matriz singular
+                }
+
+                continue;
+            }
+
+            let pivot = *trig_matrix.get(i, i)?;
+
+            if pivot.is_zero() {
+                return Err(MatrixError::FailedGauss);
+            }
 
             // The main gaussian elimination, not even I remember how
             // i did it in such a asimple way
             for x in (i + 1)..trig_matrix.rows {
-                let m = (*trig_matrix.get(x, i).unwrap() / pivot).unwrap();
+                let m = (*trig_matrix.get(x, i)? / pivot)?;
 
                 let row_x = trig_matrix.get_row(x)?;
                 let row_i = trig_matrix.get_row(i)?;
@@ -303,14 +321,14 @@ impl Matrix {
 
         let mut dummy = Fraction::from(1);
 
-        for i in 0..self.columns {
+        for i in 0..min(self.columns, self.rows) {
             let pivot_exists = new_matrix.partial_pivoting(i, &mut dummy)?;
 
             if !pivot_exists {
                 return Err(MatrixError::FailedGaussJordan);
             }
 
-            let pivot = *new_matrix.get(i, i).unwrap();
+            let pivot = *new_matrix.get(i, i)?;
 
             let new_pivot_row = new_matrix
                 .get_row(i)?
@@ -324,7 +342,7 @@ impl Matrix {
                     continue;
                 }
 
-                let factor = *new_matrix.get(r, i).unwrap();
+                let factor = *new_matrix.get(r, i)?;
 
                 if factor.is_zero() {
                     continue;
@@ -2128,4 +2146,412 @@ mod tests {
             }
         }
     }
+
+    #[test]
+    fn test_partial_pivoting_no_swap_needed() {
+        let mut m = Matrix {
+            rows: 3,
+            columns: 3,
+            data: vec![
+                Fraction::from(5),
+                Fraction::from(0),
+                Fraction::from(0),
+                Fraction::from(2),
+                Fraction::from(0),
+                Fraction::from(0),
+                Fraction::from(1),
+                Fraction::from(0),
+                Fraction::from(0),
+            ],
+        };
+
+        let mut sign = Fraction::from(1);
+
+        let res = m.partial_pivoting(0, &mut sign).unwrap();
+
+        assert!(res);
+        assert_eq!(sign, Fraction::from(1)); // no swap
+    }
+
+    #[test]
+    fn test_partial_pivoting_with_swap() {
+        let mut m = Matrix {
+            rows: 3,
+            columns: 3,
+            data: vec![
+                Fraction::from(1),
+                Fraction::from(0),
+                Fraction::from(0),
+                Fraction::from(9),
+                Fraction::from(0),
+                Fraction::from(0),
+                Fraction::from(3),
+                Fraction::from(0),
+                Fraction::from(0),
+            ],
+        };
+
+        let mut sign = Fraction::from(1);
+
+        let res = m.partial_pivoting(0, &mut sign).unwrap();
+
+        assert!(res);
+
+        // ahora la fila 0 debe ser la de valor 9
+        assert_eq!(m.get(0, 0).unwrap(), &Fraction::from(9));
+
+        assert_eq!(sign, Fraction::from(-1)); // hubo swap
+    }
+
+    #[test]
+    fn test_partial_pivoting_all_zero_column() {
+        let mut m = Matrix {
+            rows: 3,
+            columns: 3,
+            data: vec![
+                Fraction::from(0),
+                Fraction::from(1),
+                Fraction::from(2),
+                Fraction::from(0),
+                Fraction::from(3),
+                Fraction::from(4),
+                Fraction::from(0),
+                Fraction::from(5),
+                Fraction::from(6),
+            ],
+        };
+
+        let mut sign = Fraction::from(1);
+
+        let res = m.partial_pivoting(0, &mut sign).unwrap();
+
+        assert!(!res); // no pivot posible
+    }
+
+    #[test]
+    fn test_partial_pivoting_largest_abs_value() {
+        let mut m = Matrix {
+            rows: 3,
+            columns: 3,
+            data: vec![
+                Fraction::from(-3),
+                Fraction::from(0),
+                Fraction::from(0),
+                Fraction::from(7),
+                Fraction::from(0),
+                Fraction::from(0),
+                Fraction::from(-10),
+                Fraction::from(0),
+                Fraction::from(0),
+            ],
+        };
+
+        let mut sign = Fraction::from(1);
+
+        m.partial_pivoting(0, &mut sign).unwrap();
+
+        assert_eq!(
+            m.get(0, 0).unwrap(),
+            &Fraction::from(-10) // mayor en valor absoluto
+        );
+    }
+
+    #[test]
+    fn test_gaussian_already_triangular() {
+        let m = Matrix {
+            rows: 3,
+            columns: 3,
+            data: vec![
+                Fraction::from(1),
+                Fraction::from(2),
+                Fraction::from(3),
+                Fraction::from(0),
+                Fraction::from(4),
+                Fraction::from(5),
+                Fraction::from(0),
+                Fraction::from(0),
+                Fraction::from(6),
+            ],
+        };
+
+        let (res, _) = m.gaussian_elimination().unwrap();
+
+        assert_eq!(res.data, m.data);
+    }
+
+    #[test]
+    fn test_gaussian_simple() {
+        let m = Matrix {
+            rows: 2,
+            columns: 2,
+            data: vec![
+                Fraction::from(2),
+                Fraction::from(1),
+                Fraction::from(4),
+                Fraction::from(3),
+            ],
+        };
+
+        let (res, _) = m.gaussian_elimination().unwrap();
+
+        // abajo izquierda debe ser 0
+        assert!(res.get(1, 0).unwrap().is_zero());
+    }
+
+    #[test]
+    fn test_gaussian_requires_pivoting() {
+        let m = Matrix {
+            rows: 2,
+            columns: 2,
+            data: vec![
+                Fraction::from(0),
+                Fraction::from(1),
+                Fraction::from(2),
+                Fraction::from(3),
+            ],
+        };
+
+        let (res, _) = m.gaussian_elimination().unwrap();
+
+        assert!(res.get(0, 0).unwrap() != &Fraction::from(0));
+    }
+
+    #[test]
+    fn test_gaussian_singular_matrix() {
+        let m = Matrix {
+            rows: 2,
+            columns: 2,
+            data: vec![
+                Fraction::from(1),
+                Fraction::from(2),
+                Fraction::from(2),
+                Fraction::from(4),
+            ],
+        };
+
+        let res = m.gaussian_elimination();
+
+        assert!(matches!(res, Err(MatrixError::FailedGauss)));
+    }
+
+    #[test]
+    fn test_gaussian_upper_triangular_property() {
+        let m = Matrix {
+            rows: 3,
+            columns: 3,
+            data: vec![
+                Fraction::from(2),
+                Fraction::from(1),
+                Fraction::from(3),
+                Fraction::from(4),
+                Fraction::from(1),
+                Fraction::from(6),
+                Fraction::from(6),
+                Fraction::from(2),
+                Fraction::from(9),
+            ],
+        };
+
+        let res = m.gaussian_elimination();
+
+        assert!(res.is_err());
+    }
+
+    #[test]
+    fn test_gj_identity_matrix() {
+        let m = Matrix {
+            rows: 3,
+            columns: 3,
+            data: vec![
+                Fraction::from(1),
+                Fraction::from(0),
+                Fraction::from(0),
+                Fraction::from(0),
+                Fraction::from(1),
+                Fraction::from(0),
+                Fraction::from(0),
+                Fraction::from(0),
+                Fraction::from(1),
+            ],
+        };
+
+        let res = m.gauss_jordan_elimination().unwrap();
+
+        assert_eq!(res.data, m.data);
+    }
+
+    #[test]
+    fn test_gj_already_reduced() {
+        let m = Matrix {
+            rows: 2,
+            columns: 2,
+            data: vec![
+                Fraction::from(1),
+                Fraction::from(0),
+                Fraction::from(0),
+                Fraction::from(1),
+            ],
+        };
+
+        let res = m.gauss_jordan_elimination().unwrap();
+
+        assert_eq!(res.data, m.data);
+    }
+
+    #[test]
+    fn test_gj_simple() {
+        let m = Matrix {
+            rows: 2,
+            columns: 2,
+            data: vec![
+                Fraction::from(2),
+                Fraction::from(1),
+                Fraction::from(4),
+                Fraction::from(3),
+            ],
+        };
+
+        let res = m.gauss_jordan_elimination().unwrap();
+
+        // debe quedar identidad
+        assert_eq!(
+            res.data,
+            vec![
+                Fraction::from(1),
+                Fraction::from(0),
+                Fraction::from(0),
+                Fraction::from(1),
+            ]
+        );
+    }
+
+    #[test]
+    fn test_gj_requires_pivoting() {
+        let m = Matrix {
+            rows: 2,
+            columns: 2,
+            data: vec![
+                Fraction::from(0),
+                Fraction::from(1),
+                Fraction::from(2),
+                Fraction::from(3),
+            ],
+        };
+
+        let res = m.gauss_jordan_elimination().unwrap();
+
+        // diagonal debe ser 1
+        assert_eq!(*res.get(0, 0).unwrap(), Fraction::from(1));
+        assert_eq!(*res.get(1, 1).unwrap(), Fraction::from(1));
+
+        // fuera de diagonal debe ser 0
+        assert!(res.get(0, 1).unwrap().is_zero());
+        assert!(res.get(1, 0).unwrap().is_zero());
+    }
+
+    #[test]
+    fn test_gj_3x3() {
+        let m = Matrix {
+            rows: 3,
+            columns: 3,
+            data: vec![
+                Fraction::from(1),
+                Fraction::from(2),
+                Fraction::from(3),
+                Fraction::from(0),
+                Fraction::from(1),
+                Fraction::from(4),
+                Fraction::from(5),
+                Fraction::from(6),
+                Fraction::from(0),
+            ],
+        };
+
+        let res = m.gauss_jordan_elimination().unwrap();
+
+        // debe quedar identidad
+        for i in 0..3 {
+            for j in 0..3 {
+                if i == j {
+                    assert_eq!(*res.get(i, j).unwrap(), Fraction::from(1));
+                } else {
+                    assert!(res.get(i, j).unwrap().is_zero());
+                }
+            }
+        }
+    }
+
+    #[test]
+    fn test_gj_singular_matrix() {
+        let m = Matrix {
+            rows: 2,
+            columns: 2,
+            data: vec![
+                Fraction::from(1),
+                Fraction::from(2),
+                Fraction::from(2),
+                Fraction::from(4),
+            ],
+        };
+
+        let res = m.gauss_jordan_elimination();
+
+        assert!(matches!(res, Err(MatrixError::FailedGaussJordan)));
+    }
+
+    #[test]
+    fn test_gj_reduced_row_echelon_form() {
+        let m = Matrix {
+            rows: 3,
+            columns: 3,
+            data: vec![
+                Fraction::from(2),
+                Fraction::from(1),
+                Fraction::from(-1),
+                Fraction::from(-3),
+                Fraction::from(-1),
+                Fraction::from(2),
+                Fraction::from(-2),
+                Fraction::from(1),
+                Fraction::from(2),
+            ],
+        };
+
+        let res = m.gauss_jordan_elimination().unwrap();
+
+        for i in 0..3 {
+            // pivote = 1
+            assert_eq!(*res.get(i, i).unwrap(), Fraction::from(1));
+
+            // columna del pivote = 0 excepto en pivote
+            for r in 0..3 {
+                if r != i {
+                    assert!(res.get(r, i).unwrap().is_zero());
+                }
+            }
+        }
+    }
+
+    #[test]
+    fn test_gj_rectangular_matrix() {
+        let m = Matrix {
+            rows: 2,
+            columns: 3,
+            data: vec![
+                Fraction::from(1),
+                Fraction::from(2),
+                Fraction::from(3),
+                Fraction::from(4),
+                Fraction::from(5),
+                Fraction::from(6),
+            ],
+        };
+
+        let res = m.gauss_jordan_elimination().unwrap();
+
+        // al menos debe cumplir forma escalonada reducida parcial
+        for i in 0..2 {
+            assert_eq!(*res.get(i, i).unwrap(), Fraction::from(1));
+        }
+    }
 }