mathhook_core/matrices/
inverse_tests.rs

1#[cfg(test)]
2mod tests {
3    use crate::core::{Expression, Number};
4    use crate::matrices::{CoreMatrixOps, Matrix};
5    use num_traits::{One, Zero};
6
7    const EPSILON: f64 = 1e-10;
8
9    fn is_approx_zero(expr: &Expression) -> bool {
10        match expr {
11            Expression::Number(Number::Integer(n)) => *n == 0,
12            Expression::Number(Number::Float(f)) => f.abs() < EPSILON,
13            Expression::Number(Number::Rational(r)) => r.is_zero(),
14            _ => expr.is_zero(),
15        }
16    }
17
18    fn is_approx_one(expr: &Expression) -> bool {
19        match expr {
20            Expression::Number(Number::Integer(n)) => *n == 1,
21            Expression::Number(Number::Float(f)) => (*f - 1.0).abs() < EPSILON,
22            Expression::Number(Number::Rational(r)) => r.is_one(),
23            _ => false,
24        }
25    }
26
27    /// Test inverse computation for identity matrix
28    #[test]
29    fn test_identity_inverse() {
30        let identity = Matrix::identity(3);
31        let inverse = identity.inverse();
32
33        assert!(matches!(inverse, Matrix::Identity(_)));
34
35        let product = identity
36            .multiply(&inverse)
37            .expect("identity * identity should succeed");
38        assert!(matches!(product, Matrix::Identity(_)));
39    }
40
41    /// Test 2x2 matrix inverse using Gauss-Jordan elimination
42    #[test]
43    fn test_2x2_gauss_jordan_inverse() {
44        let matrix = Matrix::dense(vec![
45            vec![Expression::integer(1), Expression::integer(2)],
46            vec![Expression::integer(3), Expression::integer(4)],
47        ]);
48
49        let inverse = matrix.inverse();
50
51        let product = matrix
52            .multiply(&inverse)
53            .expect("matrix * inverse should succeed");
54
55        // A * A^(-1) should be identity - may be optimized to Identity/Diagonal
56        match product {
57            Matrix::Identity(ref id) => {
58                assert_eq!(id.size, 2, "Identity matrix should be 2x2");
59            }
60            Matrix::Diagonal(ref diag) => {
61                assert_eq!(diag.diagonal_elements.len(), 2);
62                for (i, elem) in diag.diagonal_elements.iter().enumerate() {
63                    assert!(
64                        is_approx_one(elem),
65                        "Diagonal element {} = {:?} should be ~1",
66                        i,
67                        elem
68                    );
69                }
70            }
71            Matrix::Dense(ref data) => {
72                assert_eq!(data.rows.len(), 2);
73                assert_eq!(data.rows[0].len(), 2);
74                assert!(
75                    is_approx_one(&data.rows[0][0]),
76                    "Diagonal [0,0] should be ~1"
77                );
78                assert!(
79                    is_approx_zero(&data.rows[0][1]),
80                    "Off-diagonal [0,1] should be ~0"
81                );
82                assert!(
83                    is_approx_zero(&data.rows[1][0]),
84                    "Off-diagonal [1,0] should be ~0"
85                );
86                assert!(
87                    is_approx_one(&data.rows[1][1]),
88                    "Diagonal [1,1] should be ~1"
89                );
90            }
91            _ => panic!(
92                "Expected identity, diagonal or dense matrix result, got {:?}",
93                product
94            ),
95        }
96    }
97
98    /// Test 3x3 matrix inverse using Gauss-Jordan elimination
99    #[test]
100    fn test_3x3_gauss_jordan_inverse() {
101        let matrix = Matrix::dense(vec![
102            vec![
103                Expression::integer(2),
104                Expression::integer(1),
105                Expression::integer(0),
106            ],
107            vec![
108                Expression::integer(1),
109                Expression::integer(2),
110                Expression::integer(1),
111            ],
112            vec![
113                Expression::integer(0),
114                Expression::integer(1),
115                Expression::integer(2),
116            ],
117        ]);
118
119        let inverse = matrix.inverse();
120
121        let product = matrix
122            .multiply(&inverse)
123            .expect("matrix * inverse should succeed");
124
125        // A * A^(-1) should be identity - may be optimized to Identity/Diagonal
126        match product {
127            Matrix::Identity(ref id) => {
128                assert_eq!(id.size, 3, "Identity matrix should be 3x3");
129            }
130            Matrix::Diagonal(ref diag) => {
131                assert_eq!(diag.diagonal_elements.len(), 3);
132                for (i, elem) in diag.diagonal_elements.iter().enumerate() {
133                    assert!(
134                        is_approx_one(elem),
135                        "Diagonal element {} = {:?} should be ~1",
136                        i,
137                        elem
138                    );
139                }
140            }
141            Matrix::Dense(ref data) => {
142                assert_eq!(data.rows.len(), 3);
143                assert_eq!(data.rows[0].len(), 3);
144                for i in 0..3 {
145                    for j in 0..3 {
146                        if i == j {
147                            assert!(
148                                is_approx_one(&data.rows[i][j]),
149                                "Diagonal [{},{}] = {:?} should be ~1",
150                                i,
151                                j,
152                                data.rows[i][j]
153                            );
154                        } else {
155                            assert!(
156                                is_approx_zero(&data.rows[i][j]),
157                                "Off-diagonal [{},{}] = {:?} should be ~0",
158                                i,
159                                j,
160                                data.rows[i][j]
161                            );
162                        }
163                    }
164                }
165            }
166            _ => panic!(
167                "Expected identity, diagonal or dense matrix result, got {:?}",
168                product
169            ),
170        }
171    }
172
173    /// Test singular matrix detection
174    #[test]
175    fn test_singular_matrix_detection() {
176        let singular = Matrix::dense(vec![
177            vec![Expression::integer(1), Expression::integer(2)],
178            vec![Expression::integer(2), Expression::integer(4)],
179        ]);
180
181        let inverse = singular.inverse();
182
183        if let Matrix::Dense(ref data) = inverse {
184            for row in &data.rows {
185                for elem in row {
186                    assert!(elem.is_zero());
187                }
188            }
189        }
190    }
191
192    /// Test (A^(-1))^(-1) = A
193    #[test]
194    fn test_inverse_of_inverse() {
195        let matrix = Matrix::dense(vec![
196            vec![Expression::integer(1), Expression::integer(2)],
197            vec![Expression::integer(3), Expression::integer(5)],
198        ]);
199
200        let inverse = matrix.inverse();
201        let double_inverse = inverse.inverse();
202
203        let product = matrix
204            .multiply(&double_inverse)
205            .expect("matrix * inverse should succeed");
206        let a_squared = matrix
207            .multiply(&matrix)
208            .expect("matrix * matrix should succeed");
209
210        if let (Matrix::Dense(prod_data), Matrix::Dense(a2_data)) = (&product, &a_squared) {
211            assert_eq!(prod_data.rows.len(), a2_data.rows.len());
212            assert_eq!(prod_data.rows[0].len(), a2_data.rows[0].len());
213        }
214    }
215
216    /// Test (AB)^(-1) = B^(-1)A^(-1)
217    #[test]
218    fn test_inverse_product_property() {
219        let a = Matrix::dense(vec![
220            vec![Expression::integer(1), Expression::integer(2)],
221            vec![Expression::integer(0), Expression::integer(1)],
222        ]);
223        let b = Matrix::dense(vec![
224            vec![Expression::integer(2), Expression::integer(0)],
225            vec![Expression::integer(1), Expression::integer(1)],
226        ]);
227
228        let ab = a.multiply(&b).expect("a * b should succeed");
229        let ab_inv = ab.inverse();
230
231        let a_inv = a.inverse();
232        let b_inv = b.inverse();
233        let b_inv_a_inv = b_inv
234            .multiply(&a_inv)
235            .expect("b_inv * a_inv should succeed");
236
237        let test1 = ab.multiply(&ab_inv).expect("ab * ab_inv should succeed");
238        let test2 = ab
239            .multiply(&b_inv_a_inv)
240            .expect("ab * b_inv_a_inv should succeed");
241
242        if let (Matrix::Dense(t1_data), Matrix::Dense(t2_data)) = (&test1, &test2) {
243            assert_eq!(t1_data.rows.len(), t2_data.rows.len());
244            assert_eq!(t1_data.rows[0].len(), t2_data.rows[0].len());
245        }
246    }
247
248    /// Test that matrix structure is preserved through inverse operations
249    #[test]
250    fn test_inverse_structure_preservation() {
251        let matrix = Matrix::dense(vec![
252            vec![Expression::integer(3), Expression::integer(1)],
253            vec![Expression::integer(2), Expression::integer(1)],
254        ]);
255
256        let inverse = matrix.inverse();
257
258        assert!(matches!(inverse, Matrix::Dense(_)));
259
260        let (orig_rows, orig_cols) = matrix.dimensions();
261        let (inv_rows, inv_cols) = inverse.dimensions();
262        assert_eq!(orig_rows, inv_rows);
263        assert_eq!(orig_cols, inv_cols);
264
265        let product = matrix
266            .multiply(&inverse)
267            .expect("matrix * inverse should succeed");
268        let (prod_rows, prod_cols) = product.dimensions();
269        assert_eq!(prod_rows, 2);
270        assert_eq!(prod_cols, 2);
271    }
272}
mathhook_core/matrices/inverse_tests.rs

mathhook_core/matrices/
inverse_tests.rs
