Local Cubic Hermite API

hermite_interp(x, y, dy, xq; extrap=NoExtrap(), deriv=EvalValue(), search=AutoSearch(), hint=nothing)

Cubic Hermite interpolation at a single query point using user-supplied slopes.

Returns interpolated value (or derivative, if deriv is set). C$^1$ continuous — slopes are used directly, no global spline solve.

hermite_interp(x, y, dy, x_query; extrap=NoExtrap(), deriv=EvalValue(), search=AutoSearch(), hint=nothing)

Cubic Hermite interpolation at multiple query points using user-supplied slopes. Returns Vector{Tv} of interpolated values.

hermite_interp(x, y, dy; extrap=NoExtrap(), search=AutoSearch()) -> CubicHermiteInterpolant1D

Create a callable cubic Hermite interpolant from user-supplied slopes.

Arguments

• x::AbstractVector: x-coordinates (must be sorted)
• y::AbstractVector: Function values at grid points
• dy::AbstractVector: First derivatives at grid points
• extrap::AbstractExtrap: Extrapolation mode (default: NoExtrap())
• search::AbstractSearchPolicy: Default search policy (default: AutoSearch())

Returns

CubicHermiteInterpolant1D callable object. Supports:

• Scalar: itp(0.5)
• Vector: itp(xq_vec)
• In-place: itp(out, xq_vec)
• Broadcast: itp.(xq) or @. coef * itp(rho)
• Derivative: itp(0.5; deriv=DerivOp(1))

Example

x  = range(0, 2π, 100)
y  = sin.(x)
dy = cos.(x)
itp = hermite_interp(x, y, dy)
itp(1.0)                          # ≈ sin(1.0)
itp(1.0; deriv=DerivOp(1))       # ≈ cos(1.0)

hermite_interp!(output, x, y, dy, x_query; extrap=NoExtrap(), deriv=EvalValue(), search=AutoSearch(), hint=nothing)

In-place cubic Hermite interpolation using user-supplied slopes.

pchip_interp(x, y, xq; coeffs=PreCompute(), extrap=NoExtrap(), deriv=EvalValue(), search=AutoSearch(), hint=nothing)

PCHIP (monotone-preserving) interpolation at a single query point.

Computes Fritsch-Carlson slopes internally, then evaluates via cubic Hermite kernel. C$^1$ continuous, monotonicity guaranteed for monotone input data.

Keyword Arguments

• coeffs: Coefficient strategy — PreCompute() (bulk slopes) or OnTheFly() (local slopes per cell)

pchip_interp(x, y, x_query; coeffs=PreCompute(), ...)

PCHIP interpolation at multiple query points. Returns Vector{Tv}.

pchip_interp(x, y; extrap=NoExtrap(), search=AutoSearch()) -> PchipInterpolant1D

Create a callable PCHIP interpolant with monotone-preserving slopes.

Slopes are computed once at construction via the Fritsch-Carlson algorithm. Subsequent itp(xq) calls just evaluate — zero slope recomputation.

Arguments

• x::AbstractVector: x-coordinates (must be sorted, ≥2 points)
• y::AbstractVector: y-values
• extrap::AbstractExtrap: Extrapolation mode (default: NoExtrap())
• search::AbstractSearchPolicy: Default search policy (default: AutoSearch())

Returns

PchipInterpolant1D callable object. Supports:

• Scalar: itp(0.5)
• Vector: itp(xq_vec)
• In-place: itp(out, xq_vec)
• Broadcast: itp.(xq) or @. coef * itp(rho)
• Derivative: itp(0.5; deriv=DerivOp(1))

Example

x = [0.0, 1.0, 2.0, 3.0, 4.0]
y = [0.0, 1.0, 0.5, 0.8, 0.2]  # non-monotone data
itp = pchip_interp(x, y)
itp(1.5)                          # monotone within each monotone segment
itp(1.5; deriv=DerivOp(1))       # first derivative

pchip_interp!(output, x, y, x_query; coeffs=PreCompute(), ...)

In-place PCHIP interpolation with monotone-preserving slopes.

cardinal_interp(x, y, xq; coeffs=AutoCoeffs(), tension=0.0, extrap=NoExtrap(), deriv=EvalValue(), search=AutoSearch(), hint=nothing)

