Adjoint API

cubic_adjoint(x, x_query; bc=CubicFit(), autocache=true) -> CubicAdjoint

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

Computes f̄ = Wᵀȳ where W is the forward interpolation weight matrix. Maps query-space sensitivities back to grid-space sensitivities.

Arguments

• x::AbstractVector: Grid points (must be sorted)
• x_query::AbstractVector: Query points (baked into the operator)
• bc::AbstractBC: Boundary condition (default: CubicFit())
• extrap::AbstractExtrap: Extrapolation mode (default: NoExtrap()). For WrapExtrap(), queries are wrapped to the domain before anchoring. For ClampExtrap()/FillExtrap(), queries are clamped to the boundary.
• autocache::Bool: Enable automatic caching (default: true)

Example

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

# Forward
itp = cubic_interp(x, f; bc=CubicFit())

# Adjoint
adj = cubic_adjoint(x, xq; bc=CubicFit())
f_bar = adj(y_bar)                      # value adjoint
f_bar = adj(y_bar; deriv=DerivOp(1))    # 1st derivative adjoint
adj(f_bar, y_bar; deriv=DerivOp(1))     # in-place

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

Notes

• PeriodicBC() and PeriodicBC(endpoint=:exclusive) are supported. For periodic BCs, the adjoint uses Sherman-Morrison with the same symmetric factorization.
• The adjoint is the Jacobian transpose (∂y/∂f)ᵀ, NOT an inverse operator.
• For BCs with non-zero values (e.g., Deriv1(0.5)), the forward is affine: y = W·f + c. The adjoint computes Wᵀȳ, independent of the constant c.

cubic_adjoint(grids::NTuple{N}, queries::NTuple{N}; bc=CubicFit(), autocache=true)

Construct an N-dimensional cubic adjoint operator.

Arguments

• grids: N-tuple of grid vectors, one per dimension
• queries: N-tuple of query coordinate vectors (SoA format)
• bc: Boundary condition (single or per-axis tuple)
• autocache: Whether to cache Thomas LU factorizations

Returns

CubicAdjointND operator that can be called as adj(y_bar) or adj(f_bar, y_bar).

Example

x = range(0.0, 1.0, 20)
y = range(0.0, 1.0, 15)
xq = rand(100)
yq = rand(100)

adj = cubic_adjoint((x, y), (xq, yq); bc=CubicFit())
f_bar = adj(y_bar)   # returns 20×15 matrix

AbstractAdjoint{Tg<:AbstractFloat}

Abstract supertype for adjoint (transpose) operators of interpolation.

These operators compute f̄ = Wᵀȳ where W is the forward interpolation weight matrix. They are query-baked, data-free: constructed from grid + query points, then applied to any value vector ȳ regardless of its element type.

Type Parameters

• Tg: Grid/coordinate type (Float32 or Float64)

Note: Only Tg is needed (no Tv) because adjoint operators are value-type independent. The same operator works for Float64, ComplexF64, or any custom value type.

1D Protocol (adjoint_protocol.jl)

Subtypes automatically inherit 6 callable overloads (Vector/Real/Tuple × alloc/in-place), Base.size, Base.Matrix, and exclusive periodic in-place handling by implementing:

_n_queries(adj)::Int                    — number of baked query points (required)
_adjoint_output_length(adj)::Int        — user-facing output length (required)
_adjoint_1d_apply!(f_bar, adj, y_bar, deriv)  — core scatter/solve (required)
_adjoint_internal_length(adj)::Int      — alloc size (default: output length)
_adjoint_1d_has_exclusive_periodic(adj)::Bool  — (default: false)
_adjoint_1d_finalize(f_bar, adj)        — fold+truncate (default: identity)

Subtypes

• AbstractAdjoint1D: 1D adjoint operators
• AbstractAdjointND: N-dimensional adjoint operators

AbstractAdjoint1D{Tg<:AbstractFloat} <: AbstractAdjoint{Tg}

Abstract supertype for 1-dimensional adjoint operators.

Subtypes automatically inherit 6 callable overloads (Vector/Real/Tuple × alloc/in-place), Base.size, Base.Matrix, and exclusive periodic in-place handling from adjoint_protocol.jl by implementing the required interface:

_n_queries(adj)::Int                    — number of baked query points (required)
_adjoint_output_length(adj)::Int        — user-facing output length (required)
_adjoint_1d_apply!(f_bar, adj, y_bar, deriv)  — core scatter/solve (required)
_adjoint_internal_length(adj)::Int      — alloc size (default: output length)
_adjoint_1d_has_exclusive_periodic(adj)::Bool  — (default: false)
_adjoint_1d_finalize(f_bar, adj)        — fold+truncate (default: identity)

Subtypes

• LinearAdjoint{Tg, EP}: Adjoint of linear interpolation (1D, pure scatter)
• CubicAdjoint{Tg, C, BC}: Adjoint of cubic spline interpolation (1D)

CubicAdjoint{Tg, C, BC}

Adjoint (transpose) operator for cubic spline interpolation. Computes f̄ = Wᵀȳ where W is the forward interpolation weight matrix.

Constructed from a grid and query points (query-baked, data-free). The same adjoint can be applied to any ȳ vector regardless of value type.

Type Parameters

• Tg: Grid float type (Float32 or Float64)
• C: CubicSplineCache type (reused from forward interpolation)
• BC: BCPair or PeriodicBC (normalized boundary condition)

Usage

adj = cubic_adjoint(x_grid, x_query; bc=CubicFit())

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

# Derivative adjoint: adjoint of the d-th derivative operator
f_bar = adj(y_bar; deriv=DerivOp(1))   # adjoint of first derivative
f_bar = adj(y_bar; deriv=EvalDeriv2()) # adjoint of second derivative

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

# Dimensions
size(adj)    # (n_grid, n_query)

Mathematical Background

Forward: y = W·f where W = Eᵧ + E_z·A⁻¹·R Adjoint: f̄ = Wᵀȳ = Eᵧᵀȳ + Rᵀ·A⁻ᵀ·E_zᵀȳ

Here A is the tridiagonal moment matrix, R the finite-difference RHS operator, Eᵧ and E_z the evaluation weight matrices for y-values and z-moments respectively.

PolyFit stencil coefficients and periodic q_transpose = A'^{-T}u are computed on the fly at each adj(ȳ) call (O(D) and O(n) respectively, negligible vs overall pipeline).

CubicAdjointND{Tg, N, S, C, MC, BP, MBP}

Adjoint (transpose) operator for N-dimensional cubic Hermite interpolation. Computes f̄ = Wᵀȳ where W is the forward ND cubic interpolation weight matrix.

Constructed from grids and query points (query-baked, data-free). The same adjoint can be applied to any ȳ vector.

Type Parameters

• Tg: Grid float type (Float32 or Float64)
• N: Number of dimensions
• S: Spacing tuple type
• C: Cache tuple type for user BC (per-axis CubicSplineCache)
• MC: Cache tuple type for mixed-partial BC (CubicFit internally)
• BP: User BC pairs — BCPair for non-periodic, PeriodicBC for periodic axes
• MBP: Mixed-partial BC pairs — same convention as BP

Architecture — Dual Cache/BC Structure

The adjoint retains per-axis Thomas LU caches for transpose solves at every adj(ȳ) call. The ND build pipeline processes (d, p_src) pairs, where the BC used depends on whether the partial is pure or mixed:

• Pure derivatives (p_src == 1): use the user's BC → caches + bcs
• Mixed partials (p_src > 1): internally always use CubicFit → mixed_caches + mixed_bcs

For periodic axes, _get_effective_bc(PeriodicBC(), p_src, grid) returns PeriodicBC() for ALL p_src, so caches[d] == mixed_caches[d] (same pool entry) and both bcs entries are PeriodicBC. The periodic build adjoint uses Sherman-Morrison transpose solves instead of standard Thomas.

CubicSplineCache is parameterized on BCPair{L,R} or PeriodicData{Tg}, so C and MC are generally different types (e.g., BCPair{Deriv2,Deriv2} vs BCPair{Deriv1,Deriv1}). When the user's BC is already CubicFit, C == MC and BP == MBP, so Julia optimizes the branch away entirely.

