diff --git a/figures.ipynb b/figures.ipynb
index 81f605af2470d3146df7d13d0c55d3414c146967..a296429f830293cbfcd3653e34f92ebe7d5ce0c3 100644
--- a/figures.ipynb
+++ b/figures.ipynb
@@ -652,6 +652,216 @@
"for d in funcs_1D:\n",
" print(d[\"err_hom\"][-1] / d[\"err_adaptive\"][-1])"
]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Iso surface"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from functools import lru_cache\n",
+ "import numpy as np\n",
+ "import scipy.linalg\n",
+ "import scipy.spatial\n",
+ "import kwant\n",
+ "\n",
+ "\n",
+ "@lru_cache()\n",
+ "def create_syst(unit_cell):\n",
+ " lat = kwant.lattice.Polyatomic(unit_cell, [(0, 0, 0)])\n",
+ " syst = kwant.Builder(kwant.TranslationalSymmetry(*lat.prim_vecs))\n",
+ " syst[lat.shape(lambda _: True, (0, 0, 0))] = 6\n",
+ " syst[lat.neighbors()] = -1\n",
+ " return kwant.wraparound.wraparound(syst).finalized()\n",
+ "\n",
+ "\n",
+ "def get_brillouin_zone(unit_cell):\n",
+ " syst = create_syst(unit_cell)\n",
+ " A = get_A(syst)\n",
+ " neighbours = kwant.linalg.lll.voronoi(A)\n",
+ " lattice_points = np.concatenate(([[0, 0, 0]], neighbours))\n",
+ " lattice_points = 2 * np.pi * (lattice_points @ A.T)\n",
+ " vor = scipy.spatial.Voronoi(lattice_points)\n",
+ " brillouin_zone = vor.vertices[vor.regions[vor.point_region[0]]]\n",
+ " return scipy.spatial.ConvexHull(brillouin_zone)\n",
+ "\n",
+ "\n",
+ "def momentum_to_lattice(k, syst):\n",
+ " A = get_A(syst)\n",
+ " k, residuals = scipy.linalg.lstsq(A, k)[:2]\n",
+ " if np.any(abs(residuals) > 1e-7):\n",
+ " raise RuntimeError(\n",
+ " \"Requested momentum doesn't correspond to any lattice momentum.\"\n",
+ " )\n",
+ " return k\n",
+ "\n",
+ "\n",
+ "def get_A(syst):\n",
+ " B = np.asarray(syst._wrapped_symmetry.periods).T\n",
+ " return np.linalg.pinv(B).T\n",
+ "\n",
+ "\n",
+ "def energies(k, unit_cell):\n",
+ " syst = create_syst(unit_cell)\n",
+ " k_x, k_y, k_z = momentum_to_lattice(k, syst)\n",
+ " params = {\"k_x\": k_x, \"k_y\": k_y, \"k_z\": k_z}\n",
+ " H = syst.hamiltonian_submatrix(params=params)\n",
+ " energies = np.linalg.eigvalsh(H)\n",
+ " return min(energies)\n",
+ "\n",
+ "\n",
+ "from functools import partial\n",
+ "\n",
+ "from ipywidgets import interact_manual\n",
+ "\n",
+ "\n",
+ "# Define the lattice vectors of some common unit cells\n",
+ "lattices = dict(\n",
+ " hexegonal=((0, 1, 0), (np.cos(-np.pi / 6), np.sin(-np.pi / 6), 0), (0, 0, 1)),\n",
+ " simple_cubic=((1, 0, 0), (0, 1, 0), (0, 0, 1)),\n",
+ " fcc=((0, 0.5, 0.5), (0.5, 0.5, 0), (0.5, 0, 0.5)),\n",
+ " bcc=((-0.5, 0.5, 0.5), (0.5, -0.5, 0.5), (0.5, 0.5, -0.5)),\n",
+ ")\n",
+ "\n",
+ "\n",
+ "def isoline_loss_function(y_iso, sigma, priority=1):\n",
+ " from adaptive.learner.learnerND import default_loss\n",
+ "\n",
+ " def gaussian(x, mu, sigma):\n",
+ " return np.exp(-(x - mu) ** 2 / sigma ** 2)\n",
+ "\n",
+ " def loss(simplex, ys):\n",
+ " distance = np.mean([abs(y_iso - y) for y in ys])\n",
+ " return priority * gaussian(distance, 0, sigma) + default_loss(simplex, ys)\n",
+ "\n",
+ " return loss\n",
+ "\n",
+ "\n",
+ "loss = isoline_loss_function(y_iso=0.1, sigma=1, priority=0.0)\n",
+ "\n",
+ "learners = []\n",
+ "fnames = []\n",
+ "for name, unit_cell in lattices.items():\n",
+ " hull = get_brillouin_zone(unit_cell)\n",
+ " learner = adaptive.LearnerND(\n",
+ " partial(energies, unit_cell=unit_cell), hull, loss_per_simplex=loss\n",
+ " )\n",
+ " fnames.append(name)\n",
+ " learners.append(learner)\n",
+ "\n",
+ "mapping = dict(zip(fnames, learners))\n",
+ "\n",
+ "learner = adaptive.BalancingLearner(learners, strategy=\"npoints\")\n",
+ "\n",
+ "\n",
+ "def select(name, learners, fnames):\n",
+ " return learners[fnames.index(name)]\n",
+ "\n",
+ "\n",
+ "def iso(unit_cell, level=8.5):\n",
+ " l = select(unit_cell, learners, fnames)\n",
+ " adaptive.runner.simple(l, goal=lambda l: l.npoints > 1000)\n",
+ " return l.plot_isosurface(level=level)\n",
+ "\n",
+ "interact_manual(iso, level=(-6, 9, 0.1), unit_cell=lattices.keys())"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def f(xy):\n",
+ " x, y = xy\n",
+ " return x ** 2 + y ** 3\n",
+ "\n",
+ "bounds = [(-1, 1), (-1, 1)]\n",
+ "npoints = 300\n",
+ "\n",
+ "def isoline_loss_function(y_iso, sigma, priority=1):\n",
+ " from adaptive.learner.learnerND import default_loss\n",
+ "\n",
+ " def gaussian(x, mu, sigma):\n",
+ " return np.exp(-(x - mu) ** 2 / sigma ** 2 / 2)\n",
+ "\n",
+ " def loss(simplex, values, value_scale):\n",
+ " distance = np.mean([abs(y_iso * value_scale - y) for y in values])\n",
+ " return priority * gaussian(distance, 0, sigma) + default_loss(simplex, values, value_scale)\n",
+ "\n",
+ " return loss\n",
+ "\n",
+ "\n",
+ "level = 0.1\n",
+ "loss = isoline_loss_function(level, sigma=0.4, priority=0.5)\n",
+ "\n",
+ "with_tri = True\n",
+ "fig, axs = plt.subplots(1, 2, figsize=(fig_width, fig_height))\n",
+ "plt.subplots_adjust(wspace=0.3)\n",
+ "\n",
+ "for i, ax in enumerate(axs.flatten()):\n",
+ " learner = adaptive.LearnerND(f, bounds, loss_per_simplex=loss)\n",
+ " label = \"ab\"[i]\n",
+ " ax.text(\n",
+ " 0.5,\n",
+ " 1.05,\n",
+ " f\"$\\mathrm{{({label})}}$\",\n",
+ " transform=ax.transAxes,\n",
+ " horizontalalignment=\"center\",\n",
+ " verticalalignment=\"bottom\",\n",
+ " )\n",
+ " ax.xaxis.set_ticks([])\n",
+ " ax.yaxis.set_ticks([])\n",
+ " kind = \"homogeneous\" if i == 0 else \"adaptive\"\n",
+ "\n",
+ " if kind == \"homogeneous\":\n",
+ " xs, ys = [np.linspace(*bound, int(npoints**0.5)) for bound in bounds]\n",
+ " data = {xy: f(xy) for xy in itertools.product(xs, ys)}\n",
+ " learner.data = data\n",
+ " elif kind == \"adaptive\":\n",
+ " learner.data = {\n",
+ " k: v for i, (k, v) in enumerate(learner.data.items()) if i <= npoints\n",
+ " }\n",
+ " adaptive.runner.simple(learner, goal=lambda l: l.npoints >= npoints)\n",
+ "\n",
+ " if with_tri:\n",
+ " xy = np.array([learner.tri.get_vertices(s)\n",
+ " for s in learner.tri.simplices]).reshape(-1, 2)\n",
+ " triang = mtri.Triangulation(*xy.T)\n",
+ " ax.triplot(triang, c=\"w\", lw=0.2, alpha=0.8)\n",
+ "\n",
+ " # Isolines\n",
+ " vertices, lines = learner._get_iso(level, which='line')\n",
+ " paths = np.array([[vertices[i], vertices[j]] for i, j in lines]).T\n",
+ " ax.plot(*paths, c='k')\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=(-1, 1, -1, 1),\n",
+ " interpolation=\"none\",\n",
+ " )\n",
+ " ax.set_xticks([])\n",
+ " ax.set_yticks([])\n",
+ "\n",
+ "axs[0].set_ylabel(r\"$\\textrm{homogeneous}$\")\n",
+ "axs[1].set_ylabel(r\"$\\textrm{adaptive}$\")\n",
+ "plt.savefig(\"figures/isoline.pdf\", bbox_inches=\"tight\", transparent=True)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
}
],
"metadata": {
diff --git a/figures/isoline.pdf b/figures/isoline.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..8eb65518621163015de1a5a30931091ef794cf98
Binary files /dev/null and b/figures/isoline.pdf differ
diff --git a/paper.md b/paper.md
index 093ae89a07c9a81ff9b4677c6141e88d51808ce2..e44df2a5dae5126f1beec060fe1413fb75eed837 100755
--- a/paper.md
+++ b/paper.md
@@ -230,8 +230,14 @@ Here, we see that for homogeneous sampling to get the same error as sampling wit
#### The `cquad` algorithm belongs to a class that is parallelizable.
-## isosurface sampling
-<!-- figure here -->
+## isoline and isosurface sampling
+We can find isolines or isosurfaces using a loss function that prioritizes intervals that are closer to the function values that we are interested in.
+See @fig:isoline.
+
+{#fig:isoline}
# Implementation and benchmarks