Compare commits
2 Commits
bc1b3f199d
...
e2a2bc7529
| Author | SHA1 | Date |
|---|---|---|
 Andrew Cassidy
|
e2a2bc7529 | 11 months ago |
|
Andrew Cassidy
|
2b303892f7 | 12 months ago |
@ -0,0 +1,66 @@
|
||||
use crate::Matrix;
|
||||
use num_traits::{Bounded, One, Zero};
|
||||
|
||||
// Identity
|
||||
impl<T: Copy + Zero + One, const N: usize> Matrix<T, N, N> {
|
||||
/// 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::<i32,3,3>::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::<i32,4,2>::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<T: Copy + Zero, const M: usize, const N: usize> Zero for Matrix<T, M, N> {
|
||||
fn zero() -> Self {
|
||||
Matrix::fill(T::zero())
|
||||
}
|
||||
|
||||
fn is_zero(&self) -> bool {
|
||||
self.elements().all(|e| e.is_zero())
|
||||
}
|
||||
}
|
||||
|
||||
// One
|
||||
impl<T: Copy + One, const M: usize, const N: usize> One for Matrix<T, M, N> {
|
||||
fn one() -> Self {
|
||||
Matrix::fill(T::one())
|
||||
}
|
||||
}
|
||||
|
||||
// min_value and max_value
|
||||
// LowerBounded and UpperBounded are automatically implemented from this
|
||||
impl<T: Copy + Bounded, const N: usize, const M: usize> Bounded for Matrix<T, N, M> {
|
||||
fn min_value() -> Self {
|
||||
Self::fill(T::min_value())
|
||||
}
|
||||
|
||||
fn max_value() -> Self {
|
||||
Self::fill(T::max_value())
|
||||
}
|
||||
}
|
||||
@ -1,7 +1,422 @@
|
||||
extern crate core;
|
||||
|
||||
use index::Index2D;
|
||||
use std::cmp::min;
|
||||
use std::fmt::Debug;
|
||||
use std::iter::{zip, Flatten};
|
||||
use std::ops::{Index, IndexMut};
|
||||
|
||||
pub mod decompose;
|
||||
mod matrix;
|
||||
mod identities;
|
||||
pub mod index;
|
||||
mod math;
|
||||
mod ops;
|
||||
|
||||
mod util;
|
||||
|
||||
pub use matrix::{Matrix, Vector};
|
||||
/// A 2D array of values which can be operated upon.
|
||||
///
|
||||
/// Matrices have a fixed size known at compile time
|
||||
#[derive(Debug, Copy, Clone, PartialEq)]
|
||||
pub struct Matrix<T, const M: usize, const N: usize>
|
||||
where
|
||||
T: Copy,
|
||||
{
|
||||
data: [[T; N]; M], // Row-Major order
|
||||
}
|
||||
|
||||
/// An alias for a [Matrix] with a single column
|
||||
pub type Vector<T, const N: usize> = Matrix<T, N, 1>;
|
||||
|
||||
// CONSTRUCTORS
|
||||
|
||||
// Default
|
||||
impl<T: Copy + Default, const M: usize, const N: usize> Default for Matrix<T, M, N> {
|
||||
fn default() -> Self {
|
||||
Matrix::fill(T::default())
|
||||
}
|
||||
}
|
||||
|
||||
// Matrix constructors
|
||||
impl<T: Copy, const M: usize, const N: usize> Matrix<T, M, N> {
|
||||
/// Generate a new matrix from a 2D Array
|
||||
///
|
||||
/// # Arguments
|
||||
///
|
||||
/// * `data`: A 2D array of elements to copy into the new matrix
|
||||
///
|
||||
/// returns: Matrix<T, M, N>
|
||||
///
|
||||
/// # Examples
|
||||
///
|
||||
/// ```
|
||||
/// # use vector_victor::Matrix;
|
||||
/// let a = Matrix::mat([[1,2,3,4];4]);
|
||||
/// ```
|
||||
#[must_use]
|
||||
pub fn mat(data: [[T; N]; M]) -> Self {
|
||||
assert!(M > 0, "Matrix must have at least 1 row");
|
||||
assert!(N > 0, "Matrix must have at least 1 column");
|
||||
Matrix::<T, M, N> { data }
|
||||
}
|
||||
|
||||
/// Generate a new matrix from a single scalar
|
||||
///
|
||||
/// # Arguments
|
||||
///
|
||||
/// * `scalar`: Scalar value to copy into the new matrix.
|
||||
///
|
||||
/// returns: Matrix<T, M, N>
|
||||
///
|
||||
/// # Examples
|
||||
///
|
||||
/// ```
|
||||
/// # use vector_victor::Matrix;
|
||||
/// let my_matrix = Matrix::<i32,4,4>::fill(5);
|
||||
/// // is equivalent to
|
||||
/// assert_eq!(my_matrix, Matrix::mat([[5;4];4]))
|
||||
/// ```
|
||||
#[must_use]
|
||||
pub fn fill(scalar: T) -> Matrix<T, M, N> {
|
||||
assert!(M > 0, "Matrix must have at least 1 row");
|
||||
assert!(N > 0, "Matrix must have at least 1 column");
|
||||
Matrix::<T, M, N> {
|
||||
data: [[scalar; N]; M],
|
||||
}
|
||||
}
|
||||
|
||||
/// Create a matrix from an iterator of vectors
|
||||
///
|
||||
/// # Arguments
|
||||
///
|
||||
/// * `iter`: iterator of vectors to copy into rows
|
||||
///
|
||||
/// returns: Matrix<T, M, N>
|
||||
///
|
||||
/// # Examples
|
||||
///
|
||||
/// ```
|
||||
/// # use vector_victor::Matrix;
|
||||
/// 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]]))
|
||||
/// ```
|
||||
#[must_use]
|
||||
pub fn from_rows<I>(iter: I) -> Self
|
||||
where
|
||||
I: IntoIterator<Item = Vector<T, N>>,
|
||||
Self: Default,
|
||||
{
|
||||
let mut result = Self::default();
|
||||
for (m, row) in iter.into_iter().enumerate().take(M) {
|
||||
result.set_row(m, &row)
|
||||
}
|
||||
result
|
||||
}
|
||||
|
||||
/// Create a matrix from an iterator of vectors
|
||||
///
|
||||
/// # Arguments
|
||||
///
|
||||
/// * `iter`: iterator of vectors to copy into columns
|
||||
///
|
||||
/// returns: Matrix<T, M, N>
|
||||
///
|
||||
/// # Examples
|
||||
///
|
||||
/// ```
|
||||
/// # use vector_victor::Matrix;
|
||||
/// 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]]))
|
||||
/// ```
|
||||
#[must_use]
|
||||
pub fn from_cols<I>(iter: I) -> Self
|
||||
where
|
||||
I: IntoIterator<Item = Vector<T, M>>,
|
||||
Self: Default,
|
||||
{
|
||||
let mut result = Self::default();
|
||||
for (n, col) in iter.into_iter().enumerate().take(N) {
|
||||
result.set_col(n, &col)
|
||||
}
|
||||
result
|
||||
}
|
||||
}
|
||||
|
||||
// Vector constructor
|
||||
impl<T: Copy, const N: usize> Vector<T, N> {
|
||||
/// Create a vector from a 1D array.
|
||||
/// Note that vectors are always column vectors unless explicitly instantiated as row vectors
|
||||
///
|
||||
/// # 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]]));
|
||||
/// ```
|
||||
pub fn vec(data: [T; N]) -> Self {
|
||||
assert!(N > 0, "Vector must have at least 1 element");
|
||||
return Vector::<T, N> {
|
||||
data: data.map(|e| [e]),
|
||||
};
|
||||
}
|
||||
}
|
||||
|
||||
// ACCESSORS AND MUTATORS
|
||||
impl<T: Copy, const M: usize, const N: usize> Matrix<T, M, N> {
|
||||
/// Returns an iterator over the elements of the matrix in row-major order.
|
||||
///
|
||||
/// # 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()))
|
||||
/// ```
|
||||
#[must_use]
|
||||
pub fn elements<'a>(&'a self) -> impl Iterator<Item = &'a T> + 'a {
|
||||
self.data.iter().flatten()
|
||||
}
|
||||
|
||||
/// Returns a mutable iterator over the elements of the matrix in row-major order.
|
||||
#[must_use]
|
||||
pub fn elements_mut<'a>(&'a mut self) -> impl Iterator<Item = &'a mut T> + 'a {
|
||||
self.data.iter_mut().flatten()
|
||||
}
|
||||
|
||||
/// returns an iterator over the elements along the diagonal of a matrix
|
||||
#[must_use]
|
||||
pub fn diagonals<'s>(&'s self) -> impl Iterator<Item = T> + 's {
|
||||
(0..min(N, M)).map(|n| self[(n, n)])
|
||||
}
|
||||
|
||||
/// Returns an iterator over the elements directly below the diagonal of a matrix
|
||||
#[must_use]
|
||||
pub fn subdiagonals<'s>(&'s self) -> impl Iterator<Item = T> + 's {
|
||||
(0..min(N, M) - 1).map(|n| self[(n, n + 1)])
|
||||
}
|
||||
|
||||
/// Returns a reference to the element at that position in the matrix, or `None` if out of bounds.
|
||||
///
|
||||
/// # Examples
|
||||
///
|
||||
/// ```
|
||||
/// # use vector_victor::Matrix;
|
||||
/// let my_matrix = Matrix::mat([[1,2],[3,4]]);
|
||||
///
|
||||
/// // 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[],
|
||||
/// // but returns an Option instead of panicking
|
||||
/// assert_eq!(my_matrix.get(2), Some(&my_matrix[2]));
|
||||
///
|
||||
/// // index 4 is out of range, so get(4) returns None.
|
||||
/// assert_eq!(my_matrix.get(4), None);
|
||||
/// ```
|
||||
#[inline]
|
||||
#[must_use]
|
||||
pub fn get(&self, index: impl Index2D) -> Option<&T> {
|
||||
let (m, n) = index.to_2d(M, N)?;
|
||||
Some(&self.data[m][n])
|
||||
}
|
||||
|
||||
/// Returns a mutable reference to the element at that position in the matrix, or `None` if out of bounds.
|
||||
#[inline]
|
||||
#[must_use]
|
||||
pub fn get_mut(&mut self, index: impl Index2D) -> Option<&mut T> {
|
||||
let (m, n) = index.to_2d(M, N)?;
|
||||
Some(&mut self.data[m][n])
|
||||
}
|
||||
|
||||
/// Returns a row of the matrix. or [None] if index is out of bounds
|
||||
///
|
||||
/// # Examples
|
||||
///
|
||||
/// ```
|
||||
/// # use vector_victor::{Matrix, Vector};
|
||||
/// let my_matrix = Matrix::mat([[1,2],[3,4]]);
|
||||
///
|
||||
/// // row at index 1
|
||||
/// assert_eq!(my_matrix.row(1), Vector::vec([3,4]));
|
||||
/// ```
|
||||
#[inline]
|
||||
#[must_use]
|
||||
pub fn row(&self, m: usize) -> Vector<T, N> {
|
||||
assert!(
|
||||
m < M,
|
||||
"Row index {} out of bounds for {}x{} matrix",
|
||||
m,
|
||||
M,
|
||||
N
|
||||
);
|
||||
Vector::<T, N>::vec(self.data[m])
|
||||
}
|
||||
|
||||
#[inline]
|
||||
pub fn set_row(&mut self, m: usize, val: &Vector<T, N>) {
|
||||
assert!(
|
||||
m < M,
|
||||
"Row index {} out of bounds for {}x{} matrix",
|
||||
m,
|
||||
M,
|
||||
N
|
||||
);
|
||||
for n in 0..N {
|
||||
self.data[m][n] = val.data[n][0];
|
||||
}
|
||||
}
|
||||
|
||||
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);
|
||||
}
|
||||
|
||||
#[inline]
|
||||
#[must_use]
|
||||
pub fn col(&self, n: usize) -> Vector<T, M> {
|
||||
assert!(
|
||||
n < N,
|
||||
"Column index {} out of bounds for {}x{} matrix",
|
||||
n,
|
||||
M,
|
||||
N
|
||||
);
|
||||
Vector::<T, M>::vec(self.data.map(|r| r[n]))
|
||||
}
|
||||
|
||||
#[inline]
|
||||
pub fn set_col(&mut self, n: usize, val: &Vector<T, M>) {
|
||||
assert!(
|
||||
n < N,
|
||||
"Column index {} out of bounds for {}x{} matrix",
|
||||
n,
|
||||
M,
|
||||
N
|
||||
);
|
||||
|
||||
for m in 0..M {
|
||||
self.data[m][n] = val.data[m][0];
|
||||
}
|
||||
}
|
||||
|
||||
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);
|
||||
}
|
||||
|
||||
#[must_use]
|
||||
pub fn rows<'a>(&'a self) -> impl Iterator<Item = Vector<T, N>> + 'a {
|
||||
(0..M).map(|m| self.row(m))
|
||||
}
|
||||
|
||||
#[must_use]
|
||||
pub fn cols<'a>(&'a self) -> impl Iterator<Item = Vector<T, M>> + 'a {
|
||||
(0..N).map(|n| self.col(n))
|
||||
}
|
||||
|
||||
#[must_use]
|
||||
pub fn permute_rows(&self, ms: &Vector<usize, M>) -> Self
|
||||
where
|
||||
T: Default,
|
||||
{
|
||||
Self::from_rows(ms.elements().map(|&m| self.row(m)))
|
||||
}
|
||||
|
||||
#[must_use]
|
||||
pub fn permute_cols(&self, ns: &Vector<usize, N>) -> Self
|
||||
where
|
||||
T: Default,
|
||||
{
|
||||
Self::from_cols(ns.elements().map(|&n| self.col(n)))
|
||||
}
|
||||
|
||||
pub fn transpose(&self) -> Matrix<T, N, M>
|
||||
where
|
||||
Matrix<T, N, M>: Default,
|
||||
{
|
||||
Matrix::<T, N, M>::from_rows(self.cols())
|
||||
}
|
||||
}
|
||||
|
||||
// Index
|
||||
impl<I, T, const M: usize, const N: usize> Index<I> for Matrix<T, M, N>
|
||||
where
|
||||
I: Index2D,
|
||||
T: Copy,
|
||||
{
|
||||
type Output = T;
|
||||
|
||||
#[inline(always)]
|
||||
fn index(&self, index: I) -> &Self::Output {
|
||||
self.get(index).expect(&*format!(
|
||||
"index {:?} out of range for {}x{} Matrix",
|
||||
index, M, N
|
||||
))
|
||||
}
|
||||
}
|
||||
|
||||
// IndexMut
|
||||
impl<I, T, const M: usize, const N: usize> IndexMut<I> for Matrix<T, M, N>
|
||||
where
|
||||
I: Index2D,
|
||||
T: Copy,
|
||||
{
|
||||
#[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, M, N
|
||||
))
|
||||
}
|
||||
}
|
||||
|
||||
// CONVERSIONS
|
||||
|
||||
// Convert from 2D Array (equivalent to new)
|
||||
impl<T: Copy, const M: usize, const N: usize> From<[[T; N]; M]> for Matrix<T, M, N> {
|
||||
fn from(data: [[T; N]; M]) -> Self {
|
||||
Self::mat(data)
|
||||
}
|
||||
}
|
||||
|
||||
// Convert from 1D Array (equivalent to vec)
|
||||
impl<T: Copy, const M: usize> From<[T; M]> for Vector<T, M> {
|
||||
fn from(data: [T; M]) -> Self {
|
||||
Self::vec(data)
|
||||
}
|
||||
}
|
||||
|
||||
// Convert from scalar (equivalent to fill)
|
||||
impl<T: Copy, const M: usize, const N: usize> From<T> for Matrix<T, M, N> {
|
||||
fn from(scalar: T) -> Self {
|
||||
Self::fill(scalar)
|
||||
}
|
||||
}
|
||||
|
||||
// IntoIter
|
||||
impl<T: Copy, const M: usize, const N: usize> IntoIterator for Matrix<T, M, N> {
|
||||
type Item = T;
|
||||
type IntoIter = Flatten<std::array::IntoIter<[T; N], M>>;
|
||||
|
||||
fn into_iter(self) -> Self::IntoIter {
|
||||
self.data.into_iter().flatten()
|
||||
}
|
||||
}
|
||||
|
||||
// FromIterator
|
||||
impl<T: Copy, const M: usize, const N: usize> FromIterator<T> for Matrix<T, M, N>
|
||||
where
|
||||
Self: Default,
|
||||
{
|
||||
fn from_iter<I: IntoIterator<Item = T>>(iter: I) -> Self {
|
||||
let mut result: Self = Default::default();
|
||||
for (l, r) in zip(result.elements_mut(), iter) {
|
||||
*l = r;
|
||||
}
|
||||
result
|
||||
}
|
||||
}
|
||||
|
||||
@ -0,0 +1,182 @@
|
||||
use crate::{Matrix, Vector};
|
||||
use num_traits::real::Real;
|
||||
use num_traits::{Inv, Num, NumOps, One, Pow, Signed, Zero};
|
||||
use std::iter::{zip, Product, Sum};
|
||||
use std::ops::{Add, Mul};
|
||||
|
||||
/// Operations for column vectors
|
||||
impl<T: Copy, const N: usize> Vector<T, N> {
|
||||
/// Compute the dot product of two vectors, otherwise known as the scalar product.
|
||||
/// This is the sum of the elementwise product, or in math terms
|
||||
///
|
||||
/// $vec(a) * vec(b) = sum_(i=1)^n a_i b_i = a_1 b_1 + a_2 b_2 + ... + a_n b_n$
|
||||
///
|
||||
/// for example, $[[1],[2],[3]] * [[4],[5],[6]] = (1 * 4) + (2 * 5) + (3 * 6) = 32$
|
||||
///
|
||||
/// For vectors in euclidean space, this has the property that it is equal to the magnitudes of
|
||||
/// the vectors times the cosine of the angle between them.
|
||||
///
|
||||
/// $vec(a) * vec(b) = |vec(a)| |vec(b)| cos(theta)$
|
||||
///
|
||||
/// this also gives it the special property that the dot product of a vector and itself is the
|
||||
/// square of its magnitude. You may recognize the 2D version as the
|
||||
/// [pythagorean theorem](https://en.wikipedia.org/wiki/Pythagorean_theorem).
|
||||
///
|
||||
/// see [dot product](https://en.wikipedia.org/wiki/Dot_product) on Wikipedia for more
|
||||
/// information.
|
||||
pub fn dot<R>(&self, rhs: &R) -> T
|
||||
where
|
||||
for<'s> &'s Self: Mul<&'s R, Output = Self>,
|
||||
T: Sum<T>,
|
||||
{
|
||||
(self * rhs).elements().cloned().sum()
|
||||
}
|
||||
|
||||
pub fn sqrmag(&self) -> T
|
||||
where
|
||||
for<'s> &'s Self: Mul<&'s Self, Output = Self>,
|
||||
T: Sum<T>,
|
||||
{
|
||||
self.dot(self)
|
||||
}
|
||||
|
||||
pub fn mag(&self) -> T
|
||||
where
|
||||
T: Sum<T> + Mul<T> + Real,
|
||||
{
|
||||
self.sqrmag().sqrt()
|
||||
}
|
||||
|
||||
pub fn normalized(&self) -> Option<Self>
|
||||
where
|
||||
T: Sum<T> + Mul<T> + Real,
|
||||
{
|
||||
match self.mag() {
|
||||
mag if mag.abs() < T::epsilon() => None,
|
||||
mag => Some(self / mag),
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// Cross product operations for column vectors in $RR^3$
|
||||
impl<T: Copy> Vector<T, 3> {
|
||||
pub fn cross_r<R: Copy>(&self, rhs: &Vector<R, 3>) -> Self
|
||||
where
|
||||
T: NumOps<R> + NumOps,
|
||||
{
|
||||
Self::vec([
|
||||
(self[1] * rhs[2]) - (self[2] * rhs[1]),
|
||||
(self[2] * rhs[0]) - (self[0] * rhs[2]),
|
||||
(self[0] * rhs[1]) - (self[1] * rhs[0]),
|
||||
])
|
||||
}
|
||||
|
||||
pub fn cross_l<R: Copy>(&self, rhs: &Vector<R, 3>) -> Vector<R, 3>
|
||||
where
|
||||
R: NumOps<T> + NumOps,
|
||||
{
|
||||
rhs.cross_r(self)
|
||||
}
|
||||
}
|
||||
|
||||
/// Operations for Matrices
|
||||
impl<T: Copy, const M: usize, const N: usize> Matrix<T, M, N> {
|
||||
pub fn mmul<R: Copy, const P: usize>(&self, rhs: &Matrix<R, N, P>) -> Matrix<T, M, P>
|
||||
where
|
||||
T: Default + NumOps<R> + Sum,
|
||||
{
|
||||
let mut result: Matrix<T, M, P> = Default::default();
|
||||
|
||||
for (m, a) in self.rows().enumerate() {
|
||||
for (n, b) in rhs.cols().enumerate() {
|
||||
result[(m, n)] = a.dot(&b)
|
||||
}
|
||||
}
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
/// Computes the absolute value of each element of the matrix
|
||||
pub fn abs(&self) -> Self
|
||||
where
|
||||
T: Signed + Default,
|
||||
{
|
||||
self.elements().map(|&x| x.abs()).collect()
|
||||
}
|
||||
|
||||
/// Computes the sign of each element of the matrix
|
||||
pub fn signum(&self) -> Self
|
||||
where
|
||||
T: Signed + Default,
|
||||
{
|
||||
self.elements().map(|&x| x.signum()).collect()
|
||||
}
|
||||
|
||||
/// Raises every element to the power of rhs, where rhs is either a scalar or a matrix of exponents
|
||||
pub fn pow<R, O>(self, rhs: R) -> O
|
||||
where
|
||||
Self: Pow<R, Output = O>,
|
||||
{
|
||||
Pow::pow(self, rhs)
|
||||
}
|
||||
}
|
||||
|
||||
// Sum up matrices
|
||||
impl<T: Copy, const M: usize, const N: usize> Sum for Matrix<T, M, N>
|
||||
where
|
||||
Self: Zero + Add<Output = Self>,
|
||||
{
|
||||
fn sum<I: Iterator<Item = Self>>(iter: I) -> Self {
|
||||
iter.fold(Self::zero(), Self::add)
|
||||
}
|
||||
}
|
||||
|
||||
// Product of matrices
|
||||
impl<T: Copy, const M: usize, const N: usize> Product for Matrix<T, M, N>
|
||||
where
|
||||
Self: One + Mul<Output = Self>,
|
||||
{
|
||||
fn product<I: Iterator<Item = Self>>(iter: I) -> Self {
|
||||
iter.fold(Self::one(), Self::mul)
|
||||
}
|
||||
}
|
||||
|
||||
/// Inverse trait. Note that this is the elementwise inverse, not the matrix inverse.
|
||||
/// For the inverse matrix see [`LUDecomposable::inv()`](crate::decompose::LUDecompose::inv())
|
||||
impl<T: Copy + Inv<Output = T> + Default, const M: usize, const N: usize> Inv for Matrix<T, M, N> {
|
||||
type Output = Self;
|
||||
|
||||
fn inv(self) -> Self::Output {
|
||||
self.elements().map(|t| t.inv()).collect()
|
||||
}
|
||||
}
|
||||
|
||||
/// Pow for $Matrix^{scalar}$
|
||||
impl<T, R, O, const M: usize, const N: usize> Pow<R> for Matrix<T, M, N>
|
||||
where
|
||||
T: Copy + Pow<R, Output = O>,
|
||||
R: Copy + Num,
|
||||
O: Copy + Default,
|
||||
{
|
||||
type Output = Matrix<O, M, N>;
|
||||
|
||||
fn pow(self, rhs: R) -> Self::Output {
|
||||
self.elements().map(|&x| x.pow(rhs)).collect()
|
||||
}
|
||||
}
|
||||
|
||||
/// Pow for $Matrix^{Matrix}$
|
||||
impl<T, R, O, const M: usize, const N: usize> Pow<Matrix<R, M, N>> for Matrix<T, M, N>
|
||||
where
|
||||
T: Copy + Pow<R, Output = O>,
|
||||
R: Copy,
|
||||
O: Copy + Default,
|
||||
{
|
||||
type Output = Matrix<O, M, N>;
|
||||
|
||||
fn pow(self, rhs: Matrix<R, M, N>) -> Self::Output {
|
||||
zip(self.elements(), rhs.elements())
|
||||
.map(|(x, &r)| x.pow(r))
|
||||
.collect()
|
||||
}
|
||||
}
|
||||
@ -1,533 +0,0 @@
|
||||
use index::Index2D;
|
||||
|
||||
use num_traits::{NumOps, One, Zero};
|
||||
use std::fmt::Debug;
|
||||
use std::iter::{zip, Flatten, Product, Sum};
|
||||
|
||||
use std::ops::{Add, AddAssign, Deref, DerefMut, Index, IndexMut, Mul, MulAssign, Neg};
|
||||
|
||||
pub mod ops;
|
||||
mod index;
|
||||
|
||||
/// A 2D array of values which can be operated upon.
|
||||
///
|
||||
/// Matrices have a fixed size known at compile time
|
||||
#[derive(Debug, Copy, Clone, PartialEq)]
|
||||
pub struct Matrix<T, const M: usize, const N: usize>
|
||||
where
|
||||
T: Copy,
|
||||
{
|
||||
data: [[T; N]; M], // Row-Major order
|
||||
}
|
||||
|
||||
/// An alias for a [Matrix] with a single column
|
||||
pub type Vector<T, const N: usize> = Matrix<T, N, 1>;
|
||||
|
||||
// Simple access functions that only require T be copyable
|
||||
impl<T: Copy, const M: usize, const N: usize> Matrix<T, M, N> {
|
||||
/// Generate a new matrix from a 2D Array
|
||||
///
|
||||
/// # Arguments
|
||||
///
|
||||
/// * `data`: A 2D array of elements to copy into the new matrix
|
||||
///
|
||||
/// returns: Matrix<T, M, N>
|
||||
///
|
||||
/// # Examples
|
||||
///
|
||||
/// ```
|
||||
/// # use vector_victor::Matrix;
|
||||
/// let a = Matrix::new([[1,2,3,4];4]);
|
||||
/// ```
|
||||
#[must_use]
|
||||
pub fn new(data: [[T; N]; M]) -> Self {
|
||||
assert!(M > 0, "Matrix must have at least 1 row");
|
||||
assert!(N > 0, "Matrix must have at least 1 column");
|
||||
Matrix::<T, M, N> { data }
|
||||
}
|
||||
|
||||
/// Generate a new matrix from a single scalar
|
||||
///
|
||||
/// # Arguments
|
||||
///
|
||||
/// * `scalar`: Scalar value to copy into the new matrix.
|
||||
///
|
||||
/// returns: Matrix<T, M, N>
|
||||
///
|
||||
/// # Examples
|
||||
///
|
||||
/// ```
|
||||
/// # use vector_victor::Matrix;
|
||||
/// let my_matrix = Matrix::<i32,4,4>::fill(5);
|
||||
/// // is equivalent to
|
||||
/// assert_eq!(my_matrix, Matrix::new([[5;4];4]))
|
||||
/// ```
|
||||
#[must_use]
|
||||
pub fn fill(scalar: T) -> Matrix<T, M, N> {
|
||||
assert!(M > 0, "Matrix must have at least 1 row");
|
||||
assert!(N > 0, "Matrix must have at least 1 column");
|
||||
Matrix::<T, M, N> {
|
||||
data: [[scalar; N]; M],
|
||||
}
|
||||
}
|
||||
|
||||
/// Create a matrix from an iterator of vectors
|
||||
///
|
||||
/// # Arguments
|
||||
///
|
||||
/// * `iter`: iterator of vectors to copy into rows
|
||||
///
|
||||
/// returns: Matrix<T, M, N>
|
||||
///
|
||||
/// # Examples
|
||||
///
|
||||
/// ```
|
||||
/// # use vector_victor::Matrix;
|
||||
/// let my_matrix = Matrix::new([[1,2,3],[4,5,6]]);
|
||||
/// let transpose : Matrix<_,3,2>= Matrix::from_rows(my_matrix.cols());
|
||||
/// assert_eq!(transpose, Matrix::new([[1,4],[2,5],[3,6]]))
|
||||
/// ```
|
||||
#[must_use]
|
||||
pub fn from_rows<I>(iter: I) -> Self
|
||||
where
|
||||
I: IntoIterator<Item = Vector<T, N>>,
|
||||
Self: Default,
|
||||
{
|
||||
let mut result = Self::default();
|
||||
for (m, row) in iter.into_iter().enumerate().take(M) {
|
||||
result.set_row(m, &row)
|
||||
}
|
||||
result
|
||||
}
|
||||
|
||||
/// Create a matrix from an iterator of vectors
|
||||
///
|
||||
/// # Arguments
|
||||
///
|
||||
/// * `iter`: iterator of vectors to copy into columns
|
||||
///
|
||||
/// returns: Matrix<T, M, N>
|
||||
///
|
||||
/// # Examples
|
||||
///
|
||||
/// ```
|
||||
/// # use vector_victor::Matrix;
|
||||
/// let my_matrix = Matrix::new([[1,2,3],[4,5,6]]);
|
||||
/// let transpose : Matrix<_,3,2>= Matrix::from_cols(my_matrix.rows());
|
||||
/// assert_eq!(transpose, Matrix::new([[1,4],[2,5],[3,6]]))
|
||||
/// ```
|
||||
#[must_use]
|
||||
pub fn from_cols<I>(iter: I) -> Self
|
||||
where
|
||||
I: IntoIterator<Item = Vector<T, M>>,
|
||||
Self: Default,
|
||||
{
|
||||
let mut result = Self::default();
|
||||
for (n, col) in iter.into_iter().enumerate().take(N) {
|
||||
result.set_col(n, &col)
|
||||
}
|
||||
result
|
||||
}
|
||||
|
||||
/// Returns an iterator over the elements of the matrix in row-major order.
|
||||
///
|
||||
/// # Examples
|
||||
/// ```
|
||||
/// # use vector_victor::Matrix;
|
||||
/// let my_matrix = Matrix::new([[1,2],[3,4]]);
|
||||
/// assert!(vec![1,2,3,4].iter().eq(my_matrix.elements()))
|
||||
/// ```
|
||||
#[must_use]
|
||||
pub fn elements<'a>(&'a self) -> impl Iterator<Item = &'a T> + 'a {
|
||||
self.data.iter().flatten()
|
||||
}
|
||||
|
||||
/// Returns a mutable iterator over the elements of the matrix in row-major order.
|
||||
#[must_use]
|
||||
pub fn elements_mut<'a>(&'a mut self) -> impl Iterator<Item = &'a mut T> + 'a {
|
||||
self.data.iter_mut().flatten()
|
||||
}
|
||||
|
||||
/// Returns a reference to the element at that position in the matrix, or `None` if out of bounds.
|
||||
///
|
||||
/// # Examples
|
||||
///
|
||||
/// ```
|
||||
/// # use vector_victor::Matrix;
|
||||
/// let my_matrix = Matrix::new([[1,2],[3,4]]);
|
||||
///
|
||||
/// // 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[],
|
||||
/// // but returns an Option instead of panicking
|
||||
/// assert_eq!(my_matrix.get(2), Some(&my_matrix[2]));
|
||||
///
|
||||
/// // index 4 is out of range, so get(4) returns None.
|
||||
/// assert_eq!(my_matrix.get(4), None);
|
||||
/// ```
|
||||
#[inline]
|
||||
#[must_use]
|
||||
pub fn get(&self, index: impl Index2D) -> Option<&T> {
|
||||
let (m, n) = index.to_2d(M, N)?;
|
||||
Some(&self.data[m][n])
|
||||
}
|
||||
|
||||
/// Returns a mutable reference to the element at that position in the matrix, or `None` if out of bounds.
|
||||
#[inline]
|
||||
#[must_use]
|
||||
pub fn get_mut(&mut self, index: impl Index2D) -> Option<&mut T> {
|
||||
let (m, n) = index.to_2d(M, N)?;
|
||||
Some(&mut self.data[m][n])
|
||||
}
|
||||
|
||||
/// Returns a row of the matrix. or [None] if index is out of bounds
|
||||
///
|
||||
/// # Examples
|
||||
///
|
||||
/// ```
|
||||
/// # use vector_victor::{Matrix, Vector};
|
||||
/// let my_matrix = Matrix::new([[1,2],[3,4]]);
|
||||
///
|
||||
/// // row at index 1
|
||||
/// assert_eq!(my_matrix.row(1), Vector::vec([3,4]));
|
||||
/// ```
|
||||
#[inline]
|
||||
#[must_use]
|
||||
pub fn row(&self, m: usize) -> Vector<T, N> {
|
||||
assert!(
|
||||
m < M,
|
||||
"Row index {} out of bounds for {}x{} matrix",
|
||||
m,
|
||||
M,
|
||||
N
|
||||
);
|
||||
Vector::<T, N>::vec(self.data[m])
|
||||
}
|
||||
|
||||
#[inline]
|
||||
pub fn set_row(&mut self, m: usize, val: &Vector<T, N>) {
|
||||
assert!(
|
||||
m < M,
|
||||
"Row index {} out of bounds for {}x{} matrix",
|
||||
m,
|
||||
M,
|
||||
N
|
||||
);
|
||||
for n in 0..N {
|
||||
self.data[m][n] = val.data[n][0];
|
||||
}
|
||||
}
|
||||
|
||||
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);
|
||||
}
|
||||
|
||||
#[inline]
|
||||
#[must_use]
|
||||
pub fn col(&self, n: usize) -> Vector<T, M> {
|
||||
assert!(
|
||||
n < N,
|
||||
"Column index {} out of bounds for {}x{} matrix",
|
||||
n,
|
||||
M,
|
||||
N
|
||||
);
|
||||
Vector::<T, M>::vec(self.data.map(|r| r[n]))
|
||||
}
|
||||
|
||||
#[inline]
|
||||
pub fn set_col(&mut self, n: usize, val: &Vector<T, M>) {
|
||||
assert!(
|
||||
n < N,
|
||||
"Column index {} out of bounds for {}x{} matrix",
|
||||
n,
|
||||
M,
|
||||
N
|
||||
);
|
||||
|
||||
for m in 0..M {
|
||||
self.data[m][n] = val.data[m][0];
|
||||
}
|
||||
}
|
||||
|
||||
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);
|
||||
}
|
||||
|
||||
#[must_use]
|
||||
pub fn rows<'a>(&'a self) -> impl Iterator<Item = Vector<T, N>> + 'a {
|
||||
(0..M).map(|m| self.row(m))
|
||||
}
|
||||
|
||||
#[must_use]
|
||||
pub fn cols<'a>(&'a self) -> impl Iterator<Item = Vector<T, M>> + 'a {
|
||||
(0..N).map(|n| self.col(n))
|
||||
}
|
||||
|
||||
#[must_use]
|
||||
pub fn permute_rows(&self, ms: &Vector<usize, M>) -> Self
|
||||
where
|
||||
T: Default,
|
||||
{
|
||||
Self::from_rows(ms.elements().map(|&m| self.row(m)))
|
||||
}
|
||||
|
||||
#[must_use]
|
||||
pub fn permute_cols(&self, ns: &Vector<usize, N>) -> Self
|
||||
where
|
||||
T: Default,
|
||||
{
|
||||
Self::from_cols(ns.elements().map(|&n| self.col(n)))
|
||||
}
|
||||
|
||||
pub fn transpose(&self) -> Matrix<T, N, M>
|
||||
where
|
||||
Matrix<T, N, M>: Default,
|
||||
{
|
||||
Matrix::<T, N, M>::from_rows(self.cols())
|
||||
}
|
||||
|
||||
pub fn abs(&self) -> Self
|
||||
where
|
||||
T: Default + PartialOrd + Zero + Neg<Output = T>,
|
||||
{
|
||||
self.elements()
|
||||
.map(|&x| match x > T::zero() {
|
||||
true => x,
|
||||
false => -x,
|
||||
})
|
||||
.collect()
|
||||
}
|
||||
}
|
||||
|
||||
// 1D vector implementations
|
||||
impl<T: Copy, const N: usize> Vector<T, N> {
|
||||
/// Create a vector from a 1D array.
|
||||
/// Note that vectors are always column vectors unless explicitly instantiated as row vectors
|
||||
///
|
||||
/// # Examples
|
||||
/// ```
|
||||
/// # use vector_victor::{Matrix, Vector};
|
||||
/// let my_vector = Vector::vec([1,2,3,4]);
|
||||
/// // is equivalent to
|
||||
/// assert_eq!(my_vector, Matrix::new([[1],[2],[3],[4]]));
|
||||
/// ```
|
||||
pub fn vec(data: [T; N]) -> Self {
|
||||
assert!(N > 0, "Vector must have at least 1 element");
|
||||
return Vector::<T, N> {
|
||||
data: data.map(|e| [e]),
|
||||
};
|
||||
}
|
||||
|
||||
pub fn dot<R>(&self, rhs: &R) -> T
|
||||
where
|
||||
for<'s> &'s Self: Mul<&'s R, Output = Self>,
|
||||
T: Sum<T>,
|
||||
{
|
||||
(self * rhs).elements().cloned().sum()
|
||||
}
|
||||
}
|
||||
|
||||
// Cross Product
|
||||
impl<T: Copy> Vector<T, 3> {
|
||||
pub fn cross_r<R: Copy>(&self, rhs: &Vector<R, 3>) -> Self
|
||||
where
|
||||
T: NumOps<R> + NumOps,
|
||||
{
|
||||
Self::vec([
|
||||
(self[1] * rhs[2]) - (self[2] * rhs[1]),
|
||||
(self[2] * rhs[0]) - (self[0] * rhs[2]),
|
||||
(self[0] * rhs[1]) - (self[1] * rhs[0]),
|
||||
])
|
||||
}
|
||||
|
||||
pub fn cross_l<R: Copy>(&self, rhs: &Vector<R, 3>) -> Vector<R, 3>
|
||||
where
|
||||
R: NumOps<T> + NumOps,
|
||||
{
|
||||
rhs.cross_r(self)
|
||||
}
|
||||
}
|
||||
|
||||
//Matrix Multiplication
|
||||
impl<T: Copy, const M: usize, const N: usize> Matrix<T, M, N> {
|
||||
pub fn mmul<R: Copy, const P: usize>(&self, rhs: &Matrix<R, N, P>) -> Matrix<T, M, P>
|
||||
where
|
||||
T: Default + NumOps<R> + Sum,
|
||||
{
|
||||
let mut result: Matrix<T, M, P> = Default::default();
|
||||
|
||||
for (m, a) in self.rows().enumerate() {
|
||||
for (n, b) in rhs.cols().enumerate() {
|
||||
result[(m, n)] = a.dot(&b)
|
||||
}
|
||||
}
|
||||
|
||||
return result;
|
||||
}
|
||||
}
|
||||
|
||||
// Square matrix implementations
|
||||
impl<T: Copy, const N: usize> Matrix<T, N, N> {
|
||||
/// Create an identity matrix
|
||||
#[must_use]
|
||||
pub fn identity() -> Self
|
||||
where
|
||||
T: Zero + One,
|
||||
{
|
||||
let mut result = Self::zero();
|
||||
for i in 0..N {
|
||||
result[(i, i)] = T::one();
|
||||
}
|
||||
return result;
|
||||
}
|
||||
|
||||
/// returns an iterator over the elements along the diagonal of a square matrix
|
||||
#[must_use]
|
||||
pub fn diagonals<'s>(&'s self) -> impl Iterator<Item = T> + 's {
|
||||
(0..N).map(|n| self[(n, n)])
|
||||
}
|
||||
|
||||
/// Returns an iterator over the elements directly below the diagonal of a square matrix
|
||||
#[must_use]
|
||||
pub fn subdiagonals<'s>(&'s self) -> impl Iterator<Item = T> + 's {
|
||||
(0..N - 1).map(|n| self[(n, n + 1)])
|
||||
}
|
||||
}
|
||||
|
||||
// Index
|
||||
impl<I, T, const M: usize, const N: usize> Index<I> for Matrix<T, M, N>
|
||||
where
|
||||
I: Index2D,
|
||||
T: Copy,
|
||||
{
|
||||
type Output = T;
|
||||
|
||||
#[inline(always)]
|
||||
fn index(&self, index: I) -> &Self::Output {
|
||||
self.get(index).expect(&*format!(
|
||||
"index {:?} out of range for {}x{} Matrix",
|
||||
index, M, N
|
||||
))
|
||||
}
|
||||
}
|
||||
|
||||
// IndexMut
|
||||
impl<I, T, const M: usize, const N: usize> IndexMut<I> for Matrix<T, M, N>
|
||||
where
|
||||
I: Index2D,
|
||||
T: Copy,
|
||||
{
|
||||
#[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, M, N
|
||||
))
|
||||
}
|
||||
}
|
||||
// Default
|
||||
impl<T: Copy + Default, const M: usize, const N: usize> Default for Matrix<T, M, N> {
|
||||
fn default() -> Self {
|
||||
Matrix::fill(T::default())
|
||||
}
|
||||
}
|
||||
|
||||
// Zero
|
||||
impl<T: Copy + Zero, const M: usize, const N: usize> Zero for Matrix<T, M, N> {
|
||||
fn zero() -> Self {
|
||||
Matrix::fill(T::zero())
|
||||
}
|
||||
|
||||
fn is_zero(&self) -> bool {
|
||||
self.elements().all(|e| e.is_zero())
|
||||
}
|
||||
}
|
||||
|
||||
// One
|
||||
impl<T: Copy + One, const M: usize, const N: usize> One for Matrix<T, M, N> {
|
||||
fn one() -> Self {
|
||||
Matrix::fill(T::one())
|
||||
}
|
||||
}
|
||||
|
||||
impl<T: Copy, const M: usize, const N: usize> From<[[T; N]; M]> for Matrix<T, M, N> {
|
||||
fn from(data: [[T; N]; M]) -> Self {
|
||||
Self::new(data)
|
||||
}
|
||||
}
|
||||
|
||||
impl<T: Copy, const M: usize> From<[T; M]> for Vector<T, M> {
|
||||
fn from(data: [T; M]) -> Self {
|
||||
Self::vec(data)
|
||||
}
|
||||
}
|
||||
|
||||
impl<T: Copy, const M: usize, const N: usize> From<T> for Matrix<T, M, N> {
|
||||
fn from(scalar: T) -> Self {
|
||||
Self::fill(scalar)
|
||||
}
|
||||
}
|
||||
|
||||
// deref 1x1 matrices to a scalar automatically
|
||||
impl<T: Copy> Deref for Matrix<T, 1, 1> {
|
||||
type Target = T;
|
||||
|
||||
fn deref(&self) -> &Self::Target {
|
||||
&self.data[0][0]
|
||||
}
|
||||
}
|
||||
|
||||
// deref 1x1 matrices to a mutable scalar automatically
|
||||
impl<T: Copy> DerefMut for Matrix<T, 1, 1> {
|
||||
fn deref_mut(&mut self) -> &mut Self::Target {
|
||||
&mut self.data[0][0]
|
||||
}
|
||||
}
|
||||
|
||||
// IntoIter
|
||||
impl<T: Copy, const M: usize, const N: usize> IntoIterator for Matrix<T, M, N> {
|
||||
type Item = T;
|
||||
type IntoIter = Flatten<std::array::IntoIter<[T; N], M>>;
|
||||
|
||||
fn into_iter(self) -> Self::IntoIter {
|
||||
self.data.into_iter().flatten()
|
||||
}
|
||||
}
|
||||
|
||||
// FromIterator
|
||||
impl<T: Copy, const M: usize, const N: usize> FromIterator<T> for Matrix<T, M, N>
|
||||
where
|
||||
Self: Default,
|
||||
{
|
||||
fn from_iter<I: IntoIterator<Item = T>>(iter: I) -> Self {
|
||||
let mut result: Self = Default::default();
|
||||
for (l, r) in zip(result.elements_mut(), iter) {
|
||||
*l = r;
|
||||
}
|
||||
result
|
||||
}
|
||||
}
|
||||
|
||||
impl<T: Copy + AddAssign, const M: usize, const N: usize> Sum for Matrix<T, M, N>
|
||||
where
|
||||
Self: Zero + Add<Output = Self>,
|
||||
{
|
||||
fn sum<I: Iterator<Item = Self>>(iter: I) -> Self {
|
||||
iter.fold(Self::zero(), Self::add)
|
||||
}
|
||||
}
|
||||
|
||||
impl<T: Copy + MulAssign, const M: usize, const N: usize> Product for Matrix<T, M, N>
|
||||
where
|
||||
Self: One + Mul<Output = Self>,
|
||||
{
|
||||
fn product<I: Iterator<Item = Self>>(iter: I) -> Self {
|
||||
iter.fold(Self::one(), Self::mul)
|
||||
}
|
||||
}
|
||||
@ -0,0 +1 @@
|
||||
|
||||
@ -1,4 +1,4 @@
|
||||
use crate::matrix::Matrix;
|
||||
use crate::Matrix;
|
||||
use num_traits::Num;
|
||||
|
||||
// borrowed from the auto_ops crate
|
||||
Loading…
Reference in New Issue