1use crate::result::{MathError, Result};
2
3use crate::matrix::traits::{Identity, Matrix, Parseable};
4
5#[derive(Debug, Clone)]
6pub struct MatrixF32 {
7 pub content: Vec<Vec<f32>>,
8 tolerance: f32,
9}
10
11#[macro_export]
12macro_rules! matrix_f32 {
13 ($expression:tt, $tol:expr) => {
14 MatrixF32::parse($expression, $tol)
15 };
16}
17pub use matrix_f32;
18
19impl Matrix for MatrixF32 {
20 type T = f32;
21 fn columns(&self) -> usize {
22 self.content
23 .iter()
24 .next()
25 .map(|row| row.len())
26 .expect("There is no row") }
28
29 fn rows(&self) -> usize {
30 self.content.len()
31 }
32
33 fn is_square(&self) -> bool {
34 self.columns() == self.rows()
35 }
36
37 fn is_symmetric(&self) -> bool {
38 if !self.is_square() {
39 return false;
40 }
41 for i in 0..self.rows() {
42 for j in 0..self.columns() {
43 if self.get(i, j).unwrap() != self.get(j, i).unwrap() {
44 return false;
45 }
46 }
47 }
48 true
49 }
50
51 fn get(&self, row: usize, column: usize) -> Result<&Self::T> {
52 if row > self.rows() || column > self.columns() {
53 return Err(MathError::MatrixError(format!("Cannot get element from position ({row},{column}) if the matrix has {}x{} dimensions!", self.rows(), self.columns())));
54 }
55 let Some(matrix_row) = self.content.get(row) else {
56 return Err(MathError::MatrixError(format!("Cannot get row {row} if the matrix has {}x{} dimensions!", self.rows(), self.columns())));
57 };
58 matrix_row.get(column).ok_or(MathError::MatrixError(format!(
59 "Cannot get element from position ({row},{column}) if the matrix has {}x{} dimensions!",
60 self.rows(),
61 self.columns()
62 )))
63 }
64
65 fn get_mut(&mut self, row: usize, column: usize) -> Result<&mut Self::T> {
66 if row > self.rows() || column > self.columns() {
67 return Err(MathError::MatrixError(format!("Cannot get element from position ({row},{column}) if the matrix has {}x{} dimensions!", self.rows(), self.columns())));
68 }
69 let rows = self.rows();
70 let columns = self.columns();
71 let Some(matrix_row) = self.content.get_mut(row) else {
72 return Err(MathError::MatrixError(format!("Cannot get row {row} if the matrix has {}x{} dimensions!", rows, columns)));
73 };
74 matrix_row
75 .get_mut(column)
76 .ok_or(MathError::MatrixError(format!(
77 "Cannot get element from position ({row},{column}) if the matrix has {}x{} dimensions!",
78 rows,
79 columns
80 )))
81 }
82
83 fn set(&mut self, row: usize, column: usize, value: Self::T) -> Result<()> {
84 *self.get_mut(row, column)? = value;
85 Ok(())
86 }
87
88 fn swap_rows(&mut self, row1: usize, row2: usize) -> Result<()> {
89 if row1 > self.rows() || row2 > self.rows() {
90 return Err(MathError::MatrixError(format!(
91 "Cannot swap rows {row1} and {row2} if the matrix has {} rows!",
92 self.rows()
93 )));
94 }
95 self.content.swap(row1, row2);
96 Ok(())
97 }
98
99 fn transpose(&self) -> Self {
100 let mut new_matrix = self.clone(); for i in 0..self.rows() {
102 for j in 0..self.columns() {
103 new_matrix
104 .set(j, i, self.get(i, j).unwrap().clone())
105 .unwrap();
106 }
107 }
108 new_matrix
109 }
110
111 fn gaussian_triangulation(&self) -> Result<Self> {
112 let mut new_matrix = self.clone();
113 for i in 0..self.rows() {
114 let mut max_row = i;
115 for j in i + 1..self.rows() {
116 if new_matrix.get(j, i).unwrap().abs() > new_matrix.get(max_row, i).unwrap().abs() {
117 max_row = j;
118 }
119 }
120 new_matrix.swap_rows(i, max_row)?;
121 if i < self.rows() - 1 && new_matrix.get(i, i).unwrap().abs() <= new_matrix.tolerance()
122 {
123 return Err(MathError::MatrixError(
124 "Cannot perform gaussian elimination on a matrix with zero pivot".to_string(),
125 ));
126 }
127 for j in i + 1..self.rows() {
128 let factor = new_matrix.get(j, i).unwrap() / new_matrix.get(i, i).unwrap();
129 for k in i..self.columns() {
130 let new_value =
131 new_matrix.get(j, k).unwrap() - factor * new_matrix.get(i, k).unwrap();
132 new_matrix.set(j, k, new_value).unwrap();
133 }
134 }
135 }
136 Ok(new_matrix)
137 }
138
139 fn lu_decomposition(&self) -> Result<(Self, Self)> {
140 if !self.is_square() {
141 return Err(MathError::MatrixError(
142 "Cannot perform LU decomposition on a non-square matrix".to_string(),
143 ));
144 }
145 let mut lower = Self::id(self.rows(), self.tolerance());
146 let mut upper = self.clone();
147 for i in 0..self.rows() {
148 for j in 0..self.columns() {
149 if j < i {
150 let factor = upper.get(i, j).unwrap() / upper.get(j, j).unwrap();
151 lower.set(i, j, factor).unwrap();
152 for k in j..self.columns() {
153 let new_value =
154 upper.get(i, k).unwrap() - factor * upper.get(j, k).unwrap();
155 upper.set(i, k, new_value).unwrap();
156 }
157 }
158 }
159 }
160 Ok((lower, upper))
161 }
162
163 fn determinant_using_gauss(&self) -> Option<Self::T> {
164 let gaussian_elimination_result = self.gaussian_triangulation().ok()?;
165 let mut mult = Self::T::id(0, 0.0);
166 for i in 0..gaussian_elimination_result.rows() {
167 for j in 0..gaussian_elimination_result.columns() {
168 if i == j {
169 let value = gaussian_elimination_result.get(i, j).unwrap();
170 mult = mult * *value;
171 }
172 }
173 }
174 Some(mult)
175 }
176
177 fn determinant_using_lu(&self) -> Option<Self::T> {
178 let (_, upper) = self.lu_decomposition().ok()?;
179 let mut mult = Self::T::id(0, 0.0);
180 for i in 0..self.rows() {
181 for j in 0..self.columns() {
182 if i == j {
183 let upper_value = upper.get(i, j).unwrap();
184 mult = mult * *upper_value;
185 }
186 }
187 }
188 Some(mult)
189 }
190
191 fn cholesky_decomposition(&self) -> Result<Self> {
192 todo!("To be done");
193 }
194}
195
196impl MatrixF32 {
197 pub fn new(content: Vec<Vec<f32>>, tolerance: f32) -> Result<Self> {
198 let mut content_iter = content.iter();
199 if let Some(length) = content_iter.next().map(|row| row.len()) {
200 while let Some(next_length) = content_iter.next().map(|row| row.len()) {
201 if length != next_length {
202 return Err(MathError::MatrixError(
203 "Cannot build matrix from rows with different lenght".to_string(),
204 ));
205 }
206 }
207 }
208 Ok(Self { content, tolerance })
209 }
210
211 pub fn tolerance(&self) -> f32 {
212 self.tolerance
213 }
214}
215
216impl TryFrom<&str> for MatrixF32 {
217 fn try_from(value: &str) -> Result<Self> {
219 MatrixF32::parse(value, 1e-4)
220 }
221
222 type Error = MathError;
223}
224
225#[cfg(test)]
226mod test {
227 use crate::matrix::traits::{Matrix, Parseable, CheckedAdd};
228
229 use super::{matrix_f32, MatrixF32};
230 use pretty_assertions;
231
232 const TOLERANCE: f32 = 1e-10;
233
234 #[test]
235 fn get_1() {
236 let matrix = MatrixF32::new(vec![vec![1f32, 2f32], vec![3f32, 4f32]], TOLERANCE).unwrap();
237 let item = matrix.get(0, 1).unwrap();
238 pretty_assertions::assert_eq!(item.clone(), 2f32)
239 }
240
241 #[test]
242 fn get_2() {
243 let matrix = MatrixF32::new(
244 vec![
245 vec![1.1, 2.2, 3.3],
246 vec![4.4, 5.5, 6.6],
247 vec![7.7, 8.8, 9.9],
248 ],
249 TOLERANCE,
250 )
251 .unwrap();
252 vec![((0, 0), 1.1), ((0, 2), 3.3), ((1, 1), 5.5), ((2, 0), 7.7)]
253 .into_iter()
254 .for_each(|item| {
255 pretty_assertions::assert_eq!(
256 matrix.get(item.0 .0, item.0 .1).unwrap().clone(),
257 item.1
258 )
259 })
260 }
261
262 #[test]
263 fn equal_matrices() {
264 let matrix_a = MatrixF32::new(vec![vec![1f32, 1f32], vec![2f32, 2f32]], TOLERANCE)
265 .expect("Should've been able to built this matrix");
266 let matrix_b = MatrixF32::new(vec![vec![1f32, 1f32], vec![2f32, 2f32]], TOLERANCE)
267 .expect("Should've been able to built this matrix");
268 pretty_assertions::assert_eq!(matrix_a, matrix_b)
269 }
270
271 #[test]
272 fn create_matrix_from_macro_1() {
273 let matrix = matrix_f32!("{{1,2},{3,4}}", TOLERANCE)
274 .expect("Should have been able to build matrix from macro");
275 println!("{matrix}");
276 let other = matrix_f32!("{{1,2},{3,4}}", TOLERANCE).expect("asdf");
277 pretty_assertions::assert_eq!(matrix, other)
278 }
279
280 #[test]
281 fn sum_with_macro() {
282 let result = (matrix_f32!("{{1,1},{1,1}}", TOLERANCE)
283 .expect("asdf")
284 .checked_add(&matrix_f32!("{{2,2},{2,2}}", TOLERANCE).expect("asdf")))
285 .expect("asdf");
286 println!("{result}")
287 }
288
289 #[test]
290 fn gaussian_elimination_1() {
291 let matrix = matrix_f32!("{{1,2,3},{4,5,6},{7,8,9}}", TOLERANCE).expect("asdf");
292 let gauss = matrix.gaussian_triangulation().expect("asdf");
293 pretty_assertions::assert_eq!(
294 gauss,
295 matrix_f32!(
296 "{{+7.0000000000000, +8.0000000000000, +9.0000000000000},
297 {+0.0000000000000, +0.8571428060532, +1.7142856121063 },
298 {+0.0000000000000, +0.0000000000000, +0.0000000000000}}",
299 TOLERANCE
300 )
301 .expect("asdf")
302 );
303 }
304
305 #[test]
306 fn gaussian_elimination_2() {
307 let matrix =
308 matrix_f32!("{{1,2,1,2,1},{-1,2,-3,2,1},{0,1,-3,2,1}}", TOLERANCE).expect("asdf");
309 let gauss = matrix.gaussian_triangulation().expect("asdf");
310 pretty_assertions::assert_eq!(
311 gauss,
312 matrix_f32!("{{1,2,1,2,1},{0,4,-2,4,2},{0,0,-2.5,1,0.5}}", TOLERANCE).expect("asdf")
313 );
314 }
315
316 #[test]
317 fn lu_decomposition_1() {
318 let matrix = matrix_f32!("{{1,2,3},{4,5,6},{7,8,9}}", TOLERANCE).expect("asdf");
319 let (l, u) = matrix.lu_decomposition().expect("asdf");
320 pretty_assertions::assert_eq!(
321 l,
322 matrix_f32!("{{1,0,0},{4,1,0},{7,2,1}}", TOLERANCE).expect("asdf")
323 );
324 pretty_assertions::assert_eq!(
325 u,
326 matrix_f32!("{{1,2,3},{0,-3,-6},{0,0,0}}", TOLERANCE).expect("asdf")
327 );
328 }
329
330 #[test]
331 fn determinant_1() {
332 let matrix = matrix_f32!("{{1,2,3},{4,5,6},{7,8,9}}", TOLERANCE).expect("asdf");
333 let determinant = matrix.determinant_using_gauss().expect("asdf");
334 pretty_assertions::assert_eq!(determinant, 0f32);
335 }
336
337 #[test]
338 fn determinant_using_lu_1() {
339 let matrix = matrix_f32!("{{1,2,3},{4,5,6},{7,8,1}}", TOLERANCE).expect("asdf");
340 let (lower, upper) = matrix.lu_decomposition().expect("asdf");
341 println!("{lower}");
342 println!("{upper}");
343 let determinant = matrix.determinant_using_lu().expect("asdf");
344 pretty_assertions::assert_eq!(determinant, 24f32);
345 }
346}