Kirill latest

bayes_optim/acquisition/acquisition_fun.py

@@ -4,9 +4,9 @@
 from typing import Tuple, Union
 
 import numpy as np
+from bayes_optim.mylogging import eprintf
 from numpy import exp, sqrt
 from scipy.stats import norm
-from bayes_optim.mylogging import eprintf
 
 from ..surrogate import GaussianProcess, RandomForest
 
@@ -156,34 +156,44 @@ def __call__(
     ) -> Union[np.ndarray, Tuple[np.ndarray, np.ndarray]]:
         X = self.check_X(X)
         n_sample = X.shape[0]
-        y_hat, sd = self._predict(X)
+        xi = 0.0
 
+        if return_dx:
+            mu, std, mu_grad, std_grad = self._model.predict(
+                X, return_std=True, return_mean_grad=True, return_std_grad=True
+            )
+        else:
+            mu, std = self._model.predict(X, return_std=True)
         # if the Kriging variance is to small
         # TODO: check the rationale of 1e-6 and why the ratio if intended
         if hasattr(self._model, "sigma2"):
-            if sd / np.sqrt(self._model.sigma2) < 1e-6:
+            if std / np.sqrt(self._model.sigma2) < 1e-6:
                 return (0, np.zeros((len(X[0]), 1))) if return_dx else 0
         else:
             # TODO: implement a counterpart of 'sigma2' for randomforest
             # or simply put a try...except around the code below
-            if sd < 1e-10:
+            if std < 1e-10:
                 return (0, np.zeros((len(X[0]), 1))) if return_dx else 0
         try:
-            xcr_ = self.plugin - y_hat
-            xcr = xcr_ / sd
-            xcr_prob, xcr_dens = norm.cdf(xcr), norm.pdf(xcr)
-            value = xcr_ * xcr_prob + sd * xcr_dens
+            impr = self.plugin - xi - mu
+            xcr = impr / std
+            cdf, pdf = norm.cdf(xcr), norm.pdf(xcr)
+            value = impr * cdf + std * pdf
             if n_sample == 1:
                 value = sum(value)
         except Exception:  # in case of numerical errors
             # IMPORTANT: always keep the output in the same type
             value = 0
 
         if return_dx:
-            y_dx, sd2_dx = self._gradient(X)
-            sd_dx = sd2_dx / (2.0 * sd)
             try:
-                dx = -y_dx * xcr_prob + sd_dx * xcr_dens
+                improve_grad = -mu_grad * std - std_grad * impr
+                improve_grad /= std**2
+                cdf_grad = improve_grad * pdf
+                pdf_grad = -impr * cdf_grad
+                exploit_grad = -mu_grad * cdf - pdf_grad
+                explore_grad = std_grad * pdf + pdf_grad
+                dx = exploit_grad + explore_grad
             except Exception:
                 dx = np.zeros((len(X[0]), 1))
             return value, dx

bayes_optim/acquisition/optim/__init__.py

@@ -1,5 +1,6 @@
 import logging
 import time
+from abc import ABC, abstractmethod
 from typing import Callable, List, Tuple
 
 import numpy as np
@@ -19,6 +20,7 @@
     "default_AQ_max_FEs",
     "default_AQ_n_restart",
     "default_AQ_wait_iter",
+    "OptimizationListener",
 ]
 
 
@@ -52,16 +54,24 @@ def __call__(self, x: np.ndarray) -> Tuple[float, np.ndarray]:
         return f, fg
 
 
+class OptimizationListener(ABC):
+    @abstractmethod
+    def on_optimum_found(self, fopt, xopt):
+        pass
+
+
 def argmax_restart(
     obj_func: Callable,
     search_space,
+    obj_func_: Callable = None,
     h: Callable = None,
     g: Callable = None,
     eval_budget: int = 100,
     n_restart: int = 10,
     wait_iter: int = 3,
     optimizer: str = "BFGS",
     logger: logging.Logger = None,
+    listener: OptimizationListener = None,
 ):
     # lists of the best solutions and acquisition values from each restart
     xopt, fopt = [], []