Cardinal spline interpolation at a single query point. Default tension=0 is Catmull-Rom. C$^1$ continuous.

Keyword Arguments

• coeffs: Coefficient strategy — PreCompute() (bulk slopes) or OnTheFly() (local slopes per cell)
• tension: Cardinal tension parameter (0 = Catmull-Rom, 1 = linear)

cardinal_interp(x, y, x_query; coeffs=AutoCoeffs(), tension=0.0, ...)

Cardinal spline interpolation at multiple query points. Returns Vector{Tv}.

cardinal_interp(x, y; tension=0.0, extrap=NoExtrap(), search=AutoSearch()) -> CardinalInterpolant1D

Create a callable cardinal spline interpolant.

Default tension=0 gives Catmull-Rom (central finite difference slopes). Increasing tension reduces overshoot; tension=1 gives zero slopes at knots (smooth S-curves between knots).

Example

itp = cardinal_interp(x, y)                # CatmullRom (tension=0)
itp = cardinal_interp(x, y; tension=0.5)   # tighter curve
itp(0.5)

cardinal_interp!(output, x, y, x_query; coeffs=AutoCoeffs(), tension=0.0, ...)

In-place cardinal spline interpolation.

akima_interp(x, y, xq; coeffs=PreCompute(), extrap=NoExtrap(), deriv=EvalValue(), search=AutoSearch(), hint=nothing)

Akima interpolation at a single query point.

Outlier-robust, C$^1$ continuous.

Keyword Arguments

• coeffs: Coefficient strategy — PreCompute() (bulk slopes) or OnTheFly() (local slopes per cell)

akima_interp(x, y, x_query; coeffs=PreCompute(), ...)

Akima interpolation at multiple query points. Returns Vector{Tv}.

akima_interp(x, y; extrap=NoExtrap(), search=AutoSearch()) -> AkimaInterpolant1D

Create a callable Akima interpolant with outlier-robust slopes.

Slopes are computed once at construction. Requires ≥ 2 grid points.

Example

itp = akima_interp(x, y)
itp(0.5)
itp(0.5; deriv=DerivOp(1))

akima_interp!(output, x, y, x_query; coeffs=PreCompute(), ...)

In-place Akima interpolation with outlier-robust slopes.

AbstractInterpolant1D{Tg<:AbstractFloat, Tv} <: AbstractInterpolant{Tg, Tv}

Abstract supertype for 1-dimensional interpolant objects.

All four 1D interpolant types inherit from this, enabling shared callable dispatch via the interpolant_protocol.jl interface.

Protocol (interpolant_protocol.jl)

Subtypes automatically inherit 3 callable overloads (scalar, vector-alloc, vector-inplace) by implementing the required interface:

_itp_eval_scalar(itp, xq, extrap, op, searcher)   — core scalar evaluation (required)
_itp_vector_loop!(out, itp, xq, extrap, op, searcher) — core vector loop (required)
_itp_grid(itp)                  — grid vector access (default: itp.x)
_itp_extrap(itp)                — extrapolation mode (default: itp.extrap)
_itp_search(itp)                — default search policy (default: itp.search_policy)

Subtypes

• LinearInterpolant{Tg, Tv}: Piecewise linear interpolation
• ConstantInterpolant{Tg, Tv}: Piecewise constant (step) interpolation
• QuadraticInterpolant{Tg, Tv}: C1 piecewise quadratic spline
• CubicInterpolant{Tg, Tv}: C2 cubic spline
• AbstractHermiteInterpolant1D{Tg, Tv}: Cubic Hermite family (see subtypes below)

AbstractHermiteInterpolant1D{Tg, Tv} <: AbstractInterpolant1D{Tg, Tv}

Abstract supertype for 1D cubic Hermite family interpolants.

All subtypes store (x, y, dy, spacing, extrap, search_policy) and share the same evaluation kernel (_hermite_eval_at_point) and integration kernel (_hermite_integral_kernel_1d). This enables a single dispatch point for _itp_eval_scalar, _itp_vector_loop!, and integrate.

Subtypes

