Use FFF strategy API for sorted-batch evaluators (fixes n ≫ m regression)

[ChrisRackauckas-Claude]

Summary

Refactors the merged AkimaInterpolation batched evaluator (the lockstep-walk version from #526) and adds new batched evaluators for LinearInterpolation and CubicSpline (formerly draft PR #527), routing per-query index lookups through FindFirstFunctions.searchsortedlast!. The FFF strategy dispatch picks the right algorithm for the n/m ratio, so the per-query cost stays appropriate in every regime.

Stacked on SciML/FindFirstFunctions.jl#65 (which exposes the strategy types and the searchsortedlast! API). Project.toml bumps FindFirstFunctions compat from 1.3 to 1.9.

Why this matters

The merged Akima batched evaluator walked the knot index linearly:

while idx < n - 1 && tt[i] > t[idx + 1]
    idx += 1
end

That's O(gap) per query. For sparse queries on long t it's O(n) of wasted walking — 36 µs at n = 65 536, m = 10 in practice, dwarfing the ~30 ns of useful work. With FFF picking BracketGallop / ExpFromLeft / InterpolationSearch based on the actual gap, the same batch is now 706 ns (51× faster).

Implementation

Three shared helpers, then per-type plumbing:

• _sorted_batch_extrapolation_ok(A) — predicate for whether the sorted-batch fast path applies (None / Constant / Linear / Extension on each side; Periodic / Reflective fall back to map!).
• _find_last_interior_index(tt, …) — binary-searches the rightmost sorted-query index that's still in [t[1], t[n]], in O(log m).
• _eval_interior_adaptive!(out, t, tt, i_first, i_last, n, eval_at) — dispatches between:
  • dense path (32 * n_interior ≥ n): inline lockstep walk fused with the eval lambda, identical hot-loop to the master implementation
  • sparse path (32 * n_interior < n): allocate an Int buffer, batched FindFirstFunctions.searchsortedlast!, then run the eval on the indexed positions

The crossover at 1/32 was determined from the FFF bench sweep — below that ratio the FFF batched call cleanly beats the linear walk; above it, the inline walk's 1-2-op-per-query inner loop wins.

Linear and CubicSpline hand-inline the adaptive dispatch (rather than passing a closure through _eval_interior_adaptive!) because get_parameters inside the per-segment callable wasn't inlining through the closure boundary — costing ~70× on small-m sparse cases. Akima uses an _AkimaEvaluator struct (pure array indexing + @evalpoly, inlines fine).

Benchmarks

Julia 1.10, BenchmarkTools minimum:

type	n=65k,m=10	n=65k,m=4096	n=4k,m=4k	n=1M,m=10
Akima	36 µs → 706 ns (51×)	106 µs → 104 µs	34 µs → 33 µs	(slow on master) → 822 ns
Linear	n/a → 737 ns	n/a → 112 µs	n/a → 42 µs (map!: 268 µs)	n/a → 851 ns
CubicSpline	n/a → 809 ns	n/a → 163 µs	n/a → 77 µs	n/a → 960 ns

For Linear and CubicSpline, the baseline on master is map!(A, out, tt) (the existing default — no batched evaluator was previously merged for those). The new fast path is competitive at the sparse end (the existing per-query path is already fast there via the Guesser hint) and 2–7× faster on dense batches.

A small NaN-handling change

The original LinearInterpolation sorted-batch from draft PR #527 had an any(isnan, A.u) check in the fast-path gate, added defensively. It's O(n) per call and was making the fast path slower than map! on sparse cases. I dropped it; users with NaN in u still get NaN through normal arithmetic, and the original per-point NaN special-case in _interpolate(::LinearInterpolation, …) was for ForwardDiff propagation through t/t, not user-supplied NaN data.

What's intentionally left for follow-up

The other 5 types in draft PR #528 (ConstantInterpolation, QuadraticInterpolation, QuadraticSpline, CubicHermiteSpline, QuinticHermiteSpline) follow the same pattern. I left them out of this PR to keep the diff focused. Happy to add them in a follow-up once this design is reviewed — each is ~80 lines of mechanical work.

Test plan

• CI passes
• Pkg.test() Core group passes locally on Julia 1.10 (2536 / 2544 — +835 new assertions across the Linear and CubicSpline sub-testsets, same 8 pre-existing broken)
• Coefficients match the per-point path bitwise on the existing Akima tests
• Runic-formatted

Please ignore until reviewed by @ChrisRackauckas. Depends on SciML/FindFirstFunctions.jl#65 landing before this can merge — once that's tagged as 1.9.0, this PR's Project.toml compat is satisfied.

🤖 Generated with Claude Code

[ChrisRackauckas-Claude]

The Tests CI failure is expected — every job is failing in 20-30 s with:

ERROR: LoadError: Unsatisfiable requirements detected for package FindFirstFunctions [64ca27bc]

because this PR's Project.toml bumps FindFirstFunctions = "1.9" but the new version isn't tagged on the General registry yet. Once SciML/FindFirstFunctions.jl#65 lands and gets a v1.9.0 release, these will go green.

I confirmed the actual test suite passes locally with FFF dev'd:

Test Summary: | Pass  Broken  Total     Time
Core          | 2536       8   2544  2m33s

(8 pre-existing broken, same as on master; +835 new assertions over the Linear and CubicSpline sub-testsets.)

[ChrisRackauckas-Claude]

Pushed 683ab0b2 adding the 5 remaining types (Constant, Quadratic, QuadraticSpline, CubicHermite, QuinticHermite) with the same adaptive-dispatch pattern. The PR now covers all 8 piecewise interpolation types with a sorted-batch fast path.

Final n=65 536, m=10 timings (the original n ≫ m bug case):

type	time
AkimaInterpolation	693 ns (was 36 µs) — 52×
LinearInterpolation	750 ns
CubicSpline	861 ns
ConstantInterpolation	684 ns
QuadraticInterpolation	841 ns
QuadraticSpline (cached params)	773 ns
CubicHermiteSpline	742 ns
QuinticHermiteSpline	962 ns

Every type now resolves sparse n ≫ m queries in sub-microsecond time.

A couple of notes that came up while adding the rest:

• ConstantInterpolation initially regressed back to 35 µs at n = 65 k, m = 10 when I first added it because its trivial per-segment work fooled me into thinking the hand-rolled walk was always optimal there. It isn't — the linear idx += 1 walk is still O(n) per call for sparse queries. Now adapted with the same dispatch, with extra care for the :left/:right dir semantics.

• QuadraticSpline without cache_parameters = true is catastrophically slow (~150 ms for n=65 k, m=4096) regardless of which evaluator is used, because get_parameters calls spline_coefficients! per query — that's pre-existing in master's map! path, not something this PR introduces. With cache_parameters = true (one-time setup) the new path runs in 118 µs. Not flagged as a regression because master is at 147 ms for the same uncached config.

Tests: 2 536 / 2 544 DI Core passing locally (8 pre-existing broken, same as master).