diff --git a/Cargo.lock b/Cargo.lock index ef6361d..bf40cdb 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -2,6 +2,12 @@ # It is not intended for manual editing. version = 4 +[[package]] +name = "jedvek" +version = "0.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2994ad4a94f739dff0d4eb50dc5be6a8ba738f1d0dd44bf0e553a5fb359c700f" + [[package]] name = "proc-macro2" version = "1.0.105" @@ -24,6 +30,7 @@ dependencies = [ name = "spindalis" version = "0.6.0" dependencies = [ + "jedvek", "spindalis_core", "spindalis_macros", ] diff --git a/README.md b/README.md index 93696ff..51e68eb 100644 --- a/README.md +++ b/README.md @@ -83,20 +83,6 @@ include multivariate polynomials, functions such as sin or cos, constants such a pi, fractional exponents, negative exponents, and optimistically also polynomial expansion. -### Arr2D - -The Arr2D struct provides a convenient way to use a two-dimensional matrix that is -a one-dimensional vector under the hood. -This allows for efficient memory usage while enabling standard matrix operations. - -Arr2D implements a wide range of methods and traits, including: - -- Linear Algebra: dot product, inverse, transpose. -- Arithmetic: Implementations of multiplication by matrices, vectors, and scalars and division by scalars -- Manipulation: shape, size, full, reshape, map -- Conversion: from_flat, From, TryFrom -- Utility: new, max, min, is_empty - ## Project layout This repository is split into three crates. The tests for spindalis_core and @@ -112,13 +98,12 @@ Within these crates, the following modules `spindalis::` are provid | Module | Description | | --------------- | ------------------------------------------------------------------------------------------------------ | -| `utils` | Utility functions such as `Arr2D`, `Arr2DError`, forward substitution, and back substitution | +| `utils` | Utility functions such as different types of mean and standard deviation | | `polynomials` | Parsing and evaluating simple and intermediate polynomials | | `derivatives` | Differentiating simple and intermediate polynomials | | `integrals` | Integrating simple and intermediate polynomials | | `solvers` | Solving equations and differential equations, including root-finding, extrema-finding, and ODE solvers | | `eigen` | Algorithms to solve eigenvalue and eigenvector problems | -| `decomposition` | Decomposition algorithms including LU decomposition and LU decomposition with partial pivoting | | `regressors` | Linear and non-linear regression, including least-squares, Gaussian, and polynomial regression | | `reduction` | Linear and non-linear dimensionality reduction algorithms, including PCA | @@ -149,7 +134,11 @@ Stay connected via our **[Discord Server](https://discord.gg/PdVZCtcgaH)** ## Stability -This project is in the alpha stage. APIs may change without warning until version +> [!CAUTION] +> Arr2D functionality and decomposition functions moved to [jedvek](https://github.com/lignum-vitae/jedvek)! +> Arr2D is now Matrix2D. + +This project is in the beta stage. APIs may change without warning until version 1.0.0. ## License diff --git a/spindalis/Cargo.toml b/spindalis/Cargo.toml index 1e54dcd..b4da6c3 100644 --- a/spindalis/Cargo.toml +++ b/spindalis/Cargo.toml @@ -15,5 +15,6 @@ name = "spindalis" path = "src/lib.rs" [dependencies] +jedvek = "0.1.0" spindalis_core = { path = "../spindalis_core", version = ">=0.2.0" } spindalis_macros = { path = "../spindalis_macros", version = "0.1.0" } diff --git a/spindalis/examples/README.md b/spindalis/examples/README.md index c54be09..f02a5d1 100644 --- a/spindalis/examples/README.md +++ b/spindalis/examples/README.md @@ -115,8 +115,8 @@ parsing to parse polynomials into an abstract syntax tree. ## Math Where applicable, functions accept any coefficient matrix that can be coerced -into a `Arr2D` type. That includes nested vecs of ints or floats, -nested arrays of ints or floats, and Arr2D vectors of types other than +into a `Matrix2D` type. That includes nested vecs of ints or floats, +nested arrays of ints or floats, and Matrix2D vectors of types other than `f64`. The matrix must be borrowed for conversion to work. The right hand side vector also accepts a vector containing any numerical values that can be converted into `f64`. diff --git a/spindalis/examples/eigen.rs b/spindalis/examples/eigen.rs index 1f337e1..202e14a 100644 --- a/spindalis/examples/eigen.rs +++ b/spindalis/examples/eigen.rs @@ -1,8 +1,8 @@ +use jedvek::{Matrix2D, Rounding}; use spindalis::eigen::power_method; -use spindalis::utils::{Arr2D, Rounding}; fn main() { - let matrix = Arr2D::from(&[[2.0, 8.0, 10.0], [8.0, 4.0, 5.0], [10.0, 5.0, 7.0]]); + let matrix = Matrix2D::from(&[[2.0, 8.0, 10.0], [8.0, 4.0, 5.0], [10.0, 5.0, 7.0]]); println!("Original Matrix:\n{matrix}\n"); let result = power_method(&matrix, 1e-10); match result { @@ -16,7 +16,7 @@ fn main() { Ok(inverted) => inverted, Err(e) => { eprintln!("Unable to invert matrix: {e:?}"); - Arr2D::from(&[[]]) + Matrix2D::from(&[[]]) } }; if !inverse_matrix.is_empty() { diff --git a/spindalis/examples/lu_decomposition.rs b/spindalis/examples/lu_decomposition.rs deleted file mode 100644 index 3ce30a6..0000000 --- a/spindalis/examples/lu_decomposition.rs +++ /dev/null @@ -1,17 +0,0 @@ -use spindalis::decomposition::{lu_decomposition, lu_pivot_decomposition}; -use spindalis::utils::Arr2D; - -fn main() { - let mat = Arr2D::from(&[[2, -1, -2], [-4, 6, 3], [-4, -2, 8]]); - // LU Decomposition - let (lower, upper) = lu_decomposition(&mat).unwrap(); - println!("LU Decomposition"); - println!("Lower Triangle Matrix:\n{lower}\nUpper Triange Matrix:\n{upper}\n"); - - // PLU Decomposition - let (lower, upper, permutation) = lu_pivot_decomposition(&mat).unwrap(); - println!("LU Decomposition with partial pivoting (PLU decomposition)"); - println!( - "Lower Triange Matrix:\n{lower}\nUpper Triange Matrix:\n{upper}\nPermutation matrix:\n{permutation}" - ); -} diff --git a/spindalis/src/lib.rs b/spindalis/src/lib.rs index 566c601..7d47185 100644 --- a/spindalis/src/lib.rs +++ b/spindalis/src/lib.rs @@ -46,11 +46,6 @@ pub mod integrals { pub use spindalis_core::integrals::univariate_definite::romberg_definite; } -pub mod decomposition { - pub use crate::solvers::decomposition::lu::lu_decomposition; - pub use crate::solvers::decomposition::plu::lu_pivot_decomposition; -} - pub mod eigen { pub use crate::solvers::eigen::power_method::power_method; } diff --git a/spindalis/src/reduction/dimension/linear/pca.rs b/spindalis/src/reduction/dimension/linear/pca.rs index fc9159e..6a9168e 100644 --- a/spindalis/src/reduction/dimension/linear/pca.rs +++ b/spindalis/src/reduction/dimension/linear/pca.rs @@ -1,5 +1,6 @@ use crate::reduction::dimension::{DimensionError, ReductionError}; -use crate::utils::{StdDevType, arith_mean, arr2D::Arr2D, std_dev}; +use crate::utils::{StdDevType, arith_mean, std_dev}; +use jedvek::Matrix2D; // ┌─────────────┬────────┬────────┬────────┬────────┬──────────┐ // │ │ USA │ France │ Belgium│ UK │ Czechia │ @@ -16,9 +17,9 @@ use crate::utils::{StdDevType, arith_mean, arr2D::Arr2D, std_dev}; pub fn pca() {} fn _center_data( - data: &Arr2D, + data: &Matrix2D, std_type: Option, -) -> Result, ReductionError> { +) -> Result, ReductionError> { let mut result = data.clone(); let mut std = 1_f64; for row in &mut result { @@ -76,10 +77,10 @@ fn _covariance(x_data: &[f64], y_data: &[f64]) -> Result { Ok(cov_sum / (n - 1.0)) } -fn _cov_mat(data: &Arr2D) -> Result, ReductionError> { +fn _cov_mat(data: &Matrix2D) -> Result, ReductionError> { let height = data.height; - let mut covariance_matrix = Arr2D::full(0.0, height, height); + let mut covariance_matrix = Matrix2D::full(0.0, height, height); for (i, x) in data.into_iter().enumerate() { for (j, y) in data.into_iter().enumerate() { @@ -99,10 +100,10 @@ mod tests { #[test] fn test_center_data() { - let data = Arr2D::from(&[[1., 2., 3.], [4., 5., 6.]]); + let data = Matrix2D::from(&[[1., 2., 3.], [4., 5., 6.]]); let centered = _center_data(&data, None).unwrap(); - let expected = Arr2D::from(&[[-1., 0., 1.], [-1., 0., 1.]]); + let expected = Matrix2D::from(&[[-1., 0., 1.], [-1., 0., 1.]]); assert_eq!(centered, expected); } diff --git a/spindalis/src/reduction/dimension/mod.rs b/spindalis/src/reduction/dimension/mod.rs index cc6889e..94ddbc8 100644 --- a/spindalis/src/reduction/dimension/mod.rs +++ b/spindalis/src/reduction/dimension/mod.rs @@ -1,18 +1,18 @@ pub mod linear; pub mod non_linear; -use crate::utils::Arr2DError; +use jedvek::Matrix2DError; pub use linear::pca::pca; #[derive(Debug)] pub enum ReductionError { ShapeError(DimensionError), - InvalidFlatVector(Arr2DError), + InvalidFlatVector(Matrix2DError), ZeroMean, } -impl From for ReductionError { - fn from(err: Arr2DError) -> Self { +impl From for ReductionError { + fn from(err: Matrix2DError) -> Self { ReductionError::InvalidFlatVector(err) } } diff --git a/spindalis/src/reduction/matrix/hessenberg.rs b/spindalis/src/reduction/matrix/hessenberg.rs index 7fd12b0..085597d 100644 --- a/spindalis/src/reduction/matrix/hessenberg.rs +++ b/spindalis/src/reduction/matrix/hessenberg.rs @@ -1,19 +1,21 @@ use crate::solvers::SolverError; -use crate::utils::Arr2D; +use jedvek::Matrix2D; /// Computes the Hessenberg reduction of a square matrix A. /// Returns (H, Q) such that H = Q^T * A * Q, where H is upper Hessenberg and Q is orthogonal. -pub fn hessenberg_reduction(matrix: &Arr2D) -> Result<(Arr2D, Arr2D), SolverError> { +pub fn hessenberg_reduction( + matrix: &Matrix2D, +) -> Result<(Matrix2D, Matrix2D), SolverError> { if matrix.height != matrix.width { return Err(SolverError::NonSquareMatrix); } let n = matrix.height; if n <= 2 { - return Ok((matrix.clone(), Arr2D::identity(n))); + return Ok((matrix.clone(), Matrix2D::identity(n))); } let mut h = matrix.clone(); - let mut q = Arr2D::identity(n); + let mut q = Matrix2D::identity(n); for k in 0..n - 2 { // x = h[k+1..n, k] @@ -83,23 +85,23 @@ pub fn hessenberg_reduction(matrix: &Arr2D) -> Result<(Arr2D, Arr2D(matrix: M) -> Result<(Arr2D, Arr2D), SolverError> -where - M: TryInto, Error = Arr2DError>, -{ - let matrix: Arr2D = matrix.try_into()?; - if matrix.height != matrix.width { - return Err(SolverError::NonSquareMatrix); - } - - let mut lower: Arr2D = Arr2D::full(0.0, matrix.height, matrix.width); - let mut upper: Arr2D = Arr2D::full(0.0, matrix.height, matrix.width); - - for i in 0..matrix.height { - for k in i..matrix.height { - let mut total = 0.0; - for j in 0..i { - total += lower[i][j] * upper[j][k]; - } - upper[i][k] = matrix[i][k] - total; - } - for k in i..matrix.height { - if i == k { - lower[i][i] = 1_f64; - } else { - let mut total = 0.0; - for j in 0..i { - total += lower[k][j] * upper[j][i]; - } - lower[k][i] = (matrix[k][i] - total) / upper[i][i]; - } - } - } - - Ok((lower, upper)) -} - -#[cfg(test)] -mod tests { - use super::*; - use crate::solvers::SolverError; - use crate::utils::Arr2D; - - #[test] - fn test_known_solution() { - let mat = Arr2D::from(&[[2, -1, -2], [-4, 6, 3], [-4, -2, 8]]); - let result = lu_decomposition(&mat).unwrap(); - let lower_exp = Arr2D::from(&[[1.0, 0.0, 0.0], [-2.0, 1.0, 0.0], [-2.0, -1.0, 1.0]]); - let upper_exp = Arr2D::from(&[[2.0, -1.0, -2.0], [0.0, 4.0, -1.0], [0.0, 0.0, 3.0]]); - let expected = (lower_exp, upper_exp); - - assert_eq!(result, expected); - } - - #[test] - fn test_known_solution_nested_vectors() { - let mat = vec![ - vec![2.0, -1.0, -2.0], - vec![-4.0, 6.0, 3.0], - vec![-4.0, -2.0, 8.0], - ]; - let result = lu_decomposition(mat).unwrap(); - let lower_exp = Arr2D::from(&[[1.0, 0.0, 0.0], [-2.0, 1.0, 0.0], [-2.0, -1.0, 1.0]]); - let upper_exp = Arr2D::from(&[[2.0, -1.0, -2.0], [0.0, 4.0, -1.0], [0.0, 0.0, 3.0]]); - let expected = (lower_exp, upper_exp); - - assert_eq!(result, expected); - } - - #[test] - fn test_known_solution_nested_borrowed_vectors() { - let mat = vec![vec![2, -1, -2], vec![-4, 6, 3], vec![-4, -2, 8]]; - let result = lu_decomposition(&mat).unwrap(); - let lower_exp = Arr2D::from(&[[1.0, 0.0, 0.0], [-2.0, 1.0, 0.0], [-2.0, -1.0, 1.0]]); - let upper_exp = Arr2D::from(&[[2.0, -1.0, -2.0], [0.0, 4.0, -1.0], [0.0, 0.0, 3.0]]); - let expected = (lower_exp, upper_exp); - - assert_eq!(result, expected); - } - - #[test] - fn test_non_square_matrix() { - let mat = Arr2D::from(&[[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]); // 2×3 - let result = lu_decomposition(&mat); - - assert!(matches!(result, Err(SolverError::NonSquareMatrix))); - } -} diff --git a/spindalis/src/solvers/decomposition/mod.rs b/spindalis/src/solvers/decomposition/mod.rs deleted file mode 100644 index 4d0ed11..0000000 --- a/spindalis/src/solvers/decomposition/mod.rs +++ /dev/null @@ -1,2 +0,0 @@ -pub mod lu; -pub mod plu; diff --git a/spindalis/src/solvers/decomposition/plu.rs b/spindalis/src/solvers/decomposition/plu.rs deleted file mode 100644 index 9332243..0000000 --- a/spindalis/src/solvers/decomposition/plu.rs +++ /dev/null @@ -1,129 +0,0 @@ -// Source for implementation https://www.math.hkust.edu.hk/~machas/numerical-methods-for-engineers.pdf -// Companion video on youtube https://www.youtube.com/watch?v=WgaFycuL8z0 -// In my implementation, the L and U matrices are computed simultaneously -// instead of calculating the U matrix then the L matrix as in the video -// L matrix is composed of multipliers for elimination (as stated at 6:55) -use crate::solvers::SolverError; -use crate::utils::{Arr2D, Arr2DError}; - -pub type PLUResult = (Arr2D, Arr2D, Arr2D); - -// LU Decomposition with Partial pivoting -pub fn lu_pivot_decomposition(matrix: M) -> Result -where - M: TryInto, Error = Arr2DError>, -{ - let matrix: Arr2D = matrix.try_into()?; - if matrix.height != matrix.width { - return Err(SolverError::NonSquareMatrix); - } - - let size = matrix.height; - let mut lu = matrix.clone(); - let mut permutation: Arr2D = Arr2D::identity(size); - - for i in 0..size { - // Find the best pivot row (p) for the current column (i) - let mut pivot_row = i; - let mut max_value = lu[i][i].abs(); - - // Select best pivot row based on highest scaled pivot ratio in column - for k in (i + 1)..size { - let value = lu[k][i].abs(); - if value > max_value { - max_value = value; - pivot_row = k; - } - } - - // Max value placed on diagonal (column and row index for swap are the same) - if pivot_row != i { - lu.swap_rows(pivot_row, i); - permutation.swap_rows(pivot_row, i); - } - // Check for singularity BEFORE division - if lu[i][i].abs() < f64::EPSILON { - return Err(SolverError::SingularMatrix); - } - - for k in (i + 1)..size { - // Calculate the lower triangular matrix - // Calculates multiplier needed to eliminate the element lu[k][i] - lu[k][i] /= lu[i][i]; // division by 0 stopped by singularity check - - // Elimination step - // lu[k][j] will form part of the final upper triangular matrix - for j in (i + 1)..size { - lu[k][j] -= lu[k][i] * lu[i][j]; - } - } - } - let mut lower: Arr2D = Arr2D::full(0.0, size, size); - let mut upper: Arr2D = Arr2D::full(0.0, size, size); - - for i in 0..size { - for j in 0..size { - if i == j { - // U takes the diagonal elements from the LU matrix. - upper[i][j] = lu[i][j]; - // L is a unit lower triangular matrix, so its diagonal elements are 1.0. - lower[i][j] = 1.0; - } else if i > j { - // Lower triangle (below diagonal) - lower[i][j] = lu[i][j]; - } else { - // Upper triangle (above diagonal) - upper[i][j] = lu[i][j]; - } - } - } - Ok((lower, upper, permutation)) -} - -#[cfg(test)] -mod tests { - use super::*; - use crate::utils::Arr2D; - - #[test] - fn test_known_solution() { - let matrix = Arr2D::from(&[[0, 1, -2], [1, 0, 2], [3, -2, 2]]); - let perm_exp = Arr2D::from(&[[0.0, 0.0, 1.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]]); - let lower_exp = Arr2D::from(&[ - [1.0, 0.0, 0.0], - [0.0, 1.0, 0.0], - [1.0 / 3.0, 2.0 / 3.0, 1.0], - ]); - let upper_exp = Arr2D::from(&[ - [3.0, -2.0, 2.0], - [0.0, 1.0, -2.0], - [0.0, 0.0, 2.666666666666667], - ]); - - let expected = (lower_exp, upper_exp, perm_exp); - - let result = lu_pivot_decomposition(&matrix).unwrap(); - assert_eq!(result, expected); - } - - #[test] - fn test_known_solution_2() { - let matrix = Arr2D::from(&[[-2, 2, -1], [6, -6, 7], [3, -8, 4]]); - let perm_exp = Arr2D::from(&[[0.0, 1.0, 0.0], [0.0, 0.0, 1.0], [1.0, 0.0, 0.0]]); - let lower_exp = Arr2D::from(&[ - [1.0, 0.0, 0.0], - [1.0 / 2.0, 1.0, 0.0], - [-1.0 / 3.0, -0.0, 1.0], // result returns -0.0 instead of 0.0 - ]); - let upper_exp = Arr2D::from(&[ - [6.0, -6.0, 7.0], - [0.0, -5.0, 1.0 / 2.0], - [0.0, 0.0, 1.333333333333333], // result is short by one decimal place compared to 4/3 - ]); - - let expected = (lower_exp, upper_exp, perm_exp); - - let result = lu_pivot_decomposition(&matrix).unwrap(); - assert_eq!(result, expected); - } -} diff --git a/spindalis/src/solvers/eigen/power_method.rs b/spindalis/src/solvers/eigen/power_method.rs index 9c45645..56d519a 100644 --- a/spindalis/src/solvers/eigen/power_method.rs +++ b/spindalis/src/solvers/eigen/power_method.rs @@ -1,16 +1,16 @@ -use crate::utils::{Arr2D, Arr2DError}; +use jedvek::{Matrix2D, Matrix2DError}; -pub fn power_method(matrix: M, es: f64) -> Result<(f64, Arr2D), Arr2DError> +pub fn power_method(matrix: M, es: f64) -> Result<(f64, Matrix2D), Matrix2DError> where - M: TryInto, Error = Arr2DError>, + M: TryInto, Error = Matrix2DError>, { - let matrix: Arr2D = matrix.try_into()?; + let matrix: Matrix2D = matrix.try_into()?; if matrix.height != matrix.width || matrix.height == 0 || matrix.width == 0 { - return Err(Arr2DError::NonSquareMatrix); + return Err(Matrix2DError::NonSquareMatrix); } - let initial_eigenvector = Arr2D::from(&[[1.0], [1.0], [1.0]]); + let initial_eigenvector = Matrix2D::from(&[[1.0], [1.0], [1.0]]); let mut eigenvector = &matrix * initial_eigenvector; - // Arr2D.max() only returns None if the matrix is empty + // Matrix2D.max() only returns None if the matrix is empty let mut eigenvalue = eigenvector.max().unwrap(); // Matrix won't be empty here eigenvector = eigenvector / eigenvalue; // Normalised Eigenvector loop { @@ -37,17 +37,17 @@ where #[cfg(test)] mod tests { use super::*; - use crate::utils::arr2D::Rounding; + use jedvek::Rounding; #[test] fn test_known_solution() { // Largest eigenvalue and eigenvector - let matrix = Arr2D::from(&[ + let matrix = Matrix2D::from(&[ [3.556, -1.778, 0.0], [-1.778, 3.556, -1.778], [0.0, -1.778, 3.556], ]); - let expected = (6.070, Arr2D::from(&[[1.0], [-1.414], [1.0]])); + let expected = (6.070, Matrix2D::from(&[[1.0], [-1.414], [1.0]])); let (mut eigenvalue, mut eigenvector) = power_method(&matrix, 1e-10).unwrap(); eigenvalue = (eigenvalue * 10_f64.powi(2)).round() / 10_f64.powi(2); @@ -59,7 +59,7 @@ mod tests { #[test] fn test_known_inverse_solution() { // Smallest eigenvalue and eigenvector - let matrix = Arr2D::from(&[ + let matrix = Matrix2D::from(&[ [3.556, -1.778, 0.0], [-1.778, 3.556, -1.778], [0.0, -1.778, 3.556], @@ -69,7 +69,7 @@ mod tests { ((1_f64 / 0.960127) * 10_f64.powi(5)).round() / 10_f64.powi(5); let expected = ( rounded_expected_eigenvalue, - Arr2D::from(&[[0.707], [1.0], [0.707]]), + Matrix2D::from(&[[0.707], [1.0], [0.707]]), ); let (converged_value, mut eigenvector) = power_method(&inverse_res, 1e-10).unwrap(); @@ -83,8 +83,8 @@ mod tests { #[test] fn test_known_solution_2() { // Largest eigenvalue and eigenvector - let matrix = Arr2D::from(&[[2, 8, 10], [8, 4, 5], [10, 5, 7]]); - let expected = (19.88, Arr2D::from(&[[0.9035], [0.7698], [1.0]])); + let matrix = Matrix2D::from(&[[2, 8, 10], [8, 4, 5], [10, 5, 7]]); + let expected = (19.88, Matrix2D::from(&[[0.9035], [0.7698], [1.0]])); let (mut eigenvalue, mut eigenvector) = power_method(&matrix, 1e-10).unwrap(); eigenvalue = (eigenvalue * 10_f64.powi(2)).round() / 10_f64.powi(2); @@ -96,8 +96,8 @@ mod tests { #[test] fn test_known_inverse_solution_2() { // Smallest eigenvalue and eigenvector - let matrix = Arr2D::from(&[[2, 8, 10], [8, 4, 5], [10, 5, 7]]); - let expected = (0.29, Arr2D::from(&[[0.04117], [1.0], [-0.80702]])); + let matrix = Matrix2D::from(&[[2, 8, 10], [8, 4, 5], [10, 5, 7]]); + let expected = (0.29, Matrix2D::from(&[[0.04117], [1.0], [-0.80702]])); let inverse_res = matrix.inverse().unwrap(); let (converged_value, mut eigenvector) = power_method(&inverse_res, 1e-10).unwrap(); diff --git a/spindalis/src/solvers/gaussian_elim.rs b/spindalis/src/solvers/gaussian_elim.rs index 797f0a3..fe09f99 100644 --- a/spindalis/src/solvers/gaussian_elim.rs +++ b/spindalis/src/solvers/gaussian_elim.rs @@ -1,5 +1,5 @@ use crate::solvers::SolverError; -use crate::utils::{Arr2DError, arr2D::Arr2D, back_substitution}; +use jedvek::{Matrix2D, Matrix2DError, substitution::back_substitution}; pub fn gaussian_elimination( matrix: M, @@ -7,10 +7,10 @@ pub fn gaussian_elimination( tolerance: f64, ) -> Result, SolverError> where - M: TryInto, Error = Arr2DError>, + M: TryInto, Error = Matrix2DError>, W: Into + Copy, { - let mut coeff_matrix: Arr2D = matrix.try_into()?; + let mut coeff_matrix: Matrix2D = matrix.try_into()?; if coeff_matrix.height != coeff_matrix.width { return Err(SolverError::NonSquareMatrix); } @@ -51,7 +51,7 @@ where } fn forward_elimination( - coeff_matrix: &mut Arr2D, + coeff_matrix: &mut Matrix2D, scale_factor: &mut [f64], size: usize, rhs_vector: &mut [f64], @@ -78,7 +78,7 @@ fn forward_elimination( } fn partial_pivot( - coeff_matrix: &mut Arr2D, + coeff_matrix: &mut Matrix2D, rhs_vector: &mut [f64], scale_factor: &mut [f64], size: usize, @@ -130,7 +130,7 @@ mod tests { fn test_known_elim_2() { let matrix = vec![vec![3, 6], vec![5, -8]]; - let coeff_matrix = Arr2D::try_from(matrix).unwrap(); + let coeff_matrix = Matrix2D::try_from(matrix).unwrap(); let rhs_vector = vec![12, 2]; let tol = 1e-12; let expected: Vec = vec![2.0, 1.0]; @@ -148,7 +148,7 @@ mod tests { vec![3.0, -5.0, 5.0, -4.0], ]; - let coeff_matrix: Arr2D = Arr2D::try_from(matrix).unwrap(); + let coeff_matrix: Matrix2D = Matrix2D::try_from(matrix).unwrap(); let rhs_vector = vec![0.0; 4]; let tol = 1e-12; diff --git a/spindalis/src/solvers/mod.rs b/spindalis/src/solvers/mod.rs index 38c6b7e..0d13bd3 100644 --- a/spindalis/src/solvers/mod.rs +++ b/spindalis/src/solvers/mod.rs @@ -1,5 +1,4 @@ pub mod bisection; -pub mod decomposition; pub mod eigen; pub mod gaussian_elim; pub mod nrm; @@ -9,7 +8,7 @@ pub use gaussian_elim::gaussian_elimination; pub use nrm::newton_raphson_method; use crate::polynomials::PolynomialError; -use crate::utils::Arr2DError; +use jedvek::Matrix2DError; #[derive(PartialEq)] pub enum SolveMode { @@ -24,13 +23,13 @@ pub enum SolverError { XInitOutOfBounds, NonSquareMatrix, SingularMatrix, - InvalidVector(Arr2DError), + InvalidVector(Matrix2DError), FunctionError(PolynomialError), NumArgumentsMismatch { num_rows: usize, rhs_len: usize }, } -impl From for SolverError { - fn from(err: Arr2DError) -> Self { +impl From for SolverError { + fn from(err: Matrix2DError) -> Self { SolverError::InvalidVector(err) } } diff --git a/spindalis/src/utils/arr2D.rs b/spindalis/src/utils/arr2D.rs deleted file mode 100644 index 90aca46..0000000 --- a/spindalis/src/utils/arr2D.rs +++ /dev/null @@ -1,1281 +0,0 @@ -use crate::decomposition::lu_pivot_decomposition; -use crate::utils::{Arr2DError, back_substitution, forward_substitution}; -use std::{ - any::type_name, - fmt::{self, Display}, - ops::{Div, Index, IndexMut, Mul}, -}; - -#[derive(Debug, PartialEq, Eq, Clone, Hash, Default)] -pub struct Arr2D { - inner: Vec, - pub height: usize, - pub width: usize, -} - -impl Arr2D { - pub fn new() -> Self { - Arr2D { - inner: Vec::new(), - height: 0, - width: 0, - } - } - - /// Return a tuple of (height, width) - pub fn shape(&self) -> (usize, usize) { - (self.height, self.width) - } - - /// Return the total number of items (i.e. height * width) - pub fn size(&self) -> usize { - self.inner.len() - } - - pub fn is_empty(&self) -> bool { - self.height == 0 || self.width == 0 - } - - pub fn max(&self) -> Option - where - T: std::cmp::PartialOrd + Clone, - { - if self.is_empty() { - return None; - } - Some( - self.inner - .iter() - .reduce(|a, b| if a > b { a } else { b }) - .unwrap() - .clone(), - ) - } - - pub fn min(&self) -> Option - where - T: std::cmp::PartialOrd + Clone, - { - if self.is_empty() { - return None; - } - Some( - self.inner - .iter() - .reduce(|a, b| if a < b { a } else { b }) - .unwrap() - .clone(), - ) - } - - pub fn inverse(&self) -> Result, Arr2DError> - where - Arr2D: for<'a> TryFrom<&'a Arr2D>, - { - if self.height != self.width { - return Err(Arr2DError::NonSquareMatrix); - } - let coeff_matrix: &Arr2D = - &self.try_into().map_err(|_| Arr2DError::ConversionFailed { - from: type_name::(), - to: type_name::(), - })?; - let size = self.height; - - let (l, u, p) = - lu_pivot_decomposition(coeff_matrix).map_err(|_| Arr2DError::SingularMatrix)?; - - let mut inverse_matrix = Arr2D::full(0.0, size, size); - - for j in 0..size { - let mut b_prime = vec![0.0; size]; - for i in 0..size { - // b' = P * e_j is the j-th column of P - b_prime[i] = p[i][j]; - } - // L * y = b' -> solve for y - let mut y = vec![0.0; size]; - forward_substitution(&l, size, &b_prime, &mut y); - - // U * x = y -> solve for x - let mut x_j = vec![0.0; size]; - back_substitution(&u, size, &y, &mut x_j); - - // Solution vector placed into j-th column of inverse matrix - for i in 0..size { - inverse_matrix[i][j] = x_j[i]; - } - } - Ok(inverse_matrix) - } - - /// Change the height and width - /// - /// # Errors - /// - /// This function will return an error if the size isn't divislbe by the new height - pub fn reshape(&mut self, height: usize) -> Result<(), Arr2DError> { - let size = self.height * self.width; - if !size.is_multiple_of(height) { - return Err(Arr2DError::InvalidReshape { - size, - new_height: height, - }); - } - - self.height = height; - self.width = size / height; - Ok(()) - } - - /// Element-wise map operatiorn - pub fn map(&self, f: F) -> Arr2D - where - F: Fn(&T) -> U, - { - let inner: Vec<_> = self.inner.iter().map(f).collect(); - - Arr2D { - inner, - height: self.height, - width: self.width, - } - } - - /// Create an iterator of refs to rows - pub fn rows(&self) -> Arr2DRows<'_, T> { - Arr2DRows { - data: &self.inner, - width: self.width, - remaining: self.height, - } - } - - /// Create an iterator of mut refs to rows - pub fn rows_mut(&mut self) -> Arr2DRowsMut<'_, T> { - Arr2DRowsMut { - data: self.inner.as_mut_slice(), - width: self.width, - remaining: self.height, - } - } - - // Create 2D Array from flat vector - pub fn from_flat( - inner: D, - default_val: T, - height: usize, - width: usize, - ) -> Result - where - D: AsRef<[T]>, - T: Clone, - { - let vec_len = inner.as_ref().len(); - let Arr2D_size = height * width; - if vec_len > Arr2D_size || Arr2D_size == 0 { - return Err(Arr2DError::InvalidShape { - input_size: (vec_len), - output_size: (Arr2D_size), - }); - } - - let inner = inner.as_ref().to_vec(); - if vec_len < Arr2D_size { - let mut new_inner = inner.clone(); - new_inner.resize(Arr2D_size, default_val); - return Ok(Self { - inner: new_inner, - height, - width, - }); - } - - Ok(Self { - inner, - height, - width, - }) - } - - // Dot product for scalar x matrix, vector x matrix, and matrix x matrix - pub fn dot(&self, rhs: &Self) -> Result - where - T: Copy + std::default::Default + std::ops::AddAssign + std::ops::Mul, - { - if self.height == 1 && self.width == 1 || rhs.height == 1 && rhs.width == 2 { - let mut matrix = Arr2D::new(); - let mut scalar = T::default(); - if rhs.height == 1 { - matrix = self.clone(); - scalar = rhs[0][0]; - } else if self.height == 1 { - matrix = rhs.clone(); - scalar = self[0][0]; - } - - let mut result = Arr2D::full(T::default(), matrix.height, matrix.width); - for i in 0..matrix.height { - for j in 0..matrix.width { - result[i][j] = scalar * matrix[i][j]; - } - } - return Ok(result); - } - if self.width != rhs.height { - return Err(Arr2DError::InvalidDotShape { - lhs: self.width, - rhs: rhs.height, - }); - } - let mut result = Arr2D::full(T::default(), self.height, rhs.width); - for i in 0..self.height { - for j in 0..rhs.width { - let mut sum = T::default(); - for k in 0..self.width { - sum += self[i][k] * rhs[k][j] - } - result[i][j] = sum; - } - } - Ok(result) - } -} - -// Mul implementation for &Arr2D * &Arr2D (&matrix * &matrix) -impl<'b, T> Mul<&'b Arr2D> for &Arr2D -where - T: Mul + Clone + std::default::Default + std::marker::Copy + std::ops::AddAssign, -{ - type Output = Arr2D; - - fn mul(self, rhs: &'b Arr2D) -> Arr2D { - self.dot(rhs).unwrap_or_default() - } -} - -// Mul implementation for Arr2D * Arr2D (matrix * matrix) -impl Mul> for Arr2D -where - T: Mul + Clone + std::default::Default + std::marker::Copy + std::ops::AddAssign, -{ - type Output = Arr2D; - - fn mul(self, rhs: Arr2D) -> Arr2D { - self.dot(&rhs).unwrap_or_default() - } -} - -// Mul implementation for Arr2D * &Arr2D (matrix * &matrix) -impl<'b, T> Mul<&'b Arr2D> for Arr2D -where - T: Mul + Clone + std::default::Default + std::marker::Copy + std::ops::AddAssign, -{ - type Output = Arr2D; - - fn mul(self, rhs: &'b Arr2D) -> Arr2D { - self.dot(rhs).unwrap_or_default() - } -} - -// Mul implementation for &Arr2D * Arr2D (&matrix * matrix) -impl Mul> for &Arr2D -where - T: Mul + Clone + std::default::Default + std::marker::Copy + std::ops::AddAssign, -{ - type Output = Arr2D; - - fn mul(self, rhs: Arr2D) -> Arr2D { - self.dot(&rhs).unwrap_or_default() - } -} - -// Mul implementation for Arr2D * scalar (matrix * scalar) -impl Mul for Arr2D -where - T: Mul + Clone + std::default::Default + std::marker::Copy, -{ - type Output = Arr2D; - - fn mul(self, rhs: T) -> Arr2D { - &self * rhs - } -} - -// Mul implementation for &Arr2D * scalar (&matrix * scalar) -impl Mul for &Arr2D -where - T: Mul + Clone + std::default::Default + std::marker::Copy, -{ - type Output = Arr2D; - - fn mul(self, rhs: T) -> Arr2D { - let mut result = Arr2D::full(T::default(), self.height, self.width); - for i in 0..self.height { - for j in 0..self.width { - result[i][j] = self[i][j] * rhs; - } - } - result - } -} - -// Div implementation for &Arr2D / scalar (&matrix / scalar) -impl Div for &Arr2D -where - T: Div + Clone + std::default::Default + std::marker::Copy, -{ - type Output = Arr2D; - - fn div(self, rhs: T) -> Arr2D { - let mut result = Arr2D::full(T::default(), self.height, self.width); - for i in 0..self.height { - for j in 0..self.width { - result[i][j] = self[i][j] / rhs; - } - } - result - } -} - -// Div implementation for Arr2D / scalar (matrix / scalar) -impl Div for Arr2D -where - T: Div + Clone + std::default::Default + std::marker::Copy, -{ - type Output = Arr2D; - - fn div(self, rhs: T) -> Self { - &self / rhs - } -} - -impl Arr2D { - pub fn as_scalar(&self) -> Option { - if self.height == 1 && self.width == 1 { - Some(self.inner[0]) - } else { - None - } - } - - pub fn as_scalar_unchecked(&self) -> T { - self.inner[0] - } - - pub fn transpose(&self) -> Arr2D { - let mut new_inner = Vec::with_capacity(self.inner.len()); - for col in 0..self.width { - for row in 0..self.height { - new_inner.push(self[(row, col)]); - } - } - Arr2D { - inner: new_inner, - height: self.width, - width: self.height, - } - } - pub fn transpose_mut(&mut self) { - let mut new_inner = Vec::with_capacity(self.inner.len()); - for col in 0..self.width { - for row in 0..self.height { - new_inner.push(self[(row, col)]); - } - } - self.inner = new_inner; - std::mem::swap(&mut self.width, &mut self.height); - } - - pub fn swap_rows(&mut self, mut a: usize, mut b: usize) { - if a == b { - return; - } - if a > b { - std::mem::swap(&mut a, &mut b); - } - - let w = self.width; - - let (left, right) = self.inner.split_at_mut(b * w); - let row_b = &mut right[..w]; - let row_a = &mut left[a * w..(a + 1) * w]; - - row_a.swap_with_slice(row_b); - } - - /// Create Arr2D with all elements being the specified value - pub fn full(val: T, height: usize, width: usize) -> Self { - Arr2D { - inner: vec![val; height * width], - height, - width, - } - } - - pub fn identity(size: usize) -> Self - where - T: From + Default + Copy, - { - let mut ident_mat: Arr2D = Arr2D::full(T::from(0), size, size); - for i in 0..size { - ident_mat[i][i] = T::from(1); - } - ident_mat - } -} - -/// Iterator for Arr2D -pub struct Arr2DRows<'a, T> { - data: &'a [T], - width: usize, - /// remaining row count - remaining: usize, -} - -impl<'a, T> Iterator for Arr2DRows<'a, T> { - type Item = &'a [T]; - - fn next(&mut self) -> Option { - if self.remaining == 0 { - return None; - } - - self.remaining -= 1; - - if self.width == 0 { - Some(&self.data[..0]) - } else { - let (row, rest) = self.data.split_at(self.width); - self.data = rest; - Some(row) - } - } - - fn size_hint(&self) -> (usize, Option) { - (self.remaining, Some(self.remaining)) - } -} - -/// Mut iterator for Arr2D -pub struct Arr2DRowsMut<'a, T> { - data: &'a mut [T], - width: usize, - remaining: usize, -} - -impl<'a, T> Iterator for Arr2DRowsMut<'a, T> { - type Item = &'a mut [T]; - - fn next(&mut self) -> Option { - if self.remaining == 0 { - return None; - } - - // it's important to take a temporary ownership over the mut slice here - // because `self.data` would be alive only during an invocation of this method - // (which would lead to borrow checker error). - // `std::mem::take` allows you to own the slice, leaving `self.data` with an empty slice. - let (row, rest) = std::mem::take(&mut self.data).split_at_mut(self.width); - - self.data = rest; - self.remaining -= 1; - - Some(row) - } - - fn size_hint(&self) -> (usize, Option) { - (self.remaining, Some(self.remaining)) - } -} - -// Convert a nested Vec to Arr2D -impl TryFrom>> for Arr2D { - type Error = Arr2DError; - - fn try_from(values: Vec>) -> Result { - if values.is_empty() { - return Ok(Self { - inner: vec![], - height: 0, - width: 0, - }); - } - let width = values[0].len(); - let height = values.len(); - let mut inner: Vec = Vec::with_capacity(height * width); - for row in values { - if row.len() != width { - return Err(Arr2DError::InconsistentRowLengths); - } - inner.extend(row); - } - Ok(Self { - inner, - height, - width, - }) - } -} - -// &Vec> -> Arr2D -impl> TryFrom<&Vec>> for Arr2D { - type Error = Arr2DError; - - fn try_from(values: &Vec>) -> Result { - if values.is_empty() { - return Ok(Self { - inner: vec![], - height: 0, - width: 0, - }); - } - - let width = values[0].len(); - let mut inner = Vec::with_capacity(values.len() * width); - - for row in values { - if row.len() != width { - return Err(Arr2DError::InconsistentRowLengths); - } - - for x in row { - inner.push( - U::try_from(x.clone()).map_err(|_| Arr2DError::ConversionFailed { - from: type_name::(), - to: type_name::(), - })?, - ); - } - } - - Ok(Self { - inner, - height: values.len(), - width, - }) - } -} - -// Arr2D -> Arr2D -impl> TryFrom<&Arr2D> for Arr2D { - type Error = Arr2DError; - - fn try_from(arr: &Arr2D) -> Result { - if arr.is_empty() { - return Ok(Self { - inner: vec![], - height: 0, - width: 0, - }); - } - - let mut inner = Vec::with_capacity(arr.inner.len()); - for x in &arr.inner { - inner.push( - U::try_from(x.clone()).map_err(|_| Arr2DError::ConversionFailed { - from: type_name::(), - to: type_name::(), - })?, - ); - } - - Ok(Self { - inner, - height: arr.height, - width: arr.width, - }) - } -} - -impl<'a, T> IntoIterator for &'a Arr2D { - type Item = &'a [T]; - type IntoIter = Arr2DRows<'a, T>; - - fn into_iter(self) -> Self::IntoIter { - self.rows() - } -} - -impl<'a, T> IntoIterator for &'a mut Arr2D { - type Item = &'a mut [T]; - type IntoIter = Arr2DRowsMut<'a, T>; - - fn into_iter(self) -> Self::IntoIter { - self.rows_mut() - } -} - -// Convert a nested array like `&[[1, 2], [3, 4]]` to Arr2D -impl From<&[[T; N]; M]> for Arr2D -where - T: Clone, -{ - fn from(values: &[[T; N]; M]) -> Self { - let mut inner = Vec::with_capacity(M * N); - - for row in values.iter() { - inner.extend_from_slice(row); - } - - Self { - inner, - height: M, - width: N, - } - } -} - -// Allow indexing Arr2D items like `arr[(0, 1)]` -impl Index<(usize, usize)> for Arr2D { - type Output = T; - - fn index(&self, idx: (usize, usize)) -> &Self::Output { - let (row, col) = idx; - if row >= self.height || col >= self.width { - panic!( - "Out of bound index ({row},{col}) into Arr2D of shape ({},{})", - self.height, self.width - ) - } - &self.inner[row * self.width + col] - } -} -impl IndexMut<(usize, usize)> for Arr2D { - fn index_mut(&mut self, idx: (usize, usize)) -> &mut Self::Output { - let (row, col) = idx; - if row >= self.height || col >= self.width { - panic!( - "Out of bound index ({row},{col}) into Arr2D of shape ({},{})", - self.height, self.width - ) - } - &mut self.inner[row * self.width + col] - } -} - -// Allow indexing Arr2D rows like `arr[1]` -impl Index for Arr2D { - type Output = [T]; - - fn index(&self, row: usize) -> &Self::Output { - if row >= self.height { - panic!( - "Out of bound row index {row} into Arr2D of shape ({},{})", - self.height, self.width - ) - } - &self.inner[row * self.width..(row + 1) * self.width] - } -} -impl IndexMut for Arr2D { - fn index_mut(&mut self, row: usize) -> &mut Self::Output { - if row >= self.height { - panic!( - "Out of bound row index {row} into Arr2D of shape ({},{})", - self.height, self.width - ) - } - &mut self.inner[row * self.width..(row + 1) * self.width] - } -} - -impl Display for Arr2D { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - if self.height == 0 || self.width == 0 { - writeln!(f, "[]")?; - return Ok(()); - } - - // Get max width per column - let mut col_widths = vec![0; self.width]; - for c in 0..self.width { - col_widths[c] = (0..self.height) - .map(|r| format!("{}", self[(r, c)]).len()) - .max() - .unwrap_or(0); - } - - // Print - for r in 0..self.height { - if r == 0 { - write!(f, "[[ ")?; - } else { - write!(f, " [ ")?; - } - for c in 0..self.width { - let item = &self[(r, c)]; - write!(f, "{:>width$}", *item, width = col_widths[c])?; - if c + 1 != self.width { - write!(f, ", ")?; - } - } - if r + 1 == self.height { - write!(f, " ]]")?; - } else { - writeln!(f, " ]")?; - } - } - - Ok(()) - } -} - -pub trait Rounding { - fn round_to_decimal(&self, decimals: u32) -> Self; -} - -impl Rounding for Arr2D { - fn round_to_decimal(&self, decimals: u32) -> Arr2D { - let factor = (10.0_f64).powi(decimals as i32); - - self.clone().map(|&val| (val * factor).round() / factor) - } -} - -impl PartialEq>> for Arr2D { - fn eq(&self, other: &Vec>) -> bool { - if self.height != other.len() { - return false; - } - if self.height == 0 { - return true; // both empty - } - - let width = self.width; - - if other.iter().any(|row| row.len() != width) { - return false; - } - - for r in 0..self.height { - for c in 0..self.width { - if self[r][c] != other[r][c] { - return false; - } - } - } - true - } -} - -impl PartialEq> for Vec> { - fn eq(&self, other: &Arr2D) -> bool { - other == self - } -} - -#[cfg(test)] -mod tests { - use super::*; - - // --- basic getters --- - - #[test] - fn test_getting_shape() { - let data = Arr2D::from(&[[1, 2, 3], [6, 5, 4]]); - assert_eq!(data.shape(), (2, 3)); - - let data: Arr2D = Arr2D::from(&[[1, 2]; 0]); - assert_eq!(data.shape(), (0, 2)); - - let data: Arr2D = Arr2D::from(&[[]; 0]); - assert_eq!(data.shape(), (0, 0)); - - let data: Arr2D = Arr2D::from(&[[], []]); - assert_eq!(data.shape(), (2, 0)); - } - - #[test] - fn test_getting_size() { - let data = Arr2D::from(&[[1, 2, 3], [6, 5, 4]]); - assert_eq!(data.size(), 6); - - let data: Arr2D = Arr2D::from(&[[1, 2]; 0]); - assert_eq!(data.size(), 0); - - let data: Arr2D = Arr2D::from(&[[]; 0]); - assert_eq!(data.size(), 0); - - let data: Arr2D = Arr2D::from(&[[], []]); - assert_eq!(data.size(), 0); - } - - // --- indexing --- - - #[test] - fn test_indexing_item() { - let data = Arr2D::from(&[[1, 2, 3], [6, 5, 4]]); - assert_eq!(data[(0, 0)], 1); - assert_eq!(data[(0, 1)], 2); - assert_eq!(data[(0, 2)], 3); - assert_eq!(data[(1, 0)], 6); - assert_eq!(data[(1, 1)], 5); - assert_eq!(data[(1, 2)], 4); - } - - #[test] - fn test_2D_indexing_item() { - let data = Arr2D::from(&[[1, 2, 3], [6, 5, 4]]); - assert_eq!(data[0][0], 1); - assert_eq!(data[0][1], 2); - assert_eq!(data[0][2], 3); - assert_eq!(data[1][0], 6); - assert_eq!(data[1][1], 5); - assert_eq!(data[1][2], 4); - } - - #[test] - fn test_mut_indexing_item() { - let mut data = Arr2D::from(&[[1, 2, 3], [6, 5, 4]]); - data[(1, 2)] = 10; - data[(0, 0)] = 11; - assert_eq!(data[(0, 0)], 11); - assert_eq!(data[(0, 1)], 2); - assert_eq!(data[(0, 2)], 3); - assert_eq!(data[(1, 0)], 6); - assert_eq!(data[(1, 1)], 5); - assert_eq!(data[(1, 2)], 10); - } - - #[test] - fn test_2D_mut_indexing_item() { - let mut data = Arr2D::from(&[[1, 2, 3], [6, 5, 4]]); - data[(1, 2)] = 10; - data[(0, 0)] = 11; - assert_eq!(data[0][0], 11); - assert_eq!(data[0][1], 2); - assert_eq!(data[0][2], 3); - assert_eq!(data[1][0], 6); - assert_eq!(data[1][1], 5); - assert_eq!(data[1][2], 10); - } - - #[test] - fn test_indexing_row() { - let data = Arr2D::from(&[[1, 2, 3], [6, 5, 4]]); - assert_eq!(data[0], [1, 2, 3]); - assert_eq!(data[1], [6, 5, 4]); - } - #[test] - fn test_mut_indexing_row() { - let mut data = Arr2D::from(&[[1, 2, 3], [6, 5, 4]]); - data[0][2] = 9; - data[1][0] = 10; - assert_eq!(data[0], [1, 2, 9]); - assert_eq!(data[1], [10, 5, 4]); - } - - #[test] - #[should_panic] - fn test_index_out_of_bounds_panics() { - let data = Arr2D::from(&[[1, 2, 3], [4, 5, 6]]); - - let _ = data[(2, 0)]; - } - #[test] - #[should_panic] - fn test_mut_index_out_of_bounds_panics() { - let mut data = Arr2D::from(&[[1, 2, 3], [4, 5, 6]]); - - let _ = &mut data[(2, 0)]; - } - - #[test] - #[should_panic] - fn test_row_index_out_of_bounds_panics() { - let data = Arr2D::from(&[[1, 2, 3], [4, 5, 6]]); - - let _ = &data[2]; - } - #[test] - #[should_panic] - fn test_row_mut_index_out_of_bounds_panics() { - let mut data = Arr2D::from(&[[1, 2, 3], [4, 5, 6]]); - - let _ = &mut data[2]; - } - - // --- iterating --- - #[test] - fn test_rows_iterator_returns_slices() { - let data = Arr2D::from(&[[1, 2, 3], [4, 5, 6]]); - let rows: Vec<&[i32]> = data.rows().collect(); - - assert_eq!(rows.len(), 2); - assert_eq!(rows[0], &[1, 2, 3]); - assert_eq!(rows[1], &[4, 5, 6]); - } - - #[test] - fn test_iterating_through_items() { - let data = Arr2D::from(&[[1, 2, 3], [6, 5, 4]]); - let mut expected = [1, 2, 3, 6, 5, 4].iter(); - - for row in data.rows() { - for item in row { - assert_eq!(item, expected.next().unwrap()); - } - } - } - - #[test] - fn test_rows_mut_iterator_allows_mutation() { - let mut data = Arr2D::from(&[[1, 2, 3], [4, 5, 6]]); - - for row in data.rows_mut() { - row.reverse(); - } - - let expected = Arr2D::from(&[[3, 2, 1], [6, 5, 4]]); - assert_eq!(data, expected); - } - - // --- transformation --- - - #[test] - fn test_reshape() { - let mut data = Arr2D::from(&[[1, 2, 3], [6, 5, 4]]); - data.reshape(3).unwrap(); - let expected = Arr2D::from(&[[1, 2], [3, 6], [5, 4]]); - - assert_eq!(data, expected); - } - - #[test] - fn test_reshape_invalid_height_errors() { - let mut data = Arr2D::from(&[[1, 2], [3, 4]]); - let err = data - .reshape(3) - .expect_err("reshape should fail when new height mismatches size"); - - assert!(matches!( - err, - Arr2DError::InvalidReshape { - size: 4, - new_height: 3 - } - )); - } - - #[test] - fn test_transpose() { - let mut data = Arr2D::from(&[[1, 2, 3], [6, 5, 4]]); - data.transpose_mut(); - let expected = Arr2D::from(&[[1, 6], [2, 5], [3, 4]]); - assert_eq!(data, expected); - } - - #[test] - fn test_arr_from_vec() { - let data = Arr2D::try_from(vec![vec![1, 2, 3], vec![6, 5, 4]]).unwrap(); - let expected = Arr2D::from(&[[1, 2, 3], [6, 5, 4]]); - assert_eq!(data, expected); - } - - #[test] - fn test_try_from_inconsistent_rows_returns_error() { - let err = Arr2D::try_from(vec![vec![1, 2, 3], vec![4, 5]]) - .expect_err("rows with different widths should error"); - - assert!(matches!(err, Arr2DError::InconsistentRowLengths)); - } - - #[test] - fn test_try_from_empty_vec_creates_empty_arr() { - let data = Arr2D::try_from(Vec::>::new()).unwrap(); - - assert_eq!(data.rows().count(), 0); - } - - #[test] - fn test_map_transforms_elements() { - let data = Arr2D::from(&[[1, 2], [3, 4]]); - let mapped = data.map(|value| value * 2); - let expected = Arr2D::from(&[[2, 4], [6, 8]]); - - assert_eq!(mapped, expected); - } - - #[test] - fn test_full() { - let data = Arr2D::full(10, 3, 4); - for row in data.rows() { - for item in row { - assert_eq!(*item, 10); - } - } - } - - // --- from flat vector to Arr2D --- - - #[test] - fn test_from_flat() { - let data = Arr2D::from_flat(vec![1, 2, 3, 4, 5, 6], 0, 2, 3).unwrap(); - let out = Arr2D::from(&[[1, 2, 3], [4, 5, 6]]); - - assert_eq!(data, out); - } - - #[test] - fn test_from_flat_ref() { - let data = Arr2D::from_flat(&vec![1, 2, 3, 4, 5, 6], 0, 2, 3).unwrap(); - let out = Arr2D::from(&[[1, 2, 3], [4, 5, 6]]); - - assert_eq!(data, out); - } - - #[test] - fn test_from_flat_slice() { - let data = Arr2D::from_flat(&[1, 2, 3, 4, 5, 6], 0, 2, 3).unwrap(); - let out = Arr2D::from(&[[1, 2, 3], [4, 5, 6]]); - - assert_eq!(data, out); - } - - #[test] - fn test_from_flat_with_default() { - let data = Arr2D::from_flat(vec![1, 2, 3, 4], 0, 2, 3).unwrap(); - let out = Arr2D::from(&[[1, 2, 3], [4, 0, 0]]); - - assert_eq!(data, out); - } - - #[test] - fn test_from_flat_with_default_ref() { - let data = Arr2D::from_flat(&vec![1, 2, 3, 4], 0, 2, 3).unwrap(); - let out = Arr2D::from(&[[1, 2, 3], [4, 0, 0]]); - - assert_eq!(data, out); - } - - #[test] - fn test_from_flat_full_zeros() { - let flat_data = Arr2D::from_flat(&vec![], 0, 2, 3).unwrap(); - let full_data = Arr2D::full(0, 2, 3); - - assert_eq!(flat_data, full_data); - } - - #[test] - fn test_from_flat_slice_full_zeros() { - let data = Arr2D::from_flat(&[], 0, 2, 3).unwrap(); - let out = Arr2D::from(&[[0, 0, 0], [0, 0, 0]]); - - assert_eq!(data, out); - } - - #[test] - fn test_from_flat_slice_full_zeros_no_ref() { - let data = Arr2D::from_flat([], 0, 2, 3).unwrap(); - let out = Arr2D::from(&[[0, 0, 0], [0, 0, 0]]); - - assert_eq!(data, out); - } - - // --- dot product --- - - #[test] - fn test_mat_mul_mat() { - let arr1 = Arr2D::from_flat(vec![1, 2, 3, 4, 5, 6], 0, 2, 3).unwrap(); - let arr2 = Arr2D::from_flat(vec![7, 8, 9, 10, 11, 12], 0, 3, 2).unwrap(); - - let expected = Arr2D::from_flat(vec![58, 64, 139, 154], 0, 2, 2).unwrap(); - - let res = arr1.dot(&arr2).unwrap(); - assert_eq!(res, expected); - } - - #[test] - fn test_vec_mul_mat() { - let arr1 = Arr2D::from_flat(vec![3, 4, 2], 0, 1, 3).unwrap(); - let arr2 = Arr2D::from_flat(vec![13, 9, 7, 15, 8, 7, 4, 6, 6, 4, 0, 3], 0, 3, 4).unwrap(); - - let expected = Arr2D::from_flat(vec![83, 63, 37, 75], 0, 1, 4).unwrap(); - - let res = arr1.dot(&arr2).unwrap(); - assert_eq!(res, expected); - } - - #[test] - fn test_scalar_mul_mat() { - let scal = Arr2D::from_flat(vec![2], 0, 1, 1).unwrap(); - let mat = Arr2D::from_flat(vec![4, 0, 1, -9], 0, 2, 2).unwrap(); - - let expected = Arr2D::from_flat(vec![8, 0, 2, -18], 0, 2, 2).unwrap(); - - let res = scal.dot(&mat).unwrap(); - assert_eq!(res, expected); - } - - #[test] - fn test_Mul_trait_scalar() { - let scal = 2; - let mat = Arr2D::from_flat(vec![4, 0, 1, -9], 0, 2, 2).unwrap(); - - let expected = Arr2D::from_flat(vec![8, 0, 2, -18], 0, 2, 2).unwrap(); - - let res = mat * scal; - assert_eq!(res, expected) - } - - #[test] - fn test_Mul_trait_scalar_2() { - let scal = 2; - let mat = Arr2D::from_flat(vec![4, 0, 1, -9], 0, 2, 2).unwrap(); - - let expected = Arr2D::from_flat(vec![8, 0, 2, -18], 0, 2, 2).unwrap(); - - let res = &mat * scal; - assert_eq!(res, expected) - } - - #[test] - fn test_Mul_trait_mat_mul_mat() { - let arr1 = Arr2D::from_flat(vec![1, 2, 3, 4, 5, 6], 0, 2, 3).unwrap(); - let arr2 = Arr2D::from_flat(vec![7, 8, 9, 10, 11, 12], 0, 3, 2).unwrap(); - - let expected = Arr2D::from_flat(vec![58, 64, 139, 154], 0, 2, 2).unwrap(); - - let res = arr1 * arr2; - assert_eq!(res, expected); - } - - #[test] - fn test__borrowed_Mul_trait_mat_mul_mat() { - let arr1 = Arr2D::from_flat(vec![1, 2, 3, 4, 5, 6], 0, 2, 3).unwrap(); - let arr2 = Arr2D::from_flat(vec![7, 8, 9, 10, 11, 12], 0, 3, 2).unwrap(); - - let expected = Arr2D::from_flat(vec![58, 64, 139, 154], 0, 2, 2).unwrap(); - - let res = arr1 * &arr2; - assert_eq!(res, expected); - } - - #[test] - fn test__borrowed_Mul_trait_mat_mul_mat_2() { - let arr1 = Arr2D::from_flat(vec![1, 2, 3, 4, 5, 6], 0, 2, 3).unwrap(); - let arr2 = Arr2D::from_flat(vec![7, 8, 9, 10, 11, 12], 0, 3, 2).unwrap(); - - let expected = Arr2D::from_flat(vec![58, 64, 139, 154], 0, 2, 2).unwrap(); - - let res = &arr1 * arr2; - assert_eq!(res, expected); - } - - #[test] - fn test__double_borrowed_Mul_trait_mat_mul_mat() { - let arr1 = Arr2D::from_flat(vec![1, 2, 3, 4, 5, 6], 0, 2, 3).unwrap(); - let arr2 = Arr2D::from_flat(vec![7, 8, 9, 10, 11, 12], 0, 3, 2).unwrap(); - - let expected = Arr2D::from_flat(vec![58, 64, 139, 154], 0, 2, 2).unwrap(); - - let res = &arr1 * &arr2; - assert_eq!(res, expected); - } - - #[test] - fn test_Mul_trait_vec_mul_mat() { - let arr1 = Arr2D::from_flat(vec![3, 4, 2], 0, 1, 3).unwrap(); - let arr2 = Arr2D::from_flat(vec![13, 9, 7, 15, 8, 7, 4, 6, 6, 4, 0, 3], 0, 3, 4).unwrap(); - - let expected = Arr2D::from_flat(vec![83, 63, 37, 75], 0, 1, 4).unwrap(); - - let res = arr1 * arr2; - assert_eq!(res, expected); - } - - #[test] - fn test_Div_trait_mat_div_scalar() { - let mat = Arr2D::from_flat(vec![1.778, 0.0, 1.778], 0.0, 3, 1).unwrap(); - let result = mat / 1.778; - let expected = Arr2D::from(&[[1.0], [0.0], [1.0]]); - - assert_eq!(result, expected); - } - - #[test] - fn test_borrow_Div_trait_mat_div_scalar() { - let mat = Arr2D::from_flat(vec![1.778, 0.0, 1.778], 0.0, 3, 1).unwrap(); - let result = &mat / 1.778; - let expected = Arr2D::from(&[[1.0], [0.0], [1.0]]); - - assert_eq!(result, expected); - } - - // --- misc --- - - #[test] - fn test_display() { - let data = Arr2D::from(&[[1.2, 34.5678], [789.02, 0.123]]); - let out = format!("{data}"); - let expected = r#" -[[ 1.2, 34.5678 ] - [ 789.02, 0.123 ]]"#; - assert_eq!(&out, &expected[1..]); - } - - #[test] - fn test_int_max() { - let data = Arr2D::from(&[[12], [10], [11], [20], [9]]); - let result = data.max().unwrap(); - let expected = 20; - - assert_eq!(result, expected); - } - - #[test] - fn test_float_max() { - let data = Arr2D::from(&[[23.3], [12.4], [23.4], [10.4]]); - let result = data.max().unwrap(); - let expected = 23.4; - - assert_eq!(result, expected); - } - - #[test] - fn test_int_min() { - let data = Arr2D::from(&[[12], [10], [11], [20], [9]]); - let result = data.min().unwrap(); - let expected = 9; - - assert_eq!(result, expected); - } - - #[test] - fn test_float_min() { - let data = Arr2D::from(&[[23.3], [12.4], [23.4], [10.4]]); - let result = data.min().unwrap(); - let expected = 10.4; - - assert_eq!(result, expected); - } - - #[test] - fn test_known_inverse() { - let matrix = Arr2D::from(&[[3, 0, 2], [2, 0, -2], [0, 1, 1]]); - let result = matrix.inverse().unwrap(); - let rounded_result = result.round_to_decimal(1); - let expected = Arr2D::from(&[[0.2, 0.2, 0.0], [-0.2, 0.3, 1.0], [0.2, -0.3, 0.0]]); - - assert_eq!(rounded_result, expected); - } - - #[test] - fn test_known_inverse_2() { - let matrix = Arr2D::from(&[ - [3.556, -1.778, 0.0], - [-1.778, 3.556, -1.778], - [0.0, -1.778, 3.556], - ]); - let result = matrix.inverse().unwrap(); - let rounded_result = result.round_to_decimal(3); - let expected = Arr2D::from(&[ - [0.422, 0.281, 0.141], - [0.281, 0.562, 0.281], - [0.141, 0.281, 0.422], - ]); - - assert_eq!(rounded_result, expected); - } -} diff --git a/spindalis/src/utils/mod.rs b/spindalis/src/utils/mod.rs index 2a71aab..f515d9f 100644 --- a/spindalis/src/utils/mod.rs +++ b/spindalis/src/utils/mod.rs @@ -1,39 +1,9 @@ -#[allow(non_snake_case)] -pub mod arr2D; -pub mod substitution; pub mod variation; -pub use arr2D::Arr2D; -pub use arr2D::Rounding; -pub use substitution::back_substitution; -pub use substitution::forward_substitution; pub use variation::arith_mean; pub use variation::geom_mean; pub use variation::std_dev; -#[derive(Debug)] -pub enum Arr2DError { - InconsistentRowLengths, - NonSquareMatrix, - SingularMatrix, - InvalidReshape { - size: usize, - new_height: usize, - }, - InvalidShape { - input_size: usize, - output_size: usize, - }, - InvalidDotShape { - lhs: usize, - rhs: usize, - }, - ConversionFailed { - from: &'static str, - to: &'static str, - }, -} - #[derive(Copy, Clone)] pub enum StdDevType { Poulation, diff --git a/spindalis/src/utils/substitution.rs b/spindalis/src/utils/substitution.rs deleted file mode 100644 index f9deef9..0000000 --- a/spindalis/src/utils/substitution.rs +++ /dev/null @@ -1,34 +0,0 @@ -use crate::utils::arr2D::Arr2D; - -#[allow(clippy::needless_range_loop)] -pub fn back_substitution( - coeff_matrix: &Arr2D, - size: usize, - rhs_vector: &[f64], - solution: &mut [f64], -) { - solution[size - 1] = rhs_vector[size - 1] / coeff_matrix[size - 1][size - 1]; - for i in (0..(size - 1)).rev() { - let mut sum = 0.0; - for j in (i + 1)..size { - sum += coeff_matrix[i][j] * solution[j] - } - solution[i] = (rhs_vector[i] - sum) / coeff_matrix[i][i]; - } -} - -#[allow(clippy::needless_range_loop)] -pub fn forward_substitution( - coeff_matrix: &Arr2D, - size: usize, - rhs_vector: &[f64], - solution: &mut [f64], -) { - for i in 0..size { - let mut sum = 0.0; - for j in 0..i { - sum += coeff_matrix[i][j] * solution[j] - } - solution[i] = (rhs_vector[i] - sum) / coeff_matrix[i][i]; - } -}