diff --git a/Cargo.toml b/Cargo.toml index 87a8dbbe..1a120447 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "peroxide" -version = "0.36.0" +version = "0.36.1" authors = ["axect "] edition = "2018" description = "Rust comprehensive scientific computation library contains linear algebra, numerical analysis, statistics and machine learning tools with farmiliar syntax" diff --git a/RELEASES.md b/RELEASES.md index f3df7bbd..0793cd66 100644 --- a/RELEASES.md +++ b/RELEASES.md @@ -1,3 +1,39 @@ +# Release 0.36.1 (2024-04-09) + +- Fix all warnings in peroxide +- Change redundant method + - `Vec::resize` -> `Vec::reshape` +- Error handling for concatenation + - `cbind` & `rbind` now returns `Result` +- New non-macro utils + - `column_stack(&[Vec]) -> Result` + - `row_stack(&[Vec]) -> Result` + - `rand_with_rng(usize, usize, &mut Rng) -> Matrix` +- Generic Butcher tableau trait (now for embedded Runge-Kutta methods) + + ```rust + pub trait ButcherTableau { + const C: &'static [f64]; + const A: &'static [&'static [f64]]; + const BH: &'static [f64]; + const BL: &'static [f64]; + + fn tol(&self) -> f64; + fn safety_factor(&self) -> f64; + fn max_step_size(&self) -> f64; + fn min_step_size(&self) -> f64; + fn max_step_iter(&self) -> usize; + } + ``` + + - Implement `ODEIntegrator` for `ButcherTableau` + - Just declare `ButcherTableau` then `step` is free + + - Three available embedded Runge-Kutta methods + - `RKF45`: Runge-Kutta-Fehlberg 4/5th order + - `DP45`: Dormand-Prince 4/5th order + - `TSIT45`: Tsitouras 4/5th order + # Release 0.36.0 (2024-04-08) ## Huge Update - Error handling & Whole new ODE diff --git a/example_data/lorenz_dp45.png b/example_data/lorenz_dp45.png new file mode 100644 index 00000000..e8c2af97 Binary files /dev/null and b/example_data/lorenz_dp45.png differ diff --git a/example_data/lorenz_euler.png b/example_data/lorenz_euler.png deleted file mode 100644 index 76eddc57..00000000 Binary files a/example_data/lorenz_euler.png and /dev/null differ diff --git a/example_data/lorenz_rkf45.png b/example_data/lorenz_rkf45.png index 6d9f94fe..fe7f21b1 100644 Binary files a/example_data/lorenz_rkf45.png and b/example_data/lorenz_rkf45.png differ diff --git a/example_data/lorenz_tsit45.png b/example_data/lorenz_tsit45.png new file mode 100644 index 00000000..a2c3e4c5 Binary files /dev/null and b/example_data/lorenz_tsit45.png differ diff --git a/examples/lorenz_dp45.rs b/examples/lorenz_dp45.rs new file mode 100644 index 00000000..9f3e5624 --- /dev/null +++ b/examples/lorenz_dp45.rs @@ -0,0 +1,47 @@ +use peroxide::fuga::*; + +#[allow(unused_variables)] +fn main() -> Result<(), Box> { + let dp45 = DP45::new(1e-4, 0.9, 1e-6, 1e-1, 100); + let basic_ode_solver = BasicODESolver::new(dp45); + let (_, y_vec) = basic_ode_solver.solve( + &Lorenz, + (0f64, 100f64), + 1e-2, + )?; + let y_mat = py_matrix(y_vec); + let y0 = y_mat.col(0); + let y2 = y_mat.col(2); + + #[cfg(feature = "plot")] + { + let mut plt = Plot2D::new(); + plt + .set_domain(y0) + .insert_image(y2) + .set_xlabel(r"$y_0$") + .set_ylabel(r"$y_2$") + .set_style(PlotStyle::Nature) + .tight_layout() + .set_dpi(600) + .set_path("example_data/lorenz_dp45.png") + .savefig()?; + } + + Ok(()) +} + +struct Lorenz; + +impl ODEProblem for Lorenz { + fn initial_conditions(&self) -> Vec { + vec![10f64, 1f64, 1f64] + } + + fn rhs(&self, _t: f64, y: &[f64], dy: &mut [f64]) -> Result<(), ODEError> { + dy[0] = 10f64 * (y[1] - y[0]); + dy[1] = 28f64 * y[0] - y[1] - y[0] * y[2]; + dy[2] = -8f64 / 3f64 * y[2] + y[0] * y[1]; + Ok(()) + } +} diff --git a/examples/lorenz_gl4.rs b/examples/lorenz_gl4.rs index e0937bf0..7b4477d0 100644 --- a/examples/lorenz_gl4.rs +++ b/examples/lorenz_gl4.rs @@ -1,7 +1,8 @@ use peroxide::fuga::*; +#[allow(unused_variables)] fn main() -> Result<(), Box> { - let gl4 = GL4::new(ImplicitSolver::FixedPoint, 100, 1e-6); + let gl4 = GL4::new(ImplicitSolver::FixedPoint, 1e-6, 100); let basic_ode_solver = BasicODESolver::new(gl4); let (_, y_vec) = basic_ode_solver.solve( &Lorenz, diff --git a/examples/lorenz_rk4.rs b/examples/lorenz_rk4.rs index ad016087..c4e51a28 100644 --- a/examples/lorenz_rk4.rs +++ b/examples/lorenz_rk4.rs @@ -1,5 +1,6 @@ use peroxide::fuga::*; +#[allow(unused_variables)] fn main() -> Result<(), Box> { let basic_ode_solver = BasicODESolver::new(RK4); let (_, y_vec) = basic_ode_solver.solve( @@ -36,7 +37,7 @@ impl ODEProblem for Lorenz { vec![10f64, 1f64, 1f64] } - fn rhs(&self, t: f64, y: &[f64], dy: &mut [f64]) -> Result<(), ODEError> { + fn rhs(&self, _t: f64, y: &[f64], dy: &mut [f64]) -> Result<(), ODEError> { dy[0] = 10f64 * (y[1] - y[0]); dy[1] = 28f64 * y[0] - y[1] - y[0] * y[2]; dy[2] = -8f64 / 3f64 * y[2] + y[0] * y[1]; diff --git a/examples/lorenz_rkf45.rs b/examples/lorenz_rkf45.rs index db294511..2184049c 100644 --- a/examples/lorenz_rkf45.rs +++ b/examples/lorenz_rkf45.rs @@ -1,5 +1,6 @@ use peroxide::fuga::*; +#[allow(unused_variables)] fn main() -> Result<(), Box> { let rkf45 = RKF45::new(1e-4, 0.9, 1e-6, 1e-1, 100); let basic_ode_solver = BasicODESolver::new(rkf45); @@ -37,7 +38,7 @@ impl ODEProblem for Lorenz { vec![10f64, 1f64, 1f64] } - fn rhs(&self, t: f64, y: &[f64], dy: &mut [f64]) -> Result<(), ODEError> { + fn rhs(&self, _t: f64, y: &[f64], dy: &mut [f64]) -> Result<(), ODEError> { dy[0] = 10f64 * (y[1] - y[0]); dy[1] = 28f64 * y[0] - y[1] - y[0] * y[2]; dy[2] = -8f64 / 3f64 * y[2] + y[0] * y[1]; diff --git a/examples/lorenz_tsit45.rs b/examples/lorenz_tsit45.rs new file mode 100644 index 00000000..748a732a --- /dev/null +++ b/examples/lorenz_tsit45.rs @@ -0,0 +1,47 @@ +use peroxide::fuga::*; + +#[allow(unused_variables)] +fn main() -> Result<(), Box> { + let tsit45 = TSIT45::new(1e-2, 0.9, 1e-6, 1e-1, 100); + let basic_ode_solver = BasicODESolver::new(tsit45); + let (_, y_vec) = basic_ode_solver.solve( + &Lorenz, + (0f64, 100f64), + 1e-2, + )?; + let y_mat = py_matrix(y_vec); + let y0 = y_mat.col(0); + let y2 = y_mat.col(2); + + #[cfg(feature = "plot")] + { + let mut plt = Plot2D::new(); + plt + .set_domain(y0) + .insert_image(y2) + .set_xlabel(r"$y_0$") + .set_ylabel(r"$y_2$") + .set_style(PlotStyle::Nature) + .tight_layout() + .set_dpi(600) + .set_path("example_data/lorenz_tsit45.png") + .savefig()?; + } + + Ok(()) +} + +struct Lorenz; + +impl ODEProblem for Lorenz { + fn initial_conditions(&self) -> Vec { + vec![10f64, 1f64, 1f64] + } + + fn rhs(&self, _t: f64, y: &[f64], dy: &mut [f64]) -> Result<(), ODEError> { + dy[0] = 10f64 * (y[1] - y[0]); + dy[1] = 28f64 * y[0] - y[1] - y[0] * y[2]; + dy[2] = -8f64 / 3f64 * y[2] + y[0] * y[1]; + Ok(()) + } +} diff --git a/examples/ode_env.rs b/examples/ode_env.rs index fa9c4af1..58985adc 100644 --- a/examples/ode_env.rs +++ b/examples/ode_env.rs @@ -1,5 +1,6 @@ use peroxide::fuga::*; +#[allow(unused_variables)] fn main() -> Result<(), Box> { let t = seq(0, 10, 1); let y = t diff --git a/examples/ode_test_gl4.rs b/examples/ode_test_gl4.rs index bc1f9682..a1e388f7 100644 --- a/examples/ode_test_gl4.rs +++ b/examples/ode_test_gl4.rs @@ -1,9 +1,8 @@ -#[macro_use] -extern crate peroxide; use peroxide::fuga::*; +#[allow(unused_variables)] fn main() -> Result<(), Box> { - let gl4 = GL4::new(ImplicitSolver::FixedPoint, 100, 1e-6); + let gl4 = GL4::new(ImplicitSolver::FixedPoint, 1e-6, 100); let basic_ode_solver = BasicODESolver::new(gl4); let (t_vec, y_vec) = basic_ode_solver.solve( &Test, diff --git a/examples/ode_test_rk4.rs b/examples/ode_test_rk4.rs index 2a22a0b4..879c07a3 100644 --- a/examples/ode_test_rk4.rs +++ b/examples/ode_test_rk4.rs @@ -1,5 +1,6 @@ use peroxide::fuga::*; +#[allow(unused_variables)] fn main() -> Result<(), Box> { let basic_ode_solver = BasicODESolver::new(RK4); let (t_vec, y_vec) = basic_ode_solver.solve( diff --git a/examples/ode_test_rkf45.rs b/examples/ode_test_rkf45.rs index bdd64bcd..8c13e148 100644 --- a/examples/ode_test_rkf45.rs +++ b/examples/ode_test_rkf45.rs @@ -1,5 +1,6 @@ use peroxide::fuga::*; +#[allow(unused_variables)] fn main() -> Result<(), Box> { let rkf = RKF45::new(1e-4, 0.9, 1e-6, 1e-1, 100); let basic_ode_solver = BasicODESolver::new(rkf); diff --git a/src/lib.rs b/src/lib.rs index 4fcfc138..5b8221e6 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -147,14 +147,17 @@ //! //! ## Useful tips for features //! -//! * After `0.28.0`, `dataframe` feature is replaced by `nc` feature. -//! * If you want to use *QR* or *SVD* or *Cholesky Decomposition* then should use `O3` feature (there are no implementations for these decompositions in `default`) -//! * If you want to write your numerical results, then use `parquet` or `nc` features (corresponding to `parquet` or `netcdf` format. (It is much more effective than `csv` and `json`.) -//! * After `0.23.0`, there are two options - `fuga`, `prelude`. Choose proper option for your computations. -//! * To plot, use `parquet` or `nc` feature to export data as parquet or netcdf format and use python to draw plot. -//! * `plot` feature has limited plot abilities. -//! * To read parquet file in python, use `pandas` & `pyarrow` libraries. -//! * There is a template of python code for netcdf. - [Socialst](https://github.com/Axect/Socialst/blob/master/Templates/PyPlot_Template/nc_plot.py) +//! * If you want to use _QR_, _SVD_, or _Cholesky Decomposition_, you should use the `O3` feature. These decompositions are not implemented in the `default` feature. +//! +//! * If you want to save your numerical results, consider using the `parquet` or `nc` features, which correspond to the `parquet` and `netcdf` file formats, respectively. These formats are much more efficient than `csv` and `json`. +//! +//! * For plotting, it is recommended to use the `plot` feature. However, if you require more customization, you can use the `parquet` or `nc` feature to export your data in the parquet or netcdf format and then use Python to create the plots. +//! +//! * To read parquet files in Python, you can use the `pandas` and `pyarrow` libraries. +//! +//! * A template for Python code that works with netcdf files can be found in the [Socialst](https://github.com/Axect/Socialst/blob/master/Templates/PyPlot_Template/nc_plot.py) repository. + +//! #![cfg_attr(docsrs, feature(doc_auto_cfg))] #[cfg(feature = "O3")] diff --git a/src/numerical/eigen.rs b/src/numerical/eigen.rs index ad9c379c..688764eb 100644 --- a/src/numerical/eigen.rs +++ b/src/numerical/eigen.rs @@ -2,13 +2,9 @@ //! //! * Reference : Press, William H., and William T. Vetterling. *Numerical Recipes.* Cambridge: Cambridge Univ. Press, 2007. -#[cfg(feature = "O3")] -use lapack::{dorgtr, dsteqr, dsytrd}; - pub use self::EigenMethod::*; use crate::structure::matrix::Matrix; use crate::util::non_macro::eye_shape; -use std::f64::EPSILON; #[derive(Debug, Copy, Clone)] pub enum EigenMethod { @@ -114,11 +110,11 @@ impl JacobiTemp { for ip in 0..n - 1 { for iq in ip + 1..n { let g = 100f64 * a[(ip, iq)].abs(); - if i > 4 && g <= EPSILON * d[ip].abs() && g <= EPSILON * d[iq].abs() { + if i > 4 && g <= f64::EPSILON * d[ip].abs() && g <= f64::EPSILON * d[iq].abs() { a[(ip, iq)] = 0f64; } else if a[(ip, iq)].abs() > tresh { h = d[iq] - d[ip]; - let t = if g <= EPSILON * h.abs() { + let t = if g <= f64::EPSILON * h.abs() { a[(ip, iq)] / h } else { let theta = 0.5 * h / a[(ip, iq)]; diff --git a/src/numerical/mod.rs b/src/numerical/mod.rs index 3a26b42c..e228e556 100644 --- a/src/numerical/mod.rs +++ b/src/numerical/mod.rs @@ -1,7 +1,5 @@ //! Differential equations & Numerical Analysis tools -//pub mod bdf; -//pub mod gauss_legendre; pub mod eigen; pub mod integral; pub mod interp; diff --git a/src/numerical/newton.rs b/src/numerical/newton.rs index 79ddb150..f05998bc 100644 --- a/src/numerical/newton.rs +++ b/src/numerical/newton.rs @@ -12,12 +12,12 @@ pub fn newton) -> Vec + Copy>(init_cond: Vec, f: F, rtol let mut x_next = init_cond; let mut x = x_next.clone(); update(&mut x_next, f); - let mut err = (x_next.sub_vec(&x)).norm(Norm::L2); + let mut err = x_next.sub_vec(&x).norm(Norm::L2); while err >= rtol { x = x_next.clone(); update(&mut x_next, f); - err = (x_next.sub_vec(&x)).norm(Norm::L2); + err = x_next.sub_vec(&x).norm(Norm::L2); } x_next diff --git a/src/numerical/ode.rs b/src/numerical/ode.rs index 95748a66..fbfe846e 100644 --- a/src/numerical/ode.rs +++ b/src/numerical/ode.rs @@ -2,6 +2,28 @@ //! //! This module provides traits and structs for solving ordinary differential equations (ODEs). //! +//! ## Overview +//! +//! - `ODEProblem`: Trait for defining an ODE problem. +//! - `ODEIntegrator`: Trait for ODE integrators. +//! - `ODESolver`: Trait for ODE solvers. +//! +//! ## Available integrators +//! +//! - **Explicit** +//! - Runge-Kutta 4th order (RK4) +//! - Runge-Kutta-Fehlberg 4/5th order (RKF45) +//! - Dormand-Prince 4/5th order (DP45) +//! - Tsitouras 4/5th order (TSIT45) +//! - **Implicit** +//! - Gauss-Legendre 4th order (GL4) +//! +//! ## Available solvers +//! +//! - `BasicODESolver`: A basic ODE solver using a specified integrator. +//! +//! You can implement your own ODE solver by implementing the `ODESolver` trait. +//! //! ## Example //! //! ```rust @@ -89,8 +111,6 @@ pub trait ODEIntegrator { /// Enum for ODE errors. - - #[derive(Debug, Clone, Copy, Error)] pub enum ODEError { #[error("constraint violation")] @@ -213,13 +233,102 @@ impl ODEIntegrator for RK4 { } // ┌─────────────────────────────────────────────────────────┐ -// Runge-Kutta-Fehlberg +// Embedded Runge-Kutta // └─────────────────────────────────────────────────────────┘ +/// Trait for Butcher tableau +/// +/// ```text +/// C | A +/// - - - +/// | BH (higher order) +/// | BL (lower order) +/// ``` +/// +/// # References +/// +/// - J. R. Dormand and P. J. Prince, _A family of embedded Runge-Kutta formulae_, J. Comp. Appl. Math., 6(1), 19-26, 1980. +/// - Wikipedia: [List of Runge-Kutta methods](https://en.wikipedia.org/wiki/List_of_Runge%E2%80%93Kutta_methods) +pub trait ButcherTableau { + const C: &'static [f64]; + const A: &'static [&'static [f64]]; + const BH: &'static [f64]; + const BL: &'static [f64]; + + fn tol(&self) -> f64; + fn safety_factor(&self) -> f64; + fn max_step_size(&self) -> f64; + fn min_step_size(&self) -> f64; + fn max_step_iter(&self) -> usize; +} + +impl ODEIntegrator for BT { + fn step(&self, problem: &P, t: f64, y: &mut [f64], dt: f64) -> Result { + let n = y.len(); + let mut iter_count = 0usize; + let mut dt = dt; + let n_k = Self::C.len(); + + loop { + let mut k_vec = vec![vec![0.0; n]; n_k]; + let mut y_temp = y.to_vec(); + + for i in 0 .. n_k { + for i in 0 .. n { + let mut s = 0.0; + for j in 0 .. i { + s += Self::A[i][j] * k_vec[j][i]; + } + y_temp[i] = y[i] + dt * s; + } + problem.rhs(t + dt * Self::C[i], &y_temp, &mut k_vec[i])?; + } + + let mut error = 0f64; + for i in 0 .. n { + let mut s = 0.0; + for j in 0 .. n_k { + s += (Self::BH[j] - Self::BL[j]) * k_vec[j][i]; + } + error = error.max(dt * s.abs()) + } + + if error < self.tol() { + for i in 0 .. n { + let mut s = 0.0; + for j in 0 .. n_k { + s += Self::BH[j] * k_vec[j][i]; + } + y[i] += dt * s; + } + let new_dt = self.safety_factor() * dt * (self.tol() / error).powf(0.25); + let new_dt = new_dt.max(self.min_step_size()).min(self.max_step_size()); + return Ok(new_dt); + } else { + iter_count += 1; + if iter_count >= self.max_step_iter() { + return Err(ReachedMaxStepIter); + } + let new_dt = self.safety_factor() * dt * (self.tol() / error).powf(0.2); + let new_dt = new_dt.max(self.min_step_size()).min(self.max_step_size()); + dt = new_dt; + } + } + } +} + /// Runge-Kutta-Fehlberg 4/5th order integrator. /// /// This integrator uses the Runge-Kutta-Fehlberg method, which is an adaptive step size integrator. /// It calculates six intermediate values (k1, k2, k3, k4, k5, k6) to estimate the next step solution and the error. /// The step size is automatically adjusted based on the estimated error to maintain the desired tolerance. +/// +/// # Member variables +/// +/// - `tol`: The tolerance for the estimated error. +/// - `safety_factor`: The safety factor for the step size adjustment. +/// - `min_step_size`: The minimum step size. +/// - `max_step_size`: The maximum step size. +/// - `max_step_iter`: The maximum number of iterations per step. #[derive(Debug, Clone, Copy)] #[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))] pub struct RKF45 { @@ -252,117 +361,174 @@ impl RKF45 { max_step_iter, } } +} - pub fn set_tol(&mut self, tol: f64) { - self.tol = tol; - } - - pub fn set_safety_factor(&mut self, safety_factor: f64) { - self.safety_factor = safety_factor; - } - - pub fn set_min_step_size(&mut self, min_step_size: f64) { - self.min_step_size = min_step_size; - } - - pub fn set_max_step_size(&mut self, max_step_size: f64) { - self.max_step_size = max_step_size; - } +impl ButcherTableau for RKF45 { + const C: &'static [f64] = &[0.0, 1.0 / 4.0, 3.0 / 8.0, 12.0 / 13.0, 1.0, 1.0 / 2.0]; + const A: &'static [&'static [f64]] = &[ + &[0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + &[0.25, 0.0, 0.0, 0.0, 0.0, 0.0], + &[3.0 / 32.0, 9.0 / 32.0, 0.0, 0.0, 0.0, 0.0], + &[1932.0 / 2197.0, -7200.0 / 2197.0, 7296.0 / 2197.0, 0.0, 0.0, 0.0], + &[439.0 / 216.0, -8.0, 3680.0 / 513.0, -845.0 / 4104.0, 0.0, 0.0], + &[-8.0 / 27.0, 2.0, -3544.0 / 2565.0, 1859.0 / 4104.0, -11.0 / 40.0, 0.0], + ]; + const BH: &'static [f64] = &[16.0 / 135.0, 0.0, 6656.0 / 12825.0, 28561.0 / 56430.0, -9.0 / 50.0, 2.0 / 55.0]; + const BL: &'static [f64] = &[25.0 / 216.0, 0.0, 1408.0 / 2565.0, 2197.0 / 4104.0, -1.0 / 5.0, 0.0]; + + fn tol(&self) -> f64 { self.tol } + fn safety_factor(&self) -> f64 { self.safety_factor } + fn min_step_size(&self) -> f64 { self.min_step_size } + fn max_step_size(&self) -> f64 { self.max_step_size } + fn max_step_iter(&self) -> usize { self.max_step_iter } +} - pub fn set_max_step_iter(&mut self, max_step_iter: usize) { - self.max_step_iter = max_step_iter; - } +/// Dormand-Prince 5(4) method +/// +/// This is an adaptive step size integrator based on a 5th order Runge-Kutta method with +/// 4th order embedded error estimation. +/// +/// # Member variables +/// +/// - `tol`: The tolerance for the estimated error. +/// - `safety_factor`: The safety factor for the step size adjustment. +/// - `min_step_size`: The minimum step size. +/// - `max_step_size`: The maximum step size. +/// - `max_step_iter`: The maximum number of iterations per step. +#[derive(Debug, Clone, Copy)] +#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))] +pub struct DP45 { + tol: f64, + safety_factor: f64, + min_step_size: f64, + max_step_size: f64, + max_step_iter: usize, +} - pub fn get_tol(&self) -> f64 { - self.tol +impl Default for DP45 { + fn default() -> Self { + Self { + tol: 1e-6, + safety_factor: 0.9, + min_step_size: 1e-6, + max_step_size: 1e-1, + max_step_iter: 100, + } } +} - pub fn get_safety_factor(&self) -> f64 { - self.safety_factor +impl DP45 { + pub fn new(tol: f64, safety_factor: f64, min_step_size: f64, max_step_size: f64, max_step_iter: usize) -> Self { + Self { + tol, + safety_factor, + min_step_size, + max_step_size, + max_step_iter, + } } +} - pub fn get_min_step_size(&self) -> f64 { - self.min_step_size - } +impl ButcherTableau for DP45 { + const C: &'static [f64] = &[0.0, 0.2, 0.3, 0.8, 8.0 / 9.0, 1.0, 1.0]; + const A: &'static [&'static [f64]] = &[ + &[0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + &[0.2, 0.0, 0.0, 0.0, 0.0, 0.0], + &[0.075, 0.225, 0.0, 0.0, 0.0, 0.0], + &[44.0 / 45.0, -56.0 / 15.0, 32.0 / 9.0, 0.0, 0.0, 0.0], + &[19372.0 / 6561.0, -25360.0 / 2187.0, 64448.0 / 6561.0, -212.0 / 729.0, 0.0, 0.0], + &[9017.0 / 3168.0, -355.0 / 33.0, 46732.0 / 5247.0, 49.0 / 176.0, -5103.0 / 18656.0, 0.0], + &[35.0 / 384.0, 0.0, 500.0 / 1113.0, 125.0 / 192.0, -2187.0 / 6784.0, 11.0 / 84.0], + ]; + const BH: &'static [f64] = &[35.0 / 384.0, 0.0, 500.0 / 1113.0, 125.0 / 192.0, -2187.0 / 6784.0, 11.0 / 84.0, 0.0]; + const BL: &'static [f64] = &[5179.0 / 57600.0, 0.0, 7571.0 / 16695.0, 393.0 / 640.0, -92097.0 / 339200.0, 187.0 / 2100.0, 1.0 / 40.0]; + + fn tol(&self) -> f64 { self.tol } + fn safety_factor(&self) -> f64 { self.safety_factor } + fn min_step_size(&self) -> f64 { self.min_step_size } + fn max_step_size(&self) -> f64 { self.max_step_size } + fn max_step_iter(&self) -> usize { self.max_step_iter } +} - pub fn get_max_step_size(&self) -> f64 { - self.max_step_size - } +/// Tsitouras 5(4) method +/// +/// This is an adaptive step size integrator based on a 5th order Runge-Kutta method with +/// 4th order embedded error estimation, using the coefficients from Tsitouras (2011). +/// +/// # Member variables +/// +/// - `tol`: The tolerance for the estimated error. +/// - `safety_factor`: The safety factor for the step size adjustment. +/// - `min_step_size`: The minimum step size. +/// - `max_step_size`: The maximum step size. +/// - `max_step_iter`: The maximum number of iterations per step. +/// +/// # References +/// +/// - Ch. Tsitouras, Comput. Math. Appl. 62 (2011) 770-780. +#[derive(Debug, Clone, Copy)] +#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))] +pub struct TSIT45 { + tol: f64, + safety_factor: f64, + min_step_size: f64, + max_step_size: f64, + max_step_iter: usize, +} - pub fn get_max_step_iter(&self) -> usize { - self.max_step_iter +impl Default for TSIT45 { + fn default() -> Self { + Self { + tol: 1e-6, + safety_factor: 0.9, + min_step_size: 1e-6, + max_step_size: 1e-1, + max_step_iter: 100, + } } } -impl ODEIntegrator for RKF45 { - fn step(&self, problem: &P, t: f64, y: &mut [f64], dt: f64) -> Result { - let mut iter_count = 0usize; - let mut dt = dt; - let n = y.len(); - - loop { - let mut k1 = vec![0.0; n]; - let mut k2 = vec![0.0; n]; - let mut k3 = vec![0.0; n]; - let mut k4 = vec![0.0; n]; - let mut k5 = vec![0.0; n]; - let mut k6 = vec![0.0; n]; - - problem.rhs(t, y, &mut k1)?; - let mut y_temp = y.to_vec(); - for i in 0..n { - y_temp[i] = y[i] + dt * (1.0 / 4.0) * k1[i]; - } - - problem.rhs(t + dt * (1.0 / 4.0), &y_temp, &mut k2)?; - for i in 0..n { - y_temp[i] = y[i] + dt * (3.0 / 32.0 * k1[i] + 9.0 / 32.0 * k2[i]); - } - - problem.rhs(t + dt * (3.0 / 8.0), &y_temp, &mut k3)?; - for i in 0..n { - y_temp[i] = y[i] + dt * (1932.0 / 2197.0 * k1[i] - 7200.0 / 2197.0 * k2[i] + 7296.0 / 2197.0 * k3[i]); - } - - problem.rhs(t + dt * (12.0 / 13.0), &y_temp, &mut k4)?; - for i in 0..n { - y_temp[i] = y[i] + dt * (439.0 / 216.0 * k1[i] - 8.0 * k2[i] + 3680.0 / 513.0 * k3[i] - 845.0 / 4104.0 * k4[i]); - } - - problem.rhs(t + dt, &y_temp, &mut k5)?; - for i in 0..n { - y_temp[i] = y[i] + dt * (-8.0 / 27.0 * k1[i] + 2.0 * k2[i] - 3544.0 / 2565.0 * k3[i] + 1859.0 / 4104.0 * k4[i] - 11.0 / 40.0 * k5[i]); - } - - problem.rhs(t + dt * 0.5, &y_temp, &mut k6)?; - - let mut error = 0.0; - for i in 0..n { - let y_err = dt * (1.0 / 360.0 * k1[i] - 128.0 / 4275.0 * k3[i] - 2197.0 / 75240.0 * k4[i] + 1.0 / 50.0 * k5[i] + 2.0 / 55.0 * k6[i]); - error += y_err * y_err; - } - - error = error.sqrt() / (n as f64).sqrt(); - - let new_dt = self.safety_factor * dt * (self.tol * dt / error).powf(0.25); - let new_dt = new_dt.max(self.min_step_size).min(self.max_step_size); - - if error <= self.tol { - for i in 0..n { - y[i] += dt * (16.0 / 135.0 * k1[i] + 6656.0 / 12825.0 * k3[i] + 28561.0 / 56430.0 * k4[i] - 9.0 / 50.0 * k5[i] + 2.0 / 55.0 * k6[i]); - } - return Ok(new_dt); - } else { - iter_count += 1; - if iter_count >= self.max_step_iter { - return Err(ReachedMaxStepIter); - } - dt = new_dt; - } +impl TSIT45 { + pub fn new(tol: f64, safety_factor: f64, min_step_size: f64, max_step_size: f64, max_step_iter: usize) -> Self { + Self { + tol, + safety_factor, + min_step_size, + max_step_size, + max_step_iter, } } } +impl ButcherTableau for TSIT45 { + const C: &'static [f64] = &[0.0, 0.161, 0.327, 0.9, 0.9800255409045097, 1.0, 1.0]; + const A: &'static [&'static [f64]] = &[ + &[0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + &[Self::C[1], 0.0, 0.0, 0.0, 0.0, 0.0], + &[Self::C[2] - 0.335480655492357, 0.335480655492357, 0.0, 0.0, 0.0, 0.0], + &[Self::C[3] - (-6.359448489975075 + 4.362295432869581), -6.359448489975075, 4.362295432869581, 0.0, 0.0, 0.0], + &[Self::C[4] - (-11.74888356406283 + 7.495539342889836 - 0.09249506636175525), -11.74888356406283, 7.495539342889836, -0.09249506636175525, 0.0, 0.0], + &[Self::C[5] - (-12.92096931784711 + 8.159367898576159 - 0.0715849732814010 - 0.02826905039406838), -12.92096931784711, 8.159367898576159, -0.0715849732814010, -0.02826905039406838, 0.0], + &[Self::BH[0], Self::BH[1], Self::BH[2], Self::BH[3], Self::BH[4], Self::BH[5]], + ]; + const BH: &'static [f64] = &[0.09646076681806523, 0.01, 0.4798896504144996, 1.379008574103742, -3.290069515436081, 2.324710524099774, 0.0]; + const BL: &'static [f64] = &[ + 0.001780011052226, + 0.000816434459657, + - 0.007880878010262, + 0.144711007173263, + - 0.582357165452555, + 0.458082105929187, + 1.0 / 66.0, + ]; + + fn tol(&self) -> f64 { self.tol } + fn safety_factor(&self) -> f64 { self.safety_factor } + fn min_step_size(&self) -> f64 { self.min_step_size } + fn max_step_size(&self) -> f64 { self.max_step_size } + fn max_step_iter(&self) -> usize { self.max_step_iter } +} + // ┌─────────────────────────────────────────────────────────┐ // Gauss-Legendre 4th order // └─────────────────────────────────────────────────────────┘ @@ -383,30 +549,36 @@ pub enum ImplicitSolver { /// This integrator uses the 4th order Gauss-Legendre Runge-Kutta method, which is an implicit integrator. /// It requires solving a system of nonlinear equations at each step, which is done using the specified implicit solver (e.g., fixed-point iteration). /// The Gauss-Legendre method has better stability properties compared to explicit methods, especially for stiff ODEs. +/// +/// # Member variables +/// +/// - `solver`: The implicit solver to use. +/// - `tol`: The tolerance for the implicit solver. +/// - `max_step_iter`: The maximum number of iterations for the implicit solver per step. #[derive(Debug, Clone, Copy)] #[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))] pub struct GL4 { solver: ImplicitSolver, - max_iter: usize, tol: f64, + max_step_iter: usize, } impl Default for GL4 { fn default() -> Self { GL4 { solver: ImplicitSolver::FixedPoint, - max_iter: 100, tol: 1e-6, + max_step_iter: 100, } } } impl GL4 { - pub fn new(solver: ImplicitSolver, max_iter: usize, tol: f64) -> Self { + pub fn new(solver: ImplicitSolver, tol: f64, max_step_iter: usize) -> Self { GL4 { solver, - max_iter, tol, + max_step_iter, } } } @@ -428,7 +600,7 @@ impl ODEIntegrator for GL4 { match self.solver { ImplicitSolver::FixedPoint => { // Fixed-point iteration - for _ in 0..self.max_iter { + for _ in 0..self.max_step_iter { problem.rhs(t + c * dt, &y1, &mut k1)?; problem.rhs(t + d * dt, &y2, &mut k2)?; diff --git a/src/numerical/optimize.rs b/src/numerical/optimize.rs index 38fd6139..b894e28b 100644 --- a/src/numerical/optimize.rs +++ b/src/numerical/optimize.rs @@ -287,7 +287,7 @@ where let mut chi2 = ((&y - &y_hat).t() * (&y - &y_hat))[(0, 0)]; let mut nu = 2f64; let lambda_0 = *self.hyperparams.get("lambda_init").unwrap_or(&1e-3); - let lambda_max = *self.hyperparams.get("lambda_max").unwrap_or(&std::f64::MAX.sqrt()); + let lambda_max = *self.hyperparams.get("lambda_max").unwrap_or(&f64::MAX.sqrt()); let mut lambda = lambda_0 * max(jtj.diag()); diff --git a/src/numerical/root.rs b/src/numerical/root.rs index f205ffd6..8769f1c9 100644 --- a/src/numerical/root.rs +++ b/src/numerical/root.rs @@ -435,7 +435,7 @@ pub fn bisection AD>( tol: f64, ) -> Result { let (a, b) = interval; - let mut root_finder = RootFinder::new(RootState::I(a, b), RootFind::Bisection, f)?; + let mut root_finder = RootFinder::new(I(a, b), RootFind::Bisection, f)?; root_finder.set_tol(tol); root_finder.set_times(times); root_finder.find_root() @@ -465,7 +465,7 @@ pub fn false_position AD>( tol: f64, ) -> Result { let (a, b) = interval; - let mut root_finder = RootFinder::new(RootState::I(a, b), RootFind::FalsePosition, f)?; + let mut root_finder = RootFinder::new(I(a, b), RootFind::FalsePosition, f)?; root_finder.set_tol(tol); root_finder.set_times(times); root_finder.find_root() @@ -494,7 +494,7 @@ pub fn newton AD>( times: usize, tol: f64, ) -> Result { - let mut root_finder = RootFinder::new(RootState::P(initial_guess), RootFind::Newton, f)?; + let mut root_finder = RootFinder::new(P(initial_guess), RootFind::Newton, f)?; root_finder.set_tol(tol); root_finder.set_times(times); root_finder.find_root() @@ -524,7 +524,7 @@ pub fn secant AD>( tol: f64, ) -> Result { let (a, b) = initial_guess; - let mut root_finder = RootFinder::new(RootState::I(a, b), RootFind::Secant, f)?; + let mut root_finder = RootFinder::new(I(a, b), RootFind::Secant, f)?; root_finder.set_tol(tol); root_finder.set_times(times); root_finder.find_root() diff --git a/src/numerical/spline.rs b/src/numerical/spline.rs index 16bbb4ca..c32f9fc6 100644 --- a/src/numerical/spline.rs +++ b/src/numerical/spline.rs @@ -37,7 +37,7 @@ //! ```rust //! use peroxide::fuga::*; //! -//! fn main() -> Result<(), Box> { +//! fn main() -> Result<(), Box> { //! let x = seq(0, 10, 1); //! let y = x.fmap(|t| t.sin()); //! @@ -91,7 +91,7 @@ //! ```rust //! use peroxide::fuga::*; //! -//! fn main() -> Result<(), Box> { +//! fn main() -> Result<(), Box> { //! let x = seq(0, 10, 1); //! let y = x.fmap(|t| t.sin()); //! @@ -141,7 +141,7 @@ //! use peroxide::fuga::*; //! use std::f64::consts::PI; //! -//! fn main() -> Result<(), Box> { +//! fn main() -> Result<(), Box> { //! let x = seq(0, 10, 1); //! let y = x.fmap(|t| t.sin()); //! @@ -293,7 +293,7 @@ pub trait Spline { /// extern crate peroxide; /// use peroxide::fuga::*; /// -/// fn main() -> Result<(), Box> { +/// fn main() -> Result<(), Box> { /// let x = c!(0.9, 1.3, 1.9, 2.1); /// let y = c!(1.3, 1.5, 1.85, 2.1); /// @@ -347,7 +347,7 @@ pub fn cubic_hermite_spline( /// extern crate peroxide; /// use peroxide::fuga::*; /// -/// fn main() -> Result<(), Box> { +/// fn main() -> Result<(), Box> { /// let x = c!(0.9, 1.3, 1.9, 2.1); /// let y = c!(1.3, 1.5, 1.85, 2.1); /// @@ -402,7 +402,7 @@ impl CubicSpline { /// extern crate peroxide; /// use peroxide::fuga::*; /// - /// fn main() -> Result<(), Box> { + /// fn main() -> Result<(), Box> { /// let x = c!(0.9, 1.3, 1.9, 2.1); /// let y = c!(1.3, 1.5, 1.85, 2.1); /// diff --git a/src/special/lanczos.rs b/src/special/lanczos.rs index b16dcabd..eeded0dd 100644 --- a/src/special/lanczos.rs +++ b/src/special/lanczos.rs @@ -49,7 +49,7 @@ pub fn gamma_approx(z: f64) -> f64 { /// Lanczos Approximation Coefficient pub fn tlg1(g: f64, n: usize) -> Vec { - (lanczos_coeff(g, n - 1).ox() * g.exp() / (2f64 * std::f64::consts::PI).sqrt()).red() + (lanczos_coeff(g, n - 1).ox() * g.exp() / (2f64 * PI).sqrt()).red() } fn lanczos_coeff(g: f64, n: usize) -> Vec { diff --git a/src/statistics/dist.rs b/src/statistics/dist.rs index 59dd3fa7..d0857493 100644 --- a/src/statistics/dist.rs +++ b/src/statistics/dist.rs @@ -301,7 +301,7 @@ impl WeightedUniform { /// extern crate peroxide; /// use peroxide::fuga::*; /// - /// fn main() -> Result<(), Box> { + /// fn main() -> Result<(), Box> { /// let weights = vec![1f64, 3f64, 0f64, 2f64]; /// let intervals = vec![0f64, 1f64, 2f64, 4f64, 5f64]; /// let w = WeightedUniform::new(weights, intervals)?; @@ -357,7 +357,7 @@ impl WeightedUniform { /// extern crate peroxide; /// use peroxide::fuga::*; /// - /// fn main() -> Result<(), Box> { + /// fn main() -> Result<(), Box> { /// let w = WeightedUniform::from_max_pool_1d(f, (-2f64, 3f64), 10, 1e-3)?; /// w.weights().print(); /// @@ -366,7 +366,7 @@ impl WeightedUniform { /// /// fn f(x: f64) -> f64 { /// if x.abs() < 1f64 { - /// (1f64 - x.abs()) + /// 1f64 - x.abs() /// } else { /// 0f64 /// } diff --git a/src/statistics/rand.rs b/src/statistics/rand.rs index 2db9557c..b5a3ee0c 100644 --- a/src/statistics/rand.rs +++ b/src/statistics/rand.rs @@ -23,7 +23,6 @@ extern crate rand; use self::rand::distributions::uniform::SampleUniform; use self::rand::prelude::*; -use std::u32; #[allow(unused_imports)] use crate::structure::matrix::*; @@ -69,7 +68,7 @@ pub fn stdrng_from_seed(seed: u64) -> StdRng { /// use peroxide::fuga::*; /// /// let mut rng = thread_rng(); -/// println!("{}", rand_num(&mut rng, 1, 7)); // Roll a dice +/// println!("{}", rand_num(&mut rng, 1, 7)); // Roll a die /// println!("{}", rand_num(&mut rng, 0f64, 1f64)); // Uniform [0,1) /// ``` pub fn rand_num(rng: &mut ThreadRng, start: T, end: T) -> T @@ -453,7 +452,7 @@ pub fn ziggurat(rng: &mut ThreadRng, sigma: f64) -> f64 { /// ``` /// use peroxide::fuga::*; /// -/// fn main() -> Result<(), Box> { +/// fn main() -> Result<(), Box> { /// let f = |x: f64| { /// if (0f64..=2f64).contains(&x) { /// -(x - 1f64).powi(2) + 1f64 @@ -514,7 +513,7 @@ where F: Fn(f64) -> f64 + Copy { /// ``` /// use peroxide::fuga::*; /// -/// fn main() -> Result<(), Box> { +/// fn main() -> Result<(), Box> { /// let mut rng = smallrng_from_seed(42); /// let f = |x: f64| { /// if (0f64..=2f64).contains(&x) { diff --git a/src/structure/ad.rs b/src/structure/ad.rs index fcd25985..e850c716 100644 --- a/src/structure/ad.rs +++ b/src/structure/ad.rs @@ -1209,7 +1209,7 @@ impl ADFn { pub fn grad(&self) -> Self { assert!(self.grad_level < 2, "Higher order AD is not allowed"); ADFn { - f: (self.f).clone(), + f: self.f.clone(), grad_level: self.grad_level + 1, } } @@ -1311,18 +1311,18 @@ impl FPVector for Vec { self.iter().fold(init.into(), |x, &y| f(x,y)) } - fn filter(&self, f: F) -> Self - where - F: Fn(Self::Scalar) -> bool { - self.iter().filter(|&x| f(*x)).cloned().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() } + fn filter(&self, f: F) -> Self + where + F: Fn(Self::Scalar) -> bool { + self.iter().filter(|&x| f(*x)).cloned().collect() + } + fn take(&self, n: usize) -> Self { self.iter().take(n).cloned().collect() } diff --git a/src/structure/dataframe.rs b/src/structure/dataframe.rs index 2d8feec8..51b03107 100644 --- a/src/structure/dataframe.rs +++ b/src/structure/dataframe.rs @@ -950,19 +950,19 @@ fn dtype_to_arrow(dt: DType) -> DataType { #[cfg(feature="parquet")] fn arrow_to_dtype(dt: DataType) -> DType { match dt { - DataType::Boolean => DType::Bool, - DataType::Int8 => DType::I8, - DataType::Int16 => DType::I16, - DataType::Int32 => DType::I32, - DataType::Int64 => DType::I64, - DataType::UInt8 => DType::U8, - DataType::UInt16 => DType::U16, - DataType::UInt32 => DType::U32, - DataType::UInt64 => DType::U64, + DataType::Boolean => Bool, + DataType::Int8 => I8, + DataType::Int16 => I16, + DataType::Int32 => I32, + DataType::Int64 => I64, + DataType::UInt8 => U8, + DataType::UInt16 => U16, + DataType::UInt32 => U32, + DataType::UInt64 => U64, // DataType::Float16 => DType::F16, - DataType::Float32 => DType::F32, - DataType::Float64 => DType::F64, - DataType::Utf8 => DType::Str, + DataType::Float32 => F32, + DataType::Float64 => F64, + DataType::Utf8 => Str, _ => unimplemented!() } } @@ -1192,7 +1192,7 @@ impl Vector for Series { /// assert_eq!(c, Series::new(vec![4,4,4])); /// } /// ``` - fn add_vec<'a, 'b>(&'a self, rhs: &'b Self) -> Self { + fn add_vec(&self, rhs: &Self) -> Self { assert_eq!(self.dtype, rhs.dtype, "DTypes are not same (add_vec)"); dtype_match!( N; @@ -1218,7 +1218,7 @@ impl Vector for Series { /// assert_eq!(c, Series::new(vec![3,3,3])); /// } /// ``` - fn sub_vec<'a, 'b>(&'a self, rhs: &'b Self) -> Self { + fn sub_vec(&self, rhs: &Self) -> Self { assert_eq!(self.dtype, rhs.dtype, "DTypes are not same (add_vec)"); dtype_match!( N; diff --git a/src/structure/matrix.rs b/src/structure/matrix.rs index b0cbb4f4..797792ca 100644 --- a/src/structure/matrix.rs +++ b/src/structure/matrix.rs @@ -40,7 +40,6 @@ //! * Type: `ml_matrix(&str)` //! //! ```rust -//! extern crate peroxide; //! use peroxide::fuga::*; //! //! fn main() { @@ -58,7 +57,6 @@ //! * Type: `py_matrix(Vec>) where T: std::convert::Into + Copy` //! //! ```rust -//! extern crate peroxide; //! use peroxide::fuga::*; //! //! fn main() { @@ -151,7 +149,6 @@ //! * `write_with_header(&self, file_path, header: Vec<&str>)`: Write with header //! //! ```rust -//! extern crate peroxide; //! use peroxide::fuga::*; //! //! fn main() { @@ -174,7 +171,6 @@ //! * Description: `read(file_path, is_header, delimiter)` //! //! ```no_run -//! extern crate peroxide; //! use peroxide::fuga::*; //! //! fn main() { @@ -233,27 +229,28 @@ //! * `cbind`: Concatenate two matrices by column direction. //! * `rbind`: Concatenate two matrices by row direction. //! -//! ```rust -//! extern crate peroxide; -//! use peroxide::fuga::*; +//! ```rust +//! use peroxide::fuga::*; //! -//! fn main() { -//! let a = ml_matrix("1 2;3 4"); -//! let b = ml_matrix("5 6;7 8"); +//! fn main() -> Result<(), Box> { +//! let a = ml_matrix("1 2;3 4"); +//! let b = ml_matrix("5 6;7 8"); //! -//! cbind(a.clone(), b.clone()).print(); -//! // c[0] c[1] c[2] c[3] -//! // r[0] 1 2 5 7 -//! // r[1] 3 4 6 8 +//! cbind(a.clone(), b.clone())?.print(); +//! // c[0] c[1] c[2] c[3] +//! // r[0] 1 2 5 7 +//! // r[1] 3 4 6 8 //! -//! rbind(a, b).print(); -//! // c[0] c[1] -//! // r[0] 1 2 -//! // r[1] 3 4 -//! // r[2] 5 6 -//! // r[3] 7 8 -//! } -//! ``` +//! rbind(a, b)?.print(); +//! // c[0] c[1] +//! // r[0] 1 2 +//! // r[1] 3 4 +//! // r[2] 5 6 +//! // r[3] 7 8 +//! +//! Ok(()) +//! } +//! ``` //! //! ## Matrix operations //! @@ -287,8 +284,8 @@ //! * `clone` is too annoying - We can use reference operations! //! //! ```rust -//! extern crate peroxide; //! use peroxide::fuga::*; +//! //! fn main() { //! let a = ml_matrix("1 2;3 4"); //! let b = ml_matrix("5 6;7 8"); @@ -342,9 +339,9 @@ //! * `zeros(usize, usize)`: Construct matrix which elements are all zero //! * `eye(usize)`: Identity matrix //! * `rand(usize, usize)`: Construct random uniform matrix (from 0 to 1) +//! * `rand_with_rng(usize, usize, &mut Rng)`: Construct random matrix with user-defined RNG //! //! ```rust -//! extern crate peroxide; //! use peroxide::fuga::*; //! //! fn main() { @@ -393,7 +390,6 @@ //! * The structure of `PQLU` is as follows: //! //! ```rust -//! extern crate peroxide; //! use peroxide::fuga::*; //! //! #[derive(Debug, Clone)] @@ -470,7 +466,6 @@ //! can get better performance (via memoization) //! //! ```rust -//! extern crate peroxide; //! use peroxide::fuga::*; //! //! fn main() { @@ -496,7 +491,6 @@ //! * Example //! //! ```ignore -//! extern crate peroxide; //! use peroxide::fuga::*; //! //! let a = ml_matrix(" @@ -529,7 +523,6 @@ //! * Return `SVD` structure //! //! ```no_run -//! extern crate peroxide; //! use peroxide::fuga::*; //! //! pub struct SVD { @@ -542,7 +535,6 @@ //! * Example //! //! ``` -//! extern crate peroxide; //! use peroxide::fuga::*; //! //! fn main() { @@ -618,9 +610,6 @@ extern crate lapack; use blas::{daxpy, dgemm, dgemv}; #[cfg(feature = "O3")] use lapack::{dgecon, dgeqrf, dgetrf, dgetri, dgetrs, dorgqr, dgesvd, dpotrf}; -#[cfg(feature = "O3")] -use std::f64::NAN; - #[cfg(feature = "serde")] use serde::{Deserialize, Serialize}; @@ -639,7 +628,6 @@ use crate::util::{ }; use crate::structure::dataframe::{Series, TypedVector}; use std::cmp::{max, min}; -use std::convert; pub use std::error::Error; use std::fmt; use std::ops::{Add, Div, Index, IndexMut, Mul, Neg, Sub}; @@ -655,7 +643,6 @@ pub type Perms = Vec<(usize, usize)>; /// /// # Examples /// ``` -/// extern crate peroxide; /// use peroxide::fuga::*; /// /// let a = matrix(vec![1,2,3,4], 2, 2, Row); @@ -687,7 +674,6 @@ impl fmt::Display for Shape { /// # Examples /// /// ``` -/// extern crate peroxide; /// use peroxide::fuga::*; /// /// let a = Matrix { @@ -725,7 +711,7 @@ pub struct Matrix { /// ``` pub fn matrix(v: Vec, r: usize, c: usize, shape: Shape) -> Matrix where - T: convert::Into, + T: Into, { Matrix { data: v.into_iter().map(|t| t.into()).collect::>(), @@ -738,7 +724,7 @@ where /// R-like matrix constructor (Explicit ver.) pub fn r_matrix(v: Vec, r: usize, c: usize, shape: Shape) -> Matrix where - T: convert::Into, + T: Into, { matrix(v, r, c, shape) } @@ -759,7 +745,7 @@ where /// ``` pub fn py_matrix(v: Vec>) -> Matrix where - T: convert::Into + Copy, + T: Into + Copy, { let r = v.len(); let c = v[0].len(); @@ -846,14 +832,13 @@ impl Matrix { /// /// # Examples /// ``` - /// extern crate peroxide; /// use peroxide::fuga::*; /// /// let a = matrix(vec![1,2,3,4], 2, 2, Col); /// let b = a.as_slice(); /// assert_eq!(b, &[1f64,2f64,3f64,4f64]); /// ``` - pub fn as_slice<'a>(&'a self) -> &'a [f64] { + pub fn as_slice(&self) -> &[f64] { &self.data[..] } @@ -861,7 +846,6 @@ impl Matrix { /// /// # Examples /// ``` - /// extern crate peroxide; /// use peroxide::fuga::*; /// /// let mut a = matrix(vec![1,2,3,4], 2, 2, Col); @@ -870,7 +854,7 @@ impl Matrix { /// assert_eq!(b, &[5f64, 2f64, 3f64, 4f64]); /// assert_eq!(a, matrix(vec![5,2,3,4], 2, 2, Col)); /// ``` - pub fn as_mut_slice<'a>(&'a mut self) -> &'a mut [f64] { + pub fn as_mut_slice(&mut self) -> &mut [f64] { &mut self.data[..] } @@ -880,7 +864,6 @@ impl Matrix { /// /// # Examples /// ``` - /// extern crate peroxide; /// use peroxide::fuga::*; /// /// let a = matrix(vec![1,2,3,4],2,2,Row); @@ -922,7 +905,6 @@ impl Matrix { /// /// # Examples /// ``` - /// extern crate peroxide; /// use peroxide::fuga::*; /// /// let a = matrix(vec![1,2,3,4],2,2,Row); @@ -961,7 +943,6 @@ impl Matrix { /// /// # Examples /// ``` - /// extern crate peroxide; /// use peroxide::fuga::*; /// /// let a = matrix(vec![1,2,3,4],2,2,Row); @@ -1118,7 +1099,6 @@ impl Matrix { /// /// # Examples /// ``` - /// extern crate peroxide; /// use peroxide::fuga::*; /// /// let a = matrix(vec![1,2,3,4], 2, 2, Row); @@ -1156,9 +1136,11 @@ impl Matrix { /// extern crate peroxide; /// use peroxide::fuga::*; /// - /// fn main() { + /// fn main() -> Result<(), Box> { /// let a = matrix(c!(1,2,3,3,2,1), 3, 2, Col); - /// a.write("example_data/test.csv"); + /// a.write("example_data/test.csv")?; + /// + /// Ok(()) /// } /// ``` #[cfg(feature="csv")] @@ -1185,9 +1167,11 @@ impl Matrix { /// extern crate peroxide; /// use peroxide::fuga::*; /// - /// fn main() { + /// fn main() -> Result<(), Box> { /// let a = matrix(c!(1,2,3,3,2,1), 3, 2, Col); - /// a.write_round("example_data/test.csv", 0); + /// a.write_round("example_data/test.csv", 0)?; + /// + /// Ok(()) /// } /// ``` #[cfg(feature="csv")] @@ -1260,9 +1244,9 @@ impl Matrix { /// use peroxide::fuga::*; /// use std::process; /// - /// fn main() { + /// fn main() -> Result<(), Box> { /// let a = matrix(c!(1,2,3,3,2,1), 3, 2, Col); - /// a.write_round("example_data/test.csv", 0); + /// a.write_round("example_data/test.csv", 0)?; /// /// let b = Matrix::read("example_data/test.csv", false, ','); // header = false, delimiter = ',' /// match b { @@ -1272,6 +1256,8 @@ impl Matrix { /// process::exit(1); /// } /// } + /// + /// Ok(()) /// } /// ``` #[cfg(feature="csv")] @@ -1366,7 +1352,7 @@ impl Matrix { } pub fn to_diag(&self) -> Matrix { - assert!(self.row == self.col, "Should be square matrix"); + assert_eq!(self.row, self.col, "Should be square matrix"); let mut result = matrix(vec![0f64; self.row * self.col], self.row, self.col, Row); let diag = self.diag(); for i in 0..self.row { @@ -1580,12 +1566,12 @@ impl Normed for Matrix { for i in 0..self.row { s_row += self[(i, j)].powi(p as i32); } - s += s_row.powf(q as f64 / (p as f64)); + s += s_row.powf(q / p); } - s.powf(1f64 / (q as f64)) + s.powf(1f64 / q) } Norm::L1 => { - let mut m = std::f64::MIN; + let mut m = f64::MIN; match self.shape { Row => self.change_shape().norm(Norm::L1), Col => { @@ -1600,7 +1586,7 @@ impl Normed for Matrix { } } Norm::LInf => { - let mut m = std::f64::MIN; + let mut m = f64::MIN; match self.shape { Col => self.change_shape().norm(Norm::LInf), Row => { @@ -1684,16 +1670,16 @@ impl MatrixProduct for Matrix { for j in 1..c1 { let n = self[(0, j)] * other; - result = cbind(result, n); + result = 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 = cbind(m, n); + m = cbind(m, n).unwrap(); } - result = rbind(result, m); + result = rbind(result, m).unwrap(); } result } @@ -1814,7 +1800,7 @@ impl<'a, 'b> Add<&'b Matrix> for &'a Matrix { /// ``` impl Add for Matrix where - T: convert::Into + Copy, + T: Into + Copy, { type Output = Self; fn add(self, other: T) -> Self { @@ -1837,7 +1823,7 @@ where /// Element-wise addition between &Matrix & f64 impl<'a, T> Add for &'a Matrix where - T: convert::Into + Copy, + T: Into + Copy, { type Output = Matrix; fn add(self, other: T) -> Self::Output { @@ -1927,7 +1913,7 @@ impl<'a> Add<&'a Matrix> for i32 { /// /// fn main() { /// let a = matrix!(1;4;1, 2, 2, Row); -/// assert_eq!(1 as usize + a, matrix!(2;5;1, 2, 2, Row)); +/// assert_eq!(1usize + a, matrix!(2;5;1, 2, 2, Row)); /// } /// ``` impl Add for usize { @@ -2074,7 +2060,7 @@ impl<'a, 'b> Sub<&'b Matrix> for &'a Matrix { /// Subtraction between Matrix & f64 impl Sub for Matrix where - T: convert::Into + Copy, + T: Into + Copy, { type Output = Self; @@ -2098,7 +2084,7 @@ where /// Subtraction between &Matrix & f64 impl<'a, T> Sub for &'a Matrix where - T: convert::Into + Copy, + T: Into + Copy, { type Output = Matrix; @@ -2736,7 +2722,7 @@ impl FPMatrix for Matrix { fn reduce(&self, init: T, f: F) -> f64 where F: Fn(f64, f64) -> f64, - T: convert::Into, + T: Into, { self.data.iter().fold(init.into(), |x, y| f(x, *y)) } @@ -3360,7 +3346,7 @@ impl LinearAlgebra for Matrix { () => { let opt_dgrf = lapack_dgetrf(self); match opt_dgrf { - None => NAN, + None => f64::NAN, Some(dgrf) => match dgrf.status { LAPACK_STATUS::Singular => 0f64, LAPACK_STATUS::NonSingular => { @@ -3654,7 +3640,7 @@ impl MutMatrix for Matrix { unsafe fn col_mut(&mut self, idx: usize) -> Vec<*mut f64> { assert!(idx < self.col, "Index out of range"); match self.shape { - Shape::Col => { + Col => { let mut v: Vec<*mut f64> = vec![&mut 0f64; self.row]; let start_idx = idx * self.row; let p = self.mut_ptr(); @@ -3663,7 +3649,7 @@ impl MutMatrix for Matrix { } v } - Shape::Row => { + Row => { let mut v: Vec<*mut f64> = vec![&mut 0f64; self.row]; let p = self.mut_ptr(); for i in 0..v.len() { @@ -3677,7 +3663,7 @@ impl MutMatrix for Matrix { unsafe fn row_mut(&mut self, idx: usize) -> Vec<*mut f64> { assert!(idx < self.row, "Index out of range"); match self.shape { - Shape::Row => { + Row => { let mut v: Vec<*mut f64> = vec![&mut 0f64; self.col]; let start_idx = idx * self.col; let p = self.mut_ptr(); @@ -3686,7 +3672,7 @@ impl MutMatrix for Matrix { } v } - Shape::Col => { + Col => { let mut v: Vec<*mut f64> = vec![&mut 0f64; self.col]; let p = self.mut_ptr(); for i in 0..v.len() { @@ -3699,8 +3685,8 @@ impl MutMatrix for Matrix { 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)), + Col => swap_vec_ptr(&mut self.col_mut(idx1), &mut self.col_mut(idx2)), + Row => swap_vec_ptr(&mut self.row_mut(idx1), &mut self.row_mut(idx2)), } } @@ -3777,6 +3763,14 @@ impl PowOps for Matrix { } impl TrigOps for Matrix { + fn sin_cos(&self) -> (Self, Self) { + let (sin, cos) = self.data.iter().map(|x| x.sin_cos()).unzip(); + ( + matrix(sin, self.row, self.col, self.shape), + matrix(cos, self.row, self.col, self.shape), + ) + } + fn sin(&self) -> Self { self.fmap(|x| x.sin()) } @@ -3813,14 +3807,6 @@ impl TrigOps for Matrix { self.fmap(|x| x.atan()) } - fn sin_cos(&self) -> (Self, Self) { - let (sin, cos) = self.data.iter().map(|x| x.sin_cos()).unzip(); - ( - matrix(sin, self.row, self.col, self.shape), - matrix(cos, self.row, self.col, self.shape), - ) - } - fn asinh(&self) -> Self { self.fmap(|x| x.asinh()) } @@ -4318,7 +4304,7 @@ pub fn lapack_dgetrf(mat: &Matrix) -> Option { }; let mut info = 0i32; - let mut ipiv: Vec = vec![0i32; std::cmp::min(m, n)]; + let mut ipiv: Vec = vec![0i32; min(m, n)]; unsafe { dgetrf(m_i32, n_i32, &mut a, m_i32, &mut ipiv, &mut info); @@ -4515,7 +4501,7 @@ pub fn lapack_dgesvd(mat: &Matrix) -> Option { if info == 0 { Some(DGESVD { - s: s, + s, u: matrix(u, m as usize, m as usize, Col), vt: matrix(vt, n as usize, n as usize, Col), status: SVD_STATUS::Success, @@ -4524,7 +4510,7 @@ pub fn lapack_dgesvd(mat: &Matrix) -> Option { None } else { Some(DGESVD { - s: s, + s, u: matrix(u, m as usize, m as usize, Col), vt: matrix(vt, n as usize, n as usize, Col), status: SVD_STATUS::Diverge(info), diff --git a/src/structure/polynomial.rs b/src/structure/polynomial.rs index 0487f181..07232363 100644 --- a/src/structure/polynomial.rs +++ b/src/structure/polynomial.rs @@ -9,7 +9,6 @@ use serde::{Deserialize, Serialize}; use peroxide_num::PowOps; use crate::traits::fp::FPVector; use std::cmp::{max, min}; -use std::convert; use std::fmt; use std::ops::{Add, Div, Mul, Neg, Sub}; @@ -180,7 +179,7 @@ impl Polynomial { /// ``` pub fn eval(&self, x: T) -> f64 where - T: convert::Into + Copy, + T: Into + Copy, { let x = x.into(); let l = self.coef.len() - 1; @@ -212,7 +211,7 @@ impl Polynomial { /// ``` pub fn translate_x(&self, x: X) -> Self where - X: convert::Into + Copy, + X: Into + Copy, { let d = Self::new(vec![1f64, x.into()]); @@ -299,7 +298,7 @@ impl Add for Polynomial { impl Add for Polynomial where - T: convert::Into + Copy, + T: Into + Copy, { type Output = Self; fn add(self, other: T) -> Self { @@ -318,7 +317,7 @@ impl Sub for Polynomial { impl Sub for Polynomial where - T: convert::Into + Copy, + T: Into + Copy, { type Output = Self; fn sub(self, other: T) -> Self { @@ -330,7 +329,7 @@ where impl Mul for Polynomial where - T: convert::Into + Copy, + T: Into + Copy, { type Output = Self; fn mul(self, other: T) -> Self { @@ -372,7 +371,7 @@ impl Mul for Polynomial { impl Div for Polynomial where - T: convert::Into + Copy, + T: Into + Copy, { type Output = Self; fn div(self, other: T) -> Self { @@ -454,7 +453,7 @@ impl Mul for f64 { type Output = Polynomial; fn mul(self, rhs: Polynomial) -> Self::Output { - rhs.mul(self as f64) + rhs.mul(self) } } diff --git a/src/structure/sparse.rs b/src/structure/sparse.rs index 4a026152..e640177f 100644 --- a/src/structure/sparse.rs +++ b/src/structure/sparse.rs @@ -159,6 +159,19 @@ impl LinearAlgebra for SPMatrix { self.to_dense().qr() } + fn svd(&self) -> SVD { + unimplemented!() + } + + #[cfg(feature="O3")] + fn cholesky(&self, _uplo: UPLO) -> Matrix { + unimplemented!() + } + + fn rref(&self) -> Matrix { + self.to_dense().rref() + } + fn det(&self) -> f64 { self.to_dense().det() } @@ -175,10 +188,6 @@ impl LinearAlgebra for SPMatrix { self.to_dense().pseudo_inv() } - fn rref(&self) -> Matrix { - self.to_dense().rref() - } - fn solve(&self, _b: &Vec, _sk: SolveKind) -> Vec { unimplemented!() } @@ -187,15 +196,6 @@ impl LinearAlgebra for SPMatrix { unimplemented!() } - fn svd(&self) -> SVD { - unimplemented!() - } - - #[cfg(feature="O3")] - fn cholesky(&self, uplo: UPLO) -> Matrix { - unimplemented!() - } - fn is_symmetric(&self) -> bool { unimplemented!() } diff --git a/src/structure/vector.rs b/src/structure/vector.rs index a297e8cf..9eb817bd 100644 --- a/src/structure/vector.rs +++ b/src/structure/vector.rs @@ -277,7 +277,6 @@ use crate::traits::{ pointer::{Oxide, Redox, RedoxCommon}, }; use std::cmp::min; -use std::convert; impl FPVector for Vec { type Scalar = f64; @@ -320,7 +319,7 @@ impl FPVector for Vec { fn reduce(&self, init: T, f: F) -> f64 where F: Fn(f64, f64) -> f64, - T: convert::Into, + T: Into, { self.iter().fold(init.into(), |x, &y| f(x, y)) } @@ -492,7 +491,7 @@ impl Algorithm for Vec { /// ``` fn rank(&self) -> Vec { let l = self.len(); - let idx = (1..(l + 1)).map(|x| x as usize).collect::>(); + let idx = (1..(l + 1)).collect::>(); let mut vec_tup = self .clone() @@ -571,7 +570,7 @@ impl Algorithm for Vec { //self.into_iter().enumerate().max_by(|x1, x2| x1.1.partial_cmp(&x2.1).unwrap()).unwrap().0 self.into_iter() .enumerate() - .fold((0usize, std::f64::MIN), |acc, (ics, &val)| { + .fold((0usize, f64::MIN), |acc, (ics, &val)| { if acc.1 < val { (ics, val) } else { @@ -603,7 +602,7 @@ impl Algorithm for Vec { _ => { self.into_iter() .enumerate() - .fold((0usize, std::f64::MAX), |acc, (ics, &val)| { + .fold((0usize, f64::MAX), |acc, (ics, &val)| { if acc.1 > val { (ics, val) } else { @@ -627,7 +626,7 @@ impl Algorithm for Vec { self[idx] } _ => { - self.into_iter().fold(std::f64::MIN, |acc, &val| { + self.into_iter().fold(f64::MIN, |acc, &val| { if acc < val { val } else { @@ -641,7 +640,7 @@ impl Algorithm for Vec { fn min(&self) -> f64 { match () { _ => { - self.into_iter().fold(std::f64::MAX, |acc, &val| { + self.into_iter().fold(f64::MAX, |acc, &val| { if acc > val { val } else { diff --git a/src/traits/math.rs b/src/traits/math.rs index ee7d6c14..228aad9e 100644 --- a/src/traits/math.rs +++ b/src/traits/math.rs @@ -7,8 +7,8 @@ use crate::structure::matrix::Matrix; /// And a space of the vector should closed for that operations. pub trait Vector { type Scalar; - fn add_vec<'a, 'b>(&'a self, rhs: &'b Self) -> Self; - fn sub_vec<'a, 'b>(&'a self, rhs: &'b Self) -> Self; + fn add_vec(&self, rhs: &Self) -> Self; + fn sub_vec(&self, rhs: &Self) -> Self; fn mul_scalar(&self, rhs: Self::Scalar) -> Self; } @@ -71,11 +71,11 @@ pub trait MatrixProduct { impl Vector for f64 { type Scalar = Self; - fn add_vec<'a, 'b>(&'a self, rhs: &'b Self) -> Self { + fn add_vec(&self, rhs: &Self) -> Self { self + rhs } - fn sub_vec<'a, 'b>(&'a self, rhs: &'b Self) -> Self { + fn sub_vec(&self, rhs: &Self) -> Self { self - rhs } diff --git a/src/traits/sugar.rs b/src/traits/sugar.rs index fc28431c..a51d9269 100644 --- a/src/traits/sugar.rs +++ b/src/traits/sugar.rs @@ -40,14 +40,14 @@ where Self::Scalar: Copy + Clone pub trait Scalable { type Vec; - fn resize(&self, size: (usize, usize), shape: Shape) -> Matrix; + fn reshape(&self, size: (usize, usize), shape: Shape) -> Matrix; fn add_col(&self, v: &Self::Vec) -> Matrix; fn add_row(&self, v: &Self::Vec) -> Matrix; } pub trait ScalableMut { type Vec; - fn resize_mut(&mut self, size: (usize, usize), shape: Shape); + fn reshape_mut(&mut self, size: (usize, usize), shape: Shape); fn add_col_mut(&mut self, v: &Self::Vec); fn add_row_mut(&mut self, v: &Self::Vec); } @@ -210,13 +210,13 @@ impl Scalable for Vec { /// /// fn main() { /// let a = c!(1,2,3,4,5,6); - /// let b1 = a.resize((3,2), Row); - /// let b2 = a.resize((3,2), Col); + /// let b1 = a.reshape((3,2), Row); + /// let b2 = a.reshape((3,2), Col); /// assert_eq!(b1, ml_matrix("1 2;3 4;5 6")); /// assert_eq!(b2, ml_matrix("1 4;2 5;3 6")); /// } /// ``` - fn resize(&self, (r, c): (usize, usize), shape: Shape) -> Matrix { + fn reshape(&self, (r, c): (usize, usize), shape: Shape) -> Matrix { assert_eq!(self.len(), r * c); let mut m = zeros_shape(r, c, shape); m.data = self[..].to_vec(); @@ -239,11 +239,11 @@ impl Scalable for Vec { /// assert_eq!(c2, ml_matrix("1 4;2 5;3 6")); /// } /// ``` - fn add_row(&self, v: &Self::Vec) -> Matrix { + fn add_col(&self, v: &Self::Vec) -> Matrix { assert_eq!(self.len(), v.len()); let mut x = self[..].to_vec(); x.extend_from_slice(&v[..]); - x.resize((2, self.len()), Shape::Row) + x.reshape((self.len(), 2), Shape::Col) } /// Vector + Vector = Matrix @@ -262,11 +262,11 @@ impl Scalable for Vec { /// assert_eq!(c2, ml_matrix("1 4;2 5;3 6")); /// } /// ``` - fn add_col(&self, v: &Self::Vec) -> Matrix { + fn add_row(&self, v: &Self::Vec) -> Matrix { assert_eq!(self.len(), v.len()); let mut x = self[..].to_vec(); x.extend_from_slice(&v[..]); - x.resize((self.len(), 2), Shape::Col) + x.reshape((2, self.len()), Shape::Row) } } @@ -283,20 +283,20 @@ impl Scalable for Matrix { /// /// fn main() { /// let a = ml_matrix("1 2 3;4 5 6"); // ml_matrix has shape `Col` - /// let b1 = a.resize((3, 2), Row); - /// let b2 = a.resize((3, 2), Col); + /// let b1 = a.reshape((3, 2), Row); + /// let b2 = a.reshape((3, 2), Col); /// assert_eq!(b1, ml_matrix("1 2;3 4;5 6")); /// assert_eq!(b2, ml_matrix("1 4;2 5;3 6")); /// } /// ``` - fn resize(&self, (r, c): (usize, usize), shape: Shape) -> Matrix { + fn reshape(&self, (r, c): (usize, usize), shape: Shape) -> Matrix { assert_eq!(self.row * self.col, r * c); let mut m = zeros_shape(r, c, shape); m.data = self.data[..].to_vec(); m } - /// Add row + /// Add column /// /// ``` /// #[macro_use] @@ -305,30 +305,30 @@ impl Scalable for Matrix { /// /// fn main() { /// let a = ml_matrix("1 2 3;4 5 6"); - /// let b = c!(7,8,9); - /// let c = a.add_row(&b); - /// assert_eq!(c, ml_matrix("1 2 3;4 5 6;7 8 9")); + /// let b = c!(7,8); + /// let c = a.add_col(&b); + /// assert_eq!(c, ml_matrix("1 2 3 7;4 5 6 8")); /// } /// ``` - fn add_row(&self, v: &Self::Vec) -> Matrix { - assert_eq!(self.col, v.len()); + fn add_col(&self, v: &Self::Vec) -> Matrix { + assert_eq!(self.row, v.len()); match self.shape { - Shape::Row => { + Shape::Col => { let mut m = self.clone(); m.data.extend_from_slice(&v[..]); - m.row += 1; + m.col += 1; m } - Shape::Col => { + Shape::Row => { let mut m = self.change_shape(); m.data.extend_from_slice(&v[..]); - m.row += 1; + m.col += 1; m } } } - /// Add column + /// Add row /// /// ``` /// #[macro_use] @@ -337,24 +337,24 @@ impl Scalable for Matrix { /// /// fn main() { /// let a = ml_matrix("1 2 3;4 5 6"); - /// let b = c!(7,8); - /// let c = a.add_col(&b); - /// assert_eq!(c, ml_matrix("1 2 3 7;4 5 6 8")); + /// let b = c!(7,8,9); + /// let c = a.add_row(&b); + /// assert_eq!(c, ml_matrix("1 2 3;4 5 6;7 8 9")); /// } /// ``` - fn add_col(&self, v: &Self::Vec) -> Matrix { - assert_eq!(self.row, v.len()); + fn add_row(&self, v: &Self::Vec) -> Matrix { + assert_eq!(self.col, v.len()); match self.shape { - Shape::Col => { + Shape::Row => { let mut m = self.clone(); m.data.extend_from_slice(&v[..]); - m.col += 1; + m.row += 1; m } - Shape::Row => { + Shape::Col => { let mut m = self.change_shape(); m.data.extend_from_slice(&v[..]); - m.col += 1; + m.row += 1; m } } @@ -373,20 +373,20 @@ impl ScalableMut for Matrix { /// /// fn main() { /// let mut a = ml_matrix("1 2 3;4 5 6"); // ml_matrix has shape `Row` - /// a.resize_mut((3, 2), Row); + /// a.reshape_mut((3, 2), Row); /// assert_eq!(a, ml_matrix("1 2;3 4;5 6")); - /// a.resize_mut((3, 2), Col); + /// a.reshape_mut((3, 2), Col); /// assert_eq!(a, ml_matrix("1 4;2 5;3 6")); /// } /// ``` - fn resize_mut(&mut self, (r, c): (usize, usize), shape: Shape) { + fn reshape_mut(&mut self, (r, c): (usize, usize), shape: Shape) { assert_eq!(self.row * self.col, r * c); self.row = r; self.col = c; self.shape = shape; } - /// Add row (Mutable) + /// Add column (Mutable) /// /// ``` /// #[macro_use] @@ -395,27 +395,27 @@ impl ScalableMut for Matrix { /// /// fn main() { /// let mut a = ml_matrix("1 2 3;4 5 6"); - /// let b = c!(7,8,9); - /// a.add_row_mut(&b); - /// assert_eq!(a, ml_matrix("1 2 3;4 5 6;7 8 9")); + /// let b = c!(7,8); + /// a.add_col_mut(&b); + /// assert_eq!(a, ml_matrix("1 2 3 7;4 5 6 8")); /// } /// ``` - fn add_row_mut(&mut self, v: &Self::Vec) { - assert_eq!(self.col, v.len()); + fn add_col_mut(&mut self, v: &Self::Vec) { + assert_eq!(self.row, v.len()); match self.shape { - Shape::Row => { + Shape::Col => { self.data.extend_from_slice(&v[..]); - self.row += 1; + self.col += 1; } - Shape::Col => { + Shape::Row => { self.change_shape_mut(); self.data.extend_from_slice(&v[..]); - self.row += 1; + self.col += 1; } } } - /// Add column (Mutable) + /// Add row (Mutable) /// /// ``` /// #[macro_use] @@ -424,22 +424,22 @@ impl ScalableMut for Matrix { /// /// fn main() { /// let mut a = ml_matrix("1 2 3;4 5 6"); - /// let b = c!(7,8); - /// a.add_col_mut(&b); - /// assert_eq!(a, ml_matrix("1 2 3 7;4 5 6 8")); + /// let b = c!(7,8,9); + /// a.add_row_mut(&b); + /// assert_eq!(a, ml_matrix("1 2 3;4 5 6;7 8 9")); /// } /// ``` - fn add_col_mut(&mut self, v: &Self::Vec) { - assert_eq!(self.row, v.len()); + fn add_row_mut(&mut self, v: &Self::Vec) { + assert_eq!(self.col, v.len()); match self.shape { - Shape::Col => { + Shape::Row => { self.data.extend_from_slice(&v[..]); - self.col += 1; + self.row += 1; } - Shape::Row => { + Shape::Col => { self.change_shape_mut(); self.data.extend_from_slice(&v[..]); - self.col += 1; + self.row += 1; } } } diff --git a/src/util/api.rs b/src/util/api.rs index f3c624bf..ab6c2956 100644 --- a/src/util/api.rs +++ b/src/util/api.rs @@ -1,7 +1,6 @@ //! Choose api - MATLAB, R, Python use crate::structure::matrix::*; -use std::convert; pub trait MATLAB { fn new(s: &str) -> Self; @@ -10,13 +9,13 @@ pub trait MATLAB { pub trait PYTHON { fn new(v: Vec>) -> Self where - T: convert::Into + Copy; + T: Into + Copy; } pub trait R { fn new(v: Vec, x: usize, y: usize, shape: Shape) -> Self where - T: convert::Into; + T: Into; } impl MATLAB for Matrix { @@ -47,7 +46,7 @@ impl MATLAB for Matrix { impl PYTHON for Matrix { fn new(v: Vec>) -> Self where - T: convert::Into + Copy, + T: Into + Copy, { let r = v.len(); let c = v[0].len(); @@ -65,7 +64,7 @@ impl PYTHON for Matrix { impl R for Matrix { fn new(v: Vec, x: usize, y: usize, shape: Shape) -> Self where - T: convert::Into, + T: Into, { matrix( v.into_iter().map(|t| t.into()).collect::>(), diff --git a/src/util/low_level.rs b/src/util/low_level.rs index 9151a7f8..8b2a1f52 100644 --- a/src/util/low_level.rs +++ b/src/util/low_level.rs @@ -1,17 +1,17 @@ pub unsafe fn copy_vec_ptr(dst: &mut Vec<*mut f64>, src: &Vec) { - assert!(dst.len() == src.len(), "Should use same length vectors"); + 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>) { - assert!(lhs.len() == rhs.len(), "Should use same length vectors"); + 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<'a>(pv: &'a Vec<*const f64>) -> Vec { +pub unsafe fn ptr_to_vec(pv: &Vec<*const f64>) -> Vec { pv.iter().map(|&x| *x).collect() } diff --git a/src/util/non_macro.rs b/src/util/non_macro.rs index 25ee56c0..3263da53 100644 --- a/src/util/non_macro.rs +++ b/src/util/non_macro.rs @@ -1,4 +1,33 @@ //! Macro to non macro function +//! +//! # R like non-macro functions +//! +//! - seq +//! - seq_with_precision +//! - rbind +//! - cbind +//! +//! # MATLAB like non-macro functions +//! +//! - eye +//! - eye_shape +//! - zeros +//! - zeros_shape +//! - linspace +//! - linspace_with_precision +//! - rand +//! - rand_with_rng +//! +//! # Numpy like non-macro functions +//! +//! - logspace +//! - column_stack +//! - row_stack +//! +//! # Haskell like non-macro functions +//! +//! - concat +//! - cat extern crate rand; use self::rand::prelude::*; @@ -6,13 +35,22 @@ use crate::structure::{ matrix::Shape::{Col, Row}, matrix::{matrix, Matrix, Shape}, }; +use thiserror::Error; use crate::traits::float::FloatWithPrecision; +#[derive(Debug, Copy, Clone, Error)] +pub enum ConcatenateError { + #[error("To concatenate, vectors or matrices must have the same length")] + DifferentLength, +} + +// ┌─────────────────────────────────────────────────────────┐ +// R-like non-macro functions +// └─────────────────────────────────────────────────────────┘ /// R like seq function /// /// # Example /// ``` -/// extern crate peroxide; /// use peroxide::fuga::*; /// /// let a = seq(1, 10, 2); @@ -37,144 +75,42 @@ where let l: usize = factor.floor() as usize + 1; let mut v: Vec = vec![0f64; l]; - for i in 0..l { - v[i] = s + step * (i as f64); + for (i, v) in v.iter_mut().enumerate() { + *v = s + step * (i as f64); } v } -/// MATLAB like zeros (Matrix) +/// Seq with Precision /// -/// # Examples +/// # Example /// ``` -/// extern crate peroxide; /// use peroxide::fuga::*; /// -/// let a = zeros(2, 2); -/// assert_eq!(a, matrix(vec![0f64;4], 2, 2, Row)); -/// ``` -pub fn zeros(r: usize, c: usize) -> Matrix { - matrix(vec![0f64; r * c], r, c, Row) -} - -/// Zeros with custom shape -pub fn zeros_shape(r: usize, c: usize, shape: Shape) -> Matrix { - matrix(vec![0f64; r * c], r, c, shape) -} - -pub fn concat(v1: &Vec, v2: &Vec) -> Vec { - let mut v = v1.clone(); - v.extend_from_slice(&v2[..]); - - v -} - -pub fn cat(val: T, vec: &Vec) -> Vec { - let mut v = vec![val]; - v.extend_from_slice(&vec[..]); - - v -} - -/// MATLAB like eye - Identity matrix -/// -/// # Examples -/// ``` -/// extern crate peroxide; -/// use peroxide::fuga::*; +/// let x = seq(0, 1e-2, 1e-3); +/// assert_ne!(x[9], 0.009); /// -/// let a = eye(2); -/// assert_eq!(a, MATLAB::new("1 0;0 1")); -/// ``` -pub fn eye(n: usize) -> Matrix { - let mut m = zeros(n, n); - for i in 0..n { - m[(i, i)] = 1f64; - } - m -} - -/// eye with custom shape -pub fn eye_shape(n: usize, shape: Shape) -> Matrix { - let mut m = zeros_shape(n, n, shape); - for i in 0..n { - m[(i, i)] = 1f64; - } - m -} - -/// MATLAB like linspace -/// -/// # Examples -/// ``` -/// extern crate peroxide; -/// use peroxide::fuga::*; -/// -/// fn main() { -/// let a = linspace(1, 10, 10); -/// assert_eq!(a, seq(1,10,1)); -/// assert_eq!(a.len(), 10); -/// } -/// ``` -pub fn linspace(start: S, end: T, length: usize) -> Vec -where - S: Into + Copy, - T: Into + Copy, -{ - let step: f64 = if length > 1 { - (end.into() - start.into()) / (length as f64 - 1f64) - } else { - 0f64 - }; - - let mut v = vec![0f64; length]; - v[0] = start.into(); - v[length - 1] = end.into(); - - for i in 1..length - 1 { - v[i] = v[0] + step * (i as f64); - } - v -} - -/// Numpy like logspace -/// -/// # Examples -/// ``` -/// extern crate peroxide; -/// use peroxide::fuga::*; -/// -/// fn main() { -/// let a = logspace(0, 10, 11, 2); -/// let b = vec![1f64, 2f64, 4f64, 8f64, 16f64, 32f64, 64f64, 128f64, 256f64, 512f64, 1024f64]; -/// assert_eq!(a, b); -/// -/// let single = logspace(0f64, 0f64, 1, 10); -/// assert_eq!(single, vec![1f64]); -/// } +/// let x = seq_with_precision(0, 1e-2, 1e-3, 3); +/// assert_eq!(x[9], 0.009); /// ``` -pub fn logspace(start: S, end: T, length: usize, base: U) -> Vec +pub fn seq_with_precision(start: S, end: T, step: U, precision: usize) -> Vec where S: Into + Copy, T: Into + Copy, U: Into + Copy, { - let s: f64 = start.into(); - let e: f64 = end.into(); - let b: f64 = base.into(); + let s = start.into(); + let e = end.into(); + let step = step.into(); assert!(e >= s); - let step: f64 = if length > 1 { - (e - s) / (length as f64 - 1f64) - } else { - 0f64 - }; - - let mut v: Vec = vec![0f64; length]; + let factor: f64 = (e - s) / step; + let l: usize = factor.floor() as usize + 1; + let mut v: Vec = vec![0f64; l]; - for i in 0..length { - v[i] = b.powf(s + step * (i as f64)); + for (i, v) in v.iter_mut().enumerate() { + *v = (s + step * (i as f64)).round_with_precision(precision); } v } @@ -187,14 +123,15 @@ where /// extern crate peroxide; /// use peroxide::fuga::*; /// -/// fn main() { +/// fn main() -> Result<(), Box> { /// let a = matrix!(1;4;1, 2, 2, Col); /// let b = matrix!(5;8;1, 2, 2, Col); /// let c = matrix!(1;8;1, 2, 4, Col); -/// assert_eq!(cbind(a,b), c); +/// assert_eq!(cbind(a,b)?, c); +/// Ok(()) /// } /// ``` -pub fn cbind(m1: Matrix, m2: Matrix) -> Matrix { +pub fn cbind(m1: Matrix, m2: Matrix) -> Result { let mut temp = m1; if temp.shape != Col { temp = temp.change_shape(); @@ -209,10 +146,12 @@ pub fn cbind(m1: Matrix, m2: Matrix) -> Matrix { let mut c = temp.col; let r = temp.row; - assert_eq!(r, temp2.row); + if r != temp2.row { + return Err(ConcatenateError::DifferentLength); + } v.extend_from_slice(&temp2.data[..]); c += temp2.col; - matrix(v, r, c, Col) + Ok(matrix(v, r, c, Col)) } /// R like rbind - concatenate two matrix by row direction @@ -223,14 +162,15 @@ pub fn cbind(m1: Matrix, m2: Matrix) -> Matrix { /// extern crate peroxide; /// use peroxide::fuga::*; /// -/// fn main() { +/// fn main() -> Result<(), Box> { /// let a = matrix!(1;4;1, 2, 2, Row); /// let b = matrix!(5;8;1, 2, 2, Row); /// let c = matrix!(1;8;1, 4, 2, Row); -/// assert_eq!(rbind(a,b), c); +/// assert_eq!(rbind(a,b)?, c); +/// Ok(()) /// } /// ``` -pub fn rbind(m1: Matrix, m2: Matrix) -> Matrix { +pub fn rbind(m1: Matrix, m2: Matrix) -> Result { let mut temp = m1; if temp.shape != Row { temp = temp.change_shape(); @@ -245,61 +185,88 @@ pub fn rbind(m1: Matrix, m2: Matrix) -> Matrix { let c = temp.col; let mut r = temp.row; - assert_eq!(c, temp2.col); + if c != temp2.col { + return Err(ConcatenateError::DifferentLength); + } v.extend_from_slice(&temp2.data[..]); r += temp2.row; - matrix(v, r, c, Row) + Ok(matrix(v, r, c, Row)) } -/// Rand matrix +// ┌─────────────────────────────────────────────────────────┐ +// MATLAB like non-macro functions +// └─────────────────────────────────────────────────────────┘ +/// MATLAB like zeros (Matrix) /// -/// # Description +/// # Examples +/// ``` +/// use peroxide::fuga::*; /// -/// Range = from 0 to 1 -pub fn rand(r: usize, c: usize) -> Matrix { - let mut m = zeros(r, c); - let mut rng = thread_rng(); - for i in 0..r { - for j in 0..c { - m[(i, j)] = rng.gen_range(0f64..=1f64); - } - } - m +/// let a = zeros(2, 2); +/// assert_eq!(a, matrix(vec![0f64;4], 2, 2, Row)); +/// ``` +pub fn zeros(r: usize, c: usize) -> Matrix { + matrix(vec![0f64; r * c], r, c, Row) } -/// Seq with Precision +/// Zeros with custom shape +pub fn zeros_shape(r: usize, c: usize, shape: Shape) -> Matrix { + matrix(vec![0f64; r * c], r, c, shape) +} + +/// MATLAB like eye - Identity matrix /// -/// # Example +/// # Examples /// ``` -/// extern crate peroxide; /// use peroxide::fuga::*; /// -/// fn main() { -/// let x = seq(0, 1e-2, 1e-3); -/// assert_ne!(x[9], 0.009); -/// -/// let x = seq_with_precision(0, 1e-2, 1e-3, 3); -/// assert_eq!(x[9], 0.009); -/// } +/// let a = eye(2); +/// assert_eq!(a, MATLAB::new("1 0;0 1")); /// ``` -pub fn seq_with_precision(start: S, end: T, step: U, precision: usize) -> Vec +pub fn eye(n: usize) -> Matrix { + let mut m = zeros(n, n); + for i in 0..n { + m[(i, i)] = 1f64; + } + m +} + +/// eye with custom shape +pub fn eye_shape(n: usize, shape: Shape) -> Matrix { + let mut m = zeros_shape(n, n, shape); + for i in 0..n { + m[(i, i)] = 1f64; + } + m +} + +/// MATLAB like linspace +/// +/// # Examples +/// ``` +/// use peroxide::fuga::*; +/// +/// let a = linspace(1, 10, 10); +/// assert_eq!(a, seq(1,10,1)); +/// assert_eq!(a.len(), 10); +/// ``` +pub fn linspace(start: S, end: T, length: usize) -> Vec where S: Into + Copy, T: Into + Copy, - U: Into + Copy, { - let s = start.into(); - let e = end.into(); - let step = step.into(); - - assert!(e >= s); + let step: f64 = if length > 1 { + (end.into() - start.into()) / (length as f64 - 1f64) + } else { + 0f64 + }; - let factor: f64 = (e - s) / step; - let l: usize = factor.floor() as usize + 1; - let mut v: Vec = vec![0f64; l]; + let mut v = vec![0f64; length]; + v[0] = start.into(); + v[length - 1] = end.into(); - for i in 0..l { - v[i] = (s + step * (i as f64)).round_with_precision(precision); + for i in 1..length - 1 { + v[i] = v[0] + step * (i as f64); } v } @@ -308,16 +275,13 @@ where /// /// # Example /// ``` -/// extern crate peroxide; /// use peroxide::fuga::*; /// -/// fn main() { -/// let x = linspace(0, 1e-2, 11); -/// assert_ne!(x[9], 0.009); +/// let x = linspace(0, 1e-2, 11); +/// assert_ne!(x[9], 0.009); /// -/// let x = linspace_with_precision(0, 1e-2, 11, 3); -/// assert_eq!(x[9], 0.009); -/// } +/// let x = linspace_with_precision(0, 1e-2, 11, 3); +/// assert_eq!(x[9], 0.009); /// ``` pub fn linspace_with_precision(start: S, end: T, length: usize, precision: usize) -> Vec where @@ -339,3 +303,115 @@ where } v } + +/// Rand matrix +/// +/// # Description +/// +/// Range = from 0 to 1 +pub fn rand(r: usize, c: usize) -> Matrix { + let mut m = zeros(r, c); + let mut rng = thread_rng(); + for i in 0..r { + for j in 0..c { + m[(i, j)] = rng.gen_range(0f64..=1f64); + } + } + m +} + +/// Rand matrix with specific rng +/// +/// # Description +/// +/// Range = from 0 to 1 +pub fn rand_with_rng(r: usize, c: usize, rng: &mut R) -> Matrix { + let mut m = zeros(r, c); + for i in 0..r { + for j in 0..c { + m[(i, j)] = rng.gen_range(0f64..=1f64); + } + } + m +} + +// ┌─────────────────────────────────────────────────────────┐ +// Numpy like non-macro functions +// └─────────────────────────────────────────────────────────┘ +/// Numpy like logspace +/// +/// # Examples +/// ``` +/// use peroxide::fuga::*; +/// +/// let a = logspace(0, 10, 11, 2); +/// let b = vec![1f64, 2f64, 4f64, 8f64, 16f64, 32f64, 64f64, 128f64, 256f64, 512f64, 1024f64]; +/// assert_eq!(a, b); +/// +/// let single = logspace(0f64, 0f64, 1, 10); +/// assert_eq!(single, vec![1f64]); +/// ``` +pub fn logspace(start: S, end: T, length: usize, base: U) -> Vec +where + S: Into + Copy, + T: Into + Copy, + U: Into + Copy, +{ + let s: f64 = start.into(); + let e: f64 = end.into(); + let b: f64 = base.into(); + + assert!(e >= s); + + let step: f64 = if length > 1 { + (e - s) / (length as f64 - 1f64) + } else { + 0f64 + }; + + let mut v: Vec = vec![0f64; length]; + + for (i, v) in v.iter_mut().enumerate() { + *v = b.powf(s + step * (i as f64)); + } + v +} + +/// Numpy like column_stack +pub fn column_stack(v: &[Vec]) -> Result { + let row = v[0].len(); + if v.iter().any(|x| x.len() != row) { + return Err(ConcatenateError::DifferentLength); + } + let data = v.iter().flatten().copied().collect(); + Ok(matrix(data, row, v.len(), Col)) +} + +/// Numpy like row_stack +pub fn row_stack(v: &[Vec]) -> Result { + let col = v[0].len(); + if v.iter().any(|x| x.len() != col) { + return Err(ConcatenateError::DifferentLength); + } + let data = v.iter().flatten().copied().collect(); + Ok(matrix(data, v.len(), col, Row)) +} + +// ┌─────────────────────────────────────────────────────────┐ +// Haskell like non-macro functions +// └─────────────────────────────────────────────────────────┘ +/// Concatenate two vectors into one +pub fn concat(v1: &[T], v2: &[T]) -> Vec { + let mut v = v1.to_vec(); + v.extend_from_slice(v2); + + v +} + +/// Concatenate a value and vector +pub fn cat(val: T, vec: &[T]) -> Vec { + let mut v = vec![val]; + v.extend_from_slice(vec); + + v +} diff --git a/src/util/useful.rs b/src/util/useful.rs index c288debf..b2a5df7f 100644 --- a/src/util/useful.rs +++ b/src/util/useful.rs @@ -1,7 +1,5 @@ //! Useful missing tools -use std::convert; -use std::f64::MAX; use std::ops::Range; // ============================================================================= @@ -19,14 +17,14 @@ use std::ops::Range; /// ``` pub fn nearly_eq(x: S, y: T) -> bool where - S: convert::Into, - T: convert::Into, + S: Into, + T: Into, { let mut b: bool = false; let e = 1e-7; let p: f64 = x.into().abs(); let q: f64 = y.into().abs(); - if (p - q).abs() < e || (p - q).abs() / (p + q).min(MAX) < e { + if (p - q).abs() < e || (p - q).abs() / (p + q).min(f64::MAX) < e { b = true; } b