Linear Interpolation API

linear_interp(x, y; extrap=:none, search=Binary()) -> LinearInterpolant

Create a callable interpolant for broadcast fusion and reuse.

Arguments

• x::AbstractVector: x-coordinates (must be sorted)
• y::AbstractVector: y-values (can be real or complex)
• extrap::Symbol: :none (default, throws DomainError), :constant, :extension, or :wrap
• search::AbstractSearchPolicy: Default search policy for interval lookup (default: Binary())

Type Handling

• x: Grid coordinates → converted to AbstractFloat, Range structure preserved
• y: Value type determines return type:
  • Real types → promoted to float(eltype(y))
  • Complex{Real} → promoted to Complex{Tg} where Tg = float(real(eltype(y)))
  • Already matching types → no conversion

Returns

LinearInterpolant{Tg, Tv} object where:

• Tg: Grid type (Float32 or Float64)
• Tv: Value type (Tg for real, Complex{Tg} for complex)

Can be:

• Called with scalar: itp(0.5) (uses stored search policy)
• Called with search override: itp(0.5; search=Binary()) (override stored policy)
• Broadcasted: itp.(rho) or @. coef * itp(rho)
• Reused multiple times without re-creating

Examples

# Create with default Binary() search policy
itp = linear_interp(x_data, y_data)

# Create with LinearBinary() as default policy (optimal for sorted queries)
itp = linear_interp(x_data, y_data; search=LinearBinary())

# Scalar call (uses stored policy)
val = itp(0.5)

# Scalar call with search policy override
val = itp(0.5; search=Binary())

# Vector call with hint for ODE/streaming patterns
hint = Ref(1)
for batch in batches
    vals = itp(batch; hint=hint)  # hint persists across calls
end

# Broadcast (creates array)
vals = itp.(query_points)

# Fused broadcast (optimal - no intermediate arrays)
result = @. coefficient * itp(rho) * ne / Te^2

# Wrap to domain (for periodic-like data)
itp_wrap = linear_interp(x_data, y_data; extrap=:wrap)
val_wrap = itp_wrap(2.5)  # wraps to domain

# Compare with 3-argument form (returns array immediately)
vals_direct = linear_interp(x_data, y_data, query_points)

Performance Notes

• Returns lightweight callable (~56 bytes), best for reuse and broadcast fusion
• 3-argument form returns array immediately, best for single use
• Use search=LinearBinary() for sorted query sequences
• Use hint=Ref(idx) for ODE/streaming patterns with persistent hint

linear_interp!(output, x, y, x_targets; extrap=:none, deriv=0, search=Binary())

Zero-allocation linear interpolation with automatic dispatch:

• For AbstractRange x: O(1) direct indexing
• For general AbstractVector x: Search algorithm determined by search parameter

Arguments

• output: Pre-allocated output vector (must be floating-point type)
• extrap::Symbol: :none (default, throws DomainError), :constant, :extension, or :wrap
• deriv::Int: Derivative order (0=value, 1=first derivative, 2=second derivative)
• search::AbstractSearchPolicy: Search algorithm for interval finding
  • Binary() (default): O(log n) binary search, stateless
  • HintedBinary(): O(1) if hint valid, O(log n) fallback
  • LinearBinary(linear_window=8): Linear search within window, then binary fallback

Example

rho = 0.0:0.01:1.0  # Uniform grid → fast O(1) path
y = sin.(rho)
out = Vector{Float64}(undef, 2)
linear_interp!(out, rho, y, [0.55, 0.77])  # throws error if outside domain
linear_interp!(out, rho, y, [-0.1, 1.2]; extrap=:extension)  # linear extrapolation

# Optimized for sorted queries
sorted_queries = sort(rand(1000))
output = zeros(1000)
linear_interp!(output, x_vec, y_vec, sorted_queries; search=LinearBinary(linear_window=8))

Implementation Note

• Optimized core works with AbstractFloat types (calls optimized scalar version)
• Integer/Real inputs automatically promoted via wrapper methods

LinearInterpolant{Tg,Tv,X,Y,P}

Lightweight callable interpolant for broadcast fusion optimization. Returned by linear_interp(x, y) (2-argument form).

Type Parameters

