Use fonction constraint value to determine best result

relf · relf · commit 420be04c5ed3 · 2025-02-12T11:38:20.000+01:00
diff --git a/crates/ego/src/solver/egor_impl.rs b/crates/ego/src/solver/egor_impl.rs
@@ -83,7 +83,7 @@ impl<SB: SurrogateBuilder + DeserializeOwned, C: CstrFn> EgorSolver<SB, C> {
             &cstr_tol,
             &sampling,
             None,
-            find_best_result_index(y_data, &cstr_tol),
+            find_best_result_index(y_data, &c_data, &cstr_tol),
             &fcstrs,
         );
         x_dat
@@ -375,6 +375,7 @@ where
             state.best_index.unwrap(),
             y_data.nrows() - add_count as usize,
             &y_data,
+            &c_data,
             &new_state.cstr_tol,
         );
         new_state.prev_best_index = state.best_index;
diff --git a/crates/ego/src/solver/egor_solver.rs b/crates/ego/src/solver/egor_solver.rs
@@ -236,7 +236,7 @@ where
         let c_data = self.eval_fcstrs(problem, &x_data);
 
         let mut initial_state = state
-            .data((x_data, y_data.clone(), c_data))
+            .data((x_data, y_data.clone(), c_data.clone()))
             .clusterings(clusterings)
             .theta_inits(theta_inits)
             .sampling(sampling);
@@ -251,7 +251,7 @@ where
             .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);
+        let best_index = find_best_result_index(&y_data, &c_data, &initial_state.cstr_tol);
         initial_state.best_index = Some(best_index);
         initial_state.prev_best_index = Some(best_index);
         initial_state.last_best_iter = 0;
diff --git a/crates/ego/src/solver/egor_state.rs b/crates/ego/src/solver/egor_state.rs
@@ -408,11 +408,11 @@ where
     /// assert!(state.is_best());
     /// ```
     fn update(&mut self) {
-        if let Some((x_data, y_data, _c_data)) = self.data.as_ref() {
+        if let Some((x_data, y_data, c_data)) = self.data.as_ref() {
             let best_index = self
                 .best_index
                 // TODO: use cdata in find_best_result_index
-                .unwrap_or_else(|| find_best_result_index(y_data, &self.cstr_tol));
+                .unwrap_or_else(|| find_best_result_index(y_data, c_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);
diff --git a/crates/ego/src/solver/trego.rs b/crates/ego/src/solver/trego.rs
@@ -99,6 +99,7 @@ impl<SB: SurrogateBuilder + DeserializeOwned, C: CstrFn> EgorSolver<SB, C> {
             best_index,
             y_data.nrows() - 1,
             &y_data,
+            &c_data,
             &new_state.cstr_tol,
         );
         new_state = new_state.data((x_data, y_data, c_data));
diff --git a/crates/ego/src/utils/find_result.rs b/crates/ego/src/utils/find_result.rs
@@ -1,4 +1,4 @@
-use ndarray::{Array1, Array2, ArrayBase, Axis, Data, Ix1, Ix2, Zip};
+use ndarray::{concatenate, Array1, Array2, ArrayBase, Axis, Data, Ix1, Ix2, Zip};
 use ndarray_stats::QuantileExt;
 
 use crate::utils::sort_axis::*;
@@ -46,39 +46,54 @@ pub fn find_best_result_index_from<F: Float>(
     current_index: usize, /* current best index */
     offset_index: usize,
     ydata: &ArrayBase<impl Data<Elem = F>, Ix2>, /* the whole data so far */
+    cdata: &ArrayBase<impl Data<Elem = F>, Ix2>, /* the whole function cstrs data so far */
     cstr_tol: &Array1<F>,
 ) -> usize {
-    let new_ydata = ydata.slice(s![offset_index.., ..]);
+    // let new_ydata = ydata.slice(s![offset_index.., ..]);
+
+    let alldata = concatenate![Axis(1), ydata.to_owned(), cdata.to_owned()];
+    let alltols = concatenate![Axis(0), cstr_tol.to_owned(), Array1::zeros(cdata.ncols())];
+
+    let new_ydata = alldata.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))
+            std::cmp::min_by(a, b, |(i, u), (j, v)| cstr_min((*i, u), (*j, v), &alltols))
         });
     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)  
