inverse works but rel equal doesn't

submatrixes and determinant can be used in the same method now

matrix/src/lib.rs

@@ -1,36 +1,45 @@
-#![feature(const_generics)]
-#![feature(const_evaluatable_checked)]
-#![allow(incomplete_features)]
 #[macro_use]
 extern crate approx;
 
 use structs::Tuple;
 
-use std::ops::Index;
+use std::ops::{Index, IndexMut};
 
 #[derive(Debug)]
-pub struct Matrix<const H: usize, const W: usize> {
-    matrix: [[f32; W]; H],
+pub struct Matrix {
+    matrix: Vec<Vec<f32>>,
 }
 
-impl<const H: usize, const W: usize> Matrix<H, W> {
+impl Matrix {
 
-    pub fn default() -> Self {
+    pub fn default(width: usize, height: usize) -> Self {
         Matrix {
-            matrix: [[0f32; W]; H],
+            matrix: vec![vec![0.0f32; width]; height],
         }
     }
 
-    pub fn from_array(matrix: [[f32; W]; H]) -> Matrix<H, W> {
+    pub fn from_array<const H: usize, const W: usize>(array: [[f32; W]; H]) -> Matrix {
+        let mut matrix: Vec<Vec<f32>> = Vec::with_capacity(H);
+        for r in array.iter() {
+            let mut row: Vec<f32> = Vec::with_capacity(W);
+            for v in r.iter() {
+                row.push(*v);
+            }
+            matrix.push(row);
+        }
         Matrix {
             matrix,
         }
     }
 
-    pub fn identity() -> Matrix<H, W> {
-        // I can't figure out how to assign a 2d array to matrix inside the generic
-        // so I instead create the new and then assign 1.0 to the necessary values
-        let mut m = Self::default();
+    pub fn from_vec(matrix: Vec<Vec<f32>>) -> Matrix {
+        Matrix {
+            matrix,
+        }
+    }
+
+    pub fn identity(size: usize) -> Matrix {
+        let mut m = Self::default(size, size);
         for i in 0..m.matrix.len() {
             m.matrix[i][i] = 1.0;
         }
@@ -48,43 +57,85 @@ impl<const H: usize, const W: usize> Matrix<H, W> {
     }
 
     pub fn determinant(&self) -> f32 {
-        self.matrix[0][0] * self.matrix[1][1] - self.matrix[0][1] * self.matrix[1][0]
-    }
-
-    pub fn minor(&self, row: usize, col: usize) -> f32 where
-        [(); H - 1]: ,
-        [(); W - 1]: ,
-    {
-        self.sub_matrix(row, col).determinant()
-    }
-
-    pub fn sub_matrix(&self, skip_row: usize, skip_col: usize) -> Matrix<{H - 1}, {W - 1}>
-    {
-        let mut idx_row: usize = 0;
-
-        let mut arr = [[0f32; W - 1]; H - 1];
-        for (i, row) in self.matrix.iter().enumerate().take(H) {
-            if i == skip_row { continue; }
-            let mut idx_col: usize = 0;
-            for (j, col) in row.iter().enumerate().take(W) {
-                if j == skip_col { continue; }
-                arr[idx_row][idx_col] = *col;
-                idx_col += 1;
+        if self.matrix[0].len() == 2 {
+            self.matrix[0][0] * self.matrix[1][1] - self.matrix[0][1] * self.matrix[1][0]
+        } else {
+            let mut sum = 0.0;
+            for (col, val) in self.matrix[0].iter().enumerate().take(self.matrix[0].len()) {
+                sum += val * self.cofactor(0, col);
             }
-            idx_row += 1;
+            sum
         }
-        Matrix::from_array(arr)
+    }
+
+    pub fn minor(&self, row: usize, col: usize) -> f32 {
+        let m = self.sub_matrix(row, col);
+        let det = m.determinant();
+        det
+    }
+
+    pub fn cofactor(&self, row: usize, col: usize) -> f32 {
+        let minor = self.minor(row, col);
+        if (row + col) & 0x1 == 0 {
+            minor
+        } else {
+            minor * -1.0
+        }
+    }
+
+    pub fn sub_matrix(&self, skip_row: usize, skip_col: usize) -> Matrix
+    {
+        let mut m = Vec::<Vec<f32>>::with_capacity(self.matrix.len() - 1);
+        for (i, row) in self.matrix.iter().enumerate().take(self.matrix.len()) {
+            if i == skip_row { continue; }
+            let mut r = Vec::<f32>::with_capacity(row.len() - 1);
+            for (j, col) in row.iter().enumerate().take(row.len()) {
+                if j == skip_col { continue; }
+                r.push(*col);
+            }
+            m.push(r);
+        }
+        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<const H: usize, const W: usize> Index<usize> for Matrix<H, W> {
-    type Output = [f32; W];
+impl Index<usize> for Matrix {
+    type Output = Vec<f32>;
     fn index(&self, index: usize) -> &Self::Output {
         &self.matrix[index]
     }
 }
 
-impl<const H: usize, const W: usize> PartialEq for Matrix<H, W> {
+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() {
             return false;
@@ -104,10 +155,10 @@ impl<const H: usize, const W: usize> PartialEq for Matrix<H, W> {
     }
 }
 
-impl<const H: usize, const W: usize> Matrix<H, W> {
-    fn calc_val_for_mul(&self, row: usize, rhs: &Matrix<H, W>, col: usize) -> f32 {
+impl Matrix {
+    fn calc_val_for_mul(&self, row: usize, rhs: &Matrix, col: usize) -> f32 {
         let mut sum = 0.0;
-        for i in 0..W {
+        for i in 0..self.matrix.len() {
             sum += self.matrix[row][i] * rhs.matrix[i][col];
         }
         sum
@@ -121,22 +172,25 @@ impl<const H: usize, const W: usize> Matrix<H, W> {
     }
 }
 
-impl<const H: usize, const W: usize> std::ops::Mul<Matrix<H, W>> for Matrix<H, W> {
-    type Output = Matrix<H, W>;
+impl std::ops::Mul<Matrix> for Matrix {
+    type Output = Matrix;
 
-    fn mul(self, _rhs: Matrix<H, W>) -> Matrix<H, W> {
-        let mut result = [[0f32; W]; H];
-        for (row, val) in result.iter_mut().enumerate().take(H) {
-            for (col, v) in val.iter_mut().enumerate().take(W) {
-                *v = self.calc_val_for_mul(row, &_rhs, col);
+    fn mul(self, _rhs: Matrix) -> Matrix {
+        let mut result: Vec<Vec<f32>> = Vec::with_capacity(self.matrix.len());
+        for row in 0..self.matrix.len() {
+            let width = self.matrix[row].len();
+            let mut new_col = Vec::with_capacity(width);
+            for col in 0..width {
+                new_col.push( self.calc_val_for_mul(row, &_rhs, col));
             }
+            result.push(new_col);
         }
-        Matrix::from_array(result)
+        Matrix::from_vec(result)
     }
 
 }
 
-impl<const H: usize, const W: usize> std::ops::Mul<Tuple> for Matrix<H, W> {
+impl std::ops::Mul<Tuple> for Matrix {
     type Output = Tuple;
 
     fn mul(self, _rhs: Tuple) -> Tuple {
@@ -319,7 +373,7 @@ mod tests {
             [4.0, 8.0, 16.0, 32.0,]
         ]);
 
-        assert_eq!(matrix * Matrix::identity(), expected);
+        assert_eq!(matrix * Matrix::identity(4), expected);
     }
 
     #[test]
@@ -327,7 +381,7 @@ mod tests {
         let t = Tuple::new(1.0, 2.0, 3.0, 4.0);
         let expected = Tuple::new(1.0, 2.0, 3.0, 4.0);
 
-        assert_eq!(Matrix::<4, 4>::identity() * t, expected);
+        assert_eq!(Matrix::identity(4) * t, expected);
     }
 
     #[test]
@@ -351,9 +405,9 @@ mod tests {
 
     #[test]
     fn transpose_identity() {
-        let mut m = Matrix::identity();
+        let mut m = Matrix::identity(4);
         m.transpose();
-        assert_eq!(m, Matrix::<4, 4>::identity());
+        assert_eq!(m, Matrix::identity(4));
     }
 
     #[test]
@@ -412,4 +466,98 @@ mod tests {
         assert_eq!(25.0, s.determinant());
         assert_eq!(25.0, m.minor(1, 0));
     }
+
+    #[test]
+    fn cofactor_3x3() {
+        let m = Matrix::from_array([
+            [3.0, 5.0, 0.0],
+            [2.0, -1.0, -7.0],
+            [6.0, -1.0, 5.0],
+        ]);
+        assert_eq!(-12.0, m.minor(0, 0));
+        assert_eq!(-12.0, m.cofactor(0, 0));
+        assert_eq!(25.0, m.minor(1, 0));
+        assert_eq!(-25.0, m.cofactor(1, 0));
+    }
+
+    #[test]
+    fn determinant_3x3() {
+        let m = Matrix::from_array([
+            [1.0, 2.0, 6.0],
+            [-5.0, 8.0, -4.0],
+            [2.0, 6.0, 4.0],
+        ]);
+        assert_eq!(56.0, m.cofactor(0, 0));
+        assert_eq!(12.0, m.cofactor(0, 1));
+        assert_eq!(-46.0, m.cofactor(0, 2));
+        assert_eq!(-196.0, m.determinant());
+    }
+
+    #[test]
+    fn determinant_4x4() {
+        let m = Matrix::from_array([
+            [-2.0, -8.0, 3.0, 5.0],
+            [-3.0, 1.0, 7.0, 3.0],
+            [1.0, 2.0, -9.0, 6.0],
+            [-6.0, 7.0, 7.0, -9.0],
+        ]);
+        assert_eq!(690.0, m.cofactor(0, 0));
+        assert_eq!(447.0, m.cofactor(0, 1));
+        assert_eq!(210.0, m.cofactor(0, 2));
+        assert_eq!(51.0, m.cofactor(0, 3));
+        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);
+    }
 }