• Tg<:AbstractFloat: Grid type (Float32/Float64) - coordinates and geometry
• Tv: Value type - element type of y (can be Tg, Complex{Tg}, or other Number)
• X<:AbstractVector{Tg}: Grid vector type (preserves Range for O(1) lookup)
• Y<:AbstractVector{Tv}: Values vector type
• P<:AbstractSearchPolicy: Search policy type

Fields

• x::X: x-coordinates (sorted)
• y::Y: y-values
• extrap::Val: Extrapolation mode (Val(:none), Val(:extension), Val(:constant), or Val(:wrap))
• search_policy::P: Default search policy for interval lookup

Usage

# Create interpolator (minimal allocation)
itp = linear_interp(x, y)  # default extrap=:none, search=Binary()

# Create with custom search policy (baked-in default)
itp = linear_interp(x, y; search=LinearBinary())

# Complex-valued interpolation
y_complex = exp.(2im .* x)
itp_c = linear_interp(x, y_complex)  # Works natively with Complex

# Use in broadcast (fused, no intermediate arrays)
result = @. coef * itp(rho) * other_terms

# Reuse interpolator multiple times
vals1 = itp.(query_points1)
vals2 = @. compute(itp(query_points2))

# Extrapolation options
itp_ext = linear_interp(x, y; extrap=:extension)  # linear extrap
itp_const = linear_interp(x, y; extrap=:constant)  # clamp to boundary values
itp_wrap = linear_interp(x, y; extrap=:wrap)  # wrap to domain
val = itp_wrap(2.5)  # wraps to domain

# Override search policy at call time
itp(0.5; search=Binary())  # override stored policy

LinearInterpolantND{Tg, Tv, N, G, S, E, P}

N-dimensional multilinear interpolant for tensor-product linear interpolation.

Performs bilinear interpolation in 2D, trilinear in 3D, and n-linear in N dimensions. The interpolation is exact at grid points and linearly blended between them.

Type Parameters

• Tg<:AbstractFloat: Grid/coordinate type (Float32 or Float64)
• Tv: Value type (can be Tg, Complex{Tg}, or other Number)
• N: Number of dimensions
• G<:Tuple{Vararg{AbstractVector,N}}: Grid tuple type (supports heterogeneous grids)
• S<:Tuple{Vararg{AbstractGridSpacing,N}}: Spacing tuple type
• E<:Tuple{Vararg{ExtrapVal,N}}: Extrapolation mode tuple type
• P<:Tuple{Vararg{AbstractSearchPolicy,N}}: Search policy tuple type

Fields

• grids: N-tuple of grid vectors for each dimension
• spacings: N-tuple of grid spacing info (for O(1) h lookup)
• data: N-dimensional array of values at grid points
• extraps: N-tuple of extrapolation modes
• searches: N-tuple of search policies

Performance

• Construction: O(1) - just stores references, no precomputation
• Query: O(2^N) - evaluates 2^N corner weights and sums
• Memory: Only stores reference to data array

Thread-Safety

Immutable after construction; safe for concurrent read access.

Mathematical Formulation

For a query point (x₁, x₂, ..., xₙ), compute:

1. Find cell indices (i₁, i₂, ..., iₙ) containing the point

2. Compute normalized coordinates αₖ = (xₖ - grid[iₖ]) / hₖ

3. Sum over all 2^N corners:

   f(x) = Σ_{b∈{0,1}^N} data[i₁+b₁, i₂+b₂, ..., iₙ+bₙ] × ∏ₖ wₖ(bₖ)

   where wₖ(0) = 1-αₖ and wₖ(1) = αₖ

Derivatives

• ∂f/∂xₖ: Replace wₖ with derivative weights: w'ₖ(0) = -1/hₖ, w'ₖ(1) = 1/hₖ
• ∂²f/∂xₖ² = 0 (linear interpolation has zero second derivative within cells)

Example

x = range(0.0, 1.0, 11)
y = range(0.0, 2.0, 21)
data = [sin(xi) * cos(yj) for xi in x, yj in y]

itp = linear_interp((x, y), data)  # Returns LinearInterpolantND{..., 2, ...}
itp((0.5, 1.0))                     # Bilinear interpolation at (0.5, 1.0)
itp((0.5, 1.0); deriv=Val((1,0)))   # ∂f/∂x at (0.5, 1.0)