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

add 2D figure

parent 293d05a9
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 178 additions and 17 deletions
figures.ipynb
+ 178
− 17
View file @ 3691b864
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` ```
import numpy as np import numpy as np
import matplotlib import matplotlib
matplotlib.use("agg") matplotlib.use("agg")
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
%matplotlib inline %matplotlib inline
%config InlineBackend.figure_format = 'svg' %config InlineBackend.figure_format = 'svg'
golden_mean = (np.sqrt(5) - 1) / 2 # Aesthetic ratio golden_mean = (np.sqrt(5) - 1) / 2 # Aesthetic ratio
fig_width_pt = 246.0 # Columnwidth fig_width_pt = 246.0 # Columnwidth
inches_per_pt = 1 / 72.27 # Convert pt to inches inches_per_pt = 1 / 72.27 # Convert pt to inches
fig_width = fig_width_pt * inches_per_pt fig_width = fig_width_pt * inches_per_pt
fig_height = fig_width * golden_mean # height in inches fig_height = fig_width * golden_mean # height in inches
fig_size = [fig_width, fig_height] fig_size = [fig_width, fig_height]
params = { params = {
"backend": "ps", "backend": "ps",
"axes.labelsize": 13, "axes.labelsize": 13,
"font.size": 13, "font.size": 13,
"legend.fontsize": 10, "legend.fontsize": 10,
"xtick.labelsize": 10, "xtick.labelsize": 10,
"ytick.labelsize": 10, "ytick.labelsize": 10,
"text.usetex": True, "text.usetex": True,
"figure.figsize": fig_size, "figure.figsize": fig_size,
"font.family": "serif", "font.family": "serif",
"font.serif": "Computer Modern Roman", "font.serif": "Computer Modern Roman",
"legend.frameon": True, "legend.frameon": True,
"savefig.dpi": 300, "savefig.dpi": 300,
} }
plt.rcParams.update(params) plt.rcParams.update(params)
plt.rc("text.latex", preamble=[r"\usepackage{xfrac}", r"\usepackage{siunitx}"]) plt.rc("text.latex", preamble=[r"\usepackage{xfrac}", r"\usepackage{siunitx}"])
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Fig 1. # Fig 1.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` ```
np.random.seed(1) np.random.seed(1)
xs = np.array([0.1, 0.3, 0.35, 0.45]) xs = np.array([0.1, 0.3, 0.35, 0.45])
f = lambda x: x ** 3 f = lambda x: x ** 3
ys = f(xs) ys = f(xs)
means = lambda x: np.convolve(x, np.ones(2) / 2, mode="valid") means = lambda x: np.convolve(x, np.ones(2) / 2, mode="valid")
xs_means = means(xs) xs_means = means(xs)
ys_means = means(ys) ys_means = means(ys)
fig, ax = plt.subplots(figsize=fig_size) fig, ax = plt.subplots(figsize=fig_size)
ax.scatter(xs, ys, c="k") ax.scatter(xs, ys, c="k")
ax.plot(xs, ys, c="k") ax.plot(xs, ys, c="k")
# ax.scatter() # ax.scatter()
ax.annotate( ax.annotate(
s=r"$L_{1,2} = \sqrt{\Delta x^2 + \Delta y^2}$", 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]])), xy=(np.mean([xs[0], xs[1]]), np.mean([ys[0], ys[1]])),
xytext=(xs[0] + 0.05, ys[0] - 0.05), xytext=(xs[0] + 0.05, ys[0] - 0.05),
arrowprops=dict(arrowstyle="->"), arrowprops=dict(arrowstyle="->"),
ha="center", ha="center",
zorder=10, zorder=10,
) )
for i, (x, y) in enumerate(zip(xs, ys)): for i, (x, y) in enumerate(zip(xs, ys)):
sign = [1, -1][i % 2] sign = [1, -1][i % 2]
ax.annotate( ax.annotate(
s=fr"$x_{i+1}, y_{i+1}$", s=fr"$x_{i+1}, y_{i+1}$",
xy=(x, y), xy=(x, y),
xytext=(x + 0.01, y + sign * 0.04), xytext=(x + 0.01, y + sign * 0.04),
arrowprops=dict(arrowstyle="->"), arrowprops=dict(arrowstyle="->"),
ha="center", ha="center",
) )
ax.scatter(xs, ys, c="green", s=5, zorder=5, label="existing data") ax.scatter(xs, ys, c="green", s=5, zorder=5, label="existing data")
losses = np.hypot(xs[1:] - xs[:-1], ys[1:] - ys[:-1]) losses = np.hypot(xs[1:] - xs[:-1], ys[1:] - ys[:-1])
ax.scatter( ax.scatter(
xs_means, ys_means, c="red", s=300 * losses, zorder=8, label="candidate points" xs_means, ys_means, c="red", s=300 * losses, zorder=8, label="candidate points"
) )
xs_dense = np.linspace(xs[0], xs[-1], 400) xs_dense = np.linspace(xs[0], xs[-1], 400)
ax.plot(xs_dense, f(xs_dense), alpha=0.3, zorder=7, label="function") ax.plot(xs_dense, f(xs_dense), alpha=0.3, zorder=7, label="function")
ax.legend() ax.legend()
ax.axis("off") ax.axis("off")
plt.savefig("figures/loss_1D.pdf", bbox_inches="tight", transparent=True) plt.savefig("figures/loss_1D.pdf", bbox_inches="tight", transparent=True)
plt.show() plt.show()
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Fig 2. # Fig 2.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` ```
import adaptive import adaptive
def f(x, offset=0.123): def f(x, offset=0.123):
a = 0.02 a = 0.02
return x + a ** 2 / (a ** 2 + (x - offset) ** 2) return x + a ** 2 / (a ** 2 + (x - offset) ** 2)
def g(x): def g(x):
return np.tanh(x * 40) return np.tanh(x * 40)
def h(x): def h(x):
return np.sin(100 * x) * np.exp(-x ** 2 / 0.1 ** 2) return np.sin(100 * x) * np.exp(-x ** 2 / 0.1 ** 2)
funcs = [ funcs = [
dict(function=f, bounds=(-1, 1), title="peak"), dict(function=f, bounds=(-1, 1), title="peak"),
dict(function=g, bounds=(-1, 1), title="tanh"), dict(function=g, bounds=(-1, 1), title="tanh"),
dict(function=h, bounds=(-0.3, 0.3), title="wave packet"), dict(function=h, bounds=(-0.3, 0.3), title="wave packet"),
] ]
fig, axs = plt.subplots(2, len(funcs), figsize=(fig_width, 1.5 * fig_height)) fig, axs = plt.subplots(2, len(funcs), figsize=(fig_width, 1.5 * fig_height))
n_points = 50 n_points = 50
for i, ax in enumerate(axs.T.flatten()): for i, ax in enumerate(axs.T.flatten()):
ax.xaxis.set_ticks([]) ax.xaxis.set_ticks([])
ax.yaxis.set_ticks([]) ax.yaxis.set_ticks([])
if i % 2 == 0: if i % 2 == 0:
d = funcs[i // 2] d = funcs[i // 2]
# homogeneous # homogeneous
xs = np.linspace(*d["bounds"], n_points) xs = np.linspace(*d["bounds"], n_points)
ys = d["function"](xs) ys = d["function"](xs)
ax.set_title(rf"\textrm{{{d['title']}}}") ax.set_title(rf"\textrm{{{d['title']}}}")
else: else:
d = funcs[(i - 1) // 2] d = funcs[(i - 1) // 2]
loss = adaptive.learner.learner1D.curvature_loss_function() loss = adaptive.learner.learner1D.curvature_loss_function()
learner = adaptive.Learner1D( learner = adaptive.Learner1D(
d["function"], bounds=d["bounds"], loss_per_interval=loss d["function"], bounds=d["bounds"], loss_per_interval=loss
) )
adaptive.runner.simple(learner, goal=lambda l: l.npoints >= n_points) adaptive.runner.simple(learner, goal=lambda l: l.npoints >= n_points)
# adaptive # adaptive
xs, ys = zip(*sorted(learner.data.items())) xs, ys = zip(*sorted(learner.data.items()))
xs_dense = np.linspace(*d["bounds"], 1000) xs_dense = np.linspace(*d["bounds"], 1000)
ax.plot(xs_dense, d["function"](xs_dense), c="red", alpha=0.3, lw=0.5) ax.plot(xs_dense, d["function"](xs_dense), c="red", alpha=0.3, lw=0.5)
ax.scatter(xs, ys, s=0.5, c="k") ax.scatter(xs, ys, s=0.5, c="k")
axs[0][0].set_ylabel(r"$\textrm{homogeneous}$") axs[0][0].set_ylabel(r"$\textrm{homogeneous}$")
axs[1][0].set_ylabel(r"$\textrm{adaptive}$") axs[1][0].set_ylabel(r"$\textrm{adaptive}$")
plt.savefig("figures/adaptive_vs_grid.pdf", bbox_inches="tight", transparent=True) plt.savefig("figures/adaptive_vs_grid.pdf", bbox_inches="tight", transparent=True)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Fig 3. # Fig 3.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` ```
from scipy import interpolate from scipy import interpolate
import functools import functools
import itertools import itertools
import adaptive import adaptive
import holoviews.plotting.mpl import holoviews.plotting.mpl
import matplotlib.tri as mtri import matplotlib.tri as mtri
def f(xy, offset=0.123): def f(xy, offset=0.123):
a = 0.1 a = 0.2
x, y = xy x, y = xy
return x * y + a ** 2 / (a ** 2 + (x - offset) ** 2 + (y - offset) ** 2) return x + np.exp(-(x ** 2 + y ** 2 - 0.75 ** 2) ** 2 / a ** 4)
@functools.lru_cache() @functools.lru_cache()
def g_setup(fname): def g_setup(fname):
data = adaptive.utils.load(fname) data = adaptive.utils.load(fname)
points = np.array(list(data.keys())) points = np.array(list(data.keys()))
values = np.array(list(data.values()), dtype=float) values = np.array(list(data.values()), dtype=float)
bounds = [(points[:, 0].min(), points[:, 0].max()), (points[:, 1].min(), points[:, 1].max())] bounds = [
(points[:, 0].min(), points[:, 0].max()),
(points[:, 1].min(), points[:, 1].max()),
]
ll, ur = np.reshape(bounds, (2, 2)).T ll, ur = np.reshape(bounds, (2, 2)).T
inds = np.all(np.logical_and(ll <= points, points <= ur), axis=1) inds = np.all(np.logical_and(ll <= points, points <= ur), axis=1)
points, values = points[inds], values[inds].reshape(-1, 1) points, values = points[inds], values[inds].reshape(-1, 1)
return interpolate.LinearNDInterpolator(points, values), bounds return interpolate.LinearNDInterpolator(points, values), bounds
def g(xy, fname): def g(xy, fname):
ip, _ = g_setup(fname) ip, _ = g_setup(fname)
return ip(xy) return np.round(ip(xy))
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 h(xy): def h(xy):
x, y = xy x, y = xy
return np.sin(100 * x * y) * np.exp(-x ** 2 / 0.1 ** 2 - y ** 2 / 0.4 ** 2) return density(x, y) + y
funcs = [ funcs = [
dict(function=f, bounds=[(-1, 1), (-1, 1)], title="peak", npoints=50,), dict(function=f, bounds=[(-1, 1), (-1, 1)], npoints=33),
dict( dict(
function=g, function=g,
bounds=g_setup("phase_diagram.pickle")[1], bounds=g_setup("phase_diagram.pickle")[1],
title="tanh", npoints=100,
npoints=140,
fname="phase_diagram.pickle", fname="phase_diagram.pickle",
), ),
dict( dict(function=h, bounds=[(-1, 1), (-3, 3)], npoints=50),
function=h,
bounds=[(-0.3, 0.3), (-0.3, 0.3)],
title="wave packet",
npoints=50,
),
] ]
fig, axs = plt.subplots(2, len(funcs), figsize=(fig_width, 1.5 * fig_height)) fig, axs = plt.subplots(2, len(funcs), figsize=(fig_width, 1.5 * fig_height))
plt.subplots_adjust(hspace=-0.1, wspace=0.1) plt.subplots_adjust(hspace=-0.1, wspace=0.1)
with_tri = False with_tri = False
for i, ax in enumerate(axs.T.flatten()): 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.xaxis.set_ticks([])
ax.yaxis.set_ticks([]) ax.yaxis.set_ticks([])
kind = "homogeneous" if i % 2 == 0 else "adaptive" kind = "homogeneous" if i % 2 == 0 else "adaptive"
d = funcs[i // 2] if kind == "homogeneous" else funcs[(i - 1) // 2] d = funcs[i // 2] if kind == "homogeneous" else funcs[(i - 1) // 2]
bounds = d["bounds"] bounds = d["bounds"]
npoints = d["npoints"] npoints = d["npoints"]
f = d["function"] f = d["function"]
fname = d.get("fname") fname = d.get("fname")
if fname is not None: if fname is not None:
f = functools.partial(f, fname=fname) f = functools.partial(f, fname=fname)
if kind == "homogeneous": if kind == "homogeneous":
ax.set_title(rf"\textrm{{{d['title']}}}")
xs, ys = [np.linspace(*bound, npoints) for bound in bounds] xs, ys = [np.linspace(*bound, npoints) for bound in bounds]
data = {xy: f(xy) for xy in itertools.product(xs, ys)} data = {xy: f(xy) for xy in itertools.product(xs, ys)}
learner = adaptive.Learner2D(f, bounds=bounds) learner = adaptive.Learner2D(f, bounds=bounds)
learner.data = data learner.data = data
d["learner_hom"] = learner
elif kind == "adaptive": elif kind == "adaptive":
learner = adaptive.Learner2D(f, bounds=bounds) learner = adaptive.Learner2D(f, bounds=bounds)
if fname is not None: if fname is not None:
learner.load(fname) 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) adaptive.runner.simple(learner, goal=lambda l: l.npoints >= npoints ** 2)
d["learner"] = learner
if with_tri: if with_tri:
tri = learner.ip().tri tri = learner.ip().tri
triang = mtri.Triangulation(*tri.points.T, triangles=tri.vertices) triang = mtri.Triangulation(*tri.points.T, triangles=tri.vertices)
ax.triplot(triang, c="w", lw=0.2, alpha=0.8) ax.triplot(triang, c="w", lw=0.2, alpha=0.8)
values = np.array(list(learner.data.values())) values = np.array(list(learner.data.values()))
ax.imshow(learner.plot().Image.I.data, extent=(-0.5, 0.5, -0.5, 0.5)) ax.imshow(
learner.plot(npoints if kind == "homogeneous" else None).Image.I.data,
extent=(-0.5, 0.5, -0.5, 0.5),
interpolation="none",
)
ax.set_xticks([]) ax.set_xticks([])
ax.set_yticks([]) ax.set_yticks([])
axs[0][0].set_ylabel(r"$\textrm{homogeneous}$") axs[0][0].set_ylabel(r"$\textrm{homogeneous}$")
axs[1][0].set_ylabel(r"$\textrm{adaptive}$") axs[1][0].set_ylabel(r"$\textrm{adaptive}$")
plt.savefig("figures/adaptive_2D.pdf", bbox_inches="tight", transparent=True) plt.savefig("figures/adaptive_2D.pdf", bbox_inches="tight", transparent=True)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` ```
from scipy import interpolate
import functools
import itertools
import adaptive
import holoviews.plotting.mpl
import matplotlib.tri as mtri
def f(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 g_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 g(xy, fname):
ip, _ = g_setup(fname)
return np.round(ip(xy))
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 h(xy):
x, y = xy
return density(x, y) + y
funcs = [
dict(function=f, bounds=[(-1, 1), (-1, 1)], npoints=33),
dict(
function=g,
bounds=g_setup("phase_diagram.pickle")[1],
npoints=100,
fname="phase_diagram.pickle",
),
dict(function=h, bounds=[(-1, 1), (-3, 3)], npoints=50),
]
fig, axs = plt.subplots(len(funcs), 2, figsize=(fig_width, 2 * fig_height))
plt.subplots_adjust(hspace=0.1, wspace=0.1)
with_tri = False
for i, ax in enumerate(axs.flatten()):
label = "abcdef"[i]
ax.text(
-0.03,
0.98,
f"$\mathrm{{({label})}}$",
transform=ax.transAxes,
horizontalalignment="right",
verticalalignment="top",
)
ax.xaxis.set_ticks([])
ax.yaxis.set_ticks([])
kind = "homogeneous" if i % 2 == 0 else "adaptive"
d = funcs[i // 2] if kind == "homogeneous" else funcs[(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,
extent=(-0.5, 0.5, -0.5, 0.5),
interpolation="none",
)
ax.set_xticks([])
ax.set_yticks([])
axs[0][0].set_title(r"$\textrm{homogeneous}$")
axs[0][1].set_title(r"$\textrm{adaptive}$")
plt.savefig("figures/adaptive_2D.pdf", bbox_inches="tight", transparent=True)
``` ```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment