rust-matrix 1.0.0-alpha

use crate::pow::Pow;
use crate::sqrt::Sqrt;
use std::fmt;
use std::fmt::Display;
use std::ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Sub};

pub trait MatrixElementRequiredTraits<T>:
    Add<Output = T>
    + Sub<Output = T>
    + Mul<Output = T>
    + Div<Output = T>
    + Add<Output = T>
    + AddAssign
    + MulAssign
    + DivAssign
    + Clone
    + Copy
    + Display
    + PartialOrd<T>
    + Neg<Output = T>
    + Default
    + Sqrt
    + Pow
    + From<u8>
{
}

impl<
        T: Add<Output = T>
            + Sub<Output = T>
            + Mul<Output = T>
            + Div<Output = T>
            + Add<Output = T>
            + AddAssign
            + MulAssign
            + DivAssign
            + Clone
            + Copy
            + Display
            + PartialOrd<T>
            + Neg<Output = T>
            + Default
            + Sqrt
            + Pow
            + From<u8>
            + ?Sized,
    > MatrixElementRequiredTraits<T> for T
{
}

#[derive(PartialEq, Debug, PartialOrd, Clone)]
pub struct Matrix<T: MatrixElementRequiredTraits<T>> {
    pub n: usize,
    pub m: usize,
    pub entries: Vec<T>,
}

impl<T: MatrixElementRequiredTraits<T>> Matrix<T> {
    /// Standard initializer for a Matrix
    /// Takes dimensions `m` (number of rows) and `n` (number of columns)
    /// together with a vector representing the elements. The number of elements
    /// should be equal to `m` x `n`
    /// ```
    /// use matrix::matrix_algebra::Matrix;
    /// let test_matrix = Matrix::new(2, 3, [1.0, 2.0, -1.0, 0.0, 3.0, 7.0].to_vec());
    /// println!("{test_matrix}");
    /// ```
    /// Will output  
    /// ⌜1    2    -1   ⌝  
    /// ⌞0    3    7    ⌟
    pub fn new(m: usize, n: usize, entries: Vec<T>) -> Self {
        Matrix { m, n, entries }
    }

    /// Convenience initializer to initialize a matrix of dimension m x n with all entries
    /// equal to the provided value
    /// ```
    /// use matrix::matrix_algebra::Matrix;
    /// let test_matrix = Matrix::new_constant_value(2, 3, 5.0).expect("Unable to initialize matrix");
    /// println!("{test_matrix}");
    /// ```
    /// /// Will output  
    /// /// ⌜5    5    5    ⌝  
    /// /// ⌞5    5    5    ⌟
    pub fn new_constant_value(m: usize, n: usize, value: T) -> Result<Self, &'static str> {
        if m == 0 || n == 0 {
            return Err("Matrix dimensions must each be greater than zero!");
        }

