diff --git a/figures.ipynb b/figures.ipynb
index 93d0a7807642961037a35e4b725558c75fd03bee..0ec869533be667fa64ed9e54d79ac250f934a5bb 100644
--- a/figures.ipynb
+++ b/figures.ipynb
@@ -39,7 +39,16 @@
"}\n",
"\n",
"plt.rcParams.update(params)\n",
- "plt.rc(\"text.latex\", preamble=[r\"\\usepackage{xfrac}\", r\"\\usepackage{siunitx}\"])"
+ "plt.rc(\"text.latex\", preamble=[r\"\\usepackage{xfrac}\", r\"\\usepackage{siunitx}\"])\n",
+ "\n",
+ "import adaptive\n",
+ "from scipy import interpolate\n",
+ "import functools\n",
+ "import itertools\n",
+ "import adaptive\n",
+ "import holoviews.plotting.mpl\n",
+ "import time\n",
+ "import matplotlib.tri as mtri"
]
},
{
@@ -113,9 +122,6 @@
"metadata": {},
"outputs": [],
"source": [
- "import adaptive\n",
- "\n",
- "\n",
"def f(x, offset=0.123):\n",
" a = 0.02\n",
" return x + a ** 2 / (a ** 2 + (x - offset) ** 2)\n",
@@ -129,33 +135,35 @@
" return np.sin(100 * x) * np.exp(-x ** 2 / 0.1 ** 2)\n",
"\n",
"\n",
- "funcs = [\n",
+ "funcs_1D = [\n",
" dict(function=f, bounds=(-1, 1), title=\"peak\"),\n",
" dict(function=g, bounds=(-1, 1), title=\"tanh\"),\n",
" dict(function=h, bounds=(-0.3, 0.3), title=\"wave packet\"),\n",
"]\n",
- "fig, axs = plt.subplots(2, len(funcs), figsize=(fig_width, 1.5 * fig_height))\n",
+ "fig, axs = plt.subplots(2, len(funcs_1D), figsize=(fig_width, 1.5 * fig_height))\n",
"n_points = 50\n",
"for i, ax in enumerate(axs.T.flatten()):\n",
" ax.xaxis.set_ticks([])\n",
" ax.yaxis.set_ticks([])\n",
- " if i % 2 == 0:\n",
- " d = funcs[i // 2]\n",
- " # homogeneous\n",
- " xs = np.linspace(*d[\"bounds\"], n_points)\n",
- " ys = d[\"function\"](xs)\n",
+ " kind = \"homogeneous\" if i % 2 == 0 else \"adaptive\"\n",
+ " index = i // 2 if kind == \"homogeneous\" else (i - 1) // 2\n",
+ " d = funcs_1D[index]\n",
+ " bounds = d[\"bounds\"]\n",
+ " f = d[\"function\"]\n",
+ "\n",
+ " if kind == \"homogeneous\":\n",
+ " xs = np.linspace(*bounds, n_points)\n",
+ " ys = f(xs)\n",
" ax.set_title(rf\"\\textrm{{{d['title']}}}\")\n",
- " else:\n",
- " d = funcs[(i - 1) // 2]\n",
+ " elif kind == \"adaptive\":\n",
" loss = adaptive.learner.learner1D.curvature_loss_function()\n",
" learner = adaptive.Learner1D(\n",
- " d[\"function\"], bounds=d[\"bounds\"], loss_per_interval=loss\n",
+ " f, bounds=bounds, loss_per_interval=loss\n",
" )\n",
" adaptive.runner.simple(learner, goal=lambda l: l.npoints >= n_points)\n",
- " # adaptive\n",
" xs, ys = zip(*sorted(learner.data.items()))\n",
- " xs_dense = np.linspace(*d[\"bounds\"], 1000)\n",
- " ax.plot(xs_dense, d[\"function\"](xs_dense), c=\"red\", alpha=0.3, lw=0.5)\n",
+ " xs_dense = np.linspace(*bounds, 1000)\n",
+ " ax.plot(xs_dense, f(xs_dense), c=\"red\", alpha=0.3, lw=0.5)\n",
" ax.scatter(xs, ys, s=0.5, c=\"k\")\n",
"\n",
"axs[0][0].set_ylabel(r\"$\\textrm{homogeneous}$\")\n",
@@ -176,13 +184,6 @@
"metadata": {},
"outputs": [],
"source": [
- "from scipy import interpolate\n",
- "import functools\n",
- "import itertools\n",
- "import adaptive\n",
- "import holoviews.plotting.mpl\n",
- "import matplotlib.tri as mtri\n",
- "\n",
"\n",
"def f(xy, offset=0.123):\n",
" a = 0.2\n",
@@ -310,134 +311,23 @@
"metadata": {},
"outputs": [],
"source": [
- "from scipy import interpolate\n",
- "import functools\n",
- "import itertools\n",
- "import adaptive\n",
- "import holoviews.plotting.mpl\n",
- "import matplotlib.tri as mtri\n",
- "\n",
- "\n",
- "def f(xy, offset=0.123):\n",
- " a = 0.2\n",
- " x, y = xy\n",
- " return x + np.exp(-(x ** 2 + y ** 2 - 0.75 ** 2) ** 2 / a ** 4)\n",
- "\n",
- "\n",
- "@functools.lru_cache()\n",
- "def g_setup(fname):\n",
- " data = adaptive.utils.load(fname)\n",
- " points = np.array(list(data.keys()))\n",
- " values = np.array(list(data.values()), dtype=float)\n",
- " bounds = [\n",
- " (points[:, 0].min(), points[:, 0].max()),\n",
- " (points[:, 1].min(), points[:, 1].max()),\n",
- " ]\n",
- " ll, ur = np.reshape(bounds, (2, 2)).T\n",
- " inds = np.all(np.logical_and(ll <= points, points <= ur), axis=1)\n",
- " points, values = points[inds], values[inds].reshape(-1, 1)\n",
- " return interpolate.LinearNDInterpolator(points, values), bounds\n",
- "\n",
- "\n",
- "def g(xy, fname):\n",
- " ip, _ = g_setup(fname)\n",
- " return np.round(ip(xy))\n",
- "\n",
- "\n",
- "def density(x, eps=0):\n",
- " e = [0.8, 0.2]\n",
- " delta = [0.5, 0.5, 0.5]\n",
- " c = 3\n",
- " omega = [0.02, 0.05]\n",
- "\n",
- " H = np.array(\n",
- " [\n",
- " [e[0] + 1j * omega[0], delta[0], delta[1]],\n",
- " [delta[0], e[1] + c * x + 1j * omega[1], delta[1]],\n",
- " [delta[1], delta[2], e[1] - c * x + 1j * omega[1]],\n",
- " ]\n",
- " )\n",
- " H += np.eye(3) * eps\n",
- " return np.trace(np.linalg.inv(H)).imag\n",
- "\n",
- "\n",
- "def h(xy):\n",
- " x, y = xy\n",
- " return density(x, y) + y\n",
- "\n",
- "\n",
- "funcs = [\n",
- " dict(function=f, bounds=[(-1, 1), (-1, 1)], npoints=33),\n",
- " dict(\n",
- " function=g,\n",
- " bounds=g_setup(\"phase_diagram.pickle\")[1],\n",
- " npoints=100,\n",
- " fname=\"phase_diagram.pickle\",\n",
- " ),\n",
- " dict(function=h, bounds=[(-1, 1), (-3, 3)], npoints=50),\n",
- "]\n",
- "fig, axs = plt.subplots(len(funcs), 2, figsize=(fig_width, 2 * fig_height))\n",
- "\n",
- "plt.subplots_adjust(hspace=0.1, wspace=0.1)\n",
- "\n",
- "with_tri = False\n",
- "\n",
- "for i, ax in enumerate(axs.flatten()):\n",
- " label = \"abcdef\"[i]\n",
- " ax.text(\n",
- " -0.03,\n",
- " 0.98,\n",
- " f\"$\\mathrm{{({label})}}$\",\n",
- " transform=ax.transAxes,\n",
- " horizontalalignment=\"right\",\n",
- " verticalalignment=\"top\",\n",
- " )\n",
- " ax.xaxis.set_ticks([])\n",
- " ax.yaxis.set_ticks([])\n",
- " kind = \"homogeneous\" if i % 2 == 0 else \"adaptive\"\n",
- " d = funcs[i // 2] if kind == \"homogeneous\" else funcs[(i - 1) // 2]\n",
- " bounds = d[\"bounds\"]\n",
- " npoints = d[\"npoints\"]\n",
- " f = d[\"function\"]\n",
- " fname = d.get(\"fname\")\n",
- " if fname is not None:\n",
- " f = functools.partial(f, fname=fname)\n",
- "\n",
- " if kind == \"homogeneous\":\n",
- " xs, ys = [np.linspace(*bound, npoints) for bound in bounds]\n",
- " data = {xy: f(xy) for xy in itertools.product(xs, ys)}\n",
- " learner = adaptive.Learner2D(f, bounds=bounds)\n",
- " learner.data = data\n",
- " d[\"learner_hom\"] = learner\n",
- " elif kind == \"adaptive\":\n",
- " learner = adaptive.Learner2D(f, bounds=bounds)\n",
- " if fname is not None:\n",
- " learner.load(fname)\n",
- " learner.data = {\n",
- " k: v for i, (k, v) in enumerate(learner.data.items()) if i <= npoints ** 2\n",
- " }\n",
- " adaptive.runner.simple(learner, goal=lambda l: l.npoints >= npoints ** 2)\n",
- " d[\"learner\"] = learner\n",
- "\n",
- " if with_tri:\n",
- " tri = learner.ip().tri\n",
- " triang = mtri.Triangulation(*tri.points.T, triangles=tri.vertices)\n",
- " ax.triplot(triang, c=\"w\", lw=0.2, alpha=0.8)\n",
- "\n",
- " values = np.array(list(learner.data.values()))\n",
- " ax.imshow(\n",
- " learner.plot(npoints if kind == \"homogeneous\" else None).Image.I.data,\n",
- " extent=(-0.5, 0.5, -0.5, 0.5),\n",
- " interpolation=\"none\",\n",
- " )\n",
- " ax.set_xticks([])\n",
- " ax.set_yticks([])\n",
- "\n",
- "axs[0][0].set_title(r\"$\\textrm{homogeneous}$\")\n",
- "axs[0][1].set_title(r\"$\\textrm{adaptive}$\")\n",
- "\n",
- "plt.savefig(\"figures/adaptive_2D.pdf\", bbox_inches=\"tight\", transparent=True)"
+ "learner = adaptive.Learner1D(funcs_1D[0]['function'], bounds=funcs_1D[0]['bounds'])\n",
+ "times = []\n",
+ "for i in range(10000):\n",
+ " t_start = time.time()\n",
+ " points, _ = learner.ask(1)\n",
+ " t_end = time.time()\n",
+ " times.append(t_end - t_start)\n",
+ " learner.tell(points[0], learner.function(points[0]))\n",
+ "plt.plot(np.cumsum(times))"
]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
}
],
"metadata": {
diff --git a/paper.md b/paper.md
index adc66dc304902a0744fe844724b37a02031c7425..ec4765729df596178771c0d61f135306431c0ca0 100755
--- a/paper.md
+++ b/paper.md
@@ -128,8 +128,16 @@ A more complex loss function that also takes the first neighboring intervals int
Figure @fig:adaptive_vs_grid shows a comparison between a result using this loss and a function that is sampled on a grid.
#### In general local loss functions only have a logarithmic overhead.
+<!-- Bas: not sure what to write here -->
#### With many points, due to the loss being local, parallel sampling incurs no additional cost.
+<!-- Bas: the text below doesn't really describe what's written above, but is an essential part nonetheless -->
+So far, the description of the general algorithm did not include parallelism.
+It needs to be able to suggest multiple points at the same time and remember which points it suggests.
+When a new point $\bm{x}_\textrm{new}$ with the largest loss $L_\textrm{max}$ is suggested, the interval it belongs to splits up into $N$ new intervals (here $N$ depends on the dimensionality of the function $f$.)
+A temporary loss $L_\textrm{temp} = L_\textrm{max}/N$ is assigned to these newly created intervals until $f(\bm{x})$ is calculated and the temporary loss can be replaced by the actual loss $L \equiv L((\bm{x},\bm{y})_\textrm{new}, (\bm{x},\bm{y})_\textrm{neigbors})$ of these new intervals, where $L \ge L_\textrm{temp}$.
+For a one-dimensional scalar function, this procedure is equivalent to temporarily using the function values of the neighbors of $x_\textrm{new}$ and assign the interpolated value to $y_\textrm{new}$ until it is known.
+When querying $n>1$ points, the former procedure simply repeats $n$ times.
# Loss function design