Release notes

ConvolutionInterpolations.jl

The notes of every release, as published on GitHub, newest first. Each release is archived on Zenodo under its own DOI; the citation file in the repository has the details.

v0.28.2 — 2026-10-14

ConvolutionInterpolations v0.28.2

Diff since v0.28.1

ConvolutionInterpolations.jl v0.28.2

Adds a CITATION.cff file so the package can be cited, with each release archived on Zenodo under a DOI. No code changes.

Breaking changes

None.

Migration guide

No changes needed.

v0.28.1 — 2026-09-29

ConvolutionInterpolations v0.28.1

Diff since v0.28.0

ConvolutionInterpolations.jl v0.28.1

Patch release: Relaxes the Project compat field to be less restrictive in the patch version. Furthermore, the readme has been updated with the correct new integral orders.

v0.28.0 — 2026-09-28

ConvolutionInterpolations v0.28.0

Diff since v0.27.0

ConvolutionInterpolations.jl v0.28.0

Highlights

Scattered data are now smoothed by default. fit_scattered and convolution_interpolation(points, values) now use mode = :smooth unless told otherwise. Scattered data are usually measurements, and even a little noise is harmful under exact interpolation: the high-order kernels reproduce the noise and amplify it between the points, most of all in derivatives. With smoothing as the default, noisy data give sound results without any tuning, while noise-free data still come out close to interpolation. Trusted, noise-free data can opt into exact interpolation with mode = :exact.

Faster, more reliable smoothing defaults. :smooth now fits with the :a4 kernel and roughness order q = 5 for every kernel. With noisy data the accuracy is set by the noise rather than the kernel, and :a4’s narrow stencil makes the fit 2–3× faster than before in 2D and 3–9× faster in 3D, at comparable accuracy. Wider kernels such as :b7 also become reliable in smoothing mode.

Clearer feedback. :smooth warns when the grid limit set by oversample stops the fit from resolving the chosen smoothing. It also raises a clear error when there are too few data points to smooth, which means 5^D or fewer: 25 in 2D, 125 in 3D.

Breaking changes

• The default mode of fit_scattered and of convolution_interpolation(points, values) is now :smooth instead of :exact. Calls without mode now smooth the data rather than interpolate it.
• With the new default, the third value returned by fit_scattered is the named tuple (lambda, noise) instead of the constraint violation viol.
• :smooth defaults to the :a4 kernel and q = 5, so results of mode = :smooth without an explicit kernel or q differ slightly from v0.27.0.
• :smooth raises an error for 5^D or fewer data points.

Migration guide

• For exact interpolation of trusted data, add mode = :exact: convolution_interpolation(points, values) becomes convolution_interpolation(points, values; mode = :exact), and the same applies to fit_scattered. Results are then identical to v0.27.0.
• Code that used the third return value of fit_scattered as viol should pass mode = :exact.
• For noisy data, no change is needed: the default now does the right thing. To reproduce v0.27.0’s smoothing exactly, pass kernel = :b5.
• For very small data sets, use mode = :exact, or pass a lower q.

v0.27.0 — 2026-09-27

ConvolutionInterpolations v0.27.0

Diff since v0.26.0

ConvolutionInterpolations.jl v0.27.0

Highlights

Smoothing of noisy scattered data: mode = :smooth

fit_scattered and convolution_interpolation(points, values) gain mode = :smooth for noisy scattered data. Instead of passing through every point, the fit balances closeness to the data against smoothness, and chooses that balance itself by generalized cross-validation, so there is no smoothing parameter to tune. The result is a high-order interpolant like any other: evaluate it, differentiate it, integrate it.

itp = convolution_interpolation(points, values; mode=:smooth)    # smoothing chosen automatically
knots, vals, info = fit_scattered(points, values; mode=:smooth)  # the fit, with diagnostics
info.noise                                                       # estimated noise level of the data
info.lambda                                                      # the chosen smoothing; pass it back to reproduce the fit

Why it matters. The previous options for noisy scattered data each needed hand-tuning and fell short. The scattered_to_grid + convolution_smooth pipeline left the imprint of the nearest-neighbour gridding in the result and blurred fine features. mode = :lsq fitted the noise at its default settings and only worked once its grid density was tuned by hand for each noise level and signal. In my comparisons, :smooth with its automatic choice matched or beat both of them even at their best hand-tuned settings, with lower errors in the values and derivatives. It also reports an estimate of the data’s noise level, a useful check against a known measurement error. It stays fast: the fit’s grid is refined only as far as the chosen smoothing needs.

Breaking changes

• fit_scattered(...; mode = :lsq) has been removed and now raises an ArgumentError. mode = :smooth replaces it.
• With mode = :smooth, the third value returned by fit_scattered is a named tuple (lambda, noise) instead of a number. mode = :exact still returns the constraint violation.

Deprecations

• scattered_to_grid is deprecated and will be removed in v1.0. It still works, with a one-time warning.

Migration guide

Least-squares fits of noisy scattered data:

# before
knots, vals, _ = fit_scattered(points, values; mode=:lsq, gridsize=14)
# after: the grid and the smoothing are chosen automatically
knots, vals, info = fit_scattered(points, values; mode=:smooth)

The gridding-and-smoothing pipeline:

# before
grid = scattered_to_grid(points, values, (xs, ys))
itp = convolution_interpolation((xs, ys), convolution_smooth((xs, ys), grid, 0.1))
# after: fit and smooth in one step
itp = convolution_interpolation(points, values; mode=:smooth)

The smoothed interpolant is defined on the data’s bounding box rather than on a grid you choose, so evaluate inside that box, or pass extrap. For noisy data that is already on a grid, convolution_smooth remains the tool.

v0.26.0 — 2026-09-27

ConvolutionInterpolations v0.26.0

Diff since v0.25.0

ConvolutionInterpolations.jl v0.26.0

Highlights

• Lazy mode gives exactly eager’s boundary values, in every dimension. Lazy boundary cells now compute the same ghost points as eager construction, including the :detect boundary check. Lazy mode supports real boundary ghosts and derivatives in any dimension, and is thread-safe.
• Lazy mode is chosen automatically for large grids. With the new default lazy=nothing, a uniform grid is built lazily when eager expansion would be costly in memory (more than 2²⁷ coefficients) or construction time (more than 2²⁴ ghost points), and eagerly otherwise.
• More accurate boundaries. Ghost points are computed with compensated summation directly from the data, removing a rounding error that grew on N-D edges and corners.
• Faster, simpler antiderivatives. One evaluator handles every integral order, with exactly derived weights; first-order integrals in 2D and 3D evaluate up to 2× faster, and integral construction is much faster.

Bug fixes

• convolution_interpolation(knots, values) without keywords on nonuniform knots built a uniform interpolant; it now builds the nonuniform one, as the keyword call does.
• Grids with different knot types per axis (e.g. a range and a vector) failed on the nonuniform b-kernel path.
• Lazy interpolants mixing :a0/:a1 axes with higher-order kernels failed at boundary cells.

Breaking changes

• boundary_fallback is deprecated and has no effect (a warning is shown; it will be removed in v1.0). Lazy interpolants that used boundary_fallback=true now return exact values near the boundaries instead of linear ones.
• The default for lazy is now nothing (chosen by grid size, see above) instead of false. Large grids are built lazily unless lazy=false is given; keyword-free calls in 4D are no longer lazy for small grids.
• Lazy interpolants with bc=:detect now follow eager’s boundary decisions, so their values near the boundaries can differ from v0.25.
• Values near the boundaries and antiderivative values change at rounding level.
• First-order antiderivative interpolants display as derivative=-1 instead of integral.

Migration guide

• Remove boundary_fallback from your calls.
• To keep eager construction regardless of grid size, pass lazy=false; to force lazy mode, pass lazy=true.
• Code that accessed internal fields of the interpolant types (such as left_values, the integral tail arrays, lazy_workspace or boundary_fallback) must be updated; these fields were removed.

v0.25.0 — 2026-09-24

ConvolutionInterpolations v0.25.0

Diff since v0.24.0

ConvolutionInterpolations.jl v0.25.0

Highlights

• Integrals of any order. derivative = -m gives the m-fold antiderivative of the interpolant, anchored at zero at the leftmost interior knot, up to order 2 (:a0, :a1), 4 (other :a kernels), 6 (:b5) or 8 (:b7–:b13). All orders are computed from exactly derived kernel antiderivatives and converge at the kernel’s full rate.
• In any dimension, mixed freely. Orders can be chosen per dimension and combined with interpolation or derivatives, e.g. derivative = (-2, 0) or (-3, 1, -1).
• Faster integrals. Integral construction and evaluation do substantially less work, and mixed integral/derivative evaluation is 20–30% faster.

Fixes

• Lazy mode combined with an antiderivative now raises a clear error at construction, instead of failing at the first evaluation.

Breaking changes

• The internal fields and type parameters of FastConvolutionInterpolation for integrals have changed. This only affects code that accesses these internals directly.
• Mixed integral interpolants now display as derivative=(…) instead of mixed=(…).

Migration guide

• No changes are needed for code using the public API.

Disclosure of AI assistance

This release was developed with the assistance of Claude (Anthropic).

v0.24.0 — 2026-09-24

ConvolutionInterpolations v0.24.0

Diff since v0.23.2

ConvolutionInterpolations.jl v0.24.0

Highlights

• Exact kernel evaluation. The fast path no longer uses precomputed kernel tables. Kernel weights are computed exactly from each kernel’s polynomial pieces, so values, derivatives and integrals are exact to rounding in any precision (Float32, Float64, BigFloat).
• Faster evaluation, especially in higher dimensions, e.g. 3D :b13 from ~60 μs to under 1 μs per evaluation.
• Nothing is written to disk. The kernel table artifact and the on-disk cache are gone, and so are the Artifacts, Scratch and Serialization dependencies.

Fixes

• Pure N-D integrals (derivative=-1 in all of 4 or more dimensions) skipped the trailing coefficients in every dimension and returned wrong values.
• On densely sampled grids, rounding in the grid spacing raised the error floor in proportion to the number of points. Values, derivatives, integrals and resampling now reach machine precision.

Deprecations

• The precompute and subgrid keywords no longer have any effect. They are still accepted, with a warning, and will be removed in a future release.

Breaking changes

• The type parameters of FastConvolutionInterpolation have changed. This only affects code that dispatches on or constructs that type with explicit parameters.
• Invalid subgrid/kernel combinations in convolution_resample no longer raise an error; the keyword is ignored with a deprecation warning.
• Results differ slightly from earlier versions, now being exact rather than table-interpolated.

Migration guide

• Remove precompute and subgrid from your calls.
• Optionally, delete the kernel cache written by earlier versions: ~/.julia/scratchspaces/75b5ba44-29d7-490a-960d-52c835d47405.

v0.23.2 — 2026-09-22

ConvolutionInterpolations v0.23.2

Diff since v0.23.1

ConvolutionInterpolations.jl v0.23.2

Fixes

• Fixed an error when the same depot is used by several Julia versions (e.g. via juliaup): cached kernel tables written by a newer Julia could not be read by an older one. Each Julia minor version now keeps its own cache, and unreadable cache files are regenerated automatically.
• The zero-allocation tests now pass on Julia 1.10 (LTS). The functors were already allocation-free; the tests were measuring call-site overhead on 1.10.

Other changes

• CI now tests the minimum supported Julia version, the LTS, the latest stable release and the upcoming prerelease on Linux, plus the latest stable release on Windows.

Breaking changes

