diff --git a/ego/Cargo.toml b/ego/Cargo.toml index e3378cd6..b381ab80 100644 --- a/ego/Cargo.toml +++ b/ego/Cargo.toml @@ -64,3 +64,7 @@ criterion = "0.5" approx = "0.4" argmin_testfunctions = "0.2" serial_test = "3.1.0" + +[[bench]] +name = "ego" +harness = false diff --git a/ego/benches/ego.rs b/ego/benches/ego.rs new file mode 100644 index 00000000..72b7b4f6 --- /dev/null +++ b/ego/benches/ego.rs @@ -0,0 +1,42 @@ +use criterion::{black_box, criterion_group, criterion_main, Criterion}; +use egobox_ego::{EgorBuilder, InfillStrategy}; +use egobox_moe::{CorrelationSpec, RegressionSpec}; +use ndarray::{array, Array2, ArrayView2, Zip}; + +/// Ackley test function: min f(x)=0 at x=(0, 0, 0) +fn ackley(x: &ArrayView2) -> Array2 { + let mut y: Array2 = Array2::zeros((x.nrows(), 1)); + Zip::from(y.rows_mut()) + .and(x.rows()) + .par_for_each(|mut yi, xi| yi.assign(&array![argmin_testfunctions::ackley(&xi.to_vec(),)])); + y +} + +fn criterion_ego(c: &mut Criterion) { + let xlimits = array![[-32.768, 32.768], [-32.768, 32.768], [-32.768, 32.768]]; + let mut group = c.benchmark_group("ego"); + group.bench_function("ego ackley", |b| { + b.iter(|| { + black_box( + EgorBuilder::optimize(ackley) + .configure(|config| { + config + .regression_spec(RegressionSpec::CONSTANT) + .correlation_spec(CorrelationSpec::ABSOLUTEEXPONENTIAL) + .infill_strategy(InfillStrategy::WB2S) + .max_iters(10) + .target(5e-1) + .seed(42) + }) + .min_within(&xlimits) + .run() + .expect("Minimize failure"), + ) + }); + }); + + group.finish(); +} + +criterion_group!(benches, criterion_ego); +criterion_main!(benches); diff --git a/ego/examples/g24.rs b/ego/examples/g24.rs index 1632e28f..abab4729 100644 --- a/ego/examples/g24.rs +++ b/ego/examples/g24.rs @@ -34,7 +34,7 @@ fn main() { // We use Egor optimizer as a service let egor = EgorServiceBuilder::optimize() - .configure(|config| config.n_cstr(2).random_seed(42)) + .configure(|config| config.n_cstr(2).seed(42)) .min_within(&xlimits); let mut y_doe = f_g24(&doe.view()); diff --git a/ego/src/egor.rs b/ego/src/egor.rs index 44df8ed2..c836638b 100644 --- a/ego/src/egor.rs +++ b/ego/src/egor.rs @@ -197,6 +197,7 @@ impl Egor { y_opt: result.state.get_full_best_cost().unwrap().to_owned(), x_hist: x_data, y_hist: y_data, + state: result.state, } } else { let x_data = to_discrete_space(&xtypes, &x_data.view()); @@ -214,6 +215,7 @@ impl Egor { y_opt: result.state.get_full_best_cost().unwrap().to_owned(), x_hist: x_data, y_hist: y_data, + state: result.state, } }; info!("Optim Result: min f(x)={} at x={}", res.y_opt, res.x_opt); @@ -280,6 +282,7 @@ mod tests { .max_iters(20) .regression_spec(RegressionSpec::ALL) .correlation_spec(CorrelationSpec::ALL) + .seed(1) }) .min_within(&array![[0.0, 25.0]]) .run() @@ -321,13 +324,7 @@ mod tests { let xlimits = array![[0.0, 25.0]]; let doe = Lhs::new(&xlimits).sample(10); let res = EgorBuilder::optimize(xsinx) - .configure(|config| { - config - .max_iters(15) - .doe(&doe) - .outdir(outdir) - .random_seed(42) - }) + .configure(|config| config.max_iters(15).doe(&doe).outdir(outdir).seed(42)) .min_within(&xlimits) .run() .expect("Minimize failure"); @@ -335,13 +332,7 @@ mod tests { assert_abs_diff_eq!(expected, res.x_opt, epsilon = 1e-1); let res = EgorBuilder::optimize(xsinx) - .configure(|config| { - config - .max_iters(5) - .outdir(outdir) - .hot_start(true) - .random_seed(42) - }) + .configure(|config| config.max_iters(5).outdir(outdir).hot_start(true).seed(42)) .min_within(&xlimits) .run() .expect("Egor should minimize xsinx"); @@ -375,7 +366,7 @@ mod tests { .regression_spec(RegressionSpec::ALL) .correlation_spec(CorrelationSpec::ALL) .target(1e-2) - .random_seed(42) + .seed(42) }) .min_within(&xlimits) .run() @@ -395,7 +386,7 @@ mod tests { .with_rng(Xoshiro256Plus::seed_from_u64(42)) .sample(10); let res = EgorBuilder::optimize(rosenb) - .configure(|config| config.doe(&doe).max_iters(20).random_seed(42)) + .configure(|config| config.doe(&doe).max_iters(20).seed(42)) .min_within(&xlimits) .run() .expect("Minimize failure"); @@ -445,7 +436,7 @@ mod tests { .doe(&doe) .max_iters(20) .cstr_tol(array![2e-6, 1e-6]) - .random_seed(42) + .seed(42) }) .min_within(&xlimits) .run() @@ -474,7 +465,7 @@ mod tests { .doe(&doe) .target(-5.5030) .max_iters(30) - .random_seed(42) + .seed(42) }) .min_within(&xlimits) .run() @@ -508,7 +499,7 @@ mod tests { .max_iters(max_iters) .target(-15.1) .infill_strategy(InfillStrategy::EI) - .random_seed(42) + .seed(42) }) .min_within_mixint_space(&xtypes) .run() @@ -530,7 +521,7 @@ mod tests { .max_iters(max_iters) .target(-15.1) .infill_strategy(InfillStrategy::EI) - .random_seed(42) + .seed(42) }) .min_within_mixint_space(&xtypes) .run() @@ -550,7 +541,7 @@ mod tests { .regression_spec(egobox_moe::RegressionSpec::CONSTANT) .correlation_spec(egobox_moe::CorrelationSpec::SQUAREDEXPONENTIAL) .max_iters(max_iters) - .random_seed(42) + .seed(42) }) .min_within_mixint_space(&xtypes) .run() @@ -601,7 +592,7 @@ mod tests { .regression_spec(egobox_moe::RegressionSpec::CONSTANT) .correlation_spec(egobox_moe::CorrelationSpec::SQUAREDEXPONENTIAL) .max_iters(max_iters) - .random_seed(42) + .seed(42) }) .min_within_mixint_space(&xtypes) .run() @@ -632,7 +623,7 @@ mod tests { let xlimits = as_continuous_limits::(&xtypes); EgorBuilder::optimize(mixobj) - .configure(|config| config.outdir(outdir).max_iters(1).random_seed(42)) + .configure(|config| config.outdir(outdir).max_iters(1).seed(42)) .min_within_mixint_space(&xtypes) .run() .unwrap(); @@ -644,13 +635,7 @@ mod tests { // Check that with no iteration, obj function is never called // as the DOE does not need to be evaluated! EgorBuilder::optimize(|_x| panic!("Should not call objective function!")) - .configure(|config| { - config - .outdir(outdir) - .hot_start(true) - .max_iters(0) - .random_seed(42) - }) + .configure(|config| config.outdir(outdir).hot_start(true).max_iters(0).seed(42)) .min_within_mixint_space(&xtypes) .run() .unwrap(); diff --git a/ego/src/egor_config.rs b/ego/src/egor_config.rs index a13b96c8..420aa573 100644 --- a/ego/src/egor_config.rs +++ b/ego/src/egor_config.rs @@ -240,7 +240,7 @@ impl EgorConfig { /// Allow to specify a seed for random number generator to allow /// reproducible runs. - pub fn random_seed(mut self, seed: u64) -> Self { + pub fn seed(mut self, seed: u64) -> Self { self.seed = Some(seed); self } diff --git a/ego/src/egor_service.rs b/ego/src/egor_service.rs index e6379c30..902a6c95 100644 --- a/ego/src/egor_service.rs +++ b/ego/src/egor_service.rs @@ -20,7 +20,7 @@ //! conf.regression_spec(RegressionSpec::ALL) //! .correlation_spec(CorrelationSpec::ALL) //! .infill_strategy(InfillStrategy::EI) -//! .random_seed(42) +//! .seed(42) //! }) //! .min_within(&array![[0., 25.]]); //! @@ -156,7 +156,7 @@ mod tests { conf.regression_spec(RegressionSpec::ALL) .correlation_spec(CorrelationSpec::ALL) .infill_strategy(InfillStrategy::EI) - .random_seed(42) + .seed(42) }) .min_within(&array![[0., 25.]]); diff --git a/ego/src/egor_solver.rs b/ego/src/egor_solver.rs index 3a215508..df17a2af 100644 --- a/ego/src/egor_solver.rs +++ b/ego/src/egor_solver.rs @@ -109,7 +109,7 @@ use crate::egor_config::EgorConfig; use crate::egor_state::{find_best_result_index, EgorState, MAX_POINT_ADDITION_RETRY}; use crate::errors::{EgoError, Result}; -use crate::mixint::*; +use crate::{find_best_result_index_from, mixint::*}; use crate::optimizer::*; @@ -361,7 +361,7 @@ where let no_point_added_retries = MAX_POINT_ADDITION_RETRY; let mut initial_state = state - .data((x_data, y_data)) + .data((x_data, y_data.clone())) .clusterings(clusterings) .theta_inits(theta_inits) .sampling(sampling); @@ -375,6 +375,10 @@ where .clone() .unwrap_or(Array1::from_elem(self.config.n_cstr, DEFAULT_CSTR_TOL)); initial_state.target_cost = self.config.target; + + let best_index = find_best_result_index(&y_data, &initial_state.cstr_tol); + initial_state.best_index = Some(best_index); + initial_state.last_best_iter = 0; debug!("INITIAL STATE = {:?}", initial_state); Ok((initial_state, None)) } @@ -437,7 +441,7 @@ where let (x_dat, y_dat, infill_value) = self.next_points( init, - new_state.get_iter(), + state.get_iter(), recluster, &mut clusterings, &mut theta_inits, @@ -532,16 +536,26 @@ where info!("Save doe shape {:?} in {:?}", doe.shape(), filepath); write_npy(filepath, &doe).expect("Write current doe"); } - let best_index = find_best_result_index(&y_data, &new_state.cstr_tol); + + let best_index = find_best_result_index_from( + state.best_index.unwrap(), + y_data.nrows() - add_count as usize, + &y_data, + &new_state.cstr_tol, + ); + // let best = find_best_result_index(&y_data, &new_state.cstr_tol); + // assert!(best_index == best); + new_state.best_index = Some(best_index); info!( "********* End iteration {}/{} in {:.3}s: Best fun(x)={} at x={}", - new_state.get_iter() + 1, - new_state.get_max_iters(), + state.get_iter() + 1, + state.get_max_iters(), now.elapsed().as_secs_f64(), y_data.row(best_index), x_data.row(best_index) ); new_state = new_state.data((x_data, y_data.clone())); + Ok((new_state, None)) } @@ -549,6 +563,8 @@ where debug!(">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> end iteration"); debug!("Current Cost {:?}", state.get_cost()); debug!("Best cost {:?}", state.get_best_cost()); + debug!("Best index {:?}", state.best_index); + debug!("Data {:?}", state.data.as_ref().unwrap()); TerminationStatus::NotTerminated } diff --git a/ego/src/egor_state.rs b/ego/src/egor_state.rs index 77b35438..6daf2993 100644 --- a/ego/src/egor_state.rs +++ b/ego/src/egor_state.rs @@ -4,7 +4,7 @@ use argmin::core::{ArgminFloat, Problem, State, TerminationReason, TerminationSt use egobox_doe::Lhs; use egobox_moe::Clustering; use linfa::Float; -use ndarray::{array, s, Array1, Array2, ArrayBase, Axis, Data, Ix2, Zip}; +use ndarray::{array, s, Array1, Array2, ArrayBase, Axis, Data, Ix1, Ix2, Zip}; use ndarray_stats::QuantileExt; use rand_xoshiro::Xoshiro256Plus; use serde::{Deserialize, Serialize}; @@ -15,6 +15,63 @@ use std::{collections::HashMap, iter::zip}; /// to train surrogate models modeling objective and constraints functions. pub(crate) const MAX_POINT_ADDITION_RETRY: i32 = 3; +fn cstr_sum(y: &ArrayBase, Ix1>, cstr_tol: &Array1) -> F { + y.slice(s![1..]) + .iter() + .enumerate() + .filter(|(i, &c)| c > cstr_tol[*i]) + .fold(F::zero(), |acc, (i, &c)| acc + (c - cstr_tol[i]).abs()) +} + +pub fn cstr_min( + (_, y1): (usize, &ArrayBase, Ix1>), + (_, y2): (usize, &ArrayBase, Ix1>), + cstr_tol: &Array1, +) -> std::cmp::Ordering { + if y1.len() > 1 { + let sum_c1 = cstr_sum(y1, cstr_tol); + let sum_c2 = cstr_sum(y2, cstr_tol); + if sum_c1 > F::zero() && sum_c2 > F::zero() { + sum_c1.partial_cmp(&sum_c2).unwrap() + } else if sum_c1 == F::zero() && sum_c2 == F::zero() { + y1[0].partial_cmp(&y2[0]).unwrap() + } else if sum_c1 == F::zero() { + std::cmp::Ordering::Less + } else { + std::cmp::Ordering::Greater + } + } else { + // unconstrained optimization + y1[0].partial_cmp(&y2[0]).unwrap() + } +} + +/// This method find the best result in ydata wrt cstr_min partial order +/// given the current best index of the current best and the offset index +/// the starting index of the new data to compare +/// +/// This method supposes indexes consistency: current_index < offset_index < ydata.nrows() +pub fn find_best_result_index_from( + current_index: usize, /* current best index */ + offset_index: usize, + ydata: &ArrayBase, Ix2>, /* the whole data so far */ + cstr_tol: &Array1, +) -> usize { + let new_ydata = ydata.slice(s![offset_index.., ..]); + + let best = ydata.row(current_index); + let min = new_ydata + .outer_iter() + .enumerate() + .fold((usize::MAX, best), |a, b| { + std::cmp::min_by(a, b, |(i, u), (j, v)| cstr_min((*i, u), (*j, v), cstr_tol)) + }); + match min { + (usize::MAX, _) => current_index, + (index, _) => offset_index + index, + } +} + /// Find best (eg minimal) cost value (y_data\[0\]) with valid constraints (y_data\[1..\] < cstr_tol). /// y_data containing ns samples [objective, cstr_1, ... cstr_nc] is given as a matrix (ns, nc + 1) pub fn find_best_result_index( @@ -85,6 +142,38 @@ mod tests { use super::*; use approx::assert_abs_diff_eq; + #[test] + fn test_cstr_min() { + // cstr respected, y1 > y2 + let (y1, y2) = (array![1.0, -0.15], array![-1.0, -0.01]); + let cstr_tol = Array1::from_elem(1, 0.1); + assert_eq!( + std::cmp::Ordering::Greater, + cstr_min((0, &y1), (1, &y2), &cstr_tol) + ); + // cstr respected, y1 > y2 + let (y1, y2) = (array![-1.0, -0.01], array![-2.0, -0.2]); + let cstr_tol = Array1::from_elem(1, 0.1); + assert_eq!( + std::cmp::Ordering::Greater, + cstr_min((0, &y1), (1, &y2), &cstr_tol) + ); + // cstr out of tolerance, y1 < y2 + let (y1, y2) = (array![-2.0, 1.], array![-1.0, 2.]); + let cstr_tol = Array1::from_elem(1, 0.1); + assert_eq!( + std::cmp::Ordering::Less, + cstr_min((0, &y1), (1, &y2), &cstr_tol) + ); + // cstr1 out of tolerance, y1 > y2 + let (y1, y2) = (array![-2.0, 1.], array![-1.0, 0.01]); + let cstr_tol = Array1::from_elem(1, 0.1); + assert_eq!( + std::cmp::Ordering::Greater, + cstr_min((0, &y1), (1, &y2), &cstr_tol) + ); + } + #[test] fn test_find_best_obj() { // respect constraint (0, 1, 2) and minimize obj (1) @@ -142,6 +231,36 @@ mod tests { let index = find_best_result_index(&y_data, &cstr_tol); assert_eq!(17, index); } + + #[test] + fn test_find_best_result2() { + let y_data = array![ + [-1.05044744051, -3.854157649246, 0.03068950747], + [-2.74965213562, -1.955703115787, -0.70939921583], + [-2.35364705246, 0.322760821911, -31.26001920874], + [-4.30684045535, -0.188609601161, 0.12375208631], + [-2.66585377971, -1.665992782883, -3.31489212502], + [-5.76598597442, 1.767753631322, -0.23219495778], + [-3.84677718652, -0.164470342807, -0.43935857142], + [-4.23672675117, -2.343687786724, -0.86266607911], + [-1.23999705899, -1.653209288978, -12.42363834689], + [-5.81590801801, -11.725502513342, 2.72175031293], + [-5.57379997815, -0.075893786744, 0.12260068082], + [-5.26821022904, -0.093334332384, -0.29931405911], + [-5.50558228637, -0.008847697249, 0.00015874647], + [-5.50802373110, -2.479726473358e-5, 2.46930218281e-5], + [-5.50802210236, -2.586721399788e-5, 2.28386911871e-5], + [-5.50801726964, 2.607167473023e-7, 5.50684865174e-6], + [-5.50801509642, 1.951629235996e-7, 2.48275059533e-6], + [-5.50801399313, -6.707576982734e-8, 1.03991762046e-6] + ]; + let cstr_tol = Array1::from_vec(vec![2e-6; 2]); + let index = find_best_result_index_from(11, 12, &y_data, &cstr_tol); + assert_eq!(17, index); + let cstr_tol = Array1::from_vec(vec![2e-6; 2]); + let index = find_best_result_index(&y_data, &cstr_tol); + assert_eq!(17, index); + } } /// Maintains the state from iteration to iteration of the [crate::EgorSolver]. @@ -207,6 +326,8 @@ pub struct EgorState { pub(crate) theta_inits: Option>>>, /// Historic data (params, objective and constraints) pub(crate) data: Option<(Array2, Array2)>, + /// index of best result in data + pub(crate) best_index: Option, /// Sampling method used to generate space filling samples pub(crate) sampling: Option>, } @@ -274,7 +395,7 @@ where /// # use egobox_ego::EgorState; /// # use argmin::core::{State, ArgminFloat}; /// # let state: EgorState = EgorState::new(); - /// # assert_eq!(state.max_iters, std::u64::MAX); + /// # assert_eq!(state.max_iters, u64::MAX); /// let state = state.max_iters(1000); /// # assert_eq!(state.max_iters, 1000); /// ``` @@ -439,7 +560,7 @@ where /// # assert_eq!(state.target_cost, f64::NEG_INFINITY); /// # assert_eq!(state.iter, 0); /// # assert_eq!(state.last_best_iter, 0); - /// # assert_eq!(state.max_iters, std::u64::MAX); + /// # assert_eq!(state.max_iters, u64::MAX); /// # assert_eq!(state.counts, HashMap::new()); /// # assert_eq!(state.time.unwrap(), instant::Duration::new(0, 0)); /// # assert_eq!(state.termination_status, TerminationStatus::NotTerminated); @@ -459,7 +580,7 @@ where iter: 0, last_best_iter: 0, - max_iters: std::u64::MAX, + max_iters: u64::MAX, counts: HashMap::new(), time: Some(instant::Duration::new(0, 0)), termination_status: TerminationStatus::NotTerminated, @@ -474,6 +595,7 @@ where clusterings: None, data: None, + best_index: None, sampling: None, theta_inits: None, } @@ -506,35 +628,24 @@ where /// assert!(state.is_best()); /// ``` fn update(&mut self) { - // TODO: better implementation should only track - // current and best index in data and compare just them - // without finding best in data each time - let data = self.data.as_ref(); - match data { - None => { - // should not occur data should be some - println!("Warning: update should occur after data initialization"); - } - Some((x_data, y_data)) => { - let best_index = find_best_result_index(y_data, &self.cstr_tol); - let best_iter = best_index.saturating_sub(self.doe_size) as u64 + 1; - if best_iter > self.last_best_iter { - let param = x_data.row(best_index).to_owned(); - std::mem::swap(&mut self.prev_best_param, &mut self.best_param); - self.best_param = Some(param); - - let cost = y_data.row(best_index).to_owned(); - std::mem::swap(&mut self.prev_best_cost, &mut self.best_cost); - self.best_cost = Some(cost); - if best_index > self.doe_size { - self.last_best_iter = best_iter; - } else { - // best point in doe => self.last_best_iter remains 0 - } - } - if self.best_cost.is_none() {} + if let Some((x_data, y_data)) = self.data.as_ref() { + let best_index = self + .best_index + .unwrap_or_else(|| find_best_result_index(y_data, &self.cstr_tol)); + + let param = x_data.row(best_index).to_owned(); + std::mem::swap(&mut self.prev_best_param, &mut self.best_param); + self.best_param = Some(param); + + let cost = y_data.row(best_index).to_owned(); + std::mem::swap(&mut self.prev_best_cost, &mut self.best_cost); + self.best_cost = Some(cost); + if best_index > self.doe_size { + self.last_best_iter = (best_index + 1 - self.doe_size) as u64; + } else { + // best point in doe => self.last_best_iter remains 0 } - }; + } } /// Returns a reference to the current parameter vector diff --git a/ego/src/lib.rs b/ego/src/lib.rs index 5916fcdd..ecd4b6eb 100644 --- a/ego/src/lib.rs +++ b/ego/src/lib.rs @@ -75,7 +75,7 @@ //! config.doe(&doe) // we pass the initial doe //! .max_iters(max_iters) //! .infill_strategy(InfillStrategy::EI) -//! .random_seed(42)) +//! .seed(42)) //! .min_within_mixint_space(&xtypes) // We build a mixed-integer optimizer //! .run() //! .expect("Egor minimization"); diff --git a/ego/src/types.rs b/ego/src/types.rs index 09c77fd8..2d44f1af 100644 --- a/ego/src/types.rs +++ b/ego/src/types.rs @@ -1,4 +1,4 @@ -use crate::errors::Result; +use crate::{errors::Result, EgorState}; use argmin::core::CostFunction; use egobox_moe::{Clustering, MixtureGpSurrogate, ThetaTuning}; use linfa::Float; @@ -19,6 +19,8 @@ pub struct OptimResult { pub x_hist: Array2, /// History of successive y values (e.g f(x_hist)) pub y_hist: Array2, + /// EgorSolver final state + pub state: EgorState, } /// Infill criterion used to select next promising point diff --git a/gp/src/utils.rs b/gp/src/utils.rs index edd075a7..4ed0d400 100644 --- a/gp/src/utils.rs +++ b/gp/src/utils.rs @@ -58,20 +58,17 @@ pub struct DistanceMatrix { pub d: Array2, pub d_indices: Array2, pub n_obs: usize, - pub n_features: usize, } impl DistanceMatrix { pub fn new(x: &ArrayBase, Ix2>) -> DistanceMatrix { let (d, d_indices) = Self::_cross_distances(x); let n_obs = x.nrows(); - let n_features = x.ncols(); DistanceMatrix { d: d.to_owned(), d_indices: d_indices.to_owned(), n_obs, - n_features, } } diff --git a/src/egor.rs b/src/egor.rs index b6a5f95d..683eae9e 100644 --- a/src/egor.rs +++ b/src/egor.rs @@ -468,7 +468,7 @@ impl Egor { config = config.outdir(outdir); }; if let Some(seed) = self.seed { - config = config.random_seed(seed); + config = config.seed(seed); }; config }