Skip to content
Snippets Groups Projects
Commit 7c34ff07 authored by Bas Nijholt's avatar Bas Nijholt
Browse files

consistent colors in fig

parent 5c4f2c88
No related branches found
No related tags found
No related merge requests found
Pipeline #21708 passed
Stage: test
Hide whitespace changes
Inline Side-by-side
Showing
with 56 additions and 30 deletions
figures.ipynb
+ 56
− 30
View file @ 7c34ff07
%% Cell type:markdown id: tags:
# Adaptive paper figures
%% Cell type:markdown id: tags:
## imports and function definitions
%% Cell type:code id: tags:
```
import functools
import itertools
import os
import pickle
import holoviews.plotting.mpl
import matplotlib; matplotlib.use("agg")
import matplotlib
matplotlib.use("agg")
import matplotlib.pyplot as plt
import matplotlib.tri as mtri
import numpy as np
from scipy import interpolate
import adaptive
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
golden_mean = (np.sqrt(5) - 1) / 2 # Aesthetic ratio
fig_width_pt = 246.0 # Columnwidth
inches_per_pt = 1 / 72.27 # Convert pt to inches
fig_width = fig_width_pt * inches_per_pt
fig_height = fig_width * golden_mean # height in inches
fig_size = [fig_width, fig_height]
params = {
"backend": "ps",
"axes.labelsize": 13,
"font.size": 13,
"legend.fontsize": 10,
"xtick.labelsize": 10,
"ytick.labelsize": 10,
"text.usetex": True,
"figure.figsize": fig_size,
"font.family": "serif",
"font.serif": "Computer Modern Roman",
"legend.frameon": True,
"savefig.dpi": 300,
}
plt.rcParams.update(params)
plt.rc("text.latex", preamble=[r"\usepackage{xfrac}", r"\usepackage{siunitx}"])
class HistogramNormalize(matplotlib.colors.Normalize):
def __init__(self, data, vmin=None, vmax=None, mixing_degree=1):
self.mixing_degree = mixing_degree
if vmin is not None:
data = data[data >= vmin]
if vmax is not None:
data = data[data <= vmax]
self.sorted_data = np.sort(data.flatten())
matplotlib.colors.Normalize.__init__(self, vmin, vmax)
def __call__(self, value, clip=None):
hist_norm = np.ma.masked_array(
np.searchsorted(self.sorted_data, value) / len(self.sorted_data)
)
linear_norm = super().__call__(value, clip)
return self.mixing_degree * hist_norm + (1 - self.mixing_degree) * linear_norm
```
%% Cell type:code id: tags:
```
# 1D funcs
def f(x, offset=0.123):
a = 0.02
return x + a ** 2 / (a ** 2 + (x - offset) ** 2)
def g(x):
return np.tanh(x * 40)
def h(x):
return np.sin(100 * x) * np.exp(-x ** 2 / 0.1 ** 2)
funcs_1D = [
dict(function=f, bounds=(-1, 1), title="peak"),
dict(function=g, bounds=(-1, 1), title="tanh"),
dict(function=h, bounds=(-0.3, 0.3), title="wave packet"),
]
# 2D funcs
def ring(xy, offset=0.123):
a = 0.2
x, y = xy
return x + np.exp(-(x ** 2 + y ** 2 - 0.75 ** 2) ** 2 / a ** 4)
@functools.lru_cache()
def phase_diagram_setup(fname):
data = adaptive.utils.load(fname)
points = np.array(list(data.keys()))
values = np.array(list(data.values()), dtype=float)
bounds = [
(points[:, 0].min(), points[:, 0].max()),
(points[:, 1].min(), points[:, 1].max()),
]
ll, ur = np.reshape(bounds, (2, 2)).T
inds = np.all(np.logical_and(ll <= points, points <= ur), axis=1)
points, values = points[inds], values[inds].reshape(-1, 1)
return interpolate.LinearNDInterpolator(points, values), bounds
def phase_diagram(xy, fname):
ip, _ = phase_diagram_setup(fname)
zs = ip(xy)
return np.round(zs)
def density(x, eps=0):
e = [0.8, 0.2]
delta = [0.5, 0.5, 0.5]
c = 3
omega = [0.02, 0.05]
H = np.array(
[
[e[0] + 1j * omega[0], delta[0], delta[1]],
[delta[0], e[1] + c * x + 1j * omega[1], delta[1]],
[delta[1], delta[2], e[1] - c * x + 1j * omega[1]],
]
)
H += np.eye(3) * eps
return np.trace(np.linalg.inv(H)).imag
def level_crossing(xy):
x, y = xy
return density(x, y) + y
funcs_2D = [
dict(function=ring, bounds=[(-1, 1), (-1, 1)], npoints=33, title="ring"),
dict(
function=functools.partial(phase_diagram, fname="phase_diagram.pickle"),
bounds=phase_diagram_setup("phase_diagram.pickle")[1],
npoints=100,
fname="phase_diagram.pickle",
title="phase diagram",
),
dict(function=level_crossing, bounds=[(-1, 1), (-3, 3)], npoints=50, title="level crossing"),
]
```
%% Cell type:markdown id: tags:
## Interval and loss
%% Cell type:code id: tags:
```
np.random.seed(1)
xs = np.array([0.1, 0.3, 0.35, 0.45])
f = lambda x: x ** 3
ys = f(xs)
means = lambda x: np.convolve(x, np.ones(2) / 2, mode="valid")
xs_means = means(xs)
ys_means = means(ys)
fig, ax = plt.subplots(figsize=fig_size)
ax.scatter(xs, ys, c="k")
ax.plot(xs, ys, c="k")
ax.annotate(
s=r"$L_{1,2} = \sqrt{\Delta x^2 + \Delta y^2}$",
xy=(np.mean([xs[0], xs[1]]), np.mean([ys[0], ys[1]])),
xytext=(xs[0] + 0.05, ys[0] - 0.05),
arrowprops=dict(arrowstyle="->"),
ha="center",
zorder=10,
)
for i, (x, y) in enumerate(zip(xs, ys)):
sign = [1, -1][i % 2]
ax.annotate(
s=fr"$x_{i+1}, y_{i+1}$",
xy=(x, y),
xytext=(x + 0.01, y + sign * 0.04),
arrowprops=dict(arrowstyle="->"),
ha="center",
)
ax.scatter(xs, ys, c="green", s=5, zorder=5, label="existing data")
losses = np.hypot(xs[1:] - xs[:-1], ys[1:] - ys[:-1])
ax.scatter(
xs_means, ys_means, c="red", s=300 * losses, zorder=8, label="candidate points"
)
xs_dense = np.linspace(xs[0], xs[-1], 400)
ax.plot(xs_dense, f(xs_dense), alpha=0.3, zorder=7, label="function")
ax.legend()
ax.axis("off")
plt.savefig("figures/loss_1D.pdf", bbox_inches="tight", transparent=True)
plt.show()
```
%% Cell type:markdown id: tags:
## Learner1D vs grid
%% Cell type:code id: tags:
```
fig, axs = plt.subplots(2, len(funcs_1D), figsize=(fig_width, 1.5 * fig_height))
n_points = 50
for i, ax in enumerate(axs.T.flatten()):
ax.xaxis.set_ticks([])
ax.yaxis.set_ticks([])
kind = "homogeneous" if i % 2 == 0 else "adaptive"
index = i // 2 if kind == "homogeneous" else (i - 1) // 2
d = funcs_1D[index]
bounds = d["bounds"]
f = d["function"]
if kind == "homogeneous":
xs = np.linspace(*bounds, n_points)
ys = f(xs)
ax.set_title(rf"\textrm{{{d['title']}}}")
elif kind == "adaptive":
loss = adaptive.learner.learner1D.curvature_loss_function()
learner = adaptive.Learner1D(f, bounds=bounds, loss_per_interval=loss)
adaptive.runner.simple(learner, goal=lambda l: l.npoints >= n_points)
xs, ys = zip(*sorted(learner.data.items()))
xs_dense = np.linspace(*bounds, 1000)
ax.plot(xs_dense, f(xs_dense), c="red", alpha=0.3, lw=0.5)
ax.scatter(xs, ys, s=0.5, c="k")
axs[0][0].set_ylabel(r"$\textrm{homogeneous}$")
axs[1][0].set_ylabel(r"$\textrm{adaptive}$")
plt.savefig("figures/Learner1D.pdf", bbox_inches="tight", transparent=True)
```
%% Cell type:markdown id: tags:
## Learner2D vs grid
%% Cell type:code id: tags:
```
fig, axs = plt.subplots(2, len(funcs_2D), figsize=(fig_width, 1.5 * fig_height))
plt.subplots_adjust(hspace=-0.1, wspace=0.1)
with_tri = False
for i, ax in enumerate(axs.T.flatten()):
label = "abcdef"[i]
ax.text(
0.5,
1.05,
f"$\mathrm{{({label})}}$",
transform=ax.transAxes,
horizontalalignment="center",
verticalalignment="bottom",
)
ax.xaxis.set_ticks([])
ax.yaxis.set_ticks([])
kind = "homogeneous" if i % 2 == 0 else "adaptive"
d = funcs_2D[i // 2] if kind == "homogeneous" else funcs_2D[(i - 1) // 2]
bounds = d["bounds"]
npoints = d["npoints"]
f = d["function"]
fname = d.get("fname")
if fname is not None:
f = functools.partial(f, fname=fname)
if kind == "homogeneous":
xs, ys = [np.linspace(*bound, npoints) for bound in bounds]
data = {xy: f(xy) for xy in itertools.product(xs, ys)}
learner = adaptive.Learner2D(f, bounds=bounds)
learner.data = data
d["learner_hom"] = learner
elif kind == "adaptive":
learner = adaptive.Learner2D(f, bounds=bounds)
if fname is not None:
learner.load(fname)
learner.data = {
k: v for i, (k, v) in enumerate(learner.data.items()) if i <= npoints ** 2
}
adaptive.runner.simple(learner, goal=lambda l: l.npoints >= npoints ** 2)
d["learner"] = learner
if with_tri:
tri = learner.ip().tri
triang = mtri.Triangulation(*tri.points.T, triangles=tri.vertices)
ax.triplot(triang, c="w", lw=0.2, alpha=0.8)
values = np.array(list(learner.data.values()))
ax.imshow(
learner.plot(npoints if kind == "homogeneous" else None).Image.I.data,
plot_data = learner.plot(npoints if kind == "homogeneous" else None).Image.I.data
im = ax.imshow(
plot_data,
extent=(-0.5, 0.5, -0.5, 0.5),
interpolation="none",
)
ax.set_xticks([])
ax.set_yticks([])
if i in [2, 3]:
norm = HistogramNormalize(plot_data, mixing_degree=0.6)
im.set_norm(norm)
axs[0][0].set_ylabel(r"$\textrm{homogeneous}$")
axs[1][0].set_ylabel(r"$\textrm{adaptive}$")
plt.savefig("figures/Learner2D.pdf", bbox_inches="tight", transparent=True)
```
%% Cell type:code id: tags:
```
import time
learner = adaptive.Learner1D(funcs_1D[0]["function"], bounds=funcs_1D[0]["bounds"])
times = []
for i in range(10000):
t_start = time.time()
points, _ = learner.ask(1)
t_end = time.time()
times.append(t_end - t_start)
learner.tell(points[0], learner.function(points[0]))
plt.plot(np.cumsum(times))
```
%% Cell type:markdown id: tags:
## Algorithm explaination: Learner1D step-by-step
%% Cell type:code id: tags:
```
fig, axs = plt.subplots(1, 1, figsize=(fig_width, 3 * fig_height))
def f(x, offset=0.123):
a = 0.2
return a ** 2 / (a ** 2 + (x - offset) ** 2)
learner = adaptive.Learner1D(
f, (0, 2), loss_per_interval=adaptive.learner.learner1D.curvature_loss_function()
)
learner._recompute_losses_factor = 0.1
xs_dense = np.linspace(*learner.bounds, 400)
ys_dense = f(xs_dense)
step = 0.4
for i in range(11):
offset = -i * step
x = learner.ask(1)[0][0]
y = f(x)
learner.tell(x, y)
xs, ys = map(np.array, zip(*sorted(learner.data.items())))
ys = ys + offset
if i >= 1:
axs.plot(xs_dense, ys_dense + offset, c="k", alpha=0.3, zorder=0)
axs.plot(xs, ys, zorder=1, c="k")
axs.scatter(xs, ys, alpha=1, zorder=2, c="k")
(x_left, x_right), loss = list(learner.losses.items())[
0
] # it's a ItemSortedDict
] # it's an ItemSortedDict
(y_left, y_right) = [
learner.data[x_left] + offset,
learner.data[x_right] + offset,
]
axs.scatter([x_left, x_right], [y_left, y_right], c="r", s=10, zorder=3)
x_mid = np.mean((x_left, x_right))
y_mid = np.interp(x_mid, (x_left, x_right), (y_left, y_right))
axs.scatter(x_mid, y_mid, zorder=4, marker="x", c="green")
axs.text(
-0.1,
0.5,
(r"$\mathrm{iterations}$" + "\n" + "$\longleftarrow$"),
transform=axs.transAxes,
horizontalalignment="center",
verticalalignment="center",
rotation=90,
fontsize=18,
)
# legend
import matplotlib.patches as mpatches
import matplotlib.lines as mlines
class LargestInterval:
pass
class Interval:
pass
class Function:
pass
class IntervalHandler:
def __init__(self, with_inner=True, length=20, *args, **kwargs):
super().__init__(*args, **kwargs)
self.with_inner = with_inner
self.length = length
def legend_artist(self, legend, orig_handle, fontsize, handlebox):
x0, y0 = handlebox.xdescent, handlebox.ydescent
offsets = [0, self.length]
line = mlines.Line2D((0, offsets[-1]), (0, 0), zorder=0, c="k")
handlebox.add_artist(line)
for offset in offsets:
circle1 = mpatches.Circle(
[x0 + offset, y0], 4, facecolor="k", lw=3, zorder=1
)
handlebox.add_artist(circle1)
if self.with_inner:
circle2 = mpatches.Circle(
[x0 + offset, y0], 3, facecolor="red", lw=3, zorder=1
)
handlebox.add_artist(circle2)
class FunctionHandler:
def __init__(self, xs, ys, *args, **kwargs):
super().__init__(*args, **kwargs)
self.xs = xs / xs.ptp() * 20
self.ys = ys - ys.mean()
def legend_artist(self, legend, orig_handle, fontsize, handlebox):
x0, y0 = handlebox.xdescent, handlebox.ydescent
line = mlines.Line2D(self.xs, self.ys * 10, zorder=0, c="k", alpha=0.3)
handlebox.add_artist(line)
plt.legend(
[
Function(),
mlines.Line2D([], [], marker="o", lw=0, c="k"),
LargestInterval(),
Interval(),
mlines.Line2D([], [], marker="x", lw=0, c="green"),
],
[
"original function",
"known point",
"interval",
"largest loss interval",
"next candidate point",
],
handler_map={
LargestInterval: IntervalHandler(False),
Interval: IntervalHandler(True),
Function: FunctionHandler(xs, ys),
},
bbox_to_anchor=(0.25, 0.9, 1.0, 0.0),
ncol=1,
)
# On grid, uncomment for homogeneous plot
# axs.plot(learner.bounds, [-(i + 0.5) * step, -(i + 0.5) * step], c='k', ls='--')
# xs_hom = np.linspace(*learner.bounds, i)
# ys_hom = f(xs_hom) - (i + 3) * step
# axs.plot(xs_hom, ys_hom, zorder=1, c="k")
# axs.scatter(xs_hom, ys_hom, alpha=1, zorder=2, c="k")
axs.axis("off")
plt.savefig("figures/algo.pdf", bbox_inches="tight", transparent=True)
plt.show()
```
%% Cell type:markdown id: tags:
## Line loss visualization
%% Cell type:code id: tags:
```
from matplotlib.patches import Polygon
fig, axs = plt.subplots(5, 1, figsize=(fig_width, 1.5 * fig_height))
f = lambda x: np.sin(x) ** 2
xs = np.array([0, 1.3, 3, 5, 7, 8])
ys = f(xs)
def plot(xs, ax):
ys = f(xs)
xs_dense = np.linspace(xs[0], xs[-1], 300)
ys_dense = f(xs_dense)
ax.plot(xs_dense, ys_dense, alpha=0.3, c="k")
ax.plot(xs, ys, c="k")
ax.scatter(xs, ys, zorder=10, s=14, c="k")
plot(xs, axs[0])
plot(xs, axs[1])
for i, ax in enumerate(axs):
ax.axis("off")
ax.set_ylim(-1.5, 1.5)
label = "abcde"[i]
label = "abcdefgh"[i]
ax.text(
0.5,
0.9,
f"$\mathrm{{({label})}}$",
transform=ax.transAxes,
horizontalalignment="center",
verticalalignment="bottom",
)
def plot_tri(xs, ax):
def plot_tri(xs, ax, colors):
ys = f(xs)
for i in range(len(xs)):
if i == 0 or i == len(xs) - 1:
continue
color = f"C{i}"
verts = [(xs[i - 1], ys[i - 1]), (xs[i], ys[i]), (xs[i + 1], ys[i + 1])]
poly = Polygon(verts, facecolor=color, alpha=0.4)
poly = Polygon(verts, facecolor=colors[xs[i]], alpha=0.4)
ax.add_patch(poly)
ax.scatter([xs[i]], [ys[i]], c=color, s=6, zorder=11)
ax.scatter([xs[i]], [ys[i]], c=colors[xs[i]], s=6, zorder=11)
ax.plot(
[xs[i - 1], xs[i + 1]], [ys[i - 1], ys[i + 1]], c="k", ls="--", alpha=0.3
)
plot_tri(xs, axs[1])
learner = adaptive.Learner1D(
f,
(xs[0], xs[-1]),
loss_per_interval=adaptive.learner.learner1D.curvature_loss_function(),
)
learner.tell_many(xs, ys)
for i in [2, 3, 4]:
ax = axs[i]
learner = adaptive.Learner1D(
f,
(xs[0], xs[-1]),
loss_per_interval=adaptive.learner.learner1D.curvature_loss_function(),
)
learner.tell_many(xs, ys)
x_new = learner.ask(1)[0][0]
learner.tell(x_new, f(x_new))
xs, ys = zip(*sorted(learner.data.items()))
for i, ax in enumerate(axs[1:]):
if i != 0:
x_new = learner.ask(1)[0][0]
learner.tell(x_new, f(x_new))
xs, ys = zip(*sorted(learner.data.items()))
plot(xs, ax)
plot_tri(xs, ax)
colors = {x: f"C{i}" for i, x in enumerate(learner.data.keys())}
plot_tri(xs, ax, colors=colors)
plt.savefig("figures/line_loss.pdf", bbox_inches="tight", transparent=True)
plt.show()
```
%% Cell type:markdown id: tags:
## Line loss L1-error(N)
%% Cell type:code id: tags:
```
import collections
from adaptive import Learner1D, LearnerND
from scipy import interpolate
def err_1D(xs, ys, xs_rand, f):
ip = interpolate.interp1d(xs, ys)
abserr = np.abs(ip(xs_rand) - f(xs_rand))
return np.average(abserr ** 2) ** 0.5
def get_err_1D(learner, N):
xs_rand = np.random.uniform(*learner.bounds, size=100_000)
xs = np.linspace(*learner.bounds, N)
ys = learner.function(xs)
err_lin = err_1D(xs, ys, xs_rand, learner.function)
xs, ys = zip(*learner.data.items())
xs, ys = xs[:N], ys[:N]
ip = interpolate.interp1d(xs, ys)
err_learner = err_1D(xs, ys, xs_rand, learner.function)
return err_lin, err_learner
def err_2D(zs, zs_rand):
abserr = np.abs(zs - zs_rand)
return np.average(abserr ** 2) ** 0.5
def get_err_2D(learner, N):
xys_rand = np.vstack(
[
np.random.uniform(*learner.bounds[0], size=int(100_000)),
np.random.uniform(*learner.bounds[1], size=int(100_000)),
]
)
xys_hom = np.array(
[
(x, y)
for x in np.linspace(*learner.bounds[0], int(N ** 0.5))
for y in np.linspace(*learner.bounds[1], int(N ** 0.5))
]
)
N = len(xys_hom)
try:
# Vectorized
zs_hom = learner.function(xys_hom.T)
zs_rand = learner.function(xys_rand)
except:
# Non-vectorized
zs_hom = np.array([learner.function(xy) for xy in xys_hom])
zs_rand = np.array([learner.function(xy) for xy in xys_rand.T])
ip = interpolate.LinearNDInterpolator(xys_hom, zs_hom)
zs = ip(xys_rand.T)
err_lin = err_2D(zs, zs_rand)
ip = interpolate.LinearNDInterpolator(learner.points[:N], learner.values[:N])
zs = ip(xys_rand.T)
err_learner = err_2D(zs, zs_rand)
return err_lin, err_learner
recalculate = False
N_max = 10_000
Ns = np.geomspace(4, N_max, 50).astype(int)
fname = "error_line_loss.pickle"
if not os.path.exists(fname) and not recalculate:
loss_1D = adaptive.learner.learner1D.curvature_loss_function()
loss_2D = adaptive.learner.learnerND.curvature_loss_function()
err = collections.defaultdict(dict)
for i, (funcs, loss, Learner, get_err) in enumerate(
zip(
[funcs_1D, funcs_2D],
[loss_1D, loss_2D],
[Learner1D, LearnerND],
[get_err_1D, get_err_2D],
)
):
for d in funcs:
learner = Learner(d["function"], d["bounds"], loss)
adaptive.runner.simple(learner, goal=lambda l: l.npoints >= N_max)
errs = [get_err(learner, N) for N in Ns]
err_hom, err_adaptive = zip(*errs)
err[i][d["title"]] = (err_hom, err_adaptive)
with open(fname, "wb") as f:
pickle.dump(err, f)
else:
with open(fname, "rb") as f:
err = pickle.load(f)
```
%% Cell type:code id: tags:
```
fig, axs = plt.subplots(2, 1, figsize=(fig_width, 1.6 * fig_height))
plt.subplots_adjust(hspace=0.3)
axs[1].set_xlabel("$N$")
for i, ax in enumerate(axs):
ax.set_ylabel(r"$\text{Err}_{1}(\tilde{f})$")
for j, (title, (err_hom, err_adaptive)) in enumerate(err[i].items()):
color = f"C{j}"
label = "abc"[j]
label = f"$\mathrm{{({label})}}$ {title}"
# ax.loglog(Ns, err_hom, ls="--", c=color)
# ax.loglog(Ns, err_adaptive, label=label, c=color)
error = np.array(err_hom) / np.array(err_adaptive)
if i == 0:
ax.loglog(Ns[:36], error[:36], c=color, label=label)
else:
ax.loglog(Ns, error, c=color, label=label)
ax.loglog(Ns, err_hom, ls="--", c=color)
ax.loglog(Ns, err_adaptive, label=label, c=color)
# error = np.array(err_hom) / np.array(err_adaptive)
# if i == 0:
# ax.loglog(Ns[:36], error[:36], c=color, label=label)
# else:
# ax.loglog(Ns, error, c=color, label=label)
ax.legend()
plt.savefig("figures/line_loss_error.pdf", bbox_inches="tight", transparent=True)
plt.show()
```
%% Cell type:code id: tags:
```
# Error reduction
print("1D")
for title, (err_hom, err_adaptive) in err[0].items():
print(title, err_hom[-1] / err_adaptive[-1])
print("\n2D")
for title, (err_hom, err_adaptive) in err[1].items():
print(title, err_hom[-1] / err_adaptive[-1])
```
%% Cell type:markdown id: tags:
## iso-line plots
%% Cell type:code id: tags:
```
from functools import lru_cache
import numpy as np
import scipy.linalg
import scipy.spatial
import kwant
@lru_cache()
def create_syst(unit_cell):
lat = kwant.lattice.Polyatomic(unit_cell, [(0, 0, 0)])
syst = kwant.Builder(kwant.TranslationalSymmetry(*lat.prim_vecs))
syst[lat.shape(lambda _: True, (0, 0, 0))] = 6
syst[lat.neighbors()] = -1
return kwant.wraparound.wraparound(syst).finalized()
def get_brillouin_zone(unit_cell):
syst = create_syst(unit_cell)
A = get_A(syst)
neighbours = kwant.linalg.lll.voronoi(A)
lattice_points = np.concatenate(([[0, 0, 0]], neighbours))
lattice_points = 2 * np.pi * (lattice_points @ A.T)
vor = scipy.spatial.Voronoi(lattice_points)
brillouin_zone = vor.vertices[vor.regions[vor.point_region[0]]]
return scipy.spatial.ConvexHull(brillouin_zone)
def momentum_to_lattice(k, syst):
A = get_A(syst)
k, residuals = scipy.linalg.lstsq(A, k)[:2]
if np.any(abs(residuals) > 1e-7):
raise RuntimeError(
"Requested momentum doesn't correspond to any lattice momentum."
)
return k
def get_A(syst):
B = np.asarray(syst._wrapped_symmetry.periods).T
return np.linalg.pinv(B).T
def energies(k, unit_cell):
syst = create_syst(unit_cell)
k_x, k_y, k_z = momentum_to_lattice(k, syst)
params = {"k_x": k_x, "k_y": k_y, "k_z": k_z}
H = syst.hamiltonian_submatrix(params=params)
energies = np.linalg.eigvalsh(H)
return min(energies)
from functools import partial
from ipywidgets import interact_manual
# Define the lattice vectors of some common unit cells
lattices = dict(
hexegonal=((0, 1, 0), (np.cos(-np.pi / 6), np.sin(-np.pi / 6), 0), (0, 0, 1)),
simple_cubic=((1, 0, 0), (0, 1, 0), (0, 0, 1)),
fcc=((0, 0.5, 0.5), (0.5, 0.5, 0), (0.5, 0, 0.5)),
bcc=((-0.5, 0.5, 0.5), (0.5, -0.5, 0.5), (0.5, 0.5, -0.5)),
)
def isoline_loss_function(y_iso, sigma, priority=1):
from adaptive.learner.learnerND import default_loss
def gaussian(x, mu, sigma):
return np.exp(-(x - mu) ** 2 / sigma ** 2)
def loss(simplex, ys):
distance = np.mean([abs(y_iso - y) for y in ys])
return priority * gaussian(distance, 0, sigma) + default_loss(simplex, ys)
return loss
loss = isoline_loss_function(y_iso=0.1, sigma=1, priority=0.0)
learners = []
fnames = []
for name, unit_cell in lattices.items():
hull = get_brillouin_zone(unit_cell)
learner = adaptive.LearnerND(
partial(energies, unit_cell=unit_cell), hull, loss_per_simplex=loss
)
fnames.append(name)
learners.append(learner)
mapping = dict(zip(fnames, learners))
learner = adaptive.BalancingLearner(learners, strategy="npoints")
def select(name, learners, fnames):
return learners[fnames.index(name)]
def iso(unit_cell, level=8.5):
l = select(unit_cell, learners, fnames)
adaptive.runner.simple(l, goal=lambda l: l.npoints > 1000)
return l.plot_isosurface(level=level)
interact_manual(iso, level=(-6, 9, 0.1), unit_cell=lattices.keys())
```
%% Cell type:code id: tags:
```
def f(xy):
x, y = xy
return x ** 2 + y ** 3
bounds = [(-1, 1), (-1, 1)]
npoints = 17 ** 2
def isoline_loss_function(y_iso, sigma, priority=1):
from adaptive.learner.learnerND import default_loss
def gaussian(x, mu, sigma):
return np.exp(-(x - mu) ** 2 / sigma ** 2 / 2)
def loss(simplex, values, value_scale):
distance = np.mean([abs(y_iso * value_scale - y) for y in values])
return priority * gaussian(distance, 0, sigma) + default_loss(
simplex, values, value_scale
)
return loss
level = 0.1
loss = isoline_loss_function(level, sigma=0.4, priority=0.5)
with_tri = True
fig, axs = plt.subplots(1, 2, figsize=(fig_width, fig_height))
plt.subplots_adjust(wspace=0.3)
for i, ax in enumerate(axs.flatten()):
learner = adaptive.LearnerND(f, bounds, loss_per_simplex=loss)
label = "ab"[i]
ax.text(
0.5,
1.05,
f"$\mathrm{{({label})}}$",
transform=ax.transAxes,
horizontalalignment="center",
verticalalignment="bottom",
)
ax.xaxis.set_ticks([])
ax.yaxis.set_ticks([])
kind = "homogeneous" if i == 0 else "adaptive"
if kind == "homogeneous":
xs, ys = [np.linspace(*bound, int(npoints ** 0.5)) for bound in bounds]
data = {xy: f(xy) for xy in itertools.product(xs, ys)}
learner.data = data
elif kind == "adaptive":
learner.data = {
k: v for i, (k, v) in enumerate(learner.data.items()) if i <= npoints
}
adaptive.runner.simple(learner, goal=lambda l: l.npoints >= npoints)
if with_tri:
xy = np.array([learner.tri.get_vertices(s) for s in learner.tri.simplices])
print(f"{len(xy)} triangles")
triang = mtri.Triangulation(*xy.reshape(-1, 2).T)
ax.triplot(triang, c="w", lw=0.2, alpha=0.8)
# Isolines
vertices, lines = learner._get_iso(level, which="line")
paths = np.array([[vertices[i], vertices[j]] for i, j in lines]).T
print("{} line segments".format(len(paths.T)))
ax.plot(*paths, c="k")
values = np.array(list(learner.data.values()))
ax.imshow(
learner.plot(npoints if kind == "homogeneous" else None).Image.I.data,
extent=(-1, 1, -1, 1),
interpolation="none",
)
ax.set_xticks([])
ax.set_yticks([])
axs[0].set_ylabel(r"$\textrm{homogeneous}$")
axs[1].set_ylabel(r"$\textrm{adaptive}$")
plt.savefig("figures/isoline.pdf", bbox_inches="tight", transparent=True)
```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment