|
| 1 | +""" |
| 2 | +============================================== |
| 3 | +Use different base estimators for optimization |
| 4 | +============================================== |
| 5 | +
|
| 6 | +Sigurd Carlen, September 2019. |
| 7 | +Reformatted by Holger Nahrstaedt 2020 |
| 8 | +
|
| 9 | +.. currentmodule:: skopt |
| 10 | +
|
| 11 | +
|
| 12 | +To use different base_estimator or create a regressor with different parameters, |
| 13 | +we can create a regressor object and set it as kernel. |
| 14 | +
|
| 15 | +""" |
| 16 | +print(__doc__) |
| 17 | + |
| 18 | +import numpy as np |
| 19 | +np.random.seed(1234) |
| 20 | +import matplotlib.pyplot as plt |
| 21 | + |
| 22 | + |
| 23 | +############################################################################# |
| 24 | +# Toy example |
| 25 | +# ----------- |
| 26 | +# |
| 27 | +# Let assume the following noisy function :math:`f`: |
| 28 | + |
| 29 | +noise_level = 0.1 |
| 30 | + |
| 31 | +# Our 1D toy problem, this is the function we are trying to |
| 32 | +# minimize |
| 33 | + |
| 34 | +def objective(x, noise_level=noise_level): |
| 35 | + return np.sin(5 * x[0]) * (1 - np.tanh(x[0] ** 2))\ |
| 36 | + + np.random.randn() * noise_level |
| 37 | + |
| 38 | +############################################################################# |
| 39 | + |
| 40 | +from skopt import Optimizer |
| 41 | +opt_gp = Optimizer([(-2.0, 2.0)], base_estimator="GP", n_initial_points=5, |
| 42 | + acq_optimizer="sampling", random_state=42) |
| 43 | + |
| 44 | +############################################################################# |
| 45 | + |
| 46 | +x = np.linspace(-2, 2, 400).reshape(-1, 1) |
| 47 | +fx = np.array([objective(x_i, noise_level=0.0) for x_i in x]) |
| 48 | + |
| 49 | +############################################################################# |
| 50 | + |
| 51 | +from skopt.acquisition import gaussian_ei |
| 52 | + |
| 53 | +def plot_optimizer(res, next_x, x, fx, n_iter, max_iters=5): |
| 54 | + x_gp = res.space.transform(x.tolist()) |
| 55 | + gp = res.models[-1] |
| 56 | + curr_x_iters = res.x_iters |
| 57 | + curr_func_vals = res.func_vals |
| 58 | + |
| 59 | + # Plot true function. |
| 60 | + ax = plt.subplot(max_iters, 2, 2 * n_iter + 1) |
| 61 | + plt.plot(x, fx, "r--", label="True (unknown)") |
| 62 | + plt.fill(np.concatenate([x, x[::-1]]), |
| 63 | + np.concatenate([fx - 1.9600 * noise_level, |
| 64 | + fx[::-1] + 1.9600 * noise_level]), |
| 65 | + alpha=.2, fc="r", ec="None") |
| 66 | + if n_iter < max_iters - 1: |
| 67 | + ax.get_xaxis().set_ticklabels([]) |
| 68 | + # Plot GP(x) + contours |
| 69 | + y_pred, sigma = gp.predict(x_gp, return_std=True) |
| 70 | + plt.plot(x, y_pred, "g--", label=r"$\mu_{GP}(x)$") |
| 71 | + plt.fill(np.concatenate([x, x[::-1]]), |
| 72 | + np.concatenate([y_pred - 1.9600 * sigma, |
| 73 | + (y_pred + 1.9600 * sigma)[::-1]]), |
| 74 | + alpha=.2, fc="g", ec="None") |
| 75 | + |
| 76 | + # Plot sampled points |
| 77 | + plt.plot(curr_x_iters, curr_func_vals, |
| 78 | + "r.", markersize=8, label="Observations") |
| 79 | + plt.title(r"x* = %.4f, f(x*) = %.4f" % (res.x[0], res.fun)) |
| 80 | + # Adjust plot layout |
| 81 | + plt.grid() |
| 82 | + |
| 83 | + if n_iter == 0: |
| 84 | + plt.legend(loc="best", prop={'size': 6}, numpoints=1) |
| 85 | + |
| 86 | + if n_iter != 4: |
| 87 | + plt.tick_params(axis='x', which='both', bottom='off', |
| 88 | + top='off', labelbottom='off') |
| 89 | + |
| 90 | + # Plot EI(x) |
| 91 | + ax = plt.subplot(max_iters, 2, 2 * n_iter + 2) |
| 92 | + acq = gaussian_ei(x_gp, gp, y_opt=np.min(curr_func_vals)) |
| 93 | + plt.plot(x, acq, "b", label="EI(x)") |
| 94 | + plt.fill_between(x.ravel(), -2.0, acq.ravel(), alpha=0.3, color='blue') |
| 95 | + |
| 96 | + if n_iter < max_iters - 1: |
| 97 | + ax.get_xaxis().set_ticklabels([]) |
| 98 | + |
| 99 | + next_acq = gaussian_ei(res.space.transform([next_x]), gp, |
| 100 | + y_opt=np.min(curr_func_vals)) |
| 101 | + plt.plot(next_x, next_acq, "bo", markersize=6, label="Next query point") |
| 102 | + |
| 103 | + # Adjust plot layout |
| 104 | + plt.ylim(0, 0.07) |
| 105 | + plt.grid() |
| 106 | + if n_iter == 0: |
| 107 | + plt.legend(loc="best", prop={'size': 6}, numpoints=1) |
| 108 | + |
| 109 | + if n_iter != 4: |
| 110 | + plt.tick_params(axis='x', which='both', bottom='off', |
| 111 | + top='off', labelbottom='off') |
| 112 | + |
| 113 | +############################################################################# |
| 114 | +# GP kernel |
| 115 | +# --------- |
| 116 | + |
| 117 | +fig = plt.figure() |
| 118 | +fig.suptitle("Standard GP kernel") |
| 119 | +for i in range(10): |
| 120 | + next_x = opt_gp.ask() |
| 121 | + f_val = objective(next_x) |
| 122 | + res = opt_gp.tell(next_x, f_val) |
| 123 | + if i >= 5: |
| 124 | + plot_optimizer(res, opt_gp._next_x, x, fx, n_iter=i-5, max_iters=5) |
| 125 | +plt.tight_layout(rect=[0, 0.03, 1, 0.95]) |
| 126 | +plt.plot() |
| 127 | + |
| 128 | +############################################################################# |
| 129 | +# Test different kernels |
| 130 | +# ---------------------- |
| 131 | + |
| 132 | +from skopt.learning import GaussianProcessRegressor |
| 133 | +from skopt.learning.gaussian_process.kernels import ConstantKernel, Matern |
| 134 | +# Gaussian process with Matérn kernel as surrogate model |
| 135 | + |
| 136 | +from sklearn.gaussian_process.kernels import (RBF, Matern, RationalQuadratic, |
| 137 | + ExpSineSquared, DotProduct, |
| 138 | + ConstantKernel) |
| 139 | + |
| 140 | + |
| 141 | +kernels = [1.0 * RBF(length_scale=1.0, length_scale_bounds=(1e-1, 10.0)), |
| 142 | + 1.0 * RationalQuadratic(length_scale=1.0, alpha=0.1), |
| 143 | + 1.0 * ExpSineSquared(length_scale=1.0, periodicity=3.0, |
| 144 | + length_scale_bounds=(0.1, 10.0), |
| 145 | + periodicity_bounds=(1.0, 10.0)), |
| 146 | + ConstantKernel(0.1, (0.01, 10.0)) |
| 147 | + * (DotProduct(sigma_0=1.0, sigma_0_bounds=(0.1, 10.0)) ** 2), |
| 148 | + 1.0 * Matern(length_scale=1.0, length_scale_bounds=(1e-1, 10.0), |
| 149 | + nu=2.5)] |
| 150 | + |
| 151 | +############################################################################# |
| 152 | + |
| 153 | +for kernel in kernels: |
| 154 | + gpr = GaussianProcessRegressor(kernel=kernel, alpha=noise_level ** 2, |
| 155 | + normalize_y=True, noise="gaussian", |
| 156 | + n_restarts_optimizer=2 |
| 157 | + ) |
| 158 | + opt = Optimizer([(-2.0, 2.0)], base_estimator=gpr, n_initial_points=5, |
| 159 | + acq_optimizer="sampling", random_state=42) |
| 160 | + fig = plt.figure() |
| 161 | + fig.suptitle(repr(kernel)) |
| 162 | + for i in range(10): |
| 163 | + next_x = opt.ask() |
| 164 | + f_val = objective(next_x) |
| 165 | + res = opt.tell(next_x, f_val) |
| 166 | + if i >= 5: |
| 167 | + plot_optimizer(res, opt._next_x, x, fx, n_iter=i - 5, max_iters=5) |
| 168 | + plt.tight_layout(rect=[0, 0.03, 1, 0.95]) |
| 169 | + plt.show() |
0 commit comments