Refactoring mostly for shebo

README.md

@@ -33,7 +33,7 @@ Surrogate-based 0-th Order Global Optimization for black-box problems.
 | `MaximizeEI`             | Maximize the expected improvement acquisition function for Gaussian processes. Use the dispersion-enhanced strategy from [(Müller; 2024)][Muller2024] for batch sampling.                                                                                                                                 |
 | `ParetoFront`            | Sample at the Pareto front of the multi-objective surrogate model to fill gaps in the surface [(Müller; 2017a)][Muller2017a].                                                                                                                                                                             |
 | `MinimizeMOSurrogate`    | Obtain pareto-optimal sample points for the multi-objective surrogate model [(Müller; 2017a)][Muller2017a].                                                                                                                                                                                               |
-| `GosacSample`            | Minimize a function with surrogate constraints to obtain a single new sample point [(Müller; 2017b)][Muller2017b].                                                                                                                                                                                        |
+| `GosacSample`            | Minimize a function with constraints to obtain a single new sample point [(Müller; 2017b)][Muller2017b].                                                                                                                                                                                                  |
 | `TransitionSearch`       | Weighted acquisition function that balances local and global search using a weighted score. Filters candidate points using evaluability surrogate. [(Müller & Day; 2019)][Muller2019].                                                                                                                    |
 | `MaximizeDistance`       | Maximizes the minimum distance to the set of current points. Used in `shebo()` and as a fallback in `EndPointsParetoFront` and `GosacSample` [(Müller & Day; 2019)][Muller2019].                                                                                                                          |
 

soogo/acquisition/gosac_sample.py

@@ -60,7 +60,7 @@ def optimize(
         surrogateModel: Surrogate,
         bounds,
         n: int = 1,
-        constraintTransform: callable = None,
+        constr_fun=None,
         **kwargs,
     ) -> np.ndarray:
         """Acquire 1 point.
@@ -69,10 +69,9 @@ def optimize(
         :param sequence bounds: List with the limits [x_min,x_max] of each
             direction x in the space.
         :param n: Unused.
-        :param constraintTransform: Function to transform the constraints.
-            If not provided, use the identity function. The optimizer, pymoo,
-            expects that negative and zero values are feasible, and positive
-            values are infeasible.
+        :param constr_fun: Constraint function to be applied to surrogate model
+            predictions. If none is provided, use the surrogate model as
+            the constraint function.
         :return: 1-by-dim matrix with the selected points.
         """
         dim = len(bounds)