• CubicHermiteInterpolant1D: User-supplied slopes
• PchipInterpolant1D: Fritsch-Carlson monotone-preserving slopes
• CardinalInterpolant1D: Central finite difference (+ tension parameter)
• AkimaInterpolant1D: Weighted-average outlier-robust slopes

CubicHermiteInterpolant1D{Tg, Tv, X, Y, DY, S, E, P}

Callable interpolant for cubic Hermite interpolation with user-supplied slopes. Returned by hermite_interp(x, y, dy) (interpolant form).

Stores precomputed grid spacing for O(1) h/inv_h lookup on uniform grids. Evaluation uses _hermite_kernel_1d (derivative-based Hermite basis functions), NOT _cubic_kernel (moment-based spline formulation).

Properties

• C$^1$ continuous (continuous first derivative, discontinuous second derivative at knots)
• User-supplied slopes: dy is the first derivative at grid points — no global solve, no BC
• Local: each cell depends only on its 2 endpoints — change 1 data point, only 2 cells affected

Usage

x  = range(0, 2π, 100)
y  = sin.(x)
dy = cos.(x)
itp = hermite_interp(x, y, dy)
itp(0.5)                          # scalar query
itp(xq_vec)                       # vector query (allocating)
itp(out, xq_vec)                  # vector query (in-place)
itp(0.5; deriv=DerivOp(1))       # first derivative
itp(0.5; search=BinarySearch())  # override search

PchipInterpolant1D{Tg, Tv, X, Y, DY, S, E, P}

Callable interpolant for PCHIP (Piecewise Cubic Hermite Interpolating Polynomial). Returned by pchip_interp(x, y) (2-argument form).

Slopes are computed once at construction via the Fritsch-Carlson monotone-preserving algorithm. Evaluation uses the same cubic Hermite kernel as CubicHermiteInterpolant1D.

Properties

• C$^1$ continuous (continuous first derivative, discontinuous second derivative at knots)
• Monotonicity preserving: monotone input data → monotone interpolant
• Local: each slope depends only on neighboring data — no global solve

Usage

itp = pchip_interp(x, y)
itp(0.5)                          # scalar query
itp(xq_vec)                       # vector query (allocating)
itp(out, xq_vec)                  # vector query (in-place)
itp(0.5; deriv=DerivOp(1))       # first derivative

CardinalInterpolant1D{Tg, Tv, X, Y, DY, S, E, P}

Callable interpolant for cardinal spline interpolation. Returned by cardinal_interp(x, y) (2-argument form).

Slopes are computed once at construction via: d_k = (1 - tension) * (y[k+1] - y[k-1]) / (x[k+1] - x[k-1])

When tension=0, this is the Catmull-Rom spline (simple central FD).

Properties

• C$^1$ continuous (continuous first derivative)
• Local: each slope depends only on immediate neighbors
• Tension parameter: controls overshoot (0 = CatmullRom, 1 = zero slopes at knots)

Usage

itp = cardinal_interp(x, y)                     # default: CatmullRom (tension=0)
itp = cardinal_interp(x, y; tension=0.5)        # tighter curve
itp(0.5)
itp(0.5; deriv=DerivOp(1))

AkimaInterpolant1D{Tg, Tv, X, Y, DY, S, E, P}

Callable interpolant for Akima interpolation. Returned by akima_interp(x, y) (2-argument form).

Slopes are computed once at construction via the Akima (1970) weighted-average algorithm. Evaluation uses the same cubic Hermite kernel.

Properties

• C$^1$ continuous (continuous first derivative)
• Local: each slope depends on 4 adjacent secants (5-point stencil)
• Outlier-robust: deviant secants receive less weight

Usage

itp = akima_interp(x, y)
itp(0.5)
itp(0.5; deriv=DerivOp(1))

hermite_adjoint(x, x_query; extrap=NoExtrap()) -> HermiteAdjoint1D

Create a cubic Hermite adjoint operator for user-supplied slopes (query-baked, data-free).

Computes f̄ = Wᵀȳ where W is the forward Hermite interpolation weight matrix restricted to the y-dependence. User-supplied slopes dy are treated as fixed and do not receive gradient contributions.

This adjoint satisfies the dot-product identity:

dot(hermite_interp(x, y, dy, xq), ȳ) ≈ dot(y, hermite_adjoint(x, xq)(ȳ))

