From ddc59d10516078b85575d69e81c21d419d9f71da Mon Sep 17 00:00:00 2001 From: Weslley da Silva Pereira Date: Tue, 16 Dec 2025 22:40:21 -0700 Subject: [PATCH 1/2] Fix bug introduced in last commit --- soogo/acquisition/utils.py | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/soogo/acquisition/utils.py b/soogo/acquisition/utils.py index 579f6bc5..c3edec3a 100644 --- a/soogo/acquisition/utils.py +++ b/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. From 85bc1b73b4a96983c57093c091e61de0b042bd53 Mon Sep 17 00:00:00 2001 From: Weslley da Silva Pereira Date: Tue, 16 Dec 2025 15:03:41 -0700 Subject: [PATCH 2/2] Refactoring mostly for shebo - evaluate_and_log_point and NomadProblem propagate both infs and nans - evaluate_and_log_point does not update best output values. This change will make it possible to use evaluate_and_log_point with other optimization problems in the library - WeightedAcquisition was generalized to deal with constraints and perturbation probability in the API. This almost make TransitionSearch superseeded, intentionally. We might remove TransitionSearch in future commits. - GosacSample.optimize has been simplified Changes in shebo: - Now the surrogate training sets are synced in the begining of the algorithm. - The new code works directly with the original bounds - The acquisition function falls back into MaximizeDistance if all points generated by the original shebo strategy happen to be discarded. --- README.md | 2 +- soogo/acquisition/gosac_sample.py | 20 +- soogo/acquisition/utils.py | 19 +- soogo/acquisition/weighted_acquisition.py | 27 +- soogo/integrations/nomad.py | 6 +- soogo/optimize/fsapso.py | 28 +- soogo/optimize/shebo.py | 523 +++++++++++----------- soogo/optimize/utils.py | 50 +-- 8 files changed, 346 insertions(+), 329 deletions(-) diff --git a/README.md b/README.md index f2b7ccde..4b8e3a8a 100644 --- a/README.md +++ b/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]. | diff --git a/soogo/acquisition/gosac_sample.py b/soogo/acquisition/gosac_sample.py index 0f65b681..e0177cc1 100644 --- a/soogo/acquisition/gosac_sample.py +++ b/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( diff --git a/soogo/acquisition/utils.py b/soogo/acquisition/utils.py index c3edec3a..b9862c8d 100644 --- a/soogo/acquisition/utils.py +++ b/soogo/acquisition/utils.py @@ -252,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 @@ -276,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)] diff --git a/soogo/acquisition/weighted_acquisition.py b/soogo/acquisition/weighted_acquisition.py index 6717a965..cc1dac58 100644 --- a/soogo/acquisition/weighted_acquisition.py +++ b/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) diff --git a/soogo/integrations/nomad.py b/soogo/integrations/nomad.py index dd570099..2ba9673f 100644 --- a/soogo/integrations/nomad.py +++ b/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 diff --git a/soogo/optimize/fsapso.py b/soogo/optimize/fsapso.py index 4b1887bd..b3318a82 100644 --- a/soogo/optimize/fsapso.py +++ b/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) diff --git a/soogo/optimize/shebo.py b/soogo/optimize/shebo.py index 28b743a7..5f314815 100644 --- a/soogo/optimize/shebo.py +++ b/soogo/optimize/shebo.py @@ -15,26 +15,27 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . +import time import warnings from typing import Callable, Optional import numpy as np -from scipy.spatial.distance import cdist, pdist, squareform +from scipy.spatial.distance import cdist + +from soogo.acquisition.base import Acquisition +from soogo.model.rbf import MedianLpfFilter from ..acquisition import ( AlternatedAcquisition, GosacSample, MaximizeDistance, - TransitionSearch, -) -from ..model import ( - CubicRadialBasisFunction, - LinearRadialBasisFunction, - RbfModel, - Surrogate, + WeightedAcquisition, + MultipleAcquisition, ) +from ..sampling import NormalSampler, SamplingStrategy +from ..acquisition.utils import FarEnoughSampleFilter +from ..model import LinearRadialBasisFunction, RbfModel, Surrogate from .utils import OptimizeResult, evaluate_and_log_point -from ..sampling import Sampler from ..termination import IterateNTimes from ..integrations.nomad import NomadProblem @@ -44,14 +45,41 @@ PyNomad = None +class NormalAndUniformSampler(NormalSampler): + """Sampler that generates half normal and half uniform samples.""" + + def get_sample( + self, bounds, *, iindex: tuple[int, ...] = (), **kwargs + ) -> np.ndarray: + """Generate n samples, half from normal distribution and half from + uniform distribution. + + :param bounds: List with the limits [x_min,x_max] of each direction x + in the space. + :param iindex: Indices of the dimensions to be sampled. If None, all + dimensions are sampled. + :return: Matrix with a sample point per line. + """ + _n = self.n + + self.n = _n // 2 + x_normal = super().get_sample(bounds, iindex=iindex, **kwargs) + + self.n = _n - self.n + x_uniform = self.get_uniform_sample(bounds, iindex=iindex) + + self.n = _n + return np.vstack((x_normal, x_uniform)) + + def shebo( fun, bounds, maxeval: int, *, - objSurrogate: Optional[Surrogate] = None, + objSurrogate: Optional[RbfModel] = None, evalSurrogate: Optional[Surrogate] = None, - acquisitionFunc: Optional[AlternatedAcquisition] = None, + acquisitionFunc: Optional[Acquisition] = None, disp: bool = False, callback: Optional[Callable[[OptimizeResult], None]] = None, ) -> OptimizeResult: @@ -96,288 +124,237 @@ def shebo( warnings.warn( "PyNomad package is required but not installed. Install the PyNomad package and try again." ) - return - - # Initialize parameters - weightPattern = [1, 0.95, 0.85, 0.75, 0.5, 0.35, 0.25, 0.1, 0.0] - dim = len(bounds) - nStart = 4 * (dim + 1) - bounds = np.asarray(bounds) + return OptimizeResult() - # Define function wrapper to rescale variables - # This is needed as SHEBO internally rescales all variables to [0, 1] - def rescaledFunc(x): - rescaledX = x * (bounds[:, 1] - bounds[:, 0]) + bounds[:, 0] - return fun(rescaledX) - - # Initialize output - out = OptimizeResult( - x=np.empty((0, dim)), - fx=np.inf, - nfev=0, - sample=np.zeros((maxeval, dim)), - fsample=np.zeros(maxeval), - ) + dim = len(bounds) # Dimension of the problem + assert dim > 0 + rtol = 1e-3 # Relative tolerance for distance-based operations + return_surrogate = (objSurrogate is not None) or ( + evalSurrogate is not None + ) # Whether to return the surrogate models - # Create Nomad wrapper - nomadFunction = NomadProblem(rescaledFunc, out) - - # Initialize or use provided surrogates + # Initialize optional variables if objSurrogate is None: - objSurrogate = RbfModel(CubicRadialBasisFunction()) + objSurrogate = RbfModel(filter=MedianLpfFilter()) if evalSurrogate is None: evalSurrogate = RbfModel(LinearRadialBasisFunction()) - - # Reserve space for the surrogates - objSurrogate.reserve(objSurrogate.ntrain + maxeval, dim) - evalSurrogate.reserve(evalSurrogate.ntrain + maxeval, dim) - - # Check if surrogates already trained - objTrained = objSurrogate.ntrain > 0 - evalTrained = evalSurrogate.ntrain > 0 - - # Default acquisition function if acquisitionFunc is None: - acquisitionFunc = AlternatedAcquisition( + acquisitionFunc = MultipleAcquisition( [ - GosacSample( - objSurrogate, rtol=0.001, termination=IterateNTimes(1) + AlternatedAcquisition( + [ + GosacSample( + objSurrogate, + rtol=rtol, + termination=IterateNTimes(1), + ), + WeightedAcquisition( + sampler=NormalAndUniformSampler( + 1000 * dim, + sigma=0.02, + strategy=SamplingStrategy.DDS, + ), + weightpattern=[ + 1.0, + 0.95, + 0.85, + 0.75, + 0.5, + 0.35, + 0.25, + 0.1, + 0.0, + ], + sigma_min=0.02, + sigma_max=0.02, + rtol=rtol, + termination=IterateNTimes(9), + ), + MaximizeDistance( + rtol=rtol, termination=IterateNTimes(1) + ), + ] ), - TransitionSearch(rtol=0.001, termination=IterateNTimes(9)), - MaximizeDistance(rtol=0.001, termination=IterateNTimes(1)), + MaximizeDistance(rtol=rtol), ] ) - if objTrained and evalTrained: - # Both surrogates are pre-trained, skip initial sampling - if disp: - print("Both surrogates pre-trained. Skipping initial sampling.") - out.fx = np.min(objSurrogate.Y) - out.x = objSurrogate.X[np.argmin(objSurrogate.Y)] - - elif objTrained and not evalTrained: - # Only obj surrogate trained - initialize eval surrogate - if disp: - print( - "Initializing evaluation surrogate from objective surrogate." - ) + # Create lists of points not in each surrogate + x_not_in_obj = np.empty((0, dim)) + x_not_in_eval = np.empty((0, dim)) + if evalSurrogate.ntrain == 0 and objSurrogate.ntrain > 0: + x_not_in_eval = objSurrogate.X + elif evalSurrogate.ntrain > 0 and objSurrogate.ntrain == 0: + x_not_in_obj = evalSurrogate.X[evalSurrogate.Y == 1] + elif evalSurrogate.ntrain > 0 and objSurrogate.ntrain > 0: + x_not_in_eval = FarEnoughSampleFilter(evalSurrogate.X, tol=rtol)( + objSurrogate.X + ) + x_not_in_obj = FarEnoughSampleFilter(objSurrogate.X, tol=rtol)( + evalSurrogate.X[evalSurrogate.Y == 1] + ) - evalSurrogate.update(objSurrogate.X, np.ones(len(objSurrogate.X))) - out.fx = np.min(objSurrogate.Y) - out.x = objSurrogate.X[np.argmin(objSurrogate.Y)] + # Reserve space for the surrogates + objSurrogate.reserve(objSurrogate.ntrain + maxeval, dim) + evalSurrogate.reserve( + evalSurrogate.ntrain + + maxeval + + len(x_not_in_eval) + - len(x_not_in_obj), + dim, + ) - elif not objTrained and evalTrained: - # Only eval surrogate trained - initialize obj surrogate - if disp: - print( - "Initializing objective surrogate from evaluation surrogate." - ) + # Update evalSurrogate with points from the objective surrogate + # Assumption: points in objSurrogate are enough to train evalSurrogate + if len(x_not_in_eval) > 0: + evalSurrogate.update(x_not_in_eval, np.ones(len(x_not_in_eval))) - # Extract all points in evalSurrogate with a value of 1 - evalPoints = evalSurrogate.X[evalSurrogate.Y == 1] + # Initialize output + # At this point, either + # + # - both evalSurrogate and objSurrogate are empty, or + # - evalSurrogate is initialized. + out = OptimizeResult() + out.init(fun, bounds, 1, maxeval, evalSurrogate) - # Evaluate each point - for x in evalPoints: - if out.nfev >= maxeval: - break + # Evaluate x_not_in_obj points and log results + if len(x_not_in_obj) > 0: + evaluate_and_log_point(fun, x_not_in_obj, out) - _ = evaluate_and_log_point(rescaledFunc, x, out) + # Initialize best values in out + out.init_best_values(objSurrogate) - if disp: - print("fEvals: %d" % out.nfev) - print("Best value: %f" % out.fx) + # Call the callback function with the current optimization result + if callback is not None: + callback(out) - # Check if more points are needed for the objective surrogate - if not objSurrogate.check_initial_design( - np.array( - out.sample[: out.nfev][~np.isnan(out.fsample[: out.nfev])] - ) - ): - maximizeDistance = MaximizeDistance(rtol=0.001) + # Keep adding points until there is a sufficient initial design for + # the objective surrogate + if objSurrogate.ntrain == 0: + while not objSurrogate.check_initial_design( + out.sample[: out.nfev][np.isfinite(out.fsample[: out.nfev])] + ) and (out.nfev < maxeval): if disp: print( - "Sampling additional points to initialize the surrogate model..." - ) - - while ( - not objSurrogate.check_initial_design( - np.array( - out.sample[: out.nfev][ - ~np.isnan(out.fsample[: out.nfev]) - ] - ) - ) - ) and (out.nfev < maxeval): - ## Generate new point - xNew = maximizeDistance.optimize( - evalSurrogate, [[0, 1] for _ in range(dim)], 1 - ) - - f = evaluate_and_log_point(rescaledFunc, xNew, out) - - if disp: - print("fEvals: %d" % out.nfev) - print("Best value: %f" % out.fx) - - # Update the surrogate model with the new point - evalSurrogate.update( - np.array(xNew), np.logical_not(np.isnan(f)).astype(float) + "Iteration: %d (Objective surrogate under construction)" + % out.nit ) - - else: - # Neither surrogate is trained - generate initial sample - if disp: - print("Performing initial sampling...") - - # Generate initial points using Latin Hypercube sampling - sampler = Sampler(nStart) - x0 = sampler.get_slhd_sample([[0, 1] for _ in range(dim)]) - - # Check that all points are far enough apart - distances = squareform(pdist(x0)) - keptIndices = [0] - - for i in range(1, len(x0)): - # Check if this point is far enough from all previously kept points - if all( - distances[i, kept_idx] >= 0.001 for kept_idx in keptIndices - ): - keptIndices.append(i) - - x0 = x0[np.array(keptIndices)] - - if disp: - print("Evaluating initial points...") - - # Evaluate the initial points - for x in x0: - if out.nfev >= maxeval: - break - - _ = evaluate_and_log_point(rescaledFunc, x, out) - - if disp: print("fEvals: %d" % out.nfev) - print("Best value: %f" % out.fx) - - # Build the evaluability surrogate model - evalSurrogate.update( - np.array(out.sample[0 : out.nfev, :]), - np.logical_not(np.isnan(out.fsample[0 : out.nfev])).astype(float), - ) - - # Check if initial design sufficient to initialize objective surrogate - if not objSurrogate.check_initial_design( - np.array( - out.sample[: out.nfev][~np.isnan(out.fsample[: out.nfev])] - ) - ): - maximizeDistance = MaximizeDistance(rtol=0.001) - if disp: print( - "Sampling additional points to initialize the surrogate model..." + "Number of feasible points: %d" + % np.sum(np.isfinite(out.fsample[: out.nfev])) ) - while ( - not objSurrogate.check_initial_design( - np.array( - out.sample[: out.nfev][ - ~np.isnan(out.fsample[: out.nfev]) - ] - ) - ) - ) and (out.nfev < maxeval): - ## Generate new point - xNew = maximizeDistance.optimize( - evalSurrogate, [[0, 1] for _ in range(dim)], 1 + dist = cdist(out.sample[: out.nfev], out.sample[: out.nfev]) + dist += np.eye(out.nfev) * np.max(dist) + print( + "Max distance between neighbors: %f" + % np.max(np.min(dist, axis=1)) ) + print(f"Last sampled point: {out.sample[out.nfev - 1]}") - f = evaluate_and_log_point(rescaledFunc, xNew, out) - - if disp: - print("fEvals: %d" % out.nfev) - print("Best value: %f" % out.fx) + # Acquire new sample point + xNew = MaximizeDistance(rtol=rtol).optimize( + objSurrogate, bounds, 1, points=out.sample[: out.nfev] + ) - # Update the surrogate model with the new point - evalSurrogate.update( - np.array(xNew), np.logical_not(np.isnan(f)).astype(float) - ) + # Compute f(xNew) and update out + evaluate_and_log_point(fun, xNew, out) + out.init_best_values(objSurrogate) + out.nit += 1 - # If we have run out of evaluations, we cannot continue - if out.nfev >= maxeval: - print( - "Maximum number of evaluations reached before enough points were sampled to initialize the surrogate model." - ) - return out + # Call the callback function + if callback is not None: + callback(out) - # Generate the surrogate model for the objective function - objSurrogate.update( - np.array(out.sample[: out.nfev][~np.isnan(out.fsample[: out.nfev])]), - np.array(out.fsample[: out.nfev][~np.isnan(out.fsample[: out.nfev])]), - ) + # Prepare for the main optimization loop + # + # - At this point, we have enough points to build both surrogates + nStart = out.nfev + nomadFunction = NomadProblem(fun, out) + xselected = np.array(out.sample[0 : out.nfev, :], copy=True) + ySelected = np.array(out.fsample[0 : out.nfev], copy=True) - if disp: - print( - "Objective surrogate model initialized with %d points." - % objSurrogate.ntrain - ) - print("Starting optimization search...") + # do until max number of f-evals reached or local min found + while out.nfev < maxeval: + if disp: + print("Iteration: %d" % out.nit) + print("fEvals: %d" % out.nfev) + print("Best value: %f" % out.fx) - # Call the callback function with the current optimization result - if callback is not None: - callback(out) + # Update surrogate models + t0 = time.time() + feasible_idx = np.isfinite(ySelected) + evalSurrogate.update(xselected, feasible_idx.astype(float)) + if np.any(feasible_idx): + objSurrogate.update( + xselected[feasible_idx], ySelected[feasible_idx] + ) + tf = time.time() + if disp: + print("Time to update surrogate model: %f s" % (tf - t0)) - # Main optimization loop - while out.nfev < maxeval: # Calculate the threshold for evaluability - threshold = np.log(max(1, out.nfev - nStart + 1)) / np.log( - maxeval - nStart + threshold = float( + np.log(max(1, out.nfev - nStart + 1)) / np.log(maxeval - nStart) ) - # Generate new point - xNew = acquisitionFunc.optimize( + # Set perturbation probability + if dim <= 10: + perturbProbability = 1.0 + else: + perturbProbability = np.random.uniform(0, 1) + + # Acquire new sample points + t0 = time.time() + xselected = acquisitionFunc.optimize( objSurrogate, - [[0, 1] for _ in range(dim)], - n=1, + bounds, + 1, points=evalSurrogate.X, - constraintTransform=lambda x: -x + threshold, - evaluabilitySurrogate=evalSurrogate, - evaluabilityThreshold=threshold, - scoreWeight=weightPattern[ - acquisitionFunc.acquisitionFuncArray[ - acquisitionFunc.idx - ].termination.iterationCount - ], + constr_fun=lambda x: threshold - evalSurrogate(x), + perturbation_probability=perturbProbability, ) - - # Evaluate new point - f = evaluate_and_log_point(rescaledFunc, xNew, out) - + if len(xselected) == 0: + threshold = float(np.finfo(float).eps) + tf = time.time() if disp: - print("fEvals: %d" % out.nfev) - print("Best value: %f" % out.fx) - - # Update evaluability surrogate model - evalSurrogate.update( - np.array(xNew), np.logical_not(np.isnan(f)).astype(float) - ) + print("Time to acquire new sample points: %f s" % (tf - t0)) + + # Compute f(xselected) + if len(xselected) > 0: + ySelected = evaluate_and_log_point(fun, xselected, out) + else: + ySelected = np.empty((0,)) + out.nit = out.nit + 1 + print( + "Acquisition function has failed to find a new sample! " + "Consider modifying it." + ) + break - # If successful, update the obj surrogate - if not np.isnan(f): - objSurrogate.update(np.array(xNew), np.array(f)) + # determine best one of newly sampled points + iSelectedBest = np.argmin(ySelected).item() + fxSelectedBest = ySelected[iSelectedBest] + if fxSelectedBest < out.fx: + out.x[:] = xselected[iSelectedBest, :] + out.fx = fxSelectedBest + new_best_point_found = True + else: + new_best_point_found = False # If the new point was better than current best, run NOMAD - if np.array_equiv(xNew, out.x): + if new_best_point_found: if disp: print("New best point found, running NOMAD...") nomadFunction.reset() - PyNomad.optimize( + res = PyNomad.optimize( fBB=nomadFunction, - pX0=xNew.flatten(), - pLB=np.array([0 for _ in range(dim)]), - pUB=np.array([1 for _ in range(dim)]), + pX0=out.x, + pLB=[b[0] for b in bounds], + pUB=[b[1] for b in bounds], params=[ "BB_OUTPUT_TYPE OBJ", f"MAX_BB_EVAL {min(4 * dim, maxeval - out.nfev)}", @@ -386,40 +363,54 @@ def rescaledFunc(x): ], ) + # Use the best point found by NOMAD + if res["f_best"] < out.fx: + out.x[:] = res["x_best"] + out.fx = res["f_best"] + # Get the points sampled by NOMAD - nomadSample = nomadFunction.get_x_history() - nomadFSample = nomadFunction.get_f_history() - - for i in range(len(nomadSample)): - point = nomadSample[i].reshape(1, -1) - # Check distance to all previously sampled points - if cdist(point, evalSurrogate.X).min() >= 1e-7: - fval = nomadFSample[i] - # Update surrogates - evalSurrogate.update( - point, np.logical_not(np.isnan(fval)).astype(float) - ) - if not np.isnan(fval): - objSurrogate.update(point, fval) + nomadSample = np.array(nomadFunction.get_x_history()) + nomadFSample = np.array(nomadFunction.get_f_history()) + if disp: print( f"NOMAD optimization completed. NOMAD used {len(nomadSample)} evaluations." ) - print("fEvals: %d" % out.nfev) - print("Best value: %f" % out.fx) - # Call the callback function with the current optimization result + # Filter out points that are too close to existing samples + idxes = FarEnoughSampleFilter( + np.vstack((evalSurrogate.X, xselected)), tol=rtol + ).indices(nomadSample) + xselected = np.vstack((xselected, nomadSample[idxes])) + ySelected = np.hstack((ySelected, nomadFSample[idxes])) + + # Update out.nit + out.nit = out.nit + 1 + + # Call the callback function if callback is not None: callback(out) - # Update the acquisition function + # Terminate if acquisition function has converged acquisitionFunc.update(out, objSurrogate) if acquisitionFunc.has_converged(): break - # Rescale the x and sample arrays to the original bounds - out.x = out.x * (bounds[:, 1] - bounds[:, 0]) + bounds[:, 0] - out.sample = out.sample * (bounds[:, 1] - bounds[:, 0]) + bounds[:, 0] + # Update output + out.sample.resize(out.nfev, dim) + out.fsample.resize(out.nfev) + + # Update surrogate model if it lives outside the function scope + if return_surrogate and evalSurrogate.ntrain > 0: + t0 = time.time() + feasible_idx = np.isfinite(ySelected) + evalSurrogate.update(xselected, feasible_idx.astype(float)) + if np.any(feasible_idx): + objSurrogate.update( + xselected[feasible_idx], ySelected[feasible_idx] + ) + tf = time.time() + if disp: + print("Time to update surrogate model: %f s" % (tf - t0)) - # Return OptimizeResult return out diff --git a/soogo/optimize/utils.py b/soogo/optimize/utils.py index a7a96134..9b6b5f2c 100644 --- a/soogo/optimize/utils.py +++ b/soogo/optimize/utils.py @@ -25,47 +25,35 @@ from .result import OptimizeResult -def evaluate_and_log_point(fun: Callable, x: np.ndarray, out: OptimizeResult): - """ - Evaluate a point or array of points and log the results. If the function - errors or the result is invalid (NaN or Inf), the output is logged as NaN. - If the function value is less than the current best, the current best ( - out.x, out.fx) is updated. Provided points are evaluated as a batch. This - function only supports single-objective functions. +def evaluate_and_log_point( + fun: Callable, x: np.ndarray, out: OptimizeResult +) -> np.ndarray: + """Evaluate an array of points and log the results in out. :param fun: The function to evaluate. - :param x: The point(s) to evaluate. Can be a 1D array (single point) or - 2D array (multiple points). + :param x: 2D array of points to evaluate. :param out: The output object to log the results. - :return: The function value(s) or NaN. Returns a scalar for single point, - array for multiple points. + :return: The function value(s) or NaN. """ - x = np.atleast_2d(x) + assert out.sample is not None and out.fsample is not None, ( + "Output object not initialized." + ) + # Evaluate function + n = len(x) try: - results = fun(x) - results = np.atleast_1d(results) + y = np.asarray(fun(x)) except Exception: - results = np.full(x.shape[0], np.nan) - - # Process each result individually - for i, y in enumerate(results): - if hasattr(y, "__len__") and len(y) > 0: - y = y[0] - if np.isnan(y) or np.isinf(y): - y = np.nan - results[i] = y - - out.sample[out.nfev, :] = x[i] - out.fsample[out.nfev] = y - out.nfev += 1 + shape_y = (n,) if out.fsample.ndim == 1 else (n, out.fsample.shape[1]) + y = np.full(shape_y, np.nan) - if not np.isnan(y) and (out.fx is None or y < out.fx): - out.x = x[i] - out.fx = y + # Log results + out.sample[out.nfev : out.nfev + n] = x + out.fsample[out.nfev : out.nfev + n] = y + out.nfev += n - return results[0] if len(results) == 1 else results + return y def uncertainty_score(candidates, points, fvals, k=3):