        Ok(Matrix {
            n,
            m,
            entries: vec![value; n * m],
        })
    }

    /// Getter for entry at position i,j (zero indexed)
    /// Returns a reference to the entry
    /// ```
    /// use matrix::matrix_algebra::Matrix;
    /// let test_matrix = Matrix::new(2, 3, [1.0, 2.0, -1.0, 0.0, 3.0, 7.0].to_vec());
    /// let entry = test_matrix.get_entry_ij(1, 1);
    /// println!("{entry}");
    /// ```
    /// Will output  
    /// 3  
    pub fn get_entry_ij(&self, i: usize, j: usize) -> &T {
        &self.entries[(i * self.n) + j]
    }

    /// Setter for entry at position i, j (zero indexed)
    /// ```
    /// use matrix::matrix_algebra::Matrix;
    /// let mut test_matrix = Matrix::new(2, 3, [1.0, 2.0, -1.0, 0.0, 3.0, 7.0].to_vec());
    /// let entry = test_matrix.get_entry_ij(1, 1);
    /// println!("{entry}");
    /// test_matrix.set_entry_ij(1, 1, &4.0);
    /// let entry = test_matrix.get_entry_ij(1, 1);
    /// println!("{entry}");
    /// ```
    /// Will output:  
    /// 3  
    /// 4  
    pub fn set_entry_ij(&mut self, i: usize, j: usize, new_value: &T) {
        self.entries[(i * self.n) + j] = new_value.clone();
    }

    /// Retrieves the rows of a matrix
    /// ```
    /// use matrix::matrix_algebra::Matrix;
    /// let mut test_matrix = Matrix::new(2, 3, [1.0, 2.0, -1.0, 0.0, 3.0, 7.0].to_vec());
    /// let rows = test_matrix.rows();
    /// ```
    pub fn rows(&self) -> Vec<Vec<T>> {
        let mut rows = Vec::new();

        for i in (0..(self.m * self.n)).step_by(self.n) {
            rows.push(self.entries[i..(i + self.n)].to_vec());
        }

        rows
    }

    /// The first elementary operation on a matrix
    /// Returns a new matrix with the rows at the specified
    /// indices interchanged.
    /// ```
    /// use matrix::matrix_algebra::Matrix;
    /// let test_matrix = Matrix::new(2, 3, [1.0, 2.0, 3.0, 4.0, 5.0, 6.0].to_vec());
    /// let result = test_matrix.row_interchange(0, 1);
    ///
    /// assert_eq!(
    ///    result,
    ///    Matrix::new(2, 3, [4.0, 5.0, 6.0, 1.0, 2.0, 3.0].to_vec()),
    /// );
    /// ```
    pub fn row_interchange(&self, first_row_index: usize, second_row_index: usize) -> Self {
        let rows = self.rows();
        let first_row = &rows[first_row_index];
        let second_row = &rows[second_row_index];

        let mut new_entries: Vec<T> = Vec::new();

        for row_index in 0..rows.len() {
            if row_index == first_row_index {
                for entry in second_row {
                    new_entries.push(*entry);
                }
            } else if row_index == second_row_index {
                for entry in first_row {
                    new_entries.push(*entry);
                }
            } else {
                for entry in rows[row_index].clone() {
                    new_entries.push(entry);
                }
            }
        }

        Matrix {
            m: self.m,
            n: self.n,
            entries: new_entries,
        }
    }

    /// The second elementary operation on a matrix
    /// Returns a new matrix with all elements in the
    ///  row at the specified index multiplied by the
    /// specified scalar.
    /// ```
    /// use matrix::matrix_algebra::Matrix;
    /// let test_matrix = Matrix::new(3, 3, [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0].to_vec());
    /// let result = test_matrix.multiply_row_by_scalar(2, 2.0);
    ///
    /// assert_eq!(
    ///    result,
    ///    Matrix::new(
    ///        3,
    ///        3,
    ///        [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 14.0, 16.0, 18.0].to_vec()
    ///    ),
    /// );
    /// ```
    pub fn multiply_row_by_scalar(&self, row_index: usize, scalar: T) -> Self {
        let rows = self.rows();
        let row_to_be_multiplied = &rows[row_index];

        let mut new_entries: Vec<T> = Vec::new();

        for i in 0..rows.len() {
            if i == row_index {
                for entry in row_to_be_multiplied {
                    new_entries.push(*entry * scalar);
                }
            } else {
                for entry in rows[i].clone() {
                    new_entries.push(entry);
                }
            }
        }

        Matrix {
            m: self.m,
            n: self.n,
            entries: new_entries,
        }
    }

    /// The third elementary operation on a matrix
    /// Adds a scalar multiple of a row to an existing row
    /// in a matrix
    /// ```
    /// use matrix::matrix_algebra::Matrix;
    /// let test_matrix = Matrix::new(3, 3, [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0].to_vec());
    /// let result = test_matrix.add_row_to_scalar_multiple_of_row(0, 2, 5.0);
    ///
    /// assert_eq!(
    ///     result,
    ///     Matrix::new(
    ///         3,
    ///         3,
    ///         [36.0, 42.0, 48.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0].to_vec()
    ///     ),
    /// );
    /// ```
    pub fn add_row_to_scalar_multiple_of_row(
        &self,
        target_index: usize,
        source_index: usize,
        scalar: T,
    ) -> Self {
        let rows = self.rows();
        let row_to_be_added_to = rows[target_index].clone();
        let row_to_be_added: Vec<T> = rows[source_index]
            .clone()
            .into_iter()
            .map(|entry| -> T { scalar * entry })
            .collect();

        let mut new_entries: Vec<T> = Vec::new();

        for i in 0..rows.len() {
            if i == target_index {
                let new_row = row_to_be_added.iter().zip(row_to_be_added_to.iter());
                for (_, (lhs, rhs)) in new_row.enumerate() {
                    new_entries.push(*lhs + *rhs);
                }
            } else {
                for entry in rows[i].clone() {
                    new_entries.push(entry);
                }
            }
        }

        Matrix {
            m: self.m,
            n: self.n,
            entries: new_entries,
        }
    }

    fn all_zeroes(&self) -> bool {
        for i in 0..self.m {
            for j in 0..self.n {
                if *self.get_entry_ij(i, j) != T::default() {
                    return false;
                }
            }
        }

        true
    }

    fn first_row_with_non_zero_entry_in_column(&self, column_index: usize) -> usize {
        let rows = self.rows();
        let mut first_row_with_non_zero_entry = 0;

        for row in rows {
            if row[column_index] != T::default() {
                break;
            }
            first_row_with_non_zero_entry += 1;
        }

        first_row_with_non_zero_entry
    }

    fn reduce_rows_relative_to_row(&self, column_index: usize, row_index: usize) -> Self {
        let mut reduced = self.clone();
        for j in 0..self.m {
            if j != row_index {
                let value_in_non_zero_column_of_row = self.get_entry_ij(j, column_index);
                if *value_in_non_zero_column_of_row != T::default() {
                    let scalar = -T::from(1) * *value_in_non_zero_column_of_row;
                    reduced = reduced.add_row_to_scalar_multiple_of_row(j, row_index, scalar);
                }
            }
        }

        reduced
    }

    fn reduce_first_row(&self, column_index: usize) -> (Self, T) {
        let mut coefficient = T::from(1);
        let first_row_with_non_zero_entry =
            self.first_row_with_non_zero_entry_in_column(column_index);
        let a_k_interchanged = self.row_interchange(0, first_row_with_non_zero_entry);
        let scalar = T::from(1) / *a_k_interchanged.get_entry_ij(0, column_index);
        coefficient = coefficient / scalar;
        (
            a_k_interchanged.multiply_row_by_scalar(0, scalar),
            coefficient,
        )
    }

    fn row_echolon_form_recursive_with_coefficient(&self, k: usize) -> (Self, T) {
        let mut coefficient = T::from(1);
        let mut partitioned: Vec<Matrix<T>> = Vec::new();

        if k != 0 {
            partitioned = self.partition(&[self.n], &[k, self.m - k]);
        }

        let a_k = if k == 0 { self } else { &partitioned[1] };

        if a_k.all_zeroes() {
            return (self.clone(), coefficient);
        }

        let columns = a_k.columns();
        let first_non_zero_column_index = first_non_zero_vec_index(&columns);
        let (mut a_k, new_coefficient) = a_k.reduce_first_row(first_non_zero_column_index);

        coefficient = coefficient * -new_coefficient;

        let combined_entries = if k == 0 {
            &a_k.entries
        } else {
            let _ = &partitioned[0].entries.append(&mut a_k.entries);
            &partitioned[0].entries
        };

        let combined = Matrix::new(self.m, self.n, combined_entries.to_vec());
        let reduced = combined.reduce_rows_relative_to_row(first_non_zero_column_index, k);

        if k != self.m - 1 {
            let (row_echolon, new_coefficient) =
                reduced.row_echolon_form_recursive_with_coefficient(k + 1);
            return (row_echolon, new_coefficient * coefficient);
        }

        (reduced, coefficient)
    }

    fn row_echolon_form_recursive(&self, k: usize) -> Self {
        self.row_echolon_form_recursive_with_coefficient(k).0
    }

    /// Reduces a matrix to Row Echolon Form
    /// ```
    /// use matrix::matrix_algebra::Matrix;
    /// let test_matrix = Matrix::new(
    ///     3,
    ///     4,
    ///     [
    ///         0.0, 0.0, 3.0, -1.0, 0.0, -1.0, 4.0, 7.0, 0.0, -1.0, 7.0, 6.0,
    ///     ]
    ///     .to_vec(),
    /// );
    ///
    /// let result = test_matrix.row_echolon_form();
    ///
    /// assert_eq!(
    ///     result,
    ///     Matrix::new(
    ///         3,
    ///         4,
    ///         [
    ///             0.0,
    ///             1.0,
    ///             0.0,
    ///             -25.0 / 3.0,
    ///             0.0,
    ///             0.0,
    ///             1.0,
    ///             -1.0 / 3.0,
    ///             0.0,
    ///             0.0,
    ///             0.0,
    ///             0.0
    ///         ]
    ///         .to_vec()
    ///     ),
    /// );
    /// ```
    pub fn row_echolon_form(&self) -> Self {
        self.row_echolon_form_recursive(0)
    }

    /// Reduces a matrix to column echolon form
    /// ```
    /// use matrix::matrix_algebra::Matrix;
    /// let test_matrix = Matrix::new(2, 3, [1.0, 2.0, 3.0, 2.0, 3.0, 4.0].to_vec());
    /// let result = test_matrix.column_echolon_form();
    ///
    /// assert_eq!(
    ///     result,
    ///     Matrix::new(2, 3, [1.0, 0.0, 0.0, 0.0, 1.0, 0.0,].to_vec())
    /// );
    /// let test_matrix = Matrix::new(3, 3, [1.0, 1.0, 1.0, 2.0, 2.0, 2.0, 3.0, 3.0, 3.0].to_vec());
    /// let result = test_matrix.column_echolon_form();
    /// assert_eq!(
    ///     result,
    ///     Matrix::new(3, 3, [1.0, 0.0, 0.0, 2.0, 0.0, 0.0, 3.0, 0.0, 0.0].to_vec())
    /// );
    /// ```
    pub fn column_echolon_form(&self) -> Self {
        self.transpose().row_echolon_form().transpose()
    }

    /// Convenience function that extracts the columns
    /// of a matrix
    /// ```
    /// use matrix::matrix_algebra::Matrix;
    /// let test_matrix = Matrix::new(
    /// 3,
    /// 4,
    /// [
    ///    1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0,
    /// ]
    /// .to_vec(),
    /// );
    ///
    /// let columns = test_matrix.columns();
    ///
    /// assert_eq!(columns.len(), 4);
    ///
    /// assert_eq!(columns[0], vec![1.0, 5.0, 9.0]);
    /// assert_eq!(columns[1], vec![2.0, 6.0, 10.0]);
    /// assert_eq!(columns[2], vec![3.0, 7.0, 11.0]);
    /// assert_eq!(columns[3], vec![4.0, 8.0, 12.0]);
    ///
    /// ```
    pub fn columns(&self) -> Vec<Vec<T>> {
        let mut columns = Vec::new();

        for i in 0..self.n {
            let mut new_column: Vec<T> = Vec::new();
            for j in 0..self.m {
                new_column.push(self.entries[(self.n * j) + i].clone());
            }
            columns.push(new_column.to_vec());
        }

        columns
    }

    /// Creates a transpose of a matrix
    /// ```
    /// use matrix::matrix_algebra::Matrix;
    /// let test_matrix = Matrix::new(2, 3, [1.0, 2.0, -1.0, 0.0, 3.0, 7.0].to_vec());
    ///
    /// let transpose_matrix = test_matrix.transpose();
    ///
    /// assert_eq!(transpose_matrix.m, test_matrix.n);
    /// assert_eq!(transpose_matrix.n, test_matrix.m);
    /// assert_eq!(transpose_matrix.entries, [1.0, 0.0, 2.0, 3.0, -1.0, 7.0]);
    /// ```
    pub fn transpose(&self) -> Self {
        let columns = self.columns();
        let mut transposed_entries = Vec::new();

        for column in columns {
            for entry in column {
                transposed_entries.push(entry);
            }
        }

        Matrix::<T> {
            m: self.n,
            n: self.m,
            entries: transposed_entries,
        }
    }

    /// Calculates the determinant of a matrix
    /// The following example shows that there can
    /// be some rounding issues with the calculations.
    /// ```
    /// use matrix::matrix_algebra::Matrix;
    /// use matrix::complex_number::ComplexNumber;
    /// let test_matrix = Matrix::new(
    ///     3,
    ///     3,
    ///     [
    ///         ComplexNumber::new(1.0, 0.0),
    ///         ComplexNumber::new(1.0, 0.0),
    ///         ComplexNumber::new(0.0, 1.0),
    ///         ComplexNumber::new(1.0, 1.0),
    ///         ComplexNumber::new(1.0, 1.0),
    ///         ComplexNumber::new(1.0, 0.0),
    ///         ComplexNumber::new(2.0, 3.0),
    ///         ComplexNumber::new(0.0, -1.0),
    ///         ComplexNumber::new(3.0, 0.0),
    ///     ]
    ///     .to_vec(),
    /// );
    ///
    /// let determinant = test_matrix.determinant();
    ///
    /// fn delta_real(determinant: ComplexNumber<f64>, expectation: f64) -> bool {
    ///     determinant.real - expectation < (1.0 / 1000000000000000000000000000000.0)
    /// }
    ///
    /// fn delta_complex(determinant: ComplexNumber<f64>, expectation: f64) -> bool {
    ///     determinant.complex - expectation < (1.0 / 1000000000000000000000000000000.0)
    /// }
    /// assert!(delta_real(
    ///     determinant.expect("Unable to calculate determinant"),
    ///     8.0
    /// ));
    /// assert!(delta_complex(
    ///     determinant.expect("Unable to calculate determinant"),
    ///     6.0
    /// ));
    /// ```
    pub fn determinant(&self) -> Result<T, &'static str> {
        if self.n != self.m {
            return Err("Non square matrices do not have determinants!");
        }
        let (row_echolon_form, coefficient) = self.row_echolon_form_recursive_with_coefficient(0);
        let mut determinant: Option<T> = None;

        for i in 0..row_echolon_form.n {
            for j in 0..row_echolon_form.m {
                if i == j {
                    let entry = row_echolon_form.get_entry_ij(i, j);
                    match determinant {
                        None => determinant = Some(*entry),
                        Some(det) => determinant = Some(det * *entry),
                    }
                }
            }
        }

        match determinant {
            Some(det) => return Ok(det * coefficient),
            None => Err("Unable to calculate determinant!"),
        }
    }

    /// Calculates the minor of a matrix (the determinant of
    /// a given submatrix) based on the provided column indices
    /// and row indices to be selected from the matrix
    /// ```
    /// use matrix::matrix_algebra::Matrix;
    /// let test_matrix = Matrix::new(3, 3, [1.0, 2.0, 1.0, 6.0, -1.0, 0.0, -1.0, -2.0, -1.0].to_vec());
    /// assert_eq!(test_matrix.minor(&[1, 2], &[1,2]).expect("Unable to create minor"), 1.0);
    /// assert_eq!(test_matrix.minor(&[0, 2], &[1,2]).expect("Unable to create minor"), -6.0);
    /// assert_eq!(test_matrix.minor(&[0, 1], &[0,1]).expect("Unable to create minor"), -13.0);
    /// ```
    pub fn minor(
        &self,
        column_indices: &[usize],
        row_indices: &[usize],
    ) -> Result<T, &'static str> {
        let submatrix = self.submatrix(column_indices, row_indices);

        submatrix.determinant()
    }

    /// Convenience function to identify a non-singular matrix
    pub fn is_nonsingular(&self) -> Result<bool, &'static str> {
        match self.determinant() {
            Err(error) => Err(error),
            Ok(determinant) => Ok(determinant != T::default()),
        }
    }

    /// Creates a matrix of cofactors from a matrix
    pub fn matrix_of_cofactors(&self) -> Result<Self, &'static str> {
        let mut entries: Vec<T> = Vec::new();
        let mut sign = T::from(1);
        for j in 0..self.m {
            for i in 0..self.n {
                let mut column_indices: Vec<usize> = Vec::new();
                for index in 0..self.n {
                    if index != i {
                        column_indices.push(index);
                    }
                }
                let mut row_indices: Vec<usize> = Vec::new();
                for index in 0..self.m {
                    if index != j {
                        row_indices.push(index);
                    }
                }
                let minor = self.minor(&column_indices, &row_indices);
                match minor {
                    Ok(minor) => {
                        entries.push(minor * sign);
                        sign = -sign;
                    }
                    Err(err) => return Err(err),
                }
            }
        }

        Ok(Matrix {
            n: self.n,
            m: self.m,
            entries,
        })
    }

    /// Caluculates the adjoint of a matrix
    /// ```
    /// use matrix::matrix_algebra::Matrix;
    /// let test_matrix = Matrix::new(3, 3, [3.0, 4.0, 3.0, 5.0, 7.0, 2.0, 0.0, 0.0, 1.0].to_vec());
    ///
    /// assert_eq!(
    ///     test_matrix
    ///         .adjoint()
    ///         .expect("test_adjoint: unable to calculate matrix adjoint"),
    ///     Matrix::new(
    ///         3,
    ///         3,
    ///         [
    ///             7.0,
    ///             -4.0,
    ///             -13.0,
    ///             -5.0,
    ///             3.0,
    ///             9.0,
    ///             0.0,
    ///             0.0,
    ///             1.0000000000000018
    ///         ]
    ///         .to_vec()
    ///     )
    /// );
    /// ```
    pub fn adjoint(&self) -> Result<Self, &'static str> {
        match self.matrix_of_cofactors() {
            Ok(matrix_of_cofactors) => Ok(matrix_of_cofactors.transpose()),
            Err(err) => Err(err),
        }
    }

    /// Calculates the inverse of a matrix
    /// ```
    /// use matrix::matrix_algebra::Matrix;
    /// let test_matrix = Matrix::new(
    ///     3,
    ///     3,
    ///     [2.0, 3.0, 0.0, 0.0, 3.0, -3.0, -2.0, 3.0, 3.0].to_vec(),
    /// );
    ///
    /// let inverse = test_matrix
    /// .inverse()
    /// .expect("Unable to compute matrix inverse in test_inverse");
    ///
    /// assert_eq!(inverse.m, 3);
    /// assert_eq!(inverse.n, 3);
    ///
    /// let expected = [
    ///     1.0 / 3.0,
    ///     -1.0 / 6.0,
    ///     -1.0 / 6.0,
    ///     1.0 / 9.0,
    ///     1.0 / 9.0,
    ///     1.0 / 9.0,
    ///     1.0 / 9.0,
    ///     -2.0 / 9.0,
    ///     1.0 / 9.0,
    ///     ]
    ///     .to_vec();
    ///     
    ///     for i in 0..expected.len() {
    ///         assert!(expected[i] - inverse.entries[i] < 1.0 / 1000000000000000.0);
    ///     }
    /// ```
    pub fn inverse(&self) -> Result<Self, &'static str> {
        if self.n != self.m {
            return Err("Unable to compute inverse when n is not equal to m");
        }
        let identity_matrix = new_identity_matrix::<T>(self.m);

        match identity_matrix {
            Ok(identity_matrix) => {
                let identity_matrix_rows = identity_matrix.rows();
                let rows = self.rows();

                let mut new_entries = Vec::new();

                for i in 0..rows.len() {
                    for element in &rows[i] {
                        new_entries.push(*element);
                    }
                    for element in &identity_matrix_rows[i] {
                        new_entries.push(*element);
                    }
                }
                let enlarged = Matrix::new(self.m, 2 * self.n, new_entries);
                let enlarged_ref = enlarged.row_echolon_form();

                return Ok(enlarged_ref.submatrix(
                    &Vec::from_iter(self.n..(2 * self.n)).as_slice(),
                    &Vec::from_iter(0..self.m).as_slice(),
                ));
            }
            Err(err) => return Err(err),
        }
    }

    /// Partitions a matrix based on the provided column and row indices
    /// ```
    /// use matrix::matrix_algebra::Matrix;
    /// let test_matrix = Matrix::new(
    /// 5,
    /// 5,
    /// [
    ///     1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0,
    ///     16.0, 17.0, 18.0, 19.0, 20.0, 21.0, 22.0, 23.0, 24.0, 25.0,
    /// ]
    /// .to_vec(),
    /// );
    ///
    /// let submatrices = test_matrix.partition(&[3, 2], &[2, 3]);
    ///
    /// assert_eq!(
    ///     submatrices,
    ///     [
    ///         Matrix::new(2, 3, [1.0, 2.0, 3.0, 6.0, 7.0, 8.0].to_vec()),
    ///         Matrix::new(2, 2, [4.0, 5.0, 9.0, 10.0].to_vec()),
    ///         Matrix::new(
    ///             3,
    ///             3,
    ///             [11.0, 12.0, 13.0, 16.0, 17.0, 18.0, 21.0, 22.0, 23.0].to_vec()
    ///         ),
    ///         Matrix::new(3, 2, [14.0, 15.0, 19.0, 20.0, 24.0, 25.0].to_vec())
    ///     ]
    /// );
    /// ```
    pub fn partition(
        &self,
        column_partitioning: &[usize],
        row_partitioning: &[usize],
    ) -> Vec<Self> {
        let mut partitioned_matrices: Vec<Matrix<T>> = Vec::new();
        let mut row_offset = 0;

        for row_partition in row_partitioning {
            let mut column_offset = 0;
            for column_partition in column_partitioning {
                let mut new_entries: Vec<T> = Vec::new();

                for i in row_offset..(row_offset + *row_partition) {
                    for j in column_offset..(column_offset + *column_partition) {
                        new_entries.push(*self.get_entry_ij(i, j));
                    }
                }
                partitioned_matrices.push(Matrix {
                    m: *row_partition,
                    n: *column_partition,
                    entries: new_entries,
                });
                column_offset += *column_partition;
            }
            row_offset += *row_partition;
        }

        partitioned_matrices
    }

    /// Generates a submatrix based on the provided column
    /// and row indices
    /// ```
    /// use matrix::matrix_algebra::Matrix;
    /// let test_matrix = Matrix::new(
    ///     5,
    ///     5,
    ///     [
    ///         1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0,
    ///         16.0, 17.0, 18.0, 19.0, 20.0, 21.0, 22.0, 23.0, 24.0, 25.0,
    ///     ]
    ///     .to_vec(),
    /// );
    ///
    /// let submatrix = test_matrix.submatrix(&[1, 3], &[0, 2, 4]);
    /// assert_eq!(submatrix.n, 2);
    /// assert_eq!(submatrix.m, 3);
    /// assert_eq!(submatrix.entries, [2.0, 4.0, 12.0, 14.0, 22.0, 24.0]);
    /// ```
    pub fn submatrix(&self, column_indices: &[usize], row_indices: &[usize]) -> Self {
        let rows = self.rows();
        let mut new_entries: Vec<T> = Vec::new();

        for row_index in row_indices {
            for column_index in column_indices {
                let new_entry = &rows[*row_index][*column_index];

                new_entries.push(*new_entry);
            }
        }

        Matrix {
            n: column_indices.len(),
            m: row_indices.len(),
            entries: new_entries,
        }
    }
}