None.

Migration guide

No action needed. Kernel tables are cached once more on first use after upgrading.

v0.23.1 — 2026-09-22

ConvolutionInterpolations v0.23.1

Diff since v0.23.0

ConvolutionInterpolations.jl v0.23.1

Changes

• The precomputed kernel tables used by the fast interpolation path are now shipped as a Julia artifact instead of as ~3 MB of source constants. This reduces precompilation time and source size.
• The artifact is downloaded automatically on install and verified to be bit-identical to the previous tables.
• Adds Artifacts (standard library) as a dependency.

Breaking changes

None. Results are numerically identical to v0.23.0.

Migration guide

No action needed.

v0.23.0 — 2026-09-15

ConvolutionInterpolations v0.23.0

Diff since v0.22.0

ConvolutionInterpolations.jl v0.23.0

A release focused on error handling and kernel routing. No API changes; correct existing code behaves identically.

Breaking changes

• FastConvolutionInterpolation now throws on nonuniform grids. Previously a direct call with nonuniform knots was accepted and silently interpolated on the wrong spacing (the fast kernels are precomputed on a uniform subgrid and cannot adapt their weights). Code that relied on this path was producing incorrect results; see the migration guide below.

Nothing else changes behaviour for code that previously worked correctly. Two call forms that used to error now work (:n3 on a uniform grid, and a direct ConvolutionInterpolation call with a Symbol kernel on a uniform grid).

Migration guide

If you construct FastConvolutionInterpolation directly and your knots are nonuniform in any dimension, switch to one of:

• convolution_interpolation(knots, values; ...) — the default entry point. It detects the grid and selects the fast path for uniform grids and the weight-adjusting slow path for nonuniform grids automatically. This is the recommended replacement.
• ConvolutionInterpolation(knots, values; ...) — the slow constructor, whose kernels adjust their weights to nonuniform spacing.

Calls through convolution_interpolation need no changes; it already routed nonuniform grids away from the fast path.

If you pass :a4, :a5 or :a7 on a nonuniform grid, nothing changes in behaviour — those calls errored before and still do — but the docstrings no longer claim a fallback to :n3. On nonuniform grids the supported kernels are :a0, :a1, :a3/:n3 and the full b-series (:b5, :b7, :b9, :b11, :b13); all of these work on both uniform and nonuniform grids, with the b-series giving higher order at the cost of a per-interval weight projection on nonuniform data.

:n3 and :a3 are now interchangeable

The nonuniform cubic kernel :n3 previously raised a MethodError on a fully uniform grid, because it has no uniform evaluator. It is now routed to :a3, which is mathematically identical on uniform grids and faster. Since :a3 already fell back to :n3 on nonuniform grids, the two names now describe one cubic kernel: whichever you pass, a uniform grid evaluates with the fast :a3 kernel and a nonuniform grid with the :n3 weights.

This makes :n3 (or :a3) a good default when you don’t know in advance whether the data will be uniform and want construction to stay fast either way — for example, user-supplied spectral or measurement grids. The b-series kernels also work on both grid types and give higher order, but on nonuniform data they pay a per-interval weight projection at construction; :a3/:n3 is the choice when that cost matters more than the extra order.

Bug fixes

• ConvolutionInterpolation(knots, values; kernel = :b7) — a direct call with a Symbol kernel on a uniform grid — failed with getindex(::Symbol, ::Int). The uniform builder received the raw keyword instead of the normalised kernel tuple. This was masked by convolution_interpolation, which always passes a tuple, and so never appeared in normal use.
• Docstrings corrected: only :a3 falls back to :n3 on nonuniform grids; :a4, :a5 and :a7 are uniform-only.

Tests

New tests cover :n3/:a3 equality on uniform grids, the direct-constructor Symbol path, nonuniform rejection in FastConvolutionInterpolation, and :a4/:a5/:a7 rejection on nonuniform grids. Gaussian kernel test tolerance adjusted. 1913 tests pass.

Housekeeping

• Manifest.toml removed from the repository.
• SparseArrays compat updated.

kernel-tables-v1 — 2026-09-22

v0.22.0 — 2026-08-18

ConvolutionInterpolations v0.22.0

Diff since v0.21.0

ConvolutionInterpolations.jl v0.22.0

Exact interpolation of fully scattered data — the answer to “does it handle fully scattered points?” is now yes, directly:

itp = convolution_interpolation(points, values)   # points: D x N matrix

Scattered values enter as exact linear constraints on an internal uniform kernel grid; the result is an ordinary completed grid, so the entire existing machinery, evaluation, derivatives, antiderivatives, applies unchanged.

Highlights

• Exact interpolation: passes through every data point (violations ~1e-12), ~7th-order convergence for quasi-uniform point sets, 1D-3D recommended
• Derivatives: ~6-7th-order accurate first derivatives via the standard derivative keyword
• Box integrals of scattered data: mixed antiderivatives accurate to ~1e-13 via inclusion-exclusion, to my knowledge unique among scattered-data methods
• Noisy data: mode=:lsq with a coarser gridsize gives a least-squares regression (noise suppression ~ sqrt(DOF/N)); the gridding+smoothing pipeline remains the robust choice for wild data
• Solvers: :cholesky (fast, default), :qr (lowest error floor), :nullspace (direct, small systems, also the validation oracle in the test suite)
• Per-dimension kernels, Float32/BigFloat support, selectable extrapolation, graceful least-squares fallback (with warning) for strongly clustered data

Breaking changes

None. The scattered-data entrypoint is a new method (convolution_interpolation(::AbstractMatrix, ::AbstractVector)) on a first-argument type no existing method accepts, and the new fit_scattered export introduces no conflicts. All existing calls, defaults, and results are unchanged.

Migration guide

No migration required. Existing code runs unmodified. To adopt the new feature, pass scattered points as a D x N matrix (each column one point) instead of knot vectors, see the README’s “Scattered Data Interpolation (trusted data)” section.

Other

• New dependency: SparseArrays (stdlib)
• Test suite grown to 1885 passing tests

This feature was prompted by an audience question at JuliaCon Mainz 2026 (the moment on tape) — asked on August 12, shipped in this release.

v0.21.0 — 2026-08-08

ConvolutionInterpolations v0.21.0

Diff since v0.20.1

Non-uniform path: much faster construction, rigorous convergence rates

The per-interval weight construction for non-uniform grids has been substantially optimized:

• Extended-precision arithmetic replaces exact rationals in the internal projection: BigFloat throughout, truncated to Float64 only at the final storage step. Accuracy is unchanged; construction is dramatically faster.
• Normalized local coordinates: knot spacings are scaled per interval before the projection, making the system’s conditioning independent of mesh size. Fine grids (tested to N = 10⁴ and beyond) now construct without numerical degradation.
• QR factorization for the projection solve, and parallelized per-interval computation (intervals are independent).

Together these make the non-uniform path viable at grid sizes where it was previously impractical.

The README now also states the rigorous convergence rates for non-uniform grids, verified in extended-precision arithmetic: the construction enforces polynomial reproduction up to each kernel’s full degree, so convergence order is degree + 1: 6th order for :b5, 8th for :b7, 10th for :b9, 12th for :b11, and 14th order for :b13. For :b7 and above this exceeds the uniform design order of 7. On the classical Runge benchmark, :b13 reaches Float64 machine precision at roughly 500 non-uniform samples.

Extended derivative support

Exact-arithmetic verification of the kernel piece polynomials showed the b-series kernels are smoother than previously documented: the family is C^(degree - 2), i.e. :b7 is C⁵ (previously stated C⁴), :b9 is C⁷, :b11 is C⁹, and :b13 is C¹¹. The :b5 documentation was already exact (C³).

Following this, the fast-path derivative range has been extended by one order for :b7 through :b13:

Kernel	Continuity	Derivative range (was)
:b7	C⁵	0..5 (0..4)
:b9	C⁷	0..6 (0..5)
:b11	C⁹	0..7 (0..6)
:b13	C¹¹	0..7 (0..6)

The line is drawn at the 7th derivative: although the kernels support deeper differentiation in exact arithmetic, derivative orders beyond 7 are dominated by roundoff noise amplification in Float64 (the chain-rule factor (1/h)ⁿ) and are deliberately not exposed on the fast path. The full ladder up to each kernel’s continuity class remains reachable by evaluating the exact rational kernel coefficients directly in extended precision: undocumented and unsupported, but there for those who need it.

New derivative orders are convergent on non-uniform grids (order degree + 1 − n) and were validated against analytical derivatives on randomly perturbed grids.

Breaking changes

None. The API is fully backwards compatible. One internal change may be observable: the non-uniform weight construction now uses extended-precision floating point instead of exact rational arithmetic, so stored Float64 weights may differ from v0.20.x at the last-bit level. Convergence behavior is unchanged or improved.

Migration guide

No action required. Existing code runs unchanged; the new derivative orders and faster non-uniform construction are available immediately.

Other

• Kernel-table regeneration fixes, including a corrected :b5 second-derivative table
• Subgrid selection now falls back to :linear at each kernel’s maximum derivative order (Hermite subgrids require derivative tables above the requested order)
• README: rigorous continuity classes and non-uniform convergence rates throughout; refreshed resampling benchmarks; documented the dimension-dependent default kernel (:b7 for 1–2D, :b5 for 3D, :a4 for 4–5D, :a3 for ≥6D)

Full test suite passing (1858+ tests).

v0.20.1 — 2026-08-05

ConvolutionInterpolations v0.20.1

Diff since v0.20.0

ConvolutionInterpolations.jl v0.20.1

Pure optimization patch — no API or behavior changes.

Faster wrapped evaluation

The extrapolation wrapper now caches domain bounds (lo/hi/tolerance per dimension) as plain floats at construction instead of reconstructing them from the expanded knots on every call. Wrapped evaluation is now within ~1 ns of the bare interpolant (~23 ns for 1D :b7, previously ~26 ns).

Evaluation improved modestly across all kernels and dimensions; construction of b-series kernels in 1D/2D is now indistinguishable from the low-order kernels (~2 μs).

Breaking changes

None. Domain-check tolerance semantics are identical to v0.20.0 — the tolerance is now precomputed at wrap time rather than per evaluation, with the same values.

One internals note: ConvolutionExtrapolation gained three fields (lo, hi, tol). Code constructing this struct positionally or accessing its fields directly (not part of the public API) must account for them.

Migration guide

Nothing to do — update and benefit.

Full test suite passes.

v0.20.0 — 2026-08-05

ConvolutionInterpolations v0.20.0

Diff since v0.19.1

ConvolutionInterpolations.jl v0.20.0

This release changes the default kernel selection to a dimension-dependent policy and fixes a long-standing boundary issue where evaluating at exact domain endpoints could throw under the default extrap=:throw.

Dimension-dependent default kernel

The default kernel keyword is now :auto, which selects the kernel based on dimension:

Dimension	Default kernel
N ≤ 2	:b7
N = 3	:b5
N = 4–5	:a4
N ≥ 6	:a3

The policy trades kernel width against tensor-product evaluation cost: in 1D and 2D the wider :b7 kernel improves accuracy at near-negligible evaluation cost, while in higher dimensions the (2·eqs)^N support growth favors narrower kernels. The keyword constructors and the positional fast constructors now share a single policy and cannot drift apart.

Passing an explicit kernel symbol or tuple behaves exactly as before; :auto may also be passed explicitly to opt into the policy.

Boundary tolerance fix