For PCHIP/Cardinal/Akima where slopes are computed from y, use their dedicated adjoint constructors instead (which scatter to both y and dy).

Arguments

• x::AbstractVector: Grid points (must be sorted)
• x_query::AbstractVector: Query points (baked into the operator)
• extrap::AbstractExtrap: Extrapolation mode (default: NoExtrap())

Example

using LinearAlgebra
x = collect(range(0, 1, 50))
xq = sort(rand(30))
y  = sin.(x)
dy = cos.(x)
y_bar = randn(30)

itp = hermite_interp(x, y, dy)
adj = hermite_adjoint(x, xq)
f_bar = adj(y_bar)

# Dot-product identity: ⟨W·y, ȳ⟩ = ⟨y, Wᵀȳ⟩
@assert dot(itp.(xq), y_bar) ≈ dot(y, adj(y_bar))

Notes

• No boundary conditions needed (Hermite uses supplied slopes directly).
• Pure scatter operation (no tridiagonal solve).
• 4th+ derivative adjoints always return zero (cubic polynomial).

pchip_adjoint(x, y, x_query; extrap=NoExtrap()) -> PchipAdjoint1D

Create a PCHIP adjoint operator (query-baked, linearized at data y).

Computes f_bar = W^T * y_bar where W is the forward PCHIP interpolation weight matrix, including the slope-from-data dependence.

Because PCHIP slopes involve data-dependent branching (monotonicity clamping), the adjoint is linearized at the given y. The dot-product identity holds exactly at this operating point:

dot(pchip_interp(x, y).(xq), y_bar) == dot(y, pchip_adjoint(x, y, xq)(y_bar))

Arguments

• x::AbstractVector: Grid points (must be sorted)
• y::AbstractVector: Data values (determines slope branch conditions)
• x_query::AbstractVector: Query points (baked into the operator)
• extrap::AbstractExtrap: Extrapolation mode (default: NoExtrap())

Example

using LinearAlgebra
x = collect(range(0, 1, 50))
y = randn(50)
xq = sort(rand(30))
y_bar = randn(30)

itp = pchip_interp(x, y)
adj = pchip_adjoint(x, y, xq)
f_bar = adj(y_bar)

# Dot-product identity: <W*y, y_bar> = <y, W^T*y_bar>
@assert dot(itp.(xq), y_bar) ≈ dot(y, adj(y_bar))

cardinal_adjoint(x, x_query; tension=0.0, extrap=NoExtrap()) -> CardinalAdjoint1D

Create a cardinal spline adjoint operator (query-baked, data-free).

Computes f_bar = W^T * y_bar where W is the forward cardinal interpolation weight matrix, including the slope-from-data dependence.

Arguments

• x::AbstractVector: Grid points (must be sorted)
• x_query::AbstractVector: Query points (baked into the operator)
• tension::Real: Tension parameter (0 = CatmullRom, 1 = zero slopes)
• extrap::AbstractExtrap: Extrapolation mode (default: NoExtrap())

Example

using LinearAlgebra
x = collect(range(0, 1, 50))
xq = sort(rand(30))
y = randn(50)
y_bar = randn(30)

itp = cardinal_interp(x, y; tension=0.5)
adj = cardinal_adjoint(x, xq; tension=0.5)
f_bar = adj(y_bar)

# Dot-product identity: <W*y, y_bar> = <y, W^T*y_bar>
@assert dot(itp.(xq), y_bar) ≈ dot(y, adj(y_bar))

akima_adjoint(x, y, x_query; extrap=NoExtrap()) -> AkimaAdjoint1D

Create an Akima adjoint operator (query-baked, linearized at data y).

Computes f_bar = W^T * y_bar where W is the forward Akima interpolation weight matrix, including the slope-from-data dependence.

Because Akima slopes involve data-dependent branching (absolute-value weights), the adjoint is linearized at the given y. The dot-product identity holds exactly at this operating point:

dot(akima_interp(x, y).(xq), y_bar) == dot(y, akima_adjoint(x, y, xq)(y_bar))

Arguments

• x::AbstractVector: Grid points (must be sorted)
• y::AbstractVector: Data values (determines slope weight conditions)
• x_query::AbstractVector: Query points (baked into the operator)
• extrap::AbstractExtrap: Extrapolation mode (default: NoExtrap())

