diff --git a/Cargo.toml b/Cargo.toml index 11204f7a..dd1dccfc 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -10,7 +10,14 @@ categories = ["science"] readme = "README.md" documentation = "https://axect.github.io/Peroxide_Doc" keywords = ["Numeric", "Science", "Dataframe", "Plot", "LinearAlgebra"] -exclude = ["example_data/", "src/bin/", "benches/", "example/", "test_data/", "peroxide-ad2"] +exclude = [ + "example_data/", + "src/bin/", + "benches/", + "example/", + "test_data/", + "peroxide-ad2", +] [badges] travis-ci = { repository = "axect/peroxide" } @@ -32,17 +39,26 @@ anyhow = "1.0" paste = "1.0" #num-complex = "0.3" netcdf = { version = "0.7", optional = true, default-features = false } -pyo3 = { version = "0.22", optional = true, features = ["auto-initialize", "gil-refs"] } +pyo3 = { version = "0.22", optional = true, features = [ + "auto-initialize", + "gil-refs", +] } blas = { version = "0.22", optional = true } lapack = { version = "0.19", optional = true } serde = { version = "1.0", features = ["derive"], optional = true } json = { version = "0.12", optional = true } -arrow2 = { version = "0.18", features = ["io_parquet", "io_parquet_compression"], optional = true } +arrow2 = { version = "0.18", features = [ + "io_parquet", + "io_parquet_compression", +], optional = true } num-complex = { version = "0.4", optional = true } -lambert_w = { version = "0.4.0", default-features = false, features = ["24bits", "50bits"] } +lambert_w = { version = "0.4.0", default-features = false, features = [ + "24bits", + "50bits", +] } [package.metadata.docs.rs] -rustdoc-args = [ "--html-in-header", "katex-header.html", "--cfg", "docsrs"] +rustdoc-args = ["--html-in-header", "katex-header.html", "--cfg", "docsrs"] [features] default = [] diff --git a/src/complex/integral.rs b/src/complex/integral.rs new file mode 100644 index 00000000..63fb2427 --- /dev/null +++ b/src/complex/integral.rs @@ -0,0 +1,53 @@ +use crate::complex::C64; +use crate::numerical::integral::{GKIntegrable, GLKIntegrable, NCIntegrable}; +use crate::structure::polynomial::{lagrange_polynomial, Calculus, Polynomial}; + +// Newton Cotes Quadrature for Complex Functions of one Real Variable +impl NCIntegrable for C64 { + type NodeY = (Vec, Vec); + type NCPolynomial = (Polynomial, Polynomial); + + fn compute_node_y(f: F, node_x: &[f64]) -> Self::NodeY + where + F: Fn(f64) -> Self, + { + node_x + .iter() + .map(|x| { + let z = f(*x); + (z.re, z.im) + }) + .unzip() + } + + fn compute_polynomial(node_x: &[f64], node_y: &Self::NodeY) -> Self::NCPolynomial { + ( + lagrange_polynomial(node_x.to_vec(), node_y.0.to_vec()), + lagrange_polynomial(node_x.to_vec(), node_y.1.to_vec()), + ) + } + + fn integrate_polynomial(p: &Self::NCPolynomial) -> Self::NCPolynomial { + (p.0.integral(), p.1.integral()) + } + + fn evaluate_polynomial(p: &Self::NCPolynomial, x: f64) -> Self { + p.0.eval(x) + C64::I * p.1.eval(x) + } +} + +// Gauss Lagrange and Kronrod Quadrature for Complex Functions of one Real Variable +impl GLKIntegrable for C64 { + const ZERO: Self = C64::ZERO; +} + +// Gauss Kronrod Quadrature for Complex Functions of one Real Variable +impl GKIntegrable for C64 { + fn is_finite(&self) -> bool { + C64::is_finite(*self) + } + + fn gk_norm(&self) -> f64 { + self.norm() + } +} diff --git a/src/complex/matrix.rs b/src/complex/matrix.rs new file mode 100644 index 00000000..2d6eb3ec --- /dev/null +++ b/src/complex/matrix.rs @@ -0,0 +1,2809 @@ +use std::{ + cmp::{max, min}, + fmt, + ops::{Add, Div, Index, IndexMut, Mul, Neg, Sub}, +}; + +use anyhow::{bail, Result}; +use matrixmultiply::CGemmOption; +use num_complex::Complex; +use peroxide_num::{ExpLogOps, PowOps, TrigOps}; +use rand_distr::num_traits::{One, Zero}; + +use crate::{ + structure::matrix::Shape, + traits::fp::{FPMatrix, FPVector}, + traits::general::Algorithm, + traits::math::{InnerProduct, LinearOp, MatrixProduct, Norm, Normed, Vector}, + traits::mutable::MutMatrix, + util::low_level::{copy_vec_ptr, swap_vec_ptr}, + util::non_macro::ConcatenateError, + util::useful::{nearly_eq, tab}, +}; + +/// R-like complex matrix structure +/// +/// # Examples +/// +/// ```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), +/// ], +/// row: 2, +/// col: 2, +/// shape: Row, +/// }; // [[1+1i,2+2i],[3+3i,4+4i]] +/// ``` +#[derive(Debug, Clone, Default)] +pub struct ComplexMatrix { + pub data: Vec>, + pub row: usize, + pub col: usize, + pub shape: Shape, +} + +// ============================================================================= +// Various complex matrix constructor +// ============================================================================= + +/// R-like complex matrix constructor +/// +/// # Examples +/// ```rust +/// #[macro_use] +/// extern crate peroxide; +/// use peroxide::fuga::*; +/// use num_complex::Complex64; +/// use peroxide::complex::matrix::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)], +/// 2, 2, Row +/// ); +/// a.col.print(); // Print matrix column +/// } +/// ``` +pub fn complex_matrix(v: Vec, r: usize, c: usize, shape: Shape) -> ComplexMatrix +where + T: Into>, +{ + ComplexMatrix { + data: v + .into_iter() + .map(|t| t.into()) + .collect::>>(), + row: r, + col: c, + shape, + } +} + +/// R-like complex matrix constructor (Explicit ver.) +pub fn r_complex_matrix(v: Vec, r: usize, c: usize, shape: Shape) -> ComplexMatrix +where + T: Into>, +{ + complex_matrix(v, r, c, shape) +} + +/// Python-like complex matrix constructor +/// +/// # Examples +/// ```rust +/// #[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 b = complex_matrix(vec![Complex64::new(1f64, 1f64), +/// Complex64::new(2f64, 2f64), +/// Complex64::new(3f64, 3f64), +/// Complex64::new(4f64, 4f64)], +/// 2, 2, Row +/// ); +/// assert_eq!(a, b); +/// } +/// ``` +pub fn py_complex_matrix(v: Vec>) -> ComplexMatrix +where + 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) +} + +/// Matlab-like matrix constructor +/// +/// Note that the entries to the `ml_complex_matrix` +/// needs to be in the `a+bi` format +/// without any spaces between the real and imaginary +/// parts of the Complex number. +/// +/// # Examples +/// ```rust +/// #[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 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)], +/// 2, 2, Row +/// ); +/// assert_eq!(a, b); +/// } +/// ``` +pub fn ml_complex_matrix(s: &str) -> ComplexMatrix { + let str_row = s.split(";").collect::>(); + let r = str_row.len(); + let str_data = str_row + .iter() + .map(|x| x.trim().split(" ").collect::>()) + .collect::>>(); + let c = str_data[0].len(); + let data = str_data + .iter() + .flat_map(|t| { + t.iter() + .map(|x| x.parse::>().unwrap()) + .collect::>>() + }) + .collect::>>(); + + complex_matrix(data, r, c, Shape::Row) +} + +/// Pretty Print +impl fmt::Display for ComplexMatrix { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> fmt::Result { + write!(f, "{}", self.spread()) + } +} + +/// PartialEq implements +impl PartialEq for ComplexMatrix { + fn eq(&self, other: &ComplexMatrix) -> bool { + if self.shape == other.shape { + self.data + .clone() + .into_iter() + .zip(other.data.clone()) + .all(|(x, y)| nearly_eq(x.re, y.re) && nearly_eq(x.im, y.im)) + && self.row == other.row + } else { + self.eq(&other.change_shape()) + } + } +} + +impl ComplexMatrix { + /// Raw pointer for `self.data` + pub fn ptr(&self) -> *const Complex { + &self.data[0] as *const Complex + } + + /// Raw mutable pointer for `self.data` + pub fn mut_ptr(&mut self) -> *mut Complex { + &mut self.data[0] as *mut Complex + } + + /// Slice of `self.data` + /// + /// # 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 b = a.as_slice(); + /// assert_eq!(b, &[Complex64::new(1f64, 1f64), + /// Complex64::new(2f64, 2f64), + /// Complex64::new(3f64, 3f64), + /// Complex64::new(4f64, 4f64)]); + /// ``` + pub fn as_slice(&self) -> &[Complex] { + &self.data[..] + } + + /// Mutable slice of `self.data` + /// + /// # 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 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)], + /// 2, 2, Row)); + /// ``` + pub fn as_mut_slice(&mut self) -> &mut [Complex] { + &mut self.data[..] + } + + /// Change Bindings + /// + /// `Row` -> `Col` or `Col` -> `Row` + /// + /// # 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 + /// ); + /// assert_eq!(a.shape, Row); + /// let b = a.change_shape(); + /// assert_eq!(b.shape, Col); + /// ``` + pub 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 ref_data = &self.data; + + match self.shape { + Shape::Row => { + for i in 0..l { + let s = (i * c) % l; + data[i] = ref_data[s]; + } + data[l] = ref_data[l]; + complex_matrix(data, r, c, Shape::Col) + } + Shape::Col => { + for i in 0..l { + let s = (i * r) % l; + data[i] = ref_data[s]; + } + data[l] = ref_data[l]; + complex_matrix(data, r, c, Shape::Row) + } + } + } + + /// Change Bindings Mutably + /// + /// `Row` -> `Col` or `Col` -> `Row` + /// + /// # 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 + /// ); + /// assert_eq!(a.shape, Row); + /// let b = a.change_shape_mut(); + /// assert_eq!(b.shape, Col); + /// ``` + pub fn change_shape_mut(&mut 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 ref_data = &self.data; + + match self.shape { + Shape::Row => { + for i in 0..l { + let s = (i * c) % l; + data[i] = ref_data[s]; + } + data[l] = ref_data[l]; + complex_matrix(data, r, c, Shape::Col) + } + Shape::Col => { + for i in 0..l { + let s = (i * r) % l; + data[i] = ref_data[s]; + } + data[l] = ref_data[l]; + complex_matrix(data, r, c, Shape::Row) + } + } + } + + /// Spread data(1D vector) to 2D formatted String + /// + /// # 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 + /// ); + /// println!("{}", a.spread()); // same as println!("{}", a); + /// // Result: + /// // c[0] c[1] + /// // r[0] 1+1i 3+3i + /// // r[1] 2+2i 4+4i + /// ``` + pub fn spread(&self) -> String { + assert_eq!(self.row * self.col, self.data.len()); + let r = self.row; + let c = self.col; + let mut key_row = 20usize; + let mut key_col = 20usize; + + if r > 100 || c > 100 || (r > 20 && c > 20) { + let part = if r <= 10 { + key_row = r; + key_col = 100; + self.take_col(100) + } else if c <= 10 { + key_row = 100; + key_col = c; + self.take_row(100) + } else { + self.take_row(20).take_col(20) + }; + return format!( + "Result is too large to print - {}x{}\n only print {}x{} parts:\n{}", + self.row.to_string(), + self.col.to_string(), + key_row.to_string(), + key_col.to_string(), + part.spread() + ); + } + + // Find maximum length of data + let sample = self.data.clone(); + let mut space: usize = sample + .into_iter() + .map( + |x| min(format!("{:.4}", x).len(), x.to_string().len()), // Choose minimum of approx vs normal + ) + .fold(0, |x, y| max(x, y)) + + 1; + + if space < 5 { + space = 5; + } + + let mut result = String::new(); + + result.push_str(&tab("", 5)); + for i in 0..c { + result.push_str(&tab(&format!("c[{}]", i), space)); // Header + } + result.push('\n'); + + for i in 0..r { + result.push_str(&tab(&format!("r[{}]", i), 5)); + for j in 0..c { + let st1 = format!("{:.4}", self[(i, j)]); // Round at fourth position + let st2 = self[(i, j)].to_string(); // Normal string + let mut st = st2.clone(); + + // Select more small thing + if st1.len() < st2.len() { + st = st1; + } + + result.push_str(&tab(&st, space)); + } + if i == (r - 1) { + break; + } + result.push('\n'); + } + + return result; + } + + /// Extract Column + /// + /// # Examples + /// ```rust + /// #[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)], + /// 2, 2, Row + /// ); + /// assert_eq!(a.col(0), vec![Complex64::new(1f64, 1f64), Complex64::new(3f64, 3f64)]); + /// } + /// ``` + pub fn col(&self, index: usize) -> Vec> { + assert!(index < self.col); + let mut container: Vec> = vec![Complex::zero(); self.row]; + for i in 0..self.row { + container[i] = self[(i, index)]; + } + container + } + + /// Extract Row + /// + /// # Examples + /// ```rust + /// #[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)], + /// 2, 2, Row + /// ); + /// assert_eq!(a.row(0), vec![Complex64::new(1f64, 1f64), Complex64::new(2f64, 2f64)]); + /// } + /// ``` + pub fn row(&self, index: usize) -> Vec> { + assert!(index < self.row); + let mut container: Vec> = vec![Complex::zero(); self.col]; + for i in 0..self.col { + container[i] = self[(index, i)]; + } + container + } + + /// Extract diagonal components + /// + /// # Examples + /// ```rust + /// #[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)], + /// 2, 2, Row + /// ); + /// assert_eq!(a.diag(), vec![Complex64::new(1f64, 1f64) ,Complex64::new(4f64, 4f64)]); + /// } + /// ``` + pub fn diag(&self) -> Vec> { + let mut container = vec![Complex::zero(); self.row]; + let r = self.row; + let c = self.col; + assert_eq!(r, c); + + let c2 = c + 1; + for i in 0..r { + container[i] = self.data[i * c2]; + } + container + } + + /// Transpose + /// + /// # 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.transpose(), a_t); + /// ``` + pub 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); + } + } + } + } + + /// 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)) -> ComplexMatrix + where + F: Fn(usize, usize) -> G + Copy, + G: Into>, + { + let row = size.0; + let col = size.1; + + let mut mat = complex_matrix(vec![Complex::zero(); row * col], row, col, Shape::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![Complex::zero(); self.col]; self.row]; + for i in 0..self.row { + result[i] = self.row(i); + } + result + } + + pub fn to_diag(&self) -> ComplexMatrix { + assert_eq!(self.row, self.col, "Should be square matrix"); + let mut result = complex_matrix( + vec![Complex::zero(); self.row * self.col], + self.row, + self.col, + Shape::Row, + ); + let diag = self.diag(); + for i in 0..self.row { + result[(i, i)] = diag[i]; + } + result + } + + /// Submatrix + /// + /// # Description + /// Return below elements of complex matrix to a new complex matrix + /// + /// $$ + /// \begin{pmatrix} + /// \\ddots & & & & \\\\ + /// & start & \\cdots & end.1 & \\\\ + /// & \\vdots & \\ddots & \\vdots & \\\\ + /// & end.0 & \\cdots & end & \\\\ + /// & & & & \\ddots + /// \end{pmatrix} + /// $$ + /// + /// # Examples + /// ```rust + /// #[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; + /// 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)], + /// 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 { + 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); + for i in 0..row { + for j in 0..col { + result[(i, j)] = self[(start.0 + i, start.1 + j)]; + } + } + result + } + + /// Substitute complex matrix to specific position + /// + /// # Description + /// Substitute below elements of complex 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::*; + /// 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; + /// 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)], + /// 2, 2, Row); + /// let c = ml_complex_matrix("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) { + 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)]; + } + } + } +} + +// ============================================================================= +// Mathematics for Matrix +// ============================================================================= +impl Vector for ComplexMatrix { + type Scalar = Complex; + + 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); + for i in 0..self.row { + for j in 0..self.col { + result[(i, j)] += other[(i, j)]; + } + } + result + } + + fn sub_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); + for i in 0..self.row { + for j in 0..self.col { + result[(i, j)] -= other[(i, j)]; + } + } + result + } + + fn mul_scalar(&self, other: Self::Scalar) -> Self { + let scalar = other; + self.fmap(|x| x * scalar) + } +} + +impl Normed for ComplexMatrix { + type UnsignedScalar = f64; + + fn norm(&self, kind: Norm) -> Self::UnsignedScalar { + match kind { + Norm::F => { + let mut s = Complex::zero(); + for i in 0..self.data.len() { + s += self.data[i].powi(2); + } + s.sqrt().re + } + Norm::Lpq(p, q) => { + let mut s = Complex::zero(); + for j in 0..self.col { + let mut s_row = Complex::zero(); + for i in 0..self.row { + s_row += self[(i, j)].powi(p as i32); + } + s += s_row.powf(q / p); + } + s.powf(1f64 / q).re + } + Norm::L1 => { + let mut m = Complex::zero(); + match self.shape { + Shape::Row => self.change_shape().norm(Norm::L1), + Shape::Col => { + for c in 0..self.col { + let s: Complex = self.col(c).iter().sum(); + if s.re > m.re { + m = s; + } + } + m.re + } + } + } + Norm::LInf => { + let mut m = Complex::zero(); + match self.shape { + Shape::Col => self.change_shape().norm(Norm::LInf), + Shape::Row => { + for r in 0..self.row { + let s: Complex = self.row(r).iter().sum(); + if s.re > m.re { + m = s; + } + } + m.re + } + } + } + Norm::L2 => { + unimplemented!() + } + Norm::Lp(_) => unimplemented!(), + } + } + + fn normalize(&self, _kind: Norm) -> Self + where + Self: Sized, + { + unimplemented!() + } +} + +/// Frobenius inner product +impl InnerProduct for ComplexMatrix { + fn dot(&self, rhs: &Self) -> Complex { + if self.shape == rhs.shape { + self.data.dot(&rhs.data) + } else { + self.data.dot(&rhs.change_shape().data) + } + } +} + +/// TODO: Transpose + +/// Matrix as Linear operator for Vector +#[allow(non_snake_case)] +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); + c + } +} + +/// R like cbind - concatenate two comlex matrix by column direction +pub fn complex_cbind(m1: ComplexMatrix, m2: ComplexMatrix) -> Result { + let mut temp = m1; + if temp.shape != Shape::Col { + temp = temp.change_shape(); + } + + let mut temp2 = m2; + if temp2.shape != Shape::Col { + temp2 = temp2.change_shape(); + } + + let mut v = temp.data; + let mut c = temp.col; + let r = temp.row; + + if r != temp2.row { + bail!(ConcatenateError::DifferentLength); + } + v.extend_from_slice(&temp2.data[..]); + c += temp2.col; + Ok(complex_matrix(v, r, c, Shape::Col)) +} + +/// R like rbind - concatenate two complex matrix by row direction +/// ``` +pub fn complex_rbind(m1: ComplexMatrix, m2: ComplexMatrix) -> Result { + let mut temp = m1; + if temp.shape != Shape::Row { + temp = temp.change_shape(); + } + + let mut temp2 = m2; + if temp2.shape != Shape::Row { + temp2 = temp2.change_shape(); + } + + let mut v = temp.data; + let c = temp.col; + let mut r = temp.row; + + if c != temp2.col { + bail!(ConcatenateError::DifferentLength); + } + v.extend_from_slice(&temp2.data[..]); + r += temp2.row; + Ok(complex_matrix(v, r, c, Shape::Row)) +} + +impl MatrixProduct for ComplexMatrix { + fn kronecker(&self, other: &Self) -> Self { + let r1 = self.row; + let c1 = self.col; + + let mut result = self[(0, 0)] * other; + + for j in 1..c1 { + let n = self[(0, j)] * other; + result = complex_cbind(result, n).unwrap(); + } + + for i in 1..r1 { + let mut m = self[(i, 0)] * other; + for j in 1..c1 { + let n = self[(i, j)] * other; + m = complex_cbind(m, n).unwrap(); + } + result = complex_rbind(result, m).unwrap(); + } + result + } + + fn hadamard(&self, other: &Self) -> Self { + assert_eq!(self.row, other.row); + assert_eq!(self.col, other.col); + + let r = self.row; + let c = self.col; + + let mut m = complex_matrix(vec![Complex::zero(); r * c], r, c, self.shape); + for i in 0..r { + for j in 0..c { + m[(i, j)] = self[(i, j)] * other[(i, j)] + } + } + m + } +} + +// ============================================================================= +// Common Properties of Matrix & 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> { + &self.data + } +} + +/// `Vec>` to `ComplexMatrix` +impl Into for Vec> { + fn into(self) -> ComplexMatrix { + let l = self.len(); + complex_matrix(self, l, 1, Shape::Col) + } +} + +/// `&Vec>` to `ComplexMatrix` +impl Into for &Vec> { + fn into(self) -> ComplexMatrix { + let l = self.len(); + complex_matrix(self.clone(), l, 1, Shape::Col) + } +} + +// ============================================================================= +// Standard Operation for Complex Matrix (ADD) +// ============================================================================= + +/// Element-wise addition of Complex Matrix +/// +/// # Caution +/// > You should remember ownership. +/// > If you use ComplexMatrix `a,b` then you can't use them after. +impl Add for ComplexMatrix { + type Output = Self; + + fn add(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); + for i in 0..self.row { + for j in 0..self.col { + result[(i, j)] += other[(i, j)]; + } + } + result + } +} + +impl<'a, 'b> Add<&'b ComplexMatrix> for &'a ComplexMatrix { + type Output = ComplexMatrix; + + fn add(self, rhs: &'b ComplexMatrix) -> Self::Output { + self.add_vec(rhs) + } +} + +/// Element-wise addition between Complex Matrix & Complex +/// +/// # 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; +/// 4.0+4.0i 5.0+5.0i"); +/// let a_exp = ml_complex_matrix("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); +/// } +/// ``` +impl Add for ComplexMatrix +where + T: Into> + Copy, +{ + type Output = Self; + fn add(self, other: T) -> Self { + self.fmap(|x| x + other.into()) + } +} + +/// Element-wise addition between &ComplexMatrix & Complex +impl<'a, T> Add for &'a ComplexMatrix +where + T: Into> + Copy, +{ + type Output = ComplexMatrix; + + fn add(self, other: T) -> Self::Output { + self.fmap(|x| x + other.into()) + } +} + +// Element-wise addition between Complex & ComplexMatrix +/// +/// # 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; +/// 4.0+4.0i 5.0+5.0i"); +/// let a_exp = ml_complex_matrix("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); +/// } +/// ``` +impl Add for Complex { + type Output = ComplexMatrix; + + fn add(self, other: ComplexMatrix) -> Self::Output { + other.add(self) + } +} + +/// Element-wise addition between Complex & &ComplexMatrix +impl<'a> Add<&'a ComplexMatrix> for Complex { + type Output = ComplexMatrix; + + fn add(self, other: &'a ComplexMatrix) -> Self::Output { + other.add(self) + } +} + +// ============================================================================= +// Standard Operation for Matrix (Neg) +// ============================================================================= +/// Negation of Complex Matrix +/// +/// # 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)], +/// 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)], +/// 2, 2, Row); +/// assert_eq!(-a, a_neg); +/// ``` +impl Neg for ComplexMatrix { + type Output = Self; + + fn neg(self) -> Self { + complex_matrix( + self.data + .into_iter() + .map(|x: Complex| -x) + .collect::>>(), + self.row, + self.col, + self.shape, + ) + } +} + +/// Negation of &'a Complex Matrix +impl<'a> Neg for &'a ComplexMatrix { + type Output = ComplexMatrix; + + fn neg(self) -> Self::Output { + complex_matrix( + self.data + .clone() + .into_iter() + .map(|x: Complex| -x) + .collect::>>(), + self.row, + self.col, + self.shape, + ) + } +} + +// ============================================================================= +// Standard Operation for Matrix (Sub) +// ============================================================================= +/// Subtraction between Complex Matrix +/// +/// # Examples +/// ```rust +/// #[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; +/// 40.0+40.0i 50.0+50.0i"); +/// let b = ml_complex_matrix("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; +/// 36.0+36.0i 45.0+45.0i"); +/// assert_eq!(a-b, diff); +/// } +/// ``` +impl Sub for ComplexMatrix { + type Output = Self; + + 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); + for i in 0..self.row { + for j in 0..self.col { + result[(i, j)] -= other[(i, j)]; + } + } + result + } +} + +impl<'a, 'b> Sub<&'b ComplexMatrix> for &'a ComplexMatrix { + type Output = ComplexMatrix; + + fn sub(self, rhs: &'b ComplexMatrix) -> Self::Output { + self.sub_vec(rhs) + } +} + +/// Subtraction between Complex Matrix & Complex +impl Sub for ComplexMatrix +where + T: Into> + Copy, +{ + type Output = Self; + + fn sub(self, other: T) -> Self::Output { + self.fmap(|x| x - other.into()) + } +} + +/// Subtraction between &Complex Matrix & Complex +impl<'a, T> Sub for &'a ComplexMatrix +where + T: Into> + Copy, +{ + type Output = ComplexMatrix; + + fn sub(self, other: T) -> Self::Output { + self.fmap(|x| x - other.into()) + } +} + +/// Subtraction Complex Matrix with Complex +/// +/// # 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; +/// 4.0+4.0i 5.0+5.0i"); +/// let a_exp = ml_complex_matrix("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); +/// } +/// ``` +impl Sub for Complex { + type Output = ComplexMatrix; + + fn sub(self, other: ComplexMatrix) -> Self::Output { + -other.sub(self) + } +} + +impl<'a> Sub<&'a ComplexMatrix> for f64 { + type Output = ComplexMatrix; + + fn sub(self, other: &'a ComplexMatrix) -> Self::Output { + -other.sub(self) + } +} + +// ============================================================================= +// Multiplication for Complex Matrix +// ============================================================================= +/// Element-wise multiplication between Complex Matrix vs Complex +impl Mul> for ComplexMatrix { + type Output = Self; + + fn mul(self, other: Complex) -> Self::Output { + self.fmap(|x| x * other) + } +} + +impl Mul for Complex { + type Output = ComplexMatrix; + + fn mul(self, other: ComplexMatrix) -> Self::Output { + other.mul(self) + } +} + +impl<'a> Mul<&'a ComplexMatrix> for Complex { + type Output = ComplexMatrix; + + fn mul(self, other: &'a ComplexMatrix) -> Self::Output { + other.mul_scalar(self) + } +} + +/// Matrix Multiplication +/// +/// # 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; +/// 4.0+4.0i 5.0+5.0i"); +/// let mut b = ml_complex_matrix("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; +/// 0.0+66.0i 0.0+66.0i"); +/// assert_eq!(a * b, prod); +/// } +/// ``` +impl Mul for ComplexMatrix { + type Output = Self; + + fn mul(self, other: Self) -> Self::Output { + matmul(&self, &other) + } +} + +impl<'a, 'b> Mul<&'b ComplexMatrix> for &'a ComplexMatrix { + type Output = ComplexMatrix; + + fn mul(self, other: &'b ComplexMatrix) -> Self::Output { + matmul(self, other) + } +} + +#[allow(non_snake_case)] +impl Mul>> for ComplexMatrix { + type Output = Vec>; + + 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>; + + fn mul(self, other: &'b Vec>) -> Self::Output { + self.apply(other) + } +} + +/// 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); + let mut c = vec![Complex::zero(); other.col]; + complex_gevm(Complex::one(), &self, &other, Complex::zero(), &mut c); + c + } +} + +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); + let mut c = vec![Complex::zero(); other.col]; + complex_gevm(Complex::one(), self, other, Complex::zero(), &mut c); + c + } +} + +// ============================================================================= +// Standard Operation for Matrix (DIV) +// ============================================================================= +/// Element-wise division between Complex Matrix vs Complex +impl Div> for ComplexMatrix { + type Output = Self; + + fn div(self, other: Complex) -> Self::Output { + self.fmap(|x| x / other) + } +} + +impl<'a> Div> for &'a ComplexMatrix { + type Output = ComplexMatrix; + + fn div(self, other: Complex) -> Self::Output { + self.fmap(|x| x / other) + } +} + +/// Index for Complex Matrix +/// +/// `(usize, usize) -> Complex` +/// +/// # 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)], +/// 2, 2, Row +/// ); +/// assert_eq!(a[(0,1)], Complex64::new(2f64, 2f64)); +/// ``` +impl Index<(usize, usize)> for ComplexMatrix { + type Output = Complex; + + fn index(&self, pair: (usize, usize)) -> &Complex { + let p = self.ptr(); + let i = pair.0; + let j = pair.1; + assert!(i < self.row && j < self.col, "Index out of range"); + match self.shape { + Shape::Row => unsafe { &*p.add(i * self.col + j) }, + Shape::Col => unsafe { &*p.add(i + j * self.row) }, + } + } +} + +/// IndexMut for Complex Matrix (Assign) +/// +/// `(usize, usize) -> Complex` +/// +/// # 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)], +/// 2, 2, Row +/// ); +/// assert_eq!(a[(0,1)], Complex64::new(2f64, 2f64)); +/// ``` +impl IndexMut<(usize, usize)> for ComplexMatrix { + fn index_mut(&mut self, pair: (usize, usize)) -> &mut Complex { + let i = pair.0; + let j = pair.1; + let r = self.row; + let c = self.col; + assert!(i < self.row && j < self.col, "Index out of range"); + let p = self.mut_ptr(); + match self.shape { + Shape::Row => { + let idx = i * c + j; + unsafe { &mut *p.add(idx) } + } + Shape::Col => { + let idx = i + j * r; + unsafe { &mut *p.add(idx) } + } + } + } +} + +// ============================================================================= +// Functional Programming Tools (Hand-written) +// ============================================================================= + +impl FPMatrix for ComplexMatrix { + type Scalar = Complex; + + fn take_row(&self, n: usize) -> Self { + if n >= self.row { + return self.clone(); + } + match self.shape { + Shape::Row => { + let new_data = self + .data + .clone() + .into_iter() + .take(n * self.col) + .collect::>>(); + complex_matrix(new_data, n, self.col, Shape::Row) + } + Shape::Col => { + 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) + } + } + } + + fn take_col(&self, n: usize) -> Self { + if n >= self.col { + return self.clone(); + } + match self.shape { + Shape::Col => { + let new_data = self + .data + .clone() + .into_iter() + .take(n * self.row) + .collect::>>(); + complex_matrix(new_data, self.row, n, Shape::Col) + } + Shape::Row => { + 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) + } + } + } + + 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(); + for i in n..self.row { + temp_data.extend(self.row(i)); + } + complex_matrix(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(); + for i in n..self.col { + temp_data.extend(self.col(i)); + } + complex_matrix(temp_data, self.row, self.col - n, Shape::Col) + } + + fn fmap(&self, f: F) -> Self + where + F: Fn(Complex) -> Complex, + { + let result = self + .data + .iter() + .map(|x| f(*x)) + .collect::>>(); + complex_matrix(result, self.row, self.col, self.shape) + } + + /// Column map + /// + /// # 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)], + /// 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)], + /// 2, 2, Col + /// ); + /// + /// assert_eq!(y, y_col_map); + /// } + /// ``` + fn col_map(&self, f: F) -> ComplexMatrix + where + F: Fn(Vec>) -> Vec>, + { + let mut result = complex_matrix( + vec![Complex::zero(); self.row * self.col], + self.row, + self.col, + Shape::Col, + ); + + for i in 0..self.col { + result.subs_col(i, &f(self.col(i))); + } + + result + } + + /// Row map + /// + /// # 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)], + /// 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)], + /// 2, 2, Row + /// ); + /// + /// assert_eq!(y, y_row_map); + /// } + /// ``` + fn row_map(&self, f: F) -> ComplexMatrix + where + F: Fn(Vec>) -> Vec>, + { + let mut result = complex_matrix( + vec![Complex::zero(); self.row * self.col], + self.row, + self.col, + Shape::Row, + ); + + for i in 0..self.row { + result.subs_row(i, &f(self.row(i))); + } + + result + } + + fn col_mut_map(&mut self, f: F) + where + F: Fn(Vec>) -> Vec>, + { + for i in 0..self.col { + unsafe { + let mut p = self.col_mut(i); + let fv = f(self.col(i)); + for j in 0..p.len() { + *p[j] = fv[j]; + } + } + } + } + + fn row_mut_map(&mut self, f: F) + where + F: Fn(Vec>) -> Vec>, + { + for i in 0..self.col { + unsafe { + let mut p = self.row_mut(i); + let fv = f(self.row(i)); + for j in 0..p.len() { + *p[j] = fv[j]; + } + } + } + } + + fn reduce(&self, init: T, f: F) -> Complex + where + F: Fn(Complex, Complex) -> Complex, + 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, + { + assert_eq!(self.data.len(), other.data.len()); + let mut a = other.clone(); + if self.shape != other.shape { + a = a.change_shape(); + } + let result = self + .data + .iter() + .zip(a.data.iter()) + .map(|(x, y)| f(*x, *y)) + .collect::>>(); + complex_matrix(result, self.row, self.col, self.shape) + } + + fn col_reduce(&self, f: F) -> Vec> + where + F: Fn(Vec>) -> Complex, + { + let mut v = vec![Complex::zero(); self.col]; + for i in 0..self.col { + v[i] = f(self.col(i)); + } + v + } + + fn row_reduce(&self, f: F) -> Vec> + where + F: Fn(Vec>) -> Complex, + { + let mut v = vec![Complex::zero(); self.row]; + for i in 0..self.row { + v[i] = f(self.row(i)); + } + v + } +} + +pub fn diag(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) +} + +/// 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 { + pub fn extract(&self) -> (Vec, Vec, ComplexMatrix, ComplexMatrix) { + ( + self.p.clone(), + self.q.clone(), + self.l.clone(), + self.u.clone(), + ) + } + + pub fn det(&self) -> Complex { + // sgn of perms + let mut sgn_p = 1f64; + let mut sgn_q = 1f64; + for (i, &j) in self.p.iter().enumerate() { + if i != j { + sgn_p *= -1f64; + } + } + for (i, &j) in self.q.iter().enumerate() { + if i != j { + sgn_q *= -1f64; + } + } + + self.u.diag().reduce(Complex::one(), |x, y| x * y) * sgn_p * sgn_q + } + + pub fn inv(&self) -> ComplexMatrix { + let (p, q, l, u) = self.extract(); + let mut m = complex_inv_u(u) * complex_inv_l(l); + // Q = Q1 Q2 Q3 .. + for (idx1, idx2) in q.into_iter().enumerate().rev() { + unsafe { + m.swap(idx1, idx2, Shape::Row); + } + } + // P = Pn-1 .. P3 P2 P1 + for (idx1, idx2) in p.into_iter().enumerate().rev() { + unsafe { + m.swap(idx1, idx2, Shape::Col); + } + } + m + } +} + +/// 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); + for i in 0..n { + m[(i, i)] = Complex::one(); + } + m +} + +// ============================================================================= +// 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 { + /// Backward Substitution for Upper Triangular + fn back_subs(&self, b: &Vec>) -> Vec> { + let n = self.col; + let mut y = vec![Complex::zero(); n]; + y[n - 1] = b[n - 1] / self[(n - 1, n - 1)]; + for i in (0..n - 1).rev() { + let mut s = Complex::zero(); + for j in i + 1..n { + s += self[(i, j)] * y[j]; + } + y[i] = 1f64 / self[(i, i)] * (b[i] - s); + } + y + } + + /// Forward substitution for Lower Triangular + fn forward_subs(&self, b: &Vec>) -> Vec> { + let n = self.col; + let mut y = vec![Complex::zero(); n]; + y[0] = b[0] / self[(0, 0)]; + for i in 1..n { + let mut s = Complex::zero(); + for j in 0..i { + s += self[(i, j)] * y[j]; + } + y[i] = 1f64 / self[(i, i)] * (b[i] - s); + } + y + } + + /// LU Decomposition Implements (Complete Pivot) + /// + /// # Description + /// It use complete pivoting LU decomposition. + /// You can get two permutations, and LU matrices. + /// + /// # Caution + /// It returns `Option` - You should unwrap to obtain real value. + /// `PQLU` has four field - `p`, `q`, `l`, `u`. + /// `p`, `q` are permutations. + /// `l`, `u` are matrices. + /// + /// # 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 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 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 pqlu = a.lu(); + /// 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, l_exp); + /// assert_eq!(u, u_exp); + /// } + /// ``` + 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 temp = self.clone(); + let (p, q) = gecp(&mut temp); + for i in 0..n { + for j in 0..i { + // Inverse multiplier + l[(i, j)] = -temp[(i, j)]; + } + for j in i..n { + u[(i, j)] = temp[(i, j)]; + } + } + // Pivoting L + for i in 0..n - 1 { + unsafe { + let l_i = l.col_mut(i); + for j in i + 1..l.col - 1 { + let dst = p[j]; + std::ptr::swap(l_i[j], l_i[dst]); + } + } + } + PQLU { p, q, l, u } + } + + /// 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 + /// ); + /// assert_eq!(a.det().norm(), 4f64); + /// } + /// ``` + fn det(&self) -> Complex { + assert_eq!(self.row, self.col); + self.lu().det() + } + + /// Block Partition + /// + /// # Examples + /// ```rust + /// #[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 (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")); + /// } + /// ``` + fn block(&self) -> (Self, Self, Self, Self) { + let r = self.row; + let c = self.col; + let l_r = self.row / 2; + let l_c = self.col / 2; + 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); + + for idx_row in 0..r { + for idx_col in 0..c { + match (idx_row, idx_col) { + (i, j) if (i < l_r) && (j < l_c) => { + m1[(i, j)] = self[(i, j)]; + } + (i, j) if (i < l_r) && (j >= l_c) => { + m2[(i, j - l_c)] = self[(i, j)]; + } + (i, j) if (i >= l_r) && (j < l_c) => { + m3[(i - l_r, j)] = self[(i, j)]; + } + (i, j) if (i >= l_r) && (j >= l_c) => { + m4[(i - l_r, j - l_c)] = self[(i, j)]; + } + _ => (), + } + } + } + (m1, m2, m3, m4) + } + + /// Inverse of Matrix + /// + /// # Caution + /// + /// `inv` function returns `Option` + /// Thus, you should use pattern matching or `unwrap` to obtain inverse. + /// + /// # Examples + /// ``` + /// #[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_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 + /// ); + /// assert_eq!(a.inv(), a_inv_exp); + /// } + /// ``` + fn inv(&self) -> Self { + self.lu().inv() + } + + /// 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> { + match sk { + SolveKind::LU => { + let lu = self.lu(); + let (p, q, l, u) = lu.extract(); + let mut v = b.clone(); + v.swap_with_perm(&p.into_iter().enumerate().collect()); + let z = l.forward_subs(&v); + let mut y = u.back_subs(&z); + y.swap_with_perm(&q.into_iter().enumerate().rev().collect()); + y + } + SolveKind::WAZ => { + unimplemented!() + } + } + } + + fn solve_mat(&self, m: &ComplexMatrix, sk: SolveKind) -> ComplexMatrix { + match sk { + SolveKind::LU => { + let lu = self.lu(); + let (p, q, l, u) = lu.extract(); + let mut x = complex_matrix( + vec![Complex::zero(); self.col * m.col], + self.col, + m.col, + Shape::Col, + ); + for i in 0..m.col { + let mut v = m.col(i).clone(); + for (r, &s) in p.iter().enumerate() { + v.swap(r, s); + } + let z = l.forward_subs(&v); + let mut y = u.back_subs(&z); + for (r, &s) in q.iter().enumerate() { + y.swap(r, s); + } + unsafe { + let mut c = x.col_mut(i); + copy_vec_ptr(&mut c, &y); + } + } + x + } + SolveKind::WAZ => { + unimplemented!() + } + } + } + + fn is_symmetric(&self) -> bool { + if self.row != self.col { + return false; + } + + for i in 0..self.row { + for j in i..self.col { + if (!nearly_eq(self[(i, j)].re, self[(j, i)].re)) + && (!nearly_eq(self[(i, j)].im, self[(j, i)].im)) + { + return false; + } + } + } + true + } +} + +#[allow(non_snake_case)] +pub fn solve(A: &ComplexMatrix, b: &ComplexMatrix, sk: SolveKind) -> ComplexMatrix { + A.solve_mat(b, sk) +} + +impl MutMatrix for ComplexMatrix { + type Scalar = Complex; + + unsafe fn col_mut(&mut self, idx: usize) -> Vec<*mut Complex> { + 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 start_idx = idx * self.row; + let p = self.mut_ptr(); + for (i, j) in (start_idx..start_idx + v.len()).enumerate() { + v[i] = p.add(j); + } + v + } + Shape::Row => { + let mut v: Vec<*mut Complex> = 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); + } + v + } + } + } + + unsafe fn row_mut(&mut self, idx: usize) -> Vec<*mut Complex> { + 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 start_idx = idx * self.col; + let p = self.mut_ptr(); + for (i, j) in (start_idx..start_idx + v.len()).enumerate() { + v[i] = p.add(j); + } + v + } + Shape::Col => { + let mut v: Vec<*mut Complex> = 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); + } + v + } + } + } + + unsafe fn swap(&mut self, idx1: usize, idx2: usize, shape: Shape) { + match shape { + Shape::Col => swap_vec_ptr(&mut self.col_mut(idx1), &mut self.col_mut(idx2)), + Shape::Row => swap_vec_ptr(&mut self.row_mut(idx1), &mut self.row_mut(idx2)), + } + } + + unsafe fn swap_with_perm(&mut self, p: &Vec<(usize, usize)>, shape: Shape) { + for (i, j) in p.iter() { + self.swap(*i, *j, shape) + } + } +} + +impl ExpLogOps for ComplexMatrix { + type Float = Complex; + + fn exp(&self) -> Self { + self.fmap(|x| x.exp()) + } + fn ln(&self) -> Self { + self.fmap(|x| x.ln()) + } + fn log(&self, base: Self::Float) -> Self { + self.fmap(|x| x.ln() / base.ln()) // Using `Log: change of base` formula + } + fn log2(&self) -> Self { + self.fmap(|x| x.ln() / 2.0.ln()) // Using `Log: change of base` formula + } + fn log10(&self) -> Self { + self.fmap(|x| x.ln() / 10.0.ln()) // Using `Log: change of base` formula + } +} + +impl PowOps for ComplexMatrix { + type Float = Complex; + + fn powi(&self, n: i32) -> Self { + self.fmap(|x| x.powi(n)) + } + + fn powf(&self, f: Self::Float) -> Self { + self.fmap(|x| x.powc(f)) + } + + fn pow(&self, _f: Self) -> Self { + unimplemented!() + } + + fn sqrt(&self) -> Self { + self.fmap(|x| x.sqrt()) + } +} + +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), + ) + } + + fn sin(&self) -> Self { + self.fmap(|x| x.sin()) + } + + fn cos(&self) -> Self { + self.fmap(|x| x.cos()) + } + + fn tan(&self) -> Self { + self.fmap(|x| x.tan()) + } + + fn sinh(&self) -> Self { + self.fmap(|x| x.sinh()) + } + + fn cosh(&self) -> Self { + self.fmap(|x| x.cosh()) + } + + fn tanh(&self) -> Self { + self.fmap(|x| x.tanh()) + } + + fn asin(&self) -> Self { + self.fmap(|x| x.asin()) + } + + fn acos(&self) -> Self { + self.fmap(|x| x.acos()) + } + + fn atan(&self) -> Self { + self.fmap(|x| x.atan()) + } + + fn asinh(&self) -> Self { + self.fmap(|x| x.asinh()) + } + + fn acosh(&self) -> Self { + self.fmap(|x| x.acosh()) + } + + fn atanh(&self) -> Self { + self.fmap(|x| x.atanh()) + } +} + +// ============================================================================= +// Back-end Utils +// ============================================================================= +/// Combine separated Complex Matrix to one Complex Matrix +/// +/// # 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 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)], +/// 2, 2, Row +/// ); +/// +/// assert_eq!(y, y_exp); +/// } +/// ``` +pub fn complex_combine( + m1: ComplexMatrix, + m2: ComplexMatrix, + m3: ComplexMatrix, + m4: ComplexMatrix, +) -> ComplexMatrix { + let l_r = m1.row; + let l_c = m1.col; + let c_l = m2.col; + let r_l = m3.row; + + 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); + + for idx_row in 0..r { + for idx_col in 0..c { + match (idx_row, idx_col) { + (i, j) if (i < l_r) && (j < l_c) => { + m[(i, j)] = m1[(i, j)]; + } + (i, j) if (i < l_r) && (j >= l_c) => { + m[(i, j)] = m2[(i, j - l_c)]; + } + (i, j) if (i >= l_r) && (j < l_c) => { + m[(i, j)] = m3[(i - l_r, j)]; + } + (i, j) if (i >= l_r) && (j >= l_c) => { + m[(i, j)] = m4[(i - l_r, j - l_c)]; + } + _ => (), + } + } + } + m +} + +/// Inverse of Lower matrix +/// +/// # Examples +/// ```rust +/// #[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; +/// 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)], +/// 2, 2, Row +/// ); +/// assert_eq!(complex_inv_l(a), b); +/// } +/// ``` +pub fn complex_inv_l(l: ComplexMatrix) -> ComplexMatrix { + let mut m = l.clone(); + + match l.row { + 1 => l, + 2 => { + m[(1, 0)] = -m[(1, 0)]; + m + } + _ => { + let (l1, l2, l3, l4) = l.block(); + + let m1 = complex_inv_l(l1); + let m2 = l2; + let m4 = complex_inv_l(l4); + let m3 = -(&(&m4 * &l3) * &m1); + + complex_combine(m1, m2, m3, m4) + } + } +} + +/// Inverse of upper triangular matrix +/// +/// # Examples +/// ```rust +/// #[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; +/// 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)], +/// 2, 2, Row +/// ); +/// assert_eq!(complex_inv_u(a), b); +/// } +/// ``` +pub fn complex_inv_u(u: ComplexMatrix) -> ComplexMatrix { + let mut w = u.clone(); + + match u.row { + 1 => { + w[(0, 0)] = 1f64 / w[(0, 0)]; + w + } + 2 => { + let a = w[(0, 0)]; + let b = w[(0, 1)]; + let c = w[(1, 1)]; + let d = a * c; + + w[(0, 0)] = 1f64 / a; + w[(0, 1)] = -b / d; + w[(1, 1)] = 1f64 / c; + w + } + _ => { + let (u1, u2, u3, u4) = u.block(); + let m1 = complex_inv_u(u1); + let m3 = u3; + let m4 = complex_inv_u(u4); + let m2 = -(m1.clone() * u2 * m4.clone()); + + complex_combine(m1, m2, m3, m4) + } + } +} + +/// Matrix multiply back-ends +pub fn matmul(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); + c +} + +/// GEMM wrapper for Matrixmultiply +/// +/// # Examples +/// ```rust +/// #[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; +/// 0.0+0.0i 1.0+1.0i"); +/// let b = ml_complex_matrix("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; +/// 1.0+1.0i 1.0+1.0i"); +/// let mul_val = ml_complex_matrix("-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); +/// assert_eq!(c1, mul_val); +/// } +pub fn complex_gemm( + alpha: Complex, + a: &ComplexMatrix, + b: &ComplexMatrix, + beta: Complex, + c: &mut ComplexMatrix, +) { + let m = a.row; + let k = a.col; + let n = b.col; + let (rsa, csa) = match a.shape { + Shape::Row => (a.col as isize, 1isize), + Shape::Col => (1isize, a.row as isize), + }; + let (rsb, csb) = match b.shape { + Shape::Row => (b.col as isize, 1isize), + Shape::Col => (1isize, b.row as isize), + }; + let (rsc, csc) = match c.shape { + Shape::Row => (c.col as isize, 1isize), + Shape::Col => (1isize, c.row as isize), + }; + + unsafe { + matrixmultiply::zgemm( + // Requires crate feature "cgemm" + CGemmOption::Standard, + CGemmOption::Standard, + m, + k, + n, + [alpha.re, alpha.im], + a.ptr() as *const _, + rsa, + csa, + b.ptr() as *const _, + rsb, + csb, + [beta.re, beta.im], + c.mut_ptr() as *mut _, + rsc, + csc, + ) + } +} + +/// General Matrix-Vector multiplication +pub fn complex_gemv( + alpha: Complex, + a: &ComplexMatrix, + b: &Vec>, + beta: Complex, + c: &mut Vec>, +) { + let m = a.row; + let k = a.col; + let n = 1usize; + let (rsa, csa) = match a.shape { + Shape::Row => (a.col as isize, 1isize), + Shape::Col => (1isize, a.row as isize), + }; + let (rsb, csb) = (1isize, 1isize); + let (rsc, csc) = (1isize, 1isize); + + unsafe { + matrixmultiply::zgemm( + // Requires crate feature "cgemm" + CGemmOption::Standard, + CGemmOption::Standard, + m, + k, + n, + [alpha.re, alpha.im], + a.ptr() as *const _, + rsa, + csa, + b.as_ptr() as *const _, + rsb, + csb, + [beta.re, beta.im], + c.as_mut_ptr() as *mut _, + rsc, + csc, + ) + } +} + +/// General Vector-Matrix multiplication +pub fn complex_gevm( + alpha: Complex, + a: &Vec>, + b: &ComplexMatrix, + beta: Complex, + c: &mut Vec>, +) { + let m = 1usize; + let k = a.len(); + let n = b.col; + let (rsa, csa) = (1isize, 1isize); + let (rsb, csb) = match b.shape { + Shape::Row => (b.col as isize, 1isize), + Shape::Col => (1isize, b.row as isize), + }; + let (rsc, csc) = (1isize, 1isize); + + unsafe { + matrixmultiply::zgemm( + // Requires crate feature "cgemm" + CGemmOption::Standard, + CGemmOption::Standard, + m, + k, + n, + [alpha.re, alpha.im], + a.as_ptr() as *const _, + rsa, + csa, + b.ptr() as *const _, + rsb, + csb, + [beta.re, beta.im], + c.as_mut_ptr() as *mut _, + rsc, + csc, + ) + } +} + +/// LU via Gaussian Elimination with Partial Pivoting +#[allow(dead_code)] +fn gepp(m: &mut ComplexMatrix) -> Vec { + let mut r = vec![0usize; m.col - 1]; + for k in 0..(m.col - 1) { + // Find the pivot row + let r_k = m + .col(k) + .into_iter() + .skip(k) + .enumerate() + .max_by(|x1, x2| x1.1.norm().partial_cmp(&x2.1.norm()).unwrap()) + .unwrap() + .0 + + k; + r[k] = r_k; + + // Interchange the rows r_k and k + for j in k..m.col { + unsafe { + std::ptr::swap(&mut m[(k, j)], &mut m[(r_k, j)]); + println!("Swap! k:{}, r_k:{}", k, r_k); + } + } + // Form the multipliers + for i in k + 1..m.col { + m[(i, k)] = -m[(i, k)] / m[(k, k)]; + } + // Update the entries + for i in k + 1..m.col { + for j in k + 1..m.col { + let local_m = m[(i, k)] * m[(k, j)]; + m[(i, j)] += local_m; + } + } + } + r +} + +/// LU via Gauss Elimination with Complete Pivoting +fn gecp(m: &mut ComplexMatrix) -> (Vec, Vec) { + let n = m.col; + let mut r = vec![0usize; n - 1]; + let mut s = vec![0usize; n - 1]; + for k in 0..n - 1 { + // Find pivot + let (r_k, s_k) = match m.shape { + Shape::Col => { + let mut row_ics = 0usize; + let mut col_ics = 0usize; + let mut max_val = 0f64; + for i in k..n { + let c = m + .col(i) + .into_iter() + .skip(k) + .enumerate() + .max_by(|x1, x2| x1.1.norm().partial_cmp(&x2.1.norm()).unwrap()) + .unwrap(); + let c_ics = c.0 + k; + let c_val = c.1.norm(); + if c_val > max_val { + row_ics = c_ics; + col_ics = i; + max_val = c_val; + } + } + (row_ics, col_ics) + } + Shape::Row => { + let mut row_ics = 0usize; + let mut col_ics = 0usize; + let mut max_val = 0f64; + for i in k..n { + let c = m + .row(i) + .into_iter() + .skip(k) + .enumerate() + .max_by(|x1, x2| x1.1.norm().partial_cmp(&x2.1.norm()).unwrap()) + .unwrap(); + let c_ics = c.0 + k; + let c_val = c.1.norm(); + if c_val > max_val { + col_ics = c_ics; + row_ics = i; + max_val = c_val; + } + } + (row_ics, col_ics) + } + }; + r[k] = r_k; + s[k] = s_k; + + // Interchange rows + for j in k..n { + unsafe { + std::ptr::swap(&mut m[(k, j)], &mut m[(r_k, j)]); + } + } + + // Interchange cols + for i in 0..n { + unsafe { + std::ptr::swap(&mut m[(i, k)], &mut m[(i, s_k)]); + } + } + + // Form the multipliers + for i in k + 1..n { + m[(i, k)] = -m[(i, k)] / m[(k, k)]; + for j in k + 1..n { + let local_m = m[(i, k)] * m[(k, j)]; + m[(i, j)] += local_m; + } + } + } + (r, s) +} diff --git a/src/complex/mod.rs b/src/complex/mod.rs index 21f80b41..b17fc18f 100644 --- a/src/complex/mod.rs +++ b/src/complex/mod.rs @@ -2,4 +2,6 @@ use num_complex::Complex; pub type C64 = Complex; +pub mod integral; +pub mod matrix; pub mod vector; diff --git a/src/complex/vector.rs b/src/complex/vector.rs index 2900a4c6..e5836643 100644 --- a/src/complex/vector.rs +++ b/src/complex/vector.rs @@ -1,7 +1,8 @@ -use num_complex::Complex; +use crate::fuga::Algorithm; use crate::traits::fp::FPVector; -use crate::traits::math::{Vector, Normed, Norm, InnerProduct}; +use crate::traits::math::{InnerProduct, Norm, Normed, Vector}; use crate::traits::sugar::VecOps; +use num_complex::Complex; impl Vector for Complex { type Scalar = Self; @@ -31,7 +32,8 @@ impl Normed for Complex { fn normalize(&self, kind: Norm) -> Self where - Self: Sized { + Self: Sized, + { let n = self.norm(kind); self / n } @@ -48,26 +50,33 @@ impl FPVector for Vec> { fn fmap(&self, f: F) -> Self where - F: Fn(Self::Scalar) -> Self::Scalar { + F: Fn(Self::Scalar) -> Self::Scalar, + { self.iter().map(|&x| f(x)).collect() } fn zip_with(&self, f: F, other: &Self) -> Self where - F: Fn(Self::Scalar, Self::Scalar) -> Self::Scalar { - self.iter().zip(other.iter()).map(|(&x, &y)| f(x,y)).collect() + F: Fn(Self::Scalar, Self::Scalar) -> Self::Scalar, + { + self.iter() + .zip(other.iter()) + .map(|(&x, &y)| f(x, y)) + .collect() } fn reduce(&self, init: T, f: F) -> Self::Scalar where - F: Fn(Self::Scalar, Self::Scalar) -> Self::Scalar, - T: Into { - self.iter().fold(init.into(), |x, &y| f(x,y)) + F: Fn(Self::Scalar, Self::Scalar) -> Self::Scalar, + T: Into, + { + self.iter().fold(init.into(), |x, &y| f(x, y)) } fn filter(&self, f: F) -> Self where - F: Fn(Self::Scalar) -> bool { + F: Fn(Self::Scalar) -> bool, + { self.iter().filter(|&x| f(*x)).cloned().collect() } @@ -110,13 +119,14 @@ impl Normed for Vec> { fn norm(&self, kind: Norm) -> Self::UnsignedScalar { match kind { Norm::L1 => self.iter().map(|x| Complex::::norm(*x).abs()).sum(), - _ => unimplemented!() + _ => unimplemented!(), } } - fn normalize(&self, kind: Norm) -> Self + fn normalize(&self, _kind: Norm) -> Self where - Self: Sized { + Self: Sized, + { unimplemented!() } } @@ -128,3 +138,37 @@ impl InnerProduct for Vec> { } impl VecOps for Vec> {} + +impl Algorithm for Vec> { + type Scalar = Complex; + + fn rank(&self) -> Vec { + unimplemented!() + } + + fn sign(&self) -> Complex { + unimplemented!() + } + + fn arg_max(&self) -> usize { + unimplemented!() + } + + fn arg_min(&self) -> usize { + unimplemented!() + } + + fn max(&self) -> Complex { + unimplemented!() + } + + fn min(&self) -> Complex { + unimplemented!() + } + + fn swap_with_perm(&mut self, p: &Vec<(usize, usize)>) { + for (i, j) in p.iter() { + self.swap(*i, *j); + } + } +} diff --git a/src/numerical/integral.rs b/src/numerical/integral.rs index ac2165dd..f7f669df 100644 --- a/src/numerical/integral.rs +++ b/src/numerical/integral.rs @@ -1,4 +1,4 @@ -use crate::structure::polynomial::{lagrange_polynomial, Calculus}; +use crate::structure::polynomial::{lagrange_polynomial, Calculus, Polynomial}; use crate::traits::fp::FPVector; use crate::util::non_macro::seq; @@ -132,6 +132,78 @@ impl Integral { } } +// Types that can be integrated using Newton Cotes Quadrature +pub trait NCIntegrable: std::ops::Sub + Sized { + type NodeY; + type NCPolynomial; + + fn compute_node_y(f: F, node_x: &[f64]) -> Self::NodeY + where + F: Fn(f64) -> Self; + + fn compute_polynomial(node_x: &[f64], node_y: &Self::NodeY) -> Self::NCPolynomial; + + fn integrate_polynomial(p: &Self::NCPolynomial) -> Self::NCPolynomial; + + fn evaluate_polynomial(p: &Self::NCPolynomial, x: f64) -> Self; +} + +impl NCIntegrable for f64 { + type NodeY = Vec; + type NCPolynomial = Polynomial; + + fn compute_node_y(f: F, node_x: &[f64]) -> Vec + where + F: Fn(f64) -> Self, + { + node_x.to_vec().fmap(|x| f(x)) + } + + fn compute_polynomial(node_x: &[f64], node_y: &Vec) -> Polynomial { + lagrange_polynomial(node_x.to_vec(), node_y.to_vec()) + } + + fn integrate_polynomial(p: &Polynomial) -> Polynomial { + p.integral() + } + + fn evaluate_polynomial(p: &Polynomial, x: f64) -> Self { + p.eval(x) + } +} + +// Types that can be integrated using Gauss Lagrange or Kronrod Quadrature +pub trait GLKIntegrable: + std::ops::Add + std::ops::Mul + Sized +{ + const ZERO: Self; +} + +impl GLKIntegrable for f64 { + const ZERO: Self = 0f64; +} + +// Types that can be integrated using Gauss Kronrod Quadrature +pub trait GKIntegrable: GLKIntegrable + std::ops::Sub + Sized + Clone { + fn is_finite(&self) -> bool; + + fn gk_norm(&self) -> f64; + + fn is_within_tol(&self, tol: f64) -> bool { + self.gk_norm() < tol + } +} + +impl GKIntegrable for f64 { + fn is_finite(&self) -> bool { + f64::is_finite(*self) + } + + fn gk_norm(&self) -> f64 { + self.abs() + } +} + /// Numerical Integration /// /// # Description @@ -159,29 +231,31 @@ impl Integral { /// * `G20K41R` /// * `G25K51R` /// * `G30K61R` -pub fn integrate(f: F, (a, b): (f64, f64), method: Integral) -> f64 +pub fn integrate(f: F, (a, b): (f64, f64), method: Integral) -> Y where - F: Fn(f64) -> f64 + Copy, + F: Fn(f64) -> Y + Copy, + Y: GKIntegrable + NCIntegrable, { match method { Integral::GaussLegendre(n) => gauss_legendre_quadrature(f, n, (a, b)), Integral::NewtonCotes(n) => newton_cotes_quadrature(f, n, (a, b)), - method => gauss_kronrod_quadrature(f, (a,b), method), + method => gauss_kronrod_quadrature(f, (a, b), method), } } /// Newton Cotes Quadrature -pub fn newton_cotes_quadrature(f: F, n: usize, (a, b): (f64, f64)) -> f64 +pub fn newton_cotes_quadrature(f: F, n: usize, (a, b): (f64, f64)) -> Y where - F: Fn(f64) -> f64, + F: Fn(f64) -> Y, + Y: NCIntegrable, { let h = (b - a) / (n as f64); let node_x = seq(a, b, h); - let node_y = node_x.fmap(|x| f(x)); - let p = lagrange_polynomial(node_x, node_y); - let q = p.integral(); + let node_y = Y::compute_node_y(f, &node_x); + let p = Y::compute_polynomial(&node_x, &node_y); + let q = Y::integrate_polynomial(&p); - q.eval(b) - q.eval(a) + Y::evaluate_polynomial(&q, b) - Y::evaluate_polynomial(&q, a) } /// Gauss Legendre Quadrature @@ -195,11 +269,12 @@ where /// # Reference /// * A. N. Lowan et al. (1942) /// * [Keisan Online Calculator](https://keisan.casio.com/exec/system/1329114617) -pub fn gauss_legendre_quadrature(f: F, n: usize, (a, b): (f64, f64)) -> f64 +pub fn gauss_legendre_quadrature(f: F, n: usize, (a, b): (f64, f64)) -> Y where - F: Fn(f64) -> f64, + F: Fn(f64) -> Y, + Y: GLKIntegrable, { - (b - a) / 2f64 * unit_gauss_legendre_quadrature(|x| f(x * (b - a) / 2f64 + (a + b) / 2f64), n) + unit_gauss_legendre_quadrature(|x| f(x * (b - a) / 2f64 + (a + b) / 2f64), n) * ((b - a) / 2f64) } /// Gauss Kronrod Quadrature @@ -215,16 +290,17 @@ where /// * arXiv: [1003.4629](https://arxiv.org/abs/1003.4629) /// * [Keisan Online Calculator](https://keisan.casio.com/exec/system/1329114617) #[allow(non_snake_case)] -pub fn gauss_kronrod_quadrature(f: F, (a, b): (T, S), method: Integral) -> f64 +pub fn gauss_kronrod_quadrature(f: F, (a, b): (T, S), method: Integral) -> Y where - F: Fn(f64) -> f64 + Copy, - T: Into, - S: Into, + F: Fn(f64) -> Y + Copy, + T: Into, + S: Into, + Y: GKIntegrable, { let (g, k) = method.get_gauss_kronrod_order(); let tol = method.get_tol(); let max_iter = method.get_max_iter(); - let mut I = 0f64; + let mut I = Y::ZERO; let mut S: Vec<(f64, f64, f64, u32)> = vec![]; S.push((a.into(), b.into(), tol, max_iter)); @@ -235,15 +311,15 @@ where let K = kronrod_quadrature(f, k as usize, (a, b)); let c = (a + b) / 2f64; let tol_curr = if method.is_relative() { - tol * G + tol * G.gk_norm() } else { tol }; - if (G - K).abs() < tol_curr || a == b || max_iter == 0 { - if ! G.is_finite() { + if (G.clone() - K).is_within_tol(tol_curr) || a == b || max_iter == 0 { + if !G.is_finite() { return G; } - I += G; + I = I + G; } else { S.push((a, c, tol / 2f64, max_iter - 1)); S.push((c, b, tol / 2f64, max_iter - 1)); @@ -255,24 +331,26 @@ where I } -pub fn kronrod_quadrature(f: F, n: usize, (a, b): (f64, f64)) -> f64 +pub fn kronrod_quadrature(f: F, n: usize, (a, b): (f64, f64)) -> Y where - F: Fn(f64) -> f64, + F: Fn(f64) -> Y, + Y: GLKIntegrable, { - (b - a) / 2f64 * unit_kronrod_quadrature(|x| f(x * (b-a) / 2f64 + (a + b) / 2f64), n) + unit_kronrod_quadrature(|x| f(x * (b - a) / 2f64 + (a + b) / 2f64), n) * ((b - a) / 2f64) } // ============================================================================= // Gauss Legendre Backends // ============================================================================= -fn unit_gauss_legendre_quadrature(f: F, n: usize) -> f64 +fn unit_gauss_legendre_quadrature(f: F, n: usize) -> Y where - F: Fn(f64) -> f64, + F: Fn(f64) -> Y, + Y: GLKIntegrable, { let (a, x) = gauss_legendre_table(n); - let mut s = 0f64; + let mut s = Y::ZERO; for i in 0..a.len() { - s += a[i] * f(x[i]); + s = s + f(x[i]) * a[i]; } s } @@ -378,14 +456,15 @@ fn gauss_legendre_table(n: usize) -> (Vec, Vec) { // Gauss Kronrod Backends // ============================================================================= -fn unit_kronrod_quadrature(f: F, n: usize) -> f64 +fn unit_kronrod_quadrature(f: F, n: usize) -> Y where - F: Fn(f64) -> f64, + F: Fn(f64) -> Y, + Y: GLKIntegrable, { - let (a,x) = kronrod_table(n); - let mut s = 0f64; - for i in 0 .. a.len() { - s += a[i] * f(x[i]); + let (a, x) = kronrod_table(n); + let mut s = Y::ZERO; + for i in 0..a.len() { + s = s + f(x[i]) * a[i]; } s } @@ -414,28 +493,28 @@ fn kronrod_table(n: usize) -> (Vec, Vec) { match n % 2 { 0 => { - for i in 0 .. ref_node.len() { + for i in 0..ref_node.len() { result_node[i] = ref_node[i]; result_weight[i] = ref_weight[i]; } - for i in ref_node.len() .. n { - result_node[i] = -ref_node[n-i-1]; - result_weight[i] = ref_weight[n-i-1]; + for i in ref_node.len()..n { + result_node[i] = -ref_node[n - i - 1]; + result_weight[i] = ref_weight[n - i - 1]; } } 1 => { - for i in 0 .. ref_node.len() { + for i in 0..ref_node.len() { result_node[i] = ref_node[i]; result_weight[i] = ref_weight[i]; } - for i in ref_node.len() .. n { - result_node[i] = -ref_node[n-i]; - result_weight[i] = ref_weight[n-i]; + for i in ref_node.len()..n { + result_node[i] = -ref_node[n - i]; + result_weight[i] = ref_weight[n - i]; } } - _ => unreachable!() + _ => unreachable!(), } (result_weight, result_node) } @@ -922,7 +1001,7 @@ const LEGENDRE_WEIGHT_25: [f64; 13] = [ 0.054904695975835192, 0.040939156701306313, 0.026354986615032137, - 0.011393798501026288, + 0.011393798501026288, ]; const LEGENDRE_WEIGHT_26: [f64; 13] = [ 0.118321415279262277, diff --git a/src/structure/matrix.rs b/src/structure/matrix.rs index 797792ca..0dc55e7c 100644 --- a/src/structure/matrix.rs +++ b/src/structure/matrix.rs @@ -192,7 +192,7 @@ //! To write columns or rows, `DataFrame` and `nc` feature could be the best choice. //! //! * `nc` feature should be required - `netcdf` or `libnetcdf` are prerequisites. -//! +//! //! ```no_run //! #[macro_use] //! extern crate peroxide; @@ -596,10 +596,10 @@ //! } //! ``` -#[cfg(feature="csv")] +#[cfg(feature = "csv")] extern crate csv; -#[cfg(feature="csv")] +#[cfg(feature = "csv")] use self::csv::{ReaderBuilder, StringRecord, WriterBuilder}; #[cfg(feature = "O3")] @@ -609,30 +609,30 @@ extern crate lapack; #[cfg(feature = "O3")] use blas::{daxpy, dgemm, dgemv}; #[cfg(feature = "O3")] -use lapack::{dgecon, dgeqrf, dgetrf, dgetri, dgetrs, dorgqr, dgesvd, dpotrf}; +use lapack::{dgecon, dgeqrf, dgesvd, dgetrf, dgetri, dgetrs, dorgqr, dpotrf}; #[cfg(feature = "serde")] use serde::{Deserialize, Serialize}; pub use self::Shape::{Col, Row}; use crate::numerical::eigen::{eigen, EigenMethod}; +use crate::structure::dataframe::{Series, TypedVector}; +use crate::traits::sugar::ScalableMut; use crate::traits::{ - general::Algorithm, fp::{FPMatrix, FPVector}, + general::Algorithm, math::{InnerProduct, LinearOp, MatrixProduct, Norm, Normed, Vector}, mutable::MutMatrix, }; use crate::util::{ - low_level::{swap_vec_ptr, copy_vec_ptr}, + low_level::{copy_vec_ptr, swap_vec_ptr}, non_macro::{cbind, eye, rbind, zeros}, useful::{nearly_eq, tab}, }; -use crate::structure::dataframe::{Series, TypedVector}; +use peroxide_num::{ExpLogOps, Numeric, PowOps, TrigOps}; use std::cmp::{max, min}; pub use std::error::Error; use std::fmt; use std::ops::{Add, Div, Index, IndexMut, Mul, Neg, Sub}; -use crate::traits::sugar::ScalableMut; -use peroxide_num::{ExpLogOps, PowOps, TrigOps, Numeric}; pub type Perms = Vec<(usize, usize)>; @@ -1143,7 +1143,7 @@ impl Matrix { /// Ok(()) /// } /// ``` - #[cfg(feature="csv")] + #[cfg(feature = "csv")] pub fn write(&self, file_path: &str) -> Result<(), Box> { let mut wtr = WriterBuilder::new().from_path(file_path)?; let r = self.row; @@ -1174,7 +1174,7 @@ impl Matrix { /// Ok(()) /// } /// ``` - #[cfg(feature="csv")] + #[cfg(feature = "csv")] pub fn write_round(&self, file_path: &str, round: usize) -> Result<(), Box> { let mut wtr = WriterBuilder::new().from_path(file_path)?; let r = self.row; @@ -1190,7 +1190,7 @@ impl Matrix { Ok(()) } - #[cfg(feature="csv")] + #[cfg(feature = "csv")] pub fn write_with_header( &self, file_path: &str, @@ -1212,7 +1212,7 @@ impl Matrix { Ok(()) } - #[cfg(feature="csv")] + #[cfg(feature = "csv")] pub fn write_with_header_round( &self, file_path: &str, @@ -1260,7 +1260,7 @@ impl Matrix { /// Ok(()) /// } /// ``` - #[cfg(feature="csv")] + #[cfg(feature = "csv")] pub fn read(file_path: &str, header: bool, delimiter: char) -> Result> { let mut rdr = ReaderBuilder::new() .has_headers(header) @@ -1365,7 +1365,7 @@ impl Matrix { /// /// # Description /// Return below elements of matrix to new matrix - /// + /// /// $$ /// \begin{pmatrix} /// \\ddots & & & & \\\\ @@ -1380,7 +1380,7 @@ impl Matrix { /// ``` /// 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"); @@ -1392,8 +1392,8 @@ impl 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 { + for i in 0..row { + for j in 0..col { result[(i, j)] = self[(start.0 + i, start.1 + j)]; } } @@ -1404,7 +1404,7 @@ impl Matrix { /// /// # Description /// Substitute below elements of matrix - /// + /// /// $$ /// \begin{pmatrix} /// \\ddots & & & & \\\\ @@ -1414,12 +1414,12 @@ impl Matrix { /// & & & & \\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"); @@ -1431,8 +1431,8 @@ impl Matrix { 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 { + for i in 0..row { + for j in 0..col { self[(start.0 + i, start.1 + j)] = m[(i, j)]; } } @@ -2561,6 +2561,8 @@ impl IndexMut<(usize, usize)> for Matrix { // ============================================================================= impl FPMatrix for Matrix { + type Scalar = f64; + fn take_row(&self, n: usize) -> Self { if n >= self.row { return self.clone(); @@ -2918,7 +2920,7 @@ impl SVD { pub fn s_mat(&self) -> Matrix { let mut mat = zeros(self.u.col, self.vt.row); - for i in 0 .. mat.row.min(mat.col) { + for i in 0..mat.row.min(mat.col) { mat[(i, i)] = self.s[i]; } mat @@ -2961,11 +2963,7 @@ impl SVD { } } - SVD { - s, - u, - vt, - } + SVD { s, u, vt } } } @@ -3151,7 +3149,7 @@ impl LinearAlgebra for Matrix { #[allow(non_snake_case)] fn qr(&self) -> QR { match () { - #[cfg(feature="O3")] + #[cfg(feature = "O3")] () => { let opt_dgeqrf = lapack_dgeqrf(self); match opt_dgeqrf { @@ -3159,10 +3157,7 @@ impl LinearAlgebra for Matrix { Some(dgeqrf) => { let q = dgeqrf.get_Q(); let r = dgeqrf.get_R(); - QR { - q, - r - } + QR { q, r } } } } @@ -3209,7 +3204,7 @@ impl LinearAlgebra for Matrix { /// ``` fn svd(&self) -> SVD { match () { - #[cfg(feature="O3")] + #[cfg(feature = "O3")] () => { let opt_dgesvd = lapack_dgesvd(self); match opt_dgesvd { @@ -3218,14 +3213,12 @@ impl LinearAlgebra for Matrix { SVD_STATUS::Diverge(i) => { panic!("Divergent occurs in SVD - {} iterations", i) } - SVD_STATUS::Success => { - SVD { - s: dgesvd.s, - u: dgesvd.u, - vt: dgesvd.vt, - } - } - } + SVD_STATUS::Success => SVD { + s: dgesvd.s, + u: dgesvd.u, + vt: dgesvd.vt, + }, + }, } } _ => { @@ -3465,7 +3458,9 @@ impl LinearAlgebra for Matrix { None => panic!("Singular matrix (opt_dgrf)"), Some(dgrf) => match dgrf.status { LAPACK_STATUS::NonSingular => lapack_dgetri(&dgrf).unwrap(), - LAPACK_STATUS::Singular => panic!("Singular matrix (LAPACK_STATUS Singular)"), + LAPACK_STATUS::Singular => { + panic!("Singular matrix (LAPACK_STATUS Singular)") + } }, } } @@ -3500,13 +3495,13 @@ impl LinearAlgebra for Matrix { /// ``` fn pseudo_inv(&self) -> Self { match () { - #[cfg(feature="O3")] + #[cfg(feature = "O3")] () => { let svd = self.svd(); let row = svd.vt.row; let col = svd.u.row; let mut sp = zeros(row, col); - for i in 0 .. row.min(col) { + for i in 0..row.min(col) { sp[(i, i)] = 1f64 / svd.s[i]; } svd.vt.t() * sp * svd.u.t() @@ -3620,9 +3615,9 @@ impl LinearAlgebra for Matrix { return false; } - for i in 0 .. self.row { - for j in i .. self.col { - if !nearly_eq(self[(i,j)], self[(j,i)]) { + for i in 0..self.row { + for j in i..self.col { + if !nearly_eq(self[(i, j)], self[(j, i)]) { return false; } } @@ -3637,6 +3632,8 @@ pub fn solve(A: &Matrix, b: &Matrix, sk: SolveKind) -> Matrix { } impl MutMatrix for Matrix { + type Scalar = f64; + unsafe fn col_mut(&mut self, idx: usize) -> Vec<*mut f64> { assert!(idx < self.col, "Index out of range"); match self.shape { @@ -3722,7 +3719,6 @@ impl Div for Matrix { } } - impl ExpLogOps for Matrix { type Float = f64; fn exp(&self) -> Self { @@ -4248,7 +4244,7 @@ pub enum POSITIVE_STATUS { #[derive(Debug, Copy, Clone, Eq, PartialEq)] pub enum UPLO { Upper, - Lower + Lower, } /// Temporary data structure from `dgetrf` @@ -4279,7 +4275,7 @@ pub struct DGESVD { pub struct DPOTRF { pub fact_mat: Matrix, pub uplo: UPLO, - pub status: POSITIVE_STATUS + pub status: POSITIVE_STATUS, } ///// Temporary data structure from `dgeev` @@ -4473,7 +4469,22 @@ pub fn lapack_dgesvd(mat: &Matrix) -> Option { // Workspace query unsafe { - dgesvd(jobu, jobvt, m, n, &mut A.data, lda, &mut s, &mut u, ldu, &mut vt, ldvt, &mut work, lwork, &mut info); + dgesvd( + jobu, + jobvt, + m, + n, + &mut A.data, + lda, + &mut s, + &mut u, + ldu, + &mut vt, + ldvt, + &mut work, + lwork, + &mut info, + ); } let optimal_lwork = work[0] as usize; @@ -4532,35 +4543,23 @@ pub fn lapack_dpotrf(mat: &Matrix, UPLO: UPLO) -> Option { let mut info = 0i32; let uplo = match UPLO { UPLO::Upper => b'U', - UPLO::Lower => b'L' + UPLO::Lower => b'L', }; - unsafe { - dpotrf( - uplo, - N, - &mut A.data, - lda, - &mut info, - ) - } + unsafe { dpotrf(uplo, N, &mut A.data, lda, &mut info) } if info == 0 { - Some( - DPOTRF { - fact_mat: matrix(A.data, mat.row, mat.col, Col), - uplo: UPLO, - status: POSITIVE_STATUS::Success - } - ) + Some(DPOTRF { + fact_mat: matrix(A.data, mat.row, mat.col, Col), + uplo: UPLO, + status: POSITIVE_STATUS::Success, + }) } else if info > 0 { - Some( - DPOTRF { - fact_mat: matrix(A.data, mat.row, mat.col, Col), - uplo: UPLO, - status: POSITIVE_STATUS::Failed(info) - } - ) + Some(DPOTRF { + fact_mat: matrix(A.data, mat.row, mat.col, Col), + uplo: UPLO, + status: POSITIVE_STATUS::Failed(info), + }) } else { None } @@ -4674,8 +4673,8 @@ impl DPOTRF { let mat = &self.fact_mat; let n = mat.col; let mut result = matrix(vec![0f64; n.pow(2)], n, n, mat.shape); - for i in 0 .. n { - for j in i .. n { + for i in 0..n { + for j in i..n { result[(i, j)] = mat[(i, j)]; } } @@ -4691,8 +4690,8 @@ impl DPOTRF { let n = mat.col; let mut result = matrix(vec![0f64; n.pow(2)], n, n, mat.shape); - for i in 0 .. n { - for j in 0 .. i+1 { + for i in 0..n { + for j in 0..i + 1 { result[(i, j)] = mat[(i, j)]; } } diff --git a/src/structure/vector.rs b/src/structure/vector.rs index 0e977361..9a0dcd71 100644 --- a/src/structure/vector.rs +++ b/src/structure/vector.rs @@ -473,6 +473,8 @@ impl MutFP for Vec { } impl Algorithm for Vec { + type Scalar = f64; + /// Assign rank /// /// # Examples diff --git a/src/traits/fp.rs b/src/traits/fp.rs index 41f11463..a773583d 100644 --- a/src/traits/fp.rs +++ b/src/traits/fp.rs @@ -1,5 +1,3 @@ -use crate::structure::matrix::Matrix; - /// Functional Programming tools for Vector pub trait FPVector { type Scalar; @@ -23,38 +21,40 @@ pub trait FPVector { fn prod(&self) -> Self::Scalar; } -/// Functional Programming for Matrix +/// Functional Programming for Matrix and ComplexMatrix pub trait FPMatrix { - fn take_row(&self, n: usize) -> Matrix; - fn take_col(&self, n: usize) -> Matrix; - fn skip_row(&self, n: usize) -> Matrix; - fn skip_col(&self, n: usize) -> Matrix; - fn fmap(&self, f: F) -> Matrix + type Scalar; + + fn take_row(&self, n: usize) -> Self; + fn take_col(&self, n: usize) -> Self; + fn skip_row(&self, n: usize) -> Self; + fn skip_col(&self, n: usize) -> Self; + fn fmap(&self, f: F) -> Self where - F: Fn(f64) -> f64; - fn col_map(&self, f: F) -> Matrix + F: Fn(Self::Scalar) -> Self::Scalar; + fn col_map(&self, f: F) -> Self where - F: Fn(Vec) -> Vec; - fn row_map(&self, f: F) -> Matrix + F: Fn(Vec) -> Vec; + fn row_map(&self, f: F) -> Self where - F: Fn(Vec) -> Vec; + F: Fn(Vec) -> Vec; fn col_mut_map(&mut self, f: F) where - F: Fn(Vec) -> Vec; + F: Fn(Vec) -> Vec; fn row_mut_map(&mut self, f: F) where - F: Fn(Vec) -> Vec; - fn reduce(&self, init: T, f: F) -> f64 + F: Fn(Vec) -> Vec; + fn reduce(&self, init: T, f: F) -> Self::Scalar where - F: Fn(f64, f64) -> f64, - T: Into; - fn zip_with(&self, f: F, other: &Matrix) -> Matrix + F: Fn(Self::Scalar, Self::Scalar) -> Self::Scalar, + T: Into; + fn zip_with(&self, f: F, other: &Self) -> Self where - F: Fn(f64, f64) -> f64; - fn col_reduce(&self, f: F) -> Vec + F: Fn(Self::Scalar, Self::Scalar) -> Self::Scalar; + fn col_reduce(&self, f: F) -> Vec where - F: Fn(Vec) -> f64; - fn row_reduce(&self, f: F) -> Vec + F: Fn(Vec) -> Self::Scalar; + fn row_reduce(&self, f: F) -> Vec where - F: Fn(Vec) -> f64; + F: Fn(Vec) -> Self::Scalar; } diff --git a/src/traits/general.rs b/src/traits/general.rs index 2737f97b..4712d304 100644 --- a/src/traits/general.rs +++ b/src/traits/general.rs @@ -1,10 +1,12 @@ /// Some algorithms for Vector pub trait Algorithm { + type Scalar; + fn rank(&self) -> Vec; - fn sign(&self) -> f64; + fn sign(&self) -> Self::Scalar; fn arg_max(&self) -> usize; fn arg_min(&self) -> usize; - fn max(&self) -> f64; - fn min(&self) -> f64; + fn max(&self) -> Self::Scalar; + fn min(&self) -> Self::Scalar; fn swap_with_perm(&mut self, p: &Vec<(usize, usize)>); } diff --git a/src/traits/math.rs b/src/traits/math.rs index 228aad9e..4e05e671 100644 --- a/src/traits/math.rs +++ b/src/traits/math.rs @@ -60,8 +60,8 @@ pub trait VectorProduct: Vector { /// Matrix Products pub trait MatrixProduct { - fn kronecker(&self, other: &Self) -> Matrix; - fn hadamard(&self, other: &Self) -> Matrix; + fn kronecker(&self, other: &Self) -> Self; + fn hadamard(&self, other: &Self) -> Self; } // ============================================================================= diff --git a/src/traits/mutable.rs b/src/traits/mutable.rs index 7b0932fc..f7f098c2 100644 --- a/src/traits/mutable.rs +++ b/src/traits/mutable.rs @@ -11,8 +11,10 @@ pub trait MutFP { } pub trait MutMatrix { - unsafe fn col_mut(&mut self, idx: usize) -> Vec<*mut f64>; - unsafe fn row_mut(&mut self, idx: usize) -> Vec<*mut f64>; + type Scalar; + + unsafe fn col_mut(&mut self, idx: usize) -> Vec<*mut Self::Scalar>; + unsafe fn row_mut(&mut self, idx: usize) -> Vec<*mut Self::Scalar>; unsafe fn swap(&mut self, idx1: usize, idx2: usize, shape: Shape); unsafe fn swap_with_perm(&mut self, p: &Vec<(usize, usize)>, shape: Shape); } diff --git a/src/util/low_level.rs b/src/util/low_level.rs index 8b2a1f52..53fe5016 100644 --- a/src/util/low_level.rs +++ b/src/util/low_level.rs @@ -1,17 +1,26 @@ -pub unsafe fn copy_vec_ptr(dst: &mut Vec<*mut f64>, src: &Vec) { +pub unsafe fn copy_vec_ptr(dst: &mut Vec<*mut T>, src: &Vec) +where + T: Copy, +{ assert_eq!(dst.len(), src.len(), "Should use same length vectors"); for (&mut p, &s) in dst.iter_mut().zip(src) { *p = s; } } -pub unsafe fn swap_vec_ptr(lhs: &mut Vec<*mut f64>, rhs: &mut Vec<*mut f64>) { +pub unsafe fn swap_vec_ptr(lhs: &mut Vec<*mut T>, rhs: &mut Vec<*mut T>) +where + T: Copy, +{ assert_eq!(lhs.len(), rhs.len(), "Should use same length vectors"); for (&mut l, &mut r) in lhs.iter_mut().zip(rhs.iter_mut()) { std::ptr::swap(l, r); } } -pub unsafe fn ptr_to_vec(pv: &Vec<*const f64>) -> Vec { +pub unsafe fn ptr_to_vec(pv: &Vec<*const T>) -> Vec +where + T: Copy, +{ pv.iter().map(|&x| *x).collect() } diff --git a/tests/complex_matrix.rs b/tests/complex_matrix.rs new file mode 100644 index 00000000..113c413e --- /dev/null +++ b/tests/complex_matrix.rs @@ -0,0 +1,23 @@ +extern crate peroxide; +#[allow(unused_imports)] +use peroxide::fuga::*; + +#[cfg(feature = "complex")] +#[test] +fn test_seq() { + 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), + ], + row: 2, + col: 2, + shape: Row, + }; + assert_eq!(v1.data[0], Complex64::new(1f64, 1f64)); +}