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

@@ -1,6 +1,6 @@
 [package]
 name = "peroxide"
-version = "0.37.1"
+version = "0.37.2"
 authors = ["axect <[email protected]>"]
 edition = "2018"
 description = "Rust comprehensive scientific computation library contains linear algebra, numerical analysis, statistics and machine learning tools with farmiliar syntax"
@@ -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"] }

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])`)

examples/broyden_test.rs

@@ -3,7 +3,7 @@ use peroxide::numerical::root::{Pt, Intv};
 
 fn main() -> Result<(), Box<dyn std::error::Error>> {
     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)?;
 

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)
+}

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(())
 }

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

src/numerical/root.rs

@@ -51,6 +51,47 @@
 //! - `Jaco<const R: usize, const C: usize>`: Represents the Jacobian matrix of a function. (`[[f64; C]; R]`)
 //! - `Hess<const R: usize, const C: usize>`: 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<Jaco<1, 1>> {
+                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<const N: usize> = [f64; N];
 /// Interval alias (`([f64; N], [f64; N])`)
@@ -186,6 +375,9 @@ pub type Jaco<const R: usize, const C: usize> = [[f64; C]; R];
 /// Hessian alias (`[[[f64; C]; C]; R]`)
 pub type Hess<const R: usize, const C: usize> = [[[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]));
     }
@@ -588,7 +781,7 @@ impl RootFinder<1, 1, (f64, f64)> for FalsePositionMethod {
 /// 
 /// fn main() -> Result<(), Box<dyn std::error::Error>> {
 ///     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)?;
@@ -618,9 +811,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<const I: usize, const O: usize> RootFinder<I, O, Intv<I>> for BroydenMethod {
     fn max_iter(&self) -> usize {
@@ -651,7 +844,7 @@ impl<const I: usize, const O: usize> RootFinder<I, O, Intv<I>> for BroydenMethod
                 return Ok(x1);
             }
             let dx = x1.iter().zip(x0.iter()).map(|(x1, x0)| x1 - x0).collect::<Vec<_>>();
-            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::<Vec<_>>();

src/prelude/mod.rs

@@ -218,3 +218,4 @@ pub use simpler::SimpleParquet;
 pub use crate::util::plot::*;
 
 pub use anyhow;
+pub use paste;

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)?;