From ccb4184fb33d521535531e82fdbc1e70adccc853 Mon Sep 17 00:00:00 2001
From: Jon Janzen <jonjanzen@me.com>
Date: Fri, 2 Apr 2021 17:15:04 -0600
Subject: [PATCH] inverse works but rel equal doesn't

---
 matrix/src/lib.rs | 83 ++++++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 82 insertions(+), 1 deletion(-)

diff --git a/matrix/src/lib.rs b/matrix/src/lib.rs
index 7e8bd48..5b35000 100644
--- a/matrix/src/lib.rs
+++ b/matrix/src/lib.rs
@@ -3,7 +3,7 @@ extern crate approx;
 
 use structs::Tuple;
 
-use std::ops::Index;
+use std::ops::{Index, IndexMut};
 
 #[derive(Debug)]
 pub struct Matrix {
@@ -97,6 +97,29 @@ impl Matrix {
         }
         Matrix::from_vec(m)
     }
+
+    pub fn is_invertable(&self) -> bool {
+        self.determinant() != 0.0
+    }
+
+    pub fn inverse(&self) -> Matrix {
+        // seems dangerous
+        if !self.is_invertable() {
+            panic!("We can't invert {:?}", self.matrix);
+        }
+
+        //let mut matrix: Vec<Vec<f32>> = Vec::with_capacity(self.matrix.len());
+        let mut matrix = Matrix::default(self.matrix.len(), self.matrix[0].len());
+        let det = self.determinant();
+        for (row_idx, row) in self.matrix.iter().enumerate().take(self.matrix.len()) {
+            for (col_idx, _) in row.iter().enumerate().take(row.len()) {
+                let c = self.cofactor(row_idx, col_idx);
+                let val = c / det;
+                matrix[col_idx][row_idx] = val;
+            }
+        }
+        matrix
+    }
 }
 
 impl Index<usize> for Matrix {
@@ -106,6 +129,12 @@ impl Index<usize> for Matrix {
     }
 }
 
+impl IndexMut<usize> for Matrix {
+    fn index_mut(&mut self, index: usize) -> &mut Self::Output {
+        &mut self.matrix[index]
+    }
+}
+
 impl PartialEq for Matrix {
     fn eq(&self, _rhs: &Self) -> bool {
         if self.matrix.len() != _rhs.matrix.len() {
@@ -479,4 +508,56 @@ mod tests {
         assert_eq!(-4071.0, m.determinant());
     }
 
+    #[test]
+    fn can_invert_invertable() {
+        let m = Matrix::from_array([
+            [6.0, 4.0, 4.0, 4.0],
+            [5.0, 5.0, 7.0, 6.0],
+            [4.0, -9.0, 3.0, -7.0],
+            [9.0, 1.0, 7.0, -6.0],
+        ]);
+
+        assert_eq!(-2120.0, m.determinant());
+        assert_eq!(true, m.is_invertable());
+    }
+
+    #[test]
+    fn can_invert_not_invertable() {
+        let m = Matrix::from_array([
+            [-4.0, 2.0, -2.0, -3.0],
+            [9.0, 6.0, 2.0, 6.0],
+            [0.0, -5.0, 1.0, -5.0],
+            [0.0, 0.0, 0.0, 0.0],
+        ]);
+
+        assert_eq!(0.0, m.determinant());
+        assert_eq!(false, m.is_invertable());
+
+    }
+
+    #[test]
+    fn inverse() {
+        let m = Matrix::from_array([
+            [-5.0, 2.0, 6.0, -8.0],
+            [1.0, -5.0, 1.0, 8.0],
+            [7.0, 7.0, -6.0, -7.0],
+            [1.0, -3.0, 7.0, 4.0],
+        ]);
+
+        let b = m.inverse();
+        assert_eq!(532.0, m.determinant());
+        assert_eq!(-160.0, m.cofactor(2, 3));
+        assert_eq!(-160.0/532.0, b[3][2]);
+        assert_eq!(105.0, m.cofactor(3, 2));
+        assert_eq!(105.0/532.0, b[2][3]);
+
+        let expected = Matrix::from_array([
+            [0.21805, 0.45113, 0.24060, -0.04511],
+            [-0.80827, -1.45677, -0.44361, 0.52068],
+            [-0.07895, -0.22368, -0.05263, 0.19737],
+            [-0.52256, -0.81392, -0.30075, 0.30639],
+        ]);
+
+        assert_eq!(expected, b);
+    }
 }