rust_matrix/
matrix_algebra.rs

1use crate::pow::Pow;
2use crate::sqrt::Sqrt;
3use std::fmt;
4use std::fmt::Display;
5use std::ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Sub};
6
7pub trait MatrixElementRequiredTraits<T>:
8    Add<Output = T>
9    + Sub<Output = T>
10    + Mul<Output = T>
11    + Div<Output = T>
12    + Add<Output = T>
13    + AddAssign
14    + MulAssign
15    + DivAssign
16    + Clone
17    + Copy
18    + Display
19    + PartialOrd<T>
20    + Neg<Output = T>
21    + Default
22    + Sqrt
23    + Pow
24    + From<u8>
25{
26}
27
28impl<
29        T: Add<Output = T>
30            + Sub<Output = T>
31            + Mul<Output = T>
32            + Div<Output = T>
33            + Add<Output = T>
34            + AddAssign
35            + MulAssign
36            + DivAssign
37            + Clone
38            + Copy
39            + Display
40            + PartialOrd<T>
41            + Neg<Output = T>
42            + Default
43            + Sqrt
44            + Pow
45            + From<u8>
46            + ?Sized,
47    > MatrixElementRequiredTraits<T> for T
48{
49}
50
51#[derive(PartialEq, Debug, PartialOrd, Clone)]
52pub struct Matrix<T: MatrixElementRequiredTraits<T>> {
53    pub n: usize,
54    pub m: usize,
55    pub entries: Vec<T>,
56}
57
58impl<T: MatrixElementRequiredTraits<T>> Matrix<T> {
59    /// Standard initializer for a Matrix
60    /// Takes dimensions `m` (number of rows) and `n` (number of columns)
61    /// together with a vector representing the elements. The number of elements
62    /// should be equal to `m` x `n`
63    /// ```
64    /// use matrix::matrix_algebra::Matrix;
65    /// let test_matrix = Matrix::new(2, 3, [1.0, 2.0, -1.0, 0.0, 3.0, 7.0].to_vec());
66    /// println!("{test_matrix}");
67    /// ```
68    /// Will output  
69    /// ⌜1    2    -1   ⌝  
70    /// ⌞0    3    7    ⌟
71    pub fn new(m: usize, n: usize, entries: Vec<T>) -> Self {
72        Matrix { m, n, entries }
73    }
74
75    /// Convenience initializer to initialize a matrix of dimension m x n with all entries
76    /// equal to the provided value
77    /// ```
78    /// use matrix::matrix_algebra::Matrix;
79    /// let test_matrix = Matrix::new_constant_value(2, 3, 5.0).expect("Unable to initialize matrix");
80    /// println!("{test_matrix}");
81    /// ```
82    /// /// Will output  
83    /// /// ⌜5    5    5    ⌝  
84    /// /// ⌞5    5    5    ⌟
85    pub fn new_constant_value(m: usize, n: usize, value: T) -> Result<Self, &'static str> {
86        if m == 0 || n == 0 {
87            return Err("Matrix dimensions must each be greater than zero!");
88        }
89
90        Ok(Matrix {
91            n,
92            m,
93            entries: vec![value; n * m],
94        })
95    }
96
97    /// Getter for entry at position i,j (zero indexed)
98    /// Returns a reference to the entry
99    /// ```
100    /// use matrix::matrix_algebra::Matrix;
101    /// let test_matrix = Matrix::new(2, 3, [1.0, 2.0, -1.0, 0.0, 3.0, 7.0].to_vec());
102    /// let entry = test_matrix.get_entry_ij(1, 1);
103    /// println!("{entry}");
104    /// ```
105    /// Will output  
106    /// 3  
107    pub fn get_entry_ij(&self, i: usize, j: usize) -> &T {
108        &self.entries[(i * self.n) + j]
109    }
110
111    /// Setter for entry at position i, j (zero indexed)
112    /// ```
113    /// use matrix::matrix_algebra::Matrix;
114    /// let mut test_matrix = Matrix::new(2, 3, [1.0, 2.0, -1.0, 0.0, 3.0, 7.0].to_vec());
115    /// let entry = test_matrix.get_entry_ij(1, 1);
116    /// println!("{entry}");
117    /// test_matrix.set_entry_ij(1, 1, &4.0);
118    /// let entry = test_matrix.get_entry_ij(1, 1);
119    /// println!("{entry}");
120    /// ```
121    /// Will output:  
122    /// 3  
123    /// 4  
124    pub fn set_entry_ij(&mut self, i: usize, j: usize, new_value: &T) {
125        self.entries[(i * self.n) + j] = new_value.clone();
126    }
127
128    /// Retrieves the rows of a matrix
129    /// ```
130    /// use matrix::matrix_algebra::Matrix;
131    /// let mut test_matrix = Matrix::new(2, 3, [1.0, 2.0, -1.0, 0.0, 3.0, 7.0].to_vec());
132    /// let rows = test_matrix.rows();
133    /// ```
134    pub fn rows(&self) -> Vec<Vec<T>> {
135        let mut rows = Vec::new();
136
137        for i in (0..(self.m * self.n)).step_by(self.n) {
138            rows.push(self.entries[i..(i + self.n)].to_vec());
139        }
140
141        rows
142    }
143
144    /// The first elementary operation on a matrix
145    /// Returns a new matrix with the rows at the specified
146    /// indices interchanged.
147    /// ```
148    /// use matrix::matrix_algebra::Matrix;
149    /// let test_matrix = Matrix::new(2, 3, [1.0, 2.0, 3.0, 4.0, 5.0, 6.0].to_vec());
150    /// let result = test_matrix.row_interchange(0, 1);
151    ///
152    /// assert_eq!(
153    ///    result,
154    ///    Matrix::new(2, 3, [4.0, 5.0, 6.0, 1.0, 2.0, 3.0].to_vec()),
155    /// );
156    /// ```
157    pub fn row_interchange(&self, first_row_index: usize, second_row_index: usize) -> Self {
158        let rows = self.rows();
159        let first_row = &rows[first_row_index];
160        let second_row = &rows[second_row_index];
161
162        let mut new_entries: Vec<T> = Vec::new();
163
164        for row_index in 0..rows.len() {
165            if row_index == first_row_index {
166                for entry in second_row {
167                    new_entries.push(*entry);
168                }
169            } else if row_index == second_row_index {
170                for entry in first_row {
171                    new_entries.push(*entry);
172                }
173            } else {
174                for entry in rows[row_index].clone() {
175                    new_entries.push(entry);
176                }
177            }
178        }
179
180        Matrix {
181            m: self.m,
182            n: self.n,
183            entries: new_entries,
184        }
185    }
186
187    /// The second elementary operation on a matrix
188    /// Returns a new matrix with all elements in the
189    ///  row at the specified index multiplied by the
190    /// specified scalar.
191    /// ```
192    /// use matrix::matrix_algebra::Matrix;
193    /// let test_matrix = Matrix::new(3, 3, [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0].to_vec());
194    /// let result = test_matrix.multiply_row_by_scalar(2, 2.0);
195    ///
196    /// assert_eq!(
197    ///    result,
198    ///    Matrix::new(
199    ///        3,
200    ///        3,
201    ///        [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 14.0, 16.0, 18.0].to_vec()
202    ///    ),
203    /// );
204    /// ```
205    pub fn multiply_row_by_scalar(&self, row_index: usize, scalar: T) -> Self {
206        let rows = self.rows();
207        let row_to_be_multiplied = &rows[row_index];
208
209        let mut new_entries: Vec<T> = Vec::new();
210
211        for i in 0..rows.len() {
212            if i == row_index {
213                for entry in row_to_be_multiplied {
214                    new_entries.push(*entry * scalar);
215                }
216            } else {
217                for entry in rows[i].clone() {
218                    new_entries.push(entry);
219                }
220            }
221        }
222
223        Matrix {
224            m: self.m,
225            n: self.n,
226            entries: new_entries,
227        }
228    }
229
230    /// The third elementary operation on a matrix
231    /// Adds a scalar multiple of a row to an existing row
232    /// in a matrix
233    /// ```
234    /// use matrix::matrix_algebra::Matrix;
235    /// let test_matrix = Matrix::new(3, 3, [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0].to_vec());
236    /// let result = test_matrix.add_row_to_scalar_multiple_of_row(0, 2, 5.0);
237    ///
238    /// assert_eq!(
239    ///     result,
240    ///     Matrix::new(
241    ///         3,
242    ///         3,
243    ///         [36.0, 42.0, 48.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0].to_vec()
244    ///     ),
245    /// );
246    /// ```
247    pub fn add_row_to_scalar_multiple_of_row(
248        &self,
249        target_index: usize,
250        source_index: usize,
251        scalar: T,
252    ) -> Self {
253        let rows = self.rows();
254        let row_to_be_added_to = rows[target_index].clone();
255        let row_to_be_added: Vec<T> = rows[source_index]
256            .clone()
257            .into_iter()
258            .map(|entry| -> T { scalar * entry })
259            .collect();
260
261        let mut new_entries: Vec<T> = Vec::new();
262
263        for i in 0..rows.len() {
264            if i == target_index {
265                let new_row = row_to_be_added.iter().zip(row_to_be_added_to.iter());
266                for (_, (lhs, rhs)) in new_row.enumerate() {
267                    new_entries.push(*lhs + *rhs);
268                }
269            } else {
270                for entry in rows[i].clone() {
271                    new_entries.push(entry);
272                }
273            }
274        }
275
276        Matrix {
277            m: self.m,
278            n: self.n,
279            entries: new_entries,
280        }
281    }
282
283    fn all_zeroes(&self) -> bool {
284        for i in 0..self.m {
285            for j in 0..self.n {
286                if *self.get_entry_ij(i, j) != T::default() {
287                    return false;
288                }
289            }
290        }
291
292        true
293    }
294
295    fn first_row_with_non_zero_entry_in_column(&self, column_index: usize) -> usize {
296        let rows = self.rows();
297        let mut first_row_with_non_zero_entry = 0;
298
299        for row in rows {
300            if row[column_index] != T::default() {
301                break;
302            }
303            first_row_with_non_zero_entry += 1;
304        }
305
306        first_row_with_non_zero_entry
307    }
308
309    fn reduce_rows_relative_to_row(&self, column_index: usize, row_index: usize) -> Self {
310        let mut reduced = self.clone();
311        for j in 0..self.m {
312            if j != row_index {
313                let value_in_non_zero_column_of_row = self.get_entry_ij(j, column_index);
314                if *value_in_non_zero_column_of_row != T::default() {
315                    let scalar = -T::from(1) * *value_in_non_zero_column_of_row;
316                    reduced = reduced.add_row_to_scalar_multiple_of_row(j, row_index, scalar);
317                }
318            }
319        }
320
321        reduced
322    }
323
324    fn reduce_first_row(&self, column_index: usize) -> (Self, T) {
325        let mut coefficient = T::from(1);
326        let first_row_with_non_zero_entry =
327            self.first_row_with_non_zero_entry_in_column(column_index);
328        let a_k_interchanged = self.row_interchange(0, first_row_with_non_zero_entry);
329        let scalar = T::from(1) / *a_k_interchanged.get_entry_ij(0, column_index);
330        coefficient = coefficient / scalar;
331        (
332            a_k_interchanged.multiply_row_by_scalar(0, scalar),
333            coefficient,
334        )
335    }
336
337    fn row_echolon_form_recursive_with_coefficient(&self, k: usize) -> (Self, T) {
338        let mut coefficient = T::from(1);
339        let mut partitioned: Vec<Matrix<T>> = Vec::new();
340
341        if k != 0 {
342            partitioned = self.partition(&[self.n], &[k, self.m - k]);
343        }
344
345        let a_k = if k == 0 { self } else { &partitioned[1] };
346
347        if a_k.all_zeroes() {
348            return (self.clone(), coefficient);
349        }
350
351        let columns = a_k.columns();
352        let first_non_zero_column_index = first_non_zero_vec_index(&columns);
353        let (mut a_k, new_coefficient) = a_k.reduce_first_row(first_non_zero_column_index);
354
355        coefficient = coefficient * -new_coefficient;
356
357        let combined_entries = if k == 0 {
358            &a_k.entries
359        } else {
360            let _ = &partitioned[0].entries.append(&mut a_k.entries);
361            &partitioned[0].entries
362        };
363
364        let combined = Matrix::new(self.m, self.n, combined_entries.to_vec());
365        let reduced = combined.reduce_rows_relative_to_row(first_non_zero_column_index, k);
366
367        if k != self.m - 1 {
368            let (row_echolon, new_coefficient) =
369                reduced.row_echolon_form_recursive_with_coefficient(k + 1);
370            return (row_echolon, new_coefficient * coefficient);
371        }
372
373        (reduced, coefficient)
374    }
375
376    fn row_echolon_form_recursive(&self, k: usize) -> Self {
377        self.row_echolon_form_recursive_with_coefficient(k).0
378    }
379
380    /// Reduces a matrix to Row Echolon Form
381    /// ```
382    /// use matrix::matrix_algebra::Matrix;
383    /// let test_matrix = Matrix::new(
384    ///     3,
385    ///     4,
386    ///     [
387    ///         0.0, 0.0, 3.0, -1.0, 0.0, -1.0, 4.0, 7.0, 0.0, -1.0, 7.0, 6.0,
388    ///     ]
389    ///     .to_vec(),
390    /// );
391    ///
392    /// let result = test_matrix.row_echolon_form();
393    ///
394    /// assert_eq!(
395    ///     result,
396    ///     Matrix::new(
397    ///         3,
398    ///         4,
399    ///         [
400    ///             0.0,
401    ///             1.0,
402    ///             0.0,
403    ///             -25.0 / 3.0,
404    ///             0.0,
405    ///             0.0,
406    ///             1.0,
407    ///             -1.0 / 3.0,
408    ///             0.0,
409    ///             0.0,
410    ///             0.0,
411    ///             0.0
412    ///         ]
413    ///         .to_vec()
414    ///     ),
415    /// );
416    /// ```
417    pub fn row_echolon_form(&self) -> Self {
418        self.row_echolon_form_recursive(0)
419    }
420
421    /// Reduces a matrix to column echolon form
422    /// ```
423    /// use matrix::matrix_algebra::Matrix;
424    /// let test_matrix = Matrix::new(2, 3, [1.0, 2.0, 3.0, 2.0, 3.0, 4.0].to_vec());
425    /// let result = test_matrix.column_echolon_form();
426    ///
427    /// assert_eq!(
428    ///     result,
429    ///     Matrix::new(2, 3, [1.0, 0.0, 0.0, 0.0, 1.0, 0.0,].to_vec())
430    /// );
431    /// let test_matrix = Matrix::new(3, 3, [1.0, 1.0, 1.0, 2.0, 2.0, 2.0, 3.0, 3.0, 3.0].to_vec());
432    /// let result = test_matrix.column_echolon_form();
433    /// assert_eq!(
434    ///     result,
435    ///     Matrix::new(3, 3, [1.0, 0.0, 0.0, 2.0, 0.0, 0.0, 3.0, 0.0, 0.0].to_vec())
436    /// );
437    /// ```
438    pub fn column_echolon_form(&self) -> Self {
439        self.transpose().row_echolon_form().transpose()
440    }
441
442    /// Convenience function that extracts the columns
443    /// of a matrix
444    /// ```
445    /// use matrix::matrix_algebra::Matrix;
446    /// let test_matrix = Matrix::new(
447    /// 3,
448    /// 4,
449    /// [
450    ///    1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0,
451    /// ]
452    /// .to_vec(),
453    /// );
454    ///
455    /// let columns = test_matrix.columns();
456    ///
457    /// assert_eq!(columns.len(), 4);
458    ///
459    /// assert_eq!(columns[0], vec![1.0, 5.0, 9.0]);
460    /// assert_eq!(columns[1], vec![2.0, 6.0, 10.0]);
461    /// assert_eq!(columns[2], vec![3.0, 7.0, 11.0]);
462    /// assert_eq!(columns[3], vec![4.0, 8.0, 12.0]);
463    ///
464    /// ```
465    pub fn columns(&self) -> Vec<Vec<T>> {
466        let mut columns = Vec::new();
467
468        for i in 0..self.n {
469            let mut new_column: Vec<T> = Vec::new();
470            for j in 0..self.m {
471                new_column.push(self.entries[(self.n * j) + i].clone());
472            }
473            columns.push(new_column.to_vec());
474        }
475
476        columns
477    }
478
479    /// Creates a transpose of a matrix
480    /// ```
481    /// use matrix::matrix_algebra::Matrix;
482    /// let test_matrix = Matrix::new(2, 3, [1.0, 2.0, -1.0, 0.0, 3.0, 7.0].to_vec());
483    ///
484    /// let transpose_matrix = test_matrix.transpose();
485    ///
486    /// assert_eq!(transpose_matrix.m, test_matrix.n);
487    /// assert_eq!(transpose_matrix.n, test_matrix.m);
488    /// assert_eq!(transpose_matrix.entries, [1.0, 0.0, 2.0, 3.0, -1.0, 7.0]);
489    /// ```
490    pub fn transpose(&self) -> Self {
491        let columns = self.columns();
492        let mut transposed_entries = Vec::new();
493
494        for column in columns {
495            for entry in column {
496                transposed_entries.push(entry);
497            }
498        }
499
500        Matrix::<T> {
501            m: self.n,
502            n: self.m,
503            entries: transposed_entries,
504        }
505    }
506
507    /// Calculates the determinant of a matrix
508    /// The following example shows that there can
509    /// be some rounding issues with the calculations.
510    /// ```
511    /// use matrix::matrix_algebra::Matrix;
512    /// use matrix::complex_number::ComplexNumber;
513    /// let test_matrix = Matrix::new(
514    ///     3,
515    ///     3,
516    ///     [
517    ///         ComplexNumber::new(1.0, 0.0),
518    ///         ComplexNumber::new(1.0, 0.0),
519    ///         ComplexNumber::new(0.0, 1.0),
520    ///         ComplexNumber::new(1.0, 1.0),
521    ///         ComplexNumber::new(1.0, 1.0),
522    ///         ComplexNumber::new(1.0, 0.0),
523    ///         ComplexNumber::new(2.0, 3.0),
524    ///         ComplexNumber::new(0.0, -1.0),
525    ///         ComplexNumber::new(3.0, 0.0),
526    ///     ]
527    ///     .to_vec(),
528    /// );
529    ///
530    /// let determinant = test_matrix.determinant();
531    ///
532    /// fn delta_real(determinant: ComplexNumber<f64>, expectation: f64) -> bool {
533    ///     determinant.real - expectation < (1.0 / 1000000000000000000000000000000.0)
534    /// }
535    ///
536    /// fn delta_complex(determinant: ComplexNumber<f64>, expectation: f64) -> bool {
537    ///     determinant.complex - expectation < (1.0 / 1000000000000000000000000000000.0)
538    /// }
539    /// assert!(delta_real(
540    ///     determinant.expect("Unable to calculate determinant"),
541    ///     8.0
542    /// ));
543    /// assert!(delta_complex(
544    ///     determinant.expect("Unable to calculate determinant"),
545    ///     6.0
546    /// ));
547    /// ```
548    pub fn determinant(&self) -> Result<T, &'static str> {
549        if self.n != self.m {
550            return Err("Non square matrices do not have determinants!");
551        }
552        let (row_echolon_form, coefficient) = self.row_echolon_form_recursive_with_coefficient(0);
553        let mut determinant: Option<T> = None;
554
555        for i in 0..row_echolon_form.n {
556            for j in 0..row_echolon_form.m {
557                if i == j {
558                    let entry = row_echolon_form.get_entry_ij(i, j);
559                    match determinant {
560                        None => determinant = Some(*entry),
561                        Some(det) => determinant = Some(det * *entry),
562                    }
563                }
564            }
565        }
566
567        match determinant {
568            Some(det) => return Ok(det * coefficient),
569            None => Err("Unable to calculate determinant!"),
570        }
571    }
572
573    /// Calculates the minor of a matrix (the determinant of
574    /// a given submatrix) based on the provided column indices
575    /// and row indices to be selected from the matrix
576    /// ```
577    /// use matrix::matrix_algebra::Matrix;
578    /// let test_matrix = Matrix::new(3, 3, [1.0, 2.0, 1.0, 6.0, -1.0, 0.0, -1.0, -2.0, -1.0].to_vec());
579    /// assert_eq!(test_matrix.minor(&[1, 2], &[1,2]).expect("Unable to create minor"), 1.0);
580    /// assert_eq!(test_matrix.minor(&[0, 2], &[1,2]).expect("Unable to create minor"), -6.0);
581    /// assert_eq!(test_matrix.minor(&[0, 1], &[0,1]).expect("Unable to create minor"), -13.0);
582    /// ```
583    pub fn minor(
584        &self,
585        column_indices: &[usize],
586        row_indices: &[usize],
587    ) -> Result<T, &'static str> {
588        let submatrix = self.submatrix(column_indices, row_indices);
589
590        submatrix.determinant()
591    }
592
593    /// Convenience function to identify a non-singular matrix
594    pub fn is_nonsingular(&self) -> Result<bool, &'static str> {
595        match self.determinant() {
596            Err(error) => Err(error),
597            Ok(determinant) => Ok(determinant != T::default()),
598        }
599    }
600
601    /// Creates a matrix of cofactors from a matrix
602    pub fn matrix_of_cofactors(&self) -> Result<Self, &'static str> {
603        let mut entries: Vec<T> = Vec::new();
604        let mut sign = T::from(1);
605        for j in 0..self.m {
606            for i in 0..self.n {
607                let mut column_indices: Vec<usize> = Vec::new();
608                for index in 0..self.n {
609                    if index != i {
610                        column_indices.push(index);
611                    }
612                }
613                let mut row_indices: Vec<usize> = Vec::new();
614                for index in 0..self.m {
615                    if index != j {
616                        row_indices.push(index);
617                    }
618                }
619                let minor = self.minor(&column_indices, &row_indices);
620                match minor {
621                    Ok(minor) => {
622                        entries.push(minor * sign);
623                        sign = -sign;
624                    }
625                    Err(err) => return Err(err),
626                }
627            }
628        }
629
630        Ok(Matrix {
631            n: self.n,
632            m: self.m,
633            entries,
634        })
635    }
636
637    /// Caluculates the adjoint of a matrix
638    /// ```
639    /// use matrix::matrix_algebra::Matrix;
640    /// let test_matrix = Matrix::new(3, 3, [3.0, 4.0, 3.0, 5.0, 7.0, 2.0, 0.0, 0.0, 1.0].to_vec());
641    ///
642    /// assert_eq!(
643    ///     test_matrix
644    ///         .adjoint()
645    ///         .expect("test_adjoint: unable to calculate matrix adjoint"),
646    ///     Matrix::new(
647    ///         3,
648    ///         3,
649    ///         [
650    ///             7.0,
651    ///             -4.0,
652    ///             -13.0,
653    ///             -5.0,
654    ///             3.0,
655    ///             9.0,
656    ///             0.0,
657    ///             0.0,
658    ///             1.0000000000000018
659    ///         ]
660    ///         .to_vec()
661    ///     )
662    /// );
663    /// ```
664    pub fn adjoint(&self) -> Result<Self, &'static str> {
665        match self.matrix_of_cofactors() {
666            Ok(matrix_of_cofactors) => Ok(matrix_of_cofactors.transpose()),
667            Err(err) => Err(err),
668        }
669    }
670
671    /// Calculates the inverse of a matrix
672    /// ```
673    /// use matrix::matrix_algebra::Matrix;
674    /// let test_matrix = Matrix::new(
675    ///     3,
676    ///     3,
677    ///     [2.0, 3.0, 0.0, 0.0, 3.0, -3.0, -2.0, 3.0, 3.0].to_vec(),
678    /// );
679    ///
680    /// let inverse = test_matrix
681    /// .inverse()
682    /// .expect("Unable to compute matrix inverse in test_inverse");
683    ///
684    /// assert_eq!(inverse.m, 3);
685    /// assert_eq!(inverse.n, 3);
686    ///
687    /// let expected = [
688    ///     1.0 / 3.0,
689    ///     -1.0 / 6.0,
690    ///     -1.0 / 6.0,
691    ///     1.0 / 9.0,
692    ///     1.0 / 9.0,
693    ///     1.0 / 9.0,
694    ///     1.0 / 9.0,
695    ///     -2.0 / 9.0,
696    ///     1.0 / 9.0,
697    ///     ]
698    ///     .to_vec();
699    ///     
700    ///     for i in 0..expected.len() {
701    ///         assert!(expected[i] - inverse.entries[i] < 1.0 / 1000000000000000.0);
702    ///     }
703    /// ```
704    pub fn inverse(&self) -> Result<Self, &'static str> {
705        if self.n != self.m {
706            return Err("Unable to compute inverse when n is not equal to m");
707        }
708        let identity_matrix = new_identity_matrix::<T>(self.m);
709
710        match identity_matrix {
711            Ok(identity_matrix) => {
712                let identity_matrix_rows = identity_matrix.rows();
713                let rows = self.rows();
714
715                let mut new_entries = Vec::new();
716
717                for i in 0..rows.len() {
718                    for element in &rows[i] {
719                        new_entries.push(*element);
720                    }
721                    for element in &identity_matrix_rows[i] {
722                        new_entries.push(*element);
723                    }
724                }
725                let enlarged = Matrix::new(self.m, 2 * self.n, new_entries);
726                let enlarged_ref = enlarged.row_echolon_form();
727
728                return Ok(enlarged_ref.submatrix(
729                    &Vec::from_iter(self.n..(2 * self.n)).as_slice(),
730                    &Vec::from_iter(0..self.m).as_slice(),
731                ));
732            }
733            Err(err) => return Err(err),
734        }
735    }
736
737    /// Partitions a matrix based on the provided column and row indices
738    /// ```
739    /// use matrix::matrix_algebra::Matrix;
740    /// let test_matrix = Matrix::new(
741    /// 5,
742    /// 5,
743    /// [
744    ///     1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0,
745    ///     16.0, 17.0, 18.0, 19.0, 20.0, 21.0, 22.0, 23.0, 24.0, 25.0,
746    /// ]
747    /// .to_vec(),
748    /// );
749    ///
750    /// let submatrices = test_matrix.partition(&[3, 2], &[2, 3]);
751    ///
752    /// assert_eq!(
753    ///     submatrices,
754    ///     [
755    ///         Matrix::new(2, 3, [1.0, 2.0, 3.0, 6.0, 7.0, 8.0].to_vec()),
756    ///         Matrix::new(2, 2, [4.0, 5.0, 9.0, 10.0].to_vec()),
757    ///         Matrix::new(
758    ///             3,
759    ///             3,
760    ///             [11.0, 12.0, 13.0, 16.0, 17.0, 18.0, 21.0, 22.0, 23.0].to_vec()
761    ///         ),
762    ///         Matrix::new(3, 2, [14.0, 15.0, 19.0, 20.0, 24.0, 25.0].to_vec())
763    ///     ]
764    /// );
765    /// ```
766    pub fn partition(
767        &self,
768        column_partitioning: &[usize],
769        row_partitioning: &[usize],
770    ) -> Vec<Self> {
771        let mut partitioned_matrices: Vec<Matrix<T>> = Vec::new();
772        let mut row_offset = 0;
773
774        for row_partition in row_partitioning {
775            let mut column_offset = 0;
776            for column_partition in column_partitioning {
777                let mut new_entries: Vec<T> = Vec::new();
778
779                for i in row_offset..(row_offset + *row_partition) {
780                    for j in column_offset..(column_offset + *column_partition) {
781                        new_entries.push(*self.get_entry_ij(i, j));
782                    }
783                }
784                partitioned_matrices.push(Matrix {
785                    m: *row_partition,
786                    n: *column_partition,
787                    entries: new_entries,
788                });
789                column_offset += *column_partition;
790            }
791            row_offset += *row_partition;
792        }
793
794        partitioned_matrices
795    }
796
797    /// Generates a submatrix based on the provided column
798    /// and row indices
799    /// ```
800    /// use matrix::matrix_algebra::Matrix;
801    /// let test_matrix = Matrix::new(
802    ///     5,
803    ///     5,
804    ///     [
805    ///         1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0,
806    ///         16.0, 17.0, 18.0, 19.0, 20.0, 21.0, 22.0, 23.0, 24.0, 25.0,
807    ///     ]
808    ///     .to_vec(),
809    /// );
810    ///
811    /// let submatrix = test_matrix.submatrix(&[1, 3], &[0, 2, 4]);
812    /// assert_eq!(submatrix.n, 2);
813    /// assert_eq!(submatrix.m, 3);
814    /// assert_eq!(submatrix.entries, [2.0, 4.0, 12.0, 14.0, 22.0, 24.0]);
815    /// ```
816    pub fn submatrix(&self, column_indices: &[usize], row_indices: &[usize]) -> Self {
817        let rows = self.rows();
818        let mut new_entries: Vec<T> = Vec::new();
819
820        for row_index in row_indices {
821            for column_index in column_indices {
822                let new_entry = &rows[*row_index][*column_index];
823
824                new_entries.push(*new_entry);
825            }
826        }
827
828        Matrix {
829            n: column_indices.len(),
830            m: row_indices.len(),
831            entries: new_entries,
832        }
833    }
834}
835
836/// Implementation of the Display trait for Matrix<T>
837/// which can be used for debugging purposes
838impl<T: MatrixElementRequiredTraits<T>> fmt::Display for Matrix<T> {
839    fn fmt(&self, fmt: &mut fmt::Formatter) -> fmt::Result {
840        let rows = self.rows();
841        for (i, row) in rows.clone().into_iter().enumerate() {
842            let printable_row: String = row
843                .iter()
844                .map(|&entry| {
845                    let mut entry_as_string = entry.to_string();
846
847                    while entry_as_string.len() < 4 {
848                        entry_as_string += " ";
849                    }
850
851                    entry_as_string + " "
852                })
853                .collect();
854            if i == 0 {
855                let _ = fmt.write_str("⌜");
856            } else if i == &rows.len() - 1 {
857                let _ = fmt.write_str("⌞");
858            } else {
859                let _ = fmt.write_str("⎸");
860            }
861
862            let _ = fmt.write_str(&printable_row)?;
863            if i == 0 {
864                let _ = fmt.write_str("⌝");
865            } else if i == rows.len() - 1 {
866                let _ = fmt.write_str("⌟");
867            } else {
868                let _ = fmt.write_str("⎹");
869            }
870            let _ = fmt.write_str("\n");
871        }
872        Ok(())
873    }
874}
875
876/// Convenience initialiser that creates a `Matrix<T>` using the default
877/// value of the type
878/// ```
879/// use matrix::matrix_algebra::{Matrix, new_all_default};
880/// use rand::prelude::*;
881/// let mut rng = rand::thread_rng();
882/// let n = rng.gen_range(1..=100);
883/// let m = rng.gen_range(1..=100);
884/// let test_matrix_f64 = new_all_default::<f64>(m, n)
885///     .expect("test_new_all_default: unable to create test_matrix_f64");
886///
887/// for entry in test_matrix_f64.entries {
888///     assert_eq!(entry, f64::default());
889/// }
890///
891/// let test_matrix_f32 = new_all_default::<f32>(m, n)
892///     .expect("test_new_all_default: unable to create test_matrix_f64");
893///
894/// for entry in test_matrix_f32.entries {
895///     assert_eq!(entry, f32::default());
896/// }
897/// ```
898pub fn new_all_default<T>(m: usize, n: usize) -> Result<Matrix<T>, &'static str>
899where
900    T: MatrixElementRequiredTraits<T>,
901{
902    if m == 0 || n == 0 {
903        return Err("Matrix dimensions must each be greater than zero!");
904    }
905
906    Ok(Matrix::<T> {
907        n,
908        m,
909        entries: vec![T::default(); n * m],
910    })
911}
912
913/// Creates an identity matrix of the specified dimension using whatever
914/// corresponds to unity for type T
915/// ```
916/// use matrix::matrix_algebra::{Matrix, new_identity_matrix};
917/// use rand::prelude::*;
918/// let mut rng = rand::thread_rng();
919/// let dimension = rng.gen_range(1..=100);
920///
921/// let test_matrix_f64 = new_identity_matrix::<f64>(dimension)
922///     .expect("test_new_identity_matrix: unable to create test_matrix_f64");
923///
924/// for index in 0..dimension {
925///     assert_eq!(*test_matrix_f64.get_entry_ij(index, index), 1.0);
926/// }
927///
928/// let test_matrix_f32 = new_identity_matrix::<f32>(dimension)
929///     .expect("test_new_identity_matrix: unable to create test_matrix_f64");
930///
931/// for index in 0..dimension {
932///     assert_eq!(*test_matrix_f32.get_entry_ij(index, index), 1.0);
933/// }
934/// ```
935pub fn new_identity_matrix<T>(dimension: usize) -> Result<Matrix<T>, &'static str>
936where
937    T: MatrixElementRequiredTraits<T>,
938{
939    let new_empty_matrix = new_all_default::<T>(dimension, dimension);
940
941    match new_empty_matrix {
942        Err(err) => return Err(err),
943        Ok(mut new_empty_matrix) => {
944            for index in 0..dimension {
945                new_empty_matrix.set_entry_ij(index, index, &T::from(1));
946            }
947            Ok(new_empty_matrix)
948        }
949    }
950}
951
952/// Implementation of matrix addition
953/// ```
954/// use matrix::matrix_algebra::Matrix;
955/// let test_matrix_a = Matrix::new(
956///     3,
957///     4,
958///     [
959///         1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0,
960///     ]
961///     .to_vec(),
962/// );
963/// let test_matrix_b = Matrix::new(
964///     3,
965///     4,
966///     [
967///         12.0, 11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 1.0,
968///     ]
969///     .to_vec(),
970/// );
971///
972/// let matrix_sum = test_matrix_a + test_matrix_b;
973///
974/// assert_eq!(matrix_sum.entries.len(), 12);
975/// assert_eq!(
976///     matrix_sum.entries,
977///     [13.0, 13.0, 13.0, 13.0, 13.0, 13.0, 13.0, 13.0, 13.0, 13.0, 13.0, 13.0]
978/// );
979/// ```
980fn matrix_add<T: MatrixElementRequiredTraits<T>>(
981    a: &Matrix<T>,
982    b: &Matrix<T>,
983) -> Result<Matrix<T>, &'static str>
984where
985    for<'a> &'a T: Add<Output = T>,
986{
987    if !is_additively_conformable(a, b) {
988        return Err("Matrices are not additively conformable!");
989    }
990    let mut entries: Vec<T> = Vec::new();
991
992    for i in 0..a.m {
993        for j in 0..a.n {
994            let new_value = a.get_entry_ij(i, j) + b.get_entry_ij(i, j);
995            entries.push(new_value);
996        }
997    }
998
999    Ok(Matrix::<T> {
1000        n: a.m,
1001        m: a.n,
1002        entries,
1003    })
1004}
1005
1006/// Implementation of the `Add` trait for `Matrix<T>`
1007/// Uses `matrix_add`
1008impl<T: MatrixElementRequiredTraits<T>> Add for Matrix<T>
1009where
1010    for<'a> &'a T: Add<Output = T>,
1011{
1012    type Output = Self;
1013
1014    fn add(self, other: Self) -> Self {
1015        match matrix_add(&self, &other) {
1016            Ok(result) => return result,
1017            Err(err) => panic!("{err}"),
1018        }
1019    }
1020}
1021
1022/// Implementation of the `Mul` trait for `Matrix<T>`
1023/// Uses `matrix_multiply`
1024/// ```
1025/// use matrix::matrix_algebra::Matrix;
1026/// let test_matrix_a = Matrix::new(
1027///     3,
1028///     4,
1029///     [
1030///         1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0,
1031///     ]
1032///     .to_vec(),
1033/// );
1034/// let test_matrix_b = Matrix::new(
1035///     4,
1036///     3,
1037///     [
1038///         12.0, 11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 1.0,
1039///     ]
1040///     .to_vec(),
1041/// );
1042///
1043/// let matrix_product = test_matrix_a * test_matrix_b;
1044///
1045/// assert_eq!(matrix_product.entries.len(), 9);
1046///
1047/// assert_eq!(
1048/// matrix_product.entries,
1049///     [
1050///         1.0 * 12.0 + 2.0 * 9.0 + 3.0 * 6.0 + 4.0 * 3.0,
1051///         1.0 * 11.0 + 2.0 * 8.0 + 3.0 * 5.0 + 4.0 * 2.0,
1052///         1.0 * 10.0 + 2.0 * 7.0 + 3.0 * 4.0 + 4.0 * 1.0,
1053///         5.0 * 12.0 + 6.0 * 9.0 + 7.0 * 6.0 + 8.0 * 3.0,
1054///         5.0 * 11.0 + 6.0 * 8.0 + 7.0 * 5.0 + 8.0 * 2.0,
1055///         5.0 * 10.0 + 6.0 * 7.0 + 7.0 * 4.0 + 8.0 * 1.0,
1056///         9.0 * 12.0 + 10.0 * 9.0 + 11.0 * 6.0 + 12.0 * 3.0,
1057///         9.0 * 11.0 + 10.0 * 8.0 + 11.0 * 5.0 + 12.0 * 2.0,
1058///         9.0 * 10.0 + 10.0 * 7.0 + 11.0 * 4.0 + 12.0 * 1.0
1059///     ]
1060/// );
1061/// ```
1062impl<T: MatrixElementRequiredTraits<T>> Mul for Matrix<T> {
1063    type Output = Self;
1064
1065    fn mul(self, rhs: Self) -> Self {
1066        match matrix_multiply::<T>(&self, &rhs) {
1067            Err(err) => panic!("{err}"),
1068            Ok(result) => result,
1069        }
1070    }
1071}
1072
1073/// Helper function to check that two matrices
1074/// can be multiplied
1075fn is_multiplicatively_conformable<T: MatrixElementRequiredTraits<T>>(
1076    a: &Matrix<T>,
1077    b: &Matrix<T>,
1078) -> bool {
1079    a.n == b.m
1080}
1081
1082/// Helper function to check that two matrices
1083/// can be summed
1084fn is_additively_conformable<T: MatrixElementRequiredTraits<T>>(
1085    a: &Matrix<T>,
1086    b: &Matrix<T>,
1087) -> bool {
1088    a.n == b.n && a.m == b.m
1089}
1090
1091/// Multiplication of two matrices
1092/// This is currently a naive implementation with no
1093/// optimisations.
1094fn matrix_multiply<T: MatrixElementRequiredTraits<T>>(
1095    a: &Matrix<T>,
1096    b: &Matrix<T>,
1097) -> Result<Matrix<T>, &'static str> {
1098    if !is_multiplicatively_conformable(&a, &b) {
1099        return Err("Matrices are not multiplicatively conformable!");
1100    }
1101    let mut entries: Vec<T> = Vec::new();
1102    let a_rows = a.rows();
1103    let b_columns = b.columns();
1104
1105    for i in 0..a_rows.len() {
1106        for j in 0..b_columns.len() {
1107            let mut new_entry: Option<T> = None;
1108            for k in 0..a_rows[i].len() {
1109                match new_entry {
1110                    Some(value) => new_entry = Some(value + a_rows[i][k] * b_columns[j][k]),
1111                    None => new_entry = Some(a_rows[i][k] * b_columns[j][k]),
1112                }
1113            }
1114
1115            if entries.len() == 0 {
1116                entries = vec![new_entry.expect("new_entry not initialised!")];
1117            } else {
1118                entries.push(new_entry.expect("new_entry not initialised!"));
1119            }
1120        }
1121    }
1122
1123    Ok(Matrix::<T> {
1124        m: a_rows.len(),
1125        n: b_columns.len(),
1126        entries,
1127    })
1128}
1129
1130fn first_non_zero_vec_index<T: MatrixElementRequiredTraits<T>>(input: &Vec<Vec<T>>) -> usize {
1131    let mut first_nonzero_vec_index = 0;
1132
1133    for vector in input {
1134        for element in vector {
1135            if *element != T::default() {
1136                return first_nonzero_vec_index;
1137            }
1138        }
1139        first_nonzero_vec_index += 1;
1140    }
1141
1142    first_nonzero_vec_index
1143}
1144
1145#[cfg(test)]
1146mod tests {
1147    use std::ops::Add;
1148
1149    use rand::prelude::*;
1150
1151    use crate::{complex_number::ComplexNumber, matrix_algebra::first_non_zero_vec_index};
1152
1153    use super::{new_all_default, new_identity_matrix, Matrix, MatrixElementRequiredTraits};
1154
1155    fn index_of_first_non_zero_element_in_vec<T: MatrixElementRequiredTraits<T>>(
1156        input: &Vec<T>,
1157    ) -> usize {
1158        let mut index_of_first_non_zero_element = 0;
1159        for element in input {
1160            if *element != T::default() {
1161                return index_of_first_non_zero_element;
1162            }
1163            index_of_first_non_zero_element += 1;
1164        }
1165
1166        index_of_first_non_zero_element
1167    }
1168
1169    #[test]
1170    fn test_constant_value_initialiser() {
1171        let mut rng = rand::thread_rng();
1172        let n = rng.gen_range(1..=100);
1173        let m = rng.gen_range(1..=100);
1174        let value: f64 = rng.gen_range(1.0..=100.0);
1175        let test_matrix =
1176            Matrix::new_constant_value(m, n, value).expect("Unable to create test_matrix");
1177
1178        assert_eq!(test_matrix.entries.len(), n * m);
1179
1180        for entry in test_matrix.entries {
1181            assert_eq!(entry, value);
1182        }
1183    }
1184
1185    #[test]
1186    fn test_initialiser() {
1187        let mut rng = rand::thread_rng();
1188        let n = rng.gen_range(1..=100);
1189        let m = rng.gen_range(1..=100);
1190        let mut entries: Vec<f32> = Vec::new();
1191
1192        for _i in 0..m {
1193            for _j in 0..n {
1194                entries.push(rng.gen_range(1.0..=100.0));
1195            }
1196        }
1197
1198        let test_matrix = Matrix { m, n, entries };
1199
1200        assert_eq!(test_matrix.entries.len(), n * m);
1201    }
1202
1203    #[test]
1204    fn test_all_zeroes_initialiser() {
1205        let mut rng = rand::thread_rng();
1206        let n = rng.gen_range(1..=100);
1207        let m = rng.gen_range(1..=100);
1208        let test_matrix =
1209            Matrix::new_constant_value(m, n, 0.0).expect("Unable to create test_matrix");
1210
1211        assert_eq!(test_matrix.entries.len(), n * m);
1212
1213        for entry in test_matrix.entries {
1214            assert_eq!(entry, 0.0);
1215        }
1216    }
1217
1218    #[test]
1219    fn test_zero_first_argument_to_initialiser() {
1220        let test_matrix = Matrix::new_constant_value(0, 1, 1.0);
1221        assert_eq!(
1222            test_matrix,
1223            Err("Matrix dimensions must each be greater than zero!")
1224        );
1225    }
1226
1227    #[test]
1228    fn test_zero_second_argument_to_initialiser() {
1229        let test_matrix = Matrix::new_constant_value(1, 0, 1.0);
1230        assert_eq!(
1231            test_matrix,
1232            Err("Matrix dimensions must each be greater than zero!")
1233        );
1234    }
1235
1236    #[test]
1237    #[should_panic]
1238    fn test_panic_on_non_multiplicatively_conformable_matrices() {
1239        let test_matrix_a =
1240            Matrix::new_constant_value(3, 4, 5.0).expect("Unable to create test_matrix_a");
1241        let test_matrix_b =
1242            Matrix::new_constant_value(5, 7, 4.0).expect("Unable to create test_matrix_b");
1243
1244        let _ = test_matrix_a * test_matrix_b;
1245    }
1246
1247    #[test]
1248    #[should_panic]
1249    fn test_panic_on_non_additively_conformable_matrices() {
1250        let test_matrix_a =
1251            Matrix::new_constant_value(3, 4, 5.0).expect("Unable to create test_matrix_a");
1252        let test_matrix_b =
1253            Matrix::new_constant_value(5, 7, 4.0).expect("Unable to create test_matrix_b");
1254
1255        let _ = test_matrix_a.add(test_matrix_b);
1256    }
1257
1258    #[test]
1259    fn test_matrix_add() {
1260        let test_matrix_a = Matrix::new(
1261            3,
1262            4,
1263            [
1264                1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0,
1265            ]
1266            .to_vec(),
1267        );
1268        let test_matrix_b = Matrix::new(
1269            3,
1270            4,
1271            [
1272                12.0, 11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 1.0,
1273            ]
1274            .to_vec(),
1275        );
1276
1277        let matrix_sum = test_matrix_a + test_matrix_b;
1278
1279        assert_eq!(matrix_sum.entries.len(), 12);
1280        assert_eq!(
1281            matrix_sum.entries,
1282            [13.0, 13.0, 13.0, 13.0, 13.0, 13.0, 13.0, 13.0, 13.0, 13.0, 13.0, 13.0]
1283        );
1284    }
1285
1286    #[test]
1287    fn test_matrix_add_operator() {
1288        let test_matrix_a = Matrix::new(
1289            3,
1290            4,
1291            [
1292                1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0,
1293            ]
1294            .to_vec(),
1295        );
1296        let test_matrix_b = Matrix::new(
1297            3,
1298            4,
1299            [
1300                12.0, 11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 1.0,
1301            ]
1302            .to_vec(),
1303        );
1304        let matrix_sum = test_matrix_a + test_matrix_b;
1305
1306        assert_eq!(matrix_sum.entries.len(), 12);
1307        assert_eq!(
1308            matrix_sum.entries,
1309            [13.0, 13.0, 13.0, 13.0, 13.0, 13.0, 13.0, 13.0, 13.0, 13.0, 13.0, 13.0,]
1310        );
1311    }
1312
1313    #[test]
1314    fn test_columns() {
1315        let test_matrix = Matrix::new(
1316            3,
1317            4,
1318            [
1319                1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0,
1320            ]
1321            .to_vec(),
1322        );
1323
1324        let columns = test_matrix.columns();
1325
1326        assert_eq!(columns.len(), 4);
1327
1328        assert_eq!(columns[0], vec![1.0, 5.0, 9.0]);
1329        assert_eq!(columns[1], vec![2.0, 6.0, 10.0]);
1330        assert_eq!(columns[2], vec![3.0, 7.0, 11.0]);
1331        assert_eq!(columns[3], vec![4.0, 8.0, 12.0]);
1332    }
1333
1334    #[test]
1335    fn test_rows() {
1336        let test_matrix = Matrix::new(
1337            3,
1338            4,
1339            [
1340                1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0,
1341            ]
1342            .to_vec(),
1343        );
1344
1345        let rows = test_matrix.rows();
1346
1347        assert_eq!(rows.len(), 3);
1348
1349        assert_eq!(rows[0], vec![1.0, 2.0, 3.0, 4.0,]);
1350        assert_eq!(rows[1], vec![5.0, 6.0, 7.0, 8.0,]);
1351        assert_eq!(rows[2], vec![9.0, 10.0, 11.0, 12.0,]);
1352    }
1353
1354    #[test]
1355    fn test_matrix_multiply() {
1356        let test_matrix_a = Matrix::new(
1357            3,
1358            4,
1359            [
1360                1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0,
1361            ]
1362            .to_vec(),
1363        );
1364        let test_matrix_b = Matrix::new(
1365            4,
1366            3,
1367            [
1368                12.0, 11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 1.0,
1369            ]
1370            .to_vec(),
1371        );
1372
1373        let matrix_product = test_matrix_a * test_matrix_b;
1374
1375        assert_eq!(matrix_product.entries.len(), 9);
1376
1377        assert_eq!(
1378            matrix_product.entries,
1379            [
1380                1.0 * 12.0 + 2.0 * 9.0 + 3.0 * 6.0 + 4.0 * 3.0,
1381                1.0 * 11.0 + 2.0 * 8.0 + 3.0 * 5.0 + 4.0 * 2.0,
1382                1.0 * 10.0 + 2.0 * 7.0 + 3.0 * 4.0 + 4.0 * 1.0,
1383                5.0 * 12.0 + 6.0 * 9.0 + 7.0 * 6.0 + 8.0 * 3.0,
1384                5.0 * 11.0 + 6.0 * 8.0 + 7.0 * 5.0 + 8.0 * 2.0,
1385                5.0 * 10.0 + 6.0 * 7.0 + 7.0 * 4.0 + 8.0 * 1.0,
1386                9.0 * 12.0 + 10.0 * 9.0 + 11.0 * 6.0 + 12.0 * 3.0,
1387                9.0 * 11.0 + 10.0 * 8.0 + 11.0 * 5.0 + 12.0 * 2.0,
1388                9.0 * 10.0 + 10.0 * 7.0 + 11.0 * 4.0 + 12.0 * 1.0
1389            ]
1390        );
1391    }
1392
1393    #[test]
1394    fn test_matrix_multiply_constant_value_initialiser() {
1395        let test_matrix_a =
1396            Matrix::new_constant_value(3, 4, 5.0).expect("Unable to create test_matrix_a");
1397        let test_matrix_b =
1398            Matrix::new_constant_value(4, 3, 4.0).expect("Unable to create test_matrix_b");
1399
1400        let matrix_product = test_matrix_a * test_matrix_b;
1401
1402        assert_eq!(matrix_product.entries.len(), 9);
1403
1404        assert_eq!(
1405            matrix_product.entries,
1406            [80.0, 80.0, 80.0, 80.0, 80.0, 80.0, 80.0, 80.0, 80.0,]
1407        );
1408    }
1409
1410    #[test]
1411    fn test_transpose() {
1412        let test_matrix = Matrix::new(2, 3, [1.0, 2.0, -1.0, 0.0, 3.0, 7.0].to_vec());
1413
1414        let transpose_matrix = test_matrix.transpose();
1415
1416        assert_eq!(transpose_matrix.m, test_matrix.n);
1417        assert_eq!(transpose_matrix.n, test_matrix.m);
1418        assert_eq!(transpose_matrix.entries, [1.0, 0.0, 2.0, 3.0, -1.0, 7.0]);
1419    }
1420
1421    #[test]
1422    fn test_new_all_default() {
1423        let mut rng = rand::thread_rng();
1424        let n = rng.gen_range(1..=100);
1425        let m = rng.gen_range(1..=100);
1426        let test_matrix_f64 = new_all_default::<f64>(m, n)
1427            .expect("test_new_all_default: unable to create test_matrix_f64");
1428
1429        for entry in test_matrix_f64.entries {
1430            assert_eq!(entry, f64::default());
1431        }
1432
1433        let test_matrix_f32 = new_all_default::<f32>(m, n)
1434            .expect("test_new_all_default: unable to create test_matrix_f64");
1435
1436        for entry in test_matrix_f32.entries {
1437            assert_eq!(entry, f32::default());
1438        }
1439    }
1440
1441    #[test]
1442    fn test_new_identity_matrix() {
1443        let mut rng = rand::thread_rng();
1444        let dimension = rng.gen_range(1..=100);
1445
1446        let test_matrix_f64 = new_identity_matrix::<f64>(dimension)
1447            .expect("test_new_identity_matrix: unable to create test_matrix_f64");
1448
1449        for index in 0..dimension {
1450            assert_eq!(*test_matrix_f64.get_entry_ij(index, index), 1.0);
1451        }
1452
1453        let test_matrix_f32 = new_identity_matrix::<f32>(dimension)
1454            .expect("test_new_identity_matrix: unable to create test_matrix_f64");
1455
1456        for index in 0..dimension {
1457            assert_eq!(*test_matrix_f32.get_entry_ij(index, index), 1.0);
1458        }
1459    }
1460
1461    #[test]
1462    fn test_get_entry() {
1463        let entries = [
1464            1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0,
1465            17.0, 18.0, 19.0, 20.0, 21.0, 22.0, 23.0, 24.0, 25.0,
1466        ]
1467        .to_vec();
1468        let test_matrix = Matrix::new(5, 5, entries);
1469        let mut index = 1.0;
1470
1471        for i in 0..5 {
1472            for j in 0..5 {
1473                assert_eq!(*test_matrix.get_entry_ij(i, j), index);
1474                index += 1.0;
1475            }
1476        }
1477    }
1478
1479    #[test]
1480    fn test_set_entry() {
1481        let mut test_matrix =
1482            new_all_default::<f64>(5, 5).expect("test_set_entry: unable to create test_matrix");
1483        let mut index = 1.0;
1484
1485        for i in 0..5 {
1486            for j in 0..5 {
1487                test_matrix.set_entry_ij(i, j, &index);
1488                assert_eq!(*test_matrix.get_entry_ij(i, j), index);
1489                index += 1.0;
1490            }
1491        }
1492    }
1493
1494    #[test]
1495    fn test_partition() {
1496        let test_matrix = Matrix::new(
1497            5,
1498            5,
1499            [
1500                1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0,
1501                16.0, 17.0, 18.0, 19.0, 20.0, 21.0, 22.0, 23.0, 24.0, 25.0,
1502            ]
1503            .to_vec(),
1504        );
1505
1506        let submatrices = test_matrix.partition(&[3, 2], &[2, 3]);
1507
1508        assert_eq!(
1509            submatrices,
1510            [
1511                Matrix::new(2, 3, [1.0, 2.0, 3.0, 6.0, 7.0, 8.0].to_vec()),
1512                Matrix::new(2, 2, [4.0, 5.0, 9.0, 10.0].to_vec()),
1513                Matrix::new(
1514                    3,
1515                    3,
1516                    [11.0, 12.0, 13.0, 16.0, 17.0, 18.0, 21.0, 22.0, 23.0].to_vec()
1517                ),
1518                Matrix::new(3, 2, [14.0, 15.0, 19.0, 20.0, 24.0, 25.0].to_vec())
1519            ]
1520        );
1521
1522        let test_matrix = Matrix::new(
1523            3,
1524            4,
1525            [
1526                0.0, 1.0, -4.0, -7.0, 0.0, 0.0, 3.0, -1.0, 0.0, 0.0, 3.0, -1.0,
1527            ]
1528            .to_vec(),
1529        );
1530
1531        let submatrices = test_matrix.partition(&[4], &[1, 2]);
1532
1533        assert_eq!(
1534            submatrices,
1535            [
1536                Matrix::new(1, 4, [0.0, 1.0, -4.0, -7.0].to_vec()),
1537                Matrix::new(2, 4, [0.0, 0.0, 3.0, -1.0, 0.0, 0.0, 3.0, -1.0,].to_vec()),
1538            ]
1539        );
1540    }
1541
1542    #[test]
1543    fn test_submatrix() {
1544        let test_matrix = Matrix::new(
1545            5,
1546            5,
1547            [
1548                1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0,
1549                16.0, 17.0, 18.0, 19.0, 20.0, 21.0, 22.0, 23.0, 24.0, 25.0,
1550            ]
1551            .to_vec(),
1552        );
1553
1554        let submatrix = test_matrix.submatrix(&[1, 3], &[0, 2, 4]);
1555
1556        assert_eq!(submatrix.n, 2);
1557        assert_eq!(submatrix.m, 3);
1558        assert_eq!(submatrix.entries, [2.0, 4.0, 12.0, 14.0, 22.0, 24.0]);
1559    }
1560
1561    #[test]
1562    fn test_row_interchange() {
1563        let test_matrix = Matrix::new(2, 3, [1.0, 2.0, 3.0, 4.0, 5.0, 6.0].to_vec());
1564
1565        let result = test_matrix.row_interchange(0, 1);
1566
1567        assert_eq!(
1568            result,
1569            Matrix::new(2, 3, [4.0, 5.0, 6.0, 1.0, 2.0, 3.0].to_vec()),
1570        );
1571    }
1572
1573    #[test]
1574    fn test_multiply_row_by_scalar() {
1575        let test_matrix = Matrix::new(3, 3, [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0].to_vec());
1576
1577        let result = test_matrix.multiply_row_by_scalar(2, 2.0);
1578
1579        assert_eq!(
1580            result,
1581            Matrix::new(
1582                3,
1583                3,
1584                [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 14.0, 16.0, 18.0].to_vec()
1585            ),
1586        );
1587    }
1588
1589    #[test]
1590    fn test_add_row_to_scalar_multiple_of_row() {
1591        let test_matrix = Matrix::new(3, 3, [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0].to_vec());
1592
1593        let result = test_matrix.add_row_to_scalar_multiple_of_row(0, 2, 5.0);
1594
1595        assert_eq!(
1596            result,
1597            Matrix::new(
1598                3,
1599                3,
1600                [36.0, 42.0, 48.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0].to_vec()
1601            ),
1602        );
1603    }
1604
1605    #[test]
1606    fn test_all_zeroes() {
1607        let mut entries = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0].to_vec();
1608        let test_matrix = Matrix::new(3, 4, entries.clone());
1609
1610        assert_eq!(true, test_matrix.all_zeroes());
1611
1612        let mut rng = rand::thread_rng();
1613        let random_index = rng.gen_range(0..12);
1614        entries[random_index] = 1.0;
1615
1616        let test_matrix = Matrix::new(3, 4, entries.clone());
1617
1618        assert_eq!(false, test_matrix.all_zeroes());
1619    }
1620
1621    #[test]
1622    fn test_first_non_zero_vec_index() {
1623        let mut rng = rand::thread_rng();
1624        let random_index = rng.gen_range(0..12);
1625        let mut test_vec: Vec<Vec<f64>> = vec![];
1626        let zero_vec = [0.0, 0.0, 0.0].to_vec();
1627        let non_zero_vec = [1.0, 2.0, 3.0].to_vec();
1628
1629        for i in 0..12 {
1630            if i == random_index {
1631                test_vec.push(non_zero_vec.clone());
1632            } else {
1633                test_vec.push(zero_vec.clone());
1634            }
1635        }
1636
1637        assert_eq!(random_index, first_non_zero_vec_index(&test_vec));
1638    }
1639
1640    #[test]
1641    fn test_index_of_first_non_zero_element_in_vec() {
1642        let mut rng = rand::thread_rng();
1643        let random_index = rng.gen_range(0..12);
1644        let mut test_vec: Vec<f64> = vec![];
1645
1646        for i in 0..12 {
1647            if i == random_index {
1648                test_vec.push(1.0);
1649            } else {
1650                test_vec.push(0.0);
1651            }
1652        }
1653
1654        assert_eq!(
1655            random_index,
1656            index_of_first_non_zero_element_in_vec(&test_vec)
1657        );
1658    }
1659
1660    #[test]
1661    fn test_row_echolon_form() {
1662        let test_matrix = Matrix::new(
1663            3,
1664            4,
1665            [
1666                0.0, 0.0, 3.0, -1.0, 0.0, -1.0, 4.0, 7.0, 0.0, -1.0, 7.0, 6.0,
1667            ]
1668            .to_vec(),
1669        );
1670
1671        let result = test_matrix.row_echolon_form();
1672
1673        assert_eq!(
1674            result,
1675            Matrix::new(
1676                3,
1677                4,
1678                [
1679                    0.0,
1680                    1.0,
1681                    0.0,
1682                    -25.0 / 3.0,
1683                    0.0,
1684                    0.0,
1685                    1.0,
1686                    -1.0 / 3.0,
1687                    0.0,
1688                    0.0,
1689                    0.0,
1690                    0.0
1691                ]
1692                .to_vec()
1693            ),
1694        );
1695
1696        let test_matrix = Matrix::new(
1697            3,
1698            3,
1699            [4.0, -8.0, 16.0, 1.0, -3.0, 6.0, 2.0, 1.0, 1.0].to_vec(),
1700        );
1701
1702        let result = test_matrix.row_echolon_form();
1703
1704        assert_eq!(
1705            result,
1706            new_identity_matrix(3)
1707                .expect("test_row_echolon_form: unable to create new_identity_matrix")
1708        );
1709
1710        let test_matrix = Matrix::new(
1711            3,
1712            4,
1713            [0.0, 2.0, 1.0, 4.0, 0.0, 0.0, 2.0, 6.0, 1.0, 0.0, -3.0, 2.0].to_vec(),
1714        );
1715
1716        let result = test_matrix.row_echolon_form();
1717
1718        assert_eq!(
1719            result,
1720            Matrix::new(
1721                3,
1722                4,
1723                [1.0, 0.0, 0.0, 11.0, 0.0, 1.0, 0.0, 0.5, 0.0, 0.0, 1.0, 3.0].to_vec()
1724            )
1725        );
1726
1727        let test_matrix = Matrix::new(
1728            3,
1729            3,
1730            [
1731                ComplexNumber::new(2.0, 0.0),
1732                ComplexNumber::new(8.0, 2.0),
1733                ComplexNumber::new(6.0, -6.0),
1734                ComplexNumber::new(0.0, 1.0),
1735                ComplexNumber::new(0.0, 5.0),
1736                ComplexNumber::new(3.0, 3.0),
1737                ComplexNumber::new(1.0, 2.0),
1738                ComplexNumber::new(4.0, 11.0),
1739                ComplexNumber::new(9.0, 3.0),
1740            ]
1741            .to_vec(),
1742        );
1743
1744        let result = test_matrix.row_echolon_form();
1745
1746        assert_eq!(
1747            result,
1748            Matrix::new(
1749                3,
1750                3,
1751                [
1752                    ComplexNumber::new(1.0, 0.0),
1753                    ComplexNumber::new(0.0, 0.0),
1754                    ComplexNumber::new(3.0, -3.0),
1755                    ComplexNumber::new(0.0, 0.0),
1756                    ComplexNumber::new(1.0, 0.0),
1757                    ComplexNumber::new(0.0, 0.0),
1758                    ComplexNumber::new(0.0, 0.0),
1759                    ComplexNumber::new(0.0, 0.0),
1760                    ComplexNumber::new(0.0, 0.0),
1761                ]
1762                .to_vec()
1763            )
1764        );
1765
1766        let test_matrix = Matrix::new(
1767            3,
1768            2,
1769            [
1770                ComplexNumber::new(1.0, 1.0),
1771                ComplexNumber::new(1.0, -1.0),
1772                ComplexNumber::new(2.0, 0.0),
1773                ComplexNumber::new(2.0, 0.0),
1774                ComplexNumber::new(1.0, 2.0),
1775                ComplexNumber::new(2.0, -1.0),
1776            ]
1777            .to_vec(),
1778        );
1779
1780        let result = test_matrix.row_echolon_form();
1781
1782        assert_eq!(
1783            result,
1784            Matrix::new(
1785                3,
1786                2,
1787                [
1788                    ComplexNumber::new(1.0, 0.0),
1789                    ComplexNumber::new(0.0, 0.0),
1790                    ComplexNumber::new(0.0, 0.0),
1791                    ComplexNumber::new(1.0, 0.0),
1792                    ComplexNumber::new(0.0, 0.0),
1793                    ComplexNumber::new(0.0, 0.0),
1794                ]
1795                .to_vec()
1796            )
1797        );
1798    }
1799
1800    #[test]
1801    fn test_column_echolon_form() {
1802        let test_matrix = Matrix::new(2, 3, [1.0, 2.0, 3.0, 2.0, 3.0, 4.0].to_vec());
1803
1804        let result = test_matrix.column_echolon_form();
1805
1806        assert_eq!(
1807            result,
1808            Matrix::new(2, 3, [1.0, 0.0, 0.0, 0.0, 1.0, 0.0,].to_vec())
1809        );
1810
1811        let test_matrix = Matrix::new(3, 3, [1.0, 1.0, 1.0, 2.0, 2.0, 2.0, 3.0, 3.0, 3.0].to_vec());
1812
1813        let result = test_matrix.column_echolon_form();
1814
1815        assert_eq!(
1816            result,
1817            Matrix::new(3, 3, [1.0, 0.0, 0.0, 2.0, 0.0, 0.0, 3.0, 0.0, 0.0].to_vec())
1818        );
1819    }
1820
1821    #[test]
1822    fn test_determinant() {
1823        let test_matrix = Matrix::new(
1824            4,
1825            4,
1826            [
1827                2.0, 2.0, -1.0, 9.0, 4.0, 2.0, 1.0, 17.0, 1.0, -1.0, 1.0, 2.0, 2.0, 1.0, 1.0, 9.0,
1828            ]
1829            .to_vec(),
1830        );
1831        let result = test_matrix.determinant();
1832        assert_eq!(result.expect("Unable to calculate determinant"), 0.0);
1833
1834        let test_matrix = Matrix::new(
1835            3,
1836            3,
1837            [
1838                ComplexNumber::new(1.0, 0.0),
1839                ComplexNumber::new(1.0, 0.0),
1840                ComplexNumber::new(0.0, 1.0),
1841                ComplexNumber::new(1.0, 1.0),
1842                ComplexNumber::new(1.0, 1.0),
1843                ComplexNumber::new(1.0, 0.0),
1844                ComplexNumber::new(2.0, 3.0),
1845                ComplexNumber::new(0.0, -1.0),
1846                ComplexNumber::new(3.0, 0.0),
1847            ]
1848            .to_vec(),
1849        );
1850
1851        let determinant = test_matrix.determinant();
1852
1853        fn delta_real(determinant: ComplexNumber<f64>, expectation: f64) -> bool {
1854            determinant.real - expectation < (1.0 / 1000000000000000000000000000000.0)
1855        }
1856
1857        fn delta_complex(determinant: ComplexNumber<f64>, expectation: f64) -> bool {
1858            determinant.complex - expectation < (1.0 / 1000000000000000000000000000000.0)
1859        }
1860        assert!(delta_real(
1861            determinant.expect("Unable to calculate determinant"),
1862            8.0
1863        ));
1864        assert!(delta_complex(
1865            determinant.expect("Unable to calculate determinant"),
1866            6.0
1867        ));
1868    }
1869
1870    #[test]
1871    fn test_determinant_err() {
1872        let test_matrix = new_all_default::<f64>(3, 2)
1873            .expect("test_determinant_err: unable to create tesdt_matrix");
1874
1875        let result = test_matrix.determinant();
1876
1877        assert_eq!(result, Err("Non square matrices do not have determinants!"));
1878    }
1879
1880    #[test]
1881    fn test_adjoint() {
1882        let test_matrix = Matrix::new(3, 3, [3.0, 4.0, 3.0, 5.0, 7.0, 2.0, 0.0, 0.0, 1.0].to_vec());
1883
1884        assert_eq!(
1885            test_matrix
1886                .adjoint()
1887                .expect("test_adjoint: unable to calculate matrix adjoint"),
1888            Matrix::new(
1889                3,
1890                3,
1891                [
1892                    7.0,
1893                    -4.0,
1894                    -13.0,
1895                    -5.0,
1896                    3.0,
1897                    9.0,
1898                    0.0,
1899                    0.0,
1900                    1.0000000000000018
1901                ]
1902                .to_vec()
1903            )
1904        );
1905    }
1906
1907    #[test]
1908    fn test_inverse() {
1909        let test_matrix = Matrix::new(
1910            3,
1911            3,
1912            [2.0, 3.0, 0.0, 0.0, 3.0, -3.0, -2.0, 3.0, 3.0].to_vec(),
1913        );
1914
1915        let inverse = test_matrix
1916            .inverse()
1917            .expect("Unable to compute matrix inverse in test_inverse");
1918
1919        assert_eq!(inverse.m, 3);
1920        assert_eq!(inverse.n, 3);
1921
1922        let expected = [
1923            1.0 / 3.0,
1924            -1.0 / 6.0,
1925            -1.0 / 6.0,
1926            1.0 / 9.0,
1927            1.0 / 9.0,
1928            1.0 / 9.0,
1929            1.0 / 9.0,
1930            -2.0 / 9.0,
1931            1.0 / 9.0,
1932        ]
1933        .to_vec();
1934
1935        for i in 0..expected.len() {
1936            assert!(expected[i] - inverse.entries[i] < 1.0 / 1000000000000000.0);
1937        }
1938
1939        let test_matrix = Matrix::new(
1940            5,
1941            5,
1942            [
1943                1.0, 2.0, 1.0, 2.0, 3.0, 2.0, 3.0, 1.0, 0.0, 1.0, 2.0, 2.0, 1.0, 0.0, 0.0, 1.0,
1944                1.0, 1.0, 1.0, 1.0, 0.0, -2.0, 0.0, 2.0, -2.0,
1945            ]
1946            .to_vec(),
1947        );
1948
1949        assert_eq!(
1950            test_matrix
1951                .inverse()
1952                .expect("Error computing inverse of test case 2 in test_inverse"),
1953            Matrix::new(
1954                5,
1955                5,
1956                [
1957                    2.0, -4.0, 5.0, -3.0, -0.5, -1.0, 3.0, -3.0, 1.0, 0.5, -2.0, 2.0, -3.0, 4.0,
1958                    0.0, 0.0, 1.0, -1.0, 0.0, 0.5, 1.0, -2.0, 2.0, -1.0, -0.5
1959                ]
1960                .to_vec()
1961            )
1962        );
1963    }
1964
1965    #[test]
1966    fn test_inverse_err_case() {
1967        let mut rng = rand::thread_rng();
1968        let random_m = rng.gen_range(0..100);
1969        let mut random_n: usize = random_m.clone();
1970
1971        while random_n == random_m {
1972            random_n = rng.gen_range(0..100);
1973        }
1974
1975        let test_matrix = new_all_default::<f64>(random_m, random_n);
1976
1977        assert_eq!(
1978            test_matrix
1979                .expect("Unable to construct new_all_default matrix in test_inverse_err_case")
1980                .inverse(),
1981            Err("Unable to compute inverse when n is not equal to m")
1982        );
1983    }
1984}
rust_matrix/matrix_algebra.rs

rust_matrix/
matrix_algebra.rs