/// Implementation of the Display trait for Matrix<T>
/// which can be used for debugging purposes
impl<T: MatrixElementRequiredTraits<T>> fmt::Display for Matrix<T> {
    fn fmt(&self, fmt: &mut fmt::Formatter) -> fmt::Result {
        let rows = self.rows();
        for (i, row) in rows.clone().into_iter().enumerate() {
            let printable_row: String = row
                .iter()
                .map(|&entry| {
                    let mut entry_as_string = entry.to_string();

                    while entry_as_string.len() < 4 {
                        entry_as_string += " ";
                    }

                    entry_as_string + " "
                })
                .collect();
            if i == 0 {
                let _ = fmt.write_str("⌜");
            } else if i == &rows.len() - 1 {
                let _ = fmt.write_str("⌞");
            } else {
                let _ = fmt.write_str("⎸");
            }

            let _ = fmt.write_str(&printable_row)?;
            if i == 0 {
                let _ = fmt.write_str("⌝");
            } else if i == rows.len() - 1 {
                let _ = fmt.write_str("⌟");
            } else {
                let _ = fmt.write_str("⎹");
            }
            let _ = fmt.write_str("\n");
        }
        Ok(())
    }
}

/// Convenience initialiser that creates a `Matrix<T>` using the default
/// value of the type
/// ```
/// use matrix::matrix_algebra::{Matrix, new_all_default};
/// use rand::prelude::*;
/// let mut rng = rand::thread_rng();
/// let n = rng.gen_range(1..=100);
/// let m = rng.gen_range(1..=100);
/// let test_matrix_f64 = new_all_default::<f64>(m, n)
///     .expect("test_new_all_default: unable to create test_matrix_f64");
///
/// for entry in test_matrix_f64.entries {
///     assert_eq!(entry, f64::default());
/// }
///
/// let test_matrix_f32 = new_all_default::<f32>(m, n)
///     .expect("test_new_all_default: unable to create test_matrix_f64");
///
/// for entry in test_matrix_f32.entries {
///     assert_eq!(entry, f32::default());
/// }
/// ```
pub fn new_all_default<T>(m: usize, n: usize) -> Result<Matrix<T>, &'static str>
where
    T: MatrixElementRequiredTraits<T>,
{
    if m == 0 || n == 0 {
        return Err("Matrix dimensions must each be greater than zero!");
    }

    Ok(Matrix::<T> {
        n,
        m,
        entries: vec![T::default(); n * m],
    })
}