@@ -71,16 +81,25 @@ def argmax_restart(
         optimizer = "MIES"
         logger.warning("L-BFGS-B can only be applied on continuous search space")
 
+    X = search_space.sample(int(2e2 * search_space.dim), method="uniform")
+    if obj_func_ is not None:
+        af_value = [obj_func_(x) for x in X]
+    else:
+        af_value = [obj_func(x)[0] for x in X]
+
+    idx = np.argsort(af_value)[::-1]
+    X = X[idx, :]
+
     for iteration in range(n_restart):
-        x0 = search_space.sample(N=1, method="uniform")[0]
+        x0 = X[iteration]
         if optimizer == "BFGS":
             bounds = np.array(search_space.bounds)
             # TODO: this is still subject to testing
             xopt_, fopt_, stop_dict = fmin_l_bfgs_b(
                 Penalized(obj_func, h, g),
                 x0,
-                pgtol=1e-8,
-                factr=1e6,
+                factr=1e5,
+                approx_grad=False,
                 bounds=bounds,
                 maxfun=eval_budget,
             )
@@ -90,7 +109,7 @@ def argmax_restart(
             fopt_ = -fopt_
 
             if logger is not None and stop_dict["warnflag"] != 0:
-                logger.debug("L-BFGS-B terminated abnormally with the state: %s" % stop_dict)
+                logger.info("L-BFGS-B terminated abnormally with the state: %s" % stop_dict)
 
         elif optimizer == "OnePlusOne_Cholesky_CMA":
             opt = OnePlusOne_Cholesky_CMA(
@@ -132,9 +151,8 @@ def argmax_restart(
                 wait_count = 0
 
                 if logger is not None:
-                    logger.debug(
-                        "restart : %d - funcalls : %d - Fopt : %f"
-                        % (iteration + 1, stop_dict["funcalls"], fopt_)
+                    logger.info(
+                        f"restart : {iteration + 1} - funcalls : {stop_dict['funcalls']} - Fopt : {fopt_:.8e}"
                     )
             else:
                 wait_count += 1
@@ -143,7 +161,11 @@ def argmax_restart(
             xopt.append(xopt_)
             fopt.append(fopt_)
 
-        if eval_budget <= 0 or wait_count >= wait_iter:
+            if listener is not None:
+                listener.on_optimum_found(fopt_, xopt_)
+
+        if eval_budget <= 0:
+            # or wait_count >= wait_iter:
             break
 
     if len(xopt) == 0:

bayes_optim/acquisition/optim/one_plus_one_cma_es.py

@@ -7,6 +7,8 @@
 
 from ...search_space import RealSpace, SearchSpace
 from ...utils import dynamic_penalty, get_logger, handle_box_constraint, set_bounds
+from bayes_optim.mylogging import eprintf
+import random
 
 Vector = List[float]
 Matrix = List[Vector]
@@ -182,6 +184,7 @@ def random_seed(self, seed):
             self._random_seed = int(seed)
             if self._random_seed:
                 np.random.seed(self._random_seed)
+                random.seed(self._random_seed)
 
     @property
     def C(self):
@@ -213,7 +216,7 @@ def x(self, x):
         else:
             # sample `x` u.a.r. in `[lb, ub]`
             assert all(~np.isinf(self.lb)) & all(~np.isinf(self.ub))
-            x = (self.ub - self.lb) * np.random.rand(self.dim) + self.lb
+            x = np.array([random.uniform(self.lb[i], self.ub[i]) for i in range(self.dim)])
 
         self._x = x
         y = self.evaluate(x)

bayes_optim/acquisition/optim/option.py

@@ -6,8 +6,8 @@
     "MIES": lambda dim: int(1e3 * dim),
     "BFGS": lambda dim: int(1e2 * dim),
     "OnePlusOne_Cholesky_CMA": lambda dim: int(1e3 * dim),
-    "DE": lambda dim: int(1e3 * dim)
+    "DE": lambda dim: int(1e3 * dim),
 }
 
-default_AQ_n_restart: Callable = lambda dim: int(5 * dim)
+default_AQ_n_restart: Callable = lambda dim: max(5, int(dim / 2))
 default_AQ_wait_iter: int = 3