Error and Validation Examples¶
Example Code¶
Representative script:
"""IMSL SPLEZ example: ITYPE sweep with max errors and plot."""
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 (
bspline_derivative,
bspline_interpolate,
bspline_optimal_knots,
cubic_spline_akima,
cubic_spline_derivative,
cubic_spline_evaluate,
cubic_spline_not_a_knot,
cubic_spline_shape_preserving,
least_squares_bspline_fixed_knots,
least_squares_bspline_variable_knots,
)
def _eval_pp(pp, xq: np.ndarray, der: int) -> np.ndarray:
if der == 0:
return np.asarray([cubic_spline_evaluate(pp, float(xv)) for xv in xq], dtype=float)
return np.asarray([cubic_spline_derivative(pp, float(xv), k=1) for xv in xq], dtype=float)
def _eval_bs(rep, xq: np.ndarray, der: int) -> np.ndarray:
if der == 0:
return np.asarray([rep_eval for rep_eval in [bspline_interpolate(np.array([0.0,1.0]), np.array([0.0,1.0]), order=2)]], dtype=float)
return np.asarray([0.0], dtype=float)
def run_demo_imsl_splez() -> Dict[str, Any]:
"""Run SPLEZ-like ITYPE sweep and return max-error table."""
ndata = 21
n = 2 * ndata - 1
x_data = np.linspace(0.0, 3.0, ndata)
f_data = np.sin(x_data**2)
x_vec = np.linspace(0.0, 3.0, n)
f_true = np.sin(x_vec**2)
fp_true = 2.0 * x_vec * np.cos(x_vec**2)
# ITYPE mapping using available API subset.
emax_f = np.zeros(15, dtype=float)
emax_fp = np.zeros(15, dtype=float)
# Precompute representatives.
pp_csint = cubic_spline_not_a_knot(x_data, f_data)
pp_akima = cubic_spline_akima(x_data, f_data)
pp_shape = cubic_spline_shape_preserving(x_data, f_data)
bs_by_order = {}
for order in [2, 3, 4, 5, 6]:
bs_by_order[order] = bspline_interpolate(x_data, f_data, order=order)
# Least squares reps
ls_fixed = {}
for order in [2, 3, 4]:
knots = bspline_optimal_knots(x_data, order=order, n_coeff=ndata)
ls_fixed[order] = least_squares_bspline_fixed_knots(x_data, f_data, knots, order)
ls_var = {}
for order in [2, 3, 4]:
ls_var[order] = least_squares_bspline_variable_knots(x_data, f_data, order=order, n_coeff=ndata)
for itype in range(1, 16):
if itype == 1:
y = _eval_pp(pp_csint, x_vec, 0)
yp = _eval_pp(pp_csint, x_vec, 1)
elif itype == 2:
y = _eval_pp(pp_akima, x_vec, 0)
yp = _eval_pp(pp_akima, x_vec, 1)
elif itype == 3:
y = _eval_pp(pp_shape, x_vec, 0)
yp = _eval_pp(pp_shape, x_vec, 1)
elif itype in [4, 5, 6, 7, 8]:
order = itype - 2
rep = bs_by_order[order]
y = np.asarray([bspline_derivative(rep, float(xv), k=0) for xv in x_vec], dtype=float)
yp = np.asarray([bspline_derivative(rep, float(xv), k=1) for xv in x_vec], dtype=float)
elif itype == 9:
y = _eval_pp(pp_shape, x_vec, 0)
yp = _eval_pp(pp_shape, x_vec, 1)
elif itype in [10, 11, 12]:
order = itype - 8
rep = ls_fixed[order]
y = np.asarray([bspline_derivative(rep, float(xv), k=0) for xv in x_vec], dtype=float)
yp = np.asarray([bspline_derivative(rep, float(xv), k=1) for xv in x_vec], dtype=float)
else:
order = itype - 11
rep = ls_var[order]
y = np.asarray([bspline_derivative(rep, float(xv), k=0) for xv in x_vec], dtype=float)
yp = np.asarray([bspline_derivative(rep, float(xv), k=1) for xv in x_vec], dtype=float)
emax_f[itype - 1] = float(np.max(np.abs(f_true - y)))
emax_fp[itype - 1] = float(np.max(np.abs(fp_true - yp)))
out_dir = Path("test_output")
out_dir.mkdir(parents=True, exist_ok=True)
out_path = out_dir / "demo_imsl_splez_max_errors.svg"
itypes = np.arange(1, 16, dtype=int)
fig, ax = plt.subplots(figsize=(10, 5))
ax.plot(itypes, emax_f, marker="o", linewidth=1.8, label="Max error for f")
ax.plot(itypes, emax_fp, marker="s", linewidth=1.8, label="Max error for f'")
ax.set_title("IMSL SPLEZ ITYPE Sweep (Mapped)")
ax.set_xlabel("ITYPE")
ax.set_ylabel("Maximum Error")
ax.grid(True, alpha=0.3)
ax.legend(loc="best")
fig.tight_layout()
fig.savefig(out_path, format="svg")
plt.close(fig)
return {
"itype": itypes,
"max_error_f": emax_f,
"max_error_fp": emax_fp,
"plot_path": str(out_path),
}
if __name__ == "__main__":
run_demo_imsl_splez()
Input (Console)¶
Run the validation script from the package root:
python examples/example_imsl_splez.py
Plot Output¶
Generated SVG plot:
Output Console¶
Summary console output:
script input setup representative numeric outputs
example_imsl_splez.py 15 spline type cases over 61 evaluation points max_error_f head = [0.009888476542265101, 0.02468177029579066, 0.05860109234948407, 0.08361467534116784, 0.00775744963901559]; max_error_fp head = [0.4952542339509511, 0.8977610263391025, 1.552541606680284, 1.786471229161286, 0.4262806226986982]; overall maxima = 0.08361467534116784 and 1.786471229161286
Detailed integration result tables:
SPLEZ full per-type error output
itype max_error_f max_error_fp
1 0.009888476542265101 0.4952542339509511
2 0.02468177029579066 0.8977610263391025
3 0.05860109234948407 1.552541606680284
4 0.08361467534116784 1.786471229161286
5 0.00775744963901559 0.4262806226986982
6 0.002924592999776543 0.05109599902921835
7 0.0004553465734867812 0.01035561382465056
8 0.00478972791596588 0.3503967644460637
9 0.05860109234948407 1.552541606680284
10 0.08361467534116784 1.786471229161286
11 0.007757449639016034 0.4262806226986982
12 0.002924592999779319 0.0510959990847295
13 0.08361467534116784 1.786471229161286
14 0.007757449639016034 0.4262806226986982
15 0.002924592999779319 0.0510959990847295