/// Creates an identity matrix of the specified dimension using whatever
/// corresponds to unity for type T
/// ```
/// use matrix::matrix_algebra::{Matrix, new_identity_matrix};
/// use rand::prelude::*;
/// let mut rng = rand::thread_rng();
/// let dimension = rng.gen_range(1..=100);
///
/// let test_matrix_f64 = new_identity_matrix::<f64>(dimension)
///     .expect("test_new_identity_matrix: unable to create test_matrix_f64");
///
/// for index in 0..dimension {
///     assert_eq!(*test_matrix_f64.get_entry_ij(index, index), 1.0);
/// }
///
/// let test_matrix_f32 = new_identity_matrix::<f32>(dimension)
///     .expect("test_new_identity_matrix: unable to create test_matrix_f64");
///
/// for index in 0..dimension {
///     assert_eq!(*test_matrix_f32.get_entry_ij(index, index), 1.0);
/// }
/// ```
pub fn new_identity_matrix<T>(dimension: usize) -> Result<Matrix<T>, &'static str>
where
    T: MatrixElementRequiredTraits<T>,
{
    let new_empty_matrix = new_all_default::<T>(dimension, dimension);

    match new_empty_matrix {
        Err(err) => return Err(err),
        Ok(mut new_empty_matrix) => {
            for index in 0..dimension {
                new_empty_matrix.set_entry_ij(index, index, &T::from(1));
            }
            Ok(new_empty_matrix)
        }
    }
}

