Solve a bug in the termination of the acquisition optimizers.

- The pymoo optimization now always terminate after 5000 func evals by default.
- Synced pool_size default with DEFAULT_N_MAX_EVALS_OPTIMIZER
- Don't use SpaceFillingSampler to generate candidates. It is too expensive.
- Use a low dispersion sequence (LHD) to generate candidates in the minimize_surrogate
- Improve GP kernel
- Avoid using non default rtol inside the algorithms

Author: Weslley da Silva Pereira

examples/gosac.ipynb

@@ -36,17 +36,14 @@
 
 from copy import deepcopy
 import importlib
-from math import sqrt
 from typing import Optional
 import numpy as np
 import matplotlib.pyplot as plt
 
 import soogo
 from soogo import OptimizeResult
-from soogo.sampling import dds_uniform_sample
 from soogo.acquisition import (
     CoordinatePerturbation,
-    TargetValueAcquisition,
     Acquisition,
     MinimizeSurrogate,
     MultipleAcquisition,
@@ -335,7 +332,7 @@ def main(config: int) -> list[OptimizeResult]:
         optres = read_and_run(
             data_file="datainput_BraninWithInteger",
             acquisitionFunc=CoordinatePerturbation(
-                sampler=dds_uniform_sample,
+                sampling_strategy="dds_uniform",
                 pool_size=200,
                 weightpattern=[0.3, 0.5, 0.8, 0.95],
                 sigma=BoundedParameter(0.2, 0.2 * 0.5**5, 0.2),
@@ -351,12 +348,7 @@ def main(config: int) -> list[OptimizeResult]:
     elif config == 5:
         optres = read_and_run(
             data_file="datainput_BraninWithInteger",
-            acquisitionFunc=MultipleAcquisition(
-                (
-                    TargetValueAcquisition(),
-                    MaximizeDistance(),
-                )
-            ),
+            acquisitionFunc=None,
             maxeval=100,
             Ntrials=3,
             batchSize=1,
@@ -388,10 +380,7 @@ def main(config: int) -> list[OptimizeResult]:
         optres = read_and_run(
             data_file="datainput_BraninWithInteger",
             acquisitionFunc=MultipleAcquisition(
-                (
-                    MinimizeSurrogate(100, 0.005 * sqrt(2)),
-                    MaximizeDistance(rtol=0.005 * sqrt(2)),
-                )
+                (MinimizeSurrogate(), MaximizeDistance())
             ),
             maxeval=100,
             Ntrials=3,

examples/legacy_opt/optimization_program_1.py

@@ -96,7 +96,6 @@
         acquisitionFunc=CoordinatePerturbation(
             pool_size=nCand,
             weightpattern=[0.3, 0.5],
-            rtol=1e-3,
             sigma=BoundedParameter(0.2, 0.2 * 0.5**5, 0.2),
             perturbation_strategy="fixed",
             termination=RobustCondition(UnsuccessfulImprovement(0.001), 5),

examples/legacy_opt/optimization_program_2.py

@@ -23,7 +23,6 @@
 __deprecated__ = False
 
 import os
-import numpy as np
 import pickle
 import time
 from benchmark import (
@@ -196,8 +195,8 @@ def run_optimizer(
     "optimizer": optimize.surrogate_optimization,
     "acquisition": acquisition.MultipleAcquisition(
         (
-            acquisition.MinimizeSurrogate(1, 0.005 * np.sqrt(2), seed=42),
-            acquisition.MaximizeDistance(rtol=0.005 * np.sqrt(2), seed=42),
+            acquisition.MinimizeSurrogate(seed=42),
+            acquisition.MaximizeDistance(seed=42),
         )
     ),
 }

examples/optimization.ipynb

@@ -30,6 +30,7 @@
     DefaultSingleObjectiveTermination,
     DefaultMultiObjectiveTermination,
 )
+from pymoo.util.display.multi import MultiObjectiveOutput
 
 # Local imports
 from ..model import Surrogate
@@ -89,7 +90,7 @@ class Acquisition(ABC):
     DEFAULT_RTOL = 1e-6
 
     #: Default maximum number of function evaluations for the acquisition
-    DEFAULT_N_MAX_EVALS_OPTIMIZER = 1000
+    DEFAULT_N_MAX_EVALS_OPTIMIZER = 5000
 
     def __init__(
         self,
@@ -135,39 +136,38 @@ def default_optimizer(
         :return: The default optimizer.
         """
         if not multi_objective:
+            termination = DefaultSingleObjectiveTermination(
+                xtol=self.rtol, n_max_evals=n_max_evals_optimizer
+            )
             if not mixed_integer:
-                return DE(
-                    termination=DefaultSingleObjectiveTermination(
-                        xtol=self.rtol, n_max_evals=n_max_evals_optimizer
-                    )
-                )
+                optim = DE()
+                optim.termination = termination
+                return optim
             else:
                 return MixedVariableGA(
                     eliminate_duplicates=ListDuplicateElimination(),
                     mating=MixedVariableMating(
                         eliminate_duplicates=ListDuplicateElimination()
                     ),
-                    termination=DefaultSingleObjectiveTermination(
-                        xtol=self.rtol, n_max_evals=n_max_evals_optimizer
-                    ),
+                    termination=termination,
                 )
         else:
+            termination = DefaultMultiObjectiveTermination(
+                xtol=self.rtol, n_max_evals=n_max_evals_optimizer
+            )
             if not mixed_integer:
-                return NSGA2(
-                    termination=DefaultMultiObjectiveTermination(
-                        xtol=self.rtol, n_max_evals=n_max_evals_optimizer
-                    )
-                )
+                optim = NSGA2()
+                optim.termination = termination
+                return optim
             else:
                 return MixedVariableGA(
                     eliminate_duplicates=ListDuplicateElimination(),
                     mating=MixedVariableMating(
                         eliminate_duplicates=ListDuplicateElimination()
                     ),
                     survival=RankAndCrowding(),
-                    termination=DefaultMultiObjectiveTermination(
-                        xtol=self.rtol, n_max_evals=n_max_evals_optimizer
-                    ),
+                    output=MultiObjectiveOutput(),
+                    termination=termination,
                 )
 
     @classmethod

examples/socemo.ipynb

@@ -27,8 +27,9 @@
 
 from .base import Acquisition
 from ..model import GaussianProcess
-from ..sampling import SpaceFillingSampler
+from ..sampling import random_sample
 from ..integrations.pymoo import PymooProblem
+from .utils import FarEnoughSampleFilter
 
 logger = logging.getLogger(__name__)
 
@@ -76,21 +77,18 @@ class MaximizeEI(Acquisition):
 
     def __init__(
         self,
-        sampler=None,
-        pool_size: int = 100,
+        pool_size: int = Acquisition.DEFAULT_N_MAX_EVALS_OPTIMIZER,
         avoid_clusters: bool = True,
+        n_max_evals_optimizer: Optional[int] = None,
         seed=None,
         **kwargs,
     ) -> None:
-        super().__init__(**kwargs)
+        super().__init__(
+            n_max_evals_optimizer=n_max_evals_optimizer or pool_size, **kwargs
+        )
         self.pool_size = pool_size
         self.avoid_clusters = avoid_clusters
         self.rng = np.random.default_rng(seed)
-        self.sampler = (
-            sampler
-            if sampler is not None
-            else SpaceFillingSampler(seed=self.rng)
-        )
 
     def optimize(
         self,
@@ -171,30 +169,14 @@ def optimize(
             xs = res.X
 
         # Returns xs if n == 1
-        if res.success and n == 1:
+        # print(f"MaximizeEI selected point with EI = {-res.F[0]}")
+        # print(f"At location: x = {res.X}")
+        # print(f"Success: {res.success}")
+        if n == 1:
             return np.asarray([xs])
 
         # Generate the complete pool of candidates
-        pool_size = self.pool_size
-        exclusion_set = (
-            np.vstack((exclusion_set, surrogateModel.X))
-            if exclusion_set is not None
-            else surrogateModel.X
-        )
-        if xs is not None:
-            exclusion_set = np.concatenate((exclusion_set, [xs]), axis=0)
-            pool_size -= 1
-        if xbest is not None:
-            exclusion_set = np.concatenate((exclusion_set, [xbest]), axis=0)
-            pool_size -= 1
-        x = (
-            self.sampler.generate(
-                self.pool_size, bounds, current_sample=exclusion_set
-            )
-            if pool_size > 0
-            else np.empty((0, dim))
-        )
-
+        x = random_sample(self.pool_size, bounds, seed=self.rng)
         if xs is not None:
             x = np.concatenate(([xs], x), axis=0)
         if xbest is not None:
@@ -287,4 +269,11 @@ def optimize(
             iBest[j] = np.argmax(score)
             eiCand[iBest[j]] = 0.0  # Remove this candidate expectancy
 
-        return x[iBest, :]
+        exclusion_set = (
+            np.vstack((exclusion_set, surrogateModel.X))
+            if exclusion_set is not None
+            else surrogateModel.X
+        )
+        return FarEnoughSampleFilter(exclusion_set, self.tol(bounds))(
+            x[iBest, :]
+        )

examples/vlse_benchmark/vlse_bench.py

@@ -32,10 +32,6 @@ class MinimizeMOSurrogate(Acquisition):
     """Obtain pareto-optimal sample points for the multi-objective surrogate
     model.
 
-    :param optimizer: Continuous multi-objective optimizer. If None, use
-        NSGA2 from pymoo.
-    :param mi_optimizer: Mixed-integer multi-objective optimizer. If None, use
-        MixedVariableGA from pymoo with RankAndCrowding survival strategy.
     :param seed: Seed for random number generator.
 
     .. attribute:: rng
@@ -44,12 +40,8 @@ class MinimizeMOSurrogate(Acquisition):
 
     """
 
-    def __init__(
-        self, optimizer=None, mi_optimizer=None, seed=None, **kwargs
-    ) -> None:
-        super().__init__(
-            optimizer, mi_optimizer, multi_objective=True, **kwargs
-        )
+    def __init__(self, seed=None, **kwargs) -> None:
+        super().__init__(multi_objective=True, **kwargs)
         self.rng = np.random.default_rng(seed)
 
     def optimize(

soogo/acquisition/base.py

@@ -23,6 +23,7 @@
 from scipy.spatial.distance import cdist
 from scipy.special import gamma
 from scipy.optimize import minimize
+from scipy.stats.qmc import LatinHypercube
 from typing import Optional
 
 from .base import Acquisition
@@ -76,8 +77,6 @@ class MinimizeSurrogate(Acquisition):
     collected for the new sample.
 
     :param pool_size: Number of uniform candidates generated each iteration.
-    :param rtol: Minimum distance between a candidate and already selected
-        points, relative to the domain size. Default is ``1e-3``.
     :param seed: Seed or random number generator.
 
     .. attribute:: pool_size
@@ -98,9 +97,12 @@ class MinimizeSurrogate(Acquisition):
     """
 
     def __init__(
-        self, pool_size: int, rtol: float = 1e-3, seed=None, **kwargs
+        self,
+        pool_size: int = Acquisition.DEFAULT_N_MAX_EVALS_OPTIMIZER,
+        seed=None,
+        **kwargs,
     ) -> None:
-        super().__init__(rtol=rtol, **kwargs)
+        super().__init__(**kwargs)
         self.pool_size = pool_size
         self.rng = np.random.default_rng(seed)
 
@@ -147,12 +149,15 @@ def optimize(
         )
 
         # Local space to store information
-        candidates = np.empty((self.pool_size * maxiter, dim))
-        distCandidates = np.empty(
-            (self.pool_size * maxiter, self.pool_size * maxiter)
+        candidates = random_sample(
+            self.pool_size,
+            bounds,
+            iindex=surrogateModel.iindex,
+            seed=LatinHypercube(d=dim, seed=self.rng),
         )
-        fcand = np.empty(self.pool_size * maxiter)
-        startpID = np.full((self.pool_size * maxiter,), False)
+        distCandidates = np.empty((self.pool_size, self.pool_size))
+        fcand = np.empty(self.pool_size)
+        startpID = np.full((self.pool_size,), False)
         selected = np.empty((n, dim))
 
         # Create a KDTree with the training data points
@@ -165,9 +170,10 @@ def optimize(
 
         iter = 0
         k = 0
+        iter_pool_size = int(np.ceil(self.pool_size / maxiter))
         while iter < maxiter and k < n and remevals > 0:
-            iStart = iter * self.pool_size
-            iEnd = (iter + 1) * self.pool_size
+            iStart = iter * iter_pool_size
+            iEnd = min((iter + 1) * iter_pool_size, self.pool_size)
 
             # if computational budget is exhausted, then return
             if remevals <= iEnd - iStart:
@@ -179,14 +185,6 @@ def optimize(
             # Consider only the best points to start local minimization
             counterLocalStart = iEnd // maxiter
 
-            # Choose candidate points uniformly in the search space
-            candidates[iStart:iEnd, :] = random_sample(
-                self.pool_size,
-                bounds,
-                iindex=surrogateModel.iindex,
-                seed=self.rng,
-            )
-
             # Compute the distance between the candidate points
             distCandidates[iStart:iEnd, iStart:iEnd] = cdist(
                 candidates[iStart:iEnd, :], candidates[iStart:iEnd, :]