+/// Find best (eg minimal) cost value (y_data\[0\]) with valid constraints, meaning
+/// * y_data\[1..\] < cstr_tol
+/// * c_data[..] < 0
+///
+/// y_data containing ns samples [objective, cstr_1, ... cstr_nc] is given as a matrix (ns, nc + 1)
+/// c_data containing [fcstr_1, ... fcstr1_nfc] where fcstr_i is the value of function constraints at x_i
 pub fn find_best_result_index<F: Float>(
     y_data: &ArrayBase<impl Data<Elem = F>, Ix2>,
+    c_data: &ArrayBase<impl Data<Elem = F>, Ix2>,
     cstr_tol: &Array1<F>,
 ) -> usize {
-    if y_data.ncols() > 1 {
+    if y_data.ncols() > 1 || c_data.ncols() > 0 {
+        // Merge metamodelised constraints and function constraints
+        let alldata = concatenate![Axis(1), y_data.to_owned(), c_data.to_owned()];
+        let alltols = concatenate![Axis(0), cstr_tol.to_owned(), Array1::zeros(c_data.ncols())];
+
         // Compute sum of violated constraints
-        let cstrs = y_data.slice(s![.., 1..]);
+        let cstrs = &alldata.slice(s![.., 1..]);
         let mut c_obj = Array2::zeros((y_data.nrows(), 2));
 
         Zip::from(c_obj.rows_mut())
             .and(cstrs.rows())
-            .and(y_data.slice(s![.., 0]))
+            .and(alldata.slice(s![.., 0]))
             .for_each(|mut c_obj_row, c_row, obj| {
-                let c_sum = zip(c_row, cstr_tol)
+                let c_sum = zip(c_row, &alltols)
                     .filter(|(c, ctol)| *c > ctol)
                     .fold(F::zero(), |acc, (c, ctol)| acc + (*c - *ctol).abs());
                 c_obj_row.assign(&array![c_sum, *obj]);
@@ -106,13 +121,13 @@ pub fn find_best_result_index<F: Float>(
             let mut index = 0;
 
             // sort regardoing minimal objective
-            let perm = y_data.sort_axis_by(Axis(0), |i, j| y_data[[i, 0]] < y_data[[j, 0]]);
-            let y_sort = y_data.to_owned().permute_axis(Axis(0), &perm);
+            let perm = alldata.sort_axis_by(Axis(0), |i, j| alldata[[i, 0]] < alldata[[j, 0]]);
+            let y_sort = alldata.to_owned().permute_axis(Axis(0), &perm);
 
             // Take the first one which do not violate constraints
             for (i, row) in y_sort.axis_iter(Axis(0)).enumerate() {
                 let success =
-                    zip(row.slice(s![1..]), cstr_tol).fold(true, |acc, (c, tol)| acc && c < tol);
+                    zip(row.slice(s![1..]), &alltols).fold(true, |acc, (c, tol)| acc && c < tol);
 
                 if success {
                     index = i;
@@ -170,28 +185,29 @@ mod tests {
     fn test_find_best_obj() {
         // respect constraint (0, 1, 2) and minimize obj (1)
         let ydata = array![[1.0, -0.15], [-1.0, -0.01], [2.0, -0.2], [-3.0, 2.0]];
+        let cdata = array![[], [], [], []];
         let cstr_tol = Array1::from_elem(4, 0.1);
-        assert_abs_diff_eq!(1, find_best_result_index(&ydata, &cstr_tol));
+        assert_abs_diff_eq!(1, find_best_result_index(&ydata, &cdata, &cstr_tol));
 
         // respect constraint (0, 1, 2) and minimize obj (2)
         let ydata = array![[1.0, -0.15], [-1.0, -0.01], [-2.0, -0.2], [-3.0, 2.0]];
         let cstr_tol = Array1::from_elem(4, 0.1);
-        assert_abs_diff_eq!(2, find_best_result_index(&ydata, &cstr_tol));
+        assert_abs_diff_eq!(2, find_best_result_index(&ydata, &cdata, &cstr_tol));
 
         // all out of tolerance => minimize constraint overshoot sum (0)
         let ydata = array![[1.0, 0.15], [-1.0, 0.3], [2.0, 0.2], [-3.0, 2.0]];
         let cstr_tol = Array1::from_elem(4, 0.1);
-        assert_abs_diff_eq!(0, find_best_result_index(&ydata, &cstr_tol));
+        assert_abs_diff_eq!(0, find_best_result_index(&ydata, &cdata, &cstr_tol));
 
         // all in tolerance => min obj
         let ydata = array![[1.0, 0.15], [-1.0, 0.3], [2.0, 0.2], [-3.0, 2.0]];
         let cstr_tol = Array1::from_elem(4, 3.0);
-        assert_abs_diff_eq!(3, find_best_result_index(&ydata, &cstr_tol));
+        assert_abs_diff_eq!(3, find_best_result_index(&ydata, &cdata, &cstr_tol));
 
         // unconstrained => min obj
         let ydata = array![[1.0], [-1.0], [2.0], [-3.0]];
         let cstr_tol = Array1::from_elem(4, 0.1);
-        assert_abs_diff_eq!(3, find_best_result_index(&ydata, &cstr_tol));
+        assert_abs_diff_eq!(3, find_best_result_index(&ydata, &cdata, &cstr_tol));
     }
 
     #[test]
@@ -216,11 +232,12 @@ mod tests {
             [-5.50801509642, 1.951629235996e-7, 2.48275059533e-6],
             [-5.50801399313, -6.707576982734e-8, 1.03991762046e-6]
         ];
+        let c_data = Array2::zeros((y_data.nrows(), 0));
         let cstr_tol = Array1::from_vec(vec![1e-6; 2]); // this is the default
-        let index = find_best_result_index(&y_data, &cstr_tol);
+        let index = find_best_result_index(&y_data, &c_data, &cstr_tol);
         assert_eq!(11, index);
         let cstr_tol = Array1::from_vec(vec![2e-6; 2]);
-        let index = find_best_result_index(&y_data, &cstr_tol);
+        let index = find_best_result_index(&y_data, &c_data, &cstr_tol);
         assert_eq!(17, index);
     }
 
@@ -246,11 +263,12 @@ mod tests {
             [-5.50801509642, 1.951629235996e-7, 2.48275059533e-6],
             [-5.50801399313, -6.707576982734e-8, 1.03991762046e-6]
         ];
+        let c_data = Array2::zeros((y_data.nrows(), 0));
         let cstr_tol = Array1::from_vec(vec![2e-6; 2]);
-        let index = find_best_result_index_from(11, 12, &y_data, &cstr_tol);
+        let index = find_best_result_index_from(11, 12, &y_data, &c_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);
+        let index = find_best_result_index(&y_data, &c_data, &cstr_tol);
         assert_eq!(17, index);
     }
 }
diff --git a/crates/moe/src/algorithm.rs b/crates/moe/src/algorithm.rs
@@ -103,7 +103,7 @@ impl GpMixtureValidParams<f64> {
 
         let training = if recomb == Recombination::Smooth(None) && self.n_clusters() > 1 {
             // Extract 5% of data for validation to find best heaviside factor
-            // TODO: Use cross-validation ? Performances
+            // TODO: Better use cross-validation... but performances impact?
             let (_, training_data) = extract_part(&data, 5);
             training_data
         } else {
@@ -176,7 +176,7 @@ impl GpMixtureValidParams<f64> {
 
         if recomb == Recombination::Smooth(None) && self.n_clusters() > 1 {
             // Extract 5% of data for validation to find best heaviside factor
-            // TODO: Use cross-validation ? Performances
+            // TODO: Better use cross-validation... but performances impact?
             let (test, _) = extract_part(&data, 5);
             let xtest = test.slice(s![.., ..nx]).to_owned();
             let ytest = test.slice(s![.., nx..]).to_owned().remove_axis(Axis(1));
diff --git a/python/src/egor.rs b/python/src/egor.rs
@@ -376,6 +376,7 @@ impl Egor {
 
     /// This function gives the best evaluation index given the outputs
     /// of the function (objective wrt constraints) under minimization.
+    /// Caveat: This function does not take into account function constraints values
     ///
     /// # Parameters
     ///     y_doe (array[ns, 1 + n_cstr]): ns values of objective and constraints
@@ -386,11 +387,14 @@ impl Egor {
     #[pyo3(signature = (y_doe))]
     fn get_result_index(&self, y_doe: PyReadonlyArray2<f64>) -> usize {
         let y_doe = y_doe.as_array();
-        find_best_result_index(&y_doe, &self.cstr_tol())
+        // TODO: Make c_doe an optional argument ?
+        let c_doe = Array2::zeros((y_doe.ncols(), 0));
+        find_best_result_index(&y_doe, &c_doe, &self.cstr_tol())
     }
 
     /// This function gives the best result given inputs and outputs
     /// of the function (objective wrt constraints) under minimization.
+    /// Caveat: This function does not take into account function constraints values
     ///
     /// # Parameters
     ///     x_doe (array[ns, nx]): ns samples where function has been evaluated
@@ -410,7 +414,9 @@ impl Egor {
     ) -> OptimResult {
         let x_doe = x_doe.as_array();
         let y_doe = y_doe.as_array();
-        let idx = find_best_result_index(&y_doe, &self.cstr_tol());
+        // TODO: Make c_doe an optional argument ?
+        let c_doe = Array2::zeros((y_doe.ncols(), 0));
+        let idx = find_best_result_index(&y_doe, &c_doe, &self.cstr_tol());
         let x_opt = x_doe.row(idx).to_pyarray_bound(py).into();
         let y_opt = y_doe.row(idx).to_pyarray_bound(py).into();
         let x_doe = x_doe.to_pyarray_bound(py).into();