/// Implementation of matrix addition
/// ```
/// use matrix::matrix_algebra::Matrix;
/// let test_matrix_a = Matrix::new(
///     3,
///     4,
///     [
///         1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0,
///     ]
///     .to_vec(),
/// );
/// let test_matrix_b = Matrix::new(
///     3,
///     4,
///     [
///         12.0, 11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 1.0,
///     ]
///     .to_vec(),
/// );
///
/// let matrix_sum = test_matrix_a + test_matrix_b;
///
/// assert_eq!(matrix_sum.entries.len(), 12);
/// assert_eq!(
///     matrix_sum.entries,
///     [13.0, 13.0, 13.0, 13.0, 13.0, 13.0, 13.0, 13.0, 13.0, 13.0, 13.0, 13.0]
/// );
/// ```
fn matrix_add<T: MatrixElementRequiredTraits<T>>(
    a: &Matrix<T>,
    b: &Matrix<T>,
) -> Result<Matrix<T>, &'static str>
where
    for<'a> &'a T: Add<Output = T>,
{
    if !is_additively_conformable(a, b) {
        return Err("Matrices are not additively conformable!");
    }
    let mut entries: Vec<T> = Vec::new();

    for i in 0..a.m {
        for j in 0..a.n {
            let new_value = a.get_entry_ij(i, j) + b.get_entry_ij(i, j);
            entries.push(new_value);
        }
    }

    Ok(Matrix::<T> {
        n: a.m,
        m: a.n,
        entries,
    })
}

/// Implementation of the `Add` trait for `Matrix<T>`
/// Uses `matrix_add`
impl<T: MatrixElementRequiredTraits<T>> Add for Matrix<T>
where
    for<'a> &'a T: Add<Output = T>,
{
    type Output = Self;

    fn add(self, other: Self) -> Self {
        match matrix_add(&self, &other) {
            Ok(result) => return result,
            Err(err) => panic!("{err}"),
        }
    }
}

/// Implementation of the `Mul` trait for `Matrix<T>`
/// Uses `matrix_multiply`
/// ```
/// use matrix::matrix_algebra::Matrix;
/// let test_matrix_a = Matrix::new(
///     3,
///     4,
///     [
///         1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0,
///     ]
///     .to_vec(),
/// );
/// let test_matrix_b = Matrix::new(
///     4,
///     3,
///     [
///         12.0, 11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 1.0,
///     ]
///     .to_vec(),
/// );
///
/// let matrix_product = test_matrix_a * test_matrix_b;
///
/// assert_eq!(matrix_product.entries.len(), 9);
///
/// assert_eq!(
/// matrix_product.entries,
///     [
///         1.0 * 12.0 + 2.0 * 9.0 + 3.0 * 6.0 + 4.0 * 3.0,
///         1.0 * 11.0 + 2.0 * 8.0 + 3.0 * 5.0 + 4.0 * 2.0,
///         1.0 * 10.0 + 2.0 * 7.0 + 3.0 * 4.0 + 4.0 * 1.0,
///         5.0 * 12.0 + 6.0 * 9.0 + 7.0 * 6.0 + 8.0 * 3.0,
///         5.0 * 11.0 + 6.0 * 8.0 + 7.0 * 5.0 + 8.0 * 2.0,
///         5.0 * 10.0 + 6.0 * 7.0 + 7.0 * 4.0 + 8.0 * 1.0,
///         9.0 * 12.0 + 10.0 * 9.0 + 11.0 * 6.0 + 12.0 * 3.0,
///         9.0 * 11.0 + 10.0 * 8.0 + 11.0 * 5.0 + 12.0 * 2.0,
///         9.0 * 10.0 + 10.0 * 7.0 + 11.0 * 4.0 + 12.0 * 1.0
///     ]
/// );
/// ```
impl<T: MatrixElementRequiredTraits<T>> Mul for Matrix<T> {
    type Output = Self;

    fn mul(self, rhs: Self) -> Self {
        match matrix_multiply::<T>(&self, &rhs) {
            Err(err) => panic!("{err}"),
            Ok(result) => result,
        }
    }
}

/// Helper function to check that two matrices
/// can be multiplied
fn is_multiplicatively_conformable<T: MatrixElementRequiredTraits<T>>(
    a: &Matrix<T>,
    b: &Matrix<T>,
) -> bool {
    a.n == b.m
}

/// Helper function to check that two matrices
/// can be summed
fn is_additively_conformable<T: MatrixElementRequiredTraits<T>>(
    a: &Matrix<T>,
    b: &Matrix<T>,
) -> bool {
    a.n == b.n && a.m == b.m
}

/// Multiplication of two matrices
/// This is currently a naive implementation with no
/// optimisations.
fn matrix_multiply<T: MatrixElementRequiredTraits<T>>(
    a: &Matrix<T>,
    b: &Matrix<T>,
) -> Result<Matrix<T>, &'static str> {
    if !is_multiplicatively_conformable(&a, &b) {
        return Err("Matrices are not multiplicatively conformable!");
    }
    let mut entries: Vec<T> = Vec::new();
    let a_rows = a.rows();
    let b_columns = b.columns();

    for i in 0..a_rows.len() {
        for j in 0..b_columns.len() {
            let mut new_entry: Option<T> = None;
            for k in 0..a_rows[i].len() {
                match new_entry {
                    Some(value) => new_entry = Some(value + a_rows[i][k] * b_columns[j][k]),
                    None => new_entry = Some(a_rows[i][k] * b_columns[j][k]),
                }
            }

            if entries.len() == 0 {
                entries = vec![new_entry.expect("new_entry not initialised!")];
            } else {
                entries.push(new_entry.expect("new_entry not initialised!"));
            }
        }
    }

    Ok(Matrix::<T> {
        m: a_rows.len(),
        n: b_columns.len(),
        entries,
    })
}

fn first_non_zero_vec_index<T: MatrixElementRequiredTraits<T>>(input: &Vec<Vec<T>>) -> usize {
    let mut first_nonzero_vec_index = 0;

    for vector in input {
        for element in vector {
            if *element != T::default() {
                return first_nonzero_vec_index;
            }
        }
        first_nonzero_vec_index += 1;
    }

    first_nonzero_vec_index
}

#[cfg(test)]
mod tests {
    use std::ops::Add;

    use rand::prelude::*;

    use crate::{complex_number::ComplexNumber, matrix_algebra::first_non_zero_vec_index};

    use super::{new_all_default, new_identity_matrix, Matrix, MatrixElementRequiredTraits};

    fn index_of_first_non_zero_element_in_vec<T: MatrixElementRequiredTraits<T>>(
        input: &Vec<T>,
    ) -> usize {
        let mut index_of_first_non_zero_element = 0;
        for element in input {
            if *element != T::default() {
                return index_of_first_non_zero_element;
            }
            index_of_first_non_zero_element += 1;
        }

        index_of_first_non_zero_element
    }

    #[test]
    fn test_constant_value_initialiser() {
        let mut rng = rand::thread_rng();
        let n = rng.gen_range(1..=100);
        let m = rng.gen_range(1..=100);
        let value: f64 = rng.gen_range(1.0..=100.0);
        let test_matrix =
            Matrix::new_constant_value(m, n, value).expect("Unable to create test_matrix");

        assert_eq!(test_matrix.entries.len(), n * m);

        for entry in test_matrix.entries {
            assert_eq!(entry, value);
        }
    }

    #[test]
    fn test_initialiser() {
        let mut rng = rand::thread_rng();
        let n = rng.gen_range(1..=100);
        let m = rng.gen_range(1..=100);
        let mut entries: Vec<f32> = Vec::new();

        for _i in 0..m {
            for _j in 0..n {
                entries.push(rng.gen_range(1.0..=100.0));
            }
        }

        let test_matrix = Matrix { m, n, entries };

        assert_eq!(test_matrix.entries.len(), n * m);
    }

