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

use fig from orbital field paper

parent 8e841729
No related branches found
No related tags found
No related merge requests found
Pipeline #21226 passed
%% Cell type:code id: tags:
```
import numpy as np
import matplotlib
matplotlib.use("agg")
import matplotlib.pyplot as plt
%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}"])
```
%% Cell type:markdown id: tags:
# Fig 1.
%% 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.scatter()
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:
# Fig 2.
%% Cell type:code id: tags:
```
import adaptive
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 = [
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"),
]
fig, axs = plt.subplots(2, len(funcs), 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([])
if i % 2 == 0:
d = funcs[i // 2]
# homogeneous
xs = np.linspace(*d["bounds"], n_points)
ys = d["function"](xs)
ax.set_title(rf"\textrm{{{d['title']}}}")
else:
d = funcs[(i - 1) // 2]
loss = adaptive.learner.learner1D.curvature_loss_function()
learner = adaptive.Learner1D(
d["function"], bounds=d["bounds"], loss_per_interval=loss
)
adaptive.runner.simple(learner, goal=lambda l: l.npoints >= n_points)
# adaptive
xs, ys = zip(*sorted(learner.data.items()))
xs_dense = np.linspace(*d["bounds"], 1000)
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")
axs[0][0].set_ylabel(r"$\textrm{homogeneous}$")
axs[1][0].set_ylabel(r"$\textrm{adaptive}$")
plt.savefig("figures/adaptive_vs_grid.pdf", bbox_inches="tight", transparent=True)
```
%% Cell type:markdown id: tags:
# Fig 3.
%% 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.1
x, y = xy
return x * y + a ** 2 / (a ** 2 + (x - offset) ** 2 + (y - offset) ** 2)
def g(xy):
x, y = xy
return np.tanh(x * 40) * np.tanh(y * 40)
@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 ip(xy)
def h(xy):
x, y = xy
return np.sin(100 * x * y) * np.exp(-x ** 2 / 0.1 ** 2 - y ** 2 / 0.4 ** 2)
funcs = [
dict(function=f, bounds=[(-1, 1), (-1, 1)], title="peak"),
dict(function=g, bounds=[(-1, 1), (-1, 1)], title="tanh"),
dict(function=h, bounds=[(-0.3, 0.3), (-0.3, 0.3)], title="wave packet"),
dict(function=f, bounds=[(-1, 1), (-1, 1)], title="peak", npoints=50,),
dict(
function=g,
bounds=g_setup("phase_diagram.pickle")[1],
title="tanh",
npoints=140,
fname="phase_diagram.pickle",
),
dict(
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))
plt.subplots_adjust(hspace=-0.1, wspace=0.1)
n_points = 50
with_tri = False
for i, ax in enumerate(axs.T.flatten()):
ax.xaxis.set_ticks([])
ax.yaxis.set_ticks([])
if i % 2 == 0:
d = funcs[i // 2]
# homogeneous
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":
ax.set_title(rf"\textrm{{{d['title']}}}")
x, y = d["bounds"]
xs = np.linspace(*x, n_points)
ys = np.linspace(*y, n_points)
data = {xy: d["function"](xy) for xy in itertools.product(xs, ys)}
learner = adaptive.Learner2D(d["function"], bounds=d["bounds"])
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
else:
# adaptive
d = funcs[(i - 1) // 2]
learner = adaptive.Learner2D(d["function"], bounds=d["bounds"])
adaptive.runner.simple(learner, goal=lambda l: l.npoints >= n_points ** 2)
tri = learner.ip().tri
triang = mtri.Triangulation(*tri.points.T, triangles=tri.vertices)
# ax.triplot(triang, c="w", lw=0.2, alpha=0.8)
elif kind == "adaptive":
learner = adaptive.Learner2D(f, bounds=bounds)
if fname is not None:
learner.load(fname)
adaptive.runner.simple(learner, goal=lambda l: l.npoints >= npoints ** 2)
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().Image.I.data, extent=(-0.5, 0.5, -0.5, 0.5))
ax.set_xticks([])
ax.set_yticks([])
axs[0][0].set_ylabel(r"$\textrm{homogeneous}$")
axs[1][0].set_ylabel(r"$\textrm{adaptive}$")
plt.savefig("figures/adaptive_2D.pdf", bbox_inches="tight", transparent=True)
```
%% Cell type:code id: tags:
```
```
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment