(fix + refac): Method-aware output-eltype via kernel-shape op trait

src/akima/akima_oneshot.jl

@@ -177,7 +177,7 @@ function akima_interp(
         search::AbstractSearchPolicy = AutoSearch(),
         hint::Union{Nothing, Base.RefValue{Int}} = nothing
     ) where {Tg, Tv, Tq <: Real}
-    Tr = _output_eltype(Tv, _promote_grid_float(Tg, Tv), Tq)
+    Tr = _output_eltype(_arithmetic_kernel_shape, _promote_grid_float(Tg, Tv), Tv, Tq)
     output = Vector{Tr}(undef, length(x_query))
     akima_interp!(output, x, y, x_query; bc = bc, coeffs = coeffs, extrap = extrap, deriv = deriv, search = search, hint = hint)
     return output

src/cardinal/cardinal_oneshot.jl

@@ -199,7 +199,7 @@ function cardinal_interp(
         search::AbstractSearchPolicy = AutoSearch(),
         hint::Union{Nothing, Base.RefValue{Int}} = nothing
     ) where {Tg, Tv, Tq <: Real}
-    Tr = _output_eltype(Tv, _promote_grid_float(Tg, Tv), Tq)
+    Tr = _output_eltype(_arithmetic_kernel_shape, _promote_grid_float(Tg, Tv), Tv, Tq)
     output = Vector{Tr}(undef, length(x_query))
     cardinal_interp!(output, x, y, x_query; bc = bc, coeffs = coeffs, tension = tension, extrap = extrap, deriv = deriv, search = search, hint = hint)
     return output

src/constant/constant_adjoint.jl

@@ -229,8 +229,8 @@ function constant_adjoint(
         side::AbstractSide = NearestSide(),
         extrap::AbstractExtrap = NoExtrap(),
     ) where {Tg}
-    # Selection kernel → raw eltype contract (cubic/linear/… keep using the
-    # shared `_promote_adjoint_inputs` for their Float-promotion needs).
+    # Grid stays raw `Tg` (no `_promote_adjoint_inputs` Float widening).
+    # Adjoint buffer eltype comes from the protocol's `_output_eltype`.
     x_p = x
     xq_p = _promote_query_typed(x_query, Tg)
 

src/constant/constant_anchor.jl

@@ -273,14 +273,14 @@ end
 # Canonical evaluation functions that take raw y vector + explicit params.
 # Used by both interpolant anchor dispatch AND series one-shot evaluation.
 
-# Default case (extension, wrap, inbounds): kernel with right-boundary check
+# Default case (extension, wrap, inbounds): kernel with right-boundary check.
+# Short-circuits use `* one(aq.xq)` to match the kernel's `* one(dL)` carrier
+# propagation — without it, Int y + Float xq returns Union{Int,Float}.
 @inline function _constant_eval_at_anchor(
         y::AbstractVector, x_last, aq::_ConstantAnchoredQuery,
         op::AbstractEvalOp, side_param::AbstractSide, ::AbstractExtrap
     )
-    # Right-edge short-circuit: `y[aq.idxR]` resolves to `y[end]` for non-periodic
-    # (idxR == n) and to the cyclic `y[1]` for `_ExclusivePeriodicAxis` (idxR == 1).
-    aq.xq == x_last && return (op isa EvalValue ? (@inbounds y[aq.idxR]) : 0 * first(y))
+    aq.xq == x_last && return (op isa EvalValue ? (@inbounds y[aq.idxR] * one(aq.xq)) : 0 * first(y) * one(aq.xq))
     @inbounds return _constant_kernel(op, y[aq.idxL], y[aq.idxR], aq.h, aq.dL, side_param)
 end
 
@@ -290,9 +290,7 @@ end
         op::AbstractEvalOp, side_param::AbstractSide, ::NoExtrap
     )
     aq.state != IN_DOMAIN && throw(DomainError(aq.xq, "query point outside domain"))
