From 35a3705763a3c231daf79a9b1d1ba7c2fbecf1c8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=A9r=C3=B4me?= Date: Sat, 16 Mar 2024 10:21:26 +0100 Subject: [PATCH] Add ITP solve method. --- src/main.rs | 38 ++++++++++--------- src/univariate_solvers.rs | 79 ++++++++++++++++++++++++++++++++++++++- 2 files changed, 98 insertions(+), 19 deletions(-) diff --git a/src/main.rs b/src/main.rs index 64d7af4..191fac4 100644 --- a/src/main.rs +++ b/src/main.rs @@ -119,7 +119,7 @@ fn test_univariate_solvers(verbose: bool) { 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 mut num_tests_total : u32 = 0; let x_mathematica: f64 = -3.26650043678562449167148755288;// 30 digits of precision let x_mathematica_2: f64 = -6.27133405258685307845641527902;// 30 digits of precision @@ -131,37 +131,39 @@ fn test_univariate_solvers(verbose: bool) { let x_bisection : f64 = univariate_solvers::bisection_solve(&(fct as fn(f64) -> f64), -5.0, 1.0, tol).unwrap(); let x_secant : f64 = univariate_solvers::secant_solve(&(fct as fn(f64) -> f64), -1.0, 1.0, tol, max_iter); let x_ridder : f64 = univariate_solvers::ridder_solve(&(fct as fn(f64) -> f64), -5.0, 1.0, tol, max_iter).unwrap(); - num_tests_passed += check_result(x_newton, x_mathematica, tol, "Newton's method", verbose); - num_tests_passed += check_result(x_newton_num, x_mathematica, tol, "Newton's method (num)", verbose); - num_tests_passed += check_result(x_halley, x_mathematica_2, tol, "Halley's method", verbose); - num_tests_passed += check_result(x_halley_num, x_mathematica_2, tol, "Halley's method (num)", verbose); - num_tests_passed += check_result(x_bisection, x_mathematica, tol, "Bisection method", verbose); - num_tests_passed += check_result(x_secant, x_mathematica, tol, "Secant method", verbose); - num_tests_passed += check_result(x_ridder, x_mathematica, tol, "Ridder's method", verbose); + let x_itp : f64 = univariate_solvers::itp_solve(&(fct as fn(f64) -> f64), -5.0, 1.0, 0.1, 2.0, 1.0, tol); + num_tests_passed += check_result(x_newton, x_mathematica, tol, "Newton's method", verbose); num_tests_total += 1; + num_tests_passed += check_result(x_newton_num, x_mathematica, tol, "Newton's method (num)", verbose); num_tests_total += 1; + num_tests_passed += check_result(x_halley, x_mathematica_2, tol, "Halley's method", verbose); num_tests_total += 1; + num_tests_passed += check_result(x_halley_num, x_mathematica_2, tol, "Halley's method (num)", verbose); num_tests_total += 1; + num_tests_passed += check_result(x_bisection, x_mathematica, tol, "Bisection method", verbose); num_tests_total += 1; + num_tests_passed += check_result(x_secant, x_mathematica, tol, "Secant method", verbose); num_tests_total += 1; + num_tests_passed += check_result(x_ridder, x_mathematica, tol, "Ridder's method", verbose); num_tests_total += 1; + num_tests_passed += check_result(x_itp, x_mathematica, tol, "ITP method", verbose); num_tests_total += 1; 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 max_iter : u32 = 100; + // let dx_num : f64 = 1e-6; let mut num_tests_passed : u32 = 0; - let num_tests_total : u32 = 1; + let mut num_tests_total : u32 = 0; let x_mathematica: f64 = -4.54295618675514754103476876324;// 30 digits of precision - let y_mathematica: f64 = -0.206327079359226884630654987440;// 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); + num_tests_passed += check_result(x_golden_section, x_mathematica, tol*1e2, "Golden section search", verbose); num_tests_total += 1; print_test_results(num_tests_passed, num_tests_total); } fn test_multivariate_optimizers(verbose: bool) { let tol : f64 = 1e-10; let max_iter : u32 = 1000; - let dx_num : f64 = 1e-6; + // let dx_num : f64 = 1e-6; let mut num_tests_passed : u32 = 0; - let num_tests_total : u32 = 1; + let mut num_tests_total : u32 = 0; let tol_x: f64 = 1e-4; let tol_f_x: f64 = 1e-6; @@ -170,7 +172,7 @@ fn test_multivariate_optimizers(verbose: bool) { let f_x_true: f64 = rosenbrock(&x_true); let sol_nelder_mead: (na::DVector, f64) = nelder_mead::nelder_mead_minimize(&(rosenbrock as fn(&na::DVector) -> f64), &na::DVector::from_vec(vec![2.0,-1.0]), 0.1, tol, max_iter, false); // num_tests_passed += check_result(x_nelder_mead, x_true, tol*1e2, "Nelder-Mead", verbose); - num_tests_passed += check_result_optim(&sol_nelder_mead.0, sol_nelder_mead.1, &x_true, f_x_true, tol_x, tol_f_x, "Nelder-Mead", verbose); + num_tests_passed += check_result_optim(&sol_nelder_mead.0, sol_nelder_mead.1, &x_true, f_x_true, tol_x, tol_f_x, "Nelder-Mead", verbose); num_tests_total += 1; print_test_results(num_tests_passed, num_tests_total); } @@ -187,7 +189,7 @@ fn test_particle_swarm_debug() { fn test_non_linear_lsqr_solvers(verbose: bool) { let mut num_tests_passed : u32 = 0; - let num_tests_total : u32 = 1; + let mut num_tests_total : u32 = 0; let tol: f64 = 1e-6; let dx_num: f64 = 1e-7; @@ -209,7 +211,7 @@ fn test_non_linear_lsqr_solvers(verbose: bool) { beta_gauss_newton = non_linear_least_squares::gauss_newton_lsqr(&xp, &yp, &fct_lsqr, &beta_gauss_newton, tol, n_iter_max, dx_num, false); // println!("beta_gauss_newton = {}\tf(beta_gauss_newton) - yp = {}", beta_gauss_newton, fct_lsqr(&xp, &beta_gauss_newton) - yp); - num_tests_passed += check_result_vector(&beta_gauss_newton, &beta_numpy, tol, "Gauss-Newton least squares", verbose); + num_tests_passed += check_result_vector(&beta_gauss_newton, &beta_numpy, tol, "Gauss-Newton least squares", verbose); num_tests_total += 1; print_test_results(num_tests_passed, num_tests_total); } diff --git a/src/univariate_solvers.rs b/src/univariate_solvers.rs index 4a100ea..a9049f8 100644 --- a/src/univariate_solvers.rs +++ b/src/univariate_solvers.rs @@ -225,4 +225,81 @@ where F : Fn(f64) -> f64 } } return Err("Maximum number of iterations exceeded.") -} \ No newline at end of file +} + +fn swap(a: &mut f64, b: &mut f64) { + let temp = *a; + *a = *b; + *b = temp; +} + +pub fn itp_solve( + f: &dyn Fn(f64) -> f64, + a: f64, + b: f64, + k1: f64, + k2: f64, + n0: f64, + epsilon: f64 +) -> f64 { + let mut a: f64 = a; + let mut b: f64 = b; + + let mut ya: f64 = f(a); + let mut yb: f64 = f(b); + + // Ensure that ya < yb + if ya > yb { + swap(&mut a, &mut b); + swap(&mut ya, &mut yb); + } + + // Preprocessing + let nhalf: f64 = f64::log2((b - a)/(2.0 * epsilon)); + let nmax: i32 = (nhalf + n0) as i32; + let mut j: i32 = 0; + + while f64::abs(b - a) > (2.0 * epsilon) { + // Calculating parameters + let x_half: f64 = (a + b) / 2.0; + let r: f64 = epsilon * (2.0 as f64).powi(nmax - j); + let delta: f64 = k1 * (b - a).powf(k2); + + // Interpolation + let xf: f64 = (yb * a - ya * b) / (yb - ya); + + // Truncation + let sigma: f64 = f64::signum(x_half - xf); + let xt: f64; + if delta <= f64::abs(x_half - xf) { + xt = xf + sigma * delta; + } else { + xt = x_half; + } + + // Projection + let x_itp: f64; + if f64::abs(xt - x_half) <= r { + x_itp = xt; + } else { + x_itp = x_half - sigma * r; + } + + // Updating interval + let y_itp: f64 = f(x_itp); + if y_itp > 0.0 { + b = x_itp; + yb = y_itp; + } else if y_itp < 0.0 { + a = x_itp; + ya = y_itp; + } else { + a = x_itp; + b = x_itp; + } + + j += 1; + } + + (a + b) / 2.0 +}