From 108b077d705a02e584369f5fad1645e274533617 Mon Sep 17 00:00:00 2001 From: Min-Gu Yoo Date: Mon, 23 Feb 2026 16:15:04 -0800 Subject: [PATCH 01/16] perf!: optimize _search_linear_binary! and change default window to 2 Optimize the LinearBinary search hot path: - Remove redundant hint clamp (internal hints are always valid) - Skip hint write on direct hit (ix unchanged) - Use single comparison per linear step (direction already bounds one side) - Change default linear_window from 8 to 2 (sweet spot for minimal random overhead while retaining sorted/clustered locality gains) Benchmark results (Grid=2000, vector call): - Random queries: ~2x slower than Binary (structural; branch predictor) - Sorted queries: ~5x faster than Binary - Clustered: ~4.7x faster | DenseLocal: ~5.9x faster --- src/core/search.jl | 64 +++++++++++++++++++++++++--------------------- 1 file changed, 35 insertions(+), 29 deletions(-) diff --git a/src/core/search.jl b/src/core/search.jl index ccbaa698..f0835aff 100644 --- a/src/core/search.jl +++ b/src/core/search.jl @@ -40,8 +40,8 @@ pure binary search is used. # Example ```julia -val = linear_interp(x, y, 0.5; search=Binary()) # pure binary search -val = linear_interp(x, y, 0.5) # same (default) +val = linear_interp(x, y, 0.5; search=Binary()) # explicit binary search +val = linear_interp(x, y, 0.5) # default: LinearBinary() # With hint: auto-upgrades to HintedBinary behavior hint = Ref(1) @@ -126,7 +126,7 @@ fall in adjacent or nearby intervals. # Construction Use the factory function (recommended) to construct with a curated set of values: ```julia -LinearBinary() # default MAX=8 +LinearBinary() # default MAX=2 LinearBinary(linear_window=4) # custom MAX=4 ``` @@ -147,7 +147,7 @@ struct LinearBinary{MAX} <: AbstractSearchPolicy end """ LinearBinary(linear_window::Integer) - LinearBinary(; linear_window::Integer=8) + LinearBinary(; linear_window::Integer=2) Factory constructor for `LinearBinary{MAX}` with a **curated set of `linear_window` values**. @@ -160,7 +160,7 @@ By restricting to powers of 2, we limit specialization to just 7 variants while covering the practical range of use cases. # Arguments -- `linear_window::Integer=8`: Size of the linear search window before binary fallback. +- `linear_window::Integer=2`: Size of the linear search window before binary fallback. **Allowed values**: `1, 2, 4, 8, 16, 32, 64, 128` (powers of 2 from 2⁰ to 2⁷) # Throws @@ -168,15 +168,15 @@ covering the practical range of use cases. # Example ```julia -policy = LinearBinary() # LinearBinary{8}() (default) +policy = LinearBinary() # LinearBinary{2}() (default) policy = LinearBinary(linear_window=4) # LinearBinary{4}() -policy = LinearBinary(linear_window=16) # LinearBinary{16}() +policy = LinearBinary(linear_window=8) # LinearBinary{8}() policy = LinearBinary(linear_window=3) # ERROR: ArgumentError ``` # Choosing `linear_window` -- **Small values (1–4)**: Lower overhead, but more frequent binary fallbacks -- **Medium values (8–16)**: Good balance for typical sorted query patterns +- **Small values (1–2)**: Minimal overhead, best for default usage and mixed query patterns +- **Medium values (4–16)**: Good balance for known-sorted query patterns - **Large values (32–128)**: For highly localized queries or very large datasets """ function LinearBinary(linear_window::Integer) @@ -190,7 +190,7 @@ function LinearBinary(linear_window::Integer) linear_window == 128 && return LinearBinary{128}() throw(ArgumentError("`linear_window` must be one of (1, 2, 4, 8, 16, 32, 64, 128), got $linear_window")) end -LinearBinary(; linear_window::Integer=8) = LinearBinary(linear_window) +LinearBinary(; linear_window::Integer=2) = LinearBinary(linear_window) # ---------------------------------------- # Hint Types (Internal) @@ -532,6 +532,11 @@ end Bounded linear search within MAX-sized window, then binary fallback. Optimal for monotonic query sequences. + +# Optimizations over naive implementation: +- No hint clamp: internal hints are always valid (initialized to 1, updated to valid idx) +- No hint write on direct hit: `ix` is unchanged, skip redundant `hint_ref[] = ix` +- Single comparison per linear step: direction already determines one bound """ @inline function _search_linear_binary!( x::AbstractVector{T}, @@ -542,31 +547,32 @@ Optimal for monotonic query sequences. ix = hint_ref[] n = length(x) @inbounds begin - ix = clamp(ix, 1, n - 1) - if x[ix] <= xq < x[ix + 1] - hint_ref[] = ix - return ix, x[ix], x[ix + 1] - end - if xq < x[ix] - for _ in 1:MAX + # Direct hit — most common for sorted/monotonic queries + xL = x[ix] + xR = x[ix + 1] + xL <= xq < xR && return ix, xL, xR # no hint write (ix unchanged) + + if xq < xL + # Walk left: xq < x[ix] guaranteed ⟹ after ix-=1, + # x[ix+1] = old x[ix] > xq — right bound already satisfied. + # Only need: x[ix] <= xq (single comparison per step) + lo = max(1, ix - MAX) + while ix > lo ix -= 1 - ix < 1 && break - if x[ix] <= xq < x[ix + 1] - hint_ref[] = ix - return ix, x[ix], x[ix + 1] - end + x[ix] <= xq && (hint_ref[] = ix; return ix, x[ix], x[ix + 1]) end - else - for _ in 1:MAX + else # xq >= xR + # Walk right: xq >= x[ix+1] guaranteed ⟹ after ix+=1, + # x[ix] = old x[ix+1] <= xq — left bound already satisfied. + # Only need: xq < x[ix+1] (single comparison per step) + hi = min(n - 1, ix + MAX) + while ix < hi ix += 1 - ix > n - 1 && break - if x[ix] <= xq < x[ix + 1] - hint_ref[] = ix - return ix, x[ix], x[ix + 1] - end + xq < x[ix + 1] && (hint_ref[] = ix; return ix, x[ix], x[ix + 1]) end end end + # Binary fallback — full range (narrowing saves < 1 iteration, not worth extra branches) idx, xL, xR = _search_binary(x, xq) hint_ref[] = idx return idx, xL, xR From be16a2d6107dacd4da9ab4e12e7f4e686d4ba635 Mon Sep 17 00:00:00 2001 From: Min-Gu Yoo Date: Mon, 23 Feb 2026 16:15:21 -0800 Subject: [PATCH 02/16] refactor!: change default search policy from Binary() to LinearBinary() MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Update all interpolation constructors and docstrings across linear, cubic, quadratic, and constant modules (1D and ND) to use LinearBinary() as the default search policy instead of Binary(). This is a breaking change for code that relied on the default being Binary() — most users won't notice since LinearBinary automatically falls back to binary search when the linear walk misses. Affected: 25 source files across all interpolation types. --- src/constant/constant_interpolant.jl | 6 ++-- src/constant/constant_oneshot.jl | 18 +++++----- src/constant/constant_series_interp.jl | 14 ++++---- src/constant/constant_types.jl | 2 +- src/constant/nd/constant_nd_interpolant.jl | 4 +-- src/constant/nd/constant_nd_oneshot.jl | 10 +++--- src/cubic/cubic_eval.jl | 2 +- src/cubic/cubic_interpolant.jl | 18 +++++----- src/cubic/cubic_oneshot.jl | 38 ++++++++++---------- src/cubic/cubic_series_interp.jl | 16 ++++----- src/cubic/cubic_types.jl | 2 +- src/cubic/nd/cubic_nd_interpolant.jl | 4 +-- src/cubic/nd/cubic_nd_oneshot.jl | 10 +++--- src/linear/linear_interpolant.jl | 8 ++--- src/linear/linear_oneshot.jl | 18 +++++----- src/linear/linear_series_interp.jl | 14 ++++---- src/linear/linear_types.jl | 4 +-- src/linear/nd/linear_nd_interpolant.jl | 4 +-- src/linear/nd/linear_nd_oneshot.jl | 10 +++--- src/quadratic/nd/quadratic_nd_interpolant.jl | 4 +-- src/quadratic/nd/quadratic_nd_oneshot.jl | 10 +++--- src/quadratic/quadratic_interpolant.jl | 6 ++-- src/quadratic/quadratic_oneshot.jl | 18 +++++----- src/quadratic/quadratic_series_interp.jl | 14 ++++---- src/quadratic/quadratic_types.jl | 2 +- 25 files changed, 128 insertions(+), 128 deletions(-) diff --git a/src/constant/constant_interpolant.jl b/src/constant/constant_interpolant.jl index 686d7332..c0ec865b 100644 --- a/src/constant/constant_interpolant.jl +++ b/src/constant/constant_interpolant.jl @@ -51,7 +51,7 @@ end # ======================================== """ - constant_interp(x, y; extrap=NoExtrap(), side=NearestSide(), search=Binary()) -> ConstantInterpolant + constant_interp(x, y; extrap=NoExtrap(), side=NearestSide(), search=LinearBinary()) -> ConstantInterpolant Create a callable interpolant for broadcast fusion and reuse. @@ -60,7 +60,7 @@ Create a callable interpolant for broadcast fusion and reuse. - `y::AbstractVector`: y-values (can be Real or Complex) - `extrap::AbstractExtrap`: `NoExtrap()` (default), `ConstExtrap()`, `ExtendExtrap()`, or `WrapExtrap()` - `side::AbstractSide`: Side selection (NearestSide(), LeftSide(), RightSide()) -- `search::AbstractSearchPolicy`: Default search policy (default: `Binary()`) +- `search::AbstractSearchPolicy`: Default search policy (default: `LinearBinary()`) # Returns `ConstantInterpolant{Tg, Tv}` object for scalar/broadcast evaluation. @@ -109,7 +109,7 @@ end y::AbstractVector{TY}; extrap::AbstractExtrap=NoExtrap(), side::AbstractSide=NearestSide(), - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {TX<:Real, TY} x_p, y_p = _promote_itp_inputs(x, y) return ConstantInterpolant(x_p, y_p; extrap, side, search) diff --git a/src/constant/constant_oneshot.jl b/src/constant/constant_oneshot.jl index 22cb841b..2cecc9fb 100644 --- a/src/constant/constant_oneshot.jl +++ b/src/constant/constant_oneshot.jl @@ -149,7 +149,7 @@ end # ======================================== """ - constant_interp(x, y, xi; extrap=NoExtrap(), side=NearestSide(), deriv=EvalValue(), search=Binary()) + constant_interp(x, y, xi; extrap=NoExtrap(), side=NearestSide(), deriv=EvalValue(), search=LinearBinary()) Constant (step/piecewise constant) interpolation at a single point. @@ -200,7 +200,7 @@ vals = constant_interp(x, y, sorted_queries; search=LinearBinary(linear_window=8 extrap::AbstractExtrap=NoExtrap(), side::AbstractSide=NearestSide(), deriv::DerivOp=EvalValue(), - search=Binary(), + search=LinearBinary(), hint::Union{Nothing,Base.RefValue{Int}}=nothing ) where {Tg<:AbstractFloat, Tv, Tq<:Real} @boundscheck length(y) == length(x) || throw(ArgumentError("x and y must have same length")) @@ -214,7 +214,7 @@ end # ======================================== """ - constant_interp!(output, x, y, x_targets; extrap=NoExtrap(), side=NearestSide(), deriv=EvalValue(), search=Binary()) + constant_interp!(output, x, y, x_targets; extrap=NoExtrap(), side=NearestSide(), deriv=EvalValue(), search=LinearBinary()) Zero-allocation constant interpolation for multiple query points. @@ -246,7 +246,7 @@ function constant_interp!( extrap::AbstractExtrap=NoExtrap(), side::AbstractSide=NearestSide(), deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:AbstractFloat, Tv} @assert length(y) == length(x) "x and y must have same length" @assert length(output) == length(x_targets) "output must match x_targets length" @@ -264,7 +264,7 @@ end # ======================================== """ - constant_interp(x, y, x_targets; extrap=NoExtrap(), side=NearestSide(), deriv=EvalValue(), search=Binary()) + constant_interp(x, y, x_targets; extrap=NoExtrap(), side=NearestSide(), deriv=EvalValue(), search=LinearBinary()) Constant interpolation for multiple query points (allocating version). @@ -287,7 +287,7 @@ function constant_interp( extrap::AbstractExtrap=NoExtrap(), side::AbstractSide=NearestSide(), deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:AbstractFloat, Tv} output = Vector{Tv}(undef, length(x_targets)) constant_interp!(output, x, y, x_targets; extrap, side, deriv, search) @@ -314,7 +314,7 @@ end extrap::AbstractExtrap=NoExtrap(), side::AbstractSide=NearestSide(), deriv::DerivOp=EvalValue(), - search=Binary(), + search=LinearBinary(), hint::Union{Nothing,Base.RefValue{Int}}=nothing ) where {Tg<:Real, Tv, Tq<:Real} x_typed, y_typed = _promote_itp_inputs(x, y) @@ -334,7 +334,7 @@ function constant_interp( extrap::AbstractExtrap=NoExtrap(), side::AbstractSide=NearestSide(), deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:Real, Tv, Tq<:Real} x_typed, y_typed, xq_typed = _promote_itp_inputs(x, y, x_targets) Tv_float = eltype(y_typed) @@ -355,7 +355,7 @@ function constant_interp!( extrap::AbstractExtrap=NoExtrap(), side::AbstractSide=NearestSide(), deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:Real, Tv, Tq<:Real} @assert length(y) == length(x) "x and y must have same length" @assert length(output) == length(x_targets) "output must match x_targets length" diff --git a/src/constant/constant_series_interp.jl b/src/constant/constant_series_interp.jl index c3afad5e..dfc4ce8a 100644 --- a/src/constant/constant_series_interp.jl +++ b/src/constant/constant_series_interp.jl @@ -77,7 +77,7 @@ mutable struct ConstantSeriesInterpolant{Tg<:AbstractFloat, Tv, E<:AbstractExtra y::Matrix{Tv}, extrap::E, side::SD, - search::P=Binary() + search::P=LinearBinary() ) where {Tg<:AbstractFloat, Tv, E<:AbstractExtrap, SD<:AbstractSide, P<:AbstractSearchPolicy, X<:AbstractVector{Tg}} new{Tg,Tv,E,SD,P,X}(x, y, LazyTranspose{Tv}(), extrap, side, search) end @@ -333,7 +333,7 @@ function constant_interp( ys::AbstractVector{<:AbstractVector{Tv}}; side::AbstractSide=NearestSide(), extrap::AbstractExtrap=NoExtrap(), - search::P=Binary() + search::P=LinearBinary() ) where {Tg<:AbstractFloat, Tv, P<:AbstractSearchPolicy} # Check if Tv's float base requires grid widening (not for Int types) # Int-based types (Complex{Int}) are handled by internal _value_type conversion @@ -394,7 +394,7 @@ function constant_interp( Y::AbstractMatrix{Tv}; side::AbstractSide=NearestSide(), extrap::AbstractExtrap=NoExtrap(), - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:AbstractFloat, Tv} # Check if Tv's float base requires grid widening Tv_real = _real_eltype(Tv) @@ -431,7 +431,7 @@ function constant_interp( ys::AbstractVector{<:AbstractVector{Tv}}; side::AbstractSide=NearestSide(), extrap::AbstractExtrap=NoExtrap(), - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:Real, Tv} # Compute promoted grid type (Tg may be Int, promotes to Float) Tg_float = float(promote_type(Tg, _real_eltype(Tv))) @@ -445,7 +445,7 @@ function constant_interp( Y::AbstractMatrix{Tv}; side::AbstractSide=NearestSide(), extrap::AbstractExtrap=NoExtrap(), - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:Real, Tv} Tg_float = float(promote_type(Tg, _real_eltype(Tv))) x_typed = _to_float(x, Tg_float) @@ -457,7 +457,7 @@ end # ======================================== """ - (sitp::ConstantSeriesInterpolant)(xq::Real; deriv=EvalValue(), search=Binary()) + (sitp::ConstantSeriesInterpolant)(xq::Real; deriv=EvalValue(), search=LinearBinary()) Evaluate multi-Y interpolant at scalar query point (out-of-place). @@ -483,7 +483,7 @@ function (sitp::ConstantSeriesInterpolant{Tg,Tv,P})( end """ - (sitp::ConstantSeriesInterpolant)(output::AbstractVector, xq::Real; deriv=EvalValue(), search=Binary()) + (sitp::ConstantSeriesInterpolant)(output::AbstractVector, xq::Real; deriv=EvalValue(), search=LinearBinary()) Evaluate multi-Y interpolant at scalar query point (in-place). """ diff --git a/src/constant/constant_types.jl b/src/constant/constant_types.jl index ded6a338..10e87394 100644 --- a/src/constant/constant_types.jl +++ b/src/constant/constant_types.jl @@ -78,7 +78,7 @@ end y::Y; extrap::AbstractExtrap=NoExtrap(), side::AbstractSide=NearestSide(), - search::P=Binary() + search::P=LinearBinary() ) where {Tg<:AbstractFloat, Tv, X<:AbstractVector{Tg}, Y<:AbstractVector{Tv}, P<:AbstractSearchPolicy} E = typeof(extrap) SD = typeof(side) diff --git a/src/constant/nd/constant_nd_interpolant.jl b/src/constant/nd/constant_nd_interpolant.jl index e358aa41..ec0d9c80 100644 --- a/src/constant/nd/constant_nd_interpolant.jl +++ b/src/constant/nd/constant_nd_interpolant.jl @@ -21,7 +21,7 @@ Create an N-dimensional constant interpolant with tuple-grid API. # Keyword Arguments - `side=NearestSide()`: Side selection mode (`NearestSide()`, `LeftSide()`, `RightSide()`) or per-axis tuple - `extrap=NoExtrap()`: Extrapolation mode (`NoExtrap()`, `ConstExtrap()`, `ExtendExtrap()`, `WrapExtrap()`) or per-axis tuple -- `search=Binary()`: Search policy or per-axis tuple +- `search=LinearBinary()`: Search policy or per-axis tuple # Returns - `ConstantInterpolantND{Tg,Tv,N}`: Callable interpolant object @@ -51,7 +51,7 @@ function constant_interp( data::AbstractArray{Tv_raw, N}; side::Union{AbstractSide, Tuple{Vararg{AbstractSide}}} = NearestSide(), extrap::Union{AbstractExtrap, NTuple{N, AbstractExtrap}} = NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = Binary() + search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = LinearBinary() ) where {N, Tv_raw} # Validate grid dimensions _validate_nd_grids(grids, data) diff --git a/src/constant/nd/constant_nd_oneshot.jl b/src/constant/nd/constant_nd_oneshot.jl index cf80b2c2..044ca587 100644 --- a/src/constant/nd/constant_nd_oneshot.jl +++ b/src/constant/nd/constant_nd_oneshot.jl @@ -150,7 +150,7 @@ function constant_interp( query::Tuple{Vararg{Real, N}}; side::Union{AbstractSide, Tuple{Vararg{AbstractSide}}} = NearestSide(), extrap::Union{AbstractExtrap, NTuple{N, AbstractExtrap}} = NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = Binary(), + search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = LinearBinary(), deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue() ) where {Tv, N} # Any derivative of constant interpolation is zero @@ -183,7 +183,7 @@ function constant_interp( queries::NTuple{N, AbstractVector{<:Real}}; side::Union{AbstractSide, Tuple{Vararg{AbstractSide}}} = NearestSide(), extrap::Union{AbstractExtrap, NTuple{N, AbstractExtrap}} = NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = Binary(), + search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = LinearBinary(), deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue() ) where {Tv, N} if _is_any_deriv(deriv) @@ -216,7 +216,7 @@ function constant_interp( queries::AbstractVector{<:Tuple{Vararg{Real, N}}}; side::Union{AbstractSide, Tuple{Vararg{AbstractSide}}} = NearestSide(), extrap::Union{AbstractExtrap, NTuple{N, AbstractExtrap}} = NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = Binary(), + search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = LinearBinary(), deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue() ) where {Tv, N} if _is_any_deriv(deriv) @@ -254,7 +254,7 @@ function constant_interp!( queries::NTuple{N, AbstractVector{<:Real}}; side::Union{AbstractSide, Tuple{Vararg{AbstractSide}}} = NearestSide(), extrap::Union{AbstractExtrap, NTuple{N, AbstractExtrap}} = NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = Binary(), + search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = LinearBinary(), deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue() ) where {Tv, N} if _is_any_deriv(deriv) @@ -288,7 +288,7 @@ function constant_interp!( queries::AbstractVector{<:Tuple{Vararg{Real, N}}}; side::Union{AbstractSide, Tuple{Vararg{AbstractSide}}} = NearestSide(), extrap::Union{AbstractExtrap, NTuple{N, AbstractExtrap}} = NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = Binary(), + search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = LinearBinary(), deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue() ) where {Tv, N} if _is_any_deriv(deriv) diff --git a/src/cubic/cubic_eval.jl b/src/cubic/cubic_eval.jl index ef8b737b..c6a63f13 100644 --- a/src/cubic/cubic_eval.jl +++ b/src/cubic/cubic_eval.jl @@ -273,7 +273,7 @@ Uses task-local pool for workspace allocation. x_query::Tg; extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), - search=Binary(), + search=LinearBinary(), hint::Union{Nothing,Base.RefValue{Int}}=nothing ) where {Tg<:AbstractFloat, Tv, X, F, BC, S<:AbstractGridSpacing{Tg}} @assert length(y) == length(cache.x) "y length must match cache grid" diff --git a/src/cubic/cubic_interpolant.jl b/src/cubic/cubic_interpolant.jl index b01b4b73..ecd7b79d 100644 --- a/src/cubic/cubic_interpolant.jl +++ b/src/cubic/cubic_interpolant.jl @@ -293,7 +293,7 @@ so the pool memory can be safely reused after this function returns. bc_pair::BCPair{L,R}, extrap::AbstractExtrap, autocache::Bool, - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:AbstractFloat, Tv, L<:PointBC, R<:PointBC} # Cache uses structural equivalent (PolyFit → Deriv1 via _cache_bc_pair internally) cache = _get_cubic_cache(x, bc_pair, autocache) @@ -319,7 +319,7 @@ so the pool memory can be safely reused after this function returns. y::AbstractVector{Tv}, bc::PeriodicBC, autocache::Bool, - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:AbstractFloat, Tv} x, y = _prepare_periodic(x, y, bc) _check_periodic_endpoints(y) @@ -351,7 +351,7 @@ Handles conversion of Real BC values to Complex when needed. # ======================================== """ - cubic_interp(x, y; bc=CubicFit(), extrap=NoExtrap(), autocache=true, search=Binary()) -> CubicInterpolant + cubic_interp(x, y; bc=CubicFit(), extrap=NoExtrap(), autocache=true, search=LinearBinary()) -> CubicInterpolant Create a callable interpolant for broadcast fusion and reuse. @@ -364,7 +364,7 @@ enabling true zero-allocation scalar evaluations in broadcast operations. - `bc::AbstractBC`: Boundary condition (default: `CubicFit()`) - `extrap::AbstractExtrap`: `NoExtrap()` (default), `ConstExtrap()`, `ExtendExtrap()`, or `WrapExtrap()` - `autocache::Bool`: Enable automatic caching (default: `true`) -- `search::AbstractSearchPolicy`: Default search policy (default: `Binary()`) +- `search::AbstractSearchPolicy`: Default search policy (default: `LinearBinary()`) # Example ```julia @@ -393,7 +393,7 @@ val = itp(0.5) # returns ComplexF64 bc::AbstractBC, extrap::AbstractExtrap, autocache::Bool, - search::P=Binary() + search::P=LinearBinary() ) where {Tg<:AbstractFloat, Tv, P<:AbstractSearchPolicy} if _is_periodic_bc(bc) return _build_interpolant_periodic(x, y, bc, autocache, search) @@ -410,7 +410,7 @@ function cubic_interp( bc::AbstractBC=CubicFit(), extrap::AbstractExtrap=NoExtrap(), autocache::Bool=true, - search::P=Binary() + search::P=LinearBinary() ) where {Tg<:AbstractFloat, Tv, P<:AbstractSearchPolicy} # Auto-promote x/y types (zero allocation if already compatible) x_p, y_p = _promote_itp_inputs(x, y) @@ -419,7 +419,7 @@ function cubic_interp( end """ - cubic_interp(cache, y; extrap=NoExtrap(), search=Binary()) -> CubicInterpolant + cubic_interp(cache, y; extrap=NoExtrap(), search=LinearBinary()) -> CubicInterpolant Create a callable interpolant from a pre-built cache. @@ -441,7 +441,7 @@ so the pool memory can be safely reused after this function returns. cache::CubicSplineCache{Tg}, y::AbstractVector{Tv}; extrap::AbstractExtrap=NoExtrap(), - search::P=Binary() + search::P=LinearBinary() ) where {Tg<:AbstractFloat, Tv, P<:AbstractSearchPolicy} tmp_z = similar!(pool, y) _solve_system!(tmp_z, cache, y, cache.bc_config) @@ -462,7 +462,7 @@ function cubic_interp( bc::AbstractBC=CubicFit(), extrap::AbstractExtrap=NoExtrap(), autocache::Bool=true, - search::P=Binary() + search::P=LinearBinary() ) where {TX<:Real, TY, P<:AbstractSearchPolicy} x_p, y_p = _promote_itp_inputs(x, y) bc_promoted = _promote_bc(bc, eltype(y_p)) diff --git a/src/cubic/cubic_oneshot.jl b/src/cubic/cubic_oneshot.jl index c5019b62..cf35b427 100644 --- a/src/cubic/cubic_oneshot.jl +++ b/src/cubic/cubic_oneshot.jl @@ -12,7 +12,7 @@ # ======================================== """ - cubic_interp!(output, cache, y, x_query; extrap=NoExtrap(), deriv=EvalValue(), search=Binary()) + cubic_interp!(output, cache, y, x_query; extrap=NoExtrap(), deriv=EvalValue(), search=LinearBinary()) In-place cubic spline interpolation using cached LU factorization. @@ -26,7 +26,7 @@ Thread-safe: workspaces allocated from task-local pool. - `x_query::AbstractVector{T}`: Query points - `extrap::AbstractExtrap`: `NoExtrap()` (default), `ConstExtrap()`, `ExtendExtrap()`, or `WrapExtrap()` - `deriv::DerivOp=EvalValue()`: Derivative order (0=value, 1=first derivative, 2=second derivative) -- `search::AbstractSearchPolicy=Binary()`: Search algorithm for interval finding +- `search::AbstractSearchPolicy=LinearBinary()`: Search algorithm for interval finding """ @inline @with_pool pool function cubic_interp!( output::AbstractVector{Tv}, @@ -35,7 +35,7 @@ Thread-safe: workspaces allocated from task-local pool. x_query::AbstractVector{Tg}; extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:AbstractFloat, Tv, X, F, BC} @assert length(y) == length(cache.x) "y length must match cache grid" @assert length(output) == length(x_query) "output length must match x_query" @@ -216,7 +216,7 @@ Pool-based exclusive extension: zero-alloc after warmup. end """ - cubic_interp!(output, x, y, x_query; bc=CubicFit(), extrap=NoExtrap(), autocache=true, deriv=EvalValue(), search=Binary()) + cubic_interp!(output, x, y, x_query; bc=CubicFit(), extrap=NoExtrap(), autocache=true, deriv=EvalValue(), search=LinearBinary()) In-place cubic spline interpolation with optional automatic caching. """ @@ -229,7 +229,7 @@ In-place cubic spline interpolation with optional automatic caching. extrap::AbstractExtrap=NoExtrap(), autocache::Bool=true, deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:AbstractFloat, Tv} searcher = _to_searcher(search) # Periodic BC @@ -251,7 +251,7 @@ end x_query::Tg; extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:AbstractFloat, Tv, X, F, BC} @assert length(output) >= 1 "output must have at least 1 element" output[1] = cubic_interp_scalar(cache, y, x_query; extrap=extrap, deriv=deriv, search=search) @@ -267,7 +267,7 @@ end extrap::AbstractExtrap=NoExtrap(), autocache::Bool=true, deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:AbstractFloat, Tv} @assert length(output) >= 1 "output must have at least 1 element" output[1] = cubic_interp(x, y, x_query; bc, extrap, autocache, deriv, search) @@ -279,13 +279,13 @@ end # ======================================== """ - cubic_interp(cache, y, x_query; extrap=NoExtrap(), deriv=EvalValue(), search=Binary()) -> Vector{T} + cubic_interp(cache, y, x_query; extrap=NoExtrap(), deriv=EvalValue(), search=LinearBinary()) -> Vector{T} Allocating version of cubic spline interpolation using cached LU factorization. # Arguments - `deriv::DerivOp=EvalValue()`: Derivative order (0=value, 1=first derivative, 2=second derivative) -- `search::AbstractSearchPolicy=Binary()`: Search algorithm for interval finding +- `search::AbstractSearchPolicy=LinearBinary()`: Search algorithm for interval finding # Example ```julia @@ -305,7 +305,7 @@ function cubic_interp( x_query::AbstractVector{Tg}; extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:AbstractFloat, Tv} output = Vector{Tv}(undef, length(x_query)) cubic_interp!(output, cache, y, x_query; extrap=extrap, deriv=deriv, search=search) @@ -313,13 +313,13 @@ function cubic_interp( end """ - cubic_interp(x, y, x_query; bc=CubicFit(), extrap=NoExtrap(), autocache=true, deriv=EvalValue(), search=Binary()) -> Vector{T} + cubic_interp(x, y, x_query; bc=CubicFit(), extrap=NoExtrap(), autocache=true, deriv=EvalValue(), search=LinearBinary()) -> Vector{T} Cubic spline interpolation with optional automatic caching. # Arguments - `deriv::DerivOp=EvalValue()`: Derivative order (0=value, 1=first derivative, 2=second derivative) -- `search::AbstractSearchPolicy=Binary()`: Search algorithm for interval finding +- `search::AbstractSearchPolicy=LinearBinary()`: Search algorithm for interval finding # Extrapolation Modes - `NoExtrap()` (default): Throws DomainError if query point is outside domain @@ -347,7 +347,7 @@ function cubic_interp( extrap::AbstractExtrap=NoExtrap(), autocache::Bool=true, deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:AbstractFloat, Tv} output = Vector{Tv}(undef, length(x_query)) cubic_interp!(output, x, y, x_query; bc, extrap, autocache, deriv, search) @@ -356,7 +356,7 @@ end # Scalar query - zero allocation cubic_interp(cache::CubicSplineCache{Tg}, y::AbstractVector{Tv}, - x_query::Tg; extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), search=Binary(), hint::Union{Nothing,Base.RefValue{Int}}=nothing) where {Tg<:AbstractFloat, Tv} = + x_query::Tg; extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), search=LinearBinary(), hint::Union{Nothing,Base.RefValue{Int}}=nothing) where {Tg<:AbstractFloat, Tv} = cubic_interp_scalar(cache, y, x_query; extrap=extrap, deriv=deriv, search=search, hint=hint) # Primary scalar method - AD-compatible @@ -369,7 +369,7 @@ function cubic_interp( extrap::AbstractExtrap=NoExtrap(), autocache::Bool=true, deriv::DerivOp=EvalValue(), - search=Binary(), + search=LinearBinary(), hint::Union{Nothing,Base.RefValue{Int}}=nothing ) where {Tg<:AbstractFloat, Tv, Tq<:Real} searcher = _to_searcher(search, hint) @@ -398,7 +398,7 @@ function cubic_interp( extrap::AbstractExtrap=NoExtrap(), autocache::Bool=true, deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:Real, Tv, Tq<:Real} x_typed, y_typed, xq_typed = _promote_itp_inputs(x, y, x_query) return cubic_interp(x_typed, y_typed, xq_typed; bc, extrap, autocache, deriv, search) @@ -414,7 +414,7 @@ function cubic_interp( extrap::AbstractExtrap=NoExtrap(), autocache::Bool=true, deriv::DerivOp=EvalValue(), - search=Binary(), + search=LinearBinary(), hint::Union{Nothing,Base.RefValue{Int}}=nothing ) where {Tg<:Real, Tv, Tq<:Real} x_typed, y_typed = _promote_itp_inputs(x, y) @@ -432,7 +432,7 @@ function cubic_interp!( extrap::AbstractExtrap=NoExtrap(), autocache::Bool=true, deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:Real, Tv, Tq<:Real} @assert length(y) == length(x) "x and y must have same length" @assert length(output) == length(x_query) "output must match x_query length" @@ -463,7 +463,7 @@ function cubic_interp!( extrap::AbstractExtrap=NoExtrap(), autocache::Bool=true, deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:Real, Tv, Tq<:Real} @assert length(output) >= 1 "output must have at least 1 element" diff --git a/src/cubic/cubic_series_interp.jl b/src/cubic/cubic_series_interp.jl index 49245259..be5d2749 100644 --- a/src/cubic/cubic_series_interp.jl +++ b/src/cubic/cubic_series_interp.jl @@ -99,7 +99,7 @@ mutable struct CubicSeriesInterpolant{ y::Matrix{Tv}, z::Matrix{Tv}, extrap::E, - search::P=Binary() + search::P=LinearBinary() ) where {Tg<:AbstractFloat, Tv, C<:CubicSplineCache{Tg}, B, E<:AbstractExtrap, P<:AbstractSearchPolicy} new{Tg, Tv, C, B, E, P}( cache, bc_for_solve, y, z, @@ -574,7 +574,7 @@ function cubic_interp( extrap::AbstractExtrap=NoExtrap(), autocache::Bool=true, precompute_transpose::Bool=false, - search::P=Binary() + search::P=LinearBinary() ) where {Tg<:AbstractFloat, Tv, P<:AbstractSearchPolicy} # Validate input @assert !isempty(ys) "ys must not be empty" @@ -642,7 +642,7 @@ function _build_series_periodic( n_series_count::Int, autocache::Bool, precompute_transpose::Bool, - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:AbstractFloat, Tv} # Extend data for exclusive endpoint x, y_mat = _prepare_periodic(x, y_mat, bc) @@ -706,7 +706,7 @@ function cubic_interp( extrap::AbstractExtrap=NoExtrap(), autocache::Bool=true, precompute_transpose::Bool=false, - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:AbstractFloat, Tv} n_pts = length(x) @@ -764,7 +764,7 @@ function cubic_interp( extrap::AbstractExtrap=NoExtrap(), autocache::Bool=true, precompute_transpose::Bool=false, - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:Real, Tv} # Compute promoted grid type (Tg may be Int, promotes to Float) Tg_float = float(promote_type(Tg, _real_eltype(Tv))) @@ -782,7 +782,7 @@ function cubic_interp( extrap::AbstractExtrap=NoExtrap(), autocache::Bool=true, precompute_transpose::Bool=false, - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:Real, Tv} Tg_float = float(promote_type(Tg, _real_eltype(Tv))) Tv_float = _value_type(Tv, Tg_float) @@ -797,7 +797,7 @@ end # ======================================== """ - (sitp::CubicSeriesInterpolant)(xq::Real; deriv=EvalValue(), search=Binary()) + (sitp::CubicSeriesInterpolant)(xq::Real; deriv=EvalValue(), search=LinearBinary()) Evaluate multi-Y interpolant at scalar query point (out-of-place). @@ -826,7 +826,7 @@ function (sitp::CubicSeriesInterpolant{Tg,Tv})( end """ - (sitp::CubicSeriesInterpolant)(output::AbstractVector, xq::Real; deriv=EvalValue(), search=Binary()) + (sitp::CubicSeriesInterpolant)(output::AbstractVector, xq::Real; deriv=EvalValue(), search=LinearBinary()) Evaluate multi-Y interpolant at scalar query point (in-place). diff --git a/src/cubic/cubic_types.jl b/src/cubic/cubic_types.jl index 1d2c1838..a1cd36db 100644 --- a/src/cubic/cubic_types.jl +++ b/src/cubic/cubic_types.jl @@ -141,7 +141,7 @@ struct CubicInterpolant{Tg<:AbstractFloat,Tv,C<:CubicSplineCache{Tg},E<:Abstract z::AbstractVector{Tv}, bc::BC, extrap::E, - search::P=Binary() + search::P=LinearBinary() ) where {Tg<:AbstractFloat, Tv, C<:CubicSplineCache{Tg}, E<:AbstractExtrap, P<:AbstractSearchPolicy, BC<:CubicBC} @assert length(cache.x) == length(y) "cache grid and y must have same length" @assert length(cache.x) == length(z) "z coefficients must match grid length" diff --git a/src/cubic/nd/cubic_nd_interpolant.jl b/src/cubic/nd/cubic_nd_interpolant.jl index 7e559138..a408a2b4 100644 --- a/src/cubic/nd/cubic_nd_interpolant.jl +++ b/src/cubic/nd/cubic_nd_interpolant.jl @@ -25,7 +25,7 @@ Create an N-dimensional cubic Hermite interpolant from grid vectors and data arr - `extrap=NoExtrap()`: Extrapolation mode(s). Can be: - Single `AbstractExtrap`: Applied to all axes (`NoExtrap()`, `ConstExtrap()`, `WrapExtrap()`) - `NTuple{N,AbstractExtrap}`: Per-axis modes -- `search=Binary()`: Search policy(s). Can be: +- `search=LinearBinary()`: Search policy(s). Can be: - Single `AbstractSearchPolicy`: Applied to all axes - `NTuple{N,AbstractSearchPolicy}`: Per-axis policies - `coeffs=PreCompute()`: Coefficient computation strategy @@ -62,7 +62,7 @@ function cubic_interp( data::AbstractArray{Tv_raw, N}; bc::Union{AbstractBC, NTuple{N,AbstractBC}}=CubicFit(), extrap::Union{AbstractExtrap, NTuple{N,AbstractExtrap}}=NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=Binary(), + search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=LinearBinary(), coeffs::AbstractCoeffStrategy=PreCompute() ) where {N, Tv_raw} # Zero-allocation type promotion (uses @generated function) diff --git a/src/cubic/nd/cubic_nd_oneshot.jl b/src/cubic/nd/cubic_nd_oneshot.jl index a171964e..e5da5f14 100644 --- a/src/cubic/nd/cubic_nd_oneshot.jl +++ b/src/cubic/nd/cubic_nd_oneshot.jl @@ -32,7 +32,7 @@ function cubic_interp( deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue(), bc::Union{AbstractBC, NTuple{N,AbstractBC}}=CubicFit(), extrap::Union{AbstractExtrap, NTuple{N,AbstractExtrap}}=NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=Binary(), + search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=LinearBinary(), coeffs::AbstractCoeffStrategy=PreCompute() ) where {Tv, N} # Type promotion + validation (same as constructor path) @@ -66,7 +66,7 @@ function cubic_interp( deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue(), bc::Union{AbstractBC, NTuple{N,AbstractBC}}=CubicFit(), extrap::Union{AbstractExtrap, NTuple{N,AbstractExtrap}}=NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=Binary(), + search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=LinearBinary(), coeffs::AbstractCoeffStrategy=PreCompute() ) where {Tv, N} Tg = _promote_grid_eltype(grids) @@ -90,7 +90,7 @@ function cubic_interp( deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue(), bc::Union{AbstractBC, NTuple{N,AbstractBC}}=CubicFit(), extrap::Union{AbstractExtrap, NTuple{N,AbstractExtrap}}=NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=Binary(), + search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=LinearBinary(), coeffs::AbstractCoeffStrategy=PreCompute() ) where {Tv, N} Tg = _promote_grid_eltype(grids) @@ -258,7 +258,7 @@ function cubic_interp!( deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue(), bc::Union{AbstractBC, NTuple{N,AbstractBC}}=CubicFit(), extrap::Union{AbstractExtrap, NTuple{N,AbstractExtrap}}=NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=Binary(), + search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=LinearBinary(), coeffs::AbstractCoeffStrategy=PreCompute() ) where {Tv, N} Tg = _promote_grid_eltype(grids) @@ -290,7 +290,7 @@ function cubic_interp!( deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue(), bc::Union{AbstractBC, NTuple{N,AbstractBC}}=CubicFit(), extrap::Union{AbstractExtrap, NTuple{N,AbstractExtrap}}=NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=Binary(), + search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=LinearBinary(), coeffs::AbstractCoeffStrategy=PreCompute() ) where {Tv, N} Tg = _promote_grid_eltype(grids) diff --git a/src/linear/linear_interpolant.jl b/src/linear/linear_interpolant.jl index c002ace1..3524d41b 100644 --- a/src/linear/linear_interpolant.jl +++ b/src/linear/linear_interpolant.jl @@ -61,7 +61,7 @@ end # ======================================== """ - linear_interp(x, y; extrap=NoExtrap(), search=Binary()) -> LinearInterpolant + linear_interp(x, y; extrap=NoExtrap(), search=LinearBinary()) -> LinearInterpolant Create a callable interpolant for broadcast fusion and reuse. @@ -69,7 +69,7 @@ Create a callable interpolant for broadcast fusion and reuse. - `x::AbstractVector`: x-coordinates (must be sorted) - `y::AbstractVector`: y-values (can be real or complex) - `extrap::AbstractExtrap`: `NoExtrap()` (default), `ConstExtrap()`, `ExtendExtrap()`, or `WrapExtrap()` -- `search::AbstractSearchPolicy`: Default search policy for interval lookup (default: `Binary()`) +- `search::AbstractSearchPolicy`: Default search policy for interval lookup (default: `LinearBinary()`) # Type Handling - x: Grid coordinates → converted to AbstractFloat, Range structure preserved @@ -91,7 +91,7 @@ Can be: # Examples ```julia -# Create with default Binary() search policy +# Create with default LinearBinary() search policy itp = linear_interp(x_data, y_data) # Create with LinearBinary() as default policy (optimal for sorted queries) @@ -143,7 +143,7 @@ function linear_interp end x::AbstractVector{TX}, y::AbstractVector{TY}; extrap::AbstractExtrap=NoExtrap(), - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {TX<:Real, TY} x_p, y_p = _promote_itp_inputs(x, y) return LinearInterpolant(x_p, y_p; extrap, search) diff --git a/src/linear/linear_oneshot.jl b/src/linear/linear_oneshot.jl index 79c8f908..f4150e9f 100644 --- a/src/linear/linear_oneshot.jl +++ b/src/linear/linear_oneshot.jl @@ -16,7 +16,7 @@ # ======================================== """ - linear_interp!(output, x, y, x_targets; extrap=NoExtrap(), deriv=EvalValue(), search=Binary()) + linear_interp!(output, x, y, x_targets; extrap=NoExtrap(), deriv=EvalValue(), search=LinearBinary()) Zero-allocation linear interpolation with automatic dispatch: - For `AbstractRange` x: O(1) direct indexing @@ -62,7 +62,7 @@ function linear_interp!( x_targets::AbstractVector{Tg}; extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:AbstractFloat, Tv} @assert length(y) == length(x) "x and y must have same length" @assert length(output) == length(x_targets) "output must match x_targets length" @@ -155,7 +155,7 @@ end x_targets::AbstractVector{Tg}; extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:AbstractFloat, Tv} @assert length(y) == length(x) "x and y must have same length" @assert length(output) == length(x_targets) "output must match x_targets length" @@ -170,7 +170,7 @@ end # ======================================== """ - linear_interp(x, y, xq::Real; extrap=NoExtrap(), deriv=EvalValue(), search=Binary()) -> AbstractFloat + linear_interp(x, y, xq::Real; extrap=NoExtrap(), deriv=EvalValue(), search=LinearBinary()) -> AbstractFloat Zero-allocation scalar linear interpolation with automatic dispatch: - For `AbstractRange` x: O(1) direct indexing @@ -385,7 +385,7 @@ end xq::Tq; extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), - search=Binary(), + search=LinearBinary(), hint::Union{Nothing,Base.RefValue{Int}}=nothing ) where {Tg<:AbstractFloat, Tv, Tq<:Real} @boundscheck length(y) == length(x) || throw(ArgumentError("x and y must have same length")) @@ -412,7 +412,7 @@ function linear_interp( x_targets::AbstractVector{Tg}; extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:AbstractFloat, Tv} output = Vector{Tv}(undef, length(x_targets)) linear_interp!(output, x, y, x_targets; extrap, deriv, search) @@ -433,7 +433,7 @@ function linear_interp!( x_targets::AbstractVector{Tq}; extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:Real, Tv, Tq<:Real} @assert length(y) == length(x) "x and y must have same length" @assert length(output) == length(x_targets) "output must match x_targets length" @@ -465,7 +465,7 @@ end xq::Tq; extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), - search=Binary(), + search=LinearBinary(), hint::Union{Nothing,Base.RefValue{Int}}=nothing ) where {Tg<:Real, Tv, Tq<:Real} x_typed, y_typed = _promote_itp_inputs(x, y) @@ -483,7 +483,7 @@ function linear_interp( x_targets::AbstractVector{Tq}; extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:Real, Tv, Tq<:Real} x_typed, y_typed, xq_typed = _promote_itp_inputs(x, y, x_targets) Tv_float = eltype(y_typed) diff --git a/src/linear/linear_series_interp.jl b/src/linear/linear_series_interp.jl index 714fe01f..f42bc286 100644 --- a/src/linear/linear_series_interp.jl +++ b/src/linear/linear_series_interp.jl @@ -78,7 +78,7 @@ mutable struct LinearSeriesInterpolant{Tg<:AbstractFloat, Tv, E<:AbstractExtrap, x::X, y::Matrix{Tv}, extrap::E, - search::P=Binary() + search::P=LinearBinary() ) where {Tg<:AbstractFloat, Tv, E<:AbstractExtrap, P<:AbstractSearchPolicy, X<:AbstractVector{Tg}} new{Tg,Tv,E,P,X}(x, y, LazyTranspose{Tv}(), extrap, search) end @@ -321,7 +321,7 @@ function linear_interp( x::AbstractVector{Tg}, ys::AbstractVector{<:AbstractVector{Tv}}; extrap::AbstractExtrap=NoExtrap(), - search::P=Binary() + search::P=LinearBinary() ) where {Tg<:AbstractFloat, Tv, P<:AbstractSearchPolicy} # Check if Tv's float base requires grid widening (not for Int types) # Int-based types (Complex{Int}) are handled by internal _value_type conversion @@ -385,7 +385,7 @@ function linear_interp( x::AbstractVector{Tg}, Y::AbstractMatrix{Tv}; extrap::AbstractExtrap=NoExtrap(), - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:AbstractFloat, Tv} # Check if Tv's float base requires grid widening Tv_real = _real_eltype(Tv) @@ -421,7 +421,7 @@ function linear_interp( x::AbstractVector{Tg}, ys::AbstractVector{<:AbstractVector{Tv}}; extrap::AbstractExtrap=NoExtrap(), - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:Real, Tv} # Compute promoted grid type (Tg may be Int, promotes to Float) Tg_float = float(promote_type(Tg, _real_eltype(Tv))) @@ -434,7 +434,7 @@ function linear_interp( x::AbstractVector{Tg}, Y::AbstractMatrix{Tv}; extrap::AbstractExtrap=NoExtrap(), - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:Real, Tv} Tg_float = float(promote_type(Tg, _real_eltype(Tv))) x_typed = _to_float(x, Tg_float) @@ -446,7 +446,7 @@ end # ======================================== """ - (sitp::LinearSeriesInterpolant)(xq::Real; deriv=EvalValue(), search=Binary()) + (sitp::LinearSeriesInterpolant)(xq::Real; deriv=EvalValue(), search=LinearBinary()) Evaluate multi-Y interpolant at scalar query point (out-of-place). @@ -465,7 +465,7 @@ function (sitp::LinearSeriesInterpolant{Tg,Tv,P})( end """ - (sitp::LinearSeriesInterpolant)(output::AbstractVector, xq::Real; deriv=EvalValue(), search=Binary()) + (sitp::LinearSeriesInterpolant)(output::AbstractVector, xq::Real; deriv=EvalValue(), search=LinearBinary()) Evaluate multi-Y interpolant at scalar query point (in-place). diff --git a/src/linear/linear_types.jl b/src/linear/linear_types.jl index 130dc4cb..3b298a6c 100644 --- a/src/linear/linear_types.jl +++ b/src/linear/linear_types.jl @@ -27,7 +27,7 @@ Returned by `linear_interp(x, y)` (2-argument form). # Usage ```julia # Create interpolator (minimal allocation) -itp = linear_interp(x, y) # default extrap=NoExtrap(), search=Binary() +itp = linear_interp(x, y) # default extrap=NoExtrap(), search=LinearBinary() # Create with custom search policy (baked-in default) itp = linear_interp(x, y; search=LinearBinary()) @@ -85,7 +85,7 @@ end x::X, y::Y; extrap::AbstractExtrap=NoExtrap(), - search::P=Binary() + search::P=LinearBinary() ) where {Tg<:AbstractFloat, Tv, X<:AbstractVector{Tg}, Y<:AbstractVector{Tv}, P<:AbstractSearchPolicy} E = typeof(extrap) return LinearInterpolant{Tg,Tv,X,Y,E,P}(x, y, extrap, search) diff --git a/src/linear/nd/linear_nd_interpolant.jl b/src/linear/nd/linear_nd_interpolant.jl index 5f9c8047..56bf2982 100644 --- a/src/linear/nd/linear_nd_interpolant.jl +++ b/src/linear/nd/linear_nd_interpolant.jl @@ -23,7 +23,7 @@ The interpolation is exact at grid points and linearly blended between them. # Keyword Arguments - `extrap=NoExtrap()`: Extrapolation mode (`NoExtrap()`, `ConstExtrap()`, `ExtendExtrap()`, `WrapExtrap()`) or per-axis tuple -- `search=Binary()`: Search policy or per-axis tuple +- `search=LinearBinary()`: Search policy or per-axis tuple # Returns - `LinearInterpolantND{Tg,Tv,N}`: Callable interpolant object @@ -62,7 +62,7 @@ function linear_interp( grids::NTuple{N, AbstractVector}, data::AbstractArray{Tv_raw, N}; extrap::Union{AbstractExtrap, NTuple{N, AbstractExtrap}} = NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = Binary() + search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = LinearBinary() ) where {N, Tv_raw} # Validate grid dimensions _validate_nd_grids(grids, data) diff --git a/src/linear/nd/linear_nd_oneshot.jl b/src/linear/nd/linear_nd_oneshot.jl index 37779bce..7341f02e 100644 --- a/src/linear/nd/linear_nd_oneshot.jl +++ b/src/linear/nd/linear_nd_oneshot.jl @@ -112,7 +112,7 @@ function linear_interp( data::AbstractArray{Tv, N}, query::Tuple{Vararg{Real, N}}; extrap::Union{AbstractExtrap, NTuple{N, AbstractExtrap}} = NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = Binary(), + search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = LinearBinary(), deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue() ) where {Tv, N} Tg = _promote_grid_eltype(grids) @@ -139,7 +139,7 @@ function linear_interp( data::AbstractArray{Tv, N}, queries::NTuple{N, AbstractVector{<:Real}}; extrap::Union{AbstractExtrap, NTuple{N, AbstractExtrap}} = NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = Binary(), + search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = LinearBinary(), deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue() ) where {Tv, N} Tg = _promote_grid_eltype(grids) @@ -161,7 +161,7 @@ function linear_interp( data::AbstractArray{Tv, N}, queries::AbstractVector{<:Tuple{Vararg{Real, N}}}; extrap::Union{AbstractExtrap, NTuple{N, AbstractExtrap}} = NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = Binary(), + search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = LinearBinary(), deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue() ) where {Tv, N} Tg = _promote_grid_eltype(grids) @@ -188,7 +188,7 @@ function linear_interp!( data::AbstractArray{Tv, N}, queries::NTuple{N, AbstractVector{<:Real}}; extrap::Union{AbstractExtrap, NTuple{N, AbstractExtrap}} = NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = Binary(), + search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = LinearBinary(), deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue() ) where {Tv, N} Tg = _promote_grid_eltype(grids) @@ -215,7 +215,7 @@ function linear_interp!( data::AbstractArray{Tv, N}, queries::AbstractVector{<:Tuple{Vararg{Real, N}}}; extrap::Union{AbstractExtrap, NTuple{N, AbstractExtrap}} = NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = Binary(), + search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = LinearBinary(), deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue() ) where {Tv, N} Tg = _promote_grid_eltype(grids) diff --git a/src/quadratic/nd/quadratic_nd_interpolant.jl b/src/quadratic/nd/quadratic_nd_interpolant.jl index 846de4b2..e88724a2 100644 --- a/src/quadratic/nd/quadratic_nd_interpolant.jl +++ b/src/quadratic/nd/quadratic_nd_interpolant.jl @@ -82,7 +82,7 @@ Create an N-dimensional quadratic interpolant from grid vectors and data array. - Single BC: Applied to all axes - `NTuple{N}`: Per-axis BCs - `extrap=NoExtrap()`: Extrapolation mode(s) -- `search=Binary()`: Search policy(s) +- `search=LinearBinary()`: Search policy(s) # Returns - `QuadraticInterpolantND{Tg, Tv, N, ...}`: Callable interpolant object @@ -102,7 +102,7 @@ function quadratic_interp( data::AbstractArray{Tv_raw, N}; bc::Union{AbstractBC, NTuple{N,AbstractBC}}=Left(QuadraticFit()), extrap::Union{AbstractExtrap, NTuple{N,AbstractExtrap}}=NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=Binary() + search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=LinearBinary() ) where {N, Tv_raw} # Zero-allocation type promotion Tg = _promote_grid_eltype(grids) diff --git a/src/quadratic/nd/quadratic_nd_oneshot.jl b/src/quadratic/nd/quadratic_nd_oneshot.jl index 4ad5b5df..67acf973 100644 --- a/src/quadratic/nd/quadratic_nd_oneshot.jl +++ b/src/quadratic/nd/quadratic_nd_oneshot.jl @@ -144,7 +144,7 @@ function quadratic_interp( deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue(), bc::Union{AbstractBC, NTuple{N,AbstractBC}}=Left(QuadraticFit()), extrap::Union{AbstractExtrap, NTuple{N,AbstractExtrap}}=NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=Binary() + search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=LinearBinary() ) where {Tv, N} Tg = _promote_grid_eltype(grids) Tg = Tg <: AbstractFloat ? Tg : Float64 @@ -174,7 +174,7 @@ function quadratic_interp( deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue(), bc::Union{AbstractBC, NTuple{N,AbstractBC}}=Left(QuadraticFit()), extrap::Union{AbstractExtrap, NTuple{N,AbstractExtrap}}=NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=Binary() + search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=LinearBinary() ) where {Tv, N} Tg = _promote_grid_eltype(grids) Tg = Tg <: AbstractFloat ? Tg : Float64 @@ -197,7 +197,7 @@ function quadratic_interp( deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue(), bc::Union{AbstractBC, NTuple{N,AbstractBC}}=Left(QuadraticFit()), extrap::Union{AbstractExtrap, NTuple{N,AbstractExtrap}}=NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=Binary() + search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=LinearBinary() ) where {Tv, N} Tg = _promote_grid_eltype(grids) Tg = Tg <: AbstractFloat ? Tg : Float64 @@ -225,7 +225,7 @@ function quadratic_interp!( deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue(), bc::Union{AbstractBC, NTuple{N,AbstractBC}}=Left(QuadraticFit()), extrap::Union{AbstractExtrap, NTuple{N,AbstractExtrap}}=NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=Binary() + search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=LinearBinary() ) where {Tv, N} Tg = _promote_grid_eltype(grids) Tg = Tg <: AbstractFloat ? Tg : Float64 @@ -255,7 +255,7 @@ function quadratic_interp!( deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue(), bc::Union{AbstractBC, NTuple{N,AbstractBC}}=Left(QuadraticFit()), extrap::Union{AbstractExtrap, NTuple{N,AbstractExtrap}}=NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=Binary() + search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=LinearBinary() ) where {Tv, N} Tg = _promote_grid_eltype(grids) Tg = Tg <: AbstractFloat ? Tg : Float64 diff --git a/src/quadratic/quadratic_interpolant.jl b/src/quadratic/quadratic_interpolant.jl index 561fca83..1a003eec 100644 --- a/src/quadratic/quadratic_interpolant.jl +++ b/src/quadratic/quadratic_interpolant.jl @@ -55,7 +55,7 @@ end # ======================================== """ - quadratic_interp(x, y; bc=Left(QuadraticFit()), extrap=NoExtrap(), search=Binary()) -> QuadraticInterpolant + quadratic_interp(x, y; bc=Left(QuadraticFit()), extrap=NoExtrap(), search=LinearBinary()) -> QuadraticInterpolant Create a callable interpolant for broadcast fusion and reuse. @@ -64,7 +64,7 @@ Create a callable interpolant for broadcast fusion and reuse. - `y::AbstractVector`: y-values (can be Real or Complex) - `bc`: Boundary condition (Left, Right, MinCurvFit, or Left/Right with QuadraticFit) - `extrap::AbstractExtrap`: `NoExtrap()` (default), `ConstExtrap()`, `ExtendExtrap()`, or `WrapExtrap()` -- `search::AbstractSearchPolicy`: Default search policy (default: `Binary()`) +- `search::AbstractSearchPolicy`: Default search policy (default: `LinearBinary()`) # Returns `QuadraticInterpolant{Tg, Tv}` object for scalar/broadcast evaluation. @@ -115,7 +115,7 @@ end y::AbstractVector{TY}; bc::QuadraticBC=Left(QuadraticFit()), extrap::AbstractExtrap=NoExtrap(), - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {TX<:Real, TY} x_p, y_p = _promote_itp_inputs(x, y) bc_p = _promote_bc(bc, eltype(x_p)) diff --git a/src/quadratic/quadratic_oneshot.jl b/src/quadratic/quadratic_oneshot.jl index edb0ff51..3b81d318 100644 --- a/src/quadratic/quadratic_oneshot.jl +++ b/src/quadratic/quadratic_oneshot.jl @@ -149,7 +149,7 @@ end # ======================================== """ - quadratic_interp(x, y, xi; bc=Left(QuadraticFit()), extrap=NoExtrap(), deriv=EvalValue(), search=Binary()) + quadratic_interp(x, y, xi; bc=Left(QuadraticFit()), extrap=NoExtrap(), deriv=EvalValue(), search=LinearBinary()) C1 piecewise quadratic spline interpolation at a single point. @@ -202,7 +202,7 @@ vals = quadratic_interp(x, y, sorted_queries; search=LinearBinary(linear_window= bc::QuadraticBC=Left(QuadraticFit()), extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), - search=Binary(), + search=LinearBinary(), hint::Union{Nothing,Base.RefValue{Int}}=nothing ) where {Tg<:AbstractFloat, Tv, Tq<:Real} @boundscheck length(y) == length(x) || throw(ArgumentError("x and y must have same length")) @@ -226,7 +226,7 @@ end # ======================================== """ - quadratic_interp!(output, x, y, x_targets; bc=Left(QuadraticFit()), extrap=NoExtrap(), deriv=EvalValue(), search=Binary()) + quadratic_interp!(output, x, y, x_targets; bc=Left(QuadraticFit()), extrap=NoExtrap(), deriv=EvalValue(), search=LinearBinary()) In-place quadratic spline interpolation for multiple query points. @@ -256,7 +256,7 @@ quadratic_interp!(output, x, y, sorted_queries; search=LinearBinary(linear_windo bc::QuadraticBC=Left(QuadraticFit()), extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:AbstractFloat, Tv} @assert length(y) == length(x) "x and y must have same length" @assert length(output) == length(x_targets) "output must match x_targets length" @@ -284,7 +284,7 @@ end # ======================================== """ - quadratic_interp(x, y, x_targets; bc=Left(QuadraticFit()), extrap=NoExtrap(), deriv=EvalValue(), search=Binary()) + quadratic_interp(x, y, x_targets; bc=Left(QuadraticFit()), extrap=NoExtrap(), deriv=EvalValue(), search=LinearBinary()) Quadratic spline interpolation for multiple query points (allocating version). @@ -307,7 +307,7 @@ function quadratic_interp( bc::QuadraticBC=Left(QuadraticFit()), extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:AbstractFloat, Tv} output = Vector{Tv}(undef, length(x_targets)) quadratic_interp!(output, x, y, x_targets; bc, extrap, deriv, search) @@ -362,7 +362,7 @@ end bc::QuadraticBC=Left(QuadraticFit()), extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), - search=Binary(), + search=LinearBinary(), hint::Union{Nothing,Base.RefValue{Int}}=nothing ) where {Tg<:Real, Tv, Tq<:Real} x_typed, y_typed = _promote_itp_inputs(x, y) @@ -384,7 +384,7 @@ function quadratic_interp( bc::QuadraticBC=Left(QuadraticFit()), extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:Real, Tv, Tq<:Real} x_typed, y_typed, xq_typed = _promote_itp_inputs(x, y, x_targets) Tv_float = eltype(y_typed) @@ -406,7 +406,7 @@ function quadratic_interp!( bc::QuadraticBC=Left(QuadraticFit()), extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:Real, Tv, Tq<:Real} @assert length(y) == length(x) "x and y must have same length" @assert length(output) == length(x_targets) "output must match x_targets length" diff --git a/src/quadratic/quadratic_series_interp.jl b/src/quadratic/quadratic_series_interp.jl index 1f378ac8..aec34a81 100644 --- a/src/quadratic/quadratic_series_interp.jl +++ b/src/quadratic/quadratic_series_interp.jl @@ -93,7 +93,7 @@ mutable struct QuadraticSeriesInterpolant{Tg<:AbstractFloat, Tv, E<:AbstractExtr d::Matrix{Tv}, h::Vector{Tg}, extrap::E, - search::P=Binary() + search::P=LinearBinary() ) where {Tg<:AbstractFloat, Tv, E<:AbstractExtrap, P<:AbstractSearchPolicy, X<:AbstractVector{Tg}} new{Tg,Tv,E,P,X}(x, y, a, d, h, LazyTransposeTriple{Tv}(), extrap, search) end @@ -367,7 +367,7 @@ function quadratic_interp( ys::AbstractVector{<:AbstractVector{Tv}}; bc::QuadraticBC=Left(QuadraticFit()), extrap::AbstractExtrap=NoExtrap(), - search::P=Binary() + search::P=LinearBinary() ) where {Tg<:AbstractFloat, Tv, P<:AbstractSearchPolicy} # Check if Tv's float base requires grid widening (not for Int types) # Int-based types (Complex{Int}) are handled by internal _value_type conversion @@ -449,7 +449,7 @@ function quadratic_interp( Y::AbstractMatrix{Tv}; bc::QuadraticBC=Left(QuadraticFit()), extrap::AbstractExtrap=NoExtrap(), - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:AbstractFloat, Tv} ys = [Y[:, k] for k in axes(Y, 2)] return quadratic_interp(x, ys; bc=bc, extrap=extrap, search=search) @@ -466,7 +466,7 @@ function quadratic_interp( ys::AbstractVector{<:AbstractVector{Tv}}; bc=Left(QuadraticFit()), extrap::AbstractExtrap=NoExtrap(), - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:Real, Tv} # Compute promoted grid type (Tg may be Int, promotes to Float) Tg_float = float(promote_type(Tg, _real_eltype(Tv))) @@ -481,7 +481,7 @@ function quadratic_interp( Y::AbstractMatrix{Tv}; bc=Left(QuadraticFit()), extrap::AbstractExtrap=NoExtrap(), - search::AbstractSearchPolicy=Binary() + search::AbstractSearchPolicy=LinearBinary() ) where {Tg<:Real, Tv} Tg_float = float(promote_type(Tg, _real_eltype(Tv))) x_typed = _to_float(x, Tg_float) @@ -495,7 +495,7 @@ end # Scalar evaluation (explicit implementation for deriv keyword support) """ - (sitp::QuadraticSeriesInterpolant)(xq::Real; deriv=EvalValue(), search=Binary()) + (sitp::QuadraticSeriesInterpolant)(xq::Real; deriv=EvalValue(), search=LinearBinary()) Evaluate all series at scalar query point (out-of-place). @@ -520,7 +520,7 @@ function (sitp::QuadraticSeriesInterpolant{Tg,Tv,P})( end """ - (sitp::QuadraticSeriesInterpolant)(output::AbstractVector, xq::Real; deriv=EvalValue(), search=Binary()) + (sitp::QuadraticSeriesInterpolant)(output::AbstractVector, xq::Real; deriv=EvalValue(), search=LinearBinary()) Evaluate all series at scalar query point (in-place, zero allocation). diff --git a/src/quadratic/quadratic_types.jl b/src/quadratic/quadratic_types.jl index 7e4cd9b1..074d3741 100644 --- a/src/quadratic/quadratic_types.jl +++ b/src/quadratic/quadratic_types.jl @@ -88,7 +88,7 @@ end a::Vector{Tv}, d::Vector{Tv}; extrap::AbstractExtrap=NoExtrap(), - search::P=Binary() + search::P=LinearBinary() ) where {Tg<:AbstractFloat, Tv, X<:AbstractVector{Tg}, Y<:AbstractVector{Tv}, P<:AbstractSearchPolicy} E = typeof(extrap) return QuadraticInterpolant{Tg,Tv,X,Y,E,P}(x, y, h, a, d, extrap, search) From d10948aa341c25ff6f4a6683cff981c1b0da86f6 Mon Sep 17 00:00:00 2001 From: Min-Gu Yoo Date: Mon, 23 Feb 2026 16:15:32 -0800 Subject: [PATCH 03/16] test: update expected defaults to LinearBinary{2} Update test assertions to reflect the new default search policy: - LinearBinary() now produces LinearBinary{2} (was {8}) - Default interpolant search_policy is LinearBinary{2} (was Binary) - Show/format tests updated accordingly --- test/test_search.jl | 12 ++++++------ test/test_show.jl | 4 ++-- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/test/test_search.jl b/test/test_search.jl index e9a383fa..16d159cb 100644 --- a/test/test_search.jl +++ b/test/test_search.jl @@ -711,8 +711,8 @@ using FastInterpolations: search_interval, _search_binary, _search_direct, _sear @test LinearBinary(4) isa LinearBinary{4} @test LinearBinary(16) isa LinearBinary{16} - # Default is 8 - @test LinearBinary() isa LinearBinary{8} + # Default is 2 + @test LinearBinary() isa LinearBinary{2} end @testset "Invalid linear_window Throws ArgumentError" begin @@ -767,15 +767,15 @@ using FastInterpolations: search_interval, _search_binary, _search_direct, _sear @testset "Single Interpolant: stored policy is used by default" begin # Create with LinearBinary as default itp_lb = linear_interp(x, y; search=LinearBinary()) - @test itp_lb.search_policy isa LinearBinary{8} + @test itp_lb.search_policy isa LinearBinary{2} # Create with HintedBinary as default itp_hb = linear_interp(x, y; search=HintedBinary()) @test itp_hb.search_policy isa HintedBinary - # Create with default (Binary) + # Create with default (LinearBinary) itp_bin = linear_interp(x, y) - @test itp_bin.search_policy isa Binary + @test itp_bin.search_policy isa LinearBinary{2} end @testset "Baked-in policy with hint: hint updates when policy supports it" begin @@ -846,7 +846,7 @@ using FastInterpolations: search_interval, _search_binary, _search_direct, _sear @testset "LinearSeriesInterpolant stored policy" begin sitp = linear_interp(x, [y1, y2]; search=LinearBinary()) - @test sitp.search_policy isa LinearBinary{8} + @test sitp.search_policy isa LinearBinary{2} # Scalar call uses stored policy hint = Ref(400) diff --git a/test/test_show.jl b/test/test_show.jl index d8fc132b..a893b79c 100644 --- a/test/test_show.jl +++ b/test/test_show.jl @@ -211,7 +211,7 @@ @test occursin("HintedBinary", verbose_hinted) verbose_lb = sprint(show, MIME("text/plain"), itp_linear_binary) - @test occursin("LinearBinary{8}", verbose_lb) + @test occursin("LinearBinary{2}", verbose_lb) verbose_lb16 = sprint(show, MIME("text/plain"), itp_linear_binary_16) @test occursin("LinearBinary{16}", verbose_lb16) @@ -468,7 +468,7 @@ @test FI._format_search(Linear()) == "Linear" @test FI._format_search(Binary()) == "Binary" @test FI._format_search(HintedBinary()) == "HintedBinary" - @test FI._format_search(LinearBinary()) == "LinearBinary{8}" + @test FI._format_search(LinearBinary()) == "LinearBinary{2}" @test FI._format_search(LinearBinary(linear_window=4)) == "LinearBinary{4}" # DerivativeView with unknown parent type (no .x or .cache.x) From 9f158660c9506613cb65425da410351d02de52bb Mon Sep 17 00:00:00 2001 From: Min-Gu Yoo Date: Mon, 23 Feb 2026 16:15:44 -0800 Subject: [PATCH 04/16] bench: add end-to-end Binary vs LinearBinary comparison script MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Tests all 4 interpolation types × 5 query patterns × 2 grid sizes in 3 calling modes (vector, scalar+hint, scalar no-hint). Validates the Binary→LinearBinary default change with real interpolation workloads rather than isolated search benchmarks. --- benchmark/default_search_comparison.jl | 259 +++++++++++++++++++++++++ 1 file changed, 259 insertions(+) create mode 100644 benchmark/default_search_comparison.jl diff --git a/benchmark/default_search_comparison.jl b/benchmark/default_search_comparison.jl new file mode 100644 index 00000000..614731d0 --- /dev/null +++ b/benchmark/default_search_comparison.jl @@ -0,0 +1,259 @@ +# ═══════════════════════════════════════════════════════════════════════════════ +# Default Search Policy Comparison: Binary() vs LinearBinary(2) +# ═══════════════════════════════════════════════════════════════════════════════ +# +# End-to-end benchmark measuring actual interpolation performance +# with the old default (Binary) vs new default (LinearBinary, window=2). +# +# Tests all 4 interpolation types × 5 query patterns × 2 grid sizes +# in both scalar-loop and vector-call modes. +# +# Run: julia --project benchmark/default_search_comparison.jl +# ═══════════════════════════════════════════════════════════════════════════════ + +using FastInterpolations +using BenchmarkTools +using Printf +using Random + +# ═══════════════════════════════════════════════════════════════════════════════ +# Query Pattern Generators +# ═══════════════════════════════════════════════════════════════════════════════ +# Each returns a Float64 vector of `n` queries within [lo, hi]. +# These simulate real-world usage patterns. + +const PATTERNS = [ + # Random: worst case for LinearBinary — no locality at all + ("Random", (n, lo, hi) -> rand(n) .* (hi - lo) .+ lo), + + # Sorted: typical ODE/PDE integration pattern — monotone increasing + ("Sorted", (n, lo, hi) -> sort(rand(n) .* (hi - lo) .+ lo)), + + # Clustered: physics simulation — queries bunch around a few hotspots + ("Clustered", (n, lo, hi) -> begin + ctrs = range(lo + 0.1 * (hi - lo), hi - 0.1 * (hi - lo), length=5) + w = (hi - lo) * 0.03 + clamp!(vcat([sort(randn(n ÷ 5) .* w .+ c) for c in ctrs]...), lo, hi) + end), + + # Reverse: sorted descending — tests backward walk in LinearBinary + ("Reverse", (n, lo, hi) -> reverse(sort(rand(n) .* (hi - lo) .+ lo))), + + # Dense local: extreme locality — ODE with tiny timestep + ("DenseLocal", (n, lo, hi) -> begin + center = (lo + hi) / 2 + span = (hi - lo) * 0.02 # 2% of domain + sort(rand(n) .* span .+ (center - span / 2)) + end), +] + +# ═══════════════════════════════════════════════════════════════════════════════ +# Benchmark Harness: Scalar Loop (DCE-safe) +# ═══════════════════════════════════════════════════════════════════════════════ +# Accumulates results to prevent dead-code elimination. + +function bench_scalar_loop(itp, queries::Vector{Float64}) + acc = 0.0 + @inbounds for i in eachindex(queries) + acc += itp(queries[i]) + end + return acc +end + +# With external hint — simulates persistent-hint pattern +function bench_scalar_loop_hint(itp, queries::Vector{Float64}, hint::Base.RefValue{Int}) + acc = 0.0 + @inbounds for i in eachindex(queries) + acc += itp(queries[i]; hint=hint) + end + return acc +end + +# ═══════════════════════════════════════════════════════════════════════════════ +# Benchmark Harness: Vector Call (in-place) +# ═══════════════════════════════════════════════════════════════════════════════ + +function bench_vector_call!(itp, out::Vector{Float64}, queries::Vector{Float64}) + itp(out, queries) + return out[1] # prevent DCE +end + +# ═══════════════════════════════════════════════════════════════════════════════ +# Interpolant Factory +# ═══════════════════════════════════════════════════════════════════════════════ +# Creates each interpolation type with the given search policy. + +function make_interpolants(x, y, search) + return ( + Linear = linear_interp(x, y; search), + Cubic = cubic_interp(x, y; search, autocache=false), + Quadratic = quadratic_interp(x, y; search), + Constant = constant_interp(x, y; search), + ) +end + +# ═══════════════════════════════════════════════════════════════════════════════ +# Main Benchmark +# ═══════════════════════════════════════════════════════════════════════════════ + +function main() + Random.seed!(42) + + println("=" ^ 78) + println(" End-to-End Search Policy Benchmark: Binary() vs LinearBinary(2)") + println(" Measuring actual interpolation time (search + kernel), ns per query") + println("=" ^ 78) + + grid_sizes = [500, 2000] + + for n_grid in grid_sizes + x = collect(range(0.0, 10.0, length=n_grid)) + y = sin.(x) .+ 0.1 .* cos.(3.0 .* x) # non-trivial function + nq = n_grid * 2 + + lo = x[2] # avoid exact first grid point + hi = x[end-1] # avoid exact last grid point + + # Build interpolants once per grid size + itps_bin = make_interpolants(x, y, Binary()) + itps_lb = make_interpolants(x, y, LinearBinary()) + + # ── Section Header ── + println() + println("━" ^ 78) + @printf(" Grid = %d points, Queries = %d\n", n_grid, nq) + println("━" ^ 78) + + # ──────────────────────────────────────────── + # Part 1: Vector call (itp(out, queries)) + # ──────────────────────────────────────────── + println() + println(" ▸ Vector Call: itp(out, queries)") + println(" " * "─" ^ 74) + @printf(" %-11s %-10s %9s %9s %8s\n", + "Pattern", "InterpType", "Binary", "LB{2}", "Speedup") + println(" " * "─" ^ 74) + + for (pname, pfn) in PATTERNS + queries = pfn(nq, lo, hi) + out = Vector{Float64}(undef, nq) + + for (itp_name, itp_b, itp_l) in [ + ("Linear", itps_bin.Linear, itps_lb.Linear), + ("Cubic", itps_bin.Cubic, itps_lb.Cubic), + ("Quadratic", itps_bin.Quadratic, itps_lb.Quadratic), + ("Constant", itps_bin.Constant, itps_lb.Constant), + ] + # Warmup + bench_vector_call!(itp_b, out, queries) + bench_vector_call!(itp_l, out, queries) + + t_b = median(@benchmark bench_vector_call!($itp_b, $out, $queries) samples=60 evals=10).time + t_l = median(@benchmark bench_vector_call!($itp_l, $out, $queries) samples=60 evals=10).time + + speedup = t_b / t_l + marker = speedup > 1.05 ? "▲" : speedup < 0.95 ? "▼" : "=" + + @printf(" %-11s %-10s %7.1f ns %7.1f ns %5.2fx %s\n", + itp_name == "Linear" ? pname : "", + itp_name, + t_b / nq, t_l / nq, + speedup, marker) + end + end + + # ──────────────────────────────────────────── + # Part 2: Scalar loop with persistent hint + # ──────────────────────────────────────────── + println() + println(" ▸ Scalar Loop with Persistent Hint: for q in queries; itp(q; hint=hint); end") + println(" " * "─" ^ 74) + @printf(" %-11s %-10s %9s %9s %8s\n", + "Pattern", "InterpType", "Binary", "LB{2}", "Speedup") + println(" " * "─" ^ 74) + + for (pname, pfn) in PATTERNS + queries = pfn(nq, lo, hi) + + # Only test Linear and Cubic for scalar loop (representative) + for (itp_name, itp_b, itp_l) in [ + ("Linear", itps_bin.Linear, itps_lb.Linear), + ("Cubic", itps_bin.Cubic, itps_lb.Cubic), + ] + hint_b = Ref(1) + hint_l = Ref(1) + + # Warmup + bench_scalar_loop_hint(itp_b, queries, hint_b) + bench_scalar_loop_hint(itp_l, queries, hint_l) + + t_b = median(@benchmark bench_scalar_loop_hint($itp_b, $queries, $hint_b) samples=60 evals=10).time + t_l = median(@benchmark bench_scalar_loop_hint($itp_l, $queries, $hint_l) samples=60 evals=10).time + + speedup = t_b / t_l + marker = speedup > 1.05 ? "▲" : speedup < 0.95 ? "▼" : "=" + + @printf(" %-11s %-10s %7.1f ns %7.1f ns %5.2fx %s\n", + itp_name == "Linear" ? pname : "", + itp_name, + t_b / nq, t_l / nq, + speedup, marker) + end + end + + # ──────────────────────────────────────────── + # Part 3: Scalar loop WITHOUT hint (broadcast-like) + # ──────────────────────────────────────────── + # This is the itp(q) call — no external hint. + # For Binary(), there's no hint at all. + # For LinearBinary(), a fresh internal hint is created per call. + println() + println(" ▸ Scalar Loop (no hint): for q in queries; itp(q); end") + println(" " * "─" ^ 74) + @printf(" %-11s %-10s %9s %9s %8s\n", + "Pattern", "InterpType", "Binary", "LB{2}", "Speedup") + println(" " * "─" ^ 74) + + for (pname, pfn) in PATTERNS + queries = pfn(nq, lo, hi) + + for (itp_name, itp_b, itp_l) in [ + ("Linear", itps_bin.Linear, itps_lb.Linear), + ("Cubic", itps_bin.Cubic, itps_lb.Cubic), + ] + # Warmup + bench_scalar_loop(itp_b, queries) + bench_scalar_loop(itp_l, queries) + + t_b = median(@benchmark bench_scalar_loop($itp_b, $queries) samples=60 evals=10).time + t_l = median(@benchmark bench_scalar_loop($itp_l, $queries) samples=60 evals=10).time + + speedup = t_b / t_l + marker = speedup > 1.05 ? "▲" : speedup < 0.95 ? "▼" : "=" + + @printf(" %-11s %-10s %7.1f ns %7.1f ns %5.2fx %s\n", + itp_name == "Linear" ? pname : "", + itp_name, + t_b / nq, t_l / nq, + speedup, marker) + end + end + end + + # ── Legend ── + println() + println("=" ^ 78) + println(" Legend:") + println(" ▲ = LinearBinary faster ▼ = Binary faster = = within 5%") + println(" Speedup = Binary_time / LB_time (>1 means LB wins)") + println() + println(" Expected behavior:") + println(" Random: LB ≈ 0.5x (2x slower) — structural penalty from hint walk") + println(" Sorted: LB ≈ 5-10x faster — hint almost always hits") + println(" Clustered: LB ≈ 3-8x faster — local queries exploit hint") + println(" Reverse: LB ≈ 5-10x faster — backward walk exploits hint") + println(" DenseLocal: LB ≈ 10-50x faster — nearly all queries hit hint") + println("=" ^ 78) +end + +main() From 3d9e66464842ee4aa2f246feed5ac4a9887304aa Mon Sep 17 00:00:00 2001 From: Min-Gu Yoo Date: Tue, 24 Feb 2026 12:41:56 -0800 Subject: [PATCH 05/16] Add AutoSearch sentinel type with query-type adaptive resolution MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit AutoSearch resolves at call time: scalar→Binary(), vector→LinearBinary(). Includes _resolve_search dispatch, _to_searcher fallbacks, export, and show format. --- src/FastInterpolations.jl | 2 +- src/core/search.jl | 74 +++++++++++++++++++++++++++++++++++---- src/core/show.jl | 1 + 3 files changed, 70 insertions(+), 7 deletions(-) diff --git a/src/FastInterpolations.jl b/src/FastInterpolations.jl index 767ba59a..2f1f9201 100644 --- a/src/FastInterpolations.jl +++ b/src/FastInterpolations.jl @@ -42,7 +42,7 @@ export precompute_transpose! # Pre-allocate point-contiguous layout for scalar export set_cubic_cache_size!, get_cubic_cache_size, clear_cubic_cache! # Search policy types (user-facing algorithm selection) -export AbstractSearchPolicy, Binary, HintedBinary, Linear, LinearBinary +export AbstractSearchPolicy, Binary, HintedBinary, Linear, LinearBinary, AutoSearch # Boundary condition types export AbstractBC, PointBC, Deriv1, Deriv2, Deriv3, BCPair diff --git a/src/core/search.jl b/src/core/search.jl index f0835aff..83e881f3 100644 --- a/src/core/search.jl +++ b/src/core/search.jl @@ -192,6 +192,56 @@ function LinearBinary(linear_window::Integer) end LinearBinary(; linear_window::Integer=2) = LinearBinary(linear_window) +""" + AutoSearch <: AbstractSearchPolicy + +Adaptive search policy that resolves at call time based on query type: +- **Scalar** queries (`Real`, `Tuple{Vararg{Real}}`) → `Binary()` — no hint locality to exploit +- **Vector** queries (`AbstractVector`, `Tuple{Vararg{AbstractVector}}`) → `LinearBinary()` — hint continuity benefits sorted sequences +- **Broadcast** (`itp.(xs)`) → `Binary()` per element — fresh searcher each call + +This is the default search policy for all interpolants. For known query patterns, +specify `Binary()` or `LinearBinary()` explicitly for optimal performance. + +# Example +```julia +itp = linear_interp(x, y) # stores AutoSearch() +itp(0.5) # → Binary() internally (scalar) +itp([0.1, 0.5, 0.9]) # → LinearBinary() internally (vector) +itp(0.5; search=LinearBinary()) # override: use LinearBinary explicitly +``` + +See also: [`Binary`](@ref), [`LinearBinary`](@ref) +""" +struct AutoSearch <: AbstractSearchPolicy end + +# ---------------------------------------- +# AutoSearch Resolution (query-type adaptive) +# ---------------------------------------- +# Resolves AutoSearch to a concrete policy based on query type. +# For non-AutoSearch policies, passes through unchanged (user override honored). +# Must be called BEFORE _to_searcher in all eval paths. + +# 1D scalar: no hint locality → pure binary search +@inline _resolve_search(::AutoSearch, ::Real) = Binary() + +# 1D vector: sorted locality → linear window with binary fallback +@inline _resolve_search(::AutoSearch, ::AbstractVector) = LinearBinary() + +# ND SoA batch: tuple of vectors → LinearBinary per axis +# NOTE: must be above Tuple{Vararg{Real}} to avoid ambiguity (AbstractVector <: Any) +@inline _resolve_search(::AutoSearch, ::Tuple{Vararg{AbstractVector}}) = LinearBinary() + +# ND scalar: tuple of reals → Binary per axis +@inline _resolve_search(::AutoSearch, ::Tuple) = Binary() + +# Passthrough: explicit policies are honored as-is +@inline _resolve_search(p::AbstractSearchPolicy, _) = p + +# Tuple of policies: resolve each element (for ND per-axis storage) +@inline _resolve_search(ps::Tuple{Vararg{AbstractSearchPolicy}}, query) = + map(p -> _resolve_search(p, query), ps) + # ---------------------------------------- # Hint Types (Internal) # ---------------------------------------- @@ -296,6 +346,12 @@ Creates a new RefHint for stateful policies, ensuring thread safety. @inline _to_searcher(::LinearBinary{MAX}, ::Nothing) where {MAX} = Searcher{LinearBinary{MAX},RefHint}(RefHint()) @inline _to_searcher(::LinearBinary{MAX}, hint::Base.RefValue{Int}) where {MAX} = Searcher{LinearBinary{MAX},RefHint}(RefHint(hint)) +# AutoSearch fallbacks: _resolve_search should be called first, but if any +# code path misses resolution, fall back to Binary (safe stateless default). +@inline _to_searcher(::AutoSearch) = Searcher{Binary,NoHint}(NoHint()) +@inline _to_searcher(::AutoSearch, ::Nothing) = Searcher{Binary,NoHint}(NoHint()) +@inline _to_searcher(::AutoSearch, hint::Base.RefValue{Int}) = Searcher{HintedBinary,RefHint}(RefHint(hint)) + # ---------------------------------------- # Searcher passthrough (advanced usage) # ---------------------------------------- @@ -355,6 +411,10 @@ end _search_binary(x::AbstractVector{T}, xq::T) where {T<:Real} O(log n) binary search for non-uniform grids (AbstractVector). + +Uses branchless `for` loop with precomputed iteration count via `leading_zeros` +for predictable loop exit on modern CPUs. The inner comparison uses `ifelse` to +compile to ARM64 `csel` / x86 `cmov` — fully branchless binary search body. """ @inline function _search_binary(x::AbstractVector{T}, xq::T) where {T<:Real} n = length(x) @@ -365,13 +425,15 @@ O(log n) binary search for non-uniform grids (AbstractVector). idx = n - 1 else lo, hi = 1, n - while hi - lo > 1 + # Precompute exact iteration count: ceil(log2(hi - lo)) via CLZ + # n >= 2 guaranteed (grid has ≥2 points), so hi - lo - 1 >= 0 + # Use %UInt64 (bit reinterpret) to avoid InexactError cold path in ASM + iters = 64 - leading_zeros((hi - lo - 1) % UInt64) + for _ in 1:iters mid = (lo + hi) >> 1 - if x[mid] <= xq - lo = mid - else - hi = mid - end + cond = x[mid] <= xq + lo = ifelse(cond, mid, lo) + hi = ifelse(cond, hi, mid) end idx = lo end diff --git a/src/core/show.jl b/src/core/show.jl index 362a7460..0d74cec6 100644 --- a/src/core/show.jl +++ b/src/core/show.jl @@ -122,6 +122,7 @@ _format_search(::Binary) = "Binary" _format_search(::HintedBinary) = "HintedBinary" _format_search(::Linear) = "Linear" _format_search(::LinearBinary{MAX}) where {MAX} = "LinearBinary{$MAX}" +_format_search(::AutoSearch) = "AutoSearch (scalar→Binary, vector→LinearBinary)" """Format boundary condition for display.""" _format_bc(::ZeroCurvBC) = "ZeroCurv (S''=0 at ends)" From de0e48416c6ec40c1eab5b784c7c74782fce3068 Mon Sep 17 00:00:00 2001 From: Min-Gu Yoo Date: Tue, 24 Feb 2026 12:48:04 -0800 Subject: [PATCH 06/16] Change all default search policy from LinearBinary() to AutoSearch() Updates keyword defaults in constructors, oneshot functions, series constructors, and ND constructors across all 4 interpolant types (linear, cubic, quadratic, constant). Docstring signatures updated accordingly. Explicit user examples preserved. --- src/constant/constant_interpolant.jl | 6 ++-- src/constant/constant_oneshot.jl | 18 +++++----- src/constant/constant_series_interp.jl | 14 ++++---- src/constant/constant_types.jl | 2 +- src/constant/nd/constant_nd_interpolant.jl | 4 +-- src/constant/nd/constant_nd_oneshot.jl | 10 +++--- src/cubic/cubic_eval.jl | 2 +- src/cubic/cubic_interpolant.jl | 18 +++++----- src/cubic/cubic_oneshot.jl | 38 ++++++++++---------- src/cubic/cubic_series_interp.jl | 16 ++++----- src/cubic/cubic_types.jl | 2 +- src/cubic/nd/cubic_nd_interpolant.jl | 4 +-- src/cubic/nd/cubic_nd_oneshot.jl | 10 +++--- src/linear/linear_interpolant.jl | 8 ++--- src/linear/linear_oneshot.jl | 18 +++++----- src/linear/linear_series_interp.jl | 14 ++++---- src/linear/linear_types.jl | 6 ++-- src/linear/nd/linear_nd_interpolant.jl | 4 +-- src/linear/nd/linear_nd_oneshot.jl | 10 +++--- src/quadratic/nd/quadratic_nd_interpolant.jl | 4 +-- src/quadratic/nd/quadratic_nd_oneshot.jl | 10 +++--- src/quadratic/quadratic_interpolant.jl | 6 ++-- src/quadratic/quadratic_oneshot.jl | 18 +++++----- src/quadratic/quadratic_series_interp.jl | 14 ++++---- src/quadratic/quadratic_types.jl | 2 +- 25 files changed, 129 insertions(+), 129 deletions(-) diff --git a/src/constant/constant_interpolant.jl b/src/constant/constant_interpolant.jl index c0ec865b..1ea1f1e0 100644 --- a/src/constant/constant_interpolant.jl +++ b/src/constant/constant_interpolant.jl @@ -51,7 +51,7 @@ end # ======================================== """ - constant_interp(x, y; extrap=NoExtrap(), side=NearestSide(), search=LinearBinary()) -> ConstantInterpolant + constant_interp(x, y; extrap=NoExtrap(), side=NearestSide(), search=AutoSearch()) -> ConstantInterpolant Create a callable interpolant for broadcast fusion and reuse. @@ -60,7 +60,7 @@ Create a callable interpolant for broadcast fusion and reuse. - `y::AbstractVector`: y-values (can be Real or Complex) - `extrap::AbstractExtrap`: `NoExtrap()` (default), `ConstExtrap()`, `ExtendExtrap()`, or `WrapExtrap()` - `side::AbstractSide`: Side selection (NearestSide(), LeftSide(), RightSide()) -- `search::AbstractSearchPolicy`: Default search policy (default: `LinearBinary()`) +- `search::AbstractSearchPolicy`: Default search policy (default: `AutoSearch()`) # Returns `ConstantInterpolant{Tg, Tv}` object for scalar/broadcast evaluation. @@ -109,7 +109,7 @@ end y::AbstractVector{TY}; extrap::AbstractExtrap=NoExtrap(), side::AbstractSide=NearestSide(), - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {TX<:Real, TY} x_p, y_p = _promote_itp_inputs(x, y) return ConstantInterpolant(x_p, y_p; extrap, side, search) diff --git a/src/constant/constant_oneshot.jl b/src/constant/constant_oneshot.jl index 2cecc9fb..9ea451de 100644 --- a/src/constant/constant_oneshot.jl +++ b/src/constant/constant_oneshot.jl @@ -149,7 +149,7 @@ end # ======================================== """ - constant_interp(x, y, xi; extrap=NoExtrap(), side=NearestSide(), deriv=EvalValue(), search=LinearBinary()) + constant_interp(x, y, xi; extrap=NoExtrap(), side=NearestSide(), deriv=EvalValue(), search=AutoSearch()) Constant (step/piecewise constant) interpolation at a single point. @@ -200,7 +200,7 @@ vals = constant_interp(x, y, sorted_queries; search=LinearBinary(linear_window=8 extrap::AbstractExtrap=NoExtrap(), side::AbstractSide=NearestSide(), deriv::DerivOp=EvalValue(), - search=LinearBinary(), + search=AutoSearch(), hint::Union{Nothing,Base.RefValue{Int}}=nothing ) where {Tg<:AbstractFloat, Tv, Tq<:Real} @boundscheck length(y) == length(x) || throw(ArgumentError("x and y must have same length")) @@ -214,7 +214,7 @@ end # ======================================== """ - constant_interp!(output, x, y, x_targets; extrap=NoExtrap(), side=NearestSide(), deriv=EvalValue(), search=LinearBinary()) + constant_interp!(output, x, y, x_targets; extrap=NoExtrap(), side=NearestSide(), deriv=EvalValue(), search=AutoSearch()) Zero-allocation constant interpolation for multiple query points. @@ -246,7 +246,7 @@ function constant_interp!( extrap::AbstractExtrap=NoExtrap(), side::AbstractSide=NearestSide(), deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:AbstractFloat, Tv} @assert length(y) == length(x) "x and y must have same length" @assert length(output) == length(x_targets) "output must match x_targets length" @@ -264,7 +264,7 @@ end # ======================================== """ - constant_interp(x, y, x_targets; extrap=NoExtrap(), side=NearestSide(), deriv=EvalValue(), search=LinearBinary()) + constant_interp(x, y, x_targets; extrap=NoExtrap(), side=NearestSide(), deriv=EvalValue(), search=AutoSearch()) Constant interpolation for multiple query points (allocating version). @@ -287,7 +287,7 @@ function constant_interp( extrap::AbstractExtrap=NoExtrap(), side::AbstractSide=NearestSide(), deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:AbstractFloat, Tv} output = Vector{Tv}(undef, length(x_targets)) constant_interp!(output, x, y, x_targets; extrap, side, deriv, search) @@ -314,7 +314,7 @@ end extrap::AbstractExtrap=NoExtrap(), side::AbstractSide=NearestSide(), deriv::DerivOp=EvalValue(), - search=LinearBinary(), + search=AutoSearch(), hint::Union{Nothing,Base.RefValue{Int}}=nothing ) where {Tg<:Real, Tv, Tq<:Real} x_typed, y_typed = _promote_itp_inputs(x, y) @@ -334,7 +334,7 @@ function constant_interp( extrap::AbstractExtrap=NoExtrap(), side::AbstractSide=NearestSide(), deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:Real, Tv, Tq<:Real} x_typed, y_typed, xq_typed = _promote_itp_inputs(x, y, x_targets) Tv_float = eltype(y_typed) @@ -355,7 +355,7 @@ function constant_interp!( extrap::AbstractExtrap=NoExtrap(), side::AbstractSide=NearestSide(), deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:Real, Tv, Tq<:Real} @assert length(y) == length(x) "x and y must have same length" @assert length(output) == length(x_targets) "output must match x_targets length" diff --git a/src/constant/constant_series_interp.jl b/src/constant/constant_series_interp.jl index dfc4ce8a..e765e08c 100644 --- a/src/constant/constant_series_interp.jl +++ b/src/constant/constant_series_interp.jl @@ -77,7 +77,7 @@ mutable struct ConstantSeriesInterpolant{Tg<:AbstractFloat, Tv, E<:AbstractExtra y::Matrix{Tv}, extrap::E, side::SD, - search::P=LinearBinary() + search::P=AutoSearch() ) where {Tg<:AbstractFloat, Tv, E<:AbstractExtrap, SD<:AbstractSide, P<:AbstractSearchPolicy, X<:AbstractVector{Tg}} new{Tg,Tv,E,SD,P,X}(x, y, LazyTranspose{Tv}(), extrap, side, search) end @@ -333,7 +333,7 @@ function constant_interp( ys::AbstractVector{<:AbstractVector{Tv}}; side::AbstractSide=NearestSide(), extrap::AbstractExtrap=NoExtrap(), - search::P=LinearBinary() + search::P=AutoSearch() ) where {Tg<:AbstractFloat, Tv, P<:AbstractSearchPolicy} # Check if Tv's float base requires grid widening (not for Int types) # Int-based types (Complex{Int}) are handled by internal _value_type conversion @@ -394,7 +394,7 @@ function constant_interp( Y::AbstractMatrix{Tv}; side::AbstractSide=NearestSide(), extrap::AbstractExtrap=NoExtrap(), - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:AbstractFloat, Tv} # Check if Tv's float base requires grid widening Tv_real = _real_eltype(Tv) @@ -431,7 +431,7 @@ function constant_interp( ys::AbstractVector{<:AbstractVector{Tv}}; side::AbstractSide=NearestSide(), extrap::AbstractExtrap=NoExtrap(), - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:Real, Tv} # Compute promoted grid type (Tg may be Int, promotes to Float) Tg_float = float(promote_type(Tg, _real_eltype(Tv))) @@ -445,7 +445,7 @@ function constant_interp( Y::AbstractMatrix{Tv}; side::AbstractSide=NearestSide(), extrap::AbstractExtrap=NoExtrap(), - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:Real, Tv} Tg_float = float(promote_type(Tg, _real_eltype(Tv))) x_typed = _to_float(x, Tg_float) @@ -457,7 +457,7 @@ end # ======================================== """ - (sitp::ConstantSeriesInterpolant)(xq::Real; deriv=EvalValue(), search=LinearBinary()) + (sitp::ConstantSeriesInterpolant)(xq::Real; deriv=EvalValue(), search=AutoSearch()) Evaluate multi-Y interpolant at scalar query point (out-of-place). @@ -483,7 +483,7 @@ function (sitp::ConstantSeriesInterpolant{Tg,Tv,P})( end """ - (sitp::ConstantSeriesInterpolant)(output::AbstractVector, xq::Real; deriv=EvalValue(), search=LinearBinary()) + (sitp::ConstantSeriesInterpolant)(output::AbstractVector, xq::Real; deriv=EvalValue(), search=AutoSearch()) Evaluate multi-Y interpolant at scalar query point (in-place). """ diff --git a/src/constant/constant_types.jl b/src/constant/constant_types.jl index 10e87394..b61f0aba 100644 --- a/src/constant/constant_types.jl +++ b/src/constant/constant_types.jl @@ -78,7 +78,7 @@ end y::Y; extrap::AbstractExtrap=NoExtrap(), side::AbstractSide=NearestSide(), - search::P=LinearBinary() + search::P=AutoSearch() ) where {Tg<:AbstractFloat, Tv, X<:AbstractVector{Tg}, Y<:AbstractVector{Tv}, P<:AbstractSearchPolicy} E = typeof(extrap) SD = typeof(side) diff --git a/src/constant/nd/constant_nd_interpolant.jl b/src/constant/nd/constant_nd_interpolant.jl index ec0d9c80..17c7ba99 100644 --- a/src/constant/nd/constant_nd_interpolant.jl +++ b/src/constant/nd/constant_nd_interpolant.jl @@ -21,7 +21,7 @@ Create an N-dimensional constant interpolant with tuple-grid API. # Keyword Arguments - `side=NearestSide()`: Side selection mode (`NearestSide()`, `LeftSide()`, `RightSide()`) or per-axis tuple - `extrap=NoExtrap()`: Extrapolation mode (`NoExtrap()`, `ConstExtrap()`, `ExtendExtrap()`, `WrapExtrap()`) or per-axis tuple -- `search=LinearBinary()`: Search policy or per-axis tuple +- `search=AutoSearch()`: Search policy or per-axis tuple # Returns - `ConstantInterpolantND{Tg,Tv,N}`: Callable interpolant object @@ -51,7 +51,7 @@ function constant_interp( data::AbstractArray{Tv_raw, N}; side::Union{AbstractSide, Tuple{Vararg{AbstractSide}}} = NearestSide(), extrap::Union{AbstractExtrap, NTuple{N, AbstractExtrap}} = NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = LinearBinary() + search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = AutoSearch() ) where {N, Tv_raw} # Validate grid dimensions _validate_nd_grids(grids, data) diff --git a/src/constant/nd/constant_nd_oneshot.jl b/src/constant/nd/constant_nd_oneshot.jl index 044ca587..f8286dd6 100644 --- a/src/constant/nd/constant_nd_oneshot.jl +++ b/src/constant/nd/constant_nd_oneshot.jl @@ -150,7 +150,7 @@ function constant_interp( query::Tuple{Vararg{Real, N}}; side::Union{AbstractSide, Tuple{Vararg{AbstractSide}}} = NearestSide(), extrap::Union{AbstractExtrap, NTuple{N, AbstractExtrap}} = NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = LinearBinary(), + search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = AutoSearch(), deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue() ) where {Tv, N} # Any derivative of constant interpolation is zero @@ -183,7 +183,7 @@ function constant_interp( queries::NTuple{N, AbstractVector{<:Real}}; side::Union{AbstractSide, Tuple{Vararg{AbstractSide}}} = NearestSide(), extrap::Union{AbstractExtrap, NTuple{N, AbstractExtrap}} = NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = LinearBinary(), + search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = AutoSearch(), deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue() ) where {Tv, N} if _is_any_deriv(deriv) @@ -216,7 +216,7 @@ function constant_interp( queries::AbstractVector{<:Tuple{Vararg{Real, N}}}; side::Union{AbstractSide, Tuple{Vararg{AbstractSide}}} = NearestSide(), extrap::Union{AbstractExtrap, NTuple{N, AbstractExtrap}} = NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = LinearBinary(), + search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = AutoSearch(), deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue() ) where {Tv, N} if _is_any_deriv(deriv) @@ -254,7 +254,7 @@ function constant_interp!( queries::NTuple{N, AbstractVector{<:Real}}; side::Union{AbstractSide, Tuple{Vararg{AbstractSide}}} = NearestSide(), extrap::Union{AbstractExtrap, NTuple{N, AbstractExtrap}} = NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = LinearBinary(), + search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = AutoSearch(), deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue() ) where {Tv, N} if _is_any_deriv(deriv) @@ -288,7 +288,7 @@ function constant_interp!( queries::AbstractVector{<:Tuple{Vararg{Real, N}}}; side::Union{AbstractSide, Tuple{Vararg{AbstractSide}}} = NearestSide(), extrap::Union{AbstractExtrap, NTuple{N, AbstractExtrap}} = NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = LinearBinary(), + search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = AutoSearch(), deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue() ) where {Tv, N} if _is_any_deriv(deriv) diff --git a/src/cubic/cubic_eval.jl b/src/cubic/cubic_eval.jl index c6a63f13..08d95f51 100644 --- a/src/cubic/cubic_eval.jl +++ b/src/cubic/cubic_eval.jl @@ -273,7 +273,7 @@ Uses task-local pool for workspace allocation. x_query::Tg; extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), - search=LinearBinary(), + search=AutoSearch(), hint::Union{Nothing,Base.RefValue{Int}}=nothing ) where {Tg<:AbstractFloat, Tv, X, F, BC, S<:AbstractGridSpacing{Tg}} @assert length(y) == length(cache.x) "y length must match cache grid" diff --git a/src/cubic/cubic_interpolant.jl b/src/cubic/cubic_interpolant.jl index ecd7b79d..6b2072b8 100644 --- a/src/cubic/cubic_interpolant.jl +++ b/src/cubic/cubic_interpolant.jl @@ -293,7 +293,7 @@ so the pool memory can be safely reused after this function returns. bc_pair::BCPair{L,R}, extrap::AbstractExtrap, autocache::Bool, - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:AbstractFloat, Tv, L<:PointBC, R<:PointBC} # Cache uses structural equivalent (PolyFit → Deriv1 via _cache_bc_pair internally) cache = _get_cubic_cache(x, bc_pair, autocache) @@ -319,7 +319,7 @@ so the pool memory can be safely reused after this function returns. y::AbstractVector{Tv}, bc::PeriodicBC, autocache::Bool, - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:AbstractFloat, Tv} x, y = _prepare_periodic(x, y, bc) _check_periodic_endpoints(y) @@ -351,7 +351,7 @@ Handles conversion of Real BC values to Complex when needed. # ======================================== """ - cubic_interp(x, y; bc=CubicFit(), extrap=NoExtrap(), autocache=true, search=LinearBinary()) -> CubicInterpolant + cubic_interp(x, y; bc=CubicFit(), extrap=NoExtrap(), autocache=true, search=AutoSearch()) -> CubicInterpolant Create a callable interpolant for broadcast fusion and reuse. @@ -364,7 +364,7 @@ enabling true zero-allocation scalar evaluations in broadcast operations. - `bc::AbstractBC`: Boundary condition (default: `CubicFit()`) - `extrap::AbstractExtrap`: `NoExtrap()` (default), `ConstExtrap()`, `ExtendExtrap()`, or `WrapExtrap()` - `autocache::Bool`: Enable automatic caching (default: `true`) -- `search::AbstractSearchPolicy`: Default search policy (default: `LinearBinary()`) +- `search::AbstractSearchPolicy`: Default search policy (default: `AutoSearch()`) # Example ```julia @@ -393,7 +393,7 @@ val = itp(0.5) # returns ComplexF64 bc::AbstractBC, extrap::AbstractExtrap, autocache::Bool, - search::P=LinearBinary() + search::P=AutoSearch() ) where {Tg<:AbstractFloat, Tv, P<:AbstractSearchPolicy} if _is_periodic_bc(bc) return _build_interpolant_periodic(x, y, bc, autocache, search) @@ -410,7 +410,7 @@ function cubic_interp( bc::AbstractBC=CubicFit(), extrap::AbstractExtrap=NoExtrap(), autocache::Bool=true, - search::P=LinearBinary() + search::P=AutoSearch() ) where {Tg<:AbstractFloat, Tv, P<:AbstractSearchPolicy} # Auto-promote x/y types (zero allocation if already compatible) x_p, y_p = _promote_itp_inputs(x, y) @@ -419,7 +419,7 @@ function cubic_interp( end """ - cubic_interp(cache, y; extrap=NoExtrap(), search=LinearBinary()) -> CubicInterpolant + cubic_interp(cache, y; extrap=NoExtrap(), search=AutoSearch()) -> CubicInterpolant Create a callable interpolant from a pre-built cache. @@ -441,7 +441,7 @@ so the pool memory can be safely reused after this function returns. cache::CubicSplineCache{Tg}, y::AbstractVector{Tv}; extrap::AbstractExtrap=NoExtrap(), - search::P=LinearBinary() + search::P=AutoSearch() ) where {Tg<:AbstractFloat, Tv, P<:AbstractSearchPolicy} tmp_z = similar!(pool, y) _solve_system!(tmp_z, cache, y, cache.bc_config) @@ -462,7 +462,7 @@ function cubic_interp( bc::AbstractBC=CubicFit(), extrap::AbstractExtrap=NoExtrap(), autocache::Bool=true, - search::P=LinearBinary() + search::P=AutoSearch() ) where {TX<:Real, TY, P<:AbstractSearchPolicy} x_p, y_p = _promote_itp_inputs(x, y) bc_promoted = _promote_bc(bc, eltype(y_p)) diff --git a/src/cubic/cubic_oneshot.jl b/src/cubic/cubic_oneshot.jl index cf35b427..a37fe206 100644 --- a/src/cubic/cubic_oneshot.jl +++ b/src/cubic/cubic_oneshot.jl @@ -12,7 +12,7 @@ # ======================================== """ - cubic_interp!(output, cache, y, x_query; extrap=NoExtrap(), deriv=EvalValue(), search=LinearBinary()) + cubic_interp!(output, cache, y, x_query; extrap=NoExtrap(), deriv=EvalValue(), search=AutoSearch()) In-place cubic spline interpolation using cached LU factorization. @@ -26,7 +26,7 @@ Thread-safe: workspaces allocated from task-local pool. - `x_query::AbstractVector{T}`: Query points - `extrap::AbstractExtrap`: `NoExtrap()` (default), `ConstExtrap()`, `ExtendExtrap()`, or `WrapExtrap()` - `deriv::DerivOp=EvalValue()`: Derivative order (0=value, 1=first derivative, 2=second derivative) -- `search::AbstractSearchPolicy=LinearBinary()`: Search algorithm for interval finding +- `search::AbstractSearchPolicy=AutoSearch()`: Search algorithm for interval finding """ @inline @with_pool pool function cubic_interp!( output::AbstractVector{Tv}, @@ -35,7 +35,7 @@ Thread-safe: workspaces allocated from task-local pool. x_query::AbstractVector{Tg}; extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:AbstractFloat, Tv, X, F, BC} @assert length(y) == length(cache.x) "y length must match cache grid" @assert length(output) == length(x_query) "output length must match x_query" @@ -216,7 +216,7 @@ Pool-based exclusive extension: zero-alloc after warmup. end """ - cubic_interp!(output, x, y, x_query; bc=CubicFit(), extrap=NoExtrap(), autocache=true, deriv=EvalValue(), search=LinearBinary()) + cubic_interp!(output, x, y, x_query; bc=CubicFit(), extrap=NoExtrap(), autocache=true, deriv=EvalValue(), search=AutoSearch()) In-place cubic spline interpolation with optional automatic caching. """ @@ -229,7 +229,7 @@ In-place cubic spline interpolation with optional automatic caching. extrap::AbstractExtrap=NoExtrap(), autocache::Bool=true, deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:AbstractFloat, Tv} searcher = _to_searcher(search) # Periodic BC @@ -251,7 +251,7 @@ end x_query::Tg; extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:AbstractFloat, Tv, X, F, BC} @assert length(output) >= 1 "output must have at least 1 element" output[1] = cubic_interp_scalar(cache, y, x_query; extrap=extrap, deriv=deriv, search=search) @@ -267,7 +267,7 @@ end extrap::AbstractExtrap=NoExtrap(), autocache::Bool=true, deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:AbstractFloat, Tv} @assert length(output) >= 1 "output must have at least 1 element" output[1] = cubic_interp(x, y, x_query; bc, extrap, autocache, deriv, search) @@ -279,13 +279,13 @@ end # ======================================== """ - cubic_interp(cache, y, x_query; extrap=NoExtrap(), deriv=EvalValue(), search=LinearBinary()) -> Vector{T} + cubic_interp(cache, y, x_query; extrap=NoExtrap(), deriv=EvalValue(), search=AutoSearch()) -> Vector{T} Allocating version of cubic spline interpolation using cached LU factorization. # Arguments - `deriv::DerivOp=EvalValue()`: Derivative order (0=value, 1=first derivative, 2=second derivative) -- `search::AbstractSearchPolicy=LinearBinary()`: Search algorithm for interval finding +- `search::AbstractSearchPolicy=AutoSearch()`: Search algorithm for interval finding # Example ```julia @@ -305,7 +305,7 @@ function cubic_interp( x_query::AbstractVector{Tg}; extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:AbstractFloat, Tv} output = Vector{Tv}(undef, length(x_query)) cubic_interp!(output, cache, y, x_query; extrap=extrap, deriv=deriv, search=search) @@ -313,13 +313,13 @@ function cubic_interp( end """ - cubic_interp(x, y, x_query; bc=CubicFit(), extrap=NoExtrap(), autocache=true, deriv=EvalValue(), search=LinearBinary()) -> Vector{T} + cubic_interp(x, y, x_query; bc=CubicFit(), extrap=NoExtrap(), autocache=true, deriv=EvalValue(), search=AutoSearch()) -> Vector{T} Cubic spline interpolation with optional automatic caching. # Arguments - `deriv::DerivOp=EvalValue()`: Derivative order (0=value, 1=first derivative, 2=second derivative) -- `search::AbstractSearchPolicy=LinearBinary()`: Search algorithm for interval finding +- `search::AbstractSearchPolicy=AutoSearch()`: Search algorithm for interval finding # Extrapolation Modes - `NoExtrap()` (default): Throws DomainError if query point is outside domain @@ -347,7 +347,7 @@ function cubic_interp( extrap::AbstractExtrap=NoExtrap(), autocache::Bool=true, deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:AbstractFloat, Tv} output = Vector{Tv}(undef, length(x_query)) cubic_interp!(output, x, y, x_query; bc, extrap, autocache, deriv, search) @@ -356,7 +356,7 @@ end # Scalar query - zero allocation cubic_interp(cache::CubicSplineCache{Tg}, y::AbstractVector{Tv}, - x_query::Tg; extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), search=LinearBinary(), hint::Union{Nothing,Base.RefValue{Int}}=nothing) where {Tg<:AbstractFloat, Tv} = + x_query::Tg; extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), search=AutoSearch(), hint::Union{Nothing,Base.RefValue{Int}}=nothing) where {Tg<:AbstractFloat, Tv} = cubic_interp_scalar(cache, y, x_query; extrap=extrap, deriv=deriv, search=search, hint=hint) # Primary scalar method - AD-compatible @@ -369,7 +369,7 @@ function cubic_interp( extrap::AbstractExtrap=NoExtrap(), autocache::Bool=true, deriv::DerivOp=EvalValue(), - search=LinearBinary(), + search=AutoSearch(), hint::Union{Nothing,Base.RefValue{Int}}=nothing ) where {Tg<:AbstractFloat, Tv, Tq<:Real} searcher = _to_searcher(search, hint) @@ -398,7 +398,7 @@ function cubic_interp( extrap::AbstractExtrap=NoExtrap(), autocache::Bool=true, deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:Real, Tv, Tq<:Real} x_typed, y_typed, xq_typed = _promote_itp_inputs(x, y, x_query) return cubic_interp(x_typed, y_typed, xq_typed; bc, extrap, autocache, deriv, search) @@ -414,7 +414,7 @@ function cubic_interp( extrap::AbstractExtrap=NoExtrap(), autocache::Bool=true, deriv::DerivOp=EvalValue(), - search=LinearBinary(), + search=AutoSearch(), hint::Union{Nothing,Base.RefValue{Int}}=nothing ) where {Tg<:Real, Tv, Tq<:Real} x_typed, y_typed = _promote_itp_inputs(x, y) @@ -432,7 +432,7 @@ function cubic_interp!( extrap::AbstractExtrap=NoExtrap(), autocache::Bool=true, deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:Real, Tv, Tq<:Real} @assert length(y) == length(x) "x and y must have same length" @assert length(output) == length(x_query) "output must match x_query length" @@ -463,7 +463,7 @@ function cubic_interp!( extrap::AbstractExtrap=NoExtrap(), autocache::Bool=true, deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:Real, Tv, Tq<:Real} @assert length(output) >= 1 "output must have at least 1 element" diff --git a/src/cubic/cubic_series_interp.jl b/src/cubic/cubic_series_interp.jl index be5d2749..8086eca1 100644 --- a/src/cubic/cubic_series_interp.jl +++ b/src/cubic/cubic_series_interp.jl @@ -99,7 +99,7 @@ mutable struct CubicSeriesInterpolant{ y::Matrix{Tv}, z::Matrix{Tv}, extrap::E, - search::P=LinearBinary() + search::P=AutoSearch() ) where {Tg<:AbstractFloat, Tv, C<:CubicSplineCache{Tg}, B, E<:AbstractExtrap, P<:AbstractSearchPolicy} new{Tg, Tv, C, B, E, P}( cache, bc_for_solve, y, z, @@ -574,7 +574,7 @@ function cubic_interp( extrap::AbstractExtrap=NoExtrap(), autocache::Bool=true, precompute_transpose::Bool=false, - search::P=LinearBinary() + search::P=AutoSearch() ) where {Tg<:AbstractFloat, Tv, P<:AbstractSearchPolicy} # Validate input @assert !isempty(ys) "ys must not be empty" @@ -642,7 +642,7 @@ function _build_series_periodic( n_series_count::Int, autocache::Bool, precompute_transpose::Bool, - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:AbstractFloat, Tv} # Extend data for exclusive endpoint x, y_mat = _prepare_periodic(x, y_mat, bc) @@ -706,7 +706,7 @@ function cubic_interp( extrap::AbstractExtrap=NoExtrap(), autocache::Bool=true, precompute_transpose::Bool=false, - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:AbstractFloat, Tv} n_pts = length(x) @@ -764,7 +764,7 @@ function cubic_interp( extrap::AbstractExtrap=NoExtrap(), autocache::Bool=true, precompute_transpose::Bool=false, - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:Real, Tv} # Compute promoted grid type (Tg may be Int, promotes to Float) Tg_float = float(promote_type(Tg, _real_eltype(Tv))) @@ -782,7 +782,7 @@ function cubic_interp( extrap::AbstractExtrap=NoExtrap(), autocache::Bool=true, precompute_transpose::Bool=false, - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:Real, Tv} Tg_float = float(promote_type(Tg, _real_eltype(Tv))) Tv_float = _value_type(Tv, Tg_float) @@ -797,7 +797,7 @@ end # ======================================== """ - (sitp::CubicSeriesInterpolant)(xq::Real; deriv=EvalValue(), search=LinearBinary()) + (sitp::CubicSeriesInterpolant)(xq::Real; deriv=EvalValue(), search=AutoSearch()) Evaluate multi-Y interpolant at scalar query point (out-of-place). @@ -826,7 +826,7 @@ function (sitp::CubicSeriesInterpolant{Tg,Tv})( end """ - (sitp::CubicSeriesInterpolant)(output::AbstractVector, xq::Real; deriv=EvalValue(), search=LinearBinary()) + (sitp::CubicSeriesInterpolant)(output::AbstractVector, xq::Real; deriv=EvalValue(), search=AutoSearch()) Evaluate multi-Y interpolant at scalar query point (in-place). diff --git a/src/cubic/cubic_types.jl b/src/cubic/cubic_types.jl index a1cd36db..114888b8 100644 --- a/src/cubic/cubic_types.jl +++ b/src/cubic/cubic_types.jl @@ -141,7 +141,7 @@ struct CubicInterpolant{Tg<:AbstractFloat,Tv,C<:CubicSplineCache{Tg},E<:Abstract z::AbstractVector{Tv}, bc::BC, extrap::E, - search::P=LinearBinary() + search::P=AutoSearch() ) where {Tg<:AbstractFloat, Tv, C<:CubicSplineCache{Tg}, E<:AbstractExtrap, P<:AbstractSearchPolicy, BC<:CubicBC} @assert length(cache.x) == length(y) "cache grid and y must have same length" @assert length(cache.x) == length(z) "z coefficients must match grid length" diff --git a/src/cubic/nd/cubic_nd_interpolant.jl b/src/cubic/nd/cubic_nd_interpolant.jl index a408a2b4..53a05da5 100644 --- a/src/cubic/nd/cubic_nd_interpolant.jl +++ b/src/cubic/nd/cubic_nd_interpolant.jl @@ -25,7 +25,7 @@ Create an N-dimensional cubic Hermite interpolant from grid vectors and data arr - `extrap=NoExtrap()`: Extrapolation mode(s). Can be: - Single `AbstractExtrap`: Applied to all axes (`NoExtrap()`, `ConstExtrap()`, `WrapExtrap()`) - `NTuple{N,AbstractExtrap}`: Per-axis modes -- `search=LinearBinary()`: Search policy(s). Can be: +- `search=AutoSearch()`: Search policy(s). Can be: - Single `AbstractSearchPolicy`: Applied to all axes - `NTuple{N,AbstractSearchPolicy}`: Per-axis policies - `coeffs=PreCompute()`: Coefficient computation strategy @@ -62,7 +62,7 @@ function cubic_interp( data::AbstractArray{Tv_raw, N}; bc::Union{AbstractBC, NTuple{N,AbstractBC}}=CubicFit(), extrap::Union{AbstractExtrap, NTuple{N,AbstractExtrap}}=NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=LinearBinary(), + search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=AutoSearch(), coeffs::AbstractCoeffStrategy=PreCompute() ) where {N, Tv_raw} # Zero-allocation type promotion (uses @generated function) diff --git a/src/cubic/nd/cubic_nd_oneshot.jl b/src/cubic/nd/cubic_nd_oneshot.jl index e5da5f14..8ae5bd3e 100644 --- a/src/cubic/nd/cubic_nd_oneshot.jl +++ b/src/cubic/nd/cubic_nd_oneshot.jl @@ -32,7 +32,7 @@ function cubic_interp( deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue(), bc::Union{AbstractBC, NTuple{N,AbstractBC}}=CubicFit(), extrap::Union{AbstractExtrap, NTuple{N,AbstractExtrap}}=NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=LinearBinary(), + search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=AutoSearch(), coeffs::AbstractCoeffStrategy=PreCompute() ) where {Tv, N} # Type promotion + validation (same as constructor path) @@ -66,7 +66,7 @@ function cubic_interp( deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue(), bc::Union{AbstractBC, NTuple{N,AbstractBC}}=CubicFit(), extrap::Union{AbstractExtrap, NTuple{N,AbstractExtrap}}=NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=LinearBinary(), + search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=AutoSearch(), coeffs::AbstractCoeffStrategy=PreCompute() ) where {Tv, N} Tg = _promote_grid_eltype(grids) @@ -90,7 +90,7 @@ function cubic_interp( deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue(), bc::Union{AbstractBC, NTuple{N,AbstractBC}}=CubicFit(), extrap::Union{AbstractExtrap, NTuple{N,AbstractExtrap}}=NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=LinearBinary(), + search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=AutoSearch(), coeffs::AbstractCoeffStrategy=PreCompute() ) where {Tv, N} Tg = _promote_grid_eltype(grids) @@ -258,7 +258,7 @@ function cubic_interp!( deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue(), bc::Union{AbstractBC, NTuple{N,AbstractBC}}=CubicFit(), extrap::Union{AbstractExtrap, NTuple{N,AbstractExtrap}}=NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=LinearBinary(), + search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=AutoSearch(), coeffs::AbstractCoeffStrategy=PreCompute() ) where {Tv, N} Tg = _promote_grid_eltype(grids) @@ -290,7 +290,7 @@ function cubic_interp!( deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue(), bc::Union{AbstractBC, NTuple{N,AbstractBC}}=CubicFit(), extrap::Union{AbstractExtrap, NTuple{N,AbstractExtrap}}=NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=LinearBinary(), + search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=AutoSearch(), coeffs::AbstractCoeffStrategy=PreCompute() ) where {Tv, N} Tg = _promote_grid_eltype(grids) diff --git a/src/linear/linear_interpolant.jl b/src/linear/linear_interpolant.jl index 3524d41b..a6238bf8 100644 --- a/src/linear/linear_interpolant.jl +++ b/src/linear/linear_interpolant.jl @@ -61,7 +61,7 @@ end # ======================================== """ - linear_interp(x, y; extrap=NoExtrap(), search=LinearBinary()) -> LinearInterpolant + linear_interp(x, y; extrap=NoExtrap(), search=AutoSearch()) -> LinearInterpolant Create a callable interpolant for broadcast fusion and reuse. @@ -69,7 +69,7 @@ Create a callable interpolant for broadcast fusion and reuse. - `x::AbstractVector`: x-coordinates (must be sorted) - `y::AbstractVector`: y-values (can be real or complex) - `extrap::AbstractExtrap`: `NoExtrap()` (default), `ConstExtrap()`, `ExtendExtrap()`, or `WrapExtrap()` -- `search::AbstractSearchPolicy`: Default search policy for interval lookup (default: `LinearBinary()`) +- `search::AbstractSearchPolicy`: Default search policy for interval lookup (default: `AutoSearch()`) # Type Handling - x: Grid coordinates → converted to AbstractFloat, Range structure preserved @@ -91,7 +91,7 @@ Can be: # Examples ```julia -# Create with default LinearBinary() search policy +# Create with default AutoSearch() search policy itp = linear_interp(x_data, y_data) # Create with LinearBinary() as default policy (optimal for sorted queries) @@ -143,7 +143,7 @@ function linear_interp end x::AbstractVector{TX}, y::AbstractVector{TY}; extrap::AbstractExtrap=NoExtrap(), - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {TX<:Real, TY} x_p, y_p = _promote_itp_inputs(x, y) return LinearInterpolant(x_p, y_p; extrap, search) diff --git a/src/linear/linear_oneshot.jl b/src/linear/linear_oneshot.jl index f4150e9f..532d8b26 100644 --- a/src/linear/linear_oneshot.jl +++ b/src/linear/linear_oneshot.jl @@ -16,7 +16,7 @@ # ======================================== """ - linear_interp!(output, x, y, x_targets; extrap=NoExtrap(), deriv=EvalValue(), search=LinearBinary()) + linear_interp!(output, x, y, x_targets; extrap=NoExtrap(), deriv=EvalValue(), search=AutoSearch()) Zero-allocation linear interpolation with automatic dispatch: - For `AbstractRange` x: O(1) direct indexing @@ -62,7 +62,7 @@ function linear_interp!( x_targets::AbstractVector{Tg}; extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:AbstractFloat, Tv} @assert length(y) == length(x) "x and y must have same length" @assert length(output) == length(x_targets) "output must match x_targets length" @@ -155,7 +155,7 @@ end x_targets::AbstractVector{Tg}; extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:AbstractFloat, Tv} @assert length(y) == length(x) "x and y must have same length" @assert length(output) == length(x_targets) "output must match x_targets length" @@ -170,7 +170,7 @@ end # ======================================== """ - linear_interp(x, y, xq::Real; extrap=NoExtrap(), deriv=EvalValue(), search=LinearBinary()) -> AbstractFloat + linear_interp(x, y, xq::Real; extrap=NoExtrap(), deriv=EvalValue(), search=AutoSearch()) -> AbstractFloat Zero-allocation scalar linear interpolation with automatic dispatch: - For `AbstractRange` x: O(1) direct indexing @@ -385,7 +385,7 @@ end xq::Tq; extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), - search=LinearBinary(), + search=AutoSearch(), hint::Union{Nothing,Base.RefValue{Int}}=nothing ) where {Tg<:AbstractFloat, Tv, Tq<:Real} @boundscheck length(y) == length(x) || throw(ArgumentError("x and y must have same length")) @@ -412,7 +412,7 @@ function linear_interp( x_targets::AbstractVector{Tg}; extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:AbstractFloat, Tv} output = Vector{Tv}(undef, length(x_targets)) linear_interp!(output, x, y, x_targets; extrap, deriv, search) @@ -433,7 +433,7 @@ function linear_interp!( x_targets::AbstractVector{Tq}; extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:Real, Tv, Tq<:Real} @assert length(y) == length(x) "x and y must have same length" @assert length(output) == length(x_targets) "output must match x_targets length" @@ -465,7 +465,7 @@ end xq::Tq; extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), - search=LinearBinary(), + search=AutoSearch(), hint::Union{Nothing,Base.RefValue{Int}}=nothing ) where {Tg<:Real, Tv, Tq<:Real} x_typed, y_typed = _promote_itp_inputs(x, y) @@ -483,7 +483,7 @@ function linear_interp( x_targets::AbstractVector{Tq}; extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:Real, Tv, Tq<:Real} x_typed, y_typed, xq_typed = _promote_itp_inputs(x, y, x_targets) Tv_float = eltype(y_typed) diff --git a/src/linear/linear_series_interp.jl b/src/linear/linear_series_interp.jl index f42bc286..e78cfb93 100644 --- a/src/linear/linear_series_interp.jl +++ b/src/linear/linear_series_interp.jl @@ -78,7 +78,7 @@ mutable struct LinearSeriesInterpolant{Tg<:AbstractFloat, Tv, E<:AbstractExtrap, x::X, y::Matrix{Tv}, extrap::E, - search::P=LinearBinary() + search::P=AutoSearch() ) where {Tg<:AbstractFloat, Tv, E<:AbstractExtrap, P<:AbstractSearchPolicy, X<:AbstractVector{Tg}} new{Tg,Tv,E,P,X}(x, y, LazyTranspose{Tv}(), extrap, search) end @@ -321,7 +321,7 @@ function linear_interp( x::AbstractVector{Tg}, ys::AbstractVector{<:AbstractVector{Tv}}; extrap::AbstractExtrap=NoExtrap(), - search::P=LinearBinary() + search::P=AutoSearch() ) where {Tg<:AbstractFloat, Tv, P<:AbstractSearchPolicy} # Check if Tv's float base requires grid widening (not for Int types) # Int-based types (Complex{Int}) are handled by internal _value_type conversion @@ -385,7 +385,7 @@ function linear_interp( x::AbstractVector{Tg}, Y::AbstractMatrix{Tv}; extrap::AbstractExtrap=NoExtrap(), - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:AbstractFloat, Tv} # Check if Tv's float base requires grid widening Tv_real = _real_eltype(Tv) @@ -421,7 +421,7 @@ function linear_interp( x::AbstractVector{Tg}, ys::AbstractVector{<:AbstractVector{Tv}}; extrap::AbstractExtrap=NoExtrap(), - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:Real, Tv} # Compute promoted grid type (Tg may be Int, promotes to Float) Tg_float = float(promote_type(Tg, _real_eltype(Tv))) @@ -434,7 +434,7 @@ function linear_interp( x::AbstractVector{Tg}, Y::AbstractMatrix{Tv}; extrap::AbstractExtrap=NoExtrap(), - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:Real, Tv} Tg_float = float(promote_type(Tg, _real_eltype(Tv))) x_typed = _to_float(x, Tg_float) @@ -446,7 +446,7 @@ end # ======================================== """ - (sitp::LinearSeriesInterpolant)(xq::Real; deriv=EvalValue(), search=LinearBinary()) + (sitp::LinearSeriesInterpolant)(xq::Real; deriv=EvalValue(), search=AutoSearch()) Evaluate multi-Y interpolant at scalar query point (out-of-place). @@ -465,7 +465,7 @@ function (sitp::LinearSeriesInterpolant{Tg,Tv,P})( end """ - (sitp::LinearSeriesInterpolant)(output::AbstractVector, xq::Real; deriv=EvalValue(), search=LinearBinary()) + (sitp::LinearSeriesInterpolant)(output::AbstractVector, xq::Real; deriv=EvalValue(), search=AutoSearch()) Evaluate multi-Y interpolant at scalar query point (in-place). diff --git a/src/linear/linear_types.jl b/src/linear/linear_types.jl index 3b298a6c..84064973 100644 --- a/src/linear/linear_types.jl +++ b/src/linear/linear_types.jl @@ -27,10 +27,10 @@ Returned by `linear_interp(x, y)` (2-argument form). # Usage ```julia # Create interpolator (minimal allocation) -itp = linear_interp(x, y) # default extrap=NoExtrap(), search=LinearBinary() +itp = linear_interp(x, y) # default extrap=NoExtrap(), search=AutoSearch() # Create with custom search policy (baked-in default) -itp = linear_interp(x, y; search=LinearBinary()) +itp = linear_interp(x, y; search=AutoSearch()) # Complex-valued interpolation y_complex = exp.(2im .* x) @@ -85,7 +85,7 @@ end x::X, y::Y; extrap::AbstractExtrap=NoExtrap(), - search::P=LinearBinary() + search::P=AutoSearch() ) where {Tg<:AbstractFloat, Tv, X<:AbstractVector{Tg}, Y<:AbstractVector{Tv}, P<:AbstractSearchPolicy} E = typeof(extrap) return LinearInterpolant{Tg,Tv,X,Y,E,P}(x, y, extrap, search) diff --git a/src/linear/nd/linear_nd_interpolant.jl b/src/linear/nd/linear_nd_interpolant.jl index 56bf2982..54185d84 100644 --- a/src/linear/nd/linear_nd_interpolant.jl +++ b/src/linear/nd/linear_nd_interpolant.jl @@ -23,7 +23,7 @@ The interpolation is exact at grid points and linearly blended between them. # Keyword Arguments - `extrap=NoExtrap()`: Extrapolation mode (`NoExtrap()`, `ConstExtrap()`, `ExtendExtrap()`, `WrapExtrap()`) or per-axis tuple -- `search=LinearBinary()`: Search policy or per-axis tuple +- `search=AutoSearch()`: Search policy or per-axis tuple # Returns - `LinearInterpolantND{Tg,Tv,N}`: Callable interpolant object @@ -62,7 +62,7 @@ function linear_interp( grids::NTuple{N, AbstractVector}, data::AbstractArray{Tv_raw, N}; extrap::Union{AbstractExtrap, NTuple{N, AbstractExtrap}} = NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = LinearBinary() + search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = AutoSearch() ) where {N, Tv_raw} # Validate grid dimensions _validate_nd_grids(grids, data) diff --git a/src/linear/nd/linear_nd_oneshot.jl b/src/linear/nd/linear_nd_oneshot.jl index 7341f02e..459b6390 100644 --- a/src/linear/nd/linear_nd_oneshot.jl +++ b/src/linear/nd/linear_nd_oneshot.jl @@ -112,7 +112,7 @@ function linear_interp( data::AbstractArray{Tv, N}, query::Tuple{Vararg{Real, N}}; extrap::Union{AbstractExtrap, NTuple{N, AbstractExtrap}} = NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = LinearBinary(), + search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = AutoSearch(), deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue() ) where {Tv, N} Tg = _promote_grid_eltype(grids) @@ -139,7 +139,7 @@ function linear_interp( data::AbstractArray{Tv, N}, queries::NTuple{N, AbstractVector{<:Real}}; extrap::Union{AbstractExtrap, NTuple{N, AbstractExtrap}} = NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = LinearBinary(), + search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = AutoSearch(), deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue() ) where {Tv, N} Tg = _promote_grid_eltype(grids) @@ -161,7 +161,7 @@ function linear_interp( data::AbstractArray{Tv, N}, queries::AbstractVector{<:Tuple{Vararg{Real, N}}}; extrap::Union{AbstractExtrap, NTuple{N, AbstractExtrap}} = NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = LinearBinary(), + search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = AutoSearch(), deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue() ) where {Tv, N} Tg = _promote_grid_eltype(grids) @@ -188,7 +188,7 @@ function linear_interp!( data::AbstractArray{Tv, N}, queries::NTuple{N, AbstractVector{<:Real}}; extrap::Union{AbstractExtrap, NTuple{N, AbstractExtrap}} = NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = LinearBinary(), + search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = AutoSearch(), deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue() ) where {Tv, N} Tg = _promote_grid_eltype(grids) @@ -215,7 +215,7 @@ function linear_interp!( data::AbstractArray{Tv, N}, queries::AbstractVector{<:Tuple{Vararg{Real, N}}}; extrap::Union{AbstractExtrap, NTuple{N, AbstractExtrap}} = NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = LinearBinary(), + search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = AutoSearch(), deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue() ) where {Tv, N} Tg = _promote_grid_eltype(grids) diff --git a/src/quadratic/nd/quadratic_nd_interpolant.jl b/src/quadratic/nd/quadratic_nd_interpolant.jl index e88724a2..43d645d3 100644 --- a/src/quadratic/nd/quadratic_nd_interpolant.jl +++ b/src/quadratic/nd/quadratic_nd_interpolant.jl @@ -82,7 +82,7 @@ Create an N-dimensional quadratic interpolant from grid vectors and data array. - Single BC: Applied to all axes - `NTuple{N}`: Per-axis BCs - `extrap=NoExtrap()`: Extrapolation mode(s) -- `search=LinearBinary()`: Search policy(s) +- `search=AutoSearch()`: Search policy(s) # Returns - `QuadraticInterpolantND{Tg, Tv, N, ...}`: Callable interpolant object @@ -102,7 +102,7 @@ function quadratic_interp( data::AbstractArray{Tv_raw, N}; bc::Union{AbstractBC, NTuple{N,AbstractBC}}=Left(QuadraticFit()), extrap::Union{AbstractExtrap, NTuple{N,AbstractExtrap}}=NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=LinearBinary() + search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=AutoSearch() ) where {N, Tv_raw} # Zero-allocation type promotion Tg = _promote_grid_eltype(grids) diff --git a/src/quadratic/nd/quadratic_nd_oneshot.jl b/src/quadratic/nd/quadratic_nd_oneshot.jl index 67acf973..f1a0d7ab 100644 --- a/src/quadratic/nd/quadratic_nd_oneshot.jl +++ b/src/quadratic/nd/quadratic_nd_oneshot.jl @@ -144,7 +144,7 @@ function quadratic_interp( deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue(), bc::Union{AbstractBC, NTuple{N,AbstractBC}}=Left(QuadraticFit()), extrap::Union{AbstractExtrap, NTuple{N,AbstractExtrap}}=NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=LinearBinary() + search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=AutoSearch() ) where {Tv, N} Tg = _promote_grid_eltype(grids) Tg = Tg <: AbstractFloat ? Tg : Float64 @@ -174,7 +174,7 @@ function quadratic_interp( deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue(), bc::Union{AbstractBC, NTuple{N,AbstractBC}}=Left(QuadraticFit()), extrap::Union{AbstractExtrap, NTuple{N,AbstractExtrap}}=NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=LinearBinary() + search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=AutoSearch() ) where {Tv, N} Tg = _promote_grid_eltype(grids) Tg = Tg <: AbstractFloat ? Tg : Float64 @@ -197,7 +197,7 @@ function quadratic_interp( deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue(), bc::Union{AbstractBC, NTuple{N,AbstractBC}}=Left(QuadraticFit()), extrap::Union{AbstractExtrap, NTuple{N,AbstractExtrap}}=NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=LinearBinary() + search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=AutoSearch() ) where {Tv, N} Tg = _promote_grid_eltype(grids) Tg = Tg <: AbstractFloat ? Tg : Float64 @@ -225,7 +225,7 @@ function quadratic_interp!( deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue(), bc::Union{AbstractBC, NTuple{N,AbstractBC}}=Left(QuadraticFit()), extrap::Union{AbstractExtrap, NTuple{N,AbstractExtrap}}=NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=LinearBinary() + search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=AutoSearch() ) where {Tv, N} Tg = _promote_grid_eltype(grids) Tg = Tg <: AbstractFloat ? Tg : Float64 @@ -255,7 +255,7 @@ function quadratic_interp!( deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue(), bc::Union{AbstractBC, NTuple{N,AbstractBC}}=Left(QuadraticFit()), extrap::Union{AbstractExtrap, NTuple{N,AbstractExtrap}}=NoExtrap(), - search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=LinearBinary() + search::Union{AbstractSearchPolicy, NTuple{N,AbstractSearchPolicy}}=AutoSearch() ) where {Tv, N} Tg = _promote_grid_eltype(grids) Tg = Tg <: AbstractFloat ? Tg : Float64 diff --git a/src/quadratic/quadratic_interpolant.jl b/src/quadratic/quadratic_interpolant.jl index 1a003eec..9126772a 100644 --- a/src/quadratic/quadratic_interpolant.jl +++ b/src/quadratic/quadratic_interpolant.jl @@ -55,7 +55,7 @@ end # ======================================== """ - quadratic_interp(x, y; bc=Left(QuadraticFit()), extrap=NoExtrap(), search=LinearBinary()) -> QuadraticInterpolant + quadratic_interp(x, y; bc=Left(QuadraticFit()), extrap=NoExtrap(), search=AutoSearch()) -> QuadraticInterpolant Create a callable interpolant for broadcast fusion and reuse. @@ -64,7 +64,7 @@ Create a callable interpolant for broadcast fusion and reuse. - `y::AbstractVector`: y-values (can be Real or Complex) - `bc`: Boundary condition (Left, Right, MinCurvFit, or Left/Right with QuadraticFit) - `extrap::AbstractExtrap`: `NoExtrap()` (default), `ConstExtrap()`, `ExtendExtrap()`, or `WrapExtrap()` -- `search::AbstractSearchPolicy`: Default search policy (default: `LinearBinary()`) +- `search::AbstractSearchPolicy`: Default search policy (default: `AutoSearch()`) # Returns `QuadraticInterpolant{Tg, Tv}` object for scalar/broadcast evaluation. @@ -115,7 +115,7 @@ end y::AbstractVector{TY}; bc::QuadraticBC=Left(QuadraticFit()), extrap::AbstractExtrap=NoExtrap(), - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {TX<:Real, TY} x_p, y_p = _promote_itp_inputs(x, y) bc_p = _promote_bc(bc, eltype(x_p)) diff --git a/src/quadratic/quadratic_oneshot.jl b/src/quadratic/quadratic_oneshot.jl index 3b81d318..a49392b5 100644 --- a/src/quadratic/quadratic_oneshot.jl +++ b/src/quadratic/quadratic_oneshot.jl @@ -149,7 +149,7 @@ end # ======================================== """ - quadratic_interp(x, y, xi; bc=Left(QuadraticFit()), extrap=NoExtrap(), deriv=EvalValue(), search=LinearBinary()) + quadratic_interp(x, y, xi; bc=Left(QuadraticFit()), extrap=NoExtrap(), deriv=EvalValue(), search=AutoSearch()) C1 piecewise quadratic spline interpolation at a single point. @@ -202,7 +202,7 @@ vals = quadratic_interp(x, y, sorted_queries; search=LinearBinary(linear_window= bc::QuadraticBC=Left(QuadraticFit()), extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), - search=LinearBinary(), + search=AutoSearch(), hint::Union{Nothing,Base.RefValue{Int}}=nothing ) where {Tg<:AbstractFloat, Tv, Tq<:Real} @boundscheck length(y) == length(x) || throw(ArgumentError("x and y must have same length")) @@ -226,7 +226,7 @@ end # ======================================== """ - quadratic_interp!(output, x, y, x_targets; bc=Left(QuadraticFit()), extrap=NoExtrap(), deriv=EvalValue(), search=LinearBinary()) + quadratic_interp!(output, x, y, x_targets; bc=Left(QuadraticFit()), extrap=NoExtrap(), deriv=EvalValue(), search=AutoSearch()) In-place quadratic spline interpolation for multiple query points. @@ -256,7 +256,7 @@ quadratic_interp!(output, x, y, sorted_queries; search=LinearBinary(linear_windo bc::QuadraticBC=Left(QuadraticFit()), extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:AbstractFloat, Tv} @assert length(y) == length(x) "x and y must have same length" @assert length(output) == length(x_targets) "output must match x_targets length" @@ -284,7 +284,7 @@ end # ======================================== """ - quadratic_interp(x, y, x_targets; bc=Left(QuadraticFit()), extrap=NoExtrap(), deriv=EvalValue(), search=LinearBinary()) + quadratic_interp(x, y, x_targets; bc=Left(QuadraticFit()), extrap=NoExtrap(), deriv=EvalValue(), search=AutoSearch()) Quadratic spline interpolation for multiple query points (allocating version). @@ -307,7 +307,7 @@ function quadratic_interp( bc::QuadraticBC=Left(QuadraticFit()), extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:AbstractFloat, Tv} output = Vector{Tv}(undef, length(x_targets)) quadratic_interp!(output, x, y, x_targets; bc, extrap, deriv, search) @@ -362,7 +362,7 @@ end bc::QuadraticBC=Left(QuadraticFit()), extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), - search=LinearBinary(), + search=AutoSearch(), hint::Union{Nothing,Base.RefValue{Int}}=nothing ) where {Tg<:Real, Tv, Tq<:Real} x_typed, y_typed = _promote_itp_inputs(x, y) @@ -384,7 +384,7 @@ function quadratic_interp( bc::QuadraticBC=Left(QuadraticFit()), extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:Real, Tv, Tq<:Real} x_typed, y_typed, xq_typed = _promote_itp_inputs(x, y, x_targets) Tv_float = eltype(y_typed) @@ -406,7 +406,7 @@ function quadratic_interp!( bc::QuadraticBC=Left(QuadraticFit()), extrap::AbstractExtrap=NoExtrap(), deriv::DerivOp=EvalValue(), - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:Real, Tv, Tq<:Real} @assert length(y) == length(x) "x and y must have same length" @assert length(output) == length(x_targets) "output must match x_targets length" diff --git a/src/quadratic/quadratic_series_interp.jl b/src/quadratic/quadratic_series_interp.jl index aec34a81..26f84088 100644 --- a/src/quadratic/quadratic_series_interp.jl +++ b/src/quadratic/quadratic_series_interp.jl @@ -93,7 +93,7 @@ mutable struct QuadraticSeriesInterpolant{Tg<:AbstractFloat, Tv, E<:AbstractExtr d::Matrix{Tv}, h::Vector{Tg}, extrap::E, - search::P=LinearBinary() + search::P=AutoSearch() ) where {Tg<:AbstractFloat, Tv, E<:AbstractExtrap, P<:AbstractSearchPolicy, X<:AbstractVector{Tg}} new{Tg,Tv,E,P,X}(x, y, a, d, h, LazyTransposeTriple{Tv}(), extrap, search) end @@ -367,7 +367,7 @@ function quadratic_interp( ys::AbstractVector{<:AbstractVector{Tv}}; bc::QuadraticBC=Left(QuadraticFit()), extrap::AbstractExtrap=NoExtrap(), - search::P=LinearBinary() + search::P=AutoSearch() ) where {Tg<:AbstractFloat, Tv, P<:AbstractSearchPolicy} # Check if Tv's float base requires grid widening (not for Int types) # Int-based types (Complex{Int}) are handled by internal _value_type conversion @@ -449,7 +449,7 @@ function quadratic_interp( Y::AbstractMatrix{Tv}; bc::QuadraticBC=Left(QuadraticFit()), extrap::AbstractExtrap=NoExtrap(), - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:AbstractFloat, Tv} ys = [Y[:, k] for k in axes(Y, 2)] return quadratic_interp(x, ys; bc=bc, extrap=extrap, search=search) @@ -466,7 +466,7 @@ function quadratic_interp( ys::AbstractVector{<:AbstractVector{Tv}}; bc=Left(QuadraticFit()), extrap::AbstractExtrap=NoExtrap(), - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:Real, Tv} # Compute promoted grid type (Tg may be Int, promotes to Float) Tg_float = float(promote_type(Tg, _real_eltype(Tv))) @@ -481,7 +481,7 @@ function quadratic_interp( Y::AbstractMatrix{Tv}; bc=Left(QuadraticFit()), extrap::AbstractExtrap=NoExtrap(), - search::AbstractSearchPolicy=LinearBinary() + search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:Real, Tv} Tg_float = float(promote_type(Tg, _real_eltype(Tv))) x_typed = _to_float(x, Tg_float) @@ -495,7 +495,7 @@ end # Scalar evaluation (explicit implementation for deriv keyword support) """ - (sitp::QuadraticSeriesInterpolant)(xq::Real; deriv=EvalValue(), search=LinearBinary()) + (sitp::QuadraticSeriesInterpolant)(xq::Real; deriv=EvalValue(), search=AutoSearch()) Evaluate all series at scalar query point (out-of-place). @@ -520,7 +520,7 @@ function (sitp::QuadraticSeriesInterpolant{Tg,Tv,P})( end """ - (sitp::QuadraticSeriesInterpolant)(output::AbstractVector, xq::Real; deriv=EvalValue(), search=LinearBinary()) + (sitp::QuadraticSeriesInterpolant)(output::AbstractVector, xq::Real; deriv=EvalValue(), search=AutoSearch()) Evaluate all series at scalar query point (in-place, zero allocation). diff --git a/src/quadratic/quadratic_types.jl b/src/quadratic/quadratic_types.jl index 074d3741..4aa397f2 100644 --- a/src/quadratic/quadratic_types.jl +++ b/src/quadratic/quadratic_types.jl @@ -88,7 +88,7 @@ end a::Vector{Tv}, d::Vector{Tv}; extrap::AbstractExtrap=NoExtrap(), - search::P=LinearBinary() + search::P=AutoSearch() ) where {Tg<:AbstractFloat, Tv, X<:AbstractVector{Tg}, Y<:AbstractVector{Tv}, P<:AbstractSearchPolicy} E = typeof(extrap) return QuadraticInterpolant{Tg,Tv,X,Y,E,P}(x, y, h, a, d, extrap, search) From b7640c0ad8de9554fc22147d3b1ea7b97d9b9468 Mon Sep 17 00:00:00 2001 From: Min-Gu Yoo Date: Tue, 24 Feb 2026 12:51:17 -0800 Subject: [PATCH 07/16] Insert _resolve_search at all 1D eval sites Adds AutoSearch resolution before _to_searcher in all interpolant callables, oneshot hot paths, and series eval methods. Scalar queries now resolve to Binary(), vector queries to LinearBinary(). --- src/constant/constant_interpolant.jl | 9 ++++++--- src/constant/constant_oneshot.jl | 6 ++++-- src/constant/constant_series_interp.jl | 6 ++++-- src/cubic/cubic_eval.jl | 3 ++- src/cubic/cubic_interpolant.jl | 9 ++++++--- src/cubic/cubic_oneshot.jl | 9 ++++++--- src/cubic/cubic_series_interp.jl | 9 ++++++--- src/linear/linear_interpolant.jl | 9 ++++++--- src/linear/linear_oneshot.jl | 9 ++++++--- src/linear/linear_series_interp.jl | 6 ++++-- src/quadratic/quadratic_interpolant.jl | 9 ++++++--- src/quadratic/quadratic_oneshot.jl | 6 ++++-- src/quadratic/quadratic_series_interp.jl | 11 +++++++---- 13 files changed, 67 insertions(+), 34 deletions(-) diff --git a/src/constant/constant_interpolant.jl b/src/constant/constant_interpolant.jl index 1ea1f1e0..627a07c2 100644 --- a/src/constant/constant_interpolant.jl +++ b/src/constant/constant_interpolant.jl @@ -12,7 +12,8 @@ # Type parameters: Tg = grid type, Tv = value type, Tq = query type # ───────────────────────────────────────────────────────────── @inline function (itp::ConstantInterpolant{Tg,Tv})(xq::Tq; deriv::DerivOp=EvalValue(), search=itp.search_policy, hint::Union{Nothing,Base.RefValue{Int}}=nothing) where {Tg<:AbstractFloat, Tv, Tq<:Real} - searcher = _to_searcher(search, hint) + resolved = _resolve_search(search, xq) + searcher = _to_searcher(resolved, hint) _constant_eval_at_point(itp.x, itp.y, xq, itp.extrap, itp.side, deriv, searcher) end @@ -24,7 +25,8 @@ end function (itp::ConstantInterpolant{Tg,Tv})(xq::AbstractVector{Tq}; deriv::DerivOp=EvalValue(), search=itp.search_policy, hint::Union{Nothing,Base.RefValue{Int}}=nothing) where {Tg<:AbstractFloat, Tv, Tq<:Real} T_out = promote_type(Tv, Tq) # Lossless: wider type to avoid precision loss output = Vector{T_out}(undef, length(xq)) - searcher = _to_searcher(search, hint) + resolved = _resolve_search(search, xq) + searcher = _to_searcher(resolved, hint) @boundscheck _check_domain(itp.x, xq, itp.extrap) @inbounds for i in eachindex(xq, output) output[i] = _constant_eval_at_point(itp.x, itp.y, xq[i], itp.extrap, itp.side, deriv, searcher) @@ -37,7 +39,8 @@ end # ───────────────────────────────────────────────────────────── function (itp::ConstantInterpolant{Tg,Tv})(output::AbstractVector, xq::AbstractVector{Tq}; deriv::DerivOp=EvalValue(), search=itp.search_policy, hint::Union{Nothing,Base.RefValue{Int}}=nothing) where {Tg<:AbstractFloat, Tv, Tq<:Real} @assert length(output) == length(xq) "output length must match xq length" - searcher = _to_searcher(search, hint) + resolved = _resolve_search(search, xq) + searcher = _to_searcher(resolved, hint) @boundscheck _check_domain(itp.x, xq, itp.extrap) @inbounds for i in eachindex(xq, output) output[i] = _constant_eval_at_point(itp.x, itp.y, xq[i], itp.extrap, itp.side, deriv, searcher) diff --git a/src/constant/constant_oneshot.jl b/src/constant/constant_oneshot.jl index 9ea451de..b7fcd46b 100644 --- a/src/constant/constant_oneshot.jl +++ b/src/constant/constant_oneshot.jl @@ -205,7 +205,8 @@ vals = constant_interp(x, y, sorted_queries; search=LinearBinary(linear_window=8 ) where {Tg<:AbstractFloat, Tv, Tq<:Real} @boundscheck length(y) == length(x) || throw(ArgumentError("x and y must have same length")) - searcher = _to_searcher(search, hint) + resolved = _resolve_search(search, xi) + searcher = _to_searcher(resolved, hint) _constant_eval_at_point(x, y, xi, extrap, side, deriv, searcher) end @@ -251,7 +252,8 @@ function constant_interp!( @assert length(y) == length(x) "x and y must have same length" @assert length(output) == length(x_targets) "output must match x_targets length" - searcher = _to_searcher(search) + resolved = _resolve_search(search, x_targets) + searcher = _to_searcher(resolved) @boundscheck _check_domain(x, x_targets, extrap) @inbounds for i in eachindex(x_targets, output) output[i] = _constant_eval_at_point(x, y, x_targets[i], extrap, side, deriv, searcher) diff --git a/src/constant/constant_series_interp.jl b/src/constant/constant_series_interp.jl index e765e08c..0014d13e 100644 --- a/src/constant/constant_series_interp.jl +++ b/src/constant/constant_series_interp.jl @@ -504,7 +504,8 @@ function (sitp::ConstantSeriesInterpolant{Tg,Tv,P})( xq_typed = Tg(xq_primal) # Build anchor using primal value - aq = _make_anchor(sitp, xq_typed, _to_searcher(search, hint)) + resolved = _resolve_search(search, xq) + aq = _make_anchor(sitp, xq_typed, _to_searcher(resolved, hint)) # Dispatch on derivative order - pass original xq for AD support _eval_constant_series_point!(output, sitp, aq, xq, deriv) @@ -569,8 +570,9 @@ Uses task-local pool for anchor vector to achieve zero allocation after warmup. _validate_series_outputs(outputs, n_ser, n_query) # Build anchors from pool (zero allocation after warmup) + resolved = _resolve_search(search, xq) aq_vec = acquire!(pool, _ConstantAnchoredQuery{Tg}, length(xq)) - _fill_anchors!(aq_vec, sitp.x, xq, Val(:constant); wrap=_should_wrap(sitp), searcher=_to_searcher(search, hint)) + _fill_anchors!(aq_vec, sitp.x, xq, Val(:constant); wrap=_should_wrap(sitp), searcher=_to_searcher(resolved, hint)) # Extract matrices for argument-passing pattern y = sitp.y diff --git a/src/cubic/cubic_eval.jl b/src/cubic/cubic_eval.jl index 08d95f51..f8c5582d 100644 --- a/src/cubic/cubic_eval.jl +++ b/src/cubic/cubic_eval.jl @@ -281,7 +281,8 @@ Uses task-local pool for workspace allocation. z = similar!(pool, y) _solve_system!(z, cache, y, cache.bc_config) - searcher = _to_searcher(search, hint) + resolved = _resolve_search(search, x_query) + searcher = _to_searcher(resolved, hint) @boundscheck _check_domain(cache.x, x_query, extrap) _eval_with_bc(cache, y, z, x_query, extrap, deriv, searcher) end diff --git a/src/cubic/cubic_interpolant.jl b/src/cubic/cubic_interpolant.jl index 6b2072b8..8c9dae7e 100644 --- a/src/cubic/cubic_interpolant.jl +++ b/src/cubic/cubic_interpolant.jl @@ -27,7 +27,8 @@ # Unified method: accepts any query type (Tg, Real, or Dual for AD) @inline function (itp::CubicInterpolant{Tg,Tv})(xq; deriv::DerivOp=EvalValue(), search=itp.search_policy, hint::Union{Nothing,Base.RefValue{Int}}=nothing) where {Tg<:AbstractFloat, Tv} @boundscheck _check_domain(itp.cache.x, xq, itp.extrap) - searcher = _to_searcher(search, hint) + resolved = _resolve_search(search, xq) + searcher = _to_searcher(resolved, hint) # Pass original xq to preserve Dual type for AD _eval_with_bc(itp.cache, itp.y, itp.z, xq, itp.extrap, deriv, searcher) end @@ -38,7 +39,8 @@ end function (itp::CubicInterpolant{Tg,Tv})(xq::AbstractVector{Tq}; deriv::DerivOp=EvalValue(), search=itp.search_policy, hint::Union{Nothing,Base.RefValue{Int}}=nothing) where {Tg<:AbstractFloat, Tv, Tq<:Real} T_out = promote_type(Tv, Tq) # Lossless: wider type to avoid precision loss output = Vector{T_out}(undef, length(xq)) - searcher = _to_searcher(search, hint) + resolved = _resolve_search(search, xq) + searcher = _to_searcher(resolved, hint) _cubic_vector_loop!(output, itp.cache, itp.y, itp.z, xq, itp.extrap, deriv, searcher) return output end @@ -46,7 +48,8 @@ end # In-place vector call with deriv keyword support function (itp::CubicInterpolant{Tg,Tv})(output::AbstractVector, xq::AbstractVector{Tq}; deriv::DerivOp=EvalValue(), search=itp.search_policy, hint::Union{Nothing,Base.RefValue{Int}}=nothing) where {Tg<:AbstractFloat, Tv, Tq<:Real} @assert length(output) == length(xq) "output length must match xq length" - searcher = _to_searcher(search, hint) + resolved = _resolve_search(search, xq) + searcher = _to_searcher(resolved, hint) _cubic_vector_loop!(output, itp.cache, itp.y, itp.z, xq, itp.extrap, deriv, searcher) return output end diff --git a/src/cubic/cubic_oneshot.jl b/src/cubic/cubic_oneshot.jl index a37fe206..84f8a6cb 100644 --- a/src/cubic/cubic_oneshot.jl +++ b/src/cubic/cubic_oneshot.jl @@ -43,7 +43,8 @@ Thread-safe: workspaces allocated from task-local pool. z = similar!(pool, y) _solve_system!(z, cache, y, cache.bc_config) - searcher = _to_searcher(search) + resolved = _resolve_search(search, x_query) + searcher = _to_searcher(resolved) _cubic_vector_loop!(output, cache, y, z, x_query, extrap, deriv, searcher) return output @@ -231,7 +232,8 @@ In-place cubic spline interpolation with optional automatic caching. deriv::DerivOp=EvalValue(), search::AbstractSearchPolicy=AutoSearch() ) where {Tg<:AbstractFloat, Tv} - searcher = _to_searcher(search) + resolved = _resolve_search(search, x_query) + searcher = _to_searcher(resolved) # Periodic BC if _is_periodic_bc(bc) return _cubic_interp_periodic!(output, x, y, x_query, bc, autocache, deriv, searcher) @@ -372,7 +374,8 @@ function cubic_interp( search=AutoSearch(), hint::Union{Nothing,Base.RefValue{Int}}=nothing ) where {Tg<:AbstractFloat, Tv, Tq<:Real} - searcher = _to_searcher(search, hint) + resolved = _resolve_search(search, xq) + searcher = _to_searcher(resolved, hint) if _is_periodic_bc(bc) return _cubic_interp_periodic_scalar(x, y, xq, bc, autocache, deriv, searcher) end diff --git a/src/cubic/cubic_series_interp.jl b/src/cubic/cubic_series_interp.jl index 8086eca1..1606c766 100644 --- a/src/cubic/cubic_series_interp.jl +++ b/src/cubic/cubic_series_interp.jl @@ -819,7 +819,8 @@ function (sitp::CubicSeriesInterpolant{Tg,Tv})( output = Vector{T_out}(undef, n_series(sitp)) # Build anchor preserving Dual type in xq - aq = _make_anchor(sitp, xq_promoted, _to_searcher(search, hint)) + resolved = _resolve_search(search, xq) + aq = _make_anchor(sitp, xq_promoted, _to_searcher(resolved, hint)) _eval_series_at_anchor!(output, sitp, aq, deriv) return output @@ -846,7 +847,8 @@ function (sitp::CubicSeriesInterpolant{Tg,Tv})( xq_promoted = _promote_for_anchor(xq, Tg) # Build anchor preserving Dual type in xq (for AD) - aq = _make_anchor(sitp, xq_promoted, _to_searcher(search, hint)) + resolved = _resolve_search(search, xq) + aq = _make_anchor(sitp, xq_promoted, _to_searcher(resolved, hint)) _eval_series_at_anchor!(output, sitp, aq, deriv) return output @@ -931,8 +933,9 @@ Builds anchors from original `xq` (preserving precision in weights) for scalar/v # Build anchors - pool handles both Tq===Tg and mixed-type cases # Each unique type combination gets its own pool slot + resolved = _resolve_search(search, xq) aq_vec = acquire!(pool, _CubicAnchoredQuery{Tg,Tq}, n_query) - _fill_anchors!(aq_vec, sitp.cache.x, xq, Val(:cubic); wrap=_should_wrap(sitp), searcher=_to_searcher(search, hint)) + _fill_anchors!(aq_vec, sitp.cache.x, xq, Val(:cubic); wrap=_should_wrap(sitp), searcher=_to_searcher(resolved, hint)) # Extract matrices for argument-passing pattern (series-contiguous layout) # This is faster than point-contiguous for vector queries because: diff --git a/src/linear/linear_interpolant.jl b/src/linear/linear_interpolant.jl index a6238bf8..2b9ff7c5 100644 --- a/src/linear/linear_interpolant.jl +++ b/src/linear/linear_interpolant.jl @@ -23,7 +23,8 @@ # - ForwardDiff.Dual queries (automatic differentiation) @inline function (itp::LinearInterpolant{Tg,Tv,X,Y,E,P})(xq; deriv::DerivOp=EvalValue(), search=itp.search_policy, hint::Union{Nothing,Base.RefValue{Int}}=nothing) where {Tg<:AbstractFloat, Tv, X, Y, E, P} @boundscheck _check_domain(itp.x, xq, itp.extrap) - searcher = _to_searcher(search, hint) + resolved = _resolve_search(search, xq) + searcher = _to_searcher(resolved, hint) # Pass original xq to preserve Dual type for AD _linear_with_extrap(itp.x, itp.y, xq, itp.extrap, deriv, searcher) end @@ -36,7 +37,8 @@ function (itp::LinearInterpolant{Tg,Tv,X,Y,E,P})(xq::AbstractVector{Tq}; deriv:: @boundscheck _check_domain(itp.x, xq, itp.extrap) T_out = promote_type(Tv, Tq) # Lossless: wider type to avoid precision loss output = Vector{T_out}(undef, length(xq)) - searcher = _to_searcher(search, hint) + resolved = _resolve_search(search, xq) + searcher = _to_searcher(resolved, hint) @inbounds for i in eachindex(xq, output) output[i] = _linear_with_extrap(itp.x, itp.y, xq[i], itp.extrap, deriv, searcher) end @@ -49,7 +51,8 @@ end function (itp::LinearInterpolant{Tg,Tv,X,Y,E,P})(output::AbstractVector, xq::AbstractVector{Tq}; deriv::DerivOp=EvalValue(), search=itp.search_policy, hint::Union{Nothing,Base.RefValue{Int}}=nothing) where {Tg<:AbstractFloat, Tv, X, Y, E, P, Tq<:Real} @assert length(output) == length(xq) "output length must match xq length" @boundscheck _check_domain(itp.x, xq, itp.extrap) - searcher = _to_searcher(search, hint) + resolved = _resolve_search(search, xq) + searcher = _to_searcher(resolved, hint) @inbounds for i in eachindex(xq, output) output[i] = _linear_with_extrap(itp.x, itp.y, xq[i], itp.extrap, deriv, searcher) end diff --git a/src/linear/linear_oneshot.jl b/src/linear/linear_oneshot.jl index 532d8b26..d837af17 100644 --- a/src/linear/linear_oneshot.jl +++ b/src/linear/linear_oneshot.jl @@ -67,7 +67,8 @@ function linear_interp!( @assert length(y) == length(x) "x and y must have same length" @assert length(output) == length(x_targets) "output must match x_targets length" - searcher = _to_searcher(search) + resolved = _resolve_search(search, x_targets) + searcher = _to_searcher(resolved) @boundscheck _check_domain(x, x_targets, extrap) _linear_interp_loop!(output, x, y, x_targets, extrap, deriv, searcher) end @@ -160,7 +161,8 @@ end @assert length(y) == length(x) "x and y must have same length" @assert length(output) == length(x_targets) "output must match x_targets length" - searcher = _to_searcher(search) + resolved = _resolve_search(search, x_targets) + searcher = _to_searcher(resolved) @boundscheck _check_domain(x, x_targets, extrap) _linear_interp_loop!(output, x, y, x_targets, extrap, deriv, searcher) end @@ -390,7 +392,8 @@ end ) where {Tg<:AbstractFloat, Tv, Tq<:Real} @boundscheck length(y) == length(x) || throw(ArgumentError("x and y must have same length")) - searcher = _to_searcher(search, hint) + resolved = _resolve_search(search, xq) + searcher = _to_searcher(resolved, hint) linear_interp(x, y, xq, extrap, deriv, searcher) end diff --git a/src/linear/linear_series_interp.jl b/src/linear/linear_series_interp.jl index e78cfb93..54c3aafc 100644 --- a/src/linear/linear_series_interp.jl +++ b/src/linear/linear_series_interp.jl @@ -490,7 +490,8 @@ function (sitp::LinearSeriesInterpolant{Tg,Tv,P})( xq_typed = Tg(xq_primal) # Build anchor from primal - aq = _make_anchor(sitp, xq_typed, _to_searcher(search, hint)) + resolved = _resolve_search(search, xq) + aq = _make_anchor(sitp, xq_typed, _to_searcher(resolved, hint)) # Dispatch on derivative order with Dual-aware evaluation _eval_linear_series_point!(output, sitp, aq, xq, deriv) # Pass original xq @@ -563,7 +564,8 @@ Pool handles both same-type and mixed-type cases efficiently. # Build anchors - pool handles both same-type and mixed-type cases Tq_eff = promote_type(Tq, Tg) aq_vec = acquire!(pool, _LinearAnchoredQuery{Tg, Tq_eff}, n_query) - _fill_anchors!(aq_vec, sitp.x, xq, Val(:linear); wrap=_should_wrap(sitp), searcher=_to_searcher(search, hint)) + resolved = _resolve_search(search, xq) + _fill_anchors!(aq_vec, sitp.x, xq, Val(:linear); wrap=_should_wrap(sitp), searcher=_to_searcher(resolved, hint)) # Extract matrices for argument-passing pattern y = sitp.y diff --git a/src/quadratic/quadratic_interpolant.jl b/src/quadratic/quadratic_interpolant.jl index 9126772a..c4c52d69 100644 --- a/src/quadratic/quadratic_interpolant.jl +++ b/src/quadratic/quadratic_interpolant.jl @@ -14,7 +14,8 @@ # ───────────────────────────────────────────────────────────── @inline function (itp::QuadraticInterpolant{Tg,Tv})(xq; deriv::DerivOp=EvalValue(), search=itp.search_policy, hint::Union{Nothing,Base.RefValue{Int}}=nothing) where {Tg<:AbstractFloat, Tv} @boundscheck _check_domain(itp.x, xq, itp.extrap) - searcher = _to_searcher(search, hint) + resolved = _resolve_search(search, xq) + searcher = _to_searcher(resolved, hint) # Pass original xq to preserve Dual type for AD _quadratic_eval_at_point(itp.x, itp.y, itp.h, itp.a, itp.d, xq, itp.extrap, deriv, searcher) end @@ -27,7 +28,8 @@ end function (itp::QuadraticInterpolant{Tg,Tv})(xi::AbstractVector{S}; deriv::DerivOp=EvalValue(), search=itp.search_policy, hint::Union{Nothing,Base.RefValue{Int}}=nothing) where {Tg<:AbstractFloat, Tv, S<:Real} T_out = promote_type(Tv, S) # Lossless: wider type to avoid precision loss output = Vector{T_out}(undef, length(xi)) - searcher = _to_searcher(search, hint) + resolved = _resolve_search(search, xi) + searcher = _to_searcher(resolved, hint) @boundscheck _check_domain(itp.x, xi, itp.extrap) @inbounds for i in eachindex(xi, output) output[i] = _quadratic_eval_at_point(itp.x, itp.y, itp.h, itp.a, itp.d, xi[i], itp.extrap, deriv, searcher) @@ -41,7 +43,8 @@ end # ───────────────────────────────────────────────────────────── function (itp::QuadraticInterpolant{Tg,Tv})(output::AbstractVector, xi::AbstractVector{<:Real}; deriv::DerivOp=EvalValue(), search=itp.search_policy, hint::Union{Nothing,Base.RefValue{Int}}=nothing) where {Tg<:AbstractFloat, Tv} @assert length(output) == length(xi) "output length must match xi length" - searcher = _to_searcher(search, hint) + resolved = _resolve_search(search, xi) + searcher = _to_searcher(resolved, hint) @boundscheck _check_domain(itp.x, xi, itp.extrap) @inbounds for i in eachindex(xi, output) output[i] = _quadratic_eval_at_point(itp.x, itp.y, itp.h, itp.a, itp.d, xi[i], itp.extrap, deriv, searcher) diff --git a/src/quadratic/quadratic_oneshot.jl b/src/quadratic/quadratic_oneshot.jl index a49392b5..19579928 100644 --- a/src/quadratic/quadratic_oneshot.jl +++ b/src/quadratic/quadratic_oneshot.jl @@ -217,7 +217,8 @@ vals = quadratic_interp(x, y, sorted_queries; search=LinearBinary(linear_window= bc_promoted = _promote_bc(bc, Tv) _compute_quadratic_coeffs!(h, d, a, x, y, bc_promoted) - searcher = _to_searcher(search, hint) + resolved = _resolve_search(search, xq) + searcher = _to_searcher(resolved, hint) _quadratic_eval_at_point(x, y, h, a, d, xq, extrap, deriv, searcher) end @@ -271,7 +272,8 @@ quadratic_interp!(output, x, y, sorted_queries; search=LinearBinary(linear_windo bc_promoted = _promote_bc(bc, Tv) _compute_quadratic_coeffs!(h, d, a, x, y, bc_promoted) - searcher = _to_searcher(search) + resolved = _resolve_search(search, x_targets) + searcher = _to_searcher(resolved) @boundscheck _check_domain(x, x_targets, extrap) @inbounds for i in eachindex(x_targets, output) output[i] = _quadratic_eval_at_point(x, y, h, a, d, x_targets[i], extrap, deriv, searcher) diff --git a/src/quadratic/quadratic_series_interp.jl b/src/quadratic/quadratic_series_interp.jl index 26f84088..d73dbd4a 100644 --- a/src/quadratic/quadratic_series_interp.jl +++ b/src/quadratic/quadratic_series_interp.jl @@ -512,7 +512,8 @@ function (sitp::QuadraticSeriesInterpolant{Tg,Tv,P})( # Promote for anchor: Int→Float, Int-backed Dual→Float-backed Dual (no-op for Float/Float-backed Dual) xq_promoted = _promote_for_anchor(xq, Tg) T_out = promote_type(Tv, typeof(xq_promoted)) - aq = _anchor_query(sitp.x, xq_promoted, Val(:quadratic); wrap=_should_wrap(sitp), searcher=_to_searcher(search, hint)) + resolved = _resolve_search(search, xq) + aq = _anchor_query(sitp.x, xq_promoted, Val(:quadratic); wrap=_should_wrap(sitp), searcher=_to_searcher(resolved, hint)) output = Vector{T_out}(undef, n_series(sitp)) _eval_series_at_anchor!(output, sitp, aq, deriv) @@ -540,7 +541,8 @@ function (sitp::QuadraticSeriesInterpolant{Tg,Tv,P})( # Promote for anchor: Int→Float, Int-backed Dual→Float-backed Dual xq_promoted = _promote_for_anchor(xq, Tg) - aq = _anchor_query(sitp.x, xq_promoted, Val(:quadratic); wrap=_should_wrap(sitp), searcher=_to_searcher(search, hint)) + resolved = _resolve_search(search, xq) + aq = _anchor_query(sitp.x, xq_promoted, Val(:quadratic); wrap=_should_wrap(sitp), searcher=_to_searcher(resolved, hint)) _eval_series_at_anchor!(output, sitp, aq, deriv) return output @@ -600,12 +602,13 @@ Pool handles both same-type and mixed-type cases efficiently. # Build anchors - pool handles both same-type and mixed-type cases Tq_eff = promote_type(Tq, Tg) aq_vec = acquire!(pool, _QuadraticAnchoredQuery{Tg, Tq_eff}, n_query) + resolved = _resolve_search(search, xq) if Tq === Tg - _fill_anchors!(aq_vec, sitp.x, xq, Val(:quadratic); wrap=_should_wrap(sitp), searcher=_to_searcher(search, hint)) + _fill_anchors!(aq_vec, sitp.x, xq, Val(:quadratic); wrap=_should_wrap(sitp), searcher=_to_searcher(resolved, hint)) else # Mixed type: convert query points to preserve precision xq_promoted = _promote_for_anchor.(xq, Tg) - _fill_anchors!(aq_vec, sitp.x, xq_promoted, Val(:quadratic); wrap=_should_wrap(sitp), searcher=_to_searcher(search, hint)) + _fill_anchors!(aq_vec, sitp.x, xq_promoted, Val(:quadratic); wrap=_should_wrap(sitp), searcher=_to_searcher(resolved, hint)) end # Evaluate all series - anchor already has correct dL precision From 317d93b402c0b2123f5f6347f68c2500bdbdf88b Mon Sep 17 00:00:00 2001 From: Min-Gu Yoo Date: Tue, 24 Feb 2026 12:54:41 -0800 Subject: [PATCH 08/16] Insert _resolve_search at all ND eval sites Adds per-axis AutoSearch resolution via map after _resolve_search_nd in all ND eval callables and oneshot functions. Scalar ND queries resolve to Binary(), SoA/AoS batch queries resolve to LinearBinary(). --- src/constant/nd/constant_nd_eval.jl | 3 +++ src/constant/nd/constant_nd_oneshot.jl | 5 +++++ src/cubic/nd/cubic_nd_eval.jl | 3 +++ src/cubic/nd/cubic_nd_oneshot.jl | 3 +++ src/linear/nd/linear_nd_eval.jl | 3 +++ src/linear/nd/linear_nd_oneshot.jl | 3 +++ src/quadratic/nd/quadratic_nd_eval.jl | 3 +++ src/quadratic/nd/quadratic_nd_oneshot.jl | 3 +++ 8 files changed, 26 insertions(+) diff --git a/src/constant/nd/constant_nd_eval.jl b/src/constant/nd/constant_nd_eval.jl index d205a3fc..32de296c 100644 --- a/src/constant/nd/constant_nd_eval.jl +++ b/src/constant/nd/constant_nd_eval.jl @@ -18,6 +18,7 @@ ) where {Tg, Tv, N} ops = _resolve_deriv_nd(deriv, Val(N)) search_tuple = _resolve_search_nd(search, Val(N)) + search_tuple = map(p -> _resolve_search(p, first(query)), search_tuple) return _eval_constant_nd(itp, query, ops, search_tuple, hint) end @@ -49,6 +50,7 @@ function (itp::ConstantInterpolantND{Tg,Tv,N})( end ops = _resolve_deriv_nd(deriv, Val(N)) search_tuple = _resolve_search_nd(search, Val(N)) + search_tuple = map(p -> _resolve_search(p, first(queries)), search_tuple) if _has_any_derivative(ops, Val(N)) fill!(output, zero(eltype(output))) return output @@ -76,6 +78,7 @@ function (itp::ConstantInterpolantND{Tg,Tv,N})( )) ops = _resolve_deriv_nd(deriv, Val(N)) search_tuple = _resolve_search_nd(search, Val(N)) + search_tuple = map(p -> _resolve_search(p, queries), search_tuple) if _has_any_derivative(ops, Val(N)) fill!(output, zero(eltype(output))) return output diff --git a/src/constant/nd/constant_nd_oneshot.jl b/src/constant/nd/constant_nd_oneshot.jl index f8286dd6..4b0512d5 100644 --- a/src/constant/nd/constant_nd_oneshot.jl +++ b/src/constant/nd/constant_nd_oneshot.jl @@ -165,6 +165,7 @@ function constant_interp( sides = _resolve_side_nd(side, Val(N)) searches = _resolve_search_nd(search, Val(N)) + searches = map(p -> _resolve_search(p, first(query)), searches) extraps_val = _resolve_extrap_nd(extrap, nothing, Val(N)) return _constant_interp_nd_oneshot( @@ -198,6 +199,7 @@ function constant_interp( sides = _resolve_side_nd(side, Val(N)) searches = _resolve_search_nd(search, Val(N)) + searches = map(p -> _resolve_search(p, first(queries)), searches) extraps_val = _resolve_extrap_nd(extrap, nothing, Val(N)) return _constant_interp_nd_oneshot_soa( @@ -231,6 +233,7 @@ function constant_interp( sides = _resolve_side_nd(side, Val(N)) searches = _resolve_search_nd(search, Val(N)) + searches = map(p -> _resolve_search(p, queries), searches) extraps_val = _resolve_extrap_nd(extrap, nothing, Val(N)) return _constant_interp_nd_oneshot_aos( @@ -269,6 +272,7 @@ function constant_interp!( sides = _resolve_side_nd(side, Val(N)) searches = _resolve_search_nd(search, Val(N)) + searches = map(p -> _resolve_search(p, first(queries)), searches) extraps_val = _resolve_extrap_nd(extrap, nothing, Val(N)) return _constant_interp_nd_oneshot_soa!( @@ -303,6 +307,7 @@ function constant_interp!( sides = _resolve_side_nd(side, Val(N)) searches = _resolve_search_nd(search, Val(N)) + searches = map(p -> _resolve_search(p, queries), searches) extraps_val = _resolve_extrap_nd(extrap, nothing, Val(N)) return _constant_interp_nd_oneshot_aos!( diff --git a/src/cubic/nd/cubic_nd_eval.jl b/src/cubic/nd/cubic_nd_eval.jl index 386c6c02..2def1103 100644 --- a/src/cubic/nd/cubic_nd_eval.jl +++ b/src/cubic/nd/cubic_nd_eval.jl @@ -42,6 +42,7 @@ itp((1.0, 0.5); deriv=(DerivOp(1), EvalValue())) # ∂f/∂x only # Note: Don't convert to Tg - preserve query type for AD support ops = _resolve_deriv_nd(deriv, Val(N)) search_tuple = _resolve_search_nd(search, Val(N)) + search_tuple = map(p -> _resolve_search(p, first(query)), search_tuple) return _eval_nd_hermite(itp, query, ops, search_tuple, hint) end @@ -73,6 +74,7 @@ function (itp::CubicInterpolantND{Tg, Tv, N})( end ops = _resolve_deriv_nd(deriv, Val(N)) search_tuple = _resolve_search_nd(search, Val(N)) + search_tuple = map(p -> _resolve_search(p, first(queries)), search_tuple) _batch_nd_soa!(output, itp, queries, ops, search_tuple, hint) return output end @@ -96,6 +98,7 @@ function (itp::CubicInterpolantND{Tg, Tv, N})( )) ops = _resolve_deriv_nd(deriv, Val(N)) search_tuple = _resolve_search_nd(search, Val(N)) + search_tuple = map(p -> _resolve_search(p, queries), search_tuple) _batch_nd_aos!(output, itp, queries, ops, search_tuple, hint) return output end diff --git a/src/cubic/nd/cubic_nd_oneshot.jl b/src/cubic/nd/cubic_nd_oneshot.jl index 8ae5bd3e..3f886e1d 100644 --- a/src/cubic/nd/cubic_nd_oneshot.jl +++ b/src/cubic/nd/cubic_nd_oneshot.jl @@ -44,6 +44,7 @@ function cubic_interp( bcs = _resolve_bcs_nd(bc, Val(N)) searches = _resolve_search_nd(search, Val(N)) + searches = map(p -> _resolve_search(p, first(query)), searches) # Validate BC requirements (once, before dispatch). _validate_nd_bcs!(grids_typed, bcs, data, Val(N)) @@ -268,6 +269,7 @@ function cubic_interp!( bcs = _resolve_bcs_nd(bc, Val(N)) searches = _resolve_search_nd(search, Val(N)) + searches = map(p -> _resolve_search(p, first(queries)), searches) _validate_nd_bcs!(grids_typed, bcs, data, Val(N)) @@ -300,6 +302,7 @@ function cubic_interp!( bcs = _resolve_bcs_nd(bc, Val(N)) searches = _resolve_search_nd(search, Val(N)) + searches = map(p -> _resolve_search(p, queries), searches) _validate_nd_bcs!(grids_typed, bcs, data, Val(N)) diff --git a/src/linear/nd/linear_nd_eval.jl b/src/linear/nd/linear_nd_eval.jl index 721b9c77..7f8731ad 100644 --- a/src/linear/nd/linear_nd_eval.jl +++ b/src/linear/nd/linear_nd_eval.jl @@ -22,6 +22,7 @@ ) where {Tg, Tv, N} ops = _resolve_deriv_nd(deriv, Val(N)) search_tuple = _resolve_search_nd(search, Val(N)) + search_tuple = map(p -> _resolve_search(p, first(query)), search_tuple) return _eval_linear_nd(itp, query, ops, search_tuple, hint) end @@ -53,6 +54,7 @@ function (itp::LinearInterpolantND{Tg,Tv,N})( end ops = _resolve_deriv_nd(deriv, Val(N)) search_tuple = _resolve_search_nd(search, Val(N)) + search_tuple = map(p -> _resolve_search(p, first(queries)), search_tuple) if _has_second_or_higher_derivative(ops, Val(N)) fill!(output, zero(eltype(output))) return output @@ -80,6 +82,7 @@ function (itp::LinearInterpolantND{Tg,Tv,N})( )) ops = _resolve_deriv_nd(deriv, Val(N)) search_tuple = _resolve_search_nd(search, Val(N)) + search_tuple = map(p -> _resolve_search(p, queries), search_tuple) if _has_second_or_higher_derivative(ops, Val(N)) fill!(output, zero(eltype(output))) return output diff --git a/src/linear/nd/linear_nd_oneshot.jl b/src/linear/nd/linear_nd_oneshot.jl index 459b6390..5dbaa759 100644 --- a/src/linear/nd/linear_nd_oneshot.jl +++ b/src/linear/nd/linear_nd_oneshot.jl @@ -122,6 +122,7 @@ function linear_interp( Tr = promote_type(Tv, Tg, typeof.(query)...) searches = _resolve_search_nd(search, Val(N)) + searches = map(p -> _resolve_search(p, first(query)), searches) extraps_val = _resolve_extrap_nd(extrap, nothing, Val(N)) ops = _resolve_deriv_nd(deriv, Val(N)) @@ -197,6 +198,7 @@ function linear_interp!( _validate_nd_grids(grids_typed, data) searches = _resolve_search_nd(search, Val(N)) + searches = map(p -> _resolve_search(p, first(queries)), searches) extraps_val = _resolve_extrap_nd(extrap, nothing, Val(N)) ops = _resolve_deriv_nd(deriv, Val(N)) @@ -224,6 +226,7 @@ function linear_interp!( _validate_nd_grids(grids_typed, data) searches = _resolve_search_nd(search, Val(N)) + searches = map(p -> _resolve_search(p, queries), searches) extraps_val = _resolve_extrap_nd(extrap, nothing, Val(N)) ops = _resolve_deriv_nd(deriv, Val(N)) diff --git a/src/quadratic/nd/quadratic_nd_eval.jl b/src/quadratic/nd/quadratic_nd_eval.jl index 9e641e8d..015a5c61 100644 --- a/src/quadratic/nd/quadratic_nd_eval.jl +++ b/src/quadratic/nd/quadratic_nd_eval.jl @@ -67,6 +67,7 @@ itp((1.0, 0.5); deriv=(DerivOp(1), EvalValue())) # ∂f/∂x only ) where {Tg, Tv, N} ops = _resolve_deriv_nd(deriv, Val(N)) search_tuple = _resolve_search_nd(search, Val(N)) + search_tuple = map(p -> _resolve_search(p, first(query)), search_tuple) return _eval_nd_quadratic(itp, query, ops, search_tuple, hint) end @@ -98,6 +99,7 @@ function (itp::QuadraticInterpolantND{Tg, Tv, N})( end ops = _resolve_deriv_nd(deriv, Val(N)) search_tuple = _resolve_search_nd(search, Val(N)) + search_tuple = map(p -> _resolve_search(p, first(queries)), search_tuple) _batch_nd_soa!(output, itp, queries, ops, search_tuple, hint) return output end @@ -121,6 +123,7 @@ function (itp::QuadraticInterpolantND{Tg, Tv, N})( )) ops = _resolve_deriv_nd(deriv, Val(N)) search_tuple = _resolve_search_nd(search, Val(N)) + search_tuple = map(p -> _resolve_search(p, queries), search_tuple) _batch_nd_aos!(output, itp, queries, ops, search_tuple, hint) return output end diff --git a/src/quadratic/nd/quadratic_nd_oneshot.jl b/src/quadratic/nd/quadratic_nd_oneshot.jl index f1a0d7ab..be0cae96 100644 --- a/src/quadratic/nd/quadratic_nd_oneshot.jl +++ b/src/quadratic/nd/quadratic_nd_oneshot.jl @@ -154,6 +154,7 @@ function quadratic_interp( bcs = _resolve_bcs_nd_quadratic(bc, Val(N)) searches = _resolve_search_nd(search, Val(N)) + searches = map(p -> _resolve_search(p, first(query)), searches) extraps_val = _resolve_extrap_nd(extrap, bcs, Val(N)) ops = _resolve_deriv_nd(deriv, Val(N)) @@ -234,6 +235,7 @@ function quadratic_interp!( bcs = _resolve_bcs_nd_quadratic(bc, Val(N)) searches = _resolve_search_nd(search, Val(N)) + searches = map(p -> _resolve_search(p, first(queries)), searches) extraps_val = _resolve_extrap_nd(extrap, bcs, Val(N)) ops = _resolve_deriv_nd(deriv, Val(N)) @@ -264,6 +266,7 @@ function quadratic_interp!( bcs = _resolve_bcs_nd_quadratic(bc, Val(N)) searches = _resolve_search_nd(search, Val(N)) + searches = map(p -> _resolve_search(p, queries), searches) extraps_val = _resolve_extrap_nd(extrap, bcs, Val(N)) ops = _resolve_deriv_nd(deriv, Val(N)) From 4948294aebb4b15a205216e923e52198f21b2d42 Mon Sep 17 00:00:00 2001 From: Min-Gu Yoo Date: Tue, 24 Feb 2026 13:02:14 -0800 Subject: [PATCH 09/16] Fix test failures: add Searcher passthrough and update default policy assertions - Add _resolve_search(::Searcher, _) passthrough after Searcher struct definition so pre-built Searcher objects injected via `search=` keyword skip resolution - Update test_search.jl default policy assertion from LinearBinary{2} to AutoSearch --- src/core/search.jl | 4 ++++ test/test_search.jl | 6 +++--- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/src/core/search.jl b/src/core/search.jl index 83e881f3..884dc364 100644 --- a/src/core/search.jl +++ b/src/core/search.jl @@ -307,6 +307,10 @@ struct Searcher{P<:AbstractSearchPolicy,H<:AbstractHint} hint::H end +# Searcher passthrough for _resolve_search: pre-built Searcher objects skip resolution. +# Must be defined after Searcher struct (Julia requires types to be defined before use). +@inline _resolve_search(s::Searcher, _) = s + """ DEFAULT_SEARCHER diff --git a/test/test_search.jl b/test/test_search.jl index 16d159cb..dc42b338 100644 --- a/test/test_search.jl +++ b/test/test_search.jl @@ -773,9 +773,9 @@ using FastInterpolations: search_interval, _search_binary, _search_direct, _sear itp_hb = linear_interp(x, y; search=HintedBinary()) @test itp_hb.search_policy isa HintedBinary - # Create with default (LinearBinary) - itp_bin = linear_interp(x, y) - @test itp_bin.search_policy isa LinearBinary{2} + # Create with default (AutoSearch) + itp_auto = linear_interp(x, y) + @test itp_auto.search_policy isa AutoSearch end @testset "Baked-in policy with hint: hint updates when policy supports it" begin From fb3f232367d39f458d30fd81512678366cba3513 Mon Sep 17 00:00:00 2001 From: Min-Gu Yoo Date: Tue, 24 Feb 2026 13:06:29 -0800 Subject: [PATCH 10/16] Update docstrings and show tests for AutoSearch default - Replace "uses LinearBinary() by default" with AutoSearch resolution docs - Add AutoSearch format test and default interpolant show test - Update Binary() docstring default reference --- src/constant/constant_interpolant.jl | 7 ++++--- src/constant/constant_types.jl | 9 +++++---- src/core/search.jl | 2 +- src/cubic/cubic_interpolant.jl | 9 +++++---- src/cubic/cubic_types.jl | 11 ++++++----- src/derivative_view.jl | 4 ++-- src/linear/linear_interpolant.jl | 7 ++++--- src/quadratic/quadratic_interpolant.jl | 7 ++++--- src/quadratic/quadratic_types.jl | 9 +++++---- test/test_show.jl | 6 ++++++ 10 files changed, 42 insertions(+), 29 deletions(-) diff --git a/src/constant/constant_interpolant.jl b/src/constant/constant_interpolant.jl index 627a07c2..f48f41c2 100644 --- a/src/constant/constant_interpolant.jl +++ b/src/constant/constant_interpolant.jl @@ -85,9 +85,10 @@ y = [1.0+2.0im, 3.0+4.0im, 5.0+6.0im] itp = constant_interp(x, y) itp(0.5) # 1.0+2.0im (ComplexF64) -# Create with custom search policy -itp = constant_interp(x, y; search=LinearBinary()) -val = itp(0.5) # uses LinearBinary() by default +# Search policy: AutoSearch adapts to query type (scalar→Binary, vector→LinearBinary) +itp = constant_interp(x, y) +val = itp(0.5) # AutoSearch resolves to Binary() for scalar +itp = constant_interp(x, y; search=LinearBinary()) # explicit override # Fused broadcast (optimal) result = @. coef * itp(query) diff --git a/src/constant/constant_types.jl b/src/constant/constant_types.jl index b61f0aba..aafab45c 100644 --- a/src/constant/constant_types.jl +++ b/src/constant/constant_types.jl @@ -43,10 +43,11 @@ val = itp(0.5) # returns ComplexF64 itp_left = constant_interp(x, y; side=LeftSide()) itp_wrap = constant_interp(x, y; extrap=WrapExtrap(), side=RightSide()) -# Custom search policy -itp = constant_interp(x, y; search=LinearBinary()) -val = itp(0.5) # uses LinearBinary() by default -val = itp(0.5; search=Binary()) # override with Binary() +# Search policy: AutoSearch adapts to query type (scalar→Binary, vector→LinearBinary) +itp = constant_interp(x, y) +val = itp(0.5) # AutoSearch resolves to Binary() for scalar +itp = constant_interp(x, y; search=LinearBinary()) # explicit override +val = itp(0.5; search=Binary()) # per-call override ``` """ struct ConstantInterpolant{Tg<:AbstractFloat, Tv, X<:AbstractVector{Tg}, Y<:AbstractVector{Tv}, E<:AbstractExtrap, SD<:AbstractSide, P<:AbstractSearchPolicy} <: AbstractInterpolant{Tg, Tv} diff --git a/src/core/search.jl b/src/core/search.jl index 884dc364..afbb6612 100644 --- a/src/core/search.jl +++ b/src/core/search.jl @@ -41,7 +41,7 @@ pure binary search is used. # Example ```julia val = linear_interp(x, y, 0.5; search=Binary()) # explicit binary search -val = linear_interp(x, y, 0.5) # default: LinearBinary() +val = linear_interp(x, y, 0.5) # default: AutoSearch() # With hint: auto-upgrades to HintedBinary behavior hint = Ref(1) diff --git a/src/cubic/cubic_interpolant.jl b/src/cubic/cubic_interpolant.jl index 8c9dae7e..1840b655 100644 --- a/src/cubic/cubic_interpolant.jl +++ b/src/cubic/cubic_interpolant.jl @@ -376,10 +376,11 @@ val = itp(0.5) # Scalar (zero-allocation) vals = itp.(query_points) # Broadcast result = @. coef * itp(rho) * ne # Fused broadcast -# Create with custom search policy -itp = cubic_interp(x, y; search=LinearBinary()) -val = itp(0.5) # Uses LinearBinary() by default -val = itp(0.5; search=Binary()) # Override with Binary() +# Search policy: AutoSearch adapts to query type (scalar→Binary, vector→LinearBinary) +itp = cubic_interp(x, y) +val = itp(0.5) # AutoSearch resolves to Binary() for scalar +itp = cubic_interp(x, y; search=LinearBinary()) # explicit override +val = itp(0.5; search=Binary()) # per-call override # Complex values x = [0.0, 1.0, 2.0, 3.0, 4.0] diff --git a/src/cubic/cubic_types.jl b/src/cubic/cubic_types.jl index 114888b8..1a984c9a 100644 --- a/src/cubic/cubic_types.jl +++ b/src/cubic/cubic_types.jl @@ -93,7 +93,7 @@ Returned by `cubic_interp(x, y)` (2-argument form). - `Tv`: Value type for y-values (can be Tg, Complex{Tg}, or other Number) - `C`: CubicSplineCache type (preserves grid type info for O(1) vs O(log n) lookup) - `E`: Extrapolation mode type (compile-time specialized) -- `P`: Search policy type (Binary, HintedBinary, LinearBinary, etc.) +- `P`: Search policy type (AutoSearch, Binary, HintedBinary, LinearBinary, etc.) - `BC`: Boundary condition type (BCPair or PeriodicBC) # Fields @@ -110,10 +110,11 @@ itp = cubic_interp(x, y) result = @. coef * itp(rho) * other_terms # fused, zero-allocation per call val = itp(0.5) # scalar (zero-allocation) -# Create with custom search policy -itp = cubic_interp(x, y; search=LinearBinary()) -val = itp(0.5) # uses LinearBinary() by default -val = itp(0.5; search=Binary()) # override with Binary() +# Search policy: AutoSearch adapts to query type (scalar→Binary, vector→LinearBinary) +itp = cubic_interp(x, y) +val = itp(0.5) # AutoSearch resolves to Binary() for scalar +itp = cubic_interp(x, y; search=LinearBinary()) # explicit override +val = itp(0.5; search=Binary()) # per-call override # Complex values x = [0.0, 1.0, 2.0, 3.0, 4.0] diff --git a/src/derivative_view.jl b/src/derivative_view.jl index 9c66f103..caf66de9 100644 --- a/src/derivative_view.jl +++ b/src/derivative_view.jl @@ -131,7 +131,7 @@ d1 = deriv1(itp) val = d1(0.5) # Override search policy -val = d1(0.5; search=LinearBinary()) +val = d1(0.5; search=Binary()) # Use hint for sequential access hint = Ref(1) @@ -150,7 +150,7 @@ results = @. itp(query_pts) + 0.1 * d1(query_pts) # In-place evaluation (single-series) output = similar(query_pts) d1(output, query_pts) -d1(output, query_pts; search=LinearBinary()) +d1(output, query_pts; search=LinearBinary()) # explicit: LinearBinary for sorted vectors ``` """ (deriv1, deriv2, deriv3) diff --git a/src/linear/linear_interpolant.jl b/src/linear/linear_interpolant.jl index 2b9ff7c5..8fdec84b 100644 --- a/src/linear/linear_interpolant.jl +++ b/src/linear/linear_interpolant.jl @@ -97,8 +97,8 @@ Can be: # Create with default AutoSearch() 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()) +# Default AutoSearch: scalar→Binary, vector→LinearBinary +itp = linear_interp(x_data, y_data) # Scalar call (uses stored policy) val = itp(0.5) @@ -129,7 +129,8 @@ 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 +- Default `AutoSearch()` adapts: scalar→`Binary()`, vector→`LinearBinary()` +- Use `search=LinearBinary()` to force linear-binary for all query types - Use `hint=Ref(idx)` for ODE/streaming patterns with persistent hint """ function linear_interp end diff --git a/src/quadratic/quadratic_interpolant.jl b/src/quadratic/quadratic_interpolant.jl index c4c52d69..ed666d42 100644 --- a/src/quadratic/quadratic_interpolant.jl +++ b/src/quadratic/quadratic_interpolant.jl @@ -90,9 +90,10 @@ y = [1.0+2.0im, 3.0+4.0im, 5.0+6.0im, 7.0+8.0im] itp = quadratic_interp(x, y) itp(0.5) # returns ComplexF64 -# Create with custom search policy -itp = quadratic_interp(x, y; search=LinearBinary()) -val = itp(0.5) # uses LinearBinary() by default +# Search policy: AutoSearch adapts to query type (scalar→Binary, vector→LinearBinary) +itp = quadratic_interp(x, y) +val = itp(0.5) # AutoSearch resolves to Binary() for scalar +itp = quadratic_interp(x, y; search=LinearBinary()) # explicit override # Fused broadcast (optimal) result = @. coef * itp(query) diff --git a/src/quadratic/quadratic_types.jl b/src/quadratic/quadratic_types.jl index 4aa397f2..c33600a5 100644 --- a/src/quadratic/quadratic_types.jl +++ b/src/quadratic/quadratic_types.jl @@ -49,10 +49,11 @@ y = [1.0+2.0im, 3.0+4.0im, 5.0+6.0im, 7.0+8.0im] itp = quadratic_interp(x, y) val = itp(0.5) # returns ComplexF64 -# Custom search policy -itp = quadratic_interp(x, y; search=LinearBinary()) -val = itp(0.5) # uses LinearBinary() by default -val = itp(0.5; search=Binary()) # override with Binary() +# Search policy: AutoSearch adapts to query type (scalar→Binary, vector→LinearBinary) +itp = quadratic_interp(x, y) +val = itp(0.5) # AutoSearch resolves to Binary() for scalar +itp = quadratic_interp(x, y; search=LinearBinary()) # explicit override +val = itp(0.5; search=Binary()) # per-call override ``` """ struct QuadraticInterpolant{Tg<:AbstractFloat, Tv, X<:AbstractVector{Tg}, Y<:AbstractVector{Tv}, E<:AbstractExtrap, P<:AbstractSearchPolicy} <: AbstractInterpolant{Tg, Tv} diff --git a/test/test_show.jl b/test/test_show.jl index a893b79c..89f3c664 100644 --- a/test/test_show.jl +++ b/test/test_show.jl @@ -215,6 +215,11 @@ verbose_lb16 = sprint(show, MIME("text/plain"), itp_linear_binary_16) @test occursin("LinearBinary{16}", verbose_lb16) + + # Default-constructed interpolant uses AutoSearch + itp_default = linear_interp(x_vec, y) + verbose_default = sprint(show, MIME("text/plain"), itp_default) + @test occursin("AutoSearch", verbose_default) end @testset "Extrapolation mode display" begin @@ -470,6 +475,7 @@ @test FI._format_search(HintedBinary()) == "HintedBinary" @test FI._format_search(LinearBinary()) == "LinearBinary{2}" @test FI._format_search(LinearBinary(linear_window=4)) == "LinearBinary{4}" + @test FI._format_search(AutoSearch()) == "AutoSearch (scalar→Binary, vector→LinearBinary)" # DerivativeView with unknown parent type (no .x or .cache.x) struct DummyInterpolant{T} <: FastInterpolations.AbstractInterpolant{T, T} end From 3700bcbe3ed67ba8aeb65a850506ffd83e4d49b2 Mon Sep 17 00:00:00 2001 From: Min-Gu Yoo Date: Tue, 24 Feb 2026 14:01:28 -0800 Subject: [PATCH 11/16] Refactor ND search resolution: merge 2-line pattern into 3-arg _resolve_search_nd Add _resolve_search_nd(search, Val(N), query_sample) overload that combines dimension broadcast + AutoSearch resolution in one call. Replaces 26 two-line patterns across all ND eval and oneshot files with a single-line call. Zero allocation verified for both interpolant eval and oneshot paths. --- src/constant/nd/constant_nd_eval.jl | 9 +++------ src/constant/nd/constant_nd_oneshot.jl | 15 +++++---------- src/core/nd_utils.jl | 15 ++++++++++++++- src/cubic/nd/cubic_nd_eval.jl | 9 +++------ src/cubic/nd/cubic_nd_oneshot.jl | 9 +++------ src/linear/nd/linear_nd_eval.jl | 9 +++------ src/linear/nd/linear_nd_oneshot.jl | 9 +++------ src/quadratic/nd/quadratic_nd_eval.jl | 9 +++------ src/quadratic/nd/quadratic_nd_oneshot.jl | 9 +++------ 9 files changed, 40 insertions(+), 53 deletions(-) diff --git a/src/constant/nd/constant_nd_eval.jl b/src/constant/nd/constant_nd_eval.jl index 32de296c..9468a848 100644 --- a/src/constant/nd/constant_nd_eval.jl +++ b/src/constant/nd/constant_nd_eval.jl @@ -17,8 +17,7 @@ hint::Union{Nothing, NTuple{N, Base.RefValue{Int}}} = nothing ) where {Tg, Tv, N} ops = _resolve_deriv_nd(deriv, Val(N)) - search_tuple = _resolve_search_nd(search, Val(N)) - search_tuple = map(p -> _resolve_search(p, first(query)), search_tuple) + search_tuple = _resolve_search_nd(search, Val(N), first(query)) return _eval_constant_nd(itp, query, ops, search_tuple, hint) end @@ -49,8 +48,7 @@ function (itp::ConstantInterpolantND{Tg,Tv,N})( )) end ops = _resolve_deriv_nd(deriv, Val(N)) - search_tuple = _resolve_search_nd(search, Val(N)) - search_tuple = map(p -> _resolve_search(p, first(queries)), search_tuple) + search_tuple = _resolve_search_nd(search, Val(N), first(queries)) if _has_any_derivative(ops, Val(N)) fill!(output, zero(eltype(output))) return output @@ -77,8 +75,7 @@ function (itp::ConstantInterpolantND{Tg,Tv,N})( "output length $(length(output)) must match query length $n_queries" )) ops = _resolve_deriv_nd(deriv, Val(N)) - search_tuple = _resolve_search_nd(search, Val(N)) - search_tuple = map(p -> _resolve_search(p, queries), search_tuple) + search_tuple = _resolve_search_nd(search, Val(N), queries) if _has_any_derivative(ops, Val(N)) fill!(output, zero(eltype(output))) return output diff --git a/src/constant/nd/constant_nd_oneshot.jl b/src/constant/nd/constant_nd_oneshot.jl index 4b0512d5..e613b201 100644 --- a/src/constant/nd/constant_nd_oneshot.jl +++ b/src/constant/nd/constant_nd_oneshot.jl @@ -164,8 +164,7 @@ function constant_interp( _validate_nd_grids(grids_typed, data) sides = _resolve_side_nd(side, Val(N)) - searches = _resolve_search_nd(search, Val(N)) - searches = map(p -> _resolve_search(p, first(query)), searches) + searches = _resolve_search_nd(search, Val(N), first(query)) extraps_val = _resolve_extrap_nd(extrap, nothing, Val(N)) return _constant_interp_nd_oneshot( @@ -198,8 +197,7 @@ function constant_interp( _validate_nd_grids(grids_typed, data) sides = _resolve_side_nd(side, Val(N)) - searches = _resolve_search_nd(search, Val(N)) - searches = map(p -> _resolve_search(p, first(queries)), searches) + searches = _resolve_search_nd(search, Val(N), first(queries)) extraps_val = _resolve_extrap_nd(extrap, nothing, Val(N)) return _constant_interp_nd_oneshot_soa( @@ -232,8 +230,7 @@ function constant_interp( _validate_nd_grids(grids_typed, data) sides = _resolve_side_nd(side, Val(N)) - searches = _resolve_search_nd(search, Val(N)) - searches = map(p -> _resolve_search(p, queries), searches) + searches = _resolve_search_nd(search, Val(N), queries) extraps_val = _resolve_extrap_nd(extrap, nothing, Val(N)) return _constant_interp_nd_oneshot_aos( @@ -271,8 +268,7 @@ function constant_interp!( _validate_nd_grids(grids_typed, data) sides = _resolve_side_nd(side, Val(N)) - searches = _resolve_search_nd(search, Val(N)) - searches = map(p -> _resolve_search(p, first(queries)), searches) + searches = _resolve_search_nd(search, Val(N), first(queries)) extraps_val = _resolve_extrap_nd(extrap, nothing, Val(N)) return _constant_interp_nd_oneshot_soa!( @@ -306,8 +302,7 @@ function constant_interp!( _validate_nd_grids(grids_typed, data) sides = _resolve_side_nd(side, Val(N)) - searches = _resolve_search_nd(search, Val(N)) - searches = map(p -> _resolve_search(p, queries), searches) + searches = _resolve_search_nd(search, Val(N), queries) extraps_val = _resolve_extrap_nd(extrap, nothing, Val(N)) return _constant_interp_nd_oneshot_aos!( diff --git a/src/core/nd_utils.jl b/src/core/nd_utils.jl index c5657d92..697b4cef 100644 --- a/src/core/nd_utils.jl +++ b/src/core/nd_utils.jl @@ -106,9 +106,16 @@ end """ _resolve_search_nd(search, Val(N)) -> NTuple{N, AbstractSearchPolicy} -Resolve search policy input to canonical N-tuple. +Resolve search policy input to canonical N-tuple (broadcast only, no AutoSearch resolution). - Single `AbstractSearchPolicy` → broadcast to all N axes - `NTuple{N, AbstractSearchPolicy}` → passthrough + + _resolve_search_nd(search, Val(N), query_sample) -> NTuple{N, AbstractSearchPolicy} + +Broadcast + resolve AutoSearch in one step. `query_sample` determines scalar vs vector resolution: +- `query_sample::Real` (from `first(query)`) → `Binary()` per axis +- `query_sample::AbstractVector` (from `first(queries)` or `queries`) → `LinearBinary()` per axis +- Explicit policies pass through unchanged. """ @inline _resolve_search_nd(s::AbstractSearchPolicy, ::Val{N}) where {N} = ntuple(_ -> s, Val(N)) @@ -118,6 +125,12 @@ Resolve search policy input to canonical N-tuple. throw(ArgumentError("search tuple must have $N elements to match grid dimensions, got $(length(s))")) end +# 3-arg: broadcast + resolve AutoSearch per-axis in one step +@inline function _resolve_search_nd(s, ::Val{N}, query_sample) where {N} + tuple = _resolve_search_nd(s, Val(N)) + return map(p -> _resolve_search(p, query_sample), tuple) +end + # ======================================== # Boundary Condition Resolution # ======================================== diff --git a/src/cubic/nd/cubic_nd_eval.jl b/src/cubic/nd/cubic_nd_eval.jl index 2def1103..aa34b26b 100644 --- a/src/cubic/nd/cubic_nd_eval.jl +++ b/src/cubic/nd/cubic_nd_eval.jl @@ -41,8 +41,7 @@ itp((1.0, 0.5); deriv=(DerivOp(1), EvalValue())) # ∂f/∂x only ) where {Tg, Tv, N} # Note: Don't convert to Tg - preserve query type for AD support ops = _resolve_deriv_nd(deriv, Val(N)) - search_tuple = _resolve_search_nd(search, Val(N)) - search_tuple = map(p -> _resolve_search(p, first(query)), search_tuple) + search_tuple = _resolve_search_nd(search, Val(N), first(query)) return _eval_nd_hermite(itp, query, ops, search_tuple, hint) end @@ -73,8 +72,7 @@ function (itp::CubicInterpolantND{Tg, Tv, N})( )) end ops = _resolve_deriv_nd(deriv, Val(N)) - search_tuple = _resolve_search_nd(search, Val(N)) - search_tuple = map(p -> _resolve_search(p, first(queries)), search_tuple) + search_tuple = _resolve_search_nd(search, Val(N), first(queries)) _batch_nd_soa!(output, itp, queries, ops, search_tuple, hint) return output end @@ -97,8 +95,7 @@ function (itp::CubicInterpolantND{Tg, Tv, N})( "output length $(length(output)) must match query length $n_queries" )) ops = _resolve_deriv_nd(deriv, Val(N)) - search_tuple = _resolve_search_nd(search, Val(N)) - search_tuple = map(p -> _resolve_search(p, queries), search_tuple) + search_tuple = _resolve_search_nd(search, Val(N), queries) _batch_nd_aos!(output, itp, queries, ops, search_tuple, hint) return output end diff --git a/src/cubic/nd/cubic_nd_oneshot.jl b/src/cubic/nd/cubic_nd_oneshot.jl index 3f886e1d..4828af89 100644 --- a/src/cubic/nd/cubic_nd_oneshot.jl +++ b/src/cubic/nd/cubic_nd_oneshot.jl @@ -43,8 +43,7 @@ function cubic_interp( Tr = promote_type(Tv, Tg, typeof.(query)...) bcs = _resolve_bcs_nd(bc, Val(N)) - searches = _resolve_search_nd(search, Val(N)) - searches = map(p -> _resolve_search(p, first(query)), searches) + searches = _resolve_search_nd(search, Val(N), first(query)) # Validate BC requirements (once, before dispatch). _validate_nd_bcs!(grids_typed, bcs, data, Val(N)) @@ -268,8 +267,7 @@ function cubic_interp!( _validate_nd_grids(grids_typed, data) bcs = _resolve_bcs_nd(bc, Val(N)) - searches = _resolve_search_nd(search, Val(N)) - searches = map(p -> _resolve_search(p, first(queries)), searches) + searches = _resolve_search_nd(search, Val(N), first(queries)) _validate_nd_bcs!(grids_typed, bcs, data, Val(N)) @@ -301,8 +299,7 @@ function cubic_interp!( _validate_nd_grids(grids_typed, data) bcs = _resolve_bcs_nd(bc, Val(N)) - searches = _resolve_search_nd(search, Val(N)) - searches = map(p -> _resolve_search(p, queries), searches) + searches = _resolve_search_nd(search, Val(N), queries) _validate_nd_bcs!(grids_typed, bcs, data, Val(N)) diff --git a/src/linear/nd/linear_nd_eval.jl b/src/linear/nd/linear_nd_eval.jl index 7f8731ad..2145adde 100644 --- a/src/linear/nd/linear_nd_eval.jl +++ b/src/linear/nd/linear_nd_eval.jl @@ -21,8 +21,7 @@ hint::Union{Nothing, NTuple{N, Base.RefValue{Int}}} = nothing ) where {Tg, Tv, N} ops = _resolve_deriv_nd(deriv, Val(N)) - search_tuple = _resolve_search_nd(search, Val(N)) - search_tuple = map(p -> _resolve_search(p, first(query)), search_tuple) + search_tuple = _resolve_search_nd(search, Val(N), first(query)) return _eval_linear_nd(itp, query, ops, search_tuple, hint) end @@ -53,8 +52,7 @@ function (itp::LinearInterpolantND{Tg,Tv,N})( )) end ops = _resolve_deriv_nd(deriv, Val(N)) - search_tuple = _resolve_search_nd(search, Val(N)) - search_tuple = map(p -> _resolve_search(p, first(queries)), search_tuple) + search_tuple = _resolve_search_nd(search, Val(N), first(queries)) if _has_second_or_higher_derivative(ops, Val(N)) fill!(output, zero(eltype(output))) return output @@ -81,8 +79,7 @@ function (itp::LinearInterpolantND{Tg,Tv,N})( "output length $(length(output)) must match query length $n_queries" )) ops = _resolve_deriv_nd(deriv, Val(N)) - search_tuple = _resolve_search_nd(search, Val(N)) - search_tuple = map(p -> _resolve_search(p, queries), search_tuple) + search_tuple = _resolve_search_nd(search, Val(N), queries) if _has_second_or_higher_derivative(ops, Val(N)) fill!(output, zero(eltype(output))) return output diff --git a/src/linear/nd/linear_nd_oneshot.jl b/src/linear/nd/linear_nd_oneshot.jl index 5dbaa759..3f14ca72 100644 --- a/src/linear/nd/linear_nd_oneshot.jl +++ b/src/linear/nd/linear_nd_oneshot.jl @@ -121,8 +121,7 @@ function linear_interp( _validate_nd_grids(grids_typed, data) Tr = promote_type(Tv, Tg, typeof.(query)...) - searches = _resolve_search_nd(search, Val(N)) - searches = map(p -> _resolve_search(p, first(query)), searches) + searches = _resolve_search_nd(search, Val(N), first(query)) extraps_val = _resolve_extrap_nd(extrap, nothing, Val(N)) ops = _resolve_deriv_nd(deriv, Val(N)) @@ -197,8 +196,7 @@ function linear_interp!( grids_typed = _convert_grids_typed(grids, Tg) _validate_nd_grids(grids_typed, data) - searches = _resolve_search_nd(search, Val(N)) - searches = map(p -> _resolve_search(p, first(queries)), searches) + searches = _resolve_search_nd(search, Val(N), first(queries)) extraps_val = _resolve_extrap_nd(extrap, nothing, Val(N)) ops = _resolve_deriv_nd(deriv, Val(N)) @@ -225,8 +223,7 @@ function linear_interp!( grids_typed = _convert_grids_typed(grids, Tg) _validate_nd_grids(grids_typed, data) - searches = _resolve_search_nd(search, Val(N)) - searches = map(p -> _resolve_search(p, queries), searches) + searches = _resolve_search_nd(search, Val(N), queries) extraps_val = _resolve_extrap_nd(extrap, nothing, Val(N)) ops = _resolve_deriv_nd(deriv, Val(N)) diff --git a/src/quadratic/nd/quadratic_nd_eval.jl b/src/quadratic/nd/quadratic_nd_eval.jl index 015a5c61..5bfb1167 100644 --- a/src/quadratic/nd/quadratic_nd_eval.jl +++ b/src/quadratic/nd/quadratic_nd_eval.jl @@ -66,8 +66,7 @@ itp((1.0, 0.5); deriv=(DerivOp(1), EvalValue())) # ∂f/∂x only hint::Union{Nothing, NTuple{N, Base.RefValue{Int}}}=nothing ) where {Tg, Tv, N} ops = _resolve_deriv_nd(deriv, Val(N)) - search_tuple = _resolve_search_nd(search, Val(N)) - search_tuple = map(p -> _resolve_search(p, first(query)), search_tuple) + search_tuple = _resolve_search_nd(search, Val(N), first(query)) return _eval_nd_quadratic(itp, query, ops, search_tuple, hint) end @@ -98,8 +97,7 @@ function (itp::QuadraticInterpolantND{Tg, Tv, N})( )) end ops = _resolve_deriv_nd(deriv, Val(N)) - search_tuple = _resolve_search_nd(search, Val(N)) - search_tuple = map(p -> _resolve_search(p, first(queries)), search_tuple) + search_tuple = _resolve_search_nd(search, Val(N), first(queries)) _batch_nd_soa!(output, itp, queries, ops, search_tuple, hint) return output end @@ -122,8 +120,7 @@ function (itp::QuadraticInterpolantND{Tg, Tv, N})( "output length $(length(output)) must match query length $n_queries" )) ops = _resolve_deriv_nd(deriv, Val(N)) - search_tuple = _resolve_search_nd(search, Val(N)) - search_tuple = map(p -> _resolve_search(p, queries), search_tuple) + search_tuple = _resolve_search_nd(search, Val(N), queries) _batch_nd_aos!(output, itp, queries, ops, search_tuple, hint) return output end diff --git a/src/quadratic/nd/quadratic_nd_oneshot.jl b/src/quadratic/nd/quadratic_nd_oneshot.jl index be0cae96..cc8fd07f 100644 --- a/src/quadratic/nd/quadratic_nd_oneshot.jl +++ b/src/quadratic/nd/quadratic_nd_oneshot.jl @@ -153,8 +153,7 @@ function quadratic_interp( Tr = promote_type(Tv, Tg, typeof.(query)...) bcs = _resolve_bcs_nd_quadratic(bc, Val(N)) - searches = _resolve_search_nd(search, Val(N)) - searches = map(p -> _resolve_search(p, first(query)), searches) + searches = _resolve_search_nd(search, Val(N), first(query)) extraps_val = _resolve_extrap_nd(extrap, bcs, Val(N)) ops = _resolve_deriv_nd(deriv, Val(N)) @@ -234,8 +233,7 @@ function quadratic_interp!( _validate_nd_grids(grids_typed, data) bcs = _resolve_bcs_nd_quadratic(bc, Val(N)) - searches = _resolve_search_nd(search, Val(N)) - searches = map(p -> _resolve_search(p, first(queries)), searches) + searches = _resolve_search_nd(search, Val(N), first(queries)) extraps_val = _resolve_extrap_nd(extrap, bcs, Val(N)) ops = _resolve_deriv_nd(deriv, Val(N)) @@ -265,8 +263,7 @@ function quadratic_interp!( _validate_nd_grids(grids_typed, data) bcs = _resolve_bcs_nd_quadratic(bc, Val(N)) - searches = _resolve_search_nd(search, Val(N)) - searches = map(p -> _resolve_search(p, queries), searches) + searches = _resolve_search_nd(search, Val(N), queries) extraps_val = _resolve_extrap_nd(extrap, bcs, Val(N)) ops = _resolve_deriv_nd(deriv, Val(N)) From 0db609d1dcf6b3c1dfa0a7a9f576e874bccdf97b Mon Sep 17 00:00:00 2001 From: Min-Gu Yoo Date: Tue, 24 Feb 2026 15:27:20 -0800 Subject: [PATCH 12/16] Fix AutoSearch post-impl issues: OOB clamp, missing resolution sites, weak tests - search.jl: clamp hint in _search_linear_binary! to guard against user-provided out-of-range hints (Ref(0), stale hints from shorter grids); mirrors _search_linear! - vector_calculus.jl: all 5 @generated fns (gradient/gradient!/hessian/hessian!/laplacian) now use 3-arg _resolve_search_nd with first(query) so AutoSearch resolves correctly - integrate_api.jl: all 8 integrate fns now call _resolve_search before _to_searcher, routing AutoSearch through the intended path instead of the safety-net fallback - nd_utils.jl: document AoS batch dispatch semantics in _resolve_search_nd docstring - test_search.jl: replace weak isfinite hint test with value+position assertions; add OOB clamp safety test; add gradient/hessian/laplacian AutoSearch resolution tests --- src/core/nd_utils.jl | 6 +- src/core/search.jl | 5 +- src/integral/integrate_api.jl | 24 +- src/vector_calculus.jl | 10 +- test/test_search.jl | 457 +++++++++++++++++++++++++++++++++- 5 files changed, 485 insertions(+), 17 deletions(-) diff --git a/src/core/nd_utils.jl b/src/core/nd_utils.jl index 697b4cef..4c2574f6 100644 --- a/src/core/nd_utils.jl +++ b/src/core/nd_utils.jl @@ -114,7 +114,11 @@ Resolve search policy input to canonical N-tuple (broadcast only, no AutoSearch Broadcast + resolve AutoSearch in one step. `query_sample` determines scalar vs vector resolution: - `query_sample::Real` (from `first(query)`) → `Binary()` per axis -- `query_sample::AbstractVector` (from `first(queries)` or `queries`) → `LinearBinary()` per axis +- `query_sample::AbstractVector` → `LinearBinary()` per axis + - SoA batch: pass `first(queries)` where `queries::NTuple{N,AbstractVector}` + - AoS batch: pass `queries` directly where `queries::AbstractVector{<:Tuple}`; + `AbstractVector{<:Tuple} <: AbstractVector` so resolves to `LinearBinary()` — correct + because AoS batches have the same sorted-locality property as SoA batches. - Explicit policies pass through unchanged. """ @inline _resolve_search_nd(s::AbstractSearchPolicy, ::Val{N}) where {N} = ntuple(_ -> s, Val(N)) diff --git a/src/core/search.jl b/src/core/search.jl index afbb6612..869bd04a 100644 --- a/src/core/search.jl +++ b/src/core/search.jl @@ -600,7 +600,9 @@ Bounded linear search within MAX-sized window, then binary fallback. Optimal for monotonic query sequences. # Optimizations over naive implementation: -- No hint clamp: internal hints are always valid (initialized to 1, updated to valid idx) +- Hint clamped once at start: guards against user-provided out-of-range hints (e.g. Ref(0), + stale hints from a different grid). Internal hints (initialized to 1, updated to valid idx) + are already valid, so the clamp is a no-op on the hot path. - No hint write on direct hit: `ix` is unchanged, skip redundant `hint_ref[] = ix` - Single comparison per linear step: direction already determines one bound """ @@ -612,6 +614,7 @@ Optimal for monotonic query sequences. ) where {T<:Real,MAX} ix = hint_ref[] n = length(x) + ix = clamp(ix, 1, n - 1) # guard against user-provided bad hints (e.g. Ref(0), stale) @inbounds begin # Direct hit — most common for sorted/monotonic queries xL = x[ix] diff --git a/src/integral/integrate_api.jl b/src/integral/integrate_api.jl index 423dd03c..d6496251 100644 --- a/src/integral/integrate_api.jl +++ b/src/integral/integrate_api.jl @@ -20,7 +20,8 @@ end y = itp.y z = itp.z Tout = promote_type(Tv, Tg, typeof(x0), typeof(x1)) - searcher = _to_searcher(search, hint) + resolved = _resolve_search(search, x0) + searcher = _to_searcher(resolved, hint) partial = @inline (i, xL, h, a2, b2) -> begin @inbounds _cubic_integral_kernel( @@ -49,7 +50,8 @@ end x = itp.x y = itp.y Tout = promote_type(Tv, Tg, typeof(x0), typeof(x1)) - searcher = _to_searcher(search, hint) + resolved = _resolve_search(search, x0) + searcher = _to_searcher(resolved, hint) partial = @inline (i, xL, h, a2, b2) -> begin @inbounds _linear_integral_kernel( @@ -77,7 +79,8 @@ end ) where {Tg<:AbstractFloat, Tv} x = itp.x Tout = promote_type(Tv, Tg, typeof(x0), typeof(x1)) - searcher = _to_searcher(search, hint) + resolved = _resolve_search(search, x0) + searcher = _to_searcher(resolved, hint) partial = @inline (i, xL, h, a2, b2) -> begin @inbounds _quadratic_integral_kernel( @@ -104,7 +107,8 @@ end hint::Union{Nothing,Base.RefValue{Int}}=nothing ) where {Tg<:AbstractFloat, Tv} Tout = promote_type(Tv, Tg, typeof(x0), typeof(x1)) - searcher = _to_searcher(search, hint) + resolved = _resolve_search(search, x0) + searcher = _to_searcher(resolved, hint) return _integrate_constant_1d_impl(itp.x, itp.y, itp.side, itp.extrap, x0, x1, searcher, Tg, Tout) end @@ -146,7 +150,8 @@ end y = sitp.y z = sitp.z Tout = promote_type(Tv, Tg, typeof(x0), typeof(x1)) - searcher = _to_searcher(search, hint) + resolved = _resolve_search(search, x0) + searcher = _to_searcher(resolved, hint) n = n_series(sitp) results = Vector{Tout}(undef, n) @inbounds for k in 1:n @@ -173,7 +178,8 @@ end x = sitp.x y = sitp.y Tout = promote_type(Tv, Tg, typeof(x0), typeof(x1)) - searcher = _to_searcher(search, hint) + resolved = _resolve_search(search, x0) + searcher = _to_searcher(resolved, hint) n = n_series(sitp) results = Vector{Tout}(undef, n) @inbounds for k in 1:n @@ -199,7 +205,8 @@ end ) where {Tg<:AbstractFloat, Tv} x = sitp.x Tout = promote_type(Tv, Tg, typeof(x0), typeof(x1)) - searcher = _to_searcher(search, hint) + resolved = _resolve_search(search, x0) + searcher = _to_searcher(resolved, hint) n = n_series(sitp) results = Vector{Tout}(undef, n) @inbounds for k in 1:n @@ -224,7 +231,8 @@ end hint::Union{Nothing,Base.RefValue{Int}}=nothing ) where {Tg<:AbstractFloat, Tv} Tout = promote_type(Tv, Tg, typeof(x0), typeof(x1)) - searcher = _to_searcher(search, hint) + resolved = _resolve_search(search, x0) + searcher = _to_searcher(resolved, hint) return _integrate_constant_series_1d(sitp.x, sitp.y, sitp.side, sitp.extrap, x0, x1, searcher, Tg, Tout) end diff --git a/src/vector_calculus.jl b/src/vector_calculus.jl index b2a21895..3567bdbf 100644 --- a/src/vector_calculus.jl +++ b/src/vector_calculus.jl @@ -50,7 +50,7 @@ See also: [`gradient!`](@ref), [`hessian`](@ref), [`laplacian`](@ref) end for i in 1:N] return quote - search = _resolve_search_nd(itp.searches, Val($N)) + search = _resolve_search_nd(itp.searches, Val($N), first(query)) cell = _locate_cell(itp, query, search, hint) return tuple($(deriv_calls...)) end @@ -105,7 +105,7 @@ See also: [`gradient`](@ref), [`hessian!`](@ref) @boundscheck length(G) >= $N || throw(DimensionMismatch( "gradient output vector must have at least $($N) elements, got $(length(G))" )) - search = _resolve_search_nd(itp.searches, Val($N)) + search = _resolve_search_nd(itp.searches, Val($N), first(query)) cell = _locate_cell(itp, query, search, hint) @inbounds begin $(stmts...) @@ -180,7 +180,7 @@ See also: [`gradient`](@ref), [`hessian!`](@ref), [`laplacian`](@ref) return quote Tq = promote_type(eltype(query), $Tg, $Tv) H = Matrix{Tq}(undef, $N, $N) - search = _resolve_search_nd(itp.searches, Val($N)) + search = _resolve_search_nd(itp.searches, Val($N), first(query)) cell = _locate_cell(itp, query, search, hint) @inbounds begin $(stmts...) @@ -251,7 +251,7 @@ See also: [`hessian`](@ref), [`gradient!`](@ref) @boundscheck size(H) == ($N, $N) || throw(DimensionMismatch( "Hessian output matrix must be $($N)×$($N), got $(size(H))" )) - search = _resolve_search_nd(itp.searches, Val($N)) + search = _resolve_search_nd(itp.searches, Val($N), first(query)) cell = _locate_cell(itp, query, search, hint) @inbounds begin $(stmts...) @@ -316,7 +316,7 @@ See also: [`gradient`](@ref), [`hessian`](@ref) end for i in 1:N] return quote - search = _resolve_search_nd(itp.searches, Val($N)) + search = _resolve_search_nd(itp.searches, Val($N), first(query)) cell = _locate_cell(itp, query, search, hint) return +($(deriv_calls...)) end diff --git a/test/test_search.jl b/test/test_search.jl index dc42b338..fd5ae957 100644 --- a/test/test_search.jl +++ b/test/test_search.jl @@ -1,8 +1,9 @@ using Test using FastInterpolations using FastInterpolations: search_interval, _search_binary, _search_direct, _search_interval, - Searcher, Binary, HintedBinary, Linear, LinearBinary, - NoHint, RefHint, DEFAULT_SEARCHER, ScalarSpacing, _create_spacing, _to_searcher + Searcher, Binary, HintedBinary, Linear, LinearBinary, AutoSearch, + NoHint, RefHint, DEFAULT_SEARCHER, ScalarSpacing, _create_spacing, _to_searcher, + _resolve_search @testset "Search Module" begin @@ -724,6 +725,22 @@ using FastInterpolations: search_interval, _search_binary, _search_direct, _sear end end + @testset "LinearBinary: out-of-range hint is clamped safely" begin + x = collect(range(0.0, 1.0, 101)) + y = x .^ 2 + itp = linear_interp(x, y; search=LinearBinary()) + + # Ref(0) — below valid range [1, n-1] + bad_hint = Ref(0) + @test isfinite(itp(0.5; hint=bad_hint)) + @test bad_hint[] >= 1 # hint was clamped then updated to valid index + + # Ref(n) — above valid range [1, n-1] + big_hint = Ref(200) + @test isfinite(itp(0.5; hint=big_hint)) + @test big_hint[] <= 100 # hint was clamped then updated + end + @testset "Integrated test" begin x = collect(range(0.0, 1.0, 101)) y = x.^3 @@ -1097,4 +1114,440 @@ using FastInterpolations: search_interval, _search_binary, _search_direct, _sear end end + # ============================================================================ + # AutoSearch Resolution Tests + # ============================================================================ + # Verifies that AutoSearch resolves to the correct concrete policy at every + # API layer: unit dispatch, interpolant callable, oneshot, and series. + + @testset "AutoSearch Resolution" begin + + # ======================================== + # Unit Tests: _resolve_search dispatch + # ======================================== + @testset "_resolve_search dispatch" begin + auto = AutoSearch() + + # 1D scalar → Binary + @test _resolve_search(auto, 0.5) isa Binary + @test _resolve_search(auto, 1) isa Binary # Int <: Real + @test _resolve_search(auto, Float32(0.5)) isa Binary + + # 1D vector → LinearBinary + @test _resolve_search(auto, [0.1, 0.5]) isa LinearBinary + @test _resolve_search(auto, 0.0:0.1:1.0) isa LinearBinary # Range <: AbstractVector + @test _resolve_search(auto, Float32[1.0, 2.0]) isa LinearBinary + + # ND scalar (Tuple of Reals) → Binary + @test _resolve_search(auto, (0.5, 0.3)) isa Binary + @test _resolve_search(auto, (0.1, 0.2, 0.3)) isa Binary + + # ND vector (Tuple of Vectors) → LinearBinary + @test _resolve_search(auto, ([0.1, 0.2], [0.3, 0.4])) isa LinearBinary + @test _resolve_search(auto, (0.0:0.1:1.0, [0.5])) isa LinearBinary + + # Passthrough: explicit policies are returned unchanged + @test _resolve_search(Binary(), 0.5) === Binary() + @test _resolve_search(Binary(), [0.1]) === Binary() + @test _resolve_search(LinearBinary(), 0.5) isa LinearBinary + @test _resolve_search(LinearBinary(), [0.1]) isa LinearBinary + @test _resolve_search(HintedBinary(), 0.5) === HintedBinary() + @test _resolve_search(Linear(), [0.1]) === Linear() + + # Tuple of policies: each resolved independently + mixed = (AutoSearch(), Binary(), AutoSearch()) + resolved_scalar = _resolve_search(mixed, 0.5) + @test resolved_scalar == (Binary(), Binary(), Binary()) + resolved_vec = _resolve_search(mixed, [0.1, 0.2]) + @test resolved_vec[1] isa LinearBinary + @test resolved_vec[2] === Binary() + @test resolved_vec[3] isa LinearBinary + + # Searcher passthrough: pre-built Searcher skips resolution + s = Searcher{Binary,NoHint}(NoHint()) + @test _resolve_search(s, 0.5) === s + @test _resolve_search(s, [0.1]) === s + end + + # ======================================== + # Interpolant API: AutoSearch resolves per query type + # ======================================== + @testset "Interpolant API" begin + x = collect(range(0.0, 1.0, 201)) + y = sin.(2π .* x) + + @testset "Linear" begin + itp = linear_interp(x, y) + @test itp.search_policy isa AutoSearch + + # Scalar call → correct result (AutoSearch → Binary internally) + val = itp(0.25) + @test val ≈ sin(2π * 0.25) atol=1e-3 + + # Vector call → correct results (AutoSearch → LinearBinary internally) + xq = [0.1, 0.3, 0.7] + vals = itp(xq) + @test length(vals) == 3 + @test all(vals .≈ sin.(2π .* xq) .|| abs.(vals .- sin.(2π .* xq)) .< 0.02) + + # In-place vector call + out = zeros(3) + itp(out, xq) + @test out ≈ vals + end + + @testset "Cubic" begin + itp = cubic_interp(x, y) + @test itp.search_policy isa AutoSearch + + val = itp(0.25) + @test val ≈ sin(2π * 0.25) atol=1e-5 + + xq = [0.1, 0.3, 0.7] + vals = itp(xq) + @test length(vals) == 3 + for i in 1:3 + @test vals[i] ≈ sin(2π * xq[i]) atol=1e-4 + end + end + + @testset "Quadratic" begin + itp = quadratic_interp(x, y) + @test itp.search_policy isa AutoSearch + + val = itp(0.25) + @test val ≈ sin(2π * 0.25) atol=1e-3 + + xq = [0.1, 0.3, 0.7] + vals = itp(xq) + @test length(vals) == 3 + for i in 1:3 + @test vals[i] ≈ sin(2π * xq[i]) atol=1e-2 + end + end + + @testset "Constant" begin + itp = constant_interp(x, y) + @test itp.search_policy isa AutoSearch + + # Constant interpolation snaps to nearest grid point + val = itp(0.25) + @test isfinite(val) + + xq = [0.1, 0.3, 0.7] + vals = itp(xq) + @test length(vals) == 3 + @test all(isfinite, vals) + end + end + + # ======================================== + # Interpolant API: hint behavior verifies resolved policy + # ======================================== + # LinearBinary with hint tracks sequentially through sorted queries. + # Binary (no hint) doesn't update the hint. + # AutoSearch + vector → LinearBinary → hint should track. + # AutoSearch + scalar → Binary → hint unchanged (no hint used). + @testset "Interpolant hint tracking reveals resolved policy" begin + x = collect(range(0.0, 1.0, 1001)) + y = x .^ 2 + + itp = linear_interp(x, y) # default AutoSearch + + # Vector call with hint: AutoSearch → LinearBinary → hint tracks + hint = Ref(1) + xq_sorted = collect(range(0.1, 0.9, 100)) + itp(zeros(100), xq_sorted; hint=hint) + # LinearBinary leaves hint near last query (~900) + @test hint[] >= 850 + + # Scalar call without hint: AutoSearch → Binary → correct value + # y = x^2, so itp(0.5) ≈ 0.25 + @test itp(0.5) ≈ 0.25 atol=1e-6 + + # Scalar call with hint: Binary+hint auto-upgrades to HintedBinary → hint updates + hint_scalar = Ref(1) + itp(0.5; hint=hint_scalar) + @test hint_scalar[] >= 490 && hint_scalar[] <= 510 # hint moved to ~500 + end + + # ======================================== + # Oneshot API: AutoSearch resolves per query type + # ======================================== + @testset "Oneshot API" begin + x = collect(range(0.0, 1.0, 201)) + y = sin.(2π .* x) + + @testset "linear_interp scalar oneshot (AutoSearch default)" begin + val = linear_interp(x, y, 0.25) + @test val ≈ sin(2π * 0.25) atol=1e-3 + end + + @testset "linear_interp vector oneshot (AutoSearch default)" begin + xq = collect(range(0.1, 0.9, 50)) + vals = linear_interp(x, y, xq) + @test length(vals) == 50 + for i in 1:50 + @test vals[i] ≈ sin(2π * xq[i]) atol=1e-2 + end + end + + @testset "linear_interp! vector oneshot (AutoSearch default)" begin + xq = collect(range(0.1, 0.9, 50)) + out = zeros(50) + linear_interp!(out, x, y, xq) + for i in 1:50 + @test out[i] ≈ sin(2π * xq[i]) atol=1e-2 + end + end + + @testset "cubic_interp scalar oneshot" begin + val = cubic_interp(x, y, 0.25) + @test val ≈ sin(2π * 0.25) atol=1e-5 + end + + @testset "cubic_interp! vector oneshot" begin + xq = collect(range(0.1, 0.9, 50)) + out = zeros(50) + cubic_interp!(out, x, y, xq) + for i in 1:50 + @test out[i] ≈ sin(2π * xq[i]) atol=1e-3 + end + end + + @testset "quadratic_interp scalar oneshot" begin + val = quadratic_interp(x, y, 0.25) + @test val ≈ sin(2π * 0.25) atol=1e-3 + end + + @testset "quadratic_interp! vector oneshot" begin + xq = collect(range(0.1, 0.9, 50)) + out = zeros(50) + quadratic_interp!(out, x, y, xq) + for i in 1:50 + @test out[i] ≈ sin(2π * xq[i]) atol=1e-2 + end + end + + @testset "constant_interp scalar oneshot" begin + val = constant_interp(x, y, 0.25) + @test isfinite(val) + end + + @testset "constant_interp! vector oneshot" begin + xq = collect(range(0.1, 0.9, 50)) + out = zeros(50) + constant_interp!(out, x, y, xq) + @test all(isfinite, out) + end + end + + # ======================================== + # Series Interpolant API: AutoSearch resolves per query type + # ======================================== + @testset "Series API" begin + x = collect(range(0.0, 1.0, 201)) + y1 = sin.(2π .* x) + y2 = cos.(2π .* x) + + @testset "LinearSeries" begin + sitp = linear_interp(x, [y1, y2]) + @test sitp.search_policy isa AutoSearch + + # Scalar call → AutoSearch → Binary + vals = sitp(0.25) + @test length(vals) == 2 + @test vals[1] ≈ sin(2π * 0.25) atol=1e-3 + @test vals[2] ≈ cos(2π * 0.25) atol=1e-3 + + # Vector call → AutoSearch → LinearBinary + xq = [0.1, 0.3, 0.7] + result = sitp(xq) + @test length(result) == 2 # 2 series + @test length(result[1]) == 3 # 3 query points + @test result[1][1] ≈ sin(2π * 0.1) atol=1e-2 + @test result[2][1] ≈ cos(2π * 0.1) atol=1e-2 + end + + @testset "CubicSeries" begin + sitp = cubic_interp(x, [y1, y2]) + @test sitp.search_policy isa AutoSearch + + vals = sitp(0.25) + @test vals[1] ≈ sin(2π * 0.25) atol=1e-5 + @test vals[2] ≈ cos(2π * 0.25) atol=1e-5 + + xq = [0.1, 0.3, 0.7] + result = sitp(xq) + @test length(result) == 2 + for i in 1:3 + @test result[1][i] ≈ sin(2π * xq[i]) atol=1e-3 + end + end + + @testset "QuadraticSeries" begin + sitp = quadratic_interp(x, [y1, y2]) + @test sitp.search_policy isa AutoSearch + + vals = sitp(0.25) + @test vals[1] ≈ sin(2π * 0.25) atol=1e-3 + + xq = [0.1, 0.3, 0.7] + result = sitp(xq) + @test length(result) == 2 + @test length(result[1]) == 3 + end + + @testset "ConstantSeries" begin + sitp = constant_interp(x, [y1, y2]) + @test sitp.search_policy isa AutoSearch + + vals = sitp(0.25) + @test length(vals) == 2 + @test all(isfinite, vals) + + xq = [0.1, 0.3, 0.7] + result = sitp(xq) + @test length(result) == 2 + @test length(result[1]) == 3 + end + + @testset "Series hint tracking (vector → LinearBinary)" begin + sitp = linear_interp(x, [y1, y2]) + hint = Ref(1) + xq_sorted = collect(range(0.1, 0.9, 100)) + + # Vector call: AutoSearch → LinearBinary → hint tracks + sitp(xq_sorted; hint=hint) + @test hint[] >= 170 # near end of sorted queries + end + end + + # ======================================== + # Explicit User Override: policy honored across all APIs + # ======================================== + @testset "Explicit user override" begin + x = collect(range(0.0, 1.0, 201)) + y = sin.(2π .* x) + y1 = sin.(2π .* x) + y2 = cos.(2π .* x) + + @testset "Interpolant: explicit Binary on vector call" begin + itp = linear_interp(x, y; search=Binary()) + @test itp.search_policy isa Binary + + # Vector call with Binary baked in → still works correctly + xq = [0.1, 0.5, 0.9] + vals = itp(xq) + for i in 1:3 + @test vals[i] ≈ sin(2π * xq[i]) atol=1e-2 + end + end + + @testset "Interpolant: explicit LinearBinary on scalar call" begin + itp = linear_interp(x, y; search=LinearBinary()) + @test itp.search_policy isa LinearBinary + + # Scalar call with LinearBinary baked in → still works + val = itp(0.25) + @test val ≈ sin(2π * 0.25) atol=1e-3 + end + + @testset "Interpolant: call-site override trumps stored policy" begin + itp = linear_interp(x, y) # AutoSearch stored + @test itp.search_policy isa AutoSearch + + # Override with explicit Binary at call-site + val = itp(0.25; search=Binary()) + @test val ≈ sin(2π * 0.25) atol=1e-3 + + # Override with explicit LinearBinary at call-site + vals = itp([0.1, 0.5]; search=LinearBinary()) + @test length(vals) == 2 + end + + @testset "Oneshot: explicit policy" begin + xq = collect(range(0.1, 0.9, 50)) + out = zeros(50) + + # Explicit Binary + linear_interp!(out, x, y, xq; search=Binary()) + @test all(isfinite, out) + + # Explicit LinearBinary + linear_interp!(out, x, y, xq; search=LinearBinary()) + @test all(isfinite, out) + + # Scalar oneshot with explicit policy + val = linear_interp(x, y, 0.25; search=LinearBinary()) + @test val ≈ sin(2π * 0.25) atol=1e-3 + end + + @testset "Series: explicit policy" begin + sitp = linear_interp(x, [y1, y2]; search=Binary()) + @test sitp.search_policy isa Binary + + # Scalar call + vals = sitp(0.25) + @test vals[1] ≈ sin(2π * 0.25) atol=1e-3 + + # Vector call + xq = [0.1, 0.3] + result = sitp(xq) + @test length(result) == 2 + @test length(result[1]) == 2 + end + + @testset "Series: call-site override" begin + sitp = linear_interp(x, [y1, y2]) # AutoSearch + @test sitp.search_policy isa AutoSearch + + # Override at call-site + vals = sitp(0.25; search=Binary()) + @test length(vals) == 2 + + result = sitp([0.1, 0.3]; search=LinearBinary()) + @test length(result) == 2 + end + end + + # ======================================== + # ND vector calculus: gradient/hessian/laplacian resolve AutoSearch + # ======================================== + @testset "gradient/hessian/laplacian resolve AutoSearch" begin + x = collect(range(0.0, 1.0, 21)) + y = collect(range(0.0, 1.0, 21)) + data = [xi^2 + yi^2 for xi in x, yi in y] + itp = cubic_interp((x, y), data) # AutoSearch stored on each axis + @test itp.searches[1] isa AutoSearch + @test itp.searches[2] isa AutoSearch + + q = (0.5, 0.5) + + # gradient: ∇(x²+y²) = (2x, 2y) + g = gradient(itp, q) + @test all(isfinite, g) + @test g[1] ≈ 2 * 0.5 atol=0.05 + @test g[2] ≈ 2 * 0.5 atol=0.05 + + # gradient! (in-place) + G = zeros(2) + gradient!(G, itp, q) + @test G ≈ collect(g) + + # hessian: H(x²+y²) = [[2,0],[0,2]] + h = hessian(itp, q) + @test all(isfinite, h) + @test h[1, 1] ≈ 2.0 atol=0.1 + @test h[2, 2] ≈ 2.0 atol=0.1 + @test abs(h[1, 2]) < 0.05 # off-diagonal ≈ 0 + + # laplacian: ∇²(x²+y²) = 4 + lap = laplacian(itp, q) + @test isfinite(lap) + @test lap ≈ 4.0 atol=0.1 + end + + end # @testset "AutoSearch Resolution" + end # @testset "Search Module" From 20981baa143a844a75e2108d7b99d099194984de Mon Sep 17 00:00:00 2001 From: Min-Gu Yoo Date: Tue, 24 Feb 2026 15:42:55 -0800 Subject: [PATCH 13/16] Fix AutoSearch review issues: ND integrate preamble + clarify comments + tests MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - integrate_common_nd.jl: _integrate_nd_preamble now uses 3-arg _resolve_search_nd(search, Val(N), first(lo)), routing AutoSearch through the intended resolution path instead of the _to_searcher safety-net fallback - search.jl: fix misleading dispatch ordering comment (Tuple{Vararg{Real}} never existed; the fallback is bare ::Tuple); add n>=2 precondition note on the clamp line - vector_calculus.jl: add inline comment on all 5 @generated sites explaining why first(query)::Real is correct (scalar-point-only API) - 9 ND eval/oneshot files: add inline comment at each AoS dispatch site documenting the AbstractVector{<:Tuple} <: AbstractVector trick - test_search.jl: document expected indices in hint-tracking assertions; widen series bound from >=170 to >=160 (margin 11→21); improve clamp test comment to note valid range --- src/constant/nd/constant_nd_eval.jl | 2 +- src/constant/nd/constant_nd_oneshot.jl | 4 ++-- src/core/search.jl | 4 +++- src/cubic/nd/cubic_nd_eval.jl | 2 +- src/cubic/nd/cubic_nd_oneshot.jl | 2 +- src/integral/integrate_common_nd.jl | 2 +- src/linear/nd/linear_nd_eval.jl | 2 +- src/linear/nd/linear_nd_oneshot.jl | 2 +- src/quadratic/nd/quadratic_nd_eval.jl | 2 +- src/quadratic/nd/quadratic_nd_oneshot.jl | 2 +- src/vector_calculus.jl | 10 +++++----- test/test_search.jl | 7 ++++--- 12 files changed, 22 insertions(+), 19 deletions(-) diff --git a/src/constant/nd/constant_nd_eval.jl b/src/constant/nd/constant_nd_eval.jl index 9468a848..976fed9e 100644 --- a/src/constant/nd/constant_nd_eval.jl +++ b/src/constant/nd/constant_nd_eval.jl @@ -75,7 +75,7 @@ function (itp::ConstantInterpolantND{Tg,Tv,N})( "output length $(length(output)) must match query length $n_queries" )) ops = _resolve_deriv_nd(deriv, Val(N)) - search_tuple = _resolve_search_nd(search, Val(N), queries) + search_tuple = _resolve_search_nd(search, Val(N), queries) # AoS: AbstractVector{<:Tuple} <: AbstractVector → LinearBinary if _has_any_derivative(ops, Val(N)) fill!(output, zero(eltype(output))) return output diff --git a/src/constant/nd/constant_nd_oneshot.jl b/src/constant/nd/constant_nd_oneshot.jl index e613b201..457ad9d8 100644 --- a/src/constant/nd/constant_nd_oneshot.jl +++ b/src/constant/nd/constant_nd_oneshot.jl @@ -230,7 +230,7 @@ function constant_interp( _validate_nd_grids(grids_typed, data) sides = _resolve_side_nd(side, Val(N)) - searches = _resolve_search_nd(search, Val(N), queries) + searches = _resolve_search_nd(search, Val(N), queries) # AoS: AbstractVector{<:Tuple} <: AbstractVector → LinearBinary extraps_val = _resolve_extrap_nd(extrap, nothing, Val(N)) return _constant_interp_nd_oneshot_aos( @@ -302,7 +302,7 @@ function constant_interp!( _validate_nd_grids(grids_typed, data) sides = _resolve_side_nd(side, Val(N)) - searches = _resolve_search_nd(search, Val(N), queries) + searches = _resolve_search_nd(search, Val(N), queries) # AoS: AbstractVector{<:Tuple} <: AbstractVector → LinearBinary extraps_val = _resolve_extrap_nd(extrap, nothing, Val(N)) return _constant_interp_nd_oneshot_aos!( diff --git a/src/core/search.jl b/src/core/search.jl index 869bd04a..299ef6bc 100644 --- a/src/core/search.jl +++ b/src/core/search.jl @@ -229,7 +229,8 @@ struct AutoSearch <: AbstractSearchPolicy end @inline _resolve_search(::AutoSearch, ::AbstractVector) = LinearBinary() # ND SoA batch: tuple of vectors → LinearBinary per axis -# NOTE: must be above Tuple{Vararg{Real}} to avoid ambiguity (AbstractVector <: Any) +# NOTE: must precede the bare ::Tuple fallback — Tuple{Vararg{AbstractVector}} <: Tuple, +# so Julia's specificity rules handle ordering correctly, but explicit ordering avoids confusion. @inline _resolve_search(::AutoSearch, ::Tuple{Vararg{AbstractVector}}) = LinearBinary() # ND scalar: tuple of reals → Binary per axis @@ -615,6 +616,7 @@ Optimal for monotonic query sequences. ix = hint_ref[] n = length(x) ix = clamp(ix, 1, n - 1) # guard against user-provided bad hints (e.g. Ref(0), stale) + # Precondition: n >= 2 (enforced by all interpolant constructors) @inbounds begin # Direct hit — most common for sorted/monotonic queries xL = x[ix] diff --git a/src/cubic/nd/cubic_nd_eval.jl b/src/cubic/nd/cubic_nd_eval.jl index aa34b26b..3048022e 100644 --- a/src/cubic/nd/cubic_nd_eval.jl +++ b/src/cubic/nd/cubic_nd_eval.jl @@ -95,7 +95,7 @@ function (itp::CubicInterpolantND{Tg, Tv, N})( "output length $(length(output)) must match query length $n_queries" )) ops = _resolve_deriv_nd(deriv, Val(N)) - search_tuple = _resolve_search_nd(search, Val(N), queries) + search_tuple = _resolve_search_nd(search, Val(N), queries) # AoS: AbstractVector{<:Tuple} <: AbstractVector → LinearBinary _batch_nd_aos!(output, itp, queries, ops, search_tuple, hint) return output end diff --git a/src/cubic/nd/cubic_nd_oneshot.jl b/src/cubic/nd/cubic_nd_oneshot.jl index 4828af89..80573edd 100644 --- a/src/cubic/nd/cubic_nd_oneshot.jl +++ b/src/cubic/nd/cubic_nd_oneshot.jl @@ -299,7 +299,7 @@ function cubic_interp!( _validate_nd_grids(grids_typed, data) bcs = _resolve_bcs_nd(bc, Val(N)) - searches = _resolve_search_nd(search, Val(N), queries) + searches = _resolve_search_nd(search, Val(N), queries) # AoS: AbstractVector{<:Tuple} <: AbstractVector → LinearBinary _validate_nd_bcs!(grids_typed, bcs, data, Val(N)) diff --git a/src/integral/integrate_common_nd.jl b/src/integral/integrate_common_nd.jl index 053d3d69..79d6a21a 100644 --- a/src/integral/integrate_common_nd.jl +++ b/src/integral/integrate_common_nd.jl @@ -55,7 +55,7 @@ end _check_nd_integrate_domain(grids[d], lo2[d], extraps[d]) _check_nd_integrate_domain(grids[d], hi2[d], extraps[d]) end - search_tuple = _resolve_search_nd(search, Val(N)) + search_tuple = _resolve_search_nd(search, Val(N), first(lo)) # lo[1]::Real → Binary/axis idx_lo, idx_hi = _nd_cell_ranges(grids, spacings, lo2, hi2, search_tuple, hint) return (sign, lo2, hi2, idx_lo, idx_hi) end diff --git a/src/linear/nd/linear_nd_eval.jl b/src/linear/nd/linear_nd_eval.jl index 2145adde..a01aafca 100644 --- a/src/linear/nd/linear_nd_eval.jl +++ b/src/linear/nd/linear_nd_eval.jl @@ -79,7 +79,7 @@ function (itp::LinearInterpolantND{Tg,Tv,N})( "output length $(length(output)) must match query length $n_queries" )) ops = _resolve_deriv_nd(deriv, Val(N)) - search_tuple = _resolve_search_nd(search, Val(N), queries) + search_tuple = _resolve_search_nd(search, Val(N), queries) # AoS: AbstractVector{<:Tuple} <: AbstractVector → LinearBinary if _has_second_or_higher_derivative(ops, Val(N)) fill!(output, zero(eltype(output))) return output diff --git a/src/linear/nd/linear_nd_oneshot.jl b/src/linear/nd/linear_nd_oneshot.jl index 3f14ca72..9e68c05d 100644 --- a/src/linear/nd/linear_nd_oneshot.jl +++ b/src/linear/nd/linear_nd_oneshot.jl @@ -223,7 +223,7 @@ function linear_interp!( grids_typed = _convert_grids_typed(grids, Tg) _validate_nd_grids(grids_typed, data) - searches = _resolve_search_nd(search, Val(N), queries) + searches = _resolve_search_nd(search, Val(N), queries) # AoS: AbstractVector{<:Tuple} <: AbstractVector → LinearBinary extraps_val = _resolve_extrap_nd(extrap, nothing, Val(N)) ops = _resolve_deriv_nd(deriv, Val(N)) diff --git a/src/quadratic/nd/quadratic_nd_eval.jl b/src/quadratic/nd/quadratic_nd_eval.jl index 5bfb1167..692b1d06 100644 --- a/src/quadratic/nd/quadratic_nd_eval.jl +++ b/src/quadratic/nd/quadratic_nd_eval.jl @@ -120,7 +120,7 @@ function (itp::QuadraticInterpolantND{Tg, Tv, N})( "output length $(length(output)) must match query length $n_queries" )) ops = _resolve_deriv_nd(deriv, Val(N)) - search_tuple = _resolve_search_nd(search, Val(N), queries) + search_tuple = _resolve_search_nd(search, Val(N), queries) # AoS: AbstractVector{<:Tuple} <: AbstractVector → LinearBinary _batch_nd_aos!(output, itp, queries, ops, search_tuple, hint) return output end diff --git a/src/quadratic/nd/quadratic_nd_oneshot.jl b/src/quadratic/nd/quadratic_nd_oneshot.jl index cc8fd07f..03e62179 100644 --- a/src/quadratic/nd/quadratic_nd_oneshot.jl +++ b/src/quadratic/nd/quadratic_nd_oneshot.jl @@ -263,7 +263,7 @@ function quadratic_interp!( _validate_nd_grids(grids_typed, data) bcs = _resolve_bcs_nd_quadratic(bc, Val(N)) - searches = _resolve_search_nd(search, Val(N), queries) + searches = _resolve_search_nd(search, Val(N), queries) # AoS: AbstractVector{<:Tuple} <: AbstractVector → LinearBinary extraps_val = _resolve_extrap_nd(extrap, bcs, Val(N)) ops = _resolve_deriv_nd(deriv, Val(N)) diff --git a/src/vector_calculus.jl b/src/vector_calculus.jl index 3567bdbf..3c9f7f37 100644 --- a/src/vector_calculus.jl +++ b/src/vector_calculus.jl @@ -50,7 +50,7 @@ See also: [`gradient!`](@ref), [`hessian`](@ref), [`laplacian`](@ref) end for i in 1:N] return quote - search = _resolve_search_nd(itp.searches, Val($N), first(query)) + search = _resolve_search_nd(itp.searches, Val($N), first(query)) # Real sample → Binary/axis (scalar ND point) cell = _locate_cell(itp, query, search, hint) return tuple($(deriv_calls...)) end @@ -105,7 +105,7 @@ See also: [`gradient`](@ref), [`hessian!`](@ref) @boundscheck length(G) >= $N || throw(DimensionMismatch( "gradient output vector must have at least $($N) elements, got $(length(G))" )) - search = _resolve_search_nd(itp.searches, Val($N), first(query)) + search = _resolve_search_nd(itp.searches, Val($N), first(query)) # Real sample → Binary/axis (scalar ND point) cell = _locate_cell(itp, query, search, hint) @inbounds begin $(stmts...) @@ -180,7 +180,7 @@ See also: [`gradient`](@ref), [`hessian!`](@ref), [`laplacian`](@ref) return quote Tq = promote_type(eltype(query), $Tg, $Tv) H = Matrix{Tq}(undef, $N, $N) - search = _resolve_search_nd(itp.searches, Val($N), first(query)) + search = _resolve_search_nd(itp.searches, Val($N), first(query)) # Real sample → Binary/axis (scalar ND point) cell = _locate_cell(itp, query, search, hint) @inbounds begin $(stmts...) @@ -251,7 +251,7 @@ See also: [`hessian`](@ref), [`gradient!`](@ref) @boundscheck size(H) == ($N, $N) || throw(DimensionMismatch( "Hessian output matrix must be $($N)×$($N), got $(size(H))" )) - search = _resolve_search_nd(itp.searches, Val($N), first(query)) + search = _resolve_search_nd(itp.searches, Val($N), first(query)) # Real sample → Binary/axis (scalar ND point) cell = _locate_cell(itp, query, search, hint) @inbounds begin $(stmts...) @@ -316,7 +316,7 @@ See also: [`gradient`](@ref), [`hessian`](@ref) end for i in 1:N] return quote - search = _resolve_search_nd(itp.searches, Val($N), first(query)) + search = _resolve_search_nd(itp.searches, Val($N), first(query)) # Real sample → Binary/axis (scalar ND point) cell = _locate_cell(itp, query, search, hint) return +($(deriv_calls...)) end diff --git a/test/test_search.jl b/test/test_search.jl index fd5ae957..73afdd10 100644 --- a/test/test_search.jl +++ b/test/test_search.jl @@ -738,7 +738,7 @@ using FastInterpolations: search_interval, _search_binary, _search_direct, _sear # Ref(n) — above valid range [1, n-1] big_hint = Ref(200) @test isfinite(itp(0.5; hint=big_hint)) - @test big_hint[] <= 100 # hint was clamped then updated + @test big_hint[] <= 100 # valid range [1, n-1]=[1,100]; query 0.5 → idx ≈ 50 end @testset "Integrated test" begin @@ -1258,7 +1258,7 @@ using FastInterpolations: search_interval, _search_binary, _search_direct, _sear hint = Ref(1) xq_sorted = collect(range(0.1, 0.9, 100)) itp(zeros(100), xq_sorted; hint=hint) - # LinearBinary leaves hint near last query (~900) + # 1001-point grid, last query at x=0.9 → idx ≈ 901; margin of 50 @test hint[] >= 850 # Scalar call without hint: AutoSearch → Binary → correct value @@ -1419,7 +1419,8 @@ using FastInterpolations: search_interval, _search_binary, _search_direct, _sear # Vector call: AutoSearch → LinearBinary → hint tracks sitp(xq_sorted; hint=hint) - @test hint[] >= 170 # near end of sorted queries + # 201-point grid, last query at x=0.9 → idx ≈ 181; margin of 21 + @test hint[] >= 160 end end From b4480af42fcec6db3d2b278bad778de24e01590e Mon Sep 17 00:00:00 2001 From: Min-Gu Yoo Date: Tue, 24 Feb 2026 16:03:20 -0800 Subject: [PATCH 14/16] Refactor _resolve_search_nd call sites: drop first() extractions MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Pass query containers directly to _resolve_search_nd instead of extracting a sample element with first(). Julia dispatch on the container type already gives the correct resolution: NTuple{N,Real} → Tuple arm → Binary/axis NTuple{N,AbstractVector} → Tuple{Vararg{...}} arm → LinearBinary/axis AbstractVector{<:Tuple} → AbstractVector arm → LinearBinary/axis 24 sites across 11 files updated. nd_utils.jl docstring updated to document the direct-container passing convention. --- src/constant/nd/constant_nd_eval.jl | 4 ++-- src/constant/nd/constant_nd_oneshot.jl | 6 +++--- src/core/nd_utils.jl | 11 ++++------- src/cubic/nd/cubic_nd_eval.jl | 4 ++-- src/cubic/nd/cubic_nd_oneshot.jl | 4 ++-- src/integral/integrate_common_nd.jl | 2 +- src/linear/nd/linear_nd_eval.jl | 4 ++-- src/linear/nd/linear_nd_oneshot.jl | 4 ++-- src/quadratic/nd/quadratic_nd_eval.jl | 4 ++-- src/quadratic/nd/quadratic_nd_oneshot.jl | 4 ++-- src/vector_calculus.jl | 10 +++++----- 11 files changed, 27 insertions(+), 30 deletions(-) diff --git a/src/constant/nd/constant_nd_eval.jl b/src/constant/nd/constant_nd_eval.jl index 976fed9e..0c795a67 100644 --- a/src/constant/nd/constant_nd_eval.jl +++ b/src/constant/nd/constant_nd_eval.jl @@ -17,7 +17,7 @@ hint::Union{Nothing, NTuple{N, Base.RefValue{Int}}} = nothing ) where {Tg, Tv, N} ops = _resolve_deriv_nd(deriv, Val(N)) - search_tuple = _resolve_search_nd(search, Val(N), first(query)) + search_tuple = _resolve_search_nd(search, Val(N), query) # NTuple{N,Real} <: Tuple → Binary/axis return _eval_constant_nd(itp, query, ops, search_tuple, hint) end @@ -48,7 +48,7 @@ function (itp::ConstantInterpolantND{Tg,Tv,N})( )) end ops = _resolve_deriv_nd(deriv, Val(N)) - search_tuple = _resolve_search_nd(search, Val(N), first(queries)) + search_tuple = _resolve_search_nd(search, Val(N), queries) # NTuple{N,AbstractVector} <: Tuple{Vararg{AbstractVector}} → LinearBinary/axis if _has_any_derivative(ops, Val(N)) fill!(output, zero(eltype(output))) return output diff --git a/src/constant/nd/constant_nd_oneshot.jl b/src/constant/nd/constant_nd_oneshot.jl index 457ad9d8..9ea1697e 100644 --- a/src/constant/nd/constant_nd_oneshot.jl +++ b/src/constant/nd/constant_nd_oneshot.jl @@ -164,7 +164,7 @@ function constant_interp( _validate_nd_grids(grids_typed, data) sides = _resolve_side_nd(side, Val(N)) - searches = _resolve_search_nd(search, Val(N), first(query)) + searches = _resolve_search_nd(search, Val(N), query) # NTuple{N,Real} <: Tuple → Binary/axis extraps_val = _resolve_extrap_nd(extrap, nothing, Val(N)) return _constant_interp_nd_oneshot( @@ -197,7 +197,7 @@ function constant_interp( _validate_nd_grids(grids_typed, data) sides = _resolve_side_nd(side, Val(N)) - searches = _resolve_search_nd(search, Val(N), first(queries)) + searches = _resolve_search_nd(search, Val(N), queries) # NTuple{N,AbstractVector} <: Tuple{Vararg{AbstractVector}} → LinearBinary/axis extraps_val = _resolve_extrap_nd(extrap, nothing, Val(N)) return _constant_interp_nd_oneshot_soa( @@ -268,7 +268,7 @@ function constant_interp!( _validate_nd_grids(grids_typed, data) sides = _resolve_side_nd(side, Val(N)) - searches = _resolve_search_nd(search, Val(N), first(queries)) + searches = _resolve_search_nd(search, Val(N), queries) # NTuple{N,AbstractVector} <: Tuple{Vararg{AbstractVector}} → LinearBinary/axis extraps_val = _resolve_extrap_nd(extrap, nothing, Val(N)) return _constant_interp_nd_oneshot_soa!( diff --git a/src/core/nd_utils.jl b/src/core/nd_utils.jl index 4c2574f6..99989073 100644 --- a/src/core/nd_utils.jl +++ b/src/core/nd_utils.jl @@ -112,13 +112,10 @@ Resolve search policy input to canonical N-tuple (broadcast only, no AutoSearch _resolve_search_nd(search, Val(N), query_sample) -> NTuple{N, AbstractSearchPolicy} -Broadcast + resolve AutoSearch in one step. `query_sample` determines scalar vs vector resolution: -- `query_sample::Real` (from `first(query)`) → `Binary()` per axis -- `query_sample::AbstractVector` → `LinearBinary()` per axis - - SoA batch: pass `first(queries)` where `queries::NTuple{N,AbstractVector}` - - AoS batch: pass `queries` directly where `queries::AbstractVector{<:Tuple}`; - `AbstractVector{<:Tuple} <: AbstractVector` so resolves to `LinearBinary()` — correct - because AoS batches have the same sorted-locality property as SoA batches. +Broadcast + resolve AutoSearch in one step. Pass the query container directly — no `first()` extraction needed: +- `query_sample::NTuple{N,Real}` (scalar ND query) → `Tuple` arm → `Binary()` per axis +- `query_sample::NTuple{N,AbstractVector}` (SoA batch) → `Tuple{Vararg{AbstractVector}}` arm → `LinearBinary()` per axis +- `query_sample::AbstractVector{<:Tuple}` (AoS batch) → `AbstractVector` arm → `LinearBinary()` per axis - Explicit policies pass through unchanged. """ @inline _resolve_search_nd(s::AbstractSearchPolicy, ::Val{N}) where {N} = ntuple(_ -> s, Val(N)) diff --git a/src/cubic/nd/cubic_nd_eval.jl b/src/cubic/nd/cubic_nd_eval.jl index 3048022e..e0858e2c 100644 --- a/src/cubic/nd/cubic_nd_eval.jl +++ b/src/cubic/nd/cubic_nd_eval.jl @@ -41,7 +41,7 @@ itp((1.0, 0.5); deriv=(DerivOp(1), EvalValue())) # ∂f/∂x only ) where {Tg, Tv, N} # Note: Don't convert to Tg - preserve query type for AD support ops = _resolve_deriv_nd(deriv, Val(N)) - search_tuple = _resolve_search_nd(search, Val(N), first(query)) + search_tuple = _resolve_search_nd(search, Val(N), query) # NTuple{N,Real} <: Tuple → Binary/axis return _eval_nd_hermite(itp, query, ops, search_tuple, hint) end @@ -72,7 +72,7 @@ function (itp::CubicInterpolantND{Tg, Tv, N})( )) end ops = _resolve_deriv_nd(deriv, Val(N)) - search_tuple = _resolve_search_nd(search, Val(N), first(queries)) + search_tuple = _resolve_search_nd(search, Val(N), queries) # NTuple{N,AbstractVector} <: Tuple{Vararg{AbstractVector}} → LinearBinary/axis _batch_nd_soa!(output, itp, queries, ops, search_tuple, hint) return output end diff --git a/src/cubic/nd/cubic_nd_oneshot.jl b/src/cubic/nd/cubic_nd_oneshot.jl index 80573edd..3d59e830 100644 --- a/src/cubic/nd/cubic_nd_oneshot.jl +++ b/src/cubic/nd/cubic_nd_oneshot.jl @@ -43,7 +43,7 @@ function cubic_interp( Tr = promote_type(Tv, Tg, typeof.(query)...) bcs = _resolve_bcs_nd(bc, Val(N)) - searches = _resolve_search_nd(search, Val(N), first(query)) + searches = _resolve_search_nd(search, Val(N), query) # NTuple{N,Real} <: Tuple → Binary/axis # Validate BC requirements (once, before dispatch). _validate_nd_bcs!(grids_typed, bcs, data, Val(N)) @@ -267,7 +267,7 @@ function cubic_interp!( _validate_nd_grids(grids_typed, data) bcs = _resolve_bcs_nd(bc, Val(N)) - searches = _resolve_search_nd(search, Val(N), first(queries)) + searches = _resolve_search_nd(search, Val(N), queries) # NTuple{N,AbstractVector} <: Tuple{Vararg{AbstractVector}} → LinearBinary/axis _validate_nd_bcs!(grids_typed, bcs, data, Val(N)) diff --git a/src/integral/integrate_common_nd.jl b/src/integral/integrate_common_nd.jl index 79d6a21a..c02b8537 100644 --- a/src/integral/integrate_common_nd.jl +++ b/src/integral/integrate_common_nd.jl @@ -55,7 +55,7 @@ end _check_nd_integrate_domain(grids[d], lo2[d], extraps[d]) _check_nd_integrate_domain(grids[d], hi2[d], extraps[d]) end - search_tuple = _resolve_search_nd(search, Val(N), first(lo)) # lo[1]::Real → Binary/axis + search_tuple = _resolve_search_nd(search, Val(N), lo) # NTuple{N,Real} <: Tuple → Binary/axis idx_lo, idx_hi = _nd_cell_ranges(grids, spacings, lo2, hi2, search_tuple, hint) return (sign, lo2, hi2, idx_lo, idx_hi) end diff --git a/src/linear/nd/linear_nd_eval.jl b/src/linear/nd/linear_nd_eval.jl index a01aafca..3ead09f1 100644 --- a/src/linear/nd/linear_nd_eval.jl +++ b/src/linear/nd/linear_nd_eval.jl @@ -21,7 +21,7 @@ hint::Union{Nothing, NTuple{N, Base.RefValue{Int}}} = nothing ) where {Tg, Tv, N} ops = _resolve_deriv_nd(deriv, Val(N)) - search_tuple = _resolve_search_nd(search, Val(N), first(query)) + search_tuple = _resolve_search_nd(search, Val(N), query) # NTuple{N,Real} <: Tuple → Binary/axis return _eval_linear_nd(itp, query, ops, search_tuple, hint) end @@ -52,7 +52,7 @@ function (itp::LinearInterpolantND{Tg,Tv,N})( )) end ops = _resolve_deriv_nd(deriv, Val(N)) - search_tuple = _resolve_search_nd(search, Val(N), first(queries)) + search_tuple = _resolve_search_nd(search, Val(N), queries) # NTuple{N,AbstractVector} <: Tuple{Vararg{AbstractVector}} → LinearBinary/axis if _has_second_or_higher_derivative(ops, Val(N)) fill!(output, zero(eltype(output))) return output diff --git a/src/linear/nd/linear_nd_oneshot.jl b/src/linear/nd/linear_nd_oneshot.jl index 9e68c05d..42dbaf47 100644 --- a/src/linear/nd/linear_nd_oneshot.jl +++ b/src/linear/nd/linear_nd_oneshot.jl @@ -121,7 +121,7 @@ function linear_interp( _validate_nd_grids(grids_typed, data) Tr = promote_type(Tv, Tg, typeof.(query)...) - searches = _resolve_search_nd(search, Val(N), first(query)) + searches = _resolve_search_nd(search, Val(N), query) # NTuple{N,Real} <: Tuple → Binary/axis extraps_val = _resolve_extrap_nd(extrap, nothing, Val(N)) ops = _resolve_deriv_nd(deriv, Val(N)) @@ -196,7 +196,7 @@ function linear_interp!( grids_typed = _convert_grids_typed(grids, Tg) _validate_nd_grids(grids_typed, data) - searches = _resolve_search_nd(search, Val(N), first(queries)) + searches = _resolve_search_nd(search, Val(N), queries) # NTuple{N,AbstractVector} <: Tuple{Vararg{AbstractVector}} → LinearBinary/axis extraps_val = _resolve_extrap_nd(extrap, nothing, Val(N)) ops = _resolve_deriv_nd(deriv, Val(N)) diff --git a/src/quadratic/nd/quadratic_nd_eval.jl b/src/quadratic/nd/quadratic_nd_eval.jl index 692b1d06..f74634ac 100644 --- a/src/quadratic/nd/quadratic_nd_eval.jl +++ b/src/quadratic/nd/quadratic_nd_eval.jl @@ -66,7 +66,7 @@ itp((1.0, 0.5); deriv=(DerivOp(1), EvalValue())) # ∂f/∂x only hint::Union{Nothing, NTuple{N, Base.RefValue{Int}}}=nothing ) where {Tg, Tv, N} ops = _resolve_deriv_nd(deriv, Val(N)) - search_tuple = _resolve_search_nd(search, Val(N), first(query)) + search_tuple = _resolve_search_nd(search, Val(N), query) # NTuple{N,Real} <: Tuple → Binary/axis return _eval_nd_quadratic(itp, query, ops, search_tuple, hint) end @@ -97,7 +97,7 @@ function (itp::QuadraticInterpolantND{Tg, Tv, N})( )) end ops = _resolve_deriv_nd(deriv, Val(N)) - search_tuple = _resolve_search_nd(search, Val(N), first(queries)) + search_tuple = _resolve_search_nd(search, Val(N), queries) # NTuple{N,AbstractVector} <: Tuple{Vararg{AbstractVector}} → LinearBinary/axis _batch_nd_soa!(output, itp, queries, ops, search_tuple, hint) return output end diff --git a/src/quadratic/nd/quadratic_nd_oneshot.jl b/src/quadratic/nd/quadratic_nd_oneshot.jl index 03e62179..2afb2356 100644 --- a/src/quadratic/nd/quadratic_nd_oneshot.jl +++ b/src/quadratic/nd/quadratic_nd_oneshot.jl @@ -153,7 +153,7 @@ function quadratic_interp( Tr = promote_type(Tv, Tg, typeof.(query)...) bcs = _resolve_bcs_nd_quadratic(bc, Val(N)) - searches = _resolve_search_nd(search, Val(N), first(query)) + searches = _resolve_search_nd(search, Val(N), query) # NTuple{N,Real} <: Tuple → Binary/axis extraps_val = _resolve_extrap_nd(extrap, bcs, Val(N)) ops = _resolve_deriv_nd(deriv, Val(N)) @@ -233,7 +233,7 @@ function quadratic_interp!( _validate_nd_grids(grids_typed, data) bcs = _resolve_bcs_nd_quadratic(bc, Val(N)) - searches = _resolve_search_nd(search, Val(N), first(queries)) + searches = _resolve_search_nd(search, Val(N), queries) # NTuple{N,AbstractVector} <: Tuple{Vararg{AbstractVector}} → LinearBinary/axis extraps_val = _resolve_extrap_nd(extrap, bcs, Val(N)) ops = _resolve_deriv_nd(deriv, Val(N)) diff --git a/src/vector_calculus.jl b/src/vector_calculus.jl index 3c9f7f37..847b1e53 100644 --- a/src/vector_calculus.jl +++ b/src/vector_calculus.jl @@ -50,7 +50,7 @@ See also: [`gradient!`](@ref), [`hessian`](@ref), [`laplacian`](@ref) end for i in 1:N] return quote - search = _resolve_search_nd(itp.searches, Val($N), first(query)) # Real sample → Binary/axis (scalar ND point) + search = _resolve_search_nd(itp.searches, Val($N), query) # NTuple{N,Real} <: Tuple → Binary/axis cell = _locate_cell(itp, query, search, hint) return tuple($(deriv_calls...)) end @@ -105,7 +105,7 @@ See also: [`gradient`](@ref), [`hessian!`](@ref) @boundscheck length(G) >= $N || throw(DimensionMismatch( "gradient output vector must have at least $($N) elements, got $(length(G))" )) - search = _resolve_search_nd(itp.searches, Val($N), first(query)) # Real sample → Binary/axis (scalar ND point) + search = _resolve_search_nd(itp.searches, Val($N), query) # NTuple{N,Real} <: Tuple → Binary/axis cell = _locate_cell(itp, query, search, hint) @inbounds begin $(stmts...) @@ -180,7 +180,7 @@ See also: [`gradient`](@ref), [`hessian!`](@ref), [`laplacian`](@ref) return quote Tq = promote_type(eltype(query), $Tg, $Tv) H = Matrix{Tq}(undef, $N, $N) - search = _resolve_search_nd(itp.searches, Val($N), first(query)) # Real sample → Binary/axis (scalar ND point) + search = _resolve_search_nd(itp.searches, Val($N), query) # NTuple{N,Real} <: Tuple → Binary/axis cell = _locate_cell(itp, query, search, hint) @inbounds begin $(stmts...) @@ -251,7 +251,7 @@ See also: [`hessian`](@ref), [`gradient!`](@ref) @boundscheck size(H) == ($N, $N) || throw(DimensionMismatch( "Hessian output matrix must be $($N)×$($N), got $(size(H))" )) - search = _resolve_search_nd(itp.searches, Val($N), first(query)) # Real sample → Binary/axis (scalar ND point) + search = _resolve_search_nd(itp.searches, Val($N), query) # NTuple{N,Real} <: Tuple → Binary/axis cell = _locate_cell(itp, query, search, hint) @inbounds begin $(stmts...) @@ -316,7 +316,7 @@ See also: [`gradient`](@ref), [`hessian`](@ref) end for i in 1:N] return quote - search = _resolve_search_nd(itp.searches, Val($N), first(query)) # Real sample → Binary/axis (scalar ND point) + search = _resolve_search_nd(itp.searches, Val($N), query) # NTuple{N,Real} <: Tuple → Binary/axis cell = _locate_cell(itp, query, search, hint) return +($(deriv_calls...)) end From f1d80348bfde76621369743c8e49cc607c48da01 Mon Sep 17 00:00:00 2001 From: Min-Gu Yoo Date: Tue, 24 Feb 2026 16:51:30 -0800 Subject: [PATCH 15/16] Address Copilot review: fix docstring example and simplify test MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - linear_types.jl: replace AutoSearch() example with Binary() in "custom search policy" docstring (was showing the default, not an override) - test_search.jl: simplify redundant test expression — collapse double sin.(2π .* xq) computation + redundant .|| into single atol check Note: did NOT apply search.jl suggestion — `(hi - lo - 1) % UInt64` is valid Julia (rem-based unchecked reinterpret, not UInt64() constructor). Changing to UInt64() would reintroduce InexactError cold path in ASM. --- docs/src/api/types.md | 3 ++ docs/src/guides/search/hints.md | 21 ++++---- docs/src/guides/search/overview.md | 20 ++++---- docs/src/guides/search/policies.md | 81 ++++++++++++++++++++++-------- docs/src/nd/overview.md | 2 +- src/linear/linear_types.jl | 4 +- test/test_search.jl | 2 +- 7 files changed, 89 insertions(+), 44 deletions(-) diff --git a/docs/src/api/types.md b/docs/src/api/types.md index 8c039795..b0f6e014 100644 --- a/docs/src/api/types.md +++ b/docs/src/api/types.md @@ -102,8 +102,11 @@ WrapExtrap Search policies control how the interpolant finds the correct interval for a query point. Different policies offer trade-offs between simplicity, performance for sequential queries, and thread safety. +`AutoSearch` is the **default** for all interpolants. It resolves to `Binary()` for scalar queries and `LinearBinary()` for vector queries at call time. + ```@docs AbstractSearchPolicy +AutoSearch Binary HintedBinary Linear diff --git a/docs/src/guides/search/hints.md b/docs/src/guides/search/hints.md index a5097565..591845da 100644 --- a/docs/src/guides/search/hints.md +++ b/docs/src/guides/search/hints.md @@ -37,12 +37,14 @@ hint = Ref(1) val = itp(0.5; search=Binary(), hint=hint) # auto-upgrades to HintedBinary ``` -This means you can: -- Keep `Binary()` as your default policy -- Pass a hint only when beneficial -- Get hinted behavior automatically when a hint is provided +This also works with the default `AutoSearch()` — scalar queries resolve to `Binary()`, and if a hint is provided, they auto-upgrade to `HintedBinary`: -Without a hint, pure binary search is used. +```julia +hint = Ref(1) +val = itp(0.5; hint=hint) # AutoSearch → Binary() → auto-upgrades to HintedBinary +``` + +Without a hint, binary search is used (no hint tracking). ## Thread Safety @@ -108,16 +110,17 @@ results = itp(queries) # O(n) total instead of O(n log n) ### Random Access Pattern -For random access, hints provide no benefit. Use the default `Binary()`: +For random access, hints provide no benefit. The default `AutoSearch()` already resolves to `Binary()` for vector queries on random data — or use `Binary()` explicitly: ```julia x = collect(range(0.0, 10.0, 10001)) y = sin.(2π .* x) -itp = linear_interp(x, y) # default Binary() +itp = linear_interp(x, y) # stores AutoSearch (default) -# Random queries → Binary search is optimal +# Random queries → AutoSearch resolves to LinearBinary, but explicit Binary is faster for random random_queries = rand(100_000) .* 10 -results = itp(random_queries) +results = itp(random_queries) # AutoSearch → LinearBinary +results = itp(random_queries; search=Binary()) # explicit Binary — better for random ``` ### Shared Hint Across Multiple Interpolants diff --git a/docs/src/guides/search/overview.md b/docs/src/guides/search/overview.md index b6add411..75b8a5b5 100644 --- a/docs/src/guides/search/overview.md +++ b/docs/src/guides/search/overview.md @@ -17,15 +17,15 @@ using FastInterpolations x = collect(range(0, 1.0, length=1000)) y = x.^3 -# For random queries +# AutoSearch is the default — resolves automatically per call xq = rand(1000) -linear_interp(x, y, xq; search=Binary()) # Default: optimal for random access -linear_interp(x, y, xq; search=LinearBinary()) # Could be slower +linear_interp(x, y, xq) # AutoSearch → LinearBinary() for vector +linear_interp(x, y, xq; search=Binary()) # Explicit: optimal for random access # For sorted/monotonic queries xq_sorted = sort(xq) -linear_interp(x, y, xq_sorted; search=Binary()) # Same as random -linear_interp(x, y, xq_sorted; search=LinearBinary()) # Much faster! +linear_interp(x, y, xq_sorted) # AutoSearch → LinearBinary() — already fast! +linear_interp(x, y, xq_sorted; search=LinearBinary()) # Same result, explicit nothing # hide ``` @@ -33,9 +33,10 @@ nothing # hide | Policy | Best For | Complexity | Thread Safety | |:-------|:---------|:-----------|:--------------| -| [`Binary()`](@ref search_policies) | Random access (default) | O(log n) | ✓ Stateless | +| [`AutoSearch()`](@ref search_policies) | **General use (default)** — adapts per query type | Delegates to Binary/LinearBinary | ✓ Stateless | +| [`Binary()`](@ref search_policies) | Random access (explicit) | O(log n) | ✓ Stateless | | [`HintedBinary()`](@ref search_policies) | Repeated queries in same region | O(1) hit, O(log n) miss | ✓ With hint | -| [`LinearBinary()`](@ref search_policies) | **Monotonic queries (recommended)** | O(1) local, O(log n) fallback | ✓ With hint | +| [`LinearBinary()`](@ref search_policies) | **Monotonic queries (explicit)** | O(1) local, O(log n) fallback | ✓ With hint | | [`Linear()`](@ref search_policies) | Close + monotonic queries (expert) | O(1) amortized | ✓ With hint | !!! note "Why No Hunt Algorithm?" @@ -47,9 +48,10 @@ nothing # hide **Which policy should I use?** -- **Random access / general use** → `Binary()` (default, most consistent performance) +- **General use / unknown pattern** → `AutoSearch()` ✅ **default** — adapts scalar→`Binary()`, vector→`LinearBinary()` +- **Known random access** → `Binary()` (explicit; skips AutoSearch dispatch) - **Queries cluster in same region** → `HintedBinary()` -- **Monotonic queries (sorted, streaming, ODE)** → `LinearBinary()` ✅ **recommended** +- **Known monotonic queries (sorted, streaming, ODE)** → `LinearBinary()` (explicit) - **Strictly monotonic, performance-critical** → `Linear()` (benchmark first! see below) !!! note "Benchmark Before Choosing Linear()" diff --git a/docs/src/guides/search/policies.md b/docs/src/guides/search/policies.md index 4aaafae6..4affec65 100644 --- a/docs/src/guides/search/policies.md +++ b/docs/src/guides/search/policies.md @@ -2,7 +2,41 @@ This page explains each search policy in detail, including when to use it and how it works internally. -## Binary (Default) +## AutoSearch (Default) + +Automatically selects `Binary()` or `LinearBinary()` based on query type at call time. This is the **default** for all interpolants — no configuration needed. + +**Resolution rules**: + +| Query type | Resolved policy | Rationale | +|:-----------|:----------------|:----------| +| Scalar (`Real`) | `Binary()` | Stateless; no hint locality to exploit | +| Vector (`AbstractVector`) | `LinearBinary()` | Exploits hint locality for sorted sequences | +| ND scalar (`Tuple{Vararg{Real}}`) | `Binary()` per axis | Same as 1D scalar | +| ND SoA batch (`NTuple{N, AbstractVector}`) | `LinearBinary()` per axis | Same as 1D vector | +| Broadcast (`itp.(xs)`) | `Binary()` per element | Each broadcast call is independent | + +```julia +itp = linear_interp(x, y) # stores AutoSearch (the default) +itp(0.5) # → Binary() internally (scalar query) +itp([0.1, 0.5, 0.9]) # → LinearBinary() internally (vector query) +itp.(xs) # → Binary() per element (broadcast) + +# Override when you know the pattern in advance: +itp(0.5; search=Binary()) # force Binary for all calls +itp(sorted_xs; search=LinearBinary()) # force LinearBinary for all calls +``` + +**How it works**: `AutoSearch` is resolved at the call site, not at construction time. The interpolant stores `AutoSearch()` and dispatches on the query argument's concrete type each time the interpolant is called. + +**When to override with an explicit policy**: +- You know queries are **always random** → explicit `Binary()` skips the dispatch check +- You know queries are **always sorted** → explicit `LinearBinary()` skips the dispatch check +- For most use cases, keeping `AutoSearch()` is the right choice + +--- + +## Binary Pure binary search. Stateless, thread-safe, no setup required. @@ -11,12 +45,12 @@ Pure binary search. Stateless, thread-safe, no setup required. **When to use**: - Random access patterns - One-off queries -- When you don't know the query pattern in advance +- When you want consistent O(log n) for any query order ```julia itp = linear_interp(x, y) -val = itp(0.5) # uses Binary() by default -val = itp(0.5; search=Binary()) # explicit +val = itp(0.5) # AutoSearch → Binary() for scalar queries +val = itp(0.5; search=Binary()) # explicit Binary ``` **How it works**: Standard binary search on the grid. Each query starts fresh with no memory of previous queries. @@ -108,15 +142,15 @@ vals = itp(sorted_queries) # O(1) amortized for sorted input You can tune the linear search window size before falling back to binary search: ```julia -LinearBinary() # default: linear_window=8 -LinearBinary(linear_window=4) # smaller bound for tightly spaced queries -LinearBinary(linear_window=16) # larger bound for sparser queries +LinearBinary() # default: linear_window=2 +LinearBinary(linear_window=4) # moderate bound for known-sorted sequences +LinearBinary(linear_window=16) # larger bound for sparser-spaced sorted queries ``` **Guidelines**: -- **Small `linear_window` (4)**: Better when queries are very close together -- **Large `linear_window` (16-32)**: Better when queries might skip a few intervals -- **Default (8)**: Good balance for most use cases +- **Small `linear_window` (1–2)**: Minimal overhead; best for mixed or unknown patterns. The default `LinearBinary()` uses `2`. +- **Medium `linear_window` (4–16)**: Good balance when queries are known-sorted +- **Large `linear_window` (32–128)**: For highly localized queries or very large datasets !!! note "Type Parameter Restriction" `linear_window` is restricted to powers of 2 (1, 2, 4, 8, 16, 32, 64) to prevent type parameter explosion. Each unique value creates a specialized method, so limiting choices keeps compile times reasonable. @@ -125,10 +159,12 @@ LinearBinary(linear_window=16) # larger bound for sparser queries ## Performance Summary -| Query Pattern | `Binary()` | `LinearBinary()` | `Linear()` | -|:--------------|:-----------|:-----------------|:-----------| -| **Random** | ✅ Best | ~2-3x slower | ❌ Up to 7x slower | -| **Monotonic** | Baseline | ✅ ~5x faster | ✅ ~6x faster | +| Query Pattern | `AutoSearch()` | `Binary()` | `LinearBinary()` | `Linear()` | +|:--------------|:---------------|:-----------|:-----------------|:-----------| +| **Random** | ✅ Same as `Binary()` | ✅ Best | ~2-3x slower | ❌ Up to 7x slower | +| **Monotonic** | ✅ Same as `LinearBinary()` | Baseline | ✅ ~5x faster | ✅ ~6x faster | + +`AutoSearch()` has negligible overhead — it dispatches to the same underlying code as the resolved policy, so you only pay for the type dispatch, not for any extra search work. !!! note "Results Vary" These are approximate results from a 500-point grid with 1000 queries. Actual performance depends on your **grid size** and **query spacing**. Run benchmark with your own data to find the best policy. @@ -139,23 +175,24 @@ LinearBinary(linear_window=16) # larger bound for sparser queries ### Baked-in Policy (Constructor) -Set the default policy when creating the interpolant: +Override the default `AutoSearch` when creating the interpolant. The stored policy applies to all subsequent calls: ```julia -# All queries use LinearBinary by default -itp = linear_interp(x, y; search=LinearBinary()) -val = itp(0.5) # uses LinearBinary +itp = linear_interp(x, y) # stores AutoSearch() — adapts per call +itp = linear_interp(x, y; search=Binary()) # always Binary(), regardless of query type +itp = linear_interp(x, y; search=LinearBinary()) # always LinearBinary() ``` ### Override at Call Time -Override the baked-in policy for specific queries: +Override the stored policy for a single call without changing the interpolant: ```julia -itp = linear_interp(x, y; search=LinearBinary()) +itp = linear_interp(x, y) # stores AutoSearch (default) -# Override for a single call -val = itp(0.5; search=Binary()) # uses Binary for this call only +# Call-site override — does not change itp.search_policy +val = itp(0.5; search=Binary()) # force Binary for this call only +vals = itp(sorted_xs; search=LinearBinary()) # force LinearBinary for this call only ``` This is useful when most queries benefit from one policy, but occasional queries need different behavior. diff --git a/docs/src/nd/overview.md b/docs/src/nd/overview.md index d5157084..f42d7fb3 100644 --- a/docs/src/nd/overview.md +++ b/docs/src/nd/overview.md @@ -39,7 +39,7 @@ Every 1D argument becomes a Tuple in ND. This applies uniformly: | BC | `bc=CubicFit()` | `bc=(CubicFit(), PeriodicBC())` | | Extrap | `extrap=ConstExtrap()` | `extrap=(ConstExtrap(), WrapExtrap())` | | Derivative | `deriv=DerivOp(1)` | `deriv=DerivOp(1, 0)` for ∂f/∂x | -| Search | `search=Binary()` | `search=(Binary(), LinearBinary())` | +| Search | `search=AutoSearch()` (default) | `search=AutoSearch()` (default) or `search=(Binary(), LinearBinary())` | **Broadcast rule**: A scalar value is broadcast to all axes. `bc=CubicFit()` is equivalent to `bc=(CubicFit(), CubicFit())` in 2D. diff --git a/src/linear/linear_types.jl b/src/linear/linear_types.jl index 84064973..4b40f959 100644 --- a/src/linear/linear_types.jl +++ b/src/linear/linear_types.jl @@ -29,8 +29,8 @@ Returned by `linear_interp(x, y)` (2-argument form). # Create interpolator (minimal allocation) itp = linear_interp(x, y) # default extrap=NoExtrap(), search=AutoSearch() -# Create with custom search policy (baked-in default) -itp = linear_interp(x, y; search=AutoSearch()) +# Create with custom search policy (overrides default AutoSearch) +itp = linear_interp(x, y; search=Binary()) # Complex-valued interpolation y_complex = exp.(2im .* x) diff --git a/test/test_search.jl b/test/test_search.jl index 73afdd10..fb3b4589 100644 --- a/test/test_search.jl +++ b/test/test_search.jl @@ -1188,7 +1188,7 @@ using FastInterpolations: search_interval, _search_binary, _search_direct, _sear xq = [0.1, 0.3, 0.7] vals = itp(xq) @test length(vals) == 3 - @test all(vals .≈ sin.(2π .* xq) .|| abs.(vals .- sin.(2π .* xq)) .< 0.02) + @test vals ≈ sin.(2π .* xq) atol=0.02 # In-place vector call out = zeros(3) From 4637ff66c2dfa90d590284ddf5e1288eafbe0a5d Mon Sep 17 00:00:00 2001 From: Min-Gu Yoo Date: Tue, 24 Feb 2026 17:01:21 -0800 Subject: [PATCH 16/16] Update performance numbers from fresh benchmark (500-2000pt grids) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Random LB penalty: ~2-3x → ~2.5-3x (measured: 2.4-2.8x for vector batch) - Monotonic LB gain: ~5x → ~4-6x (measured: 3.7-4.8x at 500pt, 4.5-6.1x at 2000pt) - Expand note: clarify "vector batch calls" context, mention scalar-no-hint is only ~1.2x slower (hint walk overhead absent without persistence) --- docs/release_notes/pr_autosearch.md | 48 +++++++++++++++++++++++++++++ docs/src/guides/search/policies.md | 10 ++++-- 2 files changed, 55 insertions(+), 3 deletions(-) create mode 100644 docs/release_notes/pr_autosearch.md diff --git a/docs/release_notes/pr_autosearch.md b/docs/release_notes/pr_autosearch.md new file mode 100644 index 00000000..ba0faca2 --- /dev/null +++ b/docs/release_notes/pr_autosearch.md @@ -0,0 +1,48 @@ +## Summary + +Introduce `AutoSearch` — an adaptive default search policy that resolves to `Binary()` for scalar queries and `LinearBinary()` for vector queries at call time. Also optimizes `LinearBinary` internals (branchless binary, window default 2→ optimal sweet spot) and ships comprehensive tests and documentation. + +## Motivation + +Previously, `Binary()` was the unconditional default. This meant vector queries over sorted data — the most common batch use case — silently used O(log n) binary search instead of the much faster O(1)-amortized `LinearBinary`. Users had to opt in manually, and many didn't know to. + +The new `AutoSearch` default removes this decision entirely for 95% of users: + +| Query type | Old behavior | New behavior | Effect | +|:-----------|:-------------|:-------------|:-------| +| Scalar (e.g., `itp(0.5)`) | `Binary()` (was default) | `Binary()` via AutoSearch | Unchanged | +| Vector (e.g., `itp(Vector)`) | `Binary()` (was default) | `LinearBinary()` via AutoSearch | Up to ~5× faster on sorted data | + +## Key Changes + +### `AutoSearch` (new type) + +```julia +itp = linear_interp(x, y) # stores AutoSearch — the new default +itp(0.5) # → Binary() (scalar) +itp([0.1, 0.5, 0.9]) # → LinearBinary() (vector) +itp(0.5; search=Binary()) # explicit override still works +``` + +Resolution happens once per call, at the eval entry point, with negligible overhead. + +### `LinearBinary` optimization + +- **Branchless binary core**: replaced `while` loop with `for _ in 1:iters` where `iters = 64 - leading_zeros(...)`. Precomputed trip count → constant-iteration loop + `ifelse` → ARM64 `csel` (branch-free). ~25–55% faster than the old loop on random queries. +- **Default window changed**: `LinearBinary()` now defaults to `linear_window=2` (was 8). Window=2 minimizes overhead for mixed/unknown patterns while still exploiting locality for sorted sequences. + +### Safety fix + +`_search_linear_binary!` now clamps the hint before first use (`ix = clamp(ix, 1, n-1)`), guarding against user-provided out-of-range hints (`Ref(0)`, stale hints from a different grid). + +## New Export + +```julia +export AutoSearch +``` + +## Impact + +- **Zero breaking changes** for existing code: explicit `search=Binary()` or `search=LinearBinary()` at constructor or call site is honored as-is. +- **Behavior change** for users who relied on the default being `Binary()` for vector queries. Vector queries now use `LinearBinary()` via AutoSearch. For **random** vector queries, this is ~2.5–3× slower than `Binary()` — use `search=Binary()` explicitly to restore the old behavior. +- All 1D, ND, series, and integration paths updated. 440+ new test assertions. diff --git a/docs/src/guides/search/policies.md b/docs/src/guides/search/policies.md index 4affec65..ed2c96e0 100644 --- a/docs/src/guides/search/policies.md +++ b/docs/src/guides/search/policies.md @@ -161,13 +161,17 @@ LinearBinary(linear_window=16) # larger bound for sparser-spaced sorted queries | Query Pattern | `AutoSearch()` | `Binary()` | `LinearBinary()` | `Linear()` | |:--------------|:---------------|:-----------|:-----------------|:-----------| -| **Random** | ✅ Same as `Binary()` | ✅ Best | ~2-3x slower | ❌ Up to 7x slower | -| **Monotonic** | ✅ Same as `LinearBinary()` | Baseline | ✅ ~5x faster | ✅ ~6x faster | +| **Random** | ✅ Same as `Binary()` | ✅ Best | ~2.5-3x slower | ❌ Up to 7x slower | +| **Monotonic** | ✅ Same as `LinearBinary()` | Baseline | ✅ ~4-6x faster | ✅ ~6x faster | `AutoSearch()` has negligible overhead — it dispatches to the same underlying code as the resolved policy, so you only pay for the type dispatch, not for any extra search work. !!! note "Results Vary" - These are approximate results from a 500-point grid with 1000 queries. Actual performance depends on your **grid size** and **query spacing**. Run benchmark with your own data to find the best policy. + These are approximate results from **vector batch calls** (`itp(out, queries)`) on 500–2000 point grids. + The random penalty (~2.5-3x) comes from the hint walk landing at the wrong position before falling back to + binary. The monotonic speedup grows with grid size (larger grids benefit more from hint locality). + For **scalar calls** without a persistent hint, the difference is much smaller (~1.2x). + Run benchmark with your own data to find the best policy. ---