diff --git a/src/decompose.rs b/src/decompose.rs index 679d7ed..d618708 100644 --- a/src/decompose.rs +++ b/src/decompose.rs @@ -124,7 +124,7 @@ impl LUDecomposition { /// This is equivalent to [`LUDecompose::det`] while allowing the LU decomposition /// to be reused pub fn det(&self) -> T { - self.parity * self.lu.diagonals().fold(T::one(), T::mul) + self.parity * self.lu.diagonals().fold(T::one(), |l, &r| l * r) } /// Calculate the inverse of the original matrix, such that $bbM xx bbM^{-1} = bbI$ diff --git a/src/identities.rs b/src/identities.rs deleted file mode 100644 index 1b87721..0000000 --- a/src/identities.rs +++ /dev/null @@ -1,66 +0,0 @@ -use crate::Matrix; -use num_traits::{Bounded, One, Zero}; - -// Identity -impl Matrix { - /// Create an identity matrix, a square matrix where the diagonals are 1 and all other elements - /// are 0. - /// for example, - /// - /// $bbI = [[1,0,0],[0,1,0],[0,0,1]]$ - /// - /// Matrix multiplication between a matrix and the identity matrix always results in itself - /// - /// $bbA xx bbI = bbA$ - /// - /// # Examples - /// ``` - /// # use vector_victor::Matrix; - /// let i = Matrix::::identity(); - /// assert_eq!(i, Matrix::mat([[1,0,0],[0,1,0],[0,0,1]])) - /// ``` - /// - /// Note that the identity only exists for matrices that are square, so this doesnt work: - /// ```compile_fail - /// # use vector_victor::Matrix; - /// let i = Matrix::::identity(); - /// ``` - #[must_use] - pub fn identity() -> Self { - let mut result = Self::zero(); - for i in 0..N { - result[(i, i)] = T::one(); - } - return result; - } -} - -// Zero -impl Zero for Matrix { - fn zero() -> Self { - Matrix::fill(T::zero()) - } - - fn is_zero(&self) -> bool { - self.elements().all(|e| e.is_zero()) - } -} - -// One -impl One for Matrix { - fn one() -> Self { - Matrix::fill(T::one()) - } -} - -// min_value and max_value -// LowerBounded and UpperBounded are automatically implemented from this -impl Bounded for Matrix { - fn min_value() -> Self { - Self::fill(T::min_value()) - } - - fn max_value() -> Self { - Self::fill(T::max_value()) - } -} diff --git a/src/index.rs b/src/index.rs index 32a153c..981f819 100644 --- a/src/index.rs +++ b/src/index.rs @@ -19,7 +19,7 @@ use std::fmt::Debug; /// assert_eq!(m[7], 8); /// /// let v = Vector::vec([4,8,15,16,23,42]); -/// assert_eq!(m[2], 15); // just like a std::vec +/// assert_eq!(v[2], 15); // just like a std::vec /// ``` /// /// Indexing by a `(usize,usize)` indexes by row and column diff --git a/src/lib.rs b/src/lib.rs index 3819f50..68368c5 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,13 +1,13 @@ extern crate core; use index::Index2D; +use num_traits::{Bounded, One, Zero}; use std::cmp::min; use std::fmt::Debug; use std::iter::{zip, Flatten}; use std::ops::{Index, IndexMut}; pub mod decompose; -mod identities; pub mod index; mod math; mod ops; @@ -37,6 +37,72 @@ impl Default for Matrix Zero for Matrix { + fn zero() -> Self { + Matrix::fill(T::zero()) + } + + fn is_zero(&self) -> bool { + self.elements().all(|e| e.is_zero()) + } +} + +// One +impl One for Matrix { + fn one() -> Self { + Matrix::fill(T::one()) + } +} + +// min_value and max_value +// LowerBounded and UpperBounded are automatically implemented from this +impl Bounded for Matrix { + fn min_value() -> Self { + Self::fill(T::min_value()) + } + + fn max_value() -> Self { + Self::fill(T::max_value()) + } +} + +// Identity +impl Matrix { + /// Create an identity matrix, a square matrix where the diagonals are 1 and all other elements + /// are 0. + /// for example, + /// + /// $bbI = [[1,0,0],[0,1,0],[0,0,1]]$ + /// + /// Matrix multiplication between a matrix and the identity matrix always results in itself + /// + /// $bbA xx bbI = bbA$ + /// + /// # Examples + /// ``` + /// # use vector_victor::Matrix; + /// let i = Matrix::::identity(); + /// assert_eq!(i, Matrix::mat([[1, 0, 0], + /// [0, 1, 0], + /// [0, 0, 1]])) + /// ``` + /// + /// Note that the identity only exists for matrices that are square, so this doesnt work: + /// ```compile_fail + /// # use vector_victor::Matrix; + /// let i = Matrix::::identity(); + /// ``` + #[must_use] + pub fn identity() -> Self { + let mut result = Self::zero(); + for i in 0..N { + result[(i, i)] = T::one(); + } + return result; + } +} + // Matrix constructors impl Matrix { /// Generate a new matrix from a 2D Array @@ -45,8 +111,6 @@ impl Matrix { /// /// * `data`: A 2D array of elements to copy into the new matrix /// - /// returns: Matrix - /// /// # Examples /// /// ``` @@ -66,15 +130,12 @@ impl Matrix { /// /// * `scalar`: Scalar value to copy into the new matrix. /// - /// returns: Matrix - /// /// # Examples /// /// ``` /// # use vector_victor::Matrix; - /// let my_matrix = Matrix::::fill(5); - /// // is equivalent to - /// assert_eq!(my_matrix, Matrix::mat([[5;4];4])) + /// // these are equivalent + /// assert_eq!(Matrix::::fill(5), Matrix::mat([[5;4];4])) /// ``` #[must_use] pub fn fill(scalar: T) -> Matrix { @@ -91,15 +152,19 @@ impl Matrix { /// /// * `iter`: iterator of vectors to copy into rows /// - /// returns: Matrix - /// /// # Examples /// + /// The following is another way of performing [`Matrix::transpose()`] /// ``` /// # use vector_victor::Matrix; - /// let my_matrix = Matrix::mat([[1,2,3],[4,5,6]]); + /// let my_matrix = Matrix::mat([[1, 2, 3], + /// [4, 5, 6]]); + /// /// let transpose : Matrix<_,3,2>= Matrix::from_rows(my_matrix.cols()); - /// assert_eq!(transpose, Matrix::mat([[1,4],[2,5],[3,6]])) + /// + /// assert_eq!(transpose, Matrix::mat([[1, 4], + /// [2, 5], + /// [3, 6]])) /// ``` #[must_use] pub fn from_rows(iter: I) -> Self @@ -120,15 +185,19 @@ impl Matrix { /// /// * `iter`: iterator of vectors to copy into columns /// - /// returns: Matrix - /// /// # Examples /// + /// The following is another way of performing [`Matrix::transpose()`] /// ``` /// # use vector_victor::Matrix; - /// let my_matrix = Matrix::mat([[1,2,3],[4,5,6]]); + /// let my_matrix = Matrix::mat([[1, 2, 3], + /// [4, 5, 6]]); + /// /// let transpose : Matrix<_,3,2>= Matrix::from_cols(my_matrix.rows()); - /// assert_eq!(transpose, Matrix::mat([[1,4],[2,5],[3,6]])) + /// + /// assert_eq!(transpose, Matrix::mat([[1, 4], + /// [2, 5], + /// [3, 6]])) /// ``` #[must_use] pub fn from_cols(iter: I) -> Self @@ -152,9 +221,8 @@ impl Vector { /// # Examples /// ``` /// # use vector_victor::{Matrix, Vector}; - /// let my_vector = Vector::vec([1,2,3,4]); - /// // is equivalent to - /// assert_eq!(my_vector, Matrix::mat([[1],[2],[3],[4]])); + /// // these are equivalent + /// assert_eq!(Vector::vec([1,2,3,4]), Matrix::mat([[1],[2],[3],[4]])); /// ``` pub fn vec(data: [T; N]) -> Self { assert!(N > 0, "Vector must have at least 1 element"); @@ -168,44 +236,89 @@ impl Vector { impl Matrix { /// Returns an iterator over the elements of the matrix in row-major order. /// + /// This is identical to the behavior of [`IntoIterator`](#associatedtype.IntoIter) + /// /// # Examples /// ``` /// # use vector_victor::Matrix; - /// let my_matrix = Matrix::mat([[1,2],[3,4]]); - /// assert!(vec![1,2,3,4].iter().eq(my_matrix.elements())) + /// let my_matrix = Matrix::mat([[1, 2], + /// [3, 4]]); + /// + /// itertools::assert_equal(my_matrix.elements(), [1,2,3,4].iter()) /// ``` #[must_use] - pub fn elements<'a>(&'a self) -> impl Iterator + 'a { + pub fn elements<'s>(&'s self) -> impl Iterator + 's { self.data.iter().flatten() } /// Returns a mutable iterator over the elements of the matrix in row-major order. + /// + /// # Examples + /// ``` + /// # use vector_victor::Matrix; + /// let mut my_matrix = Matrix::mat([[1, 2], + /// [3, 4]]); + /// + /// for elem in my_matrix.elements_mut() {*elem += 2;} + /// itertools::assert_equal(my_matrix.elements(), [3,4,5,6].iter()) + /// ``` #[must_use] - pub fn elements_mut<'a>(&'a mut self) -> impl Iterator + 'a { + pub fn elements_mut<'s>(&'s mut self) -> impl Iterator + 's { self.data.iter_mut().flatten() } /// returns an iterator over the elements along the diagonal of a matrix + /// + /// # Examples + /// ``` + /// # use vector_victor::Matrix; + /// let my_matrix = Matrix::mat([[1, 2, 3], + /// [4, 5, 6], + /// [7, 8, 9], + /// [10,11,12]]); + /// + /// itertools::assert_equal(my_matrix.diagonals(), [1,5,9].iter()) + /// ``` #[must_use] - pub fn diagonals<'s>(&'s self) -> impl Iterator + 's { - (0..min(N, M)).map(|n| self[(n, n)]) + pub fn diagonals<'s>(&'s self) -> impl Iterator + 's { + (0..min(N, M)).map(|n| &self[(n, n)]) } /// Returns an iterator over the elements directly below the diagonal of a matrix + /// returns an iterator over the elements along the diagonal of a matrix + /// + /// # Examples + /// ``` + /// # use vector_victor::Matrix; + /// let my_matrix = Matrix::mat([[1, 2, 3], + /// [4, 5, 6], + /// [7, 8, 9], + /// [10,11,12]]); + /// + /// itertools::assert_equal(my_matrix.subdiagonals(), [4,8,12].iter()); + /// ``` #[must_use] - pub fn subdiagonals<'s>(&'s self) -> impl Iterator + 's { - (0..min(N, M) - 1).map(|n| self[(n, n + 1)]) + pub fn subdiagonals<'s>(&'s self) -> impl Iterator + 's { + (0..min(N, M - 1)).map(|n| &self[(n + 1, n)]) } /// Returns a reference to the element at that position in the matrix, or `None` if out of bounds. /// + /// [`Index`](#impl-Index%3CI%3E-for-Matrix%3CT,+M,+N%3E) behaves similarly, + /// but will panic if the index is out of bounds instead of returning an option + /// + /// # Arguments + /// + /// * `index`: a 1D or 2D index into the matrix. See [Index2D] for more information on matrix indexing. + /// /// # Examples /// /// ``` /// # use vector_victor::Matrix; - /// let my_matrix = Matrix::mat([[1,2],[3,4]]); + /// let my_matrix = Matrix::mat([[1, 2], + /// [3, 4]]); /// - /// // element at index 2 is the same as the element at (row 1, column 0). + /// // element at index 2 is the same as the element at row 1, column 0. /// assert_eq!(my_matrix.get(2), my_matrix.get((1,0))); /// /// // my_matrix.get() is equivalent to my_matrix[], @@ -222,7 +335,29 @@ impl Matrix { Some(&self.data[m][n]) } - /// Returns a mutable reference to the element at that position in the matrix, or `None` if out of bounds. + /// Returns a mutable reference to the element at that position in the matrix, + /// or `None` if out of bounds. + /// + /// [`IndexMut`](#impl-IndexMut%3CI%3E-for-Matrix%3CT,+M,+N%3E) behaves similarly, + /// but will panic if the index is out of bounds instead of returning an option + /// + /// # Arguments + /// + /// * `index`: a 1D or 2D index into the matrix. See [Index2D] for more information + /// on matrix indexing. + /// + /// # Examples + /// + /// ``` + /// # use vector_victor::Matrix; + /// let mut my_matrix = Matrix::mat([[1, 2], + /// [3, 4]]); + /// + /// match my_matrix.get_mut(2) { + /// Some(t) => *t = 5, + /// None => panic!()}; + /// assert_eq!(my_matrix, Matrix::mat([[1,2],[5,4]])) + /// ``` #[inline] #[must_use] pub fn get_mut(&mut self, index: impl Index2D) -> Option<&mut T> { @@ -230,13 +365,18 @@ impl Matrix { Some(&mut self.data[m][n]) } - /// Returns a row of the matrix. or [None] if index is out of bounds + /// Returns a row of the matrix. + /// + /// # Panics + /// + /// Panics if row index `m` is out of bounds. /// /// # Examples /// /// ``` /// # use vector_victor::{Matrix, Vector}; - /// let my_matrix = Matrix::mat([[1,2],[3,4]]); + /// let my_matrix = Matrix::mat([[1, 2], + /// [3, 4]]); /// /// // row at index 1 /// assert_eq!(my_matrix.row(1), Vector::vec([3,4])); @@ -246,7 +386,7 @@ impl Matrix { pub fn row(&self, m: usize) -> Vector { assert!( m < M, - "Row index {} out of bounds for {}x{} matrix", + "Row index {} out of bounds for {}×{} matrix", m, M, N @@ -254,11 +394,27 @@ impl Matrix { Vector::::vec(self.data[m]) } + /// Sets a row of the matrix. + /// + /// # Panics + /// + /// Panics if row index `m` is out of bounds. + /// + /// # Examples + /// + /// ``` + /// # use vector_victor::{Matrix, Vector}; + /// let mut my_matrix = Matrix::mat([[1, 2], + /// [3, 4]]); + /// // row at index 1 + /// my_matrix.set_row(1, &Vector::vec([5,6])); + /// assert_eq!(my_matrix, Matrix::mat([[1,2],[5,6]])); + /// ``` #[inline] pub fn set_row(&mut self, m: usize, val: &Vector) { assert!( m < M, - "Row index {} out of bounds for {}x{} matrix", + "Row index {} out of bounds for {}×{} matrix", m, M, N @@ -268,18 +424,28 @@ impl Matrix { } } - pub fn pivot_row(&mut self, m1: usize, m2: usize) { - let tmp = self.row(m2); - self.set_row(m2, &self.row(m1)); - self.set_row(m1, &tmp); - } - + /// Returns a column of the matrix. + /// + /// # Panics + /// + /// Panics if column index `n` is out of bounds. + /// + /// # Examples + /// + /// ``` + /// # use vector_victor::{Matrix, Vector}; + /// let my_matrix = Matrix::mat([[1, 2], + /// [3, 4]]); + /// + /// // column at index 1 + /// assert_eq!(my_matrix.col(1), Vector::vec([2,4])); + /// ``` #[inline] #[must_use] pub fn col(&self, n: usize) -> Vector { assert!( n < N, - "Column index {} out of bounds for {}x{} matrix", + "Column index {} out of bounds for {}×{} matrix", n, M, N @@ -287,11 +453,27 @@ impl Matrix { Vector::::vec(self.data.map(|r| r[n])) } + /// Sets a column of the matrix. + /// + /// # Panics + /// + /// Panics if column index `n` is out of bounds. + /// + /// # Examples + /// + /// ``` + /// # use vector_victor::{Matrix, Vector}; + /// let mut my_matrix = Matrix::mat([[1, 2], + /// [3, 4]]); + /// // column at index 1 + /// my_matrix.set_col(1, &Vector::vec([5,6])); + /// assert_eq!(my_matrix, Matrix::mat([[1,5],[3,6]])); + /// ``` #[inline] pub fn set_col(&mut self, n: usize, val: &Vector) { assert!( n < N, - "Column index {} out of bounds for {}x{} matrix", + "Column index {} out of bounds for {}×{} matrix", n, M, N @@ -302,22 +484,64 @@ impl Matrix { } } - pub fn pivot_col(&mut self, n1: usize, n2: usize) { - let tmp = self.col(n2); - self.set_col(n2, &self.col(n1)); - self.set_col(n1, &tmp); - } - + /// Returns an iterator over the rows of the matrix, returning them as column vectors. #[must_use] pub fn rows<'a>(&'a self) -> impl Iterator> + 'a { (0..M).map(|m| self.row(m)) } + /// Returns an iterator over the columns of the matrix, returning them as column vectors. #[must_use] pub fn cols<'a>(&'a self) -> impl Iterator> + 'a { (0..N).map(|n| self.col(n)) } + /// Interchange two rows + /// + /// # Panics + /// + /// Panics if row index `m1` or `m2` are out of bounds + pub fn pivot_row(&mut self, m1: usize, m2: usize) { + let tmp = self.row(m2); + self.set_row(m2, &self.row(m1)); + self.set_row(m1, &tmp); + } + + /// Interchange two columns + /// + /// # Panics + /// + /// Panics if column index `n1` or `n2` are out of bounds + pub fn pivot_col(&mut self, n1: usize, n2: usize) { + let tmp = self.col(n2); + self.set_col(n2, &self.col(n1)); + self.set_col(n1, &tmp); + } + + /// Apply a permutation matrix to the rows of a matrix + /// + /// # Arguments + /// + /// * `ms`: a [`Vector`] of [`usize`] of length M. Each entry is the index of the row that will + /// appear in the result + /// + /// # Panics + /// + /// Panics if any of the row indices in `ms` is out of bounds + /// + /// # Examples + /// + /// ``` + /// # use vector_victor::{Matrix, Vector}; + /// let my_matrix = Matrix::mat([[1, 2, 3], + /// [4, 5, 6], + /// [7, 8, 9]]); + /// + /// let permuted = my_matrix.permute_rows(&Vector::vec([1, 0, 2])); + /// assert_eq!(permuted, Matrix::mat([[4, 5, 6], + /// [1, 2, 3], + /// [7, 8, 9]])) + /// ``` #[must_use] pub fn permute_rows(&self, ms: &Vector) -> Self where @@ -326,6 +550,16 @@ impl Matrix { Self::from_rows(ms.elements().map(|&m| self.row(m))) } + /// Apply a permutation matrix to the columns of a matrix + /// + /// # Arguments + /// + /// * `ns`: a [`Vector`] of [`usize`] of length N. Each entry is the index of the column that will + /// appear in the result + /// + /// # Panics + /// + /// Panics if any of the column indices in `ns` is out of bounds #[must_use] pub fn permute_cols(&self, ns: &Vector) -> Self where @@ -334,6 +568,20 @@ impl Matrix { Self::from_cols(ns.elements().map(|&n| self.col(n))) } + /// Returns the transpose $M^T$ of the matrix, or the matrix flipped across its diagonal. + /// + /// # Examples + /// ``` + /// # use vector_victor::Matrix; + /// let my_matrix = Matrix::mat([[1, 2, 3], + /// [4, 5, 6]]); + /// + /// assert_eq!( + /// my_matrix.transpose(), + /// Matrix::mat([[1, 4], + /// [2, 5], + /// [3, 6]])) + /// ``` pub fn transpose(&self) -> Matrix where Matrix: Default, @@ -353,7 +601,7 @@ where #[inline(always)] fn index(&self, index: I) -> &Self::Output { self.get(index).expect(&*format!( - "index {:?} out of range for {}x{} Matrix", + "index {:?} out of range for {}×{} Matrix", index, M, N )) } @@ -368,7 +616,7 @@ where #[inline(always)] fn index_mut(&mut self, index: I) -> &mut Self::Output { self.get_mut(index).expect(&*format!( - "index {:?} out of range for {}x{} Matrix", + "index {:?} out of range for {}×{} Matrix", index, M, N )) } diff --git a/src/mod.rs b/src/mod.rs deleted file mode 100644 index 8b13789..0000000 --- a/src/mod.rs +++ /dev/null @@ -1 +0,0 @@ - diff --git a/tests/decompose.rs b/tests/decompose.rs index 9c904d8..3a07d56 100644 --- a/tests/decompose.rs +++ b/tests/decompose.rs @@ -1,3 +1,5 @@ +//! Data structures and traits for decomposing and solving matrices + #[macro_use] mod common;