diff --git a/Cargo.toml b/Cargo.toml index fcd42ca..3b24c33 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -17,7 +17,6 @@ testing = [] [dependencies] nalgebra = "0.32.0" approx = "0.5.0" -num-traits = "0.2.0" thiserror = "1.0.0" log = "0.4.0" getset = "0.1.0" diff --git a/src/algo/lipo.rs b/src/algo/lipo.rs index a35ff49..b8af58d 100644 --- a/src/algo/lipo.rs +++ b/src/algo/lipo.rs @@ -15,7 +15,6 @@ use log::{debug, trace}; use nalgebra::{ convert, try_convert, ComplexField, DimName, Dyn, IsContiguous, OVector, StorageMut, Vector, U1, }; -use num_traits::{One, Zero}; use thiserror::Error; use crate::core::{Domain, Function, Optimizer, Problem, Sample, Solver, System}; @@ -100,6 +99,7 @@ impl Lipo

{ /// Initializes LIPO solver with given options. pub fn with_options(p: &P, dom: &Domain, options: LipoOptions

, rng: Rng) -> Self { let dim = Dyn(dom.dim()); + let zero = convert(0.0); let p_explore = options.p_explore.clamp(0.0, 1.0); @@ -114,8 +114,8 @@ impl Lipo

{ xs: Vec::new(), ys: Vec::new(), best: 0, - k: P::Field::zero(), - k_inf: P::Field::zero(), + k: zero, + k_inf: zero, rng, p_explore, tmp: OVector::zeros_generic(dim, U1::name()), @@ -127,11 +127,12 @@ impl Lipo

{ /// Resets the internal state of the solver. pub fn reset(&mut self) { + let zero = convert(0.0); self.xs.clear(); self.ys.clear(); self.best = 0; - self.k = P::Field::zero(); - self.k_inf = P::Field::zero(); + self.k = zero; + self.k_inf = zero; self.iter = 0; } @@ -156,6 +157,8 @@ impl Lipo

{ .. } = self; + let one: P::Field = convert(1.0); + if !xs.is_empty() { // By definition, the k estimation is done by finding a minimum k_l // from a sequence of ks, such that @@ -179,9 +182,9 @@ impl Lipo

{ debug!("|| x - xi || = {}", dist); } - let it = try_convert(((*k_inf).ln() / (P::Field::one() + alpha).ln()).ceil()) - .unwrap_or_default() as i32; - *k = (P::Field::one() + alpha).powi(it); + let it = + try_convert(((*k_inf).ln() / (one + alpha).ln()).ceil()).unwrap_or_default() as i32; + *k = (one + alpha).powi(it); if !k.is_finite() { return Err(LipoError::InfiniteLipschitzConstant); @@ -267,7 +270,7 @@ where // Exploitation mode is allowed only when there is enough points // evaluated and the Lipschitz constant is estimated. Then there is // randomness involved in choosing whether we explore or exploit. - if !initialization && *k != F::Field::zero() && rng.f64() >= *p_explore { + if !initialization && *k != convert(0.0) && rng.f64() >= *p_explore { debug!("exploitation mode"); let mut tmp_best = ys[*best]; diff --git a/src/algo/nelder_mead.rs b/src/algo/nelder_mead.rs index ba2c08a..8b163d0 100644 --- a/src/algo/nelder_mead.rs +++ b/src/algo/nelder_mead.rs @@ -28,7 +28,6 @@ use nalgebra::{ storage::{Storage, StorageMut}, ComplexField, Dim, DimName, Dyn, IsContiguous, OVector, RealField, Vector, U1, }; -use num_traits::{One, Zero}; use thiserror::Error; use crate::core::{Domain, Function, Optimizer, Problem, RealField as _, Solver, System}; @@ -103,13 +102,14 @@ impl NelderMeadOptions

{ } CoefficientsFamily::Balanced => { let n: P::Field = convert(dom.dim() as f64); - let n_inv = P::Field::one() / n; + let one: P::Field = convert(1.0); + let n_inv = one / n; *reflection_coeff = convert(-1.0); - *expansion_coeff = -(n_inv * convert(2.0) + convert(1.0)); - *outer_contraction_coeff = -(P::Field::one() - n_inv); + *expansion_coeff = -(n_inv * convert(2.0) + one); + *outer_contraction_coeff = -(one - n_inv); *inner_contraction_coeff = -*outer_contraction_coeff; - *shrink_coeff = P::Field::one() - n_inv; + *shrink_coeff = one - n_inv; } CoefficientsFamily::GoldenSection => { let alpha = 1.0 / (0.5 * (5f64.sqrt() + 1.0)); @@ -271,7 +271,7 @@ impl NelderMead { } // Calculate the centroid. - centroid.fill(F::Field::zero()); + centroid.fill(convert(0.0)); (0..n) .map(|i| &simplex[sort_perm[i]]) .for_each(|xi| *centroid += xi); diff --git a/src/algo/trust_region.rs b/src/algo/trust_region.rs index a8c22da..6257126 100644 --- a/src/algo/trust_region.rs +++ b/src/algo/trust_region.rs @@ -33,7 +33,6 @@ use nalgebra::{ convert, storage::StorageMut, ComplexField, DimName, Dyn, IsContiguous, OMatrix, OVector, RealField, Vector, U1, }; -use num_traits::{One, Zero}; use thiserror::Error; use crate::{ @@ -150,7 +149,7 @@ impl TrustRegion

{ let delta_init = match options.delta_init { DeltaInit::Fixed(fixed) => fixed, // Zero is recognized in the function `next`. - DeltaInit::Estimated => P::Field::zero(), + DeltaInit::Estimated => convert(0.0), }; let scale = dom .scale() @@ -181,7 +180,7 @@ impl TrustRegion

{ pub fn reset(&mut self) { self.delta = match self.options.delta_init { DeltaInit::Fixed(fixed) => fixed, - DeltaInit::Estimated => P::Field::zero(), + DeltaInit::Estimated => convert(0.0), }; self.mu = convert(0.5); self.iter = 1; @@ -246,6 +245,9 @@ impl Solver for TrustRegion { .. } = self; + let one = convert(1.0); + let zero = convert(0.0); + #[allow(clippy::needless_late_init)] let scaled_newton: &mut OVector; let scale_inv2: &mut OVector; @@ -266,7 +268,7 @@ impl Solver for TrustRegion { let rx_norm = rx.norm(); - let estimate_delta = *delta == R::Field::zero(); + let estimate_delta = *delta == zero; if estimate_delta { // Zero delta signifies that the initial delta is to be set // automatically and it has not been done yet. @@ -280,8 +282,8 @@ impl Solver for TrustRegion { // where K = 100. The approach is taken from GSL. for (j, col) in jac.column_iter().enumerate() { temp[j] = col.norm(); - if temp[j] == R::Field::zero() { - temp[j] = R::Field::one(); + if temp[j] == zero { + temp[j] = one; } } temp.component_mul_assign(x); @@ -289,7 +291,7 @@ impl Solver for TrustRegion { let factor = convert(100.0); *delta = temp.norm() * factor; - if *delta == R::Field::zero() { + if *delta == zero { *delta = factor; } } @@ -311,7 +313,7 @@ impl Solver for TrustRegion { "Newton step is invalid for ill-defined Jacobian (zero columns: {:?})", jac.column_iter() .enumerate() - .filter(|(_, col)| col.norm() == R::Field::zero()) + .filter(|(_, col)| col.norm() == zero) .map(|(i, _)| i) .collect::>() ); @@ -337,7 +339,7 @@ impl Solver for TrustRegion { let grad_norm = grad_neg.norm(); - if grad_norm == R::Field::zero() { + if grad_norm == zero { // Gradient is zero, it is useless to compute the dogleg step. // Instead, we take the Newton direction to the trust region // boundary. @@ -356,7 +358,7 @@ impl Solver for TrustRegion { // Compute D^(-2) (use p for storage). scale_inv2 = p; scale_inv2.copy_from(scale); - scale_inv2.apply(|s| *s = R::Field::one() / (*s * *s)); + scale_inv2.apply(|s| *s = one / (*s * *s)); // Compute g = -D^(-2) grad F, the steepest descent direction in // scaled space (use cauchy for storage). @@ -454,7 +456,7 @@ impl Solver for TrustRegion { #[allow(clippy::suspicious_operation_groupings)] let d = (b * b + a * c_neg).sqrt(); - let alpha = if b <= R::Field::zero() { + let alpha = if b <= zero { (-b + d) / a } else { c_neg / (b + d) @@ -496,10 +498,10 @@ impl Solver for TrustRegion { // from the solution (i.e., for large || r(x) ||). // Determine lambda. - let d = if rx_norm >= R::Field::one() { - R::Field::one() / rx_norm + let d = if rx_norm >= one { + one / rx_norm } else { - R::Field::one() + R::Field::one() / convert(*iter as f64) + one + one / convert(*iter as f64) }; let lambda = *mu * rx_norm.powf(d); @@ -581,13 +583,13 @@ impl Solver for TrustRegion { let deny = if allow_ascent { // If ascent is allowed, then check only for zero, which would // make the gain ratio calculation ill-defined. - predicted == R::Field::zero() + predicted == zero } else { // If ascent is not allowed, test positivity of the predicted // gain. Note that even if the actual reduction was positive, // the step would be rejected anyway because the gain ratio // would be negative. - predicted <= R::Field::zero() + predicted <= zero }; if deny { @@ -596,7 +598,7 @@ impl Solver for TrustRegion { } else { debug!("predicted gain <= 0"); } - R::Field::zero() + zero } else { let actual = rx_norm - rx_trial_norm; let gain_ratio = actual / predicted; @@ -606,7 +608,7 @@ impl Solver for TrustRegion { } } else { debug!("trial step is invalid, gain ratio = 0"); - R::Field::zero() + zero }; // Decide if the step is accepted or not. @@ -722,6 +724,9 @@ impl Optimizer for TrustRegion { .. } = self; + let zero = convert(0.0); + let one = convert(1.0); + #[allow(clippy::needless_late_init)] let scaled_newton: &mut OVector; let scale_inv2_grad: &mut OVector; @@ -740,7 +745,7 @@ impl Optimizer for TrustRegion { grad.compute(f, x, scale, fx); hes.compute(f, x, scale, fx); - let estimate_delta = *delta == F::Field::zero(); + let estimate_delta = *delta == zero; if estimate_delta { *delta = grad.norm() * convert(0.1); } @@ -763,7 +768,7 @@ impl Optimizer for TrustRegion { "Newton step is invalid for ill-defined Hessian (zero columns: {:?})", hes.column_iter() .enumerate() - .filter(|(_, col)| col.norm() == F::Field::zero()) + .filter(|(_, col)| col.norm() == zero) .map(|(i, _)| i) .collect::>() ); @@ -786,7 +791,7 @@ impl Optimizer for TrustRegion { let grad_norm = grad_neg.norm(); - if grad_norm == F::Field::zero() { + if grad_norm == zero { // Gradient is zero, it is useless to compute the dogleg step. // Instead, we take the Newton direction to the trust region // boundary. @@ -805,13 +810,13 @@ impl Optimizer for TrustRegion { // Compute || D^-1 grad || (use p for storage). scale_inv2_grad = p; scale_inv2_grad.copy_from(scale); - scale_inv2_grad.apply(|s| *s = F::Field::one() / *s); + scale_inv2_grad.apply(|s| *s = one / *s); scale_inv2_grad.component_mul_assign(grad_neg); let scale_inv_grad_norm = scale_inv2_grad.norm(); // Compute D^(-2) grad. scale_inv2_grad.copy_from(scale); - scale_inv2_grad.apply(|s| *s = F::Field::one() / (*s * *s)); + scale_inv2_grad.apply(|s| *s = one / (*s * *s)); scale_inv2_grad.component_mul_assign(grad_neg); scale_inv2_grad.neg_mut(); @@ -825,11 +830,11 @@ impl Optimizer for TrustRegion { hes.mul_to(scale_inv2_grad, temp); let quadratic_form = scale_inv2_grad.dot(temp); - let tau = if quadratic_form <= F::Field::zero() { - F::Field::one() + let tau = if quadratic_form <= zero { + one } else { // tau = min(|| D^-1 grad|| / delta * grad^T D^-2 H D^-2 grad, 1). - (scale_inv_grad_norm.powi(3) / (*delta * quadratic_form)).min(F::Field::one()) + (scale_inv_grad_norm.powi(3) / (*delta * quadratic_form)).min(one) }; let p = scale_inv2_grad; @@ -902,7 +907,7 @@ impl Optimizer for TrustRegion { #[allow(clippy::suspicious_operation_groupings)] let d = (b * b + a * c_neg).sqrt(); - let alpha = if b <= F::Field::zero() { + let alpha = if b <= zero { (-b + d) / a } else { c_neg / (b + d) @@ -952,13 +957,13 @@ impl Optimizer for TrustRegion { let deny = if allow_ascent { // If ascent is allowed, then check only for zero, which would // make the gain ratio calculation ill-defined. - predicted == F::Field::zero() + predicted == zero } else { // If ascent is not allowed, test positivity of the predicted // gain. Note that even if the actual reduction was positive, // the step would be rejected anyway because the gain ratio // would be negative. - predicted <= F::Field::zero() + predicted <= zero }; if deny { @@ -967,7 +972,7 @@ impl Optimizer for TrustRegion { } else { debug!("predicted gain <= 0"); } - F::Field::zero() + zero } else { let actual = fx - fx_trial; let gain_ratio = actual / predicted; @@ -977,7 +982,7 @@ impl Optimizer for TrustRegion { } } else { debug!("trial step is invalid, gain ratio = 0"); - F::Field::zero() + zero }; // Decide if the step is accepted or not. diff --git a/src/derivatives.rs b/src/derivatives.rs index 5717e14..91643cb 100644 --- a/src/derivatives.rs +++ b/src/derivatives.rs @@ -3,10 +3,10 @@ use std::ops::Deref; use nalgebra::{ + convert, storage::{Storage, StorageMut}, ComplexField, DimName, Dyn, IsContiguous, OMatrix, OVector, RealField, Vector, U1, }; -use num_traits::One; use crate::core::{Function, Problem, RealField as _, System}; @@ -68,6 +68,7 @@ impl Jacobian { Srx: Storage, { let eps = R::Field::EPSILON_SQRT; + let one: R::Field = convert(1.0); for (j, mut col) in self.jac.column_iter_mut().enumerate() { let xj = x[j]; @@ -88,7 +89,7 @@ impl Jacobian { // Rosenbrock test to fail for trust region algorithm. This is very // anecdotal evidence and I would like to understand why the sign of // the step is important. - let magnitude = R::Field::one() / scale[j]; + let magnitude = one / scale[j]; let step = eps * xj.abs().max(magnitude); // Update the point. @@ -170,12 +171,13 @@ impl Gradient { Sscale: Storage, { let eps = F::Field::EPSILON_SQRT; + let one: F::Field = convert(1.0); for i in 0..x.nrows() { let xi = x[i]; // See the implementation of Jacobian for details on computing step size. - let magnitude = F::Field::one() / scale[i]; + let magnitude = one / scale[i]; let step = eps * xi.abs().max(magnitude); // Update the point. @@ -260,12 +262,13 @@ impl Hessian { Sscale: Storage, { let eps = F::Field::EPSILON_CBRT; + let one: F::Field = convert(1.0); for i in 0..x.nrows() { let xi = x[i]; // See the implementation of Jacobian for details on computing step size. - let magnitude = F::Field::one() / scale[i]; + let magnitude = one / scale[i]; let step = eps * xi.abs().max(magnitude); // Store the step for Hessian calculation.