Evaluating at exact domain endpoints — for example itp(2π) on knots range(0, 2π, length=n) — could previously throw a domain error under the default extrap=:throw, because the domain bounds are reconstructed from the expanded knot arithmetic and could differ from the original endpoints by one ulp.

Domain checks now allow ulp-scale slack (relative to the bound magnitudes, with an absolute floor near zero) and clamp near-boundary queries before evaluation. Exact endpoints and ulp-adjacent points evaluate correctly; genuinely out-of-domain points still throw. The tolerance forgives floating-point rounding, not user error.

Breaking changes

• Default results change for constructions that did not specify a kernel. A 1D or 2D construction with default arguments now uses :b7 instead of :b5, producing slightly different (generally more accurate) values. Constructions with an explicit kernel= argument are unaffected.
• Near-boundary behavior under extrap=:throw changes. Queries within one ulp of the domain boundary previously threw and now evaluate. Code relying on the throw as a strict boundary detector at ulp precision will observe different behavior.

No exported function signatures were removed or renamed; the kernel keyword’s default value changed from :b5 to :auto.

Migration guide

• To reproduce v0.19.x default results exactly: pass kernel=:b5 explicitly.
• If you stored reference outputs generated with default kernels in 1D or 2D, regenerate them or pin kernel=:b5.
• If you relied on extrap=:throw for strict domain enforcement: the enforcement still holds beyond ulp-scale tolerance; no action needed unless you specifically depended on exact-endpoint queries throwing.
• Everything else — explicit kernels, extrapolation modes, boundary conditions, derivatives, lazy mode — is unchanged.

Full test suite passes.

v0.19.1 — 2026-08-04

ConvolutionInterpolations v0.19.1

Diff since v0.19.0

ConvolutionInterpolations.jl v0.19.1

Performance-focused patch release: constructors are ~2× faster with ~20× less memory churn in low dimensions, evaluation is slightly faster across the board, and one latent 2D bug is fixed. No API changes: existing code benefits automatically.

Faster construction

Precomputed kernel tables are now cached and shared between interpolants of the same kernel, derivative order, and element type, instead of being copied on every construction. Together with allocation reductions in the knots and boundary-setup paths, a 1D :b5 construction at 25 points goes from 3.1 μs / 33.8 KiB to 0.88 μs / 1.7 KiB. Construction is now faster than a comparable cubic spline setup in our benchmarks while retaining higher convergence order.

Range-typed knots (e.g. range(0, 1, 100)) are now preserved end-to-end rather than collected into vectors, eliminating copies during setup.

Faster evaluation

Evaluation functors now cache the grid origin and compute cell offsets arithmetically instead of indexing the knot vector, removing an index and a division per dimension per call. 1D :b5 evaluation drops from ~23 ns to ~21 ns, with similar small gains in higher dimensions. This also removes an evaluation penalty that range-typed knots would otherwise incur.

Bugfix

• Fixed the 2D quintic subgrid mode using the x-direction table spacing in the final y-direction Hermite step. Results were unaffected unless precompute differed between the two dimensions.
• Fixed a missing element-type conversion for :a1 antiderivative kernel tables at non-Float64 precision.

Internal

• Removed a legacy unused workspace array from the boundary-condition machinery.

Breaking changes

None in the API sense — no exported function signatures, keyword arguments, or defaults changed, and no user code needs modification.

Two of the bugfixes are strictly output-changing, which automated tooling may flag as breaking:

• Numerical results for 2D quintic subgrid with per-dimension precompute values change (they were previously wrong due to the Hermite spacing bug). Results with the default or uniform precompute are unaffected.
• Kernel table element types for :a1 antiderivative at non-Float64 precision are now correct (previously leaked Float64), which can shift results at the precision of the requested type.

Both are corrections toward the documented behavior, released as a patch per common Julia ecosystem practice for bugfixes.

Two internal representation changes are visible only to code reaching into non-API internals:

• FastConvolutionInterpolation gained an internal x0 field; code constructing this struct positionally must add it. The public constructors are unaffected.
• Knots passed as ranges are now stored as ranges rather than collected Vectors, so the knots field type parameter may differ from previous versions.

Migration guide

For users of the public API (convolution_interpolation, FastConvolutionInterpolation(knots, vs; kwargs...), ConvolutionInterpolation(knots, vs; kwargs...)): nothing to do — update and enjoy the speedup.

If you compared stored numerical outputs across versions and used 2D quintic mode with asymmetric precompute, regenerate your reference values; the new results are the correct ones.

If you dispatched on or accessed the internal struct fields directly (not part of the public API), account for the new x0::NTuple{N,T} field and the possibility of range-typed knots.

Full test suite (1858 tests) passes.

AI Assistance

I was assisted by Anthropic’s Claude in debugging the code and improving performance.

v0.19.0 — 2026-08-02

ConvolutionInterpolations v0.19.0

Diff since v0.18.0

ConvolutionInterpolations v0.19.0

Breaking changes

• bc=:detect (the default) now selects the polynomial ghost extension by a different criterion. Interpolated values near domain boundaries will differ from v0.18. Interior values are unchanged. The new guard tests the ns-th finite-difference residual of the extension against the local data range, replacing a first/second/third-difference smoothness heuristic that rejected the polynomial extension on ~43% of smooth inputs, including data the extension reproduces exactly. On 20,000 held-out random signals the rate of selecting the worse branch drops from ~8% to ~3%.
• convolution_gaussian now throws for B > 0.34. Above that value the kernel’s between-node relative error exceeds 1e-12. convolution_smooth evaluates on grid nodes only and is unaffected at any B.
• convolution_smooth and convolution_gaussian return different values near boundaries, because two bugs there were fixed (see below). Constant data now reproduces to machine precision at the boundary; previously there was a 16% deficit at B=0.1.
• bc=:poly and bc=:linear change only at roundoff level, from local rather than global mean-centring in the boundary setup.

Migration guide

• No source changes are required. The API is unchanged; only numerical results near boundaries differ.
• If you use convolution_gaussian with B > 0.34, either lower B (0.1 is the usual choice) or switch to convolution_smooth, which is exact at any B when you only need values on the grid nodes.
• If you have tests pinned to exact boundary values from v0.18, regenerate the expected values, or relax them to a tolerance.
• To keep the old behaviour, pin the previous minor version in Project.toml:

  [compat]
  ConvolutionInterpolations = "0.18"

Fixed

• convolution_smooth and convolution_gaussian allocated eqs_d ghost points but passed eqs = eqs_d to the boundary fill, which uses the convention ghosts = eqs - 1. The ghost slot adjacent to the data was never written and stayed zero, silently dropping one kernel tap at each boundary.
• The Gaussian ghost fallback subtracted the mean twice, offsetting every ghost point by the mean of the input.

Performance

• Boundary setup now reads only the ns+3 samples it needs rather than the whole line. Construction is ~1.4× faster in 2D and ~1.7× in 3D, for every kernel including the a-family. Evaluation cost is unchanged.

Tests

• 1858 passing (up from 1826). New: constant preservation, linear reproduction, offset invariance and boundary symmetry for the Gaussian path; guard exactness and allocation-free checks for :detect.

Parts of this release were developed with assistance from Claude (Anthropic); all code and results were verified by the author.

v0.18.0 — 2026-07-13

ConvolutionInterpolations v0.18.0

Diff since v0.17.1

ConvolutionInterpolations.jl v0.18.0

Breaking changes

None. No exported API was changed or removed: v0.18.0 is additive. All v0.17.x code runs unchanged.

Migration guide

No action required. Existing code is unaffected. Optionally:

• Scattered data pipelines using an external gridding package can switch to the new built-in scattered_to_grid.
• Smoothing calls tuned to small B values may benefit from the updated guidance (B ≈ 0.1); see the README.

New features

scattered_to_grid – scattered data gridding

Grid unstructured data onto a regular grid by nearest-neighbor assignment, in any number of dimensions:

grid = scattered_to_grid(points, values, (xs, ys); k=1) # nearest point (k=1, default)
smooth = convolution_smooth((xs, ys), grid, 0.1) # B=0.1 is recommended
itp = convolution_interpolation((xs, ys), smooth)

This completes the scattered-data-to-interpolant pipeline within the package: gridding, smoothing, and high-order interpolation/differentiation/integration of noisy scattered data, end to end. Nearest-neighbor gridding is deliberately simple: its error is spatially uncorrelated, which is exactly what the Gaussian smoother removes best.

Full Float32 support

Float32 in, Float32 out:

• Kernel tables are converted to the working precision at construction
• Scalar parameters (e.g. the smoothing parameter B) accept any Real
• Query points promote to the interpolant’s type at evaluation

Together with the existing BigFloat support, the package now preserves types and computes at native precision from Float32 through BigFloat.

Changed

• Smoothing guidance: B ≈ 0.1 is recommended as a default (previously B ≤ 0.1). Testing across 1D–4D shows the optimum is shallow and two-sided: smaller B over-smooths and attenuates extrema, larger B leaves noise behind.
• README: rewritten scattered data pipeline example (no external packages required), updated figures.

Dependencies

• New: NearestNeighbors.jl (KD-tree queries for scattered_to_grid)

Tests

• 1699 → 1826 passing tests, including new suites for scattered_to_grid (brute-force equivalence in 2D/3D, grid orientation, duplicate points, type handling) and Float32 (eltype preservation across all major paths, accuracy verification, mixed-type regression guards)

v0.17.1 — 2026-06-01

ConvolutionInterpolations v0.17.1

Diff since v0.17.0

v0.17.1

Bugfixes, expanded testsuite, and :natural boundary condition reintroduced

List of improvements - Fixing a few stale things in the readme (tuple vs vector in boundary conditions) - Fix several printing bugs which had passed the full test suite - These bugs only surfaced when tested directly in the REPL - Another small bugfix relating to the nonuniform :n3 kernel (it was missing in a dictionary) - Reintroduced the :natural extrapolation option, usefull in lower dimensions (also in the readme) - Include testing in the testsuite of the :natural boundary condition - 1699 tests now passing

v0.17.0 — 2026-04-13

ConvolutionInterpolations v0.17.0

Diff since v0.16.0

v0.17.0

New features

• boundary_fallback support for nonuniform lazy :n3 (1D/2D/3D/ND). Near-boundary evaluation falls back to linear interpolation, consistent with uniform lazy behavior. Avoids expensive ghost point resolution in high-dimensional corners (e.g. 3⁸ = 6,561 weighted lookups reduced to 2⁸ = 256 linear interpolation).

Bug fixes

• Error guard for slow lazy functors. Attempting lazy=true, fast=false on uniform grids now raises a clear error at construction time, instead of silently compiling and failing with MethodError at evaluation.

Tests

• Split lazy tests into separate uniform and nonuniform files.
• Added boundary_fallback tests for nonuniform lazy :n3 in 1D/2D/3D.
• 1241 tests passing.

Breaking changes

This is a non-breaking release. All existing APIs, defaults, and behaviors are unchanged.

Migration guide

No migration steps required. Existing code will work without modification. New features (boundary_fallback for nonuniform lazy :n3) are opt-in and default to previous behavior.

v0.16.0 — 2026-04-12

ConvolutionInterpolations v0.16.0

Diff since v0.15.0

v0.16.0

Breaking changes

• Removed :periodic boundary condition. This was legacy code from early development that predated the polynomial boundary condition system. It did not converge and was not suitable for production use.
• Removed dual extrapolation option. The double-pass boundary extrapolation mode has been removed.
• Removed slow lazy functors and lazy_ghost_value. These untested legacy code paths have been removed. Fast lazy evaluation (lazy=true) is unaffected.
• convolution_resample default kernel changed from :b9 to :b13. The cost difference is negligible while accuracy improves.
• convolution_resample default subgrid changed to :cubic. This eliminates precomputation for resampling and improves accuracy.

