feat(nd): unit axes for the ND solver families via dimensionless twins - #202
Merged
Conversation
…ire before float promotion
…— scaled [Y]-homogeneous store
…params, grid⁻ᵏ restoration at the cell seam
…ll params, shared reparam helpers hoisted to nd_utils
… engine run, volume-element restore
…e + composition-gap pins
The ND-internal fiber entry (_deriv_1d!) never runs the 1D builders'
payload normalize, so a Real Deriv2(0.0) beside unit ND data wrote 0.0
straight into the [Y]-typed RHS (DimensionError at _compute_rhs_first!).
_scale_bcs_reparam(bcs, grids, data) now threads a data sample so Real
structural payloads take the canonical _payload_val rehydration — on the
dimensionless t-axis the payload space IS [Y].
New testitem pins the composition surface: 3D mixed-unit (s x m x kg)
build/eval/deriv/integrate vs twin, per-axis PolyFit{4}, Real-zero
BCPair (=== against the ZeroCurvBC mint), exclusive-periodic unit axis,
and gradient! into an eltype-compatible store.
…the pool pipeline, op-aware Tr Cubic/quadratic ND one-shots (scalar, batch, in-place) now accept the same Real-or-unit axes as the persistent builders: - Gates unified: all five one-shot entries run the per-axis `_check_nd_reparam_grid` probe; `_check_nd_solver_grid` is deleted. - Pool backends (cubic/quadratic x scalar/batch) grow a Tg-folded twin arm — `_reparam_grids` + `_scale_bcs_reparam` feed the solve, the pooled partials stay [Y]-homogeneous, eval keeps reading the unit axes via `_compute_all_local_params_reparam`, and the cell seam restores grid^-k units through `_nd_fill_deriv_scale`. Real bodies are branch-preserved verbatim (type-folded conditions only). - `Tr` at the four public entries unified to linear's op-aware canonical (`_deriv_eltype_nd` over `_nd_value_eltype`) — unit derivative results and batch buffers are typed in value/coord^N space; identity on Real grids. - Batch backends drop the common-`Tg` grids pin (mixed-unit axes arrive as a heterogeneous per-axis-float tuple; the linear sibling relaxed first). - `_cache_axes_pooled` mirrors `_convert_cache_axes`' abstract-tag arm (mixed-unit axes wrap at their own eltype). - AutoCoeffs + non-Real scalar query resolves to the PreCompute pool by dispatch (the OnTheFly arm is Real-tuple-only); explicit OnTheFly rides the hetero collapse and gates via `_check_nd_hetero_grid`. Gotcha fixed along the way: `@with_pool`'s typed-checkpoint local scan only recognizes plain-symbol assignments, so an `acquire!` eltype bound by tuple destructuring is hoisted above its definition (runtime UndefVarError) — the twin arms use per-branch simple assignments. Tests: new P5b mirror testitem (scalar/batch/in-place vs persistent, deriv units, same-unit concrete-Tg, unit Range axis, exclusive periodic); F17 narrowed to the hetero families; F20 one-shot pins flipped to parity.
…h seams
Every Real-vs-reparam decision introduced on this branch now lives in a
dispatch arm instead of a call-site ternary ("dispatch, not a boolean test" —
the codebase's own gate idiom):
- `_reparam_solve_frame(grids, bcs, data)` — solve-frame selection for the
scaled-store families (persistent build x2 + one-shot pool x4); the Tz line
collapses to the build's canonical
`_promote_eltype(op, _promote_grid_eltype(grids_solve), Tv)` everywhere.
- `_compute_all_local_params` 5-arg is constrained to `Tg <: Real` (its body
is Real-width semantics) and gains a `::Type` fallback that routes to the
dimensionless collapse — the six local-params call sites call it
unconditionally with the width tag in scope.
- `_restore_nd_deriv_scale(r, grids, ops)` — the cell-seam grid^-k
restoration x6; the Real arm is an identity METHOD (branch preservation
stays load-bearing, no x1.0).
- `_integrate_reparam` ternary becomes tag-dispatch arms.
- `_oneshot_fill_eltype(Tg, Tv)` — cubic one-shot fill/extrap space.
Real arms are the exact previous computations (per-site type-identity
checked); the per-branch assignment blocks that worked around @with_pool's
destructuring blind spot disappear with the branches themselves.
…asing
Comment-only pass over the branch's additions: gate/build/eval/one-shot blocks
tightened to their load-bearing constraints, and two sentences left stale by
the dispatch refactor ("type-folded branch at the call sites", the inline
`true`-fold note) rewritten against the current seams.
…ector{BCPair} solve arms
Two sibling arms skipped the grid-aware `_normalize_bc(::BCPair, x, y)`
rehydration, so structural Real zeros beside unit data died in the solve
(`DimensionError: K m⁻¹ and 0.0 m`):
- The solve-side `_normalize_bc_array` overload for `Vector{<:BCPair}`
delegated to the type-only 3-arg path. Deleted — the general per-element
grid-aware arm handles `Vector{BCPair}` correctly as-is.
- The 3-arg catch-all dropped `(x, y)` for bare `PointBC`s, so scalar
`bc = Deriv2(0.0)` (plain and Series) never rehydrated. New one-line arm
routes through the pair rehydration: `_normalize_bc(BCPair(bc, bc), x, y)`.
Nonzero unitless payloads now get the same actionable ArgumentError on these
arms as on the scalar BCPair path (previously a deep DimensionError). Pins:
four success shapes === the ZeroCurvBC mint + two rejection-contract pins.
…on + FillExtrap surfaces F14 pinned only the Linear family. Extend it with the campaign's surfaces — cubic/quadratic persistent eval/deriv/integrate (full + bounded), one-shot scalar/PreCompute/deriv/batch, and FillExtrap in-domain/OOB/deriv-OOB — plus a 1D sibling pin for the rehydrated bare-PointBC / Series-array builds and unit FillExtrap. All 16 audited call shapes infer concretely with the expected unit-typed results today; these pins keep it that way.
…pace `_scale_payload_reparam`'s Real arm witnessed against the data sample alone ([Y]), so Real data on unit axes silently accepted a dimensionally incomplete nonzero payload (`Deriv1(0.25)` read as 0.25/axis-unit) that the 1D builders reject. The witness is now `first(data) * inv(u)` — the payload's true [Y/Xᵏ] space — rehydrated there and then scaled into the [Y] store. Structural zeros mint identically; typed payloads are untouched. Pins: nonzero Real payload beside Real data + unit axes → the 1D-actionable ArgumentError; zero-mint === parity; typed [Y/X] payloads stay accepted.
… parity)
Real axes float Int/Rational/Complex{Int} data at the entry, so FillExtrap
payloads normalize into the float value space. The unit-axis arm of
`_value_type` fell to the duck passthrough and kept promotable data raw —
`FillExtrap(0.5)` then normalized into Int and threw InexactError (every ND
family; linear included). A layered dispatch arm now floats promotable values
for non-Real tags while duck REAL grids (Dual) keep the raw-value contract.
Float32/ComplexF64 are fixpoints of `float`, so only non-float promotable data
changes — and now converts at build exactly like the Real-axis entry.
Pins: Int data + FillExtrap(0.5) across linear/cubic persistent + cubic
one-shot on unit axes, plus in-domain === parity with the Real-grid build.
…tity coordinate
`_resolve_grididx` rebuilds the wrapper at the axis eltype, but the payload was
pinned `T <: Real`, so every unit-axis GridIdx query (persistent and one-shot,
all families) died in the constructor with a TypeError. The struct and inner
ctor now admit `T <: Number`.
The promotion contract deliberately stays `S <: Real`: Unitful's generic
Quantity-vs-Real rules see the wrapper (a Real) and nest the numeric type
(`Quantity{Quantity{…}}` + convert ambiguity) when promoted against Quantity.
Unit coordinates instead flow through the one value-consuming seam — a direct
`-(g::GridIdx, x::Number) = g.val - x` method covering every dL/α read; Real
payloads fold identically.
Pins: persistent all-GridIdx / mixed tuples (cubic + linear, === parity),
one-shot PreCompute + mixed-tuple Auto mirrors, and the all-GridIdx+Auto
hetero-collapse refusal contract.
…ore serves The cubic/quadratic GriddedQuery fast path stores PHYSICAL h/inv_h/dL anchors and calls the cell kernels directly, bypassing the scaled-store restore seam — on unit axes the persistent form threw DimensionError and the unified one-shot (`interp(grids, data, gq; method = CubicInterp())`) died on a Quantity cache key. Dispatch now declines for non-Real tags: - persistent `_gridded_eval_itp_methods!` cubic/quadratic arms are constrained to `Tg <: Real`; unconstrained sibling arms return `false`, so the shared restore-aware `_nd_batch_pointwise!` core serves unit axes. - the one-shot hooks route through a tag-dispatched `_*_gridded_oneshot_try!` pair — non-Real declines to the query-protocol batch path. Real axes keep the fused path untouched (specificity, not a boolean test). Pins: persistent value === pointwise, deriv ≈ (two FP orderings), quadratic parity, and the unified one-shot mirror.
…load space The quadratic scalar and in-place entries passed raw `Tv` to the eager `_resolve_extrap`, so Int data + `FillExtrap(0.5)` threw InexactError at resolve time (even for in-domain queries) — on Real axes (pre-existing) and on the newly opened unit axes alike. `_oneshot_fill_eltype` moves to the shared ND home, parameterized by the family coeff witness (`_coeff_op2` cubic / `_coeff_op` quadratic); the in-place entry reuses `_nd_promote_grids`' promoted `Tv_p` exactly as the cubic mirror.
… refusal Hermite ND has no scaled-store/reparam seam (user partials live per-axis in [Y/Xᵈ]); unit axes died in deep MethodErrors — persistent ctor in `_pack_and_extend_nodal_derivs`, one-shot in abstract-Tv coerces — on master and this branch alike (A/B probed, not a branch regression). Both entries now run `_check_nd_hetero_grid` (the F11/F17 refusal contract), which also keeps the reparam params fallback structurally unreachable from non-solver families. Real axes and the whole Hermite ND Real suite are untouched.
…×deriv on unit axes The reparam bounded arm's `sign == 0` short-circuit (`z_t * vol`) and the GridIdx `-` seam feeding the reparam dLs had no unit-grid coverage — both probe correct today; pin the zero's [Y·X₁·X₂] type, the one-flip sign, and bit-identity of GridIdx derivs against flat coordinates (cubic + quadratic).
`_check_nd_reparam_eltype` probed `promote_op(_coeff_op, Tg, Tg)` — an op that never calls `oneunit` (and misuses `Tg` in the value slot) — so a `*`/`inv`-complete duck Number without `oneunit` passed the gate and died deep in float-promotion MethodErrors. The probe now runs `_reparam_op` (`x * _deriv_oneunit(x, DerivOp(1))`), the canonical per-element twin transform `_reparam_grids` applies, so the gate checks exactly what its message promises. Unitful/Dual populations unchanged.
The adapter also backs the scaled-store in-domain restore seam (_restore_nd_deriv_scale), not just FillExtrap OOB zeros — say so where a "fill" grep lands.
… stop allocating The scaled-store families reparameterized by materializing a twin axis per call (`x .* inv(oneunit(x))`), so unit Vector axes paid an O(n) array on every one-shot query and every `integrate` — and the twin, being a plain Vector, also lost the parent's cached h/inv_h, making bounded integrate 2.2× slower than its Real sibling. Both consumers (the fiber solve and the separable integrate engine) only ask an axis for `getindex` and `_get_h`/`_get_inv_h`, so the twin needs no array: `_ReparamAxis` is a view that applies `_reparam_op` per access and takes its widths from the parent's cache times a scalar. Range axes keep the arithmetic form (Base's range broadcast is already isbits and preserves the `_CachedRange` objectid fast path in the solve cache banks). While pinning the contract, bounded integrate still allocated on mixed-unit axes: `_normalize_bounds_nd` indexed the bounds tuple with a runtime `d`, leaking the element Union — the same class its neighbours already unroll. Measured (n=400 Vector axes): one-shot 3328 B → 0 B, integrate full 7232 B → 0 B, bounded 7424 B → 0 B and 245 µs → 114 µs (Real: 114 µs). Eval values are bit-identical; bounded integrals move ≤1 ULP because the cached physical width is scaled once instead of differencing two scaled coordinates.
`Base.:-(::GridIdx, ::Number)` was ambiguous against Base's `-(::Real, ::Complex)` arms (GridIdx <: Real) — Aqua caught 2, and disambiguating would need a method per Complex arm. The seams read the payload directly instead: `_coord_value` is identity for every non-GridIdx query, so Real paths keep their promotion behaviour (A/B: linear 10k batch 12.2 µs both sides, 0 alloc, values `===`). Applied at each family's canonical coordinate seam — `_alpha_of` (linear/constant 1D+ND, plus the exclusive-periodic wrapper's own arm), the reparam local params (cubic/quadratic ND), and the Constant ND kernel's `dL`. The Constant ND seam had no unit-axis GridIdx coverage; the pin now sweeps every family × (all-GridIdx, mixed, derivative).
…ding `oneunit(eltype)` twice plus a `first(inner)` element read (dead — only its type mattered) become a single `u`. Same codegen (3 pointer copies + two constant stores), one less way to trip on an empty axis.
A bare `GridIdx(k)` carries `val = NaN` until `_resolve_grididx` meets it with its axis. Scalar entries resolve at the door; the ND batch loops did not — so a GridIdx on an INTERPOLATING axis evaluated the kernel at NaN: linear/cubic returned NaN, constant silently returned a wrong cell, and unit axes threw a DimensionError. Existing batch coverage only put GridIdx on NoInterp axes, whose kernel never reads the coordinate. The query protocol grows a grid-aware `_extract_query_point`; the eight ND eval batch loops pass the axes they already evaluate against. Identity for non-GridIdx elements — A/B on 20k-point batches: linear 376.8→377.4 µs, constant 372.8→374.9 µs, cubic 582.5→564.7 µs, 0 B throughout.
…data eltype `constant_interp`/`cubic_interp` (1D), `ConstantInterpolantND` and the gridded constant entry promoted a `FillExtrap` payload against `eltype(y)`, so Int data + `FillExtrap(0.5)` threw `InexactError` at construction — even though those same interpolants already return Float64 for Int data (the kernels float against the grid). They now use linear's spelling, `_value_type(Tv, <float grid type>)`. Unit axes inherit it through the same call. Pins sweep 1D persistent (5 families) + one-shot + ND + gridded.
Both verified identical on the v0.4.17 tag, so neither is a regression from this branch: - 1-D batch `GridIdx` — the 1-D vector loops never resolve the wrapper, so the domain check sees `val = NaN` and throws. Loud, never silent. - SoA with a pinned scalar axis — a mixed `(vector, scalar)` tuple has no protocol arm, so `_query_length` reports the tuple arity. Loud except when the vector length equals N, where the size check passes and the kernel reads past a scalar under @inbounds. Each pin states the intended contract, so whoever fixes it gets an "Unexpectedly Pass" instead of having to rediscover the case.
Contributor
FastInterpolations.jl BenchmarksAll benchmarks (56 total, click to expand)
|
| Benchmark | Current | Previous | Imm. Ratio | Grad. Ratio | Tier |
|---|---|---|---|---|---|
5_linear_construct/g0100 |
37.77 ns |
32.5 ns |
1.164 |
1.089 |
immediate |
Thresholds: immediate > 1.1x (vs latest master), gradual > 1.1x (vs sliding window)
Runner:
znver3|4c— AMD EPYC 7763 64-Core Processor, julia 1.12.6. Times are min-merged and compared only against this same machine's history.
This comment was automatically generated by Benchmark workflow.
mgyoo86
changed the title
…s branch moved Three CI failures, all from tests that pinned behaviour this branch changed: - **Constant's fill space** was mine to revert, not the test's to update. `ConstantInterp` returns data values verbatim, so its fill belongs in `Tv` and `FillExtrap(NaN)` beside Int data is rejected at construction — Int has no NaN. Only Cubic keeps the promotion fix (it solves, so its value space is float). The sweep now says which families accept a float fill and why, and pins NaN propagation wherever the value space holds NaN (Float64/Float32). - **Duck-grid lint** — the scaled-store policy seams are dispatch PAIRS whose unconstrained sibling serves duck/unit grids, the same shape as the existing `_inv_const` allowlist class. Enumerated with a reason, ratchet updated. - **Duck BC boundary** — the bare-PointBC arm now follows #201's rule that a structural Real zero mints into the payload space, so `Deriv1(0.0)` builds on duck `y` while a nonzero payload still fails.
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## master #202 +/- ##
========================================
Coverage 96.31% 96.31%
========================================
Files 164 164
Lines 13852 13965 +113
========================================
+ Hits 13341 13450 +109
- Misses 511 515 +4
🚀 New features to boost your workflow:
|
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Extends #201's unit-native solver core to the ND solver families (Cubic/Quadratic PreCompute): build, eval, derivatives,
integrate(full + bounded), one-shot, gridded queries,GridIdxandFillExtrapall work on unit-carrying axes, mixed units included. The blocker was the shared2^Nnodal-derivative store, whose slots differ by derivative order ([Y],[Y/X₁],[Y/X₁X₂]) and so cannot live in one homogeneous array. Each axis is therefore reparameterized to its exact dimensionless twint = x·inv(oneunit(x)), the store stays[Y]-homogeneous, andgrid⁻ᵏis restored at a single cell seam.Mechanism
_reparam_op(x) = x * _deriv_oneunit(x, DerivOp(1))— the gate probes it, the twin applies it, its eltype comes from it. Eval never builds a twin: it scalesh/inv_h/dLper query and restores the grading with one scalar.getindexand_get_h/_get_inv_h, so_ReparamAxisis a view over the parent plus two scalar witnesses — no per-call array, widths straight from the parent's cache. Range axes keep the arithmetic form (already isbits, and it preserves the_CachedRangecache-bank fast path).Rides along
GridIdxin ND batches carried itsval = NaNpoison into the kernel: on an interpolating axis, linear/cubic returnedNaNand Constant a wrong cell — silently, on Real grids, since v0.4.17. Existing coverage only pinnedGridIdxonNoInterpaxes, whose kernel never reads the coordinate.FillExtrappromotion space — 1-Dconstant_interp/cubic_interp,ConstantInterpolantNDand gridded constant promoted the fill against the raw data eltype, so Int data plusFillExtrap(0.5)threwInexactErrorat construction.oneunit-less duck used to pass), Hermite ND gates non-Real axes at both entries, the gridded fused path declines them, and_normalize_bounds_ndunrolls its bounds tuple.Numbers
Real is unchanged by A/B: linear 10 k batch 12.2 µs both sides with values
===, ND 20 k batch within noise, 0 B throughout. Unit axes now allocate nothing on the hot paths (one-shot 3328 B → 0,integratefull 7232 B → 0, bounded 7424 B → 0) and bounded integration went 245 µs → 114 µs, matching Real, because the lazy twin keeps the parent's cached widths. Eval values are bit-identical; bounded integrals move ≤ 1 ULP.Deferred
Adjoints/AD, Hermite ND and the hetero families on unit axes all refuse with actionable errors. Two query shapes stay broken and are now recorded with
@test_broken, both verified identical on the v0.4.17 tag: 1-D batchGridIdx(loud — the 1-D loops never resolve the wrapper) and SoA queries with a pinned scalar axis (no protocol arm for a mixed(vector, scalar)tuple; fixing it is a design choice since theNoInterppre-slice entry uses the same shape).