-    # Right-edge short-circuit: `y[aq.idxR]` resolves to `y[end]` for non-periodic
-    # (idxR == n) and to the cyclic `y[1]` for `_ExclusivePeriodicAxis` (idxR == 1).
-    aq.xq == x_last && return (op isa EvalValue ? (@inbounds y[aq.idxR]) : 0 * first(y))
+    aq.xq == x_last && return (op isa EvalValue ? (@inbounds y[aq.idxR] * one(aq.xq)) : 0 * first(y) * one(aq.xq))
     @inbounds return _constant_kernel(op, y[aq.idxL], y[aq.idxR], aq.h, aq.dL, side_param)
 end
 
@@ -305,9 +303,7 @@ end
         y_bnd = aq.state == OOB_LEFT ? first(y) : last(y)
         return _eval_extrapolation(op, y_bnd, extrap, aq.xq)
     end
-    # Right-edge short-circuit: `y[aq.idxR]` resolves to `y[end]` for non-periodic
-    # (idxR == n) and to the cyclic `y[1]` for `_ExclusivePeriodicAxis` (idxR == 1).
-    aq.xq == x_last && return (op isa EvalValue ? (@inbounds y[aq.idxR]) : 0 * first(y))
+    aq.xq == x_last && return (op isa EvalValue ? (@inbounds y[aq.idxR] * one(aq.xq)) : 0 * first(y) * one(aq.xq))
     @inbounds return _constant_kernel(op, y[aq.idxL], y[aq.idxR], aq.h, aq.dL, side_param)
 end
 
@@ -334,10 +330,8 @@ end
         x_min, x_max = first(itp.x), last(itp.x)
         throw(DomainError(aq.xq, "query point outside domain [$x_min, $x_max]"))
     end
-    # Right-edge short-circuit: `idxR` resolves to `n` for non-periodic and to
-    # `n+1` (cyclic `y[1]`) on `:exclusive` PeriodicBC's extended-grid persistent path.
     if aq.xq == last(itp.x)
-        return op isa EvalValue ? (@inbounds itp.y[aq.idxR]) : zero(T)
+        return op isa EvalValue ? (@inbounds itp.y[aq.idxR] * one(aq.xq)) : zero(T) * one(aq.xq)
     end
     @inbounds return _constant_kernel(op, itp.y[aq.idxL], itp.y[aq.idxR], aq.h, aq.dL, itp.side)
 end
