diff --git a/README.md b/README.md index a4776e8..ab56723 100644 --- a/README.md +++ b/README.md @@ -31,5 +31,5 @@ let args = Array::from_vec(vec![3.0, -8.3]); let ans = minimizer.minimize(&function, args.view()); // Print the optimized values -println!("Final optimized arguments: {}", ans); +println!("Final optimized arguments: {:?}", ans); ``` diff --git a/examples/simple.rs b/examples/simple.rs index 6bed106..6adbfd9 100644 --- a/examples/simple.rs +++ b/examples/simple.rs @@ -3,6 +3,7 @@ extern crate optimize; use ndarray::prelude::*; +use optimize::minimizer::Minimizer; use optimize::vector::NelderMeadBuilder; pub fn main() { @@ -11,10 +12,10 @@ pub fn main() { let minimizer = NelderMeadBuilder::default() .xtol(1e-6f64) .ftol(1e-6f64) - .maxiter(50000) + .maxiter(50000) .build() .unwrap(); let args = Array::from_vec(vec![3.0, -8.3]); let ans = minimizer.minimize(&function, args.view()); - println!("Final optimized arguments: {}", ans); + println!("Final optimized arguments: {:?}", ans); } diff --git a/src/lib.rs b/src/lib.rs index f8be16d..6dde80c 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -11,6 +11,7 @@ //! # extern crate optimize; //! # use ndarray::prelude::*; //! # use optimize::vector::NelderMeadBuilder; +//! # use optimize::minimizer::Minimizer; //! // Define a function that we aim to minimize //! let function = |x: ArrayView1| (1.0 - x[0]).powi(2) + 100.0 * (x[1] - x[0].powi(2)).powi(2); //! @@ -29,7 +30,7 @@ //! let ans = minimizer.minimize(&function, args.view()); //! //! // Print the optimized values -//! println!("Final optimized arguments: {}", ans); +//! println!("Final optimized arguments: {}", ans.minimum_value.unwrap()); //! ``` #![deny(missing_docs)] @@ -43,7 +44,8 @@ extern crate derive_builder; extern crate float_cmp; extern crate num_traits; +pub mod minimizer; pub mod scalar; pub mod vector; -mod utils; \ No newline at end of file +mod utils; diff --git a/src/minimizer.rs b/src/minimizer.rs index 3ab348d..029a68d 100644 --- a/src/minimizer.rs +++ b/src/minimizer.rs @@ -1,5 +1,42 @@ - +//! This module provides the base framework for all minimizers present in this crate, such as the +//! base trait and return type. use ndarray::prelude::*; +use std::time::Duration; + +/// Minimizer states at the end of the run +#[derive(Debug, PartialEq)] +pub enum RunStatus { + /// Minimizer finished successfully. + Success, + /// Minimization was not succesful + Failure, +} + +/// Final minimization parameters that result in the smallest found value. +#[derive(Debug, PartialEq)] +pub enum MinResult { + /// Scalar parameter value where the minimum is found + Scalar(f64), + /// Vector of parameters where the minimum if found + Vector(Array1), +} + +/// A minimization result, storing various details of the run and the final results. +#[derive(Debug, PartialEq)] +pub struct OptimResult { + /// The runtime of the minimization according to the system clock. + pub runtime: Option, + /// The number of function evaluations performed. + pub f_evals: Option, + /// The number of iterations run. + pub iterations: Option, + /// The final parameter values. + pub minimum: Option, + /// The function value at the found minimum. + pub minimum_value: Option, + /// The minimizer success or failure status. + pub status: RunStatus, +} /// A general minimizer trait. pub trait Minimizer { @@ -9,5 +46,5 @@ pub trait Minimizer { &self, func: F, args: ArrayView1, - ) -> Array1; + ) -> OptimResult; } diff --git a/src/scalar/mod.rs b/src/scalar/mod.rs index 88edb23..f668006 100644 --- a/src/scalar/mod.rs +++ b/src/scalar/mod.rs @@ -1 +1 @@ -//! This module contains algorithms that search for local minima of functions along a single dimension. \ No newline at end of file +//! This module contains algorithms that search for local minima of functions along a single dimension. diff --git a/src/utils.rs b/src/utils.rs index fb56539..419dc81 100644 --- a/src/utils.rs +++ b/src/utils.rs @@ -31,7 +31,6 @@ where grad } - #[cfg(test)] mod tests { diff --git a/src/vector/mod.rs b/src/vector/mod.rs index 93975fa..8be0c13 100644 --- a/src/vector/mod.rs +++ b/src/vector/mod.rs @@ -2,4 +2,4 @@ mod nelder_mead; -pub use self::nelder_mead::NelderMeadBuilder; \ No newline at end of file +pub use self::nelder_mead::NelderMeadBuilder; diff --git a/src/vector/nelder_mead.rs b/src/vector/nelder_mead.rs index 81c7ad6..fcd6734 100644 --- a/src/vector/nelder_mead.rs +++ b/src/vector/nelder_mead.rs @@ -1,27 +1,28 @@ //! This implementation of Nelder-Mead is based on -//! -//! Gao, F and Han, L. Implementing the Nelder-Mead simplex algorithm with -//! adaptive parameters. 2012. Computational Optimization and Applications. +//! +//! Gao, F and Han, L. Implementing the Nelder-Mead simplex algorithm with +//! adaptive parameters. 2012. Computational Optimization and Applications. //! 51:1, pp 259--277 -//! +//! //! In particular, it adapts their suggestion to use adaptive step sizes, //! which depend on the dimensionality of the optimization problem. -//! +//! //! # Use case -//! +//! //! The Nelder-Mead algorithm does not require a gradient or a hessian. -//! As a tradeoff it typically requires a lot of function evaluations to -//! find a minimum. Further, there are few theoretical results on the +//! As a tradeoff it typically requires a lot of function evaluations to +//! find a minimum. Further, there are few theoretical results on the //! convergence of Nelder-Mead iterations. -//! +//! //! # Examples -//! +//! //! ``` //! # extern crate ndarray; //! # extern crate optimize; //! # use ndarray::prelude::*; //! # use optimize::vector::NelderMeadBuilder; -//! +//! # use optimize::minimizer::Minimizer; +//! //! let function = //! |x: ArrayView1| (1.0 - x[0]).powi(2) + 100.0 * (x[1] - x[0].powi(2)).powi(2); //! let minimizer = NelderMeadBuilder::default() @@ -31,16 +32,17 @@ //! .unwrap(); //! let args = Array::from_vec(vec![3.0, -8.3]); //! let res = minimizer.minimize(&function, args.view()); -//! println!("res: {}", res); +//! println!("res: {:?}", res); //! ``` use float_cmp::ApproxOrdUlps; +use minimizer::{MinResult, Minimizer, OptimResult, RunStatus}; use ndarray::prelude::*; -use ::utils::WrappedFunction; +use std::time::SystemTime; +use utils::WrappedFunction; type Simplex = Vec<(f64, Array1)>; - #[derive(Builder, Debug)] /// A minimizer for a scalar function of one or more variables using the Nelder-Mead algorithm. pub struct NelderMead { @@ -75,31 +77,22 @@ pub struct NelderMead { } impl NelderMead { - - /// Search for the value minimizing `func` given an initial guess - /// in the form of a point. The algorithm will explore the variable - /// space without constraints. - pub fn minimize(&self, func: F, x0: ArrayView1) -> Array1 - where F: Fn(ArrayView1) -> f64 { - let n = x0.len(); - let eps = 0.05; - let mut init_simplex = Array2::default((n+1, n)); - init_simplex.slice_mut(s![0, ..]).assign(&x0); - init_simplex.slice_mut(s![1.., ..]).assign(&(Array::eye(n) * eps + &x0 * (1.0-eps))); - - self.minimize_simplex(func, init_simplex) - } - /// Search for the value minimizing `func` given an initial guess /// in the form of a set of coordinates, the `init_simplex`. This algorithm /// only ever explores the space spanned by these initial vectors. - /// If you have parameter restrictions that effectively place your parameters + /// If you have parameter restrictions that effectively place your parameters /// in a subspace, you can enforce these restrictions by setting `init_simplex` /// to a basis of this subspace. - pub fn minimize_simplex(&self, func: F, init_simplex: Array2) -> Array1 - where F: Fn(ArrayView1) -> f64 { + pub fn minimize_simplex(&self, func: F, init_simplex: Array2) -> OptimResult + where + F: Fn(ArrayView1) -> f64, + { + let start_time = SystemTime::now(); let mut func = WrappedFunction { num: 0, func: func }; - let mut simplex = init_simplex.outer_iter().map(|xi| (func.call(xi), xi.to_owned())).collect::(); + let mut simplex = init_simplex + .outer_iter() + .map(|xi| (func.call(xi), xi.to_owned())) + .collect::(); self.order_simplex(&mut simplex); let mut centroid = self.centroid(&simplex); let n = simplex.len(); @@ -109,17 +102,18 @@ impl NelderMead { let mut iterations = 1; while !self.finished(&simplex, iterations, maxiter, func.num, maxfun) { - - let f_n1 = simplex[n-1].0; - let f_n = simplex[n-2].0; + let f_n1 = simplex[n - 1].0; + let f_n = simplex[n - 2].0; let f_0 = simplex[0].0; - let reflected = ¢roid + &(alpha * &(¢roid - &simplex[n-1].1)); + let reflected = ¢roid + &(alpha * &(¢roid - &simplex[n - 1].1)); let f_reflected = func.call(reflected.view()); - if f_reflected < f_n && f_reflected > f_0 { // try reflecting the worst point through the centroid + if f_reflected < f_n && f_reflected > f_0 { + // try reflecting the worst point through the centroid self.lean_update(&mut simplex, &mut centroid, reflected, f_reflected); - } else if f_reflected < f_0 { // try expanding beyond the centroid + } else if f_reflected < f_0 { + // try expanding beyond the centroid let expanded = ¢roid + &(beta * &(&reflected - ¢roid)); let f_expanded = func.call(expanded.view()); @@ -128,8 +122,9 @@ impl NelderMead { } else { self.lean_update(&mut simplex, &mut centroid, reflected, f_reflected); } - } else if f_reflected < f_n1 && f_reflected >= f_n { // try a contraction outwards - let contracted = ¢roid + &(gamma * (¢roid - &simplex[n-1].1)); + } else if f_reflected < f_n1 && f_reflected >= f_n { + // try a contraction outwards + let contracted = ¢roid + &(gamma * (¢roid - &simplex[n - 1].1)); let f_contracted = func.call(contracted.view()); if f_contracted < f_reflected { self.lean_update(&mut simplex, &mut centroid, contracted, f_contracted); @@ -137,8 +132,9 @@ impl NelderMead { // shrink self.shrink(&mut simplex, &mut func, delta, &mut centroid); } - } else { // try a contraction inwards - let contracted = ¢roid - &(gamma * (¢roid - &simplex[n-1].1)); + } else { + // try a contraction inwards + let contracted = ¢roid - &(gamma * (¢roid - &simplex[n - 1].1)); let f_contracted = func.call(contracted.view()); if f_contracted < f_reflected { @@ -150,7 +146,15 @@ impl NelderMead { } iterations += 1; } - simplex.remove(0).1 + let best = simplex.remove(0); + OptimResult { + f_evals: Some(func.num), + iterations: Some(iterations), + minimum: Some(MinResult::Vector(best.1)), + minimum_value: Some(best.0), + status: RunStatus::Success, + runtime: start_time.elapsed().ok(), + } } /// Helper function to keep the main loop clean. Resolves default values that can @@ -159,16 +163,21 @@ impl NelderMead { fn initialize_parameters(&self, n: usize) -> (usize, usize, f64, f64, f64, f64) { let maxiter = match self.maxiter { Some(x) => x, - None => 200 * n + None => 200 * n, }; let maxfun = match self.maxfun { Some(x) => x, - None => 200 * n + None => 200 * n, }; let (alpha, beta, gamma, delta) = if self.adaptive { let dim = n as f64; - (1.0, 1.0 + 2.0 / dim, 0.75 - 1.0 / (2.0 * dim), 1.0 - 1.0 / dim) + ( + 1.0, + 1.0 + 2.0 / dim, + 0.75 - 1.0 / (2.0 * dim), + 1.0 - 1.0 / dim, + ) } else { (1.0, 2.0, 0.5, 0.5) }; @@ -177,12 +186,22 @@ impl NelderMead { } #[inline] - fn finished(&self, simplex: &Simplex, iterations: usize, maxiter: usize, nfeval: usize, maxfun: usize) -> bool { + fn finished( + &self, + simplex: &Simplex, + iterations: usize, + maxiter: usize, + nfeval: usize, + maxfun: usize, + ) -> bool { let n = simplex.len(); - iterations > maxiter - || nfeval > maxfun - || ( simplex[n-1].0 - simplex[0].0 < self.ftol - && (&simplex[n-1].1 - &simplex[0].1).mapv(f64::abs).scalar_sum() < n as f64 * self.xtol ) + iterations > maxiter + || nfeval > maxfun + || (simplex[n - 1].0 - simplex[0].0 < self.ftol + && (&simplex[n - 1].1 - &simplex[0].1) + .mapv(f64::abs) + .scalar_sum() + < n as f64 * self.xtol) } /// Update the centroid effiently, knowing only one value changed. @@ -190,23 +209,36 @@ impl NelderMead { /// given that we inserted a single out-of-place value in a sorted vec. /// This update is O(n). #[inline] - fn lean_update(&self, simplex: &mut Simplex, centroid: &mut Array1, xnew: Array1, fnew: f64) { + fn lean_update( + &self, + simplex: &mut Simplex, + centroid: &mut Array1, + xnew: Array1, + fnew: f64, + ) { let n = simplex.len(); - *centroid += &(&xnew / (n-1) as f64); - simplex[n-1] = (fnew, xnew); + *centroid += &(&xnew / (n - 1) as f64); + simplex[n - 1] = (fnew, xnew); self.order_simplex(simplex); - *centroid -= &(&simplex[n-1].1 / (n-1) as f64); + *centroid -= &(&simplex[n - 1].1 / (n - 1) as f64); } /// shrink all points towards the best point. - /// Assumes the simplex is ordered. + /// Assumes the simplex is ordered. /// The centroid is updated by shrinking the centroid directly, /// Then removing the new 'worst x' and adding in the old 'worst x'. - /// This update of `centroid` is O(n). + /// This update of `centroid` is O(n). /// Shrinkage requires n function evaluations. #[inline] - fn shrink(&self, simplex: &mut Simplex, f: &mut WrappedFunction, sigma: f64, centroid: &mut Array1) - where F: Fn(ArrayView1) -> f64 { + fn shrink( + &self, + simplex: &mut Simplex, + f: &mut WrappedFunction, + sigma: f64, + centroid: &mut Array1, + ) where + F: Fn(ArrayView1) -> f64, + { { let mut iter = simplex.iter_mut(); let (_, x0) = iter.next().unwrap(); @@ -230,10 +262,10 @@ impl NelderMead { fn centroid(&self, simplex: &Simplex) -> Array1 { let n = simplex.len(); let mut centroid = Array1::zeros(simplex[0].1.len()); - for (_, xi) in simplex.iter().take(n-1) { + for (_, xi) in simplex.iter().take(n - 1) { centroid += xi; } - centroid / (n-1) as f64 + centroid / (n - 1) as f64 } /// This sorting algorithm should have a runtime of O(n) if only one new element is inserted. @@ -244,6 +276,26 @@ impl NelderMead { } } +impl Minimizer for NelderMead { + /// Search for the value minimizing `func` given an initial guess + /// in the form of a point. The algorithm will explore the variable + /// space without constraints. + fn minimize(&self, func: F, x0: ArrayView1) -> OptimResult + where + F: Fn(ArrayView1) -> f64, + { + let n = x0.len(); + let eps = 0.05; + let mut init_simplex = Array2::default((n + 1, n)); + init_simplex.slice_mut(s![0, ..]).assign(&x0); + init_simplex + .slice_mut(s![1.., ..]) + .assign(&(Array::eye(n) * eps + &x0 * (1.0 - eps))); + + self.minimize_simplex(func, init_simplex) + } +} + #[cfg(test)] mod tests { @@ -261,9 +313,17 @@ mod tests { .unwrap(); let args = Array::from_vec(vec![3.0, -8.3]); let res = minimizer.minimize(&function, args.view()); - println!("res: {}", res); - assert!(res[0].approx_eq(&1.0, 1e-4, 10)); - assert!(res[1].approx_eq(&1.0, 1e-4, 10)); + println!("res: {:?}", res); + assert!(res.status == RunStatus::Success); + + let values = match res.minimum.unwrap() { + MinResult::Vector(v) => Some(v), + _ => None, + }; + assert!(values.is_some()); + let values = values.unwrap(); + assert!(values[0].approx_eq(&1.0, 1e-4, 10)); + assert!(values[1].approx_eq(&1.0, 1e-4, 10)); } }