@@ -81,20 +80,11 @@ def optimize(
         iindex = surrogateModel.iindex
         optimizer = self.optimizer if len(iindex) == 0 else self.mi_optimizer
 
-        if constraintTransform is None:
-
-            def constraintTransform(x):
-                return x
-
-        def transformedConstraint(x):
-            surrogateOutput = surrogateModel(x)
-            return constraintTransform(surrogateOutput)
-
         cheapProblem = PymooProblem(
             self.fun,
             bounds,
             iindex,
-            gfunc=transformedConstraint,
+            gfunc=surrogateModel if constr_fun is None else constr_fun,
             n_ieq_constr=gdim,
         )
         res = pymoo_minimize(

soogo/acquisition/utils.py

@@ -20,6 +20,7 @@
 
 import numpy as np
 from scipy.spatial.distance import cdist
+from scipy.spatial import KDTree
 import networkx as nx
 from networkx.algorithms.approximation import maximum_independent_set
 
@@ -31,7 +32,7 @@ def weighted_score(
     sx_min: float = 0.0,
     sx_max: float = 1.0,
     dx_max: float = 1.0,
-) -> float:
+):
     r"""Computes the weighted score from the predicted value of the surrogate
     model at a point and minimum distance from the point to the set of
     previously selected evaluation points.
@@ -226,10 +227,22 @@ class FarEnoughSampleFilter:
     distance threshold from already sampled points.
 
     :param X: Matrix of existing sample points (n x d).
-    :param tol: Minimum distance threshold. Points closer than this are
-        filtered out.
+    :param tol: Minimum distance threshold.
+
+    .. attribute:: tree
+
+        KDTree built from the existing sample points for efficient distance
+        queries.
+
+    .. attribute:: tol
+
+        Minimum distance threshold. Points closer than this are filtered out.
     """
 
+    def __init__(self, X, tol):
+        self.tree = KDTree(X)
+        self.tol = tol
+
     def is_far_enough(self, x):
         """Check if a point is far enough from existing samples.
 
@@ -239,11 +252,11 @@ def is_far_enough(self, x):
         dist, _ = self.tree.query(x.reshape(1, -1))
         return dist[0] >= self.tol
 
-    def __call__(self, Xc):
+    def indices(self, Xc):
         """Filter candidates based on minimum distance criterion.
 
         :param Xc: Matrix of candidate points (m x d).
-        :return: Filtered matrix containing only points that are far
+        :return: Filtered indices for points that are far
             enough from existing samples.
         """
         # Discard points that are too close to X
@@ -263,6 +276,17 @@ def __call__(self, Xc):
                 if dist[i, j] < self.tol
             ]
         )
-
         idx = maximum_independent_set(g)
-        return Xc0[list(idx)]
+
+        # Recover original indices
+        original_indices = np.where(mask0)[0]
+        return original_indices[list(idx)]
+
+    def __call__(self, Xc):
+        """Filter candidates based on minimum distance criterion.
+
+        :param Xc: Matrix of candidate points (m x d).
+        :return: Filtered matrix containing only points that are far
+            enough from existing samples.
+        """
+        return Xc[self.indices(Xc)]

soogo/acquisition/weighted_acquisition.py

@@ -146,6 +146,9 @@ def optimize(
         surrogateModel: Surrogate,
         bounds,
         n: int = 1,
+        *,
+        constr_fun=None,
+        perturbation_probability=None,
         **kwargs,
     ) -> np.ndarray:
         """Generate a number of candidates using the :attr:`sampler`. Then,
@@ -186,12 +189,15 @@ def optimize(
                 coord = [i for i in range(dim)]
 
             # Compute probability in case DDS is used
-            if self.maxeval > 1 and self.neval < self.maxeval:
-                prob = min(20 / dim, 1) * (
-                    1 - (log(self.neval + 1) / log(self.maxeval))
-                )
+            if perturbation_probability is None:
+                if self.maxeval > 1 and self.neval < self.maxeval:
+                    prob = min(20 / dim, 1) * (
+                        1 - (log(self.neval + 1) / log(self.maxeval))
+                    )
+                else:
+                    prob = 1.0
             else:
-                prob = 1.0
+                prob = perturbation_probability
 
             x = self.sampler.get_sample(
                 bounds,
@@ -203,6 +209,17 @@ def optimize(
         else:
             x = self.sampler.get_sample(bounds, iindex=iindex)
 
+        if constr_fun is not None:
+            # Filter out candidates that do not satisfy the constraints
+            constr_values = constr_fun(x)
+            if constr_values.ndim == 1:
+                feasible_idx = constr_values <= 0
+            else:
+                feasible_idx = np.all(constr_values <= 0, axis=1)
+            x = x[feasible_idx]
+            if x.shape[0] == 0:
+                return np.empty((0, dim))
+
         # Evaluate candidates
         fx = surrogateModel(x)
 

soogo/integrations/nomad.py

@@ -48,10 +48,12 @@ def __call__(self, x):
         point = np.array([x.get_coord(i) for i in range(x.size())])
         self._xHistory.append(point)
 
-        f = evaluate_and_log_point(self.func, point, self.out)
+        f = evaluate_and_log_point(self.func, point.reshape(1, -1), self.out)[
+            0
+        ]
         self._fHistory.append(f)
 
-        if not np.isnan(f):
+        if np.isfinite(f):
             # Set NOMAD objective function value
             x.setBBO(str(f).encode("UTF-8"))
             return 1

soogo/optimize/fsapso.py

@@ -122,7 +122,11 @@ def fsapso(
 
         # Evaluate initial points
         for x0 in xInit:
-            _ = evaluate_and_log_point(fun, x0, out)
+            y0 = evaluate_and_log_point(fun, x0.reshape(1, -1), out)[0]
+
+            if y0 < out.fx:
+                out.x[:] = x0
+                out.fx = y0
 
             if disp:
                 print("fEvals: %d" % out.nfev)
@@ -211,7 +215,11 @@ def fsapso(
         # Check xMin is at least tol away from existing points
         if np.min(cdist(xMin.reshape(1, -1), out.sample[: out.nfev])) > tol:
             # Evaluate minimum with true objective
-            fMin = evaluate_and_log_point(fun, xMin, out)
+            fMin = evaluate_and_log_point(fun, xMin.reshape(1, -1), out)[0]
+
+            if fMin < out.fx:
+                out.x[:] = xMin
+                out.fx = fMin
 
             if disp:
                 print("fEvals: %d" % out.nfev)
@@ -252,7 +260,13 @@ def fsapso(
                 )
                 > tol
             ):
-                fBestParticle = evaluate_and_log_point(fun, xBestParticle, out)
+                fBestParticle = evaluate_and_log_point(
+                    fun, xBestParticle.reshape(1, -1), out
+                )[0]
+
+                if fBestParticle < out.fx:
+                    out.x[:] = xBestParticle
+                    out.fx = fBestParticle
 
                 if disp:
                     print("fEvals: %d" % out.nfev)
@@ -296,8 +310,12 @@ def fsapso(
                 > tol
             ):
                 fMostUncertain = evaluate_and_log_point(
-                    fun, xMostUncertain, out
-                )
+                    fun, xMostUncertain.reshape(1, -1), out
+                )[0]
+
+                if fMostUncertain < out.fx:
+                    out.x[:] = xMostUncertain
+                    out.fx = fMostUncertain
 
                 if disp:
                     print("fEvals: %d" % out.nfev)