From 2695b0e7580d8b522983528ffc4010f319d78fb1 Mon Sep 17 00:00:00 2001 From: Giorgio Comitini Date: Mon, 15 Apr 2024 16:34:03 +0200 Subject: [PATCH 1/4] Do not include legend box if there is no legend --- src/util/plot.rs | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/util/plot.rs b/src/util/plot.rs index de32b46a..90499195 100644 --- a/src/util/plot.rs +++ b/src/util/plot.rs @@ -651,10 +651,14 @@ impl Plot for Plot2D { } } + if !legends.is_empty() { + plot_string.push_str("plt.legend()\n"); + } + if self.tight { - plot_string.push_str(&format!("plt.legend()\nplt.savefig(pa, dpi={}, bbox_inches='tight')", dpi)[..]); + plot_string.push_str(&format!("plt.savefig(pa, dpi={}, bbox_inches='tight')", dpi)[..]); } else { - plot_string.push_str(&format!("plt.legend()\nplt.savefig(pa, dpi={})", dpi)[..]); + plot_string.push_str(&format!("plt.savefig(pa, dpi={})", dpi)[..]); } py.run(&plot_string[..], Some(&globals), None)?; From e0361301bf241e9bd280aaff6eea8299a9ba93fd Mon Sep 17 00:00:00 2001 From: axect Date: Tue, 16 Apr 2024 11:46:14 +0900 Subject: [PATCH 2/4] ADD: Add rtol field to Broyden method --- src/numerical/root.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/numerical/root.rs b/src/numerical/root.rs index 9594f2f2..8df0724e 100644 --- a/src/numerical/root.rs +++ b/src/numerical/root.rs @@ -618,9 +618,9 @@ impl RootFinder<1, 1, (f64, f64)> for FalsePositionMethod { pub struct BroydenMethod { pub max_iter: usize, pub tol: f64, + pub rtol: f64, } - #[allow(unused_variables, non_snake_case)] impl RootFinder> for BroydenMethod { fn max_iter(&self) -> usize { @@ -651,7 +651,7 @@ impl RootFinder> for BroydenMethod return Ok(x1); } let dx = x1.iter().zip(x0.iter()).map(|(x1, x0)| x1 - x0).collect::>(); - if dx.norm(Norm::L2) < self.tol { + if dx.norm(Norm::L2) < self.rtol { return Ok(x1); } let df = fx1.iter().zip(fx0.iter()).map(|(fx1, fx0)| fx1 - fx0).collect::>(); From 963ca3e988264e3e9b4d12cf73f86f60b9b92356 Mon Sep 17 00:00:00 2001 From: axect Date: Tue, 16 Apr 2024 14:22:02 +0900 Subject: [PATCH 3/4] IMPL: Implement high level macros for root finding --- Cargo.toml | 1 + examples/root_macro_test.rs | 28 ++++++ examples/root_test.rs | 3 + src/fuga/mod.rs | 1 + src/numerical/root.rs | 195 +++++++++++++++++++++++++++++++++++- src/prelude/mod.rs | 1 + 6 files changed, 228 insertions(+), 1 deletion(-) create mode 100644 examples/root_macro_test.rs diff --git a/Cargo.toml b/Cargo.toml index 8f8d4f49..4b38f596 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -29,6 +29,7 @@ matrixmultiply = { version = "0.3", features = ["threading"] } peroxide-ad = "0.3" peroxide-num = "0.1" anyhow = "1.0" +paste = "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"] } diff --git a/examples/root_macro_test.rs b/examples/root_macro_test.rs new file mode 100644 index 00000000..f2012939 --- /dev/null +++ b/examples/root_macro_test.rs @@ -0,0 +1,28 @@ +#[macro_use] +extern crate peroxide; +use peroxide::fuga::*; +use anyhow::Result; + +fn main() -> Result<()> { + let root_bisect = bisection!(f, (0.0, 2.0), 100, 1e-6); + let root_newton = newton!(f, 0.0, 100, 1e-6); + let root_false_pos = false_position!(f, (0.0, 2.0), 100, 1e-6); + let root_secant = secant!(f, (0.0, 2.0), 100, 1e-6); + + println!("root_bisect: {}", root_bisect); + println!("root_newton: {}", root_newton); + println!("root_false_pos: {}", root_false_pos); + println!("root_secant: {}", root_secant); + + assert!(f(root_bisect).abs() < 1e-6); + assert!(f(root_newton).abs() < 1e-6); + assert!(f(root_false_pos).abs() < 1e-6); + assert!(f(root_secant).abs() < 1e-6); + + Ok(()) +} + +#[ad_function] +fn f(x: f64) -> f64 { + (x - 1f64).powi(3) +} diff --git a/examples/root_test.rs b/examples/root_test.rs index b797a39d..06bca33e 100644 --- a/examples/root_test.rs +++ b/examples/root_test.rs @@ -6,12 +6,15 @@ fn main() -> Result<()> { 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 secant = SecantMethod { 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)?; + let result_secant = secant.find(&problem)?; println!("{:?}", result_bisect); println!("{:?}", result_newton); println!("{:?}", result_false_pos); + println!("{:?}", result_secant); Ok(()) } diff --git a/src/fuga/mod.rs b/src/fuga/mod.rs index 9fb97357..216dcbd8 100644 --- a/src/fuga/mod.rs +++ b/src/fuga/mod.rs @@ -195,6 +195,7 @@ pub use crate::ml::reg::*; pub use crate::util::plot::*; pub use anyhow; +pub use paste; // ============================================================================= // Enums diff --git a/src/numerical/root.rs b/src/numerical/root.rs index 8df0724e..5fabfa7d 100644 --- a/src/numerical/root.rs +++ b/src/numerical/root.rs @@ -51,6 +51,47 @@ //! - `Jaco`: Represents the Jacobian matrix of a function. (`[[f64; C]; R]`) //! - `Hess`: Represents the Hessian matrix of a function. (`[[[f64; C]; C]; R]`) //! +//! ## High-level macros +//! +//! Peroxide also provides high-level macros for root finding. +//! Assume `f: fn(f64) -> f64`. +//! +//! - `bisection!(f, (a,b), max_iter, tol)` +//! - `newton!(f, x0, max_iter, tol)`: (**Caution**: newton macro requires `#[ad_function]` attribute) +//! - `secant!(f, (x0, x1), max_iter, tol)` +//! - `false_position!(f, (a,b), max_iter, tol)` +//! +//! ```rust +//! #[macro_use] +//! extern crate peroxide; +//! use peroxide::fuga::*; +//! use anyhow::Result; +//! +//! fn main() -> Result<()> { +//! let root_bisect = bisection!(f, (0.0, 2.0), 100, 1e-6); +//! let root_newton = newton!(f, 0.0, 100, 1e-6); +//! let root_false_pos = false_position!(f, (0.0, 2.0), 100, 1e-6); +//! let root_secant = secant!(f, (0.0, 2.0), 100, 1e-6); +//! +//! println!("root_bisect: {}", root_bisect); +//! println!("root_newton: {}", root_newton); +//! println!("root_false_pos: {}", root_false_pos); +//! println!("root_secant: {}", root_secant); +//! +//! assert!(f(root_bisect).abs() < 1e-6); +//! assert!(f(root_newton).abs() < 1e-6); +//! assert!(f(root_false_pos).abs() < 1e-6); +//! assert!(f(root_secant).abs() < 1e-6); +//! +//! Ok(()) +//! } +//! +//! #[ad_function] +//! fn f(x: f64) -> f64 { +//! (x - 1f64).powi(3) +//! } +//! ``` +//! //! ## Examples //! //! ### Finding the root of a cubic function @@ -177,6 +218,154 @@ use crate::traits::math::{Normed, Norm, LinearOp}; use crate::traits::sugar::{ConvToMat, VecOps}; use crate::util::non_macro::zeros; +// ┌─────────────────────────────────────────────────────────┐ +// High level macro +// └─────────────────────────────────────────────────────────┘ +/// High level macro for bisection +/// +/// # Arguments +/// +/// - `f`: `fn(f64) -> f64` +/// - `(a, b)`: `(f64, f64)` +/// - `max_iter`: `usize` +/// - `tol`: `f64` +#[macro_export] +macro_rules! bisection { + ($f:ident, ($a:expr, $b:expr), $max_iter:expr, $tol:expr) => {{ + struct BisectionProblem; + + impl RootFindingProblem<1, 1, (f64, f64)> for BisectionProblem { + fn initial_guess(&self) -> (f64, f64) { + ($a, $b) + } + + fn function(&self, x: [f64; 1]) -> Result<[f64; 1]> { + Ok([$f(x[0])]) + } + } + + let problem = BisectionProblem; + let bisection = BisectionMethod { max_iter: $max_iter, tol: $tol }; + let root = bisection.find(&problem)?; + root[0] + }} +} + +/// High level macro for newton (using Automatic differentiation) +/// +/// # Requirements +/// +/// - This macro requires the function with `ad_function` +/// +/// ```rust +/// use peroxide::fuga::*; +/// +/// #[ad_function] +/// fn f(x: f64) -> f64 { +/// (x - 1f64).powi(3) +/// } +/// ``` +/// +/// # Arguments +/// +/// - `f`: `fn(f64) -> f64` +/// - `x`: `f64` +/// - `max_iter`: `usize` +/// - `tol`: `f64` +#[macro_export] +macro_rules! newton { + ($f:ident, $x:expr, $max_iter:expr, $tol:expr) => {{ + use paste::paste; + struct NewtonProblem; + + impl RootFindingProblem<1, 1, f64> for NewtonProblem { + fn initial_guess(&self) -> f64 { + $x + } + + fn function(&self, x: [f64; 1]) -> Result<[f64; 1]> { + Ok([$f(x[0])]) + } + + fn derivative(&self, x: [f64; 1]) -> Result> { + paste! { + let x_ad = AD1(x[0], 1f64); + Ok([[[<$f _ad>](x_ad).dx()]]) + } + } + } + + let problem = NewtonProblem; + let newton = NewtonMethod { max_iter: $max_iter, tol: $tol }; + let root = newton.find(&problem)?; + root[0] + }} +} + +/// High level macro for false position +/// +/// # Arguments +/// +/// - `f`: `fn(f64) -> f64` +/// - `(a, b)`: `(f64, f64)` +/// - `max_iter`: `usize` +/// - `tol`: `f64` +#[macro_export] +macro_rules! false_position { + ($f:ident, ($a:expr, $b:expr), $max_iter:expr, $tol:expr) => {{ + struct FalsePositionProblem; + + impl RootFindingProblem<1, 1, (f64, f64)> for FalsePositionProblem { + fn initial_guess(&self) -> (f64, f64) { + ($a, $b) + } + + fn function(&self, x: [f64; 1]) -> Result<[f64; 1]> { + Ok([$f(x[0])]) + } + } + + let problem = FalsePositionProblem; + let false_position = FalsePositionMethod { max_iter: $max_iter, tol: $tol }; + let root = false_position.find(&problem)?; + root[0] + }} +} + +/// High level macro for secant +/// +/// # Arguments +/// +/// - `f`: `fn(f64) -> f64` +/// - `(a, b)`: `(f64, f64)` +/// - `max_iter`: `usize` +/// - `tol`: `f64` +#[macro_export] +macro_rules! secant { + ($f:ident, ($a:expr, $b:expr), $max_iter:expr, $tol:expr) => {{ + struct SecantProblem; + + impl RootFindingProblem<1, 1, (f64, f64)> for SecantProblem { + fn initial_guess(&self) -> (f64, f64) { + ($a, $b) + } + + fn function(&self, x: [f64; 1]) -> Result<[f64; 1]> { + Ok([$f(x[0])]) + } + } + + let problem = SecantProblem; + let secant = SecantMethod { max_iter: $max_iter, tol: $tol }; + let root = secant.find(&problem)?; + root[0] + }} +} + + +// ┌─────────────────────────────────────────────────────────┐ +// Type aliases +// └─────────────────────────────────────────────────────────┘ /// Point alias (`[f64; N]`) pub type Pt = [f64; N]; /// Interval alias (`([f64; N], [f64; N])`) @@ -186,6 +375,9 @@ pub type Jaco = [[f64; C]; R]; /// Hessian alias (`[[[f64; C]; C]; R]`) pub type Hess = [[[f64; C]; C]; R]; +// ┌─────────────────────────────────────────────────────────┐ +// Traits +// └─────────────────────────────────────────────────────────┘ /// Trait to define a root finding problem /// /// # Type Parameters @@ -484,8 +676,9 @@ impl RootFinder<1, 1, (f64, f64)> for SecantMethod { bail!(RootError::ZeroSecant([x0], [x1])); } + let f0_old = f0; f0 = f1; - (x0, x1) = (x1, x1 - f1 * (x1 - x0) / (f1 - f0)) + (x0, x1) = (x1, x1 - f1 * (x1 - x0) / (f1 - f0_old)) } bail!(RootError::NotConverge([x1])); } diff --git a/src/prelude/mod.rs b/src/prelude/mod.rs index 2a5b9f22..7836b46d 100644 --- a/src/prelude/mod.rs +++ b/src/prelude/mod.rs @@ -218,3 +218,4 @@ pub use simpler::SimpleParquet; pub use crate::util::plot::*; pub use anyhow; +pub use paste; From 2bf4c6387046cd5dc2b0553d85e66279c10e3a2d Mon Sep 17 00:00:00 2001 From: axect Date: Tue, 16 Apr 2024 14:29:09 +0900 Subject: [PATCH 4/4] RLSE: Ver 0.37.2 - Do not include legend box if there is no legend (#58) - Add rtol for Broyden method - High-level macros for root finding --- Cargo.toml | 2 +- RELEASES.md | 10 ++++++++++ examples/broyden_test.rs | 2 +- src/numerical/root.rs | 2 +- 4 files changed, 13 insertions(+), 3 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 4b38f596..a1e0ef1e 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "peroxide" -version = "0.37.1" +version = "0.37.2" authors = ["axect "] edition = "2018" description = "Rust comprehensive scientific computation library contains linear algebra, numerical analysis, statistics and machine learning tools with farmiliar syntax" diff --git a/RELEASES.md b/RELEASES.md index 2bb06fd9..f7ee308d 100644 --- a/RELEASES.md +++ b/RELEASES.md @@ -1,3 +1,13 @@ +# Release 0.37.2 (2024-04-16) + +- Do not include legend box if there is no legend ([#58](https://github.com/Axect/Peroxide/pull/58)) (Thanks to [@GComitini](https://github.com/GComitini)) +- Add `rtol` field to `BroydenMethod` +- Implement high-level macros for root finding + - `bisection!(f, (a,b), max_iter, tol)` + - `newton!(f, x0, max_iter, tol)` (require `#[ad_function]` attribute) + - `secant!(f, (a,b), max_iter, tol)` + - `false_position!(f, (a,b), max_iter, tol)` + # Release 0.37.1 (2024-04-15) - Implement `BrodenMethod`: Broyden's method (`I>=1, O>=1, T=([f64; I], [f64; I])`) diff --git a/examples/broyden_test.rs b/examples/broyden_test.rs index 164f4c4b..6ddf5653 100644 --- a/examples/broyden_test.rs +++ b/examples/broyden_test.rs @@ -3,7 +3,7 @@ use peroxide::numerical::root::{Pt, Intv}; fn main() -> Result<(), Box> { let problem = Quadratic; - let broyden = BroydenMethod { max_iter: 100, tol: 1e-6 }; + let broyden = BroydenMethod { max_iter: 100, tol: 1e-6, rtol: 1e-6 }; let root = broyden.find(&problem)?; let result = problem.function(root)?; diff --git a/src/numerical/root.rs b/src/numerical/root.rs index 5fabfa7d..33a04cb1 100644 --- a/src/numerical/root.rs +++ b/src/numerical/root.rs @@ -781,7 +781,7 @@ impl RootFinder<1, 1, (f64, f64)> for FalsePositionMethod { /// /// fn main() -> Result<(), Box> { /// let problem = CircleTangentLine; -/// let broyden = BroydenMethod { max_iter: 100, tol: 1e-6 }; +/// let broyden = BroydenMethod { max_iter: 100, tol: 1e-6, rtol: 1e-6 }; /// /// let root = broyden.find(&problem)?; /// let result = problem.function(root)?;