Introduction (Himmelblau’s function)#
Let’s use bloptools to minimize Himmelblau’s function, which has four global minima:
[1]:
from bloptools.utils import prepare_re_env
%run -i $prepare_re_env.__file__ --db-type=temp
[2]:
import numpy as np
import matplotlib as mpl
from matplotlib import pyplot as plt
from bloptools.utils import functions
x1 = x2 = np.linspace(-6, 6, 1024)
X1, X2 = np.meshgrid(x1, x2)
F = functions.himmelblau(X1, X2)
plt.pcolormesh(x1, x2, F, norm=mpl.colors.LogNorm(vmin=1e-1, vmax=1e3), cmap="magma_r")
plt.colorbar()
plt.xlabel("x1")
plt.ylabel("x2")
[2]:
Text(0, 0.5, 'x2')
There are several things that our agent will need. The first ingredient is some degrees of freedom (these are always ophyd devices) which the agent will move around to different inputs within each DOF’s bounds (the second ingredient). We define these here:
[3]:
from bloptools.bayesian import DOF
dofs = [
DOF(name="x1", limits=(-6, 6)),
DOF(name="x2", limits=(-6, 6)),
]
We also need to give the agent something to do. We want our agent to look in the feedback for a variable called “himmelblau”, and try to minimize it.
[4]:
from bloptools.bayesian import Objective
objectives = [Objective(key="himmelblau", minimize=True)]
In our digestion function, we define our objective as a deterministic function of the inputs:
[5]:
def digestion(db, uid):
products = db[uid].table()
for index, entry in products.iterrows():
products.loc[index, "himmelblau"] = functions.himmelblau(entry.x1, entry.x2)
return products
We then combine these ingredients into an agent, giving it an instance of databroker so that it can see the output of the plans it runs.
[6]:
from bloptools.bayesian import Agent
agent = Agent(
dofs=dofs,
objectives=objectives,
digestion=digestion,
db=db,
)
Without any data, we can’t make any inferences about what the function looks like, and so we can’t use any non-trivial acquisition functions. Let’s start by quasi-randomly sampling the parameter space, and plotting our model of the function:
[7]:
RE(agent.learn("quasi-random", n=32))
agent.plot_objectives()
running iteration 1 / 1
Transient Scan ID: 1 Time: 2023-11-08 23:55:43
Persistent Unique Scan ID: '721def24-b2ab-47a1-a419-d4d3d5733394'
New stream: 'primary'
+-----------+------------+------------+------------+
| seq_num | time | x1 | x2 |
+-----------+------------+------------+------------+
| 1 | 23:55:43.6 | 0.403 | -1.335 |
| 2 | 23:55:43.6 | 1.633 | -2.829 |
| 3 | 23:55:43.6 | 2.888 | -3.278 |
| 4 | 23:55:43.7 | 1.119 | -4.772 |
| 5 | 23:55:43.7 | -0.550 | -5.070 |
| 6 | 23:55:43.7 | -1.593 | -3.565 |
| 7 | 23:55:43.7 | -3.160 | -5.762 |
| 8 | 23:55:43.7 | -4.929 | -4.256 |
| 9 | 23:55:43.7 | -5.549 | -1.664 |
| 10 | 23:55:43.7 | -4.639 | 0.345 |
| 11 | 23:55:43.7 | -4.319 | -0.158 |
| 12 | 23:55:43.7 | -2.939 | -2.519 |
| 13 | 23:55:43.7 | -0.982 | -1.014 |
| 14 | 23:55:43.7 | -1.975 | 0.826 |
| 15 | 23:55:43.7 | -3.408 | 1.851 |
| 16 | 23:55:43.7 | -4.125 | 4.443 |
| 17 | 23:55:43.7 | -5.894 | 5.949 |
| 18 | 23:55:43.7 | -2.505 | 4.882 |
| 19 | 23:55:43.7 | -1.462 | 3.377 |
| 20 | 23:55:43.7 | 0.154 | 3.090 |
| 21 | 23:55:43.7 | -0.018 | 2.331 |
| 22 | 23:55:43.7 | 1.313 | 2.642 |
| 23 | 23:55:43.7 | 1.924 | 4.585 |
| 24 | 23:55:43.7 | 4.587 | 5.581 |
| 25 | 23:55:43.7 | 3.544 | 4.087 |
| 26 | 23:55:43.7 | 3.976 | 2.184 |
| 27 | 23:55:43.7 | 2.543 | 1.147 |
| 28 | 23:55:43.7 | 3.013 | -0.503 |
| 29 | 23:55:43.7 | 5.933 | 0.691 |
| 30 | 23:55:43.7 | 4.970 | -1.997 |
| 31 | 23:55:43.7 | 5.499 | -3.899 |
| 32 | 23:55:43.7 | 4.456 | -5.393 |
+-----------+------------+------------+------------+
generator list_scan ['721def24'] (scan num: 1)
To decide which points to sample, the agent needs an acquisition function. The available acquisition function are here:
[8]:
agent.all_acq_funcs
EXPECTED IMPROVEMENT (identifiers: ['ei', 'expected_improvement'])
-> The expected value of max(f(x) - \nu, 0), where \nu is the current maximum.
MONTE CARLO EXPECTED IMPROVEMENT (identifiers: ['qei', 'monte_carlo_expected_improvement'])
-> The expected value of max(f(x) - \nu, 0), where \nu is the current maximum.
EXPECTED MEAN (identifiers: ['em', 'expected_mean'])
-> The expected value at each input.
MONTE CARLO EXPECTED MEAN (identifiers: ['qem', 'monte_carlo_expected_mean'])
-> The expected value at each input.
LOWER BOUND MAX VALUE ENTROPY (identifiers: ['lbmve', 'lbmes', 'gibbon', 'lower_bound_max_value_entropy'])
-> Max entropy search, basically
NOISY EXPECTED HYPERVOLUME IMPROVEMENT (identifiers: ['nehvi', 'noisy_expected_hypervolume_improvement'])
-> It's like a big box. How big is the box?
PROBABILITY OF IMPROVEMENT (identifiers: ['pi', 'probability_of_improvement'])
-> The probability that this input improves on the current maximum.
MONTE CARLO PROBABILITY OF IMPROVEMENT (identifiers: ['qpi', 'monte_carlo_probability_of_improvement'])
-> The probability that this input improves on the current maximum.
RANDOM (identifiers: ['r', 'random'])
-> Uniformly-sampled random points.
QUASI-RANDOM (identifiers: ['qr', 'quasi-random'])
-> Sobol-sampled quasi-random points.
GRID SCAN (identifiers: ['g', 'gs', 'grid'])
-> A grid scan over the parameters.
UPPER CONFIDENCE BOUND (identifiers: ['ucb', 'upper_confidence_bound'])
-> The expected value, plus some multiple of the uncertainty (typically \mu + 2\sigma).
MONTE CARLO UPPER CONFIDENCE BOUND (identifiers: ['qucb', 'monte_carlo_upper_confidence_bound'])
-> The expected value, plus some multiple of the uncertainty (typically \mu + 2\sigma).
Now we can start to learn intelligently. Using the shorthand acquisition functions shown above, we can see the output of a few different ones:
[9]:
agent.plot_acquisition(acq_func="qei")
To decide where to go, the agent will find the inputs that maximize a given acquisition function:
[10]:
agent.ask("qei", n=1)
[10]:
(array([[2.79727951, 2.50029115]]),
{'name': 'monte_carlo_expected_improvement',
'args': {},
'read_only_values': array([], shape=(1, 0), dtype=float64),
'duration': 0.3987107699999797})
We can also ask the agent for multiple points to sample and it will jointly maximize the acquisition function over all sets of inputs, and find the most efficient route between them:
[11]:
X, _ = agent.ask("qei", n=8, route=True)
agent.plot_acquisition(acq_func="qei")
plt.scatter(*X.T, marker="d", facecolor="w", edgecolor="k")
plt.plot(
*X.T,
color="r",
)
[11]:
[<matplotlib.lines.Line2D at 0x7f79be133bb0>]
All of this is automated inside the learn method, which will find a point (or points) to sample, sample them, and retrain the model and its hyperparameters with the new data. To do 4 learning iterations of 8 points each, we can run
[12]:
RE(agent.learn("qei", n=4, iterations=8))
running iteration 1 / 8
Transient Scan ID: 2 Time: 2023-11-08 23:55:51
Persistent Unique Scan ID: 'c9f8b5fb-dc26-4bde-b7a4-5f0e264642e3'
New stream: 'primary'
+-----------+------------+------------+------------+
| seq_num | time | x1 | x2 |
+-----------+------------+------------+------------+
| 1 | 23:55:51.7 | 3.314 | -1.898 |
| 2 | 23:55:51.7 | 2.801 | 2.501 |
| 3 | 23:55:51.7 | -3.035 | 3.223 |
| 4 | 23:55:51.8 | -6.000 | -6.000 |
+-----------+------------+------------+------------+
generator list_scan ['c9f8b5fb'] (scan num: 2)
running iteration 2 / 8
Transient Scan ID: 3 Time: 2023-11-08 23:55:53
Persistent Unique Scan ID: 'ffd0c9b7-a0b4-448c-b5df-78374a7f32d8'
New stream: 'primary'
+-----------+------------+------------+------------+
| seq_num | time | x1 | x2 |
+-----------+------------+------------+------------+
| 1 | 23:55:53.6 | -6.000 | -6.000 |
| 2 | 23:55:53.6 | -3.033 | 3.221 |
| 3 | 23:55:53.6 | 2.797 | 2.501 |
| 4 | 23:55:53.6 | 3.301 | -1.904 |
+-----------+------------+------------+------------+
generator list_scan ['ffd0c9b7'] (scan num: 3)
running iteration 3 / 8
Transient Scan ID: 4 Time: 2023-11-08 23:55:55
Persistent Unique Scan ID: '34a79422-4c28-47e5-afef-7b6daa08ce0d'
New stream: 'primary'
+-----------+------------+------------+------------+
| seq_num | time | x1 | x2 |
+-----------+------------+------------+------------+
| 1 | 23:55:55.1 | 3.318 | -1.906 |
| 2 | 23:55:55.1 | 2.797 | 2.501 |
| 3 | 23:55:55.1 | -3.034 | 3.222 |
| 4 | 23:55:55.1 | -6.000 | -6.000 |
+-----------+------------+------------+------------+
generator list_scan ['34a79422'] (scan num: 4)
running iteration 4 / 8
Transient Scan ID: 5 Time: 2023-11-08 23:55:56
Persistent Unique Scan ID: '18736dd6-1a4b-4f03-9bf3-65efa6735a4e'
New stream: 'primary'
+-----------+------------+------------+------------+
| seq_num | time | x1 | x2 |
+-----------+------------+------------+------------+
| 1 | 23:55:56.8 | -6.000 | -6.000 |
| 2 | 23:55:56.8 | -3.029 | 3.214 |
| 3 | 23:55:56.8 | 2.798 | 2.500 |
| 4 | 23:55:56.8 | 3.299 | -1.893 |
+-----------+------------+------------+------------+
generator list_scan ['18736dd6'] (scan num: 5)
running iteration 5 / 8
Transient Scan ID: 6 Time: 2023-11-08 23:55:58
Persistent Unique Scan ID: 'a390ecdb-f9b1-4c38-8862-443eb7826cec'
New stream: 'primary'
+-----------+------------+------------+------------+
| seq_num | time | x1 | x2 |
+-----------+------------+------------+------------+
| 1 | 23:55:58.6 | 3.299 | -1.902 |
| 2 | 23:55:58.6 | 2.797 | 2.500 |
| 3 | 23:55:58.6 | -3.033 | 3.221 |
| 4 | 23:55:58.6 | -6.000 | -6.000 |
+-----------+------------+------------+------------+
generator list_scan ['a390ecdb'] (scan num: 6)
running iteration 6 / 8
Transient Scan ID: 7 Time: 2023-11-08 23:56:00
Persistent Unique Scan ID: 'fbb45a85-add0-4087-bc9d-7e1e84de5e5c'
New stream: 'primary'
+-----------+------------+------------+------------+
| seq_num | time | x1 | x2 |
+-----------+------------+------------+------------+
| 1 | 23:56:00.3 | -6.000 | -6.000 |
| 2 | 23:56:00.3 | -3.032 | 3.221 |
| 3 | 23:56:00.3 | 2.797 | 2.501 |
| 4 | 23:56:00.3 | 3.307 | -1.896 |
+-----------+------------+------------+------------+
generator list_scan ['fbb45a85'] (scan num: 7)
running iteration 7 / 8
Transient Scan ID: 8 Time: 2023-11-08 23:56:02
Persistent Unique Scan ID: '72b0876c-afe3-423e-842d-7179f4f730fd'
New stream: 'primary'
+-----------+------------+------------+------------+
| seq_num | time | x1 | x2 |
+-----------+------------+------------+------------+
| 1 | 23:56:02.0 | 3.306 | -1.907 |
| 2 | 23:56:02.0 | 2.797 | 2.500 |
| 3 | 23:56:02.0 | -3.034 | 3.221 |
| 4 | 23:56:02.0 | -6.000 | -6.000 |
+-----------+------------+------------+------------+
generator list_scan ['72b0876c'] (scan num: 8)
running iteration 8 / 8
Transient Scan ID: 9 Time: 2023-11-08 23:56:04
Persistent Unique Scan ID: '7c7a8d53-4a74-4975-b526-4108d128b5eb'
New stream: 'primary'
+-----------+------------+------------+------------+
| seq_num | time | x1 | x2 |
+-----------+------------+------------+------------+
| 1 | 23:56:04.0 | -6.000 | -6.000 |
| 2 | 23:56:04.0 | -3.039 | 3.225 |
| 3 | 23:56:04.0 | 2.795 | 2.503 |
| 4 | 23:56:04.0 | 3.300 | -1.905 |
+-----------+------------+------------+------------+
generator list_scan ['7c7a8d53'] (scan num: 9)
[12]:
('c9f8b5fb-dc26-4bde-b7a4-5f0e264642e3',
'ffd0c9b7-a0b4-448c-b5df-78374a7f32d8',
'34a79422-4c28-47e5-afef-7b6daa08ce0d',
'18736dd6-1a4b-4f03-9bf3-65efa6735a4e',
'a390ecdb-f9b1-4c38-8862-443eb7826cec',
'fbb45a85-add0-4087-bc9d-7e1e84de5e5c',
'72b0876c-afe3-423e-842d-7179f4f730fd',
'7c7a8d53-4a74-4975-b526-4108d128b5eb')
Our agent has found all the global minima of Himmelblau’s function using Bayesian optimization, and we can ask it for the best point:
[13]:
agent.plot_objectives()
print(agent.best)
x1 -3.03913
x2 3.224808
himmelblau 2.26465
time 2023-11-08 23:56:04.029840231
acq_func monte_carlo_expected_improvement
himmelblau_fitness -2.26465
Name: 33, dtype: object
