Add parameter typing (only float and integer)

bayes_opt/bayesian_optimization.py

@@ -64,13 +64,13 @@ def dispatch(self, event):
 
 
 class BayesianOptimization(Observable):
-    def __init__(self, f, pbounds, random_state=None, verbose=2):
+    def __init__(self, f, pbounds, ptypes=None, random_state=None, verbose=2):
         """"""
         self._random_state = ensure_rng(random_state)
 
         # Data structure containing the function to be optimized, the bounds of
         # its domain, and a record of the evaluations we have done so far
-        self._space = TargetSpace(f, pbounds, random_state)
+        self._space = TargetSpace(f, pbounds, ptypes, random_state)
 
         # queue
         self._queue = Queue()
@@ -129,6 +129,7 @@ def suggest(self, utility_function):
             gp=self._gp,
             y_max=self._space.target.max(),
             bounds=self._space.bounds,
+            btypes=self._space.btypes,
             random_state=self._random_state
         )
 

bayes_opt/target_space.py

@@ -17,22 +17,26 @@ class TargetSpace(object):
     >>> def target_func(p1, p2):
     >>>     return p1 + p2
     >>> pbounds = {'p1': (0, 1), 'p2': (1, 100)}
-    >>> space = TargetSpace(target_func, pbounds, random_state=0)
+    >>> ptypes = {'p1': float, 'p2': int}
+    >>> space = TargetSpace(target_func, pbounds, ptypes, random_state=0)
     >>> x = space.random_points(1)[0]
     >>> y = space.register_point(x)
     >>> assert self.max_point()['max_val'] == y
     """
-    def __init__(self, target_func, pbounds, random_state=None):
+    def __init__(self, target_func, pbounds, ptypes=None,random_state=None):
         """
         Parameters
         ----------
         target_func : function
             Function to be maximized.
 
         pbounds : dict
-            Dictionary with parameters names as keys and a tuple with minimum
+            Dictionary with parameter names and list of the minimum and maximum boundaries
             and maximum values.
 
+        ptypes : dict
+            Dictionnary with parameter names and their type
+
         random_state : int, RandomState, or None
             optionally specify a seed for a random number generator
         """
@@ -44,10 +48,19 @@ def __init__(self, target_func, pbounds, random_state=None):
         # Get the name of the parameters
         self._keys = sorted(pbounds)
         # Create an array with parameters bounds
-        self._bounds = np.array(
-            [item[1] for item in sorted(pbounds.items(), key=lambda x: x[0])],
-            dtype=np.float
-        )
+        self._bounds = np.array([list(pbounds[item]) for item in self._keys], dtype=float)
+        # Create an array with the parameters type if declared
+        if ptypes is None:
+            self._btypes = None
+        else:
+            ## TODO: add exception if parameter names in btypes and ptypes do not have the same length and content
+            ## TODO: or store pbounds and ptypes has dictionnaries
+            try:
+                assert (len(ptypes) == len(pbounds))
+            except AssertionError:
+                raise AssertionError("ptypes and pbounds do not have same content."+\
+                                     "ptypes and pbounds must list exact same parameters")
+            self._btypes = np.array([ptypes[item] for item in self._keys], dtype=type)
 
         # preallocated memory for X and Y points
         self._params = np.empty(shape=(0, self.dim))
@@ -87,6 +100,10 @@ def keys(self):
     def bounds(self):
         return self._bounds
 
+    @property
+    def btypes(self):
+        return self._btypes
+
     def params_to_array(self, params):
         try:
             assert set(params) == set(self.keys)
@@ -142,11 +159,12 @@ def register(self, params, target):
 
         Notes
         -----
-        runs in ammortized constant time
+        runs in amortized constant time
 
         Example
         -------
         >>> pbounds = {'p1': (0, 1), 'p2': (1, 100)}
+        >>> ptypes = {'p1': float, 'p2':int}
         >>> space = TargetSpace(lambda p1, p2: p1 + p2, pbounds)
         >>> len(space)
         0
@@ -208,18 +226,26 @@ def random_sample(self):
         -------
         >>> target_func = lambda p1, p2: p1 + p2
         >>> pbounds = {'p1': (0, 1), 'p2': (1, 100)}
+        >>> ptypes = {'p1': float, 'p2':int}
         >>> space = TargetSpace(target_func, pbounds, random_state=0)
         >>> space.random_points(1)
-        array([[ 55.33253689,   0.54488318]])
+        array([[ 0.54488318, 55]])
         """
