From f32670f65f61a54e31a5d3118515cf380835d9b6 Mon Sep 17 00:00:00 2001 From: Axect Date: Thu, 24 Oct 2024 09:10:34 +0900 Subject: [PATCH 1/5] CHGE: MatrixTrait --- src/fuga/mod.rs | 8 + src/numerical/eigen.rs | 1 + src/numerical/optimize.rs | 1 + src/numerical/utils.rs | 1 + src/prelude/mod.rs | 10 ++ src/special/lanczos.rs | 1 + src/statistics/stat.rs | 2 +- src/structure/matrix.rs | 333 ++++++++++++++++++-------------------- src/traits/matrix.rs | 25 +++ src/traits/mod.rs | 1 + src/traits/pointer.rs | 1 + src/traits/sugar.rs | 1 + src/util/non_macro.rs | 1 + 13 files changed, 211 insertions(+), 175 deletions(-) create mode 100644 src/traits/matrix.rs diff --git a/src/fuga/mod.rs b/src/fuga/mod.rs index 5b4942d2..098c9587 100644 --- a/src/fuga/mod.rs +++ b/src/fuga/mod.rs @@ -157,6 +157,7 @@ pub use crate::traits::{ fp::{FPMatrix, FPVector}, general::Algorithm, math::{InnerProduct, LinearOp, MatrixProduct, Norm, Normed, Vector, VectorProduct}, + matrix::MatrixTrait, mutable::{MutFP, MutMatrix}, num::Real, pointer::{MatrixPtr, Oxide, Redox, RedoxCommon}, @@ -171,7 +172,14 @@ pub use crate::traits::{ mutable::ParallelMutFP, }; +#[cfg(feature = "complex")] +pub use crate::complex::{ + C64, + matrix::*, +}; + #[allow(unused_imports)] +#[allow(ambiguous_glob_reexports)] pub use crate::structure::{ matrix::*, polynomial::*, diff --git a/src/numerical/eigen.rs b/src/numerical/eigen.rs index 688764eb..4632ecc9 100644 --- a/src/numerical/eigen.rs +++ b/src/numerical/eigen.rs @@ -4,6 +4,7 @@ pub use self::EigenMethod::*; use crate::structure::matrix::Matrix; +use crate::traits::matrix::MatrixTrait; use crate::util::non_macro::eye_shape; #[derive(Debug, Copy, Clone)] diff --git a/src/numerical/optimize.rs b/src/numerical/optimize.rs index b894e28b..765fc5d2 100644 --- a/src/numerical/optimize.rs +++ b/src/numerical/optimize.rs @@ -112,6 +112,7 @@ pub use self::OptMethod::{GaussNewton, GradientDescent, LevenbergMarquardt}; use self::OptOption::{InitParam, MaxIter}; use crate::numerical::utils::jacobian; use crate::structure::matrix::{LinearAlgebra, Matrix}; +use crate::traits::matrix::MatrixTrait; use crate::structure::ad::{AD, ADVec}; use crate::util::useful::max; use std::collections::HashMap; diff --git a/src/numerical/utils.rs b/src/numerical/utils.rs index 158aafea..f28ca5a9 100644 --- a/src/numerical/utils.rs +++ b/src/numerical/utils.rs @@ -1,4 +1,5 @@ use crate::structure::matrix::*; +use crate::traits::matrix::MatrixTrait; use crate::structure::ad::*; use crate::structure::ad::AD::*; use crate::util::non_macro::{cat, zeros}; diff --git a/src/prelude/mod.rs b/src/prelude/mod.rs index d10ff75e..d111b760 100644 --- a/src/prelude/mod.rs +++ b/src/prelude/mod.rs @@ -153,6 +153,7 @@ pub use crate::traits::{ fp::{FPMatrix, FPVector}, general::Algorithm, math::{InnerProduct, LinearOp, MatrixProduct, Vector, VectorProduct}, + matrix::MatrixTrait, mutable::{MutFP, MutMatrix}, num::Real, pointer::{MatrixPtr, Oxide, Redox, RedoxCommon}, @@ -184,6 +185,15 @@ pub use crate::structure::dataframe::WithCSV; #[cfg(feature="nc")] pub use crate::structure::dataframe::WithNetCDF; +#[cfg(feature = "complex")] +#[allow(ambiguous_glob_reexports)] +pub use crate::complex::{ + C64, + matrix::*, + vector::*, + integral::*, +}; + pub use simpler::{solve, SimplerLinearAlgebra}; #[allow(unused_imports)] diff --git a/src/special/lanczos.rs b/src/special/lanczos.rs index eeded0dd..bf563d9d 100644 --- a/src/special/lanczos.rs +++ b/src/special/lanczos.rs @@ -2,6 +2,7 @@ use crate::statistics::ops::{double_factorial, factorial, C}; use crate::structure::matrix::Matrix; +use crate::traits::matrix::MatrixTrait; use crate::traits::pointer::{Oxide, RedoxCommon}; use crate::util::non_macro::zeros; use crate::util::useful::sgn; diff --git a/src/statistics/stat.rs b/src/statistics/stat.rs index 0c44b3bb..41fb4dec 100644 --- a/src/statistics/stat.rs +++ b/src/statistics/stat.rs @@ -151,7 +151,7 @@ use std::fmt; use self::QType::*; //use crate::structure::dataframe::*; use crate::structure::matrix::*; -use crate::traits::fp::FPVector; +use crate::traits::matrix::MatrixTrait; use order_stat::kth_by; /// Statistics Trait diff --git a/src/structure/matrix.rs b/src/structure/matrix.rs index 66d82553..2459d12e 100644 --- a/src/structure/matrix.rs +++ b/src/structure/matrix.rs @@ -624,6 +624,7 @@ use crate::traits::{ general::Algorithm, math::{InnerProduct, LinearOp, MatrixProduct, Norm, Normed, Vector}, mutable::MutMatrix, + matrix::MatrixTrait, }; #[cfg(feature = "parallel")] use crate::traits::math::{ParallelInnerProduct, ParallelNormed}; @@ -821,14 +822,16 @@ impl PartialEq for Matrix { /// Main matrix structure #[allow(dead_code)] -impl Matrix { +impl MatrixTrait for Matrix { + type Scalar = f64; + /// Raw pointer for `self.data` - pub fn ptr(&self) -> *const f64 { + fn ptr(&self) -> *const f64 { &self.data[0] as *const f64 } /// Raw mutable pointer for `self.data` - pub fn mut_ptr(&mut self) -> *mut f64 { + fn mut_ptr(&mut self) -> *mut f64 { &mut self.data[0] as *mut f64 } @@ -842,7 +845,7 @@ impl Matrix { /// let b = a.as_slice(); /// assert_eq!(b, &[1f64,2f64,3f64,4f64]); /// ``` - pub fn as_slice(&self) -> &[f64] { + fn as_slice(&self) -> &[f64] { &self.data[..] } @@ -858,7 +861,7 @@ impl Matrix { /// assert_eq!(b, &[5f64, 2f64, 3f64, 4f64]); /// assert_eq!(a, matrix(vec![5,2,3,4], 2, 2, Col)); /// ``` - pub fn as_mut_slice(&mut self) -> &mut [f64] { + fn as_mut_slice(&mut self) -> &mut [f64] { &mut self.data[..] } @@ -875,7 +878,7 @@ impl Matrix { /// let b = a.change_shape(); /// assert_eq!(b.shape, Col); /// ``` - pub fn change_shape(&self) -> Self { + fn change_shape(&self) -> Self { let r = self.row; let c = self.col; assert_eq!(r * c, self.data.len()); @@ -916,7 +919,7 @@ impl Matrix { /// let b = a.change_shape(); /// assert_eq!(b.shape, Col); /// ``` - pub fn change_shape_mut(&mut self) { + fn change_shape_mut(&mut self) { let r = self.row; let c = self.col; assert_eq!(r * c, self.data.len()); @@ -956,7 +959,7 @@ impl Matrix { /// // r[0] 1 3 /// // r[1] 2 4 /// ``` - pub fn spread(&self) -> String { + fn spread(&self) -> String { assert_eq!(self.row * self.col, self.data.len()); let r = self.row; let c = self.col; @@ -977,10 +980,10 @@ impl Matrix { }; return format!( "Result is too large to print - {}x{}\nonly print {}x{} parts:\n{}", - self.row.to_string(), - self.col.to_string(), - key_row.to_string(), - key_col.to_string(), + self.row, + self.col, + key_row, + key_col, part.spread() ); } @@ -992,7 +995,7 @@ impl Matrix { .map( |x| min(format!("{:.4}", x).len(), x.to_string().len()), // Choose minimum of approx vs normal ) - .fold(0, |x, y| max(x, y)) + .fold(0, max) + 1; if space < 5 { @@ -1027,7 +1030,7 @@ impl Matrix { result.push('\n'); } - return result; + result } /// Extract Column @@ -1043,7 +1046,7 @@ impl Matrix { /// assert_eq!(a.col(0), c!(1,3)); /// } /// ``` - pub fn col(&self, index: usize) -> Vec { + fn col(&self, index: usize) -> Vec { assert!(index < self.col); let mut container: Vec = vec![0f64; self.row]; for i in 0..self.row { @@ -1065,7 +1068,7 @@ impl Matrix { /// assert_eq!(a.row(0), c!(1,2)); /// } /// ``` - pub fn row(&self, index: usize) -> Vec { + fn row(&self, index: usize) -> Vec { assert!(index < self.row); let mut container: Vec = vec![0f64; self.col]; for i in 0..self.col { @@ -1087,7 +1090,7 @@ impl Matrix { /// assert_eq!(a.diag(), c!(1,4)); /// } /// ``` - pub fn diag(&self) -> Vec { + fn diag(&self) -> Vec { let mut container = vec![0f64; self.row]; let r = self.row; let c = self.col; @@ -1108,7 +1111,7 @@ impl Matrix { /// let a = matrix(vec![1,2,3,4], 2, 2, Row); /// println!("{}", a); // [[1,3],[2,4]] /// ``` - pub fn transpose(&self) -> Self { + fn transpose(&self) -> Self { match self.shape { Row => matrix(self.data.clone(), self.col, self.row, Col), Col => matrix(self.data.clone(), self.col, self.row, Row), @@ -1128,10 +1131,146 @@ impl Matrix { /// assert_eq!(a.transpose(), a.t()); /// } /// ``` - pub fn t(&self) -> Self { + fn t(&self) -> Self { self.transpose() } + /// Substitute Col + #[inline] + fn subs_col(&mut self, idx: usize, v: &[f64]) { + for i in 0..self.row { + self[(i, idx)] = v[i]; + } + } + + /// Substitute Row + #[inline] + fn subs_row(&mut self, idx: usize, v: &[f64]) { + for j in 0..self.col { + self[(idx, j)] = v[j]; + } + } + + /// From index operations + fn from_index(f: F, size: (usize, usize)) -> Matrix + where + F: Fn(usize, usize) -> G + Copy, + G: Into, + { + let row = size.0; + let col = size.1; + + let mut mat = matrix(vec![0f64; row * col], row, col, Row); + + for i in 0..row { + for j in 0..col { + mat[(i, j)] = f(i, j).into(); + } + } + mat + } + + /// Matrix to `Vec>` + /// + /// To send `Matrix` to `inline-python` + fn to_vec(&self) -> Vec> { + let mut result = vec![vec![0f64; self.col]; self.row]; + for i in 0..self.row { + result[i] = self.row(i); + } + result + } + + fn to_diag(&self) -> Matrix { + assert_eq!(self.row, self.col, "Should be square matrix"); + let mut result = matrix(vec![0f64; self.row * self.col], self.row, self.col, Row); + let diag = self.diag(); + for i in 0..self.row { + result[(i, i)] = diag[i]; + } + result + } + + /// Submatrix + /// + /// # Description + /// Return below elements of matrix to new matrix + /// + /// $$ + /// \begin{pmatrix} + /// \\ddots & & & & \\\\ + /// & start & \\cdots & end.1 & \\\\ + /// & \\vdots & \\ddots & \\vdots & \\\\ + /// & end.0 & \\cdots & end & \\\\ + /// & & & & \\ddots + /// \end{pmatrix} + /// $$ + /// + /// # Examples + /// ``` + /// extern crate peroxide; + /// use peroxide::fuga::*; + /// + /// fn main() { + /// let a = ml_matrix("1 2 3;4 5 6;7 8 9"); + /// let b = ml_matrix("5 6;8 9"); + /// let c = a.submat((1, 1), (2, 2)); + /// assert_eq!(b, c); + /// } + /// ``` + fn submat(&self, start: (usize, usize), end: (usize, usize)) -> Matrix { + let row = end.0 + 1 - start.0; + let col = end.1 + 1 - start.1; + let mut result = matrix(vec![0f64; row * col], row, col, self.shape); + for i in 0..row { + for j in 0..col { + result[(i, j)] = self[(start.0 + i, start.1 + j)]; + } + } + result + } + + /// Substitute matrix to specific position + /// + /// # Description + /// Substitute below elements of matrix + /// + /// $$ + /// \begin{pmatrix} + /// \\ddots & & & & \\\\ + /// & start & \\cdots & end.1 & \\\\ + /// & \\vdots & \\ddots & \\vdots & \\\\ + /// & end.0 & \\cdots & end & \\\\ + /// & & & & \\ddots + /// \end{pmatrix} + /// $$ + /// + /// # Examples + /// ``` + /// extern crate peroxide; + /// use peroxide::fuga::*; + /// + /// fn main() { + /// let mut a = ml_matrix("1 2 3;4 5 6;7 8 9"); + /// let b = ml_matrix("1 2;3 4"); + /// let c = ml_matrix("1 2 3;4 1 2;7 3 4"); + /// a.subs_mat((1,1), (2,2), &b); + /// assert_eq!(a, c); + /// } + /// ``` + fn subs_mat(&mut self, start: (usize, usize), end: (usize, usize), m: &Matrix) { + let row = end.0 - start.0 + 1; + let col = end.1 - start.1 + 1; + for i in 0..row { + for j in 0..col { + self[(start.0 + i, start.1 + j)] = m[(i, j)]; + } + } + } + +} + +impl Matrix { /// Write to CSV /// /// # Examples @@ -1289,159 +1428,6 @@ impl Matrix { Ok(m) } - - /// Should check shape - pub fn subs(&mut self, idx: usize, v: &Vec) { - let p = &mut self.mut_ptr(); - match self.shape { - Row => { - let c = self.col; - unsafe { - p.add(idx * c).copy_from(v.as_ptr(), c); - } - } - Col => { - let r = self.row; - unsafe { - p.add(idx * r).copy_from(v.as_ptr(), r); - } - } - } - } - - /// Substitute Col - #[inline] - pub fn subs_col(&mut self, idx: usize, v: &Vec) { - for i in 0..self.row { - self[(i, idx)] = v[i]; - } - } - - /// Substitute Row - #[inline] - pub fn subs_row(&mut self, idx: usize, v: &Vec) { - for j in 0..self.col { - self[(idx, j)] = v[j]; - } - } - - /// From index operations - pub fn from_index(f: F, size: (usize, usize)) -> Matrix - where - F: Fn(usize, usize) -> G + Copy, - G: Into, - { - let row = size.0; - let col = size.1; - - let mut mat = matrix(vec![0f64; row * col], row, col, Row); - - for i in 0..row { - for j in 0..col { - mat[(i, j)] = f(i, j).into(); - } - } - mat - } - - /// Matrix to `Vec>` - /// - /// To send `Matrix` to `inline-python` - pub fn to_vec(&self) -> Vec> { - let mut result = vec![vec![0f64; self.col]; self.row]; - for i in 0..self.row { - result[i] = self.row(i); - } - result - } - - pub fn to_diag(&self) -> Matrix { - assert_eq!(self.row, self.col, "Should be square matrix"); - let mut result = matrix(vec![0f64; self.row * self.col], self.row, self.col, Row); - let diag = self.diag(); - for i in 0..self.row { - result[(i, i)] = diag[i]; - } - result - } - - /// Submatrix - /// - /// # Description - /// Return below elements of matrix to new matrix - /// - /// $$ - /// \begin{pmatrix} - /// \\ddots & & & & \\\\ - /// & start & \\cdots & end.1 & \\\\ - /// & \\vdots & \\ddots & \\vdots & \\\\ - /// & end.0 & \\cdots & end & \\\\ - /// & & & & \\ddots - /// \end{pmatrix} - /// $$ - /// - /// # Examples - /// ``` - /// extern crate peroxide; - /// use peroxide::fuga::*; - /// - /// fn main() { - /// let a = ml_matrix("1 2 3;4 5 6;7 8 9"); - /// let b = ml_matrix("5 6;8 9"); - /// let c = a.submat((1, 1), (2, 2)); - /// assert_eq!(b, c); - /// } - /// ``` - pub fn submat(&self, start: (usize, usize), end: (usize, usize)) -> Matrix { - let row = end.0 - start.0 + 1; - let col = end.1 - start.1 + 1; - let mut result = matrix(vec![0f64; row * col], row, col, self.shape); - for i in 0..row { - for j in 0..col { - result[(i, j)] = self[(start.0 + i, start.1 + j)]; - } - } - result - } - - /// Substitute matrix to specific position - /// - /// # Description - /// Substitute below elements of matrix - /// - /// $$ - /// \begin{pmatrix} - /// \\ddots & & & & \\\\ - /// & start & \\cdots & end.1 & \\\\ - /// & \\vdots & \\ddots & \\vdots & \\\\ - /// & end.0 & \\cdots & end & \\\\ - /// & & & & \\ddots - /// \end{pmatrix} - /// $$ - /// - /// # Examples - /// ``` - /// extern crate peroxide; - /// use peroxide::fuga::*; - /// - /// fn main() { - /// let mut a = ml_matrix("1 2 3;4 5 6;7 8 9"); - /// let b = ml_matrix("1 2;3 4"); - /// let c = ml_matrix("1 2 3;4 1 2;7 3 4"); - /// a.subs_mat((1,1), (2,2), &b); - /// assert_eq!(a, c); - /// } - /// ``` - pub fn subs_mat(&mut self, start: (usize, usize), end: (usize, usize), m: &Matrix) { - let row = end.0 - start.0 + 1; - let col = end.1 - start.1 + 1; - for i in 0..row { - for j in 0..col { - self[(start.0 + i, start.1 + j)] = m[(i, j)]; - } - } - } - /// Matrix from series /// /// # Example @@ -1483,6 +1469,7 @@ impl Matrix { .collect::>(); matrix(data, row, col, Row) } + } // ============================================================================= @@ -1751,8 +1738,6 @@ impl ParallelInnerProduct for Matrix { } } -/// TODO: Transpose - /// Matrix as Linear operator for Vector #[allow(non_snake_case)] impl LinearOp, Vec> for Matrix { diff --git a/src/traits/matrix.rs b/src/traits/matrix.rs new file mode 100644 index 00000000..0f08a6b8 --- /dev/null +++ b/src/traits/matrix.rs @@ -0,0 +1,25 @@ +pub trait MatrixTrait: Sized { + type Scalar; + fn ptr(&self) -> *const Self::Scalar; + fn mut_ptr(&mut self) -> *mut Self::Scalar; + fn as_slice(&self) -> &[Self::Scalar]; + fn as_mut_slice(&mut self) -> &mut [Self::Scalar]; + fn change_shape(&self) -> Self; + fn change_shape_mut(&mut self); + fn spread(&self) -> String; + fn col(&self, index: usize) -> Vec; + fn row(&self, index: usize) -> Vec; + fn diag(&self) -> Vec; + fn transpose(&self) -> Self; + fn t(&self) -> Self { self.transpose() } + fn subs_col(&mut self, idx: usize, v: &[Self::Scalar]); + fn subs_row(&mut self, idx: usize, v: &[Self::Scalar]); + fn from_index(f: F, size: (usize, usize)) -> Self + where + F: Fn(usize, usize) -> G + Copy, + G: Into; + fn to_vec(&self) -> Vec>; + fn to_diag(&self) -> Self; + fn submat(&self, start: (usize, usize), end: (usize, usize)) -> Self; + fn subs_mat(&mut self, start: (usize, usize), end: (usize, usize), m: &Self); +} diff --git a/src/traits/mod.rs b/src/traits/mod.rs index 2709f2db..67a0d15b 100644 --- a/src/traits/mod.rs +++ b/src/traits/mod.rs @@ -2,6 +2,7 @@ pub mod float; pub mod fp; pub mod general; pub mod math; +pub mod matrix; pub mod mutable; pub mod num; pub mod pointer; diff --git a/src/traits/pointer.rs b/src/traits/pointer.rs index 65fab349..7551cbb7 100644 --- a/src/traits/pointer.rs +++ b/src/traits/pointer.rs @@ -49,6 +49,7 @@ use crate::structure::ad::AD; use crate::traits::{ fp::FPVector, math::{LinearOp, Vector}, + matrix::MatrixTrait, }; use std::ops::{Add, Deref, Div, Mul, Sub}; diff --git a/src/traits/sugar.rs b/src/traits/sugar.rs index e1ce7f1e..d9286c43 100644 --- a/src/traits/sugar.rs +++ b/src/traits/sugar.rs @@ -1,4 +1,5 @@ use crate::structure::matrix::{matrix, Matrix, Shape}; +use crate::traits::matrix::MatrixTrait; use crate::traits::fp::FPVector; use crate::util::non_macro::zeros_shape; use std::ops::{Add, Div, Mul, Sub}; diff --git a/src/util/non_macro.rs b/src/util/non_macro.rs index 8c6f69b2..d50bd3b9 100644 --- a/src/util/non_macro.rs +++ b/src/util/non_macro.rs @@ -38,6 +38,7 @@ use crate::structure::{ matrix::{matrix, Matrix, Shape}, }; use crate::traits::float::FloatWithPrecision; +use crate::traits::matrix::MatrixTrait; use anyhow::{Result, bail}; #[derive(Debug, Copy, Clone)] From d2a93b43e17a262decc804205e65ea07c1e81ebf Mon Sep 17 00:00:00 2001 From: Axect Date: Thu, 24 Oct 2024 09:53:33 +0900 Subject: [PATCH 2/5] CHGE: Generic LinearAlgebra --- src/structure/matrix.rs | 121 +++++----------------------------------- src/traits/matrix.rs | 97 +++++++++++++++++++++++++++++++- 2 files changed, 109 insertions(+), 109 deletions(-) diff --git a/src/structure/matrix.rs b/src/structure/matrix.rs index 2459d12e..c4dbab42 100644 --- a/src/structure/matrix.rs +++ b/src/structure/matrix.rs @@ -624,7 +624,7 @@ use crate::traits::{ general::Algorithm, math::{InnerProduct, LinearOp, MatrixProduct, Norm, Normed, Vector}, mutable::MutMatrix, - matrix::MatrixTrait, + matrix::{MatrixTrait, LinearAlgebra, PQLU, WAZD, QR, SVD, Form, SolveKind}, }; #[cfg(feature = "parallel")] use crate::traits::math::{ParallelInnerProduct, ParallelNormed}; @@ -2882,27 +2882,6 @@ impl FPMatrix for Matrix { // ============================================================================= // Linear Algebra // ============================================================================= - -/// Linear algebra trait -pub trait LinearAlgebra { - fn back_subs(&self, b: &Vec) -> Vec; - fn forward_subs(&self, b: &Vec) -> Vec; - fn lu(&self) -> PQLU; - fn waz(&self, d_form: Form) -> Option; - fn qr(&self) -> QR; - fn svd(&self) -> SVD; - #[cfg(feature = "O3")] - fn cholesky(&self, uplo: UPLO) -> Matrix; - fn rref(&self) -> Matrix; - fn det(&self) -> f64; - fn block(&self) -> (Matrix, Matrix, Matrix, Matrix); - fn inv(&self) -> Matrix; - fn pseudo_inv(&self) -> Matrix; - fn solve(&self, b: &Vec, sk: SolveKind) -> Vec; - fn solve_mat(&self, m: &Matrix, sk: SolveKind) -> Matrix; - fn is_symmetric(&self) -> bool; -} - pub fn diag(n: usize) -> Matrix { let mut v: Vec = vec![0f64; n * n]; for i in 0..n { @@ -2912,30 +2891,7 @@ pub fn diag(n: usize) -> Matrix { matrix(v, n, n, Row) } -/// Data structure for Complete Pivoting LU decomposition -/// -/// # Usage -/// ```rust -/// extern crate peroxide; -/// use peroxide::fuga::*; -/// -/// let a = ml_matrix("1 2;3 4"); -/// let pqlu = a.lu(); -/// let (p, q, l, u) = pqlu.extract(); -/// // p, q are permutations -/// // l, u are matrices -/// l.print(); // lower triangular -/// u.print(); // upper triangular -/// ``` -#[derive(Debug, Clone)] -pub struct PQLU { - pub p: Vec, - pub q: Vec, - pub l: Matrix, - pub u: Matrix, -} - -impl PQLU { +impl PQLU { pub fn extract(&self) -> (Vec, Vec, Matrix, Matrix) { ( self.p.clone(), @@ -2982,43 +2938,7 @@ impl PQLU { } } -#[derive(Debug, Clone)] -pub struct WAZD { - pub w: Matrix, - pub z: Matrix, - pub d: Matrix, -} - -#[derive(Debug, Copy, Clone)] -pub enum Form { - Diagonal, - Identity, -} - -#[derive(Debug, Clone)] -pub struct QR { - pub q: Matrix, - pub r: Matrix, -} - -impl QR { - pub fn q(&self) -> &Matrix { - &self.q - } - - pub fn r(&self) -> &Matrix { - &self.r - } -} - -#[derive(Debug, Clone)] -pub struct SVD { - pub s: Vec, - pub u: Matrix, - pub vt: Matrix, -} - -impl SVD { +impl SVD { pub fn u(&self) -> &Matrix { &self.u } @@ -3076,15 +2996,9 @@ impl SVD { } } -#[derive(Debug, Copy, Clone)] -pub enum SolveKind { - LU, - WAZ, -} - -impl LinearAlgebra for Matrix { +impl LinearAlgebra for Matrix { /// Backward Substitution for Upper Triangular - fn back_subs(&self, b: &Vec) -> Vec { + fn back_subs(&self, b: &[f64]) -> Vec { let n = self.col; let mut y = vec![0f64; n]; y[n - 1] = b[n - 1] / self[(n - 1, n - 1)]; @@ -3099,7 +3013,7 @@ impl LinearAlgebra for Matrix { } /// Forward substitution for Lower Triangular - fn forward_subs(&self, b: &Vec) -> Vec { + fn forward_subs(&self, b: &[f64]) -> Vec { let n = self.col; let mut y = vec![0f64; n]; y[0] = b[0] / self[(0, 0)]; @@ -3127,8 +3041,6 @@ impl LinearAlgebra for Matrix { /// /// # Examples /// ``` - /// #[macro_use] - /// extern crate peroxide; /// use peroxide::fuga::*; /// /// fn main() { @@ -3137,11 +3049,11 @@ impl LinearAlgebra for Matrix { /// let (p,q,l,u) = (pqlu.p, pqlu.q, pqlu.l, pqlu.u); /// assert_eq!(p, vec![1]); // swap 0 & 1 (Row) /// assert_eq!(q, vec![1]); // swap 0 & 1 (Col) - /// assert_eq!(l, matrix(c!(1,0,0.5,1),2,2,Row)); - /// assert_eq!(u, matrix(c!(4,3,0,-0.5),2,2,Row)); + /// assert_eq!(l, matrix(vec![1,0,0.5,1],2,2,Row)); + /// assert_eq!(u, matrix(vec![4,3,0,-0.5],2,2,Row)); /// } /// ``` - fn lu(&self) -> PQLU { + fn lu(&self) -> PQLU { assert_eq!(self.col, self.row); let n = self.row; let len: usize = n * n; @@ -3173,7 +3085,7 @@ impl LinearAlgebra for Matrix { PQLU { p, q, l, u } } - fn waz(&self, d_form: Form) -> Option { + fn waz(&self, d_form: Form) -> Option> { match d_form { Form::Diagonal => { let n = self.row; @@ -3241,7 +3153,6 @@ impl LinearAlgebra for Matrix { /// /// # Example /// ``` - /// extern crate peroxide; /// use peroxide::fuga::*; /// /// fn main() { @@ -3256,7 +3167,7 @@ impl LinearAlgebra for Matrix { /// } /// ``` #[allow(non_snake_case)] - fn qr(&self) -> QR { + fn qr(&self) -> QR { match () { #[cfg(feature = "O3")] () => { @@ -3298,7 +3209,6 @@ impl LinearAlgebra for Matrix { /// /// # Examples /// ``` - /// extern crate peroxide; /// use peroxide::fuga::*; /// /// fn main() { @@ -3311,7 +3221,7 @@ impl LinearAlgebra for Matrix { /// a.print(); /// } /// ``` - fn svd(&self) -> SVD { + fn svd(&self) -> SVD { match () { #[cfg(feature = "O3")] () => { @@ -3634,7 +3544,7 @@ impl LinearAlgebra for Matrix { /// /// * Biswa Nath Datta, *Numerical Linear Algebra and Applications, Second Edition* /// * Ke Chen, *Matrix Preconditioning Techniques and Applications*, Cambridge Monographs on Applied and Computational Mathematics - fn solve(&self, b: &Vec, sk: SolveKind) -> Vec { + fn solve(&self, b: &[f64], sk: SolveKind) -> Vec { match sk { #[cfg(feature = "O3")] SolveKind::LU => { @@ -3735,11 +3645,6 @@ impl LinearAlgebra for Matrix { } } -#[allow(non_snake_case)] -pub fn solve(A: &Matrix, b: &Matrix, sk: SolveKind) -> Matrix { - A.solve_mat(b, sk) -} - impl MutMatrix for Matrix { type Scalar = f64; diff --git a/src/traits/matrix.rs b/src/traits/matrix.rs index 0f08a6b8..85b7b7fd 100644 --- a/src/traits/matrix.rs +++ b/src/traits/matrix.rs @@ -1,4 +1,6 @@ -pub trait MatrixTrait: Sized { +use std::ops::{Index, IndexMut}; + +pub trait MatrixTrait: Sized + Index<(usize, usize), Output=Self::Scalar> + IndexMut<(usize, usize)> { type Scalar; fn ptr(&self) -> *const Self::Scalar; fn mut_ptr(&mut self) -> *mut Self::Scalar; @@ -23,3 +25,96 @@ pub trait MatrixTrait: Sized { fn submat(&self, start: (usize, usize), end: (usize, usize)) -> Self; fn subs_mat(&mut self, start: (usize, usize), end: (usize, usize), m: &Self); } + +// ┌─────────────────────────────────────────────────────────┐ +// For Linear Algebra +// └─────────────────────────────────────────────────────────┘ +/// Linear algebra trait +pub trait LinearAlgebra { + fn back_subs(&self, b: &[M::Scalar]) -> Vec; + fn forward_subs(&self, b: &[M::Scalar]) -> Vec; + fn lu(&self) -> PQLU; + fn waz(&self, d_form: Form) -> Option>; + fn qr(&self) -> QR; + fn svd(&self) -> SVD; + #[cfg(feature = "O3")] + fn cholesky(&self, uplo: UPLO) -> M; + fn rref(&self) -> M; + fn det(&self) -> f64; + fn block(&self) -> (M, M, M, M); + fn inv(&self) -> M; + fn pseudo_inv(&self) -> M; + fn solve(&self, b: &[M::Scalar], sk: SolveKind) -> Vec; + fn solve_mat(&self, m: &M, sk: SolveKind) -> M; + fn is_symmetric(&self) -> bool; +} + +#[allow(non_snake_case)] +pub fn solve>(A: &M, b: &M, sk: SolveKind) -> M { + A.solve_mat(b, sk) +} + +/// Data structure for Complete Pivoting LU decomposition +/// +/// # Usage +/// ```rust +/// use peroxide::fuga::*; +/// +/// let a = ml_matrix("1 2;3 4"); +/// let pqlu = a.lu(); +/// let (p, q, l, u) = pqlu.extract(); +/// // p, q are permutations +/// // l, u are matrices +/// l.print(); // lower triangular +/// u.print(); // upper triangular +/// ``` +#[derive(Debug, Clone)] +pub struct PQLU { + pub p: Vec, + pub q: Vec, + pub l: M, + pub u: M, +} + +#[derive(Debug, Clone)] +pub struct WAZD { + pub w: M, + pub z: M, + pub d: M, +} + +#[derive(Debug, Clone)] +pub struct QR { + pub q: M, + pub r: M, +} + +#[derive(Debug, Copy, Clone)] +pub enum Form { + Diagonal, + Identity, +} + +#[derive(Debug, Copy, Clone)] +pub enum SolveKind { + LU, + WAZ, +} + +impl QR { + pub fn q(&self) -> &M { + &self.q + } + + pub fn r(&self) -> &M { + &self.r + } +} + +#[derive(Debug, Clone)] +pub struct SVD { + pub s: Vec, + pub u: M, + pub vt: M, +} + From cae917715f69b3c75bec63d665598a96677af262 Mon Sep 17 00:00:00 2001 From: Axect Date: Thu, 24 Oct 2024 10:30:26 +0900 Subject: [PATCH 3/5] CHGE: Fix all errors in Generic LinearAlgebra --- src/fuga/mod.rs | 6 +-- src/numerical/newton.rs | 2 +- src/numerical/optimize.rs | 4 +- src/numerical/spline.rs | 1 + src/prelude/mod.rs | 4 +- src/prelude/simpler.rs | 79 ++++++++++++++++++++------------------- src/statistics/stat.rs | 2 +- src/structure/matrix.rs | 8 ++-- src/structure/sparse.rs | 19 +++++----- 9 files changed, 64 insertions(+), 61 deletions(-) diff --git a/src/fuga/mod.rs b/src/fuga/mod.rs index 098c9587..60dcfcb3 100644 --- a/src/fuga/mod.rs +++ b/src/fuga/mod.rs @@ -157,7 +157,7 @@ pub use crate::traits::{ fp::{FPMatrix, FPVector}, general::Algorithm, math::{InnerProduct, LinearOp, MatrixProduct, Norm, Normed, Vector, VectorProduct}, - matrix::MatrixTrait, + matrix::{MatrixTrait, LinearAlgebra}, mutable::{MutFP, MutMatrix}, num::Real, pointer::{MatrixPtr, Oxide, Redox, RedoxCommon}, @@ -234,11 +234,11 @@ pub use crate::numerical::integral::Integral::{ pub use crate::statistics::stat::QType::{ Type1, Type2, Type3, Type4, Type5, Type6, Type7, Type8, Type9, }; -pub use crate::structure::matrix::{ +pub use crate::traits::matrix::{ Form::{Diagonal, Identity}, SolveKind::{LU, WAZ}, - UPLO::{Upper, Lower} }; +pub use crate::structure::matrix::UPLO::{Upper, Lower}; pub use crate::structure::dataframe::DType::*; pub use crate::structure::ad::AD::*; pub use crate::numerical::spline::SlopeMethod::{Akima, Quadratic}; diff --git a/src/numerical/newton.rs b/src/numerical/newton.rs index f05998bc..afef9890 100644 --- a/src/numerical/newton.rs +++ b/src/numerical/newton.rs @@ -1,9 +1,9 @@ use crate::numerical::utils::jacobian; -use crate::structure::matrix::*; use crate::structure::ad::*; use crate::traits::{ math::{Norm, Normed, Vector}, mutable::MutFP, + matrix::LinearAlgebra, }; /// Newton-Raphson Method diff --git a/src/numerical/optimize.rs b/src/numerical/optimize.rs index 765fc5d2..5cf4f22c 100644 --- a/src/numerical/optimize.rs +++ b/src/numerical/optimize.rs @@ -111,8 +111,8 @@ pub use self::OptMethod::{GaussNewton, GradientDescent, LevenbergMarquardt}; use self::OptOption::{InitParam, MaxIter}; use crate::numerical::utils::jacobian; -use crate::structure::matrix::{LinearAlgebra, Matrix}; -use crate::traits::matrix::MatrixTrait; +use crate::structure::matrix::Matrix; +use crate::traits::matrix::{MatrixTrait, LinearAlgebra}; use crate::structure::ad::{AD, ADVec}; use crate::util::useful::max; use std::collections::HashMap; diff --git a/src/numerical/spline.rs b/src/numerical/spline.rs index c1da0df1..bd7c0dcd 100644 --- a/src/numerical/spline.rs +++ b/src/numerical/spline.rs @@ -278,6 +278,7 @@ use crate::structure::matrix::*; use crate::structure::polynomial::*; #[allow(unused_imports)] use crate::structure::vector::*; +use crate::traits::matrix::LinearAlgebra; #[allow(unused_imports)] use crate::util::non_macro::*; use crate::util::useful::zip_range; diff --git a/src/prelude/mod.rs b/src/prelude/mod.rs index d111b760..f1abc275 100644 --- a/src/prelude/mod.rs +++ b/src/prelude/mod.rs @@ -153,7 +153,7 @@ pub use crate::traits::{ fp::{FPMatrix, FPVector}, general::Algorithm, math::{InnerProduct, LinearOp, MatrixProduct, Vector, VectorProduct}, - matrix::MatrixTrait, + matrix::{MatrixTrait, PQLU, QR, WAZD}, mutable::{MutFP, MutMatrix}, num::Real, pointer::{MatrixPtr, Oxide, Redox, RedoxCommon}, @@ -170,7 +170,7 @@ pub use crate::structure::{ ad::AD::*, matrix::{ combine, diag, gemm, gemv, gen_householder, inv_l, inv_u, matrix, ml_matrix, py_matrix, - r_matrix, Col, Matrix, Row, Shape, PQLU, QR, WAZD, + r_matrix, Col, Matrix, Row, Shape, }, polynomial::{Polynomial,poly,Calculus,lagrange_polynomial,legendre_polynomial}, vector::*, diff --git a/src/prelude/simpler.rs b/src/prelude/simpler.rs index c2c0c7d7..2916443b 100644 --- a/src/prelude/simpler.rs +++ b/src/prelude/simpler.rs @@ -8,9 +8,10 @@ use crate::numerical::{ }; #[cfg(feature = "parquet")] use crate::structure::dataframe::{DataFrame, WithParquet}; -use crate::structure::matrix::{self, Matrix}; +use crate::structure::matrix::Matrix; use crate::structure::polynomial; use crate::traits::math::{Norm, Normed}; +use crate::traits::matrix::{MatrixTrait, LinearAlgebra, PQLU, WAZD, QR, Form, SolveKind}; #[cfg(feature = "parquet")] use arrow2::io::parquet::write::CompressionOptions; #[cfg(feature = "parquet")] @@ -28,22 +29,22 @@ pub fn integrate f64 + Copy>(f: F, (a, b): (f64, f64)) -> f64 { } /// Simple Linear algebra -pub trait SimplerLinearAlgebra { - fn back_subs(&self, b: &Vec) -> Vec; - fn forward_subs(&self, b: &Vec) -> Vec; - fn lu(&self) -> matrix::PQLU; - fn waz_diag(&self) -> Option; - fn waz(&self) -> Option; - fn qr(&self) -> matrix::QR; +pub trait SimplerLinearAlgebra { + fn back_subs(&self, b: &[f64]) -> Vec; + fn forward_subs(&self, b: &[f64]) -> Vec; + fn lu(&self) -> PQLU; + fn waz_diag(&self) -> Option>; + fn waz(&self) -> Option>; + fn qr(&self) -> QR; #[cfg(feature = "O3")] - fn cholesky(&self) -> Matrix; - fn rref(&self) -> Matrix; + fn cholesky(&self) -> M; + fn rref(&self) -> M; fn det(&self) -> f64; - fn block(&self) -> (Matrix, Matrix, Matrix, Matrix); - fn inv(&self) -> Matrix; - fn pseudo_inv(&self) -> Matrix; - fn solve(&self, b: &Vec) -> Vec; - fn solve_mat(&self, m: &Matrix) -> Matrix; + fn block(&self) -> (M, M, M, M); + fn inv(&self) -> M; + fn pseudo_inv(&self) -> M; + fn solve(&self, b: &[f64]) -> Vec; + fn solve_mat(&self, m: &M) -> M; fn is_symmetric(&self) -> bool; } @@ -74,73 +75,73 @@ impl SimpleNorm for Matrix { } } -impl SimplerLinearAlgebra for Matrix { - fn back_subs(&self, b: &Vec) -> Vec { - matrix::LinearAlgebra::back_subs(self, b) +impl SimplerLinearAlgebra for Matrix { + fn back_subs(&self, b: &[f64]) -> Vec { + LinearAlgebra::back_subs(self, b) } - fn forward_subs(&self, b: &Vec) -> Vec { - matrix::LinearAlgebra::forward_subs(self, b) + fn forward_subs(&self, b: &[f64]) -> Vec { + LinearAlgebra::forward_subs(self, b) } - fn lu(&self) -> matrix::PQLU { - matrix::LinearAlgebra::lu(self) + fn lu(&self) -> PQLU { + LinearAlgebra::lu(self) } - fn waz_diag(&self) -> Option { - matrix::LinearAlgebra::waz(self, matrix::Form::Diagonal) + fn waz_diag(&self) -> Option> { + LinearAlgebra::waz(self, Form::Diagonal) } - fn waz(&self) -> Option { - matrix::LinearAlgebra::waz(self, matrix::Form::Identity) + fn waz(&self) -> Option> { + LinearAlgebra::waz(self, Form::Identity) } - fn qr(&self) -> matrix::QR { - matrix::LinearAlgebra::qr(self) + fn qr(&self) -> QR { + LinearAlgebra::qr(self) } #[cfg(feature = "O3")] fn cholesky(&self) -> Matrix { - matrix::LinearAlgebra::cholesky(self, matrix::UPLO::Lower) + LinearAlgebra::cholesky(self, UPLO::Lower) } fn rref(&self) -> Matrix { - matrix::LinearAlgebra::rref(self) + LinearAlgebra::rref(self) } fn det(&self) -> f64 { - matrix::LinearAlgebra::det(self) + LinearAlgebra::det(self) } fn block(&self) -> (Matrix, Matrix, Matrix, Matrix) { - matrix::LinearAlgebra::block(self) + LinearAlgebra::block(self) } fn inv(&self) -> Matrix { - matrix::LinearAlgebra::inv(self) + LinearAlgebra::inv(self) } fn pseudo_inv(&self) -> Matrix { - matrix::LinearAlgebra::pseudo_inv(self) + LinearAlgebra::pseudo_inv(self) } - fn solve(&self, b: &Vec) -> Vec { - matrix::LinearAlgebra::solve(self, b, matrix::SolveKind::LU) + fn solve(&self, b: &[f64]) -> Vec { + LinearAlgebra::solve(self, b, SolveKind::LU) } fn solve_mat(&self, m: &Matrix) -> Matrix { - matrix::LinearAlgebra::solve_mat(self, m, matrix::SolveKind::LU) + LinearAlgebra::solve_mat(self, m, SolveKind::LU) } fn is_symmetric(&self) -> bool { - matrix::LinearAlgebra::is_symmetric(self) + LinearAlgebra::is_symmetric(self) } } /// Simple solve #[allow(non_snake_case)] pub fn solve(A: &Matrix, m: &Matrix) -> Matrix { - matrix::solve(A, m, matrix::SolveKind::LU) + crate::traits::matrix::solve(A, m, SolveKind::LU) } /// Simple Chebyshev Polynomial (First Kind) diff --git a/src/statistics/stat.rs b/src/statistics/stat.rs index 41fb4dec..729af67c 100644 --- a/src/statistics/stat.rs +++ b/src/statistics/stat.rs @@ -151,7 +151,7 @@ use std::fmt; use self::QType::*; //use crate::structure::dataframe::*; use crate::structure::matrix::*; -use crate::traits::matrix::MatrixTrait; +use crate::traits::matrix::{MatrixTrait, LinearAlgebra}; use order_stat::kth_by; /// Statistics Trait diff --git a/src/structure/matrix.rs b/src/structure/matrix.rs index c4dbab42..23feb1f7 100644 --- a/src/structure/matrix.rs +++ b/src/structure/matrix.rs @@ -3049,8 +3049,8 @@ impl LinearAlgebra for Matrix { /// let (p,q,l,u) = (pqlu.p, pqlu.q, pqlu.l, pqlu.u); /// assert_eq!(p, vec![1]); // swap 0 & 1 (Row) /// assert_eq!(q, vec![1]); // swap 0 & 1 (Col) - /// assert_eq!(l, matrix(vec![1,0,0.5,1],2,2,Row)); - /// assert_eq!(u, matrix(vec![4,3,0,-0.5],2,2,Row)); + /// assert_eq!(l, matrix(vec![1.0,0.0,0.5,1.0],2,2,Row)); + /// assert_eq!(u, matrix(vec![4.0,3.0,0.0,-0.5],2,2,Row)); /// } /// ``` fn lu(&self) -> PQLU { @@ -3563,7 +3563,7 @@ impl LinearAlgebra for Matrix { SolveKind::LU => { let lu = self.lu(); let (p, q, l, u) = lu.extract(); - let mut v = b.clone(); + let mut v = b.to_vec(); v.swap_with_perm(&p.into_iter().enumerate().collect()); let z = l.forward_subs(&v); let mut y = u.back_subs(&z); @@ -3575,7 +3575,7 @@ impl LinearAlgebra for Matrix { None => panic!("Can't solve by WAZ with Singular matrix!"), Some(obj) => obj, }; - let x = &wazd.w.t() * b; + let x = &wazd.w.t() * &b.to_vec(); let x = &wazd.z * &x; x } diff --git a/src/structure/sparse.rs b/src/structure/sparse.rs index e640177f..c523953c 100644 --- a/src/structure/sparse.rs +++ b/src/structure/sparse.rs @@ -2,7 +2,8 @@ //! //! * Reference : Press, William H., and William T. Vetterling. *Numerical Recipes.* Cambridge: Cambridge Univ. Press, 2007. -use crate::structure::matrix::{Form, LinearAlgebra, Matrix, SolveKind, PQLU, QR, WAZD, SVD}; +use crate::structure::matrix::Matrix; +use crate::traits::matrix::{Form, LinearAlgebra, SolveKind, PQLU, QR, WAZD, SVD}; use crate::traits::math::LinearOp; //use crate::traits::math::{InnerProduct, LinearOp, Norm, Normed, Vector}; use crate::util::non_macro::zeros; @@ -138,28 +139,28 @@ impl LinearOp, Vec> for SPMatrix { /// Linear algebra for sparse matrix /// /// **Caution** : In every ops in this trait, there is converting process to dense matrix -impl LinearAlgebra for SPMatrix { - fn back_subs(&self, _b: &Vec) -> Vec { +impl LinearAlgebra for SPMatrix { + fn back_subs(&self, _b: &[f64]) -> Vec { unimplemented!() } - fn forward_subs(&self, _b: &Vec) -> Vec { + fn forward_subs(&self, _b: &[f64]) -> Vec { unimplemented!() } - fn lu(&self) -> PQLU { + fn lu(&self) -> PQLU { self.to_dense().lu() } - fn waz(&self, _d_form: Form) -> Option { + fn waz(&self, _d_form: Form) -> Option> { unimplemented!() } - fn qr(&self) -> QR { + fn qr(&self) -> QR { self.to_dense().qr() } - fn svd(&self) -> SVD { + fn svd(&self) -> SVD { unimplemented!() } @@ -188,7 +189,7 @@ impl LinearAlgebra for SPMatrix { self.to_dense().pseudo_inv() } - fn solve(&self, _b: &Vec, _sk: SolveKind) -> Vec { + fn solve(&self, _b: &[f64], _sk: SolveKind) -> Vec { unimplemented!() } From 34b6945deebd007d214345e70c5495a4bc118dbe Mon Sep 17 00:00:00 2001 From: Axect Date: Thu, 24 Oct 2024 11:30:00 +0900 Subject: [PATCH 4/5] FIX: Fix wrong implemented functions and apply updates for matrix trait --- src/complex/matrix.rs | 1022 +++++++++++++++++++---------------------- src/traits/matrix.rs | 2 +- 2 files changed, 473 insertions(+), 551 deletions(-) diff --git a/src/complex/matrix.rs b/src/complex/matrix.rs index 2d6eb3ec..a564c7bd 100644 --- a/src/complex/matrix.rs +++ b/src/complex/matrix.rs @@ -15,10 +15,12 @@ use crate::{ traits::fp::{FPMatrix, FPVector}, traits::general::Algorithm, traits::math::{InnerProduct, LinearOp, MatrixProduct, Norm, Normed, Vector}, + traits::matrix::{MatrixTrait, LinearAlgebra, Form, SolveKind, PQLU, WAZD, QR, SVD}, traits::mutable::MutMatrix, util::low_level::{copy_vec_ptr, swap_vec_ptr}, util::non_macro::ConcatenateError, util::useful::{nearly_eq, tab}, + complex::C64, }; /// R-like complex matrix structure @@ -27,15 +29,14 @@ use crate::{ /// /// ```rust /// use peroxide::fuga::*; -/// use num_complex::Complex64; /// use peroxide::complex::matrix::ComplexMatrix; /// /// let v1 = ComplexMatrix { /// data: vec![ -/// Complex64::new(1f64, 1f64), -/// Complex64::new(2f64, 2f64), -/// Complex64::new(3f64, 3f64), -/// Complex64::new(4f64, 4f64), +/// C64::new(1f64, 1f64), +/// C64::new(2f64, 2f64), +/// C64::new(3f64, 3f64), +/// C64::new(4f64, 4f64), /// ], /// row: 2, /// col: 2, @@ -44,7 +45,7 @@ use crate::{ /// ``` #[derive(Debug, Clone, Default)] pub struct ComplexMatrix { - pub data: Vec>, + pub data: Vec, pub row: usize, pub col: usize, pub shape: Shape, @@ -61,28 +62,27 @@ pub struct ComplexMatrix { /// #[macro_use] /// extern crate peroxide; /// use peroxide::fuga::*; -/// use num_complex::Complex64; -/// use peroxide::complex::matrix::complex_matrix; +/// use peroxide::complex::matrix::cmatrix; /// /// fn main() { -/// let a = complex_matrix(vec![Complex64::new(1f64, 1f64), -/// Complex64::new(2f64, 2f64), -/// Complex64::new(3f64, 3f64), -/// Complex64::new(4f64, 4f64)], +/// let a = cmatrix(vec![C64::new(1f64, 1f64), +/// C64::new(2f64, 2f64), +/// C64::new(3f64, 3f64), +/// C64::new(4f64, 4f64)], /// 2, 2, Row /// ); /// a.col.print(); // Print matrix column /// } /// ``` -pub fn complex_matrix(v: Vec, r: usize, c: usize, shape: Shape) -> ComplexMatrix +pub fn cmatrix(v: Vec, r: usize, c: usize, shape: Shape) -> ComplexMatrix where - T: Into>, + T: Into, { ComplexMatrix { data: v .into_iter() .map(|t| t.into()) - .collect::>>(), + .collect::>(), row: r, col: c, shape, @@ -90,11 +90,11 @@ where } /// R-like complex matrix constructor (Explicit ver.) -pub fn r_complex_matrix(v: Vec, r: usize, c: usize, shape: Shape) -> ComplexMatrix +pub fn r_cmatrix(v: Vec, r: usize, c: usize, shape: Shape) -> ComplexMatrix where - T: Into>, + T: Into, { - complex_matrix(v, r, c, shape) + cmatrix(v, r, c, shape) } /// Python-like complex matrix constructor @@ -104,37 +104,36 @@ where /// #[macro_use] /// extern crate peroxide; /// use peroxide::fuga::*; -/// use num_complex::Complex64; /// use peroxide::complex::matrix::*; /// /// fn main() { -/// let a = py_complex_matrix(vec![vec![Complex64::new(1f64, 1f64), -/// Complex64::new(2f64, 2f64)], -/// vec![Complex64::new(3f64, 3f64), -/// Complex64::new(4f64, 4f64)] +/// let a = py_cmatrix(vec![vec![C64::new(1f64, 1f64), +/// C64::new(2f64, 2f64)], +/// vec![C64::new(3f64, 3f64), +/// C64::new(4f64, 4f64)] /// ]); -/// let b = complex_matrix(vec![Complex64::new(1f64, 1f64), -/// Complex64::new(2f64, 2f64), -/// Complex64::new(3f64, 3f64), -/// Complex64::new(4f64, 4f64)], +/// let b = cmatrix(vec![C64::new(1f64, 1f64), +/// C64::new(2f64, 2f64), +/// C64::new(3f64, 3f64), +/// C64::new(4f64, 4f64)], /// 2, 2, Row /// ); /// assert_eq!(a, b); /// } /// ``` -pub fn py_complex_matrix(v: Vec>) -> ComplexMatrix +pub fn py_cmatrix(v: Vec>) -> ComplexMatrix where - T: Into> + Copy, + T: Into + Copy, { let r = v.len(); let c = v[0].len(); let data: Vec = v.into_iter().flatten().collect(); - complex_matrix(data, r, c, Shape::Row) + cmatrix(data, r, c, Shape::Row) } /// Matlab-like matrix constructor /// -/// Note that the entries to the `ml_complex_matrix` +/// Note that the entries to the `ml_cmatrix` /// needs to be in the `a+bi` format /// without any spaces between the real and imaginary /// parts of the Complex number. @@ -144,22 +143,21 @@ where /// #[macro_use] /// extern crate peroxide; /// use peroxide::fuga::*; -/// use num_complex::Complex64; /// use peroxide::complex::matrix::*; /// /// fn main() { -/// let a = ml_complex_matrix("1.0+1.0i 2.0+2.0i; +/// let a = ml_cmatrix("1.0+1.0i 2.0+2.0i; /// 3.0+3.0i 4.0+4.0i"); -/// let b = complex_matrix(vec![Complex64::new(1f64, 1f64), -/// Complex64::new(2f64, 2f64), -/// Complex64::new(3f64, 3f64), -/// Complex64::new(4f64, 4f64)], +/// let b = cmatrix(vec![C64::new(1f64, 1f64), +/// C64::new(2f64, 2f64), +/// C64::new(3f64, 3f64), +/// C64::new(4f64, 4f64)], /// 2, 2, Row /// ); /// assert_eq!(a, b); /// } /// ``` -pub fn ml_complex_matrix(s: &str) -> ComplexMatrix { +pub fn ml_cmatrix(s: &str) -> ComplexMatrix { let str_row = s.split(";").collect::>(); let r = str_row.len(); let str_data = str_row @@ -171,12 +169,12 @@ pub fn ml_complex_matrix(s: &str) -> ComplexMatrix { .iter() .flat_map(|t| { t.iter() - .map(|x| x.parse::>().unwrap()) - .collect::>>() + .map(|x| x.parse::().unwrap()) + .collect::>() }) - .collect::>>(); + .collect::>(); - complex_matrix(data, r, c, Shape::Row) + cmatrix(data, r, c, Shape::Row) } /// Pretty Print @@ -202,15 +200,17 @@ impl PartialEq for ComplexMatrix { } } -impl ComplexMatrix { +impl MatrixTrait for ComplexMatrix { + type Scalar = C64; + /// Raw pointer for `self.data` - pub fn ptr(&self) -> *const Complex { - &self.data[0] as *const Complex + fn ptr(&self) -> *const C64 { + &self.data[0] as *const C64 } /// Raw mutable pointer for `self.data` - pub fn mut_ptr(&mut self) -> *mut Complex { - &mut self.data[0] as *mut Complex + fn mut_ptr(&mut self) -> *mut C64 { + &mut self.data[0] as *mut C64 } /// Slice of `self.data` @@ -218,22 +218,21 @@ impl ComplexMatrix { /// # Examples /// ```rust /// use peroxide::fuga::*; - /// use num_complex::Complex64; /// use peroxide::complex::matrix::*; /// - /// let a = complex_matrix(vec![Complex64::new(1f64, 1f64), - /// Complex64::new(2f64, 2f64), - /// Complex64::new(3f64, 3f64), - /// Complex64::new(4f64, 4f64)], + /// let a = cmatrix(vec![C64::new(1f64, 1f64), + /// C64::new(2f64, 2f64), + /// C64::new(3f64, 3f64), + /// C64::new(4f64, 4f64)], /// 2, 2, Row /// ); /// let b = a.as_slice(); - /// assert_eq!(b, &[Complex64::new(1f64, 1f64), - /// Complex64::new(2f64, 2f64), - /// Complex64::new(3f64, 3f64), - /// Complex64::new(4f64, 4f64)]); + /// assert_eq!(b, &[C64::new(1f64, 1f64), + /// C64::new(2f64, 2f64), + /// C64::new(3f64, 3f64), + /// C64::new(4f64, 4f64)]); /// ``` - pub fn as_slice(&self) -> &[Complex] { + fn as_slice(&self) -> &[C64] { &self.data[..] } @@ -242,28 +241,27 @@ impl ComplexMatrix { /// # Examples /// ```rust /// use peroxide::fuga::*; - /// use num_complex::Complex64; /// use peroxide::complex::matrix::*; /// - /// let mut a = complex_matrix(vec![Complex64::new(1f64, 1f64), - /// Complex64::new(2f64, 2f64), - /// Complex64::new(3f64, 3f64), - /// Complex64::new(4f64, 4f64)], + /// let mut a = cmatrix(vec![C64::new(1f64, 1f64), + /// C64::new(2f64, 2f64), + /// C64::new(3f64, 3f64), + /// C64::new(4f64, 4f64)], /// 2, 2, Row /// ); /// let mut b = a.as_mut_slice(); - /// b[1] = Complex64::new(5f64, 5f64); - /// assert_eq!(b, &[Complex64::new(1f64, 1f64), - /// Complex64::new(5f64, 5f64), - /// Complex64::new(3f64, 3f64), - /// Complex64::new(4f64, 4f64)]); - /// assert_eq!(a, complex_matrix(vec![Complex64::new(1f64, 1f64), - /// Complex64::new(5f64, 5f64), - /// Complex64::new(3f64, 3f64), - /// Complex64::new(4f64, 4f64)], + /// b[1] = C64::new(5f64, 5f64); + /// assert_eq!(b, &[C64::new(1f64, 1f64), + /// C64::new(5f64, 5f64), + /// C64::new(3f64, 3f64), + /// C64::new(4f64, 4f64)]); + /// assert_eq!(a, cmatrix(vec![C64::new(1f64, 1f64), + /// C64::new(5f64, 5f64), + /// C64::new(3f64, 3f64), + /// C64::new(4f64, 4f64)], /// 2, 2, Row)); /// ``` - pub fn as_mut_slice(&mut self) -> &mut [Complex] { + fn as_mut_slice(&mut self) -> &mut [C64] { &mut self.data[..] } @@ -274,25 +272,24 @@ impl ComplexMatrix { /// # Examples /// ```rust /// use peroxide::fuga::*; - /// use num_complex::Complex64; /// use peroxide::complex::matrix::*; /// - /// let mut a = complex_matrix(vec![Complex64::new(1f64, 1f64), - /// Complex64::new(2f64, 2f64), - /// Complex64::new(3f64, 3f64), - /// Complex64::new(4f64, 4f64)], + /// let mut a = cmatrix(vec![C64::new(1f64, 1f64), + /// C64::new(2f64, 2f64), + /// C64::new(3f64, 3f64), + /// C64::new(4f64, 4f64)], /// 2, 2, Row /// ); /// assert_eq!(a.shape, Row); /// let b = a.change_shape(); /// assert_eq!(b.shape, Col); /// ``` - pub fn change_shape(&self) -> Self { + fn change_shape(&self) -> Self { let r = self.row; let c = self.col; assert_eq!(r * c, self.data.len()); let l = r * c - 1; - let mut data: Vec> = self.data.clone(); + let mut data: Vec = self.data.clone(); let ref_data = &self.data; match self.shape { @@ -302,7 +299,7 @@ impl ComplexMatrix { data[i] = ref_data[s]; } data[l] = ref_data[l]; - complex_matrix(data, r, c, Shape::Col) + cmatrix(data, r, c, Shape::Col) } Shape::Col => { for i in 0..l { @@ -310,7 +307,7 @@ impl ComplexMatrix { data[i] = ref_data[s]; } data[l] = ref_data[l]; - complex_matrix(data, r, c, Shape::Row) + cmatrix(data, r, c, Shape::Row) } } } @@ -322,43 +319,43 @@ impl ComplexMatrix { /// # Examples /// ```rust /// use peroxide::fuga::*; - /// use num_complex::Complex64; /// use peroxide::complex::matrix::*; /// - /// let mut a = complex_matrix(vec![Complex64::new(1f64, 1f64), - /// Complex64::new(2f64, 2f64), - /// Complex64::new(3f64, 3f64), - /// Complex64::new(4f64, 4f64)], - /// 2, 2, Row - /// ); + /// let mut a = cmatrix(vec![ + /// C64::new(1f64, 1f64), + /// C64::new(2f64, 2f64), + /// C64::new(3f64, 3f64), + /// C64::new(4f64, 4f64) + /// ], + /// 2, 2, Row + /// ); /// assert_eq!(a.shape, Row); - /// let b = a.change_shape_mut(); - /// assert_eq!(b.shape, Col); + /// a.change_shape_mut(); + /// assert_eq!(a.shape, Col); /// ``` - pub fn change_shape_mut(&mut self) -> Self { + fn change_shape_mut(&mut self) { let r = self.row; let c = self.col; assert_eq!(r * c, self.data.len()); let l = r * c - 1; - let mut data: Vec> = self.data.clone(); - let ref_data = &self.data; + let ref_data = self.data.clone(); match self.shape { Shape::Row => { for i in 0..l { let s = (i * c) % l; - data[i] = ref_data[s]; + self.data[i] = ref_data[s]; } - data[l] = ref_data[l]; - complex_matrix(data, r, c, Shape::Col) + self.data[l] = ref_data[l]; + self.shape = Shape::Col; } Shape::Col => { for i in 0..l { let s = (i * r) % l; - data[i] = ref_data[s]; + self.data[i] = ref_data[s]; } - data[l] = ref_data[l]; - complex_matrix(data, r, c, Shape::Row) + self.data[l] = ref_data[l]; + self.shape = Shape::Row; } } } @@ -368,13 +365,12 @@ impl ComplexMatrix { /// # Examples /// ```rust /// use peroxide::fuga::*; - /// use num_complex::Complex64; /// use peroxide::complex::matrix::*; /// - /// let a = complex_matrix(vec![Complex64::new(1f64, 1f64), - /// Complex64::new(2f64, 2f64), - /// Complex64::new(3f64, 3f64), - /// Complex64::new(4f64, 4f64)], + /// let a = cmatrix(vec![C64::new(1f64, 1f64), + /// C64::new(2f64, 2f64), + /// C64::new(3f64, 3f64), + /// C64::new(4f64, 4f64)], /// 2, 2, Row /// ); /// println!("{}", a.spread()); // same as println!("{}", a); @@ -383,7 +379,7 @@ impl ComplexMatrix { /// // r[0] 1+1i 3+3i /// // r[1] 2+2i 4+4i /// ``` - pub fn spread(&self) -> String { + fn spread(&self) -> String { assert_eq!(self.row * self.col, self.data.len()); let r = self.row; let c = self.col; @@ -464,22 +460,21 @@ impl ComplexMatrix { /// #[macro_use] /// extern crate peroxide; /// use peroxide::fuga::*; - /// use num_complex::Complex64; /// use peroxide::complex::matrix::*; /// /// fn main() { - /// let a = complex_matrix(vec![Complex64::new(1f64, 1f64), - /// Complex64::new(2f64, 2f64), - /// Complex64::new(3f64, 3f64), - /// Complex64::new(4f64, 4f64)], + /// let a = cmatrix(vec![C64::new(1f64, 1f64), + /// C64::new(2f64, 2f64), + /// C64::new(3f64, 3f64), + /// C64::new(4f64, 4f64)], /// 2, 2, Row /// ); - /// assert_eq!(a.col(0), vec![Complex64::new(1f64, 1f64), Complex64::new(3f64, 3f64)]); + /// assert_eq!(a.col(0), vec![C64::new(1f64, 1f64), C64::new(3f64, 3f64)]); /// } /// ``` - pub fn col(&self, index: usize) -> Vec> { + fn col(&self, index: usize) -> Vec { assert!(index < self.col); - let mut container: Vec> = vec![Complex::zero(); self.row]; + let mut container: Vec = vec![Complex::zero(); self.row]; for i in 0..self.row { container[i] = self[(i, index)]; } @@ -493,22 +488,21 @@ impl ComplexMatrix { /// #[macro_use] /// extern crate peroxide; /// use peroxide::fuga::*; - /// use num_complex::Complex64; /// use peroxide::complex::matrix::*; /// /// fn main() { - /// let a = complex_matrix(vec![Complex64::new(1f64, 1f64), - /// Complex64::new(2f64, 2f64), - /// Complex64::new(3f64, 3f64), - /// Complex64::new(4f64, 4f64)], + /// let a = cmatrix(vec![C64::new(1f64, 1f64), + /// C64::new(2f64, 2f64), + /// C64::new(3f64, 3f64), + /// C64::new(4f64, 4f64)], /// 2, 2, Row /// ); - /// assert_eq!(a.row(0), vec![Complex64::new(1f64, 1f64), Complex64::new(2f64, 2f64)]); + /// assert_eq!(a.row(0), vec![C64::new(1f64, 1f64), C64::new(2f64, 2f64)]); /// } /// ``` - pub fn row(&self, index: usize) -> Vec> { + fn row(&self, index: usize) -> Vec { assert!(index < self.row); - let mut container: Vec> = vec![Complex::zero(); self.col]; + let mut container: Vec = vec![Complex::zero(); self.col]; for i in 0..self.col { container[i] = self[(index, i)]; } @@ -522,20 +516,19 @@ impl ComplexMatrix { /// #[macro_use] /// extern crate peroxide; /// use peroxide::fuga::*; - /// use num_complex::Complex64; /// use peroxide::complex::matrix::*; /// /// fn main() { - /// let a = complex_matrix(vec![Complex64::new(1f64, 1f64), - /// Complex64::new(2f64, 2f64), - /// Complex64::new(3f64, 3f64), - /// Complex64::new(4f64, 4f64)], + /// let a = cmatrix(vec![C64::new(1f64, 1f64), + /// C64::new(2f64, 2f64), + /// C64::new(3f64, 3f64), + /// C64::new(4f64, 4f64)], /// 2, 2, Row /// ); - /// assert_eq!(a.diag(), vec![Complex64::new(1f64, 1f64) ,Complex64::new(4f64, 4f64)]); + /// assert_eq!(a.diag(), vec![C64::new(1f64, 1f64) ,C64::new(4f64, 4f64)]); /// } /// ``` - pub fn diag(&self) -> Vec> { + fn diag(&self) -> Vec { let mut container = vec![Complex::zero(); self.row]; let r = self.row; let c = self.col; @@ -553,80 +546,33 @@ impl ComplexMatrix { /// # Examples /// ```rust /// use peroxide::fuga::*; - /// use num_complex::Complex64; /// use peroxide::complex::matrix::*; /// - /// let a = complex_matrix(vec![Complex64::new(1f64, 1f64), - /// Complex64::new(2f64, 2f64), - /// Complex64::new(3f64, 3f64), - /// Complex64::new(4f64, 4f64)], + /// let a = cmatrix(vec![C64::new(1f64, 1f64), + /// C64::new(2f64, 2f64), + /// C64::new(3f64, 3f64), + /// C64::new(4f64, 4f64)], /// 2, 2, Row /// ); - /// let a_t = complex_matrix(vec![Complex64::new(1f64, 1f64), - /// Complex64::new(2f64, 2f64), - /// Complex64::new(3f64, 3f64), - /// Complex64::new(4f64, 4f64)], + /// let a_t = cmatrix(vec![C64::new(1f64, 1f64), + /// C64::new(2f64, 2f64), + /// C64::new(3f64, 3f64), + /// C64::new(4f64, 4f64)], /// 2, 2, Col /// ); /// /// assert_eq!(a.transpose(), a_t); /// ``` - pub fn transpose(&self) -> Self { + fn transpose(&self) -> Self { match self.shape { - Shape::Row => complex_matrix(self.data.clone(), self.col, self.row, Shape::Col), - Shape::Col => complex_matrix(self.data.clone(), self.col, self.row, Shape::Row), - } - } - - /// R-like transpose function - /// - /// # Examples - /// ```rust - /// use peroxide::fuga::*; - /// use num_complex::Complex64; - /// use peroxide::complex::matrix::*; - /// - /// let a = complex_matrix(vec![Complex64::new(1f64, 1f64), - /// Complex64::new(2f64, 2f64), - /// Complex64::new(3f64, 3f64), - /// Complex64::new(4f64, 4f64)], - /// 2, 2, Row - /// ); - /// let a_t = complex_matrix(vec![Complex64::new(1f64, 1f64), - /// Complex64::new(2f64, 2f64), - /// Complex64::new(3f64, 3f64), - /// Complex64::new(4f64, 4f64)], - /// 2, 2, Col - /// ); - /// - /// assert_eq!(a.t(), a_t); - /// ``` - pub fn t(&self) -> Self { - self.transpose() - } - - /// Should check shape - pub fn subs(&mut self, idx: usize, v: &Vec>) { - let p = &mut self.mut_ptr(); - match self.shape { - Shape::Row => { - let c = self.col; - unsafe { - p.add(idx * c).copy_from(v.as_ptr(), c); - } - } - Shape::Col => { - let r = self.row; - unsafe { - p.add(idx * r).copy_from(v.as_ptr(), r); - } - } + Shape::Row => cmatrix(self.data.clone(), self.col, self.row, Shape::Col), + Shape::Col => cmatrix(self.data.clone(), self.col, self.row, Shape::Row), } } /// Substitute Col #[inline] - pub fn subs_col(&mut self, idx: usize, v: &Vec>) { + fn subs_col(&mut self, idx: usize, v: &[C64]) { for i in 0..self.row { self[(i, idx)] = v[i]; } @@ -634,22 +580,22 @@ impl ComplexMatrix { /// Substitute Row #[inline] - pub fn subs_row(&mut self, idx: usize, v: &Vec>) { + fn subs_row(&mut self, idx: usize, v: &[C64]) { for j in 0..self.col { self[(idx, j)] = v[j]; } } /// From index operations - pub fn from_index(f: F, size: (usize, usize)) -> ComplexMatrix + fn from_index(f: F, size: (usize, usize)) -> ComplexMatrix where F: Fn(usize, usize) -> G + Copy, - G: Into>, + G: Into, { let row = size.0; let col = size.1; - let mut mat = complex_matrix(vec![Complex::zero(); row * col], row, col, Shape::Row); + let mut mat = cmatrix(vec![Complex::zero(); row * col], row, col, Shape::Row); for i in 0..row { for j in 0..col { @@ -659,10 +605,10 @@ impl ComplexMatrix { mat } - /// Matrix to `Vec>>` + /// Matrix to `Vec>` /// /// To send `Matrix` to `inline-python` - pub fn to_vec(&self) -> Vec>> { + fn to_vec(&self) -> Vec> { let mut result = vec![vec![Complex::zero(); self.col]; self.row]; for i in 0..self.row { result[i] = self.row(i); @@ -670,9 +616,9 @@ impl ComplexMatrix { result } - pub fn to_diag(&self) -> ComplexMatrix { + fn to_diag(&self) -> ComplexMatrix { assert_eq!(self.row, self.col, "Should be square matrix"); - let mut result = complex_matrix( + let mut result = cmatrix( vec![Complex::zero(); self.row * self.col], self.row, self.col, @@ -705,27 +651,26 @@ impl ComplexMatrix { /// #[macro_use] /// extern crate peroxide; /// use peroxide::fuga::*; - /// use num_complex::Complex64; /// use peroxide::complex::matrix::*; /// /// fn main() { - /// let a = ml_complex_matrix("1.0+1.0i 2.0+2.0i 3.0+3.0i; + /// let a = ml_cmatrix("1.0+1.0i 2.0+2.0i 3.0+3.0i; /// 4.0+4.0i 5.0+5.0i 6.0+6.0i; /// 7.0+7.0i 8.0+8.0i 9.0+9.0i"); - /// let b = complex_matrix(vec![Complex64::new(5f64, 5f64), - /// Complex64::new(6f64, 6f64), - /// Complex64::new(8f64, 8f64), - /// Complex64::new(9f64, 9f64)], + /// let b = cmatrix(vec![C64::new(5f64, 5f64), + /// C64::new(6f64, 6f64), + /// C64::new(8f64, 8f64), + /// C64::new(9f64, 9f64)], /// 2, 2, Row /// ); /// let c = a.submat((1, 1), (2, 2)); /// assert_eq!(b, c); /// } /// ``` - pub fn submat(&self, start: (usize, usize), end: (usize, usize)) -> ComplexMatrix { + fn submat(&self, start: (usize, usize), end: (usize, usize)) -> ComplexMatrix { let row = end.0 - start.0 + 1; let col = end.1 - start.1 + 1; - let mut result = complex_matrix(vec![Complex::zero(); row * col], row, col, self.shape); + let mut result = cmatrix(vec![Complex::zero(); row * col], row, col, self.shape); for i in 0..row { for j in 0..col { result[(i, j)] = self[(start.0 + i, start.1 + j)]; @@ -753,26 +698,25 @@ impl ComplexMatrix { /// ``` /// extern crate peroxide; /// use peroxide::fuga::*; - /// use num_complex::Complex64; /// use peroxide::complex::matrix::*; /// /// fn main() { - /// let mut a = ml_complex_matrix("1.0+1.0i 2.0+2.0i 3.0+3.0i; + /// let mut a = ml_cmatrix("1.0+1.0i 2.0+2.0i 3.0+3.0i; /// 4.0+4.0i 5.0+5.0i 6.0+6.0i; /// 7.0+7.0i 8.0+8.0i 9.0+9.0i"); - /// let b = complex_matrix(vec![Complex64::new(1f64, 1f64), - /// Complex64::new(2f64, 2f64), - /// Complex64::new(3f64, 3f64), - /// Complex64::new(4f64, 4f64)], + /// let b = cmatrix(vec![C64::new(1f64, 1f64), + /// C64::new(2f64, 2f64), + /// C64::new(3f64, 3f64), + /// C64::new(4f64, 4f64)], /// 2, 2, Row); - /// let c = ml_complex_matrix("1.0+1.0i 2.0+2.0i 3.0+3.0i; + /// let c = ml_cmatrix("1.0+1.0i 2.0+2.0i 3.0+3.0i; /// 4.0+4.0i 1.0+1.0i 2.0+2.0i; /// 7.0+7.0i 3.0+3.0i 4.0+4.0i"); /// a.subs_mat((1,1), (2,2), &b); /// assert_eq!(a, c); /// } /// ``` - pub fn subs_mat(&mut self, start: (usize, usize), end: (usize, usize), m: &ComplexMatrix) { + fn subs_mat(&mut self, start: (usize, usize), end: (usize, usize), m: &ComplexMatrix) { let row = end.0 - start.0 + 1; let col = end.1 - start.1 + 1; for i in 0..row { @@ -787,13 +731,13 @@ impl ComplexMatrix { // Mathematics for Matrix // ============================================================================= impl Vector for ComplexMatrix { - type Scalar = Complex; + type Scalar = C64; fn add_vec(&self, other: &Self) -> Self { assert_eq!(self.row, other.row); assert_eq!(self.col, other.col); - let mut result = complex_matrix(self.data.clone(), self.row, self.col, self.shape); + let mut result = cmatrix(self.data.clone(), self.row, self.col, self.shape); for i in 0..self.row { for j in 0..self.col { result[(i, j)] += other[(i, j)]; @@ -806,7 +750,7 @@ impl Vector for ComplexMatrix { assert_eq!(self.row, other.row); assert_eq!(self.col, other.col); - let mut result = complex_matrix(self.data.clone(), self.row, self.col, self.shape); + let mut result = cmatrix(self.data.clone(), self.row, self.col, self.shape); for i in 0..self.row { for j in 0..self.col { result[(i, j)] -= other[(i, j)]; @@ -850,7 +794,7 @@ impl Normed for ComplexMatrix { Shape::Row => self.change_shape().norm(Norm::L1), Shape::Col => { for c in 0..self.col { - let s: Complex = self.col(c).iter().sum(); + let s: C64 = self.col(c).iter().sum(); if s.re > m.re { m = s; } @@ -865,7 +809,7 @@ impl Normed for ComplexMatrix { Shape::Col => self.change_shape().norm(Norm::LInf), Shape::Row => { for r in 0..self.row { - let s: Complex = self.row(r).iter().sum(); + let s: C64 = self.row(r).iter().sum(); if s.re > m.re { m = s; } @@ -891,7 +835,7 @@ impl Normed for ComplexMatrix { /// Frobenius inner product impl InnerProduct for ComplexMatrix { - fn dot(&self, rhs: &Self) -> Complex { + fn dot(&self, rhs: &Self) -> C64 { if self.shape == rhs.shape { self.data.dot(&rhs.data) } else { @@ -904,11 +848,11 @@ impl InnerProduct for ComplexMatrix { /// Matrix as Linear operator for Vector #[allow(non_snake_case)] -impl LinearOp>, Vec>> for ComplexMatrix { - fn apply(&self, other: &Vec>) -> Vec> { +impl LinearOp, Vec> for ComplexMatrix { + fn apply(&self, other: &Vec) -> Vec { assert_eq!(self.col, other.len()); let mut c = vec![Complex::zero(); self.row]; - complex_gemv(Complex::one(), self, other, Complex::zero(), &mut c); + cgemv(Complex::one(), self, other, Complex::zero(), &mut c); c } } @@ -934,7 +878,7 @@ pub fn complex_cbind(m1: ComplexMatrix, m2: ComplexMatrix) -> Result Result // ============================================================================= -/// `Complex Matrix` to `Vec>` -impl Into>> for ComplexMatrix { - fn into(self) -> Vec> { +/// `Complex Matrix` to `Vec` +impl Into> for ComplexMatrix { + fn into(self) -> Vec { self.data } } -/// `&ComplexMatrix` to `&Vec>` -impl<'a> Into<&'a Vec>> for &'a ComplexMatrix { - fn into(self) -> &'a Vec> { +/// `&ComplexMatrix` to `&Vec` +impl<'a> Into<&'a Vec> for &'a ComplexMatrix { + fn into(self) -> &'a Vec { &self.data } } -/// `Vec>` to `ComplexMatrix` -impl Into for Vec> { +/// `Vec` to `ComplexMatrix` +impl Into for Vec { fn into(self) -> ComplexMatrix { let l = self.len(); - complex_matrix(self, l, 1, Shape::Col) + cmatrix(self, l, 1, Shape::Col) } } -/// `&Vec>` to `ComplexMatrix` -impl Into for &Vec> { +/// `&Vec` to `ComplexMatrix` +impl Into for &Vec { fn into(self) -> ComplexMatrix { let l = self.len(); - complex_matrix(self.clone(), l, 1, Shape::Col) + cmatrix(self.clone(), l, 1, Shape::Col) } } @@ -1051,7 +995,7 @@ impl Add for ComplexMatrix { assert_eq!(&self.row, &other.row); assert_eq!(&self.col, &other.col); - let mut result = complex_matrix(self.data.clone(), self.row, self.col, self.shape); + let mut result = cmatrix(self.data.clone(), self.row, self.col, self.shape); for i in 0..self.row { for j in 0..self.col { result[(i, j)] += other[(i, j)]; @@ -1069,27 +1013,26 @@ impl<'a, 'b> Add<&'b ComplexMatrix> for &'a ComplexMatrix { } } -/// Element-wise addition between Complex Matrix & Complex +/// Element-wise addition between Complex Matrix & C64 /// /// # Examples /// ```rust /// #[macro_use] /// extern crate peroxide; /// use peroxide::fuga::*; -/// use num_complex::Complex64; /// use peroxide::complex::matrix::*; /// /// fn main() { -/// let mut a = ml_complex_matrix("1.0+1.0i 2.0+2.0i; +/// let mut a = ml_cmatrix("1.0+1.0i 2.0+2.0i; /// 4.0+4.0i 5.0+5.0i"); -/// let a_exp = ml_complex_matrix("2.0+2.0i 3.0+3.0i; +/// let a_exp = ml_cmatrix("2.0+2.0i 3.0+3.0i; /// 5.0+5.0i 6.0+6.0i"); -/// assert_eq!(a + Complex64::new(1_f64, 1_f64), a_exp); +/// assert_eq!(a + C64::new(1_f64, 1_f64), a_exp); /// } /// ``` impl Add for ComplexMatrix where - T: Into> + Copy, + T: Into + Copy, { type Output = Self; fn add(self, other: T) -> Self { @@ -1097,10 +1040,10 @@ where } } -/// Element-wise addition between &ComplexMatrix & Complex +/// Element-wise addition between &ComplexMatrix & C64 impl<'a, T> Add for &'a ComplexMatrix where - T: Into> + Copy, + T: Into + Copy, { type Output = ComplexMatrix; @@ -1109,7 +1052,7 @@ where } } -// Element-wise addition between Complex & ComplexMatrix +// Element-wise addition between C64 & ComplexMatrix /// /// # Examples /// @@ -1117,18 +1060,17 @@ where /// #[macro_use] /// extern crate peroxide; /// use peroxide::fuga::*; -/// use num_complex::Complex64; /// use peroxide::complex::matrix::*; /// /// fn main() { -/// let mut a = ml_complex_matrix("1.0+1.0i 2.0+2.0i; +/// let mut a = ml_cmatrix("1.0+1.0i 2.0+2.0i; /// 4.0+4.0i 5.0+5.0i"); -/// let a_exp = ml_complex_matrix("2.0+2.0i 3.0+3.0i; +/// let a_exp = ml_cmatrix("2.0+2.0i 3.0+3.0i; /// 5.0+5.0i 6.0+6.0i"); -/// assert_eq!(Complex64::new(1_f64, 1_f64) + a, a_exp); +/// assert_eq!(C64::new(1_f64, 1_f64) + a, a_exp); /// } /// ``` -impl Add for Complex { +impl Add for C64 { type Output = ComplexMatrix; fn add(self, other: ComplexMatrix) -> Self::Output { @@ -1136,8 +1078,8 @@ impl Add for Complex { } } -/// Element-wise addition between Complex & &ComplexMatrix -impl<'a> Add<&'a ComplexMatrix> for Complex { +/// Element-wise addition between C64 & &ComplexMatrix +impl<'a> Add<&'a ComplexMatrix> for C64 { type Output = ComplexMatrix; fn add(self, other: &'a ComplexMatrix) -> Self::Output { @@ -1154,18 +1096,17 @@ impl<'a> Add<&'a ComplexMatrix> for Complex { /// ```rust /// extern crate peroxide; /// use peroxide::fuga::*; -/// use num_complex::Complex64; /// use peroxide::complex::matrix::*; /// -/// let a = complex_matrix(vec![Complex64::new(1f64, 1f64), -/// Complex64::new(2f64, 2f64), -/// Complex64::new(3f64, 3f64), -/// Complex64::new(4f64, 4f64)], +/// let a = cmatrix(vec![C64::new(1f64, 1f64), +/// C64::new(2f64, 2f64), +/// C64::new(3f64, 3f64), +/// C64::new(4f64, 4f64)], /// 2, 2, Row); -/// let a_neg = complex_matrix(vec![Complex64::new(-1f64, -1f64), -/// Complex64::new(-2f64, -2f64), -/// Complex64::new(-3f64, -3f64), -/// Complex64::new(-4f64, -4f64)], +/// let a_neg = cmatrix(vec![C64::new(-1f64, -1f64), +/// C64::new(-2f64, -2f64), +/// C64::new(-3f64, -3f64), +/// C64::new(-4f64, -4f64)], /// 2, 2, Row); /// assert_eq!(-a, a_neg); /// ``` @@ -1173,11 +1114,11 @@ impl Neg for ComplexMatrix { type Output = Self; fn neg(self) -> Self { - complex_matrix( + cmatrix( self.data .into_iter() - .map(|x: Complex| -x) - .collect::>>(), + .map(|x: C64| -x) + .collect::>(), self.row, self.col, self.shape, @@ -1190,12 +1131,12 @@ impl<'a> Neg for &'a ComplexMatrix { type Output = ComplexMatrix; fn neg(self) -> Self::Output { - complex_matrix( + cmatrix( self.data .clone() .into_iter() - .map(|x: Complex| -x) - .collect::>>(), + .map(|x: C64| -x) + .collect::>(), self.row, self.col, self.shape, @@ -1213,15 +1154,14 @@ impl<'a> Neg for &'a ComplexMatrix { /// #[macro_use] /// extern crate peroxide; /// use peroxide::fuga::*; -/// use num_complex::Complex64; /// use peroxide::complex::matrix::*; /// /// fn main() { -/// let a = ml_complex_matrix("10.0+10.0i 20.0+20.0i; +/// let a = ml_cmatrix("10.0+10.0i 20.0+20.0i; /// 40.0+40.0i 50.0+50.0i"); -/// let b = ml_complex_matrix("1.0+1.0i 2.0+2.0i; +/// let b = ml_cmatrix("1.0+1.0i 2.0+2.0i; /// 4.0+4.0i 5.0+5.0i"); -/// let diff = ml_complex_matrix("9.0+9.0i 18.0+18.0i; +/// let diff = ml_cmatrix("9.0+9.0i 18.0+18.0i; /// 36.0+36.0i 45.0+45.0i"); /// assert_eq!(a-b, diff); /// } @@ -1232,7 +1172,7 @@ impl Sub for ComplexMatrix { fn sub(self, other: Self) -> Self::Output { assert_eq!(&self.row, &other.row); assert_eq!(&self.col, &other.col); - let mut result = complex_matrix(self.data.clone(), self.row, self.col, self.shape); + let mut result = cmatrix(self.data.clone(), self.row, self.col, self.shape); for i in 0..self.row { for j in 0..self.col { result[(i, j)] -= other[(i, j)]; @@ -1250,10 +1190,10 @@ impl<'a, 'b> Sub<&'b ComplexMatrix> for &'a ComplexMatrix { } } -/// Subtraction between Complex Matrix & Complex +/// Subtraction between Complex Matrix & C64 impl Sub for ComplexMatrix where - T: Into> + Copy, + T: Into + Copy, { type Output = Self; @@ -1262,10 +1202,10 @@ where } } -/// Subtraction between &Complex Matrix & Complex +/// Subtraction between &Complex Matrix & C64 impl<'a, T> Sub for &'a ComplexMatrix where - T: Into> + Copy, + T: Into + Copy, { type Output = ComplexMatrix; @@ -1274,25 +1214,24 @@ where } } -/// Subtraction Complex Matrix with Complex +/// Subtraction Complex Matrix with C64 /// /// # Examples /// ```rust /// #[macro_use] /// extern crate peroxide; /// use peroxide::fuga::*; -/// use num_complex::Complex64; /// use peroxide::complex::matrix::*; /// /// fn main() { -/// let mut a = ml_complex_matrix("1.0+1.0i 2.0+2.0i; +/// let mut a = ml_cmatrix("1.0+1.0i 2.0+2.0i; /// 4.0+4.0i 5.0+5.0i"); -/// let a_exp = ml_complex_matrix("0.0+0.0i 1.0+1.0i; +/// let a_exp = ml_cmatrix("0.0+0.0i 1.0+1.0i; /// 3.0+3.0i 4.0+4.0i"); -/// assert_eq!(a - Complex64::new(1_f64, 1_f64), a_exp); +/// assert_eq!(a - C64::new(1_f64, 1_f64), a_exp); /// } /// ``` -impl Sub for Complex { +impl Sub for C64 { type Output = ComplexMatrix; fn sub(self, other: ComplexMatrix) -> Self::Output { @@ -1311,16 +1250,16 @@ impl<'a> Sub<&'a ComplexMatrix> for f64 { // ============================================================================= // Multiplication for Complex Matrix // ============================================================================= -/// Element-wise multiplication between Complex Matrix vs Complex -impl Mul> for ComplexMatrix { +/// Element-wise multiplication between Complex Matrix vs C64 +impl Mul for ComplexMatrix { type Output = Self; - fn mul(self, other: Complex) -> Self::Output { + fn mul(self, other: C64) -> Self::Output { self.fmap(|x| x * other) } } -impl Mul for Complex { +impl Mul for C64 { type Output = ComplexMatrix; fn mul(self, other: ComplexMatrix) -> Self::Output { @@ -1328,7 +1267,7 @@ impl Mul for Complex { } } -impl<'a> Mul<&'a ComplexMatrix> for Complex { +impl<'a> Mul<&'a ComplexMatrix> for C64 { type Output = ComplexMatrix; fn mul(self, other: &'a ComplexMatrix) -> Self::Output { @@ -1343,15 +1282,14 @@ impl<'a> Mul<&'a ComplexMatrix> for Complex { /// #[macro_use] /// extern crate peroxide; /// use peroxide::fuga::*; -/// use num_complex::Complex64; /// use peroxide::complex::matrix::*; /// /// fn main() { -/// let mut a = ml_complex_matrix("1.0+1.0i 2.0+2.0i; +/// let mut a = ml_cmatrix("1.0+1.0i 2.0+2.0i; /// 4.0+4.0i 5.0+5.0i"); -/// let mut b = ml_complex_matrix("2.0+2.0i 2.0+2.0i; +/// let mut b = ml_cmatrix("2.0+2.0i 2.0+2.0i; /// 5.0+5.0i 5.0+5.0i"); -/// let prod = ml_complex_matrix("0.0+24.0i 0.0+24.0i; +/// let prod = ml_cmatrix("0.0+24.0i 0.0+24.0i; /// 0.0+66.0i 0.0+66.0i"); /// assert_eq!(a * b, prod); /// } @@ -1360,7 +1298,7 @@ impl Mul for ComplexMatrix { type Output = Self; fn mul(self, other: Self) -> Self::Output { - matmul(&self, &other) + cmatmul(&self, &other) } } @@ -1368,31 +1306,31 @@ impl<'a, 'b> Mul<&'b ComplexMatrix> for &'a ComplexMatrix { type Output = ComplexMatrix; fn mul(self, other: &'b ComplexMatrix) -> Self::Output { - matmul(self, other) + cmatmul(self, other) } } #[allow(non_snake_case)] -impl Mul>> for ComplexMatrix { - type Output = Vec>; +impl Mul> for ComplexMatrix { + type Output = Vec; - fn mul(self, other: Vec>) -> Self::Output { + fn mul(self, other: Vec) -> Self::Output { self.apply(&other) } } #[allow(non_snake_case)] -impl<'a, 'b> Mul<&'b Vec>> for &'a ComplexMatrix { - type Output = Vec>; +impl<'a, 'b> Mul<&'b Vec> for &'a ComplexMatrix { + type Output = Vec; - fn mul(self, other: &'b Vec>) -> Self::Output { + fn mul(self, other: &'b Vec) -> Self::Output { self.apply(other) } } -/// Matrix multiplication for `Vec>` vs `ComplexMatrix` -impl Mul for Vec> { - type Output = Vec>; +/// Matrix multiplication for `Vec` vs `ComplexMatrix` +impl Mul for Vec { + type Output = Vec; fn mul(self, other: ComplexMatrix) -> Self::Output { assert_eq!(self.len(), other.row); @@ -1402,8 +1340,8 @@ impl Mul for Vec> { } } -impl<'a, 'b> Mul<&'b ComplexMatrix> for &'a Vec> { - type Output = Vec>; +impl<'a, 'b> Mul<&'b ComplexMatrix> for &'a Vec { + type Output = Vec; fn mul(self, other: &'b ComplexMatrix) -> Self::Output { assert_eq!(self.len(), other.row); @@ -1416,46 +1354,45 @@ impl<'a, 'b> Mul<&'b ComplexMatrix> for &'a Vec> { // ============================================================================= // Standard Operation for Matrix (DIV) // ============================================================================= -/// Element-wise division between Complex Matrix vs Complex -impl Div> for ComplexMatrix { +/// Element-wise division between Complex Matrix vs C64 +impl Div for ComplexMatrix { type Output = Self; - fn div(self, other: Complex) -> Self::Output { + fn div(self, other: C64) -> Self::Output { self.fmap(|x| x / other) } } -impl<'a> Div> for &'a ComplexMatrix { +impl<'a> Div for &'a ComplexMatrix { type Output = ComplexMatrix; - fn div(self, other: Complex) -> Self::Output { + fn div(self, other: C64) -> Self::Output { self.fmap(|x| x / other) } } /// Index for Complex Matrix /// -/// `(usize, usize) -> Complex` +/// `(usize, usize) -> C64` /// /// # Examples /// ```rust /// extern crate peroxide; /// use peroxide::fuga::*; -/// use num_complex::Complex64; /// use peroxide::complex::matrix::*; /// -/// let a = complex_matrix(vec![Complex64::new(1f64, 1f64), -/// Complex64::new(2f64, 2f64), -/// Complex64::new(3f64, 3f64), -/// Complex64::new(4f64, 4f64)], +/// let a = cmatrix(vec![C64::new(1f64, 1f64), +/// C64::new(2f64, 2f64), +/// C64::new(3f64, 3f64), +/// C64::new(4f64, 4f64)], /// 2, 2, Row /// ); -/// assert_eq!(a[(0,1)], Complex64::new(2f64, 2f64)); +/// assert_eq!(a[(0,1)], C64::new(2f64, 2f64)); /// ``` impl Index<(usize, usize)> for ComplexMatrix { - type Output = Complex; + type Output = C64; - fn index(&self, pair: (usize, usize)) -> &Complex { + fn index(&self, pair: (usize, usize)) -> &C64 { let p = self.ptr(); let i = pair.0; let j = pair.1; @@ -1469,25 +1406,24 @@ impl Index<(usize, usize)> for ComplexMatrix { /// IndexMut for Complex Matrix (Assign) /// -/// `(usize, usize) -> Complex` +/// `(usize, usize) -> C64` /// /// # Examples /// ```rust /// extern crate peroxide; /// use peroxide::fuga::*; -/// use num_complex::Complex64; /// use peroxide::complex::matrix::*; /// -/// let mut a = complex_matrix(vec![Complex64::new(1f64, 1f64), -/// Complex64::new(2f64, 2f64), -/// Complex64::new(3f64, 3f64), -/// Complex64::new(4f64, 4f64)], +/// let mut a = cmatrix(vec![C64::new(1f64, 1f64), +/// C64::new(2f64, 2f64), +/// C64::new(3f64, 3f64), +/// C64::new(4f64, 4f64)], /// 2, 2, Row /// ); -/// assert_eq!(a[(0,1)], Complex64::new(2f64, 2f64)); +/// assert_eq!(a[(0,1)], C64::new(2f64, 2f64)); /// ``` impl IndexMut<(usize, usize)> for ComplexMatrix { - fn index_mut(&mut self, pair: (usize, usize)) -> &mut Complex { + fn index_mut(&mut self, pair: (usize, usize)) -> &mut C64 { let i = pair.0; let j = pair.1; let r = self.row; @@ -1512,7 +1448,7 @@ impl IndexMut<(usize, usize)> for ComplexMatrix { // ============================================================================= impl FPMatrix for ComplexMatrix { - type Scalar = Complex; + type Scalar = C64; fn take_row(&self, n: usize) -> Self { if n >= self.row { @@ -1525,15 +1461,15 @@ impl FPMatrix for ComplexMatrix { .clone() .into_iter() .take(n * self.col) - .collect::>>(); - complex_matrix(new_data, n, self.col, Shape::Row) + .collect::>(); + cmatrix(new_data, n, self.col, Shape::Row) } Shape::Col => { - let mut temp_data: Vec> = Vec::new(); + let mut temp_data: Vec = Vec::new(); for i in 0..n { temp_data.extend(self.row(i)); } - complex_matrix(temp_data, n, self.col, Shape::Row) + cmatrix(temp_data, n, self.col, Shape::Row) } } } @@ -1549,15 +1485,15 @@ impl FPMatrix for ComplexMatrix { .clone() .into_iter() .take(n * self.row) - .collect::>>(); - complex_matrix(new_data, self.row, n, Shape::Col) + .collect::>(); + cmatrix(new_data, self.row, n, Shape::Col) } Shape::Row => { - let mut temp_data: Vec> = Vec::new(); + let mut temp_data: Vec = Vec::new(); for i in 0..n { temp_data.extend(self.col(i)); } - complex_matrix(temp_data, self.row, n, Shape::Col) + cmatrix(temp_data, self.row, n, Shape::Col) } } } @@ -1565,33 +1501,33 @@ impl FPMatrix for ComplexMatrix { fn skip_row(&self, n: usize) -> Self { assert!(n < self.row, "Skip range is larger than row of matrix"); - let mut temp_data: Vec> = Vec::new(); + let mut temp_data: Vec = Vec::new(); for i in n..self.row { temp_data.extend(self.row(i)); } - complex_matrix(temp_data, self.row - n, self.col, Shape::Row) + cmatrix(temp_data, self.row - n, self.col, Shape::Row) } fn skip_col(&self, n: usize) -> Self { assert!(n < self.col, "Skip range is larger than col of matrix"); - let mut temp_data: Vec> = Vec::new(); + let mut temp_data: Vec = Vec::new(); for i in n..self.col { temp_data.extend(self.col(i)); } - complex_matrix(temp_data, self.row, self.col - n, Shape::Col) + cmatrix(temp_data, self.row, self.col - n, Shape::Col) } fn fmap(&self, f: F) -> Self where - F: Fn(Complex) -> Complex, + F: Fn(C64) -> C64, { let result = self .data .iter() .map(|x| f(*x)) - .collect::>>(); - complex_matrix(result, self.row, self.col, self.shape) + .collect::>(); + cmatrix(result, self.row, self.col, self.shape) } /// Column map @@ -1599,23 +1535,22 @@ impl FPMatrix for ComplexMatrix { /// # Example /// ```rust /// use peroxide::fuga::*; - /// use num_complex::Complex64; /// use peroxide::complex::matrix::*; /// use peroxide::traits::fp::FPMatrix; /// /// fn main() { - /// let x = complex_matrix(vec![Complex64::new(1f64, 1f64), - /// Complex64::new(2f64, 2f64), - /// Complex64::new(3f64, 3f64), - /// Complex64::new(4f64, 4f64)], + /// let x = cmatrix(vec![C64::new(1f64, 1f64), + /// C64::new(2f64, 2f64), + /// C64::new(3f64, 3f64), + /// C64::new(4f64, 4f64)], /// 2, 2, Row /// ); /// let y = x.col_map(|r| r.fmap(|t| t + r[0])); /// - /// let y_col_map = complex_matrix(vec![Complex64::new(2f64, 2f64), - /// Complex64::new(4f64, 4f64), - /// Complex64::new(4f64, 4f64), - /// Complex64::new(6f64, 6f64)], + /// let y_col_map = cmatrix(vec![C64::new(2f64, 2f64), + /// C64::new(4f64, 4f64), + /// C64::new(4f64, 4f64), + /// C64::new(6f64, 6f64)], /// 2, 2, Col /// ); /// @@ -1624,9 +1559,9 @@ impl FPMatrix for ComplexMatrix { /// ``` fn col_map(&self, f: F) -> ComplexMatrix where - F: Fn(Vec>) -> Vec>, + F: Fn(Vec) -> Vec, { - let mut result = complex_matrix( + let mut result = cmatrix( vec![Complex::zero(); self.row * self.col], self.row, self.col, @@ -1645,23 +1580,22 @@ impl FPMatrix for ComplexMatrix { /// # Example /// ```rust /// use peroxide::fuga::*; - /// use num_complex::Complex64; /// use peroxide::complex::matrix::*; /// use peroxide::traits::fp::FPMatrix; /// /// fn main() { - /// let x = complex_matrix(vec![Complex64::new(1f64, 1f64), - /// Complex64::new(2f64, 2f64), - /// Complex64::new(3f64, 3f64), - /// Complex64::new(4f64, 4f64)], + /// let x = cmatrix(vec![C64::new(1f64, 1f64), + /// C64::new(2f64, 2f64), + /// C64::new(3f64, 3f64), + /// C64::new(4f64, 4f64)], /// 2, 2, Row /// ); /// let y = x.row_map(|r| r.fmap(|t| t + r[0])); /// - /// let y_row_map = complex_matrix(vec![Complex64::new(2f64, 2f64), - /// Complex64::new(3f64, 3f64), - /// Complex64::new(6f64, 6f64), - /// Complex64::new(7f64, 7f64)], + /// let y_row_map = cmatrix(vec![C64::new(2f64, 2f64), + /// C64::new(3f64, 3f64), + /// C64::new(6f64, 6f64), + /// C64::new(7f64, 7f64)], /// 2, 2, Row /// ); /// @@ -1670,9 +1604,9 @@ impl FPMatrix for ComplexMatrix { /// ``` fn row_map(&self, f: F) -> ComplexMatrix where - F: Fn(Vec>) -> Vec>, + F: Fn(Vec) -> Vec, { - let mut result = complex_matrix( + let mut result = cmatrix( vec![Complex::zero(); self.row * self.col], self.row, self.col, @@ -1688,7 +1622,7 @@ impl FPMatrix for ComplexMatrix { fn col_mut_map(&mut self, f: F) where - F: Fn(Vec>) -> Vec>, + F: Fn(Vec) -> Vec, { for i in 0..self.col { unsafe { @@ -1703,7 +1637,7 @@ impl FPMatrix for ComplexMatrix { fn row_mut_map(&mut self, f: F) where - F: Fn(Vec>) -> Vec>, + F: Fn(Vec) -> Vec, { for i in 0..self.col { unsafe { @@ -1716,17 +1650,17 @@ impl FPMatrix for ComplexMatrix { } } - fn reduce(&self, init: T, f: F) -> Complex + fn reduce(&self, init: T, f: F) -> C64 where - F: Fn(Complex, Complex) -> Complex, - T: Into>, + F: Fn(C64, C64) -> C64, + T: Into, { self.data.iter().fold(init.into(), |x, y| f(x, *y)) } fn zip_with(&self, f: F, other: &ComplexMatrix) -> Self where - F: Fn(Complex, Complex) -> Complex, + F: Fn(C64, C64) -> C64, { assert_eq!(self.data.len(), other.data.len()); let mut a = other.clone(); @@ -1738,13 +1672,13 @@ impl FPMatrix for ComplexMatrix { .iter() .zip(a.data.iter()) .map(|(x, y)| f(*x, *y)) - .collect::>>(); - complex_matrix(result, self.row, self.col, self.shape) + .collect::>(); + cmatrix(result, self.row, self.col, self.shape) } - fn col_reduce(&self, f: F) -> Vec> + fn col_reduce(&self, f: F) -> Vec where - F: Fn(Vec>) -> Complex, + F: Fn(Vec) -> C64, { let mut v = vec![Complex::zero(); self.col]; for i in 0..self.col { @@ -1753,9 +1687,9 @@ impl FPMatrix for ComplexMatrix { v } - fn row_reduce(&self, f: F) -> Vec> + fn row_reduce(&self, f: F) -> Vec where - F: Fn(Vec>) -> Complex, + F: Fn(Vec) -> C64, { let mut v = vec![Complex::zero(); self.row]; for i in 0..self.row { @@ -1765,47 +1699,36 @@ impl FPMatrix for ComplexMatrix { } } -pub fn diag(n: usize) -> ComplexMatrix { - let mut v: Vec> = vec![Complex::zero(); n * n]; +pub fn cdiag(n: usize) -> ComplexMatrix { + let mut v: Vec = vec![Complex::zero(); n * n]; for i in 0..n { let idx = i * (n + 1); v[idx] = Complex::one(); } - complex_matrix(v, n, n, Shape::Row) + cmatrix(v, n, n, Shape::Row) } -/// Data structure for Complete Pivoting LU decomposition -/// -/// # Usage -/// ```rust -/// extern crate peroxide; -/// use peroxide::fuga::*; -/// use num_complex::Complex64; -/// use peroxide::complex::matrix::*; -/// use peroxide::complex::matrix::LinearAlgebra; -/// -/// let a = complex_matrix(vec![Complex64::new(1f64, 1f64), -/// Complex64::new(2f64, 2f64), -/// Complex64::new(3f64, 3f64), -/// Complex64::new(4f64, 4f64)], -/// 2, 2, Row -/// ); -/// let pqlu = a.lu(); -/// let (p, q, l, u) = pqlu.extract(); -/// // p, q are permutations -/// // l, u are matrices -/// println!("{}", l); // lower triangular -/// println!("{}", u); // upper triangular -/// ``` -#[derive(Debug, Clone)] -pub struct PQLU { - pub p: Vec, - pub q: Vec, - pub l: ComplexMatrix, - pub u: ComplexMatrix, -} - -impl PQLU { +impl PQLU { + /// Extract PQLU + /// + /// # Usage + /// ```rust + /// extern crate peroxide; + /// use peroxide::fuga::*; + /// + /// let a = cmatrix(vec![C64::new(1f64, 1f64), + /// C64::new(2f64, 2f64), + /// C64::new(3f64, 3f64), + /// C64::new(4f64, 4f64)], + /// 2, 2, Row + /// ); + /// let pqlu = a.lu(); + /// let (p, q, l, u) = pqlu.extract(); + /// // p, q are permutations + /// // l, u are matrices + /// println!("{}", l); // lower triangular + /// println!("{}", u); // upper triangular + /// ``` pub fn extract(&self) -> (Vec, Vec, ComplexMatrix, ComplexMatrix) { ( self.p.clone(), @@ -1815,7 +1738,7 @@ impl PQLU { ) } - pub fn det(&self) -> Complex { + pub fn det(&self) -> C64 { // sgn of perms let mut sgn_p = 1f64; let mut sgn_q = 1f64; @@ -1853,8 +1776,8 @@ impl PQLU { } /// MATLAB like eye - Identity matrix -pub fn eye(n: usize) -> ComplexMatrix { - let mut m = complex_matrix(vec![Complex::zero(); n * n], n, n, Shape::Row); +pub fn ceye(n: usize) -> ComplexMatrix { + let mut m = cmatrix(vec![Complex::zero(); n * n], n, n, Shape::Row); for i in 0..n { m[(i, i)] = Complex::one(); } @@ -1865,29 +1788,9 @@ pub fn eye(n: usize) -> ComplexMatrix { // Linear Algebra // ============================================================================= -#[derive(Debug, Copy, Clone)] -pub enum SolveKind { - LU, - WAZ, -} - -/// Linear algebra trait -pub trait LinearAlgebra { - fn back_subs(&self, b: &Vec>) -> Vec>; - fn forward_subs(&self, b: &Vec>) -> Vec>; - fn lu(&self) -> PQLU; - fn det(&self) -> Complex; - fn block(&self) -> (ComplexMatrix, ComplexMatrix, ComplexMatrix, ComplexMatrix); - fn inv(&self) -> ComplexMatrix; - fn solve(&self, b: &Vec>, sk: SolveKind) -> Vec>; - fn solve_mat(&self, m: &ComplexMatrix, sk: SolveKind) -> ComplexMatrix; - fn is_symmetric(&self) -> bool; - // ToDo: Add other fn of this trait from src/structure/matrix.rs -} - -impl LinearAlgebra for ComplexMatrix { +impl LinearAlgebra for ComplexMatrix { /// Backward Substitution for Upper Triangular - fn back_subs(&self, b: &Vec>) -> Vec> { + fn back_subs(&self, b: &[C64]) -> Vec { let n = self.col; let mut y = vec![Complex::zero(); n]; y[n - 1] = b[n - 1] / self[(n - 1, n - 1)]; @@ -1902,7 +1805,7 @@ impl LinearAlgebra for ComplexMatrix { } /// Forward substitution for Lower Triangular - fn forward_subs(&self, b: &Vec>) -> Vec> { + fn forward_subs(&self, b: &[C64]) -> Vec { let n = self.col; let mut y = vec![Complex::zero(); n]; y[0] = b[0] / self[(0, 0)]; @@ -1932,30 +1835,34 @@ impl LinearAlgebra for ComplexMatrix { /// ``` /// #[macro_use] /// use peroxide::fuga::*; - /// use num_complex::Complex64; /// use peroxide::complex::matrix::*; - /// use peroxide::complex::matrix::LinearAlgebra; /// /// fn main() { - /// let a = complex_matrix(vec![Complex64::new(1f64, 1f64), - /// Complex64::new(2f64, 2f64), - /// Complex64::new(3f64, 3f64), - /// Complex64::new(4f64, 4f64)], - /// 2, 2, Row + /// let a = cmatrix(vec![ + /// C64::new(1f64, 1f64), + /// C64::new(2f64, 2f64), + /// C64::new(3f64, 3f64), + /// C64::new(4f64, 4f64) + /// ], + /// 2, 2, Row /// ); /// - /// let l_exp = complex_matrix(vec![Complex64::new(1f64, 0f64), - /// Complex64::new(0f64, 0f64), - /// Complex64::new(0.5f64, -0.0f64), - /// Complex64::new(1f64, 0f64)], - /// 2, 2, Row + /// let l_exp = cmatrix(vec![ + /// C64::new(1f64, 0f64), + /// C64::new(0f64, 0f64), + /// C64::new(0.5f64, -0.0f64), + /// C64::new(1f64, 0f64) + /// ], + /// 2, 2, Row /// ); /// - /// let u_exp = complex_matrix(vec![Complex64::new(4f64, 4f64), - /// Complex64::new(3f64, 3f64), - /// Complex64::new(0f64, 0f64), - /// Complex64::new(-0.5f64, -0.5f64)], - /// 2, 2, Row + /// let u_exp = cmatrix(vec![ + /// C64::new(4f64, 4f64), + /// C64::new(3f64, 3f64), + /// C64::new(0f64, 0f64), + /// C64::new(-0.5f64, -0.5f64) + /// ], + /// 2, 2, Row /// ); /// let pqlu = a.lu(); /// let (p,q,l,u) = (pqlu.p, pqlu.q, pqlu.l, pqlu.u); @@ -1965,13 +1872,13 @@ impl LinearAlgebra for ComplexMatrix { /// assert_eq!(u, u_exp); /// } /// ``` - fn lu(&self) -> PQLU { + fn lu(&self) -> PQLU { assert_eq!(self.col, self.row); let n = self.row; let len: usize = n * n; - let mut l = eye(n); - let mut u = complex_matrix(vec![Complex::zero(); len], n, n, self.shape); + let mut l = ceye(n); + let mut u = cmatrix(vec![Complex::zero(); len], n, n, self.shape); let mut temp = self.clone(); let (p, q) = gecp(&mut temp); @@ -1997,27 +1904,40 @@ impl LinearAlgebra for ComplexMatrix { PQLU { p, q, l, u } } + fn waz(&self, d_form: Form) -> Option> { + unimplemented!() + } + + fn qr(&self) -> QR { unimplemented!() } + + fn svd(&self) -> SVD { unimplemented!() } + + #[cfg(feature = "O3")] + fn cholesky(&self, uplo: UPLO) -> ComplexMatrix { unimplemented!() } + + fn rref(&self) -> ComplexMatrix { unimplemented!() } + /// Determinant /// /// # Examples /// ``` /// #[macro_use] /// use peroxide::fuga::*; - /// use num_complex::Complex64; /// use peroxide::complex::matrix::*; - /// use peroxide::complex::matrix::LinearAlgebra; /// /// fn main() { - /// let a = complex_matrix(vec![Complex64::new(1f64, 1f64), - /// Complex64::new(2f64, 2f64), - /// Complex64::new(3f64, 3f64), - /// Complex64::new(4f64, 4f64)], - /// 2, 2, Row + /// let a = cmatrix(vec![ + /// C64::new(1f64, 1f64), + /// C64::new(2f64, 2f64), + /// C64::new(3f64, 3f64), + /// C64::new(4f64, 4f64) + /// ], + /// 2, 2, Row /// ); /// assert_eq!(a.det().norm(), 4f64); /// } /// ``` - fn det(&self) -> Complex { + fn det(&self) -> C64 { assert_eq!(self.row, self.col); self.lu().det() } @@ -2029,22 +1949,22 @@ impl LinearAlgebra for ComplexMatrix { /// #[macro_use] /// extern crate peroxide; /// use peroxide::fuga::*; - /// use num_complex::Complex64; /// use peroxide::complex::matrix::*; - /// use peroxide::complex::matrix::LinearAlgebra; /// /// fn main() { - /// let a = complex_matrix(vec![Complex64::new(1f64, 1f64), - /// Complex64::new(2f64, 2f64), - /// Complex64::new(3f64, 3f64), - /// Complex64::new(4f64, 4f64)], - /// 2, 2, Row + /// let a = cmatrix(vec![ + /// C64::new(1f64, 1f64), + /// C64::new(2f64, 2f64), + /// C64::new(3f64, 3f64), + /// C64::new(4f64, 4f64) + /// ], + /// 2, 2, Row /// ); /// let (m1, m2, m3, m4) = a.block(); - /// assert_eq!(m1, ml_complex_matrix("1.0+1.0i")); - /// assert_eq!(m2, ml_complex_matrix("2.0+2.0i")); - /// assert_eq!(m3, ml_complex_matrix("3.0+3.0i")); - /// assert_eq!(m4, ml_complex_matrix("4.0+4.0i")); + /// assert_eq!(m1, ml_cmatrix("1.0+1.0i")); + /// assert_eq!(m2, ml_cmatrix("2.0+2.0i")); + /// assert_eq!(m3, ml_cmatrix("3.0+3.0i")); + /// assert_eq!(m4, ml_cmatrix("4.0+4.0i")); /// } /// ``` fn block(&self) -> (Self, Self, Self, Self) { @@ -2055,10 +1975,10 @@ impl LinearAlgebra for ComplexMatrix { let r_l = r - l_r; let c_l = c - l_c; - let mut m1 = complex_matrix(vec![Complex::zero(); l_r * l_c], l_r, l_c, self.shape); - let mut m2 = complex_matrix(vec![Complex::zero(); l_r * c_l], l_r, c_l, self.shape); - let mut m3 = complex_matrix(vec![Complex::zero(); r_l * l_c], r_l, l_c, self.shape); - let mut m4 = complex_matrix(vec![Complex::zero(); r_l * c_l], r_l, c_l, self.shape); + let mut m1 = cmatrix(vec![Complex::zero(); l_r * l_c], l_r, l_c, self.shape); + let mut m2 = cmatrix(vec![Complex::zero(); l_r * c_l], l_r, c_l, self.shape); + let mut m3 = cmatrix(vec![Complex::zero(); r_l * l_c], r_l, l_c, self.shape); + let mut m4 = cmatrix(vec![Complex::zero(); r_l * c_l], r_l, c_l, self.shape); for idx_row in 0..r { for idx_col in 0..c { @@ -2094,24 +2014,26 @@ impl LinearAlgebra for ComplexMatrix { /// #[macro_use] /// extern crate peroxide; /// use peroxide::fuga::*; - /// use num_complex::Complex64; /// use peroxide::complex::matrix::*; - /// use peroxide::complex::matrix::LinearAlgebra; /// /// fn main() { /// // Non-singular - /// let a = complex_matrix(vec![Complex64::new(1f64, 1f64), - /// Complex64::new(2f64, 2f64), - /// Complex64::new(3f64, 3f64), - /// Complex64::new(4f64, 4f64)], - /// 2, 2, Row + /// let a = cmatrix(vec![ + /// C64::new(1f64, 1f64), + /// C64::new(2f64, 2f64), + /// C64::new(3f64, 3f64), + /// C64::new(4f64, 4f64) + /// ], + /// 2, 2, Row /// ); /// - /// let a_inv_exp = complex_matrix(vec![Complex64::new(-1.0f64, 1f64), - /// Complex64::new(0.5f64, -0.5f64), - /// Complex64::new(0.75f64, -0.75f64), - /// Complex64::new(-0.25f64, 0.25f64)], - /// 2, 2, Row + /// let a_inv_exp = cmatrix(vec![ + /// C64::new(-1.0f64, 1f64), + /// C64::new(0.5f64, -0.5f64), + /// C64::new(0.75f64, -0.75f64), + /// C64::new(-0.25f64, 0.25f64) + /// ], + /// 2, 2, Row /// ); /// assert_eq!(a.inv(), a_inv_exp); /// } @@ -2120,18 +2042,22 @@ impl LinearAlgebra for ComplexMatrix { self.lu().inv() } + fn pseudo_inv(&self) -> ComplexMatrix { + unimplemented!() + } + /// Solve with Vector /// /// # Solve options /// /// * LU: Gaussian elimination with Complete pivoting LU (GECP) /// * WAZ: Solve with WAZ decomposition - fn solve(&self, b: &Vec>, sk: SolveKind) -> Vec> { + fn solve(&self, b: &[C64], sk: SolveKind) -> Vec { match sk { SolveKind::LU => { let lu = self.lu(); let (p, q, l, u) = lu.extract(); - let mut v = b.clone(); + let mut v = b.to_vec(); v.swap_with_perm(&p.into_iter().enumerate().collect()); let z = l.forward_subs(&v); let mut y = u.back_subs(&z); @@ -2149,7 +2075,7 @@ impl LinearAlgebra for ComplexMatrix { SolveKind::LU => { let lu = self.lu(); let (p, q, l, u) = lu.extract(); - let mut x = complex_matrix( + let mut x = cmatrix( vec![Complex::zero(); self.col * m.col], self.col, m.col, @@ -2197,18 +2123,18 @@ impl LinearAlgebra for ComplexMatrix { } #[allow(non_snake_case)] -pub fn solve(A: &ComplexMatrix, b: &ComplexMatrix, sk: SolveKind) -> ComplexMatrix { +pub fn csolve(A: &ComplexMatrix, b: &ComplexMatrix, sk: SolveKind) -> ComplexMatrix { A.solve_mat(b, sk) } impl MutMatrix for ComplexMatrix { - type Scalar = Complex; + type Scalar = C64; - unsafe fn col_mut(&mut self, idx: usize) -> Vec<*mut Complex> { + unsafe fn col_mut(&mut self, idx: usize) -> Vec<*mut C64> { assert!(idx < self.col, "Index out of range"); match self.shape { Shape::Col => { - let mut v: Vec<*mut Complex> = vec![&mut Complex::zero(); self.row]; + let mut v: Vec<*mut C64> = vec![&mut Complex::zero(); self.row]; let start_idx = idx * self.row; let p = self.mut_ptr(); for (i, j) in (start_idx..start_idx + v.len()).enumerate() { @@ -2217,7 +2143,7 @@ impl MutMatrix for ComplexMatrix { v } Shape::Row => { - let mut v: Vec<*mut Complex> = vec![&mut Complex::zero(); self.row]; + let mut v: Vec<*mut C64> = vec![&mut Complex::zero(); self.row]; let p = self.mut_ptr(); for i in 0..v.len() { v[i] = p.add(idx + i * self.col); @@ -2227,11 +2153,11 @@ impl MutMatrix for ComplexMatrix { } } - unsafe fn row_mut(&mut self, idx: usize) -> Vec<*mut Complex> { + unsafe fn row_mut(&mut self, idx: usize) -> Vec<*mut C64> { assert!(idx < self.row, "Index out of range"); match self.shape { Shape::Row => { - let mut v: Vec<*mut Complex> = vec![&mut Complex::zero(); self.col]; + let mut v: Vec<*mut C64> = vec![&mut Complex::zero(); self.col]; let start_idx = idx * self.col; let p = self.mut_ptr(); for (i, j) in (start_idx..start_idx + v.len()).enumerate() { @@ -2240,7 +2166,7 @@ impl MutMatrix for ComplexMatrix { v } Shape::Col => { - let mut v: Vec<*mut Complex> = vec![&mut Complex::zero(); self.col]; + let mut v: Vec<*mut C64> = vec![&mut Complex::zero(); self.col]; let p = self.mut_ptr(); for i in 0..v.len() { v[i] = p.add(idx + i * self.row); @@ -2265,7 +2191,7 @@ impl MutMatrix for ComplexMatrix { } impl ExpLogOps for ComplexMatrix { - type Float = Complex; + type Float = C64; fn exp(&self) -> Self { self.fmap(|x| x.exp()) @@ -2285,7 +2211,7 @@ impl ExpLogOps for ComplexMatrix { } impl PowOps for ComplexMatrix { - type Float = Complex; + type Float = C64; fn powi(&self, n: i32) -> Self { self.fmap(|x| x.powi(n)) @@ -2308,8 +2234,8 @@ impl TrigOps for ComplexMatrix { fn sin_cos(&self) -> (Self, Self) { let (sin, cos) = self.data.iter().map(|x| (x.sin(), x.cos())).unzip(); ( - complex_matrix(sin, self.row, self.col, self.shape), - complex_matrix(cos, self.row, self.col, self.shape), + cmatrix(sin, self.row, self.col, self.shape), + cmatrix(cos, self.row, self.col, self.shape), ) } @@ -2370,22 +2296,21 @@ impl TrigOps for ComplexMatrix { /// # Examples /// ```rust /// use peroxide::fuga::*; -/// use num_complex::Complex64; /// use peroxide::complex::matrix::*; /// use peroxide::traits::fp::FPMatrix; /// /// fn main() { -/// let x1 = complex_matrix(vec![Complex64::new(1f64, 1f64)], 1, 1, Row); -/// let x2 = complex_matrix(vec![Complex64::new(2f64, 2f64)], 1, 1, Row); -/// let x3 = complex_matrix(vec![Complex64::new(3f64, 3f64)], 1, 1, Row); -/// let x4 = complex_matrix(vec![Complex64::new(4f64, 4f64)], 1, 1, Row); +/// let x1 = cmatrix(vec![C64::new(1f64, 1f64)], 1, 1, Row); +/// let x2 = cmatrix(vec![C64::new(2f64, 2f64)], 1, 1, Row); +/// let x3 = cmatrix(vec![C64::new(3f64, 3f64)], 1, 1, Row); +/// let x4 = cmatrix(vec![C64::new(4f64, 4f64)], 1, 1, Row); /// /// let y = complex_combine(x1, x2, x3, x4); /// -/// let y_exp = complex_matrix(vec![Complex64::new(1f64, 1f64), -/// Complex64::new(2f64, 2f64), -/// Complex64::new(3f64, 3f64), -/// Complex64::new(4f64, 4f64)], +/// let y_exp = cmatrix(vec![C64::new(1f64, 1f64), +/// C64::new(2f64, 2f64), +/// C64::new(3f64, 3f64), +/// C64::new(4f64, 4f64)], /// 2, 2, Row /// ); /// @@ -2406,7 +2331,7 @@ pub fn complex_combine( let r = l_r + r_l; let c = l_c + c_l; - let mut m = complex_matrix(vec![Complex::zero(); r * c], r, c, m1.shape); + let mut m = cmatrix(vec![Complex::zero(); r * c], r, c, m1.shape); for idx_row in 0..r { for idx_col in 0..c { @@ -2437,16 +2362,15 @@ pub fn complex_combine( /// #[macro_use] /// extern crate peroxide; /// use peroxide::fuga::*; -/// use num_complex::Complex64; /// use peroxide::complex::matrix::*; /// /// fn main() { -/// let a = ml_complex_matrix("2.0+2.0i 0.0+0.0i; +/// let a = ml_cmatrix("2.0+2.0i 0.0+0.0i; /// 2.0+2.0i 1.0+1.0i"); -/// let b = complex_matrix(vec![Complex64::new(2f64, 2f64), -/// Complex64::new(0f64, 0f64), -/// Complex64::new(-2f64, -2f64), -/// Complex64::new(1f64, 1f64)], +/// let b = cmatrix(vec![C64::new(2f64, 2f64), +/// C64::new(0f64, 0f64), +/// C64::new(-2f64, -2f64), +/// C64::new(1f64, 1f64)], /// 2, 2, Row /// ); /// assert_eq!(complex_inv_l(a), b); @@ -2481,16 +2405,15 @@ pub fn complex_inv_l(l: ComplexMatrix) -> ComplexMatrix { /// #[macro_use] /// extern crate peroxide; /// use peroxide::fuga::*; -/// use num_complex::Complex64; /// use peroxide::complex::matrix::*; /// /// fn main() { -/// let a = ml_complex_matrix("2.0+2.0i 2.0+2.0i; +/// let a = ml_cmatrix("2.0+2.0i 2.0+2.0i; /// 0.0+0.0i 1.0+1.0i"); -/// let b = complex_matrix(vec![Complex64::new(0.25f64, -0.25f64), -/// Complex64::new(-0.5f64, 0.5f64), -/// Complex64::new(0.0f64, 0.0f64), -/// Complex64::new(0.5f64, -0.5f64)], +/// let b = cmatrix(vec![C64::new(0.25f64, -0.25f64), +/// C64::new(-0.5f64, 0.5f64), +/// C64::new(0.0f64, 0.0f64), +/// C64::new(0.5f64, -0.5f64)], /// 2, 2, Row /// ); /// assert_eq!(complex_inv_u(a), b); @@ -2528,10 +2451,10 @@ pub fn complex_inv_u(u: ComplexMatrix) -> ComplexMatrix { } /// Matrix multiply back-ends -pub fn matmul(a: &ComplexMatrix, b: &ComplexMatrix) -> ComplexMatrix { +pub fn cmatmul(a: &ComplexMatrix, b: &ComplexMatrix) -> ComplexMatrix { assert_eq!(a.col, b.row); - let mut c = complex_matrix(vec![Complex::zero(); a.row * b.col], a.row, b.col, a.shape); - complex_gemm(Complex::one(), a, b, Complex::zero(), &mut c); + let mut c = cmatrix(vec![Complex::zero(); a.row * b.col], a.row, b.col, a.shape); + cgemm(Complex::one(), a, b, Complex::zero(), &mut c); c } @@ -2542,28 +2465,27 @@ pub fn matmul(a: &ComplexMatrix, b: &ComplexMatrix) -> ComplexMatrix { /// #[macro_use] /// extern crate peroxide; /// use peroxide::fuga::*; -/// use num_complex::Complex64; /// /// use peroxide::complex::matrix::*; /// /// fn main() { -/// let a = ml_complex_matrix("1.0+1.0i 2.0+2.0i; +/// let a = ml_cmatrix("1.0+1.0i 2.0+2.0i; /// 0.0+0.0i 1.0+1.0i"); -/// let b = ml_complex_matrix("1.0+1.0i 0.0+0.0i; +/// let b = ml_cmatrix("1.0+1.0i 0.0+0.0i; /// 2.0+2.0i 1.0+1.0i"); -/// let mut c1 = ml_complex_matrix("1.0+1.0i 1.0+1.0i; +/// let mut c1 = ml_cmatrix("1.0+1.0i 1.0+1.0i; /// 1.0+1.0i 1.0+1.0i"); -/// let mul_val = ml_complex_matrix("-10.0+10.0i -4.0+4.0i; +/// let mul_val = ml_cmatrix("-10.0+10.0i -4.0+4.0i; /// -4.0+4.0i -2.0+2.0i"); /// -/// complex_gemm(Complex64::new(1.0, 1.0), &a, &b, Complex64::new(0.0, 0.0), &mut c1); +/// cgemm(C64::new(1.0, 1.0), &a, &b, C64::new(0.0, 0.0), &mut c1); /// assert_eq!(c1, mul_val); /// } -pub fn complex_gemm( - alpha: Complex, +pub fn cgemm( + alpha: C64, a: &ComplexMatrix, b: &ComplexMatrix, - beta: Complex, + beta: C64, c: &mut ComplexMatrix, ) { let m = a.row; @@ -2606,12 +2528,12 @@ pub fn complex_gemm( } /// General Matrix-Vector multiplication -pub fn complex_gemv( - alpha: Complex, +pub fn cgemv( + alpha: C64, a: &ComplexMatrix, - b: &Vec>, - beta: Complex, - c: &mut Vec>, + b: &Vec, + beta: C64, + c: &mut Vec, ) { let m = a.row; let k = a.col; @@ -2648,11 +2570,11 @@ pub fn complex_gemv( /// General Vector-Matrix multiplication pub fn complex_gevm( - alpha: Complex, - a: &Vec>, + alpha: C64, + a: &Vec, b: &ComplexMatrix, - beta: Complex, - c: &mut Vec>, + beta: C64, + c: &mut Vec, ) { let m = 1usize; let k = a.len(); diff --git a/src/traits/matrix.rs b/src/traits/matrix.rs index 85b7b7fd..25d2979a 100644 --- a/src/traits/matrix.rs +++ b/src/traits/matrix.rs @@ -40,7 +40,7 @@ pub trait LinearAlgebra { #[cfg(feature = "O3")] fn cholesky(&self, uplo: UPLO) -> M; fn rref(&self) -> M; - fn det(&self) -> f64; + fn det(&self) -> M::Scalar; fn block(&self) -> (M, M, M, M); fn inv(&self) -> M; fn pseudo_inv(&self) -> M; From 206f46dbc36cdb2fb6b9bf8640055e60c35216c0 Mon Sep 17 00:00:00 2001 From: Axect Date: Fri, 25 Oct 2024 12:11:16 +0900 Subject: [PATCH 5/5] FIX: Fix old feature name --- tests/blas_test.rs | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/blas_test.rs b/tests/blas_test.rs index 5401eeac..59c2767e 100644 --- a/tests/blas_test.rs +++ b/tests/blas_test.rs @@ -3,7 +3,7 @@ extern crate peroxide; use peroxide::fuga::*; #[test] -#[cfg(feature = "native")] +#[cfg(feature = "O3")] fn daxpy_test() { let a = ml_matrix("1 2; 3 4"); let b = matrix(vec![1, 3, 2, 4], 2, 2, Col); @@ -14,7 +14,7 @@ fn daxpy_test() { } #[test] -#[cfg(feature = "native")] +#[cfg(feature = "O3")] fn dgemv_test() { let a = ml_matrix("1 2;3 4"); let b = c![1, 2]; @@ -24,7 +24,7 @@ fn dgemv_test() { } #[test] -#[cfg(feature = "native")] +#[cfg(feature = "O3")] fn dgemm_test() { // 2x2 let a = ml_matrix("1 2;3 4");