From e8301c26b019665eea64ad10a33c9eefbba1e607 Mon Sep 17 00:00:00 2001 From: Thomas Morris Date: Thu, 18 Apr 2024 11:13:28 -0700 Subject: [PATCH] added pareto example to the docs --- .../{himmelblau.ipynb => introduction.ipynb} | 0 docs/source/tutorials/pareto-fronts.ipynb | 159 ++++++++++++++++++ src/blop/agent.py | 11 +- src/blop/bayesian/transforms.py | 79 --------- src/blop/plotting.py | 2 +- src/blop/tests/test_agent.py | 3 +- src/blop/utils/functions.py | 6 +- 7 files changed, 175 insertions(+), 85 deletions(-) rename docs/source/tutorials/{himmelblau.ipynb => introduction.ipynb} (100%) create mode 100644 docs/source/tutorials/pareto-fronts.ipynb delete mode 100644 src/blop/bayesian/transforms.py diff --git a/docs/source/tutorials/himmelblau.ipynb b/docs/source/tutorials/introduction.ipynb similarity index 100% rename from docs/source/tutorials/himmelblau.ipynb rename to docs/source/tutorials/introduction.ipynb diff --git a/docs/source/tutorials/pareto-fronts.ipynb b/docs/source/tutorials/pareto-fronts.ipynb new file mode 100644 index 0000000..4b8a13c --- /dev/null +++ b/docs/source/tutorials/pareto-fronts.ipynb @@ -0,0 +1,159 @@ +{ + "cells": [ + { + "attachments": {}, + "cell_type": "markdown", + "id": "e7b5e13a-c059-441d-8d4f-fff080d52054", + "metadata": {}, + "source": [ + "# Multiobjective optimization with Pareto front mapping\n", + "\n", + "One way to do multiobjective optimization is with Pareto optimization, which explores the set of Pareto-efficient points. A point is Pareto-efficient if there are no other valid points that are better at every objective: it shows the \"trade-off\" between several objectives. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cac0177b-576c-4f01-b306-2e9e8544a05c", + "metadata": {}, + "outputs": [], + "source": [ + "from blop.utils import prepare_re_env\n", + "\n", + "%run -i $prepare_re_env.__file__ --db-type=temp" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "120812d8-5de2-4efa-8fc6-2d8e4cd0693e", + "metadata": {}, + "outputs": [], + "source": [ + "from blop import DOF, Objective, Agent\n", + "from blop.utils import functions\n", + "from blop.dofs import BrownianMotion\n", + "\n", + "import numpy as np\n", + "\n", + "\n", + "def digestion(db, uid):\n", + " products = db[uid].table()\n", + "\n", + " for index, entry in products.iterrows():\n", + " x1, x2 = entry.x1, entry.x2\n", + "\n", + " products.loc[index, \"f1\"] = (x1 - 2) ** 2 + (x2 - 1) + 2\n", + " products.loc[index, \"f2\"] = 9 * x1 - (x2 - 1) + 2\n", + " products.loc[index, \"c1\"] = x1**2 + x2**2\n", + " products.loc[index, \"c2\"] = x1 - 3 * x2 + 10\n", + "\n", + " return products\n", + "\n", + "\n", + "dofs = [\n", + " DOF(name=\"x1\", search_domain=(-20, 20)),\n", + " DOF(name=\"x2\", search_domain=(-20, 20)),\n", + "]\n", + "\n", + "\n", + "objectives = [\n", + " Objective(name=\"f1\", target=\"min\"),\n", + " Objective(name=\"f2\", target=\"min\"),\n", + " Objective(name=\"c1\", target=(-np.inf, 225)),\n", + " Objective(name=\"c2\", target=(-np.inf, 0)),\n", + "]\n", + "\n", + "agent = Agent(\n", + " dofs=dofs,\n", + " objectives=objectives,\n", + " digestion=digestion,\n", + " db=db,\n", + ")\n", + "\n", + "(uid,) = RE(agent.learn(\"qr\", n=64))" + ] + }, + { + "cell_type": "markdown", + "id": "d81c6af2-0a4c-4d31-9b7c-cc3065550b98", + "metadata": {}, + "source": [ + "We can plot our fitness and constraint objectives to see their models:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c6c08f53-e96b-4987-ba3d-a93c84468b0b", + "metadata": {}, + "outputs": [], + "source": [ + "agent.plot_objectives()" + ] + }, + { + "cell_type": "markdown", + "id": "48b976e6-048a-4e16-9a1c-500203a2e195", + "metadata": {}, + "source": [ + "We can plot the Pareto front (the set of all Pareto-efficient points), which shows the trade-off between the two fitnesses. The points in blue comprise the Pareto front, while the points in red are either not Pareto efficient or are invalidated by one of the constraints." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "990a877e-f533-419c-bf5d-569ad7e72c6b", + "metadata": {}, + "outputs": [], + "source": [ + "agent.plot_pareto_front()" + ] + }, + { + "cell_type": "markdown", + "id": "29d7f4fb-8f25-4b57-982f-28737fad2a7c", + "metadata": {}, + "source": [ + "We can explore the Pareto front by choosing a random point on the Pareto front and computing the expected improvement in the hypervolume of all fitness objectives with respect to that point (called the \"reference point\"). All this is done automatically with the `qnehvi` acquisition function:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b49e6233-b228-43a3-9d8a-722a82e93443", + "metadata": {}, + "outputs": [], + "source": [ + "# this is broken now but is fixed in the next PR\n", + "# RE(agent.learn(\"qnehvi\", n=4))" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.0" + }, + "vscode": { + "interpreter": { + "hash": "b0fa6594d8f4cbf19f97940f81e996739fb7646882a419484c72d19e05852a7e" + } + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/src/blop/agent.py b/src/blop/agent.py index a8dd2d9..7d4fe37 100644 --- a/src/blop/agent.py +++ b/src/blop/agent.py @@ -601,8 +601,12 @@ def best_f(self): @property def pareto_front_mask(self): - # a point is on the Pareto front if it is Pareto dominant - # a point is Pareto dominant if it is there is no other point that is better at every objective + """ + A mask for all data points that is true when the point is on the Pareto front. + A point is on the Pareto front if it is Pareto dominant + A point is Pareto dominant if it is there is no other point that is better at every objective + Points that violate any constraint are excluded. + """ y = self.train_targets()[:, self.objectives.type == "fitness"] in_pareto_front = ~(y.unsqueeze(1) > y.unsqueeze(0)).all(axis=-1).any(axis=0) all_constraints_satisfied = self.evaluated_constraints.all(axis=-1) @@ -610,6 +614,9 @@ def pareto_front_mask(self): @property def pareto_front(self): + """ + A subset of the data table containing only points on the Pareto front. + """ return self.table.loc[self.pareto_front_mask.numpy()] @property diff --git a/src/blop/bayesian/transforms.py b/src/blop/bayesian/transforms.py deleted file mode 100644 index 0999784..0000000 --- a/src/blop/bayesian/transforms.py +++ /dev/null @@ -1,79 +0,0 @@ -from typing import Union - -import botorch -from botorch.acquisition.objective import PosteriorTransform -from botorch.posteriors.gpytorch import GPyTorchPosterior -from botorch.posteriors.posterior_list import PosteriorList -from torch import Tensor - -# This file is deprecated for now, but we may need it later - - -def targeting_transform(y, target): - if target == "max": - return y - if target == "min": - return -y - elif not isinstance(target, tuple): - return -(y - target).abs() - else: - return -((y - 0.5 * (target[1] + target[0])).abs() - 0.5 * (target[1] - target[0])).clamp(min=0) - - -class TargetingPosteriorTransform(PosteriorTransform): - r"""An affine posterior transform for scalarizing multi-output posteriors.""" - - scalarize: bool = True - - def __init__(self, weights: Tensor, targets: Tensor) -> None: - r""" - Args: - weights: A one-dimensional tensor with `m` elements representing the - linear weights on the outputs. - offset: An offset to be added to posterior mean. - """ - super().__init__() - self.targets = targets - self.register_buffer("weights", weights) - - def sample_transform(self, y): - for i, target in enumerate(self.targets): - y[..., i] = targeting_transform(y[..., i], target) - return y @ self.weights.unsqueeze(-1) - - def mean_transform(self, mean, var): - for i, target in enumerate(self.targets): - mean[..., i] = targeting_transform(mean[..., i], target) - return mean @ self.weights.unsqueeze(-1) - - def variance_transform(self, mean, var): - return var @ self.weights.unsqueeze(-1) - - def evaluate(self, Y: Tensor) -> Tensor: - r"""Evaluate the transform on a set of outcomes. - - Args: - Y: A `batch_shape x q x m`-dim tensor of outcomes. - - Returns: - A `batch_shape x q`-dim tensor of transformed outcomes. - """ - return self.sample_transform(Y) - - def forward(self, posterior: Union[GPyTorchPosterior, PosteriorList]) -> GPyTorchPosterior: - r"""Compute the posterior of the affine transformation. - - Args: - posterior: A posterior with the same number of outputs as the - elements in `self.weights`. - - Returns: - A single-output posterior. - """ - - return botorch.posteriors.transformed.TransformedPosterior( - posterior, - sample_transform=self.sample_transform, - mean_transform=self.mean_transform, - variance_transform=self.variance_transform, - ) diff --git a/src/blop/plotting.py b/src/blop/plotting.py index 80659aa..2c9e66b 100644 --- a/src/blop/plotting.py +++ b/src/blop/plotting.py @@ -185,7 +185,7 @@ def _plot_objs_many_dofs(agent, axes=(0, 1), shading="nearest", cmap=DEFAULT_COL ) val_cbar = agent.obj_fig.colorbar(val_ax, ax=agent.obj_axes[obj_index, 0], location="bottom", aspect=32, shrink=0.8) - val_cbar.set_label(f"{obj.units if obj.units else ''}") + val_cbar.set_label(f"{obj.units or ''}") if obj.is_fitness: _ = agent.obj_fig.colorbar(fitness_ax, ax=agent.obj_axes[obj_index, 1], location="bottom", aspect=32, shrink=0.8) diff --git a/src/blop/tests/test_agent.py b/src/blop/tests/test_agent.py index 7b6f117..a04b432 100644 --- a/src/blop/tests/test_agent.py +++ b/src/blop/tests/test_agent.py @@ -4,7 +4,8 @@ def test_agent(agent, RE): RE(agent.learn("qr", n=4)) - agent.best + assert [dof.name in agent.best for dof in agent.dofs] + assert [obj.name in agent.best for obj in agent.objectives] def test_forget(agent, RE): diff --git a/src/blop/utils/functions.py b/src/blop/utils/functions.py index 745e316..74f6f42 100644 --- a/src/blop/utils/functions.py +++ b/src/blop/utils/functions.py @@ -3,8 +3,10 @@ def approximate_erf(x): - # we want to compute erf(x) to compute the definite integral of the Gaussian PDF - # we instead use an approximation with tanh(x) which is faster and better-conditioned near +/- infinity + """ + An approximation of erf(x), to compute the definite integral of the Gaussian PDF + This is faster and better-conditioned near +/- infinity + """ return torch.tanh(1.20278247 * x)