pub mod error;
mod trait_impls;
pub mod ffi;
use core::slice;
use std::{fmt::Debug, fmt::Display};
use std::ops::{AddAssign, MulAssign, Neg, SubAssign};
use std::mem::swap;
use anyhow::{Error, Result};
use error::*;
use num_traits::Num;
pub use trait_impls::*;
pub trait Element<T>: Num + Copy + Debug + Display + AddAssign + MulAssign + SubAssign + Neg<Output = T> {}
impl Element<f32> for f32 {}
impl Element<f64> for f64 {}
impl Element<i8> for i8 {}
impl Element<i16> for i16 {}
impl Element<i32> for i32 {}
impl Element<i64> for i64 {}
#[derive(Clone, Debug, PartialEq)]
pub struct Matrix<T>
where T: Element<T>
{
rows: usize,
cols: usize,
vals: Vec<T>,
}
impl <T> Matrix<T>
where T: Element<T>
{
pub fn new(rows: usize, cols: usize) -> Matrix<T>
{
let mut a = Matrix
{
rows,
cols,
vals: Vec::with_capacity(rows * cols),
};
for _ in 0..a.vals.capacity()
{
a.vals.push(T::zero());
}
a
}
pub fn new_identity(n: usize) -> Matrix<T>
{
let mut a = Matrix::new(n, n);
for i in 0..n
{
a[(i, i)] = T::one();
}
a
}
pub fn from_vec(cols: usize, vec: Vec<T>) -> Result<Matrix<T>>
{
if vec.len() % cols != 0
{
return Err(MatrixFromVecError.into())
}
Ok(Matrix {
rows: vec.len() / cols,
cols,
vals: vec,
})
}
pub fn from_row_vec(vec: Vec<T>) -> Matrix<T>
{
Matrix
{
rows: 1,
cols: vec.len(),
vals: vec,
}
}
pub fn from_col_vec(vec: Vec<T>) -> Matrix<T>
{
Matrix
{
rows: vec.len(),
cols: 1,
vals: vec,
}
}
pub fn get_rows(&self) -> usize
{
self.rows
}
pub fn get_cols(&self) -> usize
{
self.cols
}
pub fn iter(&self) -> slice::Iter<'_, T>
{
self.vals.iter()
}
pub fn inplace_row_swap(&mut self, r1: usize, r2: usize)
{
let mut storage: T;
for i in 0..self.cols
{
storage = self[(i, r1)];
self[(i, r1)] = self[(i, r2)];
self[(i, r2)] = storage;
}
}
pub fn inplace_row_scale(&mut self, row: usize, scalar: T)
{
for i in 0..self.cols
{
self[(row, i)] *= scalar;
}
}
pub fn inplace_scale(&mut self, scalar: T)
{
for i in 0..self.rows
{
self.inplace_row_scale(i, scalar);
}
}
pub fn inplace_row_add(&mut self, r1: usize, r2: usize)
{
for i in 0..self.cols
{
let addend = self[(r2, i)];
self[(r1, i)] += addend;
}
}
pub fn inplace_scaled_row_add(&mut self, r1: usize, r2: usize, scalar: T)
{
for i in 0..self.cols
{
let addend = self[(r2, i)] * scalar;
self[(r1, i)] += addend;
}
}
pub fn multiply_matrix(&self, a: &Matrix<T>) -> Result<Matrix<T>>
{
if self.cols != a.rows
{
return Err(MatrixMultiplicationError.into())
}
let n = self.cols;
let mut result = Matrix::new(self.rows, a.cols);
for i in 0..self.rows
{
for j in 0..a.cols
{
for x in 0..n
{
result[(i, j)] += self[(i, x)] * a[(x, j)]
}
}
}
Ok(result)
}
pub fn augment_with(&self, a: &Matrix<T>) -> Result<Matrix<T>>
{
if a.get_rows() != self.rows
{
return Err(MatrixAugmentationError { a: self.rows, b: a.rows }.into())
}
let mut b: Matrix<T> = Matrix::new(self.rows, self.cols + a.cols);
for i in 0..b.rows
{
for j in 0..b.cols
{
if j < self.cols
{
b[(i, j)] = self[(i, j)];
}
else
{
b[(i, j)] = a[(i, j - self.cols)];
}
}
}
Ok(b)
}
pub fn subset(&self, r1: usize, c1: usize, r2: usize, c2: usize) -> Matrix<T>
{
let mut b = Matrix::new(r2 - r1 + 1, c2 - c1 + 1);
for i in r1..r2+1
{
for j in c1..c2+1
{
b[(i-r1, j-c1)] = self[(i, j)];
}
}
b
}
pub fn trace(&self) -> Result<T>
{
if self.rows != self.cols
{
return Err(NonSquareMatrixError.into())
}
let mut total: T = T::zero();
for i in 0..self.rows
{
total += self[(i, i)];
}
Ok(total)
}
pub fn transpose(&self) -> Matrix<T>
{
let mut tspose = self.clone();
swap(&mut tspose.rows, &mut tspose.cols);
for i in 0..self.rows
{
for j in 0..self.cols
{
tspose[(j, i)] = self[(i, j)];
}
}
tspose
}
fn try_inplace_invert_2(&mut self) -> Result<()>
{
let a11 = self[(0, 0)];
let a12 = self[(0, 1)];
let a21 = self[(1, 0)];
let a22 = self[(1, 1)];
let det = a11*a22 - a12*a21;
if det == T::zero()
{
return Err(MatrixInversionError::DeterminantWasZero.into())
}
self[(0, 0)] = a22 / det;
self[(1, 0)] = - a21 / det;
self[(0, 1)] = - a12 / det;
self[(1, 1)] = a11 / det;
Ok(())
}
fn try_inplace_invert_3(&mut self) -> Result<()>
{
let a11 = self[(0, 0)];
let a12 = self[(1, 0)];
let a13 = self[(2, 0)];
let a21 = self[(0, 1)];
let a22 = self[(1, 1)];
let a23 = self[(2, 1)];
let a31 = self[(0, 2)];
let a32 = self[(1, 2)];
let a33 = self[(2, 2)];
let det = a11*a22*a33 + a21*a32*a13 + a31*a12*a23
- a11*a32*a23 - a31*a22*a13 - a21*a12*a33;
if det == T::zero()
{
return Err(MatrixInversionError::DeterminantWasZero.into())
}
self[(0, 0)] = (a22 * a33 - a23 * a32) / det;
self[(1, 0)] = (a23 * a31 - a21 * a33) / det;
self[(2, 0)] = (a21 * a32 - a22 * a31) / det;
self[(0, 1)] = (a13 * a32 - a12 * a33) / det;
self[(1, 1)] = (a11 * a33 - a13 * a31) / det;
self[(2, 1)] = (a12 * a31 - a11 * a32) / det;
self[(0, 2)] = (a12 * a23 - a13 * a22) / det;
self[(1, 2)] = (a13 * a21 - a11 * a23) / det;
self[(2, 2)] = (a11 * a22 - a12 * a21) / det;
Ok(())
}
fn try_inplace_invert_4(&mut self) -> Result<()>
{
let a11 = self[(0, 0)];
let a12 = self[(1, 0)];
let a13 = self[(2, 0)];
let a14 = self[(3, 0)];
let a21 = self[(0, 1)];
let a22 = self[(1, 1)];
let a23 = self[(2, 1)];
let a24 = self[(3, 1)];
let a31 = self[(0, 2)];
let a32 = self[(1, 2)];
let a33 = self[(2, 2)];
let a34 = self[(3, 2)];
let a41 = self[(0, 3)];
let a42 = self[(1, 3)];
let a43 = self[(2, 3)];
let a44 = self[(3, 3)];
let det = a11*a22*a33*a44 + a11*a23*a34*a42 + a11*a24*a32*a43 +
a12*a21*a34*a43 + a12*a23*a31*a44 + a12*a24*a33*a41 +
a13*a21*a32*a44 + a13*a22*a34*a41 + a13*a24*a31*a42 +
a14*a21*a33*a42 + a14*a22*a34*a43 + a14*a23*a32*a41 -
a11*a22*a34*a43 - a11*a23*a32*a44 - a11*a24*a33*a42 -
a12*a21*a33*a44 - a12*a23*a34*a41 - a12*a24*a31*a43 -
a13*a21*a34*a42 - a13*a22*a31*a44 - a13*a24*a32*a41 -
a14*a21*a32*a43 - a14*a22*a33*a41 - a14*a23*a31*a42;
if det == T::zero()
{
return Err(MatrixInversionError::DeterminantWasZero.into());
}
self[(0, 0)] = (a22*a33*a44 + a23*a34*a42 + a24*a32*a43 - a22*a34*a43 - a23*a32*a44 - a24*a33*a42) / det;
self[(1, 0)] = (a12*a34*a43 + a13*a32*a44 + a14*a33*a42 - a12*a33*a44 - a13*a34*a42 - a14*a32*a43) / det;
self[(2, 0)] = (a12*a23*a44 + a13*a24*a42 + a14*a22*a43 - a12*a24*a43 - a13*a22*a44 - a14*a23*a42) / det;
self[(3, 0)] = (a12*a24*a33 + a13*a22*a34 + a14*a23*a32 - a12*a23*a34 - a13*a24*a32 - a14*a22*a33) / det;
self[(0, 1)] = (a21*a34*a43 + a23*a31*a44 + a24*a33*a41 - a21*a33*a44 - a23*a34*a41 - a24*a31*a43) / det;
self[(1, 1)] = (a11*a33*a44 + a13*a34*a41 + a14*a31*a43 - a11*a34*a43 - a13*a31*a44 - a14*a33*a41) / det;
self[(2, 1)] = (a11*a24*a43 + a13*a21*a44 + a14*a23*a41 - a11*a23*a44 - a13*a24*a41 - a14*a21*a43) / det;
self[(3, 1)] = (a11*a23*a34 + a13*a24*a31 + a14*a21*a33 - a11*a24*a33 - a13*a21*a34 - a14*a23*a31) / det;
self[(0, 2)] = (a21*a32*a44 + a22*a34*a41 + a24*a31*a42 - a21*a34*a42 - a22*a31*a44 - a24*a32*a41) / det;
self[(1, 2)] = (a11*a34*a42 + a12*a31*a44 + a14*a32*a41 - a11*a32*a44 - a12*a34*a41 - a14*a31*a42) / det;
self[(2, 2)] = (a11*a22*a44 + a12*a24*a41 + a14*a21*a42 - a11*a24*a42 - a12*a21*a44 - a14*a22*a41) / det;
self[(3, 2)] = (a11*a24*a32 + a12*a21*a34 + a14*a22*a31 - a11*a22*a34 - a12*a24*a31 - a14*a21*a32) / det;
self[(0, 3)] = (a21*a33*a42 + a22*a31*a43 + a23*a32*a41 - a21*a32*a43 - a22*a33*a41 - a23*a31*a42) / det;
self[(1, 3)] = (a11*a32*a43 + a12*a33*a41 + a13*a31*a42 - a11*a33*a42 - a12*a31*a43 - a13*a32*a41) / det;
self[(2, 3)] = (a11*a23*a42 + a12*a21*a43 + a13*a22*a41 - a11*a22*a43 - a12*a23*a41 - a13*a21*a42) / det;
self[(3, 3)] = (a11*a22*a33 + a12*a23*a31 + a13*a21*a32 - a11*a23*a32 - a12*a21*a33 - a13*a22*a31) / det;
Ok(())
}
fn try_inplace_invert_n(&mut self) -> Result<()>
{
let n = self.rows;
let mut inv: Matrix<T> = Matrix::new_identity(n);
for j in 0..n
{
for i in 0..n
{
if i == j
{
continue;
}
else
{
if self[(j, j)] == T::zero()
{
return Err(MatrixInversionError::ZeroDuringInversion.into())
}
let scalar = self[(i, j)] / self[(j, j)];
self.inplace_scaled_row_add(i, j, -scalar);
inv.inplace_scaled_row_add(i, j, -scalar);
}
}
}
for i in 0..n
{
let scalar: T = T::one() / self[(i, i)];
self.inplace_row_scale(i, scalar);
inv.inplace_row_scale(i, scalar);
}
*self = inv;
Ok(())
}
pub fn try_inplace_invert(&mut self) -> Result<()>
{
if self.rows != self.cols
{
return Err(NonSquareMatrixError.into())
}
if self.rows == 1 && self.vals[0] == T::zero()
{
return Err(Error::new(MatrixInversionError::SingularValueWasZero))
}
match self.rows {
1 => self.vals[0] = T::one() / self.vals[0],
2 => self.try_inplace_invert_2()?,
3 => self.try_inplace_invert_3()?,
4 => self.try_inplace_invert_4()?,
_ => self.try_inplace_invert_n()?,
};
Ok(())
}
}
#[macro_export]
macro_rules! row_vec {
($($e:expr),+ $(,)?) => {
Matrix::from_row_vec(
vec![$($e),+]
)
};
}
#[macro_export]
macro_rules! col_vec {
($($e:expr),+ $(,)?) => {
Matrix::from_col_vec(
vec![$($e),+]
)
};
}