PolyFit stencil coefficients are computed on the fly from BCPair + caches[d].x (O(D) per axis, negligible), matching the forward compute_rhs! pattern. Grids are NOT stored separately — each per-axis CubicSplineCache already owns a mutation-safe copy of its grid via cache.x (inner constructor copy() + typeof(xc) rebinding).

Usage

adj = cubic_adjoint((x, y), (xq, yq); bc=CubicFit())
f_bar = adj(y_bar)              # allocating
adj(f_bar, y_bar)               # in-place (zero-allocation)

# Periodic
adj = cubic_adjoint((x, y), (xq, yq); bc=PeriodicBC())
f_bar = adj(y_bar)              # inclusive periodic

# Mixed: periodic × non-periodic
adj = cubic_adjoint((x, y), (xq, yq); bc=(PeriodicBC(), CubicFit()))

LinearAdjoint{Tg, EP}

Adjoint (transpose) operator for 1D linear interpolation. Computes f̄ = Wᵀȳ where W is the forward linear interpolation weight matrix.

Constructed from grid and query points (query-baked, data-free). The same adjoint can be applied to any ȳ vector.

Type Parameters

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

Fields

• anchors: Pre-computed _LinearAnchoredQuery per query point
• grid_size: Number of grid points (for output allocation)
• extrap: Extrapolation policy (for compile-time OOB dispatch in scatter)

Usage

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

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

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

LinearAdjointND{Tg, N, G, S, EP}

Adjoint (transpose) operator for N-dimensional linear interpolation. Computes f̄ = Wᵀȳ where W is the forward ND linear interpolation weight matrix.

Constructed from grids and query points (query-baked, data-free). The same adjoint can be applied to any ȳ vector.

Type Parameters

• Tg: Grid float type (Float32 or Float64)
• N: Number of dimensions
• G: Grid tuple type
• S: Spacing tuple type
• EP: Extrapolation tuple type

Architecture

Pure scatter operator — no caches, no boundary conditions, no tridiagonal solve. The only per-axis state is pre-baked weights in the anchors.

Usage

adj = linear_adjoint((x, y), (xq, yq))
f_bar = adj(y_bar)              # allocating
adj(f_bar, y_bar)               # in-place (zero-allocation)

linear_adjoint(x, x_query; extrap=NoExtrap()) -> LinearAdjoint

Create a linear interpolation adjoint operator (query-baked, data-free).

Computes f̄ = Wᵀȳ where W is the forward linear interpolation weight matrix. Maps query-space sensitivities back to grid-space sensitivities.

Arguments

• x::AbstractVector: Grid points (must be sorted)
• x_query::AbstractVector: Query points (baked into the operator)
• extrap::AbstractExtrap: Extrapolation mode (default: NoExtrap()). For WrapExtrap(), queries are wrapped to the domain before anchoring. For ClampExtrap()/FillExtrap(), OOB queries use boundary intervals with correct side flags for adjoint scatter skipping.

Example

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

itp = linear_interp(x, f)
adj = linear_adjoint(x, xq)
f_bar = adj(y_bar)

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

Notes

• No boundary conditions needed (linear interpolation has none).
• No caching infrastructure (no tridiagonal solve).
• 2nd+ derivative adjoints always return zero.
• NoExtrap validates all queries are in-domain at construction time.

linear_adjoint(grids::NTuple{N}, queries; extrap=NoExtrap())

Construct an N-dimensional linear adjoint operator.

Arguments

• grids: N-tuple of grid vectors, one per dimension
• queries: Query points (SoA tuple of vectors, AoS, single tuple, SVector, etc.)
• extrap: Extrapolation mode (single or per-axis tuple)

Returns

LinearAdjointND operator that can be called as adj(y_bar) or adj(f_bar, y_bar).

Example

x = range(0.0, 1.0, 20)
y = range(0.0, 1.0, 15)
xq = rand(100)
yq = rand(100)

adj = linear_adjoint((x, y), (xq, yq))
f_bar = adj(y_bar)   # returns 20×15 matrix