Add ITP solve method.
This commit is contained in:
parent
528f7900bf
commit
35a3705763
2 changed files with 98 additions and 19 deletions
38
src/main.rs
38
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>, f64) = nelder_mead::nelder_mead_minimize(&(rosenbrock as fn(&na::DVector<f64>) -> 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);
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -225,4 +225,81 @@ where F : Fn(f64) -> f64
|
|||
}
|
||||
}
|
||||
return Err("Maximum number of iterations exceeded.")
|
||||
}
|
||||
}
|
||||
|
||||
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
|
||||
}
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue