Store function constraint values

crates/ego/src/egor.rs

@@ -250,7 +250,7 @@ impl<O: GroupFunc, C: CstrFn, SB: SurrogateBuilder + DeserializeOwned> Egor<O, C
         };
 
         info!("{}", result);
-        let (x_data, y_data) = result.state().clone().take_data().unwrap();
+        let (x_data, y_data, _c_data) = result.state().clone().take_data().unwrap();
 
         let res = if !self.solver.config.discrete() {
             info!("Data: \n{}", concatenate![Axis(1), x_data, y_data]);

crates/ego/src/solver/egor_impl.rs

@@ -55,7 +55,6 @@ impl<SB: SurrogateBuilder + DeserializeOwned, C: CstrFn> EgorSolver<SB, C> {
         &self,
         x_data: &ArrayBase<impl Data<Elem = f64>, Ix2>,
         y_data: &ArrayBase<impl Data<Elem = f64>, Ix2>,
-        _cstr_funcs: &[&(dyn ObjFn<InfillObjData<f64>> + Sync)],
     ) -> Array2<f64> {
         let rng = self.rng.clone();
         let sampling = Lhs::new(&self.xlimits).with_rng(rng).kind(LhsKind::Maximin);
@@ -66,14 +65,19 @@ impl<SB: SurrogateBuilder + DeserializeOwned, C: CstrFn> EgorSolver<SB, C> {
             .cstr_tol
             .clone()
             .unwrap_or(Array1::from_elem(self.config.n_cstr, DEFAULT_CSTR_TOL));
-        let (x_dat, _, _, _) = self.next_points(
+
+        // TODO: c_data has to be passed as argument or computed using fcstrs
+        let c_data = Array2::zeros((x_data.nrows(), 0));
+
+        let (x_dat, _, _, _, _) = self.next_points(
             true,
             0,
             false, // done anyway
             &mut clusterings,
             &mut theta_tunings,
             x_data,
             y_data,
+            &c_data,
             &cstr_tol,
             &sampling,
             None,
@@ -249,7 +253,7 @@ where
             PotentialBug,
             "EgorSolver: No sampling!"
         ))?;
-        let (mut x_data, mut y_data) = new_state
+        let (mut x_data, mut y_data, mut c_data) = new_state
             .take_data()
             .ok_or_else(argmin_error_closure!(PotentialBug, "EgorSolver: No data!"))?;
 
@@ -270,14 +274,15 @@ where
 
             let problem = fobj.take_problem().unwrap();
             let fcstrs = problem.fn_constraints();
-            let (x_dat, y_dat, infill_value, infill_data) = self.next_points(
+            let (x_dat, y_dat, c_dat, infill_value, infill_data) = self.next_points(
                 init,
                 state.get_iter(),
                 recluster,
                 &mut clusterings,
                 &mut theta_inits,
                 &x_data,
                 &y_data,
+                &c_data,
                 &state.cstr_tol,
                 &sampling,
                 lhs_optim_seed,
@@ -287,12 +292,19 @@ where
             fobj.problem = Some(problem);
 
             debug!("Try adding {}", x_dat);
-            let added_indices = update_data(&mut x_data, &mut y_data, &x_dat, &y_dat);
+            let added_indices = update_data(
+                &mut x_data,
+                &mut y_data,
+                &mut c_data,
+                &x_dat,
+                &y_dat,
+                &c_dat,
+            );
 
             new_state = new_state
                 .clusterings(clusterings.clone())
                 .theta_inits(theta_inits.clone())
-                .data((x_data.clone(), y_data.clone()))
+                .data((x_data.clone(), y_data.clone(), c_data.clone()))
                 .infill_value(infill_value)
                 .sampling(sampling.clone())
                 .param(x_dat.row(0).to_owned())
@@ -365,7 +377,7 @@ where
         );
         new_state.prev_best_index = state.best_index;
         new_state.best_index = Some(best_index);
-        new_state = new_state.data((x_data.clone(), y_data.clone()));
+        new_state = new_state.data((x_data.clone(), y_data.clone(), c_data.clone()));
 
         Ok((new_state, infill_data, best_index))
     }
@@ -384,15 +396,23 @@ where
         theta_inits: &mut [Option<Array2<f64>>],
         x_data: &ArrayBase<impl Data<Elem = f64>, Ix2>,
         y_data: &ArrayBase<impl Data<Elem = f64>, Ix2>,
+        c_data: &ArrayBase<impl Data<Elem = f64>, Ix2>,
         cstr_tol: &Array1<f64>,
         sampling: &Lhs<f64, Xoshiro256Plus>,
         lhs_optim: Option<u64>,
         best_index: usize,
         cstr_funcs: &[impl CstrFn],
-    ) -> (Array2<f64>, Array2<f64>, f64, InfillObjData<f64>) {
+    ) -> (
+        Array2<f64>,
+        Array2<f64>,
+        Array2<f64>,
+        f64,
+        InfillObjData<f64>,
+    ) {
         debug!("Make surrogate with {}", x_data);
         let mut x_dat = Array2::zeros((0, x_data.ncols()));
         let mut y_dat = Array2::zeros((0, y_data.ncols()));
+        let mut c_dat = Array2::zeros((0, c_data.ncols()));
         let mut infill_val = f64::INFINITY;
         let mut infill_data = Default::default();
         for i in 0..self.config.q_points {
@@ -458,14 +478,42 @@ where
                 scale_wb2,
             };
 
+            let cstr_funcs = cstr_funcs
+                .iter()
+                .map(|cstr| {
+                    |x: &[f64],
+                     gradient: Option<&mut [f64]>,
+                     params: &mut InfillObjData<f64>|
+                     -> f64 {
+                        let x = if self.config.discrete() {
+                            let xary = Array2::from_shape_vec((1, x.len()), x.to_vec()).unwrap();
+                            // We have to cast x to folded space as EgorSolver
+                            // works internally in the continuous space while
+                            // the constraint function expects discrete variable in folded space
+                            to_discrete_space(&self.config.xtypes, &xary)
+                                .row(0)
+                                .into_owned();
+                            &xary.into_iter().collect::<Vec<_>>()
+                        } else {
+                            x
+                        };
+                        cstr(x, gradient, params)
+                    }
+                })
+                .collect::<Vec<_>>();
+            let cstr_funcs = cstr_funcs
+                .iter()
+                .map(|cstr| cstr as &(dyn ObjFn<InfillObjData<f64>> + Sync))
+                .collect::<Vec<_>>();
+
             match self.find_best_point(
                 sampling,
                 obj_model.as_ref(),
                 cstr_models,
                 cstr_tol,
                 lhs_optim,
                 &infill_data,
-                cstr_funcs,
+                &cstr_funcs,
             ) {
                 Ok((infill_obj, xk)) => {
                     match self.get_virtual_point(&xk, y_data, obj_model.as_ref(), cstr_models) {
@@ -475,7 +523,19 @@ where
                                 y_dat,
                                 Array2::from_shape_vec((1, 1 + self.config.n_cstr), yk).unwrap()
                             ];
+
+                            let ck = cstr_funcs
+                                .iter()
+                                .map(|cstr| cstr(&xk.to_vec(), None, &mut infill_data))
+                                .collect::<Vec<_>>();
+                            c_dat = concatenate![
+                                Axis(0),
+                                c_dat,
+                                Array2::from_shape_vec((1, cstr_funcs.len()), ck).unwrap()
+                            ];
+
                             x_dat = concatenate![Axis(0), x_dat, xk.insert_axis(Axis(0))];
+
                             // infill objective was minimized while infill criterion itself
                             // is expected to be maximized hence the negative sign here
                             infill_val = -infill_obj;
@@ -494,7 +554,7 @@ where
                 }
             }
         }
-        (x_dat, y_dat, infill_val, infill_data)
+        (x_dat, y_dat, c_dat, infill_val, infill_data)
     }
 
     pub(crate) fn compute_scaling(
@@ -550,7 +610,7 @@ where
         cstr_tol: &Array1<f64>,
         lhs_optim_seed: Option<u64>,
         infill_data: &InfillObjData<f64>,
-        cstr_funcs: &[impl CstrFn],
+        cstr_funcs: &[&(dyn ObjFn<InfillObjData<f64>> + Sync)],
     ) -> Result<(f64, Array1<f64>)> {
         let fmin = infill_data.fmin;
 
@@ -626,30 +686,8 @@ where
             })
             .collect();
         let mut cstr_refs: Vec<_> = cstrs.iter().map(|c| c.as_ref()).collect();
-        let cstr_funcs = cstr_funcs
-            .iter()
-            .map(|cstr| {
-                |x: &[f64], gradient: Option<&mut [f64]>, params: &mut InfillObjData<f64>| -> f64 {
-                    let x = if self.config.discrete() {
-                        let xary = Array2::from_shape_vec((1, x.len()), x.to_vec()).unwrap();
-                        // We have to cast x to folded space as EgorSolver
-                        // works internally in the continuous space while
-                        // the constraint function expects discrete variable in folded space
-                        to_discrete_space(&self.config.xtypes, &xary)
-                            .row(0)
-                            .into_owned();
-                        &xary.into_iter().collect::<Vec<_>>()
-                    } else {
-                        x
-                    };
-                    cstr(x, gradient, params)
-                }
-            })
-            .collect::<Vec<_>>();
-        let cstr_funcs = cstr_funcs
-            .iter()
-            .map(|cstr| cstr as &(dyn ObjFn<InfillObjData<f64>> + Sync));
         cstr_refs.extend(cstr_funcs);
+
         info!("Optimize infill criterion...");
         while !success && n_optim <= n_max_optim {
             let x_start = sampling.sample(self.config.n_start);
@@ -818,15 +856,52 @@ where
         pb: &mut Problem<O>,
         x: &Array2<f64>,
     ) -> Array2<f64> {
-        let params = if self.config.discrete() {
+        let x = if self.config.discrete() {
             // We have to cast x to folded space as EgorSolver
             // works internally in the continuous space while
             // the objective function expects discrete variable in folded space
             to_discrete_space(&self.config.xtypes, x)
         } else {
             x.to_owned()
         };
-        pb.problem("cost_count", |problem| problem.cost(&params))
+        pb.problem("cost_count", |problem| problem.cost(&x))
             .expect("Objective evaluation")
     }
+
+    pub fn eval_fcstrs<O: DomainConstraints<C>>(
+        &self,
+        pb: &mut Problem<O>,
+        x: &Array2<f64>,
+    ) -> Array2<f64> {
+        let problem = pb.take_problem().unwrap();
+        let fcstrs = problem.fn_constraints();
+
+        let mut unused = InfillObjData::default();
+
+        let mut res = Array2::zeros((x.nrows(), fcstrs.len()));
+        Zip::from(res.rows_mut()).for_each(|mut r| {
+            r.assign(
+                &fcstrs
+                    .iter()
+                    .map(|cstr| {
+                        let x = if self.config.discrete() {
+                            let xary = x.to_owned();
+                            // We have to cast x to folded space as EgorSolver
+                            // works internally in the continuous space while
+                            // the constraint function expects discrete variable in folded space
+                            to_discrete_space(&self.config.xtypes, &xary)
+                                .row(0)
+                                .into_owned();
+                            &xary.into_iter().collect::<Vec<_>>()
+                        } else {
+                            x.as_slice().unwrap()
+                        };
+                        cstr(x, None, &mut unused)
+                    })
+                    .collect::<Array1<_>>(),
+            )
+        });
+        pb.problem = Some(problem);
+        res
+    }
 }