Migration guide

• Replace bc=:periodic with bc=:poly (recommended) or bc=:detect (default). Both offer superior accuracy with proven convergence.
• If you relied on the dual extrapolation option, use standard single-pass extrapolation instead (extrap=Line(), extrap=Flat(), or extrap=Throw()).
• If you need the previous resampling defaults, pass kernel=:b9, subgrid=:linear explicitly to convolution_resample.
• If 3D evaluation speed is critical, pass subgrid=:linear to recover the previous linear subgrid performance. The default :cubic subgrid is slower per evaluation but significantly more accurate.

New features

• 3D cubic Hermite subgrid. 3D evaluation now uses cubic Hermite corner accumulation (64 accumulators, 21 Hermite calls) for both eager and lazy functors. This provides substantially better accuracy and convergence compared to linear subgrid.
• Zero-caching default path. With :cubic subgrid as default, precalculation is only triggered for BigFloat or top-derivative evaluation. New users will not encounter precomputation on first use in typical workflows.
• Subgrid hierarchy. Cubic Hermite is now the default for 1D–3D. 4D and above default to linear subgrid where the cost of higher-order evaluation grows too quickly.
• Shared subgrid helpers. New eval_sg_k and eval_sg_kt functions centralize subgrid kernel evaluation, reducing code duplication across functor paths.
• 5D mixed integral tests. New test coverage for the N-integral-dims functor, bringing the total to 1184 tests.

Parts of this release were developed with assistance from Claude (Anthropic). All code, methods, and scientific content have been verified and validated by the author.

v0.15.0 — 2026-04-06

ConvolutionInterpolations v0.15.0

Diff since v0.14.0

v0.15.0 Release Notes

⚠️ Breaking Changes

1. Boundary condition :auto renamed to :detect

The default boundary condition symbol has been renamed from :auto to :detect to better reflect its behavior, it detects the best boundary condition automatically based on the data.

Migration: Replace any explicit bc=:auto with bc=:detect. Users relying on the default behavior (never setting bc explicitly) are unaffected.

# Before
itp = convolution_interpolation(x, y; bc=:auto)
# After
itp = convolution_interpolation(x, y; bc=:detect)

2. Gaussian smoothing split into dedicated functions

The B parameter previously accepted by convolution_interpolation has been removed. Gaussian smoothing is now available through two dedicated functions that provide a cleaner interface and significantly better performance.

Migration:

# Before
itp = convolution_interpolation(x, y; B=0.05)

# After, for an interpolant at arbitrary locations:
itp = convolution_gaussian(x, y, 0.05)

# After, for a smoothed array on the same grid:
y_smooth = convolution_smooth(x, y, 0.05)

---

🆕 New Features

convolution_resample, High-order grid-to-grid resampling

Resample data from one uniform grid to another using high-order separable convolution interpolation. Significantly faster than constructing a full interpolant and evaluating at each output point, with the advantage growing with dimension:

Grid	convolution_resample	construct + eval	Speedup
1D 500→1000	22 μs	29 μs	1.3×
2D 100²→200²	1.9 ms	7.0 ms	3.6×
3D 30³→60³	8.9 ms	127 ms	14×

• Default kernel :b9, higher accuracy than convolution_interpolation’s default :b5
• O(N) allocations independent of grid size
• 1D to 10D type-stable default wrappers
• Supports per-dimension kernels and derivative resampling
• Antiderivatives not supported, use convolution_interpolation with derivative=-1

See Grid Resampling in the README.

convolution_smooth, Separable Gaussian smoothing

Apply Gaussian smoothing to recover clean signals from noisy data. Returns an array of the same size as the input.

Grid	B=0.02	B=0.05	B=0.1
1D n=1000	22 μs	13 μs	9.6 μs
2D 200²	1.7 ms	0.99 ms	0.72 ms
3D 50³	9.3 ms	5.7 ms	4.4 ms
4D 40⁴	273 ms	206 ms	169 ms

• 1000× speedup over previous Gaussian path via exp weight hoisting, weights precomputed once per call
• Adaptive eqs÷3 boundary fit, wider kernels use more boundary points automatically
• O(N) allocations independent of grid size
• B ≤ 0.1 recommended for effective noise reduction

See Smoothing Noisy Data in the README.

convolution_gaussian — Standalone Gaussian interpolant

Construct a Gaussian smoothing interpolant that supports evaluation at arbitrary locations, just like convolution_interpolation.

itp = convolution_gaussian(x, y_noisy, 0.05)
itp(1.5)  # evaluate anywhere in the domain

• Allocation-free evaluation via Val{B} and Val{NT} type encoding of kernel parameters
• Dedicated 1D, 2D, 3D and ND evaluation functors

Scattered Data Pipeline

The two new functions combine naturally with RBF packages for a complete scattered data workflow — see Scattered Data Pipeline in the README.

---

🚀 Performance Improvements

Full type stability throughout

• Ref{Array{T,N}} ping-pong buffers eliminate Core.Box boxing of loop variables
• Val{UG} branch elimination removes allocating empty matrices in boundary condition paths
• BoundaryWorkspace{T,1} via Val(1) ensures fully concrete workspace types
• Type-stable 1D to 10D default convenience wrappers, verified with @code_warntype

Lazy evaluation functors reworked

The 1D, 2D, and 3D lazy evaluation functors have been fully reworked for type stability and performance. Furthermore: lazy construction speedup for 4D: 1320ms → 79ms.

ND pure antiderivative allocation-free

The N-dimensional pure antiderivative path is now allocation-free and supports linear/cubic/quintic subgrid selection, achieved via nexprs/nloops macros.

---

🧪 Tests

• 1178 tests passing (was 1086, +92 new tests)
• Comprehensive test suite for convolution_resample including derivative scaling
• Comprehensive test suite for convolution_smooth including boundary behavior
• Lazy mode test suite substantially expanded, 1D/2D/3D with boundary_fallback=false, 4D+ with boundary_fallback=true, nonuniform :n3, construction speed

---

📖 Documentation

README updated with: - New sections: Smoothing Noisy Data, Grid Resampling, Scattered Data Pipeline - All performance tables based on updated benchmarks - B-value comparison figure showing effect of smoothing parameter - Updated lazy construction timings reflecting actual improvements - Updated Performance Guidelines with recommendations for new functions

---

Parts of this release were developed with assistance from Claude (Anthropic). All code, methods, and scientific content have been verified and validated by the author.

v0.14.0 — 2026-03-28

ConvolutionInterpolations v0.14.0

Diff since v0.13.0

v0.14.0: Fast Mixed Integral/Derivative, Zero-Allocation Evaluation

Highlights

The headline feature of v0.14.0 is a new class of fast mixed integral/derivative operators that combine antiderivatives, derivatives, and interpolation across different dimensions in a single allocation-free evaluation call. This extends the fast antiderivative path introduced in v0.13.0, which accelerated pure antiderivatives using precomputed tail sums, to the fully general mixed case.

New Features

Fast mixed integral/derivative operators

In v0.13.0, derivative=(-1,-1) in 2D used tail sums to achieve O(stencil^2) evaluation independent of grid size for pure antiderivatives. v0.14.0 extends this to arbitrary mixtures of integrals, derivatives, and interpolation across dimensions. For example:

# integrate in x, differentiate in y, interpolate in z, single O(stencil^3) call
itp = convolution_interpolation((xs,ys,zs), vs; derivative=(-1, 1, 0))

# definite integral via differences cancels the unknown constant:
itp(x2, y, z) - itp(x1, y, z)  # ≈ ∫ₓ₁ˣ² f dx · ∂f/∂y · f(z)

The fast path uses precomputed tail sums to evaluate contributions from coefficients outside the local stencil in O(1) per integral dimension, reducing the total cost to O(stencil^N) where N is the number of dimensions, the same asymptotic cost as pure interpolation. This is achieved by decomposing the evaluation domain, each part handled by a dedicated lookup into precomputed tail sum arrays built at construction time.

Mixed operators are supported for any combination of integral (-1), interpolation (0), and derivative (≥1) orders per dimension, in any number of dimensions. Up to 3 integral dimensions use the specialized fast path; 4 or more integral dimensions fall back to a general accurate path.

The direct (non-fast) path for mixed operators is also available via fast=false, and agreement between fast and direct paths is verified in the test suite.

End-to-end zero-allocation evaluation

All non-lazy evaluation paths, including the extrapolation wrapper, are now verified allocation-free via @allocated tests. This covers all kernels (:a0 through :b13), all derivative/integral/mixed combinations, and 1D through 4D. The key enabler was replacing runtime fields and symbol-based dispatch with compile-time singleton structs and Val-based type parameters throughout.

BigFloat precision preservation

All kernels correctly propagate BigFloat arithmetic through construction and evaluation without silent downcast to Float64. B-series kernels on dense grids (n=2000) achieve errors below 1e-16, well below Float64 precision, verified by new precision tests. This means the package can be used directly for extended precision computation by simply passing BigFloat grids and values.

Expanded test suite

Test count grows from 696 to 941, adding: - Mixed integral/derivative convergence in 2D, 3D, 4D (fast and direct paths) - Fast vs direct agreement tests in 2D, 3D, 4D - Zero allocation tests across all functors and derivative combinations - BigFloat type preservation and machine precision tests

Breaking Changes

Low-level extrapolation API: symbols → structs

The low-level ConvolutionExtrapolation constructor now requires singleton structs instead of symbols:

# v0.13 (no longer supported at low level)
ConvolutionExtrapolation(itp, :throw)

# v0.14
ConvolutionExtrapolation(itp, Throw())
ConvolutionExtrapolation(itp, Line())
ConvolutionExtrapolation(itp, Flat())

This enables compile-time dispatch and is required for zero-allocation evaluation.

Migration Guide

For the convenience constructor (convolution_interpolation): no changes needed. Both symbols and structs are accepted:

# both still work
itp = convolution_interpolation(x, y; extrap=:line)
itp = convolution_interpolation(x, y; extrap=Line())

For the low-level ConvolutionExtrapolation: replace symbols with the corresponding struct:

Symbol	Struct
:throw	Throw()
:line	Line()
:flat	Flat()

Declaration of AI Assistance

The development of v0.14.0 was carried out with assistance from Claude (Anthropic), including design of the fast mixed integral/derivative functors, type system cleanup, dispatch debugging, test suite development, and documentation. All code, methods, and scientific content have been verified and validated by the author.

v0.13.0 — 2026-03-22

ConvolutionInterpolations v0.13.0

Diff since v0.12.0

ConvolutionInterpolations.jl v0.13.0 - 3D Fast Antiderivative Support

New Features

3D fast antiderivative evaluation The 3D antiderivative functor now uses a 27-region decomposition (1 center + 6 faces + 12 edges + 8 corners) with precomputed prefix/suffix sum arrays, achieving ~1000× speedup over the slow path on a 100×100×100 grid. Accuracy reaches 1e-15 with cubic subgrid interpolation.

Cubic and quintic subgrid modes are now enabled for 1D, 2D, and 3D integral evaluation (previously forced to :linear), giving significantly better accuracy at minimal extra cost.

Breaking Changes

:nearest subgrid removed The :nearest subgrid option has been removed. It was undocumented, only partially implemented for 1D, and offered poor accuracy relative to :linear at negligible speed gain.