-        # TODO: support integer, category, and basic scipy.optimize constraints
+        # TODO: support category, and basic scipy.optimize constraints
         data = np.empty((1, self.dim))
-        for col, (lower, upper) in enumerate(self._bounds):
-            data.T[col] = self.random_state.uniform(lower, upper, size=1)
+        if self.btypes is None:
+            for col, (lower, upper) in enumerate(self._bounds):
+                data.T[col] = self.random_state.uniform(lower, upper, size=1)
+        else:
+            for col, (lower, upper) in enumerate(self._bounds):
+                if self.btypes[col] != int:
+                    data.T[col] = self.random_state.uniform(lower, upper, size=1)
+                if self.btypes[col] == int:
+                    data.T[col] = self.random_state.randint(int(lower), int(upper), size=1)
         return data.ravel()
 
     def max(self):
-        """Get maximum target value found and corresponding parametes."""
+        """Get maximum target value found and corresponding parameters."""
         try:
             res = {
                 'target': self.target.max(),
@@ -248,7 +274,19 @@ def set_bounds(self, new_bounds):
         ----------
         new_bounds : dict
             A dictionary with the parameter name and its new bounds
+
+        Returns
+        ----------
+        if type of modified parameter is int, then return rounded integer value
+        Example : new_bounds = {"p1", (1.2, 8.7)} and "p1" is integer
+        then new_bounds are (1,9)
         """
         for row, key in enumerate(self.keys):
             if key in new_bounds:
-                self._bounds[row] = new_bounds[key]
+                if self._btypes is not None:
+                    if self._btypes[row] == int:
+                        lbound = self._btypes[row](np.round(new_bounds[key][0], 0))
+                        ubound = self._btypes[row](np.round(new_bounds[key][1], 0))
+                        new_bounds[key] = (lbound, ubound)
+                    self._bounds[row] = list(new_bounds[key])
+                self._bounds[row] = list(new_bounds[key])

bayes_opt/util.py

@@ -3,8 +3,39 @@
 from scipy.stats import norm
 from scipy.optimize import minimize
 
+def generate_trials(n_events, bounds, btypes, random_state):
+    """A function to generate set of events under several constrains
 
-def acq_max(ac, gp, y_max, bounds, random_state, n_warmup=100000, n_iter=250):
+    Parameters
+    ----------
+    :param n_events:
+    The number of events to generate
+
+    :param bounds:
+    The variables bounds to limit the search of the acq max.
+
+    :param btypes:
+    The types of the variables.
+
+    :param random_state:
+    Instance of np.RandomState random number generator
+    """
+    x_trials = np.empty((n_events, bounds.shape[0]))
+    if btypes is None:
+        x_trials = random_state.uniform(bounds[:, 0], bounds[:, 1],
+                                        size=(n_events, bounds.shape[0]))
+    else:
+        for col, name in enumerate(bounds):
+            # print(col, name)
+            lower, upper = name
+            if btypes[col] != int:
+                x_trials[:, col] = random_state.uniform(lower, upper, size=n_events)
+            if btypes[col] == int:
+                x_trials[:, col] = random_state.randint(int(lower), int(upper), size=n_events)
+    return x_trials
+
+
+def acq_max(ac, gp, y_max, bounds, random_state, btypes=None, n_warmup=100000, n_iter=250):
     """
     A function to find the maximum of the acquisition function
 
@@ -26,6 +57,9 @@ def acq_max(ac, gp, y_max, bounds, random_state, n_warmup=100000, n_iter=250):
     :param bounds:
         The variables bounds to limit the search of the acq max.
 
+    :param btypes:
+        The types of the variables.
+
     :param random_state:
         instance of np.RandomState random number generator
 
@@ -39,20 +73,18 @@ def acq_max(ac, gp, y_max, bounds, random_state, n_warmup=100000, n_iter=250):
     -------
     :return: x_max, The arg max of the acquisition function.
     """
-
     # Warm up with random points
-    x_tries = random_state.uniform(bounds[:, 0], bounds[:, 1],
-                                   size=(n_warmup, bounds.shape[0]))
+    x_tries = generate_trials(n_warmup, bounds, btypes, random_state)
     ys = ac(x_tries, gp=gp, y_max=y_max)
     x_max = x_tries[ys.argmax()]
     max_acq = ys.max()
 
     # Explore the parameter space more throughly
-    x_seeds = random_state.uniform(bounds[:, 0], bounds[:, 1],
-                                   size=(n_iter, bounds.shape[0]))
+    x_seeds = generate_trials(n_iter, bounds, btypes, random_state)
     for x_try in x_seeds:
         # Find the minimum of minus the acquisition function
-        res = minimize(lambda x: -ac(x.reshape(1, -1), gp=gp, y_max=y_max),
+        ac_op = lambda x: -ac(x.reshape(1, -1), gp=gp, y_max=y_max)
+        res = minimize(ac_op,
                        x_try.reshape(1, -1),
                        bounds=bounds,
                        method="L-BFGS-B")
@@ -61,10 +93,36 @@ def acq_max(ac, gp, y_max, bounds, random_state, n_warmup=100000, n_iter=250):
         if not res.success:
             continue
 
+        # If integer in list of bounds
+        # search minimum between surroundings integers of the detected extremal point
+        if btypes is not None:
+            if int in btypes:
+                x_inf = res.x.copy()
+                x_sup = res.x.copy()
+                for i, (val, t) in enumerate(zip(res.x, btypes)):
+                    x_inf[i] = t(val)
+                    x_sup[i] = t(val + 1) if t == int else t(val)
+                # Store it if better than previous minimum(maximum).
+                x_ext = [x_inf, x_sup]
+                if max_acq is None or -res.fun[0] >= max_acq:
+                    max_acq = -1*np.minimum(ac_op(x_inf), ac_op(x_sup))
+                    x_argmax = np.argmin((ac_op(x_inf), ac_op(x_sup)))
+                    x_max = x_ext[x_argmax]
+            else:
+                # If only float in bounds
+                # store it if better than previous minimum(maximum).
+                if max_acq is None or -res.fun[0] >= max_acq:
+                    x_max = res.x
+                    max_acq = -res.fun[0]
+        else:
+            if max_acq is None or -res.fun[0] >= max_acq:
+                x_max = res.x
+                max_acq = -res.fun[0]
+
         # Store it if better than previous minimum(maximum).
-        if max_acq is None or -res.fun[0] >= max_acq:
-            x_max = res.x
-            max_acq = -res.fun[0]
+        # if max_acq is None or -res.fun[0] >= max_acq:
+        #     x_max = res.x
+        #     max_acq = -res.fun[0]
 
     # Clip output to make sure it lies within the bounds. Due to floating
     # point technicalities this is not always the case.

examples/basic-tour.ipynb

@@ -1,20 +1,5 @@
 {
  "cells": [
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# Basic tour of the Bayesian Optimization package\n",
-    "\n",
-    "This is a constrained global optimization package built upon bayesian inference and gaussian process, that attempts to find the maximum value of an unknown function in as few iterations as possible. This technique is particularly suited for optimization of high cost functions, situations where the balance between exploration and exploitation is important.\n",
-    "\n",
-    "Bayesian optimization works by constructing a posterior distribution of functions (gaussian process) that best describes the function you want to optimize. As the number of observations grows, the posterior distribution improves, and the algorithm becomes more certain of which regions in parameter space are worth exploring and which are not, as seen in the picture below.\n",
-    "\n",
-    "As you iterate over and over, the algorithm balances its needs of exploration and exploitation taking into account what it knows about the target function. At each step a Gaussian Process is fitted to the known samples (points previously explored), and the posterior distribution, combined with a exploration strategy (such as UCB (Upper Confidence Bound), or EI (Expected Improvement)), are used to determine the next point that should be explored (see the gif below).\n",
-    "\n",
-    "This process is designed to minimize the number of steps required to find a combination of parameters that are close to the optimal combination. To do so, this method uses a proxy optimization problem (finding the maximum of the acquisition function) that, albeit still a hard problem, is cheaper (in the computational sense) and common tools can be employed. Therefore Bayesian Optimization is most adequate for situations where sampling the function to be optimized is a very expensive endeavor. See the references for a proper discussion of this method."
-   ]
-  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -42,6 +27,21 @@
     "    return -x ** 2 - (y - 1) ** 2 + 1"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Basic tour of the Bayesian Optimization package\n",
+    "\n",
+    "This is a constrained global optimization package built upon bayesian inference and gaussian process, that attempts to find the maximum value of an unknown function in as few iterations as possible. This technique is particularly suited for optimization of high cost functions, situations where the balance between exploration and exploitation is important.\n",
+    "\n",
+    "Bayesian optimization works by constructing a posterior distribution of functions (gaussian process) that best describes the function you want to optimize. As the number of observations grows, the posterior distribution improves, and the algorithm becomes more certain of which regions in parameter space are worth exploring and which are not, as seen in the picture below.\n",
+    "\n",
+    "As you iterate over and over, the algorithm balances its needs of exploration and exploitation taking into account what it knows about the target function. At each step a Gaussian Process is fitted to the known samples (points previously explored), and the posterior distribution, combined with a exploration strategy (such as UCB (Upper Confidence Bound), or EI (Expected Improvement)), are used to determine the next point that should be explored (see the gif below).\n",
+    "\n",
+    "This process is designed to minimize the number of steps required to find a combination of parameters that are close to the optimal combination. To do so, this method uses a proxy optimization problem (finding the maximum of the acquisition function) that, albeit still a hard problem, is cheaper (in the computational sense) and common tools can be employed. Therefore Bayesian Optimization is most adequate for situations where sampling the function to be optimized is a very expensive endeavor. See the references for a proper discussion of this method."
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},

examples/bo_parameterTyping_example.py

@@ -0,0 +1,20 @@
+from bayes_opt import BayesianOptimization
+
+# function to be maximized - must find (x=0;y=10)
+targetFunction = lambda x, y: -(x-0.5) ** 2 - (y - 10) ** 2 + 1
+
+# define parameters bounds
+bounds = {'y': (5, 15), 'x': (-3, 3)}
+btypes = {'y':int, 'x':float}
+bo = BayesianOptimization(targetFunction, bounds) #, ptypes=btypes)
+
+bo.probe({"x": 1.4, "y": 6})
+bo.probe({"x": 2.4, "y": 12})
+bo.probe({"x": -2.4, "y": 13})
+
+bo.maximize(init_points=5, n_iter=5, kappa=2)
+
+# print results
+print(f'Estimated position of the maximum: {bo.max}')
+print(f'List of tested positions:\n{bo.res}')
+

examples/sklearn_example.py

@@ -70,7 +70,8 @@ def svc_crossval(expC, expGamma):
 
     optimizer = BayesianOptimization(
         f=svc_crossval,
-        pbounds={"expC": (-3, 2), "expGamma": (-4, -1)},
+        pbounds={"expC": [float, (-3, 2)], "expGamma": [float, (-4, -1)]},
+        ptypes={"expC": float, "expGamma":float},
         random_state=1234,
         verbose=2
     )
@@ -90,8 +91,8 @@ def rfc_crossval(n_estimators, min_samples_split, max_features):
         accordingly.
         """
         return rfc_cv(
-            n_estimators=int(n_estimators),
-            min_samples_split=int(min_samples_split),
+            n_estimators=n_estimators,
+            min_samples_split=min_samples_split,
             max_features=max(min(max_features, 0.999), 1e-3),
             data=data,
             targets=targets,
@@ -102,7 +103,11 @@ def rfc_crossval(n_estimators, min_samples_split, max_features):
         pbounds={
             "n_estimators": (10, 250),
             "min_samples_split": (2, 25),
-            "max_features": (0.1, 0.999),
+            "max_features": (0.1, 0.999)
+        },
+        ptypes={"n_estimators": int,
+            "min_samples_split": int,
+            "max_features": float
         },
         random_state=1234,
         verbose=2