optimization/
problems.rs

1//! Common optimization problems for testing purposes.
2//!
3//! Currently, the following [optimization test functions] are implemented.
4//!
5//! ## Bowl-Shaped
6//!
7//! * [`Sphere`](http://www.sfu.ca/~ssurjano/spheref.html)
8//!
9//! ## Valley-Shaped
10//!
11//! * [`Rosenbrock`](http://www.sfu.ca/~ssurjano/rosen.html)
12//!
13//! [optimization test functions]: http://www.sfu.ca/~ssurjano/optimization.html
14
15use rand::random;
16use std::f64::INFINITY;
17use std::ops::Add;
18
19use types::{Function, Function1};
20
21
22/// Specifies a well known optimization problem.
23pub trait Problem: Function + Default {
24    /// Returns the dimensionality of the input domain.
25    fn dimensions(&self) -> usize;
26
27    /// Returns the input domain of the function in terms of upper and lower,
28    /// respectively, for each input dimension.
29    fn domain(&self) -> Vec<(f64, f64)>;
30
31    /// Returns the position as well as the value of the global minimum.
32    fn minimum(&self) -> (Vec<f64>, f64);
33
34    /// Generates a random and **feasible** position to start a minimization.
35    fn random_start(&self) -> Vec<f64>;
36
37    /// Tests whether the supplied position is legal for this function.
38    fn is_legal_position(&self, position: &[f64]) -> bool {
39        position.len() == self.dimensions() &&
40        position.iter().zip(self.domain()).all(|(&x, (lower, upper))| {
41            lower < x && x < upper
42        })
43    }
44}
45
46
47macro_rules! define_problem {
48    ( $name:ident: $this:ident,
49        default: $def:expr,
50        dimensions: $dims:expr,
51        domain: $domain:expr,
52        minimum: $miny:expr,
53        at: $minx:expr,
54        start: $start:expr,
55        value: $x1:ident => $value:expr,
56        gradient: $x2:ident => $gradient:expr ) =>
57    {
58        impl Default for $name {
59            fn default() -> Self {
60                $def
61            }
62        }
63
64        impl Function for $name {
65            fn value(&$this, $x1: &[f64]) -> f64 {
66                assert!($this.is_legal_position($x1));
67
68                $value
69            }
70        }
71
72        impl Function1 for $name {
73            fn gradient(&$this, $x2: &[f64]) -> Vec<f64> {
74                assert!($this.is_legal_position($x2));
75
76                $gradient
77            }
78        }
79
80        impl Problem for $name {
81            fn dimensions(&$this) -> usize {
82                $dims
83            }
84
85            fn domain(&$this) -> Vec<(f64, f64)> {
86                $domain
87            }
88
89            fn minimum(&$this) -> (Vec<f64>, f64) {
90                ($minx, $miny)
91            }
92
93            fn random_start(&$this) -> Vec<f64> {
94                $start
95            }
96        }
97    };
98}
99
100
101/// n-dimensional Sphere function.
102///
103/// It is continuous, convex and unimodal:
104///
105/// > f(x) = ∑ᵢ xᵢ²
106///
107/// *Global minimum*: `f(0,...,0) = 0`
108#[derive(Debug, Copy, Clone)]
109pub struct Sphere {
110    dimensions: usize
111}
112
113impl Sphere {
114    pub fn new(dimensions: usize) -> Sphere {
115        assert!(dimensions > 0, "dimensions must be larger than 1");
116
117        Sphere {
118            dimensions
119        }
120    }
121}
122
123define_problem!{Sphere: self,
124    default: Sphere::new(2),
125    dimensions: self.dimensions,
126    domain: (0..self.dimensions).map(|_| (-INFINITY, INFINITY)).collect(),
127    minimum: 0.0,
128    at: (0..self.dimensions).map(|_| 0.0).collect(),
129    start: (0..self.dimensions).map(|_| random::<f64>() * 10.24 - 5.12).collect(),
130    value: x => x.iter().map(|x| x.powi(2)).fold(0.0, Add::add),
131    gradient: x => x.iter().map(|x| 2.0 * x).collect()
132}
133
134
135/// Two-dimensional Rosenbrock function.
136///
137/// A non-convex function with its global minimum inside a long, narrow, parabolic
138/// shaped flat valley:
139///
140/// > f(x, y) = (a - x)² + b (y - x²)²
141///
142/// *Global minimum*: `f(a, a²) = 0`
143#[derive(Debug, Copy, Clone)]
144pub struct Rosenbrock {
145    a: f64,
146    b: f64
147}
148
149impl Rosenbrock {
150    /// Creates a new `Rosenbrock` function given `a` and `b`, commonly definied
151    /// with 1 and 100, respectively, which also corresponds to the `default`.
152    pub fn new(a: f64, b: f64) -> Rosenbrock {
153        Rosenbrock {
154            a,
155            b
156        }
157    }
158}
159
160define_problem!{Rosenbrock: self,
161    default: Rosenbrock::new(1.0, 100.0),
162    dimensions: 2,
163    domain: vec![(-INFINITY, INFINITY), (-INFINITY, INFINITY)],
164    minimum: 0.0,
165    at: vec![self.a, self.a * self.a],
166    start: (0..2).map(|_| random::<f64>() * 4.096 - 2.048).collect(),
167    value: x => (self.a - x[0]).powi(2) + self.b * (x[1] - x[0].powi(2)).powi(2),
168    gradient: x => vec![-2.0 * self.a + 4.0 * self.b * x[0].powi(3) - 4.0 * self.b * x[0] * x[1] + 2.0 * x[0],
169                        2.0 * self.b * (x[1] - x[0].powi(2))]
170}
171
172
173/*
174pub struct McCormick;
175
176impl McCormick {
177    pub fn new() -> McCormick {
178        McCormick
179    }
180}
181
182define_problem!{McCormick: self,
183    default: McCormick::new(),
184    dimensions: 2,
185    domain: vec![(-INFINITY, INFINITY), (-INFINITY, INFINITY)],
186    minimum: -1.9133,
187    at: vec![-0.54719, -1.54719],
188    start: vec![random::<f64>() * 5.5 - 1.5, random::<f64>() * 7.0 - 3.0],
189    value: x => (x[0] + x[1]).sin() + (x[0] - x[1]).powi(2) - 1.5 * x[0] + 2.5 * x[1] + 1.0,
190    gradient: x => vec![(x[0] + x[1]).cos() + 2.0 * (x[0] - x[1]) - 1.5,
191                        (x[0] + x[1]).cos() - 2.0 * (x[0] - x[1]) + 2.5]
192}
193*/
194
195
196#[cfg(test)]
197macro_rules! test_minimizer {
198    ( $minimizer:expr, $( $name:ident => $problem:expr ),* ) => {
199        $(
200            #[test]
201            fn $name() {
202                let minimizer = $minimizer;
203                let problem = $problem;
204
205                for _ in 0..100 {
206                    let position = $crate::problems::Problem::random_start(&problem);
207
208                    let solution = $crate::Minimizer::minimize(&minimizer,
209                        &problem, position);
210
211                    let distance = $crate::Evaluation::position(&solution).iter()
212                        .zip($crate::problems::Problem::minimum(&problem).0)
213                        .map(|(a, b)| (a - b).powi(2))
214                        .fold(0.0, ::std::ops::Add::add)
215                        .sqrt();
216
217                    assert!(distance < 1.0e-2);
218                }
219            }
220        )*
221    };
222}