Migration Guide

:nearest subgrid If you were using subgrid=:nearest, switch to subgrid=:linear:

# Before
itp = convolution_interpolation(x, vs; subgrid=:nearest)
# After
itp = convolution_interpolation(x, vs; subgrid=:linear)

No other API changes are required. All other keyword arguments, kernel names, derivative orders, and boundary condition symbols are unchanged.

Internal Changes

These changes do not affect the public API. - Constructors simplified and unified to normalize all inputs to per-dim tuples - Functors reorganized for clarity and dispatch correctness - New abstract types added to prevent wrong functor dispatch - Improved type stability of 1D and 2D integral tail sums - Test suite reorganized, trimming redundant tests and adding new coverage for nonuniform per-dim kernel selection, per-dim derivative combinations, boundary conditions, extrapolation, subgrid downgrade, and Gaussian smoothing - All 696 tests pass and run in 4 minutes

v0.12.0 — 2026-03-18

ConvolutionInterpolations v0.12.0

Diff since v0.11.0

ConvolutionInterpolations.jl v0.12.0

New Features

:a0 antiderivative support

The nearest-neighbor kernel now supports derivative=-1. The antiderivative of a piecewise constant function is piecewise linear, and integrates exactly. A delta-like spike integrates to a perfect Heaviside step:

n = 1001
x = range(0.0, 2π, length=n)
h = x[2] - x[1]
y = [abs(xi - π) ≈ 0.0 ? 1.0/h : 0.0 for xi in x]

itp = convolution_interpolation(x, y; kernel=:a0, derivative=-1)
itp(π + 0.1)  # = 1.0
itp(π - 0.1)  # = 0.0

This also works in multiple dimensions and with per-dimension kernel tuples, enabling for example a sharp delta selector in one dimension combined with smooth high-order differentiation in another:

itp = convolution_interpolation((x, y), vals; kernel=(:a0, :b7), derivative=(-1, 1))

Breaking Changes

• bc=:polynomial renamed to bc=:poly. This is a one-word change that makes the API more concise and consistent with the symbol-based convention used throughout.

Migration Guide

# Before
itp = convolution_interpolation(x, y; bc=:polynomial)

# After
itp = convolution_interpolation(x, y; bc=:poly)

Improvements

• Multidimensional integration convergence claims updated: order is empirically flat (~7th) for general smooth integrands, but increases dramatically for symmetric integrands (from ~7th in 1D to ~12th in 5D for :b7). The mechanism is not yet fully understood.
• Technical background section corrected: left_values is precomputed for all dimensions; prefix sums for O(1) tail lookups are available in 1D and 2D only.
• Delta spike → Heaviside example added to README.
• Kernel reference table updated: :a0 derivative range now -1..0.

Bug Fixes

• Fixed per-dimension constructor crash when :a0 or :a1 appeared in a kernel tuple alongside higher-order kernels (e.g. kernel=(:a0, :b7)). Ghost point computation now correctly skips dimensions with eqs=1.

Tests

• New :a0 antiderivative testset: constant function, delta spike, anchor condition, fast vs direct agreement.
• Full test suite passing.

Co-developed with assistance from Claude (Anthropic).

v0.11.0 — 2026-03-17

ConvolutionInterpolations v0.11.0

Diff since v0.10.1

ConvolutionInterpolations.jl v0.11.0 — API Stabilization

This release completes the API cleanup in preparation for v1.0.0. All breaking changes are simple renames, motivated by consistency and a gradual departure from conventions inherited from Interpolations.jl that didn’t fit this package’s design.

Breaking changes

Before	After
degree=:b5	kernel=:b5
kernel_bc=:auto	bc=:auto
extrapolation_bc=Throw()	extrap=:throw
Throw(), Line(), Flat(), Natural()	:throw, :line, :flat, :natural
itp.deg	itp.kernel_sym
itp.kernel_bc	itp.bc

Periodic and Reflect extrapolation modes removed — use bc=:periodic at construction for smooth periodicity instead, ideally combined with extrap=:natural.

lazy and boundary_fallback are now fully independent — high-dimensional lazy mode no longer silently forces boundary_fallback=true.

Julia compat bumped to 1.8.

Migration guide

# Before
itp = convolution_interpolation(x, y;
    degree=:b7,
    kernel_bc=:polynomial,
    extrapolation_bc=Line(),
    derivative=1)

# After
itp = convolution_interpolation(x, y;
    kernel=:b7,
    bc=:polynomial,
    extrap=:line,
    derivative=1)

Improvements

• Test suite tightened throughout — most tolerances now 1e-6, no boundary avoidance anywhere
• README: scattered data workflow, 2D integration convergence plot, :natural extrapolation description, updated lazy/boundary_fallback documentation
• show improved for nonuniform kernel display and multi-dimension BC condensation

Co-developed with assistance from Claude (Anthropic).

v0.10.1 — 2026-03-16

ConvolutionInterpolations v0.10.1

Diff since v0.10.0

v0.10.1 — Major Bugfix: All 2D/3D Operations Now at Machine Precision

Bug Fixed

A fundamental orientation error in the boundary condition ghost point computation has been corrected. When extracting 1D slices for right-boundary polynomial extrapolation in N>1 dimensions, the slice was ordered boundary->interior rather than interior->boundary. This caused the polynomial coefficients to extrapolate in the wrong direction, corrupting all ghost points on right boundaries in 2D and higher dimensions.

The fix is a single-line correction in the slice extraction loop:

# Before
workspace.slice[j] = c[idx - workspace.slice_offset[j]]
# After  
workspace.slice[j] = c[idx - workspace.slice_offset[n_dim - j + 1]]

Impact

All multidimensional operations are affected: interpolation, differentiation up to 6th order, antiderivatives, and any mixed combination thereof, in 2D and higher. All now converge to machine precision as expected from the tensor product theory.

• 1D: Unaffected: already at machine precision
• 2D/3D interpolation, differentiation, antiderivatives, and mixed operations: Now all at machine precision

Convergence Results

The fix unlocks the full power of the tensor product approach. For integration specifically, an unexpected dimensional enhancement is observed:

Dimensions	Kernel order (1D)	Observed convergence order
1D	7th	7th
2D	7th	~9th
3D	7th	~11th
4D	7th	~13th

A 10×10×10 grid already achieves ~1e-4 accuracy for smooth 3D integrands, with machine precision around n=100 per dimension.

Note

Despite being a patch release, this fix unlocks the true strength of ConvolutionInterpolations.jl: high-order machine precision evaluation of any combination of interpolation, differentiation, and integration in any number of dimensions, from a single set of sample points.

v0.10.0 — 2026-03-15

ConvolutionInterpolations v0.10.0

Diff since v0.9.0

ConvolutionInterpolations.jl v0.10.0: Fast 1D/2D Antiderivative, Boundary Condition Fix

Breaking Change: :detect Kernel Boundary Condition Removed

The :detect kernel boundary condition has been removed. This option attempted to automatically detect periodic signals and apply a mixed linear-prediction boundary condition per slice. In practice, this per-slice detection was unreliable, different slices of the same 2D or 3D array could receive different boundary conditions, creating asymmetric coefficient arrays and degrading interpolation accuracy, particularly for wide kernels (:b9, :b11, :b13) with few grid points.

Migration guide: Replace kernel_bc=:detect with one of: - kernel_bc=:polynomial, optimal for most cases, preserves polynomial reproduction (default when sufficient points) - kernel_bc=:auto, uses :polynomial when sufficient points, falls back to :linear automatically - kernel_bc=:linear, safe, stable fallback for any function type - kernel_bc=:quadratic, second-order extrapolation (use with caution on wide kernels) - kernel_bc=:periodic, for known periodic signals with insufficient points for :polynomial

The default fallback for insufficient grid points has also changed from :detect to :linear. The :linear kernel boundary condition is unconditionally stable, it never produces the unbounded ghost point growth that the :quadratic kernel boundary condition can exhibit on small grids with wide kernels.

---

Fast Antiderivative: O(eqs) in 1D, O(eqs^2) in 2D

The fast path for derivative=-1 has been fundamentally redesigned. Previously, despite being called “fast”, the antiderivative evaluation looped over the entire coefficient array, O(N) per eval in 1D, O(N^2) in 2D, giving no advantage over the direct path for large grids.

What changed: The kernel antiderivative K̃(s) saturates to exactly +/- ½ outside the stencil support. This means the contribution of all coefficients outside the local stencil can be precomputed at construction time as a cumulative sum, the “tail”. At evaluation, the total integral decomposes into:

• Local stencil (O(eqs) terms): per-j kernel lookup, unchanged from before
• Left/right tails (O(1)): precomputed prefix sums, single array lookup

In 2D, the domain decomposes into 9 quadrants around the stencil. The four corner quadrants (both dimensions in tail) are O(1) via 2D cross-sum arrays. The four edge quadrants (one dimension in stencil, one in tail) are O(eqs) via 1D prefix lookups. Only the center quadrant (both in stencil) requires the O(eqs^2) double loop.

Speedup at N=200:

Kernel	Fast vs direct
:a3	115x
:b5	44x
:b7	1,513x
:b11	1,957x
:b13	2,117x

At N=1000, speedups reach 11,000x for :b13. Evaluation time is now flat (O(1) in N) for all grid sizes.

---

Boundary Condition Improvements

Polynomial path (right boundary sign fix): The polynomial kernel boundary condition (:polynomial) matrix is applied to the reversed signal at right boundaries. A curvature continuity check now correctly determines whether the ghost points need to be negated, handling even functions (cosine-like) and odd functions (sine-like) correctly in all dimensions. Previously, a hardcoded dimension check (dim == 1) was used, which failed for even functions in the second dimension.

Recursive path: (:linear,:quadratic,:periodic) The signal is now reversed before recursive extrapolation at right boundaries, consistent with the polynomial path. This fixes asymmetric coefficient arrays that were visible when plotting itp.coefs for 2D even functions.

---

Tests

All 1046 tests pass (up from 1038). Several tolerances tightened.

---

Co-developed with assistance from Claude (Anthropic).

v0.9.0 — 2026-03-12

ConvolutionInterpolations v0.9.0

Diff since v0.8.0

ConvolutionInterpolations.jl v0.9.0

New Features

Per-dimension kernels and derivatives

Both the uniform fast path and nonuniform b-kernel path now support per-dimension kernel degrees and derivative orders via tuples:

# Different kernel per dimension
itp = convolution_interpolation((x, y), vals; degree=(:b5, :b9))

# Different derivative order per dimension
itp = convolution_interpolation((x, y), vals; degree=(:b7, :b7), derivative=(1, 0))

# Works on nonuniform grids too
itp = convolution_interpolation((x_nu, y_nu), vals; degree=(:b5, :b9), derivative=(0, 1))

On nonuniform grids each dimension is padded by its own stencil radius, with no wasted memory. Mixed uniform/nonuniform grids are handled per-dimension automatically.

Mixed antiderivative/derivative orders (uniform grids)

Pass a tuple with -1 entries to integrate in some dimensions and differentiate in others. This enables Leibniz-rule-style operations in a single interpolation object:

# Integrate in x, differentiate in y
itp = convolution_interpolation((x, y), vals; degree=:b7, derivative=(-1, 1))
itp(1.0, 2.0)  # ≈ ∂/∂y ∫ f dx

Both slow and fast evaluation paths are supported. Not available on nonuniform grids.

Pretty printing

