use ::af;
use super::{Matrix, Axes};
use ::vector::Vector;
use ::norm::Norm;
use ::std::vec::{IntoIter};
use ::std::slice::{Iter};
impl Matrix {
pub fn new<T>(rows: usize, cols: usize, data: T) -> Matrix where T: Into<Vec<f64>> {
let data: Vec<f64> = data.into();
assert!(rows > 0 && cols > 0, "Cannot create Matrix with dimension of zero: {} x {}", rows, cols);
assert!(data.len() == rows * cols, "The amount of data does not match the dimensions. Length of data: {} != rows: {} * cols: {}",
data.len(), rows, cols);
let backend = af::Array::new(&data, af::Dim4::new(&[rows as u64, cols as u64, 1, 1]));
Matrix {rows, cols, backend, data}
}
pub fn from_fn<F>(rows: usize, cols: usize, f: F) -> Matrix where F: Fn(usize, usize) -> f64 {
let mut data: Vec<f64> = Vec::with_capacity(rows * cols);
for row in 0..rows {
for col in 0..cols {
data.push(f(rows, cols));
}
}
Matrix::new(rows, cols, data)
}
pub fn zeros(rows: usize, cols: usize) -> Matrix {
Matrix::new(rows, cols, vec![0f64; rows * cols])
}
pub fn ones(rows: usize, cols: usize) -> Matrix {
Matrix::new(rows, cols, vec![1f64; rows * cols])
}
pub fn from_diag<T>(diag: T) -> Matrix where T: Into<Vec<f64>> {
let diag: Vec<f64> = diag.into();
let size = diag.len();
let mut data: Vec<f64> = vec![0f64; size * size];
for (i, item) in diag.into_iter().enumerate().take(size) {
data[i * (size + 1)] = item;
}
Matrix::new(size, size, data)
}
pub fn identity(size: usize) -> Matrix {
let mut data: Vec<f64> = vec![0f64; size * size];
for i in 0..size {
data[i * (size + 1)] = 1f64;
}
Matrix::new(size, size, data)
}
pub fn rows(&self) -> usize {
self.rows
}
pub fn cols(&self) -> usize {
self.cols
}
pub fn size(&self) -> usize {self.rows() * self.cols()}
pub fn is_empty(&self) -> bool {af::Array::is_empty(&self.backend)}
pub fn data(&self) -> Vec<f64> {
self.data.clone()
}
pub fn into_vec(self) -> Vec<f64> {
self.data
}
pub fn get(&self, row: usize, col: usize) -> Option<f64> {
self.data.get(row * self.cols() + col).cloned()
}
pub fn get_row(&self, row_index: usize) -> Matrix {
assert!(row_index < self.rows(), "Tried to get a row index that does not exist: {} >= {}", row_index, self.rows());
let offset = row_index * self.cols();
let length = self.cols();
let data = self.data[offset..(length+offset)].to_vec();
Matrix::new(1, length, data)
}
pub fn get_col(&self, col_index: usize) -> Matrix {
assert!(col_index < self.cols(), "Tried to get a column index that does not exist: {} >= {}", col_index, self.cols());
let transpose = self.transpose();
transpose.get_row(col_index)
}
pub fn transpose(&self) -> Matrix {
Matrix::from(af::transpose(&self.backend, false))
}
pub fn apply<F>(&self, f: F) -> Matrix where F: Fn(f64) -> f64 {
let mut data: Vec<f64> = Vec::with_capacity(self.size());
for val in self {
data.push(f(val));
}
Matrix::new(self.rows(), self.cols(), data)
}
pub fn apply_self<F>(&mut self, f: F) where F: Fn(f64) -> f64 {
let mut data = self.data();
for val in &mut data {
*val = f(*val);
}
*self = Matrix::new(self.rows(), self.cols(), data);
}
pub fn max_all(&self) -> f64 { af::max_all(&self.backend).0 }
pub fn max(&self, axis: Axes) -> Vector {
match axis {
Axes::Col => Vector::from(af::max(&self.backend, 0)),
Axes::Row => Vector::from(af::max(&self.backend, 1)),
}
}
pub fn min_all(&self) -> f64 { af::min_all(&self.backend).0 }
pub fn min(&self, axis: Axes) -> Vector {
match axis {
Axes::Col => Vector::from(af::min(&self.backend, 0)),
Axes::Row => Vector::from(af::min(&self.backend, 1)),
}
}
pub fn select(&self, new_rows: usize, new_cols: usize, indxs: &[(usize, usize)]) -> Result<Matrix, String> {
let mut data: Vec<f64> = Vec::with_capacity(new_rows * new_cols);
for &(row, col) in indxs {
let o = self.get(row, col);
match o {
Some(x) => data.push(x),
None => return Err(format!("Tried to access index: {} x {} but dimensions are: {} x {}", row, col, self.rows(), self.cols()))
}
}
Ok(Matrix::new(new_rows, new_cols, data))
}
pub fn sum_row_all(&self) -> Vector {
Vector::from(af::sum(&self.backend, 1))
}
pub fn sum_col_all(&self) -> Vector {
Vector::from(af::sum(&self.backend, 0))
}
pub fn sum_row(&self, row_index: usize) -> f64 {
af::sum_all(&self.get_row(row_index).backend).0
}
pub fn sum_col(&self, col_index: usize) -> f64 {
af::sum_all(&self.get_col(col_index).backend).0
}
pub fn mean(&self, axis: Axes) -> Vector {
match axis {
Axes::Row => Vector::from(af::mean(&self.backend, 1)),
Axes::Col => Vector::from(af::mean(&self.backend, 0)),
}
}
pub fn mean_all(&self) -> f64 {
af::mean_all(&self.backend).0
}
pub fn variance(&self, axis: Axes) -> Result<Vector, String> {
let mean = self.mean(axis);
let n: usize;
let m: usize;
match axis {
Axes::Row => {
n = self.rows();
m = self.cols();
},
Axes::Col => {
n = self.cols();
m = self.rows();
}
}
if n < 2 {
return Err(format!("There must be at least two working rows or columns"));
}
let mut variance = Vector::zeros(m);
for i in 0..n {
let mut t: Vec<f64> = Vec::with_capacity(m);
unsafe { t.set_len(m) };
for j in 0..m {
t[j] = match axis {
Axes::Row => self.data.get(i * m + j).cloned().unwrap(),
Axes::Col => self.data.get(j * n + i).cloned().unwrap()
}
}
let v = Vector::new(t);
variance += (&v - &mean) * (&v - &mean);
}
let var_size = n - 1;
Ok(variance / (var_size as f64))
}
pub fn row_iter(&self) -> IntoIter<Vector> {
let mut rows: Vec<Vector> = Vec::with_capacity(self.rows());
for r in 0..self.rows() {
rows.push(Vector::new(self.get_row(r).into_vec()));
}
rows.into_iter()
}
pub fn col_iter(&self) -> IntoIter<Vector> {
let mut cols: Vec<Vector> = Vec::with_capacity(self.cols());
for c in 0..self.cols() {
cols.push(Vector::new(self.get_col(c).into_vec()));
}
cols.into_iter()
}
pub fn is_diag(&self) -> bool {
if self.rows() != self.cols() {
false
} else {
let mut next_diag: usize = 0;
self.data().iter().enumerate().all(|(i, data)| if i == next_diag {
next_diag += self.cols() + 1;
true
} else {
*data == 0.0
})
}
}
pub fn get_diag(&self) -> Vec<f64> {
let cols = self.cols();
let rows = self.rows();
let size = self.size();
let minor_axis = if cols < rows {
cols
} else {
rows
};
let mut diag: Vec<f64> = Vec::with_capacity(size);
for i in 0..minor_axis {
diag.push(self.get(i, i).unwrap());
}
diag
}
pub fn diag_iter(&self) -> IntoIter<f64> {
self.get_diag().into_iter()
}
pub fn norm(&self, norm_type: Norm) -> f64 {
match norm_type {
Norm::Taxicab => af::norm(&self.backend, af::NormType::MATRIX_1, 0.0, 0.0),
Norm::Euclidean => {
let mut norm = 0f64;
for row in self.row_iter() {
norm += row.dot(&row).unwrap();
}
norm.sqrt()
},
Norm::Infinity => af::norm(&self.backend, af::NormType::MATRIX_INF, 0.0, 0.0),
Norm::LP(p) => af::norm(&self.backend, af::NormType::MATRIX_L_PQ, p, 0.0),
}
}
pub fn metric(&self, rhs: &Matrix) -> Result<f64, String> {
if let Err(s) = Matrix::check_size(self, rhs) {
return Err(s);
}
let diff: Matrix = Matrix::from(&self.backend - &rhs.backend);
Ok(diff.norm(Norm::Euclidean))
}
pub fn select_cols(&self, indxs: &[usize]) -> Matrix {
let mut cols: Vec<Matrix> = Vec::with_capacity(indxs.len());
for c in indxs {
cols.push(self.get_col(*c));
}
Matrix::hcat_many(&cols)
}
pub fn select_rows(&self, indxs: &[usize]) -> Matrix {
let mut rows: Vec<Matrix> = Vec::with_capacity(indxs.len());
for r in indxs {
rows.push(self.get_row(*r));
}
Matrix::vcat_many(&rows)
}
pub fn select_block(&self, rows: &[usize], cols: &[usize]) -> Result<Matrix, String> {
let mut data: Vec<f64> = Vec::with_capacity(rows.len() * cols.len());
for r in rows {
for c in cols {
let val = self.get(*r, *c);
match val {
Some(v) => data.push(v),
None => return Err(format!("Row or Column index exceeded maximum. Matrix size is: {} x {}, Selected Rows are: {:?}, and Selected Cols are: {:?}",
self.rows(), self.cols(), rows, cols)),
}
}
}
Ok(Matrix::new(rows.len(), cols.len(), data))
}
pub fn hcat(&self, rhs: &Matrix) -> Matrix {
Matrix::from(af::join(1, &self.backend, &rhs.backend))
}
pub fn hcat_many(matrices: &Vec<Matrix>) -> Matrix {
let backends: Vec<&af::Array> = matrices.into_iter().map(|m| &m.backend).collect();
Matrix::from(af::join_many(1, backends))
}
pub fn vcat(&self, bottom: &Matrix) -> Matrix {
Matrix::from(af::join(0, &self.backend, &bottom.backend))
}
pub fn vcat_many(matrices: &Vec<Matrix>) -> Matrix {
let backends: Vec<&af::Array> = matrices.into_iter().map(|m| &m.backend).collect();
Matrix::from(af::join_many(0, backends))
}
fn check_size(lhs: &Matrix, rhs: &Matrix) -> Result<(), String> {
if lhs.size() == rhs.size() {
Ok(())
} else {
Err(format!("The two Matrices must have the same size. The left MAtrix has size: {} x {} and the right Matrix has size: {} x {}",
lhs.rows(), lhs.cols(), rhs.rows(), rhs.cols()))
}
}
pub fn back_substitution(&self, v: &Vector) -> Result<Vector, String> {
if self.is_empty() {
Err(format!("Matrix is empty"))
} else if self.rows() != self.cols() {
Err(format!("Matrix must be square. Is: {} x {}", self.rows(), self.cols()))
} else if self.rows() != v.size() {
Err(format!("Matrix and Vector must be dimensionally compatible. Matrix rows: {} vs Vector size: {}", self.rows(), v.size()))
} else {
let mut results: Vec<f64> = v.data();
for i in (0..self.rows()).rev() {
let row = self.get_row(i);
let diagonal = self.get(i, i).unwrap();
let dot = {
let row_part: Vec<f64> = row.into_vec()[(i + 1) .. self.rows()].to_vec();
let aug_part: Vec<f64> = results[(i + 1) .. self.rows()].to_vec();
println!("Row part: {:?}", row_part);
println!("Aug part: {:?}", aug_part);
Vector::dot(&Vector::new(row_part), &Vector::new(aug_part)).unwrap()
};
results[i] = (results[i] - dot) / diagonal;
}
Ok(Vector::new(results))
}
}
}
use std::fmt;
impl fmt::Display for Matrix {
fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
write!(f, "[");
for (j, row) in self.row_iter().enumerate() {
let values = row.into_vec();
for (i, v) in values.iter().enumerate() {
if i == values.len()-1 {
write!(f, "{}", v);
} else {
write!(f, "{}\t", v);
}
}
if j != self.rows()-1 {
write!(f, "\n ");
}
}
write!(f, "]")
}
}
impl<'a> From<&'a Matrix> for Matrix {
fn from(matrix: &Matrix) -> Matrix {
matrix.clone()
}
}
impl From<af::Array> for Matrix {
fn from(array: af::Array) -> Matrix {
let dims = array.dims();
let max = ::std::usize::MAX as u64;
assert!(dims.get()[0] < max, "Internal array representation contained more than `usize` rows: {} > {}", dims.get()[0], max);
assert!(dims.get()[1] < max, "Internal array representation contained more than `usize` columns: {} > {}", dims.get()[1], max);
assert!(dims.get()[2] == 1, "Internal array representation contained data in 3D space");
assert!(dims.get()[3] == 1, "Internal array representation contained data in 4D space");
assert!(dims.get()[0] * dims.get()[1] < max, "Internal array representation contained more than `usize` total elements: {} > {}", dims.get()[0] * dims.get()[1], max);
let rows = dims[0] as usize;
let cols = dims[1] as usize;
let mut data: Vec<f64> = vec![0.; rows * cols];
array.host(&mut data);
Matrix::new(rows, cols, data)
}
}
impl From<f64> for Matrix {
fn from(value: f64) -> Matrix {
Matrix::new(1, 1, vec![value])
}
}
impl<'a> From<&'a f64> for Matrix {
fn from(value: &f64) -> Matrix {
Matrix::new(1, 1, vec![*value])
}
}
impl From<Vector> for Matrix {
fn from(vector: Vector) -> Matrix {
Matrix::new(vector.size(), 1, vector.data())
}
}
impl<'a> From<&'a Vector> for Matrix {
fn from(vector: &Vector) -> Matrix {
Matrix::new(vector.size(), 1, vector.data())
}
}
use std::iter::{IntoIterator};
impl IntoIterator for Matrix {
type Item = f64;
type IntoIter = ::std::vec::IntoIter<f64>;
fn into_iter(self) -> Self::IntoIter {
self.data().into_iter()
}
}
impl<'a> IntoIterator for &'a Matrix {
type Item = f64;
type IntoIter = ::std::vec::IntoIter<f64>;
fn into_iter(self) -> Self::IntoIter {
self.data().into_iter()
}
}