@@ -362,11 +356,12 @@ end
 Evaluate constant interpolant at multiple anchored query points.
 Returns newly allocated vector.
 """
-function (itp::ConstantInterpolant{T})(
-        aq_vec::AbstractVector{<:_ConstantAnchoredQuery{T}};
+function (itp::ConstantInterpolant{Tg, Tv})(
+        aq_vec::AbstractVector{<:_ConstantAnchoredQuery{Tg, Tq}};
         deriv::DerivOp = EvalValue()
-    ) where {T}
-    output = Vector{T}(undef, length(aq_vec))
+    ) where {Tg, Tv, Tq <: Real}
+    T_out = _output_eltype(_constant_kernel_shape, Tg, Tv, Tq)
+    output = Vector{T_out}(undef, length(aq_vec))
     @inbounds for i in eachindex(aq_vec)
         output[i] = _constant_eval_with_anchor(itp, aq_vec[i], deriv)
     end

src/constant/constant_interpolant.jl

@@ -20,15 +20,15 @@ end
     return _constant_vector_loop!(output, itp.x, itp.y, xq, extrap, itp.side, op, searcher)
 end
 
-# ========================================
-# Selection-kernel output eltype trait
-# ========================================
-# Override the shared `_output_eltype(itp, Tq)` trait so the protocol's scalar
-# + batch callables keep raw `Tv` for plain numeric queries (selection kernel)
-# and only widen for duck-typed queries (Dual, …) so AD carriers round-trip.
-# Single source of truth — no callable overrides needed.
+# Constant declares its kernel shape — selection (`y * one(dL)`, no division).
+# Args mirror the real kernel: `(xL, yv, xq)` so `xq - xL` exposes the actual
+# `dL` carrier (e.g. `Dual` grid + `Float` xq → `Dual` dL). Trait infers the
+# exact return type via `promote_op`, so scalar/batch agree (Int×Int×Int → Int;
+# SVector × Dual → SVector{Dual}; Float y × Dual grid → Dual; etc.).
+@inline _constant_kernel_shape(xL, yv, xq) = yv * one(xq - xL)
+
 @inline _output_eltype(::ConstantInterpolant{Tg, Tv}, ::Type{Tq}) where {Tg, Tv, Tq} =
-    Tq <: _PromotableValue ? Tv : promote_type(Tv, Tq)
+    _output_eltype(_constant_kernel_shape, Tg, Tv, Tq)
 
 # ─────────────────────────────────────────────────────────────
 # Vector loop (function barrier)
@@ -110,8 +110,9 @@ end
 # ========================================
 # Generic Constructor (User API)
 # ========================================
-# Selection kernel → raw eltype contract (Int in → Int out). Signature
-# parametrized directly on `{Tg, Tv}` — no `_promote_grid_float` indirection.
+# Storage parametrized on `{Tg, Tv}` directly — no `_promote_grid_float`
+# indirection. Return type widens via the kernel's `* one(dL)` carrier
+# propagation (handled per-callable, not at construction).
 @inline function constant_interp(
         x::AbstractVector{Tg},
         y::AbstractVector{Tv};

src/constant/constant_kernels.jl

@@ -17,7 +17,11 @@
 # AD Support:
 # - dL can be ForwardDiff.Dual for automatic differentiation
 # - Comparisons use _extract_primal(dL) to get Float value
-# - Output is always Tv (no AD propagation through constant interp - derivative is 0)
+# - Multiplying the selection by `one(dL)` propagates `Tq`'s carrier (Dual,
+#   Float, …) into the result while keeping the value unchanged — aligns
+#   in-domain type with `_promote_extrap_val`'s OOB output. Derivative
+#   branches use `0 * y_left * one(dL)` so they only require `*(::Int, ::Tv)`
+#   and `*(::Tv, ::Real)` (no `zero(::Tv)` assumption beyond master's).
 #
 # Grid point behavior: When dL == 0 (exactly at grid point),
 # all side modes return y_left (the value at that grid point).
@@ -30,7 +34,7 @@ Always returns the left boundary value `y_left`.
 dL can be any Real (including ForwardDiff.Dual for AD).
 """
 @inline function _constant_kernel(::EvalValue, y_left::Tv, ::Tv, ::Tg, dL::Td, ::LeftSide) where {Tv, Tg, Td <: Real}
-    return y_left
+    return y_left * one(dL)
 end
 
 """
@@ -43,7 +47,8 @@ dL can be any Real (including ForwardDiff.Dual for AD).
 @inline function _constant_kernel(::EvalValue, y_left::Tv, y_right::Tv, ::Tg, dL::Td, ::RightSide) where {Tv, Tg, Td <: Real}
     # Use primal value for comparison (supports ForwardDiff.Dual)
     dL_primal = _extract_primal(dL)
-    return iszero(dL_primal) ? y_left : y_right
+    selected = iszero(dL_primal) ? y_left : y_right
+    return selected * one(dL)
 end
 
 """
@@ -56,7 +61,8 @@ dL can be any Real (including ForwardDiff.Dual for AD).
 @inline function _constant_kernel(::EvalValue, y_left::Tv, y_right::Tv, h::Tg, dL::Td, ::NearestSide) where {Tv, Tg, Td <: Real}
     # Use primal value for comparison (supports ForwardDiff.Dual)
     dL_primal = _extract_primal(dL)
-    return dL_primal <= h / 2 ? y_left : y_right
+    selected = dL_primal <= h / 2 ? y_left : y_right
+    return selected * one(dL)
 end
 
 """
@@ -67,7 +73,7 @@ Always returns zero (constant function has no slope).
 Uses `0 * y_left` for duck-typing support and NaN propagation.
 """
 @inline function _constant_kernel(::EvalDeriv1, y_left::Tv, ::Tv, ::Tg, dL::Td, ::AbstractSide) where {Tv, Tg, Td <: Real}