ConvolutionInterpolation, FastConvolutionInterpolation, and ConvolutionExtrapolation now display concise, human-readable summaries instead of the raw type parameter wall:

FastConvolutionInterpolation(Float64, 2D, coefs=38×42, degree=(:b5, :b9), derivative=(1, 0))
ConvolutionExtrapolation(ConvolutionInterpolation(Float64, 2D, coefs=40×44, degree=(:b5, :b9), derivative=(0, 0), nonuniform), bc=Throw)

Tests

• 1038 tests passing (up from 754)
• Added test_nonuniform_perdim_derivatives.jl, test_perdim_derivatives.jl, and test_mixed_integral.jl which test combinations for per-dim uniform/nonuniform integrals/derivatives and scalar/tuple consistency

Breaking Changes

None. All existing code continues to work without modification. Scalar degree and derivative arguments behave identically to before.

Migration Guide

No migration needed. This release is fully backwards compatible.

The new tuple syntax is purely additive, existing scalar usage is unchanged:

# These still work exactly as before
itp = convolution_interpolation(x, y; degree=:b7)
itp = convolution_interpolation((x, y), vals; degree=:b7, derivative=1)

v0.8.0 — 2026-03-08

ConvolutionInterpolations v0.8.0

Diff since v0.7.0

ConvolutionInterpolations.jl v0.8.0 Release Notes

Breaking changes

None. All existing API defaults are preserved. All existing code works without modification.

New features

Antiderivative support (derivative=-1)

Pass derivative=-1 to compute the indefinite integral of the interpolant:

x = range(0.0, 2π, length=100)
y = sin.(x)

itp = convolution_interpolation(x, y; derivative=-1)
itp(1.0)   # ≈ 0.4597  (= 1 - cos(1), anchored at x=0)

The result is zero-anchored at the leftmost interior knot (F(anchor) = 0 exactly). In N dimensions the result is the mixed partial antiderivative, consistent with the tensor product structure used throughout the package.

Supported kernels (uniform grids only): :a1, :a3, :a4, :a5, :a7, :b5, :b7, :b9, :b11, :b13. Not yet supported on non-uniform grids.

Implementation. Each kernel K is analytically integrated to produce K̃ with coefficients stored as exact rational numbers. The fast path (default) precomputes the anchor-side sum Σⱼ cⱼ · K̃((anchor − xⱼ)/h) once at construction and stores it as left_values — one Vector{Float64} per dimension. Evaluation then only computes the K̃(x) half, giving zero-allocation O(stencil) performance identical in character to derivative=0. For b-series kernels the anchor-side sum uses exact Rational{BigInt} arithmetic before conversion to Float64, preserving the kernel’s full precision.

Dispatch type cleanup

IntegralOrder and DerivativeOrder{DO} are now the canonical types for the derivative slot of the interpolant type. This replaces the previous bare Val{0} usage and resolves a dispatch ambiguity in nonuniform functors that could produce incorrect method resolution at runtime.

Tests

754 tests (up from 671), including 83 new antiderivative tests covering 1D/2D/3D convergence for fast and direct paths, fast vs. direct agreement, anchor condition (F(anchor) = 0), and all supported kernels.

Migration guide v0.7.0 → v0.8.0

No changes required. New functionality is opt-in via derivative=-1.

# New in v0.8.0
x = range(0.0, 2π, length=100)
y = sin.(x)
itp = convolution_interpolation(x, y; derivative=-1)  # antiderivative

Declaration of AI assistance

Benchmarking, release notes, and documentation review performed with assistance from Claude (Anthropic).

v0.7.0 — 2026-03-06

ConvolutionInterpolations v0.7.0

Diff since v0.6.0

ConvolutionInterpolations.jl v0.7.0 Release Notes

Breaking Changes

None. All existing API defaults are preserved. The default kernel remains :b5, lazy defaults to false, and all existing keyword arguments retain their previous behavior. This release adds new features and keyword arguments with backward-compatible defaults only.

New Features

Lazy evaluation system

New lazy keyword argument for all constructors. When lazy=true, ghost point expansion is skipped at construction time — boundary ghost values are computed on the fly during evaluation only when the stencil reaches near a domain boundary. Interior evaluation is identical to eager mode.

Construction speedup scales with grid size and dimensionality: up to 349× faster for :b5 on 4D 30⁴ grids. Lazy construction is constant-time (~0.12 ms) regardless of grid size.

Implemented for both uniform and nonuniform grids: - Uniform lazy: lazy_ghost_value computes ghost points per stencil point - Nonuniform lazy: _nonuniform_lazy_coef handles per-stencil boundary logic

When lazy=true, boundary_fallback and Line() extrapolation are automatically enabled to prevent expensive tensor-product ghost computation in high dimensions.

Dimension-independent linear extrapolation

Line() extrapolation now uses a 2-point directional derivative instead of per-dimension 4th-order finite differences. Cost is 2 interpolation evaluations regardless of dimensionality — up to 9× faster than the previous approach in 4D. The old calculate_gradient and calculate_component_gradient functions have been removed.

Boundary fallback protection

New boundary_fallback keyword argument. When lazy=true, this is automatically enabled together with Line() extrapolation, so near-boundary evaluations are linearly extrapolated from the interior instead of computing expensive tensor-product ghost points. This prevents the Julia REPL from hanging on expensive corner computations in high dimensions.

Nonuniform :a0 and :a1

Nearest-neighbor (:a0) and linear (:a1) interpolation now work natively on nonuniform grids without requiring ghost points. These kernels are automatically detected and routed to the nonuniform evaluation path.

A-series derivatives

Derivatives are now supported for all a-series kernels (:a3, :a4, :a5, :a7). The top derivative order automatically switches to linear interpolation of the kernel to match the kernel’s continuity class, via the max_smooth_derivative lookup and validate_subgrid_compatibility.

Specialized dimension dispatch

Optimized evaluation paths for 1D, 2D, and 3D via Val(N), with a separate generic HigherDimension(Val(N)) path for 4D+. This enables dimension-specific optimizations without affecting the general ND fallback.

Per-dimension kernel boundary conditions

kernel_bc now accepts a tuple of (left, right) symbol pairs for per-dimension, per-side control of boundary conditions. The constructor validates that each dimension has enough grid points for the selected boundary condition; when insufficient, :polynomial silently falls back to :detect for that dimension only, preserving optimal boundary handling in all other dimensions.

Improvements

Performance

• Shipped kernel tables: Pre-shipped tables at precompute=101 mean zero computation at construction time for the default settings (already present in v0.6.0, now documented more prominently)
• Constant-time lazy construction: ~0.12 ms regardless of grid size or dimension when lazy=true
• 2-point directional extrapolation: O(1) in dimensionality, replacing O(N) finite difference gradient

Bug fixes

• :n3 now forces fast=false: Fixes a bug where explicitly passing degree=:n3 to the low-level FastConvolutionInterpolation constructor could incorrectly use precomputed tables not designed for nonuniform kernels

Tests

• Low-level constructor tests: New test coverage for ConvolutionInterpolation, FastConvolutionInterpolation, and ConvolutionExtrapolation constructors directly, across all kernels
• Tightened convergence tests: Test thresholds now match actual kernel convergence orders (e.g. error ratio ≥ 64 per grid doubling for b-series) instead of the previous conservative threshold of 2. Separate thresholds for fast vs direct evaluation paths, and per-derivative loosening accounting for ~1 order loss per derivative
• Lazy correctness tests: Verify eager vs lazy agreement at interior and boundary points for uniform and nonuniform grids in 1D, 2D, and 3D
• Lazy boundary fallback tests: Verify auto-enabling of boundary_fallback in 4D+ and correct Line() extrapolation behavior
• Nonuniform :a0/:a1 tests: Grid point reproduction, midpoint averaging, exact linear/bilinear/trilinear reproduction, and lazy field verification
• Construction speed test: Verify lazy construction of 50³ :b5 completes in under 100ms
• Test suite: 671 tests total, covering uniform/nonuniform grids, lazy/eager modes, derivatives up to 2nd order, convergence rates, extrapolation boundary conditions, and boundary fallback across 1D–5D
• Test documentation: New TESTING.md documenting the full test suite structure, strategies, and parameters

Documentation

• Updated docstrings: Comprehensive docstrings for convolution_interpolation, ConvolutionInterpolation, FastConvolutionInterpolation, ConvolutionExtrapolation, GaussianConvolutionKernel, _domain_bounds, and individual per-kernel docstrings for all a-series and b-series kernels
• Updated README: New high-dimensional interpolation section with benchmark tables, lazy mode documentation, updated kernel reference tables with nonuniform :a0/:a1, three evaluation paths in Technical Background, updated Speed section referencing lazy mode, new performance guidelines, improved Quick Start examples

Internal Changes

• Three evaluation paths: direct, fast (precomputed tables), and binomial expansion
• _domain_bounds helper accounts for lazy mode and boundary fallback
• extrapolate_point Line branch simplified from O(N) to O(1) evaluations
• validate_subgrid_compatibility enforces subgrid/derivative/kernel compatibility
• max_smooth_derivative dict maps each kernel to its maximum smooth derivative order
• Dead code removed: calculate_gradient, calculate_component_gradient

v0.6.0 — 2026-03-01

ConvolutionInterpolations v0.6.0

Diff since v0.5.1

ConvolutionInterpolations.jl v0.6.0 Release Notes

Highlights

Constructor performance improved by up to 44× through pre-shipped kernel tables and allocation-free boundary computation. The default 1D/2D setup path now requires zero disk I/O. Evaluation performance is unchanged.

New features

Pre-shipped kernel tables. Kernel tables for all b-series kernels are now shipped as package constants computed at exact rational precision. The default constructor loads these from memory, i.e. no Scratch.jl caching or disk I/O needed. Higher resolutions, BigFloat, or linear subgrid still compute on demand and cache to disk as before.

Automatic linear subgrid for 3D+. Since cubic/quintic subgrid is only implemented for 1D and 2D, the constructor now explicitly falls back to subgrid=:linear with precompute=max(precompute, 10_000) for N≥3, ensuring adequate table resolution.

Allocation-free boundary computation. The BoundaryWorkspace now achieves fully allocation-free boundary condition computation.

Function barrier pattern. Internal helpers with concrete return types eliminate type instability from the constructor.

Performance

Representative 1D benchmarks (b5 kernel, 100-point grid):

Metric	v0.5	v0.6	Speedup
Construction	95 μs, 618 allocs	5 μs, 40 allocs	19×
Evaluation	12 ns	20 ns	—

Evaluation is slightly higher than previously reported because benchmarks now include the extrapolation wrapper and the default :cubic subgrid. Bypassing the extrapolation (itp.itp(x)) and using :linear subgrid recovers lower timings. Construction speedups increase with dimensionality (up to 47× for 4D b13).

Migration guide

No breaking changes. All existing code works without modification.

Shipped tables. The default precompute=101 with :cubic subgrid now loads pre-shipped constants for 1D/2D.

3D+ automatic fallback. For N≥3, the constructor switches to subgrid=:linear, precompute=10_000. This uses the Scratch.jl cache.

POLYNOMIAL_GHOST_COEFFS stored as Float64. Previously Matrix{Int}. The integer values are exactly representable so there is no precision loss.

Internal changes

