diff --git a/src/akima/akima_oneshot.jl b/src/akima/akima_oneshot.jl index 78f81615d..353a615a3 100644 --- a/src/akima/akima_oneshot.jl +++ b/src/akima/akima_oneshot.jl @@ -24,9 +24,11 @@ # Grid pre-normalized by the public `akima_interp` API via `_resolve_axis(x)` # before dispatching here; `_periodic_extend_1d` preserves the normalization. x_eff, y_ext, bc_eff, extrap_eff = _periodic_extend_1d(x, y, bc, extrap) - Tdy = _promote_eltype(_coeff_op, eltype(x_eff), Tv) + # Value-matched width: dy buffer + slope arithmetic run at `Tw` — see pchip_oneshot.jl. + Tw = _promote_grid_float(eltype(x_eff), Tv) + Tdy = _promote_eltype(_coeff_op, Tw, Tv) dy = acquire!(pool, Tdy, length(y_ext)) - _akima_slopes!(dy, x_eff, y_ext; bc = bc_eff) + _akima_slopes!(dy, x_eff, y_ext, Tw; bc = bc_eff) searcher = _resolve_search(x_eff, xq, search, hint) return _hermite_eval_at_point(x_eff, y_ext, dy, xq, extrap_eff, deriv, searcher) end @@ -47,9 +49,11 @@ end @boundscheck length(output) == length(x_query) || _throw_length_mismatch(length(x_query), length(output), "x_query", "output") x_eff, y_ext, bc_eff, extrap_eff = _periodic_extend_1d(x, y, bc, extrap) - Tdy = _promote_eltype(_coeff_op, eltype(x_eff), Tv) + # Value-matched width: dy buffer + slope arithmetic run at `Tw` — see pchip_oneshot.jl. + Tw = _promote_grid_float(eltype(x_eff), Tv) + Tdy = _promote_eltype(_coeff_op, Tw, Tv) dy = acquire!(pool, Tdy, length(y_ext)) - _akima_slopes!(dy, x_eff, y_ext; bc = bc_eff) + _akima_slopes!(dy, x_eff, y_ext, Tw; bc = bc_eff) searcher = _resolve_search(x_eff, x_query, search, hint) return _hermite_vector_loop!(output, x_eff, y_ext, dy, x_query, extrap_eff, deriv, searcher) end @@ -126,7 +130,8 @@ Outlier-robust, C\$^1\$ continuous. search::AbstractSearchPolicy = AutoSearch(), hint::Union{Nothing, Base.RefValue{Int}} = nothing ) where {Tg, Tv, Tq <: Real} - x = _resolve_axis(x) + # Value-matched Tg: Int/OneTo grid + Float32 data → Float32 axis. + x = _resolve_axis(x, _promote_grid_float(Tg, Tv)) extrap_eff = _resolve_extrap(extrap, bc, x, y) resolved = _resolve_coeffs(coeffs, x, xq) if resolved isa OnTheFly @@ -152,7 +157,7 @@ In-place Akima interpolation with outlier-robust slopes. search::AbstractSearchPolicy = AutoSearch(), hint::Union{Nothing, Base.RefValue{Int}} = nothing ) where {Tg, Tv, Tq <: Real} - x = _resolve_axis(x) + x = _resolve_axis(x, _promote_grid_float(Tg, Tv)) extrap_eff = _resolve_extrap(extrap, bc, x, y) resolved = _resolve_coeffs(coeffs, x, x_query) if resolved isa OnTheFly diff --git a/src/akima/akima_slopes.jl b/src/akima/akima_slopes.jl index 699a21ba7..83461a089 100644 --- a/src/akima/akima_slopes.jl +++ b/src/akima/akima_slopes.jl @@ -40,12 +40,18 @@ secant sequence. # Complexity O(n), single pass, zero allocation (writes into `dy`). """ +# Width-less form: delegate with the axis's own eltype — bit-identical to the +# historic raw behavior. One-shot PreCompute backends pass the value-matched `Tw`. +_akima_slopes!(dy::AbstractVector, x::AbstractVector, y::AbstractVector; bc::AbstractBC = NoBC()) = + _akima_slopes!(dy, x, y, eltype(x); bc) + function _akima_slopes!( dy::AbstractVector, - x::AbstractVector{Tg}, - y::AbstractVector; + x::AbstractVector, + y::AbstractVector, + ::Type{Tw}; bc::AbstractBC = NoBC() - ) where {Tg} + ) where {Tw} n = length(x) @assert n >= 2 "Akima requires at least 2 points" @assert length(y) == n "y length must match x" @@ -55,12 +61,12 @@ function _akima_slopes!( # 4-secant helper (cycle=2 for `:exclusive` yields a 2-secant alternation). if n == 2 if bc isa PeriodicBC - @inbounds dy[1] = _akima_local_4secant_periodic(x, y, 1, n, bc) - @inbounds dy[2] = _akima_local_4secant_periodic(x, y, 2, n, bc) + @inbounds dy[1] = _akima_local_4secant_periodic(Tw, x, y, 1, n, bc) + @inbounds dy[2] = _akima_local_4secant_periodic(Tw, x, y, 2, n, bc) return dy end @inbounds begin - δ = _forward_secant(x, y, 1) + δ = _forward_secant(Tw, x, y, 1) dy[1] = δ dy[2] = δ end @@ -73,13 +79,13 @@ function _akima_slopes!( # Wrap-aware path: every index is within K=5 stencil reach of the # join, so use the closed-cycle 4-secant formula at all 3 points. @inbounds for i in 1:3 - dy[i] = _akima_local_4secant_periodic(x, y, i, n, bc) + dy[i] = _akima_local_4secant_periodic(Tw, x, y, i, n, bc) end return dy end @inbounds begin - m1 = _forward_secant(x, y, 1) - m2 = _forward_secant(x, y, 2) + m1 = _forward_secant(Tw, x, y, 1) + m2 = _forward_secant(Tw, x, y, 2) dy[1] = m1 dy[2] = (m1 + m2) / 2 dy[3] = m2 @@ -101,9 +107,9 @@ function _akima_slopes!( # Simplification: extrapolate m sequence linearly. # Compute all n-1 secant slopes - @inbounds m1 = _forward_secant(x, y, 1) - @inbounds m2 = _forward_secant(x, y, 2) - @inbounds m3 = _forward_secant(x, y, 3) + @inbounds m1 = _forward_secant(Tw, x, y, 1) + @inbounds m2 = _forward_secant(Tw, x, y, 2) + @inbounds m3 = _forward_secant(Tw, x, y, 3) # Boundary virtual / wrapped secants for the LEFT side of the domain. # NoBC: linear extrapolation (Akima's original). @@ -112,8 +118,8 @@ function _akima_slopes!( # (closed cycle on n cells, with virtual seam cell). # The `_periodic_secant` abstraction absorbs both PeriodicBC variants. if bc isa PeriodicBC - @inbounds m_0 = _periodic_secant(x, y, 0, n, bc) # m[0] - @inbounds m_neg1 = _periodic_secant(x, y, -1, n, bc) # m[-1] + @inbounds m_0 = _periodic_secant(Tw, x, y, 0, n, bc) # m[0] + @inbounds m_neg1 = _periodic_secant(Tw, x, y, -1, n, bc) # m[-1] else @inbounds m_0 = 2 * m1 - m2 # virtual (NoBC) @inbounds m_neg1 = 3 * m1 - 2 * m2 # virtual (NoBC) @@ -133,7 +139,7 @@ function _akima_slopes!( m_k = m3 @inbounds for k in 3:(n - 2) - m_kp1 = _forward_secant(x, y, k + 1) + m_kp1 = _forward_secant(Tw, x, y, k + 1) dy[k] = _akima_weighted_slope(m_km2, m_km1, m_k, m_kp1) m_km2 = m_km1 m_km1 = m_k @@ -147,8 +153,8 @@ function _akima_slopes!( # PeriodicBC{:inclusive}: m[n]=m[1], m[n+1]=m[2] (closed cycle on n-1 cells). # PeriodicBC{:exclusive}: m[n]=seam, m[n+1]=m[1] (closed cycle on n cells). if bc isa PeriodicBC - @inbounds m_np1 = _periodic_secant(x, y, n, n, bc) # m[n] - @inbounds m_np2 = _periodic_secant(x, y, n + 1, n, bc) # m[n+1] + @inbounds m_np1 = _periodic_secant(Tw, x, y, n, n, bc) # m[n] + @inbounds m_np2 = _periodic_secant(Tw, x, y, n + 1, n, bc) # m[n+1] else m_np1 = 2 * m_k - m_km1 # virtual (NoBC) m_np2 = 3 * m_k - 2 * m_km1 # virtual (NoBC) diff --git a/src/cardinal/cardinal_oneshot.jl b/src/cardinal/cardinal_oneshot.jl index fe0bbef2e..8b4fdad5c 100644 --- a/src/cardinal/cardinal_oneshot.jl +++ b/src/cardinal/cardinal_oneshot.jl @@ -30,9 +30,12 @@ # `_CachedRange`/Vector) — the public `cardinal_interp` API pre-resolved # via `_resolve_axis(x)` before dispatching here, so no extra prep needed. x_eff, y_ext, bc_eff, extrap_eff = _periodic_extend_1d(x, y, bc, extrap) - Tdy = _promote_eltype(_coeff_op, eltype(x_eff), Tv) + # Value-matched width: dy buffer + slope arithmetic (incl. the `1 - tension` + # scale) run at `Tw` — see pchip_oneshot.jl. + Tw = _promote_grid_float(eltype(x_eff), Tv) + Tdy = _promote_eltype(_coeff_op, Tw, Tv) dy = acquire!(pool, Tdy, length(y_ext)) - _cardinal_slopes!(dy, x_eff, y_ext, tension; bc = bc_eff) + _cardinal_slopes!(dy, x_eff, y_ext, tension, Tw; bc = bc_eff) searcher = _resolve_search(x_eff, xq, search, hint) return _hermite_eval_at_point(x_eff, y_ext, dy, xq, extrap_eff, deriv, searcher) end @@ -54,9 +57,12 @@ end @boundscheck length(output) == length(x_query) || _throw_length_mismatch(length(x_query), length(output), "x_query", "output") x_eff, y_ext, bc_eff, extrap_eff = _periodic_extend_1d(x, y, bc, extrap) - Tdy = _promote_eltype(_coeff_op, eltype(x_eff), Tv) + # Value-matched width: dy buffer + slope arithmetic (incl. the `1 - tension` + # scale) run at `Tw` — see pchip_oneshot.jl. + Tw = _promote_grid_float(eltype(x_eff), Tv) + Tdy = _promote_eltype(_coeff_op, Tw, Tv) dy = acquire!(pool, Tdy, length(y_ext)) - _cardinal_slopes!(dy, x_eff, y_ext, tension; bc = bc_eff) + _cardinal_slopes!(dy, x_eff, y_ext, tension, Tw; bc = bc_eff) searcher = _resolve_search(x_eff, x_query, search, hint) return _hermite_vector_loop!(output, x_eff, y_ext, dy, x_query, extrap_eff, deriv, searcher) end @@ -144,7 +150,8 @@ Default `tension=0` is Catmull-Rom. C\$^1\$ continuous. search::AbstractSearchPolicy = AutoSearch(), hint::Union{Nothing, Base.RefValue{Int}} = nothing ) where {Tg, Tv, Tq <: Real} - x = _resolve_axis(x) + # Value-matched Tg: Int/OneTo grid + Float32 data → Float32 axis (tension follows). + x = _resolve_axis(x, _promote_grid_float(Tg, Tv)) tension_f = float(eltype(x))(tension) extrap_eff = _resolve_extrap(extrap, bc, x, y) resolved = _resolve_coeffs(coeffs, x, xq) @@ -172,7 +179,7 @@ In-place cardinal spline interpolation. search::AbstractSearchPolicy = AutoSearch(), hint::Union{Nothing, Base.RefValue{Int}} = nothing ) where {Tg, Tv, Tq <: Real} - x = _resolve_axis(x) + x = _resolve_axis(x, _promote_grid_float(Tg, Tv)) tension_f = float(eltype(x))(tension) extrap_eff = _resolve_extrap(extrap, bc, x, y) resolved = _resolve_coeffs(coeffs, x, x_query) diff --git a/src/cardinal/cardinal_slopes.jl b/src/cardinal/cardinal_slopes.jl index 222b0b6dd..dfc2bb8d0 100644 --- a/src/cardinal/cardinal_slopes.jl +++ b/src/cardinal/cardinal_slopes.jl @@ -24,30 +24,41 @@ Compute cardinal spline slopes in-place. # Complexity O(n), single pass, zero allocation (writes into `dy`). """ +# Width-less form: delegate with the axis's own eltype (historic behavior; the +# `float(one(Tw))` scale slot reproduces the old `one(Int) - tension::Float64` +# promotion for Int axes). One-shot PreCompute backends pass the value-matched +# `Tw` so slopes — including the `(1 - tension)` scale — are born at the value width. +_cardinal_slopes!(dy::AbstractVector, x::AbstractVector, y::AbstractVector, tension; bc::AbstractBC = NoBC()) = + _cardinal_slopes!(dy, x, y, tension, eltype(x); bc) + function _cardinal_slopes!( dy::AbstractVector, - x::AbstractVector{Tg}, + x::AbstractVector, y::AbstractVector, - tension; + tension, + ::Type{Tw}; bc::AbstractBC = NoBC() - ) where {Tg} + ) where {Tw} n = length(x) @assert n >= 2 "Cardinal spline requires at least 2 points" @assert length(y) == n "y length must match x" @assert length(dy) == n "dy length must match x" - scale = one(Tg) - tension + # Dimensionless scale at the value-matched width (`one(quantity)` strips + # units; `float` keeps Int axes on the historic Float64 arithmetic). + Tsc = typeof(float(one(Tw))) + scale = one(Tsc) - convert(Tsc, tension) # Special case: 2 points. PeriodicBC routes through the wrap-aware central # FD helper (see PCHIP n=2 note for rationale). if n == 2 if bc isa PeriodicBC - @inbounds dy[1] = _cardinal_boundary_slope(x, y, 1, n, scale, bc) - @inbounds dy[2] = _cardinal_boundary_slope(x, y, 2, n, scale, bc) + @inbounds dy[1] = _cardinal_boundary_slope(Tw, x, y, 1, n, scale, bc) + @inbounds dy[2] = _cardinal_boundary_slope(Tw, x, y, 2, n, scale, bc) return dy end @inbounds begin - δ = _forward_secant(x, y, 1) + δ = _forward_secant(Tw, x, y, 1) dy[1] = scale * δ dy[2] = scale * δ end @@ -55,18 +66,18 @@ function _cardinal_slopes!( end # Left endpoint: bc-dispatched helper. - @inbounds dy[1] = _cardinal_boundary_slope(x, y, 1, n, scale, bc) + @inbounds dy[1] = _cardinal_boundary_slope(Tw, x, y, 1, n, scale, bc) # Interior: central finite difference (K=3, no wrap needed). @inbounds for k in 2:(n - 1) - dy[k] = scale * _centered_secant(x, y, k) + dy[k] = scale * _centered_secant(Tw, x, y, k) end # Right endpoint: same bc-dispatched helper. PeriodicBC{:inclusive} yields # dy[n] == dy[1] automatically (closed-cycle symmetry); :exclusive yields # a different value using the seam secant — the helper handles both via # `_periodic_secant`/`_periodic_cell_width`. - @inbounds dy[n] = _cardinal_boundary_slope(x, y, n, n, scale, bc) + @inbounds dy[n] = _cardinal_boundary_slope(Tw, x, y, n, n, scale, bc) return dy end diff --git a/src/constant/constant_oneshot.jl b/src/constant/constant_oneshot.jl index bcf2f235d..fcc7382d6 100644 --- a/src/constant/constant_oneshot.jl +++ b/src/constant/constant_oneshot.jl @@ -213,8 +213,10 @@ vals = constant_interp(x, y, sorted_queries; search=LinearBinarySearch(linear_wi @boundscheck length(y) == length(x) || throw(ArgumentError("x and y must have same length")) # Surface-level BC-aware resolvers (zero-alloc reference wrapping). BC info - # lives in axis type after resolution → searcher uses `NoBC()`. - x_eff = _resolve_axis(x, bc) + # lives in axis type after resolution → searcher uses `NoBC()`. Raw Tg: + # the selection kernel keeps natural promotion (no float forcing) — an Int + # range stays `_CachedRange{Int}`, mirroring the ND constant scalar rule. + x_eff = _resolve_axis(x, bc, Tg) y_eff = _resolve_data(y, bc) extrap_eff = _resolve_extrap(extrap, bc, x_eff, y_eff) searcher = _resolve_search(x_eff, xi, search, hint) @@ -266,8 +268,9 @@ function constant_interp!( @assert length(y) == length(x) "x and y must have same length" @assert length(output) == length(x_targets) "output must match x_targets length" - # Surface-level BC-aware resolvers (same template as Linear oneshot). - x_eff = _resolve_axis(x, bc) + # Surface-level BC-aware resolvers (same template as Linear oneshot); + # raw Tg — selection kernel keeps natural promotion (no float forcing). + x_eff = _resolve_axis(x, bc, eltype(x)) y_eff = _resolve_data(y, bc) extrap_eff = _resolve_extrap(extrap, bc, x_eff, y_eff) searcher = _resolve_search(x_eff, x_targets, search, nothing) diff --git a/src/constant/nd/constant_nd_adjoint_types.jl b/src/constant/nd/constant_nd_adjoint_types.jl index fed64f3e0..26f61dec9 100644 --- a/src/constant/nd/constant_nd_adjoint_types.jl +++ b/src/constant/nd/constant_nd_adjoint_types.jl @@ -74,7 +74,7 @@ struct ConstantAdjointND{ Tg, N, B <: NTuple{N, AbstractBC}, EP <: Tuple{Vararg{AbstractExtrap, N}}, SD <: Tuple{Vararg{AbstractSide, N}}, Tq, } - grids_c = map((g, bc) -> _convert_copy(_cache_axis(g, bc, Tg), Tg), grids, bcs) + grids_c = map((g, bc, T) -> _convert_copy(_cache_axis(g, bc, T), T), grids, bcs, ntuple(_ -> Tg, Val(N))) return new{Tg, N, typeof(grids_c), B, EP, SD, Tq}(grids_c, bcs, extraps, sides, anchors, grid_size) end end diff --git a/src/constant/nd/constant_nd_oneshot.jl b/src/constant/nd/constant_nd_oneshot.jl index 2ac88b868..02d5ed91c 100644 --- a/src/constant/nd/constant_nd_oneshot.jl +++ b/src/constant/nd/constant_nd_oneshot.jl @@ -29,7 +29,11 @@ function _constant_interp_nd_oneshot( ops::NTuple{N, AbstractEvalOp}, hints = nothing ) where {Tv, N} - grids_eff = map(_resolve_axis, grids, bcs) + # Selection kernel: no x·y arithmetic → RAW grid eltype (no float forcing), mirroring + # `_nd_promote_grids_raw`/batch/persistent. All-Int stays Int; the output follows the + # natural promote_type(grid, data, query) — e.g. Int grid + Float32 query → Float32. + Tg = _promote_grid_eltype(grids) + grids_eff = _resolve_axes(grids, bcs, Tg) # @generated static-Tg unroll (no Type-captured closure) # Bare GridIdx(k).val is NaN → resolve to the grid coordinate for the value kernel (search still uses .idx). query = map(_resolve_grididx, query, grids_eff) # Validate AND promote per axis: an in-domain NoExtrap axis becomes InBounds for the lean diff --git a/src/constant/nd/constant_nd_types.jl b/src/constant/nd/constant_nd_types.jl index 28b8a214d..2955024de 100644 --- a/src/constant/nd/constant_nd_types.jl +++ b/src/constant/nd/constant_nd_types.jl @@ -72,7 +72,7 @@ struct ConstantInterpolantND{ bcs::NTuple{N, AbstractBC} = ntuple(_ -> NoBC(), Val(N)), store::StorePolicy = StorePolicy() ) where {Tg, Tv, N} - grids_c = map((g, bc) -> _own_or_ref_axis(_cache_axis(g, bc, Tg), Tg, store), grids, bcs) + grids_c = map((g, bc, T) -> _own_or_ref_axis(_cache_axis(g, bc, T), T, store), grids, bcs, ntuple(_ -> Tg, Val(N))) data_c = _own_or_ref_data(data, store) return new{Tg, Tv, N, typeof(grids_c), typeof(extraps), typeof(sides), typeof(searches), typeof(data_c)}( grids_c, data_c, extraps, sides, searches diff --git a/src/core/cached_range.jl b/src/core/cached_range.jl index 423a5d15f..ac2c57b8d 100644 --- a/src/core/cached_range.jl +++ b/src/core/cached_range.jl @@ -201,19 +201,58 @@ end @inline _get_h(x::AbstractRange, ::Int) = step(x) @inline _get_inv_h(x::AbstractRange, ::Int) = inv(step(x)) +# ── Width-first forms (Range hierarchy) — see cached_vector.jl for the contract. +# `_CachedRange` reuses the cached reciprocal (convert is a no-op once the axis +# is value-matched, e.g. `_CachedRange{Tw}` from `_resolve_axis(x, Tw)`); a raw +# range differences via `step()` in its own eltype, converts the span, then divides. +@inline _get_inv_h(::Type{Tw}, x::_CachedRange, i::Int) where {Tw} = + convert(_promote_eltype(_inv_op, Tw), _get_inv_h(x, i)) +@inline _get_inv_2cell(::Type{Tw}, x::_CachedRange, i::Int) where {Tw} = + convert(_promote_eltype(_inv_op, Tw), _get_inv_2cell(x, i)) +@inline _get_inv_h(::Type{Tw}, x::AbstractRange, ::Int) where {Tw} = + inv(convert(Tw, step(x))) +# Search-result form: endpoints ignored — the cached (or step-derived) reciprocal wins. +@inline _get_inv_h(::Type{Tw}, x::_CachedRange, i::Int, ::Real, ::Real) where {Tw} = + _get_inv_h(Tw, x, i) + # ======================================== # `_resolve_axis` — one-shot Range wrapping # ======================================== @inline _resolve_axis(x::AbstractRange) = _to_float(x, float(eltype(x))) @inline _resolve_axis(x::AbstractRange, ::AbstractBC) = _to_float(x, float(eltype(x))) @inline _resolve_axis(c::_CachedRange) = c +# Tg-typed 2-arg (no BC): value-matched one-shot normalization — an Int/OneTo grid +# beside Float32 data resolves to `_CachedRange{Float32}`, not the blind Float64. +@inline _resolve_axis(x::AbstractRange, ::Type{Tg}) where {Tg} = _to_float(x, Tg) +@inline _resolve_axis(c::_CachedRange, ::Type{Tg}) where {Tg} = _convert_copy(c, Tg) -# `:exclusive` raw-input one-shot path — wrap into `_ExclusivePeriodicAxis`. -@inline function _resolve_axis(x::AbstractRange, bc::PeriodicBC{:exclusive}) - bc_resolved = _resolve_bc_period(x, bc) - return _ExclusivePeriodicAxis(_to_float(x, float(eltype(x))), bc_resolved.period) +# Shared convert-first wrapper for EVERY `:exclusive` axis site (one-shot + persistent, Range + +# Vector — see cached_vector.jl). Resolving the period AGAINST the already-converted inner +# (`_resolve_bc_period` normalizes it to the inner's eltype) is what stops a wider period literal +# from re-widening a value-matched Tg axis through the `_ExclusivePeriodicAxis` ctor's seam +# promote. Callers pass the inner at the width they want (`_to_float(x, Tg)` / raw / `_CachedVector`). +@inline function _wrap_exclusive(inner, bc::PeriodicBC{:exclusive}) + bc_resolved = _resolve_bc_period(inner, bc) + return _ExclusivePeriodicAxis(inner, bc_resolved.period) end +# `:exclusive` raw-input one-shot path — wrap into `_ExclusivePeriodicAxis` (natural float width). +@inline _resolve_axis(x::AbstractRange, bc::PeriodicBC{:exclusive}) = + _wrap_exclusive(_to_float(x, float(eltype(x))), bc) + +# 3-arg Tg-aware one-shot resolution — value-matched grid float so an Int/OneTo grid beside +# Float32 data floats to Float32 (not the blind `float(eltype)`=Float64), matching the persistent +# path and the natural `promote_type(grid, data, query)` output. Mirrors `_cache_axis(x, bc, Tg)`. +@inline _resolve_axis(x::AbstractRange, ::AbstractBC, ::Type{Tg}) where {Tg} = _to_float(x, Tg) +@inline _resolve_axis(c::_CachedRange, ::AbstractBC, ::Type{Tg}) where {Tg} = _convert_copy(c, Tg) +# `:exclusive` — convert-first to Tg, then `_wrap_exclusive` resolves the period against the +# Tg-typed inner. Diagonal (`_CachedRange` × `:exclusive`) is load-bearing against ambiguity +# between the two arms above; `_convert_copy` keeps the Tg-matched wrap allocation-free. +@inline _resolve_axis(x::AbstractRange, bc::PeriodicBC{:exclusive}, ::Type{Tg}) where {Tg} = + _wrap_exclusive(_to_float(x, Tg), bc) +@inline _resolve_axis(c::_CachedRange, bc::PeriodicBC{:exclusive}, ::Type{Tg}) where {Tg} = + _wrap_exclusive(_convert_copy(c, Tg), bc) + # ======================================== # `_cache_axis` — persistent-path Range wrapping # ======================================== @@ -224,25 +263,20 @@ end @inline _cache_axis(c::_CachedRange) = c @inline _cache_axis(c::_CachedRange, ::AbstractBC) = c -# `:exclusive` 2-arg variants — produce `_ExclusivePeriodicAxis`. -@inline function _cache_axis(x::AbstractRange, bc::PeriodicBC{:exclusive}) - bc_resolved = _resolve_bc_period(x, bc) - return _ExclusivePeriodicAxis(_to_float(x, float(eltype(x))), bc_resolved.period) -end -@inline function _cache_axis(c::_CachedRange, bc::PeriodicBC{:exclusive}) - bc_resolved = _resolve_bc_period(c, bc) - return _ExclusivePeriodicAxis(c, bc_resolved.period) -end +# `:exclusive` 2-arg variants — produce `_ExclusivePeriodicAxis` (natural float width). +@inline _cache_axis(x::AbstractRange, bc::PeriodicBC{:exclusive}) = + _wrap_exclusive(_to_float(x, float(eltype(x))), bc) +@inline _cache_axis(c::_CachedRange, bc::PeriodicBC{:exclusive}) = _wrap_exclusive(c, bc) # 3-arg Tg-aware. Raw Range respects Tg via `_to_float`. Pre-wrapped passes # through — downstream `_convert_copy(_, Tg)` enforces Tg (intentional contract; # see DISPATCH TABLE in `periodic_axis.jl`). @inline _cache_axis(x::AbstractRange, ::AbstractBC, ::Type{Tg}) where {Tg} = _to_float(x, Tg) @inline _cache_axis(c::_CachedRange, bc::AbstractBC, ::Type{Tg}) where {Tg} = _cache_axis(c, bc) -@inline function _cache_axis(x::AbstractRange, bc::PeriodicBC{:exclusive}, ::Type{Tg}) where {Tg} - bc_resolved = _resolve_bc_period(x, bc) - return _ExclusivePeriodicAxis(_to_float(x, Tg), bc_resolved.period) -end +# Convert-first delegate (mirrors the Vector-side 3-arg): the 2-arg `_CachedRange` +# arm resolves the period against the Tg-typed axis, so the period follows Tg. +@inline _cache_axis(x::AbstractRange, bc::PeriodicBC{:exclusive}, ::Type{Tg}) where {Tg} = + _cache_axis(_to_float(x, Tg), bc) @inline _cache_axis(c::_CachedRange, bc::PeriodicBC{:exclusive}, ::Type{Tg}) where {Tg} = _cache_axis(c, bc) # ======================================== @@ -251,3 +285,8 @@ end # Range types have stack-only `_CachedRange` — no pool buffer needed. @inline _cache_axis_pooled(_, x::AbstractRange) = _to_float(x, float(eltype(x))) @inline _cache_axis_pooled(_, x::_CachedRange) = x +# 3-arg Tg-aware. Also load-bearing for dispatch: `AbstractRange <: AbstractVector`, so +# without these a Range would fall into the pooled-Vector 3-arg overload (cached_vector.jl) +# and lose its stack `_CachedRange` form. +@inline _cache_axis_pooled(_, x::AbstractRange, ::Type{Tg}) where {Tg} = _to_float(x, Tg) +@inline _cache_axis_pooled(_, x::_CachedRange, ::Type{Tg}) where {Tg} = _convert_copy(x, Tg) diff --git a/src/core/cached_vector.jl b/src/core/cached_vector.jl index 5b031288d..13fda3736 100644 --- a/src/core/cached_vector.jl +++ b/src/core/cached_vector.jl @@ -113,10 +113,10 @@ end # x type + Tg → result (notes) # ───────────────────────────────────────────────────────────────────────── # Vector{Tg} + Tg → pool-backed `_CachedVector{Tg}` -# (`_to_float` identity → 2-arg pool wrap) -# Vector{S} + Tg, S ≠ Tg → pool-backed `_CachedVector{Tg}` ⚠️ warns once -# (`_to_float` heap-allocates `Tg.(x)` then -# 2-arg builds the wrapper) +# (same-eltype overload → 2-arg pool wrap) +# Vector{S} + Tg, S ≠ Tg → pool-backed `_CachedVector{Tg}`, zero heap +# (converting `copyto!` into a POOLED buffer, +# then 2-arg builds the wrapper) # AbstractRange + Tg → `_CachedRange{Tg}` (zero pool, immutable # struct; eltype always Tg) # _CachedRange{Tg} + Tg → identity passthrough @@ -132,10 +132,21 @@ end # # Bottom line: `_cache_axis_pooled(pool, x, Tg)` *always* returns a wrapper # with eltype `Tg`. The fast paths (same-eltype Vector/Range/wrapped) are -# zero-alloc; eltype mismatches incur a one-time conversion alloc. +# zero-alloc; wrapped-type eltype mismatches incur a one-time conversion alloc. @inline _cache_axis_pooled(pool, x, ::Type{Tg}) where {Tg} = _cache_axis_pooled(pool, _to_float(x, Tg)) +# Vector overloads (beat the catch-all): same-eltype → plain 2-arg pool wrap; +# mismatched eltype (e.g. Int grid under a value-matched Tg) → converting +# `copyto!` into a POOLED buffer, then 2-arg wrap — warm one-shots stay zero-alloc. +@inline _cache_axis_pooled(pool, x::AbstractVector{Tg}, ::Type{Tg}) where {Tg} = + _cache_axis_pooled(pool, x) +@inline function _cache_axis_pooled(pool, x::AbstractVector, ::Type{Tg}) where {Tg} + buf = acquire!(pool, Tg, length(x)) + copyto!(buf, x) + return _cache_axis_pooled(pool, buf) +end + # ======================================== # _get_h / _get_inv_h accessors (Vector hierarchy) # ======================================== @@ -177,6 +188,37 @@ end @inline _get_h(::AbstractVector, ::Int, xL::Real, xR::Real) = xR - xL @inline _get_inv_h(::AbstractVector, ::Int, xL::Real, xR::Real) = inv(xR - xL) +# ── Width-first forms: `_get_*(Tw, x, i)` ──────────────────────────────────── +# `Tw` = value-matched coordinate width, computed ONCE at the caller's surface +# (`_promote_grid_float(Tg, Tv)`). Raw axes difference in their OWN eltype first +# (Int spans are exact; coordinates may exceed `Tw`'s ulp, so convert the SPAN, +# never the endpoints), then divide — the reciprocal is BORN at `Tw`: no +# `inv(Int)::Float64` minting beside narrower data. Wrapped axes (rows below and +# in cached_range.jl) reuse their cached reciprocal; the convert is a no-op once +# the axis is value-matched. `_get_h` needs no raw/wrapped split — converting +# the span is exact either way, and the inner call dispatches to the cache. +@inline Base.@propagate_inbounds _get_h(::Type{Tw}, x::AbstractVector, i::Int) where {Tw} = + convert(Tw, _get_h(x, i)) +@inline Base.@propagate_inbounds _get_inv_h(::Type{Tw}, x::AbstractVector, i::Int) where {Tw} = + inv(_get_h(Tw, x, i)) +@inline Base.@propagate_inbounds _get_inv_2cell(::Type{Tw}, x::AbstractVector, i::Int) where {Tw} = + @inbounds inv(convert(Tw, x[i + 1] - x[i - 1])) + +# Dispatch shields — `_CachedVector <: AbstractVector` must keep its cached +# reciprocal (width-convert the cached value instead of re-dividing). +@inline Base.@propagate_inbounds _get_inv_h(::Type{Tw}, x::_CachedVector, i::Int) where {Tw} = + convert(_promote_eltype(_inv_op, Tw), _get_inv_h(x, i)) +@inline Base.@propagate_inbounds _get_inv_2cell(::Type{Tw}, x::_CachedVector, i::Int) where {Tw} = + @inbounds inv(convert(Tw, x.h[i - 1] + x.h[i])) + +# Width-first SEARCH-RESULT form `(Tw, x, idx, xL, xR)` — same contract as above, +# but the span comes from the search endpoints (raw axes never re-index). The +# `_CachedVector` shield ignores the endpoints and reuses its cached reciprocal. +@inline _get_inv_h(::Type{Tw}, ::AbstractVector, ::Int, xL::Real, xR::Real) where {Tw} = + inv(convert(Tw, xR - xL)) +@inline Base.@propagate_inbounds _get_inv_h(::Type{Tw}, x::_CachedVector, idx::Int, ::Real, ::Real) where {Tw} = + _get_inv_h(Tw, x, idx) + # Persistent axis wrapping is split into two stages (see `periodic_axis.jl`): # - outer surface API: `_cache_axis(x, bc)` — bc-aware wrap, zero-copy # of buffer (Vector → `_CachedVector`, Range → `_CachedRange`, @@ -195,12 +237,27 @@ end @inline _resolve_axis(x::AbstractVector) = x @inline _resolve_axis(x::AbstractVector, ::AbstractBC) = x @inline _resolve_axis(c::_CachedVector) = c +# Tg-typed 2-arg (no BC): vectors stay raw — Tg only steers Range axes (the search +# promote-compares raw vectors; eager conversion would allocate). +@inline _resolve_axis(x::AbstractVector, ::Type{Tg}) where {Tg} = x +@inline _resolve_axis(c::_CachedVector, ::Type{Tg}) where {Tg} = c -# `:exclusive` raw-input one-shot path — wrap into `_ExclusivePeriodicAxis`. -@inline function _resolve_axis(x::AbstractVector, bc::PeriodicBC{:exclusive}) - bc_resolved = _resolve_bc_period(x, bc) - return _ExclusivePeriodicAxis(x, bc_resolved.period) -end +# `:exclusive` raw-input one-shot path — raw inner (2-arg has no Tg; the ctor widens against the +# period). `_wrap_exclusive` is the shared convert-first wrapper defined in cached_range.jl. +@inline _resolve_axis(x::AbstractVector, bc::PeriodicBC{:exclusive}) = _wrap_exclusive(x, bc) + +# 3-arg Tg-aware one-shot resolution — a raw Vector stays raw for NON-periodic axes (the search +# promote-compares it, no eager conversion), so `Tg` only steers Range axes there. +@inline _resolve_axis(x::AbstractVector, ::AbstractBC, ::Type{Tg}) where {Tg} = x +@inline _resolve_axis(c::_CachedVector, ::AbstractBC, ::Type{Tg}) where {Tg} = c +# `:exclusive` is the exception: the wrapped axis eltype is `promote(inner, period)`, so a raw Int +# inner + `float(Int)`=Float64 period forces a Float64 axis past a value-matched Float32 witness +# (crash). Convert-first to Tg — mirrors the Range arm + persistent `_cache_axis` — so the period +# follows Tg. Diagonals are load-bearing against ambiguity between these two and the arms above. +@inline _resolve_axis(x::AbstractVector, bc::PeriodicBC{:exclusive}, ::Type{Tg}) where {Tg} = + _wrap_exclusive(_to_float(x, Tg), bc) +@inline _resolve_axis(c::_CachedVector, bc::PeriodicBC{:exclusive}, ::Type{Tg}) where {Tg} = + _wrap_exclusive(_convert_copy(c, Tg), bc) # ======================================== # `_cache_axis` — persistent-path Vector wrapping @@ -213,14 +270,9 @@ end @inline _cache_axis(c::_CachedVector, ::AbstractBC) = c # `:exclusive` 2-arg variants — produce `_ExclusivePeriodicAxis(_CachedVector, ·)`. -@inline function _cache_axis(x::AbstractVector, bc::PeriodicBC{:exclusive}) - bc_resolved = _resolve_bc_period(x, bc) - return _ExclusivePeriodicAxis(_CachedVector(x), bc_resolved.period) -end -@inline function _cache_axis(c::_CachedVector, bc::PeriodicBC{:exclusive}) - bc_resolved = _resolve_bc_period(c, bc) - return _ExclusivePeriodicAxis(c, bc_resolved.period) -end +@inline _cache_axis(x::AbstractVector, bc::PeriodicBC{:exclusive}) = + _wrap_exclusive(_CachedVector(x), bc) +@inline _cache_axis(c::_CachedVector, bc::PeriodicBC{:exclusive}) = _wrap_exclusive(c, bc) # 3-arg Tg-aware. Raw Vector respects Tg via `_to_float` then wraps. # Pre-wrapped `_CachedVector` passes through — downstream `_convert_copy(_, Tg)` diff --git a/src/core/nd_utils.jl b/src/core/nd_utils.jl index 1a619941d..efe5cb4d4 100644 --- a/src/core/nd_utils.jl +++ b/src/core/nd_utils.jl @@ -1138,6 +1138,20 @@ Compute local cell parameters for all axes via `_get_h(grid, idx)` / # (N=0 edge: `float(promote_type())` = `float(Union{})` throws; the ntuples are # empty there, so this placeholder is never used.) Tg = N == 0 ? Float64 : float(_promote_grid_eltype(grids)) + return _compute_all_local_params(q_evals, grids, indices, Ls, Tg) +end + +# Data-aware form: the caller supplies the width type `Tg` (value-matched, e.g. +# `_promote_grid_float(grid eltype, Tv)` — Int grid + Float32 data → Float32), so raw +# Int axes don't widen the eval to Float64 via `inv(Int)`. `dLs` deliberately keep +# their natural `q - L` promotion — converting them would strip Dual-query partials. +@inline function _compute_all_local_params( + q_evals::Tuple{Vararg{Real, N}}, + grids::Tuple{Vararg{AbstractVector, N}}, + indices::NTuple{N, Int}, + Ls::Tuple{Vararg{Real, N}}, + ::Type{Tg}, + ) where {N, Tg} hs = ntuple(Val(N)) do d @inbounds convert(Tg, _get_h(grids[d], indices[d])) end @@ -1227,6 +1241,35 @@ Generates unrolled `(_convert_grid(grids[1], Tg), _convert_grid(grids[2], Tg), . return :(($(exprs...),)) end +# ── Static-Tg tuple maps for the one-shot hot paths (@generated) ───────────── +# `map(f, grids, ntuple(_ -> Tg, Val(N)))` re-captures the Type witness in the +# ntuple closure: under a degraded-inference context (a long-lived test worker) +# the tuple elements decay to `DataType` and every per-axis call goes through +# dynamic dispatch — 32 B/axis on Julia 1.12 CI, 60-90 KB downstream on LTS. +# These unroll at codegen with `Tg` as a STATIC signature parameter: no closure +# and no runtime `Type` value exist on any Julia version. + +# Pooled value-matched wrap (cubic/quadratic PreCompute scalar backends). +@generated function _cache_axes_pooled(pool, grids::NTuple{N, AbstractVector}, ::Type{Tg}) where {N, Tg} + exprs = [:(_cache_axis_pooled(pool, grids[$i], Tg)) for i in 1:N] + return :(($(exprs...),)) +end + +# BC-aware one-shot resolve (linear/constant/hetero OnTheFly surfaces). +@generated function _resolve_axes(grids::NTuple{N, AbstractVector}, bcs, ::Type{Tg}) where {N, Tg} + exprs = [:(_resolve_axis(grids[$i], bcs[$i], Tg)) for i in 1:N] + return :(($(exprs...),)) +end + +# Width-typed reciprocal spans from search results (linear ND scalar one-shot). +# Span-first via the width-first 5-arg `_get_inv_h` rows: raw axes difference in +# their own eltype, convert the span once, divide at `Tg` — the reciprocal is +# born at `Tg` (an Int axis would otherwise mint Float64 via `inv(Int)`). +@generated function _typed_inv_hs(grids::NTuple{N, AbstractVector}, idxs, Ls, Rs, ::Type{Tg}) where {N, Tg} + exprs = [:(_get_inv_h(Tg, grids[$i], idxs[$i], Ls[$i], Rs[$i])) for i in 1:N] + return :(($(exprs...),)) +end + """ _nd_promote_grids(grids, data) -> (grids_typed, Tg, Tv, Tz) @@ -1245,15 +1288,49 @@ grids_typed, _, _, _ = _nd_promote_grids(grids, data) # grid-only (batch dispa grids_typed, Tg, Tv, Tz = _nd_promote_grids(grids, data) # full (oneshot/build) ``` """ +@inline _nd_promote_grids( + grids::NTuple{N, AbstractVector}, + data::AbstractArray{Tv_raw, N} +) where {Tv_raw, N} = _nd_promote_grids(grids, data, Tv_raw) + +# 3-arg form: `Tv_extra` widens the value space beyond `eltype(data)` BEFORE the grid +# value-match — Hermite's value space is data ∪ partials (Float32 data + Float64 +# partials must give a Float64 grid, matching the one-shot rule). 2-arg delegates +# with `Tv_extra = eltype(data)` (neutral). @inline function _nd_promote_grids( + grids::NTuple{N, AbstractVector}, + data::AbstractArray{Tv_raw, N}, + ::Type{Tv_extra} + ) where {Tv_raw, Tv_extra, N} + Tv_all = promote_type(Tv_raw, Tv_extra) + # Value-matched grid float (1D rule): Int/OneTo grid + Float32 data → Float32 grid, so the + # cheap grid converts and the O(nᴺ) data aliases under copy=false. The old grid-eltype-only + # `float(...)` gave Float64 and dragged Tv (and the data, via `Tv.(data)`) up with it. + Tg = _promote_grid_float(_promote_grid_eltype(grids), Tv_all) + grids_typed = _convert_grids_typed(grids, Tg) + Tv = _value_type(Tv_all, Tg) + Tz = _promote_eltype(Tv, Tg) + return grids_typed, Tg, Tv, Tz +end + +""" + _nd_promote_types(grids, data) -> (Tg, Tv, Tz) + +Type-only counterpart of [`_nd_promote_grids`](@ref): computes the same +`(Tg, Tv, Tz)` from the grid/data *types* alone (via the `@generated` +`_promote_grid_eltype`), skipping the `_convert_grids_typed` allocation. Callers +that value-match each axis downstream (the OnTheFly hetero one-shot, whose inner +1D one-shots resolve their own axes) use this to avoid an eager grid convert on +raw Int/Rational axes. +""" +@inline function _nd_promote_types( grids::NTuple{N, AbstractVector}, data::AbstractArray{Tv_raw, N} ) where {Tv_raw, N} - Tg = float(_promote_grid_eltype(grids)) - grids_typed = _convert_grids_typed(grids, Tg) + Tg = _promote_grid_float(_promote_grid_eltype(grids), Tv_raw) Tv = _value_type(Tv_raw, Tg) Tz = _promote_eltype(Tv, Tg) - return grids_typed, Tg, Tv, Tz + return Tg, Tv, Tz end """ diff --git a/src/core/periodic_axis.jl b/src/core/periodic_axis.jl index 0f7915cba..9d5130e3a 100644 --- a/src/core/periodic_axis.jl +++ b/src/core/periodic_axis.jl @@ -223,6 +223,14 @@ end @inline Base.@propagate_inbounds _get_inv_h(g::_ExclusivePeriodicAxis, idx::Int) = idx < length(g.inner) ? _get_inv_h(g.inner, idx) : @inbounds(inv(g._x_max - g.inner[idx])) +# Width-first shield (`_ExclusivePeriodicAxis <: AbstractVector` — must not fall +# into the raw-span row in cached_vector.jl): interior cells thread `Tw` into the +# inner's own width-first row; the seam converts its span first (span-then-convert, +# see cached_vector.jl), so the reciprocal is born at `Tw`. +@inline Base.@propagate_inbounds _get_inv_h(::Type{Tw}, g::_ExclusivePeriodicAxis, idx::Int) where {Tw} = + idx < length(g.inner) ? _get_inv_h(Tw, g.inner, idx) : + @inbounds(inv(convert(Tw, g._x_max - g.inner[idx]))) + # No `_CachedRange`-specific specialization for `_get_h` / `_get_inv_h`: the # generic wrapper overload above already does the right thing for both # `_CachedRange` and `_CachedVector` / `Vector` inners — interior cells @@ -261,6 +269,20 @@ end @inline Base.@propagate_inbounds _get_inv_h(g::_ExclusivePeriodicAxis, idx::Int, xL::Real, xR::Real) = idx < length(g.inner) ? _get_inv_h(g.inner, idx, xL, xR) : inv(xR - xL) +# Width-first search-result form: interior delegates to the inner axis's row +# (cached reciprocal when wrapped); the seam cell has no stored width — span-first +# from the search endpoints (`xR == g._x_max`), reciprocal born at `Tw`. +@inline Base.@propagate_inbounds function _get_inv_h( + ::Type{Tw}, + g::_ExclusivePeriodicAxis, + idx::Int, + xL::Real, + xR::Real, + ) where {Tw} + return idx < length(g.inner) ? _get_inv_h(Tw, g.inner, idx, xL, xR) : + inv(convert(Tw, xR - xL)) +end + # `_alpha_of` for the wrapper: seam-aware so the value computation shares a # denominator with the 4-arg `_get_inv_h` at the seam cell. # - Interior cell (`R != g._x_max`): defer to inner so `_CachedRange` uses @@ -405,6 +427,12 @@ end # `_CachedRange`/`_CachedVector` → wrapper) live in their owner files # (`cached_range.jl` / `cached_vector.jl`). @inline _resolve_axis(g::_ExclusivePeriodicAxis) = g +# 2-arg/3-arg `:exclusive`: a pre-wrapped axis passes through — do NOT re-wrap. Without these, +# `_ExclusivePeriodicAxis <: AbstractVector` sends it into the raw-Vector `:exclusive` arms +# (`cached_vector.jl`), nesting toward length (n+1)+1 and throwing in the ctor. Mirrors the +# `_cache_axis` passthroughs below; non-exclusive 2-/3-arg forms already fall to the raw `= x` arm. +@inline _resolve_axis(g::_ExclusivePeriodicAxis, ::PeriodicBC{:exclusive}) = g +@inline _resolve_axis(g::_ExclusivePeriodicAxis, ::PeriodicBC{:exclusive}, ::Type{Tg}) where {Tg} = g @inline _cache_axis(g::_ExclusivePeriodicAxis) = g @inline _cache_axis(g::_ExclusivePeriodicAxis, ::AbstractBC) = g @inline _cache_axis(g::_ExclusivePeriodicAxis, ::PeriodicBC{:exclusive}) = g diff --git a/src/core/utils.jl b/src/core/utils.jl index 6dfc34e82..f428d7700 100644 --- a/src/core/utils.jl +++ b/src/core/utils.jl @@ -188,6 +188,12 @@ end # coefficients (cubic `z`, quadratic `a`/`d`, hermite `dy`) are solved before any query. @inline _coeff_op(h::Tg, yv::Tv) where {Tg, Tv} = yv + yv * inv(h) +# `_inv_op` (1-arg): reciprocal-spacing eltype — `inv(h)` for an axis already at the +# value-matched width (compose: `_promote_eltype(_inv_op, _promote_grid_float(Tg, Tv))`). +# Floats Int (`inv(Int)::Float64` never survives a narrow value space upstream), keeps +# Unitful inverse units and duck carriers. +@inline _inv_op(h) = inv(h) + # `_integrate_op` (3-arg): the definite-integral element type — value × spacing. # ∫ ≈ Σ yᵢ·hᵢ is dimensionally distinct from the eval witnesses (which weight the value # by the dimensionless offset `dL/h`). `span` is the integration length (`b2 - xL` for a @@ -217,20 +223,28 @@ end # `_CachedVector`/`_CachedRange` axis uses the cached reciprocal (no division) and a # `_UnitStep` range folds the multiply to identity. Raw `AbstractVector` computes # `inv(h)` on the fly (≤1 ULP vs `/h`, perf-neutral: the extra multiply runs free in -# the divider's shadow). `Tc` matches every call site. -@inline function _forward_secant(x, y, i) - Tc = _promote_eltype(_coeff_op, eltype(x), eltype(y)) - return @inbounds _fielddiff(Tc, y[i + 1], y[i]) * _get_inv_h(x, i) +# the divider's shadow). +# +# Width-first form: `Tw` is the value-matched coordinate width from the caller's +# surface (`_promote_grid_float(Tg, Tv)`) — the reciprocal is born at `Tw`, so a raw +# Int axis stops minting `inv(Int)::Float64` beside narrower data. The width-less +# forms delegate with `Tw = eltype(x)`: bit-identical to the historic raw behavior. +@inline function _forward_secant(::Type{Tw}, x, y, i) where {Tw} + Tc = _promote_eltype(_coeff_op, Tw, eltype(y)) + return @inbounds _fielddiff(Tc, y[i + 1], y[i]) * _get_inv_h(Tw, x, i) end +@inline _forward_secant(x, y, i) = _forward_secant(eltype(x), x, y, i) # Backward secant at i is the forward secant of the previous cell (denominator h_{i-1}). -@inline _backward_secant(x, y, i) = _forward_secant(x, y, i - 1) +@inline _backward_secant(::Type{Tw}, x, y, i) where {Tw} = _forward_secant(Tw, x, y, i - 1) +@inline _backward_secant(x, y, i) = _forward_secant(eltype(x), x, y, i - 1) # Centered (2-cell-span) secant (y[i+1]-y[i-1]) / (x[i+1]-x[i-1]) via `_get_inv_2cell`. -@inline function _centered_secant(x, y, i) - Tc = _promote_eltype(_coeff_op, eltype(x), eltype(y)) - return @inbounds _fielddiff(Tc, y[i + 1], y[i - 1]) * _get_inv_2cell(x, i) +@inline function _centered_secant(::Type{Tw}, x, y, i) where {Tw} + Tc = _promote_eltype(_coeff_op, Tw, eltype(y)) + return @inbounds _fielddiff(Tc, y[i + 1], y[i - 1]) * _get_inv_2cell(Tw, x, i) end +@inline _centered_secant(x, y, i) = _centered_secant(eltype(x), x, y, i) """ _promote_query_eltype(::Type{Tv}, q::Tuple) -> Type diff --git a/src/cubic/cubic_cache_pool.jl b/src/cubic/cubic_cache_pool.jl index 59035c81d..a18144c0a 100644 --- a/src/cubic/cubic_cache_pool.jl +++ b/src/cubic/cubic_cache_pool.jl @@ -691,6 +691,55 @@ end end end +# --------------------------------------------------------------- +# Explicit target float type `Tg` (data-aware bank; one-shot value-match) +# --------------------------------------------------------------- +# The 3-arg forms above derive the cache float from the grid alone +# (`_cache_float_type(eltype(x))`) — an Int grid always lands in the Float64 +# bank. These 4-arg forms take the value-matched float `Tg` (computed at the +# one-shot surface as `_promote_grid_float(grid_eltype, data_eltype)`) so an +# Int grid beside Float32 data is cached at Float32 (`cache.x` eltype = `Tg`). +# Only the Integer/Rational bank routing changes; Float/Range/Dual grids already +# carry their own float and ignore `Tg` (see the impl overloads below). +@inline function _get_cubic_cache( + x::AbstractVector, + bc::BCPair{L, R}, + autocache::Bool, + ::Type{Tg} + ) where {L <: PointBC, R <: PointBC, Tg} + bc_cache = _cache_bc_pair(bc, Tg) + x_norm = _resolve_axis(x) + if autocache + return _get_derivative_cache_impl(x_norm, bc_cache, Tg) + else + return _build_derivative_bc_cache(_to_float(x_norm, Tg), bc_cache.left, bc_cache.right) + end +end + +@inline function _get_cubic_cache( + x::AbstractVector, + bc::AbstractBC, + autocache::Bool, + ::Type{Tg} + ) where {Tg} + x_norm = _resolve_axis(x) + return if autocache + if bc isa PeriodicBC + return _get_periodic_cache_impl(x_norm, bc, Tg) + else + return _get_derivative_cache_impl(x_norm, _normalize_bc(bc), Tg) + end + else + x_typed = _to_float(x_norm, Tg) + if _is_periodic_bc(bc) + return _build_periodic_cache(x_typed, bc) + else + bc_normalized = _normalize_bc(bc) + return _build_derivative_bc_cache(x_typed, bc_normalized.left, bc_normalized.right) + end + end +end + # =============================================================== # Internal: Cache Implementation @@ -743,6 +792,22 @@ end return _build_derivative_bc_cache(_to_float(x, FT), bc_pair.left, bc_pair.right) end +# ---- Explicit target float `Tg` (data-aware bank) ---- +# Integer/Rational grids route to the `Vector{Tg}` bank (value-matched) instead +# of the blind `Vector{float(T)}`. All other grids (Float/Range/Dual) already +# carry their own float or are ephemeral, so `Tg` is redundant — delegate to the +# default 3-arg routing. +@inline function _get_derivative_cache_impl(x::AbstractVector{T}, bc_pair::BCPair{L, R}, ::Type{Tg}) where {T <: Integer, L <: PointBC, R <: PointBC, Tg} + bank = _get_derivative_bank(Vector{Tg}, bc_pair) + return _lookup_or_insert!(bank, x, bc_pair) +end +@inline function _get_derivative_cache_impl(x::AbstractVector{T}, bc_pair::BCPair{L, R}, ::Type{Tg}) where {T <: Rational, L <: PointBC, R <: PointBC, Tg} + bank = _get_derivative_bank(Vector{Tg}, bc_pair) + return _lookup_or_insert!(bank, x, bc_pair) +end +@inline _get_derivative_cache_impl(x::AbstractVector, bc_pair::BCPair{L, R}, ::Type{Tg}) where {L <: PointBC, R <: PointBC, Tg} = + _get_derivative_cache_impl(x, bc_pair) + """ Internal implementation for periodic BC cache lookup. `bc` is threaded through so `_get_periodic_bank` selects the right E-variant bank — `:inclusive`, @@ -781,3 +846,17 @@ end @inline function _get_periodic_cache_impl(x::AbstractVector, bc::PeriodicBC) return _build_periodic_cache(_to_float(x, _cache_float_type(eltype(x))), bc) end + +# ---- Explicit target float `Tg` (data-aware bank) ---- +# Same contract as the derivative impl: Integer/Rational route to `Vector{Tg}`; +# other grids ignore `Tg` and delegate to the default routing. +@inline function _get_periodic_cache_impl(x::AbstractVector{T}, bc::PeriodicBC, ::Type{Tg}) where {T <: Integer, Tg} + bank = _get_periodic_bank(Vector{Tg}, bc) + return _lookup_or_insert!(bank, x, bc) +end +@inline function _get_periodic_cache_impl(x::AbstractVector{T}, bc::PeriodicBC, ::Type{Tg}) where {T <: Rational, Tg} + bank = _get_periodic_bank(Vector{Tg}, bc) + return _lookup_or_insert!(bank, x, bc) +end +@inline _get_periodic_cache_impl(x::AbstractVector, bc::PeriodicBC, ::Type{Tg}) where {Tg} = + _get_periodic_cache_impl(x, bc) diff --git a/src/cubic/cubic_oneshot.jl b/src/cubic/cubic_oneshot.jl index f653806fe..a4a811b20 100644 --- a/src/cubic/cubic_oneshot.jl +++ b/src/cubic/cubic_oneshot.jl @@ -85,8 +85,11 @@ Type-Free design: handles both concrete (Deriv1{T}) and lazy (PolyFit{D}) types. @assert length(y) == length(x) "y length must match x" @assert length(output) == length(x_query) "output length must match x_query" - # Cache uses structural equivalent (PolyFit → Deriv1 via _cache_bc_pair internally) - cache = _get_cubic_cache(x, bc, _effective_autocache(autocache, Tg)) + # Cache uses structural equivalent (PolyFit → Deriv1 via _cache_bc_pair internally). + # Value-matched `Tg_eff` (Int grid + Float32 data → Float32) selects the data-aware + # cache bank so `cache.x` — and thus the solve — is at the value width. + Tg_eff = _promote_grid_float(Tg, eltype(y)) + cache = _get_cubic_cache(x, bc, _effective_autocache(autocache, Tg), Tg_eff) Tz = _promote_eltype(_coeff_op, eltype(cache.x), eltype(y)) z = acquire!(pool, Tz, length(y)) # Solve uses original BC for proper RHS materialization @@ -116,8 +119,11 @@ AD-compatible: xq is unconstrained to support ForwardDiff.Dual types. op::O, searcher::S ) where {Tg, Tv, Tq <: Real, L <: PointBC, R <: PointBC, O <: AbstractEvalOp, S <: Searcher} - # Cache uses structural equivalent (PolyFit → Deriv1 via _cache_bc_pair internally) - cache = _get_cubic_cache(x, bc, _effective_autocache(autocache, Tg)) + # Cache uses structural equivalent (PolyFit → Deriv1 via _cache_bc_pair internally). + # Value-matched `Tg_eff` selects the data-aware cache bank (Int grid + Float32 + # data → Float32 cache), so scalar ≡ batch ≡ persistent at the value width. + Tg_eff = _promote_grid_float(Tg, Tv) + cache = _get_cubic_cache(x, bc, _effective_autocache(autocache, Tg), Tg_eff) Tz = _promote_eltype(_coeff_op, eltype(cache.x), Tv) tmp_z = acquire!(pool, Tz, length(y)) # Solve uses original BC for proper RHS materialization @@ -157,7 +163,10 @@ lifetime). `y_eff` returned for caller convenience — same object as `y` # Build cache on the user's grid (BC-aware: `:inclusive` user length n+1 OR # `:exclusive` user length n → wrapped axis virtual n+1). Zero-copy. - cache = _get_cubic_cache(x, bc, _effective_autocache(autocache, Tg)) + # Value-matched `Tg_eff` selects the data-aware cache bank (Int grid + Float32 + # data → Float32 cache). + Tg_eff = _promote_grid_float(Tg, Tv) + cache = _get_cubic_cache(x, bc, _effective_autocache(autocache, Tg), Tg_eff) # `_resolve_data` handles the per-bc endpoint validation (`:inclusive` # checks `y[1] ≈ y[end]`; `:exclusive` is a no-op wrap to length n+1). y_eff = _resolve_data(y, bc) @@ -238,7 +247,10 @@ In-place cubic spline interpolation with optional automatic caching. deriv::DerivOp = EvalValue(), search::AbstractSearchPolicy = AutoSearch() ) where {Tg, Tv} - x = _resolve_axis(x) + # Value-matched Tg: Int/OneTo grid + Float32 data → Float32 axis, so the spline + # cache builds (and memoises — `_CachedRange` is isbits, objectid-deterministic) + # at the value width instead of the blind Float64. + x = _resolve_axis(x, _promote_grid_float(Tg, Tv)) # No BC on Searcher: seam handled by axis-level dispatch on `cache.x` at eval. searcher = _resolve_search(x, x_query, search, nothing) # Periodic BC @@ -353,7 +365,9 @@ function cubic_interp( search::AbstractSearchPolicy = AutoSearch(), hint::Union{Nothing, Base.RefValue{Int}} = nothing ) where {Tg, Tv, Tq <: Real} - x = _resolve_axis(x) + # Value-matched Tg (see the in-place form above): Ranges resolve to the value + # width; raw Vectors pass through (identity-keyed cache — legacy width there). + x = _resolve_axis(x, _promote_grid_float(Tg, Tv)) # No BC on Searcher: seam handled by axis-level dispatch on `cache.x` at eval. searcher = _resolve_search(x, xq, search, hint) if _is_periodic_bc(bc) diff --git a/src/cubic/nd/cubic_nd_oneshot.jl b/src/cubic/nd/cubic_nd_oneshot.jl index c23a229eb..b8378a42d 100644 --- a/src/cubic/nd/cubic_nd_oneshot.jl +++ b/src/cubic/nd/cubic_nd_oneshot.jl @@ -37,12 +37,13 @@ function cubic_interp( hint::Union{Nothing, NTuple{N, Base.RefValue{Int}}} = nothing ) where {Tv, N} # Scalar one-shot: raw grids — a stable grid id lets `_get_cubic_cache` memoise - # (a per-call `Tg.(x)` copy would miss every time + alloc). Output types from - # op-shape inference (`inv(h)`/`dL/h` float Int). Batch keeps eager-convert. - Tg = _promote_grid_eltype(grids) + # (a per-call copy would miss every time + alloc). `Tg` is value-matched (Int/OneTo grid + + # Float32 data → Float32), so the OnTheFly eval + witness `Tr` agree. Batch keeps eager-convert. + Tg = _promote_grid_float(_promote_grid_eltype(grids), Tv) Tv_p = _promote_eltype(_coeff_op, Tg, Tv) _validate_nd_grids(grids, data) - Tr = _promote_eltype(_interp_op, Tg, Tv, promote_type(typeof.(query)...)) + Tq = promote_type(typeof.(query)...) + Tr = _promote_eltype(_interp_op, Tg, Tv, Tq) bcs = _resolve_bcs_nd(bc, Val(N)) searches = _resolve_search_nd(search, Val(N), query) # NTuple{N,Real} <: Tuple → BinarySearch/axis @@ -120,9 +121,12 @@ Zero-allocation after warmup (pool reuse). ops::NTuple{N, AbstractEvalOp}, hints = nothing ) where {Tv, N} - # Raw grids: per-axis partials + `_compute_all_local_params` (promotes the cell - # width) accept a raw/heterogeneous axis. `Tg` only types the pooled buffer. - Tg = _promote_grid_eltype(grids) + # Value-matched pooled wrap, symmetric with the quadratic scalar backend: Ranges + # → isbits `_CachedRange{Tg}` (the per-axis spline caches still memoise via + # value-deterministic objectid); a mismatched Vector converts into a POOL buffer + # (warm one-shots stay zero-alloc), so the whole solve pipeline runs at `Tg`. + Tg = _promote_grid_float(_promote_grid_eltype(grids), Tv) + grids = _cache_axes_pooled(pool, grids, Tg) # @generated static-Tg unroll (no Type-captured closure) # Bare GridIdx(k).val is NaN → resolve to the grid coordinate for the value kernel (search still uses .idx). query = map(_resolve_grididx, query, grids) # 0. Validate (NoExtrap throw must precede FillExtrap short-circuit) AND promote per axis: @@ -153,10 +157,11 @@ Zero-allocation after warmup (pool reuse). # 4. Eval pipeline (all standalone functions, no Interpolant needed). # Axis-only forms — `grids_p` axes carry `h`/`inv_h` directly via `_get_h`/ - # `_get_inv_h` (cached lookup for wrapped axes, on-the-fly diff for raw Vector). + # `_get_inv_h` (cached lookup for wrapped axes, on-the-fly diff for raw Vector); + # the data-aware form width-types hs/inv_hs at `Tg` (raw Int-Vector axes included). q_evals = _handle_all_extraps(query, grids_p, extraps_eff) indices, Ls, _ = _search_all_intervals(q_evals, grids_p, searches, hints, extraps_eff) - hs, inv_hs, dLs = _compute_all_local_params(q_evals, grids_p, indices, Ls) + hs, inv_hs, dLs = _compute_all_local_params(q_evals, grids_p, indices, Ls, Tg) # 6. Tensor-product kernel evaluation return _eval_nd_cell(partials, indices, hs, inv_hs, dLs, ops) diff --git a/src/cubic/nd/cubic_nd_types.jl b/src/cubic/nd/cubic_nd_types.jl index 1436fde61..b73ea57db 100644 --- a/src/cubic/nd/cubic_nd_types.jl +++ b/src/cubic/nd/cubic_nd_types.jl @@ -105,7 +105,7 @@ struct CubicInterpolantND{ # When the outer API has already wrapped + extended the grids, this # is just a per-axis `copy` of the wrapper (no buffer duplication # for Range-backed wrappers). - grids_c = map((g, bc) -> _convert_copy(_cache_axis(g, bc, Tg), Tg), grids, bcs) + grids_c = map((g, bc, T) -> _convert_copy(_cache_axis(g, bc, T), T), grids, bcs, ntuple(_ -> Tg, Val(N))) return new{Tg, Tv, N, NP1, typeof(grids_c), typeof(bcs), typeof(extraps), typeof(searches)}( grids_c, nodal_derivs, bcs, extraps, searches ) diff --git a/src/hermite/hermite_eval.jl b/src/hermite/hermite_eval.jl index db2988964..c832e0a64 100644 --- a/src/hermite/hermite_eval.jl +++ b/src/hermite/hermite_eval.jl @@ -23,11 +23,16 @@ op::O, searcher::S ) where {Tg, Tv, Tq, O <: AbstractEvalOp, S <: Searcher} + # Compile-time cell-geometry types at the value-matched width (y ∪ dy): the + # converts are no-ops for float/duck/Unitful axes and float an Int axis's + # `inv(h)::Float64` to the value width. dL keeps query blood (Dual partials). + Tw = _hermite_grid_float(Tg, Tv, eltype(dy)) + Tinv = _promote_eltype(_inv_op, Tw) xq = _resolve_grididx(xq, x) idx, idx_R, xL, _ = search_interval(searcher, x, xq, e) dL = xq - xL - h = _get_h(x, idx) - inv_h = _get_inv_h(x, idx) + h = convert(Tw, _get_h(x, idx)) + inv_h = convert(Tinv, _get_inv_h(x, idx)) @inbounds return _hermite_kernel_1d(op, y[idx], y[idx_R], dy[idx], dy[idx_R], h, inv_h, dL) end @@ -46,11 +51,13 @@ end # NoExtrap → InBounds for the search once the domain check passes (lean search). # ExtendExtrap passes through: it may arrive OOB → standard two-sided-clamp search (not the # lean InBounds one, whose one-sided clamp would give idx ≤ 0 OOB-left); boundary cell extrapolates. + Tw = _hermite_grid_float(Tg, Tv, eltype(dy)) + Tinv = _promote_eltype(_inv_op, Tw) extrap_eff = _check_domain(x, xq, extrap) idx, idx_R, xL, _ = search_interval(searcher, x, xq, extrap_eff) dL = xq - xL - h = _get_h(x, idx) - inv_h = _get_inv_h(x, idx) + h = convert(Tw, _get_h(x, idx)) + inv_h = convert(Tinv, _get_inv_h(x, idx)) @inbounds return _hermite_kernel_1d(op, y[idx], y[idx_R], dy[idx], dy[idx_R], h, inv_h, dL) end @@ -140,14 +147,19 @@ end op::O, searcher::S ) where {Tg, Tv, Tq, O <: AbstractEvalOp, S <: Searcher} + # Compile-time value-matched width (value space = y — slopes derive from y): + # local slopes AND the kernel cell geometry all run at `Tw`, so a raw Int + # axis stops minting `inv(Int)::Float64` beside narrower data. dL keeps + # query blood (Dual-query safety). + Tw = _promote_grid_float(Tg, Tv) xq = _resolve_grididx(xq, x) idx, idx_R, xL, _ = search_interval(searcher, x, xq, e) n = _data_length(x) - dyL = _local_slope(sm, x, y, idx, n) - dyR = _local_slope(sm, x, y, idx_R, n) + dyL = _local_slope(Tw, sm, x, y, idx, n) + dyR = _local_slope(Tw, sm, x, y, idx_R, n) dL = xq - xL - h = _get_h(x, idx) - inv_h = _get_inv_h(x, idx) + h = _get_h(Tw, x, idx) + inv_h = _get_inv_h(Tw, x, idx) yr = _raw(y) @inbounds return _hermite_kernel_1d(op, yr[idx], yr[idx_R], dyL, dyR, h, inv_h, dL) end @@ -165,14 +177,16 @@ end # NoExtrap → InBounds for the search once the domain check passes (lean search). # ExtendExtrap passes through: it may arrive OOB → standard two-sided-clamp search (not the # lean InBounds one, whose one-sided clamp would give idx ≤ 0 OOB-left); boundary cell extrapolates. + # Value-matched width for slopes + cell geometry — see the InBounds arm above. + Tw = _promote_grid_float(Tg, Tv) extrap_eff = _check_domain(x, xq, extrap) idx, idx_R, xL, _ = search_interval(searcher, x, xq, extrap_eff) n = _data_length(x) - dyL = _local_slope(sm, x, y, idx, n) - dyR = _local_slope(sm, x, y, idx_R, n) + dyL = _local_slope(Tw, sm, x, y, idx, n) + dyR = _local_slope(Tw, sm, x, y, idx_R, n) dL = xq - xL - h = _get_h(x, idx) - inv_h = _get_inv_h(x, idx) + h = _get_h(Tw, x, idx) + inv_h = _get_inv_h(Tw, x, idx) yr = _raw(y) @inbounds return _hermite_kernel_1d(op, yr[idx], yr[idx_R], dyL, dyR, h, inv_h, dL) end diff --git a/src/hermite/hermite_local_slopes.jl b/src/hermite/hermite_local_slopes.jl index afce92717..357fa1673 100644 --- a/src/hermite/hermite_local_slopes.jl +++ b/src/hermite/hermite_local_slopes.jl @@ -30,7 +30,41 @@ # Boundary (NoBC): _pchip_endpoint_slope (3-point one-sided FD + monotonicity clamping) # Boundary (PeriodicBC): closed-cycle interior formula via abstraction +# Two parallel forms per method: the width-less historic form (integrate / +# adjoint callers — raw-eltype semantics, incl. Rational exactness and the +# float-lifted `:exclusive` seam) and the width-first `(Tw, sm, …)` form the +# eval paths use (hermite_eval.jl sm arms) with the value-matched `Tw`. The +# width-first seam converts period/span AT `Tw`, so it must receive a real +# (float/duck) width — never delegate it a bare Int eltype. @inline function _local_slope(sm::PchipSlopes, x::AbstractVector{Tg}, y::AbstractVector{Tv}, i::Int, n::Int) where {Tg, Tv} + # See the width-first body below for the `_raw` rationale. + xr = _raw(x) + yr = _raw(y) + if n == 2 + if sm.bc isa PeriodicBC + return _pchip_boundary_slope(xr, yr, i, n, sm.bc) + end + return _forward_secant(xr, yr, 1) + end + if i == 1 || i == n + return _pchip_boundary_slope(xr, yr, i, n, sm.bc) + end + @inbounds begin + h_prev = _get_h(xr, i - 1) + h_curr = _get_h(xr, i) + δ_prev = _backward_secant(xr, yr, i) + δ_curr = _forward_secant(xr, yr, i) + end + if sign(δ_prev) != sign(δ_curr) + return zero(_promote_eltype(_coeff_op, Tg, Tv)) + else + w1 = 2 * h_curr + h_prev + w2 = h_curr + 2 * h_prev + return _pchip_harmonic_mean(w1, w2, δ_prev, δ_curr) + end +end + +@inline function _local_slope(::Type{Tw}, sm::PchipSlopes, x::AbstractVector, y::AbstractVector{Tv}, i::Int, n::Int) where {Tw, Tv} # Strip wrappers — all internal accesses are guaranteed in [1, n_raw] for # both NoBC and PeriodicBC paths (boundary helpers route via # `_periodic_secant`/`_periodic_cell_width` which special-case the seam @@ -46,24 +80,26 @@ # diverge from the persistent path's wrap-aware result. if n == 2 if sm.bc isa PeriodicBC - return _pchip_boundary_slope(xr, yr, i, n, sm.bc) + return _pchip_boundary_slope(Tw, xr, yr, i, n, sm.bc) end - return _forward_secant(xr, yr, 1) + return _forward_secant(Tw, xr, yr, 1) end if i == 1 || i == n - return _pchip_boundary_slope(xr, yr, i, n, sm.bc) + return _pchip_boundary_slope(Tw, xr, yr, i, n, sm.bc) end - # Interior: weighted harmonic mean (Fritsch-Carlson) + # Interior: weighted harmonic mean (Fritsch-Carlson). Raw cell widths are + # width-neutral WEIGHTS (Int × Tc-secant stays Tc — no `inv`); only the + # secants carry the value-matched width. @inbounds begin h_prev = _get_h(xr, i - 1) h_curr = _get_h(xr, i) - δ_prev = _backward_secant(xr, yr, i) - δ_curr = _forward_secant(xr, yr, i) + δ_prev = _backward_secant(Tw, xr, yr, i) + δ_curr = _forward_secant(Tw, xr, yr, i) end if sign(δ_prev) != sign(δ_curr) - return zero(_promote_eltype(_coeff_op, Tg, Tv)) + return zero(_promote_eltype(_coeff_op, Tw, Tv)) else w1 = 2 * h_curr + h_prev w2 = h_curr + 2 * h_prev @@ -117,6 +153,42 @@ end end end +# ── Width-first boundary variants — secants (and the `:exclusive` seam width, +# whose resolved period would otherwise carry Float64) run at `Tw`. +@inline function _pchip_boundary_slope(::Type{Tw}, x::AbstractVector, y::AbstractVector, i, n, ::NoBC) where {Tw} + if i == 1 + @inbounds begin + h1 = _get_h(x, 1) + h2 = _get_h(x, 2) + δ1 = _forward_secant(Tw, x, y, 1) + δ2 = _forward_secant(Tw, x, y, 2) + end + return _pchip_endpoint_slope(h1, h2, δ1, δ2) + else # i == n + @inbounds begin + h_last = _get_h(x, n - 1) + h_prev = _get_h(x, n - 2) + δ_last = _backward_secant(Tw, x, y, n) + δ_prev = _backward_secant(Tw, x, y, n - 1) + end + return _pchip_endpoint_slope(h_last, h_prev, δ_last, δ_prev) + end +end + +@inline function _pchip_boundary_slope(::Type{Tw}, x::AbstractVector, y::AbstractVector{Tv}, i, n, bc::PeriodicBC) where {Tw, Tv} + δ_prev = _periodic_secant(Tw, x, y, i - 1, n, bc) + δ_curr = _periodic_secant(Tw, x, y, i, n, bc) + h_prev = _periodic_cell_width(Tw, x, i - 1, n, bc) + h_curr = _periodic_cell_width(Tw, x, i, n, bc) + if sign(δ_prev) != sign(δ_curr) + return zero(_promote_eltype(_coeff_op, Tw, Tv)) + else + w1 = 2 * h_curr + h_prev + w2 = h_curr + 2 * h_prev + return _pchip_harmonic_mean(w1, w2, δ_prev, δ_curr) + end +end + # ======================================== # Cardinal Local Slope # ======================================== @@ -147,6 +219,27 @@ end @inbounds return scale * _centered_secant(xr, yr, i) end +@inline function _local_slope(::Type{Tw}, sm::CardinalSlopes, x::AbstractVector, y::AbstractVector, i::Int, n::Int) where {Tw} + xr = _raw(x) + yr = _raw(y) + # Dimensionless scale at the value-matched width: `float(one(Tw))` is + # Unitful-safe (one(quantity) is dimensionless) and keeps a Float64 default + # `tension` from re-widening narrower computations. + Tsc = typeof(float(one(Tw))) + scale = one(Tsc) - convert(Tsc, sm.tension) + + if n == 2 + if sm.bc isa PeriodicBC + return _cardinal_boundary_slope(Tw, xr, yr, i, n, scale, sm.bc) + end + @inbounds return scale * _forward_secant(Tw, xr, yr, 1) + end + if i == 1 || i == n + return _cardinal_boundary_slope(Tw, xr, yr, i, n, scale, sm.bc) + end + @inbounds return scale * _centered_secant(Tw, xr, yr, i) +end + @inline function _cardinal_boundary_slope(x, y, i, n, scale, ::NoBC) return if i == 1 @inbounds scale * _forward_secant(x, y, 1) @@ -167,6 +260,25 @@ end return scale * (m_prev * h_prev + m_curr * h_curr) / (h_prev + h_curr) end +# ── Width-first boundary variants — secants and the cell widths (this one +# DIVIDES by them, and the `:exclusive` seam width carries the resolved +# period) all run at `Tw`. +@inline function _cardinal_boundary_slope(::Type{Tw}, x, y, i, n, scale, ::NoBC) where {Tw} + return if i == 1 + @inbounds scale * _forward_secant(Tw, x, y, 1) + else + @inbounds scale * _backward_secant(Tw, x, y, n) + end +end + +@inline function _cardinal_boundary_slope(::Type{Tw}, x, y, i, n, scale, bc::PeriodicBC) where {Tw} + h_prev = _periodic_cell_width(Tw, x, i - 1, n, bc) + h_curr = _periodic_cell_width(Tw, x, i, n, bc) + m_prev = _periodic_secant(Tw, x, y, i - 1, n, bc) + m_curr = _periodic_secant(Tw, x, y, i, n, bc) + return scale * (m_prev * h_prev + m_curr * h_curr) / (h_prev + h_curr) +end + # ======================================== # Akima Local Slope # ======================================== @@ -222,6 +334,33 @@ end return _akima_local_4secant(xr, yr, i, n) end +@inline function _local_slope(::Type{Tw}, sm::AkimaSlopes, x::AbstractVector, y::AbstractVector, i::Int, n::Int) where {Tw} + xr = _raw(x) + yr = _raw(y) + if n == 2 + if sm.bc isa PeriodicBC + return _akima_local_4secant_periodic(Tw, xr, yr, i, n, sm.bc) + end + return _forward_secant(Tw, xr, yr, 1) + end + if n == 3 + if sm.bc isa PeriodicBC + return _akima_local_4secant_periodic(Tw, xr, yr, i, n, sm.bc) + end + @inbounds begin + m1 = _forward_secant(Tw, xr, yr, 1) + m2 = _forward_secant(Tw, xr, yr, 2) + end + i == 1 && return m1 + i == 3 && return m2 + return (m1 + m2) / 2 # i == 2 + end + if sm.bc isa PeriodicBC && (i <= 2 || i >= n - 1) + return _akima_local_4secant_periodic(Tw, xr, yr, i, n, sm.bc) + end + return _akima_local_4secant(Tw, xr, yr, i, n) +end + # Akima 4-secant computation with wrap-aware secant access. Dispatches via # `_periodic_secant` so inclusive (m_n = m_1 wrap, mod1 j n-1) and exclusive # (m_n = seam virtual, mod1 j n) both work with one body. @@ -230,6 +369,11 @@ end return _akima_weighted_slope(_ws(i - 2), _ws(i - 1), _ws(i), _ws(i + 1)) end +@inline function _akima_local_4secant_periodic(::Type{Tw}, x, y, i::Int, n::Int, bc::PeriodicBC) where {Tw} + @inline _ws(j) = _periodic_secant(Tw, x, y, j, n, bc) + return _akima_weighted_slope(_ws(i - 2), _ws(i - 1), _ws(i), _ws(i + 1)) +end + # Helper: compute the 4 secants needed for Akima slope at index i (NoBC). # Secant at interval j → j+1 is m_j = (y[j+1] - y[j]) / (x[j+1] - x[j]) # Real secants exist for j = 1..n-1. Virtual secants extend the sequence linearly. @@ -273,3 +417,34 @@ end _safe_secant(i - 2), _safe_secant(i - 1), _safe_secant(i), _safe_secant(i + 1) ) end + +# Width-first NoBC variant — same virtual-secant extrapolation, secants at `Tw`. +@inline function _akima_local_4secant(::Type{Tw}, x::AbstractVector, y::AbstractVector, i::Int, n::Int) where {Tw} + @inline _secant(j) = _forward_secant(Tw, x, y, j) + + nm1 = n - 1 # last valid secant index + + @inline function _safe_secant(j) + if 1 <= j <= nm1 + return _secant(j) + elseif j == 0 + return 2 * _secant(1) - _secant(2) + elseif j == -1 + m1 = _secant(1) + m2 = _secant(2) + return 3 * m1 - 2 * m2 + elseif j == nm1 + 1 # n + return 2 * _secant(nm1) - _secant(nm1 - 1) + elseif j == nm1 + 2 # n+1 + m_last = _secant(nm1) + m_prev = _secant(nm1 - 1) + return 3 * m_last - 2 * m_prev + else + error("_safe_secant: j=$j out of expected range [-1, $(nm1 + 2)] for n=$(nm1 + 1)") + end + end + + return _akima_weighted_slope( + _safe_secant(i - 2), _safe_secant(i - 1), _safe_secant(i), _safe_secant(i + 1) + ) +end diff --git a/src/hermite/hermite_oneshot.jl b/src/hermite/hermite_oneshot.jl index a7f5dc4c1..3935f8f2e 100644 --- a/src/hermite/hermite_oneshot.jl +++ b/src/hermite/hermite_oneshot.jl @@ -30,7 +30,8 @@ C\$^1\$ continuous — slopes are used directly, no global spline solve. search::AbstractSearchPolicy = AutoSearch(), hint::Union{Nothing, Base.RefValue{Int}} = nothing, ) where {Tg, Tv, Tq <: Real} - x = _resolve_axis(x) + # Value space = y ∪ dy: the axis floats against both widths (matching ND partials). + x = _resolve_axis(x, _hermite_grid_float(Tg, Tv, eltype(dy))) @boundscheck length(y) == length(x) || _throw_length_mismatch(length(x), length(y)) @boundscheck length(dy) == length(x) || _throw_length_mismatch(length(x), length(dy), "x", "dy") @boundscheck length(x) >= 2 || throw(ArgumentError("Hermite interpolation requires at least 2 points, got $(length(x))")) @@ -60,7 +61,7 @@ function hermite_interp!( search::AbstractSearchPolicy = AutoSearch(), hint::Union{Nothing, Base.RefValue{Int}} = nothing, ) where {Tg, Tv, Tq <: Real} - x = _resolve_axis(x) + x = _resolve_axis(x, _hermite_grid_float(Tg, Tv, eltype(dy))) @boundscheck length(y) == length(x) || _throw_length_mismatch(length(x), length(y)) @boundscheck length(dy) == length(x) || _throw_length_mismatch(length(x), length(dy), "x", "dy") @boundscheck length(output) == length(x_query) || _throw_length_mismatch(length(x_query), length(output), "x_query", "output") @@ -108,6 +109,13 @@ end # ║ INPUT PROMOTION HELPER ║ # ╚═══════════════════════════════════════════════════════════════════════════╝ +# Hermite grid float: the value space is y ∪ dy, so the axis width must see both +# (x::Float32 + y::Float32 + dy::Float64 → Float64; duck dy leaves the grid alone). +@inline function _hermite_grid_float(::Type{TX}, ::Type{TY}, ::Type{TDY}) where {TX, TY, TDY} + Tg_y = _promote_grid_float(TX, TY) + return TDY <: _PromotableValue ? promote_type(Tg_y, float(_real_eltype(TDY))) : Tg_y +end + # Joint promotion of (x, y, dy) — grid type Tg considers all three inputs, # so e.g. x::Float32 + y::Float32 + dy::Float64 → Tg=Float64 (no precision loss). @inline function _promote_hermite_inputs( @@ -115,8 +123,7 @@ end y::AbstractVector{TY}, dy::AbstractVector{TDY}, ) where {TX, TY, TDY} - Tg_y = _promote_grid_float(TX, TY) - Tg = TDY <: _PromotableValue ? promote_type(Tg_y, float(_real_eltype(TDY))) : Tg_y + Tg = _hermite_grid_float(TX, TY, TDY) x_p = _to_float(x, Tg) # Only promote values when Tg is a standard float type — duck-typed Tg (Dual etc.) # leaves y/dy as-is; Julia arithmetic promotion handles the rest in kernels. diff --git a/src/hermite/hermite_periodic_slopes.jl b/src/hermite/hermite_periodic_slopes.jl index 0a5795162..3a9e18ebd 100644 --- a/src/hermite/hermite_periodic_slopes.jl +++ b/src/hermite/hermite_periodic_slopes.jl @@ -75,6 +75,19 @@ end return _forward_secant(x, y, jw) end +# ── Width-first forms — `Tw` = value-matched coordinate width (see utils.jl +# secants). Real-cell secants/widths thread `Tw` straight through; the +# `:exclusive` seam resolves its period AT `Tw` and differences span-first, so +# neither the F64 period lift nor `inv(Int)` can widen narrower data. +@inline function _periodic_secant(::Type{Tw}, x::AbstractVector, y::AbstractVector, j::Int, n::Int, ::PeriodicBC{:inclusive}) where {Tw} + jw = mod1(j, n - 1) + return _forward_secant(Tw, x, y, jw) +end +@inline function _periodic_secant(::Type{Tw}, x::AbstractVector, y::AbstractVector, j::Int, n::Int, ::PeriodicBC{:extended}) where {Tw} + jw = mod1(j, n - 1) + return _forward_secant(Tw, x, y, jw) +end + # Cast the resolved exclusive period to the grid's promoted-float type so # seam-cell arithmetic (`seam_h`, seam secant) stays in the grid eltype. # Without this cast, a `Float64` user-supplied period combined with a @@ -100,6 +113,16 @@ end return _forward_secant(x, y, jw) end +@inline function _periodic_secant(::Type{Tw}, x::AbstractVector, y::AbstractVector, j::Int, n::Int, bc::PeriodicBC{:exclusive}) where {Tw} + jw = mod1(j, n) + if jw == n + seam_h = _periodic_cell_width(Tw, x, n, n, bc) + Tc = _promote_eltype(_coeff_op, Tw, eltype(y)) + @inbounds return _fielddiff(Tc, y[1], y[n]) / seam_h + end + return _forward_secant(Tw, x, y, jw) +end + """ _periodic_cell_width(x, j, n, bc) -> Tg @@ -129,5 +152,22 @@ end return _get_h(x, jw) end +# Width-first cell widths: real cells convert the span (exact for Int axes); the +# `:exclusive` seam resolves period and span both AT `Tw`, so the seam width — +# which feeds PCHIP's harmonic-mean weights and Cardinal's boundary divisor — +# cannot re-widen a value-matched computation. +@inline _periodic_cell_width(::Type{Tw}, x::AbstractVector, j::Int, n::Int, ::PeriodicBC{:inclusive}) where {Tw} = + _get_h(Tw, x, mod1(j, n - 1)) +@inline _periodic_cell_width(::Type{Tw}, x::AbstractVector, j::Int, n::Int, ::PeriodicBC{:extended}) where {Tw} = + _get_h(Tw, x, mod1(j, n - 1)) +@inline function _periodic_cell_width(::Type{Tw}, x::AbstractVector, j::Int, n::Int, bc::PeriodicBC{:exclusive}) where {Tw} + jw = mod1(j, n) + if jw == n + period = convert(Tw, _resolve_exclusive_period(x, bc)) + return period - convert(Tw, @inbounds x[n] - x[1]) + end + return _get_h(Tw, x, jw) +end + # `_bc_after_extend` lives in `src/core/periodic.jl` next to `_periodic_extend_1d` # — it is generic post-extension BC normalization, not Hermite-specific. diff --git a/src/hermite/nd/hermite_nd_build.jl b/src/hermite/nd/hermite_nd_build.jl index c77beee76..48dca0308 100644 --- a/src/hermite/nd/hermite_nd_build.jl +++ b/src/hermite/nd/hermite_nd_build.jl @@ -317,11 +317,11 @@ end data::AbstractArray{Tv, N}, partials::HermitePartials{N, Tv, K}, bcs::Tuple{Vararg{AbstractBC, N}}, - ) where {Tv, N, K} - # `Tg` only feeds the `:exclusive` virtual-endpoint extension below; raw grids - # may be Int/heterogeneous, so float it — a non-integer period on an Int grid - # must not be forced through `Int` (mirrors the cubic/quad twin in periodic.jl). - Tg = float(_promote_grid_eltype(grids)) + # `Tg` only feeds the `:exclusive` virtual-endpoint extension below — value-matched + # by the callers (Int grid + Float32 data → Float32) so an extended axis cannot + # reintroduce Float64 into `dL = q - xL`; still floats a non-integer period. + ::Type{Tg} = float(_promote_grid_eltype(grids)), + ) where {Tv, N, K, Tg} extended = ntuple(d -> bcs[d] isa PeriodicBC{:exclusive}, Val(N)) n_orig = size(data) n_ext = ntuple(d -> extended[d] ? n_orig[d] + 1 : n_orig[d], Val(N)) diff --git a/src/hermite/nd/hermite_nd_interpolant.jl b/src/hermite/nd/hermite_nd_interpolant.jl index 87e7c1f91..6357f4d54 100644 --- a/src/hermite/nd/hermite_nd_interpolant.jl +++ b/src/hermite/nd/hermite_nd_interpolant.jl @@ -37,9 +37,9 @@ function CubicHermiteInterpolantND( # already enforces this, so this guards direct struct construction. K == (1 << N) - 1 || _throw_partials_not_full_mixed(N, K) - # Promote across (grid, data, partials) to a single Tv. - grids_typed, _, Tv_promoted, _ = _nd_promote_grids(grids, data) - Tv = promote_type(Tv_promoted, Tv_part) + # Promote across (grid, data, partials) to a single Tv — the grid value-match must + # see the partials' width too (value space = data ∪ partials, matching one-shot). + grids_typed, _, Tv, _ = _nd_promote_grids(grids, data, Tv_part) data_typed = _coerce_data_eltype(data, Tv, Val(N)) partials_typed = _coerce_partials_eltype(partials, Tv, Val(N)) diff --git a/src/hermite/nd/hermite_nd_oneshot.jl b/src/hermite/nd/hermite_nd_oneshot.jl index 2b2932b55..6d4c59cf2 100644 --- a/src/hermite/nd/hermite_nd_oneshot.jl +++ b/src/hermite/nd/hermite_nd_oneshot.jl @@ -61,7 +61,8 @@ function hermite_interp( search::Union{AbstractSearchPolicy, NTuple{N, AbstractSearchPolicy}} = AutoSearch(), hint::Union{Nothing, NTuple{N, Base.RefValue{Int}}} = nothing, ) where {Tv, N, Tv_part} - Tg = _promote_grid_eltype(grids) + # `Tg` value-matched as in `_hermite_oneshot_prepare`, so `Tr` agrees with the eval. + Tg = _promote_grid_float(_promote_grid_eltype(grids), promote_type(Tv, Tv_part)) Tq = _query_eltype(queries) Tr = _promote_eltype(_interp_op, Tg, promote_type(Tv, Tv_part), Tq) output = Vector{Tr}(undef, _query_length(queries)) @@ -112,7 +113,9 @@ end K == (1 << N) - 1 || _throw_partials_not_full_mixed(N, K) # Raw grids: the pack + cell-eval float the cell width, so no eager `Tg.(x)` copy. - Tg = _promote_grid_eltype(grids) + # `Tg` value-matched to data∪partials (Int grid + Float32 → Float32) keeps `Tv` narrow, + # so `_coerce_*_eltype` pass through instead of copying data + K partials per call. + Tg = _promote_grid_float(_promote_grid_eltype(grids), promote_type(eltype(data), Tv_part)) Tv_promoted = _promote_eltype(_coeff_op, Tg, eltype(data)) Tv = promote_type(Tv_promoted, Tv_part) data_typed = _coerce_data_eltype(data, Tv, Val(N)) @@ -158,8 +161,10 @@ end oob_result !== nothing && return oob_result # `bcs_post` (3rd return) is unused: no partial-solve step here, and - # extrap is materialized from `extraps_val` + `grids_p`. - grids_p, buf, _ = _pack_and_extend_nodal_derivs_pooled(pool, grids, data, partials, bcs) + # extrap is materialized from `extraps_val` + `grids_p`. `Tg` value-matches the + # width/extension type (Int grid + Float32 data → Float32) — grids stay raw. + Tg = _promote_grid_float(_promote_grid_eltype(grids), Tv) + grids_p, buf, _ = _pack_and_extend_nodal_derivs_pooled(pool, grids, data, partials, bcs, Tg) extraps_eff = map(_resolve_extrap, extraps_val, grids_p) # `mono` is all-true for a single point (no monotone-hint optimization). @@ -168,7 +173,7 @@ end q_evals = _handle_all_extraps(query, grids_p, extraps_eff) indices, Ls, _ = _search_all_intervals(q_evals, grids_p, searches, hints, mono, extraps_eff) - hs, inv_hs, dLs = _compute_all_local_params(q_evals, grids_p, indices, Ls) + hs, inv_hs, dLs = _compute_all_local_params(q_evals, grids_p, indices, Ls, Tg) return _eval_nd_cell(buf, indices, hs, inv_hs, dLs, ops) end @@ -196,7 +201,8 @@ end _query_validate(queries) # `bcs_post` (3rd return) is unused here — see the scalar path note. - grids_p, buf, _ = _pack_and_extend_nodal_derivs_pooled(pool, grids, data, partials, bcs) + Tg = _promote_grid_float(_promote_grid_eltype(grids), Tv) + grids_p, buf, _ = _pack_and_extend_nodal_derivs_pooled(pool, grids, data, partials, bcs, Tg) extraps_eff = map(_resolve_extrap, extraps_val, grids_p) extraps_eff = _validate_nd_domain(grids_p, queries, extraps_eff) @@ -209,7 +215,7 @@ end end q_evals = _handle_all_extraps(query_k, grids_p, extraps_eff) indices, Ls, _ = _search_all_intervals(q_evals, grids_p, policies, hints, extraps_eff) - hs, inv_hs, dLs = _compute_all_local_params(q_evals, grids_p, indices, Ls) + hs, inv_hs, dLs = _compute_all_local_params(q_evals, grids_p, indices, Ls, Tg) output[k] = _eval_nd_cell(buf, indices, hs, inv_hs, dLs, ops) end return output diff --git a/src/hermite/nd/hermite_nd_types.jl b/src/hermite/nd/hermite_nd_types.jl index 20d8fba44..3eb663b8d 100644 --- a/src/hermite/nd/hermite_nd_types.jl +++ b/src/hermite/nd/hermite_nd_types.jl @@ -108,7 +108,7 @@ struct CubicHermiteInterpolantND{ searches::Tuple{Vararg{AbstractSearchPolicy, N}}, ) where {Tg, Tv, N, NP1} NP1 == N + 1 || throw(ArgumentError("NP1 must equal N+1")) - grids_c = map((g, bc) -> _convert_copy(_cache_axis(g, bc, Tg), Tg), grids, bcs) + grids_c = map((g, bc, T) -> _convert_copy(_cache_axis(g, bc, T), T), grids, bcs, ntuple(_ -> Tg, Val(N))) return new{Tg, Tv, N, NP1, typeof(grids_c), typeof(bcs), typeof(extraps), typeof(searches)}( grids_c, nodal_derivs, bcs, extraps, searches, ) diff --git a/src/hetero/hetero_oneshot.jl b/src/hetero/hetero_oneshot.jl index fa4755c5e..490b19381 100644 --- a/src/hetero/hetero_oneshot.jl +++ b/src/hetero/hetero_oneshot.jl @@ -146,7 +146,12 @@ end # `_wrap_to_domain` / `search_interval` without the raw n-length Vector's # `last - first ≠ period` mismatch (Linear/Constant ND mirror this pattern). bcs = map(_bc_for_periodic_check, methods) - grids_eff = map(_resolve_axis, grids, bcs) + # Value-matched grid float (Int grid + Float32 data → Float32) — output matches the caller's witness. + Tg = _promote_grid_float(_promote_grid_eltype(grids), Tv) + # @generated static-Tg unroll: a Type captured in a closure (or an + # `ntuple(_ -> Tg, …)` element) de-optimizes under weak const-prop — + # LTS per-fiber heap, and per-axis dynamic dispatch on 1.12 CI workers. + grids_eff = _resolve_axes(grids, bcs, Tg) # NOTE: inclusive PeriodicBC slice validation is NOT performed here — it is # hoisted to the callers (`_interp_nd_oneshot_dispatch` and the OnTheFly # branch of `_interp_nd_oneshot_batch_dispatch!`) so the batch path pays the @@ -154,8 +159,7 @@ end extraps_eff = map(_resolve_extrap, extraps_val, bcs, grids_eff) q_eval = _handle_all_extraps(query, grids_eff, extraps_eff) # Tr promotes data with grid + query eltypes → Dual-safe pool buffers for AD. - # Raw grids may be Int/heterogeneous; `_promote_eltype` floats Int internally. - Tg = _promote_grid_eltype(grids) + # `Tg` (value-matched above) already floats Int; reuse it here. Tr = _promote_eltype(Tv, Tg, typeof.(q_eval)...) # GridIdx safety gate: same reason as the persistent path — a GridIdx on a @@ -197,7 +201,11 @@ end ) end - # Pure global-solve path: no pre-search, full windows, bit-for-bit pre-Phase-3 behavior. + # Global-solve path: grids straight through. The caller + # (`_interp_nd_oneshot_dispatch`) promotes types only — the axes arrive raw, + # and each inner 1D one-shot value-matches its own axis via the data-aware + # cache / width-first geometry (and wraps `:exclusive` axes itself — user + # length n, not the wrapped virtual n+1). `_collapse_dims` emits `Tr` directly. full_windows = map(Base.OneTo, size(data)) return _collapse_dims(Tr, data, grids, methods, extraps_eff, q_eval, ops, searches, hints, full_windows) end @@ -252,8 +260,13 @@ function _interp_nd_oneshot_dispatch( methods::Tuple{Vararg{AbstractInterpMethod, N}}, deriv, extrap, search, hints, coeffs, ) where {N} - grids_typed, Tg, Tv, _ = _nd_promote_grids(grids, data) - _validate_nd_grids(grids_typed, data) + # Type-only promotion: no eager grid convert. The OnTheFly path value-matches + # each raw axis at its own inner 1D surface (data-aware caches / width-first + # geometry), so it takes the raw `grids` straight through. The PreCompute path + # (`_interp_nd_hetero_oneshot`) requires homogeneous-eltype `Tg` grids and uses + # them directly, so it still gets the converted grids. + Tg, Tv, _ = _nd_promote_types(grids, data) + _validate_nd_grids(grids, data) Tr = _promote_eltype(eltype(data), Tg, typeof.(query)...) # bc-aware extrap: NoExtrap → WrapExtrap on PeriodicBC axes. @@ -266,12 +279,12 @@ function _interp_nd_oneshot_dispatch( extraps_val = _resolve_extrap(extrap, bcs, Val(N), Tv) searches = _resolve_search_nd(search, Val(N), query) ops = _resolve_deriv_nd(deriv, Val(N)) - _validate_axis_methods(grids_typed, methods, extraps_val) + _validate_axis_methods(grids, methods, extraps_val) if coeffs isa OnTheFly - return _interp_nd_oneshot_onthefly(grids_typed, data, query, methods, extraps_val, searches, ops, hints)::Tr + return _interp_nd_oneshot_onthefly(grids, data, query, methods, extraps_val, searches, ops, hints)::Tr end - return _interp_nd_hetero_oneshot(grids_typed, data, query, methods, extraps_val, searches, ops, hints)::Tr + return _interp_nd_hetero_oneshot(_convert_grids_typed(grids, Tg), data, query, methods, extraps_val, searches, ops, hints)::Tr end # ======================================== @@ -326,15 +339,17 @@ end methods::Tuple{Vararg{AbstractInterpMethod, N}}, deriv, extrap, search, hints, coeffs, ) where {N} - grids_typed, _, Tv, _ = _nd_promote_grids(grids, data) - _validate_nd_grids(grids_typed, data) + # Type-only promotion (see the scalar dispatch): OnTheFly takes raw `grids` + # and value-matches per axis; PreCompute gets the converted grids. + Tg, Tv, _ = _nd_promote_types(grids, data) + _validate_nd_grids(grids, data) _query_check_ndims(queries, Val(N)) # bc-aware extrap (matches scalar dispatch). bcs = map(_bc_for_periodic_check, methods) extraps_val = _resolve_extrap(extrap, bcs, Val(N), Tv) ops = _resolve_deriv_nd(deriv, Val(N)) - _validate_axis_methods(grids_typed, methods, extraps_val) + _validate_axis_methods(grids, methods, extraps_val) if coeffs isa OnTheFly nq = _query_length(queries) @@ -352,12 +367,12 @@ end searches = _resolve_search_nd(search, Val(N)) @inbounds for k in 1:nq query_k = _extract_query_point(queries, k, Val(N)) - output[k] = _interp_nd_oneshot_onthefly(grids_typed, data, query_k, methods, extraps_val, searches, ops, hints) + output[k] = _interp_nd_oneshot_onthefly(grids, data, query_k, methods, extraps_val, searches, ops, hints) end return output end - return _interp_nd_hetero_batch_dispatch!(output, grids_typed, data, queries, methods, extraps_val, ops, search, hints) + return _interp_nd_hetero_batch_dispatch!(output, _convert_grids_typed(grids, Tg), data, queries, methods, extraps_val, ops, search, hints) end # ======================================== diff --git a/src/hetero/hetero_types.jl b/src/hetero/hetero_types.jl index 928748ca9..4be476dde 100644 --- a/src/hetero/hetero_types.jl +++ b/src/hetero/hetero_types.jl @@ -82,7 +82,7 @@ struct HeteroInterpolantND{ store::StorePolicy = StorePolicy() ) where {Tg, N} Tv = eltype(data) - grids_c = map((g, bc, m) -> _own_or_ref_axis(_cache_axis_for_method(g, bc, Tg, m), Tg, store), grids, bcs, methods) + grids_c = map((g, bc, m, T) -> _own_or_ref_axis(_cache_axis_for_method(g, bc, T, m), T, store), grids, bcs, methods, ntuple(_ -> Tg, Val(N))) return new{Tg, Tv, N, typeof(grids_c), typeof(methods), typeof(extraps), typeof(searches), typeof(data)}( grids_c, data, methods, extraps, searches ) diff --git a/src/linear/linear_oneshot.jl b/src/linear/linear_oneshot.jl index 80e0a05cb..cda4c2e94 100644 --- a/src/linear/linear_oneshot.jl +++ b/src/linear/linear_oneshot.jl @@ -77,7 +77,8 @@ function linear_interp!( # for `:exclusive`). # BC info lives in the axis type after resolution → searcher uses `NoBC()`, # the seam is handled by the wrapper (or by the naturally-extended Range). - x_eff = _resolve_axis(x, bc) + # Value-matched Tg keeps the loop interior bit-consistent with the scalar path. + x_eff = _resolve_axis(x, bc, _promote_grid_float(eltype(x), eltype(y))) y_eff = _resolve_data(y, bc) extrap_eff = _resolve_extrap(extrap, bc, x_eff, y_eff) searcher = _resolve_search(x_eff, x_targets, search, nothing) @@ -204,11 +205,11 @@ For ForwardDiff compatibility, `xq` can be a Dual type: # ======================================== # Core eval: extrap dispatch → search → kernel (no intermediate layers) # ======================================== -# Oneshot path (no spacing): α via direct (q-L)/(R-L) on plain Vector grid -# (`_alpha_of(q, L, R, grid)`); inv_h recomputed per query. -# Persistent path (with spacing): α via cached `inv_h * (q-L)`; mirrors the ND -# `_locate_cell` design — exact at knots is sacrificed for query-time speed -# (1 ULP off possible at right knots; oneshot path remains exact). +# Cell geometry enters at the value-matched width: `inv_h` converts to the +# compile-time `typeof(inv(oneunit(Tw)))` (the `_CachedRange` Tinv idiom) — a no-op +# for float/duck/Unitful axes, and the value-width float for an Int axis whose +# `inv(h::Int)` would otherwise widen everything to Float64. α = (q−L)·inv_h keeps +# query blood (Dual partials) and mirrors the persistent + ND `_locate_cell` form. # Core in-bounds path: search + kernel. All non-InBounds extrap overloads # delegate here after their preprocessing — see cubic_eval.jl for the same @@ -222,16 +223,13 @@ For ForwardDiff compatibility, `xq` can be a Dual type: op::O, searcher::S ) where {Tg, Tv, Tq, O <: AbstractEvalOp, S <: Searcher} + # Compile-time cell-geometry type: reciprocal spacing at the value-matched width. + Tinv = _promote_eltype(_inv_op, _promote_grid_float(Tg, Tv)) xq = _resolve_grididx(xq, x) - idx, idx_R, xL, xR = search_interval(searcher, x, xq, e) - # Independent computation of `α` and `inv_h`. The kernel uses only one - # (EvalValue → α, `DerivOp(1)` → inv_h, `DerivOp(2)` → neither), so the - # unused branch's fdiv is dead-code-eliminated by LLVM. `_alpha_of` - # dispatches on grid type: - # `_CachedRange` → `(q-L) * x.inv_h` (cached fmul, no fdiv) - # raw `AbstractVector` → `(q-L) / float(R-L)` (single fdiv) - α = _alpha_of(xq, xL, xR, x) - @inbounds return _linear_kernel(op, y[idx], y[idx_R], _get_inv_h(x, idx), α) + idx, idx_R, xL, _ = search_interval(searcher, x, xq, e) + inv_h = convert(Tinv, _get_inv_h(x, idx)) + α = _alpha_of(xq, xL, inv_h) + @inbounds return _linear_kernel(op, y[idx], y[idx_R], inv_h, α) end # NoExtrap / ExtendExtrap / others matching AbstractExtrap: domain check (NoExtrap throws @@ -251,10 +249,12 @@ end xq = _resolve_grididx(xq, x) # NoExtrap → InBounds for the search once the domain check passes (lean search); # ExtendExtrap passes through and keeps the two-sided-clamp search (it may arrive OOB). + Tinv = _promote_eltype(_inv_op, _promote_grid_float(Tg, Tv)) extrap_eff = _check_domain(x, xq, extrap) - idx, idx_R, xL, xR = search_interval(searcher, x, xq, extrap_eff) - α = _alpha_of(xq, xL, xR, x) - @inbounds return _linear_kernel(op, y[idx], y[idx_R], _get_inv_h(x, idx), α) + idx, idx_R, xL, _ = search_interval(searcher, x, xq, extrap_eff) + inv_h = convert(Tinv, _get_inv_h(x, idx)) + α = _alpha_of(xq, xL, inv_h) + @inbounds return _linear_kernel(op, y[idx], y[idx_R], inv_h, α) end # ClampExtrap / FillExtrap (all ops): boundary check → extrap value or delegate. @@ -295,11 +295,13 @@ end op::O, searcher::S ) where {Tg, Tv, Tq, O <: AbstractEvalOp, S <: Searcher} + Tinv = _promote_eltype(_inv_op, _promote_grid_float(Tg, Tv)) xq_wrapped = _wrap_to_domain(_resolve_grididx(xq, x), x) - idx, idx_R, xL, xR = search_interval(searcher, x, xq_wrapped) - α = _alpha_of(xq_wrapped, xL, xR, x) + idx, idx_R, xL, _ = search_interval(searcher, x, xq_wrapped) + inv_h = convert(Tinv, _get_inv_h(x, idx)) + α = _alpha_of(xq_wrapped, xL, inv_h) yi = _raw(y) - @inbounds return _linear_kernel(op, yi[idx], yi[idx_R], _get_inv_h(x, idx), α) + @inbounds return _linear_kernel(op, yi[idx], yi[idx_R], inv_h, α) end # Public scalar one-shot API. @@ -320,8 +322,9 @@ end # Same surface-level resolution as the in-place vector form. Zero-alloc: # `_resolve_axis` returns either `x` (Vector passthrough), a stack-allocated # `_CachedRange`, or an `_ExclusivePeriodicAxis` reference wrapper. - # `_resolve_data` is reference-only. - x_eff = _resolve_axis(x, bc) + # `_resolve_data` is reference-only. The Tg is value-matched: an Int/OneTo + # grid beside Float32 data floats to Float32, not the blind Float64. + x_eff = _resolve_axis(x, bc, _promote_grid_float(Tg, Tv)) y_eff = _resolve_data(y, bc) extrap_eff = _resolve_extrap(extrap, bc, x_eff, y_eff) searcher = _resolve_search(x_eff, xq, search, hint) diff --git a/src/linear/nd/linear_nd_adjoint_types.jl b/src/linear/nd/linear_nd_adjoint_types.jl index e5b9a47bc..7e6cd915e 100644 --- a/src/linear/nd/linear_nd_adjoint_types.jl +++ b/src/linear/nd/linear_nd_adjoint_types.jl @@ -96,7 +96,7 @@ struct LinearAdjointND{ Tg, N, B <: NTuple{N, AbstractBC}, EP <: Tuple{Vararg{AbstractExtrap, N}}, } - grids_c = map((g, bc) -> _convert_copy(_cache_axis(g, bc, Tg), Tg), grids, bcs) + grids_c = map((g, bc, T) -> _convert_copy(_cache_axis(g, bc, T), T), grids, bcs, ntuple(_ -> Tg, Val(N))) return new{Tg, N, typeof(grids_c), B, EP}(grids_c, bcs, extraps, anchors, grid_size) end end diff --git a/src/linear/nd/linear_nd_oneshot.jl b/src/linear/nd/linear_nd_oneshot.jl index cc370896c..1ee16c1d8 100644 --- a/src/linear/nd/linear_nd_oneshot.jl +++ b/src/linear/nd/linear_nd_oneshot.jl @@ -42,7 +42,9 @@ function _linear_interp_nd_oneshot( # seam) so kernels read raw `data[..., idx_R, ...]` directly. With BC info # encoded in the axis type, the searcher carries `NoBC()` — the wrapper's # dispatch handles seam regardless. - grids_eff = map(_resolve_axis, grids, bcs) + # Value-matched grid float (Int grid + Float32 data → Float32) — eval matches the `Tr` witness. + Tg = _promote_grid_float(_promote_grid_eltype(grids), Tv) + grids_eff = _resolve_axes(grids, bcs, Tg) # @generated static-Tg unroll (no Type-captured closure) # Bare GridIdx(k).val is NaN → resolve to the grid coordinate for the value kernel (search still uses .idx). query = map(_resolve_grididx, query, grids_eff) # Validate (NoExtrap throw must precede the FillExtrap short-circuit) AND promote per axis: @@ -56,12 +58,14 @@ function _linear_interp_nd_oneshot( extraps_eff = _resolve_extrap(extraps_val, bcs, grids_eff, data, Val(N)) q_eval = _handle_all_extraps(query, grids_eff, extraps_eff) stencils, Ls, Rs = _search_all_intervals_stencil(q_eval, grids_eff, searches, hints, extraps_eff) - αs = map(_alpha_of, q_eval, Ls, Rs, grids_eff) - # 4-arg `_get_inv_h(g, idx, xL, xR)` — `_CachedVector`/`_CachedRange` use - # cached fields (idx-indexed or scalar); raw `Vector` falls back to - # `inv(xR - xL)`. `first(stencil)` = idxL for K=2 (linear) stencils. + # Width-first `_get_inv_h(Tg, g, idx, xL, xR)` — `_CachedVector`/`_CachedRange` + # use cached fields; a raw `Vector` spans `xR - xL` in its own eltype and divides + # at the value-matched `Tg` (span-first: an Int axis would otherwise mint Float64 + # via `inv(Int)`). α derives from the typed inv_h — query-blood promotion preserved + # (Dual query ⇒ Dual α) and the seam cell shares the deriv path's denominator. idxLs = map(first, stencils) - inv_hs = map(_get_inv_h, grids_eff, idxLs, Ls, Rs) + inv_hs = _typed_inv_hs(grids_eff, idxLs, Ls, Rs, Tg) + αs = map(_alpha_of, q_eval, Ls, inv_hs) return _multilinear_sum(data, stencils, inv_hs, αs, ops, Val(N)) end @@ -133,10 +137,9 @@ function linear_interp( deriv::Union{DerivOp, Tuple{Vararg{DerivOp, N}}} = EvalValue(), hint::Union{Nothing, NTuple{N, Base.RefValue{Int}}} = nothing ) where {Tv, N} - # Scalar one-shot: raw grids — the kernel's `map(_resolve_axis, …)` shapes each - # axis and the search promote-compares, so no eager `Tg.(x)` copy. `Tr` from - # op-shape inference (`dL/h` floats Int). Batch keeps eager-convert (amortised). - Tg = _promote_grid_eltype(grids) + # Scalar one-shot: raw grids shaped per axis by `_resolve_axis(…, Tg)` at the value-matched + # float (Int grid + Float32 data → Float32) — no eager copy, witness `Tr` matches the eval. + Tg = _promote_grid_float(_promote_grid_eltype(grids), Tv) _validate_nd_grids(grids, data) Tr = _promote_eltype(_interp_op, Tg, Tv, promote_type(typeof.(query)...)) diff --git a/src/linear/nd/linear_nd_types.jl b/src/linear/nd/linear_nd_types.jl index f5e85581f..a4461329e 100644 --- a/src/linear/nd/linear_nd_types.jl +++ b/src/linear/nd/linear_nd_types.jl @@ -92,7 +92,7 @@ struct LinearInterpolantND{ bcs::NTuple{N, AbstractBC} = ntuple(_ -> NoBC(), Val(N)), store::StorePolicy = StorePolicy() ) where {Tg, Tv, N} - grids_c = map((g, bc) -> _own_or_ref_axis(_cache_axis(g, bc, Tg), Tg, store), grids, bcs) + grids_c = map((g, bc, T) -> _own_or_ref_axis(_cache_axis(g, bc, T), T, store), grids, bcs, ntuple(_ -> Tg, Val(N))) data_c = _own_or_ref_data(data, store) return new{Tg, Tv, N, typeof(grids_c), typeof(extraps), typeof(searches), typeof(data_c)}( grids_c, data_c, extraps, searches diff --git a/src/pchip/pchip_oneshot.jl b/src/pchip/pchip_oneshot.jl index c8ee1f868..c932b5490 100644 --- a/src/pchip/pchip_oneshot.jl +++ b/src/pchip/pchip_oneshot.jl @@ -27,9 +27,12 @@ # Grid pre-normalized by the public `pchip_interp` API via `_resolve_axis(x)` # before dispatching here; `_periodic_extend_1d` preserves the normalization. x_eff, y_ext, bc_eff, extrap_eff = _periodic_extend_1d(x, y, bc, extrap) - Tdy = _promote_eltype(_coeff_op, eltype(x_eff), Tv) + # Value-matched width: the dy buffer AND the slope arithmetic inside the + # filler both run at `Tw` (raw Int axes stop minting `inv(Int)::Float64`). + Tw = _promote_grid_float(eltype(x_eff), Tv) + Tdy = _promote_eltype(_coeff_op, Tw, Tv) dy = acquire!(pool, Tdy, length(y_ext)) - _pchip_slopes!(dy, x_eff, y_ext; bc = bc_eff) + _pchip_slopes!(dy, x_eff, y_ext, Tw; bc = bc_eff) searcher = _resolve_search(x_eff, xq, search, hint) return _hermite_eval_at_point(x_eff, y_ext, dy, xq, extrap_eff, deriv, searcher) end @@ -50,9 +53,12 @@ end @boundscheck length(output) == length(x_query) || _throw_length_mismatch(length(x_query), length(output), "x_query", "output") x_eff, y_ext, bc_eff, extrap_eff = _periodic_extend_1d(x, y, bc, extrap) - Tdy = _promote_eltype(_coeff_op, eltype(x_eff), Tv) + # Value-matched width: the dy buffer AND the slope arithmetic inside the + # filler both run at `Tw` (raw Int axes stop minting `inv(Int)::Float64`). + Tw = _promote_grid_float(eltype(x_eff), Tv) + Tdy = _promote_eltype(_coeff_op, Tw, Tv) dy = acquire!(pool, Tdy, length(y_ext)) - _pchip_slopes!(dy, x_eff, y_ext; bc = bc_eff) + _pchip_slopes!(dy, x_eff, y_ext, Tw; bc = bc_eff) searcher = _resolve_search(x_eff, x_query, search, hint) return _hermite_vector_loop!(output, x_eff, y_ext, dy, x_query, extrap_eff, deriv, searcher) end @@ -134,8 +140,9 @@ C\$^1\$ continuous, monotonicity guaranteed for monotone input data. # Unified entry — bc flows through the BC-aware extrap/search resolvers and # into the slope routines. No `_is_periodic_bc` branch, no extension copy: # `_resolve_search`'s seam dispatch + bc-aware slope formulas handle the - # closed-cycle on the user's n-length grid (Linear pattern). - x = _resolve_axis(x) + # closed-cycle on the user's n-length grid (Linear pattern). Value-matched Tg: + # Int/OneTo grid + Float32 data → Float32 axis. + x = _resolve_axis(x, _promote_grid_float(Tg, Tv)) extrap_eff = _resolve_extrap(extrap, bc, x, y) resolved = _resolve_coeffs(coeffs, x, xq) if resolved isa OnTheFly @@ -161,7 +168,7 @@ In-place PCHIP interpolation with monotone-preserving slopes. search::AbstractSearchPolicy = AutoSearch(), hint::Union{Nothing, Base.RefValue{Int}} = nothing ) where {Tg, Tv, Tq <: Real} - x = _resolve_axis(x) + x = _resolve_axis(x, _promote_grid_float(Tg, Tv)) extrap_eff = _resolve_extrap(extrap, bc, x, y) resolved = _resolve_coeffs(coeffs, x, x_query) if resolved isa OnTheFly diff --git a/src/pchip/pchip_slopes.jl b/src/pchip/pchip_slopes.jl index 97a48e6bf..3d371175b 100644 --- a/src/pchip/pchip_slopes.jl +++ b/src/pchip/pchip_slopes.jl @@ -73,12 +73,19 @@ Compute PCHIP (Fritsch-Carlson) monotone-preserving slopes in-place. # Complexity O(n), single pass, zero allocation (writes into `dy`). """ +# Width-less form: delegate with the axis's own eltype — bit-identical to the +# historic raw behavior. One-shot PreCompute backends pass the value-matched +# `Tw` explicitly so slopes are born at the value width (raw Int axes included). +_pchip_slopes!(dy::AbstractVector, x::AbstractVector, y::AbstractVector; bc::AbstractBC = NoBC()) = + _pchip_slopes!(dy, x, y, eltype(x); bc) + function _pchip_slopes!( dy::AbstractVector, - x::AbstractVector{Tg}, - y::AbstractVector; + x::AbstractVector, + y::AbstractVector, + ::Type{Tw}; bc::AbstractBC = NoBC() - ) where {Tg} + ) where {Tw} n = length(x) @assert n >= 2 "PCHIP requires at least 2 points" @assert length(y) == n "y length must match x" @@ -90,12 +97,12 @@ function _pchip_slopes!( # n=2 grids (where the seam secant can have opposite sign of cell-1). if n == 2 if bc isa PeriodicBC - @inbounds dy[1] = _pchip_boundary_slope(x, y, 1, n, bc) - @inbounds dy[2] = _pchip_boundary_slope(x, y, 2, n, bc) + @inbounds dy[1] = _pchip_boundary_slope(Tw, x, y, 1, n, bc) + @inbounds dy[2] = _pchip_boundary_slope(Tw, x, y, 2, n, bc) return dy end @inbounds begin - δ = _forward_secant(x, y, 1) + δ = _forward_secant(Tw, x, y, 1) dy[1] = δ dy[2] = δ end @@ -104,17 +111,19 @@ function _pchip_slopes!( Tc = eltype(dy) - # Compute secant slopes for first two intervals (needed for first interior) + # Compute secant slopes for first two intervals (needed for first interior). + # Cell widths stay raw — they are width-neutral WEIGHTS (Int × secant keeps + # the secant's type); only the secants carry the value-matched width. @inbounds h_prev = _get_h(x, 1) - @inbounds δ_prev = _forward_secant(x, y, 1) + @inbounds δ_prev = _forward_secant(Tw, x, y, 1) @inbounds h_curr = _get_h(x, 2) - @inbounds δ_curr = _forward_secant(x, y, 2) + @inbounds δ_curr = _forward_secant(Tw, x, y, 2) # Left endpoint: bc-dispatched helper. # NoBC: one-sided 3-point FD with monotonicity clamping. # PeriodicBC: closed-cycle interior formula via wrap-aware abstraction. - @inbounds dy[1] = _pchip_boundary_slope(x, y, 1, n, bc) + @inbounds dy[1] = _pchip_boundary_slope(Tw, x, y, 1, n, bc) # Interior slopes (k = 2:n-1) — unchanged. K=3 stencil never crosses join. @inbounds for k in 2:(n - 1) @@ -133,7 +142,7 @@ function _pchip_slopes!( h_prev = h_curr δ_prev = δ_curr h_curr = _get_h(x, k + 1) - δ_curr = _forward_secant(x, y, k + 1) + δ_curr = _forward_secant(Tw, x, y, k + 1) end end @@ -146,7 +155,7 @@ function _pchip_slopes!( # differs from i=1's (m_0=seam, m_1) — dy[1] ≠ dy[n] in general, both # correctly wrap-aware via the seam secant. # - NoBC: helper falls back to the original one-sided FD. - @inbounds dy[n] = _pchip_boundary_slope(x, y, n, n, bc) + @inbounds dy[n] = _pchip_boundary_slope(Tw, x, y, n, n, bc) return dy end diff --git a/src/quadratic/nd/quadratic_nd_adjoint_types.jl b/src/quadratic/nd/quadratic_nd_adjoint_types.jl index e3215d777..119718f64 100644 --- a/src/quadratic/nd/quadratic_nd_adjoint_types.jl +++ b/src/quadratic/nd/quadratic_nd_adjoint_types.jl @@ -87,7 +87,7 @@ struct QuadraticAdjointND{ ) where { Tg, N, BP <: NTuple{N, AbstractBC}, } - grids_c = map((g, bc) -> _convert_copy(_cache_axis(g, bc, Tg), Tg), grids, bcs) + grids_c = map((g, bc, T) -> _convert_copy(_cache_axis(g, bc, T), T), grids, bcs, ntuple(_ -> Tg, Val(N))) return new{Tg, N, typeof(grids_c), BP}(grids_c, bcs, anchors, grid_size, mincurv_Cs) end end diff --git a/src/quadratic/nd/quadratic_nd_interpolant.jl b/src/quadratic/nd/quadratic_nd_interpolant.jl index 178eb4284..688cce5f9 100644 --- a/src/quadratic/nd/quadratic_nd_interpolant.jl +++ b/src/quadratic/nd/quadratic_nd_interpolant.jl @@ -116,7 +116,7 @@ function _build_nd_quadratic_interpolant( # Cache axes for the build phase — inner ctor of `QuadraticInterpolantND` # handles the owned `_convert_copy` separately, so we only wrap (no copy) # here. Already-cached axes pass through idempotently in the ctor. - grids_cached = map((g, bc) -> _cache_axis(g, bc, Tg), grids, bcs) + grids_cached = map(_cache_axis, grids, bcs, ntuple(_ -> Tg, Val(N))) nodal_derivs = _build_nd_coeffs_quadratic(grids_cached, data, bcs) return QuadraticInterpolantND(grids_cached, nodal_derivs, bcs, extraps_val, searches) diff --git a/src/quadratic/nd/quadratic_nd_oneshot.jl b/src/quadratic/nd/quadratic_nd_oneshot.jl index 578090152..ff024366c 100644 --- a/src/quadratic/nd/quadratic_nd_oneshot.jl +++ b/src/quadratic/nd/quadratic_nd_oneshot.jl @@ -39,14 +39,14 @@ Zero-allocation after warmup (pool reuse). oob_result = _try_fill_oob(query, grids, extraps_val, ops, @inbounds first(data)) oob_result !== nothing && return oob_result - # 1. Pool-backed per-axis cache — build phase reuses h/inv_h across slices. - # `map` over the tuple dispatches per-element on the concrete axis type - # (Range vs Vector); `ntuple(d -> grids[d], ...)` would Union-box. - grids_c = map(g -> _cache_axis_pooled(pool, g), grids) - - # 2. Pool-allocate partials array (THE KEY: pool instead of heap) - # Tz widens Tv with Tg (Dual grid → Dual derivs). `Tg` raw; `_coeff_op` floats Int. - Tg = _promote_grid_eltype(grids) + # 1. Value-matched grid float (Int/OneTo grid + Float32 data → Float32) shared by the pooled + # per-axis cache and the partials type, so the pool buffer + output follow natural promote_type. + # `map` dispatches per-element on the concrete axis type (Range vs Vector). + Tg = _promote_grid_float(_promote_grid_eltype(grids), Tv) + grids_c = _cache_axes_pooled(pool, grids, Tg) # @generated static-Tg unroll (no Type-captured closure) + + # 2. Pool-allocate partials array (THE KEY: pool instead of heap). Tz widens Tv with Tg + # (Dual grid → Dual derivs); `_coeff_op` floats Int. Tz = _promote_eltype(_coeff_op, Tg, Tv) n_partials = 1 << N partials = acquire!(pool, Tz, (n_partials, size(data)...)) @@ -135,21 +135,13 @@ One-shot ND quadratic interpolation at a single point. Zero-allocation after warmup. # Strategy selection (`coeffs`) -- `AutoCoeffs()` (default): resolves to `PreCompute()` for both scalar and - batch quadratic queries. This matches the interpolant constructor and AD - rules **bit-exactly**, at the cost of always building the `2^N` partial - array. Cubic `AutoCoeffs()` differs (it routes scalar queries to `OnTheFly`) - because cubic OnTheFly↔PreCompute equivalence holds at `≈ rtol=1e-10` and - AD seed agreement does not require strict equality there. -- `PreCompute()`: explicitly build all nodal partial derivatives first, then - evaluate. The user `bc` is applied uniformly to every partial (pure and - mixed), so the stored mixed partials match OnTheFly's composition formula. -- `OnTheFly()`: sequential 1D interpolation per fiber via `_collapse_dims`. - Applies the user-specified `bc` uniformly at every 1D step. - -`PreCompute()` and `OnTheFly()` are equivalent up to a few ULPs of FP -reordering noise for every user BC, but **not** bit-identical — that is why -quadratic `AutoCoeffs` defaults to `PreCompute`. +- `AutoCoeffs()` (default): shared `_resolve_coeffs_nd_oneshot` policy, same as cubic — + scalar query → `OnTheFly()` (cell-local collapse, skips the `2^N` partial build), + batch → `PreCompute()` (amortized). +- `PreCompute()`: build all nodal partials first — bit-identical to the persistent + interpolant (`quadratic_interp(grids, data)(q)`); opt in when `===` parity matters. +- `OnTheFly()`: sequential 1D collapse; matches PreCompute to a few ULPs of + FP-reordering noise (not bit-exact) for every user `bc`. """ function quadratic_interp( grids::NTuple{N, AbstractVector}, @@ -163,8 +155,9 @@ function quadratic_interp( hint::Union{Nothing, NTuple{N, Base.RefValue{Int}}} = nothing ) where {Tv, N} # Scalar one-shot: raw grids — both PreCompute and OnTheFly accept them - # (`_compute_all_local_params` floats the cell width). Batch keeps eager-convert. - Tg = _promote_grid_eltype(grids) + # (`_compute_all_local_params` floats the cell width). `Tg` value-matched (Int/OneTo grid + + # Float32 data → Float32) so eval + witness `Tr` agree. Batch keeps eager-convert. + Tg = _promote_grid_float(_promote_grid_eltype(grids), Tv) _validate_nd_grids(grids, data) Tr = _promote_eltype(_interp_op, Tg, Tv, promote_type(typeof.(query)...)) @@ -174,9 +167,9 @@ function quadratic_interp( extraps_val = _resolve_extrap(extrap, bcs, Val(N), Tv) ops = _resolve_deriv_nd(deriv, Val(N)) - # Keep AutoCoeffs on the specialized PreCompute path so scalar one-shot - # calls match the interpolant constructor and AD rules. - coeffs_resolved = coeffs isa AutoCoeffs ? PreCompute() : coeffs + # Central policy (same as cubic): scalar AutoCoeffs → OnTheFly. Bit-exact parity + # with the persistent interpolant stays available via explicit coeffs=PreCompute(). + coeffs_resolved = _resolve_coeffs_nd_oneshot(coeffs, query, ntuple(_ -> QuadraticInterp(), Val(N))) if coeffs_resolved isa OnTheFly sample = @inbounds first(data) methods = map(bc_i -> QuadraticInterp(_to_quadratic_bc(bc_i, sample)), bcs) diff --git a/src/quadratic/nd/quadratic_nd_types.jl b/src/quadratic/nd/quadratic_nd_types.jl index 72f07ef72..9f07855a5 100644 --- a/src/quadratic/nd/quadratic_nd_types.jl +++ b/src/quadratic/nd/quadratic_nd_types.jl @@ -90,7 +90,7 @@ struct QuadraticInterpolantND{ searches::Tuple{Vararg{AbstractSearchPolicy, N}} ) where {Tg, Tv, N, NP1} NP1 == N + 1 || throw(ArgumentError("NP1 must equal N+1")) - grids_c = map((g, bc) -> _convert_copy(_cache_axis(g, bc, Tg), Tg), grids, bcs) + grids_c = map((g, bc, T) -> _convert_copy(_cache_axis(g, bc, T), T), grids, bcs, ntuple(_ -> Tg, Val(N))) return new{Tg, Tv, N, NP1, typeof(grids_c), typeof(bcs), typeof(extraps), typeof(searches)}( grids_c, nodal_derivs, bcs, extraps, searches ) diff --git a/src/quadratic/quadratic_oneshot.jl b/src/quadratic/quadratic_oneshot.jl index d67ebba0e..ae5afcb77 100644 --- a/src/quadratic/quadratic_oneshot.jl +++ b/src/quadratic/quadratic_oneshot.jl @@ -169,7 +169,9 @@ vals = quadratic_interp(x, y, sorted_queries; search=LinearBinarySearch(linear_w @boundscheck length(y) == length(x) || throw(ArgumentError("x and y must have same length")) @boundscheck length(x) >= 2 || throw(ArgumentError("x must have at least 2 elements")) - x = _cache_axis_pooled(pool, x) + # Value-matched pooled wrap: Int/OneTo grid + Float32 data → Float32 axis + # (vector conversion lands in a pool buffer — warm one-shots stay zero-alloc). + x = _cache_axis_pooled(pool, x, _promote_grid_float(Tg, Tv)) # Compute coefficients using temporary arrays from pool. The grid `x` # carries cached `h`/`inv_h` when wrapped, or computes them on the fly. nx = length(x) @@ -226,7 +228,8 @@ quadratic_interp!(output, x, y, sorted_queries; search=LinearBinarySearch(linear @assert length(output) == length(x_targets) "output must match x_targets length" @assert length(x) >= 2 "x must have at least 2 elements" - x = _cache_axis_pooled(pool, x) + # Value-matched pooled wrap — keeps the batch interior consistent with scalar. + x = _cache_axis_pooled(pool, x, _promote_grid_float(Tg, Tv)) # Compute coefficients using temporary arrays from pool. The grid `x` # carries cached `h`/`inv_h` when wrapped, or computes them on the fly. nx = length(x) diff --git a/test/setup.jl b/test/setup.jl index 420416148..69a6e909a 100644 --- a/test/setup.jl +++ b/test/setup.jl @@ -10,6 +10,7 @@ using TestItemRunner # accessing FastInterpolations.AdaptiveArrayPools.RUNTIME_CHECK is safe. @testsnippet AllocConstants begin const AAP_RUNTIME_CHECK = FastInterpolations.AdaptiveArrayPools.RUNTIME_CHECK + # LTS keeps a small 240-byte margin for genuine warm-path noise; 1.12+ is strict (0). const ALLOC_THRESHOLD = VERSION >= v"1.12" ? 0 : (2 * AAP_RUNTIME_CHECK + 1) * 240 const ND_ALLOC_THRESHOLD = VERSION >= v"1.12" ? 0 : (2 * AAP_RUNTIME_CHECK + 1) * 240 end diff --git a/test/test_1d_output_type_promotion.jl b/test/test_1d_output_type_promotion.jl new file mode 100644 index 000000000..41d88b71d --- /dev/null +++ b/test/test_1d_output_type_promotion.jl @@ -0,0 +1,152 @@ +# Value-match contract for 1D ONE-SHOT paths, mirroring test_nd_output_type_promotion.jl. +# Arithmetic kernels: output = promote_type(value-matched grid float, data, query), +# where the value-matched grid float = float(promote(grid eltype, data eltype)) — an +# Int/OneTo grid beside Float32 data floats to Float32, not the blind Float64. +# Selection kernel (constant): pure natural promote_type(grid, data, query). +# Persistent 1D already conforms (via the 3-arg `_cache_axis`) — guarded here too. + +@testitem "1D scalar one-shot output types: arithmetic value-match, @inferred" begin + # Direct function-value calls (NOT closures): a closure capturing testitem-level + # locals boxes them as `Any`, which poisons `@inferred` for every cell. + data32 = Float32[sin(0.4i) + 0.5 for i in 1:7] + dy32 = Float32[0.4cos(0.4i) for i in 1:7] + + gridspecs = ( + ("IntRange", Int, () -> 1:7), + ("IntOneTo", Int, () -> Base.OneTo(7)), + ("IntVec", Int, () -> collect(1:7)), + ("F32range", Float32, () -> 1.0f0:1.0f0:7.0f0), + ("F32vec", Float32, () -> collect(Float32, 1:7)), + ) + fns = (linear_interp, quadratic_interp, cubic_interp, pchip_interp, akima_interp, cardinal_interp) + + for (gname, Te, gbuild) in gridspecs, f in fns, q in (2.4f0, 2.4) + Tq = typeof(q) + Tg = float(promote_type(Te, Float32)) + expected = promote_type(Tg, Float32, Tq) + g = gbuild() + @testset "$(nameof(f)) $gname Tq=$Tq → $expected" begin + # Every 1D one-shot value-matches: an Int-VECTOR axis beside Float32 data + # returns Float32. Cubic's spline cache is now data-aware (the eltype-aware + # bank routes Int + Float32 → a Float32 cache), so scalar ≡ batch ≡ persistent. + @test f(g, data32, q) isa expected + @test (@inferred f(g, data32, q)) isa expected + end + end + + # hermite carries user slopes — same value-match rule with value space y ∪ dy. + for (gname, Te, gbuild) in gridspecs, q in (2.4f0, 2.4) + Tq = typeof(q) + expected = promote_type(float(promote_type(Te, Float32)), Float32, Tq) + g = gbuild() + @testset "hermite $gname Tq=$Tq → $expected" begin + @test hermite_interp(g, data32, dy32, q) isa expected + @test (@inferred hermite_interp(g, data32, dy32, q)) isa expected + end + end +end + +@testitem "1D Series one-shot output eltypes: value-matched (guard)" begin + # The Series wrappers already value-match their axes at the entry (arith families + # via `_promote_grid_float(Tg, _series_eltype(s))`; constant via raw Tg) — pin it. + data32a = Float32[sin(0.4i) for i in 1:7] + data32b = Float32[cos(0.4i) for i in 1:7] + for g in (1:7, collect(1:7)) + s = Series(data32a, data32b) + @test eltype(linear_interp(g, s, 2.4f0)) === Float32 + @test eltype(quadratic_interp(g, s, 2.4f0)) === Float32 + @test eltype(cubic_interp(g, s, 2.4f0)) === Float32 + @test eltype(constant_interp(g, s, 2.4f0)) === Float32 + end +end + +@testitem "1D persistent parity guard: same value-matched output types" begin + data32 = Float32[sin(0.4i) + 0.5 for i in 1:7] + dy32 = Float32[0.4cos(0.4i) for i in 1:7] + + gridspecs = ( + ("IntRange", Int, () -> 1:7), + ("IntVec", Int, () -> collect(1:7)), + ("F32vec", Float32, () -> collect(Float32, 1:7)), + ) + builders = ( + ("linear", g -> linear_interp(g, data32)), + ("quadratic", g -> quadratic_interp(g, data32)), + ("cubic", g -> cubic_interp(g, data32)), + ("pchip", g -> pchip_interp(g, data32)), + ("akima", g -> akima_interp(g, data32)), + ("cardinal", g -> cardinal_interp(g, data32)), + ("hermite", g -> hermite_interp(g, data32, dy32)), + ) + + for (gname, Te, gbuild) in gridspecs, (mname, build) in builders, q in (2.4f0, 2.4) + Tq = typeof(q) + expected = promote_type(float(promote_type(Te, Float32)), Float32, Tq) + itp = build(gbuild()) + @testset "$mname $gname Tq=$Tq → $expected" begin + @test itp(q) isa expected + end + end +end + +@testitem "1D batch one-shot output eltypes: value-matched" begin + data32 = Float32[sin(0.4i) + 0.5 for i in 1:7] + dy32 = Float32[0.4cos(0.4i) for i in 1:7] + qs32 = Float32[2.4, 3.1] + + fns = ( + ("linear", g -> linear_interp(g, data32, qs32)), + ("quadratic", g -> quadratic_interp(g, data32, qs32)), + ("cubic", g -> cubic_interp(g, data32, qs32)), + ("pchip", g -> pchip_interp(g, data32, qs32)), + ("akima", g -> akima_interp(g, data32, qs32)), + ("cardinal", g -> cardinal_interp(g, data32, qs32)), + ("hermite", g -> hermite_interp(g, data32, dy32, qs32)), + ) + + for (gname, g) in (("IntRange", 1:7), ("IntVec", collect(1:7))), (mname, f) in fns + @testset "$mname $gname → Vector{Float32}" begin + @test eltype(f(g)) === Float32 + end + end +end + +@testitem "1D constant: selection kernel stays natural (no float forcing)" begin + data32 = Float32[sin(0.4i) + 0.5 for i in 1:7] + dataI = [2i for i in 1:7] + + for (gname, g) in (("IntRange", 1:7), ("IntOneTo", Base.OneTo(7)), ("IntVec", collect(1:7))) + @testset "constant $gname" begin + # F32 data: natural promote(Int, F32, Tq) + @test constant_interp(g, data32, 2.4f0) isa Float32 + @test constant_interp(g, data32, 2.4) isa Float64 + @test constant_interp(g, data32)(2.4f0) isa Float32 # persistent parity + # all-Int: nearest-neighbor selection returns Int for Int query + @test constant_interp(g, dataI, 3) isa Int + @test constant_interp(g, dataI, 2.4f0) isa Float32 + @test constant_interp(g, dataI)(3) isa Int # persistent parity + # batch eltype follows the same natural rule + @test eltype(constant_interp(g, data32, Float32[2.4, 3.1])) === Float32 + end + end +end + +@testitem "1D hermite mixed y/dy widths: value space = y ∪ dy" begin + data32 = Float32[sin(0.4i) + 0.5 for i in 1:7] + dy64 = [0.4cos(0.4i) for i in 1:7] # Float64 slopes beside Float32 data + + for (gname, g) in (("IntRange", 1:7), ("F32vec", collect(Float32, 1:7))) + @testset "hermite $gname F32 data + F64 dy → Float64" begin + @test hermite_interp(g, data32, dy64, 2.4f0) isa Float64 + end + end +end + +@testitem "1D Int-data arithmetic: float-forces to Float64 (unchanged legacy)" begin + dataI = [2i for i in 1:7] + # All-Int value space: float(promote(Int, Int)) = Float64 — arithmetic kernels divide. + for f in (linear_interp, cubic_interp, pchip_interp) + @test f(1:7, dataI, 2.4f0) isa Float64 + @test f(1:7, dataI, 3) isa Float64 + end +end diff --git a/test/test_axis_data_resolvers.jl b/test/test_axis_data_resolvers.jl index 013c20248..59a8d948f 100644 --- a/test/test_axis_data_resolvers.jl +++ b/test/test_axis_data_resolvers.jl @@ -45,6 +45,100 @@ end end +@testitem "_resolve_axis 3-arg Tg dispatch table (incl. pre-wrapped × :exclusive diagonal)" begin + using FastInterpolations: + _resolve_axis, _CachedRange, _CachedVector, _ExclusivePeriodicAxis, + _to_float, NoBC, PeriodicBC + + # 3-arg `_resolve_axis(x, bc, Tg)` dispatch table. The pre-wrapped × `:exclusive` + # DIAGONAL methods are load-bearing against ambiguity: `_Cached* <: Abstract*`, so + # (specific container × generic BC) and (generic container × specific BC) cross + # without them — Aqua ambiguity RED + runtime MethodError in hetero one-shot. + bc_no = NoBC() + bc_excl = PeriodicBC(endpoint = :exclusive, period = 4.0) + + @testset "raw Range + Tg → _CachedRange{Tg}" begin + @test _resolve_axis(0:1:3, bc_no, Float32) isa _CachedRange{Float32} + @test _resolve_axis(0.0:1.0:3.0, bc_no, Float64) isa _CachedRange{Float64} + end + + @testset "_CachedRange + non-exclusive + Tg → _convert_copy (same-type identity)" begin + cr64 = _to_float(0.0:1.0:3.0, Float64) + @test _resolve_axis(cr64, bc_no, Float64) === cr64 + @test _resolve_axis(cr64, bc_no, Float32) isa _CachedRange{Float32} + end + + @testset "raw Vector / _CachedVector + non-exclusive + Tg → passthrough (Tg ignored)" begin + x32 = Float32[0.0, 1.0, 2.0, 3.0] + cv64 = _CachedVector([0.0, 1.0, 2.0, 3.0]) + @test _resolve_axis(x32, bc_no, Float64) === x32 + @test _resolve_axis(cv64, bc_no, Float32) === cv64 + end + + @testset "raw Range/Vector + :exclusive + Tg → _ExclusivePeriodicAxis" begin + ax_r = _resolve_axis(0:1:3, bc_excl, Float32) + @test ax_r isa _ExclusivePeriodicAxis + @test ax_r.inner isa _CachedRange{Float32} # Tg respected on raw Range + # Period follows the Tg-typed axis (`_resolve_bc_period` normalizes it): a + # Float64 period literal must NOT re-widen a value-matched Float32 axis. + @test ax_r.period isa Float32 + @test eltype(ax_r) === Float32 + x = [0.0, 1.0, 2.0, 3.0] + ax_v = _resolve_axis(x, bc_excl, Float64) + @test ax_v isa _ExclusivePeriodicAxis + @test ax_v.inner === x # raw Vector never converts + end + + @testset "DIAGONAL: _CachedRange + :exclusive + Tg" begin + cr64 = _to_float(0.0:1.0:3.0, Float64) + ax = _resolve_axis(cr64, bc_excl, Float64) + @test ax isa _ExclusivePeriodicAxis + @test ax.inner === cr64 # same-type `_convert_copy` identity + @test ax.period == 4.0 + @test length(ax) == 5 # virtual n+1 + # Cross-eltype: the axis converts to Tg first, and the period is resolved + # against the CONVERTED axis — so it follows Tg instead of re-widening. + ax32 = _resolve_axis(cr64, bc_excl, Float32) + @test ax32.inner isa _CachedRange{Float32} + @test ax32.period isa Float32 + end + + @testset "DIAGONAL: _CachedVector + :exclusive + Tg" begin + cv64 = _CachedVector([0.0, 1.0, 2.0, 3.0]) + ax = _resolve_axis(cv64, bc_excl, Float64) + @test ax isa _ExclusivePeriodicAxis + @test ax.inner isa _CachedVector{Float64} # Tg respected (same-width converted copy) + @test ax.period == 4.0 + # `:exclusive` converts to Tg — UNLIKE a non-periodic vector, the wrapped axis eltype flows + # into the value/witness path via the period seam, so it must follow Tg (matching the + # `_CachedRange` diagonal above and the raw-Vector arm). The old "never convert" contract + # produced a Float64 axis past a value-matched Float32 witness → crash. + ax32 = _resolve_axis(cv64, bc_excl, Float32) + @test ax32.inner isa _CachedVector{Float32} + @test ax32.period isa Float32 + end + + @testset "diagonals are type-stable" begin + cr64 = _to_float(0.0:1.0:3.0, Float64) + cv64 = _CachedVector([0.0, 1.0, 2.0, 3.0]) + f(x, bc, ::Type{T}) where {T} = _resolve_axis(x, bc, T) + @test (@inferred f(cr64, bc_excl, Float64)) isa _ExclusivePeriodicAxis + @test (@inferred f(cv64, bc_excl, Float64)) isa _ExclusivePeriodicAxis + end + + @testset "Tg-typed 2-arg (no BC): Ranges value-match, Vectors pass through" begin + # Used by the 1D one-shot entries (no bc concept at the normalize point). + @test _resolve_axis(0:1:3, Float32) isa _CachedRange{Float32} + cr64 = _to_float(0.0:1.0:3.0, Float64) + @test _resolve_axis(cr64, Float64) === cr64 # same-type identity + @test _resolve_axis(cr64, Float32) isa _CachedRange{Float32} + x32 = Float32[0.0, 1.0, 2.0, 3.0] + @test _resolve_axis(x32, Float64) === x32 # vectors never convert + cv64 = _CachedVector([0.0, 1.0, 2.0, 3.0]) + @test _resolve_axis(cv64, Float32) === cv64 + end +end + @testitem "_resolve_data — type correctness" begin using FastInterpolations: _resolve_data, _ExclusivePeriodicData, NoBC, PeriodicBC @@ -531,6 +625,16 @@ end @test ax isa _ExclusivePeriodicAxis @test ax.inner isa _CachedRange{Float64} # inner promoted to Tg end + + @testset "Int range + :exclusive + Float64 period + Tg=Float32 stays Float32" begin + # Convert-first contract: the period resolves against the Tg-typed axis, + # so a Float64 period literal cannot re-widen the value-matched axis. + r = 0:1:3 + ax = _cache_axis(r, bc_excl, Float32) + @test ax isa _ExclusivePeriodicAxis + @test ax.inner isa _CachedRange{Float32} + @test ax.period isa Float32 + end end @testitem "_cache_axis 3-arg Tg INTENTIONALLY IGNORED: pre-wrapped inputs" begin @@ -828,3 +932,115 @@ end @test (@allocated _measure_view(vw)) <= WRAPPER_OVERHEAD_LIMIT end end + +# RED PIN (#7 root / #1): the one-shot `_resolve_axis` :exclusive Vector arm must respect the +# value-matched `Tg` exactly as the persistent `_cache_axis` does. Today `_resolve_axis` resolves +# the period against the RAW grid (→ `float(Int)` = Float64) while `_cache_axis` converts-first +# (→ Float32), so an Int Vector + Tg=Float32 diverges: Float64 axis vs Float32 axis. That divergence +# IS the linear one-shot crash and the 2-arg/3-arg period-timing split — a `_wrap_exclusive` +# helper shared by both families (always convert-first) would close it. +@testitem "_resolve_axis :exclusive Vector respects Tg like _cache_axis (narrow-float root)" begin + using FastInterpolations: _resolve_axis, _cache_axis, PeriodicBC + + bc = PeriodicBC(endpoint = :exclusive, period = 2.5) + x = collect(0:2) # Int Vector; values exact in Float32 + + axr = _resolve_axis(x, bc, Float32) # one-shot arm + axc = _cache_axis(x, bc, Float32) # persistent arm (reference; already narrow) + + @test eltype(axc) === Float32 # sanity: persistent narrows correctly + @test eltype(axr) === Float32 # RED: one-shot currently Float64 + @test axr.period isa Float32 # RED: one-shot currently Float64 + @test eltype(axr) === eltype(axc) # one-shot ≡ persistent axis width +end + +# RED PIN (Copilot #182): `_resolve_axis` only has the 1-arg `_ExclusivePeriodicAxis` passthrough +# (periodic_axis.jl), so the 2-arg/3-arg `:exclusive` forms fall through to the AbstractVector arms +# and RE-WRAP an already-wrapped axis (`_ExclusivePeriodicAxis <: AbstractVector`) — nesting toward +# length (n+1)+1, which actually throws in the ctor (`inner[end] < x_max` fails). `_cache_axis` +# already defends this with a full passthrough set; `_resolve_axis` must mirror it. +@testitem "_resolve_axis :exclusive passes through a pre-wrapped axis (no re-wrap)" begin + using FastInterpolations: _resolve_axis, _wrap_exclusive, _ExclusivePeriodicAxis, PeriodicBC + + bc = PeriodicBC(endpoint = :exclusive, period = 4.0) + wrapped = _wrap_exclusive(collect(0.0:3.0), bc) # _ExclusivePeriodicAxis, length n+1 = 5 + @test wrapped isa _ExclusivePeriodicAxis + @test length(wrapped) == 5 + + # Feeding the already-wrapped axis back must PASS THROUGH (identity), not nest/throw. + @test _resolve_axis(wrapped, bc) === wrapped # 2-arg + @test _resolve_axis(wrapped, bc, Float64) === wrapped # 3-arg Tg-aware + @test length(_resolve_axis(wrapped, bc, Float64)) == 5 # not (n+1)+1 +end + +# Width-first geometry primitives: `_get_*(Tw, x, i)` — the value-matched coordinate +# width `Tw` comes from the caller's surface (`_promote_grid_float(Tg, Tv)`). Raw axes +# difference in their OWN eltype first (Int spans are exact), convert the SPAN once, +# then divide — the reciprocal is BORN at `Tw` (no `inv(Int)::Float64` minting), and +# coordinates beyond `Tw`'s ulp cannot cancel (span-first, never endpoint-convert). +# Wrapped axes reuse the cached reciprocal (convert is a no-op once value-matched). +# Width-less forms keep the historic raw-eltype behavior via delegation. +@testitem "width-first _get_h/_get_inv_h/_get_inv_2cell + secants" begin + using FastInterpolations: _get_h, _get_inv_h, _get_inv_2cell, + _forward_secant, _backward_secant, _centered_secant, + _cache_axis, _resolve_axis, NoBC + + # raw Int vector: reciprocal born at Tw + x = [1, 3, 6] + @test _get_h(Float32, x, 1) === 2.0f0 + @test _get_inv_h(Float32, x, 1) === 0.5f0 + @test _get_inv_2cell(Float32, x, 2) === inv(5.0f0) + @test _get_inv_h(Float64, x, 2) === inv(3.0) + + # span-first precision: Int coords beyond Float32's ulp — endpoint-convert + # would cancel to 0 (inv → Inf); the span itself is small and exact. + xb = [16_777_216, 16_777_218] + @test _get_inv_h(Float32, xb, 1) === 0.5f0 + + # wrapped axes: cached reciprocal reused; convert no-op when value-matched + c = _cache_axis(collect(Float32, 1:5), NoBC()) + @test _get_inv_h(Float32, c, 2) === _get_inv_h(c, 2) + r = _resolve_axis(1:5, Float32) # _CachedRange{Float32} + @test _get_inv_h(Float32, r, 1) === 1.0f0 + @test _get_h(Float32, r, 1) === 1.0f0 + + # width-first secants: Int axis + F32 data → Float32 end to end + y32 = Float32[1, 2, 4] + @test _forward_secant(Float32, x, y32, 1) === 0.5f0 + @test _backward_secant(Float32, x, y32, 2) === 0.5f0 + @test _centered_secant(Float32, x, y32, 2) === 0.6f0 + # width-less forms keep the historic raw semantics (Int axis → Float64) + @test _forward_secant(x, y32, 1) isa Float64 +end + +# Width-first SEARCH-RESULT form `_get_inv_h(Tw, g, idx, xL, xR)` — the ND linear +# scalar one-shot derives inv_h from the search endpoints. Same span-first doctrine +# as the 3-arg family: raw axes difference in their OWN eltype, convert the span +# once, divide at `Tw` (no `inv(Int)::Float64` minting). Cached/wrapped axes keep +# their cached reciprocal (endpoint args ignored; convert no-op once value-matched). +@testitem "width-first search-result _get_inv_h(Tw, g, idx, xL, xR)" begin + using FastInterpolations: _get_inv_h, _cache_axis, _resolve_axis, + NoBC, PeriodicBC + + # raw Int vector: reciprocal born at Tw from the search endpoints + x = [1, 3, 6] + @test _get_inv_h(Float32, x, 1, 1, 3) === 0.5f0 + @test _get_inv_h(Float64, x, 2, 3, 6) === inv(3.0) + # Int coords beyond Float32's ulp: the span is small and exact — + # endpoint-convert would cancel to 0 (inv → Inf) + xb = [16_777_216, 16_777_218] + @test _get_inv_h(Float32, xb, 1, 16_777_216, 16_777_218) === 0.5f0 + + # wrapped axes: cached reciprocal reused; endpoints ignored + c = _cache_axis(collect(Float32, 1:5), NoBC()) + @test _get_inv_h(Float32, c, 2, 2.0f0, 3.0f0) === _get_inv_h(c, 2) + r = _resolve_axis(1:5, Float32) # _CachedRange{Float32} + @test _get_inv_h(Float32, r, 1, 1.0f0, 2.0f0) === 1.0f0 + + # :exclusive axis: interior delegates to the wrapped inner; the seam cell has + # no stored width — span-first from the search endpoints (wrap domain edge) + bc = PeriodicBC(endpoint = :exclusive, period = 8.0) + g = _resolve_axis([0, 2, 4, 6], bc, Float32) # inner n=4, seam idx=4 + @test _get_inv_h(Float32, g, 2, 2.0f0, 4.0f0) === 0.5f0 + @test _get_inv_h(Float32, g, 4, 6.0f0, 8.0f0) === 0.5f0 +end diff --git a/test/test_cubic_autocache.jl b/test/test_cubic_autocache.jl index 55172abf9..cdccbfed7 100644 --- a/test/test_cubic_autocache.jl +++ b/test/test_cubic_autocache.jl @@ -1015,3 +1015,67 @@ end @test cache3 !== cache @test cache3.x isa FastInterpolations._CachedRange{Float64} end + +# ============================================================================= +# Eltype-aware bank: an explicit target float type routes an Int/Rational grid +# to the value-matched cache bank (data-aware), instead of the grid's blind +# natural float. The default (no target type) is unchanged. +# ============================================================================= +@testitem "Cubic Cache: eltype-aware bank (data-aware float type)" begin + import FastInterpolations: _get_cubic_cache + + clear_cubic_cache!() + xi = [0, 1, 2, 3, 4, 5, 6, 7] # Vector{Int} + bc = BCPair(Deriv1(0.0), Deriv1(0.0)) # the pre-normalized form the one-shot passes + + # Default 3-arg: an Int grid caches at its natural float (Float64). + c_default = _get_cubic_cache(xi, bc, true) + @test eltype(c_default.x) === Float64 + + # NEW 4-arg: an explicit target float type routes the Int grid to that bank. + c_f32 = _get_cubic_cache(xi, bc, true, Float32) + c_f64 = _get_cubic_cache(xi, bc, true, Float64) + @test eltype(c_f32.x) === Float32 + @test eltype(c_f64.x) === Float64 + + # Distinct banks (Float32 vs Float64) — no objectid collision on one Int grid. + @test c_f32 !== c_f64 + # Explicit-Float64 lands in the SAME bank as the default (grid natural float), + # so the common Int + Float64 case never fragments the cache. + @test c_f64 === c_default + # Warm re-lookup: same Int grid + same target → same cache object. + @test _get_cubic_cache(xi, bc, true, Float32) === c_f32 +end + +# ============================================================================= +# 1D one-shot value-match: Int Vector + Float32 data returns Float32 and agrees +# bit-for-bit across scalar / batch / persistent (all F32-solve). Warm scalar +# one-shot on the raw Int grid stays zero-alloc (cache hit via the Int objectid). +# ============================================================================= +@testitem "Cubic 1D one-shot: Int Vector + Float32 value-matches (scalar ≡ batch ≡ persistent)" setup = [AllocConstants] begin + xi = [0, 1, 2, 3, 4, 5, 6, 7] + y = Float32[sin(1.3f0 * i) + 0.4f0 * i for i in xi] + q = 3.4f0 + + # value-matched: Int axis + Float32 data → Float32 out (not the blind Float64). + v_scalar = cubic_interp(xi, y, q) + @test v_scalar isa Float32 + + # scalar ≡ batch ≡ persistent (all F32-solve). + v_batch = cubic_interp(xi, y, [q])[1] + v_persist = cubic_interp(xi, y)(q) + @test v_scalar === v_batch + @test v_scalar === v_persist + + # warm zero-alloc on the raw Int grid (cache hit via the Int grid's objectid). + function _alloc_cubic_scalar_int_f32() + x = [0, 1, 2, 3, 4, 5, 6, 7] + yy = Float32[sin(1.3f0 * i) for i in x] + qq = 3.4f0 + for _ in 1:3 + cubic_interp(x, yy, qq) + end + @allocated cubic_interp(x, yy, qq) + end + @test _alloc_cubic_scalar_int_f32() <= ALLOC_THRESHOLD +end diff --git a/test/test_hermite_nd_partials.jl b/test/test_hermite_nd_partials.jl index ea71ce12b..220cb58e7 100644 --- a/test/test_hermite_nd_partials.jl +++ b/test/test_hermite_nd_partials.jl @@ -805,3 +805,36 @@ end end # @testitem + +@testitem "Mixed data/partial widths: grid promotion sees partials (persistent ≡ one-shot)" begin + # The persistent build's grid value-match must see `data ∪ partials` — matching Float32 + # data alone would store Int coordinates as Float32, and above 2^24 those quantize: + # the persistent result then drifts from the one-shot path (which already promotes + # its grid width against data ∪ partials). + n = 6 + base = 16_777_220 # 2^24 = 16_777_216: Float32 integer limit + gx = base .+ collect(1:n) + gy = collect(1:n) + f(x, y) = 0.001x + 2.0y + data = Float32[f(x, y) for x in gx, y in gy] + parts = HermitePartials( + (1, 0) => [0.001 for _ in 1:n, _ in 1:n], + (0, 1) => [2.0 for _ in 1:n, _ in 1:n], + (1, 1) => [0.0 for _ in 1:n, _ in 1:n], + ) + + itp = hermite_interp((gx, gy), data, parts) + @test map(eltype, itp.grids) === (Float64, Float64) # data(F32) ∪ partials(F64) → F64 + q = (Float64(base) + 2.3, 3.7) + @test itp(q) ≈ hermite_interp((gx, gy), data, parts, q) rtol = 1.0e-12 + + # Guard the other direction: an all-Float32 value space must NOT over-widen the grids. + data32 = Float32[f(x, y) for x in 1:n, y in 1:n] + parts32 = HermitePartials( + (1, 0) => Float32[0.001 for _ in 1:n, _ in 1:n], + (0, 1) => Float32[2.0 for _ in 1:n, _ in 1:n], + (1, 1) => Float32[0.0 for _ in 1:n, _ in 1:n], + ) + itp32 = hermite_interp((collect(1:n), collect(1:n)), data32, parts32) + @test map(eltype, itp32.grids) === (Float32, Float32) +end diff --git a/test/test_nd_exclusive_period_oneshot.jl b/test/test_nd_exclusive_period_oneshot.jl index 4fa1c0f25..3ec1d8b12 100644 --- a/test/test_nd_exclusive_period_oneshot.jl +++ b/test/test_nd_exclusive_period_oneshot.jl @@ -103,3 +103,47 @@ end end end end + +# RED PIN (#1): narrow-float (Float32) exclusive-periodic on an Int Vector grid. The arithmetic +# kernels (linear/cubic) value-match `Tg` to the DATA width (Int grid + Float32 data → Float32), but +# the one-shot `_resolve_axis` exclusive-Vector arm resolved the period against the RAW Int grid → +# Float64 axis → the Float32 witness `::Tr` threw. Persistent already returned Float32; the one-shot +# must match. Range grids already worked via the convert-first arm; only the Vector arms regressed. +# (Hermite uses a separate `_pack_and_extend` path — pinned in its own testitem above.) +@testitem "ND :exclusive one-shot narrow-float (Float32) Int-vec grid — value-matches persistent" begin + using FastInterpolations + const FI = FastInterpolations + + x = collect(0:2) + y = collect(0:2) # Int Vector grids + data = Float32[1 2 3; 4 5 6; 7 8 9] # narrow float → value-matched Tg = Float32 + bc = FI.PeriodicBC(endpoint = :exclusive, period = 2.5) + xf = Float32.(x) + yf = Float32.(y) + + # Arithmetic kernels: Int grid + Float32 data → Float32 one-shot, equal to the persistent + # interpolant AND to the Float32-grid one-shot (the `_ExclusivePeriodicAxis` now follows Tg). + @testset "$name value-matches Float32" for (name, ctor) in + (("linear", FI.linear_interp), ("cubic", FI.cubic_interp)) + itp = ctor((x, y), data; bc = bc) # persistent reference → Float32 + @test itp(0.5f0, 0.5f0) isa Float32 + for q in ((0.5f0, 0.5f0), (1.5f0, 0.5f0), (2.25f0, 1.75f0), (0.0f0, 2.4f0)) + r = ctor((x, y), data, q; bc = bc) # pre-fix: TypeError (Float64 axis vs Float32 witness) + @test r isa Float32 # value-matched output type + @test isapprox(r, itp(q...); atol = 1.0f-5) # == persistent + @test isapprox(r, ctor((xf, yf), data, q; bc = bc); atol = 1.0f-5) # == Float32 grid + end + end + + # Selection kernel (constant): the period floats Int→Float64, so its exclusive axis (and output) + # are Float64 for BOTH one-shot and persistent — a deliberate, consistent contract (NOT the + # arithmetic value-match). Pin the consistency so a future change can't silently split them. + @testset "constant stays consistent (Int→Float64, one-shot ≡ persistent)" begin + itpc = FI.constant_interp((x, y), data; bc = bc) + for q in ((0.5f0, 0.5f0), (2.25f0, 1.75f0), (0.0f0, 2.4f0)) + rc = FI.constant_interp((x, y), data, q; bc = bc) + @test typeof(rc) === typeof(itpc(q...)) # both Float64 — must not diverge + @test rc == itpc(q...) # selection ⇒ exact match + end + end +end diff --git a/test/test_nd_oneshot_onthefly.jl b/test/test_nd_oneshot_onthefly.jl index 02eca9772..b43c37b15 100644 --- a/test/test_nd_oneshot_onthefly.jl +++ b/test/test_nd_oneshot_onthefly.jl @@ -133,18 +133,24 @@ end @test val_auto ≈ val_pre rtol = 1.0e-10 end - @testset "AutoCoeffs: quadratic AD seed matches PreCompute exactly" begin + @testset "AutoCoeffs: quadratic scalar IS the OnTheFly path (seed identity)" begin xq = range(0.0, 2.0, 15) yq = range(0.0, 2.0, 15) data_q = [sin(xi) * cos(yj) for xi in xq, yj in yq] q = (0.7, 0.9) - @test quadratic_interp((xq, yq), data_q, q) == - quadratic_interp((xq, yq), data_q, q; coeffs = PreCompute()) - @test quadratic_interp((xq, yq), data_q, q; deriv = (DerivOp(1), EvalValue())) == - quadratic_interp((xq, yq), data_q, q; deriv = (DerivOp(1), EvalValue()), coeffs = PreCompute()) - @test quadratic_interp((xq, yq), data_q, q; deriv = (EvalValue(), DerivOp(1))) == - quadratic_interp((xq, yq), data_q, q; deriv = (EvalValue(), DerivOp(1)), coeffs = PreCompute()) + # Default == explicit OnTheFly bit-exactly (same code path); PreCompute + # tracks within FP-reordering noise (bit-exact parity is opt-in via + # coeffs=PreCompute(), pinned in test_nd_raw_grid_oneshot.jl). + for deriv in ( + (EvalValue(), EvalValue()), + (DerivOp(1), EvalValue()), + (EvalValue(), DerivOp(1)), + ) + val_auto = quadratic_interp((xq, yq), data_q, q; deriv) + @test val_auto == quadratic_interp((xq, yq), data_q, q; deriv, coeffs = OnTheFly()) + @test val_auto ≈ quadratic_interp((xq, yq), data_q, q; deriv, coeffs = PreCompute()) rtol = 1.0e-12 + end end @testset "AutoCoeffs: interp scalar matches PreCompute" begin diff --git a/test/test_nd_output_type_promotion.jl b/test/test_nd_output_type_promotion.jl new file mode 100644 index 000000000..43dd6ab96 --- /dev/null +++ b/test/test_nd_output_type_promotion.jl @@ -0,0 +1,169 @@ +# ND output-type contract: the result eltype must follow the natural +# `promote_type(grid, data, query)` and be `@inferred`-stable, for every method × +# grid container × (data, query) float combination. An Int/OneTo grid value-matches +# the data float (Int grid + Float32 data → Float32 grid), so e.g. +# F32 grid+data+query → Float32, while a Float64 query still promotes → Float64. + +@testitem "ND persistent eval output type = promote_type(grid, data, query), @inferred" begin + FI = FastInterpolations + + builders = ( + ("linear", (g, d) -> linear_interp(g, d)), + ("cubic", (g, d) -> cubic_interp(g, d)), + ("quadratic", (g, d) -> quadratic_interp(g, d)), + ("constant", (g, d) -> constant_interp(g, d)), + ) + gridspecs = ( + ("F32vec", Float32, () -> (collect(Float32, 1:7), collect(Float32, 1:7))), + ("F64vec", Float64, () -> (collect(Float64, 1:7), collect(Float64, 1:7))), + ("IntOneTo", :match, () -> (Base.OneTo(7), Base.OneTo(7))), + ("IntVec", :match, () -> (collect(1:7), collect(1:7))), + ) + + for (mname, build) in builders, (gname, gTg, gbuild) in gridspecs, + Tv in (Float32, Float64), Tq in (Float32, Float64) + + gx, gy = gbuild() + xs = gx isa Base.OneTo ? (1:7) : gx + ys = gy isa Base.OneTo ? (1:7) : gy + data = Tv.([sin(3 * float(x)) + cos(2 * float(y)) for x in xs, y in ys]) + q = (Tq(2.4), Tq(3.6)) + Tg = gTg === :match ? promote_type(Float32, Tv) : gTg # Int/OneTo grid value-matches data float + expected = promote_type(Tg, Tv, Tq) + itp = build((gx, gy), data) + + @testset "$mname $gname Tv=$Tv Tq=$Tq → $expected" begin + @test itp(q) isa expected + @test (@inferred itp(q)) isa expected + end + end +end + +@testitem "ND one-shot output type = promote_type(grid, data, query), @inferred" begin + FI = FastInterpolations + + methods = ( + ("linear", (g, d, q) -> linear_interp(g, d, q)), + ("cubic", (g, d, q) -> cubic_interp(g, d, q)), + ("quadratic", (g, d, q) -> quadratic_interp(g, d, q)), + ("constant", (g, d, q) -> constant_interp(g, d, q)), + ) + gridspecs = ( + ("F32vec", Float32, () -> (collect(Float32, 1:7), collect(Float32, 1:7))), + ("F64vec", Float64, () -> (collect(Float64, 1:7), collect(Float64, 1:7))), + ("IntOneTo", :match, () -> (Base.OneTo(7), Base.OneTo(7))), + ("IntVec", :match, () -> (collect(1:7), collect(1:7))), + ) + + for (mname, mfn) in methods, (gname, gTg, gbuild) in gridspecs, + Tv in (Float32, Float64), Tq in (Float32, Float64) + + gx, gy = gbuild() + xs = gx isa Base.OneTo ? (1:7) : gx + ys = gy isa Base.OneTo ? (1:7) : gy + data = Tv.([sin(3 * float(x)) + cos(2 * float(y)) for x in xs, y in ys]) + q = (Tq(2.4), Tq(3.6)) + Tg = gTg === :match ? promote_type(Float32, Tv) : gTg # Int/OneTo grid value-matches data float + expected = promote_type(Tg, Tv, Tq) + + @testset "$mname $gname Tv=$Tv Tq=$Tq → $expected" begin + @test mfn((gx, gy), data, q) isa expected + @test (@inferred mfn((gx, gy), data, q)) isa expected + end + end +end + +@testitem "Hermite ND output type = promote_type(grid, data∪partials, query), @inferred" begin + FI = FastInterpolations + + gridspecs = ( + ("F32vec", Float32, () -> (collect(Float32, 1:5), collect(Float32, 1:5))), + ("F64vec", Float64, () -> (collect(Float64, 1:5), collect(Float64, 1:5))), + ("IntOneTo", :match, () -> (Base.OneTo(5), Base.OneTo(5))), + ("IntVec", :match, () -> (collect(1:5), collect(1:5))), + ) + + for (gname, gTg, gbuild) in gridspecs, Tv in (Float32, Float64), Tq in (Float32, Float64) + gx, gy = gbuild() + xs, ys = 1:5, 1:5 + data = Tv.([sin(0.5a) * cos(0.5b) for a in xs, b in ys]) + p = HermitePartials( + (1, 0) => Tv.([0.5cos(0.5a) * cos(0.5b) for a in xs, b in ys]), + (0, 1) => Tv.([-0.5sin(0.5a) * sin(0.5b) for a in xs, b in ys]), + (1, 1) => Tv.([-0.25cos(0.5a) * sin(0.5b) for a in xs, b in ys]), + ) + q = (Tq(2.4), Tq(3.6)) + Tg = gTg === :match ? promote_type(Float32, Tv) : gTg # Int/OneTo grid value-matches data float + expected = promote_type(Tg, Tv, Tq) + + @testset "oneshot $gname Tv=$Tv Tq=$Tq → $expected" begin + @test hermite_interp((gx, gy), data, p, q) isa expected + @test (@inferred hermite_interp((gx, gy), data, p, q)) isa expected + end + itp = hermite_interp((gx, gy), data, p) + @testset "persistent $gname Tv=$Tv Tq=$Tq → $expected" begin + @test itp(q) isa expected + @test (@inferred itp(q)) isa expected + end + end +end + +@testitem "cubic explicit PreCompute value-matches all axis containers" begin + # The scalar PreCompute backend promotes grids exactly like the batch path + # (`_nd_promote_grids`): Ranges → isbits `_CachedRange{Tg}` (spline caches memoise + # via value-deterministic objectid), same-eltype Vectors pass by identity, and a + # mismatched Vector converts (the caches match it by content). So Int/OneTo/Vector + # grids + Float32 data all return Float32, matching the OnTheFly default. + g = Base.OneTo(7) + data = Float32.([sin(3.0x) + cos(2.0y) for x in 1:7, y in 1:7]) + q = (2.4f0, 3.6f0) + v_otf = cubic_interp((g, g), data, q) + @test v_otf isa Float32 # default (OnTheFly) + v_pc = cubic_interp((g, g), data, q; coeffs = PreCompute()) + @test v_pc isa Float32 + @test v_pc ≈ v_otf rtol = 1.0f-5 + + gv = collect(1:7) + @test cubic_interp((gv, gv), data, q; coeffs = PreCompute()) isa Float32 + @test cubic_interp((g, gv), data, q; coeffs = PreCompute()) isa Float32 # mixed containers +end + +@testitem "ND Int-data output types: arithmetic float-forces, selection stays natural" begin + FI = FastInterpolations + dataI = [x + 2y for x in 1:7, y in 1:7] + + gridspecs = ( + ("IntOneTo", Int, () -> (Base.OneTo(7), Base.OneTo(7))), + ("IntVec", Int, () -> (collect(1:7), collect(1:7))), + ("F32vec", Float32, () -> (collect(Float32, 1:7), collect(Float32, 1:7))), + ("F64vec", Float64, () -> (collect(Float64, 1:7), collect(Float64, 1:7))), + ) + arith = ( + ("linear", linear_interp), + ("cubic", cubic_interp), + ("quadratic", quadratic_interp), + ) + + for (gname, Te, gbuild) in gridspecs, q in ((2.4f0, 3.6f0), (2.4, 3.6), (2, 3)) + Tq = typeof(q[1]) + gx, gy = gbuild() + + # Arithmetic kernels divide → Int floats: Tg = float(promote(grid, Int)). + Tg = float(promote_type(Te, Int)) + expected = promote_type(Tg, Tq) + for (mname, f) in arith + @testset "$mname $gname Int data Tq=$Tq → $expected" begin + @test f((gx, gy), dataI, q) isa expected + @test f((gx, gy), dataI)(q) isa expected # one-shot ≡ persistent type + end + end + + # Selection kernel (constant): no x·y arithmetic → pure natural promotion, + # NO float forcing (all-Int in → Int out); one-shot must match persistent. + expected_sel = promote_type(Te, Int, Tq) + @testset "constant $gname Int data Tq=$Tq → $expected_sel" begin + @test constant_interp((gx, gy), dataI, q) isa expected_sel + @test constant_interp((gx, gy), dataI)(q) isa expected_sel + end + end +end diff --git a/test/test_nd_raw_grid_oneshot.jl b/test/test_nd_raw_grid_oneshot.jl index 24b596f2d..70cc399fb 100644 --- a/test/test_nd_raw_grid_oneshot.jl +++ b/test/test_nd_raw_grid_oneshot.jl @@ -159,9 +159,10 @@ end end @allocated cubic_interp((x, y), data, q) end - # PreCompute (the practically-used cubic strategy) is now also raw: each axis's - # spline cache memoises by id and the cell width is promoted, so warm Int-grid - # one-shot is zero-alloc here too — no internal `Tg.(x)` convert. + # PreCompute (the practically-used cubic strategy): axes go through the POOLED + # value-matched wrap (`_cache_axes_pooled` — acquire + copyto!, no heap `Tg.(x)`), + # and the spline autocache hits via its content-match pass (pooled buffers don't + # keep a stable objectid) — so the warm Int-grid one-shot stays zero-alloc. function _alloc_cubic_nd_int_pre_2d() x = [0, 1, 2, 3, 4, 5, 6, 7] y = [0, 1, 2, 3, 4, 5] @@ -241,66 +242,149 @@ end end end +# ============================================================================ +# Int Vector grids + Float32 data — the value-matched (Tg = Float32) arms +# ============================================================================ +# +# `float(Int) === Float64`, so the Int-axis + Float64-data pins above never +# widen (axes stay raw). Float32 data forces `Tg = Float32`, and every arm is +# now zero-heap: PreCompute / linear go through POOLED wraps, the OnTheFly cubic +# collapse value-matches the raw Int axes via the data-aware cache, and the +# MIXED-method hetero OnTheFly fallback passes raw grids straight through +# (type-only promotion — no eager `_convert_grids_typed`), so the inner 1D +# one-shots value-match each axis themselves. + +@testitem "Int Vector grids + Float32 data (Tg = Float32 arms)" setup = [AllocConstants] begin + function _alloc_linear_int_f32_2d() + x = [0, 1, 2, 3, 4, 5, 6, 7] + y = [0, 1, 2, 3, 4, 5] + data = Float32[1.0f0 * a + 2.0f0 * b for a in x, b in y] + q = (3.4f0, 2.6f0) + for _ in 1:3 + linear_interp((x, y), data, q) + end + @allocated linear_interp((x, y), data, q) + end + function _alloc_cubic_int_pre_f32_2d() + x = [0, 1, 2, 3, 4, 5, 6, 7] + y = [0, 1, 2, 3, 4, 5] + data = Float32[sin(1.0f0 * a) + cos(1.0f0 * b) for a in x, b in y] + q = (3.4f0, 2.6f0) + for _ in 1:3 + cubic_interp((x, y), data, q; coeffs = PreCompute()) + end + @allocated cubic_interp((x, y), data, q; coeffs = PreCompute()) + end + function _alloc_cubic_int_otf_f32_2d() + x = [0, 1, 2, 3, 4, 5, 6, 7] + y = [0, 1, 2, 3, 4, 5] + data = Float32[sin(1.0f0 * a) + cos(1.0f0 * b) for a in x, b in y] + q = (3.4f0, 2.6f0) + for _ in 1:3 + cubic_interp((x, y), data, q) + end + @allocated cubic_interp((x, y), data, q) + end + function _alloc_hetero_mixed_int_f32_2d() + x = [0, 1, 2, 3, 4, 5, 6, 7] + y = [0, 1, 2, 3, 4, 5] + data = Float32[sin(1.0f0 * a) + 2.0f0 * b for a in x, b in y] + q = (3.4f0, 2.6f0) + m = (CubicInterp(), LinearInterp()) + for _ in 1:3 + interp((x, y), data, q; method = m) + end + @allocated interp((x, y), data, q; method = m) + end + + # value-match sanity: Int axes + Float32 data → Float32 out on every family + @testset "results are Float32" begin + x = [0, 1, 2, 3, 4] + y = [0, 1, 2, 3] + data = Float32[1.0f0 * a + 2.0f0 * b for a in x, b in y] + q = (1.4f0, 0.6f0) + @test linear_interp((x, y), data, q) isa Float32 + @test cubic_interp((x, y), data, q) isa Float32 + @test cubic_interp((x, y), data, q; coeffs = PreCompute()) isa Float32 + @test interp((x, y), data, q; method = (CubicInterp(), LinearInterp())) isa Float32 + end + + # Every arm is zero-alloc: pooled PreCompute/linear wraps, the OnTheFly cubic + # collapse (inner 1D value-matches raw Int axes via the data-aware cache), and + # the mixed-method hetero OnTheFly fallback (raw grids passed straight through + # by type-only promotion — no eager `_convert_grids_typed`). + @testset "warm zero-alloc (all raw-grid arms)" begin + @test _alloc_linear_int_f32_2d() <= ND_ALLOC_THRESHOLD + @test _alloc_cubic_int_pre_f32_2d() <= ND_ALLOC_THRESHOLD + @test _alloc_cubic_int_otf_f32_2d() <= ND_ALLOC_THRESHOLD + @test _alloc_hetero_mixed_int_f32_2d() <= ND_ALLOC_THRESHOLD + end +end + # ============================================================================ # Quadratic ND one-shot — raw Int grid, default PreCompute path # ============================================================================ # -# Quad `AutoCoeffs` defaults to PreCompute (for bit-exact AD-rule matching), so the -# DEFAULT scalar path is the PreCompute cell-eval. Passing raw grids + floating the -# cell width in `_compute_all_local_params` makes the default warm scalar one-shot -# on an Int Vector grid zero-alloc (the kernel wraps axes via `_cache_axis_pooled` -# into pool buffers, so the Int grid needs no `Tg.(x)` copy). Batch keeps eager-convert. +# Quad scalar `AutoCoeffs` routes through the shared `_resolve_coeffs_nd_oneshot` +# policy (scalar → OnTheFly, batch → PreCompute), same as cubic. Both scalar +# strategies stay zero-alloc on raw Int Vector grids: OnTheFly via pooled collapse, +# explicit PreCompute via pooled value-matched axis conversion. Batch keeps eager-convert. -@testitem "Quadratic ND one-shot raw-grid (default PreCompute, no eager convert)" setup = [AllocConstants] begin +@testitem "Quadratic ND one-shot raw-grid (no eager convert)" setup = [AllocConstants] begin using ForwardDiff - function _alloc_quad_nd_int_2d() + function _alloc_quad_nd_int_2d(coeffs) x = [0, 1, 2, 3, 4, 5, 6] y = [0, 1, 2, 3, 4] data = [sin(1.0 * a) + cos(1.0 * b) for a in x, b in y] q = (3.4, 2.6) for _ in 1:3 - quadratic_interp((x, y), data, q) + quadratic_interp((x, y), data, q; coeffs) end - @allocated quadratic_interp((x, y), data, q) + @allocated quadratic_interp((x, y), data, q; coeffs) end - function _alloc_quad_nd_int_3d() + function _alloc_quad_nd_int_3d(coeffs) x = [0, 1, 2, 3, 4] y = [0, 1, 2, 3] z = [0, 1, 2, 3, 4, 5] data = [a + 0.5b + 0.25c for a in x, b in y, c in z] q = (2.4, 1.6, 3.8) for _ in 1:3 - quadratic_interp((x, y, z), data, q) + quadratic_interp((x, y, z), data, q; coeffs) end - @allocated quadratic_interp((x, y, z), data, q) + @allocated quadratic_interp((x, y, z), data, q; coeffs) end - @testset "warm zero-alloc scalar one-shot on Int Vector grids (default PreCompute)" begin - @test _alloc_quad_nd_int_2d() <= ND_ALLOC_THRESHOLD - @test _alloc_quad_nd_int_3d() <= ND_ALLOC_THRESHOLD + @testset "warm zero-alloc scalar one-shot on Int Vector grids ($label)" for (label, coeffs) in + (("default → OnTheFly", AutoCoeffs()), ("explicit PreCompute", PreCompute())) + @test _alloc_quad_nd_int_2d(coeffs) <= ND_ALLOC_THRESHOLD + @test _alloc_quad_nd_int_3d(coeffs) <= ND_ALLOC_THRESHOLD end - # ---- bit-identical: raw Int grid === Float64 grid (default PreCompute) ---- - @testset "Int Vector grid === Float64 grid" begin + # ---- bit-identical: raw Int grid === Float64 grid (both scalar strategies) ---- + @testset "Int Vector grid === Float64 grid ($label)" for (label, coeffs) in + (("default → OnTheFly", AutoCoeffs()), ("explicit PreCompute", PreCompute())) x = [0, 1, 2, 3, 4, 5, 6] y = [0, 1, 2, 3, 4] data = [sin(1.0 * a) + cos(1.0 * b) for a in x, b in y] xf = Float64.(x) yf = Float64.(y) for q in [(3.4, 2.6), (0.3, 0.7), (5.9, 3.1)] - @test quadratic_interp((x, y), data, q) === quadratic_interp((xf, yf), data, q) + @test quadratic_interp((x, y), data, q; coeffs) === + quadratic_interp((xf, yf), data, q; coeffs) end end - # ---- one-shot === persistent QuadraticInterpolantND (both PreCompute) ---- - @testset "one-shot === persistent QuadraticInterpolantND" begin + # ---- persistent parity: bit-exact is opt-in via coeffs=PreCompute(); + # the default (OnTheFly) matches to FP-reordering noise ---- + @testset "one-shot vs persistent QuadraticInterpolantND" begin x = [0, 1, 2, 3, 4, 5, 6] y = [0, 1, 2, 3, 4] data = [sin(1.0 * a) + cos(1.0 * b) for a in x, b in y] itp = quadratic_interp((x, y), data) for q in [(3.4, 2.6), (0.3, 0.7), (5.9, 3.1)] - @test quadratic_interp((x, y), data, q) === itp(q) + @test quadratic_interp((x, y), data, q; coeffs = PreCompute()) === itp(q) + @test quadratic_interp((x, y), data, q) ≈ itp(q) rtol = 1.0e-12 end end diff --git a/test/test_nd_store_copyfalse_narrow_float.jl b/test/test_nd_store_copyfalse_narrow_float.jl new file mode 100644 index 000000000..c27f77015 --- /dev/null +++ b/test/test_nd_store_copyfalse_narrow_float.jl @@ -0,0 +1,46 @@ +# `StorePolicy(copy = false)` must be a true zero-copy contract for narrow-float ND data +# (Float32/Float16/ComplexF32): the stored array IS the caller's array (`itp.data === A`), +# never widened. Pins the fix for the ND builders' grid-eltype-only promotion, which ran +# `Tv.(data)` (an O(n²) Float64 copy) before the store policy was consulted. + +@testitem "ND store copy=false is zero-copy for narrow-float data" begin + FI = FastInterpolations + build(A) = FI.linear_interp( + axes(A), A; + extrap = FI.ClampExtrap(), + store = FI.StorePolicy(copy = false), + ) + + # The broken set: bare numerics narrower than Float64 (+ Complex), which the + # legacy promotion widens. `itp.data === A` is the exact zero-copy witness. + @testset "$T: copy=false aliases the input (no preemptive widening)" for T in + (Float32, Float16, ComplexF32) + A = rand(T, 8, 8) + itp = build(A) + @test eltype(itp.data) === T # value array not widened to Float64/ComplexF64 + @test itp.data === A # stored array IS the caller's array → zero copy + end + + # Invariant guards — these already alias today; they must not regress under the fix. + @testset "invariant: already-aliasing types stay aliased" begin + for A in (rand(Float64, 8, 8), rand(ComplexF64, 8, 8)) + @test build(A).data === A + end + end +end + +@testitem "ND store copy=false: narrow-float build allocates no O(n²) copy" begin + FI = FastInterpolations + _build_alloc(A) = @allocated FI.linear_interp( + axes(A), A; + extrap = FI.ClampExtrap(), + store = FI.StorePolicy(copy = false), + ) + @testset "$T build is not a full-array copy" for T in (Float32, Float16, ComplexF32) + A = rand(T, 256, 256) + _build_alloc(A) # warm up (compilation) + data_bytes = sizeof(T) * length(A) + # An alias allocates ~0; a widened copy allocates ≥ one data array. + @test _build_alloc(A) < data_bytes ÷ 4 + end +end