From 7b45c5f994507f3bdd009eb696d843c5d8cf9dae Mon Sep 17 00:00:00 2001 From: saudzahirr Date: Mon, 18 May 2026 00:00:57 +0500 Subject: [PATCH] Refactor --- README.md | 2 +- docs/quickstart.md | 253 +++++++++++++++++++++++++++----------- docs/usage/derivatives.md | 106 +++++++++++++++- docs/usage/integration.md | 41 +++++- docs/usage/performance.md | 80 +++++++++++- pyproject.toml | 3 - 6 files changed, 393 insertions(+), 92 deletions(-) diff --git a/README.md b/README.md index c5c9bab..88ba9d8 100644 --- a/README.md +++ b/README.md @@ -67,7 +67,7 @@ Requires gfortran, Meson, and f2py (`numpy`). - [Theory](https://eggzec.github.io/spinterp/theory/) — hierarchical basis, Smolyak construction, algorithms - [Quickstart](https://eggzec.github.io/spinterp/quickstart/) — runnable examples -- [API Reference](https://eggzec.github.io/spinterp/api/) — function signatures and arguments +- [API Reference](https://eggzec.github.io/spinterp/usage/) — function signatures and arguments ## License diff --git a/docs/quickstart.md b/docs/quickstart.md index 7b0fa7f..89ace91 100644 --- a/docs/quickstart.md +++ b/docs/quickstart.md @@ -10,111 +10,214 @@ f(x, y) = \sin(x) + \cos(y) over the domain $[0, \pi] \times [0, \pi]$ using the default **Clenshaw-Curtis** sparse grid. -### Step 1 — build the hierarchical surpluses - -The sparse grid interpolant is assembled level by level. At each level $k$, new grid points -$\mathbf{x}_k$ are generated, the function is evaluated there, and the hierarchical surplus -(the correction to the interpolant from the previous level) is computed. +The complete script below builds the interpolant level by level, evaluates it at random +points, computes gradient vectors, integrates over the domain, visualises the grid and +the error surface, and demonstrates alternative grid types: ```python import numpy as np +import matplotlib.pyplot as plt import spinterp -d = 2 # dimension -scale = np.array([np.pi, np.pi]) # domain [0,pi]^2 +# ============================================================ +# Problem setup +# ============================================================ + +d = 2 +max_level = 4 +scale = np.array([np.pi, np.pi]) # domain: [0, pi]^2 def f(x, y): return np.sin(x) + np.cos(y) -all_seq, all_surp = [], [] +# ============================================================ +# STEP 1 — Build hierarchical surpluses (Clenshaw-Curtis) +# ============================================================ + +all_seq = [] +all_surp = [] + +for k in range(max_level + 1): -for k in range(5): # levels 0 .. 4 nl = spinterp.spnlevels(k, d) - seq = spinterp.spgetseq(k, d, nl) # multi-index set, shape (nl, d) - tp = spinterp.spdim_cc(seq) - x_k = spinterp.spgrid_cc(seq, tp) # grid points in [0,1]^d, shape (tp, d) + seq = spinterp.spgetseq(k, d, nl) # multi-index set, shape (nl, d) + tp = spinterp.spdim_cc(seq) # total grid points at this level + x_k = spinterp.spgrid_cc(seq, tp) # sparse-grid points in [0,1]^d - # Evaluate f at the scaled grid points - fvals = np.array([f(*(x_k[i] * scale)) for i in range(tp)]) + x_scaled = x_k * scale + fvals = np.array([ + f(x_scaled[i, 0], x_scaled[i, 1]) for i in range(tp) + ]) if k == 0: surp_k = fvals.copy() else: - z_prev = np.concatenate(all_surp) - seq_prev = np.vstack(all_seq) - interp = spinterp.spcmpvals_cc(z_prev, x_k, seq, seq_prev) - surp_k = fvals - interp + z_prev = np.concatenate(all_surp) + seq_prev = np.vstack(all_seq) + interp_prev = spinterp.spcmpvals_cc(z_prev, x_k, seq, seq_prev) + surp_k = fvals - interp_prev all_seq.append(seq) all_surp.append(surp_k) z = np.concatenate(all_surp) # flat surplus array -seq = np.vstack(all_seq) # combined level-index matrix -``` +seq = np.vstack(all_seq) # combined level-index matrix -### Step 2 — evaluate the interpolant +print("Finished building sparse-grid interpolant.") +print("Total surpluses:", len(z)) -```python -rng = np.random.default_rng(42) -pts = rng.random((5, d)) * scale # 5 random points in [0,pi]^2 -pts_unit = pts / scale # normalise to [0,1]^2 - -ip = spinterp.spinterp_cc(z, pts_unit, seq) -exact = f(pts[:, 0], pts[:, 1]) +# ============================================================ +# STEP 2 — Evaluate interpolant at random points +# ============================================================ -print("Interpolated:", ip) -print("Exact: ", exact) -print("Max error: ", np.max(np.abs(ip - exact))) -``` +rng = np.random.default_rng(42) +pts = rng.random((10, d)) * scale # random points in [0, pi]^2 +pts_unit = pts / scale # normalise to [0,1]^2 -Example output: +ip = spinterp.spinterp_cc(z, pts_unit, seq) +exact = f(pts[:, 0], pts[:, 1]) +err = np.abs(ip - exact) + +print("\nInterpolation Test") +print("-" * 50) +for i in range(len(pts)): + print(f"Point {i+1}") + print(f" x = {pts[i]}") + print(f" interpolated = {ip[i]: .8f}") + print(f" exact = {exact[i]: .8f}") + print(f" error = {err[i]: .3e}") +print("\nMax error:", np.max(err)) + +# ============================================================ +# STEP 3 — Compute gradient +# ============================================================ + +ip2, grad = spinterp.spderiv_cc(z, pts_unit, seq) +# grad.shape == (npoints, d) +print("\nGradient Test") +print("-" * 50) +for i in range(len(pts)): + print(f"Point {i+1}") + print(f" df/dx = {grad[i, 0]: .8f}") + print(f" df/dy = {grad[i, 1]: .8f}") + +# ============================================================ +# STEP 4 — Quadrature / Integration over [0, pi]^2 +# ============================================================ + +tp_all = sum(spinterp.spdim_cc(s) for s in all_seq) +weights = spinterp.spquadw_cc(seq, tp_all) +# spquadw_cc integrates over [0,1]^d; multiply by the domain volume to +# recover the integral over [0,pi]^2 +integral = float(np.dot(weights, z)) * np.prod(scale) + +# Exact: ∫∫ [sin(x)+cos(y)] dx dy over [0,pi]^2 = 2*pi +exact_integral = 2 * np.pi + +print("\nQuadrature Test") +print("-" * 50) +print("Approx integral =", integral) +print("Exact integral =", exact_integral) +print("Absolute error =", abs(integral - exact_integral)) + +# ============================================================ +# STEP 5 — Visualise sparse-grid points +# ============================================================ + +tp_plot = spinterp.spdim_cc(seq) +grid_plot = spinterp.spgrid_cc(seq, tp_plot) +grid_scaled = grid_plot * scale + +plt.figure(figsize=(6, 6)) +plt.scatter(grid_scaled[:, 0], grid_scaled[:, 1], s=30) +plt.xlabel("x") +plt.ylabel("y") +plt.title(f"Clenshaw-Curtis Sparse Grid (level={max_level}, d={d})") +plt.grid(True) + +# ============================================================ +# STEP 6 — Compare exact function vs interpolant on a fine grid +# ============================================================ + +N = 80 +xv = np.linspace(0, np.pi, N) +yv = np.linspace(0, np.pi, N) +X, Y = np.meshgrid(xv, yv) + +Z_exact = np.sin(X) + np.cos(Y) + +pts_eval = np.column_stack([X.ravel(), Y.ravel()]) +pts_eval_unit = pts_eval / scale +Z_interp = spinterp.spinterp_cc(z, pts_eval_unit, seq).reshape(N, N) + +fig = plt.figure(figsize=(10, 7)) +ax = fig.add_subplot(111, projection='3d') +ax.plot_surface(X, Y, Z_exact) +ax.set_title("Exact Function: sin(x) + cos(y)") +ax.set_xlabel("x") +ax.set_ylabel("y") + +fig2 = plt.figure(figsize=(10, 7)) +ax2 = fig2.add_subplot(111, projection='3d') +ax2.plot_surface(X, Y, Z_interp) +ax2.set_title("Sparse Grid Interpolant") +ax2.set_xlabel("x") +ax2.set_ylabel("y") + +fig3 = plt.figure(figsize=(10, 7)) +ax3 = fig3.add_subplot(111, projection='3d') +ax3.plot_surface(X, Y, np.abs(Z_exact - Z_interp)) +ax3.set_title("Absolute Interpolation Error") +ax3.set_xlabel("x") +ax3.set_ylabel("y") + +plt.show() + +# ============================================================ +# OPTIONAL — Alternative grid types +# ============================================================ + +print("\nAlternative grid examples") + +# Chebyshev polynomial sparse grid — same point count as CC +tp_cb = spinterp.spdim_cc(seq) +x_cb = spinterp.spgrid_cb(seq, tp_cb) +interp_cb = spinterp.spinterp_cb(z, pts_unit, seq) +print("Chebyshev interpolant computed.") + +# Gauss-Patterson sparse grid — same point count as NoBoundary +tp_gp = spinterp.spdim_m(seq, boundary=0) +x_gp = spinterp.spgrid_gp(seq, tp_gp) +interp_gp = spinterp.spinterp_gp(z, pts_unit, seq) +print("Gauss-Patterson interpolant computed.") + +print("\nDone.") ``` -Interpolated: [0.641 1.765 0.278 0.832 1.113] -Exact: [0.641 1.765 0.278 0.832 1.113] -Max error: 0.0047 -``` - -### Step 3 — visualise the sparse grid - -The sparse grid at level 4 in 2-D with Clenshaw-Curtis nodes: - -![Clenshaw-Curtis sparse grid, level 4, d=2](_static/ex_firstexample_01.png) - -And comparing the true function with the sparse grid interpolant: - -![f(x,y)=sin(x)+cos(y) vs sparse grid interpolant](_static/ex_firstexample_02.png) --- ## Grid types -Choose a different grid by swapping the `spgrid_*`, `spcmpvals_*`, and `spinterp_*` calls: +Choose a different grid by swapping the `spgrid_*`, `spcmpvals_*`, and `spinterp_*` calls +and using the matching point-count function: -```python -# Chebyshev polynomial sparse grid -x_k = spinterp.spgrid_cb(seq, tp) -interp = spinterp.spcmpvals_cb(z_prev, x_k, seq, seq_prev) -ip = spinterp.spinterp_cb(z, pts_unit, seq) - -# Gauss-Patterson -x_k = spinterp.spgrid_gp(seq, tp) -interp = spinterp.spcmpvals_gp(z_prev, x_k, seq, seq_prev) -ip = spinterp.spinterp_gp(z, pts_unit, seq) -``` +| Grid type | Point-count function | `spgrid_*` | `spinterp_*` | +|---|---|---|---| +| Clenshaw-Curtis | `spdim_cc(seq)` | `spgrid_cc` | `spinterp_cc` | +| Chebyshev | `spdim_cc(seq)` | `spgrid_cb` | `spinterp_cb` | +| Gauss-Patterson | `spdim_m(seq, boundary=0)` | `spgrid_gp` | `spinterp_gp` | +| Maximum | `spdim_m(seq, boundary=1)` | `spgrid_m` | `spinterp_m` | +| NoBoundary | `spdim_m(seq, boundary=0)` | `spgrid_nb` | `spinterp_nb` | + +See the **OPTIONAL** section of the script above for Chebyshev and Gauss-Patterson examples. --- ## Derivatives -Request gradient vectors by calling `spderiv_cc` instead of `spinterp_cc`: - -```python -ip, grad = spinterp.spderiv_cc(z, pts_unit, seq) -# grad.shape == (npoints, d) -print("df/dx1:", grad[:, 0]) -print("df/dx2:", grad[:, 1]) -``` +Call `spderiv_cc` instead of `spinterp_cc` to get both the interpolated value and the full +gradient simultaneously — see **STEP 3** of the script above. See [Computing Derivatives](usage/derivatives.md) for the full discussion. @@ -122,11 +225,13 @@ See [Computing Derivatives](usage/derivatives.md) for the full discussion. ## Quadrature -Integrate $f$ over $[0,1]^d$ by dotting the quadrature weights with the surpluses: +Integrate $f$ over the domain by dotting the quadrature weights with the surpluses — see +**STEP 4** of the script above. -```python -tp_all = sum(spinterp.spdim_cc(s) for s in all_seq) -w = spinterp.spquadw_cc(seq, tp_all) -integral = float(np.dot(w, z)) -print("Integral ≈", integral) -``` +Clenshaw-Curtis sparse grid at level 4: + +![Clenshaw-Curtis sparse grid, level 4, d=2](_static/ex_firstexample_01.png) + +Exact function vs interpolant: + +![f(x,y)=sin(x)+cos(y) vs sparse grid interpolant](_static/ex_firstexample_02.png) diff --git a/docs/usage/derivatives.md b/docs/usage/derivatives.md index d312a3b..895da63 100644 --- a/docs/usage/derivatives.md +++ b/docs/usage/derivatives.md @@ -12,11 +12,42 @@ Call `spderiv_cc` (or `spderiv_cb` for the Chebyshev grid) instead of `spinterp_ function returns both the interpolated values **and** the full gradient vector: ```python +import numpy as np import spinterp +d = 2 + +def f(x, y): + return np.sin(x) + np.cos(y) + +# Build hierarchical surpluses (Clenshaw-Curtis grid) +all_seq, all_surp = [], [] +for k in range(5): + nl = spinterp.spnlevels(k, d) + seq = spinterp.spgetseq(k, d, nl) + tp = spinterp.spdim_cc(seq) + x_k = spinterp.spgrid_cc(seq, tp) + fvals = np.array([f(*x_k[i]) for i in range(tp)]) + if k == 0: + surp_k = fvals.copy() + else: + z_prev = np.concatenate(all_surp) + seq_prev = np.vstack(all_seq) + surp_k = fvals - spinterp.spcmpvals_cc(z_prev, x_k, seq, seq_prev) + all_seq.append(seq) + all_surp.append(surp_k) + +z = np.concatenate(all_surp) +seq = np.vstack(all_seq) + +pts_unit = np.array([[0.25, 0.4], [0.6, 0.75], [0.1, 0.9]]) + # ip.shape = (npoints,) # grad.shape = (npoints, d) ip, grad = spinterp.spderiv_cc(z, pts_unit, seq) +print("Values: ", ip) +print("df/dx1: ", grad[:, 0]) +print("df/dx2: ", grad[:, 1]) ``` The procedure for building the surpluses `z` and level-index matrix `seq` is identical @@ -83,17 +114,50 @@ where $\Delta$ is the cell width $1/2^{\ell_{\max}}$. Use `spcont_deriv_cc` followed by `pp_deriv` to obtain continuous derivatives: ```python -maxlev = int(seq[:, 0].max()) -ip, ipder, ipder2 = spinterp.spcont_deriv_cc(z, pts_unit, seq, maxlev) -# pp_deriv post-processes ipder in-place using ipder2 import numpy as np +import spinterp + +d = 2 + +def f(x, y): + return np.sin(x) + np.cos(y) + +# Build hierarchical surpluses (Clenshaw-Curtis grid) +all_seq, all_surp = [], [] +for k in range(5): + nl = spinterp.spnlevels(k, d) + seq = spinterp.spgetseq(k, d, nl) + tp = spinterp.spdim_cc(seq) + x_k = spinterp.spgrid_cc(seq, tp) + fvals = np.array([f(*x_k[i]) for i in range(tp)]) + if k == 0: + surp_k = fvals.copy() + else: + z_prev = np.concatenate(all_surp) + seq_prev = np.vstack(all_seq) + surp_k = fvals - spinterp.spcmpvals_cc(z_prev, x_k, seq, seq_prev) + all_seq.append(seq) + all_surp.append(surp_k) + +z = np.concatenate(all_surp) +seq = np.vstack(all_seq) + +pts_unit = np.array([[0.25, 0.4], [0.6, 0.75], [0.1, 0.9]]) +maxlev = int(seq.max()) maxlevvec = np.full(d, maxlev, dtype=np.int32) + +ip, ipder, ipder2 = spinterp.spcont_deriv_cc(z, pts_unit, seq, maxlev) + +# pp_deriv post-processes ipder in-place using ipder2 spinterp.pp_deriv( np.asfortranarray(ipder), np.asfortranarray(ipder2), maxlevvec, - np.asfortranarray(pts_unit) + np.asfortranarray(pts_unit), ) +print("Values: ", ip) +print("df/dx1: ", ipder[:, 0]) +print("df/dx2: ", ipder[:, 1]) ``` The figure below shows the same derivative after the continuity post-processing: @@ -109,7 +173,41 @@ The derivatives are computed via the **discrete cosine transform (DCT)**, using `spderiv_cb` function: ```python +import numpy as np +import spinterp + +d = 2 + +def f(x, y): + return np.sin(x) + np.cos(y) + +# Build hierarchical surpluses (Chebyshev grid) +# CB uses the same point count as CC +all_seq, all_surp = [], [] +for k in range(5): + nl = spinterp.spnlevels(k, d) + seq = spinterp.spgetseq(k, d, nl) + tp = spinterp.spdim_cc(seq) + x_k = spinterp.spgrid_cb(seq, tp) + fvals = np.array([f(*x_k[i]) for i in range(tp)]) + if k == 0: + surp_k = fvals.copy() + else: + z_prev = np.concatenate(all_surp) + seq_prev = np.vstack(all_seq) + surp_k = fvals - spinterp.spcmpvals_cb(z_prev, x_k, seq, seq_prev) + all_seq.append(seq) + all_surp.append(surp_k) + +z_cb = np.concatenate(all_surp) +seq = np.vstack(all_seq) + +pts_unit = np.array([[0.25, 0.4], [0.6, 0.75], [0.1, 0.9]]) + ip, grad = spinterp.spderiv_cb(z_cb, pts_unit, seq) +print("Values: ", ip) +print("df/dx1: ", grad[:, 0]) +print("df/dx2: ", grad[:, 1]) ``` The resulting derivatives are infinitely smooth: diff --git a/docs/usage/integration.md b/docs/usage/integration.md index 6dca92b..8f6d244 100644 --- a/docs/usage/integration.md +++ b/docs/usage/integration.md @@ -16,13 +16,46 @@ where $\mathbf{w}$ is the vector of quadrature weights and $\mathbf{z}$ is the f surplus array. ```python -import spinterp, numpy as np +import numpy as np +import spinterp + +d = 2 +max_level = 5 + +def f(x, y): + return np.sin(x) + np.cos(y) + +# Build hierarchical surpluses (Clenshaw-Curtis grid) +all_seq, all_surp = [], [] +for k in range(max_level + 1): + nl = spinterp.spnlevels(k, d) + seq = spinterp.spgetseq(k, d, nl) + tp = spinterp.spdim_cc(seq) + x_k = spinterp.spgrid_cc(seq, tp) + fvals = np.array([f(*x_k[i]) for i in range(tp)]) + if k == 0: + surp_k = fvals.copy() + else: + z_prev = np.concatenate(all_surp) + seq_prev = np.vstack(all_seq) + surp_k = fvals - spinterp.spcmpvals_cc(z_prev, x_k, seq, seq_prev) + all_seq.append(seq) + all_surp.append(surp_k) + +z = np.concatenate(all_surp) +seq = np.vstack(all_seq) # Build quadrature weights for the CC grid -# (seq and z come from the standard surplus construction) -tp_all = sum(spinterp.spdim_cc(s) for s in all_seq) -w = spinterp.spquadw_cc(seq, tp_all) +# tp_all must match the total size of the combined seq / z +tp_all = sum(spinterp.spdim_cc(s) for s in all_seq) +w = spinterp.spquadw_cc(seq, tp_all) integral = float(np.dot(w, z)) + +# Exact: ∫∫ [sin(x)+cos(y)] dx dy over [0,1]^2 = (1-cos(1)) + sin(1) +exact = (1 - np.cos(1)) + np.sin(1) +print(f"Approximate: {integral:.10f}") +print(f"Exact: {exact:.10f}") +print(f"Error: {abs(integral - exact):.3e}") ``` --- diff --git a/docs/usage/performance.md b/docs/usage/performance.md index ba26d0c..14f61af 100644 --- a/docs/usage/performance.md +++ b/docs/usage/performance.md @@ -50,7 +50,11 @@ The recommended pattern is to build the interpolant in a loop, checking the esti after each level: ```python -import numpy as np, spinterp +import numpy as np +import spinterp + +def f_vec(x1, x2): + return (x1 * x2) ** 2 d, n_max = 2, 8 all_seq, all_surp = [], [] @@ -103,14 +107,44 @@ A simple Python equivalent: keep only subgrids where $|\alpha_{\mathbf{l}}|_\inf a threshold: ```python +import numpy as np +import spinterp + +d = 2 +max_level = 6 + +def f(x, y): + return np.sin(x) + np.cos(y) + +# Build hierarchical surpluses (Clenshaw-Curtis grid) +all_seq, all_surp = [], [] +for k in range(max_level + 1): + nl = spinterp.spnlevels(k, d) + seq = spinterp.spgetseq(k, d, nl) + tp = spinterp.spdim_cc(seq) + x_k = spinterp.spgrid_cc(seq, tp) + fvals = np.array([f(*x_k[i]) for i in range(tp)]) + if k == 0: + surp_k = fvals.copy() + else: + z_prev = np.concatenate(all_surp) + seq_prev = np.vstack(all_seq) + surp_k = fvals - spinterp.spcmpvals_cc(z_prev, x_k, seq, seq_prev) + all_seq.append(seq) + all_surp.append(surp_k) + +# Purge subgrids whose surpluses are all below the drop tolerance drop_tol = 1e-5 keep = [ - (seq, surp) - for seq, surp in zip(all_seq, all_surp) + (s, surp) + for s, surp in zip(all_seq, all_surp) if np.max(np.abs(surp)) > drop_tol ] -seq_pruned = np.vstack([s for s, _ in keep]) +seq_pruned = np.vstack([s for s, _ in keep]) z_pruned = np.concatenate([surp for _, surp in keep]) + +print(f"Original subgrids: {len(all_seq)}, after purge: {len(keep)}") +print(f"Original surpluses: {sum(len(s) for s in all_surp)}, after purge: {len(z_pruned)}") ``` The trade-off: smaller `drop_tol` → higher accuracy, slower evaluation. Typical savings @@ -128,11 +162,45 @@ all query points as a single 2-D array is far more efficient than one-point-at-a calls. ```python +import numpy as np +import spinterp + +d = 2 +max_level = 4 + +def f(x, y): + return np.sin(x) + np.cos(y) + +# Build hierarchical surpluses (Clenshaw-Curtis grid) +all_seq, all_surp = [], [] +for k in range(max_level + 1): + nl = spinterp.spnlevels(k, d) + seq = spinterp.spgetseq(k, d, nl) + tp = spinterp.spdim_cc(seq) + x_k = spinterp.spgrid_cc(seq, tp) + fvals = np.array([f(*x_k[i]) for i in range(tp)]) + if k == 0: + surp_k = fvals.copy() + else: + z_prev = np.concatenate(all_surp) + seq_prev = np.vstack(all_seq) + surp_k = fvals - spinterp.spcmpvals_cc(z_prev, x_k, seq, seq_prev) + all_seq.append(seq) + all_surp.append(surp_k) + +z = np.concatenate(all_surp) +seq = np.vstack(all_seq) + +rng = np.random.default_rng(0) +pts = rng.random((1000, d)) + # Slow — Python loop over 1000 points -results = [spinterp.spinterp_cc(z, pts[i:i+1], seq)[0] for i in range(1000)] +results_slow = [spinterp.spinterp_cc(z, pts[i:i+1], seq)[0] for i in range(1000)] # Fast — single batched call -results = spinterp.spinterp_cc(z, pts, seq) # pts.shape = (1000, d) +results_fast = spinterp.spinterp_cc(z, pts, seq) # pts.shape = (1000, d) + +print(f"Max diff slow vs fast: {np.max(np.abs(np.array(results_slow) - results_fast)):.2e}") ``` Typical speed-up: **50–100×** for 1000 points, due to reduced Python overhead and better diff --git a/pyproject.toml b/pyproject.toml index 9b73eea..94e0bb6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -30,9 +30,6 @@ keywords = [ "Chebyshev", "hierarchical surplus", "high-dimensional approximation", - "Fortran", - "f2py", - "scientific computing", "numerical methods", ] classifiers = [