Example

using LinearAlgebra
x = collect(range(0, 1, 50))
y = randn(50)
xq = sort(rand(30))
y_bar = randn(30)

itp = akima_interp(x, y)
adj = akima_adjoint(x, y, xq)
f_bar = adj(y_bar)

# Dot-product identity: <W*y, y_bar> = <y, W^T*y_bar>
@assert dot(itp.(xq), y_bar) ≈ dot(y, adj(y_bar))

HermiteAdjoint1D{Tg, EP}

Adjoint (transpose) operator for cubic Hermite interpolation with user-supplied slopes. Computes f̄ = Wᵀȳ where W is the forward Hermite interpolation weight matrix, considering only the y-dependence (slopes dy are treated as fixed).

Constructed from a grid and query points (query-baked, data-free).

Type Parameters

• Tg: Grid float type (Float32 or Float64)
• EP: Extrapolation policy type (NoExtrap, ExtendExtrap, ClampExtrap, FillExtrap, WrapExtrap)

Usage

adj = hermite_adjoint(x, xq)

# Value adjoint (default)
f_bar = adj(y_bar)

# Derivative adjoint
f_bar = adj(y_bar; deriv=DerivOp(1))

# In-place
adj(f_bar, y_bar; deriv=DerivOp(1))

# Dot-product identity: ⟨W·y, ȳ⟩ = ⟨y, Wᵀȳ⟩
@assert dot(itp.(xq), y_bar) ≈ dot(y, adj(y_bar))

PchipAdjoint1D{Tg, Tv, EP}

Adjoint (transpose) operator for PCHIP (Fritsch-Carlson) interpolation. Computes f_bar = W^T * y_bar where W is the forward PCHIP interpolation weight matrix, including the slope-from-data dependence.

Because PCHIP slopes involve data-dependent branching (sign checks for monotonicity clamping), this adjoint is constructed at a specific y operating point. The dot-product identity holds exactly at that y:

dot(pchip_interp(x, y).(xq), y_bar) == dot(y, adj(y_bar))

Type Parameters

• Tg: Grid float type (Float32 or Float64)
• Tv: Value type (must support sign, abs, zero, division)
• EP: Extrapolation policy type

Usage

adj = pchip_adjoint(x, y, xq)

f_bar = adj(y_bar)                      # value adjoint
f_bar = adj(y_bar; deriv=DerivOp(1))    # derivative adjoint
adj(f_bar, y_bar)                       # in-place

CardinalAdjoint1D{Tg, EP}

Adjoint (transpose) operator for cardinal spline interpolation. Computes f_bar = W^T * y_bar where W is the forward cardinal interpolation weight matrix, including the slope-from-data dependence.

Because cardinal slopes are linear in y, the full forward operation is linear and the dot-product identity holds exactly:

dot(cardinal_interp(x, y; tension=t).(xq), y_bar) == dot(y, adj(y_bar))

Type Parameters

• Tg: Grid float type (Float32 or Float64)
• EP: Extrapolation policy type

Usage

adj = cardinal_adjoint(x, xq; tension=0.5)

f_bar = adj(y_bar)                      # value adjoint
f_bar = adj(y_bar; deriv=DerivOp(1))    # derivative adjoint
adj(f_bar, y_bar)                       # in-place

AkimaAdjoint1D{Tg, Tv, EP}

Adjoint (transpose) operator for Akima interpolation. Computes f_bar = W^T * y_bar where W is the forward Akima interpolation weight matrix, including the slope-from-data dependence.

Because Akima slopes involve data-dependent branching (absolute-value weights on secant differences), this adjoint is constructed at a specific y operating point. The dot-product identity holds exactly at that y:

dot(akima_interp(x, y).(xq), y_bar) == dot(y, adj(y_bar))

Type Parameters

• Tg: Grid float type (Float32 or Float64)
• Tv: Value type (must support sign, abs, zero, division)
• EP: Extrapolation policy type

Usage

adj = akima_adjoint(x, y, xq)

f_bar = adj(y_bar)                      # value adjoint
f_bar = adj(y_bar; deriv=DerivOp(1))    # derivative adjoint
adj(f_bar, y_bar)                       # in-place