aoclib_rs/
matrix.rs

1use std::{
2    cmp::Ordering,
3    fmt,
4    fmt::Formatter,
5    ops::{Deref, DerefMut, Index, IndexMut, Mul, MulAssign},
6};
7
8use crate::fold_while;
9
10use num_rational::Rational64 as R64;
11
12/// A row vector (in the linear algebra sense) of rational numbers.
13#[derive(Clone, Debug, PartialEq)]
14pub struct RowVec(Vec<R64>);
15
16impl RowVec {
17    // TODO: test initializers
18    pub fn new(v: Vec<R64>) -> Self {
19        Self(v)
20    }
21
22    pub fn zeros(len: usize) -> Self {
23        Self(vec![R64::ZERO; len])
24    }
25
26    pub fn from_int_vec(v: Vec<i64>) -> Self {
27        Self(v.iter().map(|i| R64::from_integer(*i)).collect())
28    }
29
30    pub fn empty() -> Self {
31        Self(Vec::new())
32    }
33
34    /// Implements a "+=" operation (without actually defining the operator) which returns an error
35    /// if the 2 `RowVec`s are of different sizes.
36    pub fn add_assign(&mut self, rhs: &Self) -> anyhow::Result<()> {
37        if self.0.len() != rhs.0.len() {
38            anyhow::bail!(
39                "addition of RowVecs of different sizes: {} vs {}",
40                self.0.len(),
41                rhs.0.len()
42            );
43        }
44
45        self.0
46            .iter_mut()
47            .enumerate()
48            .for_each(|(i, e)| *e += rhs.0[i]);
49
50        Ok(())
51    }
52
53    /// Implements a "+" operation (without actually defining the operator) which returns an error
54    /// if the 2 `RowVec`s are of different sizes.
55    pub fn add(&self, rhs: &Self) -> anyhow::Result<Self> {
56        let mut out = self.clone();
57        out.add_assign(rhs)?;
58        Ok(out)
59    }
60
61    /// Divides the `RowVec` by the `leader` (ie, the first non-zero entry), to make the leader `1`.
62    /// Useful in matrix row reduction.
63    pub fn normalize(&mut self) {
64        let Some(leader_col) = self.leader_col() else {
65            return;
66        };
67
68        if self.0[leader_col] != R64::ONE {
69            let factor = self.0[leader_col].recip();
70            *self *= factor;
71        }
72    }
73
74    pub fn is_zeros(&self) -> bool {
75        fold_while(self.0.iter(), true, |_, v| {
76            let r = *v == R64::ZERO;
77            (r, r)
78        })
79    }
80
81    /// Returns the column (or index) of the first non-zero entry, or `None` if the `RowVec` is empty
82    /// or all zero.
83    pub fn leader_col(&self) -> Option<usize> {
84        self.0.iter().position(|&e| e != R64::ZERO)
85    }
86}
87
88impl Deref for RowVec {
89    type Target = Vec<R64>;
90
91    fn deref(&self) -> &Self::Target {
92        &self.0
93    }
94}
95
96impl DerefMut for RowVec {
97    fn deref_mut(&mut self) -> &mut Self::Target {
98        &mut self.0
99    }
100}
101
102impl MulAssign<R64> for RowVec {
103    fn mul_assign(&mut self, rhs: R64) {
104        self.0.iter_mut().for_each(|cell| *cell *= rhs);
105    }
106}
107
108impl Mul<R64> for RowVec {
109    type Output = RowVec;
110
111    fn mul(mut self, rhs: R64) -> Self::Output {
112        self *= rhs;
113        self
114    }
115}
116
117impl fmt::Display for RowVec {
118    fn fmt(&self, f: &mut Formatter) -> fmt::Result {
119        write!(f, "[ ")?;
120        for v in &self.0 {
121            write!(f, "{} ", *v)?;
122        }
123        write!(f, "]")?;
124        Ok(())
125    }
126}
127
128/// A matrix (in the linear algebra sense) of rational numbers.
129#[derive(Clone, Debug, PartialEq)]
130pub struct Matrix(Vec<RowVec>);
131
132impl Matrix {
133    /// The caller is responsible for making sure the input is valid, meaning that no rows are
134    /// empty, and each row is the same size. Failing to do so could lead to undefined behaviour or
135    /// panics later on.
136    pub fn new(m: Vec<Vec<R64>>) -> Self {
137        Self(m.into_iter().map(RowVec).collect())
138    }
139
140    /// The caller is responsible for making sure the input is valid, meaning that no rows are
141    /// empty, and each row is the same size. Failing to do so could lead to undefined behaviour or
142    /// panics later on.
143    pub fn from_row_vecs(m: Vec<RowVec>) -> Self {
144        Self(m)
145    }
146
147    // TODO: test
148    /// The caller is responsible for making sure the input is valid, meaning that no rows are
149    /// empty, and each row is the same size. Failing to do so could lead to undefined behaviour or
150    /// panics later on.
151    pub fn from_int_vecs(m: Vec<Vec<i64>>) -> Self {
152        Self(m.into_iter().map(RowVec::from_int_vec).collect())
153    }
154
155    pub fn zeros(rows: usize, cols: usize) -> Self {
156        if rows == 0 || cols == 0 {
157            Self::empty()
158        } else {
159            Self(vec![RowVec::zeros(cols); rows])
160        }
161    }
162
163    pub fn empty() -> Self {
164        Self(Vec::new())
165    }
166
167    /// Puts the matrix into Reduced Row Echelon Form.
168    pub fn rref(&mut self) {
169        self.r#ref();
170        for row in 0..self.0.len() {
171            self.eliminate_above_leader(row);
172        }
173        self.leader_sort();
174        self.normalize();
175    }
176
177    // TODO(?) technically does more than necessary without the sort after each eliminate?
178    /// Puts the matrix into (unreduced) Row Echelon Form.
179    pub fn r#ref(&mut self) {
180        self.leader_sort();
181        for row in 0..self.0.len() {
182            self.eliminate_below_leader(row);
183        }
184        self.leader_sort();
185    }
186
187    /// Sort the rows of the matrix such the ones with the leftmost leaders (ie, leftmost non-zero
188    /// entries) come first. All-zero rows are at the bottom.
189    ///
190    /// ```
191    /// use aoclib_rs::matrix::Matrix;
192    ///
193    /// let mut m = Matrix::from_int_vecs(vec![
194    ///     vec![0, 2, 3, 4],
195    ///     vec![2, 3, 6, 3],
196    ///     vec![0, 0, 4, 5],
197    ///     vec![0, 0, 0, 3],
198    ///     vec![0, 5, 6, 3],
199    ///     vec![0, 0, 0, 0],
200    ///     vec![3, 4, 2, 6],
201    /// ]);
202    /// m.leader_sort();
203    ///
204    /// assert_eq!(
205    ///     m,
206    ///     Matrix::from_int_vecs(vec![
207    ///         vec![2, 3, 6, 3],
208    ///         vec![3, 4, 2, 6],
209    ///         vec![0, 2, 3, 4],
210    ///         vec![0, 5, 6, 3],
211    ///         vec![0, 0, 4, 5],
212    ///         vec![0, 0, 0, 3],
213    ///         vec![0, 0, 0, 0],
214    ///     ])
215    /// );
216    /// ```
217    pub fn leader_sort(&mut self) {
218        self.0.sort_by(|a, b| {
219            for (i, ai) in a.iter().enumerate() {
220                let bi = b[i];
221                if *ai == R64::ZERO && bi != R64::ZERO {
222                    return Ordering::Greater;
223                } else if bi == R64::ZERO && *ai != R64::ZERO {
224                    return Ordering::Less;
225                }
226            }
227            Ordering::Equal
228        });
229    }
230
231    /// Perform a step in row reduction by eliminating the entries in the same column as `row`'s
232    /// leader for all rows below `row`.
233    pub fn eliminate_below_leader(&mut self, row: usize) {
234        let Some(leader_col) = self.0[row].leader_col() else {
235            return;
236        };
237
238        for i in (row + 1)..self.0.len() {
239            self.eliminate(row, i, leader_col);
240        }
241    }
242
243    /// Perform a step in row reduction by eliminating the `other_row` entry corresponding to
244    /// `leader_col` (which is the column of the leader in `selected_row`).
245    ///
246    /// ```
247    /// use aoclib_rs::matrix::Matrix;
248    ///
249    /// let mut m = Matrix::from_int_vecs(vec![vec![0, 2, 3, 4], vec![2, 4, 6, 3]]);
250    /// m.eliminate(0, 1, 1);
251    ///
252    /// assert_eq!(
253    ///     m,
254    ///     Matrix::from_int_vecs(vec![vec![0, 2, 3, 4], vec![2, 0, 0, -5]])
255    /// );
256    /// ```
257    pub fn eliminate(&mut self, selected_row: usize, other_row: usize, leader_col: usize) {
258        if self.0[other_row][leader_col] == R64::ZERO {
259            return;
260        }
261
262        let factor = -self.0[other_row][leader_col] / self.0[selected_row][leader_col];
263        let term = self.0[selected_row].clone() * factor;
264        self.0[other_row].add_assign(&term).unwrap();
265    }
266
267    /// Divides rows to ensure that each leader is `1`.
268    pub fn normalize(&mut self) {
269        self.0.iter_mut().for_each(|row| row.normalize());
270    }
271
272    /// Perform a step in row reduction by eliminating the entries in the same column as `row`'s
273    /// leader for all rows above `row`.
274    pub fn eliminate_above_leader(&mut self, row: usize) {
275        let Some(leader_col) = self.0[row].leader_col() else {
276            return;
277        };
278
279        for i in 0..row {
280            self.eliminate(row, i, leader_col);
281        }
282    }
283
284    /// Implements a "+=" operation (without actually defining the operator) which returns an error
285    /// if the 2 `Matrix`es are of different sizes.
286    pub fn add_assign(&mut self, rhs: &Self) -> anyhow::Result<()> {
287        if self.0.len() != rhs.0.len() {
288            anyhow::bail!(
289                "addition of matrices of different heights: {} vs {}",
290                self.0.len(),
291                rhs.0.len()
292            );
293        }
294
295        if self.0.is_empty() {
296            return Ok(());
297        }
298
299        if self.0[0].len() != rhs.0[0].len() {
300            anyhow::bail!(
301                "addition of matrices of different widths: {} vs {}",
302                self.0[0].len(),
303                rhs.0[0].len()
304            );
305        }
306
307        if self.0[0].is_empty() {
308            return Ok(());
309        }
310
311        for i in 0..self.0.len() {
312            self.0[i].add_assign(&rhs.0[i])?;
313        }
314
315        Ok(())
316    }
317
318    /// Implements a "+" operation (without actually defining the operator) which returns an error
319    /// if the 2 `Matrix`es are of different sizes.
320    pub fn add(&self, rhs: &Self) -> anyhow::Result<Self> {
321        let mut out = self.clone();
322        out.add_assign(rhs)?;
323        Ok(out)
324    }
325
326    // TODO: test
327    pub fn is_empty(&self) -> bool {
328        self.0.is_empty()
329    }
330
331    /// Implements matrix multiplication. Returns errors if the two matrices are incorrect sizes
332    /// for multiplication.
333    pub fn matrix_mul(&self, rhs: &Self) -> anyhow::Result<Self> {
334        if self.is_empty() && rhs.is_empty() {
335            return Ok(Self(Vec::new()));
336        }
337
338        if self.is_empty() {
339            anyhow::bail!("multiplication of empty matrix with non-empty matrix");
340        }
341
342        if self.0[0].len() != rhs.0.len() {
343            anyhow::bail!(
344                "multiplication of incompatible matrices: lhs width {} vs rhs height {}",
345                self.0[0].len(),
346                rhs.0.len()
347            );
348        }
349
350        let mut new = Self::zeros(self.0.len(), rhs.0[0].len());
351        for i in 0..rhs.0[0].len() {
352            for j in 0..self.0.len() {
353                for k in 0..self.0[0].len() {
354                    new.0[j][i] += self.0[j][k] * rhs.0[k][i];
355                }
356            }
357        }
358        Ok(new)
359    }
360
361    /// Appends a `RowVec` to the bottom of the matrix. Note that the caller is responsible for
362    /// ensuring the RowVec is the correct size. Appending `RowVec`s of the wrong size could lead to
363    /// undefined behaviour or panics later on.
364    pub fn append_row(&mut self, r: RowVec) {
365        self.0.push(r);
366    }
367
368    /// Remove row `i` from the `Matrix`.
369    pub fn remove_row(&mut self, i: usize) {
370        self.0.remove(i);
371    }
372
373    pub fn height(&self) -> usize {
374        self.0.len()
375    }
376
377    pub fn width(&self) -> usize {
378        if self.0.is_empty() {
379            0
380        } else {
381            self.0[0].len()
382        }
383    }
384
385    /// Returns an iterator over rows of the `Matrix`.
386    pub fn iter(&self) -> impl Iterator<Item = &RowVec> {
387        self.0.iter()
388    }
389
390    pub fn is_zeros(&self) -> bool {
391        fold_while(self.0.iter(), true, |_, v| {
392            let r = v.is_zeros();
393            (r, r)
394        })
395    }
396
397    /// Returns a new matrix with a single column whose contents are a copy of a column `c` from
398    /// the original matrix.
399    pub fn get_column_copy(&self, c: usize) -> Self {
400        Self::new(self.iter().map(|r| vec![r[c]]).collect())
401    }
402}
403
404impl MulAssign<R64> for Matrix {
405    fn mul_assign(&mut self, rhs: R64) {
406        self.0.iter_mut().for_each(|row| *row *= rhs);
407    }
408}
409
410impl Mul<R64> for Matrix {
411    type Output = Self;
412
413    fn mul(mut self, rhs: R64) -> Self::Output {
414        self *= rhs;
415        self
416    }
417}
418
419impl Index<usize> for Matrix {
420    type Output = RowVec;
421
422    fn index(&self, index: usize) -> &Self::Output {
423        &self.0[index]
424    }
425}
426
427impl IndexMut<usize> for Matrix {
428    fn index_mut(&mut self, index: usize) -> &mut Self::Output {
429        &mut self.0[index]
430    }
431}
432
433impl Index<(usize, usize)> for Matrix {
434    type Output = R64;
435
436    fn index(&self, index: (usize, usize)) -> &Self::Output {
437        &self.0[index.0][index.1]
438    }
439}
440
441impl IndexMut<(usize, usize)> for Matrix {
442    fn index_mut(&mut self, index: (usize, usize)) -> &mut Self::Output {
443        &mut self.0[index.0][index.1]
444    }
445}
446
447impl IntoIterator for Matrix {
448    type Item = RowVec;
449    type IntoIter = std::vec::IntoIter<Self::Item>;
450
451    fn into_iter(self) -> Self::IntoIter {
452        self.0.into_iter()
453    }
454}
455
456impl fmt::Display for Matrix {
457    fn fmt(&self, f: &mut Formatter) -> fmt::Result {
458        write!(f, "[ ")?;
459        let mut first = true;
460        for v in &self.0 {
461            if !first {
462                write!(f, "\n  ")?;
463            }
464            first = false;
465            write!(f, "{} ", *v)?;
466        }
467        write!(f, "]")?;
468        Ok(())
469    }
470}
471
472#[cfg(test)]
473mod tests {
474    use super::*;
475
476    #[test]
477    fn test_row_vec_add_assign_good() {
478        let mut rv1 = RowVec::from_int_vec(vec![1, 2, 3]);
479        let rv2 = RowVec::from_int_vec(vec![4, 7, 9]);
480        let r = rv1.add_assign(&rv2);
481        assert!(r.is_ok());
482        assert_eq!(rv1, RowVec::from_int_vec(vec![5, 9, 12]));
483    }
484
485    #[test]
486    fn test_row_vec_add_assign_bad() {
487        let mut rv1 = RowVec::from_int_vec(vec![1, 2, 3]);
488        let rv2 = RowVec::from_int_vec(vec![4, 7, 9, 10]);
489        let r = rv1.add_assign(&rv2);
490        assert!(r.is_err());
491        assert_eq!(rv1, RowVec::from_int_vec(vec![1, 2, 3]));
492    }
493
494    #[test]
495    fn test_row_vec_add_assign_empty() {
496        let mut rv1 = RowVec::empty();
497        let rv2 = RowVec::empty();
498        let r = rv1.add_assign(&rv2);
499        assert!(r.is_ok());
500        assert_eq!(rv1, RowVec::empty());
501    }
502
503    #[test]
504    fn test_row_vec_add_good() {
505        let rv1 = RowVec::from_int_vec(vec![1, 2, 3]);
506        let rv2 = RowVec::from_int_vec(vec![4, 7, 9]);
507        let r = rv1.add(&rv2);
508        assert!(r.is_ok());
509        assert_eq!(r.unwrap(), RowVec::from_int_vec(vec![5, 9, 12]));
510    }
511
512    #[test]
513    fn test_row_vec_add_bad() {
514        let rv1 = RowVec::from_int_vec(vec![1, 2, 3]);
515        let rv2 = RowVec::from_int_vec(vec![4, 7, 9, 10]);
516        let r = rv1.add(&rv2);
517        assert!(r.is_err());
518    }
519
520    #[test]
521    fn test_row_vec_add_empty() {
522        let rv1 = RowVec::empty();
523        let rv2 = RowVec::empty();
524        let r = rv1.add(&rv2);
525        assert!(r.is_ok());
526        assert_eq!(r.unwrap(), RowVec(vec![]));
527    }
528
529    #[test]
530    fn test_row_vec_normalize_good() {
531        let mut rv = RowVec::from_int_vec(vec![7, 5, 9]);
532        rv.normalize();
533        assert_eq!(rv, RowVec(vec![R64::ONE, R64::new(5, 7), R64::new(9, 7)]));
534    }
535
536    #[test]
537    fn test_row_vec_normalize_leading_zeros() {
538        let mut rv = RowVec::from_int_vec(vec![0, 0, 7, 5, 9]);
539        rv.normalize();
540        assert_eq!(
541            rv,
542            RowVec(vec![
543                R64::ZERO,
544                R64::ZERO,
545                R64::ONE,
546                R64::new(5, 7),
547                R64::new(9, 7)
548            ])
549        );
550    }
551
552    #[test]
553    fn test_row_vec_normalize_leader_is_one() {
554        let mut rv = RowVec::from_int_vec(vec![1, 2, 3]);
555        rv.normalize();
556        assert_eq!(rv, RowVec::from_int_vec(vec![1, 2, 3]));
557    }
558
559    #[test]
560    fn test_row_vec_normalize_all_zeros() {
561        let mut rv = RowVec::from_int_vec(vec![0, 0, 0]);
562        rv.normalize();
563        assert_eq!(rv, RowVec::from_int_vec(vec![0, 0, 0]));
564    }
565
566    #[test]
567    fn test_row_vec_mul_assign_good() {
568        let mut rv = RowVec::from_int_vec(vec![1, 2, 3]);
569        rv *= R64::from_integer(4);
570        assert_eq!(rv, RowVec::from_int_vec(vec![4, 8, 12]));
571    }
572
573    #[test]
574    fn test_row_vec_mul_assign_empty() {
575        let mut rv = RowVec::empty();
576        rv *= R64::from_integer(4);
577        assert_eq!(rv, RowVec::empty());
578    }
579
580    #[test]
581    fn test_row_vec_mul_good() {
582        let rv = RowVec::from_int_vec(vec![1, 2, 3]);
583        assert_eq!(
584            rv * R64::from_integer(4),
585            RowVec::from_int_vec(vec![4, 8, 12])
586        );
587    }
588
589    #[test]
590    fn test_row_vec_mul_empty() {
591        let rv = RowVec::empty();
592        assert_eq!(rv * R64::from_integer(4), RowVec::from_int_vec(vec![]));
593    }
594
595    #[test]
596    fn test_row_vec_is_zeros_false() {
597        let rv = RowVec::from_int_vec(vec![0, 2, 0, 3, 0]);
598        assert!(!rv.is_zeros());
599    }
600
601    #[test]
602    fn test_row_vec_is_zeros_true() {
603        let rv = RowVec::from_int_vec(vec![0, 0, 0, 0, 0]);
604        assert!(rv.is_zeros());
605    }
606
607    #[test]
608    fn test_row_vec_is_zeros_empty() {
609        let rv = RowVec::empty();
610        assert!(rv.is_zeros());
611    }
612
613    #[test]
614    fn test_row_vec_leader_col_good() {
615        let rv = RowVec::from_int_vec(vec![0, 0, 1, 2, 3]);
616        assert_eq!(rv.leader_col(), Some(2));
617    }
618
619    #[test]
620    fn test_row_vec_leader_col_first() {
621        let rv = RowVec::from_int_vec(vec![4, 0, 1, 2, 3]);
622        assert_eq!(rv.leader_col(), Some(0));
623    }
624
625    #[test]
626    fn test_row_vec_leader_col_last() {
627        let rv = RowVec::from_int_vec(vec![0, 0, 0, 0, 3]);
628        assert_eq!(rv.leader_col(), Some(4));
629    }
630
631    #[test]
632    fn test_row_vec_leader_col_zeros() {
633        let rv = RowVec::from_int_vec(vec![0, 0, 0, 0, 0]);
634        assert_eq!(rv.leader_col(), None);
635    }
636
637    #[test]
638    fn test_row_vec_leader_col_empty() {
639        let rv = RowVec::empty();
640        assert_eq!(rv.leader_col(), None);
641    }
642
643    #[test]
644    fn test_matrix_zeros() {
645        let m = Matrix::zeros(4, 5);
646        assert_eq!(
647            m,
648            Matrix::from_int_vecs(vec![
649                vec![0, 0, 0, 0, 0],
650                vec![0, 0, 0, 0, 0],
651                vec![0, 0, 0, 0, 0],
652                vec![0, 0, 0, 0, 0]
653            ])
654        )
655    }
656
657    #[test]
658    fn test_matrix_zeros_no_rows() {
659        let m = Matrix::zeros(0, 5);
660        assert_eq!(m, Matrix::new(Vec::new()));
661    }
662
663    #[test]
664    fn test_matrix_zeros_no_cols() {
665        let m = Matrix::zeros(3, 0);
666        assert_eq!(m, Matrix::new(Vec::new()));
667    }
668
669    #[test]
670    fn test_matrix_zeros_empty() {
671        let m = Matrix::zeros(0, 0);
672        assert_eq!(m, Matrix::new(Vec::new()));
673    }
674
675    #[test]
676    fn test_matrix_leader_sort_good() {
677        let mut m = Matrix::from_int_vecs(vec![
678            vec![0, 0, 1, 10, 0],
679            vec![5, 0, 0, 0, 0],
680            vec![0, 0, 0, 0, 0],
681            vec![0, -1, 0, 0, 0],
682            vec![0, 1, 0, 0, 0],
683        ]);
684        m.leader_sort();
685        assert_eq!(
686            m,
687            Matrix::from_int_vecs(vec![
688                vec![5, 0, 0, 0, 0],
689                vec![0, -1, 0, 0, 0],
690                vec![0, 1, 0, 0, 0],
691                vec![0, 0, 1, 10, 0],
692                vec![0, 0, 0, 0, 0],
693            ]),
694        );
695    }
696
697    #[test]
698    fn test_matrix_leader_sort_empty() {
699        let mut m = Matrix::empty();
700        m.leader_sort();
701        assert_eq!(m, Matrix::empty());
702    }
703
704    #[test]
705    fn test_matrix_eliminate_good() {
706        let mut m = Matrix::from_int_vecs(vec![
707            vec![1, 2, 3, 4, 5],
708            vec![6, 7, 8, 9, 10],
709            vec![11, 12, 13, 14, 15],
710            vec![16, 17, 18, 19, 20],
711            vec![21, 22, 23, 24, 25],
712        ]);
713
714        m.eliminate(2, 4, 0);
715        assert_eq!(
716            m,
717            Matrix::from_row_vecs(vec![
718                RowVec::from_int_vec(vec![1, 2, 3, 4, 5]),
719                RowVec::from_int_vec(vec![6, 7, 8, 9, 10]),
720                RowVec::from_int_vec(vec![11, 12, 13, 14, 15]),
721                RowVec::from_int_vec(vec![16, 17, 18, 19, 20]),
722                RowVec::new(vec![
723                    R64::ZERO,
724                    R64::new(-10, 11),
725                    R64::new(-20, 11),
726                    R64::new(-30, 11),
727                    R64::new(-40, 11)
728                ]),
729            ])
730        );
731
732        // TODO: verify
733        m.eliminate(4, 0, 1);
734        assert_eq!(
735            m,
736            Matrix::from_row_vecs(vec![
737                RowVec::from_int_vec(vec![1, 0, -1, -2, -3]),
738                RowVec::from_int_vec(vec![6, 7, 8, 9, 10]),
739                RowVec::from_int_vec(vec![11, 12, 13, 14, 15]),
740                RowVec::from_int_vec(vec![16, 17, 18, 19, 20]),
741                RowVec::new(vec![
742                    R64::ZERO,
743                    R64::new(-10, 11),
744                    R64::new(-20, 11),
745                    R64::new(-30, 11),
746                    R64::new(-40, 11)
747                ]),
748            ])
749        );
750
751        // TODO: verify
752        m.eliminate(2, 3, 2);
753        assert_eq!(
754            m,
755            Matrix::from_row_vecs(vec![
756                RowVec::from_int_vec(vec![1, 0, -1, -2, -3]),
757                RowVec::from_int_vec(vec![6, 7, 8, 9, 10]),
758                RowVec::from_int_vec(vec![11, 12, 13, 14, 15]),
759                RowVec::new(vec![
760                    R64::new(10, 13),
761                    R64::new(5, 13),
762                    R64::ZERO,
763                    R64::new(-5, 13),
764                    R64::new(-10, 13)
765                ]),
766                RowVec::new(vec![
767                    R64::ZERO,
768                    R64::new(-10, 11),
769                    R64::new(-20, 11),
770                    R64::new(-30, 11),
771                    R64::new(-40, 11)
772                ]),
773            ])
774        );
775    }
776
777    #[test]
778    fn test_matrix_eliminate_other_already_eliminated() {
779        let mut m = Matrix::from_int_vecs(vec![vec![0, 1, 2], vec![3, 0, 5]]);
780        let m2 = m.clone();
781        m.eliminate(0, 1, 1);
782        assert_eq!(m, m2);
783    }
784
785    #[test]
786    fn test_matrix_eliminate_self() {
787        let mut m = Matrix::from_int_vecs(vec![vec![0, 1, 2]]);
788        m.eliminate(0, 0, 1);
789        assert_eq!(m, Matrix::zeros(1, 3));
790    }
791
792    #[test]
793    #[should_panic]
794    fn test_matrix_eliminate_leader_is_zero() {
795        let mut m = Matrix::from_int_vecs(vec![vec![0, 1, 2], vec![3, 4, 5]]);
796        m.eliminate(0, 1, 0);
797    }
798
799    #[test]
800    fn test_matrix_normalize() {
801        let mut m = Matrix::from_int_vecs(vec![
802            vec![7, 5, 9],
803            vec![0, 0, 7, 5, 9],
804            vec![1, 2, 3],
805            vec![0, 0, 0],
806        ]);
807        m.normalize();
808        assert_eq!(
809            m,
810            Matrix::from_row_vecs(vec![
811                RowVec::new(vec![R64::ONE, R64::new(5, 7), R64::new(9, 7)]),
812                RowVec::new(vec![
813                    R64::ZERO,
814                    R64::ZERO,
815                    R64::ONE,
816                    R64::new(5, 7),
817                    R64::new(9, 7)
818                ]),
819                RowVec::from_int_vec(vec![1, 2, 3]),
820                RowVec::from_int_vec(vec![0, 0, 0]),
821            ])
822        );
823    }
824
825    #[test]
826    fn test_matrix_add_assign_good() -> Result<(), Box<dyn std::error::Error>> {
827        let mut m1 = Matrix::from_int_vecs(vec![vec![3, 5, 7], vec![8, 2, 3], vec![9, 2, 3]]);
828        let m2 = Matrix::from_int_vecs(vec![vec![8, 2, 4], vec![1, 6, 7], vec![2, 3, 5]]);
829        m1.add_assign(&m2)?;
830        assert_eq!(
831            m1,
832            Matrix::from_int_vecs(vec![vec![11, 7, 11], vec![9, 8, 10], vec![11, 5, 8]])
833        );
834        Ok(())
835    }
836
837    #[test]
838    fn test_matrix_add_assign_both_empty() -> Result<(), Box<dyn std::error::Error>> {
839        let mut m1 = Matrix::empty();
840        let m2 = Matrix::empty();
841        m1.add_assign(&m2)?;
842        assert_eq!(m1, Matrix::empty());
843        Ok(())
844    }
845
846    #[test]
847    fn test_matrix_add_assign_left_empty() {
848        let mut m1 = Matrix::empty();
849        let m1_copy = m1.clone();
850        let m2 = Matrix::from_int_vecs(vec![vec![8, 2, 4], vec![1, 6, 7]]);
851        assert!(m1.add_assign(&m2).is_err());
852        assert_eq!(m1, m1_copy);
853    }
854
855    #[test]
856    fn test_matrix_add_assign_right_empty() {
857        let mut m1 = Matrix::from_int_vecs(vec![vec![3, 5, 7], vec![8, 2, 3], vec![9, 2, 3]]);
858        let m1_copy = m1.clone();
859        let m2 = Matrix::empty();
860        assert!(m1.add_assign(&m2).is_err());
861        assert_eq!(m1, m1_copy);
862    }
863
864    #[test]
865    fn test_matrix_add_assign_row_mismatch() {
866        let mut m1 = Matrix::from_int_vecs(vec![vec![3, 5, 7], vec![8, 2, 3], vec![9, 2, 3]]);
867        let m1_copy = m1.clone();
868        let m2 = Matrix::from_int_vecs(vec![vec![8, 2, 4], vec![1, 6, 7]]);
869        assert!(m1.add_assign(&m2).is_err());
870        assert_eq!(m1, m1_copy);
871    }
872
873    #[test]
874    fn test_matrix_add_assign_col_mismatch() {
875        let mut m1 = Matrix::from_int_vecs(vec![vec![3, 5], vec![8, 2]]);
876        let m1_copy = m1.clone();
877        let m2 = Matrix::from_int_vecs(vec![vec![8, 2, 4], vec![1, 6, 7]]);
878        assert!(m1.add_assign(&m2).is_err());
879        assert_eq!(m1, m1_copy);
880    }
881
882    #[test]
883    fn test_matrix_add_good() -> Result<(), Box<dyn std::error::Error>> {
884        let m1 = Matrix::from_int_vecs(vec![vec![3, 5, 7], vec![8, 2, 3], vec![9, 2, 3]]);
885        let m2 = Matrix::from_int_vecs(vec![vec![8, 2, 4], vec![1, 6, 7], vec![2, 3, 5]]);
886        let result = m1.add(&m2)?;
887        assert_eq!(
888            result,
889            Matrix::from_int_vecs(vec![vec![11, 7, 11], vec![9, 8, 10], vec![11, 5, 8]])
890        );
891        Ok(())
892    }
893
894    #[test]
895    fn test_matrix_add_empty() -> Result<(), Box<dyn std::error::Error>> {
896        let m1 = Matrix::empty();
897        let m2 = Matrix::empty();
898        let result = m1.add(&m2)?;
899        assert_eq!(result, Matrix::empty());
900        Ok(())
901    }
902
903    #[test]
904    fn test_matrix_add_row_mismatch() {
905        let m1 = Matrix::from_int_vecs(vec![vec![3, 5, 7], vec![8, 2, 3], vec![9, 2, 3]]);
906        let m2 = Matrix::from_int_vecs(vec![vec![8, 2, 4], vec![1, 6, 7]]);
907        assert!(m1.add(&m2).is_err());
908    }
909
910    #[test]
911    fn test_matrix_add_col_mismatch() {
912        let m1 = Matrix::from_int_vecs(vec![vec![3, 5], vec![8, 2]]);
913        let m2 = Matrix::from_int_vecs(vec![vec![8, 2, 4], vec![1, 6, 7]]);
914        assert!(m1.add(&m2).is_err());
915    }
916
917    #[test]
918    fn test_matrix_matrix_mul_good() -> Result<(), Box<dyn std::error::Error>> {
919        let m1 = Matrix::from_int_vecs(vec![vec![3, 4], vec![7, 2], vec![1, 5]]);
920        let m2 = Matrix::from_int_vecs(vec![vec![5, 6, 2, 3], vec![8, 3, 9, 7]]);
921        let result = m1.matrix_mul(&m2)?;
922        assert_eq!(
923            result,
924            Matrix::from_int_vecs(vec![
925                vec![47, 30, 42, 37],
926                vec![51, 48, 32, 35],
927                vec![45, 21, 47, 38],
928            ])
929        );
930        Ok(())
931    }
932
933    #[test]
934    fn test_matrix_matrix_mul_both_empty() -> Result<(), Box<dyn std::error::Error>> {
935        let m1 = Matrix::empty();
936        let m2 = Matrix::empty();
937        let result = m1.matrix_mul(&m2)?;
938        assert_eq!(result, Matrix::empty());
939        Ok(())
940    }
941
942    #[test]
943    fn test_matrix_matrix_mul_left_empty() {
944        let m1 = Matrix::empty();
945        let m2 = Matrix::zeros(3, 4);
946        assert!(m1.matrix_mul(&m2).is_err());
947    }
948
949    #[test]
950    fn test_matrix_matrix_mul_right_empty() {
951        let m1 = Matrix::zeros(2, 3);
952        let m2 = Matrix::empty();
953        assert!(m1.matrix_mul(&m2).is_err());
954    }
955
956    #[test]
957    fn test_matrix_matrix_mul_mismatch() {
958        let m1 = Matrix::zeros(3, 3);
959        let m2 = Matrix::zeros(4, 4);
960        assert!(m1.matrix_mul(&m2).is_err());
961    }
962
963    #[test]
964    fn test_matrix_is_empty_true() {
965        assert!(Matrix::empty().is_empty());
966    }
967
968    #[test]
969    fn test_matrix_is_empty_false() {
970        assert!(!Matrix::zeros(1, 1).is_empty());
971    }
972
973    #[test]
974    fn test_matrix_height_good() {
975        assert_eq!(Matrix::zeros(2, 3).height(), 2);
976    }
977
978    #[test]
979    fn test_matrix_height_empty() {
980        assert_eq!(Matrix::empty().height(), 0);
981    }
982
983    #[test]
984    fn test_matrix_width_good() {
985        assert_eq!(Matrix::zeros(2, 3).width(), 3);
986    }
987
988    #[test]
989    fn test_matrix_width_empty() {
990        assert_eq!(Matrix::empty().width(), 0);
991    }
992
993    #[test]
994    fn test_matrix_append_row_good() {
995        let mut m = Matrix::zeros(2, 3);
996        let v = RowVec::zeros(3);
997        m.append_row(v);
998        assert_eq!(m, Matrix::zeros(3, 3));
999    }
1000
1001    #[test]
1002    fn test_matrix_append_row_empty() {
1003        let mut m = Matrix::empty();
1004        let v = RowVec::zeros(3);
1005        m.append_row(v);
1006        assert_eq!(m, Matrix::zeros(1, 3));
1007    }
1008
1009    #[test]
1010    fn test_matrix_remove_row_good() {
1011        let mut m = Matrix::from_int_vecs(vec![vec![1, 2, 3], vec![4, 5, 6], vec![7, 8, 9]]);
1012        m.remove_row(1);
1013        assert_eq!(m, Matrix::from_int_vecs(vec![vec![1, 2, 3], vec![7, 8, 9]]));
1014    }
1015
1016    #[test]
1017    fn test_matrix_remove_row_first() {
1018        let mut m = Matrix::from_int_vecs(vec![vec![1, 2, 3], vec![4, 5, 6], vec![7, 8, 9]]);
1019        m.remove_row(0);
1020        assert_eq!(m, Matrix::from_int_vecs(vec![vec![4, 5, 6], vec![7, 8, 9]]));
1021    }
1022
1023    #[test]
1024    fn test_matrix_remove_row_last() {
1025        let mut m = Matrix::from_int_vecs(vec![vec![1, 2, 3], vec![4, 5, 6], vec![7, 8, 9]]);
1026        m.remove_row(2);
1027        assert_eq!(m, Matrix::from_int_vecs(vec![vec![1, 2, 3], vec![4, 5, 6]]));
1028    }
1029
1030    #[test]
1031    fn test_matrix_remove_row_only() {
1032        let mut m = Matrix::from_int_vecs(vec![vec![1, 2, 3]]);
1033        m.remove_row(0);
1034        assert_eq!(m, Matrix::empty());
1035    }
1036
1037    #[test]
1038    fn test_matrix_is_zeros_true() {
1039        assert!(Matrix::zeros(2, 3).is_zeros());
1040    }
1041
1042    #[test]
1043    fn test_matrix_is_zeros_false() {
1044        assert!(!Matrix::from_int_vecs(vec![vec![1]]).is_zeros());
1045    }
1046
1047    #[test]
1048    fn test_matrix_is_zeros_empty() {
1049        assert!(Matrix::empty().is_zeros());
1050    }
1051
1052    #[test]
1053    fn test_matrix_get_column_copy_good() {
1054        let m = Matrix::from_int_vecs(vec![vec![1, 2, 3], vec![4, 5, 6], vec![7, 8, 9]]);
1055        assert_eq!(
1056            m.get_column_copy(1),
1057            Matrix::from_int_vecs(vec![vec![2], vec![5], vec![8]])
1058        );
1059    }
1060
1061    #[test]
1062    fn test_matrix_get_column_copy_first() {
1063        let m = Matrix::from_int_vecs(vec![vec![1, 2, 3], vec![4, 5, 6], vec![7, 8, 9]]);
1064        assert_eq!(
1065            m.get_column_copy(0),
1066            Matrix::from_int_vecs(vec![vec![1], vec![4], vec![7]])
1067        );
1068    }
1069
1070    #[test]
1071    fn test_matrix_get_column_copy_last() {
1072        let m = Matrix::from_int_vecs(vec![vec![1, 2, 3], vec![4, 5, 6], vec![7, 8, 9]]);
1073        assert_eq!(
1074            m.get_column_copy(2),
1075            Matrix::from_int_vecs(vec![vec![3], vec![6], vec![9]])
1076        );
1077    }
1078
1079    #[test]
1080    fn test_matrix_get_column_copy_only() {
1081        let m = Matrix::from_int_vecs(vec![vec![3], vec![6], vec![9]]);
1082        assert_eq!(
1083            m.get_column_copy(0),
1084            Matrix::from_int_vecs(vec![vec![3], vec![6], vec![9]])
1085        );
1086    }
1087
1088    #[test]
1089    fn test_matrix_eliminate_below_leader_good() {
1090        let mut m = Matrix::from_int_vecs(vec![vec![1, 2, 3], vec![4, 5, 6], vec![7, 8, 9]]);
1091        m.eliminate_below_leader(1);
1092
1093        // TODO: verify
1094        assert_eq!(
1095            m,
1096            Matrix::from_row_vecs(vec![
1097                RowVec::from_int_vec(vec![1, 2, 3]),
1098                RowVec::from_int_vec(vec![4, 5, 6]),
1099                RowVec::new(vec![R64::ZERO, R64::new(-3, 4), R64::new(-3, 2)]),
1100            ])
1101        );
1102    }
1103
1104    #[test]
1105    fn test_matrix_eliminate_below_leader_first() {
1106        let mut m = Matrix::from_int_vecs(vec![vec![1, 2, 3], vec![4, 5, 6], vec![7, 8, 9]]);
1107        m.eliminate_below_leader(0);
1108
1109        // TODO: verify
1110        assert_eq!(
1111            m,
1112            Matrix::from_int_vecs(vec![vec![1, 2, 3], vec![0, -3, -6], vec![0, -6, -12]])
1113        );
1114    }
1115
1116    #[test]
1117    fn test_matrix_eliminate_below_leader_last() {
1118        let mut m = Matrix::from_int_vecs(vec![vec![1, 2, 3], vec![4, 5, 6], vec![7, 8, 9]]);
1119        let m_copy = m.clone();
1120        m.eliminate_below_leader(2);
1121        assert_eq!(m, m_copy);
1122    }
1123
1124    #[test]
1125    fn test_matrix_eliminate_above_leader_good() {
1126        let mut m = Matrix::from_int_vecs(vec![vec![1, 2, 3], vec![4, 5, 6], vec![7, 8, 9]]);
1127        m.eliminate_above_leader(1);
1128
1129        // TODO: verify
1130        assert_eq!(
1131            m,
1132            Matrix::from_row_vecs(vec![
1133                RowVec::new(vec![R64::ZERO, R64::new(3, 4), R64::new(3, 2)]),
1134                RowVec::from_int_vec(vec![4, 5, 6]),
1135                RowVec::from_int_vec(vec![7, 8, 9]),
1136            ])
1137        );
1138    }
1139
1140    #[test]
1141    fn test_matrix_eliminate_above_leader_first() {
1142        let mut m = Matrix::from_int_vecs(vec![vec![1, 2, 3], vec![4, 5, 6], vec![7, 8, 9]]);
1143        let m_copy = m.clone();
1144        m.eliminate_above_leader(0);
1145        assert_eq!(m, m_copy);
1146    }
1147
1148    #[test]
1149    fn test_matrix_eliminate_above_leader_last() {
1150        let mut m = Matrix::from_int_vecs(vec![vec![1, 2, 3], vec![4, 5, 6], vec![7, 8, 9]]);
1151        m.eliminate_above_leader(2);
1152
1153        // TODO: verify
1154        assert_eq!(
1155            m,
1156            Matrix::from_row_vecs(vec![
1157                RowVec::new(vec![R64::ZERO, R64::new(6, 7), R64::new(12, 7)]),
1158                RowVec::new(vec![R64::ZERO, R64::new(3, 7), R64::new(6, 7)]),
1159                RowVec::from_int_vec(vec![7, 8, 9]),
1160            ])
1161        );
1162    }
1163
1164    #[test]
1165    fn test_matrix_ref_good() {
1166        let mut m = Matrix::from_int_vecs(vec![vec![1, 2, 3], vec![4, 5, 6], vec![7, 8, 9]]);
1167        m.r#ref();
1168        assert_eq!(
1169            m,
1170            Matrix::from_int_vecs(vec![vec![1, 2, 3], vec![0, -3, -6], vec![0, 0, 0]])
1171        );
1172    }
1173
1174    #[test]
1175    fn test_matrix_ref_empty() {
1176        let mut m = Matrix::empty();
1177        m.r#ref();
1178        assert_eq!(m, Matrix::empty());
1179    }
1180
1181    // TODO: more ref tests?
1182
1183    #[test]
1184    fn test_matrix_rref_good() {
1185        let mut m = Matrix::from_int_vecs(vec![vec![1, 2, 3], vec![4, 5, 6], vec![7, 8, 9]]);
1186        m.rref();
1187        assert_eq!(
1188            m,
1189            Matrix::from_int_vecs(vec![vec![1, 0, -1], vec![0, 1, 2], vec![0, 0, 0]])
1190        );
1191    }
1192
1193    #[test]
1194    fn test_matrix_rref_empty() {
1195        let mut m = Matrix::empty();
1196        m.rref();
1197        assert_eq!(m, Matrix::empty());
1198    }
1199
1200    // TODO: more rref tests?
1201
1202    #[test]
1203    fn test_matrix_mul_assign_good() {
1204        let mut m = Matrix::from_int_vecs(vec![vec![1, 2, 3], vec![4, 5, 6], vec![7, 8, 9]]);
1205        m *= R64::from_integer(5);
1206        assert_eq!(
1207            m,
1208            Matrix::from_int_vecs(vec![vec![5, 10, 15], vec![20, 25, 30], vec![35, 40, 45]])
1209        );
1210    }
1211
1212    #[test]
1213    fn test_matrix_mul_assign_empty() {
1214        let mut m = Matrix::empty();
1215        m *= R64::from_integer(5);
1216        assert_eq!(m, Matrix::empty());
1217    }
1218
1219    #[test]
1220    fn test_matrix_mul_good() {
1221        let m = Matrix::from_int_vecs(vec![vec![1, 2, 3], vec![4, 5, 6], vec![7, 8, 9]]);
1222        let result = m * R64::from_integer(5);
1223        assert_eq!(
1224            result,
1225            Matrix::from_int_vecs(vec![vec![5, 10, 15], vec![20, 25, 30], vec![35, 40, 45]])
1226        );
1227    }
1228
1229    #[test]
1230    fn test_matrix_mul_empty() {
1231        let m = Matrix::empty();
1232        let result = m * R64::from_integer(5);
1233        assert_eq!(result, Matrix::empty());
1234    }
1235
1236    #[test]
1237    fn test_matrix_index() {
1238        assert_eq!(
1239            Matrix::from_int_vecs(vec![vec![1, 2, 3], vec![4, 5, 6], vec![7, 8, 9]])[(1, 2)],
1240            R64::from_integer(6)
1241        );
1242    }
1243
1244    #[test]
1245    fn test_matrix_index_mut() {
1246        let mut m = Matrix::from_int_vecs(vec![vec![1, 2, 3], vec![4, 5, 6], vec![7, 8, 9]]);
1247        m[(1, 2)] = R64::ZERO;
1248        assert_eq!(
1249            m,
1250            Matrix::from_int_vecs(vec![vec![1, 2, 3], vec![4, 5, 0], vec![7, 8, 9]])
1251        );
1252    }
1253}
aoclib_rs/matrix.rs

aoclib_rs/
matrix.rs
