1use rand::{distributions::Uniform, Rng};
2
3pub struct AlgebrustMatrix {
4 matrix: Vec<Vec<f64>>
5}
6
7impl AlgebrustMatrix {
8 pub fn new(matrix: &[&[f64]]) -> Self {
9 let matrix: Vec<Vec<f64>> = matrix
10 .iter()
11 .map(|vector| vector.to_vec())
12 .collect();
13 AlgebrustMatrix { matrix }
14 }
15
16 pub fn new_rand(size: (usize, usize), low: f64, high: f64) -> Self {
17 let mut rng = rand::thread_rng();
18 let range = Uniform::new(low, high);
19 let matrix: Vec<Vec<f64>> = (0..size.0)
20 .map(|_| (0..size.1).map(|_| rng.sample(&range)).collect())
21 .collect();
22 AlgebrustMatrix { matrix }
23 }
24
25 pub fn new_zeros(size: (usize, usize)) -> Self {
26 let matrix: Vec<Vec<f64>> = (0..size.0)
27 .map(|_| vec![0.0; size.1])
28 .collect();
29 AlgebrustMatrix { matrix }
30 }
31
32 pub fn new_identity(n: usize) -> Self {
33 let matrix = Self::get_identity_matrix(n);
34 AlgebrustMatrix { matrix }
35 }
36
37 pub fn size(&self) -> (usize, usize) {
38 (self.matrix.len(), self.matrix[0].len())
39 }
40
41 pub fn matrix_ref(&self) -> &Vec<Vec<f64>> {
42 &self.matrix
43 }
44
45 fn get_column(&self, column_index: usize, other: &AlgebrustMatrix) -> Vec<f64> {
46 other.matrix.iter().map(|row| row[column_index]).collect()
47 }
48
49 fn get_identity_matrix(n: usize) -> Vec<Vec<f64>> {
50 let mut matrix: Vec<Vec<f64>> = vec![vec![0.0; n]; n];
51 for i in 0..n { matrix[i][i] = 1.0; }
52 matrix
53 }
54
55 fn get_diagonal_product(&self, matrix: &Vec<Vec<f64>>, n: usize) -> f64 {
56 (0..n)
57 .map(|i| matrix[i][i])
58 .product()
59 }
60
61 fn forward_substitution(&self, l: &Vec<Vec<f64>>, b: &Vec<f64>) -> Vec<f64> {
62 let n = l.len();
63 let mut x = vec![0.0; n];
64 for i in 0..n {
65 let summation: f64 = (0..i)
66 .map(|j| l[i][j] * x[j])
67 .sum();
68 x[i] = (b[i] - summation) / l[i][i];
69 }
70
71 x
72 }
73
74 fn backward_substitution(&self, u: &Vec<Vec<f64>>, x: &Vec<f64>) -> Vec<f64> {
75 let n = u.len();
76 let mut y = vec![0.0; n];
77 for i in (0..n).rev() {
78 let summation: f64 = ((i + 1)..n)
79 .map(|j| u[i][j] * y[j])
80 .sum();
81 y[i] = (x[i] - summation) / u[i][i];
82 }
83
84 y
85 }
86
87 pub fn addition(&self, other: &AlgebrustMatrix) -> AlgebrustMatrix {
88 assert_eq!(self.size(), other.size());
89 let matrix: Vec<Vec<f64>> = self.matrix
90 .iter()
91 .zip(&other.matrix)
92 .map(|(vector1, vector2)| vector1
93 .iter()
94 .zip(vector2)
95 .map(|(a, b)| a + b)
96 .collect())
97 .collect();
98 AlgebrustMatrix { matrix }
99 }
100
101 pub fn subtraction(&self, other: &AlgebrustMatrix) -> AlgebrustMatrix {
102 assert_eq!(self.size(), other.size());
103 let matrix: Vec<Vec<f64>> = self.matrix
104 .iter()
105 .zip(&other.matrix)
106 .map(|(vector1, vector2)| vector1
107 .iter()
108 .zip(vector2)
109 .map(|(a, b)| a - b)
110 .collect())
111 .collect();
112 AlgebrustMatrix { matrix }
113 }
114
115 pub fn multiplication(&self, other: &AlgebrustMatrix) -> AlgebrustMatrix {
116 let (n, o) = self.size();
117 let (p, m) = other.size();
118 assert_eq!(o, p);
119 let mut matrix = vec![vec![0.0; m]; n];
120 for i in 0..n {
121 for j in 0..m {
122 matrix[i][j] = self.get_column(j, other)
123 .iter()
124 .zip(&self.matrix[i])
125 .map(|(a, b)| a * b)
126 .sum();
127 }
128 }
129 AlgebrustMatrix { matrix }
130 }
131
132 pub fn scalar_multiplication(&self, scaler: f64) -> AlgebrustMatrix {
133 let matrix: Vec<Vec<f64>> = self.matrix
134 .iter()
135 .map(|vector| vector
136 .iter()
137 .map(|a| a * scaler)
138 .collect())
139 .collect();
140 AlgebrustMatrix { matrix }
141 }
142
143 pub fn transpose(&self) -> AlgebrustMatrix {
144 let (n, m) = self.size();
145 let mut matrix: Vec<Vec<f64>> = vec![vec![0.0; n]; m];
146 for i in 0..n {
147 for j in 0..m {
148 matrix[j][i] = self.matrix[i][j];
149 }
150 }
151 AlgebrustMatrix { matrix }
152 }
153
154 pub fn lu_decomposition(&self) -> (AlgebrustMatrix, AlgebrustMatrix) {
155 let (n, m) = self.size();
156 assert_eq!(n, m);
157 let mut l: Vec<Vec<f64>> = Self::get_identity_matrix(n);
158 let mut u: Vec<Vec<f64>> = self.matrix.clone();
159 for j in 0..n {
160 for i in j+1..n {
161 l[i][j] = u[i][j] / u[j][j];
162 for k in j..n {
163 u[i][k] = u[i][k] - l[i][j] * u[j][k];
164 }
165 }
166 }
167 (AlgebrustMatrix{matrix: l}, AlgebrustMatrix{matrix: u})
168 }
169
170 pub fn determinant(&self) -> f64 {
171 let (n, m) = self.size();
172 assert_eq!(n, m);
173 let (l, u) = self.lu_decomposition();
174 let determinant_l: f64 = self.get_diagonal_product(l.matrix_ref(), n);
175 let determinant_u: f64 = self.get_diagonal_product(u.matrix_ref(), n);
176 determinant_l * determinant_u
177 }
178
179 pub fn inverse(&self) -> AlgebrustMatrix {
180 let (n, m) = self.size();
181 assert_eq!(n, m);
182 let (l, u) = self.lu_decomposition();
183 let identity_matrix = AlgebrustMatrix { matrix: Self::get_identity_matrix(n) };
184 let mut inverse_matrix = vec![vec![0.0; n]; n];
185 for i in 0..n {
186 let e_i = self.get_column(i, &identity_matrix);
187 let x_i = self.forward_substitution(l.matrix_ref(), &e_i);
188 let y_i = self.backward_substitution(u.matrix_ref(), &x_i);
189 for j in 0..n {
190 inverse_matrix[j][i] = y_i[j];
191 }
192 }
193 AlgebrustMatrix { matrix: inverse_matrix }
194 }
195}