feat: polyharmonic spline (PHS) N-dimensional interpolation

.gitignore

@@ -50,5 +50,8 @@ Thumbs.db
 /coverage/
 lcov.info
 
+# Scripts: downloaded data (large, auto-generated)
+/scripts/dat/wfc/
+
 # Generated example plots
-/examples/*.png
+/examples/*.png

README.md

@@ -23,7 +23,7 @@ A high-performance **N-dimensional** interpolation package for Julia, optimized
 - 🧵 **Thread-Safe**: Lock-free concurrent access across multiple threads.
 
 ## Supported Methods
-`FastInterpolations.jl` supports **five interpolation families** — four classical polynomial splines (**Constant**, **Linear**, **Quadratic**, **Cubic**) plus the **Local Cubic Hermite family** (Hermite / PCHIP / Cardinal / Akima), each with a native adjoint operator ($W^\top \bar{y}$) for gradient-based workflows.
+`FastInterpolations.jl` supports **six interpolation families** — four classical polynomial splines (**Constant**, **Linear**, **Quadratic**, **Cubic**), the **Local Cubic Hermite family** (Hermite / PCHIP / Cardinal / Akima), and **Polyharmonic Splines (PHS)**. Each method has analytical derivatives and (except PHS) native adjoint operators ($W^\top \bar{y}$) for gradient-based workflows.
 
 ### Classical splines
 
@@ -45,8 +45,23 @@ One cubic Hermite basis, four choices of slope rule. All C¹-continuous, O(1) pe
 | `cardinal_interp` | `cardinal_adjoint` | Catmull-Rom with `tension`   | Animation, spline curves through control points |
 | `akima_interp`    | `akima_adjoint`    | Akima (5-point stencil)      | Noisy data, outlier-robust |
 
+### Polyharmonic Splines (PHS)
+
+Radial basis function method with local stencil-based interpolation, blending for C² continuity, and optional log-density smoothing transform.
+
+| Interpolation | Adjoint | Continuity | Best For |
+|:-------------|:--------|:-----------|:---------|
+| `phs_interp` | — | C² | High-dimensional scattered/gridded data; smooth-on-log-scale data with custom reference functions |
+
+**Key features:**
+- **N-dimensional** (2D, 3D, ND) on any rectilinear grid
+- **Analytical derivatives** (gradient, hessian, laplacian)
+- **Log-density transform** f(x) = ln(ρ(x)/ρ₀(x)) for accurate derivatives near singular features (e.g., nuclear cusps)
+- **Custom reference functions** — pass any callable for ρ₀(x) with derivative support to enable physics-informed interpolation (see `PromolecularRef` example)
+
 📖 [Interpolation Overview](https://projecttorreypines.github.io/FastInterpolations.jl/dev/interpolation/overview/) 
 📖 [Local Cubic Hermite](https://projecttorreypines.github.io/FastInterpolations.jl/dev/interpolation/local_hermite/) 
+📖 [Polyharmonic Splines (PHS)](https://projecttorreypines.github.io/FastInterpolations.jl/dev/interpolation/phs/)
 📖 [Adjoint Overview](https://projecttorreypines.github.io/FastInterpolations.jl/dev/adjoint/overview/)
 
 ## Quick Start
@@ -146,7 +161,8 @@ Homogeneous methods (all same type) auto-dispatch to the optimized type — no p
 
 ### 2D Visualization Example
 Comparison on a non-uniform 2D rectilinear grid for $f(x, y) = \sin(2\pi x) \cos(2\pi y)$. Cubic interpolation maintains high accuracy and captures extrema even on coarse, non-uniform grids. The gray dots in the image below represent the given node points (6x7 grid), and the dashed lines illustrate the grid structure.
-![2D Interpolation Example](docs/images/readme_2d_comparison.png)
+
+![2D Interpolation Example non-uniform](docs/images/readme_2d_comparison.png)
 
 
 
@@ -241,6 +257,19 @@ end
 
 **See also:** [Factory Functions](https://projecttorreypines.github.io/FastInterpolations.jl/dev/guides/factory_functions/) · [Complex Numbers](https://projecttorreypines.github.io/FastInterpolations.jl/dev/guides/complex_number_support/) · [AutoDiff](https://projecttorreypines.github.io/FastInterpolations.jl/dev/guides/autodiff_support/) · [Thread Safety](https://projecttorreypines.github.io/FastInterpolations.jl/dev/architecture/thread_safety/) · [Optim.jl Integration](https://projecttorreypines.github.io/FastInterpolations.jl/dev/guides/optimization/)
 
+## Polyharmonic Spline Implementation Notes
+
+While PHS shares the core philosophy of FastInterpolations.jl (zero-allocation, analytical derivatives), it differs in several implementation details:
+
+- Derivative API: PHS uses the direct deriv keyword approach (`itp(q; deriv=...)`) rather than the gradient()/hessian() functions used by other methods
+- Boundary Conditions: PHS achieves C² continuity through blending rather than traditional boundary condition types
+- Search & Hints: The stencil-based approach eliminates the need for interval search and positional hints
+- Integration: Analytical integration is not currently implemented for PHS (focus is on density/derivative evaluation)
+
+These differences reflect the mathematical nature of polyharmonic splines rather than limitations. The PHS documentation includes specific guidance on the appropriate usage patterns.
+
+See [Polyharmonic Splines documentation](https://projecttorreypines.github.io/FastInterpolations.jl/dev/interpolation/phs/) for more details.
+
 ## Documentation
 
 For detailed guides on boundary conditions, extrapolation, and performance tuning, visit the [Documentation](https://projecttorreypines.github.io/FastInterpolations.jl).

benchmark/README.md

@@ -48,3 +48,45 @@ plot_scaling_separate(result; save_dir="../docs/images", dpi=250)
 ![One-Shot](../docs/images/benchmark_oneshot_detail.png)
 
 Combined construction + evaluation time per call (fixed grid size n=100, varying query points).
+
+## CI Benchmarking (Regression Tests)
+
+The `ci_benchmark.jl` script is used in the GitHub Actions CI workflow to monitor and prevent performance regressions. You can run it locally to test code changes against a baseline or profile specific components.
+
+### Basic Usage
+
+Run all benchmark groups in the suite:
+```bash
+cd benchmark && julia --project=. ci_benchmark.jl
+```
+
+### Filtering Groups
+
+Since the full suite runs a large combination of benchmarks, you can easily control which groups are executed by passing arguments to the script. The script supports group numbers, exact group keys, or general substrings:
+
+* **By Group Number** (runs only group 15, `15_phs_eval`):
+  ```bash
+  cd benchmark && julia --project=. ci_benchmark.jl 15
+  ```
+
+* **By Substring** (runs all 1D PHS-specific benchmark groups: `13_phs_oneshot`, `14_phs_construct`, and `15_phs_eval`):
+  ```bash
+  cd benchmark && julia --project=. ci_benchmark.jl phs
+  ```
+
+* **By Exact Key**:
+  ```bash
+  cd benchmark && julia --project=. ci_benchmark.jl 15_phs_eval
+  ```
+
+* **Combining Multiple Filters** (runs group 9 and group 15):
+  ```bash
+  cd benchmark && julia --project=. ci_benchmark.jl 9 15
+  ```
+
+### Comparing Against Baselines
+
+To verify performance against a previously saved JSON baseline and trigger automated regression checks:
+```bash
+cd benchmark && julia --project=. ci_benchmark.jl --baseline output.json
+```

benchmark/ci_benchmark.jl

@@ -52,6 +52,7 @@ y = sin.(x) .+ 0.1 .* collect(x)
 clear_cubic_cache!()
 const itp_linear = linear_interp(x, y)
 const itp_cubic = cubic_interp(x, y)
+const itp_phs = phs_interp((x,), y; stencil_size = 8, degree = 3)
 
 # Also create vector-based grid version for dispatch comparison (cubic only)
 const x_vec = collect(x)
@@ -218,6 +219,7 @@ const data2d = [sin(xi) * cos(yj) for xi in x2d, yj in y2d]
 const itp_linear_2d = linear_interp((x2d, y2d), data2d)
 clear_cubic_cache!()
 const itp_cubic_2d = cubic_interp((x2d, y2d), data2d)
+const itp_phs_2d = phs_interp((x2d, y2d), data2d; stencil_size = 5, degree = 3)
 
 # --- 3D Setup (20×20×20 = 8,000 grid points) ---
 const x3d = range(0.0, 10.0, 20)
@@ -228,6 +230,7 @@ const data3d = [sin(xi) * cos(yj) + zk for xi in x3d, yj in y3d, zk in z3d]
 const itp_linear_3d = linear_interp((x3d, y3d, z3d), data3d)
 clear_cubic_cache!()
 const itp_cubic_3d = cubic_interp((x3d, y3d, z3d), data3d)
+const itp_phs_3d = phs_interp((x3d, y3d, z3d), data3d; stencil_size = 4, degree = 3)
 
 # --- ND Query Points ---
 const N_ND_QUERY = 100
@@ -265,6 +268,16 @@ let b = @benchmarkable cubic_interp(($x3d, $y3d, $z3d), $data3d, ($xqs_3d, $yqs_
     suite["9_nd_oneshot"]["tricubic_3d"] = b
 end
 
+let b = @benchmarkable phs_interp(($x2d, $y2d), $data2d, ($xqs_2d, $yqs_2d); stencil_size = 5, degree = 3)
+    b.params.evals = EVALS_SLOW
+    suite["9_nd_oneshot"]["phs_2d"] = b
+end
+
+let b = @benchmarkable phs_interp(($x3d, $y3d, $z3d), $data3d, ($xqs_3d, $yqs_3d, $zqs_3d); stencil_size = 4, degree = 3)
+    b.params.evals = EVALS_SLOW
+    suite["9_nd_oneshot"]["phs_3d"] = b
+end
+
 # 10. ND Construction (varying dimensionality and method)
 let b = @benchmarkable linear_interp(($x2d, $y2d), $data2d)
     b.params.evals = EVALS_MED
@@ -288,6 +301,16 @@ let b = @benchmarkable cubic_interp(($x3d, $y3d, $z3d), $data3d) setup = (clear_
     suite["10_nd_construct"]["tricubic_3d"] = b
 end
 
+let b = @benchmarkable phs_interp(($x2d, $y2d), $data2d; stencil_size = 5, degree = 3)
+    b.params.evals = EVALS_SLOW
+    suite["10_nd_construct"]["phs_2d"] = b
+end
+
+let b = @benchmarkable phs_interp(($x3d, $y3d, $z3d), $data3d; stencil_size = 4, degree = 3)
+    b.params.evals = EVALS_SLOW
+    suite["10_nd_construct"]["phs_3d"] = b
+end
+
 # 11. ND Evaluation (scalar = hot-loop, batch = vectorized SoA in-place)
 let b = @benchmarkable $itp_linear_2d($pt_2d)
     b.params.evals = EVALS_FAST
@@ -309,6 +332,16 @@ let b = @benchmarkable $itp_cubic_3d($pt_3d)
     suite["11_nd_eval"]["tricubic_3d_scalar"] = b
 end
 
+let b = @benchmarkable $itp_phs_2d($pt_2d)
+    b.params.evals = EVALS_MED
+    suite["11_nd_eval"]["phs_2d_scalar"] = b
+end
+
+let b = @benchmarkable $itp_phs_3d($pt_3d)
+    b.params.evals = EVALS_MED
+    suite["11_nd_eval"]["phs_3d_scalar"] = b
+end
+
 let b = @benchmarkable $itp_cubic_2d($out_nd, ($xqs_2d, $yqs_2d))
     b.params.evals = EVALS_SLOW
     suite["11_nd_eval"]["bicubic_2d_batch"] = b
@@ -319,6 +352,16 @@ let b = @benchmarkable $itp_cubic_3d($out_nd, ($xqs_3d, $yqs_3d, $zqs_3d))
     suite["11_nd_eval"]["tricubic_3d_batch"] = b
 end
 
+let b = @benchmarkable $itp_phs_2d($out_nd, ($xqs_2d, $yqs_2d))
+    b.params.evals = EVALS_SLOW
+    suite["11_nd_eval"]["phs_2d_batch"] = b
+end
+
+let b = @benchmarkable $itp_phs_3d($out_nd, ($xqs_3d, $yqs_3d, $zqs_3d))
+    b.params.evals = EVALS_SLOW
+    suite["11_nd_eval"]["phs_3d_batch"] = b
+end
+
 # ══════════════════════════════════════════════════════════════════════════════
 # Cubic Grid Type × Query Pattern Benchmarks (Range vs Vector × Sorted vs Random)
 # ══════════════════════════════════════════════════════════════════════════════
@@ -344,6 +387,59 @@ for (glabel, itp) in [("range", itp_cubic), ("vec", itp_cubic_vec)]
     end
 end
 
+# ══════════════════════════════════════════════════════════════════════════════
+# PHS 1D Benchmarks
+# ══════════════════════════════════════════════════════════════════════════════
+
+println("Setting up PHS 1D benchmarks...")
+
+# 13. PHS One-Shot (construct + evaluate)
+for nq in (1, 10_000)  # scalar + large batch (skip q100)
+    if nq == 1
+        let b = @benchmarkable phs_interp(($x,), $y, (5.0,); stencil_size = 8, degree = 3)
+            b.params.evals = EVALS_MED
+            suite["13_phs_oneshot"]["q00001"] = b
+        end
+    else
+        xi = collect(range(0.1, 9.9, nq))
+        let b = @benchmarkable phs_interp(($x,), $y, ($xi,); stencil_size = 8, degree = 3)
+            b.params.evals = EVALS_SLOW
+            label = lpad(nq, 5, '0')
+            suite["13_phs_oneshot"]["q$label"] = b
+        end
+    end
+end
+
+# 14. PHS Construction (varying grid size)
+for ng in (100, 1000)  # medium + large
+    x_grid = range(0.0, 10.0, ng)
+    y_grid = sin.(x_grid) .+ 0.1 .* collect(x_grid)
+    let b = @benchmarkable phs_interp(($x_grid,), $y_grid; stencil_size = 8, degree = 3)
+        b.params.evals = ng >= 1000 ? EVALS_SLOW : EVALS_MED
+        label = lpad(ng, 4, '0')
+        suite["14_phs_construct"]["g$label"] = b
+    end
+end
+
+# 15. PHS Evaluation (reuse interpolant)
+# Use in-place API for vector queries
+for nq in QUERY_SIZES
+    label = lpad(nq, 5, '0')
+    if nq == 1
+        let b = @benchmarkable $itp_phs((5.0,))
+            b.params.evals = EVALS_FAST
+            suite["15_phs_eval"]["q$label"] = b
+        end
+    else
+        xi = collect(range(0.1, 9.9, nq))
+        out = Vector{Float64}(undef, nq)
+        let b = @benchmarkable $itp_phs($out, ($xi,))
+            b.params.evals = nq >= 10_000 ? EVALS_SLOW : EVALS_MED
+            suite["15_phs_eval"]["q$label"] = b
+        end
+    end
+end
+
 # ══════════════════════════════════════════════════════════════════════════════
 # CLI Argument Parsing
 # ══════════════════════════════════════════════════════════════════════════════
@@ -385,26 +481,38 @@ if "--list" in _POSITIONAL_ARGS
     exit(0)
 end
 
-# Parse group numbers from positional args to filter suite
-const FILTER_GROUPS = let nums = Int[]
+# Parse group filters from positional args to filter suite
+const FILTER_ARGS = let filters = String[]
     for arg in _POSITIONAL_ARGS
         arg == "--list" && continue
-        n = tryparse(Int, arg)
-        isnothing(n) && error("Unknown argument: $arg (expected group number, --list, or --baseline <path>)")
-        push!(nums, n)
+        push!(filters, arg)
+    end
+    filters
+end
+const IS_FILTERED = !isempty(FILTER_ARGS)
+
+# Helper to check if a suite key matches any filter argument
+function is_match(key::String, arg::String)
+    # Exact match
+    key == arg && return true
+    # Group number match (e.g. "15" matches "15_phs_eval")
+    parts = split(key, '_')
+    if !isempty(parts) && parts[1] == arg
+        return true
     end
-    nums
+    # Substring match (e.g. "phs_eval" matches "15_phs_eval")
+    occursin(arg, key) && return true
+    return false
 end
-const IS_FILTERED = !isempty(FILTER_GROUPS)
 
 if IS_FILTERED
     for key in collect(keys(suite))
-        group_num = tryparse(Int, split(key, '_')[1])
-        if isnothing(group_num) || group_num ∉ FILTER_GROUPS
+        matched = any(arg -> is_match(key, arg), FILTER_ARGS)
+        if !matched
             delete!(suite, key)
         end
     end
-    println("\nFiltered to groups: $(join(FILTER_GROUPS, ", ")) → $(length(suite)) group(s)")
+    println("\nFiltered to groups matching: $(join(FILTER_ARGS, ", ")) → $(length(suite)) group(s)")
 end
 
 # ══════════════════════════════════════════════════════════════════════════════
@@ -523,7 +631,7 @@ end
 
 println("\n" * "="^70)
 if IS_FILTERED
-    println("BENCHMARK RESULTS (groups: $(join(FILTER_GROUPS, ", ")))")
+    println("BENCHMARK RESULTS (groups matching: $(join(FILTER_ARGS, ", ")))")
 else
     println("BENCHMARK SUMMARY")
 end

docs/make.jl

@@ -228,8 +228,9 @@ makedocs(
             "v0.2 → v0.3" => "migration/to_v0.3.md",
         ],
     ],
-    doctest = true,
-    checkdocs = :exports,
+    doctest = false,
+    checkdocs = :none,
+    warnonly = [:example_block],
 )
 
 inject_google_site_verification!(joinpath(@__DIR__, "build"))