import time # wall-clock timing
from scipy.interpolate import make_interp_spline, RegularGridInterpolator, RectBivariateSpline
def ns_per_point(evaluate, n_points, repeats=5):
# the best of `repeats` timings of one vectorized call, in nanoseconds per point
best = np.inf # fastest run so far
for _ in range(repeats):
start = time.perf_counter() # high-resolution clock
evaluate() # evaluate at all points at once
best = min(best, time.perf_counter() - start) # keep the fastest run
return 1e9 * best / n_points # seconds per call → ns per point
def timing_1d(n_knots=1000, n_points=1_000_000):
# evaluation time per point in 1D, at random (unsorted) points
rng = np.random.default_rng(0) # reproducible points
x = np.linspace(0.0, 1.0, n_knots) # uniform grid
y = np.sin(8.0 * x) + 0.3 * x # smooth data on it
p = rng.uniform(0.0, 1.0, n_points) # the evaluation points
convinterp, scipy = {}, {} # ns per point, by method
for k in ["a1", "a3", "a4", "b5", "b7", "b9", "b11", "b13"]:
itp = convolution_interpolation(x, y, kernel=k) # built once, outside the timing
convinterp[f'"{k}"'] = ns_per_point(lambda: itp(p), n_points)
scipy["np.interp (linear)"] = ns_per_point(lambda: np.interp(p, x, y), n_points)
spline = CubicSpline(x, y) # cubic spline
scipy["CubicSpline"] = ns_per_point(lambda: spline(p), n_points)
for degree in (5, 7):
bspline = make_interp_spline(x, y, k=degree) # interpolating B-spline
scipy[f"B-spline, degree {degree}"] = ns_per_point(lambda: bspline(p), n_points)
return convinterp, scipy
def timing_2d(n_knots=200, n_points=200_000):
# evaluation time per point in 2D, at random points on a 200 × 200 grid
rng = np.random.default_rng(1) # reproducible points
x = np.linspace(0.0, 1.0, n_knots) # knots of axis 0
y = np.linspace(0.0, 1.0, n_knots) # knots of axis 1
data = np.sin(8.0 * x)[:, None] * np.cos(5.0 * y)[None, :] # smooth data on the grid
px = rng.uniform(0.0, 1.0, n_points) # x coordinates of the points
py = rng.uniform(0.0, 1.0, n_points) # y coordinates of the points
pxy = np.column_stack([px, py]) # the same points as rows, for SciPy
convinterp, scipy = {}, {} # ns per point, by method
for k in ["a1", "a3", "a4", "b5", "b7"]:
itp = convolution_interpolation((x, y), data, kernel=k) # built once, outside the timing
convinterp[f'"{k}"'] = ns_per_point(lambda: itp(px, py), n_points)
for method in ["linear", "cubic", "quintic"]:
rgi = RegularGridInterpolator((x, y), data, method=method) # SciPy's grid interpolator
scipy[f"RegularGridInterpolator, {method}"] = ns_per_point(lambda: rgi(pxy), n_points)
for degree in (3, 5):
rbs = RectBivariateSpline(x, y, data, kx=degree, ky=degree) # tensor-product spline
scipy[f"RectBivariateSpline, degree {degree}"] = ns_per_point(
lambda: rbs(px, py, grid=False), n_points) # grid=False: evaluate point by point
return convinterp, scipy
def bar_panel(ax, convinterp, scipy, title):
# horizontal bars: convinterp kernels on top, SciPy methods below, with the times written out
labels = [f"convinterp {k}" for k in convinterp] + list(scipy) # one bar per method
times = list(convinterp.values()) + list(scipy.values()) # ns per point, same order
colors = ["tab:blue"] * len(convinterp) + ["tab:orange"] * len(scipy)
rows = np.arange(len(labels)) # one row per bar
ax.barh(rows, times, color=colors)
ax.set_yticks(rows, labels)
ax.invert_yaxis() # first entry at the top
for row, t in zip(rows, times):
ax.text(t, row, f" {t:.0f}", va="center", fontsize=9) # the time next to each bar
ax.set_xlim(0, 1.15 * max(times)) # room for the labels
ax.set_xlabel("Evaluation time per point [ns]")
ax.set_title(title)
def timing_figure():
c1, s1 = timing_1d() # 1D timings
c2, s2 = timing_2d() # 2D timings
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 12)) # 1D and 2D, stacked
bar_panel(ax1, c1, s1, "1D: 1000 knots, 10⁶ random points")
bar_panel(ax2, c2, s2, "2D: 200 × 200 knots, 2·10⁵ random points")
fig.tight_layout()
return fig
fig = timing_figure()
plt.show()