• FastConvolutionInterpolation constructor forces subgrid=:linear with precompute≥10_000 for N≥3
• apply_boundary_conditions_for_dim! refactored: ghost matrix fetched once, in-place mean-centering, polynomial/recursive path determined per side
• get_boundary_indices uses ntuple(Val(N)) for zero-allocation index construction
• fill_ghost_points_polynomial! and fill_ghost_points_recursive! accept AbstractVector{T}
• get_recursive_coefs extracted as standalone helper; old boundary_coefs removed
• Old Scratch-only get_precomputed_kernel_and_range removed
• Comprehensive docstring update across all source files

Declaration of AI assistance

Constructor optimization and docstring review performed with assistance from Claude (Anthropic).

v0.5.1 — 2026-02-28

ConvolutionInterpolations v0.5.1

Diff since v0.5.0

ConvolutionInterpolations.jl v0.5.1

Breaking changes

No breaking changes. All changes are internal bugfixes and optimizations.

Bugfixes

• Fixed polynomial boundary condition sign error for right-side boundaries in 2D+ arrays, which caused incorrect ghost point extrapolation in multidimensional nonuniform interpolation.
• Tightened derivative convergence test tolerances, enabled by the boundary fix.

Performance

• Uniform dimensions in mixed uniform/nonuniform grids now skip the Vandermonde projection during weight precomputation, using direct binomial expansion instead. This halves setup time for 2D grids with one uniform and one non-uniform dimension, with larger savings as the number of uniform dimensions increases.

v0.5.0 — 2026-02-26

ConvolutionInterpolations v0.5.0

Diff since v0.4.1

Breaking changes

No breaking changes. The API is unchanged and all changes are internal. Existing code works without modification. Non-uniform grids are detected automatically.

How it works

On uniform grids, kernel weights are the same for every interval. On non-uniform grids, each interval needs its own weights adapted to the local grid geometry. This is achieved by expanding the kernel in a binomial series (using Rational{BigInt} arithmetic to avoid floating-point contamination in the weight generation) around each interval’s local coordinate, then projecting through a Vandermonde system to enforce polynomial reproduction up to the kernel’s design order. The result is a compact set of polynomial coefficients per interval, evaluated via Horner’s method at query time.

Ghost point values near domain boundaries are also computed via Vandermonde solves using Rational{BigInt} arithmetic.

New features

• Non-uniform b-kernels: :b5 through :b13 now work on non-uniform grids with per-interval weight polynomials precomputed via exact Rational{BigInt} arithmetic
• Non-uniform derivatives: Derivatives d1 through d6 (depending on kernel continuity) on non-uniform grids. Each derivative order loses approximately one order of convergence
• :n3 kernel: Dedicated non-uniform cubic kernel equivalent to non-uniform Catmull-Rom splines. All :a kernels fall back to :n3 on non-uniform grids
• Zero-allocation evaluation: ~40ns per point via Horner’s method

Tests

414 tests (up from ~300), including non-uniform grid point reproduction (1D/2D/3D), derivative accuracy (d1/d2 in 1D/2D), convergence rate verification (d0/d1/d2), and uniform convergence tests for both fast and direct evaluation paths.

Migration guide v0.4.x → v0.5.0

Before (v0.4.x) — non-uniform grids only supported cubic interpolation:

x = [0.0, 0.3, 1.0, 1.8, 3.0, 5.0]
y = sin.(x)
itp = convolution_interpolation(x, y)  # falls back to :n3 (cubic)

After (v0.5.0) — same API, now with high-order b-kernels and derivatives:

x = [0.0, 0.3, 1.0, 1.8, 3.0, 5.0, 5.5, 6.2]
y = sin.(x)
itp = convolution_interpolation(x, y; degree=:b5)
itp_d1 = convolution_interpolation(x, y; degree=:b5, derivative=1)

Uniform grid code is completely unaffected.

v0.4.1 — 2026-02-23

ConvolutionInterpolations v0.4.1

Diff since v0.4.0

v0.4.1

Improved

• Nonuniform grid boundary handling: Ghost points at domain boundaries are now computed via closed-form Vandermonde-derived coefficients that reproduce polynomials up to degree 2, replacing the previous clamp-based approach. Recovers the uniform integer coefficients [3, -3, 1] in the uniform limit. Clean O(h³) convergence on nonuniform grids, consistent with interior accuracy.

v0.4.0 — 2026-02-22

ConvolutionInterpolations v0.4.0

Diff since v0.3.1

v0.4.0 — Nonuniform Grid Support

New Features

Nonuniform grid interpolation in any dimension. Nonuniform grids are detected automatically — the existing API works without changes:

x = [0.0, 0.3, 1.0, 1.8, 3.0, 4.0, 4.5, 5.0]
itp = convolution_interpolation(x, sin.(x))
itp(2.5)  # just works

This extends to N-dimensional tensor product grids where each dimension can have its own nonuniform spacing:

kx = [0.0, 0.3, 1.0, 1.8, 3.0]
ky = [0.0, 0.5, 1.2, 2.0, 3.5]
data = [sin(x) * cos(y) for x in kx, y in ky]
itp = convolution_interpolation((kx, ky), data)

Properties: - Cubic polynomial weights with 3rd order convergence (O(h³)) - Exact reproduction of polynomials up to degree 2 - Allocation-free evaluation (~12 ns in 1D) - Recovers Keys’ cubic kernel (a = −1/2) in the uniform limit - Mathematically equivalent to nonuniform Catmull-Rom interpolation

Adaptive refinement. Unlike Chebyshev interpolation (which requires recalculating all samples) or uniform grids (which require roughly doubling samples), nonuniform grids allow adding a single sample exactly where needed — all existing samples are reused.

Migration Guide

No breaking changes. The existing API is unchanged. Uniform grids continue to use the same optimized fast path. The nonuniform path activates automatically only when non-equal spacings are detected. No migration required. All existing code works without changes.

v0.3.1 — 2026-02-22

ConvolutionInterpolations v0.3.1

Diff since v0.3.0

v0.3.1 Release Notes

Bugfix: Pre-range type promotion causing allocations in ND evaluation

The precomputed kernel range (pre_range) was stored as Vector{Rational{BigInt}} instead of Vector{Float64}, causing type promotion to BigFloat during evaluation. This resulted in massive allocations and degraded performance in 2D+ interpolation.

Impact: 2D evaluation ~50× slower than expected, 3D-4D progressively worse. 1D unaffected due to specialized dispatch.

Fix: pre_range is now converted to the target float type before caching. Kernel precomputation still uses exact rational arithmetic internally for precision.

Action required

Clear the cached kernels to remove stale Rational{BigInt} data:

using Scratch
Scratch.clear_scratchspaces!()

Kernels will regenerate automatically on next use.

v0.3.0 — 2026-02-12

ConvolutionInterpolations v0.3.0

Diff since v0.2.0

ConvolutionInterpolations.jl v0.3.0 Release Notes

Breaking

• Renamed b3 → a4: The classic cubic 4-tap kernel (Keys 1981) now lives in the a-series where it belongs. The b prefix is reserved exclusively for the novel kernel family discovered in this package.

New features

• Derivative interpolation: Compute derivatives through analytically predifferentiated convolution kernels — no finite differences, no step-size tuning. Supports up to 6th order depending on kernel continuity class.
• Subgrid interpolation modes: Cubic and quintic Hermite subgrid interpolation using predifferentiated kernels. Achieves high accuracy with as few as 100 precomputed points, dramatically reducing memory while improving quality over linear subgrid.

Corrections

• Continuity classes rigorously verified: b5 (C³), b7 (C⁴), b9 (C⁵), b11 (C⁶), b13 (C⁶). Previous claims of higher continuity did not hold up under verification. C⁶ is still far beyond the C¹ that literature kernels actually achieve.

---

Migration guide: v0.2.0 → v0.3.0

Kernel rename: b3 → a4

This is the only breaking change. If you used the b3 kernel, replace it with a4:

# Before (v0.2.0)
itp = convolution_interpolation(x, y; degree=:b3)

# After (v0.3.0)
itp = convolution_interpolation(x, y; degree=:a4)

This rename reflects the kernel’s origin: a4 is the 4-tap cubic kernel from the literature (Keys 1981), not a member of the novel b-series family. The kernel itself is unchanged — only the name is different.

No other breaking changes

All other API remains the same. The new derivative and subgrid keyword arguments are additive and default to their v0.2.0 behavior when not specified:

# These calls work identically in v0.2.0 and v0.3.0
itp = convolution_interpolation(x, y)
itp = convolution_interpolation(x, y; degree=:b5)
itp = convolution_interpolation((xs, ys), data; kernel_bc=:polynomial)

New features available after upgrading

# Derivatives (new in v0.3.0)
itp_d1 = convolution_interpolation(x, y; derivative=1)   # first derivative
itp_d2 = convolution_interpolation(x, y; derivative=2)   # second derivative

# Subgrid modes (new in v0.3.0)
itp = convolution_interpolation(x, y; subgrid=:cubic)    # default, uses 1st derivative
itp = convolution_interpolation(x, y; subgrid=:quintic)  # uses 1st and 2nd derivatives
itp = convolution_interpolation(x, y; subgrid=:linear)   # v0.2.0 behavior

v0.2.0 — 2025-11-20

ConvolutionInterpolations v0.2.0

Diff since v0.1.4

ConvolutionInterpolations.jl - Significant Release

This release represents a significant evolution of ConvolutionInterpolations.jl, featuring a breakthrough in boundary condition handling that unlocks the full potential of the existing b-series kernels, along with substantial performance optimizations. While the main interface remains largely unchanged, several defaults and internal behaviors have been updated for optimal performance.

🎯 Highlights

• Novel polynomial boundary conditions: Vandermonde-based ghost point computation enabling b-series kernels to achieve their theoretical 7th order convergence all the way to domain edges (previously these kernels existed but couldn’t reach their design accuracy)
• Default kernel changed: Now defaults to :b5 (quintic reproduction, C³ continuity) instead of :a3
• True O(1) performance: Elimination of binary search operations for constant-time interpolation
• Persistent kernel caching: First-use generation with disk caching via Scratch.jl for instant subsequent loads
• Extended precision support: Kernel coefficients stored as exact rationals enabling BigFloat arithmetic
• Allocation-free evaluation: All fast-mode evaluations allocate zero memory
• Kernel discovery tools: New gen/ folder with scripts for deriving kernels and boundary conditions

🔄 Breaking Changes

Default Kernel Changed

Before: Default was :a3 (cubic, C¹ continuous, 4th order accurate) After: Default is :b5 (quintic, C³ continuous, 7th order accurate)

Migration: If you require the old behavior, explicitly specify degree=:a3:

# Old code (implicit default)
itp = convolution_interpolation(x, y)

# New code for same behavior
itp = convolution_interpolation(x, y; degree=:a3)

Rationale: The b5 kernel provides dramatically better accuracy (7th order vs 4th order) with comparable performance. With the new polynomial boundary conditions, it can now achieve its designed accuracy, making it a superior default for most use cases.

Boundary Condition Default Changed

Before: Default was kernel_bc=:detect After: Default is kernel_bc=:auto

Migration: No action needed in most cases. The :auto option intelligently: - Prioritizes :polynomial boundary conditions when sufficient grid points exist - Falls back to :detect for sparse grids or compatibility

Rationale: The new :polynomial boundary condition (computed via Vandermonde systems) preserves polynomial reproduction properties at boundaries, ensuring 7th order convergence extends to domain edges.

Precompute Default Increased

Before: Default was precompute=1_000 (1,000 points) After: Default is precompute=100_000 (100,000 points)

Migration: For memory-constrained applications, explicitly reduce:

itp = convolution_interpolation(x, y; precompute=10_000)

Rationale: The default now provides maximum accuracy with negligible memory cost (~1.6 MB per kernel). Since kernel lookups are O(1) regardless of precompute resolution, there’s no performance penalty for the higher default.

Kernel Coefficients Now Type-Generic

Before: Kernel coefficients were stored as Float64 After: Coefficients stored as exact Rational{Int} and converted to evaluation type

Impact: This is transparent for Float64 usage but enables extended precision:

# Now supports arbitrary precision
x_big = big.(x)
y_big = big.(y)
itp = convolution_interpolation(x_big, y_big)  # Uses BigFloat precision

Migration: No changes needed. Existing Float64 code works identically.

New Dependencies

Before: Zero external dependencies After: Minimal dependencies added for caching functionality - LinearAlgebra: Standard library (for Vandermonde systems) - Serialization: Standard library (for kernel caching) - Scratch.jl: External package (for persistent disk cache management)

Migration: No code changes needed. Scratch.jl will be installed automatically as a dependency.

Rationale: These minimal dependencies enable persistent kernel caching (eliminating redundant computation) and the new polynomial boundary condition method.

✨ New Features

Breakthrough: Polynomial Boundary Conditions

The major advancement in this release: a Vandermonde-based boundary condition method that enables the existing b-series kernels to achieve their theoretical convergence properties.

Background: The b-series kernels (b5, b7, b9, b11, b13) have existed in the package since their discovery via systematic analytical search. However, they suffered from poor boundary handling that prevented them from achieving their designed 7th order convergence.

The Breakthrough: Novel :polynomial boundary condition method: - Computes optimal ghost point values by solving Vandermonde systems - Preserves each kernel’s polynomial reproduction properties at domain boundaries - Ensures 7th order convergence extends all the way to domain edges - Automatically selected by :auto when sufficient grid points available

Impact: The b-series kernels can now achieve their designed accuracy on the challenging Runge function:

Kernel	Polynomial Reproduction	Continuity	Order	Runge Accuracy (1000 pts, Float64)
:b5	Quintic (degree 5)	C³	7th	~10⁻¹⁴
:b7	Septic (degree 7)	C⁵	7th	~10⁻¹²
:b9	Septic (degree 7)	C⁷	7th	~10⁻¹⁰
:b11	Septic (degree 7)	C⁹	7th	~10⁻⁹
:b13	Septic (degree 7)	C¹¹	7th	~10⁻⁸

Note: Higher-degree kernels are more affected by Float64 machine precision limitations. All kernels continue converging beyond Float64 precision when using BigFloat arithmetic.

# Explicitly use polynomial boundary conditions
itp = convolution_interpolation(x, y; kernel_bc=:polynomial)

Kernel Discovery Tools

New gen/ folder containing scripts for kernel and boundary condition discovery:

1. generalized_kernel_derivation.py: Systematic kernel discovery via symbolic boundary value problems
   • Generalizes Keys (1981) method to arbitrary polynomial degrees
   • Exhaustive search through boundary condition combinations
   • Taylor series verification of convergence order
   • Frequency domain optimization via FFT analysis
   • Can reproduce and verify published kernels (Keys cubic, Meijering quintic/septic)
2. compute_ghost_coefficients.jl: Boundary condition computation
   • Solves Vandermonde systems for optimal ghost point coefficients
   • Ensures polynomial reproduction at domain boundaries
   • Generates coefficients for any kernel degree
   • Outputs exact rational coefficients for arbitrary precision

These tools enable researchers to: - Discover new convolution kernels for specific applications - Verify published kernel properties - Compute optimal boundary conditions for custom kernels - Extend the method to specialized interpolation requirements

Persistent Kernel Caching

Discretized kernels are computed once and cached to disk using Scratch.jl: - First use: 0.1-10 seconds to generate (depends on kernel order and precompute size) - Subsequent uses: <1ms to load from disk cache - Precision-aware: Separate caches for Float64, Float32, BigFloat(256), etc. - Automatic: No user intervention required

Extended Precision Support

All kernels support arbitrary precision arithmetic:

# Use BigFloat for accuracy beyond Float64 machine precision
setprecision(256)  # 256-bit precision
x = range(big(0.0), big(2π), length=100)
y = sin.(x)
itp = convolution_interpolation(x, y)  # Automatically uses BigFloat

Kernel coefficients stored as exact rationals enable convergence well beyond standard machine precision.

🚀 Performance Improvements

True O(1) Interpolation

• Eliminated all binary search operations
• Direct index calculation on uniform grids
• Interpolation time independent of grid size (verified: same ~15ns from 10 to 10,000 points)

Allocation-Free Evaluation

All fast-mode kernel evaluations are allocation-free:

# Zero allocations per evaluation
@allocated itp(0.5)  # Returns 0

Kernel Discretization with Linear Interpolation

Novel implementation approach: - Discretizes kernels at high resolution (default 100,000 points) - Uses linear interpolation between precomputed values - More numerically stable than direct polynomial evaluation for practical kernel degrees - Faster than Horner’s method while maintaining accuracy

Benchmark Results (1D)

Kernel	Initialization	Evaluation	Accuracy (Runge, 1000 pts, Float64)
:a1	~75 μs	4 ns	~10⁻³
:b5	~1.9 ms	15 ns	~10⁻¹⁴
:b7	~2.4 ms	16 ns	~10⁻¹²
:b9	~2.8 ms	17 ns	~10⁻¹⁰
:b11	~3.2 ms	18 ns	~10⁻⁹
:b13	~3.5 ms	19 ns	~10⁻⁸

Note: Higher-degree kernels (b9-b13) are progressively more affected by Float64 machine precision limitations. For applications requiring their full accuracy potential, use BigFloat arithmetic - all kernels continue converging beyond Float64 precision.

📚 API Additions

Kernel Options (Clarification)

The following kernels were already available in the package:

Low-order kernels:

:a0, :a1, :a3, :b3

Higher-order a-series (from literature):

:a5, :a7

B-series kernels (7th order convergence):

:b5, :b7, :b9, :b11, :b13

Note: These kernels existed in previous versions but could not achieve their designed accuracy due to boundary condition limitations. With the new polynomial boundary conditions, they now work as intended.

Gaussian smoothing:

B=1.0  # Existing

New Boundary Condition Options

# Automatic (prioritizes :polynomial, falls back to :detect)
kernel_bc=:auto  # ⭐ New default

# Polynomial reproduction (optimal)
kernel_bc=:polynomial  # ⭐ New

# Existing options
kernel_bc=:detect, :linear, :quadratic, :periodic

🔧 Migration Guide

Minimal Migration (Keep Existing Behavior)

If you want to maintain exact previous behavior:

# Change from:
itp = convolution_interpolation(x, y)

# To:
itp = convolution_interpolation(x, y; degree=:a3, kernel_bc=:detect, precompute=1_000)

Recommended Migration (Adopt New Defaults)

Most users should simply update their code to use the new defaults, which provide better accuracy and performance:

# Old code
itp = convolution_interpolation(x, y)  # Was :a3

# New code - no changes needed, now uses :b5!
itp = convolution_interpolation(x, y)

The b5 kernel with polynomial boundary conditions provides superior accuracy (7th order vs 4th order) with comparable performance to the old cubic default.

Memory-Constrained Environments

If memory is extremely limited (embedded systems, etc.):

# Reduce precompute resolution (degrades accuracy slightly)
itp = convolution_interpolation(x, y; precompute=10_000)

# Or disable fast mode entirely (slower but minimal memory)
itp = convolution_interpolation(x, y; fast=false)

High-Precision Applications

For applications requiring extended precision:

# Use BigFloat with desired precision
setprecision(512)  # 512-bit precision
x = big.(your_data_x)
y = big.(your_data_y)
itp = convolution_interpolation(x, y)  # Automatically uses BigFloat

📊 Comparison with Previous Version

Aspect	Previous	Current	Improvement
Default kernel	:a3 (cubic)	:b5 (quintic)	500× better accuracy
B-kernel functionality	Existed but inaccurate	Fully functional	Achieves designed 7th order
Convergence order	4th	7th	1.75× order improvement
Boundary handling	Detect-based	Polynomial-optimal	Maintains accuracy at edges
Evaluation performance	Slower, some type instabilities	~4-19ns (1D), O(1) guaranteed	Major performance upgrade
Memory allocations	Variable, not allocation-free	Zero (fast mode)	Allocation-free
Initialization	Per-call overhead	Cached after first use	~1000× faster subsequent uses
Precision support	Float64 only	Any (Float32, Float64, BigFloat)	Unlimited precision
Dependencies	Zero	Minimal (stdlib + Scratch.jl)	Enables caching & boundaries

🐛 Bug Fixes

• Fixed binary search scaling in higher dimensions
• Corrected boundary coefficient computation for periodic signals
• Improved numerical stability of kernel evaluation near boundaries
• Fixed type promotion issues with mixed-precision data

🔬 Technical Details

Kernel Discovery Method

The b-series kernels were discovered through systematic analytical search using SymPy: 1. Formulate boundary value problem as symbolic equations 2. Exhaustively search boundary condition combinations 3. Verify polynomial reproduction via Taylor series expansion 4. Validate frequency domain characteristics via FFT analysis 5. Select kernels with optimal accuracy/smoothness balance

This automated discovery process generalizes the approach of Keys (1981) to arbitrary polynomial degrees, enabling kernels that would be intractable to derive by hand.

Exact Rational Storage

Kernel coefficients stored as Rational{Int}:

# Example b5 kernel coefficient
-1228311859//656862720  # Exact rational, not Float64

This enables: - Arbitrary precision evaluation - No accumulated rounding errors in coefficient storage - Reproducible results across platforms and precision levels

📖 References

Novel Contributions (This Work)

• Polynomial boundary condition method via Vandermonde systems (major breakthrough)
• Kernel discretization with linear interpolation approach
• O(1) performance via direct indexing on uniform grids
• Persistent kernel caching via Scratch.jl
• B-series kernel discovery via systematic analytical search (kernels existed in package previously, but discovery method is novel)

Foundational Work

• R.G. Keys, “Cubic Convolution Interpolation for Digital Image Processing,” IEEE Trans. ASSP, 1981
• E.H.W. Meijering et al., “Image Reconstruction by Convolution with Symmetrical Piecewise nth-Order Polynomial Kernels,” IEEE Trans. IP, 1999

🙏 Acknowledgments

Special thanks to the Interpolations.jl maintainers for API design inspiration and to the Julia community for feedback during development.

📝 License

MIT License - see LICENSE file for details

---

Questions or Issues? - GitHub Issues: [ConvolutionInterpolations.jl/issues] - Documentation: [Full package documentation] - Examples: See README.md for comprehensive examples

Upgrading from v0.x? Most code will work without changes. For exact previous behavior, see the “Minimal Migration” section above.

Merged pull requests: - Bump actions/checkout from 4 to 5 (#5) (@dependabot [bot])

v0.1.4 — 2025-08-03

ConvolutionInterpolations v0.1.4

Diff since v0.1.3

v0.1.3 — 2025-08-02

ConvolutionInterpolations v0.1.3

Diff since v0.1.2

v0.1.2 — 2025-05-22

ConvolutionInterpolations v0.1.2

Diff since v0.1.1

Closed issues: - Support for non-uniform grid? (#3)

v0.1.1 — 2025-04-11

ConvolutionInterpolations v0.1.1

Diff since v0.1.0

v0.1.0 — 2025-04-10

ConvolutionInterpolations v0.1.0