|
| 1 | +//! Egor optimizer implements EGO algorithm with basic handling of constraints. |
| 2 | +//! |
| 3 | +//! ```no_run |
| 4 | +//! # use ndarray::{array, Array2, ArrayView1, ArrayView2, Zip}; |
| 5 | +//! # use egobox_doe::{Lhs, SamplingMethod}; |
| 6 | +//! # use egobox_ego::{EgorBuilder, InfillStrategy, InfillOptimizer}; |
| 7 | +//! |
| 8 | +//! # use rand_xoshiro::Xoshiro256Plus; |
| 9 | +//! # use ndarray_rand::rand::SeedableRng; |
| 10 | +//! use argmin_testfunctions::rosenbrock; |
| 11 | +//! |
| 12 | +//! // Rosenbrock test function: minimum y_opt = 0 at x_opt = (1, 1) |
| 13 | +//! fn rosenb(x: &ArrayView2<f64>) -> Array2<f64> { |
| 14 | +//! let mut y: Array2<f64> = Array2::zeros((x.nrows(), 1)); |
| 15 | +//! Zip::from(y.rows_mut()) |
| 16 | +//! .and(x.rows()) |
| 17 | +//! .par_for_each(|mut yi, xi| yi.assign(&array![rosenbrock(&xi.to_vec(), 1., 100.)])); |
| 18 | +//! y |
| 19 | +//! } |
| 20 | +//! |
| 21 | +//! let xlimits = array![[-2., 2.], [-2., 2.]]; |
| 22 | +//! let res = EgorBuilder::optimize(rosenb) |
| 23 | +//! .min_within(&xlimits) |
| 24 | +//! .infill_strategy(InfillStrategy::EI) |
| 25 | +//! .n_doe(10) |
| 26 | +//! .target(1e-1) |
| 27 | +//! .n_iter(30) |
| 28 | +//! .run() |
| 29 | +//! .expect("Rosenbrock minimization"); |
| 30 | +//! println!("Rosenbrock min result = {:?}", res); |
| 31 | +//! ``` |
| 32 | +//! |
| 33 | +//! Constraints are expected to be evaluated with the objective function |
| 34 | +//! meaning that the function passed to the optimizer has to return |
| 35 | +//! a vector consisting of [obj, cstr_1, ..., cstr_n] and the cstr values |
| 36 | +//! are intended to be negative at the end of the optimization. |
| 37 | +//! Constraint number should be declared with `n_cstr` setter. |
| 38 | +//! A tolerance can be adjust with `cstr_tol` setter for relaxing constraint violation |
| 39 | +//! if specified cstr values should be < `cstr_tol` (instead of < 0) |
| 40 | +//! |
| 41 | +//! ```no_run |
| 42 | +//! # use ndarray::{array, Array2, ArrayView1, ArrayView2, Zip}; |
| 43 | +//! # use egobox_doe::{Lhs, SamplingMethod}; |
| 44 | +//! # use egobox_ego::{EgorBuilder, InfillStrategy, InfillOptimizer}; |
| 45 | +//! # use rand_xoshiro::Xoshiro256Plus; |
| 46 | +//! # use ndarray_rand::rand::SeedableRng; |
| 47 | +//! |
| 48 | +//! // Function G24: 1 global optimum y_opt = -5.5080 at x_opt =(2.3295, 3.1785) |
| 49 | +//! fn g24(x: &ArrayView1<f64>) -> f64 { |
| 50 | +//! -x[0] - x[1] |
| 51 | +//! } |
| 52 | +//! |
| 53 | +//! // Constraints < 0 |
| 54 | +//! fn g24_c1(x: &ArrayView1<f64>) -> f64 { |
| 55 | +//! -2.0 * x[0].powf(4.0) + 8.0 * x[0].powf(3.0) - 8.0 * x[0].powf(2.0) + x[1] - 2.0 |
| 56 | +//! } |
| 57 | +//! |
| 58 | +//! fn g24_c2(x: &ArrayView1<f64>) -> f64 { |
| 59 | +//! -4.0 * x[0].powf(4.0) + 32.0 * x[0].powf(3.0) |
| 60 | +//! - 88.0 * x[0].powf(2.0) + 96.0 * x[0] + x[1] |
| 61 | +//! - 36.0 |
| 62 | +//! } |
| 63 | +//! |
| 64 | +//! // Gouped function : objective + constraints |
| 65 | +//! fn f_g24(x: &ArrayView2<f64>) -> Array2<f64> { |
| 66 | +//! let mut y = Array2::zeros((x.nrows(), 3)); |
| 67 | +//! Zip::from(y.rows_mut()) |
| 68 | +//! .and(x.rows()) |
| 69 | +//! .for_each(|mut yi, xi| { |
| 70 | +//! yi.assign(&array![g24(&xi), g24_c1(&xi), g24_c2(&xi)]); |
| 71 | +//! }); |
| 72 | +//! y |
| 73 | +//! } |
| 74 | +//! |
| 75 | +//! let xlimits = array![[0., 3.], [0., 4.]]; |
| 76 | +//! let doe = Lhs::new(&xlimits).sample(10); |
| 77 | +//! let res = EgorBuilder::optimize(f_g24) |
| 78 | +//! .min_within(&xlimits) |
| 79 | +//! .n_cstr(2) |
| 80 | +//! .infill_strategy(InfillStrategy::EI) |
| 81 | +//! .infill_optimizer(InfillOptimizer::Cobyla) |
| 82 | +//! .doe(&doe) |
| 83 | +//! .n_iter(40) |
| 84 | +//! .target(-5.5080) |
| 85 | +//! .run() |
| 86 | +//! .expect("g24 minimized"); |
| 87 | +//! println!("G24 min result = {:?}", res); |
| 88 | +//! ``` |
| 89 | +//! |
| 90 | +use crate::egor_solver::*; |
| 91 | +use crate::mixint::*; |
| 92 | +use crate::types::*; |
| 93 | + |
| 94 | +use egobox_moe::{CorrelationSpec, MoeParams, RegressionSpec}; |
| 95 | +use ndarray::Array1; |
| 96 | +use ndarray::{Array2, ArrayBase, Data, Ix2}; |
| 97 | +use ndarray_rand::rand::SeedableRng; |
| 98 | +use rand_xoshiro::Xoshiro256Plus; |
| 99 | + |
| 100 | +/// EGO optimizer service builder allowing to use Egor optimizer |
| 101 | +/// with an ask-and-tell interface. |
| 102 | +/// |
| 103 | +pub struct EgorServiceBuilder { |
| 104 | + seed: Option<u64>, |
| 105 | +} |
| 106 | + |
| 107 | +impl EgorServiceBuilder { |
| 108 | + /// Function to be minimized domain should be basically R^nx -> R^ny |
| 109 | + /// where nx is the dimension of input x and ny the output dimension |
| 110 | + /// equal to 1 (obj) + n (cstrs). |
| 111 | + /// But function has to be able to evaluate several points in one go |
| 112 | + /// hence take an (p, nx) matrix and return an (p, ny) matrix |
| 113 | + pub fn optimize() -> Self { |
| 114 | + EgorServiceBuilder { seed: None } |
| 115 | + } |
| 116 | + |
| 117 | + /// Allow to specify a seed for random number generator to allow |
| 118 | + /// reproducible runs. |
| 119 | + pub fn random_seed(mut self, seed: u64) -> Self { |
| 120 | + self.seed = Some(seed); |
| 121 | + self |
| 122 | + } |
| 123 | + |
| 124 | + /// Build an Egor optimizer to minimize the function within |
| 125 | + /// the continuous `xlimits` specified as [[lower, upper], ...] array where the |
| 126 | + /// number of rows gives the dimension of the inputs (continuous optimization) |
| 127 | + /// and the ith row is the interval of the ith component of the input x. |
| 128 | + pub fn min_within( |
| 129 | + self, |
| 130 | + xlimits: &ArrayBase<impl Data<Elem = f64>, Ix2>, |
| 131 | + ) -> EgorService<MoeParams<f64, Xoshiro256Plus>> { |
| 132 | + let rng = if let Some(seed) = self.seed { |
| 133 | + Xoshiro256Plus::seed_from_u64(seed) |
| 134 | + } else { |
| 135 | + Xoshiro256Plus::from_entropy() |
| 136 | + }; |
| 137 | + EgorService { |
| 138 | + solver: EgorSolver::new(xlimits, rng), |
| 139 | + } |
| 140 | + } |
| 141 | + |
| 142 | + /// Build an Egor optimizer to minimize the function R^n -> R^p taking |
| 143 | + /// inputs specified with given xtypes where some of components may be |
| 144 | + /// discrete variables (mixed-integer optimization). |
| 145 | + pub fn min_within_mixint_space(self, xtypes: &[XType]) -> EgorService<MixintMoeParams> { |
| 146 | + let rng = if let Some(seed) = self.seed { |
| 147 | + Xoshiro256Plus::seed_from_u64(seed) |
| 148 | + } else { |
| 149 | + Xoshiro256Plus::from_entropy() |
| 150 | + }; |
| 151 | + EgorService { |
| 152 | + solver: EgorSolver::new_with_xtypes(xtypes, rng), |
| 153 | + } |
| 154 | + } |
| 155 | +} |
| 156 | + |
| 157 | +/// Egor optimizer structure used to parameterize the underlying `argmin::Solver` |
| 158 | +/// and trigger the optimization using `argmin::Executor`. |
| 159 | +#[derive(Clone)] |
| 160 | +pub struct EgorService<SB: SurrogateBuilder> { |
| 161 | + solver: EgorSolver<SB>, |
| 162 | +} |
| 163 | + |
| 164 | +impl<SB: SurrogateBuilder> EgorService<SB> { |
| 165 | + /// Sets allowed number of evaluation of the function under optimization |
| 166 | + pub fn n_iter(mut self, n_iter: usize) -> Self { |
| 167 | + self.solver = self.solver.n_iter(n_iter); |
| 168 | + self |
| 169 | + } |
| 170 | + |
| 171 | + /// Sets the number of runs of infill strategy optimizations (best result taken) |
| 172 | + pub fn n_start(mut self, n_start: usize) -> Self { |
| 173 | + self.solver = self.solver.n_start(n_start); |
| 174 | + self |
| 175 | + } |
| 176 | + |
| 177 | + /// Sets Number of parallel evaluations of the function under optimization |
| 178 | + pub fn q_points(mut self, q_points: usize) -> Self { |
| 179 | + self.solver = self.solver.q_points(q_points); |
| 180 | + self |
| 181 | + } |
| 182 | + |
| 183 | + /// Number of samples of initial LHS sampling (used when DOE not provided by the user) |
| 184 | + /// |
| 185 | + /// When 0 a number of points is computed automatically regarding the number of input variables |
| 186 | + /// of the function under optimization. |
| 187 | + pub fn n_doe(mut self, n_doe: usize) -> Self { |
| 188 | + self.solver = self.solver.n_doe(n_doe); |
| 189 | + self |
| 190 | + } |
| 191 | + |
| 192 | + /// Sets the number of constraint functions |
| 193 | + pub fn n_cstr(mut self, n_cstr: usize) -> Self { |
| 194 | + self.solver = self.solver.n_cstr(n_cstr); |
| 195 | + self |
| 196 | + } |
| 197 | + |
| 198 | + /// Sets the tolerance on constraints violation (cstr < tol) |
| 199 | + pub fn cstr_tol(mut self, tol: &Array1<f64>) -> Self { |
| 200 | + self.solver = self.solver.cstr_tol(tol); |
| 201 | + self |
| 202 | + } |
| 203 | + |
| 204 | + /// Sets an initial DOE \['ns', `nt`\] containing `ns` samples. |
| 205 | + /// |
| 206 | + /// Either `nt` = `nx` then only `x` input values are specified and `ns` evals are done to get y ouput doe values, |
| 207 | + /// or `nt = nx + ny` then `x = doe\[:, :nx\]` and `y = doe\[:, nx:\]` are specified |
| 208 | + pub fn doe(mut self, doe: &Array2<f64>) -> Self { |
| 209 | + self.solver = self.solver.doe(doe); |
| 210 | + self |
| 211 | + } |
| 212 | + |
| 213 | + /// Sets the parallel infill strategy |
| 214 | + /// |
| 215 | + /// Parallel infill criterion to get virtual next promising points in order to allow |
| 216 | + /// n parallel evaluations of the function under optimization. |
| 217 | + pub fn qei_strategy(mut self, q_ei: QEiStrategy) -> Self { |
| 218 | + self.solver = self.solver.qei_strategy(q_ei); |
| 219 | + self |
| 220 | + } |
| 221 | + |
| 222 | + /// Sets the infill strategy |
| 223 | + pub fn infill_strategy(mut self, infill: InfillStrategy) -> Self { |
| 224 | + self.solver = self.solver.infill_strategy(infill); |
| 225 | + self |
| 226 | + } |
| 227 | + |
| 228 | + /// Sets the infill optimizer |
| 229 | + pub fn infill_optimizer(mut self, optimizer: InfillOptimizer) -> Self { |
| 230 | + self.solver = self.solver.infill_optimizer(optimizer); |
| 231 | + self |
| 232 | + } |
| 233 | + |
| 234 | + /// Sets the allowed regression models used in gaussian processes. |
| 235 | + pub fn regression_spec(mut self, regression_spec: RegressionSpec) -> Self { |
| 236 | + self.solver = self.solver.regression_spec(regression_spec); |
| 237 | + self |
| 238 | + } |
| 239 | + |
| 240 | + /// Sets the allowed correlation models used in gaussian processes. |
| 241 | + pub fn correlation_spec(mut self, correlation_spec: CorrelationSpec) -> Self { |
| 242 | + self.solver = self.solver.correlation_spec(correlation_spec); |
| 243 | + self |
| 244 | + } |
| 245 | + |
| 246 | + /// Sets the number of components to be used specifiying PLS projection is used (a.k.a KPLS method). |
| 247 | + /// |
| 248 | + /// This is used to address high-dimensional problems typically when `nx` > 9 wher `nx` is the dimension of `x`. |
| 249 | + pub fn kpls_dim(mut self, kpls_dim: usize) -> Self { |
| 250 | + self.solver = self.solver.kpls_dim(kpls_dim); |
| 251 | + self |
| 252 | + } |
| 253 | + |
| 254 | + /// Sets the number of clusters used by the mixture of surrogate experts. |
| 255 | + /// |
| 256 | + /// When set to 0, the number of clusters is determined automatically |
| 257 | + /// (warning in this case the optimizer runs slower) |
| 258 | + pub fn n_clusters(mut self, n_clusters: usize) -> Self { |
| 259 | + self.solver = self.solver.n_clusters(n_clusters); |
| 260 | + self |
| 261 | + } |
| 262 | + |
| 263 | + /// Sets a known target minimum to be used as a stopping criterion. |
| 264 | + pub fn target(mut self, target: f64) -> Self { |
| 265 | + self.solver = self.solver.target(target); |
| 266 | + self |
| 267 | + } |
| 268 | + |
| 269 | + /// Sets a directory to write optimization history and used as search path for hot start doe |
| 270 | + pub fn outdir(mut self, outdir: impl Into<String>) -> Self { |
| 271 | + self.solver = self.solver.outdir(outdir); |
| 272 | + self |
| 273 | + } |
| 274 | + |
| 275 | + /// Whether we start by loading last DOE saved in `outdir` as initial DOE |
| 276 | + pub fn hot_start(mut self, hot_start: bool) -> Self { |
| 277 | + self.solver = self.solver.hot_start(hot_start); |
| 278 | + self |
| 279 | + } |
| 280 | + |
| 281 | + /// Given an evaluated doe (x, y) data, return the next promising x point |
| 282 | + /// where optimum may occurs regarding the infill criterium. |
| 283 | + /// This function inverse the control of the optimization and can used |
| 284 | + /// ask-and-tell interface to the EGO optimizer. |
| 285 | + pub fn suggest( |
| 286 | + &self, |
| 287 | + x_data: &ArrayBase<impl Data<Elem = f64>, Ix2>, |
| 288 | + y_data: &ArrayBase<impl Data<Elem = f64>, Ix2>, |
| 289 | + ) -> Array2<f64> { |
| 290 | + self.solver.suggest(x_data, y_data) |
| 291 | + } |
| 292 | +} |
| 293 | + |
| 294 | +#[cfg(test)] |
| 295 | +mod tests { |
| 296 | + use super::*; |
| 297 | + use approx::assert_abs_diff_eq; |
| 298 | + use ndarray::{array, concatenate, ArrayView2, Axis}; |
| 299 | + |
| 300 | + use ndarray_stats::QuantileExt; |
| 301 | + |
| 302 | + use serial_test::serial; |
| 303 | + |
| 304 | + fn xsinx(x: &ArrayView2<f64>) -> Array2<f64> { |
| 305 | + (x - 3.5) * ((x - 3.5) / std::f64::consts::PI).mapv(|v| v.sin()) |
| 306 | + } |
| 307 | + |
| 308 | + #[test] |
| 309 | + #[serial] |
| 310 | + fn test_xsinx_egor_builder() { |
| 311 | + let ego = EgorServiceBuilder::optimize() |
| 312 | + .random_seed(42) |
| 313 | + .min_within(&array![[0., 25.]]) |
| 314 | + .regression_spec(RegressionSpec::ALL) |
| 315 | + .correlation_spec(CorrelationSpec::ALL) |
| 316 | + .infill_strategy(InfillStrategy::EI); |
| 317 | + |
| 318 | + let mut doe = array![[0.], [7.], [20.], [25.]]; |
| 319 | + let mut y_doe = xsinx(&doe.view()); |
| 320 | + for _i in 0..10 { |
| 321 | + let x_suggested = ego.suggest(&doe, &y_doe); |
| 322 | + |
| 323 | + doe = concatenate![Axis(0), doe, x_suggested]; |
| 324 | + y_doe = xsinx(&doe.view()); |
| 325 | + } |
| 326 | + |
| 327 | + let expected = -15.1; |
| 328 | + let y_opt = y_doe.min().unwrap(); |
| 329 | + assert_abs_diff_eq!(expected, *y_opt, epsilon = 1e-1); |
| 330 | + } |
| 331 | +} |
0 commit comments