    #[test]
    fn test_all_zeroes_initialiser() {
        let mut rng = rand::thread_rng();
        let n = rng.gen_range(1..=100);
        let m = rng.gen_range(1..=100);
        let test_matrix =
            Matrix::new_constant_value(m, n, 0.0).expect("Unable to create test_matrix");

        assert_eq!(test_matrix.entries.len(), n * m);

        for entry in test_matrix.entries {
            assert_eq!(entry, 0.0);
        }
    }

    #[test]
    fn test_zero_first_argument_to_initialiser() {
        let test_matrix = Matrix::new_constant_value(0, 1, 1.0);
        assert_eq!(
            test_matrix,
            Err("Matrix dimensions must each be greater than zero!")
        );
    }

    #[test]
    fn test_zero_second_argument_to_initialiser() {
        let test_matrix = Matrix::new_constant_value(1, 0, 1.0);
        assert_eq!(
            test_matrix,
            Err("Matrix dimensions must each be greater than zero!")
        );
    }

    #[test]
    #[should_panic]
    fn test_panic_on_non_multiplicatively_conformable_matrices() {
        let test_matrix_a =
            Matrix::new_constant_value(3, 4, 5.0).expect("Unable to create test_matrix_a");
        let test_matrix_b =
            Matrix::new_constant_value(5, 7, 4.0).expect("Unable to create test_matrix_b");

        let _ = test_matrix_a * test_matrix_b;
    }

    #[test]
    #[should_panic]
    fn test_panic_on_non_additively_conformable_matrices() {
        let test_matrix_a =
            Matrix::new_constant_value(3, 4, 5.0).expect("Unable to create test_matrix_a");
        let test_matrix_b =
            Matrix::new_constant_value(5, 7, 4.0).expect("Unable to create test_matrix_b");

        let _ = test_matrix_a.add(test_matrix_b);
    }

    #[test]
    fn test_matrix_add() {
        let test_matrix_a = Matrix::new(
            3,
            4,
            [
                1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0,
            ]
            .to_vec(),
        );
        let test_matrix_b = Matrix::new(
            3,
            4,
            [
                12.0, 11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 1.0,
            ]
            .to_vec(),
        );

        let matrix_sum = test_matrix_a + test_matrix_b;

        assert_eq!(matrix_sum.entries.len(), 12);
        assert_eq!(
            matrix_sum.entries,
            [13.0, 13.0, 13.0, 13.0, 13.0, 13.0, 13.0, 13.0, 13.0, 13.0, 13.0, 13.0]
        );
    }

    #[test]
    fn test_matrix_add_operator() {
        let test_matrix_a = Matrix::new(
            3,
            4,
            [
                1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0,
            ]
            .to_vec(),
        );
        let test_matrix_b = Matrix::new(
            3,
            4,
            [
                12.0, 11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 1.0,
            ]
            .to_vec(),
        );
        let matrix_sum = test_matrix_a + test_matrix_b;

        assert_eq!(matrix_sum.entries.len(), 12);
        assert_eq!(
            matrix_sum.entries,
            [13.0, 13.0, 13.0, 13.0, 13.0, 13.0, 13.0, 13.0, 13.0, 13.0, 13.0, 13.0,]
        );
    }

    #[test]
    fn test_columns() {
        let test_matrix = Matrix::new(
            3,
            4,
            [
                1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0,
            ]
            .to_vec(),
        );

        let columns = test_matrix.columns();

        assert_eq!(columns.len(), 4);

        assert_eq!(columns[0], vec![1.0, 5.0, 9.0]);
        assert_eq!(columns[1], vec![2.0, 6.0, 10.0]);
        assert_eq!(columns[2], vec![3.0, 7.0, 11.0]);
        assert_eq!(columns[3], vec![4.0, 8.0, 12.0]);
    }

    #[test]
    fn test_rows() {
        let test_matrix = Matrix::new(
            3,
            4,
            [
                1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0,
            ]
            .to_vec(),
        );

        let rows = test_matrix.rows();

        assert_eq!(rows.len(), 3);

        assert_eq!(rows[0], vec![1.0, 2.0, 3.0, 4.0,]);
        assert_eq!(rows[1], vec![5.0, 6.0, 7.0, 8.0,]);
        assert_eq!(rows[2], vec![9.0, 10.0, 11.0, 12.0,]);
    }

    #[test]
    fn test_matrix_multiply() {
        let test_matrix_a = Matrix::new(
            3,
            4,
            [
                1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0,
            ]
            .to_vec(),
        );
        let test_matrix_b = Matrix::new(
            4,
            3,
            [
                12.0, 11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 1.0,
            ]
            .to_vec(),
        );

        let matrix_product = test_matrix_a * test_matrix_b;

        assert_eq!(matrix_product.entries.len(), 9);

        assert_eq!(
            matrix_product.entries,
            [
                1.0 * 12.0 + 2.0 * 9.0 + 3.0 * 6.0 + 4.0 * 3.0,
                1.0 * 11.0 + 2.0 * 8.0 + 3.0 * 5.0 + 4.0 * 2.0,
                1.0 * 10.0 + 2.0 * 7.0 + 3.0 * 4.0 + 4.0 * 1.0,
                5.0 * 12.0 + 6.0 * 9.0 + 7.0 * 6.0 + 8.0 * 3.0,
                5.0 * 11.0 + 6.0 * 8.0 + 7.0 * 5.0 + 8.0 * 2.0,
                5.0 * 10.0 + 6.0 * 7.0 + 7.0 * 4.0 + 8.0 * 1.0,
                9.0 * 12.0 + 10.0 * 9.0 + 11.0 * 6.0 + 12.0 * 3.0,
                9.0 * 11.0 + 10.0 * 8.0 + 11.0 * 5.0 + 12.0 * 2.0,
                9.0 * 10.0 + 10.0 * 7.0 + 11.0 * 4.0 + 12.0 * 1.0
            ]
        );
    }

    #[test]
    fn test_matrix_multiply_constant_value_initialiser() {
        let test_matrix_a =
            Matrix::new_constant_value(3, 4, 5.0).expect("Unable to create test_matrix_a");
        let test_matrix_b =
            Matrix::new_constant_value(4, 3, 4.0).expect("Unable to create test_matrix_b");

        let matrix_product = test_matrix_a * test_matrix_b;

        assert_eq!(matrix_product.entries.len(), 9);

        assert_eq!(
            matrix_product.entries,
            [80.0, 80.0, 80.0, 80.0, 80.0, 80.0, 80.0, 80.0, 80.0,]
        );
    }

    #[test]
    fn test_transpose() {
        let test_matrix = Matrix::new(2, 3, [1.0, 2.0, -1.0, 0.0, 3.0, 7.0].to_vec());

        let transpose_matrix = test_matrix.transpose();

        assert_eq!(transpose_matrix.m, test_matrix.n);
        assert_eq!(transpose_matrix.n, test_matrix.m);
        assert_eq!(transpose_matrix.entries, [1.0, 0.0, 2.0, 3.0, -1.0, 7.0]);
    }

    #[test]
    fn test_new_all_default() {
        let mut rng = rand::thread_rng();
        let n = rng.gen_range(1..=100);
        let m = rng.gen_range(1..=100);
        let test_matrix_f64 = new_all_default::<f64>(m, n)
            .expect("test_new_all_default: unable to create test_matrix_f64");

        for entry in test_matrix_f64.entries {
            assert_eq!(entry, f64::default());
        }

        let test_matrix_f32 = new_all_default::<f32>(m, n)
            .expect("test_new_all_default: unable to create test_matrix_f64");

        for entry in test_matrix_f32.entries {
            assert_eq!(entry, f32::default());
        }
    }

    #[test]
    fn test_new_identity_matrix() {
        let mut rng = rand::thread_rng();
        let dimension = rng.gen_range(1..=100);

        let test_matrix_f64 = new_identity_matrix::<f64>(dimension)
            .expect("test_new_identity_matrix: unable to create test_matrix_f64");

        for index in 0..dimension {
            assert_eq!(*test_matrix_f64.get_entry_ij(index, index), 1.0);
        }

        let test_matrix_f32 = new_identity_matrix::<f32>(dimension)
            .expect("test_new_identity_matrix: unable to create test_matrix_f64");

        for index in 0..dimension {
            assert_eq!(*test_matrix_f32.get_entry_ij(index, index), 1.0);
        }
    }

