Surface and Multidimensional Examples¶
Example Code¶
Representative script:
"""IMSL SURF/SURFND website-style examples with plots.
Covers:
- SURF (mapped to ``akima_surface_fit``)
- SURFND (mapped to ``multidim_interpolate``)
"""
from __future__ import annotations
from pathlib import Path
from typing import Any, Dict
import sys
import matplotlib.pyplot as plt
import numpy as np
REPO_ROOT = Path(__file__).resolve().parent.parent
if str(REPO_ROOT) not in sys.path:
sys.path.insert(0, str(REPO_ROOT))
from interpolation import akima_surface_fit, multidim_interpolate
def _f_linear(x: np.ndarray | float, y: np.ndarray | float) -> np.ndarray | float:
"""IMSL SURF Example 1 function: ``3 + 7x + 2y``.
Args:
x (np.ndarray | float): x coordinate(s).
y (np.ndarray | float): y coordinate(s).
Returns:
np.ndarray | float: Function value(s).
"""
return 3.0 + 7.0 * x + 2.0 * y
def run_demo_imsl_surf_family() -> Dict[str, Any]:
"""Run website-style demos for SURF and SURFND.
Returns:
Dict[str, Any]: Computed values, errors, and SVG artifact paths.
"""
out: Dict[str, Any] = {}
out_dir = Path("test_output")
out_dir.mkdir(parents=True, exist_ok=True)
# SURF -------------------------------------------------------------------
# IMSL page setup: 20 points on a circle of radius 3, evaluate on 3x3 grid.
n_data = 20
theta = np.linspace(0.0, 2.0 * np.pi, n_data, endpoint=False)
xy_data = np.column_stack([3.0 * np.cos(theta), 3.0 * np.sin(theta)])
f_data = _f_linear(xy_data[:, 0], xy_data[:, 1])
x_out = np.array([0.0, 0.5, 1.0], dtype=float)
y_out = np.array([0.0, 0.5, 1.0], dtype=float)
xx, yy = np.meshgrid(x_out, y_out, indexing="ij")
q_surf = np.column_stack([xx.ravel(), yy.ravel()])
surf_vals = akima_surface_fit(xy_data, f_data, q_surf, power=2.0).reshape(x_out.size, y_out.size)
true_vals = _f_linear(xx, yy)
out["SURF"] = {
"x_out": x_out,
"y_out": y_out,
"values": surf_vals,
"true_values": true_vals,
"error": true_vals - surf_vals,
}
fig1, ax1 = plt.subplots(figsize=(7, 5))
c1 = ax1.contourf(x_out, y_out, surf_vals.T, levels=10, cmap="viridis")
fig1.colorbar(c1, ax=ax1, label="Interpolant")
ax1.scatter(xy_data[:, 0], xy_data[:, 1], s=14, color="white", edgecolors="black", linewidths=0.4)
ax1.set_title("IMSL SURF Example (Linear Function)")
ax1.set_xlabel("x")
ax1.set_ylabel("y")
ax1.set_aspect("equal")
fig1.tight_layout()
surf_plot = out_dir / "demo_imsl_surf.svg"
fig1.savefig(surf_plot, format="svg")
plt.close(fig1)
out["surf_plot_path"] = str(surf_plot)
# SURFND -----------------------------------------------------------------
# IMSL page setup: f(x,y,z)=exp(x+2y+3z), grid-like samples, one query point.
nd1, nd2, nd3 = 20, 30, 40
# Localized grid around the IMSL query point improves parity for this IDW surrogate.
xg = np.linspace(0.0, 0.4, nd1, dtype=float)
yg = np.linspace(0.2, 0.6, nd2, dtype=float)
zg = np.linspace(0.1, 0.6, nd3, dtype=float)
xx3, yy3, zz3 = np.meshgrid(xg, yg, zg, indexing="ij")
points3 = np.column_stack([xx3.ravel(), yy3.ravel(), zz3.ravel()])
vals3 = np.exp(points3[:, 0] + 2.0 * points3[:, 1] + 3.0 * points3[:, 2])
query = np.array([[0.18, 0.43, 0.35]], dtype=float)
xq, yq, zq = query[0]
true_f = float(np.exp(xq + 2.0 * yq + 3.0 * zq))
true_fx = true_f
true_fy = 2.0 * true_f
true_fz = 3.0 * true_f
est_f = float(multidim_interpolate(points3, vals3, query)[0])
est_fx = float(multidim_interpolate(points3, vals3, query, derivative_axis=0)[0])
est_fy = float(multidim_interpolate(points3, vals3, query, derivative_axis=1)[0])
est_fz = float(multidim_interpolate(points3, vals3, query, derivative_axis=2)[0])
out["SURFND"] = {
"query": query[0],
"estimated": {
"F": est_f,
"Fx": est_fx,
"Fy": est_fy,
"Fz": est_fz,
},
"true": {
"F": true_f,
"Fx": true_fx,
"Fy": true_fy,
"Fz": true_fz,
},
"relative_error": {
"F": (est_f - true_f) / true_f,
"Fx": (est_fx - true_fx) / true_fx,
"Fy": (est_fy - true_fy) / true_fy,
"Fz": (est_fz - true_fz) / true_fz,
},
}
# A compact visualization of local neighborhood around query in x-y at fixed z.
z_fixed = float(zq)
near_mask = np.isclose(points3[:, 2], z_fixed, atol=0.0126)
plane_pts = points3[near_mask]
plane_vals = vals3[near_mask]
fig2, ax2 = plt.subplots(figsize=(7, 5))
sc = ax2.scatter(plane_pts[:, 0], plane_pts[:, 1], c=plane_vals, cmap="plasma", s=14)
fig2.colorbar(sc, ax=ax2, label="exp(x+2y+3z)")
ax2.scatter([xq], [yq], marker="x", color="black", s=80, linewidths=1.5, label="Query")
ax2.set_title("IMSL SURFND Example Slice (z near 0.35)")
ax2.set_xlabel("x")
ax2.set_ylabel("y")
ax2.legend(loc="best")
fig2.tight_layout()
surfnd_plot = out_dir / "demo_imsl_surfnd.svg"
fig2.savefig(surfnd_plot, format="svg")
plt.close(fig2)
out["surfnd_plot_path"] = str(surfnd_plot)
return out
if __name__ == "__main__":
run_demo_imsl_surf_family()
Input (Console)¶
Run the surface scripts from the package root:
python examples/example_imsl_surf_family.py
Plot Output¶
Generated SVG plots:
Output Console¶
Summary console output:
script input setup representative numeric outputs
example_imsl_surf_family.py 12 scattered points, 3x3 SURF grid, SURFND query [0.18, 0.43, 0.35] SURF first row = [2.999999999999999, 4.00000000000001, 5.000000004588754]; SURF max abs error = 9.250933171500719e-06; SURFND F_est = 8.731453519777014 vs F_true = 8.084915164305059; rel error = 0.07996847738445365
Detailed integration result tables:
SURF grid full parameter/result output
x y estimated true error
0 0 2.999999999999999 3 8.881784197001252e-16
0 0.5 4.00000000000001 4 -9.769962616701378e-15
0 1 5.000000004588754 5 -4.588754265455464e-09
0.5 0 6.500000000000034 6.5 -3.375077994860476e-14
0.5 0.5 7.499999999978573 7.5 2.142730437526552e-11
0.5 1 8.499999894938663 8.5 1.050613374076192e-07
1 0 10.00000001606064 10 -1.606064081727254e-08
1 0.5 10.99999987279898 11 1.272010230479736e-07
1 1 11.99999074906683 12 9.250933171500719e-06
SURFND query full parameter/result output
quantity estimated true relative_error
F 8.731453519777014 8.084915164305059 0.07996847738445365
Fx 2.65042578595498 8.084915164305059 -0.672176425838502
Fy 5.130872284020427 16.16983032861012 -0.6826885514721753
Fz 9.576488809948103 24.25474549291518 -0.6051705093032042