added pareto example to the docs

docs/source/tutorials.rst

@@ -4,6 +4,7 @@ Tutorials
 .. toctree::
    :maxdepth: 2
 
-   tutorials/himmelblau.ipynb
+   tutorials/introduction.ipynb
    tutorials/hyperparameters.ipynb
+   tutorials/pareto-fronts.ipynb
    tutorials/passive-dofs.ipynb

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
+}

src/blop/agent.py

@@ -601,15 +601,22 @@ 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)
         return in_pareto_front & all_constraints_satisfied
 
     @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

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)

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):

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)