From 1c5f21bb76c799fa687856121791c2b2a6644f41 Mon Sep 17 00:00:00 2001 From: axect Date: Thu, 11 Apr 2024 23:11:57 +0900 Subject: [PATCH 01/15] CHGE: Start to change whole structure of root module - Remove all of old structures - Define basic traits - Implement Bisection & Newton --- src/grave/root.rs | 556 +++++++++++++++++++++++++++++++++++++ src/numerical/root.rs | 623 +++++++----------------------------------- 2 files changed, 648 insertions(+), 531 deletions(-) create mode 100644 src/grave/root.rs diff --git a/src/grave/root.rs b/src/grave/root.rs new file mode 100644 index 00000000..8769f1c9 --- /dev/null +++ b/src/grave/root.rs @@ -0,0 +1,556 @@ +//! Root Finding +//! +//! # Implemented Algorithms +//! +//! * Bisection +//! * False Position (Regula Falsi) +//! * Secant +//! * Newton +//! +//! # Low-level API +//! +//! ## Members +//! +//! * `RootState` +//! * `P(f64)` : For point-like initial guess +//! * `I(f64, f64)` : For interval-like initial guess +//! * `fn get_point(&self) -> Option` : Returns the point +//! * `fn get_interval(&self) -> Option<(f64, f64)>` : Returns the interval +//! * `RootFind` : Algorithms for root finding +//! * `Bisection` +//! * `FalsePosition` +//! * `Newton` +//! * `Secant` +//! * `RootError` : Error for root finding +//! * `MismatchedState` : Mismatched state and method (ex: Point vs Bisection) +//! * `TimesUp` : No root until `self.times` +//! * `NaNRoot` : NaN +//! * `RootFinder` : Main structure for root finding +//! * `fn new(RootState, RootFind, f)` : Create RootFinder (times: 100, tol: 1e-10) +//! * `fn condition_number(&self) -> f64` : Compute condition number +//! * `fn set_tol(&mut self, f64) -> &mut Self` : Set tolerance +//! * `fn set_times(&mut self, usize) -> &mut Self` : Set max iteration times +//! * `fn get_tol(&self) -> f64` : Get tolerance +//! * `fn get_times(&self) -> usize` : Get max iteration times +//! * `fn get_method(&self) -> RootFind` : Get method +//! * `fn get_curr(&self) -> &RootState` : Get current state +//! * `fn check_find(&self) -> RootBool` : Check if find root +//! * `fn update(&mut self)` : Update one step +//! * `fn find_root(&mut self) -> Result` : Find root +//! +//! ## Usage +//! +//! ``` +//! extern crate peroxide; +//! use peroxide::fuga::*; +//! +//! fn main() -> Result<(), RootError> { +//! let init = RootState::I(1f64, 4f64); +//! let mut rf = RootFinder::new(init, Bisection, f)?; +//! rf.set_tol(1e-15) // Default: 1e-10 +//! .set_times(200); // Default: 100 +//! let root = rf.find_root()?; +//! root.print(); +//! Ok(()) +//! } +//! +//! fn f(x: AD) -> AD { +//! x.sin() +//! } +//! ``` +//! +//! # High-level API +//! +//! ## Members +//! +//! All output type is `Result` +//! +//! * `bisection(f, interval: (f64, f64), times: usize, tol: f64)` +//! * `false_position(f, interval: (f64, f64), times: usize, tol: f64)` +//! * `secant(f, initial_guess: (f64, f64), times: usize, tol: f64)` +//! * `newton(f, initial_guess: f64, times: usize, tol: f64)` +//! +//! ## Usage +//! +//! ``` +//! extern crate peroxide; +//! use peroxide::fuga::*; +//! +//! fn main() -> Result<(), RootError> { +//! let root = bisection(f, (1f64, 4f64), 100, 1e-15)?; +//! root.print(); +//! Ok(()) +//! } +//! +//! fn f(x: AD) -> AD { +//! x.sin() +//! } +//! ``` +//! +//! # Reference +//! +//! * Walter Gautschi, *Numerical Analysis*, Springer (2012) + +use crate::structure::ad::{AD, AD::AD1, ADFn}; +use std::fmt::Display; +use RootState::{I, P}; +use crate::traits::stable::StableFn; +//use std::collections::HashMap; +//use std::marker::PhantomData; + +#[derive(Debug, Copy, Clone)] +pub enum RootState { + P(f64), + I(f64, f64), +} + +#[derive(Debug, Copy, Clone, PartialEq, Eq)] +pub enum RootFind { + Bisection, + Newton, + Secant, + FalsePosition, +} + +/// Structure for Root finding +#[derive(Debug, Clone)] +#[allow(dead_code)] +pub struct RootFinder AD> { + init: RootState, + curr: RootState, + method: RootFind, + f: Box, + find: RootBool, + times: usize, + tol: f64, + root: f64, +} + +#[derive(Debug, Copy, Clone)] +pub enum RootError { + MismatchedState, + TimesUp, + NaNRoot, +} + +#[derive(Debug, Copy, Clone, PartialEq, Eq)] +pub enum RootBool { + Find, + NotYet, + Error, +} + +impl Display for RootError { + fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result { + match self { + RootError::MismatchedState => writeln!(f, "Mismatched RootState with RootFind method"), + RootError::TimesUp => writeln!(f, "No root until set times - More times required"), + RootError::NaNRoot => writeln!(f, "Root is NaN - Should modify initial states"), + } + } +} + +impl std::error::Error for RootError {} + +impl RootState { + pub fn get_point(&self) -> Option { + match self { + P(x) => Some(*x), + _ => None, + } + } + + pub fn get_interval(&self) -> Option<(f64, f64)> { + match self { + I(x, y) => Some((*x, *y)), + _ => None, + } + } +} + +impl AD> RootFinder { + /// Create RootFinder + /// + /// # Default Options + /// + /// * Times: 100 + /// * Tolerance: 1e-10 + /// + /// # Usage + /// ``` + /// extern crate peroxide; + /// use peroxide::fuga::*; + /// + /// fn main() -> Result<(), RootError> { + /// let init = RootState::I(1f64, 3f64); + /// let mut a = RootFinder::new(init, Bisection, f)?; + /// let x = a.find_root()?; + /// x.print(); + /// Ok(()) + /// } + /// + /// fn f(x: AD) -> AD { + /// x.powi(2) - 4f64 + /// } + /// ``` + pub fn new(init: RootState, method: RootFind, f: F) -> Result { + match method { + RootFind::Bisection => match init { + P(_) => Err(RootError::MismatchedState), + _ => Ok(RootFinder { + init, + curr: init, + method, + f: Box::new(f), + find: RootBool::NotYet, + times: 100, + tol: 1e-10, + root: 0f64, + }), + }, + RootFind::Newton => match init { + P(_) => Ok(RootFinder { + init, + curr: init, + method, + f: Box::new(f), + find: RootBool::NotYet, + times: 100, + tol: 1e-10, + root: 0f64, + }), + _ => Err(RootError::MismatchedState), + }, + RootFind::Secant => match init { + P(_) => Err(RootError::MismatchedState), + _ => Ok(RootFinder { + init, + curr: init, + method, + f: Box::new(f), + find: RootBool::NotYet, + times: 100, + tol: 1e-10, + root: 0f64, + }), + }, + RootFind::FalsePosition => match init { + I(_, _) => Ok(RootFinder { + init, + curr: init, + method, + f: Box::new(f), + find: RootBool::NotYet, + times: 100, + tol: 1e-10, + root: 0f64, + }), + _ => Err(RootError::MismatchedState), + }, + } + } + + pub fn f(&self, x: AD) -> AD { + (self.f)(x) + } + + /// Condition number + pub fn condition_number(&self) -> f64 { + match self.curr { + P(p) => { + let z = AD1(p, 1f64); + let fz = self.f(z); + p * fz.dx() / fz.x() + } + I(a, b) => { + let p = (a + b) / 2f64; + let z = AD1(p, 1f64); + let fz = self.f(z); + p * fz.dx() / fz.x() + } + } + } + + /// Set max iteration times + pub fn set_times(&mut self, times: usize) -> &mut Self { + self.times = times; + self + } + + /// Set tolerance + pub fn set_tol(&mut self, tol: f64) -> &mut Self { + self.tol = tol; + self + } + + /// Get max iteration times + pub fn get_times(&self) -> usize { + self.times + } + + /// Get tolerance + pub fn get_tol(&self) -> f64 { + self.tol + } + + /// Get method + pub fn get_method(&self) -> RootFind { + self.method + } + + /// Get current RootState + pub fn get_curr(&self) -> &RootState { + &self.curr + } + + /// Check if find root + pub fn check_find(&self) -> RootBool { + self.find + } + + /// Update RootFinder + #[inline] + pub fn update(&mut self) { + match self.method { + RootFind::Bisection => match self.curr { + I(a, b) => { + let f_ad = ADFn::new(|x| self.f(x)); + let x = 0.5 * (a + b); + let fa = f_ad.call_stable(a); + let fx = f_ad.call_stable(x); + let fb = f_ad.call_stable(b); + if (a - b).abs() <= self.tol { + self.find = RootBool::Find; + self.root = x; + } else { + if fa * fx < 0f64 { + self.curr = I(a, x); + } else if fx * fb < 0f64 { + self.curr = I(x, b); + } else if fx == 0f64 { + self.find = RootBool::Find; + self.root = x; + } else { + self.find = RootBool::Error; + } + } + } + _ => unreachable!(), + }, + RootFind::FalsePosition => match self.curr { + I(a, b) => { + let f_ad = ADFn::new(|x| self.f(x)); + let fa = f_ad.call_stable(a); + let fb = f_ad.call_stable(b); + let x = (a * fb - b * fa) / (fb - fa); + let fx = f_ad.call_stable(x); + if (a - b).abs() <= self.tol || fx.abs() <= self.tol { + self.find = RootBool::Find; + self.root = x; + } else { + if fx * fa < 0f64 { + self.curr = I(a, x); + } else if fx * fb < 0f64 { + self.curr = I(x, b); + } else if fx == 0f64 { + self.find = RootBool::Find; + self.root = x; + } else { + self.find = RootBool::Error; + } + } + } + _ => unreachable!(), + }, + RootFind::Newton => match self.curr { + P(xn) => { + let z = AD1(xn, 1f64); + let fz = self.f(z); + let x = xn - (fz.x() / fz.dx()); + if (x - xn).abs() <= self.tol { + self.find = RootBool::Find; + self.root = x; + } + self.curr = P(x); + } + _ => unreachable!(), + }, + RootFind::Secant => match self.curr { + I(xn_1, xn) => { + let f_ad = ADFn::new(|x| self.f(x)); + let fxn_1 = f_ad.call_stable(xn_1); + let fxn = f_ad.call_stable(xn); + let x = xn - (xn - xn_1) / (fxn - fxn_1) * fxn; + if (x - xn).abs() <= self.tol { + self.find = RootBool::Find; + self.root = x; + } + self.curr = I(xn, x); + } + _ => unreachable!(), + }, + } + } + + /// Find Root + pub fn find_root(&mut self) -> Result { + for _i in 0..self.times { + self.update(); + match self.find { + RootBool::Find => { + //println!("{}", i+1); + return Ok(self.root); + } + RootBool::Error => { + return Err(RootError::NaNRoot); + } + _ => (), + } + } + Err(RootError::TimesUp) + } +} + +/// Bisection method to find root +/// +/// # Usage +/// ``` +/// extern crate peroxide; +/// use peroxide::fuga::*; +/// +/// fn main() -> Result<(), RootError> { +/// let x = bisection(f, (0f64, 4f64), 100, 1e-15)?; +/// assert!((x - 3f64).abs() < 1e-15); +/// Ok(()) +/// } +/// +/// fn f(x: AD) -> AD { +/// x.powi(2) - x * 2f64 - 3f64 +/// } +/// ``` +pub fn bisection AD>( + f: F, + interval: (f64, f64), + times: usize, + tol: f64, +) -> Result { + let (a, b) = interval; + 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() +} + +/// False position method to find root +/// +/// # Usage +/// ``` +/// extern crate peroxide; +/// use peroxide::fuga::*; +/// +/// fn main() -> Result<(), RootError> { +/// let x = false_position(f, (0f64, 4f64), 1000, 1e-15)?; +/// assert!((x - 3f64).abs() < 1e-15); +/// Ok(()) +/// } +/// +/// fn f(x: AD) -> AD { +/// x.powi(2) - x * 2f64 - 3f64 +/// } +/// ``` +pub fn false_position AD>( + f: F, + interval: (f64, f64), + times: usize, + tol: f64, +) -> Result { + let (a, b) = interval; + 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() +} + +/// Newton method to find root +/// +/// # Usage +/// ``` +/// extern crate peroxide; +/// use peroxide::fuga::*; +/// +/// fn main() -> Result<(), RootError> { +/// let x = newton(f, 2f64, 100, 1e-15)?; +/// assert!((x - 3f64).abs() < 1e-15); +/// Ok(()) +/// } +/// +/// fn f(x: AD) -> AD { +/// x.powi(2) - x * 2f64 - 3f64 +/// } +/// ``` +pub fn newton AD>( + f: F, + initial_guess: f64, + times: usize, + tol: f64, +) -> Result { + 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() +} + +/// Secant method to find root +/// +/// # Usage +/// ``` +/// extern crate peroxide; +/// use peroxide::fuga::*; +/// +/// fn main() -> Result<(), RootError> { +/// let x = secant(f, (2f64, 2.1f64), 100, 1e-15)?; +/// assert!((x - 3f64).abs() < 1e-15); +/// Ok(()) +/// } +/// +/// fn f(x: AD) -> AD { +/// x.powi(2) - x * 2f64 - 3f64 +/// } +/// ``` +pub fn secant AD>( + f: F, + initial_guess: (f64, f64), + times: usize, + tol: f64, +) -> Result { + let (a, b) = initial_guess; + 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() +} + +//pub trait RootFinder { +// type InitialState; +// type RootType; +// +// fn mut_update(&mut self); +// fn find_root(&mut self) -> Option; +// fn set_initial_condition(&mut self, init: Self::InitialState) -> &mut Self; +// fn set_times(&mut self) -> &mut Self; +// fn check_enough(&self) -> bool; +//} +// +//#[derive(Debug, Clone, Copy, Hash, PartialOrd, PartialEq, Eq)] +//enum RootOptions { +// InitCondition, +// Times, +//} +// +//pub struct Bisection f64> { +// init: (f64, f64), +// curr: (f64, f64), +// func: Box, +// times: usize, +// options: HashMap, +//} diff --git a/src/numerical/root.rs b/src/numerical/root.rs index 8769f1c9..1462c81a 100644 --- a/src/numerical/root.rs +++ b/src/numerical/root.rs @@ -1,556 +1,117 @@ -//! Root Finding -//! -//! # Implemented Algorithms -//! -//! * Bisection -//! * False Position (Regula Falsi) -//! * Secant -//! * Newton -//! -//! # Low-level API -//! -//! ## Members -//! -//! * `RootState` -//! * `P(f64)` : For point-like initial guess -//! * `I(f64, f64)` : For interval-like initial guess -//! * `fn get_point(&self) -> Option` : Returns the point -//! * `fn get_interval(&self) -> Option<(f64, f64)>` : Returns the interval -//! * `RootFind` : Algorithms for root finding -//! * `Bisection` -//! * `FalsePosition` -//! * `Newton` -//! * `Secant` -//! * `RootError` : Error for root finding -//! * `MismatchedState` : Mismatched state and method (ex: Point vs Bisection) -//! * `TimesUp` : No root until `self.times` -//! * `NaNRoot` : NaN -//! * `RootFinder` : Main structure for root finding -//! * `fn new(RootState, RootFind, f)` : Create RootFinder (times: 100, tol: 1e-10) -//! * `fn condition_number(&self) -> f64` : Compute condition number -//! * `fn set_tol(&mut self, f64) -> &mut Self` : Set tolerance -//! * `fn set_times(&mut self, usize) -> &mut Self` : Set max iteration times -//! * `fn get_tol(&self) -> f64` : Get tolerance -//! * `fn get_times(&self) -> usize` : Get max iteration times -//! * `fn get_method(&self) -> RootFind` : Get method -//! * `fn get_curr(&self) -> &RootState` : Get current state -//! * `fn check_find(&self) -> RootBool` : Check if find root -//! * `fn update(&mut self)` : Update one step -//! * `fn find_root(&mut self) -> Result` : Find root -//! -//! ## Usage -//! -//! ``` -//! extern crate peroxide; -//! use peroxide::fuga::*; -//! -//! fn main() -> Result<(), RootError> { -//! let init = RootState::I(1f64, 4f64); -//! let mut rf = RootFinder::new(init, Bisection, f)?; -//! rf.set_tol(1e-15) // Default: 1e-10 -//! .set_times(200); // Default: 100 -//! let root = rf.find_root()?; -//! root.print(); -//! Ok(()) -//! } -//! -//! fn f(x: AD) -> AD { -//! x.sin() -//! } -//! ``` -//! -//! # High-level API -//! -//! ## Members -//! -//! All output type is `Result` -//! -//! * `bisection(f, interval: (f64, f64), times: usize, tol: f64)` -//! * `false_position(f, interval: (f64, f64), times: usize, tol: f64)` -//! * `secant(f, initial_guess: (f64, f64), times: usize, tol: f64)` -//! * `newton(f, initial_guess: f64, times: usize, tol: f64)` -//! -//! ## Usage -//! -//! ``` -//! extern crate peroxide; -//! use peroxide::fuga::*; -//! -//! fn main() -> Result<(), RootError> { -//! let root = bisection(f, (1f64, 4f64), 100, 1e-15)?; -//! root.print(); -//! Ok(()) -//! } -//! -//! fn f(x: AD) -> AD { -//! x.sin() -//! } -//! ``` -//! -//! # Reference -//! -//! * Walter Gautschi, *Numerical Analysis*, Springer (2012) - -use crate::structure::ad::{AD, AD::AD1, ADFn}; -use std::fmt::Display; -use RootState::{I, P}; -use crate::traits::stable::StableFn; -//use std::collections::HashMap; -//use std::marker::PhantomData; - -#[derive(Debug, Copy, Clone)] -pub enum RootState { - P(f64), - I(f64, f64), +use thiserror::Error; + +#[derive(Error, Debug, Clone)] +pub enum RootFindingError { + #[error("No available root")] + NoAvailableRoot, + #[error("Zero derivative at {0}")] + ZeroDerivative(f64), + #[error("Constrain violated at {0}")] + ConstraintViolation(f64), } -#[derive(Debug, Copy, Clone, PartialEq, Eq)] -pub enum RootFind { - Bisection, - Newton, - Secant, - FalsePosition, +pub trait RootFindingProblem { + fn function(&self, x: f64) -> Result; + fn initial_guess(&self) -> T; + fn derivative(&self, x: f64) -> f64 { + unimplemented!() + } + fn hessian(&self, x: f64) -> f64 { + unimplemented!() + } } -/// Structure for Root finding -#[derive(Debug, Clone)] -#[allow(dead_code)] -pub struct RootFinder AD> { - init: RootState, - curr: RootState, - method: RootFind, - f: Box, - find: RootBool, - times: usize, - tol: f64, - root: f64, +pub trait RootFindingMethod { + fn step>(&self, problem: &P, state: T) -> Result; } -#[derive(Debug, Copy, Clone)] -pub enum RootError { - MismatchedState, - TimesUp, - NaNRoot, +pub trait RootSolver { + fn solve(&self, problem: &P, finder: &F) -> Result + where + P: RootFindingProblem, + F: RootFindingMethod; } -#[derive(Debug, Copy, Clone, PartialEq, Eq)] -pub enum RootBool { - Find, - NotYet, - Error, +// ┌─────────────────────────────────────────────────────────┐ +// Bisection method +// └─────────────────────────────────────────────────────────┘ +pub struct BisectionMethod { + pub max_iterations: usize, + pub tolerance: f64, } -impl Display for RootError { - fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result { - match self { - RootError::MismatchedState => writeln!(f, "Mismatched RootState with RootFind method"), - RootError::TimesUp => writeln!(f, "No root until set times - More times required"), - RootError::NaNRoot => writeln!(f, "Root is NaN - Should modify initial states"), +impl RootFindingMethod<(f64, f64)> for BisectionMethod { + fn step>( + &self, + problem: &P, + state: (f64, f64), + ) -> Result<(f64, f64), RootFindingError> { + let (a, b) = state; + let c = (a + b) / 2.0; + + let fa = problem.function(a)?; + let fb = problem.function(b)?; + let fc = problem.function(c)?; + + if fa * fc < 0.0 { + Ok((a, c)) + } else if fb * fc < 0.0 { + Ok((c, b)) + } else if fc == 0.0 { + Ok((c, c)) + } else { + Err(RootFindingError::NoAvailableRoot) } } } -impl std::error::Error for RootError {} - -impl RootState { - pub fn get_point(&self) -> Option { - match self { - P(x) => Some(*x), - _ => None, - } - } - - pub fn get_interval(&self) -> Option<(f64, f64)> { - match self { - I(x, y) => Some((*x, *y)), - _ => None, - } - } +// ┌─────────────────────────────────────────────────────────┐ +// Newton method +// └─────────────────────────────────────────────────────────┘ +pub struct NewtonMethod { + pub max_iterations: usize, + pub tolerance: f64, } -impl AD> RootFinder { - /// Create RootFinder - /// - /// # Default Options - /// - /// * Times: 100 - /// * Tolerance: 1e-10 - /// - /// # Usage - /// ``` - /// extern crate peroxide; - /// use peroxide::fuga::*; - /// - /// fn main() -> Result<(), RootError> { - /// let init = RootState::I(1f64, 3f64); - /// let mut a = RootFinder::new(init, Bisection, f)?; - /// let x = a.find_root()?; - /// x.print(); - /// Ok(()) - /// } - /// - /// fn f(x: AD) -> AD { - /// x.powi(2) - 4f64 - /// } - /// ``` - pub fn new(init: RootState, method: RootFind, f: F) -> Result { - match method { - RootFind::Bisection => match init { - P(_) => Err(RootError::MismatchedState), - _ => Ok(RootFinder { - init, - curr: init, - method, - f: Box::new(f), - find: RootBool::NotYet, - times: 100, - tol: 1e-10, - root: 0f64, - }), - }, - RootFind::Newton => match init { - P(_) => Ok(RootFinder { - init, - curr: init, - method, - f: Box::new(f), - find: RootBool::NotYet, - times: 100, - tol: 1e-10, - root: 0f64, - }), - _ => Err(RootError::MismatchedState), - }, - RootFind::Secant => match init { - P(_) => Err(RootError::MismatchedState), - _ => Ok(RootFinder { - init, - curr: init, - method, - f: Box::new(f), - find: RootBool::NotYet, - times: 100, - tol: 1e-10, - root: 0f64, - }), - }, - RootFind::FalsePosition => match init { - I(_, _) => Ok(RootFinder { - init, - curr: init, - method, - f: Box::new(f), - find: RootBool::NotYet, - times: 100, - tol: 1e-10, - root: 0f64, - }), - _ => Err(RootError::MismatchedState), - }, - } - } - - pub fn f(&self, x: AD) -> AD { - (self.f)(x) - } - - /// Condition number - pub fn condition_number(&self) -> f64 { - match self.curr { - P(p) => { - let z = AD1(p, 1f64); - let fz = self.f(z); - p * fz.dx() / fz.x() - } - I(a, b) => { - let p = (a + b) / 2f64; - let z = AD1(p, 1f64); - let fz = self.f(z); - p * fz.dx() / fz.x() - } - } - } - - /// Set max iteration times - pub fn set_times(&mut self, times: usize) -> &mut Self { - self.times = times; - self - } - - /// Set tolerance - pub fn set_tol(&mut self, tol: f64) -> &mut Self { - self.tol = tol; - self - } - - /// Get max iteration times - pub fn get_times(&self) -> usize { - self.times - } - - /// Get tolerance - pub fn get_tol(&self) -> f64 { - self.tol - } - - /// Get method - pub fn get_method(&self) -> RootFind { - self.method - } - - /// Get current RootState - pub fn get_curr(&self) -> &RootState { - &self.curr - } - - /// Check if find root - pub fn check_find(&self) -> RootBool { - self.find - } - - /// Update RootFinder - #[inline] - pub fn update(&mut self) { - match self.method { - RootFind::Bisection => match self.curr { - I(a, b) => { - let f_ad = ADFn::new(|x| self.f(x)); - let x = 0.5 * (a + b); - let fa = f_ad.call_stable(a); - let fx = f_ad.call_stable(x); - let fb = f_ad.call_stable(b); - if (a - b).abs() <= self.tol { - self.find = RootBool::Find; - self.root = x; - } else { - if fa * fx < 0f64 { - self.curr = I(a, x); - } else if fx * fb < 0f64 { - self.curr = I(x, b); - } else if fx == 0f64 { - self.find = RootBool::Find; - self.root = x; - } else { - self.find = RootBool::Error; - } - } - } - _ => unreachable!(), - }, - RootFind::FalsePosition => match self.curr { - I(a, b) => { - let f_ad = ADFn::new(|x| self.f(x)); - let fa = f_ad.call_stable(a); - let fb = f_ad.call_stable(b); - let x = (a * fb - b * fa) / (fb - fa); - let fx = f_ad.call_stable(x); - if (a - b).abs() <= self.tol || fx.abs() <= self.tol { - self.find = RootBool::Find; - self.root = x; - } else { - if fx * fa < 0f64 { - self.curr = I(a, x); - } else if fx * fb < 0f64 { - self.curr = I(x, b); - } else if fx == 0f64 { - self.find = RootBool::Find; - self.root = x; - } else { - self.find = RootBool::Error; - } - } - } - _ => unreachable!(), - }, - RootFind::Newton => match self.curr { - P(xn) => { - let z = AD1(xn, 1f64); - let fz = self.f(z); - let x = xn - (fz.x() / fz.dx()); - if (x - xn).abs() <= self.tol { - self.find = RootBool::Find; - self.root = x; - } - self.curr = P(x); - } - _ => unreachable!(), - }, - RootFind::Secant => match self.curr { - I(xn_1, xn) => { - let f_ad = ADFn::new(|x| self.f(x)); - let fxn_1 = f_ad.call_stable(xn_1); - let fxn = f_ad.call_stable(xn); - let x = xn - (xn - xn_1) / (fxn - fxn_1) * fxn; - if (x - xn).abs() <= self.tol { - self.find = RootBool::Find; - self.root = x; - } - self.curr = I(xn, x); - } - _ => unreachable!(), - }, +impl RootFindingMethod for NewtonMethod { + fn step>( + &self, + problem: &P, + state: f64, + ) -> Result { + let f = problem.function(state)?; + let df = problem.derivative(state); + + if df == 0.0 { + return Err(RootFindingError::ZeroDerivative(state)); } - } - /// Find Root - pub fn find_root(&mut self) -> Result { - for _i in 0..self.times { - self.update(); - match self.find { - RootBool::Find => { - //println!("{}", i+1); - return Ok(self.root); - } - RootBool::Error => { - return Err(RootError::NaNRoot); - } - _ => (), - } - } - Err(RootError::TimesUp) + Ok(state - f / df) } } -/// Bisection method to find root -/// -/// # Usage -/// ``` -/// extern crate peroxide; -/// use peroxide::fuga::*; -/// -/// fn main() -> Result<(), RootError> { -/// let x = bisection(f, (0f64, 4f64), 100, 1e-15)?; -/// assert!((x - 3f64).abs() < 1e-15); -/// Ok(()) -/// } -/// -/// fn f(x: AD) -> AD { -/// x.powi(2) - x * 2f64 - 3f64 -/// } -/// ``` -pub fn bisection AD>( - f: F, - interval: (f64, f64), - times: usize, - tol: f64, -) -> Result { - let (a, b) = interval; - 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() -} - -/// False position method to find root -/// -/// # Usage -/// ``` -/// extern crate peroxide; -/// use peroxide::fuga::*; -/// -/// fn main() -> Result<(), RootError> { -/// let x = false_position(f, (0f64, 4f64), 1000, 1e-15)?; -/// assert!((x - 3f64).abs() < 1e-15); -/// Ok(()) -/// } -/// -/// fn f(x: AD) -> AD { -/// x.powi(2) - x * 2f64 - 3f64 -/// } -/// ``` -pub fn false_position AD>( - f: F, - interval: (f64, f64), - times: usize, - tol: f64, -) -> Result { - let (a, b) = interval; - 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() +// ┌─────────────────────────────────────────────────────────┐ +// Secant method +// └─────────────────────────────────────────────────────────┘ +pub struct SecantMethod { + pub max_iterations: usize, + pub tolerance: f64, } -/// Newton method to find root -/// -/// # Usage -/// ``` -/// extern crate peroxide; -/// use peroxide::fuga::*; -/// -/// fn main() -> Result<(), RootError> { -/// let x = newton(f, 2f64, 100, 1e-15)?; -/// assert!((x - 3f64).abs() < 1e-15); -/// Ok(()) -/// } -/// -/// fn f(x: AD) -> AD { -/// x.powi(2) - x * 2f64 - 3f64 -/// } -/// ``` -pub fn newton AD>( - f: F, - initial_guess: f64, - times: usize, - tol: f64, -) -> Result { - 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() -} +impl RootFindingMethod<(f64, f64)> for SecantMethod { + fn step>( + &self, + problem: &P, + state: (f64, f64), + ) -> Result<(f64, f64), RootFindingError> { + let (x0, x1) = state; + let f0 = problem.function(x0)?; + let f1 = problem.function(x1)?; + + if f0 == f1 { + return Err(RootFindingError::ZeroDerivative(x0)); + } -/// Secant method to find root -/// -/// # Usage -/// ``` -/// extern crate peroxide; -/// use peroxide::fuga::*; -/// -/// fn main() -> Result<(), RootError> { -/// let x = secant(f, (2f64, 2.1f64), 100, 1e-15)?; -/// assert!((x - 3f64).abs() < 1e-15); -/// Ok(()) -/// } -/// -/// fn f(x: AD) -> AD { -/// x.powi(2) - x * 2f64 - 3f64 -/// } -/// ``` -pub fn secant AD>( - f: F, - initial_guess: (f64, f64), - times: usize, - tol: f64, -) -> Result { - let (a, b) = initial_guess; - 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() + unimplemented!() + } } - -//pub trait RootFinder { -// type InitialState; -// type RootType; -// -// fn mut_update(&mut self); -// fn find_root(&mut self) -> Option; -// fn set_initial_condition(&mut self, init: Self::InitialState) -> &mut Self; -// fn set_times(&mut self) -> &mut Self; -// fn check_enough(&self) -> bool; -//} -// -//#[derive(Debug, Clone, Copy, Hash, PartialOrd, PartialEq, Eq)] -//enum RootOptions { -// InitCondition, -// Times, -//} -// -//pub struct Bisection f64> { -// init: (f64, f64), -// curr: (f64, f64), -// func: Box, -// times: usize, -// options: HashMap, -//} From ca621ebea5dcf4bba4d7d63dc036d84f972b461b Mon Sep 17 00:00:00 2001 From: axect Date: Thu, 11 Apr 2024 23:14:34 +0900 Subject: [PATCH 02/15] DEP: Add anyhow --- Cargo.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Cargo.toml b/Cargo.toml index 567d134d..7f301446 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -37,6 +37,7 @@ 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 } +anyhow = "1.0.82" [package.metadata.docs.rs] rustdoc-args = [ "--html-in-header", "katex-header.html", "--cfg", "docsrs"] From 5c0657b7737684ac325d975846a766821b5d76e9 Mon Sep 17 00:00:00 2001 From: axect Date: Thu, 11 Apr 2024 23:17:46 +0900 Subject: [PATCH 03/15] CHGE: Change all RootFindingError to anyhow --- src/numerical/root.rs | 32 +++++++++++--------------------- 1 file changed, 11 insertions(+), 21 deletions(-) diff --git a/src/numerical/root.rs b/src/numerical/root.rs index 1462c81a..22db61c3 100644 --- a/src/numerical/root.rs +++ b/src/numerical/root.rs @@ -1,17 +1,7 @@ -use thiserror::Error; - -#[derive(Error, Debug, Clone)] -pub enum RootFindingError { - #[error("No available root")] - NoAvailableRoot, - #[error("Zero derivative at {0}")] - ZeroDerivative(f64), - #[error("Constrain violated at {0}")] - ConstraintViolation(f64), -} +use anyhow::{Result, bail}; pub trait RootFindingProblem { - fn function(&self, x: f64) -> Result; + fn function(&self, x: f64) -> Result; fn initial_guess(&self) -> T; fn derivative(&self, x: f64) -> f64 { unimplemented!() @@ -22,11 +12,11 @@ pub trait RootFindingProblem { } pub trait RootFindingMethod { - fn step>(&self, problem: &P, state: T) -> Result; + fn step>(&self, problem: &P, state: T) -> Result; } pub trait RootSolver { - fn solve(&self, problem: &P, finder: &F) -> Result + fn solve(&self, problem: &P, finder: &F) -> Result where P: RootFindingProblem, F: RootFindingMethod; @@ -45,7 +35,7 @@ impl RootFindingMethod<(f64, f64)> for BisectionMethod { &self, problem: &P, state: (f64, f64), - ) -> Result<(f64, f64), RootFindingError> { + ) -> Result<(f64, f64)> { let (a, b) = state; let c = (a + b) / 2.0; @@ -60,7 +50,7 @@ impl RootFindingMethod<(f64, f64)> for BisectionMethod { } else if fc == 0.0 { Ok((c, c)) } else { - Err(RootFindingError::NoAvailableRoot) + bail!("There is no root in the interval [{}, {}]", a, b); } } } @@ -78,12 +68,12 @@ impl RootFindingMethod for NewtonMethod { &self, problem: &P, state: f64, - ) -> Result { + ) -> Result { let f = problem.function(state)?; let df = problem.derivative(state); if df == 0.0 { - return Err(RootFindingError::ZeroDerivative(state)); + bail!("Zero derivative at x = {}", state); } Ok(state - f / df) @@ -103,15 +93,15 @@ impl RootFindingMethod<(f64, f64)> for SecantMethod { &self, problem: &P, state: (f64, f64), - ) -> Result<(f64, f64), RootFindingError> { + ) -> Result<(f64, f64)> { let (x0, x1) = state; let f0 = problem.function(x0)?; let f1 = problem.function(x1)?; if f0 == f1 { - return Err(RootFindingError::ZeroDerivative(x0)); + bail!("Zero secant at ({}, {})", x0, x1); } - unimplemented!() + Ok((x1, x1 - f1 * (x1 - x0) / (f1 - f0))) } } From 46e6b033d8ca44fb73bb5e9133e80d28ad86c85a Mon Sep 17 00:00:00 2001 From: axect Date: Fri, 12 Apr 2024 00:01:17 +0900 Subject: [PATCH 04/15] DOCS: Change README & CITATION.cff --- CITATION.cff | 65 ++++++-------- README.md | 245 +++++++++++++++++++++++++++------------------------ 2 files changed, 157 insertions(+), 153 deletions(-) diff --git a/CITATION.cff b/CITATION.cff index fddd8fbf..f3a911cc 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -1,40 +1,31 @@ -# This CITATION.cff file was generated with cffinit. -# Visit https://bit.ly/cffinit to generate yours today! - +abstract: Peroxide is a comprehensive numeric library written in Rust, designed to + cater to the needs of scientists, engineers, mathematicians, and anyone who desires + high performance numerical computation. The library provides robust and efficient + functionality for linear algebra, numerical analysis, statistics, and more. Peroxide + leverages the Rust language's safety, concurrency, and performance capabilities + to provide an interface that is both user-friendly and highly performant. It is + designed with simplicity in mind, aiming to offer the ease-of-use found in high-level + languages without sacrificing the speed and precision demanded by complex numerical + tasks. Peroxide stands as an indispensable tool for those who seek a performant, + safe and robust numeric computation solution. +authors: +- affiliation: Yonsei University + family-names: Kim + given-names: Tae-Geun + orcid: 0009-0000-4229-2935 cff-version: 1.2.0 +date-released: '2024-04-11' +doi: 10.5281/zenodo.10960095 +keywords: +- Rust +- Numeric +- Integration +- Linear algebra +- Differential equation +license: +- apache-2.0 +- mit +repository-code: https://github.com/Axect/Peroxide/tree/v0.36.4 title: Peroxide -message: >- - If you use this software, please cite it using the - metadata from this file. type: software -authors: - - given-names: Tae-Geun - family-names: Kim - email: ax2ct@outlook.com - affiliation: Yonsei University - orcid: 'https://orcid.org/0009-0000-4229-2935' -repository-code: 'https://github.com/Axect/Peroxide' -url: 'https://crates.io/crates/peroxide' -abstract: >- - Peroxide is a comprehensive numeric library written in - Rust, designed to cater to the needs of scientists, - engineers, mathematicians, and anyone who desires high - performance numerical computation. The library provides - robust and efficient functionality for linear algebra, - numerical analysis, statistics, and more. Peroxide - leverages the Rust language's safety, concurrency, and - performance capabilities to provide an interface that is - both user-friendly and highly performant. It is designed - with simplicity in mind, aiming to offer the ease-of-use - found in high-level languages without sacrificing the - speed and precision demanded by complex numerical tasks. - Peroxide stands as an indispensable tool for those who - seek a performant, safe and robust numeric computation - solution. -keywords: - - Rust - - Numeric - - Integration - - Linear algebra - - Differential equation -license: Apache-2.0 +version: v0.36.4 diff --git a/README.md b/README.md index f1d9f1d3..2f0ba09d 100644 --- a/README.md +++ b/README.md @@ -2,8 +2,8 @@ [![On crates.io](https://img.shields.io/crates/v/peroxide.svg)](https://crates.io/crates/peroxide) [![On docs.rs](https://docs.rs/peroxide/badge.svg)](https://axect.github.io/Peroxide_Doc) -[![builds.sr.ht status](https://builds.sr.ht/~axect/Peroxide/.build.yml.svg)](https://builds.sr.ht/~axect/Peroxide/.build.yml?) -![github](https://github.com/Axect/Peroxide/workflows/Github/badge.svg) +[![DOI](https://zenodo.org/badge/130400565.svg)](https://zenodo.org/doi/10.5281/zenodo.10815823) +![github](https://github.com/Axect/Peroxide/workflows/Github/badge.svg) ![maintenance](https://img.shields.io/badge/maintenance-actively--developed-brightgreen.svg) @@ -12,13 +12,13 @@ Rust numeric library contains linear algebra, numerical analysis, statistics and ## Table of Contents - [Why Peroxide?](#why-peroxide) - - [1. Customize features](#1-customize-features) - - [2. Easy to optimize](#2-easy-to-optimize) - - [3. Friendly syntax](#3-friendly-syntax) - - [4. Can choose two different coding styles](#4-can-choose-two-different-coding-styles) - - [5. Batteries included](#5-batteries-included) - - [6. Compatible with Mathematics](#6-compatible-with-mathematics) - - [7. Written in Rust](#7-written-in-rust) + - [1. Customize features](#1-customize-features) + - [2. Easy to optimize](#2-easy-to-optimize) + - [3. Friendly syntax](#3-friendly-syntax) + - [4. Can choose two different coding styles](#4-can-choose-two-different-coding-styles) + - [5. Batteries included](#5-batteries-included) + - [6. Compatible with Mathematics](#6-compatible-with-mathematics) + - [7. Written in Rust](#7-written-in-rust) - [Latest README version](#latest-readme-version) - [Pre-requisite](#pre-requisite) - [Install](#install) @@ -26,7 +26,7 @@ Rust numeric library contains linear algebra, numerical analysis, statistics and - [Module Structure](#module-structure) - [Documentation](#documentation) - [Examples](#examples) -- [Version Info](#version-info) +- [Release Info](#release-info) - [Contributes Guide](#contributes-guide) - [LICENSE](#license) - [TODO](#todo) @@ -37,13 +37,13 @@ Rust numeric library contains linear algebra, numerical analysis, statistics and Peroxide provides various features. -* `default` - Pure Rust (No dependencies of architecture - Perfect cross compilation) -* `O3` - BLAS & LAPACK (Perfect performance but little bit hard to set-up - Strongly recommend to look [Peroxide with BLAS](https://github.com/Axect/Peroxide_BLAS)) -* `plot` - With matplotlib of python, we can draw any plots. -* `nc` - To handle netcdf file format with DataFrame -* `csv` - To handle csv file format with Matrix or DataFrame -* `parquet` - To handle parquet file format with DataFrame -* `serde` - serialization with [Serde](https://serde.rs/). +- `default` - Pure Rust (No dependencies of architecture - Perfect cross compilation) +- `O3` - BLAS & LAPACK (Perfect performance but little bit hard to set-up - Strongly recommend to look [Peroxide with BLAS](https://github.com/Axect/Peroxide_BLAS)) +- `plot` - With matplotlib of python, we can draw any plots. +- `nc` - To handle netcdf file format with DataFrame +- `csv` - To handle csv file format with Matrix or DataFrame +- `parquet` - To handle parquet file format with DataFrame +- `serde` - serialization with [Serde](https://serde.rs/). If you want to do high performance computation and more linear algebra, then choose `O3` feature. If you don't want to depend C/C++ or Fortran libraries, then choose `default` feature. @@ -84,7 +84,7 @@ fn main() { z[(0,1)] = 2.0; z[(1,0)] = 3.0; z[(1,1)] = 4.0; - + // Simple but effective operations let c = a * b; // Matrix multiplication (BLAS integrated) @@ -106,8 +106,8 @@ fn main() { In peroxide, there are two different options. -* `prelude`: To simple use. -* `fuga`: To choose numerical algorithms explicitly. +- `prelude`: To simple use. +- `fuga`: To choose numerical algorithms explicitly. For examples, let's see norm. @@ -120,7 +120,7 @@ use peroxide::prelude::*; fn main() { let a = c!(1, 2, 3); let l2 = a.norm(); // L2 is default vector norm - + assert_eq!(l2, 14f64.sqrt()); } ``` @@ -130,7 +130,7 @@ In `fuga`, use various norms. But you should write a little bit longer than `pre #[macro_use] extern crate peroxide; use peroxide::fuga::*; - + fn main() { let a = c!(1, 2, 3); let l1 = a.norm(Norm::L1); @@ -146,86 +146,86 @@ fn main() { Peroxide can do many things. -* Linear Algebra - * Effective Matrix structure - * Transpose, Determinant, Diagonal - * LU Decomposition, Inverse matrix, Block partitioning - * QR Decomposition (`O3` feature) - * Singular Value Decomposition (SVD) (`O3` feature) - * Cholesky Decomposition (`O3` feature) - * Reduced Row Echelon form - * Column, Row operations - * Eigenvalue, Eigenvector -* Functional Programming - * Easier functional programming with `Vec` - * For matrix, there are three maps - * `fmap` : map for all elements - * `col_map` : map for column vectors - * `row_map` : map for row vectors -* Automatic Differentiation - * Taylor mode Forward AD - for nth order AD - * Exact jacobian - * `Real` trait to constrain for `f64` and `AD` (for ODE) -* Numerical Analysis - * Lagrange interpolation - * Splines - * Cubic Spline - * Cubic Hermite Spline - * Estimate slope via Akima - * Estimate slope via Quadratic interpolation - * Non-linear regression - * Gradient Descent - * Levenberg Marquardt - * Ordinary Differential Equation - * Runge-Kutta 4th order (`RK4`) - * Runge-Kutta-Fehlberg 4/5th order (`RKF45`) - * Gauss-Legendre 4th order (`GL4`; Implicit) - * Numerical Integration - * Newton-Cotes Quadrature - * Gauss-Legendre Quadrature (up to 30 order) - * Gauss-Kronrod Quadrature (Adaptive) - * G7K15, G10K21, G15K31, G20K41, G25K51, G30K61 - * Gauss-Kronrod Quadrature (Relative tolerance) - * G7K15R, G10K21R, G15K31R, G20K41R, G25K51R, G30K61R - * Root Finding - * Bisection - * False Position (Regula Falsi) - * Secant - * Newton -* Statistics - * More easy random with `rand` crate - * Ordered Statistics - * Median - * Quantile (Matched with R quantile) - * Probability Distributions - * Bernoulli - * Uniform - * Binomial - * Normal - * Gamma - * Beta - * Student's-t - * Weighted Uniform - * RNG algorithms - * Acceptance Rejection - * Marsaglia Polar - * Ziggurat - * Wrapper for `rand-dist` crate - * Confusion Matrix & Metrics -* Special functions - * Wrapper for `puruspe` crate (pure rust) -* Utils - * R-like macro & functions - * Matlab-like macro & functions - * Numpy-like macro & functions - * Julia-like macro & functions -* Plotting - * With `pyo3` & `matplotlib` -* DataFrame - * Support various types simultaneously - * Read & Write `csv` files (`csv` feature) - * Read & Write `netcdf` files (`nc` feature) - * Read & Write `parquet` files (`parquet` feature) +- Linear Algebra + - Effective Matrix structure + - Transpose, Determinant, Diagonal + - LU Decomposition, Inverse matrix, Block partitioning + - QR Decomposition (`O3` feature) + - Singular Value Decomposition (SVD) (`O3` feature) + - Cholesky Decomposition (`O3` feature) + - Reduced Row Echelon form + - Column, Row operations + - Eigenvalue, Eigenvector +- Functional Programming + - Easier functional programming with `Vec` + - For matrix, there are three maps + - `fmap` : map for all elements + - `col_map` : map for column vectors + - `row_map` : map for row vectors +- Automatic Differentiation + - Taylor mode Forward AD - for nth order AD + - Exact jacobian + - `Real` trait to constrain for `f64` and `AD` (for ODE) +- Numerical Analysis + - Lagrange interpolation + - Splines + - Cubic Spline + - Cubic Hermite Spline + - Estimate slope via Akima + - Estimate slope via Quadratic interpolation + - Non-linear regression + - Gradient Descent + - Levenberg Marquardt + - Ordinary Differential Equation + - Runge-Kutta 4th order (`RK4`) + - Runge-Kutta-Fehlberg 4/5th order (`RKF45`) + - Gauss-Legendre 4th order (`GL4`; Implicit) + - Numerical Integration + - Newton-Cotes Quadrature + - Gauss-Legendre Quadrature (up to 30 order) + - Gauss-Kronrod Quadrature (Adaptive) + - G7K15, G10K21, G15K31, G20K41, G25K51, G30K61 + - Gauss-Kronrod Quadrature (Relative tolerance) + - G7K15R, G10K21R, G15K31R, G20K41R, G25K51R, G30K61R + - Root Finding + - Bisection + - False Position (Regula Falsi) + - Secant + - Newton +- Statistics + - More easy random with `rand` crate + - Ordered Statistics + - Median + - Quantile (Matched with R quantile) + - Probability Distributions + - Bernoulli + - Uniform + - Binomial + - Normal + - Gamma + - Beta + - Student's-t + - Weighted Uniform + - RNG algorithms + - Acceptance Rejection + - Marsaglia Polar + - Ziggurat + - Wrapper for `rand-dist` crate + - Confusion Matrix & Metrics +- Special functions + - Wrapper for `puruspe` crate (pure rust) +- Utils + - R-like macro & functions + - Matlab-like macro & functions + - Numpy-like macro & functions + - Julia-like macro & functions +- Plotting + - With `pyo3` & `matplotlib` +- DataFrame + - Support various types simultaneously + - Read & Write `csv` files (`csv` feature) + - Read & Write `netcdf` files (`nc` feature) + - Read & Write `parquet` files (`parquet` feature) ### 6. Compatible with Mathematics @@ -301,60 +301,73 @@ Corresponding to `0.36.0` ## Pre-requisite -* For `O3` feature - Need `OpenBLAS` -* For `plot` feature - Need `matplotlib` and optional `scienceplots` (for publication quality) -* For `nc` feature - Need `netcdf` +- For `O3` feature - Need `OpenBLAS` +- For `plot` feature - Need `matplotlib` and optional `scienceplots` (for publication quality) +- For `nc` feature - Need `netcdf` ## Install -* Run below commands in your project directory +- Run below commands in your project directory 1. Default + ```bash cargo add peroxide ``` + 2. OpenBLAS + ```bash cargo add peroxide --features O3 ``` + 3. Plot + ```bash cargo add peroxide --features plot ``` + 4. NetCDF dependency for DataFrame + ```bash cargo add peroxide --features nc ``` + 5. CSV dependency for DataFrame + ```bash cargo add peroxide --features csv ``` + 6. Parquet dependency for DataFrame + ```bash cargo add peroxide --features parquet ``` 7. Serialize or Deserialize with Matrix or polynomial + ```bash cargo add peroxide --features serde ``` 8. All features + ```bash cargo add peroxide --features "O3 plot nc csv parquet serde" ``` ## Useful tips for features -* 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 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`. +- 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. +- 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. + - 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. + - 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. ## Module Structure @@ -368,7 +381,7 @@ Corresponding to `0.36.0` - [matlab_macro.rs](src/macros/matlab_macro.rs) : MATLAB like macro - [mod.rs](src/macros/mod.rs) - [r_macro.rs](src/macros/r_macro.rs) : R like macro - - __ml__ : For machine learning (*Beta*) + - __ml__ : For machine learning (_Beta_) - [mod.rs](src/ml/mod.rs) - [reg.rs](src/ml/reg.rs) : Regression tools - __numerical__ : To do numerical things @@ -400,7 +413,7 @@ Corresponding to `0.36.0` - [dataframe.rs](src/structure/dataframe.rs) : Dataframe - [matrix.rs](src/structure/matrix.rs) : Matrix - [polynomial.rs](src/structure/polynomial.rs) : Polynomial - - [sparse.rs](src/structure/sparse.rs) : For sparse structure (*Beta*) + - [sparse.rs](src/structure/sparse.rs) : For sparse structure (_Beta_) - [vector.rs](src/structure/vector.rs) : Extra tools for `Vec` - __traits__ - [mod.rs](src/traits/mod.rs) @@ -425,7 +438,7 @@ Corresponding to `0.36.0` ## Documentation -* [![On docs.rs](https://docs.rs/peroxide/badge.svg)](https://axect.github.io/Peroxide_Doc) +- [![On docs.rs](https://docs.rs/peroxide/badge.svg)](https://axect.github.io/Peroxide_Doc) ## Examples @@ -433,9 +446,9 @@ Corresponding to `0.36.0` - More examples are in [Peroxide Gallery](https://github.com/Axect/Peroxide_Gallery). -## Version Info +## Release Info -To see [RELEASES.md](./RELEASES.md) +To see [RELEASES.md](RELEASES.md) ## Contributes Guide From 6e7b13c41851a2283d1f9a279d5fdc8a0460f27a Mon Sep 17 00:00:00 2001 From: axect Date: Fri, 12 Apr 2024 09:57:41 +0900 Subject: [PATCH 05/15] IMPL: Completely generic RootFinding traits & useful type alias & macros --- src/numerical/root.rs | 166 +++++++++++++++++++++++++++++++++++------- 1 file changed, 141 insertions(+), 25 deletions(-) diff --git a/src/numerical/root.rs b/src/numerical/root.rs index 22db61c3..2cb20fed 100644 --- a/src/numerical/root.rs +++ b/src/numerical/root.rs @@ -1,37 +1,102 @@ use anyhow::{Result, bail}; -pub trait RootFindingProblem { - fn function(&self, x: f64) -> Result; +type PT = [f64; N]; +type INTV = (PT, PT); +type J = [[f64; C]; R]; +type H = [[[f64; C]; C]; R]; + +/// Trait to define a root finding problem +/// +/// # Type Parameters +/// +/// - `I`: Input type (e.g. `f64`, `[f64; N]`, or etc.) +/// - `O`: Output type (e.g. `f64`, `[f64; N]`, or etc.) +/// - `T`: State type (e.g. `f64`, `(f64, f64)`, or etc.) +pub trait RootFindingProblem { + fn function(&self, x: PT) -> Result>; fn initial_guess(&self) -> T; - fn derivative(&self, x: f64) -> f64 { + fn derivative(&self, x: PT) -> Result> { unimplemented!() } - fn hessian(&self, x: f64) -> f64 { + fn hessian(&self, x: PT) -> Result> { unimplemented!() } } -pub trait RootFindingMethod { - fn step>(&self, problem: &P, state: T) -> Result; +pub trait RootFindingMethod { + fn step>(&self, problem: &P, state: T) -> Result; } pub trait RootSolver { - fn solve(&self, problem: &P, finder: &F) -> Result + fn solve(&self, problem: &P, finder: &F) -> Result<[f64; I]> where - P: RootFindingProblem, - F: RootFindingMethod; + P: RootFindingProblem, + F: RootFindingMethod; +} + +/// Macro for single function +/// +/// # Description +/// +/// For I=1, O=1, it is bother to write below code. +/// +/// ```no_run +/// let fx = problem.function([x])?[0]; +/// ``` +/// +/// This macro solve this problem as follows. +/// +/// ```no_run +/// let fx = single_function!(problem, x); +/// ``` +macro_rules! single_function { + ($problem:expr, $x:expr) => {{ + $problem.function([$x])?[0] + }} +} + +/// Macro for single derivative +/// +/// # Description +/// +/// For I=1, O=1, it is bother to write below code. +/// +/// ```no_run +/// let fx = problem.derivative([x])?[0][0]; +/// ``` +/// +/// This macro solve this problem as follows. +/// +/// ```no_run +/// let fx = single_derivative!(problem, x); +/// ``` +macro_rules! single_derivative { + ($problem:expr, $x:expr) => {{ + $problem.derivative([$x])?[0][0] + }} } // ┌─────────────────────────────────────────────────────────┐ // Bisection method // └─────────────────────────────────────────────────────────┘ +/// Bisection method +/// +/// # Arguments +/// +/// - `max_iterations`: Maximum number of iterations +/// - `tolerance`: Tolerance +/// +/// # Caution +/// +/// - The function should be continuous +/// - The function should have a root in the initial interval pub struct BisectionMethod { pub max_iterations: usize, pub tolerance: f64, } -impl RootFindingMethod<(f64, f64)> for BisectionMethod { - fn step>( +impl RootFindingMethod<1, 1, (f64, f64)> for BisectionMethod { + fn step>( &self, problem: &P, state: (f64, f64), @@ -39,9 +104,9 @@ impl RootFindingMethod<(f64, f64)> for BisectionMethod { let (a, b) = state; let c = (a + b) / 2.0; - let fa = problem.function(a)?; - let fb = problem.function(b)?; - let fc = problem.function(c)?; + let fa = single_function!(problem, a); + let fb = single_function!(problem, b); + let fc = single_function!(problem, c); if fa * fc < 0.0 { Ok((a, c)) @@ -58,45 +123,67 @@ impl RootFindingMethod<(f64, f64)> for BisectionMethod { // ┌─────────────────────────────────────────────────────────┐ // Newton method // └─────────────────────────────────────────────────────────┘ +/// Newton method +/// +/// # Arguments +/// +/// - `max_iterations`: Maximum number of iterations +/// - `tolerance`: Tolerance +/// +/// # Caution +/// +/// - The function should be differentiable +/// - This method highly depends on the initial guess +/// - This method is not guaranteed to converge pub struct NewtonMethod { pub max_iterations: usize, pub tolerance: f64, } -impl RootFindingMethod for NewtonMethod { - fn step>( +impl RootFindingMethod<1, 1, f64> for NewtonMethod { + fn step>( &self, problem: &P, - state: f64, + x: f64, ) -> Result { - let f = problem.function(state)?; - let df = problem.derivative(state); + let f = single_function!(problem, x); + let df = single_derivative!(problem, x); if df == 0.0 { - bail!("Zero derivative at x = {}", state); + bail!("Zero derivative at x = {}", x); } - Ok(state - f / df) + Ok(x - f / df) } } // ┌─────────────────────────────────────────────────────────┐ // Secant method // └─────────────────────────────────────────────────────────┘ +/// Secant method +/// +/// # Arguments +/// +/// - `max_iterations`: Maximum number of iterations +/// - `tolerance`: Tolerance +/// +/// # Caution +/// +/// - The function should be differentiable pub struct SecantMethod { pub max_iterations: usize, pub tolerance: f64, } -impl RootFindingMethod<(f64, f64)> for SecantMethod { - fn step>( +impl RootFindingMethod<1, 1, (f64, f64)> for SecantMethod { + fn step>( &self, problem: &P, state: (f64, f64), ) -> Result<(f64, f64)> { let (x0, x1) = state; - let f0 = problem.function(x0)?; - let f1 = problem.function(x1)?; + let f0 = single_function!(problem, x0); + let f1 = single_function!(problem, x1); if f0 == f1 { bail!("Zero secant at ({}, {})", x0, x1); @@ -105,3 +192,32 @@ impl RootFindingMethod<(f64, f64)> for SecantMethod { Ok((x1, x1 - f1 * (x1 - x0) / (f1 - f0))) } } + +// ┌─────────────────────────────────────────────────────────┐ +// Broyden method +// └─────────────────────────────────────────────────────────┘ +/// Broyden method +/// +/// # Arguments +/// +/// - `max_iterations`: Maximum number of iterations +/// - `tolerance`: Tolerance +/// +/// # Caution +/// +/// - The function should be differentiable +pub struct BroydenMethod { + pub max_iterations: usize, + pub tolerance: f64, +} + + +impl RootFindingMethod> for BroydenMethod { + fn step>>( + &self, + problem: &P, + state: INTV, + ) -> Result> { + unimplemented!() + } +} From f9ab4a68d705359bbd5818eb37442d54cd70fcfe Mon Sep 17 00:00:00 2001 From: axect Date: Sat, 13 Apr 2024 19:19:45 +0900 Subject: [PATCH 06/15] IMPL: Refine and stablize RootFinding traits --- src/numerical/root.rs | 303 +++++++++++++++++++++++++++++++++--------- 1 file changed, 239 insertions(+), 64 deletions(-) diff --git a/src/numerical/root.rs b/src/numerical/root.rs index 2cb20fed..f782534b 100644 --- a/src/numerical/root.rs +++ b/src/numerical/root.rs @@ -23,15 +23,27 @@ pub trait RootFindingProblem { } } -pub trait RootFindingMethod { - fn step>(&self, problem: &P, state: T) -> Result; +pub trait RootFinder { + fn max_iter(&self) -> usize; + fn tol(&self) -> f64; + fn find>(&self, problem: &P) -> Result>; } -pub trait RootSolver { - fn solve(&self, problem: &P, finder: &F) -> Result<[f64; I]> - where - P: RootFindingProblem, - F: RootFindingMethod; +#[derive(Debug, Copy, Clone)] +pub enum RootError { + NotConverge(PT), + NoRoot, + ZeroDerivative(PT), +} + +impl std::fmt::Display for RootError<1> { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + RootError::NoRoot => write!(f, "There is no root in the interval"), + RootError::NotConverge(a) => write!(f, "Not yet converge. Our guess is {:?}", a), + RootError::ZeroDerivative(a) => write!(f, "Zero derivative in {:?}", a), + } + } } /// Macro for single function @@ -83,40 +95,58 @@ macro_rules! single_derivative { /// /// # Arguments /// -/// - `max_iterations`: Maximum number of iterations -/// - `tolerance`: Tolerance +/// - `max_iter`: Maximum number of iterations +/// - `tol`: tol /// /// # Caution /// /// - The function should be continuous /// - The function should have a root in the initial interval pub struct BisectionMethod { - pub max_iterations: usize, - pub tolerance: f64, + pub max_iter: usize, + pub tol: f64, } -impl RootFindingMethod<1, 1, (f64, f64)> for BisectionMethod { - fn step>( +impl RootFinder<1, 1, (f64, f64)> for BisectionMethod { + fn max_iter(&self) -> usize { + self.max_iter + } + + fn tol(&self) -> f64 { + self.tol + } + + fn find>( &self, problem: &P, - state: (f64, f64), - ) -> Result<(f64, f64)> { - let (a, b) = state; - let c = (a + b) / 2.0; + ) -> Result<[f64; 1]> { + let state = problem.initial_guess(); + let (mut a, mut b) = state; + let mut fa = single_function!(problem, a); + let mut fb = single_function!(problem, b); + + if fa * fb > 0.0 { + bail!(RootError::NoRoot); + } + + for _ in 0..self.max_iter { + let c = (a + b) / 2.0; + let fc = single_function!(problem, c); - let fa = single_function!(problem, a); - let fb = single_function!(problem, b); - let fc = single_function!(problem, c); - - if fa * fc < 0.0 { - Ok((a, c)) - } else if fb * fc < 0.0 { - Ok((c, b)) - } else if fc == 0.0 { - Ok((c, c)) - } else { - bail!("There is no root in the interval [{}, {}]", a, b); + if fc.abs() < self.tol { + return Ok([c]); + } else if fa * fc < 0.0 { + b = c; + fb = fc; + } else if fb * fc < 0.0 { + a = c; + fa = fc; + } else { + bail!(RootError::NoRoot); + } } + let c = (a + b) / 2.0; + bail!(RootError::NotConverge([c])); } } @@ -127,8 +157,8 @@ impl RootFindingMethod<1, 1, (f64, f64)> for BisectionMethod { /// /// # Arguments /// -/// - `max_iterations`: Maximum number of iterations -/// - `tolerance`: Tolerance +/// - `max_iter`: Maximum number of iterations +/// - `tol`: Absolute tolerance /// /// # Caution /// @@ -136,24 +166,55 @@ impl RootFindingMethod<1, 1, (f64, f64)> for BisectionMethod { /// - This method highly depends on the initial guess /// - This method is not guaranteed to converge pub struct NewtonMethod { - pub max_iterations: usize, - pub tolerance: f64, + pub max_iter: usize, + pub tol: f64, +} + +#[derive(Debug, Copy, Clone)] +pub enum NewtonError { + ZeroDerivative(f64), + NotConverge(f64), +} + +impl std::fmt::Display for NewtonError { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + NewtonError::ZeroDerivative(x) => { + write!(f, "Zero derivative at x = {}", x) + } + NewtonError::NotConverge(x) => { + write!(f, "Not converge at x = {}", x) + } + } + } } -impl RootFindingMethod<1, 1, f64> for NewtonMethod { - fn step>( +impl RootFinder<1, 1, f64> for NewtonMethod { + fn max_iter(&self) -> usize { + self.max_iter + } + fn tol(&self) -> f64 { + self.tol + } + fn find>( &self, problem: &P, - x: f64, - ) -> Result { - let f = single_function!(problem, x); - let df = single_derivative!(problem, x); + ) -> Result<[f64; 1]> { + let mut x = problem.initial_guess(); - if df == 0.0 { - bail!("Zero derivative at x = {}", x); - } + for _ in 0..self.max_iter { + let f = single_function!(problem, x); + let df = single_derivative!(problem, x); - Ok(x - f / df) + if f.abs() < self.tol { + return Ok([x]); + } else if df == 0.0 { + bail!(NewtonError::ZeroDerivative(x)); + } else { + x -= f / df; + } + } + bail!(NewtonError::NotConverge(x)); } } @@ -164,32 +225,141 @@ impl RootFindingMethod<1, 1, f64> for NewtonMethod { /// /// # Arguments /// -/// - `max_iterations`: Maximum number of iterations -/// - `tolerance`: Tolerance +/// - `max_iter`: Maximum number of iterations +/// - `tol`: Absolute tolerance /// /// # Caution /// /// - The function should be differentiable pub struct SecantMethod { - pub max_iterations: usize, - pub tolerance: f64, + pub max_iter: usize, + pub tol: f64, +} + +#[derive(Debug, Copy, Clone)] +pub enum SecantError { + ZeroSecant(f64, f64), + NotConverge(f64, f64), +} + +impl std::fmt::Display for SecantError { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + SecantError::ZeroSecant(x0, x1) => { + write!(f, "Zero secant at ({}, {})", x0, x1) + } + SecantError::NotConverge(x0, x1) => { + write!(f, "Not converge at ({}, {})", x0, x1) + } + } + } } -impl RootFindingMethod<1, 1, (f64, f64)> for SecantMethod { - fn step>( +impl RootFinder<1, 1, (f64, f64)> for SecantMethod { + fn max_iter(&self) -> usize { + self.max_iter + } + fn tol(&self) -> f64 { + self.tol + } + fn find>( &self, problem: &P, - state: (f64, f64), - ) -> Result<(f64, f64)> { + ) -> Result<[f64; 1]> { + let mut state = problem.initial_guess(); + + for _ in 0..self.max_iter { + let (x0, x1) = state; + let f0 = single_function!(problem, x0); + let f1 = single_function!(problem, x1); + + if f1.abs() < self.tol { + return Ok([x1]); + } + + if f0 == f1 { + bail!(SecantError::ZeroSecant(x0, x1)); + } + + state = (x1, x1 - f1 * (x1 - x0) / (f1 - f0)) + } let (x0, x1) = state; - let f0 = single_function!(problem, x0); - let f1 = single_function!(problem, x1); + bail!(SecantError::NotConverge(x0, x1)); + } +} + +// ┌─────────────────────────────────────────────────────────┐ +// False position method +// └─────────────────────────────────────────────────────────┘ +/// False position method +/// +/// # Arguments +/// +/// - `max_iter`: Maximum number of iterations +/// - `tol`: Absolute tolerance +/// +/// # Caution +/// +/// - The function should be differentiable +pub struct FalsePositionMethod { + pub max_iter: usize, + pub tol: f64, +} + +#[derive(Debug, Copy, Clone)] +pub enum FalsePositionError { + NoRoot, + NotConverge(f64, f64), +} + +impl std::fmt::Display for FalsePositionError { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + FalsePositionError::NoRoot => write!(f, "There is no root in the interval"), + FalsePositionError::NotConverge(a, b) => { + write!(f, "Not converge in [{}, {}]", a, b) + } + } + } +} + +impl RootFinder<1, 1, (f64, f64)> for FalsePositionMethod { + fn max_iter(&self) -> usize { + self.max_iter + } + fn tol(&self) -> f64 { + self.tol + } + fn find>( + &self, + problem: &P, + ) -> Result<[f64; 1]> { + let state = problem.initial_guess(); + let (mut a, mut b) = state; + let mut fa = single_function!(problem, a); + let mut fb = single_function!(problem, b); - if f0 == f1 { - bail!("Zero secant at ({}, {})", x0, x1); + if fa * fb > 0.0 { + bail!(FalsePositionError::NoRoot); } - Ok((x1, x1 - f1 * (x1 - x0) / (f1 - f0))) + for _ in 0..self.max_iter { + let c = (a * fb - b * fa) / (fb - fa); + let fc = single_function!(problem, c); + + if fc.abs() < self.tol { + return Ok([c]); + } else if fa * fc < 0.0 { + b = c; + fb = fc; + } else if fb * fc < 0.0 { + a = c; + fa = fc; + } else { + bail!(FalsePositionError::NoRoot); + } + } + bail!(FalsePositionError::NotConverge(a, b)); } } @@ -200,24 +370,29 @@ impl RootFindingMethod<1, 1, (f64, f64)> for SecantMethod { /// /// # Arguments /// -/// - `max_iterations`: Maximum number of iterations -/// - `tolerance`: Tolerance +/// - `max_iter`: Maximum number of iterations +/// - `tol`: Absolute tolerance /// /// # Caution /// /// - The function should be differentiable pub struct BroydenMethod { - pub max_iterations: usize, - pub tolerance: f64, + pub max_iter: usize, + pub tol: f64, } -impl RootFindingMethod> for BroydenMethod { - fn step>>( +impl RootFinder> for BroydenMethod { + fn max_iter(&self) -> usize { + self.max_iter + } + fn tol(&self) -> f64 { + self.tol + } + fn find>>( &self, problem: &P, - state: INTV, - ) -> Result> { + ) -> Result> { unimplemented!() } } From 59456287eb37f8ba726803a843be5cfc6f208ca4 Mon Sep 17 00:00:00 2001 From: axect Date: Sat, 13 Apr 2024 23:50:33 +0900 Subject: [PATCH 07/15] CHGE: Change test and error type for root finding --- examples/root_test.rs | 46 +++++++++ src/fuga/mod.rs | 1 - src/numerical/root.rs | 137 ++++++++++--------------- src/prelude/mod.rs | 4 +- tests/root.rs | 232 +++++++++++++++++++++--------------------- 5 files changed, 213 insertions(+), 207 deletions(-) create mode 100644 examples/root_test.rs diff --git a/examples/root_test.rs b/examples/root_test.rs new file mode 100644 index 00000000..b797a39d --- /dev/null +++ b/examples/root_test.rs @@ -0,0 +1,46 @@ +use peroxide::fuga::*; +use anyhow::Result; + +fn main() -> Result<()> { + let problem = Cubic; + let bisect = BisectionMethod { max_iter: 100, tol: 1e-6 }; + let newton = NewtonMethod { max_iter: 100, tol: 1e-6 }; + let false_pos = FalsePositionMethod { max_iter: 100, tol: 1e-6 }; + let result_bisect = bisect.find(&problem)?; + let result_newton = newton.find(&problem)?; + let result_false_pos = false_pos.find(&problem)?; + println!("{:?}", result_bisect); + println!("{:?}", result_newton); + println!("{:?}", result_false_pos); + + Ok(()) +} + +struct Cubic; + +impl Cubic { + fn eval(&self, x: Pt<1>) -> Result> { + Ok([(x[0] - 1f64).powi(3)]) + } +} + +impl RootFindingProblem<1, 1, (f64, f64)> for Cubic { + fn function(&self, x: Pt<1>) -> Result> { + self.eval(x) + } + fn initial_guess(&self) -> (f64, f64) { + (0.0, 2.0) + } +} + +impl RootFindingProblem<1, 1, f64> for Cubic { + fn function(&self, x: Pt<1>) -> Result> { + self.eval(x) + } + fn initial_guess(&self) -> f64 { + 0.0 + } + fn derivative(&self, x: Pt<1>) -> Result> { + Ok([[3.0 * (x[0] - 1f64).powi(2)]]) + } +} diff --git a/src/fuga/mod.rs b/src/fuga/mod.rs index f91274b6..006212b7 100644 --- a/src/fuga/mod.rs +++ b/src/fuga/mod.rs @@ -213,7 +213,6 @@ pub use crate::numerical::integral::Integral::{ G25K51R, G30K61R, }; -pub use crate::numerical::root::RootFind::{Bisection, FalsePosition, Newton, Secant}; pub use crate::statistics::stat::QType::{ Type1, Type2, Type3, Type4, Type5, Type6, Type7, Type8, Type9, }; diff --git a/src/numerical/root.rs b/src/numerical/root.rs index f782534b..6a4fb61a 100644 --- a/src/numerical/root.rs +++ b/src/numerical/root.rs @@ -1,9 +1,9 @@ use anyhow::{Result, bail}; -type PT = [f64; N]; -type INTV = (PT, PT); -type J = [[f64; C]; R]; -type H = [[[f64; C]; C]; R]; +pub type Pt = [f64; N]; +pub type Intv = (Pt, Pt); +pub type Jaco = [[f64; C]; R]; +pub type Hess = [[[f64; C]; C]; R]; /// Trait to define a root finding problem /// @@ -13,12 +13,14 @@ type H = [[[f64; C]; C]; R]; /// - `O`: Output type (e.g. `f64`, `[f64; N]`, or etc.) /// - `T`: State type (e.g. `f64`, `(f64, f64)`, or etc.) pub trait RootFindingProblem { - fn function(&self, x: PT) -> Result>; + fn function(&self, x: Pt) -> Result>; fn initial_guess(&self) -> T; - fn derivative(&self, x: PT) -> Result> { + #[allow(unused_variables)] + fn derivative(&self, x: Pt) -> Result> { unimplemented!() } - fn hessian(&self, x: PT) -> Result> { + #[allow(unused_variables)] + fn hessian(&self, x: Pt) -> Result> { unimplemented!() } } @@ -26,14 +28,15 @@ pub trait RootFindingProblem { pub trait RootFinder { fn max_iter(&self) -> usize; fn tol(&self) -> f64; - fn find>(&self, problem: &P) -> Result>; + fn find>(&self, problem: &P) -> Result>; } #[derive(Debug, Copy, Clone)] pub enum RootError { - NotConverge(PT), + NotConverge(Pt), NoRoot, - ZeroDerivative(PT), + ZeroDerivative(Pt), + ZeroSecant(Pt, Pt), } impl std::fmt::Display for RootError<1> { @@ -42,6 +45,7 @@ impl std::fmt::Display for RootError<1> { RootError::NoRoot => write!(f, "There is no root in the interval"), RootError::NotConverge(a) => write!(f, "Not yet converge. Our guess is {:?}", a), RootError::ZeroDerivative(a) => write!(f, "Zero derivative in {:?}", a), + RootError::ZeroSecant(a,b) => write!(f, "Zero secant in ({:?}, {:?})", a, b), } } } @@ -52,15 +56,16 @@ impl std::fmt::Display for RootError<1> { /// /// For I=1, O=1, it is bother to write below code. /// -/// ```no_run +/// ```ignore /// let fx = problem.function([x])?[0]; /// ``` /// /// This macro solve this problem as follows. /// -/// ```no_run +/// ```ignore /// let fx = single_function!(problem, x); /// ``` +#[macro_export] macro_rules! single_function { ($problem:expr, $x:expr) => {{ $problem.function([$x])?[0] @@ -73,15 +78,16 @@ macro_rules! single_function { /// /// For I=1, O=1, it is bother to write below code. /// -/// ```no_run +/// ```ignore /// let fx = problem.derivative([x])?[0][0]; /// ``` /// /// This macro solve this problem as follows. /// -/// ```no_run +/// ```ignore /// let fx = single_derivative!(problem, x); /// ``` +#[macro_export] macro_rules! single_derivative { ($problem:expr, $x:expr) => {{ $problem.derivative([$x])?[0][0] @@ -125,7 +131,11 @@ impl RootFinder<1, 1, (f64, f64)> for BisectionMethod { let mut fa = single_function!(problem, a); let mut fb = single_function!(problem, b); - if fa * fb > 0.0 { + if fa.abs() < self.tol { + return Ok([a]); + } else if fb.abs() < self.tol { + return Ok([b]); + } else if fa * fb > 0.0 { bail!(RootError::NoRoot); } @@ -170,25 +180,6 @@ pub struct NewtonMethod { pub tol: f64, } -#[derive(Debug, Copy, Clone)] -pub enum NewtonError { - ZeroDerivative(f64), - NotConverge(f64), -} - -impl std::fmt::Display for NewtonError { - fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { - match self { - NewtonError::ZeroDerivative(x) => { - write!(f, "Zero derivative at x = {}", x) - } - NewtonError::NotConverge(x) => { - write!(f, "Not converge at x = {}", x) - } - } - } -} - impl RootFinder<1, 1, f64> for NewtonMethod { fn max_iter(&self) -> usize { self.max_iter @@ -209,12 +200,12 @@ impl RootFinder<1, 1, f64> for NewtonMethod { if f.abs() < self.tol { return Ok([x]); } else if df == 0.0 { - bail!(NewtonError::ZeroDerivative(x)); + bail!(RootError::ZeroDerivative([x])); } else { x -= f / df; } } - bail!(NewtonError::NotConverge(x)); + bail!(RootError::NotConverge([x])); } } @@ -236,25 +227,6 @@ pub struct SecantMethod { pub tol: f64, } -#[derive(Debug, Copy, Clone)] -pub enum SecantError { - ZeroSecant(f64, f64), - NotConverge(f64, f64), -} - -impl std::fmt::Display for SecantError { - fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { - match self { - SecantError::ZeroSecant(x0, x1) => { - write!(f, "Zero secant at ({}, {})", x0, x1) - } - SecantError::NotConverge(x0, x1) => { - write!(f, "Not converge at ({}, {})", x0, x1) - } - } - } -} - impl RootFinder<1, 1, (f64, f64)> for SecantMethod { fn max_iter(&self) -> usize { self.max_iter @@ -266,11 +238,15 @@ impl RootFinder<1, 1, (f64, f64)> for SecantMethod { &self, problem: &P, ) -> Result<[f64; 1]> { - let mut state = problem.initial_guess(); + let state = problem.initial_guess(); + let (mut x0, mut x1) = state; + let mut f0 = single_function!(problem, x0); + + if f0.abs() < self.tol { + return Ok([x0]); + } for _ in 0..self.max_iter { - let (x0, x1) = state; - let f0 = single_function!(problem, x0); let f1 = single_function!(problem, x1); if f1.abs() < self.tol { @@ -278,13 +254,13 @@ impl RootFinder<1, 1, (f64, f64)> for SecantMethod { } if f0 == f1 { - bail!(SecantError::ZeroSecant(x0, x1)); + bail!(RootError::ZeroSecant([x0], [x1])); } - state = (x1, x1 - f1 * (x1 - x0) / (f1 - f0)) + f0 = f1; + (x0, x1) = (x1, x1 - f1 * (x1 - x0) / (f1 - f0)) } - let (x0, x1) = state; - bail!(SecantError::NotConverge(x0, x1)); + bail!(RootError::NotConverge([x1])); } } @@ -306,23 +282,6 @@ pub struct FalsePositionMethod { pub tol: f64, } -#[derive(Debug, Copy, Clone)] -pub enum FalsePositionError { - NoRoot, - NotConverge(f64, f64), -} - -impl std::fmt::Display for FalsePositionError { - fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { - match self { - FalsePositionError::NoRoot => write!(f, "There is no root in the interval"), - FalsePositionError::NotConverge(a, b) => { - write!(f, "Not converge in [{}, {}]", a, b) - } - } - } -} - impl RootFinder<1, 1, (f64, f64)> for FalsePositionMethod { fn max_iter(&self) -> usize { self.max_iter @@ -339,8 +298,12 @@ impl RootFinder<1, 1, (f64, f64)> for FalsePositionMethod { let mut fa = single_function!(problem, a); let mut fb = single_function!(problem, b); - if fa * fb > 0.0 { - bail!(FalsePositionError::NoRoot); + if fa.abs() < self.tol { + return Ok([a]); + } else if fb.abs() < self.tol { + return Ok([b]); + } else if fa * fb > 0.0 { + bail!(RootError::NoRoot); } for _ in 0..self.max_iter { @@ -356,10 +319,11 @@ impl RootFinder<1, 1, (f64, f64)> for FalsePositionMethod { a = c; fa = fc; } else { - bail!(FalsePositionError::NoRoot); + bail!(RootError::NoRoot); } } - bail!(FalsePositionError::NotConverge(a, b)); + let c = (a * fb - b * fa) / (fb - fa); + bail!(RootError::NotConverge([c])); } } @@ -382,17 +346,18 @@ pub struct BroydenMethod { } -impl RootFinder> for BroydenMethod { +#[allow(unused_variables)] +impl RootFinder> for BroydenMethod { fn max_iter(&self) -> usize { self.max_iter } fn tol(&self) -> f64 { self.tol } - fn find>>( + fn find>>( &self, problem: &P, - ) -> Result> { + ) -> Result> { unimplemented!() } } diff --git a/src/prelude/mod.rs b/src/prelude/mod.rs index 7c0d8a39..77ef606e 100644 --- a/src/prelude/mod.rs +++ b/src/prelude/mod.rs @@ -201,7 +201,7 @@ pub use crate::numerical::{ interp::*, ode::*, optimize::*, - root::{bisection, false_position, newton, secant}, + root::*, spline::{cubic_spline, CubicSpline, CubicHermiteSpline, Spline}, utils::*, }; @@ -215,4 +215,4 @@ pub use crate::statistics::stat::Metric::*; pub use simpler::SimpleParquet; #[cfg(feature="plot")] -pub use crate::util::plot::*; \ No newline at end of file +pub use crate::util::plot::*; diff --git a/tests/root.rs b/tests/root.rs index ec2e7b39..b0d8128f 100644 --- a/tests/root.rs +++ b/tests/root.rs @@ -1,146 +1,142 @@ -// #[macro_use] -extern crate peroxide; use peroxide::fuga::*; -// use peroxide::numerical::root::*; +use anyhow::Result; #[test] -fn test_root_finder() -> Result<(), RootError> { - let init = RootState::I(0.1f64, 5f64); - let init_p = RootState::P(0.1f64); - let init_s = RootState::I(0.1f64, 0.2f64); - // Bisection - let mut a1 = RootFinder::new(init, RootFind::Bisection, f_exp)?; - let mut b1 = RootFinder::new(init, RootFind::Bisection, f_ln)?; - let mut c1 = RootFinder::new(init, RootFind::Bisection, f_sqrt)?; - a1.set_tol(1e-15); - b1.set_tol(1e-15); - c1.set_tol(1e-15); - // FalsePosition - let mut a2 = RootFinder::new(init, RootFind::FalsePosition, f_exp)?; - let mut b2 = RootFinder::new(init, RootFind::FalsePosition, f_ln)?; - let mut c2 = RootFinder::new(init, RootFind::FalsePosition, f_sqrt)?; - a2.set_times(10000); - b2.set_times(10000); - c2.set_times(10000); - // Secant - let mut a3 = RootFinder::new(init_s, RootFind::Secant, f_exp)?; - let mut b3 = RootFinder::new(init_s, RootFind::Secant, f_ln)?; - let mut c3 = RootFinder::new(init_s, RootFind::Secant, f_sqrt)?; - a3.set_times(10000); - b3.set_times(10000); - c3.set_times(10000); - // Newton - let mut a4 = RootFinder::new(init_p, RootFind::Newton, f_exp)?; - let mut b4 = RootFinder::new(init_p, RootFind::Newton, f_ln)?; - let mut c4 = RootFinder::new(init_p, RootFind::Newton, f_sqrt)?; - a4.set_tol(1e-15); - b4.set_tol(1e-15); - c4.set_tol(1e-15); - - let x1 = a1.find_root()?; - let x2 = b1.find_root()?; - let x3 = c1.find_root()?; - let y1 = a2.find_root()?; - let y2 = b2.find_root()?; - let y3 = c2.find_root()?; - let z1 = a3.find_root()?; - let z2 = b3.find_root()?; - let z3 = c3.find_root()?; - let w1 = a4.find_root()?; - let w2 = b4.find_root()?; - let w3 = c4.find_root()?; - x1.print(); - x2.print(); - x3.print(); - y1.print(); - y2.print(); - y3.print(); - z1.print(); - z2.print(); - z3.print(); - w1.print(); - w2.print(); - w3.print(); +fn test_cubic_root() -> Result<()> { + let problem = Cubic; + let bisect = BisectionMethod { max_iter: 100, tol: 1e-6 }; + let newton = NewtonMethod { max_iter: 100, tol: 1e-6 }; + let false_pos = FalsePositionMethod { max_iter: 100, tol: 1e-6 }; + let root_bisect = bisect.find(&problem)?; + let root_newton = newton.find(&problem)?; + let root_false_pos = false_pos.find(&problem)?; + + let result_bisect = problem.eval(root_bisect)?[0]; + let result_newton = problem.eval(root_newton)?[0]; + let result_false_pos = problem.eval(root_false_pos)?[0]; + + assert!(result_bisect.abs() < 1e-6); + assert!(result_newton.abs() < 1e-6); + assert!(result_false_pos.abs() < 1e-6); Ok(()) } -#[test] -fn test_bisection() -> Result<(), RootError> { - let x1 = bisection(f_exp, (0f64, 5f64), 100, 1e-15)?; - let x2 = bisection(f_ln, (0.1f64, 5f64), 100, 1e-15)?; - let x3 = bisection(f_sqrt, (0f64, 5f64), 100, 1e-15)?; - let x4 = bisection(f_sin, (0f64, 5f64), 100, 1e-15)?; - - println!("Bisection: "); - x1.print(); - x2.print(); - x3.print(); - x4.print(); +struct Cubic; - Ok(()) +impl Cubic { + fn eval(&self, x: [f64; 1]) -> Result<[f64; 1]> { + Ok([(x[0] - 1f64).powi(3)]) + } } -#[test] -fn test_secant() -> Result<(), RootError> { - let x1 = secant(f_exp, (0f64, 0.1f64), 200, 1e-15)?; - let x2 = secant(f_ln, (0.1f64, 0.2f64), 200, 1e-15)?; - let x3 = secant(f_sqrt, (0f64, 0.1f64), 200, 1e-15)?; - let x4 = secant(f_sin, (2f64, 2.1f64), 200, 1e-15)?; - - println!("Secant: "); - x1.print(); - x2.print(); - x3.print(); - x4.print(); +impl RootFindingProblem<1, 1, (f64, f64)> for Cubic { + fn function(&self, x: [f64; 1]) -> Result<[f64; 1]> { + self.eval(x) + } + fn initial_guess(&self) -> (f64, f64) { + (0.0, 2.0) + } +} - Ok(()) +impl RootFindingProblem<1, 1, f64> for Cubic { + fn function(&self, x: [f64; 1]) -> Result<[f64; 1]> { + self.eval(x) + } + fn initial_guess(&self) -> f64 { + 0.0 + } + fn derivative(&self, x: [f64; 1]) -> Result> { + Ok([[3.0 * (x[0] - 1f64).powi(2)]]) + } } #[test] -fn test_false_position() -> Result<(), RootError> { - println!("False Position: "); - let x1 = false_position(f_exp, (0f64, 5f64), 1000, 1e-10)?; - let x2 = false_position(f_ln, (0.1f64, 5f64), 1000, 1e-10)?; - let x3 = false_position(f_sqrt, (0f64, 5f64), 1000, 1e-10)?; - let x4 = false_position(f_sin, (1f64, 5f64), 1000, 1e-10)?; - - x1.print(); - x2.print(); - x3.print(); - x4.print(); +fn test_sine_root() -> Result<()> { + let problem = Sine; + let bisect = BisectionMethod { max_iter: 100, tol: 1e-6 }; + let newton = NewtonMethod { max_iter: 100, tol: 1e-6 }; + let false_pos = FalsePositionMethod { max_iter: 100, tol: 1e-6 }; + let root_bisect = bisect.find(&problem)?; + let root_newton = newton.find(&problem)?; + let root_false_pos = false_pos.find(&problem)?; + + let result_bisect = problem.eval(root_bisect)?[0]; + let result_newton = problem.eval(root_newton)?[0]; + let result_false_pos = problem.eval(root_false_pos)?[0]; + + assert!(result_bisect.abs() < 1e-6); + assert!(result_newton.abs() < 1e-6); + assert!(result_false_pos.abs() < 1e-6); Ok(()) } -#[test] -fn test_newton() -> Result<(), RootError> { - let x1 = newton(f_exp, 0.1f64, 100, 1e-15)?; - let x2 = newton(f_ln, 0.1f64, 100, 1e-15)?; - let x3 = newton(f_sqrt, 0.1f64, 100, 1e-15)?; - let x4 = newton(f_sin, 2f64, 100, 1e-15)?; - - println!("Newton: "); - x1.print(); - x2.print(); - x3.print(); - x4.print(); +struct Sine; - Ok(()) +impl Sine { + fn eval(&self, x: [f64; 1]) -> Result<[f64; 1]> { + Ok([x[0].sin()]) + } } -fn f_exp(x: AD) -> AD { - x.exp() - 2f64 +impl RootFindingProblem<1, 1, (f64, f64)> for Sine { + fn function(&self, x: [f64; 1]) -> Result<[f64; 1]> { + self.eval(x) + } + fn initial_guess(&self) -> (f64, f64) { + (0.0, 2.0) + } } -fn f_ln(x: AD) -> AD { - x.ln() +impl RootFindingProblem<1, 1, f64> for Sine { + fn function(&self, x: [f64; 1]) -> Result<[f64; 1]> { + self.eval(x) + } + fn initial_guess(&self) -> f64 { + 1.0 + } + fn derivative(&self, x: [f64; 1]) -> Result> { + Ok([[x[0].cos()]]) + } } -fn f_sqrt(x: AD) -> AD { - x.sqrt() - 2f64 +#[test] +fn test_cosine_root() -> Result<()> { + let problem = Cosine; + let newton = NewtonMethod { max_iter: 100, tol: 1e-6 }; + let root_newton = match newton.find(&problem) { + Ok(x) => x, + Err(e) => { + println!("{:?}", e); + match e.downcast::>() { + Ok(RootError::ZeroDerivative(x)) => x, + Ok(e) => panic!("ok but {:?}", e), + Err(e) => panic!("err {:?}", e), + } + } + }; + assert_eq!(root_newton[0], 0.0); + + Ok(()) +} + +struct Cosine; + +impl Cosine { + fn eval(&self, x: [f64; 1]) -> Result<[f64; 1]> { + Ok([x[0].cos()]) + } } -fn f_sin(x: AD) -> AD { - x.sin() +impl RootFindingProblem<1, 1, f64> for Cosine { + fn function(&self, x: [f64; 1]) -> Result<[f64; 1]> { + self.eval(x) + } + fn initial_guess(&self) -> f64 { + 0.0 // should fail in newton (derivative is 0) + } + fn derivative(&self, x: [f64; 1]) -> Result> { + Ok([[-x[0].sin()]]) + } } From a702e6293f0c24c89e1f0a7628388aec5afb9d95 Mon Sep 17 00:00:00 2001 From: axect Date: Sat, 13 Apr 2024 23:58:36 +0900 Subject: [PATCH 08/15] DOCS: Add explanation for type parameters in root.rs --- src/numerical/root.rs | 64 ++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 63 insertions(+), 1 deletion(-) diff --git a/src/numerical/root.rs b/src/numerical/root.rs index 6a4fb61a..95fbf682 100644 --- a/src/numerical/root.rs +++ b/src/numerical/root.rs @@ -1,8 +1,12 @@ use anyhow::{Result, bail}; +/// Point alias (`[f64; N]`) pub type Pt = [f64; N]; +/// Interval alias (`([f64; N], [f64; N])`) pub type Intv = (Pt, Pt); +/// Jacobian alias (`[[f64; C]; R]`) pub type Jaco = [[f64; C]; R]; +/// Hessian alias (`[[[f64; C]; C]; R]`) pub type Hess = [[[f64; C]; C]; R]; /// Trait to define a root finding problem @@ -12,6 +16,13 @@ pub type Hess = [[[f64; C]; C]; R]; /// - `I`: Input type (e.g. `f64`, `[f64; N]`, or etc.) /// - `O`: Output type (e.g. `f64`, `[f64; N]`, or etc.) /// - `T`: State type (e.g. `f64`, `(f64, f64)`, or etc.) +/// +/// # Methods +/// +/// - `function`: Function +/// - `initial_guess`: Initial guess +/// - `derivative`: Derivative (optional) +/// - `hessian`: Hessian (optional) pub trait RootFindingProblem { fn function(&self, x: Pt) -> Result>; fn initial_guess(&self) -> T; @@ -25,6 +36,26 @@ pub trait RootFindingProblem { } } +/// Trait to define a root finder +/// +/// # Type Parameters +/// +/// - `I`: Input type (e.g. `f64`, `[f64; N]`, or etc.) +/// - `O`: Output type (e.g. `f64`, `[f64; N]`, or etc.) +/// - `T`: State type (e.g. `f64`, `(f64, f64)`, or etc.) +/// +/// # Methods +/// +/// - `max_iter`: Maximum number of iterations +/// - `tol`: Absolute tolerance +/// - `find`: Find root +/// +/// # Available root finders +/// +/// - `BisectionMethod`: `I=1, O=1, T=(f64, f64)` +/// - `FalsePositionMethod`: `I=1, O=1, T=(f64, f64)` +/// - `NewtonMethod`: `I=1, O=1, T=f64` +/// - `SecantMethod`: `I=1, O=1, T=(f64, f64)` pub trait RootFinder { fn max_iter(&self) -> usize; fn tol(&self) -> f64; @@ -99,6 +130,12 @@ macro_rules! single_derivative { // └─────────────────────────────────────────────────────────┘ /// Bisection method /// +/// # Type for `RootFinder` +/// +/// - `I`: 1 +/// - `O`: 1 +/// - `T`: `(f64, f64)` +/// /// # Arguments /// /// - `max_iter`: Maximum number of iterations @@ -165,6 +202,12 @@ impl RootFinder<1, 1, (f64, f64)> for BisectionMethod { // └─────────────────────────────────────────────────────────┘ /// Newton method /// +/// # Type for `RootFinder` +/// +/// - `I`: 1 +/// - `O`: 1 +/// - `T`: `f64` +/// /// # Arguments /// /// - `max_iter`: Maximum number of iterations @@ -214,6 +257,12 @@ impl RootFinder<1, 1, f64> for NewtonMethod { // └─────────────────────────────────────────────────────────┘ /// Secant method /// +/// # Type for `RootFinder` +/// +/// - `I`: 1 +/// - `O`: 1 +/// - `T`: `(f64, f64)` +/// /// # Arguments /// /// - `max_iter`: Maximum number of iterations @@ -222,6 +271,7 @@ impl RootFinder<1, 1, f64> for NewtonMethod { /// # Caution /// /// - The function should be differentiable +/// - This method is not guaranteed to converge pub struct SecantMethod { pub max_iter: usize, pub tol: f64, @@ -269,6 +319,12 @@ impl RootFinder<1, 1, (f64, f64)> for SecantMethod { // └─────────────────────────────────────────────────────────┘ /// False position method /// +/// # Type for `RootFinder` +/// +/// - `I`: 1 +/// - `O`: 1 +/// - `T`: `(f64, f64)` +/// /// # Arguments /// /// - `max_iter`: Maximum number of iterations @@ -276,7 +332,7 @@ impl RootFinder<1, 1, (f64, f64)> for SecantMethod { /// /// # Caution /// -/// - The function should be differentiable +/// - The function should be continuous pub struct FalsePositionMethod { pub max_iter: usize, pub tol: f64, @@ -332,6 +388,12 @@ impl RootFinder<1, 1, (f64, f64)> for FalsePositionMethod { // └─────────────────────────────────────────────────────────┘ /// Broyden method /// +/// # Type for `RootFinder` +/// +/// - `I`: free +/// - `O`: free +/// - `T`: `Intv` (=`([f64; I], [f64; I])`) +/// /// # Arguments /// /// - `max_iter`: Maximum number of iterations From 02a1eae3238944f79affca2403c1fa0a374d45bd Mon Sep 17 00:00:00 2001 From: axect Date: Sun, 14 Apr 2024 00:04:08 +0900 Subject: [PATCH 09/15] DOCS: Update README for root finding --- README.md | 20 ++++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index 2f0ba09d..a9c01ce1 100644 --- a/README.md +++ b/README.md @@ -177,9 +177,19 @@ Peroxide can do many things. - Gradient Descent - Levenberg Marquardt - Ordinary Differential Equation - - Runge-Kutta 4th order (`RK4`) - - Runge-Kutta-Fehlberg 4/5th order (`RKF45`) - - Gauss-Legendre 4th order (`GL4`; Implicit) + - Trait based ODE solver (after `v0.36.0`) + - Explicit integrator + - Ralston's 3rd order + - Runge-Kutta 4th order + - Ralston's 4th order + - Runge-Kutta 5th order + - Embedded integrator + - Bogacki-Shampine 3(2) + - Runge-Kutta-Fehlberg 4(5) + - Dormand-Prince 5(4) + - Tsitouras 5(4) + - Implicit integrator + - Gauss-Legendre 4th order - Numerical Integration - Newton-Cotes Quadrature - Gauss-Legendre Quadrature (up to 30 order) @@ -188,8 +198,9 @@ Peroxide can do many things. - Gauss-Kronrod Quadrature (Relative tolerance) - G7K15R, G10K21R, G15K31R, G20K41R, G25K51R, G30K61R - Root Finding + - Trait based root finding (after `v0.37.0`) - Bisection - - False Position (Regula Falsi) + - False Position - Secant - Newton - Statistics @@ -211,6 +222,7 @@ Peroxide can do many things. - Marsaglia Polar - Ziggurat - Wrapper for `rand-dist` crate + - Piecewise Rejection Sampling - Confusion Matrix & Metrics - Special functions - Wrapper for `puruspe` crate (pure rust) From 990baf497c6700e164f36e520bd7ce83530be462 Mon Sep 17 00:00:00 2001 From: axect Date: Sun, 14 Apr 2024 00:17:12 +0900 Subject: [PATCH 10/15] CHGE: Change thiserror -> anyhow in ODE --- Cargo.toml | 3 +-- examples/lorenz_dp45.rs | 2 +- examples/lorenz_gl4.rs | 2 +- examples/lorenz_rk4.rs | 2 +- examples/lorenz_rkf45.rs | 2 +- examples/lorenz_tsit45.rs | 2 +- examples/ode_env.rs | 2 +- examples/ode_test_gl4.rs | 2 +- examples/ode_test_rk4.rs | 2 +- examples/ode_test_rkf45.rs | 2 +- src/fuga/mod.rs | 2 ++ src/numerical/ode.rs | 39 ++++++++++++++++++++++---------------- src/prelude/mod.rs | 2 ++ 13 files changed, 37 insertions(+), 27 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 7f301446..772dca9e 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -28,7 +28,7 @@ puruspe = "0.2" matrixmultiply = { version = "0.3", features = ["threading"] } peroxide-ad = "0.3" peroxide-num = "0.1" -thiserror = "1.0" +anyhow = "1.0" #num-complex = "0.3" netcdf = { version = "0.7", optional = true, default_features = false } pyo3 = { version = "0.21", optional = true, features = ["auto-initialize", "gil-refs"] } @@ -37,7 +37,6 @@ 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 } -anyhow = "1.0.82" [package.metadata.docs.rs] rustdoc-args = [ "--html-in-header", "katex-header.html", "--cfg", "docsrs"] diff --git a/examples/lorenz_dp45.rs b/examples/lorenz_dp45.rs index 9f3e5624..da1de787 100644 --- a/examples/lorenz_dp45.rs +++ b/examples/lorenz_dp45.rs @@ -38,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]) -> anyhow::Result<()> { 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_gl4.rs b/examples/lorenz_gl4.rs index 7b4477d0..ba17bcc4 100644 --- a/examples/lorenz_gl4.rs +++ b/examples/lorenz_gl4.rs @@ -38,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]) -> anyhow::Result<()> { 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_rk4.rs b/examples/lorenz_rk4.rs index c4e51a28..79fc3b9a 100644 --- a/examples/lorenz_rk4.rs +++ b/examples/lorenz_rk4.rs @@ -37,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]) -> anyhow::Result<()> { 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 2184049c..98c3ca65 100644 --- a/examples/lorenz_rkf45.rs +++ b/examples/lorenz_rkf45.rs @@ -38,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]) -> anyhow::Result<()> { 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 index 748a732a..f923e1c5 100644 --- a/examples/lorenz_tsit45.rs +++ b/examples/lorenz_tsit45.rs @@ -38,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]) -> anyhow::Result<()> { 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/ode_env.rs b/examples/ode_env.rs index 58985adc..e627e8bd 100644 --- a/examples/ode_env.rs +++ b/examples/ode_env.rs @@ -49,7 +49,7 @@ impl ODEProblem for Test { vec![1f64] } - fn rhs(&self, t: f64, _y: &[f64], dy: &mut [f64]) -> Result<(), ODEError> { + fn rhs(&self, t: f64, _y: &[f64], dy: &mut [f64]) -> anyhow::Result<()> { Ok(dy[0] = self.cs.eval(t)) } } diff --git a/examples/ode_test_gl4.rs b/examples/ode_test_gl4.rs index a1e388f7..1d1791e5 100644 --- a/examples/ode_test_gl4.rs +++ b/examples/ode_test_gl4.rs @@ -36,7 +36,7 @@ impl ODEProblem for Test { vec![1f64] } - fn rhs(&self, t: f64, y: &[f64], dy: &mut [f64]) -> Result<(), ODEError> { + fn rhs(&self, t: f64, y: &[f64], dy: &mut [f64]) -> anyhow::Result<()> { Ok(dy[0] = (5f64 * t.powi(2) - y[0]) / (t + y[0]).exp()) } } diff --git a/examples/ode_test_rk4.rs b/examples/ode_test_rk4.rs index 879c07a3..a531aa3d 100644 --- a/examples/ode_test_rk4.rs +++ b/examples/ode_test_rk4.rs @@ -35,7 +35,7 @@ impl ODEProblem for Test { vec![1f64] } - fn rhs(&self, t: f64, y: &[f64], dy: &mut [f64]) -> Result<(), ODEError> { + fn rhs(&self, t: f64, y: &[f64], dy: &mut [f64]) -> anyhow::Result<()> { Ok(dy[0] = (5f64 * t.powi(2) - y[0]) / (t + y[0]).exp()) } } diff --git a/examples/ode_test_rkf45.rs b/examples/ode_test_rkf45.rs index 8c13e148..7441c3a1 100644 --- a/examples/ode_test_rkf45.rs +++ b/examples/ode_test_rkf45.rs @@ -36,7 +36,7 @@ impl ODEProblem for Test { vec![1f64] } - fn rhs(&self, t: f64, y: &[f64], dy: &mut [f64]) -> Result<(), ODEError> { + fn rhs(&self, t: f64, y: &[f64], dy: &mut [f64]) -> anyhow::Result<()> { Ok(dy[0] = (5f64 * t.powi(2) - y[0]) / (t + y[0]).exp()) } } diff --git a/src/fuga/mod.rs b/src/fuga/mod.rs index 006212b7..9fb97357 100644 --- a/src/fuga/mod.rs +++ b/src/fuga/mod.rs @@ -194,6 +194,8 @@ pub use crate::ml::reg::*; #[cfg(feature = "plot")] pub use crate::util::plot::*; +pub use anyhow; + // ============================================================================= // Enums // ============================================================================= diff --git a/src/numerical/ode.rs b/src/numerical/ode.rs index f961e5cf..933ae548 100644 --- a/src/numerical/ode.rs +++ b/src/numerical/ode.rs @@ -79,8 +79,7 @@ //! } //! ``` -use thiserror::Error; -use ODEError::*; +use anyhow::{Result, bail}; /// Trait for defining an ODE problem. /// @@ -107,7 +106,7 @@ use ODEError::*; /// ``` pub trait ODEProblem { fn initial_conditions(&self) -> Vec; - fn rhs(&self, t: f64, y: &[f64], dy: &mut [f64]) -> Result<(), ODEError>; + fn rhs(&self, t: f64, y: &[f64], dy: &mut [f64]) -> Result<()>; } @@ -115,7 +114,7 @@ pub trait ODEProblem { /// /// Implement this trait to define your own ODE integrator. pub trait ODEIntegrator { - fn step(&self, problem: &P, t: f64, y: &mut [f64], dt: f64) -> Result; + fn step(&self, problem: &P, t: f64, y: &mut [f64], dt: f64) -> Result; } @@ -132,6 +131,7 @@ pub trait ODEIntegrator { /// /// ```no_run /// use peroxide::fuga::*; +/// use anyhow::{Result, bail}; /// /// struct ConstrainedProblem { /// y_constraint: f64 @@ -139,9 +139,9 @@ pub trait ODEIntegrator { /// /// impl ODEProblem for ConstrainedProblem { /// fn initial_conditions(&self) -> Vec { vec![0.0] } // y_0 = 0 -/// fn rhs(&self, t: f64, y: &[f64], dy: &mut [f64]) -> Result<(), ODEError> { +/// fn rhs(&self, t: f64, y: &[f64], dy: &mut [f64]) -> Result<()> { /// if y[0] < self.y_constraint { -/// return Err(ODEError::ConstraintViolation(t, y.to_vec(), dy.to_vec())); +/// bail!(ODEError::ConstraintViolation(t, y.to_vec(), dy.to_vec())); /// } else { /// // some function /// Ok(()) @@ -149,19 +149,26 @@ pub trait ODEIntegrator { /// } /// } /// ``` -#[derive(Debug, Clone, Error)] +#[derive(Debug, Clone)] pub enum ODEError { - #[error("constraint violation")] ConstraintViolation(f64, Vec, Vec), // t, y, dy - #[error("reached maximum number of iterations per step")] ReachedMaxStepIter, } +impl std::fmt::Display for ODEError { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + ODEError::ConstraintViolation(t, y, dy) => write!(f, "Constraint violation at t = {}, y = {:?}, dy = {:?}", t, y, dy), + ODEError::ReachedMaxStepIter => write!(f, "Reached maximum number of steps per step"), + } + } +} + /// Trait for ODE solvers. /// /// Implement this trait to define your own ODE solver. pub trait ODESolver { - fn solve(&self, problem: &P, t_span: (f64, f64), dt: f64) -> Result<(Vec, Vec>), ODEError>; + fn solve(&self, problem: &P, t_span: (f64, f64), dt: f64) -> Result<(Vec, Vec>)>; } /// A basic ODE solver using a specified integrator. @@ -171,7 +178,7 @@ pub trait ODESolver { /// ``` /// use peroxide::fuga::*; /// -/// fn main() -> Result<(), Box> { +/// fn main() -> std::results::Result<(), Box> { /// let rkf = RKF45::new(1e-4, 0.9, 1e-6, 1e-1, 100); /// let basic_ode_solver = BasicODESolver::new(rkf); /// let (t_vec, y_vec) = basic_ode_solver.solve( @@ -191,7 +198,7 @@ pub trait ODESolver { /// vec![1f64] /// } /// -/// fn rhs(&self, t: f64, y: &[f64], dy: &mut [f64]) -> Result<(), ODEError> { +/// fn rhs(&self, t: f64, y: &[f64], dy: &mut [f64]) -> anyhow::Result<()> { /// dy[0] = (5f64 * t.powi(2) - y[0]) / (t + y[0]).exp(); /// Ok(()) /// } @@ -208,7 +215,7 @@ impl BasicODESolver { } impl ODESolver for BasicODESolver { - fn solve(&self, problem: &P, t_span: (f64, f64), dt: f64) -> Result<(Vec, Vec>), ODEError> { + fn solve(&self, problem: &P, t_span: (f64, f64), dt: f64) -> Result<(Vec, Vec>)> { let mut t = t_span.0; let mut dt = dt; let mut y = problem.initial_conditions(); @@ -271,7 +278,7 @@ pub trait ButcherTableau { } impl ODEIntegrator for BU { - fn step(&self, problem: &P, t: f64, y: &mut [f64], dt: f64) -> Result { + 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; @@ -320,7 +327,7 @@ impl ODEIntegrator for BU { } else { iter_count += 1; if iter_count >= self.max_step_iter() { - return Err(ReachedMaxStepIter); + bail!(ODEError::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()); @@ -742,7 +749,7 @@ impl GL4 { impl ODEIntegrator for GL4 { #[inline] - fn step(&self, problem: &P, t: f64, y: &mut [f64], dt: f64) -> Result { + fn step(&self, problem: &P, t: f64, y: &mut [f64], dt: f64) -> Result { let n = y.len(); let sqrt3 = 3.0_f64.sqrt(); let c = 0.5 * (3.0 - sqrt3) / 6.0; diff --git a/src/prelude/mod.rs b/src/prelude/mod.rs index 77ef606e..2a5b9f22 100644 --- a/src/prelude/mod.rs +++ b/src/prelude/mod.rs @@ -216,3 +216,5 @@ pub use simpler::SimpleParquet; #[cfg(feature="plot")] pub use crate::util::plot::*; + +pub use anyhow; From 70d3a1e8f52ccb31a6c8f64ceb51896b87610d5f Mon Sep 17 00:00:00 2001 From: axect Date: Sun, 14 Apr 2024 00:22:07 +0900 Subject: [PATCH 11/15] CHGE: Change thiserror -> anyhow in spline --- src/numerical/spline.rs | 53 +++++++++++++++++++++++------------------ 1 file changed, 30 insertions(+), 23 deletions(-) diff --git a/src/numerical/spline.rs b/src/numerical/spline.rs index c32f9fc6..85fe4a4e 100644 --- a/src/numerical/spline.rs +++ b/src/numerical/spline.rs @@ -214,7 +214,7 @@ use serde::{Deserialize, Serialize}; use std::cmp::{max, min}; use std::convert::From; use std::ops::{Index, Range}; -use thiserror::Error; +use anyhow::{Result, bail}; /// Trait for spline interpolation /// @@ -317,7 +317,7 @@ pub trait Spline { /// Ok(()) /// } /// ``` -pub fn cubic_spline(node_x: &[f64], node_y: &[f64]) -> Result { +pub fn cubic_spline(node_x: &[f64], node_y: &[f64]) -> Result { CubicSpline::from_nodes(node_x, node_y) } @@ -325,7 +325,7 @@ pub fn cubic_hermite_spline( node_x: &[f64], node_y: &[f64], slope_method: SlopeMethod, -) -> Result { +) -> Result { CubicHermiteSpline::from_nodes(node_x, node_y, slope_method) } @@ -383,18 +383,25 @@ impl Spline for CubicSpline { } } -#[derive(Debug, Copy, Clone, Error)] +#[derive(Debug, Copy, Clone)] pub enum SplineError { - #[error("node_x has less than 3 elements")] NotEnoughNodes, - #[error("node_x and node_y have different lengths")] NotEqualNodes, - #[error("nodes and slopes have different lengths")] NotEqualSlopes, - #[error("there are redundant nodes in node_x")] RedundantNodeX, } +impl std::fmt::Display for SplineError { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + SplineError::NotEnoughNodes => write!(f, "node_x has less than 3 elements"), + SplineError::NotEqualNodes => write!(f, "node_x and node_y have different lengths"), + SplineError::NotEqualSlopes => write!(f, "nodes and slopes have different lengths"), + SplineError::RedundantNodeX => write!(f, "there are redundant nodes in node_x"), + } + } +} + impl CubicSpline { /// # Examples /// ``` @@ -415,14 +422,14 @@ impl CubicSpline { /// Ok(()) /// } /// ``` - pub fn from_nodes(node_x: &[f64], node_y: &[f64]) -> Result { + pub fn from_nodes(node_x: &[f64], node_y: &[f64]) -> Result { let polynomials = CubicSpline::cubic_spline(node_x, node_y)?; Ok(CubicSpline { polynomials: zip_range(node_x, &polynomials), }) } - fn cubic_spline(node_x: &[f64], node_y: &[f64]) -> Result, SplineError> { + fn cubic_spline(node_x: &[f64], node_y: &[f64]) -> Result> { //! Pre calculated variables //! node_x: n+1 //! node_y: n+1 @@ -433,10 +440,10 @@ impl CubicSpline { //! z : n+1 let n = node_x.len() - 1; if n < 2 { - return Err(NotEnoughNodes); + bail!(NotEnoughNodes); } if n != node_y.len() - 1 { - return Err(NotEqualNodes); + bail!(NotEqualNodes); } // Pre-calculations @@ -504,7 +511,7 @@ impl CubicSpline { /// The method ensures that the transition between each polynomial is smooth and that the spline /// interpolation of the new nodes is calculated around `x = 0` in order to avoid that /// successive spline extensions with large x values become inaccurate. - pub fn extend_with_nodes(&mut self, node_x: Vec, node_y: Vec) -> Result<(), SplineError> { + pub fn extend_with_nodes(&mut self, node_x: Vec, node_y: Vec) -> Result<()> { let mut ext_node_x = Vec::with_capacity(node_x.len() + 1); let mut ext_node_y = Vec::with_capacity(node_x.len() + 1); @@ -631,16 +638,16 @@ impl CubicHermiteSpline { node_x: &[f64], node_y: &[f64], m: &[f64], - ) -> Result { + ) -> Result { let n = node_x.len(); if n < 3 { - return Err(NotEnoughNodes); + bail!(NotEnoughNodes); } if n != node_y.len() { - return Err(NotEqualNodes); + bail!(NotEqualNodes); } if n != m.len() { - return Err(NotEqualSlopes); + bail!(NotEqualSlopes); } let mut r = vec![Range::default(); node_x.len() - 1]; @@ -673,7 +680,7 @@ impl CubicHermiteSpline { node_x: &[f64], node_y: &[f64], slope_method: SlopeMethod, - ) -> Result { + ) -> Result { match slope_method { SlopeMethod::Akima => CubicHermiteSpline::from_nodes_with_slopes( node_x, @@ -772,9 +779,9 @@ pub enum SlopeMethod { Quadratic, } -fn akima_slopes(x: &[f64], y: &[f64]) -> Result, SplineError> { +fn akima_slopes(x: &[f64], y: &[f64]) -> Result> { if x.len() < 3 { - return Err(NotEnoughNodes); + bail!(NotEnoughNodes); } let mut m = vec![0f64; x.len()]; @@ -799,7 +806,7 @@ fn akima_slopes(x: &[f64], y: &[f64]) -> Result, SplineError> { for i in 0..new_x.len() - 1 { let dx = new_x[i + 1] - new_x[i]; if dx == 0f64 { - return Err(RedundantNodeX); + bail!(RedundantNodeX); } s[i] = (new_y[i + 1] - new_y[i]) / dx; } @@ -818,9 +825,9 @@ fn akima_slopes(x: &[f64], y: &[f64]) -> Result, SplineError> { Ok(m) } -fn quadratic_slopes(x: &[f64], y: &[f64]) -> Result, SplineError> { +fn quadratic_slopes(x: &[f64], y: &[f64]) -> Result> { if x.len() < 3 { - return Err(NotEnoughNodes); + bail!(NotEnoughNodes); } let mut m = vec![0f64; x.len()]; let q_i = lagrange_polynomial(x[0..3].to_vec(), y[0..3].to_vec()); From 140e2e8366ca54b6f6f2d92adb1ff6f30003fe90 Mon Sep 17 00:00:00 2001 From: axect Date: Sun, 14 Apr 2024 00:26:25 +0900 Subject: [PATCH 12/15] CHGE: Change thiserror -> anyhow in simpler & doc tests --- src/numerical/ode.rs | 9 ++++----- src/numerical/spline.rs | 8 ++++---- src/prelude/simpler.rs | 4 ++-- 3 files changed, 10 insertions(+), 11 deletions(-) diff --git a/src/numerical/ode.rs b/src/numerical/ode.rs index 933ae548..e550f858 100644 --- a/src/numerical/ode.rs +++ b/src/numerical/ode.rs @@ -73,7 +73,7 @@ //! vec![1f64] //! } //! -//! fn rhs(&self, t: f64, y: &[f64], dy: &mut [f64]) -> Result<(), ODEError> { +//! fn rhs(&self, t: f64, y: &[f64], dy: &mut [f64]) -> anyhow::Result<()> { //! Ok(dy[0] = (5f64 * t.powi(2) - y[0]) / (t + y[0]).exp()) //! } //! } @@ -97,7 +97,7 @@ use anyhow::{Result, bail}; /// vec![1.0, 2.0] /// } /// -/// fn rhs(&self, t: f64, y: &[f64], dy: &mut [f64]) -> Result<(), ODEError> { +/// fn rhs(&self, t: f64, y: &[f64], dy: &mut [f64]) -> anyhow::Result<()> { /// dy[0] = -0.5 * y[0]; /// dy[1] = y[0] - y[1]; /// Ok(()) @@ -131,7 +131,6 @@ pub trait ODEIntegrator { /// /// ```no_run /// use peroxide::fuga::*; -/// use anyhow::{Result, bail}; /// /// struct ConstrainedProblem { /// y_constraint: f64 @@ -139,9 +138,9 @@ pub trait ODEIntegrator { /// /// impl ODEProblem for ConstrainedProblem { /// fn initial_conditions(&self) -> Vec { vec![0.0] } // y_0 = 0 -/// fn rhs(&self, t: f64, y: &[f64], dy: &mut [f64]) -> Result<()> { +/// fn rhs(&self, t: f64, y: &[f64], dy: &mut [f64]) -> anyhow::Result<()> { /// if y[0] < self.y_constraint { -/// bail!(ODEError::ConstraintViolation(t, y.to_vec(), dy.to_vec())); +/// anyhow::bail!(ODEError::ConstraintViolation(t, y.to_vec(), dy.to_vec())); /// } else { /// // some function /// Ok(()) diff --git a/src/numerical/spline.rs b/src/numerical/spline.rs index 85fe4a4e..0a0be9b9 100644 --- a/src/numerical/spline.rs +++ b/src/numerical/spline.rs @@ -23,11 +23,11 @@ //! ## Members //! //! * `CubicSpline`: Structure for cubic spline -//! * `fn from_nodes(node_x: &[f64], node_y: &[f64]) -> Self` : Create a cubic spline from nodes -//! * `fn extend_with_nodes(&mut self, node_x: Vec, node_y: Vec)` : Extend the spline with nodes +//! * `fn from_nodes(node_x: &[f64], node_y: &[f64]) -> Result` : Create a cubic spline from nodes +//! * `fn extend_with_nodes(&mut self, node_x: Vec, node_y: Vec) -> Result<()>` : Extend the spline with nodes //! * `CubicHermiteSpline`: Structure for cubic Hermite spline -//! * `fn from_nodes_with_slopes(node_x: &[f64], node_y: &[f64], m: &[f64]) -> Self` : Create a Cubic Hermite spline from nodes with slopes -//! * `fn from_nodes(node_x: &[f64], node_y: &[f64], slope_method: SlopeMethod) -> Self` : Create a Cubic Hermite spline from nodes with slope estimation methods +//! * `fn from_nodes_with_slopes(node_x: &[f64], node_y: &[f64], m: &[f64]) -> Result` : Create a Cubic Hermite spline from nodes with slopes +//! * `fn from_nodes(node_x: &[f64], node_y: &[f64], slope_method: SlopeMethod) -> Result` : Create a Cubic Hermite spline from nodes with slope estimation methods //! * `SlopeMethod`: Enum for slope estimation methods //! * `Akima`: Akima's method to estimate slopes ([Akima (1970)](https://dl.acm.org/doi/abs/10.1145/321607.321609)) //! * `Quadratic`: Using quadratic interpolation to estimate slopes diff --git a/src/prelude/simpler.rs b/src/prelude/simpler.rs index 5436ea5c..2b1851f9 100644 --- a/src/prelude/simpler.rs +++ b/src/prelude/simpler.rs @@ -6,7 +6,7 @@ use crate::numerical::{ integral, integral::Integral::G7K15R, spline, - spline::{CubicHermiteSpline, SlopeMethod::Quadratic, SplineError}, + spline::{CubicHermiteSpline, SlopeMethod::Quadratic}, }; use crate::structure::matrix::{self, Matrix}; use crate::structure::polynomial; @@ -148,7 +148,7 @@ pub fn chebyshev_polynomial(n: usize) -> polynomial::Polynomial { polynomial::chebyshev_polynomial(n, polynomial::SpecialKind::First) } -pub fn cubic_hermite_spline(node_x: &[f64], node_y: &[f64]) -> Result { +pub fn cubic_hermite_spline(node_x: &[f64], node_y: &[f64]) -> anyhow::Result { spline::cubic_hermite_spline(node_x, node_y, Quadratic) } From 078476955c435581b720dfe868cc5f166817e4de Mon Sep 17 00:00:00 2001 From: axect Date: Sun, 14 Apr 2024 00:35:13 +0900 Subject: [PATCH 13/15] CHGE: Change thiserror -> anyhow in rand, dist, nonmacro, ode --- src/numerical/ode.rs | 2 +- src/statistics/dist.rs | 40 +++++++++++++++++++--------------------- src/statistics/rand.rs | 13 ++++++------- src/util/non_macro.rs | 29 ++++++++++++++++++----------- 4 files changed, 44 insertions(+), 40 deletions(-) diff --git a/src/numerical/ode.rs b/src/numerical/ode.rs index e550f858..09430feb 100644 --- a/src/numerical/ode.rs +++ b/src/numerical/ode.rs @@ -177,7 +177,7 @@ pub trait ODESolver { /// ``` /// use peroxide::fuga::*; /// -/// fn main() -> std::results::Result<(), Box> { +/// fn main() -> Result<(), Box> { /// let rkf = RKF45::new(1e-4, 0.9, 1e-6, 1e-1, 100); /// let basic_ode_solver = BasicODESolver::new(rkf); /// let (t_vec, y_vec) = basic_ode_solver.solve( diff --git a/src/statistics/dist.rs b/src/statistics/dist.rs index d0857493..b409c197 100644 --- a/src/statistics/dist.rs +++ b/src/statistics/dist.rs @@ -56,7 +56,6 @@ //! * Usage is very simple //! //! ```rust -//! extern crate peroxide; //! use peroxide::fuga::*; //! //! fn main() { @@ -82,7 +81,6 @@ //! * **Caution**: `Uniform(T, T)` generates `T` type samples (only for `Uniform`) //! //! ```rust -//! extern crate peroxide; //! use peroxide::fuga::*; //! //! fn main() { @@ -110,7 +108,6 @@ //! * In peroxide, main traits is Ziggurat - most efficient traits to generate random normal samples. //! * Code is based on a [C implementation](https://www.seehuhn.de/pages/ziggurat.html) by Jochen Voss. //! ```rust -//! extern crate peroxide; //! use peroxide::fuga::*; //! //! fn main() { @@ -135,7 +132,6 @@ //! * To generate beta random samples, Peroxide uses the `rand_distr::Beta` distribution from the `rand_distr` crate. //! //! ```rust -//! extern crate peroxide; //! use peroxide::fuga::*; //! //! fn main() { @@ -159,7 +155,6 @@ //! * To generate gamma random samples, Peroxide uses the `rand_distr::Gamma` distribution from the `rand_distr` crate. //! //! ```rust -//! extern crate peroxide; //! use peroxide::fuga::*; //! //! fn main() { @@ -183,7 +178,6 @@ //! * To generate binomial random samples, Peroxide uses the `rand_distr::Binomial` distribution from the `rand_distr` crate. //! //! ```rust -//! extern crate peroxide; //! use peroxide::fuga::*; //! //! fn main() { @@ -207,7 +201,6 @@ //! * To generate Student's t random samples, Peroxide uses the `rand_distr::StudentT` distribution from the `rand_distr` crate. //! //! ```rust -//! extern crate peroxide; //! use peroxide::fuga::*; //! //! fn main() { @@ -248,8 +241,8 @@ use crate::statistics::{ops::C, stat::Statistics}; use crate::util::non_macro::{linspace, seq}; use crate::util::useful::{auto_zip, find_interval}; use std::f64::consts::E; -use thiserror::Error; use self::WeightedUniformError::*; +use anyhow::{Result, bail}; /// One parameter distribution /// @@ -281,24 +274,30 @@ pub struct WeightedUniform> { intervals: Vec<(T, T)> } -#[derive(Debug, Clone, Copy, Error)] +#[derive(Debug, Clone, Copy)] pub enum WeightedUniformError { - #[error("all weights are zero")] AllZeroWeightError, - #[error("weights and intervals have different length")] LengthMismatchError, - #[error("no non-zero interval found")] NoNonZeroIntervalError, - #[error("weights are empty")] EmptyWeightError, } +impl std::fmt::Display for WeightedUniformError { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + WeightedUniformError::AllZeroWeightError => write!(f, "all weights are zero"), + WeightedUniformError::LengthMismatchError => write!(f, "weights and intervals have different length"), + WeightedUniformError::NoNonZeroIntervalError => write!(f, "no non-zero interval found"), + WeightedUniformError::EmptyWeightError => write!(f, "weights are empty"), + } + } +} + impl WeightedUniform { /// Create a new weighted uniform distribution /// /// # Examples /// ``` - /// extern crate peroxide; /// use peroxide::fuga::*; /// /// fn main() -> Result<(), Box> { @@ -311,17 +310,17 @@ impl WeightedUniform { /// Ok(()) /// } /// ``` - pub fn new(weights: Vec, intervals: Vec) -> Result { + pub fn new(weights: Vec, intervals: Vec) -> Result { let mut weights = weights; if weights.len() == 0 { - return Err(EmptyWeightError); + bail!(EmptyWeightError); } if weights.iter().all(|&x| x == 0f64) { - return Err(AllZeroWeightError); + bail!(AllZeroWeightError); } let mut intervals = auto_zip(&intervals); if weights.len() != intervals.len() { - return Err(LengthMismatchError); + bail!(LengthMismatchError); } // Remove zero weights & corresponding intervals @@ -354,7 +353,6 @@ impl WeightedUniform { /// /// # Examples /// ``` - /// extern crate peroxide; /// use peroxide::fuga::*; /// /// fn main() -> Result<(), Box> { @@ -372,7 +370,7 @@ impl WeightedUniform { /// } /// } /// ``` - pub fn from_max_pool_1d(f: F, (a, b): (f64, f64), n: usize, eps: f64) -> Result + pub fn from_max_pool_1d(f: F, (a, b): (f64, f64), n: usize, eps: f64) -> Result where F: Fn(f64) -> f64 + Copy { // Find non-zero intervals let mut a = a; @@ -393,7 +391,7 @@ impl WeightedUniform { } } if a >= b { - return Err(NoNonZeroIntervalError); + bail!(NoNonZeroIntervalError); } let domain = linspace(a, b, n+1); diff --git a/src/statistics/rand.rs b/src/statistics/rand.rs index b5a3ee0c..36590085 100644 --- a/src/statistics/rand.rs +++ b/src/statistics/rand.rs @@ -15,10 +15,10 @@ //! ``` //! //! * To want more detailed explanation, see [`rand` crate](https://crates.io/crates/rand) -//! +//! //! ## Piece-wise Rejection Sampling -//! -//! +//! +//! extern crate rand; use self::rand::distributions::uniform::SampleUniform; @@ -26,7 +26,7 @@ use self::rand::prelude::*; #[allow(unused_imports)] use crate::structure::matrix::*; -use crate::statistics::dist::{RNG, WeightedUniform, WeightedUniformError}; +use crate::statistics::dist::{RNG, WeightedUniform}; /// Small random number generator from seed /// @@ -64,7 +64,6 @@ pub fn stdrng_from_seed(seed: u64) -> StdRng { /// /// # Examples /// ``` -/// extern crate peroxide; /// use peroxide::fuga::*; /// /// let mut rng = thread_rng(); @@ -466,7 +465,7 @@ pub fn ziggurat(rng: &mut ThreadRng, sigma: f64) -> f64 { /// /// Ok(()) /// } -pub fn prs(f: F, n: usize, (a, b): (f64, f64), m: usize, eps: f64) -> Result, WeightedUniformError> +pub fn prs(f: F, n: usize, (a, b): (f64, f64), m: usize, eps: f64) -> anyhow::Result> where F: Fn(f64) -> f64 + Copy { let mut rng = thread_rng(); @@ -528,7 +527,7 @@ where F: Fn(f64) -> f64 + Copy { /// /// Ok(()) /// } -pub fn prs_with_rng(f: F, n: usize, (a, b): (f64, f64), m: usize, eps: f64, rng: &mut R) -> Result, WeightedUniformError> +pub fn prs_with_rng(f: F, n: usize, (a, b): (f64, f64), m: usize, eps: f64, rng: &mut R) -> anyhow::Result> where F: Fn(f64) -> f64 + Copy { let mut result = vec![0f64; n]; diff --git a/src/util/non_macro.rs b/src/util/non_macro.rs index 3263da53..1656d33f 100644 --- a/src/util/non_macro.rs +++ b/src/util/non_macro.rs @@ -35,15 +35,22 @@ use crate::structure::{ matrix::Shape::{Col, Row}, matrix::{matrix, Matrix, Shape}, }; -use thiserror::Error; use crate::traits::float::FloatWithPrecision; +use anyhow::{Result, bail}; -#[derive(Debug, Copy, Clone, Error)] +#[derive(Debug, Copy, Clone)] pub enum ConcatenateError { - #[error("To concatenate, vectors or matrices must have the same length")] DifferentLength, } +impl std::fmt::Display for ConcatenateError { + fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result { + match *self { + ConcatenateError::DifferentLength => write!(f, "To concatenate, vectors or matrices must have the same length"), + } + } +} + // ┌─────────────────────────────────────────────────────────┐ // R-like non-macro functions // └─────────────────────────────────────────────────────────┘ @@ -131,7 +138,7 @@ where /// Ok(()) /// } /// ``` -pub fn cbind(m1: Matrix, m2: Matrix) -> Result { +pub fn cbind(m1: Matrix, m2: Matrix) -> Result { let mut temp = m1; if temp.shape != Col { temp = temp.change_shape(); @@ -147,7 +154,7 @@ pub fn cbind(m1: Matrix, m2: Matrix) -> Result { let r = temp.row; if r != temp2.row { - return Err(ConcatenateError::DifferentLength); + bail!(ConcatenateError::DifferentLength); } v.extend_from_slice(&temp2.data[..]); c += temp2.col; @@ -170,7 +177,7 @@ pub fn cbind(m1: Matrix, m2: Matrix) -> Result { /// Ok(()) /// } /// ``` -pub fn rbind(m1: Matrix, m2: Matrix) -> Result { +pub fn rbind(m1: Matrix, m2: Matrix) -> Result { let mut temp = m1; if temp.shape != Row { temp = temp.change_shape(); @@ -186,7 +193,7 @@ pub fn rbind(m1: Matrix, m2: Matrix) -> Result { let mut r = temp.row; if c != temp2.col { - return Err(ConcatenateError::DifferentLength); + bail!(ConcatenateError::DifferentLength); } v.extend_from_slice(&temp2.data[..]); r += temp2.row; @@ -378,20 +385,20 @@ where } /// Numpy like column_stack -pub fn column_stack(v: &[Vec]) -> Result { +pub fn column_stack(v: &[Vec]) -> Result { let row = v[0].len(); if v.iter().any(|x| x.len() != row) { - return Err(ConcatenateError::DifferentLength); + bail!(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 { +pub fn row_stack(v: &[Vec]) -> Result { let col = v[0].len(); if v.iter().any(|x| x.len() != col) { - return Err(ConcatenateError::DifferentLength); + bail!(ConcatenateError::DifferentLength); } let data = v.iter().flatten().copied().collect(); Ok(matrix(data, v.len(), col, Row)) From ad2341d9c8dfd5a5014a37557d0be572a26106e2 Mon Sep 17 00:00:00 2001 From: axect Date: Sun, 14 Apr 2024 00:50:49 +0900 Subject: [PATCH 14/15] DOCS: Add docs for root --- src/numerical/ode.rs | 1 + src/numerical/root.rs | 172 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 173 insertions(+) diff --git a/src/numerical/ode.rs b/src/numerical/ode.rs index 09430feb..83be8d75 100644 --- a/src/numerical/ode.rs +++ b/src/numerical/ode.rs @@ -10,6 +10,7 @@ //! - `ODEError`: Enum for ODE errors. //! - `ReachedMaxStepIter`: Reached maximum number of steps per step. (internal error) //! - `ConstraintViolation(f64, Vec, Vec)`: Constraint violation. (user-defined error) +//! - ODE uses `anyhow` for error handling. So, you can customize your errors. //! //! ## Available integrators //! diff --git a/src/numerical/root.rs b/src/numerical/root.rs index 95fbf682..99cb2237 100644 --- a/src/numerical/root.rs +++ b/src/numerical/root.rs @@ -1,3 +1,175 @@ +//! # Root Finding Methods +//! +//! This module provides a collection of root finding algorithms for solving nonlinear equations. +//! It defines traits for representing root finding problems and root finding methods, and provides implementations of several common algorithms. +//! +//! ## Traits +//! +//! - `RootFindingProblem`: Defines the interface for a root finding problem. +//! It requires implementing the `function` method to evaluate the function at a given point, and the `initial_guess` method to provide an initial guess for the root. +//! Optionally, the `derivative` and `hessian` methods can be implemented to provide the derivative and Hessian of the function, respectively. +//! +//! - `I`: Input dimension +//! - `O`: Output dimension +//! - `T`: State type (e.g. `f64`, `(f64, f64)`, or etc.) +//! +//! - `RootFinder`: Defines the interface for a root finding method. +//! It requires implementing the `find` method, which takes a `RootFindingProblem` and returns the root of the function. +//! The `max_iter` and `tol` methods provide the maximum number of iterations and the tolerance for the root finding algorithm. +//! +//! ## Root Finding Methods +//! +//! - `BisectionMethod`: Implements the bisection method for finding roots of continuous functions. +//! It requires an initial interval that brackets the root. +//! +//! - Type Parameters: `I=1, O=1, T=(f64, f64)` +//! +//! - `NewtonMethod`: Implements Newton's method for finding roots of differentiable functions. +//! It requires an initial guess for the root and the derivative of the function. +//! +//! - Type Parameters: `I=1, O=1, T=f64` +//! +//! - `SecantMethod`: Implements the secant method for finding roots of differentiable functions. +//! It requires two initial guesses for the root. +//! +//! - Type Parameters: `I=1, O=1, T=f64` +//! +//! - `FalsePositionMethod`: Implements the false position method for finding roots of continuous functions. +//! It requires an initial interval that brackets the root. +//! +//! - Type Parameters: `I=1, O=1, T=(f64, f64)` +//! +//! - `BroydenMethod` (unimplemented): Implements Broyden's method for finding roots of systems of nonlinear equations. +//! +//! - Type Parameters: `I>=1, O>=1, T=([f64; I], [f64; I])` +//! +//! ## Convenient type aliases +//! +//! - `Pt`: Represents a point in N-dimensional space. (`[f64; N]`) +//! - `Intv`: Represents an interval in I-dimensional space. (`([f64; N], [f64; N])`) +//! - `Jaco`: Represents the Jacobian matrix of a function. (`[[f64; C]; R]`) +//! - `Hess`: Represents the Hessian matrix of a function. (`[[[f64; C]; C]; R]`) +//! +//! ## Examples +//! +//! ### Finding the root of a cubic function +//! +//! ```rust +//! use peroxide::fuga::*; +//! use anyhow::Result; +//! +//! fn main() -> Result<()> { +//! let problem = Cubic; +//! +//! let bisect = BisectionMethod { max_iter: 100, tol: 1e-6 }; +//! let newton = NewtonMethod { max_iter: 100, tol: 1e-6 }; +//! let false_pos = FalsePositionMethod { max_iter: 100, tol: 1e-6 }; +//! +//! let root_bisect = bisect.find(&problem)?; +//! let root_newton = newton.find(&problem)?; +//! let root_false_pos = false_pos.find(&problem)?; +//! +//! let result_bisect = problem.eval(root_bisect)?[0]; +//! let result_newton = problem.eval(root_newton)?[0]; +//! let result_false_pos = problem.eval(root_false_pos)?[0]; +//! +//! assert!(result_bisect.abs() < 1e-6); +//! assert!(result_newton.abs() < 1e-6); +//! assert!(result_false_pos.abs() < 1e-6); +//! +//! Ok(()) +//! } +//! +//! struct Cubic; +//! +//! impl Cubic { +//! fn eval(&self, x: [f64; 1]) -> Result<[f64; 1]> { +//! Ok([(x[0] - 1f64).powi(3)]) +//! } +//! } +//! +//! impl RootFindingProblem<1, 1, (f64, f64)> for Cubic { +//! fn function(&self, x: [f64; 1]) -> Result<[f64; 1]> { +//! self.eval(x) +//! } +//! +//! fn initial_guess(&self) -> (f64, f64) { +//! (0.0, 2.0) +//! } +//! } +//! +//! impl RootFindingProblem<1, 1, f64> for Cubic { +//! fn function(&self, x: [f64; 1]) -> Result<[f64; 1]> { +//! self.eval(x) +//! } +//! +//! fn initial_guess(&self) -> f64 { +//! 0.0 +//! } +//! +//! fn derivative(&self, x: [f64; 1]) -> Result> { +//! Ok([[3.0 * (x[0] - 1f64).powi(2)]]) +//! } +//! } +//! ``` +//! +//! This example demonstrates how to find the root of a cubic function `(x - 1)^3` using various root finding methods. +//! The `Cubic` struct implements the `RootFindingProblem` trait for both `(f64, f64)` and `f64` initial guess types, allowing the use of different root finding methods. +//! +//! ### Finding the root of the cosine function (error handling example) +//! +//! ```rust +//! use peroxide::fuga::*; +//! use anyhow::Result; +//! +//! fn main() -> Result<()> { +//! let problem = Cosine; +//! let newton = NewtonMethod { max_iter: 100, tol: 1e-6 }; +//! +//! let root_newton = match newton.find(&problem) { +//! Ok(x) => x, +//! Err(e) => { +//! println!("{:?}", e); +//! match e.downcast::>() { +//! Ok(RootError::ZeroDerivative(x)) => x, +//! Ok(e) => panic!("ok but {:?}", e), +//! Err(e) => panic!("err {:?}", e), +//! } +//! } +//! }; +//! +//! assert_eq!(root_newton[0], 0.0); +//! +//! Ok(()) +//! } +//! +//! struct Cosine; +//! +//! impl Cosine { +//! fn eval(&self, x: [f64; 1]) -> Result<[f64; 1]> { +//! Ok([x[0].cos()]) +//! } +//! } +//! +//! impl RootFindingProblem<1, 1, f64> for Cosine { +//! fn function(&self, x: [f64; 1]) -> Result<[f64; 1]> { +//! self.eval(x) +//! } +//! +//! fn initial_guess(&self) -> f64 { +//! 0.0 // should fail in newton (derivative is 0) +//! } +//! +//! fn derivative(&self, x: [f64; 1]) -> Result> { +//! Ok([[-x[0].sin()]]) +//! } +//! } +//! ``` +//! +//! This example shows how to find the root of the cosine function using Newton's method. +//! The `Cosine` struct implements the `RootFindingProblem` trait for the `f64` initial guess type. +//! The initial guess is set to `0.0`, which is a point where the derivative of the cosine function is 0. +//! This leads to the `NewtonMethod` returning a `RootError::ZeroDerivative` error, which is handled in the example. use anyhow::{Result, bail}; /// Point alias (`[f64; N]`) From e2e5852218ce4ba9120130d8c24badd8eaeea696 Mon Sep 17 00:00:00 2001 From: axect Date: Sun, 14 Apr 2024 01:08:40 +0900 Subject: [PATCH 15/15] RLSE: Ver 0.37.0 - Whole new root finding - from thiserror to anyhow - Add zenodo badge --- Cargo.toml | 2 +- README.md | 21 +++++++++++++++++++-- TODO.md | 45 +++++++++++++++++++++++---------------------- 3 files changed, 43 insertions(+), 25 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 772dca9e..eb0cc3e3 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "peroxide" -version = "0.36.4" +version = "0.37.0" 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/README.md b/README.md index a9c01ce1..020e90a1 100644 --- a/README.md +++ b/README.md @@ -309,7 +309,7 @@ How's that? Let me know if there's anything else you'd like me to improve! ## Latest README version -Corresponding to `0.36.0` +Corresponding to `0.37.0` ## Pre-requisite @@ -456,11 +456,13 @@ Corresponding to `0.36.0` - In [examples](./examples) directory, there are some examples. +- In [tests](./tests) directory, there are some useful tests. + - More examples are in [Peroxide Gallery](https://github.com/Axect/Peroxide_Gallery). ## Release Info -To see [RELEASES.md](RELEASES.md) +To see [RELEASES.md](./RELEASES.md) ## Contributes Guide @@ -473,3 +475,18 @@ Peroxide is licensed under dual licenses - Apache License 2.0 and MIT License. ## TODO To see [TODO.md](./TODO.md) + +## Cite Peroxide + +Hey there! +If you're using Peroxide in your research or project, you're not required to cite us. +But if you do, we'd be really grateful! 😊 + +To make citing Peroxide easy, we've created a DOI through Zenodo. Just click on this badge: + +[![DOI](https://zenodo.org/badge/130400565.svg)](https://zenodo.org/doi/10.5281/zenodo.10815823) + +This will take you to the Zenodo page for Peroxide. +At the bottom, you'll find the citation information in various formats like BibTeX, RIS, and APA. + +So, if you want to acknowledge the work we've put into Peroxide, citing us would be a great way to do it! Thanks for considering it, we appreciate your support! 👍 diff --git a/TODO.md b/TODO.md index 60a2de21..4842dcf4 100644 --- a/TODO.md +++ b/TODO.md @@ -1,20 +1,16 @@ # TODO -## 2021.01.14 +## 2024.04.14 ### Primary +- [ ] Whole new Optimize (Trait based) + - [ ] Pure Rust implementation of Linear Algebra - [x] LU (Completely Pivoting) - [x] LU (Partial Pivoting) - [ ] QR - [ ] SVD -- [ ] Add more IO options for DataFrame - - [x] CSV (`csv` feature) - - [x] NetCDF (`nc` feature) - - [ ] JSON - - [ ] Arrow IPC - - [ ] Parquet ### Subs @@ -30,30 +26,16 @@ - [x] Student's t - [x] Uniform - [ ] Wishart + - [x] Weighted Uniform - [ ] Implement special polynomial - [x] Legendre - [ ] Bessel - [ ] Hermite - [ ] Implement convenient structure of Neural Network -- [ ] Documentized - - [x] Vector - - [x] Matrix - - [x] Linear Algebra - - [x] Functional Programming - - [x] Statistics - - [ ] Interpolation & Spline - - [x] ODE - - [ ] Macros - - [ ] Machine Learning - - [x] Optimize - - [x] Automatic Differentiation - - [x] DataFrame - [ ] Add Statistical regression - [ ] Gaussian Kernel - [ ] Logistic Kernel -- [ ] Make or Use pure Rust plot library - [ ] Implement more Eigenvalue algorithms -- [ ] Implement more spline algorithms - [ ] Complex matrix ## Complete @@ -75,3 +57,22 @@ - [x] Replace `proc_macro` for `AD` with ordinary macro or Enum - [x] Make `csv` optional - [x] Remove `dual`, `hyperdual` and modify `Real`, `Number` (How to bind `f64` & `AD` effectively?) +- [x] Add more IO options for DataFrame + - [x] CSV (`csv` feature) + - [x] NetCDF (`nc` feature) + - [x] Parquet +- [x] Documentized + - [x] Vector + - [x] Matrix + - [x] Linear Algebra + - [x] Functional Programming + - [x] Statistics + - [x] Interpolation & Spline + - [x] ODE + - [x] Macros + - [x] Optimize + - [x] Automatic Differentiation + - [x] DataFrame +- [x] Implement more spline algorithms +- [x] Whole new ODE (trait based) +- [x] Whole new root finding (trait based)