crates/ego/src/solver/egor_service.rs

@@ -139,7 +139,7 @@ impl<SB: SurrogateBuilder + DeserializeOwned, C: CstrFn> EgorServiceApi<SB, C> {
         let xtypes = &self.solver.config.xtypes;
         let x_data = to_continuous_space(xtypes, x_data);
         // TODO: cstr_funcs not managed
-        let x = self.solver.suggest(&x_data, y_data, &[]);
+        let x = self.solver.suggest(&x_data, y_data);
         to_discrete_space(xtypes, &x).to_owned()
     }
 }

crates/ego/src/solver/egor_solver.rs

@@ -233,8 +233,10 @@ where
         let theta_inits = vec![None; self.config.n_cstr + 1];
         let no_point_added_retries = MAX_POINT_ADDITION_RETRY;
 
+        let c_data = self.eval_fcstrs(problem, &x_data);
+
         let mut initial_state = state
-            .data((x_data, y_data.clone()))
+            .data((x_data, y_data.clone(), c_data))
             .clusterings(clusterings)
             .theta_inits(theta_inits)
             .sampling(sampling);
@@ -274,7 +276,7 @@ where
         } else {
             self.ego_iteration(fobj, state)?
         };
-        let (x_data, y_data) = res.0.data.clone().unwrap();
+        let (x_data, y_data, _c_data) = res.0.data.clone().unwrap();
 
         if self.config.outdir.is_some() {
             let doe = concatenate![Axis(1), x_data, y_data];
@@ -340,7 +342,7 @@ where
         state: EgorState<f64>,
     ) -> std::result::Result<(EgorState<f64>, Option<KV>), argmin::core::Error> {
         let rho = |sigma| sigma * sigma;
-        let (_, y_data) = state.data.as_ref().unwrap(); // initialized in init
+        let (_, y_data, _) = state.data.as_ref().unwrap(); // initialized in init
         let best = state.best_index.unwrap(); // initialized in init
         let prev_best = state.prev_best_index.unwrap(); // initialized in init
 

crates/ego/src/solver/egor_state.rs

@@ -76,8 +76,8 @@ pub struct EgorState<F: Float> {
     /// ThetaTunings controlled by n_optmod configuration triggering
     /// GP surrogate models hyperparameters optimization or reusing previous ones
     pub theta_inits: Option<Vec<Option<Array2<F>>>>,
-    /// Historic data (params, objective and constraints)
-    pub data: Option<(Array2<F>, Array2<F>)>,
+    /// Historic data (params, objective and constraints values, function constraints)
+    pub data: Option<(Array2<F>, Array2<F>, Array2<F>)>,
     /// Previous index of best result in data
     pub prev_best_index: Option<usize>,
     /// index of best result in data
@@ -216,15 +216,17 @@ where
 
     /// Set the current data points as training points for the surrogate models
     /// These points are gradually selected by the EGO algorithm regarding an infill criterion.
-    /// Data is expressed as a couple (xdata, ydata) where xdata is a (p, nx matrix)
-    /// and ydata is a (p, 1 + nb of cstr) matrix and ydata_i = fcost(xdata_i) for i in [1, p].  
-    pub fn data(mut self, data: (Array2<F>, Array2<F>)) -> Self {
+    /// Data is expressed as a triple (xdata, ydata, cdata) where :
+    /// * xdata is a (p, nx matrix),
+    /// * ydata is a (p, 1 + nb of cstr) matrix and ydata_i = fcost(xdata_i) for i in [1, p],  
+    /// * cdata is a (p, nb of fcstr) matrix and cdata_i = fcstr_j(xdata_i) for i in [1, p], j in [1, nb of fcstr]
+    pub fn data(mut self, data: (Array2<F>, Array2<F>, Array2<F>)) -> Self {
         self.data = Some(data);
         self
     }
 
     /// Moves the current data out and replaces it internally with `None`.
-    pub fn take_data(&mut self) -> Option<(Array2<F>, Array2<F>)> {
+    pub fn take_data(&mut self) -> Option<(Array2<F>, Array2<F>, Array2<F>)> {
         self.data.take()
     }
 
@@ -390,7 +392,7 @@ where
     /// let mut state: EgorState<f64> = EgorState::new();
     ///
     /// // Simulating a new, better parameter vector
-    /// let mut state = state.data((array![[1.0f64], [2.0f64], [3.0]], array![[10.0], [5.0], [0.5]]));
+    /// let mut state = state.data((array![[1.0f64], [2.0f64], [3.0]], array![[10.0], [5.0], [0.5]], array![[], [], []]));
     /// state.iter = 2;
     /// state.prev_best_index = Some(0);
     /// state.best_index = Some(2);
@@ -406,9 +408,10 @@ where
     /// assert!(state.is_best());
     /// ```
     fn update(&mut self) {
-        if let Some((x_data, y_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));
 
             let param = x_data.row(best_index).to_owned();