-    return 0 * y_left
+    return 0 * y_left * one(dL)
 end
 
 """
@@ -77,7 +83,7 @@ Second derivative of constant interpolation.
 Always returns zero (constant function has no curvature).
 """
 @inline function _constant_kernel(::EvalDeriv2, y_left::Tv, ::Tv, ::Tg, dL::Td, ::AbstractSide) where {Tv, Tg, Td <: Real}
-    return 0 * y_left
+    return 0 * y_left * one(dL)
 end
 
 """
@@ -86,12 +92,12 @@ end
 Third derivative of constant interpolation is always zero.
 """
 @inline function _constant_kernel(::EvalDeriv3, y_left::Tv, ::Tv, ::Tg, dL::Td, ::AbstractSide) where {Tv, Tg, Td <: Real}
-    return 0 * y_left
+    return 0 * y_left * one(dL)
 end
 
 """Generic fallback: N-th derivative of degree-0 (constant) is zero for N ≥ 1."""
-@inline function _constant_kernel(::DerivOp{N}, y_left::Tv, ::Tv, ::Tg, ::Td, ::AbstractSide) where {N, Tv, Tg, Td <: Real}
-    return 0 * y_left
+@inline function _constant_kernel(::DerivOp{N}, y_left::Tv, ::Tv, ::Tg, dL::Td, ::AbstractSide) where {N, Tv, Tg, Td <: Real}
+    return 0 * y_left * one(dL)
 end
 
 # ========================================

src/constant/constant_oneshot.jl

@@ -30,12 +30,10 @@
         searcher::S
     ) where {Tg, Tv, Tq <: Real, S <: Searcher}
     if _extract_primal(xi) == _extract_primal(last(x))
-        # `last(x)` for `_ExclusivePeriodicAxis` is the *virtual* `inner[1] + period`
-        # (the seam right endpoint). The corresponding y at that virtual slot is
-        # `last(y) = inner[1]` for `_ExclusivePeriodicData` — cyclic via the data
-        # wrapper. For raw vectors both `last`s are the user's last entry.
-        # Single uniform `last(y)` handles both cases — no `_resolve_idx` needed.
-        return op isa EvalValue ? last(y) : 0 * first(y)
+        # `last(y)` covers both raw vectors and `_ExclusivePeriodicData` (cyclic
+        # `inner[1]`). `* one(xi)` propagates Tq carrier to match the kernel
+        # path below (without it, Int y + Float xq returns Union{Int,Float}).
+        return op isa EvalValue ? last(y) * one(xi) : 0 * first(y) * one(xi)
     end
     idx, idx_R, xL, xR = search_interval(searcher, x, xi)
     dL = xi - xL
@@ -108,7 +106,7 @@ end
     # the exact boundary. `last(_ExclusivePeriodicData) = inner[1]` so `:exclusive`
     # cyclic wrap is preserved; raw Vector yields `y[n]`.
     _extract_primal(xi_wrapped) == _extract_primal(last(x)) &&
-        return op isa EvalValue ? last(y) : 0 * first(y)
+        return op isa EvalValue ? last(y) * one(xi_wrapped) : 0 * first(y) * one(xi_wrapped)
     idx, idx_R, xL, xR = search_interval(searcher, x, xi_wrapped)
     dL = xi_wrapped - xL
     # Unwrap data once: `search_interval` already resolved the seam (idx_R = 1
@@ -157,8 +155,8 @@ Constant (step/piecewise constant) interpolation at a single point.
   - `LinearBinarySearch(linear_window=8)`: Linear search within window, then binary fallback
 
 # Returns
