Skip to content

Commit

Permalink
refactor: do not use num_traits dependency
Browse files Browse the repository at this point in the history
  • Loading branch information
@pnevyk
pnevyk committed Nov 15, 2023
1 parent 24791d1 commit b69bb8d
Show file tree
Hide file tree
Showing 5 changed files with 61 additions and 51 deletions.
1 change: 0 additions & 1 deletion Cargo.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down
21 changes: 12 additions & 9 deletions src/algo/lipo.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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};
Expand Down Expand Up @@ -100,6 +99,7 @@ impl<P: Problem> Lipo<P> {
/// Initializes LIPO solver with given options.
pub fn with_options(p: &P, dom: &Domain<P::Field>, options: LipoOptions<P>, rng: Rng) -> Self {
let dim = Dyn(dom.dim());
let zero = convert(0.0);

let p_explore = options.p_explore.clamp(0.0, 1.0);

Expand All @@ -114,8 +114,8 @@ impl<P: Problem> Lipo<P> {
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()),
Expand All @@ -127,11 +127,12 @@ impl<P: Problem> Lipo<P> {

/// 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;
}

Expand All @@ -156,6 +157,8 @@ impl<P: Problem> Lipo<P> {
..
} = 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
Expand All @@ -179,9 +182,9 @@ impl<P: Problem> Lipo<P> {
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);
Expand Down Expand Up @@ -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];
Expand Down
12 changes: 6 additions & 6 deletions src/algo/nelder_mead.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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};
Expand Down Expand Up @@ -103,13 +102,14 @@ impl<P: Problem> NelderMeadOptions<P> {
}
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));
Expand Down Expand Up @@ -271,7 +271,7 @@ impl<F: Function> NelderMead<F> {
}

// 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);
Expand Down
67 changes: 36 additions & 31 deletions src/algo/trust_region.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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::{
Expand Down Expand Up @@ -150,7 +149,7 @@ impl<P: Problem> TrustRegion<P> {
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()
Expand Down Expand Up @@ -181,7 +180,7 @@ impl<P: Problem> TrustRegion<P> {
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;
Expand Down Expand Up @@ -246,6 +245,9 @@ impl<R: System> Solver<R> for TrustRegion<R> {
..
} = self;

let one = convert(1.0);
let zero = convert(0.0);

#[allow(clippy::needless_late_init)]
let scaled_newton: &mut OVector<R::Field, Dyn>;
let scale_inv2: &mut OVector<R::Field, Dyn>;
Expand All @@ -266,7 +268,7 @@ impl<R: System> Solver<R> for TrustRegion<R> {

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.
Expand All @@ -280,16 +282,16 @@ impl<R: System> Solver<R> for TrustRegion<R> {
// 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);

let factor = convert(100.0);
*delta = temp.norm() * factor;

if *delta == R::Field::zero() {
if *delta == zero {
*delta = factor;
}
}
Expand All @@ -311,7 +313,7 @@ impl<R: System> Solver<R> for TrustRegion<R> {
"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::<Vec<_>>()
);
Expand All @@ -337,7 +339,7 @@ impl<R: System> Solver<R> for TrustRegion<R> {

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.
Expand All @@ -356,7 +358,7 @@ impl<R: System> Solver<R> for TrustRegion<R> {
// 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).
Expand Down Expand Up @@ -454,7 +456,7 @@ impl<R: System> Solver<R> for TrustRegion<R> {

#[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)
Expand Down Expand Up @@ -496,10 +498,10 @@ impl<R: System> Solver<R> for TrustRegion<R> {
// 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);
Expand Down Expand Up @@ -581,13 +583,13 @@ impl<R: System> Solver<R> for TrustRegion<R> {
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 {
Expand All @@ -596,7 +598,7 @@ impl<R: System> Solver<R> for TrustRegion<R> {
} else {
debug!("predicted gain <= 0");
}
R::Field::zero()
zero
} else {
let actual = rx_norm - rx_trial_norm;
let gain_ratio = actual / predicted;
Expand All @@ -606,7 +608,7 @@ impl<R: System> Solver<R> for TrustRegion<R> {
}
} else {
debug!("trial step is invalid, gain ratio = 0");
R::Field::zero()
zero
};

// Decide if the step is accepted or not.
Expand Down Expand Up @@ -722,6 +724,9 @@ impl<F: Function> Optimizer<F> for TrustRegion<F> {
..
} = self;

let zero = convert(0.0);
let one = convert(1.0);

#[allow(clippy::needless_late_init)]
let scaled_newton: &mut OVector<F::Field, Dyn>;
let scale_inv2_grad: &mut OVector<F::Field, Dyn>;
Expand All @@ -740,7 +745,7 @@ impl<F: Function> Optimizer<F> for TrustRegion<F> {
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);
}
Expand All @@ -763,7 +768,7 @@ impl<F: Function> Optimizer<F> for TrustRegion<F> {
"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::<Vec<_>>()
);
Expand All @@ -786,7 +791,7 @@ impl<F: Function> Optimizer<F> for TrustRegion<F> {

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.
Expand All @@ -805,13 +810,13 @@ impl<F: Function> Optimizer<F> for TrustRegion<F> {
// 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();

Expand All @@ -825,11 +830,11 @@ impl<F: Function> Optimizer<F> for TrustRegion<F> {
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;
Expand Down Expand Up @@ -902,7 +907,7 @@ impl<F: Function> Optimizer<F> for TrustRegion<F> {

#[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)
Expand Down Expand Up @@ -952,13 +957,13 @@ impl<F: Function> Optimizer<F> for TrustRegion<F> {
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 {
Expand All @@ -967,7 +972,7 @@ impl<F: Function> Optimizer<F> for TrustRegion<F> {
} else {
debug!("predicted gain <= 0");
}
F::Field::zero()
zero
} else {
let actual = fx - fx_trial;
let gain_ratio = actual / predicted;
Expand All @@ -977,7 +982,7 @@ impl<F: Function> Optimizer<F> for TrustRegion<F> {
}
} else {
debug!("trial step is invalid, gain ratio = 0");
F::Field::zero()
zero
};

// Decide if the step is accepted or not.
Expand Down
Loading

0 comments on commit b69bb8d

Please sign in to comment.