mathhook_core/core/matrix/
numeric_matrix.rs

1use crate::error::MathError;
2
3mod arithmetic;
4mod conversion;
5mod decomposition;
6mod display;
7mod multiply;
8mod solve;
9
10pub use decomposition::LUResult;
11
12const EPSILON: f64 = 1e-10;
13
14#[derive(Debug, Clone, PartialEq)]
15pub struct NumericMatrix {
16    rows: usize,
17    cols: usize,
18    data: Vec<f64>,
19}
20
21impl NumericMatrix {
22    pub fn zeros(rows: usize, cols: usize) -> Result<Self, MathError> {
23        if rows == 0 || cols == 0 {
24            return Err(MathError::DomainError {
25                operation: "NumericMatrix::zeros".to_string(),
26                value: crate::Expression::integer(0),
27                reason: "Matrix dimensions must be positive".to_string(),
28            });
29        }
30        Ok(Self {
31            rows,
32            cols,
33            data: vec![0.0; rows * cols],
34        })
35    }
36
37    pub fn identity(n: usize) -> Result<Self, MathError> {
38        if n == 0 {
39            return Err(MathError::DomainError {
40                operation: "NumericMatrix::identity".to_string(),
41                value: crate::Expression::integer(0),
42                reason: "Matrix dimension must be positive".to_string(),
43            });
44        }
45        let mut data = vec![0.0; n * n];
46        for i in 0..n {
47            data[i * n + i] = 1.0;
48        }
49        Ok(Self {
50            rows: n,
51            cols: n,
52            data,
53        })
54    }
55
56    pub fn from_flat(rows: usize, cols: usize, data: Vec<f64>) -> Result<Self, MathError> {
57        if rows == 0 || cols == 0 {
58            return Err(MathError::DomainError {
59                operation: "NumericMatrix::from_flat".to_string(),
60                value: crate::Expression::integer(0),
61                reason: "Matrix dimensions must be positive".to_string(),
62            });
63        }
64        if data.len() != rows * cols {
65            return Err(MathError::DomainError {
66                operation: "NumericMatrix::from_flat".to_string(),
67                value: crate::Expression::integer(data.len() as i64),
68                reason: format!(
69                    "Data length {} does not match dimensions {}x{}",
70                    data.len(),
71                    rows,
72                    cols
73                ),
74            });
75        }
76        Ok(Self { rows, cols, data })
77    }
78
79    pub fn from_fn<F>(rows: usize, cols: usize, mut f: F) -> Result<Self, MathError>
80    where
81        F: FnMut(usize, usize) -> f64,
82    {
83        if rows == 0 || cols == 0 {
84            return Err(MathError::DomainError {
85                operation: "NumericMatrix::from_fn".to_string(),
86                value: crate::Expression::integer(0),
87                reason: "Matrix dimensions must be positive".to_string(),
88            });
89        }
90        let mut data = Vec::with_capacity(rows * cols);
91        for i in 0..rows {
92            for j in 0..cols {
93                data.push(f(i, j));
94            }
95        }
96        Ok(Self { rows, cols, data })
97    }
98
99    pub fn dimensions(&self) -> (usize, usize) {
100        (self.rows, self.cols)
101    }
102
103    pub fn get(&self, row: usize, col: usize) -> Result<f64, MathError> {
104        if row >= self.rows || col >= self.cols {
105            return Err(MathError::DomainError {
106                operation: "NumericMatrix::get".to_string(),
107                value: crate::Expression::integer((row * self.cols + col) as i64),
108                reason: format!(
109                    "Index ({}, {}) out of bounds for {}x{} matrix",
110                    row, col, self.rows, self.cols
111                ),
112            });
113        }
114        Ok(self.data[row * self.cols + col])
115    }
116
117    pub fn set(&mut self, row: usize, col: usize, value: f64) -> Result<(), MathError> {
118        if row >= self.rows || col >= self.cols {
119            return Err(MathError::DomainError {
120                operation: "NumericMatrix::set".to_string(),
121                value: crate::Expression::integer((row * self.cols + col) as i64),
122                reason: format!(
123                    "Index ({}, {}) out of bounds for {}x{} matrix",
124                    row, col, self.rows, self.cols
125                ),
126            });
127        }
128        self.data[row * self.cols + col] = value;
129        Ok(())
130    }
131
132    pub fn is_square(&self) -> bool {
133        self.rows == self.cols
134    }
135
136    pub fn is_symmetric(&self) -> bool {
137        if !self.is_square() {
138            return false;
139        }
140        for i in 0..self.rows {
141            for j in (i + 1)..self.cols {
142                let diff = (self.data[i * self.cols + j] - self.data[j * self.cols + i]).abs();
143                if diff > EPSILON {
144                    return false;
145                }
146            }
147        }
148        true
149    }
150
151    pub fn transpose(&self) -> Self {
152        let mut data = vec![0.0; self.rows * self.cols];
153        for i in 0..self.rows {
154            for j in 0..self.cols {
155                data[j * self.rows + i] = self.data[i * self.cols + j];
156            }
157        }
158        Self {
159            rows: self.cols,
160            cols: self.rows,
161            data,
162        }
163    }
164}
165
166#[cfg(test)]
167mod tests {
168    use super::*;
169
170    #[test]
171    fn test_zeros() {
172        let m = NumericMatrix::zeros(2, 3).unwrap();
173        assert_eq!(m.dimensions(), (2, 3));
174        assert_eq!(m.get(0, 0).unwrap(), 0.0);
175        assert_eq!(m.get(1, 2).unwrap(), 0.0);
176    }
177
178    #[test]
179    fn test_identity() {
180        let m = NumericMatrix::identity(3).unwrap();
181        assert_eq!(m.dimensions(), (3, 3));
182        assert_eq!(m.get(0, 0).unwrap(), 1.0);
183        assert_eq!(m.get(1, 1).unwrap(), 1.0);
184        assert_eq!(m.get(2, 2).unwrap(), 1.0);
185        assert_eq!(m.get(0, 1).unwrap(), 0.0);
186    }
187
188    #[test]
189    fn test_from_flat() {
190        let data = vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
191        let m = NumericMatrix::from_flat(2, 3, data).unwrap();
192        assert_eq!(m.dimensions(), (2, 3));
193        assert_eq!(m.get(0, 0).unwrap(), 1.0);
194        assert_eq!(m.get(1, 2).unwrap(), 6.0);
195    }
196
197    #[test]
198    fn test_from_fn() {
199        let m = NumericMatrix::from_fn(2, 2, |i, j| (i * 2 + j) as f64).unwrap();
200        assert_eq!(m.get(0, 0).unwrap(), 0.0);
201        assert_eq!(m.get(0, 1).unwrap(), 1.0);
202        assert_eq!(m.get(1, 0).unwrap(), 2.0);
203        assert_eq!(m.get(1, 1).unwrap(), 3.0);
204    }
205
206    #[test]
207    fn test_transpose() {
208        let data = vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
209        let m = NumericMatrix::from_flat(2, 3, data).unwrap();
210        let t = m.transpose();
211        assert_eq!(t.dimensions(), (3, 2));
212        assert_eq!(t.get(0, 0).unwrap(), 1.0);
213        assert_eq!(t.get(1, 0).unwrap(), 2.0);
214        assert_eq!(t.get(2, 1).unwrap(), 6.0);
215    }
216
217    #[test]
218    fn test_is_symmetric() {
219        let sym_data = vec![1.0, 2.0, 2.0, 1.0];
220        let m = NumericMatrix::from_flat(2, 2, sym_data).unwrap();
221        assert!(m.is_symmetric());
222
223        let nonsym_data = vec![1.0, 2.0, 3.0, 1.0];
224        let m = NumericMatrix::from_flat(2, 2, nonsym_data).unwrap();
225        assert!(!m.is_symmetric());
226    }
227
228    #[test]
229    fn test_invalid_dimensions() {
230        assert!(NumericMatrix::zeros(0, 3).is_err());
231        assert!(NumericMatrix::identity(0).is_err());
232        assert!(NumericMatrix::from_flat(2, 3, vec![1.0, 2.0]).is_err());
233    }
234
235    #[test]
236    fn test_out_of_bounds() {
237        let m = NumericMatrix::zeros(2, 2).unwrap();
238        assert!(m.get(2, 0).is_err());
239        assert!(m.get(0, 2).is_err());
240    }
241}
mathhook_core/core/matrix/numeric_matrix.rs

mathhook_core/core/matrix/
numeric_matrix.rs
