From d3ddec19d59ce83650e53e8f8bd85e488375d01b Mon Sep 17 00:00:00 2001 From: Jan Rittig <63640094+jangerit@users.noreply.github.com> Date: Mon, 8 Jun 2020 02:34:36 +0200 Subject: [PATCH] Add single-objective Bayesian Optimization (#38) Co-authored-by: Kobi Felton --- poetry.lock | 39 ++-- pyproject.toml | 1 + summit/benchmarks/test_functions.py | 1 - summit/strategies/__init__.py | 1 + summit/strategies/sobo.py | 286 ++++++++++++++++++++++++++++ tests/test_strategies.py | 56 ++++++ 6 files changed, 365 insertions(+), 19 deletions(-) create mode 100644 summit/strategies/sobo.py diff --git a/poetry.lock b/poetry.lock index 7e5e523b..acb706bb 100644 --- a/poetry.lock +++ b/poetry.lock @@ -119,6 +119,23 @@ notebook = ["jupyter-client (>=4.0.6)", "ipywidgets (>=4.0.3)", "ipykernel (>=4. optional = ["mpi4py", "ipython (>=4.0.0)"] plotting = ["matplotlib (>=3.0)", "plotly (>=1.8.6)"] +[[package]] +category = "main" +description = "The Bayesian Optimization Toolbox" +name = "gpyopt" +optional = false +python-versions = "*" +version = "1.2.6" + +[package.dependencies] +GPy = ">=1.8" +numpy = ">=1.7" +scipy = ">=0.16" + +[package.extras] +docs = ["matplotlib (>=1.3)", "sphinx", "ipython"] +optimizer = ["direct", "cma", "pydoe", "sobol-seq", "emcee"] + [[package]] category = "dev" description = "Internationalized Domain Names in Applications (IDNA)" @@ -729,22 +746,6 @@ threadpoolctl = ">=2.0.0" [package.extras] alldeps = ["numpy (>=1.13.3)", "scipy (>=0.19.1)"] -[[package]] -category = "main" -description = "Integrator for Python-based quantum computing software" -name = "scikit-quant" -optional = false -python-versions = "*" -version = "0.6.1" - -[package.dependencies] -Py-BOBYQA = ">=1.1" - -[package.source] -reference = "84e158b0f902eb96167e93a6603103feb4f58d52" -type = "git" -url = "https://github.com/sustainable-processes/scikit-quant.git" - [[package]] category = "main" description = "SciPy: Scientific Library for Python" @@ -905,7 +906,7 @@ docs = ["sphinx", "jaraco.packaging (>=3.2)", "rst.linker (>=1.9)"] testing = ["jaraco.itertools", "func-timeout"] [metadata] -content-hash = "285ad44b8cc6b564d14261ce6363b7cf99a9277c6760d90f4d8ba5c5b8cdc34b" +content-hash = "8cc7cb847e0e7b7adaa7078fc6d50540ea078d65853beac6d16d563696e59611" python-versions = "^3.6" [metadata.files] @@ -968,6 +969,9 @@ gpy = [ {file = "GPy-1.9.9.win-amd64-py3.6.exe", hash = "sha256:9231368ce5f79aadbf4e06a805742c1376cdd68326f613091f9ebca7aa113fbb"}, {file = "GPy-1.9.9.win-amd64-py3.7.exe", hash = "sha256:138ac1574beabe321a2bf84615a3c00c2d4e0085a4ed74afb9ef8c086804d890"}, ] +gpyopt = [ + {file = "GPyOpt-1.2.6.tar.gz", hash = "sha256:e714daa035bb529a6db23c53665a762a4ab3456b9329c19ad3b03983f94c9b2a"}, +] idna = [ {file = "idna-2.9-py2.py3-none-any.whl", hash = "sha256:a068a21ceac8a4d63dbfd964670474107f541babbd2250d61922f029858365fa"}, {file = "idna-2.9.tar.gz", hash = "sha256:7588d1c14ae4c77d74036e8c22ff447b26d0fde8f007354fd48a7814db15b7cb"}, @@ -1274,7 +1278,6 @@ scikit-learn = [ {file = "scikit_learn-0.23.1-cp38-cp38-win32.whl", hash = "sha256:5bcea4d6ee431c814261117281363208408aa4e665633655895feb059021aca6"}, {file = "scikit_learn-0.23.1-cp38-cp38-win_amd64.whl", hash = "sha256:16feae4361be6b299d4d08df5a30956b4bfc8eadf173fe9258f6d59630f851d4"}, ] -scikit-quant = [] scipy = [ {file = "scipy-1.4.1-cp35-cp35m-macosx_10_6_intel.whl", hash = "sha256:c5cac0c0387272ee0e789e94a570ac51deb01c796b37fb2aad1fb13f85e2f97d"}, {file = "scipy-1.4.1-cp35-cp35m-manylinux1_i686.whl", hash = "sha256:a144811318853a23d32a07bc7fd5561ff0cac5da643d96ed94a4ffe967d89672"}, diff --git a/pyproject.toml b/pyproject.toml index 139c12e0..47c67c2e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -17,6 +17,7 @@ numpy = "1.16.0" platypus-opt = "^1.0" sklearn = "^0.0.0" SQSnobFit = "^0.4.3" +gpyopt = "^1.2.6" [tool.poetry.dev-dependencies] pytest = "^3.0" diff --git a/summit/benchmarks/test_functions.py b/summit/benchmarks/test_functions.py index f1e92fe1..de79bd58 100644 --- a/summit/benchmarks/test_functions.py +++ b/summit/benchmarks/test_functions.py @@ -292,4 +292,3 @@ def plot(self, **kwargs): plt.show() plt.close() - diff --git a/summit/strategies/__init__.py b/summit/strategies/__init__.py index e97f8161..ac5a9886 100644 --- a/summit/strategies/__init__.py +++ b/summit/strategies/__init__.py @@ -3,3 +3,4 @@ from .tsemo import TSEMO2 from .neldermead import NelderMead from .snobfit import SNOBFIT +from .sobo import SOBO diff --git a/summit/strategies/sobo.py b/summit/strategies/sobo.py new file mode 100644 index 00000000..039be39f --- /dev/null +++ b/summit/strategies/sobo.py @@ -0,0 +1,286 @@ +from .base import Strategy +from .random import LHS +from summit.domain import Domain, DomainError +from summit.utils.dataset import DataSet + +import GPy +import GPyOpt + +import numpy as np +import pandas as pd +from abc import ABC, abstractmethod + +class SOBO(Strategy): + ''' Single-objective Bayesian Optimization (SOBO) + + Parameters + ---------- + domain: summit.domain.Domain + The Summit domain describing the optimization problem. + gp_model_type: string, optional + The GPy Gaussian Process model type. + By default, gaussian processes with the Matern 5.2 kernel will be used. + acquisition_type: string, optional + The acquisition function type from GPyOpt. + By default, Excpected Improvement (EI). + optimizer_type: string, optional + The internal optimizer used in GPyOpt for maximization of the acquisition function. + By default, lfbgs will be used. + evaluator_type: string, optional + The evaluator type used for batch mode (how multiple points are chosen in one iteration). + By default, thompson sampling will be used. + kernel: GPy kernel object, optional + The kernel used in the GP. + By default a Matern 5.2 kernel (GPy object) will be used. + exact_feval: boolean, optional + Whether the function evaluations are exact (True) or noisy (False). + By default: False. + ard: boolean, optional + Whether automatic relevance determination should be applied (True). + By default: True. + standardize_outputs: boolean, optional + Whether the outputs should be standardized (True). + By default: True. + + Notes + ---------- + This implementation uses the python package GPyOpt provided by + the Machine Learning Group of the University of Sheffield. + + Copyright (c) 2016, the GPyOpt Authors. + All rights reserved. + + Github repository: https://github.com/SheffieldML/GPyOpt + Homepage: http://sheffieldml.github.io/GPyOpt/ + + Please cite their work, when using this strategy. + + Examples + -------- + >>> from summit.domain import Domain, ContinuousVariable, DiscreteVariable + >>> from summit.strategies import SOBO + >>> import numpy as np + >>> domain = Domain() + >>> domain += ContinuousVariable(name='temperature', description='reaction temperature in celsius', bounds=[50, 100]) + >>> domain += DiscreteVariable(name='flowrate_a', description='flow of reactant a in mL/min', levels=[1,2,3,4,5]) + >>> domain += ContinuousVariable(name='flowrate_b', description='flow of reactant b in mL/min', bounds=[0.1, 0.5]) + >>> strategy = SOBO(domain) + >>> result, xbest, fbest, param = strategy.suggest_experiments(5) + + ''' + + def __init__(self, domain, gp_model_type=None, acquisition_type=None, optimizer_type=None, evaluator_type=None, **kwargs): + Strategy.__init__(self, domain) + + # TODO: notation - discrete in our model (e.g., catalyst type) = categorical? + self.input_domain = [] + for v in self.domain.variables: + if not v.is_objective: + if v.variable_type == 'continuous': + self.input_domain.append( + {'name': v.name, + 'type': v.variable_type, + 'domain': (v.bounds[0], v.bounds[1])}) + elif v.variable_type == 'discrete': + self.input_domain.append( + {'name': v.name, + 'type': 'discrete', + 'domain': tuple(v.levels)}) + # TODO: GPyOpt currently does not support mixed-domains w/ bandit inputs, there is a PR for this though + elif v.variable_type == 'descriptors': + self.input_domain.append({'name': v.name, + 'type': 'bandit', + 'domain': [tuple(t) for t in v.ds.data_to_numpy().tolist()]}) + + ''' possible workaround for mixed-type variable problems: treat descriptor as categorical variables + self.input_domain.append({'name': v.name, + 'type': 'categorical', + 'domain': tuple(np.arange(v.ds.data_to_numpy().shape[0]).tolist())}) + ''' + else: + raise TypeError('Unknown variable type.') + + # TODO: how to handle equality constraints? Could we remove '==' from constraint types as each equality + # constraint reduces the degrees of freedom? + if self.domain.constraints is not None: + constraints = self.constr_wrapper(self.domain) + self.constraints = [{'name': 'constr_' + str(i), + 'constraint': c[0] if c[1] in ['<=', '<'] else '(' + c[0] + ')*(-1)'} + for i,c in enumerate(constraints) if not (c[1] == '==')] + else: + self.constraints = None + + self.input_dim = self.domain.num_continuous_dimensions() + self.domain.num_discrete_variables() + + """ + Gaussian Process (GP) model + GP: standard Gaussian Process + GP_MCMC: Gaussian Process with prior in hyperparameters + sparseGP: sparse Gaussian Process + warpedGP: warped Gaussian Process + InputWarpedGP: input warped Gaussian Process + RF: random forest (scikit-learn) + """ + + if gp_model_type in ['GP', 'GP_MCMC', 'sparseGP', 'warpedGP', 'InputWarpedGP', 'RF']: + self.gp_model_type = gp_model_type + else: + self.gp_model_type = 'GP' # default model type is a standard Gaussian Process (from GPy package) + + """ + Acquisition function type + EI: expected improvement + EI_MCMC: integrated expected improvement (requires GP_MCMC model) (https://dash.harvard.edu/bitstream/handle/1/11708816/snoek-bayesopt-nips-2012.pdf?sequence%3D1) + LCB: lower confidence bound + LCB_MCMC: integrated GP-Lower confidence bound (requires GP_MCMC model) + MPI: maximum probability of improvement + MPI_MCMC: maximum probability of improvement (requires GP_MCMC model) + LP: local penalization + ES: entropy search + """ + if acquisition_type in ['EI', 'EI_MCMC', 'LCB', 'LCB_MCMC', 'MPI', 'MPI_MCMC', 'LP', 'ES']: + self.acquisition_type = acquisition_type + else: + self.acquisition_type = 'EI' # default acquisition function is expected utility improvement + + """ + Method for optimization of acquisition function + lbfgs: Limited-memory Broyden–Fletcher–Goldfarb–Shanno, + DIRECT: Dividing Rectangles, + CMA: covariance matrix adaption + """ + if optimizer_type in ['lbfgs', 'DIRECT', 'CMA']: + self.optimizer_type = optimizer_type + else: + self.optimizer_type = 'lbfgs' # default optimizer: lbfgs + + if evaluator_type in ['sequential', 'random', 'local_penalization', 'thompson_sampling']: + self.evaluator_type = evaluator_type + else: + self.evaluator_type = 'random' + + # specify GPy kernel: # https://gpy.readthedocs.io/en/deploy/GPy.kern.html#subpackages + self.kernel = kwargs.get('kernel', GPy.kern.Matern52(self.input_dim)) + # Are function values exact (w/o noise)? + self.exact_feval = kwargs.get('exact_feval', False) + # automatic relevance determination + self.ARD = kwargs.get('ARD', True) + # Standardization of outputs? + self.standardize_outputs = kwargs.get('standardize_outputs', True) + + + def suggest_experiments(self, num_experiments, + prev_res: DataSet=None, prev_param=None): + """ Suggest experiments using GPyOpt single-objective Bayesian Optimization + + Parameters + ---------- + num_experiments: int + The number of experiments (i.e., samples) to generate + prev_res: summit.utils.data.DataSet, optional + Dataset with data from previous experiments of previous iteration. + If no data is passed, then random sampling will + be used to suggest an initial design. + prev_param: array-like, optional TODO: how to handle this? + File with results from previous iterations of SOBO algorithm. + If no data is passed, only results from prev_res will be used. + + Returns + ------- + next_experiments + A `Dataset` object with points to be evaluated next + xbest + Best point from all iterations. + fbest + Objective value at best point from all iterations. + param + A list containing all evaluated X and corresponding Y values. + """ + + param = None + xbest = None + fbest = float("inf") + + # Suggest random initial design + if prev_res is None: + '''lhs design does not consider constraints + lhs = LHS(self.domain) + next_experiments = lhs.suggest_experiments((num_experiments)) + return next_experiments, None, float("inf"), None + ''' + feasible_region = GPyOpt.Design_space(space=self.input_domain, constraints=self.constraints) + request = GPyOpt.experiment_design.initial_design('random', feasible_region, num_experiments) + + else: + # Get inputs and outputs + inputs, outputs = self.transform.transform_inputs_outputs(prev_res) + + # Set up maximization and minimization by converting maximization to minimization problem + for v in self.domain.variables: + if v.is_objective and v.maximize: + outputs[v.name] = -1 * outputs[v.name] + + inputs = inputs.to_numpy() + outputs = outputs.to_numpy() + + if prev_param is not None: + X_step = prev_param[0] + Y_step = prev_param[1] + + X_step = np.vstack((X_step, inputs)) + Y_step = np.vstack((Y_step, outputs)) + + else: + X_step = inputs + Y_step = outputs + + sobo_model = GPyOpt.methods.BayesianOptimization(f=None, + domain=self.input_domain, + constraints=self.constraints, + model_type=self.gp_model_type, + kernel=self.kernel, + acquisition_type=self.acquisition_type, + acquisition_optimizer_type=self.optimizer_type, + normalize_Y=self.standardize_outputs, + batch_size=num_experiments, + evaluator_type=self.evaluator_type, + maximize=False, + ARD=self.ARD, + exact_feval=self.exact_feval, + X=X_step, + Y=Y_step) + request = sobo_model.suggest_next_locations() + + + # Store parameters (history of suggested points and function evaluations) + param = [X_step, Y_step] + + fbest = np.min(Y_step) + xbest = X_step[np.argmin(Y_step)] + + + # Generate DataSet object with variable values of next + next_experiments = None + if request is not None and len(request)!=0: + next_experiments = {} + i_inp = 0 + for v in self.domain.variables: + if not v.is_objective: + next_experiments[v.name] = request[:, i_inp] + i_inp += 1 + next_experiments = DataSet.from_df(pd.DataFrame(data=next_experiments)) + next_experiments[('strategy', 'METADATA')] = 'Single-objective BayOpt' + + return next_experiments, xbest, fbest, param + + + def constr_wrapper(self, summit_domain): + v_input_names = [v.name for v in summit_domain.variables if not v.is_objective] + gpyopt_constraints = [] + for c in summit_domain.constraints: + tmp_c = c.lhs + for v_input_index, v_input_name in enumerate(v_input_names): + v_gpyopt_name = 'x[:,'+str(v_input_index)+']' + tmp_c = tmp_c.replace(v_input_name, v_gpyopt_name) + gpyopt_constraints.append([tmp_c, c.constraint_type]) + return gpyopt_constraints diff --git a/tests/test_strategies.py b/tests/test_strategies.py index 89769cad..bfb17fa7 100644 --- a/tests/test_strategies.py +++ b/tests/test_strategies.py @@ -419,3 +419,59 @@ def test_nm3D(maximize,x_start,constraint, plot=False): # (xbest[2] >= 0.851 and xbest[2] <= 0.853) and (fbest <= -3.85 and fbest >= -3.87) if plot: hartmann3D.plot(polygons=polygons_points) + + +@pytest.mark.parametrize('num_experiments', [1, 2, 4]) +@pytest.mark.parametrize('maximize', [True, False]) +@pytest.mark.parametrize('constraint', [True, False]) +def test_sobo(num_experiments, maximize, constraint, plot=False): + + hartmann3D = test_functions.Hartmann3D(maximize=maximize, constraints=constraint) + strategy = SOBO(domain=hartmann3D.domain) + + # Uncomment to start algorithm with pre-defined initial experiments + initial_exp = None + # Uncomment to create test case which results in reduction dimension and dimension recovery + #initial_exp = pd.DataFrame(data={'x_1': [0.1,0.1,0.4,0.3], 'x_2': [0.6,0.2,0.4,0.5], 'x_3': [1,1,1,0.3]}) # initial experimental points + #initial_exp = DataSet.from_df(initial_exp) + #initial_exp = hartmann3D.run_experiments(initial_exp) + + # run SOBO loop for fixed number of iteration + num_iter = 200//num_experiments # maximum number of iterations + max_stop = 80//num_experiments # allowed number of consecutive iterations w/o improvement + nstop = 0 + fbestold = float("inf") + + if initial_exp is not None: + next_experiments = initial_exp + else: + next_experiments = None + + param = None + for i in range(num_iter): + next_experiments, xbest, fbest, param = \ + strategy.suggest_experiments(num_experiments=num_experiments, prev_res=next_experiments, prev_param=param) + + # This is the part where experiments take place + next_experiments = hartmann3D.run_experiments(next_experiments) + + if fbest < fbestold: + fbestold = fbest + nstop = 0 + else: + nstop += 1 + if nstop >= max_stop: + print("No improvement in last " + str(max_stop) + " iterations.") + break + + print(next_experiments) # show next experiments + print("\n") + + xbest = np.around(xbest, decimals=3) + fbest = np.around(fbest, decimals=3) + print("Optimal setting: " + str(xbest) + " with outcome: " + str(fbest)) + # Extrema of test function without constraint: glob_min = -3.86 at (0.114,0.556,0.853) + assert (fbest <= -3.85 and fbest >= -3.87) + + if plot: + hartmann3D.plot()