    #[test]
    fn test_get_entry() {
        let entries = [
            1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0,
            17.0, 18.0, 19.0, 20.0, 21.0, 22.0, 23.0, 24.0, 25.0,
        ]
        .to_vec();
        let test_matrix = Matrix::new(5, 5, entries);
        let mut index = 1.0;

        for i in 0..5 {
            for j in 0..5 {
                assert_eq!(*test_matrix.get_entry_ij(i, j), index);
                index += 1.0;
            }
        }
    }

    #[test]
    fn test_set_entry() {
        let mut test_matrix =
            new_all_default::<f64>(5, 5).expect("test_set_entry: unable to create test_matrix");
        let mut index = 1.0;

        for i in 0..5 {
            for j in 0..5 {
                test_matrix.set_entry_ij(i, j, &index);
                assert_eq!(*test_matrix.get_entry_ij(i, j), index);
                index += 1.0;
            }
        }
    }

    #[test]
    fn test_partition() {
        let test_matrix = Matrix::new(
            5,
            5,
            [
                1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0,
                16.0, 17.0, 18.0, 19.0, 20.0, 21.0, 22.0, 23.0, 24.0, 25.0,
            ]
            .to_vec(),
        );

        let submatrices = test_matrix.partition(&[3, 2], &[2, 3]);

        assert_eq!(
            submatrices,
            [
                Matrix::new(2, 3, [1.0, 2.0, 3.0, 6.0, 7.0, 8.0].to_vec()),
                Matrix::new(2, 2, [4.0, 5.0, 9.0, 10.0].to_vec()),
                Matrix::new(
                    3,
                    3,
                    [11.0, 12.0, 13.0, 16.0, 17.0, 18.0, 21.0, 22.0, 23.0].to_vec()
                ),
                Matrix::new(3, 2, [14.0, 15.0, 19.0, 20.0, 24.0, 25.0].to_vec())
            ]
        );

        let test_matrix = Matrix::new(
            3,
            4,
            [
                0.0, 1.0, -4.0, -7.0, 0.0, 0.0, 3.0, -1.0, 0.0, 0.0, 3.0, -1.0,
            ]
            .to_vec(),
        );

        let submatrices = test_matrix.partition(&[4], &[1, 2]);

        assert_eq!(
            submatrices,
            [
                Matrix::new(1, 4, [0.0, 1.0, -4.0, -7.0].to_vec()),
                Matrix::new(2, 4, [0.0, 0.0, 3.0, -1.0, 0.0, 0.0, 3.0, -1.0,].to_vec()),
            ]
        );
    }

    #[test]
    fn test_submatrix() {
        let test_matrix = Matrix::new(
            5,
            5,
            [
                1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0,
                16.0, 17.0, 18.0, 19.0, 20.0, 21.0, 22.0, 23.0, 24.0, 25.0,
            ]
            .to_vec(),
        );

        let submatrix = test_matrix.submatrix(&[1, 3], &[0, 2, 4]);

        assert_eq!(submatrix.n, 2);
        assert_eq!(submatrix.m, 3);
        assert_eq!(submatrix.entries, [2.0, 4.0, 12.0, 14.0, 22.0, 24.0]);
    }

    #[test]
    fn test_row_interchange() {
        let test_matrix = Matrix::new(2, 3, [1.0, 2.0, 3.0, 4.0, 5.0, 6.0].to_vec());

        let result = test_matrix.row_interchange(0, 1);

        assert_eq!(
            result,
            Matrix::new(2, 3, [4.0, 5.0, 6.0, 1.0, 2.0, 3.0].to_vec()),
        );
    }

    #[test]
    fn test_multiply_row_by_scalar() {
        let test_matrix = Matrix::new(3, 3, [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0].to_vec());

        let result = test_matrix.multiply_row_by_scalar(2, 2.0);

        assert_eq!(
            result,
            Matrix::new(
                3,
                3,
                [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 14.0, 16.0, 18.0].to_vec()
            ),
        );
    }

    #[test]
    fn test_add_row_to_scalar_multiple_of_row() {
        let test_matrix = Matrix::new(3, 3, [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0].to_vec());

        let result = test_matrix.add_row_to_scalar_multiple_of_row(0, 2, 5.0);

        assert_eq!(
            result,
            Matrix::new(
                3,
                3,
                [36.0, 42.0, 48.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0].to_vec()
            ),
        );
    }

    #[test]
    fn test_all_zeroes() {
        let mut entries = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0].to_vec();
        let test_matrix = Matrix::new(3, 4, entries.clone());

        assert_eq!(true, test_matrix.all_zeroes());

        let mut rng = rand::thread_rng();
        let random_index = rng.gen_range(0..12);
        entries[random_index] = 1.0;

        let test_matrix = Matrix::new(3, 4, entries.clone());