-- Interpolated value, eltype `eltype(y)` (raw Tv; widens to `promote_type(Tv, Tq)`
-  only for duck-typed queries).
+- Interpolated value, eltype `promote_type(eltype(y), eltype(xq))` (the
+  kernel's `* one(dL)` carrier propagation; fully-Int chain preserves Int).
 
 # Example
 ```julia
@@ -276,9 +274,10 @@ sorted_queries = sort(rand(1000))
 vals = constant_interp(x, y, sorted_queries; search=LinearBinarySearch(linear_window=8))
 ```
 """
-# Unified allocating vector one-shot. Output eltype = `eltype(y)` for plain
-# numeric queries (raw Tv contract); widens to `promote_type(Tv, Tq)` for
-# duck-typed queries (Dual, Measurement, …) so AD carriers aren't stripped.
+# Buffer eltype via Constant's kernel shape — Julia infers the return type
+# from `_constant_kernel_shape(xL, yv, xq) = yv * one(xq - xL)`, matching the
+# actual kernel reality (Int×Int×Int → Int; SVector × Dual → SVector{Dual};
+# Float y × Dual grid → Dual carrier via `xq - xL`).
 function constant_interp(
         x::AbstractVector,
         y::AbstractVector,
@@ -289,10 +288,9 @@ function constant_interp(
         deriv::DerivOp = EvalValue(),
         search::AbstractSearchPolicy = AutoSearch()
     )
-    Tv = eltype(y)
-    Tq = eltype(x_targets)
-    T_out = Tq <: _PromotableValue ? Tv : promote_type(Tv, Tq)
-    output = Vector{T_out}(undef, length(x_targets))
+    output = Vector{_output_eltype(_constant_kernel_shape, eltype(x), eltype(y), eltype(x_targets))}(
+        undef, length(x_targets)
+    )
     constant_interp!(output, x, y, x_targets; bc, extrap, side, deriv, search)
     return output
 end

src/constant/constant_oneshot_series.jl

@@ -78,9 +78,7 @@ end
     x = _to_float(x, Tg)
     K = n_series(s)
     Tv = _series_eltype(s)
-    # Duck-typed queries (Dual, …) widen to keep AD carrier; plain queries
-    # preserve raw Tv. Mirrors ConstantSeriesInterpolant scalar callable.
-    Tv_out = Tq <: _PromotableValue ? Tv : promote_type(Tv, Tq)
+    Tv_out = _output_eltype(_constant_kernel_shape, Tg, Tv, Tq)
     output = Vector{Tv_out}(undef, K)
     if _is_periodic_bc(bc)
         # Helper wraps `x` via `_resolve_axis(x, bc)` and searches against the
@@ -286,8 +284,7 @@ function constant_interp(
     ) where {Tg, Tq <: Real}
     K = n_series(s)
     Tv = _series_eltype(s)
-    # Same Tv_out rule as the scalar series oneshot — duck queries widen.
-    Tv_out = Tq <: _PromotableValue ? Tv : promote_type(Tv, Tq)
+    Tv_out = _output_eltype(_constant_kernel_shape, Tg, Tv, Tq)
     outputs = [Vector{Tv_out}(undef, length(xqs)) for _ in 1:K]
     constant_interp!(outputs, x, s, xqs; bc, side, extrap, deriv, search)
     return outputs

src/constant/constant_series_interp.jl

@@ -393,21 +393,20 @@ Returns a vector of values, one per y-series.
 - `deriv=DerivOp(1),DerivOp(2)`: Returns zeros (step function derivative is zero everywhere)
 
 # AD Support
-Plain queries (`Tq <: _PromotableValue`) → output eltype `Tv` unchanged.
-Duck-typed queries (e.g. `ForwardDiff.Dual`) widen to `promote_type(Tv, Tq)`.
+Output eltype routes through the kernel-shape trait
+`_output_eltype(_constant_kernel_shape, Tg, Tv, Tq)`. `Base.promote_op` on
+`yv * one(xq - xL)` jointly considers grid (Tg), value (Tv) and query (Tq),
+so duck carriers on either Tg (e.g. `Dual` grid + `Float` xq) or Tq (`Float`
+grid + `Dual` xq) propagate uniformly. Carriers like `SVector × Dual`
+resolve correctly instead of collapsing to `Vector{Any}`.
 """
 function (sitp::ConstantSeriesInterpolant{Tg, Tv, P})(
         xq::Tq;
         deriv::DerivOp = EvalValue(),
         search::AbstractSearchPolicy = sitp.search_policy,
         hint::Union{Nothing, Base.RefValue{Int}} = nothing
     ) where {Tg, Tv, P, Tq <: Real}
-    # Selection kernel returns `y[idx]::Tv` — Tv flows through unchanged for
-    # plain numerical queries. Duck-typed queries (Dual, …) get the legacy
-    # `promote_type(Tv, Tq)` widening so AD callers keep their input-Dual →
-    # output-Dual API contract (the partial wrt xq is 0 for constant, but the
-    # carrier type must match).
-    T_out = Tq <: _PromotableValue ? Tv : promote_type(Tv, Tq)
+    T_out = _output_eltype(_constant_kernel_shape, Tg, Tv, Tq)
     out = Vector{T_out}(undef, n_series(sitp))
     return sitp(out, xq; deriv = deriv, search = search, hint = hint)
 end
@@ -459,10 +458,10 @@ function (sitp::ConstantSeriesInterpolant{Tg, Tv, P})(
     n_query = length(xq_typed)
     n_ser = n_series(sitp)
 
-    # Explicit Vector{Vector{Tv}} for type stability on Julia LTS
-    outputs = Vector{Vector{Tv}}(undef, n_ser)
+    T_out = _output_eltype(_constant_kernel_shape, Tg, Tv, Tq)
+    outputs = Vector{Vector{T_out}}(undef, n_ser)
     @inbounds for k in 1:n_ser
-        outputs[k] = Vector{Tv}(undef, n_query)
+        outputs[k] = Vector{T_out}(undef, n_query)
     end
     sitp(outputs, xq_typed; deriv = deriv, search = search, hint = hint)
 
@@ -484,7 +483,7 @@ query on the stack and reused for all K series, staying in registers across
 the K evals. No pool, no `aq_vec` scratch.
 """
 function (sitp::ConstantSeriesInterpolant{Tg, Tv, P})(
-        outputs::AbstractVector{<:AbstractVector{Tv}},
+        outputs::AbstractVector{<:AbstractVector},
         xq::AbstractVector{<:Real};
         deriv::DerivOp = EvalValue(),
         search::AbstractSearchPolicy = sitp.search_policy,

src/constant/constant_types.jl

@@ -61,9 +61,8 @@ struct ConstantInterpolant{Tg, Tv, X <: AbstractVector{Tg}, Y <: AbstractVector{
     search_policy::P  # Default search policy (immutable, thread-safe)
 
     # Inner: `_cache_axis` (insurance) then `_convert_copy` for ownership.
-    # `{Tg, Tv}` parametrized directly — selection kernel → raw eltype
-    # contract (Int in → Int out), unlike arithmetic methods that thread
-    # `_promote_grid_float(TX, TY)` through here.
+    # `{Tg, Tv}` parametrized directly — no `_promote_grid_float(TX, TY)`
+    # indirection (storage stays raw; the kernel handles return-type widening).
     function ConstantInterpolant(
             x::AbstractVector{Tg}, y::AbstractVector{Tv}, ev::E, sv::SD, search::P;
             bc::AbstractBC = NoBC()

src/constant/nd/constant_nd_adjoint.jl

@@ -240,7 +240,8 @@ end
 @inline function _constant_adjoint_dispatch(
         grids::NTuple{N, AbstractVector}, queries, bc, side, extrap
     ) where {N}
-    # Selection-kernel adjoint: raw eltype contract (no `float()` widening).
+    # Grid stays raw (no `float()` widening); adjoint buffer eltype comes
+    # from the protocol's `_output_eltype`.
     Tg = _promote_grid_eltype(grids)
     grids_typed = _convert_grids_typed(grids, Tg)