diff --git a/README.md b/README.md index ab190214..96f40c5e 100644 --- a/README.md +++ b/README.md @@ -139,9 +139,9 @@ cubic_interp(x, y, 5.0; deriv=2) # 2nd derivative at x=5.0 cubic_interp(x, y, 5.0; deriv=3) # 3rd derivative at x=5.0 # Constant interpolation — choose which side to sample -constant_interp(x, y, xq; side=:nearest) # nearest neighbor (default) -constant_interp(x, y, xq; side=:left) # left-continuous -constant_interp(x, y, xq; side=:right) # right-continuous +constant_interp(x, y, xq; side=NearestSide()) # nearest neighbor (default) +constant_interp(x, y, xq; side=LeftSide()) # left-continuous +constant_interp(x, y, xq; side=RightSide()) # right-continuous # Quadratic boundary condition — single endpoint constraint quadratic_interp(x, y, xq; bc=Left(Deriv1(0.0))) # S'(left) = 0 @@ -154,9 +154,9 @@ cubic_interp(x, y, xq; bc=BCPair(Deriv1(2.0), Deriv2(-5.0))) # custom (left, ri cubic_interp(x, y, xq; bc=CubicFit()) # Estimate derivatives using 4-point fit at both ends # Extrapolation modes — all methods support these -linear_interp(x, y, xq; extrap=:constant) # clamp to boundary values -quadratic_interp(x, y, xq; extrap=:wrap) # wrap around (periodic data) -cubic_interp(x, y, xq; extrap=:extension) # extend boundary polynomial +linear_interp(x, y, xq; extrap=ConstExtrap()) # clamp to boundary values +quadratic_interp(x, y, xq; extrap=WrapExtrap()) # wrap around (periodic data) +cubic_interp(x, y, xq; extrap=ExtendExtrap()) # extend boundary polynomial ``` ## Documentation diff --git a/docs/src/api/constant.md b/docs/src/api/constant.md index 23e4161b..ad9cfd46 100644 --- a/docs/src/api/constant.md +++ b/docs/src/api/constant.md @@ -7,7 +7,7 @@ | Function | Description | |----------|-------------| | `constant_interp(x, y, xq)` | Constant interpolation at point(s) `xq` | -| `constant_interp(x, y, xq; side=:left)` | With side mode (`:nearest`, `:left`, `:right`) | +| `constant_interp(x, y, xq; side=LeftSide())` | With side mode (`NearestSide()`, `LeftSide()`, `RightSide()`) | | `constant_interp!(out, x, y, xq)` | In-place constant interpolation | | `constant_interp!(out, x, y, xq; side)` | In-place with side mode | @@ -16,7 +16,7 @@ | Function | Description | |----------|-------------| | `itp = constant_interp(x, y)` | Create constant interpolant | -| `itp = constant_interp(x, y; side=:left)` | Create with side mode | +| `itp = constant_interp(x, y; side=LeftSide())` | Create with side mode | | `itp(xq)` | Evaluate at point(s) `xq` | | `itp(out, xq)` | Evaluate at `xq`, store result in `out` | @@ -46,6 +46,15 @@ ConstantInterpolant ConstantInterpolantND ``` +## Side Selection Types + +```@docs +AbstractSide +NearestSide +LeftSide +RightSide +``` + ## Derivative Views See [Derivatives](../interpolation/derivatives.md) for `deriv1`, `deriv2`, `deriv3` API reference. diff --git a/docs/src/interpolation/constant.md b/docs/src/interpolation/constant.md index 860d9511..bbc3ae66 100644 --- a/docs/src/interpolation/constant.md +++ b/docs/src/interpolation/constant.md @@ -2,13 +2,13 @@ Step function interpolation that returns values from neighboring data points without blending. -**Key Feature**: The `side` keyword controls which neighbor to sample. +**Key Feature**: The `side` keyword controls which neighbor to sample using [`AbstractSide`](@ref) types. | Mode | Behavior | |------|----------| -| `:nearest` | Nearest neighbor — **default** | -| `:left` | Left-continuous (use left neighbor) | -| `:right` | Right-continuous (use right neighbor) | +| `NearestSide()` | Nearest neighbor — **default** | +| `LeftSide()` | Left-continuous (use left neighbor) | +| `RightSide()` | Right-continuous (use right neighbor) | --- @@ -22,11 +22,11 @@ y = [1.0, 3.0, 2.0, 4.0, 1.5] # One-shot evaluation constant_interp(x, y, 1.5) # → 3.0 (nearest neighbor) -constant_interp(x, y, 1.5; side=:left) # → 3.0 (left neighbor) -constant_interp(x, y, 1.5; side=:right) # → 2.0 (right neighbor) +constant_interp(x, y, 1.5; side=LeftSide()) # → 3.0 (left neighbor) +constant_interp(x, y, 1.5; side=RightSide()) # → 2.0 (right neighbor) # Create reusable interpolant -itp = constant_interp(x, y; side=:nearest) +itp = constant_interp(x, y; side=NearestSide()) itp(1.5) # single point itp([1.0, 2.0]) # multiple points @@ -55,12 +55,14 @@ x = [0.0, 1.0, 2.5, 4.0, 5.0] y = [1.0, 3.0, 2.0, 4.0, 1.5] xq = range(x[1], x[end], 100) +sides = [NearestSide(), LeftSide(), RightSide()] +labels = ["NearestSide()", "LeftSide()", "RightSide()"] p = plot(layout=(1,3), size=(900, 280), legend=:topright) -for (i, side) in enumerate([:nearest, :left, :right]) +for (i, (side, lbl)) in enumerate(zip(sides, labels)) yq = constant_interp(x, y, xq; side=side) plot!(p[i], xq, yq, label="spline", linewidth=2) scatter!(p[i], x, y, label="data", markersize=7, color=:black) - title!(p[i], "side=:$side") + title!(p[i], "side=$lbl") end p ``` diff --git a/src/FastInterpolations.jl b/src/FastInterpolations.jl index 442e10e1..7e456f53 100644 --- a/src/FastInterpolations.jl +++ b/src/FastInterpolations.jl @@ -59,6 +59,9 @@ export AbstractEvalOp, EvalValue, EvalDeriv1, EvalDeriv2, EvalDeriv3 # Extrapolation mode types (typed API for ND extrap) export AbstractExtrap, NoExtrap, ConstExtrap, ExtendExtrap, WrapExtrap +# Side selection types (constant interpolation) +export AbstractSide, NearestSide, LeftSide, RightSide + # Integration export integrate, cumulative_integrate diff --git a/src/constant/constant_interpolant.jl b/src/constant/constant_interpolant.jl index 564dd296..230c246c 100644 --- a/src/constant/constant_interpolant.jl +++ b/src/constant/constant_interpolant.jl @@ -57,7 +57,7 @@ end # ======================================== """ - constant_interp(x, y; extrap=NoExtrap(), side=:nearest, search=Binary()) -> ConstantInterpolant + constant_interp(x, y; extrap=NoExtrap(), side=NearestSide(), search=Binary()) -> ConstantInterpolant Create a callable interpolant for broadcast fusion and reuse. @@ -65,7 +65,7 @@ Create a callable interpolant for broadcast fusion and reuse. - `x::AbstractVector`: x-coordinates (sorted, length ≥ 2) - `y::AbstractVector`: y-values (can be Real or Complex) - `extrap::AbstractExtrap`: `NoExtrap()` (default), `ConstExtrap()`, `ExtendExtrap()`, or `WrapExtrap()` -- `side::Symbol`: Side selection +- `side::AbstractSide`: Side selection (NearestSide(), LeftSide(), RightSide()) - `search::AbstractSearchPolicy`: Default search policy (default: `Binary()`) # Returns @@ -114,7 +114,7 @@ end x::AbstractVector{TX}, y::AbstractVector{TY}; extrap::AbstractExtrap=NoExtrap(), - side::Symbol=:nearest, + side::AbstractSide=NearestSide(), search::AbstractSearchPolicy=Binary() ) where {TX<:Real, TY} x_p, y_p = _promote_itp_inputs(x, y) diff --git a/src/constant/constant_kernels.jl b/src/constant/constant_kernels.jl index 7274ffdb..60f424d4 100644 --- a/src/constant/constant_kernels.jl +++ b/src/constant/constant_kernels.jl @@ -8,7 +8,7 @@ # - h = x_{i+1} - x_i (interval width, Tg) # - dL = xq - x_i (offset from left boundary, can be Dual for AD) # - y_left, y_right = values (Tv, can be Complex) -# - side = Val(:nearest) | Val(:left) | Val(:right) +# - side = NearestSide() | LeftSide() | RightSide() # # Type parameters: # - Tg<:AbstractFloat: Grid type (geometry) @@ -23,37 +23,37 @@ # all side modes return y_left (the value at that grid point). """ - _constant_kernel(::EvalValue, y_left, y_right, h, dL, ::Val{:left}) + _constant_kernel(::EvalValue, y_left, y_right, h, dL, ::LeftSide) Constant interpolation with left-continuous (floor) convention. Always returns the left boundary value `y_left`. dL can be any Real (including ForwardDiff.Dual for AD). """ -@inline function _constant_kernel(::EvalValue, y_left::Tv, ::Tv, ::Tg, dL::Td, ::Val{:left}) where {Tv, Tg<:AbstractFloat, Td<:Real} +@inline function _constant_kernel(::EvalValue, y_left::Tv, ::Tv, ::Tg, dL::Td, ::LeftSide) where {Tv, Tg<:AbstractFloat, Td<:Real} return y_left end """ - _constant_kernel(::EvalValue, y_left, y_right, h, dL, ::Val{:right}) + _constant_kernel(::EvalValue, y_left, y_right, h, dL, ::RightSide) Constant interpolation with right-continuous (ceiling) convention. Returns `y_left` at grid point (dL == 0), `y_right` otherwise. dL can be any Real (including ForwardDiff.Dual for AD). """ -@inline function _constant_kernel(::EvalValue, y_left::Tv, y_right::Tv, ::Tg, dL::Td, ::Val{:right}) where {Tv, Tg<:AbstractFloat, Td<:Real} +@inline function _constant_kernel(::EvalValue, y_left::Tv, y_right::Tv, ::Tg, dL::Td, ::RightSide) where {Tv, Tg<:AbstractFloat, Td<:Real} # Use primal value for comparison (supports ForwardDiff.Dual) dL_primal = _extract_primal(dL) return iszero(dL_primal) ? y_left : y_right end """ - _constant_kernel(::EvalValue, y_left, y_right, h, dL, ::Val{:nearest}) + _constant_kernel(::EvalValue, y_left, y_right, h, dL, ::NearestSide) Constant interpolation with nearest-neighbor convention and left tie-breaking. Returns `y_left` if dL <= h/2 (including midpoint), `y_right` otherwise. dL can be any Real (including ForwardDiff.Dual for AD). """ -@inline function _constant_kernel(::EvalValue, y_left::Tv, y_right::Tv, h::Tg, dL::Td, ::Val{:nearest}) where {Tv, Tg<:AbstractFloat, Td<:Real} +@inline function _constant_kernel(::EvalValue, y_left::Tv, y_right::Tv, h::Tg, dL::Td, ::NearestSide) where {Tv, Tg<:AbstractFloat, Td<:Real} # Use primal value for comparison (supports ForwardDiff.Dual) dL_primal = _extract_primal(dL) return dL_primal <= h / 2 ? y_left : y_right @@ -67,7 +67,7 @@ Always returns zero (constant function has no slope). Returns `zero(Tv)` to preserve Complex type when applicable. dL can be any Real (including ForwardDiff.Dual for AD). """ -@inline function _constant_kernel(::EvalDeriv1, y_left::Tv, ::Tv, ::Tg, dL::Td, ::SideVal) where {Tv, Tg<:AbstractFloat, Td<:Real} +@inline function _constant_kernel(::EvalDeriv1, y_left::Tv, ::Tv, ::Tg, dL::Td, ::AbstractSide) where {Tv, Tg<:AbstractFloat, Td<:Real} return zero(Tv) end @@ -79,7 +79,7 @@ Always returns zero (constant function has no curvature). Returns `zero(Tv)` to preserve Complex type when applicable. dL can be any Real (including ForwardDiff.Dual for AD). """ -@inline function _constant_kernel(::EvalDeriv2, y_left::Tv, ::Tv, ::Tg, dL::Td, ::SideVal) where {Tv, Tg<:AbstractFloat, Td<:Real} +@inline function _constant_kernel(::EvalDeriv2, y_left::Tv, ::Tv, ::Tg, dL::Td, ::AbstractSide) where {Tv, Tg<:AbstractFloat, Td<:Real} return zero(Tv) end @@ -91,6 +91,6 @@ Constant functions have all derivatives equal to zero. Returns `zero(Tv)` to preserve Complex type when applicable. dL can be any Real (including ForwardDiff.Dual for AD). """ -@inline function _constant_kernel(::EvalDeriv3, y_left::Tv, ::Tv, ::Tg, dL::Td, ::SideVal) where {Tv, Tg<:AbstractFloat, Td<:Real} +@inline function _constant_kernel(::EvalDeriv3, y_left::Tv, ::Tv, ::Tg, dL::Td, ::AbstractSide) where {Tv, Tg<:AbstractFloat, Td<:Real} return zero(Tv) end diff --git a/src/constant/constant_oneshot.jl b/src/constant/constant_oneshot.jl index 43be6b2b..0f402bd1 100644 --- a/src/constant/constant_oneshot.jl +++ b/src/constant/constant_oneshot.jl @@ -31,7 +31,7 @@ Type parameters: """ @inline function _constant_eval_extrap( y::AbstractVector{Tv}, xi::Tg, x_min::Tg, x_max::Tg, - ::ConstExtrap, ::SideVal, ::EvalValue + ::ConstExtrap, ::AbstractSide, ::EvalValue ) where {Tg<:AbstractFloat, Tv} if xi < x_min return @inbounds y[1] @@ -42,21 +42,21 @@ end @inline function _constant_eval_extrap( y::AbstractVector{Tv}, ::Tg, ::Tg, ::Tg, - ::ConstExtrap, ::SideVal, ::EvalDeriv1 + ::ConstExtrap, ::AbstractSide, ::EvalDeriv1 ) where {Tg<:AbstractFloat, Tv} return zero(Tv) end @inline function _constant_eval_extrap( y::AbstractVector{Tv}, ::Tg, ::Tg, ::Tg, - ::ConstExtrap, ::SideVal, ::EvalDeriv2 + ::ConstExtrap, ::AbstractSide, ::EvalDeriv2 ) where {Tg<:AbstractFloat, Tv} return zero(Tv) end @inline function _constant_eval_extrap( y::AbstractVector{Tv}, ::Tg, ::Tg, ::Tg, - ::ConstExtrap, ::SideVal, ::EvalDeriv3 + ::ConstExtrap, ::AbstractSide, ::EvalDeriv3 ) where {Tg<:AbstractFloat, Tv} return zero(Tv) end @@ -64,7 +64,7 @@ end # ExtendExtrap delegates to ConstExtrap (slope=0 for constant function) @inline function _constant_eval_extrap( y::AbstractVector{Tv}, xi::Tg, x_min::Tg, x_max::Tg, - ::ExtendExtrap, side::SideVal, op::AbstractEvalOp + ::ExtendExtrap, side::AbstractSide, op::AbstractEvalOp ) where {Tg<:AbstractFloat, Tv} return _constant_eval_extrap(y, xi, x_min, x_max, ConstExtrap(), side, op) end @@ -96,7 +96,7 @@ AD Support: y::AbstractVector{Tv}, xi::Tq, extrap::AbstractExtrap, - side::SideVal, + side::AbstractSide, op::AbstractEvalOp, searcher::S ) where {Tg<:AbstractFloat, Tv, Tq<:Real, S<:Searcher} @@ -149,7 +149,7 @@ end # ======================================== """ - constant_interp(x, y, xi; extrap=NoExtrap(), side=:nearest, deriv=0, search=Binary()) + constant_interp(x, y, xi; extrap=NoExtrap(), side=NearestSide(), deriv=0, search=Binary()) Constant (step/piecewise constant) interpolation at a single point. @@ -162,10 +162,10 @@ Constant (step/piecewise constant) interpolation at a single point. - `ConstExtrap()`: clamp to boundary values - `ExtendExtrap()`: same as ConstExtrap (slope=0) - `WrapExtrap()`: wrap to [x_min, x_max) -- `side::Symbol`: Side selection - - `:nearest` (default): nearest neighbor (left tie-breaking at midpoint) - - `:left`: always use left value - - `:right`: use right value (except at grid points) +- `side::AbstractSide`: Side selection + - `NearestSide()` (default): nearest neighbor (left tie-breaking at midpoint) + - `LeftSide()`: always use left value + - `RightSide()`: use right value (except at grid points) - `deriv::Int`: Derivative order (0, 1, or 2). Derivatives are always 0. - `search::AbstractSearchPolicy`: Search algorithm for interval finding - `Binary()` (default): O(log n) binary search, stateless @@ -181,8 +181,8 @@ x = [0.0, 1.0, 2.0, 3.0] y = [10.0, 20.0, 30.0, 40.0] constant_interp(x, y, 0.5) # 10.0 (nearest to left) -constant_interp(x, y, 0.5; side=:left) # 10.0 -constant_interp(x, y, 0.5; side=:right) # 20.0 +constant_interp(x, y, 0.5; side=LeftSide()) # 10.0 +constant_interp(x, y, 0.5; side=RightSide()) # 20.0 constant_interp(x, y, 1.0) # 20.0 (grid point) constant_interp(x, y, -1.0; extrap=ConstExtrap()) # 10.0 (clamped) @@ -198,7 +198,7 @@ vals = constant_interp(x, y, sorted_queries; search=LinearBinary(linear_window=8 y::AbstractVector{Tv}, xi::Tq; extrap::AbstractExtrap=NoExtrap(), - side::Symbol=:nearest, + side::AbstractSide=NearestSide(), deriv::Int=0, search=Binary(), hint::Union{Nothing,Base.RefValue{Int}}=nothing @@ -207,9 +207,7 @@ vals = constant_interp(x, y, sorted_queries; search=LinearBinary(linear_window=8 searcher = _to_searcher(search, hint) @_dispatch_deriv deriv => op begin - @_dispatch_side side => sv begin - _constant_eval_at_point(x, y, xi, extrap, sv, op, searcher) - end + _constant_eval_at_point(x, y, xi, extrap, side, op, searcher) end end @@ -218,7 +216,7 @@ end # ======================================== """ - constant_interp!(output, x, y, x_targets; extrap=NoExtrap(), side=:nearest, deriv=0, search=Binary()) + constant_interp!(output, x, y, x_targets; extrap=NoExtrap(), side=NearestSide(), deriv=0, search=Binary()) Zero-allocation constant interpolation for multiple query points. @@ -248,7 +246,7 @@ function constant_interp!( y::AbstractVector{Tv}, x_targets::AbstractVector{Tg}; extrap::AbstractExtrap=NoExtrap(), - side::Symbol=:nearest, + side::AbstractSide=NearestSide(), deriv::Int=0, search::AbstractSearchPolicy=Binary() ) where {Tg<:AbstractFloat, Tv} @@ -257,11 +255,9 @@ function constant_interp!( searcher = _to_searcher(search) @_dispatch_deriv deriv => op begin - @_dispatch_side side => sv begin - @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, sv, op, searcher) - end + @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, op, searcher) end end return output @@ -272,7 +268,7 @@ end # ======================================== """ - constant_interp(x, y, x_targets; extrap=NoExtrap(), side=:nearest, deriv=0, search=Binary()) + constant_interp(x, y, x_targets; extrap=NoExtrap(), side=NearestSide(), deriv=0, search=Binary()) Constant interpolation for multiple query points (allocating version). @@ -293,7 +289,7 @@ function constant_interp( y::AbstractVector{Tv}, x_targets::AbstractVector{Tg}; extrap::AbstractExtrap=NoExtrap(), - side::Symbol=:nearest, + side::AbstractSide=NearestSide(), deriv::Int=0, search::AbstractSearchPolicy=Binary() ) where {Tg<:AbstractFloat, Tv} @@ -320,7 +316,7 @@ end y::AbstractVector{Tv}, xi::Tq; extrap::AbstractExtrap=NoExtrap(), - side::Symbol=:nearest, + side::AbstractSide=NearestSide(), deriv::Int=0, search=Binary(), hint::Union{Nothing,Base.RefValue{Int}}=nothing @@ -340,7 +336,7 @@ function constant_interp( y::AbstractVector{Tv}, x_targets::AbstractVector{Tq}; extrap::AbstractExtrap=NoExtrap(), - side::Symbol=:nearest, + side::AbstractSide=NearestSide(), deriv::Int=0, search::AbstractSearchPolicy=Binary() ) where {Tg<:Real, Tv, Tq<:Real} @@ -361,7 +357,7 @@ function constant_interp!( y::AbstractVector{Tv}, x_targets::AbstractVector{Tq}; extrap::AbstractExtrap=NoExtrap(), - side::Symbol=:nearest, + side::AbstractSide=NearestSide(), deriv::Int=0, search::AbstractSearchPolicy=Binary() ) where {Tg<:Real, Tv, Tq<:Real} diff --git a/src/constant/constant_series_interp.jl b/src/constant/constant_series_interp.jl index d54b4925..21d73188 100644 --- a/src/constant/constant_series_interp.jl +++ b/src/constant/constant_series_interp.jl @@ -30,7 +30,7 @@ Shares a single x-grid across N y-series for efficient batch evaluation. - `y::Matrix{Tv}`: Function values (n_points × n_series) series-contiguous - `_transpose::LazyTranspose{Tv}`: Lazy point-contiguous layout for scalar SIMD - `extrap::AbstractExtrap`: Extrapolation mode -- `side::SideVal`: Side selection (:nearest, :left, :right) +- `side::SD`: Side selection (NearestSide(), LeftSide(), RightSide()) # Memory Layout Primary storage is series-contiguous (n_points × n_series): @@ -64,22 +64,22 @@ sitp_complex = constant_interp(x, y_complex) This type uses `mutable struct` with all `const` fields (Julia 1.8+) instead of plain `struct` for performance reasons. See CubicSeriesInterpolant for details. """ -mutable struct ConstantSeriesInterpolant{Tg<:AbstractFloat, Tv, E<:AbstractExtrap, P<:AbstractSearchPolicy, X<:AbstractVector{Tg}} <: AbstractSeriesInterpolant{Tg, Tv} +mutable struct ConstantSeriesInterpolant{Tg<:AbstractFloat, Tv, E<:AbstractExtrap, SD<:AbstractSide, P<:AbstractSearchPolicy, X<:AbstractVector{Tg}} <: AbstractSeriesInterpolant{Tg, Tv} const x::X # Shared x-grid (Range or Vector) const y::Matrix{Tv} # Series-contiguous y (n_points × n_series) const _transpose::LazyTranspose{Tv} # Lazy point-contiguous layout const extrap::E # Extrapolation mode (compile-time specialized) - const side::SideVal # Side selection + const side::SD # Side selection (compile-time specialized) const search_policy::P # Default search policy function ConstantSeriesInterpolant( x::X, y::Matrix{Tv}, extrap::E, - side::SideVal, + side::SD, search::P=Binary() - ) where {Tg<:AbstractFloat, Tv, E<:AbstractExtrap, P<:AbstractSearchPolicy, X<:AbstractVector{Tg}} - new{Tg,Tv,E,P,X}(x, y, LazyTranspose{Tv}(), extrap, side, search) + ) 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 end @@ -244,7 +244,7 @@ end x_max::Tg, aq::_ConstantAnchoredQuery{Tg}, ::NoExtrap, - ::SideVal, + ::AbstractSide, ::AbstractEvalOp, ::UInt8 ) where {Tg<:AbstractFloat, Tv} @@ -261,7 +261,7 @@ end ::Tg, ::_ConstantAnchoredQuery{Tg}, ::ConstExtrap, - ::SideVal, + ::AbstractSide, op::AbstractEvalOp, side::UInt8 ) where {Tg<:AbstractFloat, Tv} @@ -278,7 +278,7 @@ end ::Tg, aq::_ConstantAnchoredQuery{Tg}, ::ExtendExtrap, - side_val::SideVal, + side_val::AbstractSide, op::AbstractEvalOp, ::UInt8 ) where {Tg<:AbstractFloat, Tv} @@ -303,14 +303,14 @@ end # ======================================== """ - constant_interp(x, ys::AbstractVector{<:AbstractVector}; side=:nearest, extrap=NoExtrap()) + constant_interp(x, ys::AbstractVector{<:AbstractVector}; side=NearestSide(), extrap=NoExtrap()) Create a multi-Y constant interpolant for multiple y-data series sharing the same x-grid. # Arguments - `x::AbstractVector`: x-coordinates (sorted, length ≥ 2) - `ys`: Vector of y-value vectors (all same length as x) -- `side`: Side for discontinuities (:left, :right, :nearest) +- `side`: Side for discontinuities (LeftSide(), RightSide(), NearestSide()) - `extrap::AbstractExtrap`: `NoExtrap()`, `ConstExtrap()`, `ExtendExtrap()`, or `WrapExtrap()` # Returns @@ -331,7 +331,7 @@ vals = sitp(0.5) function constant_interp( x::AbstractVector{Tg}, ys::AbstractVector{<:AbstractVector{Tv}}; - side::Symbol=:nearest, + side::AbstractSide=NearestSide(), extrap::AbstractExtrap=NoExtrap(), search::P=Binary() ) where {Tg<:AbstractFloat, Tv, P<:AbstractSearchPolicy} @@ -367,14 +367,12 @@ function constant_interp( y_mat[:, k] .= Tv_out.(ys[k]) end - @_dispatch_side side => side_val begin - return ConstantSeriesInterpolant(x, y_mat, extrap, side_val, search) - end + return ConstantSeriesInterpolant(x, y_mat, extrap, side, search) end # Matrix input: columns as y-series """ - constant_interp(x, Y::AbstractMatrix; side=:nearest, extrap=NoExtrap()) + constant_interp(x, Y::AbstractMatrix; side=NearestSide(), extrap=NoExtrap()) Create a multi-Y constant interpolant from a matrix where each column is a y-series. @@ -394,7 +392,7 @@ sitp = constant_interp(x, Y) function constant_interp( x::AbstractVector{Tg}, Y::AbstractMatrix{Tv}; - side::Symbol=:nearest, + side::AbstractSide=NearestSide(), extrap::AbstractExtrap=NoExtrap(), search::AbstractSearchPolicy=Binary() ) where {Tg<:AbstractFloat, Tv} @@ -419,9 +417,7 @@ function constant_interp( Tv_out = _value_type(Tv, Tg) y_mat = Tv_out === Tv ? copy(Y) : Tv_out.(Y) - @_dispatch_side side => side_val begin - return ConstantSeriesInterpolant(x, y_mat, extrap, side_val, search) - end + return ConstantSeriesInterpolant(x, y_mat, extrap, side, search) end # ======================================== @@ -433,7 +429,7 @@ end function constant_interp( x::AbstractVector{Tg}, ys::AbstractVector{<:AbstractVector{Tv}}; - side::Symbol=:nearest, + side::AbstractSide=NearestSide(), extrap::AbstractExtrap=NoExtrap(), search::AbstractSearchPolicy=Binary() ) where {Tg<:Real, Tv} @@ -447,7 +443,7 @@ end function constant_interp( x::AbstractVector{Tg}, Y::AbstractMatrix{Tv}; - side::Symbol=:nearest, + side::AbstractSide=NearestSide(), extrap::AbstractExtrap=NoExtrap(), search::AbstractSearchPolicy=Binary() ) where {Tg<:Real, Tv} @@ -622,7 +618,7 @@ Uses argument-passing pattern for optimal performance. k::Int, aq_vec::AbstractVector{<:_ConstantAnchoredQuery{Tg}}, extrap::AbstractExtrap, - side_val::SideVal, + side_val::AbstractSide, op::AbstractEvalOp ) where {Tg<:AbstractFloat, Tv} @inbounds for j in eachindex(out, aq_vec) @@ -643,7 +639,7 @@ Internal: Evaluate single series at single query point with extrapolation handli k::Int, aq::_ConstantAnchoredQuery{Tg}, extrap::AbstractExtrap, - side_val::SideVal, + side_val::AbstractSide, op::AbstractEvalOp ) where {Tg<:AbstractFloat, Tv} # Special case: at right boundary (MUST be preserved!) @@ -677,7 +673,7 @@ Internal: Core constant evaluation for series k at anchored query point. y::Matrix{Tv}, k::Int, aq::_ConstantAnchoredQuery{Tg}, - side_val::SideVal, + side_val::AbstractSide, op::AbstractEvalOp ) where {Tg<:AbstractFloat, Tv} # Derivatives of constant (step) function are zero diff --git a/src/constant/constant_types.jl b/src/constant/constant_types.jl index 4670af11..ded6a338 100644 --- a/src/constant/constant_types.jl +++ b/src/constant/constant_types.jl @@ -6,7 +6,7 @@ # Internal evaluation functions are in constant_oneshot.jl. """ - ConstantInterpolant{Tg,Tv,X,Y,E,P} + ConstantInterpolant{Tg,Tv,X,Y,E,SD,P} Lightweight callable interpolant for constant (step) interpolation. Returned by `constant_interp(x, y)` (2-argument form). @@ -17,18 +17,19 @@ Returned by `constant_interp(x, y)` (2-argument form). - `X<:AbstractVector{Tg}`: Type of x-coordinates - `Y<:AbstractVector{Tv}`: Type of y-values - `E<:AbstractExtrap`: Extrapolation mode type (compile-time specialized) +- `SD<:AbstractSide`: Side selection type (compile-time specialized) - `P<:AbstractSearchPolicy`: Search policy type # Fields - `x::X`: x-coordinates (sorted) - `y::Y`: y-values - `extrap::E`: Extrapolation mode (NoExtrap(), ExtendExtrap(), ConstExtrap(), or WrapExtrap()) -- `side::SideVal`: Side selection (:nearest, :left, :right) +- `side::SD`: Side selection (NearestSide(), LeftSide(), RightSide()) - `search_policy::P`: Default search policy for interval lookup # Usage ```julia -itp = constant_interp(x, y) # default: extrap=NoExtrap(), side=:nearest +itp = constant_interp(x, y) # default: extrap=NoExtrap(), side=NearestSide() val = itp(0.5) # scalar evaluation vals = itp.(query_points) # broadcast @@ -39,8 +40,8 @@ itp = constant_interp(x, y) val = itp(0.5) # returns ComplexF64 # With options -itp_left = constant_interp(x, y; side=:left) -itp_wrap = constant_interp(x, y; extrap=WrapExtrap(), side=:right) +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()) @@ -48,20 +49,20 @@ val = itp(0.5) # uses LinearBinary() by default val = itp(0.5; search=Binary()) # override with Binary() ``` """ -struct ConstantInterpolant{Tg<:AbstractFloat, Tv, X<:AbstractVector{Tg}, Y<:AbstractVector{Tv}, E<:AbstractExtrap, P<:AbstractSearchPolicy} <: AbstractInterpolant{Tg, Tv} +struct ConstantInterpolant{Tg<:AbstractFloat, Tv, X<:AbstractVector{Tg}, Y<:AbstractVector{Tv}, E<:AbstractExtrap, SD<:AbstractSide, P<:AbstractSearchPolicy} <: AbstractInterpolant{Tg, Tv} x::X y::Y extrap::E # Extrapolation mode (compile-time specialized) - side::SideVal + side::SD # Side selection (compile-time specialized) search_policy::P # Default search policy (immutable, thread-safe) # Inner constructor: parametric, only calls new (handles validation only) - function ConstantInterpolant{Tg,Tv,X,Y,E,P}( - x::X, y::Y, ev::E, sv::SideVal, search::P - ) where {Tg<:AbstractFloat, Tv, X<:AbstractVector{Tg}, Y<:AbstractVector{Tv}, E<:AbstractExtrap, P<:AbstractSearchPolicy} + function ConstantInterpolant{Tg,Tv,X,Y,E,SD,P}( + x::X, y::Y, ev::E, sv::SD, search::P + ) where {Tg<:AbstractFloat, Tv, X<:AbstractVector{Tg}, Y<:AbstractVector{Tv}, E<:AbstractExtrap, SD<:AbstractSide, P<:AbstractSearchPolicy} @assert length(x) == length(y) "x and y must have same length" @assert length(x) >= 2 "x must have at least 2 elements" - new{Tg,Tv,X,Y,E,P}(x, y, ev, sv, search) + new{Tg,Tv,X,Y,E,SD,P}(x, y, ev, sv, search) end end @@ -76,11 +77,10 @@ end x::X, y::Y; extrap::AbstractExtrap=NoExtrap(), - side::Symbol=:nearest, + side::AbstractSide=NearestSide(), search::P=Binary() ) where {Tg<:AbstractFloat, Tv, X<:AbstractVector{Tg}, Y<:AbstractVector{Tv}, P<:AbstractSearchPolicy} E = typeof(extrap) - @_dispatch_side side => sv begin - return ConstantInterpolant{Tg,Tv,X,Y,E,P}(x, y, extrap, sv, search) - end + SD = typeof(side) + return ConstantInterpolant{Tg,Tv,X,Y,E,SD,P}(x, y, extrap, side, search) end diff --git a/src/constant/nd/constant_nd_eval.jl b/src/constant/nd/constant_nd_eval.jl index cccefc7d..62327300 100644 --- a/src/constant/nd/constant_nd_eval.jl +++ b/src/constant/nd/constant_nd_eval.jl @@ -235,7 +235,7 @@ Computes cell widths, distances from left edge, side-based offsets, and returns @generated function _constant_nd_kernel( data::AbstractArray{Tv, N}, spacings::NTuple{N, AbstractGridSpacing}, - sides::NTuple{N, SideVal}, + sides::Tuple{Vararg{AbstractSide, N}}, indices::NTuple{N, Int}, q_eval::Tuple{Vararg{Real, N}}, Ls::Tuple{Vararg{Real, N}} @@ -282,16 +282,16 @@ end # Side Offset Computation (dispatch helpers for @generated kernel) # ======================================== -@inline function _compute_single_offset(::Val{:left}, h, dL) +@inline function _compute_single_offset(::LeftSide, h, dL) return 0 end -@inline function _compute_single_offset(::Val{:right}, h, dL) +@inline function _compute_single_offset(::RightSide, h, dL) dL_primal = _extract_primal(dL) return iszero(dL_primal) ? 0 : 1 end -@inline function _compute_single_offset(::Val{:nearest}, h, dL) +@inline function _compute_single_offset(::NearestSide, h, dL) dL_primal = _extract_primal(dL) return dL_primal <= h / 2 ? 0 : 1 end diff --git a/src/constant/nd/constant_nd_interpolant.jl b/src/constant/nd/constant_nd_interpolant.jl index 3390e4cc..e358aa41 100644 --- a/src/constant/nd/constant_nd_interpolant.jl +++ b/src/constant/nd/constant_nd_interpolant.jl @@ -19,7 +19,7 @@ Create an N-dimensional constant interpolant with tuple-grid API. - `data`: N-dimensional data array where `size(data, d) == length(grids[d])` # Keyword Arguments -- `side=:nearest`: Side selection mode (`:nearest`, `:left`, `:right`) or per-axis tuple +- `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 @@ -40,7 +40,7 @@ vals = itp(points) # Batch AoS # Per-axis configuration itp = constant_interp((x, y), data; - side=(:left, :right), + side=(LeftSide(), RightSide()), extrap=(NoExtrap(), WrapExtrap()), search=(Binary(), LinearBinary()) ) @@ -49,7 +49,7 @@ itp = constant_interp((x, y), data; function constant_interp( grids::NTuple{N, AbstractVector}, data::AbstractArray{Tv_raw, N}; - side::Union{Symbol, NTuple{N, Symbol}} = :nearest, + side::Union{AbstractSide, Tuple{Vararg{AbstractSide}}} = NearestSide(), extrap::Union{AbstractExtrap, NTuple{N, AbstractExtrap}} = NoExtrap(), search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = Binary() ) where {N, Tv_raw} @@ -75,20 +75,8 @@ function constant_interp( searches = _resolve_search_nd(search, Val(N)) extrap_vals = _resolve_extrap_nd(extrap, nothing, Val(N)) - @_dispatch_side_nd sides => side_vals begin - return ConstantInterpolantND{Tg, Tv, N, - typeof(grids_typed), typeof(spacings), typeof(extrap_vals), typeof(side_vals), typeof(searches)}( - grids_typed, spacings, Array(data_typed), extrap_vals, side_vals, searches - ) - end + return ConstantInterpolantND{Tg, Tv, N, + typeof(grids_typed), typeof(spacings), typeof(extrap_vals), typeof(sides), typeof(searches)}( + grids_typed, spacings, Array(data_typed), extrap_vals, sides, searches + ) end - -# ======================================== -# Val Type Conversion -# ======================================== - -@inline function _to_side_vals(sides::NTuple{N, Symbol}) where {N} - return ntuple(i -> _to_side_val(sides[i]), Val(N)) -end - -@inline _to_side_val(s::Symbol) = Val(s) diff --git a/src/constant/nd/constant_nd_oneshot.jl b/src/constant/nd/constant_nd_oneshot.jl index 2562536b..d2330567 100644 --- a/src/constant/nd/constant_nd_oneshot.jl +++ b/src/constant/nd/constant_nd_oneshot.jl @@ -23,7 +23,7 @@ Evaluates directly from grids + data without constructing a ConstantInterpolantN data::AbstractArray{Tv, N}, query::Tuple{Vararg{Real, N}}, extraps_val::Tuple{Vararg{AbstractExtrap, N}}, - side_vals::NTuple{N, SideVal}, + side_vals::Tuple{Vararg{AbstractSide, N}}, searches::NTuple{N, AbstractSearchPolicy} ) where {Tg<:AbstractFloat, Tv, N} spacings = _create_spacings_pooled(pool, grids) @@ -44,7 +44,7 @@ Writes results into `output`. No heap allocation beyond spacings. data::AbstractArray{Tv, N}, queries::Tuple{Vararg{AbstractVector{<:Real}, N}}, extraps_val::Tuple{Vararg{AbstractExtrap, N}}, - side_vals::NTuple{N, SideVal}, + side_vals::Tuple{Vararg{AbstractSide, N}}, searches::NTuple{N, AbstractSearchPolicy} ) where {Tg<:AbstractFloat, Tv, N} n_queries = length(queries[1]) @@ -75,7 +75,7 @@ function _constant_interp_nd_oneshot_soa( data::AbstractArray{Tv, N}, queries::Tuple{Vararg{AbstractVector{<:Real}, N}}, extraps_val::Tuple{Vararg{AbstractExtrap, N}}, - side_vals::NTuple{N, SideVal}, + side_vals::Tuple{Vararg{AbstractSide, N}}, searches::NTuple{N, AbstractSearchPolicy} ) where {Tg<:AbstractFloat, Tv, N} output = Vector{Tv}(undef, length(queries[1])) @@ -94,7 +94,7 @@ Writes results into `output`. No heap allocation beyond spacings. data::AbstractArray{Tv, N}, queries::AbstractVector{<:Tuple{Vararg{Real, N}}}, extraps_val::Tuple{Vararg{AbstractExtrap, N}}, - side_vals::NTuple{N, SideVal}, + side_vals::Tuple{Vararg{AbstractSide, N}}, searches::NTuple{N, AbstractSearchPolicy} ) where {Tg<:AbstractFloat, Tv, N} n_queries = length(queries) @@ -120,7 +120,7 @@ function _constant_interp_nd_oneshot_aos( data::AbstractArray{Tv, N}, queries::AbstractVector{<:Tuple{Vararg{Real, N}}}, extraps_val::Tuple{Vararg{AbstractExtrap, N}}, - side_vals::NTuple{N, SideVal}, + side_vals::Tuple{Vararg{AbstractSide, N}}, searches::NTuple{N, AbstractSearchPolicy} ) where {Tg<:AbstractFloat, Tv, N} output = Vector{Tv}(undef, length(queries)) @@ -149,7 +149,7 @@ function constant_interp( grids::NTuple{N, AbstractVector}, data::AbstractArray{Tv, N}, query::Tuple{Vararg{Real, N}}; - side::Union{Symbol, NTuple{N, Symbol}} = :nearest, + side::Union{AbstractSide, Tuple{Vararg{AbstractSide}}} = NearestSide(), extrap::Union{AbstractExtrap, NTuple{N, AbstractExtrap}} = NoExtrap(), search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = Binary(), deriv::Union{Int, Val, NTuple{N,Int}} = 0 @@ -168,10 +168,8 @@ function constant_interp( searches = _resolve_search_nd(search, Val(N)) extraps_val = _resolve_extrap_nd(extrap, nothing, Val(N)) - @_dispatch_side_nd sides => side_vals begin - return _constant_interp_nd_oneshot( - grids_typed, data, query, extraps_val, side_vals, searches)::Tv - end + return _constant_interp_nd_oneshot( + grids_typed, data, query, extraps_val, sides, searches)::Tv end """ @@ -184,7 +182,7 @@ function constant_interp( grids::NTuple{N, AbstractVector}, data::AbstractArray{Tv, N}, queries::NTuple{N, AbstractVector{<:Real}}; - side::Union{Symbol, NTuple{N, Symbol}} = :nearest, + side::Union{AbstractSide, Tuple{Vararg{AbstractSide}}} = NearestSide(), extrap::Union{AbstractExtrap, NTuple{N, AbstractExtrap}} = NoExtrap(), search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = Binary(), deriv::Union{Int, Val, NTuple{N,Int}} = 0 @@ -203,10 +201,8 @@ function constant_interp( searches = _resolve_search_nd(search, Val(N)) extraps_val = _resolve_extrap_nd(extrap, nothing, Val(N)) - @_dispatch_side_nd sides => side_vals begin - return _constant_interp_nd_oneshot_soa( - grids_typed, data, queries, extraps_val, side_vals, searches)::Vector{Tv} - end + return _constant_interp_nd_oneshot_soa( + grids_typed, data, queries, extraps_val, sides, searches)::Vector{Tv} end """ @@ -219,7 +215,7 @@ function constant_interp( grids::NTuple{N, AbstractVector}, data::AbstractArray{Tv, N}, queries::AbstractVector{<:Tuple{Vararg{Real, N}}}; - side::Union{Symbol, NTuple{N, Symbol}} = :nearest, + side::Union{AbstractSide, Tuple{Vararg{AbstractSide}}} = NearestSide(), extrap::Union{AbstractExtrap, NTuple{N, AbstractExtrap}} = NoExtrap(), search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = Binary(), deriv::Union{Int, Val, NTuple{N,Int}} = 0 @@ -238,10 +234,8 @@ function constant_interp( searches = _resolve_search_nd(search, Val(N)) extraps_val = _resolve_extrap_nd(extrap, nothing, Val(N)) - @_dispatch_side_nd sides => side_vals begin - return _constant_interp_nd_oneshot_aos( - grids_typed, data, queries, extraps_val, side_vals, searches)::Vector{Tv} - end + return _constant_interp_nd_oneshot_aos( + grids_typed, data, queries, extraps_val, sides, searches)::Vector{Tv} end # ======================================== @@ -259,7 +253,7 @@ function constant_interp!( grids::NTuple{N, AbstractVector}, data::AbstractArray{Tv, N}, queries::NTuple{N, AbstractVector{<:Real}}; - side::Union{Symbol, NTuple{N, Symbol}} = :nearest, + side::Union{AbstractSide, Tuple{Vararg{AbstractSide}}} = NearestSide(), extrap::Union{AbstractExtrap, NTuple{N, AbstractExtrap}} = NoExtrap(), search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = Binary(), deriv::Union{Int, Val, NTuple{N,Int}} = 0 @@ -278,10 +272,8 @@ function constant_interp!( searches = _resolve_search_nd(search, Val(N)) extraps_val = _resolve_extrap_nd(extrap, nothing, Val(N)) - @_dispatch_side_nd sides => side_vals begin - return _constant_interp_nd_oneshot_soa!( - output, grids_typed, data, queries, extraps_val, side_vals, searches) - end + return _constant_interp_nd_oneshot_soa!( + output, grids_typed, data, queries, extraps_val, sides, searches) end """ @@ -295,7 +287,7 @@ function constant_interp!( grids::NTuple{N, AbstractVector}, data::AbstractArray{Tv, N}, queries::AbstractVector{<:Tuple{Vararg{Real, N}}}; - side::Union{Symbol, NTuple{N, Symbol}} = :nearest, + side::Union{AbstractSide, Tuple{Vararg{AbstractSide}}} = NearestSide(), extrap::Union{AbstractExtrap, NTuple{N, AbstractExtrap}} = NoExtrap(), search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = Binary(), deriv::Union{Int, Val, NTuple{N,Int}} = 0 @@ -314,8 +306,6 @@ function constant_interp!( searches = _resolve_search_nd(search, Val(N)) extraps_val = _resolve_extrap_nd(extrap, nothing, Val(N)) - @_dispatch_side_nd sides => side_vals begin - return _constant_interp_nd_oneshot_aos!( - output, grids_typed, data, queries, extraps_val, side_vals, searches) - end + return _constant_interp_nd_oneshot_aos!( + output, grids_typed, data, queries, extraps_val, sides, searches) end diff --git a/src/constant/nd/constant_nd_types.jl b/src/constant/nd/constant_nd_types.jl index 779469b5..dea93321 100644 --- a/src/constant/nd/constant_nd_types.jl +++ b/src/constant/nd/constant_nd_types.jl @@ -17,7 +17,7 @@ N-dimensional constant (step) interpolation with per-axis configuration. - `G<:NTuple{N, AbstractVector{Tg}}`: Grid tuple type - `S<:NTuple{N, AbstractGridSpacing{Tg}}`: Spacing tuple type - `E<:Tuple{Vararg{AbstractExtrap, N}}`: Extrapolation mode tuple type -- `SD<:NTuple{N, SideVal}`: Side selection tuple type +- `SD<:Tuple{Vararg{AbstractSide, N}}`: Side selection tuple type - `P<:NTuple{N, AbstractSearchPolicy}`: Search policy tuple type # Fields @@ -25,7 +25,7 @@ N-dimensional constant (step) interpolation with per-axis configuration. - `spacings`: Tuple of spacing objects for efficient interval lookup - `data`: N-dimensional data array - `extraps`: Per-axis extrapolation modes -- `sides`: Per-axis side selection (:nearest, :left, :right) +- `sides`: Per-axis side selection (NearestSide(), LeftSide(), RightSide()) - `searches`: Per-axis search policies # Usage @@ -35,14 +35,14 @@ x = [0.0, 1.0, 2.0] y = [0.0, 1.0, 2.0, 3.0] data = rand(3, 4) -itp = constant_interp((x, y), data) # Default: side=:nearest, extrap=NoExtrap() +itp = constant_interp((x, y), data) # Default: side=NearestSide(), extrap=NoExtrap() val = itp((0.5, 1.5)) # Single query # With configuration -itp = constant_interp((x, y), data; side=:left, extrap=ConstExtrap()) +itp = constant_interp((x, y), data; side=LeftSide(), extrap=ConstExtrap()) # Per-axis configuration -itp = constant_interp((x, y), data; side=(:left, :right), extrap=(NoExtrap(), WrapExtrap())) +itp = constant_interp((x, y), data; side=(LeftSide(), RightSide()), extrap=(NoExtrap(), WrapExtrap())) ``` """ struct ConstantInterpolantND{ @@ -52,7 +52,7 @@ struct ConstantInterpolantND{ G<:NTuple{N, AbstractVector{Tg}}, S<:NTuple{N, AbstractGridSpacing{Tg}}, E<:Tuple{Vararg{AbstractExtrap, N}}, - SD<:NTuple{N, SideVal}, + SD<:Tuple{Vararg{AbstractSide, N}}, P<:NTuple{N, AbstractSearchPolicy}, } <: AbstractInterpolantND{Tg, Tv, N} grids::G @@ -66,7 +66,7 @@ struct ConstantInterpolantND{ grids::G, spacings::S, data::Array{Tv,N}, extraps::E, sides::SD, searches::P ) where {Tg<:AbstractFloat, Tv, N, G<:NTuple{N,AbstractVector{Tg}}, S<:NTuple{N,AbstractGridSpacing{Tg}}, E, - SD<:NTuple{N,SideVal}, P<:NTuple{N,AbstractSearchPolicy}} + SD<:Tuple{Vararg{AbstractSide, N}}, P<:NTuple{N,AbstractSearchPolicy}} new{Tg,Tv,N,G,S,E,SD,P}(grids, spacings, data, extraps, sides, searches) end end diff --git a/src/core/bc_types.jl b/src/core/bc_types.jl index 896db019..d2d48352 100644 --- a/src/core/bc_types.jl +++ b/src/core/bc_types.jl @@ -723,7 +723,7 @@ end # ======================================== """ - materialize_bc(bc::PolyFit{D}, xs, ys, endpoint::Val{:left/:right}) -> Deriv1{Tv} + materialize_bc(bc::PolyFit{D}, xs, ys, endpoint::LeftSide/RightSide) -> Deriv1{Tv} Convert a polynomial-fit BC to a concrete `Deriv1` by estimating the derivative from data. @@ -734,7 +734,7 @@ allowing all existing `Deriv1` code paths to work unchanged. - `bc::PolyFit{D}`: Polynomial fit BC to materialize (lazy, no value type) - `xs::AbstractVector{Tg}`: Grid coordinates (real-valued) - `ys::AbstractVector{Tv}`: Function values at grid points (can be Complex) -- `endpoint::Val{:left}` or `Val{:right}`: Which endpoint to estimate +- `endpoint::LeftSide` or `::RightSide`: Which endpoint to estimate # Returns - `Deriv1{Tv}(estimated_value)`: Concrete first derivative BC with value type from ys @@ -742,26 +742,26 @@ allowing all existing `Deriv1` code paths to work unchanged. # Example ```julia bc = QuadraticFit() # = PolyFit{2} -concrete_bc = materialize_bc(bc, xs, ys, Val(:left)) # → Deriv1{Float64}(computed_value) +concrete_bc = materialize_bc(bc, xs, ys, LeftSide()) # → Deriv1{Float64}(computed_value) # Complex values y_complex = [1.0+2.0im, 3.0+4.0im, ...] -concrete_bc = materialize_bc(bc, xs, y_complex, Val(:left)) # → Deriv1{ComplexF64}(...) +concrete_bc = materialize_bc(bc, xs, y_complex, LeftSide()) # → Deriv1{ComplexF64}(...) ``` See also: [`PolyFit`](@ref), [`_estimate_endpoint_derivative`](@ref) """ @inline function materialize_bc( - ::PolyFit{D}, xs::AbstractVector{Tg}, ys::AbstractVector{Tv}, endpoint::Val + ::PolyFit{D}, xs::AbstractVector{Tg}, ys::AbstractVector{Tv}, endpoint::AbstractSide ) where {D, Tg<:AbstractFloat, Tv} val = _estimate_endpoint_derivative(xs, ys, endpoint, PolyFit{D}()) return Deriv1{Tv}(val) # Natural Deriv1{Tv} return - no workaround needed! end # Passthrough for already-concrete BCs (no materialization needed) -@inline materialize_bc(bc::Deriv1, ::AbstractVector, ::AbstractVector, ::Val) = bc -@inline materialize_bc(bc::Deriv2, ::AbstractVector, ::AbstractVector, ::Val) = bc -@inline materialize_bc(bc::Deriv3, ::AbstractVector, ::AbstractVector, ::Val) = bc +@inline materialize_bc(bc::Deriv1, ::AbstractVector, ::AbstractVector, ::AbstractSide) = bc +@inline materialize_bc(bc::Deriv2, ::AbstractVector, ::AbstractVector, ::AbstractSide) = bc +@inline materialize_bc(bc::Deriv3, ::AbstractVector, ::AbstractVector, ::AbstractSide) = bc # ======================================== diff --git a/src/core/eval_ops.jl b/src/core/eval_ops.jl index 439ff3a4..cc84973e 100644 --- a/src/core/eval_ops.jl +++ b/src/core/eval_ops.jl @@ -135,12 +135,69 @@ itp = cubic_interp((x, y), data; extrap=WrapExtrap()) """ struct WrapExtrap <: AbstractExtrap end +# ======================================== +# Typed Side Selection Tags (Constant Interpolation) +# ======================================== +# +# Compile-time type tags for side mode selection. +# Zero-cost type dispatch — no macros needed with proper types. +# +# 1D: Structs store SD<:AbstractSide, dispatch on concrete subtypes +# ND: Structs store Tuple{Vararg{AbstractSide, N}}, resolved via _resolve_side_nd + +""" + AbstractSide + +Abstract type for side selection mode in constant interpolation. +Determines which neighbor value to use at non-grid-point locations. + +# Concrete subtypes +- [`NearestSide`](@ref): Nearest neighbor with left tie-breaking at midpoint +- [`LeftSide`](@ref): Always use left (floor) value +- [`RightSide`](@ref): Use right (ceiling) value, except at grid points + +!!! info "Type dispatch mechanism" + Interpolant structs store `side` as a type parameter `SD<:AbstractSide`, + so dispatch is fully monomorphized at compile time — zero overhead. + For oneshot paths where `side` appears in a Union, Julia union-splits + automatically (3 concrete subtypes < 4 limit). +""" +abstract type AbstractSide end + """ - SideVal + NearestSide <: AbstractSide -Union type for side selection mode values (constant interpolation). -Using concrete Union enables Julia's union-splitting optimization. +Nearest-neighbor side selection with left tie-breaking at midpoint. +Returns left value if distance <= h/2, right value otherwise. -Valid values: `Val(:nearest)`, `Val(:left)`, `Val(:right)` +# Example +```julia +itp = constant_interp(x, y; side=NearestSide()) +``` +""" +struct NearestSide <: AbstractSide end + +""" + LeftSide <: AbstractSide + +Left-continuous (floor) side selection. Always returns the left boundary value. + +# Example +```julia +itp = constant_interp(x, y; side=LeftSide()) +``` +""" +struct LeftSide <: AbstractSide end + +""" + RightSide <: AbstractSide + +Right-continuous (ceiling) side selection. +Returns right value except at grid points (where it returns left value). + +# Example +```julia +itp = constant_interp(x, y; side=RightSide()) +``` """ -const SideVal = Union{Val{:nearest}, Val{:left}, Val{:right}} +struct RightSide <: AbstractSide end diff --git a/src/core/nd_utils.jl b/src/core/nd_utils.jl index 17da35fd..a920aa90 100644 --- a/src/core/nd_utils.jl +++ b/src/core/nd_utils.jl @@ -11,21 +11,6 @@ # - NTuple{N} → passthrough (with optional validation) # - Wrong-sized tuple → ArgumentError -# ======================================== -# Side Validation (for Constant ND) -# ======================================== - -""" - _validate_side(side::Symbol) -> Nothing - -Validate side selection symbol. Throws `ArgumentError` if invalid. - -Valid options: `:nearest`, `:left`, `:right` -""" -@inline function _validate_side(side::Symbol) - side in (:nearest, :left, :right) && return nothing - throw(ArgumentError("`side` must be :nearest, :left, or :right, got :$side")) -end # ======================================== # Extrapolation Resolution @@ -157,27 +142,17 @@ end # ======================================== """ - _resolve_side_nd(side, Val(N)) -> NTuple{N, Symbol} + _resolve_side_nd(side, Val(N)) -> NTuple{N, AbstractSide} Resolve side selection to canonical N-tuple. -- Single `Symbol` → broadcast to all N axes (validated) -- `NTuple{N, Symbol}` → validate each and passthrough - -Valid side options: `:nearest`, `:left`, `:right` +- Single `AbstractSide` → broadcast to all N axes +- `NTuple{N, AbstractSide}` → passthrough """ -@inline function _resolve_side_nd(side::Symbol, ::Val{N}) where {N} - _validate_side(side) - return ntuple(_ -> side, Val(N)) -end +@inline _resolve_side_nd(side::AbstractSide, ::Val{N}) where {N} = ntuple(_ -> side, Val(N)) -@inline function _resolve_side_nd(side::NTuple{N,Symbol}, ::Val{N}) where {N} - @inbounds for i in 1:N - _validate_side(side[i]) - end - return side -end +@inline _resolve_side_nd(side::NTuple{N,AbstractSide}, ::Val{N}) where {N} = side -@inline function _resolve_side_nd(side::Tuple{Vararg{Symbol}}, ::Val{N}) where {N} +@inline function _resolve_side_nd(side::Tuple{Vararg{AbstractSide}}, ::Val{N}) where {N} throw(ArgumentError("side tuple must have $N elements to match grid dimensions, got $(length(side))")) end diff --git a/src/core/polyfit_kernels.jl b/src/core/polyfit_kernels.jl index ac55d3c8..665a7d99 100644 --- a/src/core/polyfit_kernels.jl +++ b/src/core/polyfit_kernels.jl @@ -10,7 +10,7 @@ # # Design Philosophy: # - Single Source of Truth: All Lagrange derivative kernels live here -# - Type-Dispatched API: PolyFit{D} for degree, Val{:left/:right} for endpoint +# - Type-Dispatched API: PolyFit{D} for degree, LeftSide/RightSide for endpoint # - Two-Stage Pattern: Separate coefficient computation from application # Stage 1: Compute coefficients from grid (once per grid) # Stage 2: Apply coefficients to data (repeated for batch operations) @@ -19,7 +19,7 @@ # ======================================== # Type-Dispatched Kernel API # ======================================== -# Unified kernel interface using PolyFit{D} and Val{:left/:right} dispatch. +# Unified kernel interface using PolyFit{D} and LeftSide/RightSide dispatch. # These are the primary internal kernels - named for what they DO, not WHERE. # # API: @@ -97,7 +97,7 @@ end # _compute_deriv1: Uniform grid direct computation # ---------------------------------------- # Computes first derivative directly using known stencil coefficients. -# Dispatches on PolyFit{D} (degree) and Val{:left/:right} (endpoint). +# Dispatches on PolyFit{D} (degree) and LeftSide/RightSide (endpoint). """ _compute_deriv1(::PolyFit{D}, ::Val{Side}, f::NTuple, inv_h) -> T @@ -106,7 +106,7 @@ Compute first derivative on uniform grid using D+1 point stencil. # Arguments - `::PolyFit{D}`: Polynomial degree (D=1,2,3 supported) -- `::Val{:left}` or `::Val{:right}`: Which endpoint +- `::LeftSide` or `::RightSide`: Which endpoint - `f::NTuple{D+1,T}`: Function values at stencil points - `inv_h::T`: Inverse grid spacing (1/h) @@ -116,35 +116,35 @@ Compute first derivative on uniform grid using D+1 point stencil. - PolyFit{3} (CubicFit): 4 points, O(h³) accuracy """ # PolyFit{1} (LinearFit) - 2 points, O(h) -@inline function _compute_deriv1(::PolyFit{1}, ::Val{:left}, f::NTuple{2,T}, inv_h::T) where {T<:AbstractFloat} +@inline function _compute_deriv1(::PolyFit{1}, ::LeftSide, f::NTuple{2,T}, inv_h::T) where {T<:AbstractFloat} return (f[2] - f[1]) * inv_h end -@inline function _compute_deriv1(::PolyFit{1}, ::Val{:right}, f::NTuple{2,T}, inv_h::T) where {T<:AbstractFloat} +@inline function _compute_deriv1(::PolyFit{1}, ::RightSide, f::NTuple{2,T}, inv_h::T) where {T<:AbstractFloat} return (f[2] - f[1]) * inv_h # Same as left for linear end # PolyFit{2} (QuadraticFit) - 3 points, O(h²) -@inline function _compute_deriv1(::PolyFit{2}, ::Val{:left}, f::NTuple{3,T}, inv_h::T) where {T<:AbstractFloat} +@inline function _compute_deriv1(::PolyFit{2}, ::LeftSide, f::NTuple{3,T}, inv_h::T) where {T<:AbstractFloat} # Coefficients: (-3, 4, -1) / 2 coeff = inv_h / 2 return muladd(T(-3), f[1], muladd(T(4), f[2], -f[3])) * coeff end -@inline function _compute_deriv1(::PolyFit{2}, ::Val{:right}, f::NTuple{3,T}, inv_h::T) where {T<:AbstractFloat} +@inline function _compute_deriv1(::PolyFit{2}, ::RightSide, f::NTuple{3,T}, inv_h::T) where {T<:AbstractFloat} # Coefficients: (1, -4, 3) / 2 coeff = inv_h / 2 return muladd(T(1), f[1], muladd(T(-4), f[2], T(3) * f[3])) * coeff end # PolyFit{3} (CubicFit) - 4 points, O(h³) -@inline function _compute_deriv1(::PolyFit{3}, ::Val{:left}, f::NTuple{4,T}, inv_h::T) where {T<:AbstractFloat} +@inline function _compute_deriv1(::PolyFit{3}, ::LeftSide, f::NTuple{4,T}, inv_h::T) where {T<:AbstractFloat} # Coefficients: (-11, 18, -9, 2) / 6 coeff = inv_h / 6 return muladd(T(-11), f[1], muladd(T(18), f[2], muladd(T(-9), f[3], T(2) * f[4]))) * coeff end -@inline function _compute_deriv1(::PolyFit{3}, ::Val{:right}, f::NTuple{4,T}, inv_h::T) where {T<:AbstractFloat} +@inline function _compute_deriv1(::PolyFit{3}, ::RightSide, f::NTuple{4,T}, inv_h::T) where {T<:AbstractFloat} # Coefficients: (-2, 9, -18, 11) / 6 coeff = inv_h / 6 return muladd(T(-2), f[1], muladd(T(9), f[2], muladd(T(-18), f[3], T(11) * f[4]))) * coeff @@ -158,35 +158,35 @@ end # ---------------------------------------- # PolyFit{1} (LinearFit) - 2 points, O(h) - Mixed type -@inline function _compute_deriv1(::PolyFit{1}, ::Val{:left}, f::NTuple{2,Tv}, inv_h::Tg) where {Tv, Tg<:AbstractFloat} +@inline function _compute_deriv1(::PolyFit{1}, ::LeftSide, f::NTuple{2,Tv}, inv_h::Tg) where {Tv, Tg<:AbstractFloat} return (f[2] - f[1]) * inv_h # Tv * Tg → Tv end -@inline function _compute_deriv1(::PolyFit{1}, ::Val{:right}, f::NTuple{2,Tv}, inv_h::Tg) where {Tv, Tg<:AbstractFloat} +@inline function _compute_deriv1(::PolyFit{1}, ::RightSide, f::NTuple{2,Tv}, inv_h::Tg) where {Tv, Tg<:AbstractFloat} return (f[2] - f[1]) * inv_h end # PolyFit{2} (QuadraticFit) - 3 points, O(h²) - Mixed type -@inline function _compute_deriv1(::PolyFit{2}, ::Val{:left}, f::NTuple{3,Tv}, inv_h::Tg) where {Tv, Tg<:AbstractFloat} +@inline function _compute_deriv1(::PolyFit{2}, ::LeftSide, f::NTuple{3,Tv}, inv_h::Tg) where {Tv, Tg<:AbstractFloat} # Coefficients: (-3, 4, -1) / 2 coeff = inv_h / 2 return muladd(-3, f[1], muladd(4, f[2], -f[3])) * coeff # Int * Tv → Tv, Tv * Tg → Tv end -@inline function _compute_deriv1(::PolyFit{2}, ::Val{:right}, f::NTuple{3,Tv}, inv_h::Tg) where {Tv, Tg<:AbstractFloat} +@inline function _compute_deriv1(::PolyFit{2}, ::RightSide, f::NTuple{3,Tv}, inv_h::Tg) where {Tv, Tg<:AbstractFloat} # Coefficients: (1, -4, 3) / 2 coeff = inv_h / 2 return muladd(1, f[1], muladd(-4, f[2], 3 * f[3])) * coeff end # PolyFit{3} (CubicFit) - 4 points, O(h³) - Mixed type -@inline function _compute_deriv1(::PolyFit{3}, ::Val{:left}, f::NTuple{4,Tv}, inv_h::Tg) where {Tv, Tg<:AbstractFloat} +@inline function _compute_deriv1(::PolyFit{3}, ::LeftSide, f::NTuple{4,Tv}, inv_h::Tg) where {Tv, Tg<:AbstractFloat} # Coefficients: (-11, 18, -9, 2) / 6 coeff = inv_h / 6 return muladd(-11, f[1], muladd(18, f[2], muladd(-9, f[3], 2 * f[4]))) * coeff end -@inline function _compute_deriv1(::PolyFit{3}, ::Val{:right}, f::NTuple{4,Tv}, inv_h::Tg) where {Tv, Tg<:AbstractFloat} +@inline function _compute_deriv1(::PolyFit{3}, ::RightSide, f::NTuple{4,Tv}, inv_h::Tg) where {Tv, Tg<:AbstractFloat} # Coefficients: (-2, 9, -18, 11) / 6 coeff = inv_h / 6 return muladd(-2, f[1], muladd(9, f[2], muladd(-18, f[3], 11 * f[4]))) * coeff @@ -195,8 +195,8 @@ end # Generic fallback for D > 3 (uses barycentric differentiation) # Julia dispatch ensures D=1,2,3 use the specialized methods above. @inline @with_pool pool function _compute_deriv1( - pf::PolyFit{D}, side::Val{S}, f::NTuple{N,T}, inv_h::T -) where {D, S, N, T<:AbstractFloat} + pf::PolyFit{D}, side::AbstractSide, f::NTuple{N,T}, inv_h::T +) where {D, N, T<:AbstractFloat} # Compute coefficients on reference grid t = 0, 1, ..., D coeffs = acquire!(pool, T, N) β = acquire!(pool, T, N) @@ -228,22 +228,22 @@ Use with `_weighted_sum(coeffs, f_values)` to compute the derivative. # Arguments - `::PolyFit{D}`: Polynomial degree -- `::Val{:left}` or `::Val{:right}`: Which endpoint +- `::LeftSide` or `::RightSide`: Which endpoint - `x::NTuple{D+1,T}`: Grid coordinates at stencil points """ # PolyFit{1} (LinearFit) - 2 points -@inline function _compute_deriv1_coeffs(::PolyFit{1}, ::Val{:left}, x::NTuple{2,T}) where {T<:AbstractFloat} +@inline function _compute_deriv1_coeffs(::PolyFit{1}, ::LeftSide, x::NTuple{2,T}) where {T<:AbstractFloat} inv_dx = inv(x[2] - x[1]) return (-inv_dx, inv_dx) end -@inline function _compute_deriv1_coeffs(::PolyFit{1}, ::Val{:right}, x::NTuple{2,T}) where {T<:AbstractFloat} +@inline function _compute_deriv1_coeffs(::PolyFit{1}, ::RightSide, x::NTuple{2,T}) where {T<:AbstractFloat} inv_dx = inv(x[2] - x[1]) return (-inv_dx, inv_dx) # Same as left for linear end # PolyFit{2} (QuadraticFit) - 3 points -@inline function _compute_deriv1_coeffs(::PolyFit{2}, ::Val{:left}, x::NTuple{3,T}) where {T<:AbstractFloat} +@inline function _compute_deriv1_coeffs(::PolyFit{2}, ::LeftSide, x::NTuple{3,T}) where {T<:AbstractFloat} x1, x2, x3 = x d12, d13 = x1 - x2, x1 - x3 d23 = x2 - x3 @@ -253,7 +253,7 @@ end return (c1, c2, c3) end -@inline function _compute_deriv1_coeffs(::PolyFit{2}, ::Val{:right}, x::NTuple{3,T}) where {T<:AbstractFloat} +@inline function _compute_deriv1_coeffs(::PolyFit{2}, ::RightSide, x::NTuple{3,T}) where {T<:AbstractFloat} x1, x2, x3 = x d12, d13, d23 = x1 - x2, x1 - x3, x2 - x3 c1 = (-d23) / (d12 * d13) @@ -263,7 +263,7 @@ end end # PolyFit{3} (CubicFit) - 4 points -@inline function _compute_deriv1_coeffs(::PolyFit{3}, ::Val{:left}, x::NTuple{4,T}) where {T<:AbstractFloat} +@inline function _compute_deriv1_coeffs(::PolyFit{3}, ::LeftSide, x::NTuple{4,T}) where {T<:AbstractFloat} x1, x2, x3, x4 = x d12, d13, d14 = x1 - x2, x1 - x3, x1 - x4 d23, d24, d34 = x2 - x3, x2 - x4, x3 - x4 @@ -275,7 +275,7 @@ end return (c1, c2, c3, c4) end -@inline function _compute_deriv1_coeffs(::PolyFit{3}, ::Val{:right}, x::NTuple{4,T}) where {T<:AbstractFloat} +@inline function _compute_deriv1_coeffs(::PolyFit{3}, ::RightSide, x::NTuple{4,T}) where {T<:AbstractFloat} x1, x2, x3, x4 = x d12, d13, d14 = x1 - x2, x1 - x3, x1 - x4 d23, d24, d34 = x2 - x3, x2 - x4, x3 - x4 @@ -291,8 +291,8 @@ end # Returns NTuple for compatibility with _weighted_sum. # Julia dispatch ensures D=1,2,3 use the specialized methods above. @inline @with_pool pool function _compute_deriv1_coeffs( - pf::PolyFit{D}, side::Val{S}, x::NTuple{N,T} -) where {D, S, N, T<:AbstractFloat} + pf::PolyFit{D}, side::Union{LeftSide,RightSide}, x::NTuple{N,T} +) where {D, N, T<:AbstractFloat} c = acquire!(pool, T, N) β = acquire!(pool, T, N) _compute_deriv1_coeffs!(c, β, pf, side, x) @@ -409,7 +409,7 @@ This is the generic fallback for D > 3. For D = 1, 2, 3, the specialized - `coeffs::AbstractVector{T}`: Output buffer for coefficients (length ≥ D+1) - `β::AbstractVector{T}`: Workspace buffer for barycentric weights (length ≥ D+1) - `::PolyFit{D}`: Polynomial degree -- `side::Val{:left}` or `Val{:right}`: Which endpoint +- `side::LeftSide` or `::RightSide`: Which endpoint - `x::NTuple{N,T}`: Grid coordinates (length ≥ D+1) # Returns @@ -417,10 +417,9 @@ This is the generic fallback for D > 3. For D = 1, 2, 3, the specialized """ @inline function _compute_deriv1_coeffs!( coeffs::AbstractVector{T}, β::AbstractVector{T}, - ::PolyFit{D}, side::Val{S}, x::NTuple{N,T} -) where {D,S,N,T<:AbstractFloat} - @assert S === :left || S === :right "Invalid side: Val(:$S). Must be Val(:left) or Val(:right)." - k = (S === :left) ? 1 : N + ::PolyFit{D}, side::Union{LeftSide,RightSide}, x::NTuple{N,T} +) where {D,N,T<:AbstractFloat} + k = side isa LeftSide ? 1 : N _d1_coeffs_at_node!(coeffs, β, x, k, Val(N)) end @@ -428,10 +427,10 @@ end # ======================================== # Unified Endpoint Derivative Estimation API # ======================================== -# Higher-level API using Val{:left}/Val{:right} and PolyFit{D} dispatch. +# Higher-level API using LeftSide/RightSide and PolyFit{D} dispatch. # Handles both uniform (Range) and non-uniform (Vector) grids. # -# API signature: _estimate_endpoint_derivative(xs, ys, Val(:left/:right), PolyFit{D}()) +# API signature: _estimate_endpoint_derivative(xs, ys, LeftSide()/RightSide(), PolyFit{D}()) # # Supported polynomial degrees: # - PolyFit{1} (LinearFit): 2-point, O(h) accuracy @@ -444,16 +443,16 @@ end # ---------------------------------------- """ - _extract_stencil_values(v::AbstractVector, ::Val{:left}, ::Val{N}) -> NTuple{N} - _extract_stencil_values(v::AbstractVector, ::Val{:right}, ::Val{N}) -> NTuple{N} + _extract_stencil_values(v::AbstractVector, ::LeftSide, ::Val{N}) -> NTuple{N} + _extract_stencil_values(v::AbstractVector, ::RightSide, ::Val{N}) -> NTuple{N} Extract N values from left or right endpoint as NTuple for kernel dispatch. """ -@inline function _extract_stencil_values(v::AbstractVector{T}, ::Val{:left}, ::Val{N}) where {T, N} +@inline function _extract_stencil_values(v::AbstractVector{T}, ::LeftSide, ::Val{N}) where {T, N} @inbounds ntuple(i -> v[i], Val(N)) end -@inline function _extract_stencil_values(v::AbstractVector{T}, ::Val{:right}, ::Val{N}) where {T, N} +@inline function _extract_stencil_values(v::AbstractVector{T}, ::RightSide, ::Val{N}) where {T, N} n = length(v) @inbounds ntuple(i -> v[n - N + i], Val(N)) end @@ -522,7 +521,7 @@ Estimate first derivative at endpoint using D+1 point polynomial fit. # Arguments - `xs`: Grid coordinates (AbstractRange for uniform, AbstractVector for non-uniform) - `ys::AbstractVector{T}`: Function values (must have ≥ D+1 elements) -- `side`: `Val(:left)` or `Val(:right)` for endpoint selection +- `side`: `LeftSide()` or `RightSide()` for endpoint selection - `::PolyFit{D}`: Polynomial degree # Supported Degrees @@ -536,8 +535,8 @@ Estimate first derivative at endpoint using D+1 point polynomial fit. - D > 3: Dispatches to generic barycentric kernels (Vector-based) """ @inline function _estimate_endpoint_derivative( - xs::AbstractRange{T}, ys::AbstractVector{T}, side::Val{S}, pf::PolyFit{D} -) where {T<:AbstractFloat, S, D} + xs::AbstractRange{T}, ys::AbstractVector{T}, side::AbstractSide, pf::PolyFit{D} +) where {T<:AbstractFloat, D} _check_polyfit_requirements(D, length(ys)) @inbounds begin f = _extract_stencil_values(ys, side, Val(D + 1)) @@ -547,8 +546,8 @@ Estimate first derivative at endpoint using D+1 point polynomial fit. end @inline function _estimate_endpoint_derivative( - xs::AbstractVector{T}, ys::AbstractVector{T}, side::Val{S}, pf::PolyFit{D} -) where {T<:AbstractFloat, S, D} + xs::AbstractVector{T}, ys::AbstractVector{T}, side::AbstractSide, pf::PolyFit{D} +) where {T<:AbstractFloat, D} @assert length(xs) == length(ys) "xs and ys must have same length" _check_polyfit_requirements(D, length(ys)) @inbounds begin @@ -566,8 +565,8 @@ end # Returns Tv (same as value type) # ---------------------------------------- @inline function _estimate_endpoint_derivative( - xs::AbstractRange{Tg}, ys::AbstractVector{Tv}, side::Val{S}, pf::PolyFit{D} -) where {Tg<:AbstractFloat, Tv, S, D} + xs::AbstractRange{Tg}, ys::AbstractVector{Tv}, side::AbstractSide, pf::PolyFit{D} +) where {Tg<:AbstractFloat, Tv, D} _check_polyfit_requirements(D, length(ys)) @inbounds begin f = _extract_stencil_values(ys, side, Val(D + 1)) @@ -577,8 +576,8 @@ end end @inline function _estimate_endpoint_derivative( - xs::AbstractVector{Tg}, ys::AbstractVector{Tv}, side::Val{S}, pf::PolyFit{D} -) where {Tg<:AbstractFloat, Tv, S, D} + xs::AbstractVector{Tg}, ys::AbstractVector{Tv}, side::AbstractSide, pf::PolyFit{D} +) where {Tg<:AbstractFloat, Tv, D} @assert length(xs) == length(ys) "xs and ys must have same length" _check_polyfit_requirements(D, length(ys)) @inbounds begin diff --git a/src/core/show.jl b/src/core/show.jl index ff7b7299..ce9e4da9 100644 --- a/src/core/show.jl +++ b/src/core/show.jl @@ -111,13 +111,11 @@ function _format_extrap(mode) return "unknown" end -"""Format side selection from SideVal.""" -function _format_side(side) - side === Val(:nearest) && return ":nearest" - side === Val(:left) && return ":left" - side === Val(:right) && return ":right" - return "unknown" -end +"""Format side selection from AbstractSide.""" +_format_side(::NearestSide) = "NearestSide" +_format_side(::LeftSide) = "LeftSide" +_format_side(::RightSide) = "RightSide" +_format_side(side::AbstractSide) = string(typeof(side)) """Format search policy name.""" _format_search(::Binary) = "Binary" diff --git a/src/core/utils.jl b/src/core/utils.jl index 083cc7ce..2972742c 100644 --- a/src/core/utils.jl +++ b/src/core/utils.jl @@ -377,107 +377,4 @@ macro _dispatch_deriv(pair, body) end end -""" - @_dispatch_side sym => varname body - -Dispatch on runtime side symbol, executing body with concrete Val type. - -Converts `side::Symbol` (`:nearest`, `:left`, `:right`) to compile-time constant -`Val(:nearest)`, `Val(:left)`, or `Val(:right)`. This creates a function barrier -ensuring type stability and enabling union-splitting optimization. - -# Arguments -- `sym => varname`: Pair of symbol variable and binding name for Val type -- `body`: Expression to execute with `varname` bound to concrete Val - -# Example -```julia -@_dispatch_side side => sv begin - _constant_kernel(op, y_left, y_right, h, dL, sv) -end -``` - -Expands to: -```julia -let _side = side - if _side === :nearest - sv = Val(:nearest) - _constant_kernel(op, y_left, y_right, h, dL, sv) - elseif _side === :left - ... - end -end -``` -""" -macro _dispatch_side(pair, body) - # Parse pair: side => sv becomes Expr(:call, :(=>), :side, :sv) - pair.head === :call && pair.args[1] === :(=>) || - error("@_dispatch_side expects `sym => varname`, got: $pair") - sym = pair.args[2] - varname = pair.args[3] - svs = esc(varname) - quote - let _side = $(esc(sym)) - if _side === :nearest - $svs = Val(:nearest) - $(esc(body)) - elseif _side === :left - $svs = Val(:left) - $(esc(body)) - elseif _side === :right - $svs = Val(:right) - $(esc(body)) - else - throw(ArgumentError("`side` must be :nearest, :left, or :right, got :$_side")) - end - end - end -end - -# ======================================== -# ND Side Dispatch -# ======================================== - -@inline _is_uniform_side(sides::NTuple{N, Symbol}) where {N} = all(==(sides[1]), sides) - -""" - @_dispatch_side_nd sides => side_vals begin ... end - -Compile-time dispatch for ND side selection. Converts `NTuple{N, Symbol}` into -concrete `NTuple{N, Val{:...}}` to avoid Union boxing at function boundaries. - -Handles the common uniform case (all sides equal) with 3 branches. -Falls back to `_to_side_vals` for the rare mixed-side case. -""" -macro _dispatch_side_nd(pair, body) - pair.head === :call && pair.args[1] === :(=>) || - error("@_dispatch_side_nd expects `sides => varname`, got: $pair") - sides_expr = pair.args[2] - varname = pair.args[3] - svs = esc(varname) - - quote - let _sides_nd = $(esc(sides_expr)), _valn_side = Val(length($(esc(sides_expr)))) - if _is_uniform_side(_sides_nd) - if _sides_nd[1] === :nearest - let $svs = ntuple(_ -> Val(:nearest), _valn_side) - $(esc(body)) - end - elseif _sides_nd[1] === :left - let $svs = ntuple(_ -> Val(:left), _valn_side) - $(esc(body)) - end - else - let $svs = ntuple(_ -> Val(:right), _valn_side) - $(esc(body)) - end - end - else - let $svs = _to_side_vals(_sides_nd) - $(esc(body)) - end - end - end - end -end diff --git a/src/cubic/cubic_solver.jl b/src/cubic/cubic_solver.jl index 0de65fca..22bc1ac9 100644 --- a/src/cubic/cubic_solver.jl +++ b/src/cubic/cubic_solver.jl @@ -278,7 +278,7 @@ end d::AbstractVector{Tv}, bc::PolyFit{D}, y::AbstractVector{Tv}, x::AbstractVector{Tg}, spacing::AbstractGridSpacing{Tg} ) where {D, Tg<:AbstractFloat, Tv} # Materialize PolyFit{D} → Deriv1 using estimated derivative - concrete_bc = materialize_bc(bc, x, y, Val(:left)) + concrete_bc = materialize_bc(bc, x, y, LeftSide()) # Delegate to existing Deriv1 code path _compute_rhs_first!(d, concrete_bc, y, x, spacing) return nothing @@ -289,7 +289,7 @@ end d::AbstractVector{Tv}, bc::PolyFit{D}, y::AbstractVector{Tv}, x::AbstractVector{Tg}, spacing::AbstractGridSpacing{Tg} ) where {D, Tg<:AbstractFloat, Tv} # Materialize PolyFit{D} → Deriv1 using estimated derivative - concrete_bc = materialize_bc(bc, x, y, Val(:right)) + concrete_bc = materialize_bc(bc, x, y, RightSide()) # Delegate to existing Deriv1 code path _compute_rhs_last!(d, concrete_bc, y, x, spacing) return nothing diff --git a/src/cubic/nd/cubic_nd_math.jl b/src/cubic/nd/cubic_nd_math.jl index 69f3a515..dad3c540 100644 --- a/src/cubic/nd/cubic_nd_math.jl +++ b/src/cubic/nd/cubic_nd_math.jl @@ -319,8 +319,8 @@ end n = length(values) @assert n >= 4 "Need at least 4 points for CubicFit" - deriv_left = _estimate_endpoint_derivative(grid, values, Val(:left), CubicFit()) - deriv_right = _estimate_endpoint_derivative(grid, values, Val(:right), CubicFit()) + deriv_left = _estimate_endpoint_derivative(grid, values, LeftSide(), CubicFit()) + deriv_right = _estimate_endpoint_derivative(grid, values, RightSide(), CubicFit()) # BC values are Tv type (can be Complex) bc = BCPair(Deriv1(Tv(deriv_left)), Deriv1(Tv(deriv_right))) diff --git a/src/integral/integrate_api.jl b/src/integral/integrate_api.jl index 58f78e08..423dd03c 100644 --- a/src/integral/integrate_api.jl +++ b/src/integral/integrate_api.jl @@ -105,24 +105,15 @@ end ) where {Tg<:AbstractFloat, Tv} Tout = promote_type(Tv, Tg, typeof(x0), typeof(x1)) searcher = _to_searcher(search, hint) - # Manual if-else on side to avoid >4 union combinations. - # Each branch has concrete side → _integrate_constant_1d_impl is fully typed. - side = itp.side - if side === Val(:left) - return _integrate_constant_1d_impl(itp.x, itp.y, Val(:left), itp.extrap, x0, x1, searcher, Tg, Tout) - elseif side === Val(:right) - return _integrate_constant_1d_impl(itp.x, itp.y, Val(:right), itp.extrap, x0, x1, searcher, Tg, Tout) - else - return _integrate_constant_1d_impl(itp.x, itp.y, Val(:nearest), itp.extrap, x0, x1, searcher, Tg, Tout) - end + return _integrate_constant_1d_impl(itp.x, itp.y, itp.side, itp.extrap, x0, x1, searcher, Tg, Tout) end -# Receives concrete side Val{S} + extrap mode. +# Receives concrete side (parametric from struct) + extrap mode. # Uses the generic _integrate_1d_cellwise path — side is already concrete here. @inline function _integrate_constant_1d_impl( - x::AbstractVector, y::AbstractVector, side::Val{S}, extrap::AbstractExtrap, + x::AbstractVector, y::AbstractVector, side::AbstractSide, extrap::AbstractExtrap, x0::Real, x1::Real, searcher::SR, ::Type{Tg}, ::Type{Tout} -) where {S, SR<:Searcher, Tg, Tout} +) where {SR<:Searcher, Tg, Tout} partial = @inline (i, xL, h, a2, b2) -> begin @inbounds _constant_integral_kernel( _EvalIntegralPartial(), y[i], y[i+1], h, a2 - xL, b2 - xL, side @@ -234,20 +225,13 @@ end ) where {Tg<:AbstractFloat, Tv} Tout = promote_type(Tv, Tg, typeof(x0), typeof(x1)) searcher = _to_searcher(search, hint) - side = sitp.side - if side === Val(:left) - return _integrate_constant_series_1d(sitp.x, sitp.y, Val(:left), sitp.extrap, x0, x1, searcher, Tg, Tout) - elseif side === Val(:right) - return _integrate_constant_series_1d(sitp.x, sitp.y, Val(:right), sitp.extrap, x0, x1, searcher, Tg, Tout) - else - return _integrate_constant_series_1d(sitp.x, sitp.y, Val(:nearest), sitp.extrap, x0, x1, searcher, Tg, Tout) - end + return _integrate_constant_series_1d(sitp.x, sitp.y, sitp.side, sitp.extrap, x0, x1, searcher, Tg, Tout) end @inline function _integrate_constant_series_1d( - x::AbstractVector, y::AbstractMatrix, side::Val{S}, extrap::AbstractExtrap, + x::AbstractVector, y::AbstractMatrix, side::AbstractSide, extrap::AbstractExtrap, x0::Real, x1::Real, searcher::SR, ::Type{Tg}, ::Type{Tout} -) where {S, SR<:Searcher, Tg, Tout} +) where {SR<:Searcher, Tg, Tout} n = size(y, 2) results = Vector{Tout}(undef, n) @inbounds for k in 1:n diff --git a/src/integral/integrate_fulldomain.jl b/src/integral/integrate_fulldomain.jl index dcb88eec..e2496533 100644 --- a/src/integral/integrate_fulldomain.jl +++ b/src/integral/integrate_fulldomain.jl @@ -26,7 +26,7 @@ end @inline (i, h) -> @inbounds _quadratic_integral_kernel(_EvalIntegralCell(), a[i], d[i], y[i], h) end -@inline function _full_cell_fn(itp::ConstantInterpolant{Tg}, side::Val{S}) where {Tg, S} +@inline function _full_cell_fn(itp::ConstantInterpolant{Tg}, side::AbstractSide) where {Tg} y = itp.y @inline (i, h) -> @inbounds _constant_integral_kernel(_EvalIntegralPartial(), y[i], y[i+1], h, zero(Tg), h, side) end @@ -48,7 +48,7 @@ end @inline (i, h) -> @inbounds _quadratic_integral_kernel(_EvalIntegralCell(), a[i,k], d[i,k], y[i,k], h) end -@inline function _full_cell_fn(sitp::ConstantSeriesInterpolant{Tg}, k::Int, side::Val{S}) where {Tg, S} +@inline function _full_cell_fn(sitp::ConstantSeriesInterpolant{Tg}, k::Int, side::AbstractSide) where {Tg} y = sitp.y @inline (i, h) -> @inbounds _constant_integral_kernel(_EvalIntegralPartial(), y[i,k], y[i+1,k], h, zero(Tg), h, side) end @@ -67,21 +67,14 @@ end return _integrate_1d_fulldomain(x, _full_cell_fn(itp), Tout) end -# Constant override: resolve side union before entering loop +# Constant override: side is parametric → compiler knows concrete type @inline function integrate( itp::ConstantInterpolant{Tg,Tv}; search=nothing, hint=nothing ) where {Tg<:AbstractFloat, Tv} x = itp.x Tout = promote_type(Tv, Tg) - side = itp.side - if side === Val(:left) - return _integrate_1d_fulldomain(x, _full_cell_fn(itp, Val(:left)), Tout) - elseif side === Val(:right) - return _integrate_1d_fulldomain(x, _full_cell_fn(itp, Val(:right)), Tout) - else - return _integrate_1d_fulldomain(x, _full_cell_fn(itp, Val(:nearest)), Tout) - end + return _integrate_1d_fulldomain(x, _full_cell_fn(itp, itp.side), Tout) end # ═══════════════════════════════════════════════════════════════ @@ -103,7 +96,7 @@ end return results end -# Constant Series override: resolve side union +# Constant Series override: side is parametric → compiler knows concrete type @inline function integrate( sitp::ConstantSeriesInterpolant{Tg,Tv}; search=nothing, hint=nothing @@ -112,19 +105,8 @@ end Tout = promote_type(Tv, Tg) n = n_series(sitp) results = Vector{Tout}(undef, n) - side = sitp.side - if side === Val(:left) - @inbounds for k in 1:n - results[k] = _integrate_1d_fulldomain(x, _full_cell_fn(sitp, k, Val(:left)), Tout) - end - elseif side === Val(:right) - @inbounds for k in 1:n - results[k] = _integrate_1d_fulldomain(x, _full_cell_fn(sitp, k, Val(:right)), Tout) - end - else - @inbounds for k in 1:n - results[k] = _integrate_1d_fulldomain(x, _full_cell_fn(sitp, k, Val(:nearest)), Tout) - end + @inbounds for k in 1:n + results[k] = _integrate_1d_fulldomain(x, _full_cell_fn(sitp, k, sitp.side), Tout) end return results end @@ -195,20 +177,13 @@ function cumulative_integrate( return _cumulative_integrate_1d(x, _full_cell_fn(itp), Tout) end -# Constant override: resolve side union +# Constant override: side is parametric → compiler knows concrete type function cumulative_integrate( itp::ConstantInterpolant{Tg,Tv} ) where {Tg<:AbstractFloat, Tv} x = itp.x Tout = promote_type(Tv, Tg) - side = itp.side - if side === Val(:left) - return _cumulative_integrate_1d(x, _full_cell_fn(itp, Val(:left)), Tout) - elseif side === Val(:right) - return _cumulative_integrate_1d(x, _full_cell_fn(itp, Val(:right)), Tout) - else - return _cumulative_integrate_1d(x, _full_cell_fn(itp, Val(:nearest)), Tout) - end + return _cumulative_integrate_1d(x, _full_cell_fn(itp, itp.side), Tout) end # Generic Series: catches Cubic, Linear, Quadratic series @@ -226,7 +201,7 @@ function cumulative_integrate( return result end -# Constant Series override: resolve side union +# Constant Series override: side is parametric → compiler knows concrete type function cumulative_integrate( sitp::ConstantSeriesInterpolant{Tg,Tv} ) where {Tg<:AbstractFloat, Tv} @@ -235,16 +210,8 @@ function cumulative_integrate( n_pts = length(x) n_ser = n_series(sitp) result = Matrix{Tout}(undef, n_pts, n_ser) - side = sitp.side - if side === Val(:left) - _side = Val(:left) - elseif side === Val(:right) - _side = Val(:right) - else - _side = Val(:nearest) - end @inbounds for k in 1:n_ser - _cumulative_integrate_1d!(@view(result[:, k]), x, _full_cell_fn(sitp, k, _side)) + _cumulative_integrate_1d!(@view(result[:, k]), x, _full_cell_fn(sitp, k, sitp.side)) end return result end diff --git a/src/integral/integrate_kernels.jl b/src/integral/integrate_kernels.jl index 6cc73c99..b76a0da3 100644 --- a/src/integral/integrate_kernels.jl +++ b/src/integral/integrate_kernels.jl @@ -87,29 +87,29 @@ end # ═══════════════════════════════════════════════════════════════ # Constant integration kernels -# Piecewise constant: value depends on side mode (:left, :right, :nearest) +# Piecewise constant: value depends on side mode (LeftSide, RightSide, NearestSide) # ═══════════════════════════════════════════════════════════════ -# --- :left side — always use left value yL --- +# --- LeftSide — always use left value yL --- @inline function _constant_integral_kernel( ::_EvalIntegralPartial, - yL::Tv, yR::Tv, h::Tg, u0::Td, u1::Td, ::Val{:left} + yL::Tv, yR::Tv, h::Tg, u0::Td, u1::Td, ::LeftSide ) where {Tv, Tg<:AbstractFloat, Td<:Real} return yL * (u1 - u0) end -# --- :right side — always use right value yR --- +# --- RightSide — always use right value yR --- @inline function _constant_integral_kernel( ::_EvalIntegralPartial, - yL::Tv, yR::Tv, h::Tg, u0::Td, u1::Td, ::Val{:right} + yL::Tv, yR::Tv, h::Tg, u0::Td, u1::Td, ::RightSide ) where {Tv, Tg<:AbstractFloat, Td<:Real} return yR * (u1 - u0) end -# --- :nearest side — split at midpoint h/2 --- +# --- NearestSide — split at midpoint h/2 --- @inline function _constant_integral_kernel( ::_EvalIntegralPartial, - yL::Tv, yR::Tv, h::Tg, u0::Td, u1::Td, ::Val{:nearest} + yL::Tv, yR::Tv, h::Tg, u0::Td, u1::Td, ::NearestSide ) where {Tv, Tg<:AbstractFloat, Td<:Real} mid = h / 2 if u1 <= mid @@ -372,17 +372,17 @@ end # ND Constant integration kernel # # Per-axis weight functions for side-dependent integration: -# :left → all weight to left corner -# :right → all weight to right corner -# :nearest → split at midpoint h/2 +# LeftSide → all weight to left corner +# RightSide → all weight to right corner +# NearestSide → split at midpoint h/2 # ═══════════════════════════════════════════════════════════════ -@inline _cw0(u0, u1, h, ::Val{:left}) = u1 - u0 -@inline _cw1(u0, u1, h, ::Val{:left}) = zero(u1 - u0) -@inline _cw0(u0, u1, h, ::Val{:right}) = zero(u1 - u0) -@inline _cw1(u0, u1, h, ::Val{:right}) = u1 - u0 -@inline _cw0(u0, u1, h, ::Val{:nearest}) = max(zero(u0), min(u1, h / 2) - u0) -@inline _cw1(u0, u1, h, ::Val{:nearest}) = max(zero(u0), u1 - max(u0, h / 2)) +@inline _cw0(u0, u1, h, ::LeftSide) = u1 - u0 +@inline _cw1(u0, u1, h, ::LeftSide) = zero(u1 - u0) +@inline _cw0(u0, u1, h, ::RightSide) = zero(u1 - u0) +@inline _cw1(u0, u1, h, ::RightSide) = u1 - u0 +@inline _cw0(u0, u1, h, ::NearestSide) = max(zero(u0), min(u1, h / 2) - u0) +@inline _cw1(u0, u1, h, ::NearestSide) = max(zero(u0), u1 - max(u0, h / 2)) @inline @generated function _integrate_constant_nd_cell( data::Array{Tv, N}, @@ -390,7 +390,7 @@ end hs::NTuple{N}, ulo::NTuple{N}, uhi::NTuple{N}, - sides::NTuple{N, SideVal} + sides::Tuple{Vararg{AbstractSide, N}} ) where {Tv, N} nc = 1 << N terms = Expr[] diff --git a/src/quadratic/quadratic_solver.jl b/src/quadratic/quadratic_solver.jl index 71612d33..c6d7a343 100644 --- a/src/quadratic/quadratic_solver.jl +++ b/src/quadratic/quadratic_solver.jl @@ -237,7 +237,7 @@ For Complex y values, materialize_bc returns Deriv1{ComplexF64} naturally. @inline function _fill_slopes!(d::AbstractVector{Tv}, s::AbstractVector{Tv}, h::AbstractVector{Tg}, bc::Left{PolyFit{D}}, x::AbstractVector{Tg}, y::AbstractVector{Tv}) where {D, Tv, Tg<:AbstractFloat} # Materialize PolyFit{D} → Deriv1{Tv} using estimated derivative - concrete_bc = materialize_bc(bc.bc, x, y, Val(:left)) + concrete_bc = materialize_bc(bc.bc, x, y, LeftSide()) d1 = concrete_bc.val # Already Tv type from polynomial fit on y values _forward_recurrence!(d, s, d1) end @@ -253,7 +253,7 @@ estimated derivative directly. @inline function _fill_slopes!(d::AbstractVector{Tv}, s::AbstractVector{Tv}, h::AbstractVector{Tg}, bc::Right{PolyFit{D}}, x::AbstractVector{Tg}, y::AbstractVector{Tv}) where {D, Tv, Tg<:AbstractFloat} # Materialize PolyFit{D} → Deriv1{Tv} using estimated derivative - concrete_bc = materialize_bc(bc.bc, x, y, Val(:right)) + concrete_bc = materialize_bc(bc.bc, x, y, RightSide()) dn = concrete_bc.val # Already Tv type from polynomial fit on y values _backward_recurrence!(d, s, dn) end diff --git a/test/test_autodiff_Enzyme.jl b/test/test_autodiff_Enzyme.jl index 2e99d00b..086c4232 100644 --- a/test/test_autodiff_Enzyme.jl +++ b/test/test_autodiff_Enzyme.jl @@ -128,7 +128,7 @@ else y = [10.0, 20.0, 30.0, 40.0, 50.0, 60.0] @testset "Single Interpolant" begin - itp = constant_interp(x, y; side=:left, extrap=ExtendExtrap()) + itp = constant_interp(x, y; side=LeftSide(), extrap=ExtendExtrap()) f(xq) = itp(xq) result = Enzyme.autodiff(Enzyme.Reverse, f, Enzyme.Active, Enzyme.Active(2.5)) diff --git a/test/test_autodiff_ForwardDiff.jl b/test/test_autodiff_ForwardDiff.jl index 72f0489b..eeffaaae 100644 --- a/test/test_autodiff_ForwardDiff.jl +++ b/test/test_autodiff_ForwardDiff.jl @@ -279,7 +279,7 @@ const FI = FastInterpolations y = [10.0, 20.0, 30.0, 40.0, 50.0, 60.0] @testset "interpolant callable AD - derivative is zero" begin - itp = constant_interp(x, y; side=:left, extrap=ExtendExtrap()) + itp = constant_interp(x, y; side=LeftSide(), extrap=ExtendExtrap()) # Constant interpolation derivative should be 0 (step function) for xq in [0.5, 1.5, 2.5, 3.5] @@ -293,12 +293,12 @@ const FI = FastInterpolations @testset "one-shot API AD" begin # Verify one-shot API works with ForwardDiff xq = 2.5 - fd_deriv = ForwardDiff.derivative(q -> constant_interp(x, y, q; side=:left), xq) + fd_deriv = ForwardDiff.derivative(q -> constant_interp(x, y, q; side=LeftSide()), xq) @test fd_deriv ≈ 0.0 atol=1e-10 end @testset "value is preserved" begin - itp = constant_interp(x, y; side=:left, extrap=ExtendExtrap()) + itp = constant_interp(x, y; side=LeftSide(), extrap=ExtendExtrap()) test_points = [0.5, 1.5, 2.5, 3.5] for xq in test_points @@ -310,7 +310,7 @@ const FI = FastInterpolations @testset "side modes" begin # Test all side modes work with AD - for side_mode in [:left, :right, :nearest] + for side_mode in [LeftSide(), RightSide(), NearestSide()] itp = constant_interp(x, y; side=side_mode, extrap=ExtendExtrap()) fd_deriv = ForwardDiff.derivative(itp, 2.5) @test fd_deriv ≈ 0.0 atol=1e-10 @@ -326,7 +326,7 @@ const FI = FastInterpolations x = collect(0.0:1.0:5.0) y1 = [10.0, 20.0, 30.0, 40.0, 50.0, 60.0] y2 = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0] - sitp = constant_interp(x, [y1, y2]; side=:left, extrap=ExtendExtrap()) + sitp = constant_interp(x, [y1, y2]; side=LeftSide(), extrap=ExtendExtrap()) @testset "series derivative is zero" begin # ForwardDiff.derivative on series returns vector of derivatives @@ -356,7 +356,7 @@ const FI = FastInterpolations x = collect(0.0:1.0:5.0) y_complex = [10.0+1.0im, 20.0+2.0im, 30.0+3.0im, 40.0+4.0im, 50.0+5.0im, 60.0+6.0im] - itp = constant_interp(x, y_complex; side=:left, extrap=ExtendExtrap()) + itp = constant_interp(x, y_complex; side=LeftSide(), extrap=ExtendExtrap()) @testset "complex derivative" begin xq = 2.5 @@ -378,7 +378,7 @@ const FI = FastInterpolations @testset "Constant gradient composition" begin x = collect(0.0:1.0:5.0) y = [10.0, 20.0, 30.0, 40.0, 50.0, 60.0] - itp = constant_interp(x, y; side=:left, extrap=ExtendExtrap()) + itp = constant_interp(x, y; side=LeftSide(), extrap=ExtendExtrap()) @testset "gradient through loss function" begin function loss(params) diff --git a/test/test_autodiff_Zygote.jl b/test/test_autodiff_Zygote.jl index 35923407..3c862fb4 100644 --- a/test/test_autodiff_Zygote.jl +++ b/test/test_autodiff_Zygote.jl @@ -144,7 +144,7 @@ end @testset "Single Interpolant - gradient is nothing (not 0)" begin # Zygote returns `nothing` for constant functions, not 0 - itp = constant_interp(x, y; side=:left, extrap=ExtendExtrap()) + itp = constant_interp(x, y; side=LeftSide(), extrap=ExtendExtrap()) for xq in [0.5, 1.5, 2.5, 3.5] zy_grad = Zygote.gradient(itp, xq)[1] @@ -154,7 +154,7 @@ end end @testset "All side modes" begin - for side_mode in [:left, :right, :nearest] + for side_mode in [LeftSide(), RightSide(), NearestSide()] itp = constant_interp(x, y; side=side_mode, extrap=ExtendExtrap()) zy_grad = Zygote.gradient(itp, 2.5)[1] @test (zy_grad === nothing) || isapprox(zy_grad, 0.0; atol=1e-10) diff --git a/test/test_coeffs.jl b/test/test_coeffs.jl index 426b310c..f280bf7a 100644 --- a/test/test_coeffs.jl +++ b/test/test_coeffs.jl @@ -214,8 +214,8 @@ using FastInterpolations # Constant coeffs # ──────────────────────────────────────────── @testset "constant coeffs" begin - @testset "side=:left" begin - itp = constant_interp(x_vec, y_const; side=:left) + @testset "side=LeftSide()" begin + itp = constant_interp(x_vec, y_const; side=LeftSide()) for xq in [0.05, 0.5, 1.5, 2.5] cell = coeffs(itp, xq) @test cell(xq) ≈ itp(xq) atol=1e-12 @@ -223,16 +223,16 @@ using FastInterpolations end end - @testset "side=:right" begin - itp = constant_interp(x_vec, y_const; side=:right) + @testset "side=RightSide()" begin + itp = constant_interp(x_vec, y_const; side=RightSide()) for xq in [0.05, 0.5, 1.5, 2.5] cell = coeffs(itp, xq) @test cell(xq) ≈ itp(xq) atol=1e-12 end end - @testset "side=:nearest" begin - itp = constant_interp(x_vec, y_const; side=:nearest) + @testset "side=NearestSide()" begin + itp = constant_interp(x_vec, y_const; side=NearestSide()) for xq in [0.05, 0.5, 1.5, 2.5] cell = coeffs(itp, xq) @test cell(xq) ≈ itp(xq) atol=1e-12 diff --git a/test/test_complex_constant.jl b/test/test_complex_constant.jl index 6fa06d05..da62caa3 100644 --- a/test/test_complex_constant.jl +++ b/test/test_complex_constant.jl @@ -90,17 +90,17 @@ using FastInterpolations y = [1.0+1.0im, 2.0+2.0im, 3.0+3.0im, 4.0+4.0im] # :left mode - itp_left = constant_interp(x, y; side=:left) + itp_left = constant_interp(x, y; side=LeftSide()) @test itp_left(0.5) == 1.0+1.0im @test itp_left(1.5) == 2.0+2.0im # :right mode - itp_right = constant_interp(x, y; side=:right) + itp_right = constant_interp(x, y; side=RightSide()) @test itp_right(0.5) == 2.0+2.0im @test itp_right(1.0) == 2.0+2.0im # At grid point, still gets left value # :nearest mode - itp_nearest = constant_interp(x, y; side=:nearest) + itp_nearest = constant_interp(x, y; side=NearestSide()) @test itp_nearest(0.4) == 1.0+1.0im # Closer to left @test itp_nearest(0.6) == 2.0+2.0im # Closer to right @test itp_nearest(0.5) == 1.0+1.0im # Midpoint → left (tie-breaking) @@ -156,7 +156,7 @@ using FastInterpolations x = [0.0, 1.0, 2.0, 3.0] y = [1.0+1.0im, 2.0+2.0im, 3.0+3.0im, 4.0+4.0im] - itp = constant_interp(x, y; side=:left) + itp = constant_interp(x, y; side=LeftSide()) xq = [0.5, 1.5, 2.5] vals = itp(xq) @@ -175,7 +175,7 @@ using FastInterpolations x = [0.0, 1.0, 2.0, 3.0] y = [1.0+1.0im, 2.0+2.0im, 3.0+3.0im, 4.0+4.0im] - itp = constant_interp(x, y; side=:left) + itp = constant_interp(x, y; side=LeftSide()) xq = [0.5, 1.5, 2.5] vals = itp.(xq) @@ -224,7 +224,7 @@ using FastInterpolations x = [0.0, 1.0, 2.0, 3.0] y = [1.0+1.0im, 2.0+2.0im, 3.0+3.0im, 4.0+4.0im] - itp = constant_interp(x, y; side=:left) + itp = constant_interp(x, y; side=LeftSide()) xq = [0.5, 1.5, 2.5] output = Vector{ComplexF64}(undef, 3) diff --git a/test/test_complex_constant_series.jl b/test/test_complex_constant_series.jl index b23c7420..1a4aa08f 100644 --- a/test/test_complex_constant_series.jl +++ b/test/test_complex_constant_series.jl @@ -68,7 +68,7 @@ using FastInterpolations @test vals isa Vector{ComplexF64} # Constant interpolation returns step value (nearest by default) - # At 5.5 with :nearest, rounds to index 6 (x=5) → value 5+10i + # At 5.5 with NearestSide(), rounds to index 6 (x=5) → value 5+10i @test isapprox(vals[1], 5.0 + 10.0im, rtol=1e-10) end @@ -105,21 +105,21 @@ using FastInterpolations end # ======================================== - # Side Options (:nearest, :left, :right) + # Side Options (NearestSide, LeftSide, RightSide) # ======================================== @testset "Side options with Complex values" begin x = collect(0.0:1.0:5.0) # [0, 1, 2, 3, 4, 5] y1 = ComplexF64[1+1im, 2+2im, 3+3im, 4+4im, 5+5im, 6+6im] - # Test :left side - sitp_left = constant_interp(x, [y1]; side=:left) + # Test LeftSide() + sitp_left = constant_interp(x, [y1]; side=LeftSide()) @test sitp_left isa ConstantSeriesInterpolant{Float64, ComplexF64} vals_left = sitp_left(2.5) # Between 2 and 3 @test vals_left isa Vector{ComplexF64} @test isapprox(vals_left[1], 3.0 + 3.0im, rtol=1e-10) # Left value at x=2 - # Test :right side - sitp_right = constant_interp(x, [y1]; side=:right) + # Test RightSide() + sitp_right = constant_interp(x, [y1]; side=RightSide()) vals_right = sitp_right(2.5) @test isapprox(vals_right[1], 4.0 + 4.0im, rtol=1e-10) # Right value at x=3 end diff --git a/test/test_constant.jl b/test/test_constant.jl index f848a003..b1260889 100644 --- a/test/test_constant.jl +++ b/test/test_constant.jl @@ -1,88 +1,34 @@ # Tests for constant (step/piecewise constant) interpolation # Import internal items for testing -import FastInterpolations: @_dispatch_side, _constant_kernel, EvalValue, EvalDeriv1, EvalDeriv2 +import FastInterpolations: _constant_kernel, EvalValue, EvalDeriv1, EvalDeriv2 # ============================================================================ # Group 1: Infrastructure Tests (Phase 1) # ============================================================================ @testset "Constant Interpolation - Infrastructure" begin - @testset "@_dispatch_side macro" begin - @testset "dispatch :nearest" begin - result = @_dispatch_side :nearest => sv begin - sv - end - @test result === Val(:nearest) + @testset "AbstractSide type hierarchy" begin + @testset "AbstractSide is abstract" begin + @test isabstracttype(FastInterpolations.AbstractSide) end - @testset "dispatch :left" begin - result = @_dispatch_side :left => sv begin - sv - end - @test result === Val(:left) + @testset "concrete subtypes" begin + @test NearestSide() isa FastInterpolations.AbstractSide + @test LeftSide() isa FastInterpolations.AbstractSide + @test RightSide() isa FastInterpolations.AbstractSide end - @testset "dispatch :right" begin - result = @_dispatch_side :right => sv begin - sv - end - @test result === Val(:right) + @testset "singleton identity" begin + @test NearestSide() === NearestSide() + @test LeftSide() === LeftSide() + @test RightSide() === RightSide() end - @testset "dispatch from variable" begin - for sym in (:nearest, :left, :right) - result = @_dispatch_side sym => sv begin - sv - end - @test result === Val(sym) - end - end - - @testset "invalid symbol throws ArgumentError" begin - @test_throws ArgumentError begin - @_dispatch_side :invalid => sv begin - sv - end - end - end - - @testset "body execution with binding" begin - x = 10 - result = @_dispatch_side :left => sv begin - (sv, x * 2) - end - @test result === (Val(:left), 20) - end - - @testset "type stability" begin - function test_dispatch(side::Symbol) - @_dispatch_side side => sv begin - sv - end - end - @test test_dispatch(:nearest) === Val(:nearest) - @test test_dispatch(:left) === Val(:left) - @test test_dispatch(:right) === Val(:right) - end - end - - @testset "SideVal type" begin - @testset "SideVal definition" begin - @test isdefined(FastInterpolations, :SideVal) - @test FastInterpolations.SideVal isa Union - end - - @testset "SideVal members" begin - @test Val(:nearest) isa FastInterpolations.SideVal - @test Val(:left) isa FastInterpolations.SideVal - @test Val(:right) isa FastInterpolations.SideVal - end - - @testset "SideVal exclusions" begin - @test !(Val(:none) isa FastInterpolations.SideVal) - @test !(Val(:constant) isa FastInterpolations.SideVal) - @test !(Val(:other) isa FastInterpolations.SideVal) + @testset "type distinctness" begin + @test typeof(NearestSide()) !== typeof(LeftSide()) + @test typeof(LeftSide()) !== typeof(RightSide()) + @test typeof(NearestSide()) !== typeof(RightSide()) end end @@ -98,27 +44,27 @@ end y_right = 20.0 h = 1.0 - @testset "EvalValue - side=:left" begin + @testset "EvalValue - side=LeftSide()" begin op = EvalValue() - sv = Val(:left) + sv = LeftSide() @test _constant_kernel(op, y_left, y_right, h, 0.0, sv) == 10.0 @test _constant_kernel(op, y_left, y_right, h, 0.3, sv) == 10.0 @test _constant_kernel(op, y_left, y_right, h, 0.7, sv) == 10.0 @test _constant_kernel(op, y_left, y_right, h, 0.99, sv) == 10.0 end - @testset "EvalValue - side=:right" begin + @testset "EvalValue - side=RightSide()" begin op = EvalValue() - sv = Val(:right) + sv = RightSide() @test _constant_kernel(op, y_left, y_right, h, 0.0, sv) == 10.0 # grid point @test _constant_kernel(op, y_left, y_right, h, 0.3, sv) == 20.0 @test _constant_kernel(op, y_left, y_right, h, 0.7, sv) == 20.0 @test _constant_kernel(op, y_left, y_right, h, 0.99, sv) == 20.0 end - @testset "EvalValue - side=:nearest" begin + @testset "EvalValue - side=NearestSide()" begin op = EvalValue() - sv = Val(:nearest) + sv = NearestSide() @test _constant_kernel(op, y_left, y_right, h, 0.0, sv) == 10.0 @test _constant_kernel(op, y_left, y_right, h, 0.4, sv) == 10.0 @test _constant_kernel(op, y_left, y_right, h, 0.5, sv) == 10.0 # tie-breaking @@ -128,30 +74,30 @@ end @testset "EvalDeriv1 - all sides return zero" begin op = EvalDeriv1() - @test _constant_kernel(op, y_left, y_right, h, 0.5, Val(:nearest)) == 0.0 - @test _constant_kernel(op, y_left, y_right, h, 0.5, Val(:left)) == 0.0 - @test _constant_kernel(op, y_left, y_right, h, 0.5, Val(:right)) == 0.0 + @test _constant_kernel(op, y_left, y_right, h, 0.5, NearestSide()) == 0.0 + @test _constant_kernel(op, y_left, y_right, h, 0.5, LeftSide()) == 0.0 + @test _constant_kernel(op, y_left, y_right, h, 0.5, RightSide()) == 0.0 end @testset "EvalDeriv2 - all sides return zero" begin op = EvalDeriv2() - @test _constant_kernel(op, y_left, y_right, h, 0.5, Val(:nearest)) == 0.0 - @test _constant_kernel(op, y_left, y_right, h, 0.5, Val(:left)) == 0.0 - @test _constant_kernel(op, y_left, y_right, h, 0.5, Val(:right)) == 0.0 + @test _constant_kernel(op, y_left, y_right, h, 0.5, NearestSide()) == 0.0 + @test _constant_kernel(op, y_left, y_right, h, 0.5, LeftSide()) == 0.0 + @test _constant_kernel(op, y_left, y_right, h, 0.5, RightSide()) == 0.0 end @testset "Type preservation" begin - @test _constant_kernel(EvalValue(), 1.0, 2.0, 1.0, 0.5, Val(:left)) isa Float64 - @test _constant_kernel(EvalDeriv1(), 1.0, 2.0, 1.0, 0.5, Val(:left)) isa Float64 - @test _constant_kernel(EvalValue(), 1.0f0, 2.0f0, 1.0f0, 0.5f0, Val(:nearest)) isa Float32 - @test _constant_kernel(EvalDeriv1(), 1.0f0, 2.0f0, 1.0f0, 0.5f0, Val(:nearest)) isa Float32 + @test _constant_kernel(EvalValue(), 1.0, 2.0, 1.0, 0.5, LeftSide()) isa Float64 + @test _constant_kernel(EvalDeriv1(), 1.0, 2.0, 1.0, 0.5, LeftSide()) isa Float64 + @test _constant_kernel(EvalValue(), 1.0f0, 2.0f0, 1.0f0, 0.5f0, NearestSide()) isa Float32 + @test _constant_kernel(EvalDeriv1(), 1.0f0, 2.0f0, 1.0f0, 0.5f0, NearestSide()) isa Float32 end @testset "Non-uniform grid (h != 1)" begin op = EvalValue() - @test _constant_kernel(op, 100.0, 200.0, 2.0, 0.9, Val(:nearest)) == 100.0 - @test _constant_kernel(op, 100.0, 200.0, 2.0, 1.0, Val(:nearest)) == 100.0 # tie - @test _constant_kernel(op, 100.0, 200.0, 2.0, 1.1, Val(:nearest)) == 200.0 + @test _constant_kernel(op, 100.0, 200.0, 2.0, 0.9, NearestSide()) == 100.0 + @test _constant_kernel(op, 100.0, 200.0, 2.0, 1.0, NearestSide()) == 100.0 # tie + @test _constant_kernel(op, 100.0, 200.0, 2.0, 1.1, NearestSide()) == 200.0 end end @@ -169,15 +115,15 @@ end @test constant_interp(x, y, 1.5) == 20.0 @test constant_interp(x, y, 2.5) == 30.0 - @test constant_interp(x, y, 0.9; side=:left) == 10.0 - @test constant_interp(x, y, 1.1; side=:left) == 20.0 + @test constant_interp(x, y, 0.9; side=LeftSide()) == 10.0 + @test constant_interp(x, y, 1.1; side=LeftSide()) == 20.0 - @test constant_interp(x, y, 0.0; side=:right) == 10.0 - @test constant_interp(x, y, 0.1; side=:right) == 20.0 + @test constant_interp(x, y, 0.0; side=RightSide()) == 10.0 + @test constant_interp(x, y, 0.1; side=RightSide()) == 20.0 - @test constant_interp(x, y, 1.0; side=:nearest) == 20.0 - @test constant_interp(x, y, 1.0; side=:left) == 20.0 - @test constant_interp(x, y, 1.0; side=:right) == 20.0 + @test constant_interp(x, y, 1.0; side=NearestSide()) == 20.0 + @test constant_interp(x, y, 1.0; side=LeftSide()) == 20.0 + @test constant_interp(x, y, 1.0; side=RightSide()) == 20.0 @test constant_interp(x, y, 4.0) == 50.0 @@ -189,7 +135,7 @@ end result = constant_interp(x, y, [0.5, 1.5, 2.5]) @test result ≈ [10.0, 20.0, 30.0] - result_left = constant_interp(x, y, [0.5, 1.5]; side=:left) + result_left = constant_interp(x, y, [0.5, 1.5]; side=LeftSide()) @test result_left ≈ [10.0, 20.0] end @@ -262,8 +208,8 @@ end @test itp_int(1.5) == 20.0 # Test options pass through the wrappers correctly - @test constant_interp(x_int, y_int, 0; side=:right) == 10.0 - @test constant_interp(x_int, y_int, 0; side=:left) == 10.0 + @test constant_interp(x_int, y_int, 0; side=RightSide()) == 10.0 + @test constant_interp(x_int, y_int, 0; side=LeftSide()) == 10.0 @test constant_interp(x_int, y_int, -1; extrap=ConstExtrap()) == 10.0 @test constant_interp(x_int, y_int, 5; extrap=ConstExtrap()) == 50.0 @@ -273,7 +219,7 @@ end @test out3 ≈ [10.0, 50.0] # Test 2-arg callable with options - itp_left = constant_interp(x_int, y_int; side=:left) + itp_left = constant_interp(x_int, y_int; side=LeftSide()) @test itp_left(0.9) == 10.0 itp_wrap = constant_interp(x_int, y_int; extrap=WrapExtrap()) @test itp_wrap(4.0) == 10.0 @@ -318,10 +264,10 @@ end @test itp_const(-1.0) == 10.0 @test itp_const(5.0) == 50.0 - itp_left = constant_interp(x, y; side=:left) + itp_left = constant_interp(x, y; side=LeftSide()) @test itp_left(0.9) == 10.0 - itp_right = constant_interp(x, y; side=:right) + itp_right = constant_interp(x, y; side=RightSide()) @test itp_right(0.1) == 20.0 end @@ -329,16 +275,16 @@ end itp = constant_interp(x, y) @test @inferred(itp(0.5)) isa Float64 @test @inferred(constant_interp(x, y, 0.5)) isa Float64 - @test @inferred(constant_interp(x, y, 0.5; side=:left)) isa Float64 + @test @inferred(constant_interp(x, y, 0.5; side=LeftSide())) isa Float64 @test @inferred(constant_interp(x, y, 0.5; extrap=ConstExtrap())) isa Float64 @test @inferred(constant_interp(x, y, 0.5; deriv=1)) isa Float64 end - @testset "Invalid side argument throws ArgumentError" begin - @test_throws ArgumentError constant_interp(x, y, 0.5; side=:invalid) - @test_throws ArgumentError constant_interp(x, y, [0.5, 1.5]; side=:foo) + @testset "Invalid side argument throws TypeError" begin + @test_throws TypeError constant_interp(x, y, 0.5; side=:invalid) + @test_throws TypeError constant_interp(x, y, [0.5, 1.5]; side=:foo) out = zeros(2) - @test_throws ArgumentError constant_interp!(out, x, y, [0.5, 1.5]; side=:bar) + @test_throws TypeError constant_interp!(out, x, y, [0.5, 1.5]; side=:bar) end end @@ -362,8 +308,8 @@ end @testset "constant_interp! with side option" begin out = zeros(5) xq = [0.1, 0.3, 0.5, 0.7, 0.9] - constant_interp!(out, x, y, xq; side=:left) - allocs = @allocated constant_interp!(out, x, y, xq; side=:left) + constant_interp!(out, x, y, xq; side=LeftSide()) + allocs = @allocated constant_interp!(out, x, y, xq; side=LeftSide()) @test allocs <= ALLOC_THRESHOLD end @@ -409,8 +355,8 @@ end @test allocs <= ALLOC_THRESHOLD end - @testset "ConstantInterpolant with side=:left in-place" begin - itp = constant_interp(x, y; side=:left) + @testset "ConstantInterpolant with side=LeftSide() in-place" begin + itp = constant_interp(x, y; side=LeftSide()) out = zeros(5) xq = [0.1, 0.3, 0.5, 0.7, 0.9] itp(out, xq) @@ -513,15 +459,15 @@ end d2 = deriv2(itp_wrap) @test d2(4.5) == 0.0 - itp_left = constant_interp(x, y; side=:left) + itp_left = constant_interp(x, y; side=LeftSide()) @test deriv1(itp_left)(0.5) == 0.0 - itp_right = constant_interp(x, y; side=:right) + itp_right = constant_interp(x, y; side=RightSide()) @test deriv2(itp_right)(0.5) == 0.0 end @testset "Edge: Grid point behavior" begin - for side in (:nearest, :left, :right) + for side in (NearestSide(), LeftSide(), RightSide()) @test constant_interp(x, y, 0.0; side=side) == 10.0 @test constant_interp(x, y, 1.0; side=side) == 20.0 @test constant_interp(x, y, 2.0; side=side) == 30.0 @@ -538,11 +484,11 @@ end end @testset "Edge: Wrap boundary cases" begin - @test constant_interp(x, y, 4.0; extrap=WrapExtrap(), side=:left) == 10.0 - @test constant_interp(x, y, 4.0; extrap=WrapExtrap(), side=:right) == 10.0 - @test constant_interp(x, y, 4.0; extrap=WrapExtrap(), side=:nearest) == 10.0 - @test constant_interp(x, y, 4.5; extrap=WrapExtrap(), side=:nearest) == 10.0 - @test constant_interp(x, y, 5.0; extrap=WrapExtrap(), side=:nearest) == 20.0 + @test constant_interp(x, y, 4.0; extrap=WrapExtrap(), side=LeftSide()) == 10.0 + @test constant_interp(x, y, 4.0; extrap=WrapExtrap(), side=RightSide()) == 10.0 + @test constant_interp(x, y, 4.0; extrap=WrapExtrap(), side=NearestSide()) == 10.0 + @test constant_interp(x, y, 4.5; extrap=WrapExtrap(), side=NearestSide()) == 10.0 + @test constant_interp(x, y, 5.0; extrap=WrapExtrap(), side=NearestSide()) == 20.0 end @testset "Edge: Extrapolation + deriv" begin @@ -554,17 +500,17 @@ end end @testset "Edge: Midpoint tie-breaking" begin - @test constant_interp(x, y, 0.5; side=:nearest) == 10.0 - @test constant_interp(x, y, 1.5; side=:nearest) == 20.0 - @test constant_interp(x, y, 2.5; side=:nearest) == 30.0 + @test constant_interp(x, y, 0.5; side=NearestSide()) == 10.0 + @test constant_interp(x, y, 1.5; side=NearestSide()) == 20.0 + @test constant_interp(x, y, 2.5; side=NearestSide()) == 30.0 end @testset "Edge: Non-uniform grid" begin x_nu = [0.0, 0.5, 2.0, 2.5, 4.0] y_nu = [10.0, 20.0, 30.0, 40.0, 50.0] - @test constant_interp(x_nu, y_nu, 0.25; side=:nearest) == 10.0 - @test constant_interp(x_nu, y_nu, 1.0; side=:nearest) == 20.0 - @test constant_interp(x_nu, y_nu, 1.5; side=:nearest) == 30.0 + @test constant_interp(x_nu, y_nu, 0.25; side=NearestSide()) == 10.0 + @test constant_interp(x_nu, y_nu, 1.0; side=NearestSide()) == 20.0 + @test constant_interp(x_nu, y_nu, 1.5; side=NearestSide()) == 30.0 end @testset "Edge: Float32 type preservation" begin @@ -583,9 +529,9 @@ end @testset "Edge: Minimum grid size (2 points)" begin x_min = [0.0, 1.0] y_min = [10.0, 20.0] - @test constant_interp(x_min, y_min, 0.5; side=:nearest) == 10.0 - @test constant_interp(x_min, y_min, 0.5; side=:left) == 10.0 - @test constant_interp(x_min, y_min, 0.5; side=:right) == 20.0 + @test constant_interp(x_min, y_min, 0.5; side=NearestSide()) == 10.0 + @test constant_interp(x_min, y_min, 0.5; side=LeftSide()) == 10.0 + @test constant_interp(x_min, y_min, 0.5; side=RightSide()) == 20.0 @test constant_interp(x_min, y_min, 0.0) == 10.0 @test constant_interp(x_min, y_min, 1.0) == 20.0 end @@ -593,8 +539,8 @@ end @testset "Edge: Range input (O(1) path)" begin x_range = 0.0:0.5:4.0 y_range = collect(10.0:10.0:90.0) - @test constant_interp(x_range, y_range, 0.25; side=:nearest) == 10.0 - @test constant_interp(x_range, y_range, 0.75; side=:nearest) == 20.0 + @test constant_interp(x_range, y_range, 0.25; side=NearestSide()) == 10.0 + @test constant_interp(x_range, y_range, 0.75; side=NearestSide()) == 20.0 itp_range = constant_interp(x_range, y_range) @test itp_range(0.25) == 10.0 diff --git a/test/test_constant_anchor.jl b/test/test_constant_anchor.jl index 0a3a840d..1b613577 100644 --- a/test/test_constant_anchor.jl +++ b/test/test_constant_anchor.jl @@ -74,12 +74,12 @@ using FastInterpolations end # ======================================== - # Evaluation Tests - :nearest mode + # Evaluation Tests - NearestSide() mode # ======================================== - @testset "itp(aq) for :nearest mode" begin + @testset "itp(aq) for NearestSide() mode" begin x = collect(range(0.0, 1.0, 11)) y = collect(1.0:11.0) - itp = constant_interp(x, y; side=:nearest, extrap=ExtendExtrap()) + itp = constant_interp(x, y; side=NearestSide(), extrap=ExtendExtrap()) xq_points = [0.05, 0.15, 0.35, 0.65, 0.95] @@ -90,12 +90,12 @@ using FastInterpolations end # ======================================== - # Evaluation Tests - :left mode + # Evaluation Tests - LeftSide() mode # ======================================== - @testset "itp(aq) for :left mode" begin + @testset "itp(aq) for LeftSide() mode" begin x = collect(range(0.0, 1.0, 11)) y = collect(1.0:11.0) - itp = constant_interp(x, y; side=:left, extrap=ExtendExtrap()) + itp = constant_interp(x, y; side=LeftSide(), extrap=ExtendExtrap()) xq_points = [0.05, 0.15, 0.35, 0.65, 0.95] @@ -106,12 +106,12 @@ using FastInterpolations end # ======================================== - # Evaluation Tests - :right mode + # Evaluation Tests - RightSide() mode # ======================================== - @testset "itp(aq) for :right mode" begin + @testset "itp(aq) for RightSide() mode" begin x = collect(range(0.0, 1.0, 11)) y = collect(1.0:11.0) - itp = constant_interp(x, y; side=:right, extrap=ExtendExtrap()) + itp = constant_interp(x, y; side=RightSide(), extrap=ExtendExtrap()) xq_points = [0.05, 0.15, 0.35, 0.65, 0.95] diff --git a/test/test_constant_series_interp.jl b/test/test_constant_series_interp.jl index c302bbcd..a2254c34 100644 --- a/test/test_constant_series_interp.jl +++ b/test/test_constant_series_interp.jl @@ -47,9 +47,9 @@ const FI = FastInterpolations end @testset "side parameter preserved" begin - for side_mode in (:nearest, :left, :right) + for side_mode in (NearestSide(), LeftSide(), RightSide()) sitp = constant_interp(x, ys; side=side_mode) - @test sitp.side === Val(side_mode) + @test sitp.side === side_mode end end end @@ -78,7 +78,7 @@ const FI = FastInterpolations x = [0.0, 1.0, 2.0, 3.0] y1 = [1.0, 2.0, 3.0, 4.0] # step function y2 = [5.0, 5.0, 5.0, 5.0] # constant - sitp = constant_interp(x, [y1, y2]; side=:left) + sitp = constant_interp(x, [y1, y2]; side=LeftSide()) @testset "at grid points" begin for i in 1:(length(x)-1) @@ -88,14 +88,14 @@ const FI = FastInterpolations end end - @testset "midpoints with side=:left" begin + @testset "midpoints with side=LeftSide()" begin result = sitp(0.5) @test result[1] == 1.0 # Left value in [0,1) @test result[2] == 5.0 # Constant end - @testset "midpoints with side=:nearest" begin - sitp_nearest = constant_interp(x, [y1, y2]; side=:nearest) + @testset "midpoints with side=NearestSide()" begin + sitp_nearest = constant_interp(x, [y1, y2]; side=NearestSide()) # At 0.4, closer to 0 -> y[1]=1.0 result = sitp_nearest(0.4) @test result[1] == 1.0 @@ -276,7 +276,7 @@ const FI = FastInterpolations ys = [y1, y2] # Both should give same results - sitp = constant_interp(x, ys; side=:nearest) + sitp = constant_interp(x, ys; side=NearestSide()) # Compare at multiple points xq = [0.15, 0.35, 0.55, 0.75, 0.95] diff --git a/test/test_cumulative_integrate.jl b/test/test_cumulative_integrate.jl index 9155d33a..e12c939f 100644 --- a/test/test_cumulative_integrate.jl +++ b/test/test_cumulative_integrate.jl @@ -28,7 +28,7 @@ using FastInterpolations @testset "constant" begin y_const = collect(1.0:length(x)) - for side in (:left, :right, :nearest) + for side in (LeftSide(), RightSide(), NearestSide()) itp_k = constant_interp(x, y_const; side=side, extrap=NoExtrap()) cum = cumulative_integrate(itp_k) @test cum[1] == 0.0 @@ -80,7 +80,7 @@ using FastInterpolations @testset "constant series" begin y_c1 = collect(1.0:length(x)) y_c2 = collect(length(x):-1.0:1.0) - for side in (:left, :right, :nearest) + for side in (LeftSide(), RightSide(), NearestSide()) sitp = constant_interp(x, [y_c1, y_c2]; side=side) itp1 = constant_interp(x, y_c1; side=side) itp2 = constant_interp(x, y_c2; side=side) diff --git a/test/test_derivatives.jl b/test/test_derivatives.jl index 93b0e924..aee6f5cf 100644 --- a/test/test_derivatives.jl +++ b/test/test_derivatives.jl @@ -1824,7 +1824,7 @@ end # DerivativeView Wrapper # Constant kernel @test FastInterpolations._constant_kernel( - FastInterpolations.EvalDeriv3(), 5.0, 5.0, 0.5, 0.2, Val(:left) + FastInterpolations.EvalDeriv3(), 5.0, 5.0, 0.5, 0.2, LeftSide() ) === zero(Float64) end diff --git a/test/test_integral_1d.jl b/test/test_integral_1d.jl index fe510211..b3922208 100644 --- a/test/test_integral_1d.jl +++ b/test/test_integral_1d.jl @@ -22,7 +22,7 @@ using FastInterpolations @testset "constant finite by side mode" begin y = collect(1.0:length(x)) - for side in (:left, :right, :nearest) + for side in (LeftSide(), RightSide(), NearestSide()) itp = constant_interp(x, y; side=side, extrap=NoExtrap()) I = integrate(itp, 0.2, 0.7) @test isfinite(I) diff --git a/test/test_integral_allocation.jl b/test/test_integral_allocation.jl index 0dfaa603..40bc0ce7 100644 --- a/test/test_integral_allocation.jl +++ b/test/test_integral_allocation.jl @@ -48,7 +48,7 @@ using FastInterpolations @testset "constant 1D: zero allocation" begin x = collect(range(0.0, 1.0, length=21)) y = collect(1.0:length(x)) - for side in (:left, :right, :nearest) + for side in (LeftSide(), RightSide(), NearestSide()) itp = constant_interp(x, y; side=side, extrap=NoExtrap()) a, b = 0.2, 0.7 diff --git a/test/test_integral_extrap.jl b/test/test_integral_extrap.jl index 31debb34..1323578e 100644 --- a/test/test_integral_extrap.jl +++ b/test/test_integral_extrap.jl @@ -32,29 +32,29 @@ using FastInterpolations @test integrate(itp, 1.0, 1.3) ≈ y[end] * 0.3 atol=1e-12 end - @testset ":constant tails (constant interp, side=:left)" begin + @testset ":constant tails (constant interp, side=LeftSide())" begin # y[1]=1.0, y[2]≈1.001, y[end-1]≈1.999, y[end]=2.0 - # With side=:left, extrap=ConstExtrap() must use y[1] (left) and y[end] (right) + # With side=LeftSide(), extrap=ConstExtrap() must use y[1] (left) and y[end] (right) y = @. x^2 + 1.0 - itp = constant_interp(x, y; side=:left, extrap=ConstExtrap()) + itp = constant_interp(x, y; side=LeftSide(), extrap=ConstExtrap()) # pure left tail @test integrate(itp, -0.4, 0.0) ≈ y[1] * 0.4 atol=1e-12 # pure right tail @test integrate(itp, 1.0, 1.6) ≈ y[end] * 0.6 atol=1e-12 # mixed: left tail + in-domain - in_part = integrate(constant_interp(x, y; side=:left, extrap=NoExtrap()), 0.0, 0.5) + in_part = integrate(constant_interp(x, y; side=LeftSide(), extrap=NoExtrap()), 0.0, 0.5) @test integrate(itp, -0.3, 0.5) ≈ y[1] * 0.3 + in_part atol=1e-12 end - @testset ":constant tails (constant interp, side=:right)" begin + @testset ":constant tails (constant interp, side=RightSide())" begin y = @. x^2 + 1.0 - itp = constant_interp(x, y; side=:right, extrap=ConstExtrap()) + itp = constant_interp(x, y; side=RightSide(), extrap=ConstExtrap()) # pure left tail — must use y[1], NOT y[2] @test integrate(itp, -0.4, 0.0) ≈ y[1] * 0.4 atol=1e-12 # pure right tail — must use y[end], NOT y[end-1] @test integrate(itp, 1.0, 1.6) ≈ y[end] * 0.6 atol=1e-12 # mixed: right tail + in-domain - in_part = integrate(constant_interp(x, y; side=:right, extrap=NoExtrap()), 0.5, 1.0) + in_part = integrate(constant_interp(x, y; side=RightSide(), extrap=NoExtrap()), 0.5, 1.0) @test integrate(itp, 0.5, 1.3) ≈ in_part + y[end] * 0.3 atol=1e-12 end diff --git a/test/test_integral_fulldomain.jl b/test/test_integral_fulldomain.jl index 410c6889..119758d5 100644 --- a/test/test_integral_fulldomain.jl +++ b/test/test_integral_fulldomain.jl @@ -18,7 +18,7 @@ using FastInterpolations @test integrate(itp_l) ≈ integrate(itp_l, first(x), last(x)) atol=1e-14 @test integrate(itp_q) ≈ integrate(itp_q, first(x), last(x)) atol=1e-14 - for side in (:left, :right, :nearest) + for side in (LeftSide(), RightSide(), NearestSide()) itp_k = constant_interp(x, y_const; side=side, extrap=NoExtrap()) @test integrate(itp_k) ≈ integrate(itp_k, first(x), last(x)) atol=1e-14 end @@ -59,7 +59,7 @@ using FastInterpolations @testset "constant series" begin y_c = hcat(collect(1.0:length(x)), collect(length(x):-1.0:1.0)) - for side in (:left, :right, :nearest) + for side in (LeftSide(), RightSide(), NearestSide()) sitp = constant_interp(x, [y_c[:, 1], y_c[:, 2]]; side=side) itp1 = constant_interp(x, y_c[:, 1]; side=side) itp2 = constant_interp(x, y_c[:, 2]; side=side) @@ -97,7 +97,7 @@ using FastInterpolations end @testset "constant ND" begin - for side in ((:left, :left), (:right, :right), (:nearest, :nearest)) + for side in ((LeftSide(), LeftSide()), (RightSide(), RightSide()), (NearestSide(), NearestSide())) itp = constant_interp((xg, yg), data_2d; side=side, extrap=(NoExtrap(), NoExtrap())) lo = (first(xg), first(yg)) hi = (last(xg), last(yg)) diff --git a/test/test_integral_nd.jl b/test/test_integral_nd.jl index 330188af..771b9081 100644 --- a/test/test_integral_nd.jl +++ b/test/test_integral_nd.jl @@ -28,7 +28,7 @@ using FastInterpolations @testset "constant nd finite by side mode" begin data = [sin(xi) + cos(yj) for xi in x, yj in y] - for side in ((:left, :left), (:right, :right), (:nearest, :nearest)) + for side in ((LeftSide(), LeftSide()), (RightSide(), RightSide()), (NearestSide(), NearestSide())) itp = constant_interp((x, y), data; side=side, extrap=(NoExtrap(), NoExtrap())) @test isfinite(integrate(itp, (0.2, 0.3), (0.8, 1.7))) end diff --git a/test/test_integral_nd_exactness.jl b/test/test_integral_nd_exactness.jl index eaa3c3b2..531a7d14 100644 --- a/test/test_integral_nd_exactness.jl +++ b/test/test_integral_nd_exactness.jl @@ -265,7 +265,7 @@ using FastInterpolations y = collect(range(0.0, 3.0, length=9)) z = collect(range(0.0, 1.0, length=7)) data = fill(7.0, length(x), length(y), length(z)) - for side in (:left, :right, :nearest) + for side in (LeftSide(), RightSide(), NearestSide()) itp = constant_interp((x, y, z), data; side=side, extrap=NoExtrap()) lo, hi = (0.3, 0.5, 0.1), (1.7, 2.4, 0.8) expected = 7.0 * (hi[1]-lo[1]) * (hi[2]-lo[2]) * (hi[3]-lo[3]) @@ -278,7 +278,7 @@ using FastInterpolations y = collect(range(0.0, 3.0, length=9)) z = collect(range(0.0, 1.0, length=7)) data = [sin(xi) + cos(yj) + zk for xi in x, yj in y, zk in z] - for side in (:left, :right, :nearest) + for side in (LeftSide(), RightSide(), NearestSide()) itp = constant_interp((x, y, z), data; side=side, extrap=NoExtrap()) @test isfinite(integrate(itp, (0.3, 0.5, 0.1), (1.7, 2.4, 0.8))) end diff --git a/test/test_integral_series.jl b/test/test_integral_series.jl index 391109f7..f5edba40 100644 --- a/test/test_integral_series.jl +++ b/test/test_integral_series.jl @@ -87,7 +87,7 @@ using FastInterpolations end @testset "ConstantSeriesInterpolant" begin - for side in (:left, :right, :nearest) + for side in (LeftSide(), RightSide(), NearestSide()) @testset "side=$side" begin sitp = constant_interp(x, [y1, y2]; side=side) itp1 = constant_interp(x, y1; side=side) diff --git a/test/test_nd_constant.jl b/test/test_nd_constant.jl index 3226e3e3..0df6c82b 100644 --- a/test/test_nd_constant.jl +++ b/test/test_nd_constant.jl @@ -31,8 +31,8 @@ end 21.0 22.0 23.0 24.0; 31.0 32.0 33.0 34.0] - @testset "side=:left (default)" begin - itp = constant_interp((x, y), data; side=:left) + @testset "side=LeftSide() (default)" begin + itp = constant_interp((x, y), data; side=LeftSide()) # At grid points, returns left value of the interval containing the point # At (0,0): interval idx=(1,1), data[1,1]=11 @@ -40,30 +40,30 @@ end # At (1,2): x in interval 2, y in interval 3, data[2,3]=23 @test itp((1.0, 2.0)) == 23.0 # Interior point # At boundary (2,3): last intervals, data[2,3]=23 (not 34!) - # Because side=:left always returns left corner of interval + # Because side=LeftSide() always returns left corner of interval @test itp((2.0, 3.0)) == 23.0 # Far corner (in last interval) - # Between grid points with :left, always select left neighbor + # Between grid points with LeftSide(), always select left neighbor @test itp((0.5, 0.5)) == 11.0 # Cell (1,1) - left corner @test itp((1.5, 2.5)) == 23.0 # Cell (2,3) - left corner @test itp((0.9, 0.9)) == 11.0 # Still left corner end - @testset "side=:right" begin - itp = constant_interp((x, y), data; side=:right) + @testset "side=RightSide()" begin + itp = constant_interp((x, y), data; side=RightSide()) # At grid points, still returns left value (dL == 0) @test itp((0.0, 0.0)) == 11.0 @test itp((1.0, 2.0)) == 23.0 - # Between grid points with :right, select right neighbor + # Between grid points with RightSide(), select right neighbor @test itp((0.5, 0.5)) == 22.0 # Right corner of cell (1,1) → data[2,2] @test itp((0.1, 0.1)) == 22.0 # Any offset → right @test itp((1.5, 2.5)) == 34.0 # Cell (2,3) → data[3,4] end - @testset "side=:nearest" begin - itp = constant_interp((x, y), data; side=:nearest) + @testset "side=NearestSide()" begin + itp = constant_interp((x, y), data; side=NearestSide()) # At grid points @test itp((0.0, 0.0)) == 11.0 @@ -79,8 +79,8 @@ end end @testset "per-axis side configuration" begin - # side=(:left, :right) → left on x-axis, right on y-axis - itp = constant_interp((x, y), data; side=(:left, :right)) + # side=(LeftSide(), RightSide()) → left on x-axis, right on y-axis + itp = constant_interp((x, y), data; side=(LeftSide(), RightSide())) @test itp((0.5, 0.5)) == 12.0 # x: left (idx=1), y: right (idx=2) → data[1,2] @test itp((1.5, 0.5)) == 22.0 # x: left (idx=2), y: right (idx=2) → data[2,2] @@ -100,14 +100,14 @@ end data[i, j, k] = 100.0 * i + 10.0 * j + k end - itp = constant_interp((x, y, z), data; side=:left) + itp = constant_interp((x, y, z), data; side=LeftSide()) # At origin: all indices = 1 @test itp((0.0, 0.0, 0.0)) == 111.0 # data[1,1,1] - # At boundary (1,1,1): still in interval [0,1], side=:left → data[1,1,1] + # At boundary (1,1,1): still in interval [0,1], side=LeftSide() → data[1,1,1] @test itp((1.0, 1.0, 1.0)) == 111.0 # data[1,1,1] (left of each interval) - # Interior points with :left + # Interior points with LeftSide() @test itp((0.5, 0.5, 0.5)) == 111.0 # All left → data[1,1,1] @test itp((0.9, 0.9, 0.9)) == 111.0 # Still all left end @@ -147,7 +147,7 @@ end y = [0.0, 1.0, 2.0] data = [1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0] - itp = constant_interp((x, y), data; side=:left) + itp = constant_interp((x, y), data; side=LeftSide()) # Query with vector instead of tuple @test itp([0.5, 0.5]) == itp((0.5, 0.5)) @@ -162,7 +162,7 @@ end y = [0.0, 1.0, 2.0] data = [1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0] - itp = constant_interp((x, y), data; side=:left) + itp = constant_interp((x, y), data; side=LeftSide()) # Batch query: tuple of vectors xs = [0.5, 1.5, 0.5] @@ -183,7 +183,7 @@ end y = [0.0, 1.0, 2.0] data = [1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0] - itp = constant_interp((x, y), data; side=:left) + itp = constant_interp((x, y), data; side=LeftSide()) # Batch query: vector of tuples points = [(0.5, 0.5), (1.5, 0.5), (0.5, 1.5)] @@ -217,7 +217,7 @@ end end @testset "extrap=ConstExtrap()" begin - itp = constant_interp((x, y), data; extrap=ConstExtrap(), side=:left) + itp = constant_interp((x, y), data; extrap=ConstExtrap(), side=LeftSide()) # In domain @test itp((0.5, 0.5)) == 1.0 @@ -226,15 +226,15 @@ end @test itp((-0.5, 0.5)) == 1.0 # Clamp x to 0 → data[1,1] = 1 @test itp((0.5, -0.5)) == 1.0 # Clamp y to 0 → data[1,1] = 1 # Clamp x to 2.0, which is at last interval boundary - # With side=:left, we get data[2, 1] = 4.0 + # With side=LeftSide(), we get data[2, 1] = 4.0 @test itp((2.5, 0.5)) == 4.0 # Clamp x to 2 → interval 2, data[2,1] # Clamp y to 2.0, which is at last interval boundary - # With side=:left, we get data[1, 2] = 2.0 + # With side=LeftSide(), we get data[1, 2] = 2.0 @test itp((0.5, 2.5)) == 2.0 # Clamp y to 2 → interval 2, data[1,2] end @testset "extrap=WrapExtrap()" begin - itp = constant_interp((x, y), data; extrap=WrapExtrap(), side=:left) + itp = constant_interp((x, y), data; extrap=WrapExtrap(), side=LeftSide()) # In domain @test itp((0.5, 0.5)) == 1.0 @@ -246,7 +246,7 @@ end @testset "per-axis extrap configuration" begin # extrap=(NoExtrap(), ConstExtrap()) → strict on x, clamp on y - itp = constant_interp((x, y), data; extrap=(NoExtrap(), ConstExtrap()), side=:left) + itp = constant_interp((x, y), data; extrap=(NoExtrap(), ConstExtrap()), side=LeftSide()) # y clamped to 2.0 → data[1, 2] = 2.0 @test itp((0.5, 2.5)) == 2.0 # y clamped to last interval @@ -262,7 +262,7 @@ end y = range(0.0, 2.0, 3) # [0, 1, 2] data = [1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0] - itp = constant_interp((x, y), data; side=:left) + itp = constant_interp((x, y), data; side=LeftSide()) @test itp((0.5, 0.5)) == 1.0 @test itp((1.5, 0.5)) == 4.0 @@ -276,7 +276,7 @@ end y = [0.0, 0.5, 2.0] # Non-uniform Vector data = [1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0] - itp = constant_interp((x, y), data; side=:left) + itp = constant_interp((x, y), data; side=LeftSide()) @test itp((0.5, 0.3)) == 1.0 # First cell @test itp((1.5, 1.0)) == 5.0 # y in [0.5, 2.0) → idx 2 @@ -312,13 +312,13 @@ end data = [1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0] # Scalar one-shot - result = constant_interp((x, y), data, (0.5, 0.5); side=:left) + result = constant_interp((x, y), data, (0.5, 0.5); side=LeftSide()) @test result == 1.0 # Batch one-shot (SoA) xs = [0.5, 1.5] ys = [0.5, 0.5] - results = constant_interp((x, y), data, (xs, ys); side=:left) + results = constant_interp((x, y), data, (xs, ys); side=LeftSide()) @test results[1] == 1.0 @test results[2] == 4.0 end @@ -346,7 +346,7 @@ end y = [0.0, 1.0, 2.0] data = ComplexF64[1+1im 2+2im 3+3im; 4+4im 5+5im 6+6im; 7+7im 8+8im 9+9im] - itp = constant_interp((x, y), data; side=:left) + itp = constant_interp((x, y), data; side=LeftSide()) @test itp((0.5, 0.5)) == 1.0 + 1.0im @test itp((1.5, 1.5)) == 5.0 + 5.0im @@ -409,9 +409,9 @@ end y = range(0.0, 1.0, 6) data = [xi + yj for xi in x, yj in y] query = (1.0, 0.5) - constant_interp((x, y), data, query; side=:left) - constant_interp((x, y), data, query; side=:left) - @allocated constant_interp((x, y), data, query; side=:left) + constant_interp((x, y), data, query; side=LeftSide()) + constant_interp((x, y), data, query; side=LeftSide()) + @allocated constant_interp((x, y), data, query; side=LeftSide()) end function _alloc_test_constant_right() @@ -419,9 +419,9 @@ end y = range(0.0, 1.0, 6) data = [xi + yj for xi in x, yj in y] query = (1.0, 0.5) - constant_interp((x, y), data, query; side=:right) - constant_interp((x, y), data, query; side=:right) - @allocated constant_interp((x, y), data, query; side=:right) + constant_interp((x, y), data, query; side=RightSide()) + constant_interp((x, y), data, query; side=RightSide()) + @allocated constant_interp((x, y), data, query; side=RightSide()) end function _alloc_test_constant_extrap_constant() @@ -470,11 +470,11 @@ end @test _alloc_test_constant_default() <= ND_ALLOC_THRESHOLD end - @testset "zero-alloc scalar (Range grids, side=:left)" begin + @testset "zero-alloc scalar (Range grids, side=LeftSide())" begin @test _alloc_test_constant_left() <= ND_ALLOC_THRESHOLD end - @testset "zero-alloc scalar (Range grids, side=:right)" begin + @testset "zero-alloc scalar (Range grids, side=RightSide())" begin @test _alloc_test_constant_right() <= ND_ALLOC_THRESHOLD end @@ -554,9 +554,9 @@ end y = collect(range(0.0, 1.0, 15)) data = [xi + yj for xi in x, yj in y] query = (1.0, 0.5) - constant_interp((x, y), data, query; side=:left) - constant_interp((x, y), data, query; side=:left) - @allocated constant_interp((x, y), data, query; side=:left) + constant_interp((x, y), data, query; side=LeftSide()) + constant_interp((x, y), data, query; side=LeftSide()) + @allocated constant_interp((x, y), data, query; side=LeftSide()) end function _alloc_test_constant_vector_3d() @@ -575,7 +575,7 @@ end @test _alloc_test_constant_vector_default() <= ND_ALLOC_THRESHOLD end - @testset "zero-alloc scalar (Vector grids, side=:left)" begin + @testset "zero-alloc scalar (Vector grids, side=LeftSide())" begin @test _alloc_test_constant_vector_left() <= ND_ALLOC_THRESHOLD end diff --git a/test/test_nd_coverage.jl b/test/test_nd_coverage.jl index 3345fff7..5a781dcf 100644 --- a/test/test_nd_coverage.jl +++ b/test/test_nd_coverage.jl @@ -892,15 +892,14 @@ end end end -@testset "core/utils.jl — @_dispatch_side_nd mixed-sides fallback" begin - # Mixed side=(:left, :right) triggers the _to_side_vals fallback in @_dispatch_side_nd - # (the `else` branch when _is_uniform_side returns false) +@testset "mixed-sides fallback" begin + # Mixed side=(LeftSide(), RightSide()) grids = (collect(range(0.0, 2.0, 4)), collect(range(0.0, 3.0, 5))) data = [Float64(10i + j) for i in 1:4, j in 1:5] - result_left = constant_interp(grids, data, (0.5, 0.5); side=:left) - result_right = constant_interp(grids, data, (0.5, 0.5); side=:right) - result_mixed = constant_interp(grids, data, (0.5, 0.5); side=(:left, :right)) + result_left = constant_interp(grids, data, (0.5, 0.5); side=LeftSide()) + result_right = constant_interp(grids, data, (0.5, 0.5); side=RightSide()) + result_mixed = constant_interp(grids, data, (0.5, 0.5); side=(LeftSide(), RightSide())) @test result_mixed isa Float64 @test result_mixed != result_left # x picks left, y picks right corner diff --git a/test/test_nd_utils_shared.jl b/test/test_nd_utils_shared.jl index 856a918f..82c0070d 100644 --- a/test/test_nd_utils_shared.jl +++ b/test/test_nd_utils_shared.jl @@ -85,31 +85,31 @@ import FastInterpolations: _resolve_extrap_nd, _resolve_search_nd, _resolve_bcs_ # _resolve_side_nd (NEW - for ConstantInterpolantND) # ======================================== @testset "_resolve_side_nd" begin - @testset "broadcast single symbol to N-tuple" begin - result = _resolve_side_nd(:nearest, Val(3)) - @test result === (:nearest, :nearest, :nearest) + @testset "broadcast single AbstractSide to N-tuple" begin + result = _resolve_side_nd(NearestSide(), Val(3)) + @test result === (NearestSide(), NearestSide(), NearestSide()) - result = _resolve_side_nd(:left, Val(2)) - @test result === (:left, :left) + result = _resolve_side_nd(LeftSide(), Val(2)) + @test result === (LeftSide(), LeftSide()) - result = _resolve_side_nd(:right, Val(4)) - @test result === (:right, :right, :right, :right) + result = _resolve_side_nd(RightSide(), Val(4)) + @test result === (RightSide(), RightSide(), RightSide(), RightSide()) end @testset "passthrough matching tuple" begin - result = _resolve_side_nd((:nearest, :left, :right), Val(3)) - @test result === (:nearest, :left, :right) + result = _resolve_side_nd((NearestSide(), LeftSide(), RightSide()), Val(3)) + @test result === (NearestSide(), LeftSide(), RightSide()) end @testset "reject wrong-length tuple" begin - @test_throws ArgumentError _resolve_side_nd((:nearest, :left), Val(3)) - @test_throws ArgumentError _resolve_side_nd((:nearest, :left, :right, :nearest), Val(3)) + @test_throws ArgumentError _resolve_side_nd((NearestSide(), LeftSide()), Val(3)) + @test_throws ArgumentError _resolve_side_nd((NearestSide(), LeftSide(), RightSide(), NearestSide()), Val(3)) end - @testset "reject invalid symbol" begin - @test_throws ArgumentError _resolve_side_nd(:invalid, Val(2)) - @test_throws ArgumentError _resolve_side_nd((:nearest, :invalid), Val(2)) - @test_throws ArgumentError _resolve_side_nd(:center, Val(2)) + @testset "reject non-AbstractSide via MethodError" begin + @test_throws MethodError _resolve_side_nd(:invalid, Val(2)) + @test_throws MethodError _resolve_side_nd(:nearest, Val(2)) + @test_throws MethodError _resolve_side_nd(:center, Val(2)) end end diff --git a/test/test_polyfit_bc.jl b/test/test_polyfit_bc.jl index c19a2983..1dd06892 100644 --- a/test/test_polyfit_bc.jl +++ b/test/test_polyfit_bc.jl @@ -291,7 +291,7 @@ end # First 4 points: x = 0.0, 0.1, 0.2, 0.3 fvals = (f(0.0), f(0.1), f(0.2), f(0.3)) - d_left = FastInterpolations._compute_deriv1(PolyFit{3}(), Val(:left), fvals, inv_h) + d_left = FastInterpolations._compute_deriv1(PolyFit{3}(), LeftSide(), fvals, inv_h) @test d_left ≈ 0.0 atol=1e-10 # f'(0) = 0 end @@ -305,7 +305,7 @@ end # Last 4 points: x = 0.7, 0.8, 0.9, 1.0 fvals = (f(0.7), f(0.8), f(0.9), f(1.0)) - d_right = FastInterpolations._compute_deriv1(PolyFit{3}(), Val(:right), fvals, inv_h) + d_right = FastInterpolations._compute_deriv1(PolyFit{3}(), RightSide(), fvals, inv_h) @test d_right ≈ 3.0 atol=1e-10 # f'(1) = 3 end @@ -317,7 +317,7 @@ end f = x -> x^2 - 2x + 1 fvals = (f(0.0), f(0.25), f(0.5), f(0.75)) - d_left = FastInterpolations._compute_deriv1(PolyFit{3}(), Val(:left), fvals, inv_h) + d_left = FastInterpolations._compute_deriv1(PolyFit{3}(), LeftSide(), fvals, inv_h) @test d_left ≈ -2.0 atol=1e-10 end @@ -328,7 +328,7 @@ end f = x -> x^2 - 2x + 1 fvals = (f(1.25), f(1.5), f(1.75), f(2.0)) - d_right = FastInterpolations._compute_deriv1(PolyFit{3}(), Val(:right), fvals, inv_h) + d_right = FastInterpolations._compute_deriv1(PolyFit{3}(), RightSide(), fvals, inv_h) @test d_right ≈ 2.0 atol=1e-10 end @@ -339,11 +339,11 @@ end f = x -> 3x + 2 fvals_left = (f(0.0), f(0.5), f(1.0), f(1.5)) - d_left = FastInterpolations._compute_deriv1(PolyFit{3}(), Val(:left), fvals_left, inv_h) + d_left = FastInterpolations._compute_deriv1(PolyFit{3}(), LeftSide(), fvals_left, inv_h) @test d_left ≈ 3.0 atol=1e-14 # Exact for linear fvals_right = (f(1.5), f(2.0), f(2.5), f(3.0)) - d_right = FastInterpolations._compute_deriv1(PolyFit{3}(), Val(:right), fvals_right, inv_h) + d_right = FastInterpolations._compute_deriv1(PolyFit{3}(), RightSide(), fvals_right, inv_h) @test d_right ≈ 3.0 atol=1e-14 end @@ -353,7 +353,7 @@ end f = x -> x^2 fvals = NTuple{4, Float32}(Float32.(f.([0.0, 0.1, 0.2, 0.3]))) - d_left = FastInterpolations._compute_deriv1(PolyFit{3}(), Val(:left), fvals, inv_h) + d_left = FastInterpolations._compute_deriv1(PolyFit{3}(), LeftSide(), fvals, inv_h) @test d_left isa Float32 @test d_left ≈ Float32(0.0) atol=1e-5 # f'(0) = 0 end @@ -369,7 +369,7 @@ end @testset "Coefficient Precomputation: Left Endpoint" begin # Non-uniform grid: [0, 0.1, 0.3, 0.6] xvals = (0.0, 0.1, 0.3, 0.6) - coeffs = FastInterpolations._compute_deriv1_coeffs(PolyFit{3}(), Val(:left), xvals) + coeffs = FastInterpolations._compute_deriv1_coeffs(PolyFit{3}(), LeftSide(), xvals) # Test: coefficients should be finite @test all(isfinite, coeffs) @@ -384,7 +384,7 @@ end @testset "Coefficient Precomputation: Right Endpoint" begin # Non-uniform grid: [0.4, 0.7, 0.9, 1.0] xvals = (0.4, 0.7, 0.9, 1.0) - coeffs = FastInterpolations._compute_deriv1_coeffs(PolyFit{3}(), Val(:right), xvals) + coeffs = FastInterpolations._compute_deriv1_coeffs(PolyFit{3}(), RightSide(), xvals) # Test: coefficients should be finite @test all(isfinite, coeffs) @@ -403,13 +403,13 @@ end f_lin(x) = 2x + 3 # Left endpoint - c_left = FastInterpolations._compute_deriv1_coeffs(PolyFit{3}(), Val(:left), x_left) + c_left = FastInterpolations._compute_deriv1_coeffs(PolyFit{3}(), LeftSide(), x_left) f_left = NTuple{4}(f_lin.(collect(x_left))) d_left = FastInterpolations._weighted_sum(c_left, f_left) @test d_left ≈ 2.0 atol=1e-13 # Right endpoint - c_right = FastInterpolations._compute_deriv1_coeffs(PolyFit{3}(), Val(:right), x_right) + c_right = FastInterpolations._compute_deriv1_coeffs(PolyFit{3}(), RightSide(), x_right) f_right = NTuple{4}(f_lin.(collect(x_right))) d_right = FastInterpolations._weighted_sum(c_right, f_right) @test d_right ≈ 2.0 atol=1e-13 @@ -423,13 +423,13 @@ end f_quad(x) = x^2 - 3x + 2 # Left endpoint: f'(0) = -3 - c_left = FastInterpolations._compute_deriv1_coeffs(PolyFit{3}(), Val(:left), x_left) + c_left = FastInterpolations._compute_deriv1_coeffs(PolyFit{3}(), LeftSide(), x_left) f_left = NTuple{4}(f_quad.(collect(x_left))) d_left = FastInterpolations._weighted_sum(c_left, f_left) @test d_left ≈ -3.0 atol=1e-12 # Right endpoint: f'(1) = -1 - c_right = FastInterpolations._compute_deriv1_coeffs(PolyFit{3}(), Val(:right), x_right) + c_right = FastInterpolations._compute_deriv1_coeffs(PolyFit{3}(), RightSide(), x_right) f_right = NTuple{4}(f_quad.(collect(x_right))) d_right = FastInterpolations._weighted_sum(c_right, f_right) @test d_right ≈ -1.0 atol=1e-12 @@ -443,13 +443,13 @@ end f_cub(x) = x^3 # Left endpoint: f'(0) = 0 - c_left = FastInterpolations._compute_deriv1_coeffs(PolyFit{3}(), Val(:left), x_left) + c_left = FastInterpolations._compute_deriv1_coeffs(PolyFit{3}(), LeftSide(), x_left) f_left = NTuple{4}(f_cub.(collect(x_left))) d_left = FastInterpolations._weighted_sum(c_left, f_left) @test d_left ≈ 0.0 atol=1e-12 # Right endpoint: f'(1) = 3 - c_right = FastInterpolations._compute_deriv1_coeffs(PolyFit{3}(), Val(:right), x_right) + c_right = FastInterpolations._compute_deriv1_coeffs(PolyFit{3}(), RightSide(), x_right) f_right = NTuple{4}(f_cub.(collect(x_right))) d_right = FastInterpolations._weighted_sum(c_right, f_right) @test d_right ≈ 3.0 atol=1e-11 @@ -463,22 +463,22 @@ end # Left endpoint using precomputed coefficients x_left = (xs[1], xs[2], xs[3], xs[4]) - c_left = FastInterpolations._compute_deriv1_coeffs(PolyFit{3}(), Val(:left), x_left) + c_left = FastInterpolations._compute_deriv1_coeffs(PolyFit{3}(), LeftSide(), x_left) f_left = (ys[1], ys[2], ys[3], ys[4]) d_left_precomp = FastInterpolations._weighted_sum(c_left, f_left) # Left endpoint using unified API (4-arg with PolyFit{3}) - d_left_onfly = FastInterpolations._estimate_endpoint_derivative(xs, ys, Val(:left), PolyFit{3}()) + d_left_onfly = FastInterpolations._estimate_endpoint_derivative(xs, ys, LeftSide(), PolyFit{3}()) @test d_left_precomp ≈ d_left_onfly rtol=1e-14 # Right endpoint n = length(xs) x_right = (xs[n-3], xs[n-2], xs[n-1], xs[n]) - c_right = FastInterpolations._compute_deriv1_coeffs(PolyFit{3}(), Val(:right), x_right) + c_right = FastInterpolations._compute_deriv1_coeffs(PolyFit{3}(), RightSide(), x_right) f_right = (ys[n-3], ys[n-2], ys[n-1], ys[n]) d_right_precomp = FastInterpolations._weighted_sum(c_right, f_right) - d_right_onfly = FastInterpolations._estimate_endpoint_derivative(xs, ys, Val(:right), PolyFit{3}()) + d_right_onfly = FastInterpolations._estimate_endpoint_derivative(xs, ys, RightSide(), PolyFit{3}()) @test d_right_precomp ≈ d_right_onfly rtol=1e-14 end @@ -492,20 +492,20 @@ end # Precomputed (non-uniform kernel with uniform data) x_left = (xs_uniform[1], xs_uniform[2], xs_uniform[3], xs_uniform[4]) - c_left = FastInterpolations._compute_deriv1_coeffs(PolyFit{3}(), Val(:left), x_left) + c_left = FastInterpolations._compute_deriv1_coeffs(PolyFit{3}(), LeftSide(), x_left) f_left = (ys[1], ys[2], ys[3], ys[4]) d_precomp = FastInterpolations._weighted_sum(c_left, f_left) # Direct uniform kernel inv_h = 1 / h - d_direct = FastInterpolations._compute_deriv1(PolyFit{3}(), Val(:left), f_left, inv_h) + d_direct = FastInterpolations._compute_deriv1(PolyFit{3}(), LeftSide(), f_left, inv_h) @test d_precomp ≈ d_direct rtol=1e-13 end @testset "Float32 Precision" begin xvals = NTuple{4, Float32}((0.0f0, 0.1f0, 0.3f0, 0.6f0)) - coeffs = FastInterpolations._compute_deriv1_coeffs(PolyFit{3}(), Val(:left), xvals) + coeffs = FastInterpolations._compute_deriv1_coeffs(PolyFit{3}(), LeftSide(), xvals) @test coeffs[1] isa Float32 @test coeffs[2] isa Float32 @@ -786,12 +786,12 @@ end # Left endpoint: f'(0) = -2 fvals_left = (y[1], y[2], y[3]) - d_left = FastInterpolations._compute_deriv1(PolyFit{2}(), Val(:left), fvals_left, inv_h) + d_left = FastInterpolations._compute_deriv1(PolyFit{2}(), LeftSide(), fvals_left, inv_h) @test d_left ≈ df_quad(x[1]) rtol=1e-12 # Right endpoint: f'(2) = 2 fvals_right = (y[end-2], y[end-1], y[end]) - d_right = FastInterpolations._compute_deriv1(PolyFit{2}(), Val(:right), fvals_right, inv_h) + d_right = FastInterpolations._compute_deriv1(PolyFit{2}(), RightSide(), fvals_right, inv_h) @test d_right ≈ df_quad(x[end]) rtol=1e-12 end @@ -806,8 +806,8 @@ end fvals_left = (y[1], y[2], y[3]) fvals_right = (y[end-2], y[end-1], y[end]) - d_left = FastInterpolations._compute_deriv1(PolyFit{2}(), Val(:left), fvals_left, inv_h) - d_right = FastInterpolations._compute_deriv1(PolyFit{2}(), Val(:right), fvals_right, inv_h) + d_left = FastInterpolations._compute_deriv1(PolyFit{2}(), LeftSide(), fvals_left, inv_h) + d_right = FastInterpolations._compute_deriv1(PolyFit{2}(), RightSide(), fvals_right, inv_h) @test d_left ≈ df_linear rtol=1e-12 @test d_right ≈ df_linear rtol=1e-12 @@ -824,14 +824,14 @@ end # Left endpoint via coefficient kernels x_left = (x[1], x[2], x[3]) f_left = (y[1], y[2], y[3]) - coeffs_left = FastInterpolations._compute_deriv1_coeffs(PolyFit{2}(), Val(:left), x_left) + coeffs_left = FastInterpolations._compute_deriv1_coeffs(PolyFit{2}(), LeftSide(), x_left) d_left = FastInterpolations._weighted_sum(coeffs_left, f_left) @test d_left ≈ df_quad(x[1]) rtol=1e-12 # Right endpoint via coefficient kernels x_right = (x[end-2], x[end-1], x[end]) f_right = (y[end-2], y[end-1], y[end]) - coeffs_right = FastInterpolations._compute_deriv1_coeffs(PolyFit{2}(), Val(:right), x_right) + coeffs_right = FastInterpolations._compute_deriv1_coeffs(PolyFit{2}(), RightSide(), x_right) d_right = FastInterpolations._weighted_sum(coeffs_right, f_right) @test d_right ≈ df_quad(x[end]) rtol=1e-12 end @@ -844,8 +844,8 @@ end x_range = range(0.0, 2.0, 11) y_range = collect(f_quad.(x_range)) - d_left_range = FastInterpolations._estimate_endpoint_derivative(x_range, y_range, Val(:left), PolyFit{2}()) - d_right_range = FastInterpolations._estimate_endpoint_derivative(x_range, y_range, Val(:right), PolyFit{2}()) + d_left_range = FastInterpolations._estimate_endpoint_derivative(x_range, y_range, LeftSide(), PolyFit{2}()) + d_right_range = FastInterpolations._estimate_endpoint_derivative(x_range, y_range, RightSide(), PolyFit{2}()) @test d_left_range ≈ df_quad(x_range[1]) rtol=1e-12 @test d_right_range ≈ df_quad(x_range[end]) rtol=1e-12 @@ -854,8 +854,8 @@ end x_vec = [0.0, 0.3, 0.8, 1.5, 2.0] y_vec = f_quad.(x_vec) - d_left_vec = FastInterpolations._estimate_endpoint_derivative(x_vec, y_vec, Val(:left), PolyFit{2}()) - d_right_vec = FastInterpolations._estimate_endpoint_derivative(x_vec, y_vec, Val(:right), PolyFit{2}()) + d_left_vec = FastInterpolations._estimate_endpoint_derivative(x_vec, y_vec, LeftSide(), PolyFit{2}()) + d_right_vec = FastInterpolations._estimate_endpoint_derivative(x_vec, y_vec, RightSide(), PolyFit{2}()) @test d_left_vec ≈ df_quad(x_vec[1]) rtol=1e-12 @test d_right_vec ≈ df_quad(x_vec[end]) rtol=1e-12 @@ -869,8 +869,8 @@ end x_range = range(0.0, 2.0, 21) y_range = collect(f_cubic.(x_range)) - d_left = FastInterpolations._estimate_endpoint_derivative(x_range, y_range, Val(:left), PolyFit{3}()) - d_right = FastInterpolations._estimate_endpoint_derivative(x_range, y_range, Val(:right), PolyFit{3}()) + d_left = FastInterpolations._estimate_endpoint_derivative(x_range, y_range, LeftSide(), PolyFit{3}()) + d_right = FastInterpolations._estimate_endpoint_derivative(x_range, y_range, RightSide(), PolyFit{3}()) @test d_left ≈ df_cubic(x_range[1]) rtol=1e-10 @test d_right ≈ df_cubic(x_range[end]) rtol=1e-10 @@ -884,8 +884,8 @@ end x = range(0.0, 1.0, 11) y = collect(f_cubic.(x)) - d2_left = FastInterpolations._estimate_endpoint_derivative(x, y, Val(:left), PolyFit{2}()) - d3_left = FastInterpolations._estimate_endpoint_derivative(x, y, Val(:left), PolyFit{3}()) + d2_left = FastInterpolations._estimate_endpoint_derivative(x, y, LeftSide(), PolyFit{2}()) + d3_left = FastInterpolations._estimate_endpoint_derivative(x, y, LeftSide(), PolyFit{3}()) exact = df_cubic(x[1]) # = 0 @@ -903,12 +903,12 @@ end inv_h = inv(step(x)) fvals = (y[1], y[2], y[3]) - d_left = FastInterpolations._compute_deriv1(PolyFit{2}(), Val(:left), fvals, inv_h) + d_left = FastInterpolations._compute_deriv1(PolyFit{2}(), LeftSide(), fvals, inv_h) @test d_left isa Float32 @test isfinite(d_left) # Unified API - d_api = FastInterpolations._estimate_endpoint_derivative(x, y, Val(:left), PolyFit{2}()) + d_api = FastInterpolations._estimate_endpoint_derivative(x, y, LeftSide(), PolyFit{2}()) @test d_api isa Float32 @test isfinite(d_api) end @@ -933,12 +933,12 @@ end # Left endpoint: f'(0) = 3 fvals_left = (y[1], y[2]) - d_left = FastInterpolations._compute_deriv1(PolyFit{1}(), Val(:left), fvals_left, inv_h) + d_left = FastInterpolations._compute_deriv1(PolyFit{1}(), LeftSide(), fvals_left, inv_h) @test d_left ≈ df_linear rtol=1e-14 # Right endpoint: f'(2) = 3 fvals_right = (y[end-1], y[end]) - d_right = FastInterpolations._compute_deriv1(PolyFit{1}(), Val(:right), fvals_right, inv_h) + d_right = FastInterpolations._compute_deriv1(PolyFit{1}(), RightSide(), fvals_right, inv_h) @test d_right ≈ df_linear rtol=1e-14 end @@ -955,13 +955,13 @@ end # Left endpoint: exact f'(0) = 0, but LinearFit will have O(h) error fvals_left = (y[1], y[2]) - d_left = FastInterpolations._compute_deriv1(PolyFit{1}(), Val(:left), fvals_left, inv_h) + d_left = FastInterpolations._compute_deriv1(PolyFit{1}(), LeftSide(), fvals_left, inv_h) # Forward difference: (f(h) - f(0)) / h = h ≈ 0.1 @test abs(d_left - df_quad(x[1])) ≈ h atol=1e-12 # Right endpoint: exact f'(1) = 2 fvals_right = (y[end-1], y[end]) - d_right = FastInterpolations._compute_deriv1(PolyFit{1}(), Val(:right), fvals_right, inv_h) + d_right = FastInterpolations._compute_deriv1(PolyFit{1}(), RightSide(), fvals_right, inv_h) # Backward difference: (f(1) - f(1-h)) / h = (1 - (1-h)²) / h = 2 - h @test d_right ≈ df_quad(x[end]) - h atol=1e-12 end @@ -978,14 +978,14 @@ end # Left endpoint via coefficient kernels x_left = (x[1], x[2]) f_left = (y[1], y[2]) - coeffs_left = FastInterpolations._compute_deriv1_coeffs(PolyFit{1}(), Val(:left), x_left) + coeffs_left = FastInterpolations._compute_deriv1_coeffs(PolyFit{1}(), LeftSide(), x_left) d_left = FastInterpolations._weighted_sum(coeffs_left, f_left) @test d_left ≈ df_linear rtol=1e-13 # Right endpoint via coefficient kernels x_right = (x[end-1], x[end]) f_right = (y[end-1], y[end]) - coeffs_right = FastInterpolations._compute_deriv1_coeffs(PolyFit{1}(), Val(:right), x_right) + coeffs_right = FastInterpolations._compute_deriv1_coeffs(PolyFit{1}(), RightSide(), x_right) d_right = FastInterpolations._weighted_sum(coeffs_right, f_right) @test d_right ≈ df_linear rtol=1e-13 end @@ -997,13 +997,13 @@ end x_left = (x[1], x[2]) f_left = (y[1], y[2]) - coeffs_left = FastInterpolations._compute_deriv1_coeffs(PolyFit{1}(), Val(:left), x_left) + coeffs_left = FastInterpolations._compute_deriv1_coeffs(PolyFit{1}(), LeftSide(), x_left) d_left = FastInterpolations._weighted_sum(coeffs_left, f_left) @test d_left ≈ 0.0 atol=1e-14 x_right = (x[end-1], x[end]) f_right = (y[end-1], y[end]) - coeffs_right = FastInterpolations._compute_deriv1_coeffs(PolyFit{1}(), Val(:right), x_right) + coeffs_right = FastInterpolations._compute_deriv1_coeffs(PolyFit{1}(), RightSide(), x_right) d_right = FastInterpolations._weighted_sum(coeffs_right, f_right) @test d_right ≈ 0.0 atol=1e-14 end @@ -1016,8 +1016,8 @@ end x_range = range(0.0, 2.0, 11) y_range = collect(f_linear.(x_range)) - d_left_range = FastInterpolations._estimate_endpoint_derivative(x_range, y_range, Val(:left), PolyFit{1}()) - d_right_range = FastInterpolations._estimate_endpoint_derivative(x_range, y_range, Val(:right), PolyFit{1}()) + d_left_range = FastInterpolations._estimate_endpoint_derivative(x_range, y_range, LeftSide(), PolyFit{1}()) + d_right_range = FastInterpolations._estimate_endpoint_derivative(x_range, y_range, RightSide(), PolyFit{1}()) @test d_left_range ≈ df_linear rtol=1e-14 @test d_right_range ≈ df_linear rtol=1e-14 @@ -1026,8 +1026,8 @@ end x_vec = [0.0, 0.3, 0.8, 1.5, 2.0] y_vec = f_linear.(x_vec) - d_left_vec = FastInterpolations._estimate_endpoint_derivative(x_vec, y_vec, Val(:left), PolyFit{1}()) - d_right_vec = FastInterpolations._estimate_endpoint_derivative(x_vec, y_vec, Val(:right), PolyFit{1}()) + d_left_vec = FastInterpolations._estimate_endpoint_derivative(x_vec, y_vec, LeftSide(), PolyFit{1}()) + d_right_vec = FastInterpolations._estimate_endpoint_derivative(x_vec, y_vec, RightSide(), PolyFit{1}()) @test d_left_vec ≈ df_linear rtol=1e-13 @test d_right_vec ≈ df_linear rtol=1e-13 @@ -1042,9 +1042,9 @@ end y = collect(f_quad.(x)) h = step(x) - d1_left = FastInterpolations._estimate_endpoint_derivative(x, y, Val(:left), PolyFit{1}()) - d2_left = FastInterpolations._estimate_endpoint_derivative(x, y, Val(:left), PolyFit{2}()) - d3_left = FastInterpolations._estimate_endpoint_derivative(x, y, Val(:left), PolyFit{3}()) + d1_left = FastInterpolations._estimate_endpoint_derivative(x, y, LeftSide(), PolyFit{1}()) + d2_left = FastInterpolations._estimate_endpoint_derivative(x, y, LeftSide(), PolyFit{2}()) + d3_left = FastInterpolations._estimate_endpoint_derivative(x, y, LeftSide(), PolyFit{3}()) exact = df_quad(x[1]) # = 0 @@ -1062,12 +1062,12 @@ end inv_h = inv(step(x)) fvals = (y[1], y[2]) - d_left = FastInterpolations._compute_deriv1(PolyFit{1}(), Val(:left), fvals, inv_h) + d_left = FastInterpolations._compute_deriv1(PolyFit{1}(), LeftSide(), fvals, inv_h) @test d_left isa Float32 @test d_left ≈ 3.0f0 rtol=1e-6 # Unified API - d_api = FastInterpolations._estimate_endpoint_derivative(x, y, Val(:left), PolyFit{1}()) + d_api = FastInterpolations._estimate_endpoint_derivative(x, y, LeftSide(), PolyFit{1}()) @test d_api isa Float32 @test d_api ≈ 3.0f0 rtol=1e-6 end @@ -1299,27 +1299,27 @@ end @testset "PolyFit{D} → Deriv1 Conversion" begin # LinearFit - bc_p1 = FastInterpolations.materialize_bc(LinearFit(), x, y, Val(:left)) + bc_p1 = FastInterpolations.materialize_bc(LinearFit(), x, y, LeftSide()) @test bc_p1 isa Deriv1{Float64} @test isfinite(bc_p1.val) # QuadraticFit - bc_p2 = FastInterpolations.materialize_bc(QuadraticFit(), x, y, Val(:left)) + bc_p2 = FastInterpolations.materialize_bc(QuadraticFit(), x, y, LeftSide()) @test bc_p2 isa Deriv1{Float64} @test isfinite(bc_p2.val) # CubicFit - bc_p3 = FastInterpolations.materialize_bc(CubicFit(), x, y, Val(:left)) + bc_p3 = FastInterpolations.materialize_bc(CubicFit(), x, y, LeftSide()) @test bc_p3 isa Deriv1{Float64} @test isfinite(bc_p3.val) # Right endpoint - bc_right = FastInterpolations.materialize_bc(QuadraticFit(), x, y, Val(:right)) + bc_right = FastInterpolations.materialize_bc(QuadraticFit(), x, y, RightSide()) @test bc_right isa Deriv1{Float64} @test isfinite(bc_right.val) # LinearFit right endpoint - bc_linear_right = FastInterpolations.materialize_bc(LinearFit(), x, y, Val(:right)) + bc_linear_right = FastInterpolations.materialize_bc(LinearFit(), x, y, RightSide()) @test bc_linear_right isa Deriv1{Float64} @test isfinite(bc_linear_right.val) end @@ -1327,19 +1327,19 @@ end @testset "Passthrough for Concrete BCs" begin # Deriv1 should pass through unchanged bc_d1 = Deriv1(1.5) - result_d1 = FastInterpolations.materialize_bc(bc_d1, x, y, Val(:left)) + result_d1 = FastInterpolations.materialize_bc(bc_d1, x, y, LeftSide()) @test result_d1 === bc_d1 # Deriv2 should pass through unchanged bc_d2 = Deriv2(0.0) - result_d2 = FastInterpolations.materialize_bc(bc_d2, x, y, Val(:left)) + result_d2 = FastInterpolations.materialize_bc(bc_d2, x, y, LeftSide()) @test result_d2 === bc_d2 end @testset "Estimated Values Match Direct API" begin # materialize_bc should produce same values as _estimate_endpoint_derivative - direct_left = FastInterpolations._estimate_endpoint_derivative(x, y, Val(:left), PolyFit{2}()) - materialized = FastInterpolations.materialize_bc(QuadraticFit(), x, y, Val(:left)) + direct_left = FastInterpolations._estimate_endpoint_derivative(x, y, LeftSide(), PolyFit{2}()) + materialized = FastInterpolations.materialize_bc(QuadraticFit(), x, y, LeftSide()) @test materialized.val ≈ direct_left rtol=1e-14 end @@ -1405,7 +1405,7 @@ end # Compute coefficients using in-place API (x as NTuple) c_left = Vector{Float64}(undef, 5) β_left = Vector{Float64}(undef, 5) - FastInterpolations._compute_deriv1_coeffs!(c_left, β_left, PolyFit{4}(), Val(:left), x) + FastInterpolations._compute_deriv1_coeffs!(c_left, β_left, PolyFit{4}(), LeftSide(), x) for i in 1:5 @test c_left[i] ≈ expected_left[i] rtol=1e-12 @@ -1416,7 +1416,7 @@ end expected_right = (3.0, -16.0, 36.0, -48.0, 25.0) ./ 12.0 .* inv_h c_right = Vector{Float64}(undef, 5) β_right = Vector{Float64}(undef, 5) - FastInterpolations._compute_deriv1_coeffs!(c_right, β_right, PolyFit{4}(), Val(:right), x) + FastInterpolations._compute_deriv1_coeffs!(c_right, β_right, PolyFit{4}(), RightSide(), x) for i in 1:5 @test c_right[i] ≈ expected_right[i] rtol=1e-12 @@ -1433,7 +1433,7 @@ end c_left = Vector{Float64}(undef, 6) β_left = Vector{Float64}(undef, 6) - FastInterpolations._compute_deriv1_coeffs!(c_left, β_left, PolyFit{5}(), Val(:left), x) + FastInterpolations._compute_deriv1_coeffs!(c_left, β_left, PolyFit{5}(), LeftSide(), x) for i in 1:6 @test c_left[i] ≈ expected_left[i] rtol=1e-12 @@ -1453,7 +1453,7 @@ end # Left endpoint: f'(0) = 0 deriv_left = FastInterpolations._estimate_endpoint_derivative( - xs, ys, Val(:left), PolyFit{D}() + xs, ys, LeftSide(), PolyFit{D}() ) @test deriv_left ≈ 0.0 atol=1e-10 @@ -1461,7 +1461,7 @@ end x_end = xs[end] expected_right = D * x_end^(D - 1) deriv_right = FastInterpolations._estimate_endpoint_derivative( - xs, ys, Val(:right), PolyFit{D}() + xs, ys, RightSide(), PolyFit{D}() ) @test deriv_right ≈ expected_right rtol=1e-10 end @@ -1476,7 +1476,7 @@ end # f'(0) = 0 for D > 1 deriv_left = FastInterpolations._estimate_endpoint_derivative( - xs, ys, Val(:left), PolyFit{D}() + xs, ys, LeftSide(), PolyFit{D}() ) @test deriv_left ≈ 0.0 atol=1e-10 @@ -1484,7 +1484,7 @@ end x_end = xs[end] expected_right = D * x_end^(D - 1) deriv_right = FastInterpolations._estimate_endpoint_derivative( - xs, ys, Val(:right), PolyFit{D}() + xs, ys, RightSide(), PolyFit{D}() ) @test deriv_right ≈ expected_right rtol=1e-9 end @@ -1520,21 +1520,21 @@ end inv_h = inv(0.1) expected_d1_left = (ys[2] - ys[1]) * inv_h result_d1_left = FastInterpolations._estimate_endpoint_derivative( - xs, ys, Val(:left), PolyFit{1}() + xs, ys, LeftSide(), PolyFit{1}() ) @test result_d1_left ≈ expected_d1_left rtol=1e-14 # D=2: (-3f₁ + 4f₂ - f₃) / (2h) expected_d2_left = (-3 * ys[1] + 4 * ys[2] - ys[3]) / (2 * 0.1) result_d2_left = FastInterpolations._estimate_endpoint_derivative( - xs, ys, Val(:left), PolyFit{2}() + xs, ys, LeftSide(), PolyFit{2}() ) @test result_d2_left ≈ expected_d2_left rtol=1e-14 # D=3: (-11f₁ + 18f₂ - 9f₃ + 2f₄) / (6h) expected_d3_left = (-11 * ys[1] + 18 * ys[2] - 9 * ys[3] + 2 * ys[4]) / (6 * 0.1) result_d3_left = FastInterpolations._estimate_endpoint_derivative( - xs, ys, Val(:left), PolyFit{3}() + xs, ys, LeftSide(), PolyFit{3}() ) @test result_d3_left ≈ expected_d3_left rtol=1e-14 end @@ -1554,16 +1554,16 @@ end D = length(x_tuple) - 1 N = D + 1 - for side in (:left, :right) + for side in (LeftSide(), RightSide()) # Specialized path (NTuple-returning, for D=1,2,3) c_specialized = FastInterpolations._compute_deriv1_coeffs( - PolyFit{D}(), Val(side), x_tuple + PolyFit{D}(), side, x_tuple ) # Generic barycentric path (in-place, x as NTuple) c_barycentric = Vector{Float64}(undef, N) β_barycentric = Vector{Float64}(undef, N) - k = (side === :left) ? 1 : N + k = (side isa LeftSide) ? 1 : N FastInterpolations._d1_coeffs_at_node!(c_barycentric, β_barycentric, x_tuple, k, Val(N)) # Compare element-wise (NTuple vs Vector) @@ -1631,59 +1631,59 @@ end @testset "_compute_deriv1 (Uniform Grid)" begin # Warmup calls (trigger compilation) - FastInterpolations._compute_deriv1(PolyFit{1}(), Val(:left), f2, inv_h) - FastInterpolations._compute_deriv1(PolyFit{2}(), Val(:left), f3, inv_h) - FastInterpolations._compute_deriv1(PolyFit{3}(), Val(:left), f4, inv_h) + FastInterpolations._compute_deriv1(PolyFit{1}(), LeftSide(), f2, inv_h) + FastInterpolations._compute_deriv1(PolyFit{2}(), LeftSide(), f3, inv_h) + FastInterpolations._compute_deriv1(PolyFit{3}(), LeftSide(), f4, inv_h) # D=1 (LinearFit) - alloc_d1_left = @allocated FastInterpolations._compute_deriv1(PolyFit{1}(), Val(:left), f2, inv_h) - alloc_d1_right = @allocated FastInterpolations._compute_deriv1(PolyFit{1}(), Val(:right), f2, inv_h) + alloc_d1_left = @allocated FastInterpolations._compute_deriv1(PolyFit{1}(), LeftSide(), f2, inv_h) + alloc_d1_right = @allocated FastInterpolations._compute_deriv1(PolyFit{1}(), RightSide(), f2, inv_h) @test alloc_d1_left <= ALLOC_THRESHOLD @test alloc_d1_right <= ALLOC_THRESHOLD # D=2 (QuadraticFit) - alloc_d2_left = @allocated FastInterpolations._compute_deriv1(PolyFit{2}(), Val(:left), f3, inv_h) - alloc_d2_right = @allocated FastInterpolations._compute_deriv1(PolyFit{2}(), Val(:right), f3, inv_h) + alloc_d2_left = @allocated FastInterpolations._compute_deriv1(PolyFit{2}(), LeftSide(), f3, inv_h) + alloc_d2_right = @allocated FastInterpolations._compute_deriv1(PolyFit{2}(), RightSide(), f3, inv_h) @test alloc_d2_left <= ALLOC_THRESHOLD @test alloc_d2_right <= ALLOC_THRESHOLD # D=3 (CubicFit) - alloc_d3_left = @allocated FastInterpolations._compute_deriv1(PolyFit{3}(), Val(:left), f4, inv_h) - alloc_d3_right = @allocated FastInterpolations._compute_deriv1(PolyFit{3}(), Val(:right), f4, inv_h) + alloc_d3_left = @allocated FastInterpolations._compute_deriv1(PolyFit{3}(), LeftSide(), f4, inv_h) + alloc_d3_right = @allocated FastInterpolations._compute_deriv1(PolyFit{3}(), RightSide(), f4, inv_h) @test alloc_d3_left <= ALLOC_THRESHOLD @test alloc_d3_right <= ALLOC_THRESHOLD end @testset "_compute_deriv1_coeffs (Non-Uniform Grid)" begin # Warmup calls - FastInterpolations._compute_deriv1_coeffs(PolyFit{1}(), Val(:left), x2) - FastInterpolations._compute_deriv1_coeffs(PolyFit{2}(), Val(:left), x3) - FastInterpolations._compute_deriv1_coeffs(PolyFit{3}(), Val(:left), x4) + FastInterpolations._compute_deriv1_coeffs(PolyFit{1}(), LeftSide(), x2) + FastInterpolations._compute_deriv1_coeffs(PolyFit{2}(), LeftSide(), x3) + FastInterpolations._compute_deriv1_coeffs(PolyFit{3}(), LeftSide(), x4) # D=1 (LinearFit) - alloc_d1_left = @allocated FastInterpolations._compute_deriv1_coeffs(PolyFit{1}(), Val(:left), x2) - alloc_d1_right = @allocated FastInterpolations._compute_deriv1_coeffs(PolyFit{1}(), Val(:right), x2) + alloc_d1_left = @allocated FastInterpolations._compute_deriv1_coeffs(PolyFit{1}(), LeftSide(), x2) + alloc_d1_right = @allocated FastInterpolations._compute_deriv1_coeffs(PolyFit{1}(), RightSide(), x2) @test alloc_d1_left <= ALLOC_THRESHOLD @test alloc_d1_right <= ALLOC_THRESHOLD # D=2 (QuadraticFit) - alloc_d2_left = @allocated FastInterpolations._compute_deriv1_coeffs(PolyFit{2}(), Val(:left), x3) - alloc_d2_right = @allocated FastInterpolations._compute_deriv1_coeffs(PolyFit{2}(), Val(:right), x3) + alloc_d2_left = @allocated FastInterpolations._compute_deriv1_coeffs(PolyFit{2}(), LeftSide(), x3) + alloc_d2_right = @allocated FastInterpolations._compute_deriv1_coeffs(PolyFit{2}(), RightSide(), x3) @test alloc_d2_left <= ALLOC_THRESHOLD @test alloc_d2_right <= ALLOC_THRESHOLD # D=3 (CubicFit) - alloc_d3_left = @allocated FastInterpolations._compute_deriv1_coeffs(PolyFit{3}(), Val(:left), x4) - alloc_d3_right = @allocated FastInterpolations._compute_deriv1_coeffs(PolyFit{3}(), Val(:right), x4) + alloc_d3_left = @allocated FastInterpolations._compute_deriv1_coeffs(PolyFit{3}(), LeftSide(), x4) + alloc_d3_right = @allocated FastInterpolations._compute_deriv1_coeffs(PolyFit{3}(), RightSide(), x4) @test alloc_d3_left <= ALLOC_THRESHOLD @test alloc_d3_right <= ALLOC_THRESHOLD end @testset "_weighted_sum (NTuple)" begin # Precompute coefficients for weighted_sum tests - c2 = FastInterpolations._compute_deriv1_coeffs(PolyFit{1}(), Val(:left), x2) - c3 = FastInterpolations._compute_deriv1_coeffs(PolyFit{2}(), Val(:left), x3) - c4 = FastInterpolations._compute_deriv1_coeffs(PolyFit{3}(), Val(:left), x4) + c2 = FastInterpolations._compute_deriv1_coeffs(PolyFit{1}(), LeftSide(), x2) + c3 = FastInterpolations._compute_deriv1_coeffs(PolyFit{2}(), LeftSide(), x3) + c4 = FastInterpolations._compute_deriv1_coeffs(PolyFit{3}(), LeftSide(), x4) # Warmup calls FastInterpolations._weighted_sum(c2, f2) @@ -1706,49 +1706,49 @@ end @testset "_estimate_endpoint_derivative (Unified API, D≤3)" begin # Test with uniform grid (Range) # Warmup - FastInterpolations._estimate_endpoint_derivative(xs_uniform, ys, Val(:left), PolyFit{1}()) - FastInterpolations._estimate_endpoint_derivative(xs_uniform, ys, Val(:left), PolyFit{2}()) - FastInterpolations._estimate_endpoint_derivative(xs_uniform, ys, Val(:left), PolyFit{3}()) + FastInterpolations._estimate_endpoint_derivative(xs_uniform, ys, LeftSide(), PolyFit{1}()) + FastInterpolations._estimate_endpoint_derivative(xs_uniform, ys, LeftSide(), PolyFit{2}()) + FastInterpolations._estimate_endpoint_derivative(xs_uniform, ys, LeftSide(), PolyFit{3}()) # D=1 uniform alloc_d1_uniform = @allocated FastInterpolations._estimate_endpoint_derivative( - xs_uniform, ys, Val(:left), PolyFit{1}() + xs_uniform, ys, LeftSide(), PolyFit{1}() ) @test alloc_d1_uniform <= ALLOC_THRESHOLD # D=2 uniform alloc_d2_uniform = @allocated FastInterpolations._estimate_endpoint_derivative( - xs_uniform, ys, Val(:left), PolyFit{2}() + xs_uniform, ys, LeftSide(), PolyFit{2}() ) @test alloc_d2_uniform <= ALLOC_THRESHOLD # D=3 uniform alloc_d3_uniform = @allocated FastInterpolations._estimate_endpoint_derivative( - xs_uniform, ys, Val(:left), PolyFit{3}() + xs_uniform, ys, LeftSide(), PolyFit{3}() ) @test alloc_d3_uniform <= ALLOC_THRESHOLD # Test with non-uniform grid (Vector) - also should be zero for D≤3 # Warmup - FastInterpolations._estimate_endpoint_derivative(xs_nonuniform, ys_nonuniform, Val(:left), PolyFit{1}()) - FastInterpolations._estimate_endpoint_derivative(xs_nonuniform, ys_nonuniform, Val(:left), PolyFit{2}()) - FastInterpolations._estimate_endpoint_derivative(xs_nonuniform, ys_nonuniform, Val(:left), PolyFit{3}()) + FastInterpolations._estimate_endpoint_derivative(xs_nonuniform, ys_nonuniform, LeftSide(), PolyFit{1}()) + FastInterpolations._estimate_endpoint_derivative(xs_nonuniform, ys_nonuniform, LeftSide(), PolyFit{2}()) + FastInterpolations._estimate_endpoint_derivative(xs_nonuniform, ys_nonuniform, LeftSide(), PolyFit{3}()) # D=1 non-uniform alloc_d1_nonuniform = @allocated FastInterpolations._estimate_endpoint_derivative( - xs_nonuniform, ys_nonuniform, Val(:left), PolyFit{1}() + xs_nonuniform, ys_nonuniform, LeftSide(), PolyFit{1}() ) @test alloc_d1_nonuniform <= ALLOC_THRESHOLD # D=2 non-uniform alloc_d2_nonuniform = @allocated FastInterpolations._estimate_endpoint_derivative( - xs_nonuniform, ys_nonuniform, Val(:left), PolyFit{2}() + xs_nonuniform, ys_nonuniform, LeftSide(), PolyFit{2}() ) @test alloc_d2_nonuniform <= ALLOC_THRESHOLD # D=3 non-uniform alloc_d3_nonuniform = @allocated FastInterpolations._estimate_endpoint_derivative( - xs_nonuniform, ys_nonuniform, Val(:left), PolyFit{3}() + xs_nonuniform, ys_nonuniform, LeftSide(), PolyFit{3}() ) @test alloc_d3_nonuniform <= ALLOC_THRESHOLD end @@ -1766,22 +1766,22 @@ end @testset "_compute_deriv1 (Uniform Grid, D>3)" begin # D=4 (QuarticFit) - verify correct result - result_d4 = FastInterpolations._compute_deriv1(PolyFit{4}(), Val(:left), f5, inv_h) + result_d4 = FastInterpolations._compute_deriv1(PolyFit{4}(), LeftSide(), f5, inv_h) @test isfinite(result_d4) # D=5 (QuinticFit) - verify correct result - result_d5 = FastInterpolations._compute_deriv1(PolyFit{5}(), Val(:left), f6, inv_h) + result_d5 = FastInterpolations._compute_deriv1(PolyFit{5}(), LeftSide(), f6, inv_h) @test isfinite(result_d5) end @testset "_compute_deriv1_coeffs (Non-Uniform Grid, D>3)" begin # D=4 - verify coefficients are finite - coeffs_d4 = FastInterpolations._compute_deriv1_coeffs(PolyFit{4}(), Val(:left), x5) + coeffs_d4 = FastInterpolations._compute_deriv1_coeffs(PolyFit{4}(), LeftSide(), x5) @test all(isfinite, coeffs_d4) @test length(coeffs_d4) == 5 # D=5 - verify coefficients are finite - coeffs_d5 = FastInterpolations._compute_deriv1_coeffs(PolyFit{5}(), Val(:left), x6) + coeffs_d5 = FastInterpolations._compute_deriv1_coeffs(PolyFit{5}(), LeftSide(), x6) @test all(isfinite, coeffs_d5) @test length(coeffs_d5) == 6 end @@ -1793,13 +1793,13 @@ end # D=4 uniform - verify correctness result_d4 = FastInterpolations._estimate_endpoint_derivative( - xs_large, ys_large, Val(:right), PolyFit{4}() + xs_large, ys_large, RightSide(), PolyFit{4}() ) @test result_d4 ≈ 2.0 atol=1e-10 # f'(0) = 0 # D=5 uniform - verify correctness result_d5 = FastInterpolations._estimate_endpoint_derivative( - xs_large, ys_large, Val(:right), PolyFit{5}() + xs_large, ys_large, RightSide(), PolyFit{5}() ) @test result_d5 ≈ 2.0 atol=1e-10 @@ -1809,10 +1809,10 @@ end end # Report allocations - alloc_d4 = alloc_est( xs_large, ys_large, Val(:right), PolyFit{4}()) - alloc_d4 = alloc_est( xs_large, ys_large, Val(:right), PolyFit{4}()) - alloc_d5 = alloc_est( xs_large, ys_large, Val(:right), PolyFit{5}()) - alloc_d5 = alloc_est( xs_large, ys_large, Val(:right), PolyFit{5}()) + alloc_d4 = alloc_est( xs_large, ys_large, RightSide(), PolyFit{4}()) + alloc_d4 = alloc_est( xs_large, ys_large, RightSide(), PolyFit{4}()) + alloc_d5 = alloc_est( xs_large, ys_large, RightSide(), PolyFit{5}()) + alloc_d5 = alloc_est( xs_large, ys_large, RightSide(), PolyFit{5}()) @test alloc_d4 <= ALLOC_THRESHOLD @test alloc_d5 <= ALLOC_THRESHOLD @@ -1823,11 +1823,11 @@ end ys_nonuniform_large = xs_nonuniform_large .^ 2 result_d4_nu = FastInterpolations._estimate_endpoint_derivative( - xs_nonuniform_large, ys_nonuniform_large, Val(:right), PolyFit{4}() + xs_nonuniform_large, ys_nonuniform_large, RightSide(), PolyFit{4}() ) @test result_d4_nu ≈ 2.0 atol=1e-10 - alloc_d4_nu = alloc_est( xs_nonuniform_large, ys_nonuniform_large, Val(:right), PolyFit{4}()) + alloc_d4_nu = alloc_est( xs_nonuniform_large, ys_nonuniform_large, RightSide(), PolyFit{4}()) @test alloc_d4_nu <= ALLOC_THRESHOLD end @@ -1912,7 +1912,7 @@ end y = [0.0, 1.0] bc_d3 = Deriv3(1.0) # Should passthrough unchanged - @test FastInterpolations.materialize_bc(bc_d3, x, y, Val(:left)) === bc_d3 + @test FastInterpolations.materialize_bc(bc_d3, x, y, LeftSide()) === bc_d3 # 5. Type promotion helpers @test FastInterpolations._promote_pointbc(Deriv1(1), Float64) isa Deriv1{Float64} @@ -1957,12 +1957,12 @@ end inv_h = 1.0 # Left endpoint - result_left = FastInterpolations._compute_deriv1(PolyFit{1}(), Val(:left), f2, inv_h) + result_left = FastInterpolations._compute_deriv1(PolyFit{1}(), LeftSide(), f2, inv_h) @test result_left isa ComplexF64 @test result_left ≈ (f2[2] - f2[1]) * inv_h # Right endpoint - result_right = FastInterpolations._compute_deriv1(PolyFit{1}(), Val(:right), f2, inv_h) + result_right = FastInterpolations._compute_deriv1(PolyFit{1}(), RightSide(), f2, inv_h) @test result_right isa ComplexF64 @test result_right ≈ (f2[2] - f2[1]) * inv_h end @@ -1973,7 +1973,7 @@ end inv_h = 1.0 # Right endpoint: (1, -4, 3) / 2 - result_right = FastInterpolations._compute_deriv1(PolyFit{2}(), Val(:right), f3, inv_h) + result_right = FastInterpolations._compute_deriv1(PolyFit{2}(), RightSide(), f3, inv_h) @test result_right isa ComplexF64 expected = (1 * f3[1] + (-4) * f3[2] + 3 * f3[3]) * (inv_h / 2) @test result_right ≈ expected @@ -1985,13 +1985,13 @@ end inv_h = 1.0 # Left endpoint: (-11, 18, -9, 2) / 6 - result_left = FastInterpolations._compute_deriv1(PolyFit{3}(), Val(:left), f4, inv_h) + result_left = FastInterpolations._compute_deriv1(PolyFit{3}(), LeftSide(), f4, inv_h) @test result_left isa ComplexF64 expected_left = (-11 * f4[1] + 18 * f4[2] + (-9) * f4[3] + 2 * f4[4]) * (inv_h / 6) @test result_left ≈ expected_left # Right endpoint: (-2, 9, -18, 11) / 6 - result_right = FastInterpolations._compute_deriv1(PolyFit{3}(), Val(:right), f4, inv_h) + result_right = FastInterpolations._compute_deriv1(PolyFit{3}(), RightSide(), f4, inv_h) @test result_right isa ComplexF64 expected_right = (-2 * f4[1] + 9 * f4[2] + (-18) * f4[3] + 11 * f4[4]) * (inv_h / 6) @test result_right ≈ expected_right @@ -2065,37 +2065,37 @@ end # LinearFit (D=1) - uniform + Complex d1_left = FastInterpolations._estimate_endpoint_derivative( - xs, ys, Val(:left), PolyFit{1}() + xs, ys, LeftSide(), PolyFit{1}() ) @test d1_left isa ComplexF64 d1_right = FastInterpolations._estimate_endpoint_derivative( - xs, ys, Val(:right), PolyFit{1}() + xs, ys, RightSide(), PolyFit{1}() ) @test d1_right isa ComplexF64 # QuadraticFit (D=2) should be exact for quadratic d2_left = FastInterpolations._estimate_endpoint_derivative( - xs, ys, Val(:left), PolyFit{2}() + xs, ys, LeftSide(), PolyFit{2}() ) @test d2_left isa ComplexF64 @test d2_left ≈ f_deriv(xs[1]) rtol=1e-10 d2_right = FastInterpolations._estimate_endpoint_derivative( - xs, ys, Val(:right), PolyFit{2}() + xs, ys, RightSide(), PolyFit{2}() ) @test d2_right isa ComplexF64 @test d2_right ≈ f_deriv(xs[end]) rtol=1e-10 # CubicFit (D=3) - uniform + Complex d3_left = FastInterpolations._estimate_endpoint_derivative( - xs, ys, Val(:left), PolyFit{3}() + xs, ys, LeftSide(), PolyFit{3}() ) @test d3_left isa ComplexF64 @test d3_left ≈ f_deriv(xs[1]) rtol=1e-8 d3_right = FastInterpolations._estimate_endpoint_derivative( - xs, ys, Val(:right), PolyFit{3}() + xs, ys, RightSide(), PolyFit{3}() ) @test d3_right isa ComplexF64 @test d3_right ≈ f_deriv(xs[end]) rtol=1e-8 @@ -2112,39 +2112,39 @@ end # LinearFit (D=1) - non-uniform + Complex (covers N=2 _weighted_sum) d1_left = FastInterpolations._estimate_endpoint_derivative( - xs, ys, Val(:left), PolyFit{1}() + xs, ys, LeftSide(), PolyFit{1}() ) @test d1_left isa ComplexF64 # QuadraticFit (D=2) should be exact for quadratic d2_left = FastInterpolations._estimate_endpoint_derivative( - xs, ys, Val(:left), PolyFit{2}() + xs, ys, LeftSide(), PolyFit{2}() ) @test d2_left isa ComplexF64 @test d2_left ≈ f_deriv(xs[1]) rtol=1e-10 # CubicFit (D=3) - non-uniform + Complex (covers N=4 _weighted_sum) d3_left = FastInterpolations._estimate_endpoint_derivative( - xs, ys, Val(:left), PolyFit{3}() + xs, ys, LeftSide(), PolyFit{3}() ) @test d3_left isa ComplexF64 @test d3_left ≈ f_deriv(xs[1]) rtol=1e-8 d3_right = FastInterpolations._estimate_endpoint_derivative( - xs, ys, Val(:right), PolyFit{3}() + xs, ys, RightSide(), PolyFit{3}() ) @test d3_right isa ComplexF64 @test d3_right ≈ f_deriv(xs[end]) rtol=1e-8 # PolyFit{4} - non-uniform + Complex (covers N=5 _weighted_sum generic) d4_left = FastInterpolations._estimate_endpoint_derivative( - xs, ys, Val(:left), PolyFit{4}() + xs, ys, LeftSide(), PolyFit{4}() ) @test d4_left isa ComplexF64 # PolyFit{5} - non-uniform + Complex (covers N=6 _weighted_sum generic) d5_left = FastInterpolations._estimate_endpoint_derivative( - xs, ys, Val(:left), PolyFit{5}() + xs, ys, LeftSide(), PolyFit{5}() ) @test d5_left isa ComplexF64 end diff --git a/test/test_show.jl b/test/test_show.jl index 8b4896a1..36a7a121 100644 --- a/test/test_show.jl +++ b/test/test_show.jl @@ -284,18 +284,18 @@ @test occursin("[1, 10]", verbose) # Integer bounds end - @testset "Side options (:left, :right)" begin + @testset "Side options (LeftSide, RightSide)" begin x_vec = collect(range(0.0, 1.0, 11)) y_short = sin.(x_vec) - itp_left = constant_interp(x_vec, y_short; side=:left) - itp_right = constant_interp(x_vec, y_short; side=:right) + itp_left = constant_interp(x_vec, y_short; side=LeftSide()) + itp_right = constant_interp(x_vec, y_short; side=RightSide()) verbose_left = sprint(show, MIME("text/plain"), itp_left) verbose_right = sprint(show, MIME("text/plain"), itp_right) - @test occursin(":left", verbose_left) - @test occursin(":right", verbose_right) + @test occursin("LeftSide", verbose_left) + @test occursin("RightSide", verbose_right) end @testset "Search policy formatting (direct)" begin @@ -463,7 +463,6 @@ # Access internal formatting functions directly to ensure full coverage @test FI._format_extrap(Val(:unknown_mode)) == "unknown" - @test FI._format_side(Val(:unknown_side)) == "unknown" @test FI._format_deriv_order(4) == "4th" @test FI._format_search(Linear()) == "Linear" @@ -819,7 +818,7 @@ compact_str = sprint(show, itp) @test occursin("ConstantInterpolantND", compact_str) @test occursin("11×15", compact_str) - @test occursin(":nearest", compact_str) + @test occursin("NearestSide", compact_str) # Verbose show (Range grids → no Search row) verbose_str = sprint(show, MIME("text/plain"), itp) @@ -849,16 +848,16 @@ x2 = range(0.0, 2.0, 15) data = [Float64(i + j) for i in 1:11, j in 1:15] - itp = constant_interp((x1, x2), data; side=(:left, :right)) + itp = constant_interp((x1, x2), data; side=(LeftSide(), RightSide())) # Compact: heterogeneous sides shown as comma-separated compact_str = sprint(show, itp) - @test occursin(":left,:right", compact_str) + @test occursin("LeftSide,RightSide", compact_str) # Verbose: should show tuple format for sides verbose_str = sprint(show, MIME("text/plain"), itp) - @test occursin(":left", verbose_str) - @test occursin(":right", verbose_str) + @test occursin("LeftSide", verbose_str) + @test occursin("RightSide", verbose_str) end @testset "ConstantInterpolantND show with heterogeneous extrapolation" begin @@ -916,20 +915,20 @@ FI = FastInterpolations # All same sides → single value - sides_same = (Val(:nearest), Val(:nearest)) - @test FI._short_side_name_nd(sides_same) == ":nearest" + sides_same = (NearestSide(), NearestSide()) + @test FI._short_side_name_nd(sides_same) == "NearestSide" # Different sides → comma-joined - sides_diff = (Val(:left), Val(:right)) - @test FI._short_side_name_nd(sides_diff) == ":left,:right" + sides_diff = (LeftSide(), RightSide()) + @test FI._short_side_name_nd(sides_diff) == "LeftSide,RightSide" # Three axes, all same - sides_3d_same = (Val(:left), Val(:left), Val(:left)) - @test FI._short_side_name_nd(sides_3d_same) == ":left" + sides_3d_same = (LeftSide(), LeftSide(), LeftSide()) + @test FI._short_side_name_nd(sides_3d_same) == "LeftSide" # Three axes, mixed - sides_3d_mixed = (Val(:left), Val(:nearest), Val(:right)) - @test FI._short_side_name_nd(sides_3d_mixed) == ":left,:nearest,:right" + sides_3d_mixed = (LeftSide(), NearestSide(), RightSide()) + @test FI._short_side_name_nd(sides_3d_mixed) == "LeftSide,NearestSide,RightSide" end @testset "LinearInterpolantND show with Vector grids (Search row)" begin diff --git a/test/test_type_stability.jl b/test/test_type_stability.jl index be682872..669f5631 100644 --- a/test/test_type_stability.jl +++ b/test/test_type_stability.jl @@ -573,7 +573,6 @@ using FastInterpolations: PolyFit, LinearFit, QuadraticFit, CubicFit end @testset "ND typed extrap — constant_interp" begin - # constant_interp has nested @_dispatch_side_nd, so constructor isn't @inferred # (same limitation as Symbol path). Test construction + eval separately. itp_none = constant_interp((x_nd, y_nd), data2d; extrap=NoExtrap()) @test itp_none isa ConstantInterpolantND