From 613455947112d6ebf1dd2fdbaf376a5ec03013c7 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=A9r=C3=B4me?= <furious.jet@gmail.com>
Date: Sun, 26 Mar 2023 22:21:30 +0200
Subject: [PATCH] Added golden section search

---
 src/main.rs                  | 36 ++++++++++++++++++++------
 src/univariate_minimizers.rs | 49 ++++++++++++++++++++++++++++++++++++
 2 files changed, 78 insertions(+), 7 deletions(-)
 create mode 100644 src/univariate_minimizers.rs

diff --git a/src/main.rs b/src/main.rs
index 4a854d5..971e0ed 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -1,6 +1,7 @@
 extern crate colored;
 use colored::Colorize;
 mod univariate_solvers;
+mod univariate_minimizers;
 
 ///  x   sin(x)
 /// e  + ──────
@@ -43,7 +44,7 @@ fn check_result(x : f64, x_true : f64, tol : f64, test_name: &str, verbose: bool
         return 1;
     } else {
         if verbose {
-            println!("{}\t: x = {}\tf(x) = {}\t{} (expected {})", test_name_padded, format!("{:<20}", x), format!("{:<40}", fct(x)), "failed".red(), x_true);
+            println!("{}\t: x = {}\tf(x) = {}\t{} (expected {}, delta = {})", test_name_padded, format!("{:<20}", x), format!("{:<40}", fct(x)), "failed".red(), x_true, (x - x_true));
         } else {
             println!("{} {} : expected {}, got {}", test_name_padded, "failed".red(), x_true, x);
         }
@@ -60,18 +61,18 @@ fn print_test_results(num_tests_passed: u32, num_tests_total: u32) {
     }
 }
 
-fn main() {
-    println!("Testing Rust numerical solvers.");
-    let mut num_tests_passed : u32 = 0;
-    let num_tests_total : u32 = 7;
-    let verbose : bool = true;
+fn test_univariate_solvers(verbose: bool) {
+    println!("Testing univariate numerical solvers.");
     let x0 : f64 = 1.0;
     let tol : f64 = 1e-10;
     let max_iter : u32 = 100;
     let dx_num : f64 = 1e-6;
+    let mut num_tests_passed : u32 = 0;
+    let num_tests_total : u32 = 7;
+
     let x_mathematica: f64 = -3.26650043678562449167148755288;// 30 digits of precision
     let x_mathematica_2: f64 = -6.27133405258685307845641527902;// 30 digits of precision
-    let x_mathematica_3: f64 = -9.42553801930504142668603949182;// 30 digits of precision
+    // let x_mathematica_3: f64 = -9.42553801930504142668603949182;// 30 digits of precision
     let x_newton: f64 = univariate_solvers::newton_solve(&(fct as fn(f64) -> f64), &(dfct as fn(f64) -> f64), x0, tol, max_iter);
     let x_newton_num: f64 = univariate_solvers::newton_solve_num(&(fct as fn(f64) -> f64), x0, tol, dx_num, max_iter);
     let x_halley: f64 = univariate_solvers::halley_solve(&(fct as fn(f64) -> f64), &(dfct as fn(f64) -> f64), &(ddfct as fn(f64) -> f64), x0, tol, max_iter, false).unwrap();
@@ -89,3 +90,24 @@ fn main() {
 
     print_test_results(num_tests_passed, num_tests_total);
 }
+
+fn test_univariate_optimizers(verbose: bool) {
+    let tol : f64 = 1e-10;
+    let max_iter : u32 = 100;
+    let dx_num : f64 = 1e-6;
+    let mut num_tests_passed : u32 = 0;
+    let num_tests_total : u32 = 1;
+
+    let x_mathematica: f64 = -4.54295618675514754103476876324;// 30 digits of precision
+    let y_mathematica: f64 = -0.206327079359226884630654987440;// 30 digits of precision
+    let x_golden_section: f64 = univariate_minimizers::golden_section_minimize(&(fct as fn(f64) -> f64), -7.0, -1.0, tol);
+    num_tests_passed += check_result(x_golden_section, x_mathematica, tol*1e2, "Golden section search", verbose);
+    print_test_results(num_tests_passed, num_tests_total);
+}
+
+fn main() {
+    println!("Testing Rust numerical solvers.");
+    let verbose : bool = true;
+    test_univariate_solvers(verbose);
+    test_univariate_optimizers(verbose);
+}
diff --git a/src/univariate_minimizers.rs b/src/univariate_minimizers.rs
new file mode 100644
index 0000000..c4227cf
--- /dev/null
+++ b/src/univariate_minimizers.rs
@@ -0,0 +1,49 @@
+/// Golden section search for minimizing a function f(x)
+/// @param f function to minimize
+/// @param a left bracket
+/// @param b right bracket
+/// @param tol tolerance
+/// @return solution
+/// @note The interval [a, b] must bracket the minimum, and the function must have f''(x) > 0 over the interval [a, b] to garantee convergence.
+pub fn golden_section_minimize<F>(f : F, mut a: f64, mut b: f64, tol: f64) -> f64
+where F : Fn(f64) -> f64
+{
+    let invphi: f64 = (f64::sqrt(5.0) - 1.0) / 2.0;  // 1 / phi
+    let invphi2: f64 = (3.0 - f64::sqrt(5.0)) / 2.0;  // 1 / phi^2
+
+    a = f64::min(a, b);
+    b = f64::max(a, b);
+    let mut h: f64 = b - a;
+    if h <= tol {
+        return (a + b)/2.0;
+    }
+
+    // Required steps to achieve tolerance
+    let n: u32 = (f64::ceil(f64::ln(tol / h) / f64::ln(invphi))) as u32;
+
+    let mut c: f64 = a + invphi2 * h;
+    let mut d: f64 = a + invphi * h;
+    let mut yc: f64 = f(c);
+    let mut yd: f64 = f(d);
+
+    for _ in 0..n {
+        if yc < yd {  // yc > yd to find the maximum
+            b = d;
+            d = c;
+            yd = yc;
+            h = invphi * h;
+            c = a + invphi2 * h;
+            yc = f(c);
+        } else {
+            a = c;
+            c = d;
+            yc = yd;
+            h = invphi * h;
+            d = a + invphi * h;
+            yd = f(d);
+        }
+    }
+
+    if yc < yd { return (a + d)/2.0; }
+    else { return (c + b)/2.0; }
+}