        assert_eq!(false, test_matrix.all_zeroes());
    }

    #[test]
    fn test_first_non_zero_vec_index() {
        let mut rng = rand::thread_rng();
        let random_index = rng.gen_range(0..12);
        let mut test_vec: Vec<Vec<f64>> = vec![];
        let zero_vec = [0.0, 0.0, 0.0].to_vec();
        let non_zero_vec = [1.0, 2.0, 3.0].to_vec();

        for i in 0..12 {
            if i == random_index {
                test_vec.push(non_zero_vec.clone());
            } else {
                test_vec.push(zero_vec.clone());
            }
        }

        assert_eq!(random_index, first_non_zero_vec_index(&test_vec));
    }

    #[test]
    fn test_index_of_first_non_zero_element_in_vec() {
        let mut rng = rand::thread_rng();
        let random_index = rng.gen_range(0..12);
        let mut test_vec: Vec<f64> = vec![];

        for i in 0..12 {
            if i == random_index {
                test_vec.push(1.0);
            } else {
                test_vec.push(0.0);
            }
        }

        assert_eq!(
            random_index,
            index_of_first_non_zero_element_in_vec(&test_vec)
        );
    }

    #[test]
    fn test_row_echolon_form() {
        let test_matrix = Matrix::new(
            3,
            4,
            [
                0.0, 0.0, 3.0, -1.0, 0.0, -1.0, 4.0, 7.0, 0.0, -1.0, 7.0, 6.0,
            ]
            .to_vec(),
        );

        let result = test_matrix.row_echolon_form();

        assert_eq!(
            result,
            Matrix::new(
                3,
                4,
                [
                    0.0,
                    1.0,
                    0.0,
                    -25.0 / 3.0,
                    0.0,
                    0.0,
                    1.0,
                    -1.0 / 3.0,
                    0.0,
                    0.0,
                    0.0,
                    0.0
                ]
                .to_vec()
            ),
        );

        let test_matrix = Matrix::new(
            3,
            3,
            [4.0, -8.0, 16.0, 1.0, -3.0, 6.0, 2.0, 1.0, 1.0].to_vec(),
        );

        let result = test_matrix.row_echolon_form();

        assert_eq!(
            result,
            new_identity_matrix(3)
                .expect("test_row_echolon_form: unable to create new_identity_matrix")
        );

        let test_matrix = Matrix::new(
            3,
            4,
            [0.0, 2.0, 1.0, 4.0, 0.0, 0.0, 2.0, 6.0, 1.0, 0.0, -3.0, 2.0].to_vec(),
        );

        let result = test_matrix.row_echolon_form();

        assert_eq!(
            result,
            Matrix::new(
                3,
                4,
                [1.0, 0.0, 0.0, 11.0, 0.0, 1.0, 0.0, 0.5, 0.0, 0.0, 1.0, 3.0].to_vec()
            )
        );

        let test_matrix = Matrix::new(
            3,
            3,
            [
                ComplexNumber::new(2.0, 0.0),
                ComplexNumber::new(8.0, 2.0),
                ComplexNumber::new(6.0, -6.0),
                ComplexNumber::new(0.0, 1.0),
                ComplexNumber::new(0.0, 5.0),
                ComplexNumber::new(3.0, 3.0),
                ComplexNumber::new(1.0, 2.0),
                ComplexNumber::new(4.0, 11.0),
                ComplexNumber::new(9.0, 3.0),
            ]
            .to_vec(),
        );

        let result = test_matrix.row_echolon_form();

        assert_eq!(
            result,
            Matrix::new(
                3,
                3,
                [
                    ComplexNumber::new(1.0, 0.0),
                    ComplexNumber::new(0.0, 0.0),
                    ComplexNumber::new(3.0, -3.0),
                    ComplexNumber::new(0.0, 0.0),
                    ComplexNumber::new(1.0, 0.0),
                    ComplexNumber::new(0.0, 0.0),
                    ComplexNumber::new(0.0, 0.0),
                    ComplexNumber::new(0.0, 0.0),
                    ComplexNumber::new(0.0, 0.0),
                ]
                .to_vec()
            )
        );

        let test_matrix = Matrix::new(
            3,
            2,
            [
                ComplexNumber::new(1.0, 1.0),
                ComplexNumber::new(1.0, -1.0),
                ComplexNumber::new(2.0, 0.0),
                ComplexNumber::new(2.0, 0.0),
                ComplexNumber::new(1.0, 2.0),
                ComplexNumber::new(2.0, -1.0),
            ]
            .to_vec(),
        );

        let result = test_matrix.row_echolon_form();

        assert_eq!(
            result,
            Matrix::new(
                3,
                2,
                [
                    ComplexNumber::new(1.0, 0.0),
                    ComplexNumber::new(0.0, 0.0),
                    ComplexNumber::new(0.0, 0.0),
                    ComplexNumber::new(1.0, 0.0),
                    ComplexNumber::new(0.0, 0.0),
                    ComplexNumber::new(0.0, 0.0),
                ]
                .to_vec()
            )
        );
    }

    #[test]
    fn test_column_echolon_form() {
        let test_matrix = Matrix::new(2, 3, [1.0, 2.0, 3.0, 2.0, 3.0, 4.0].to_vec());

        let result = test_matrix.column_echolon_form();

        assert_eq!(
            result,
            Matrix::new(2, 3, [1.0, 0.0, 0.0, 0.0, 1.0, 0.0,].to_vec())
        );

        let test_matrix = Matrix::new(3, 3, [1.0, 1.0, 1.0, 2.0, 2.0, 2.0, 3.0, 3.0, 3.0].to_vec());

        let result = test_matrix.column_echolon_form();

        assert_eq!(
            result,
            Matrix::new(3, 3, [1.0, 0.0, 0.0, 2.0, 0.0, 0.0, 3.0, 0.0, 0.0].to_vec())
        );
    }

    #[test]
    fn test_determinant() {
        let test_matrix = Matrix::new(
            4,
            4,
            [
                2.0, 2.0, -1.0, 9.0, 4.0, 2.0, 1.0, 17.0, 1.0, -1.0, 1.0, 2.0, 2.0, 1.0, 1.0, 9.0,
            ]
            .to_vec(),
        );
        let result = test_matrix.determinant();
        assert_eq!(result.expect("Unable to calculate determinant"), 0.0);

        let test_matrix = Matrix::new(
            3,
            3,
            [
                ComplexNumber::new(1.0, 0.0),
                ComplexNumber::new(1.0, 0.0),
                ComplexNumber::new(0.0, 1.0),
                ComplexNumber::new(1.0, 1.0),
                ComplexNumber::new(1.0, 1.0),
                ComplexNumber::new(1.0, 0.0),
                ComplexNumber::new(2.0, 3.0),
                ComplexNumber::new(0.0, -1.0),
                ComplexNumber::new(3.0, 0.0),
            ]
            .to_vec(),
        );

        let determinant = test_matrix.determinant();

        fn delta_real(determinant: ComplexNumber<f64>, expectation: f64) -> bool {
            determinant.real - expectation < (1.0 / 1000000000000000000000000000000.0)
        }

        fn delta_complex(determinant: ComplexNumber<f64>, expectation: f64) -> bool {
            determinant.complex - expectation < (1.0 / 1000000000000000000000000000000.0)
        }
        assert!(delta_real(
            determinant.expect("Unable to calculate determinant"),
            8.0
        ));
        assert!(delta_complex(
            determinant.expect("Unable to calculate determinant"),
            6.0
        ));
    }

    #[test]
    fn test_determinant_err() {
        let test_matrix = new_all_default::<f64>(3, 2)
            .expect("test_determinant_err: unable to create tesdt_matrix");

        let result = test_matrix.determinant();

        assert_eq!(result, Err("Non square matrices do not have determinants!"));
    }

    #[test]
    fn test_adjoint() {
        let test_matrix = Matrix::new(3, 3, [3.0, 4.0, 3.0, 5.0, 7.0, 2.0, 0.0, 0.0, 1.0].to_vec());

        assert_eq!(
            test_matrix
                .adjoint()
                .expect("test_adjoint: unable to calculate matrix adjoint"),
            Matrix::new(
                3,
                3,
                [
                    7.0,
                    -4.0,
                    -13.0,
                    -5.0,
                    3.0,
                    9.0,
                    0.0,
                    0.0,
                    1.0000000000000018
                ]
                .to_vec()
            )
        );
    }

    #[test]
    fn test_inverse() {
        let test_matrix = Matrix::new(
            3,
            3,
            [2.0, 3.0, 0.0, 0.0, 3.0, -3.0, -2.0, 3.0, 3.0].to_vec(),
        );

        let inverse = test_matrix
            .inverse()
            .expect("Unable to compute matrix inverse in test_inverse");

        assert_eq!(inverse.m, 3);
        assert_eq!(inverse.n, 3);

        let expected = [
            1.0 / 3.0,
            -1.0 / 6.0,
            -1.0 / 6.0,
            1.0 / 9.0,
            1.0 / 9.0,
            1.0 / 9.0,
            1.0 / 9.0,
            -2.0 / 9.0,
            1.0 / 9.0,
        ]
        .to_vec();

        for i in 0..expected.len() {
            assert!(expected[i] - inverse.entries[i] < 1.0 / 1000000000000000.0);
        }

        let test_matrix = Matrix::new(
            5,
            5,
            [
                1.0, 2.0, 1.0, 2.0, 3.0, 2.0, 3.0, 1.0, 0.0, 1.0, 2.0, 2.0, 1.0, 0.0, 0.0, 1.0,
                1.0, 1.0, 1.0, 1.0, 0.0, -2.0, 0.0, 2.0, -2.0,
            ]
            .to_vec(),
        );

        assert_eq!(
            test_matrix
                .inverse()
                .expect("Error computing inverse of test case 2 in test_inverse"),
            Matrix::new(
                5,
                5,
                [
                    2.0, -4.0, 5.0, -3.0, -0.5, -1.0, 3.0, -3.0, 1.0, 0.5, -2.0, 2.0, -3.0, 4.0,
                    0.0, 0.0, 1.0, -1.0, 0.0, 0.5, 1.0, -2.0, 2.0, -1.0, -0.5
                ]
                .to_vec()
            )
        );
    }

    #[test]
    fn test_inverse_err_case() {
        let mut rng = rand::thread_rng();
        let random_m = rng.gen_range(0..100);
        let mut random_n: usize = random_m.clone();

        while random_n == random_m {
            random_n = rng.gen_range(0..100);
        }

        let test_matrix = new_all_default::<f64>(random_m, random_n);

        assert_eq!(
            test_matrix
                .expect("Unable to construct new_all_default matrix in test_inverse_err_case")
                .inverse(),
            Err("Unable to compute inverse when n is not equal to m")
        );
    }
}