Reject non-one-based out/tt in sorted-batch itp(out, ts) kernels

[ChrisRackauckas-Claude]

What changed and why

The sorted-batch kernels behind itp(out, ts) (LinearInterpolation,
QuadraticInterpolation, AkimaInterpolation, ConstantInterpolation,
QuadraticSpline, CubicSpline, CubicHermiteSpline, QuinticHermiteSpline)
only check length(out) == length(tt) before indexing both arrays as if
they were 1-based. An OffsetArray out or tt is therefore read/written
out of bounds, returning silently wrong values, and a prior probe
segfaulted the GC when all 10 interpolation types were exercised this way
in one process.

Update: the first version of this PR added
Base.require_one_based_indexing(out, tt) unconditionally right after the
length check. Independent review (Cursor, round 1) correctly flagged that
as a regression: the unsorted path already falls through to
map!(A, out, tt), which handles arbitrary axes correctly on master
(confirmed empirically below), so unconditionally throwing broke calls
that previously worked.

Fixed by guarding only the unsafe @inbounds sorted-batch kernel:

if !Base.has_offset_axes(out, tt) && _sorted_batch_extrapolation_ok(A) && issorted(tt)
    _linear_eval_sorted!(out, A, tt)
else
    map!(A, out, tt)
end

Non-one-based out/tt (sorted or not) now always routes to map!,
which is correct for any axes. One-based, sorted input still takes the
fast @inbounds path. No input that worked on master now errors; the
only behavior change is that the two inputs that were previously
silently wrong (sorted tt/out with offset axes) are now correct.

Failing before / passing after

Before (unmodified master, commit b98db80), reproduced in a clean
temp environment with Pkg.develop pointed at the unmodified checkout:

=== candidate 1: itp(out, ts) with OffsetArray ts, 1-based out ===
correct = [0.8753842058167891, 0.5252087174427744]
master out = [0.5252087174427744, 2.4456522643573106]  (wrong if != correct)

LinearInterpolation(sin.(1:10), 1:10)(zeros(2), OffsetArray([1.5, 2.5], 0:1))
silently returns [0.5252, 2.4457] on master where [0.8754, 0.5252] is
correct — i.e. it's reading/writing through the wrong 0-based offsets.
That input is sorted, so it hits the broken fast path.

Confirmed separately that master's unsorted path already gives the
correct answer via map!, which is exactly what this PR's fix now also
guarantees for the sorted case:

master-equivalent unsorted OffsetArray tt, 1-based out:
  result = [-0.30784124362403054, 0.8753842058167891, 0.5252087174427744]  match=true
unsorted 1-based tt, OffsetArray out:
  result = [-0.30784124362403054, 0.8753842058167891, 0.5252087174427744]  match=true

After this fix, the sorted case is also correct:

julia> A(zeros(2), OffsetArray([1.5, 2.5], 0:1))
[0.8753842058167891, 0.5252087174427744]   # correct, via map! fallback

Updated test in test/interpolation_tests.jl ("OffsetArray batch safety")
exercises all 8 covered types with offset out, offset tt, and both
offset together, for both sorted and unsorted query order, and asserts
the result matches the scalar call in every case — including the cases
that already worked on master (unsorted, offset) and the case that was
silently wrong on master (sorted, offset). This test fails against
unmodified master (wrong values on the sorted/offset case) and passes
with the fix; it also fails against the first version of this PR (throws
on the unsorted/offset case) and passes with the current fix.

What was not verified

• Did not re-run the GC segfault repro from the original probe (all 10
  types evaluated with offset arrays in one process); the out-of-bounds
  write this PR closes is the same code path that produced it.
• GPU/CUDA paths, downstream packages, and the pre Julia channel were not
  exercised locally.
• Core and QA groups were run on Julia 1.12 and 1.10 (lts) after this
  fix; pre is listed in test_groups.toml but was not run.

Anything a reviewer should push back on

Base.has_offset_axes is the right public-ish Base check for "does this
need the map! fallback", matching the semantics map!/copyto! already
use internally, but worth double-checking there isn't a cheaper
AbstractUnitRange-based check the package would prefer for consistency
with _cumulative_integral's use of Base.require_one_based_indexing
elsewhere (that call site has no working fallback, so throwing there is
still correct and unrelated to this change).

Please ignore until reviewed by @ChrisRackauckas.

Risk assessment

• Risk: low
• Blast radius: only the 8 sorted-batch itp(out, ts) methods in
  interpolation_methods.jl. One-based out/tt (the overwhelming common
  case) is unaffected; non-one-based input now always gives a correct
  result via the existing map! fallback instead of either corrupting
  memory (master) or erroring (this PR's first version). No public API
  signature changes.
• Evidence: failing-before/passing-after shown above, including the
  regression the first version introduced and this version fixes.
• Independent review: Cursor reviewed round 1 and returned CHANGES
  (regression on the unsorted offset path); addressed in this revision.
  A fresh independent pass is pending on the current head.
• Merge: needs human review — safety fix touching 8 hot-path methods;
  wants an independent pass before merge despite being low-risk.

---

🤖 Generated with Claude Code

https://claude.ai/code/session_01LPHREnnonfLg1VcE1EJovv

Independent review: Devin Fusion (fusion-claude-opus-5-5-high-sidekick-swe-2-medium) rated it low, verdict MERGE: #622 (comment)

[ChrisRackauckas]

🤖 Automated comment from an AI agent running as @ChrisRackauckas — not written or reviewed by Chris.

Independent review (Devin CLI 3000.11.3, model fusion-claude-opus-5-5-high-sidekick-swe-2-medium): MERGE, risk low.

Full review

VERDICT: MERGE
RISK: low

Review of #622 (head 0aa37ba, base origin/master = b98db80). The fix is correct. Every change to the source is the same one-term guard, !Base.has_offset_axes(out, tt) && …, added in 8 places. When out and tt are one-based this guard is a no-op, so the fast path behaves exactly as it did before. When either has offset axes, the call falls back to map!, which handles any axes. I found no wrong results, no loosened tests and no swallowed errors. Everything below is cosmetic or a follow-up.

Blocking findings

None.

Non-blocking findings

1. The PR title says the opposite of what the code does. (Confirmed by reading the diff and the PR.) The title is "Reject non-one-based out/tt in sorted-batch itp(out, ts) kernels", which describes the first commit f78f25e (require_one_based_indexing, i.e. throw). The head commit 0aa37ba routes these inputs to map! and throws nothing. A squash merge would put the wrong description into history. Retitle it, e.g. "Route offset-axis out/tt to map! in sorted-batch itp(out, ts)".

2. The test comment includes history and a reference that leads nowhere. (Confirmed by reading.) See test/interpolation_tests.jl:2163-2168:

   • "(which already handles arbitrary axes correctly, as it does on master for unsorted tt)" is retrospective wording. CLAUDE.md says to describe what exists, not what came before.
   • "(https://github.com/SciML/DataInterpolations.jl, candidate 1)" is the repo root plus the name of an internal probe. It is not an issue link and means nothing to a reader.
   • There are about 9 comment lines out of about 60 added (~15%), above the ~10% guideline.

   Suggested fix: cut it to one line saying non-one-based out/tt must give the same result as the scalar call.

3. The comment on the last test block is not checked. (Confirmed by reading.) test/interpolation_tests.jl:2200-2201 says the one-based sorted input "takes the sorted fast path", but nothing asserts that. The assertion is only out ≈ A.(q_sorted). Either drop the parenthetical or add a check. I checked the fast path myself with timings (see "What I ran"), so this is about the comment's wording, not a bug.

4. Offset inputs hit a performance cliff. (Confirmed by running.) Sorted 1e4-point CubicSpline batch:

   • one-based input, PR: 52.5 µs (master: 52.4 µs, so no change)
   • offset-tt input, PR: 410 µs, the same as the unsorted map! path
   • offset-tt input, master: 52.6 µs, but the output there is garbage

   This trade is acceptable: correctness first, and offset inputs are rare. Writing the _*_eval_sorted! kernels in terms of firstindex/eachindex would keep the fast path for them. That belongs in a separate follow-up, not this PR.

5. PR body attribution does not follow CLAUDE.md. (Confirmed by reading.) The body lacks the required first-line 🚧 UNREVIEWED banner and ends with "Generated with Claude Code" rather than the footer format (harness version and model). The commit trailers say Agent-Harness: Claude Code with no version. This has no effect on the code.

6. Out of scope, and already present on master: an interpolation built from offset-axis u/t fails with BoundsError instead of a clean error. (Confirmed by running.) LinearInterpolation and ConstantInterpolation reject offset u/t with a clean ArgumentError. CubicSpline(OffsetArray(u,-1), OffsetArray(t,-1)) throws a BoundsError during construction instead. The behaviour is the same on master and this PR does not touch it. Noted only because it is the next neighbour of the same bug class.

What I ran

All runs used TMPDIR/scratch under the review directory. Master was a pristine copy from git archive origin/master, developed into its own env.

• The PR's own example (LinearInterpolation(sin.(1:10), 1:10)(zeros(2), OffsetArray([1.5,2.5],0:1)), Julia 1.12, default bounds checking) reproduced exactly:

  • master: [0.5252…, 2.4457…] (wrong)
  • PR: [0.8754…, 0.5252…], which matches the scalar result

• My adversarial probe scratch/probe.jl covered:

  • 11 types: the 8 the PR names, plus PCHIPInterpolation, LagrangeInterpolation and BSplineInterpolation
  • extrapolation set to Extension
  • 5 query sets: sorted small, sorted with extrapolation on both sides, sorted dense (400 points, which reaches the large-n interior kernel), unsorted, and a single point
  • offsets −1, 1, 7 and −1000
  • offset tt, offset out, both offset, and both offset with different offsets, plus one-based controls

  1100 cases in total. Every result was compared against the scalar A.(q).

  • PR, Julia 1.12, --check-bounds=yes: 0 wrong, 0 errors
  • PR, Julia 1.12, default bounds: 0 wrong, 0 errors
  • PR, Julia 1.10.12 (lts), --check-bounds=yes: 0 wrong, 0 errors. So map!'s positional iteration with mismatched axes also holds on LTS.
  • master, Julia 1.12, --check-bounds=yes: 576 BoundsErrors, all in the sorted cases with offset axes. That is 9 types (PCHIP is a CubicHermiteSpline) × 4 sorted query sets × 4 offsets × 4 offset variants. Lagrange and BSpline are unaffected, as expected, because they use the generic map! method.

• The new testset alone (extracted to scratch/newtest.jl):

  • PR, Julia 1.12 with --check-bounds=yes, and PR, Julia 1.10: 56/56 pass
  • master, --check-bounds=yes: 8 errors (BoundsError)
  • master, default bounds: 8 errors (out-of-bounds reads of tt returned garbage below t[1], which raised LeftExtrapolationError)

  So the test fails before the fix and passes after it. Its reference values come from the scalar call, not from the batch output being tested.

• Full suite: GROUP=Core julia +1.12 --project -e 'using Pkg; Pkg.test()' passed. interpolation_tests.jl had 17627 pass and 5 broken; the broken count is unchanged because the diff adds no @test_broken. extrapolation_tests, interface, abstract_interface and parameter_tests all passed.

• Formatting and spelling:

  • julia +1.12 --project=@runic -m Runic --check src/interpolation_methods.jl test/interpolation_tests.jl: exit 0
  • git diff origin/master...HEAD | typos -: clean

• Public API check: Base.ispublic(Base, :has_offset_axes) is true on 1.12.7 and 1.13. It is false on 1.11, and ispublic does not exist on 1.10. QA runs on lts and 1 (= 1.12), so ExplicitImports' public-access check should pass. Base.require_one_based_indexing, which the package already uses, has the same status.

• Matrix-valued u: a sanity check of the generic (out::AbstractArray, t) method gave correct results. That method is unchanged.

• Benchmarks: scratch/probe2.jl, using BenchmarkTools; numbers are in finding 4.

What I did not verify

• GROUP=QA (Aqua/ExplicitImports/JET) was not run, because it loads Makie, Mooncake and Symbolics and the host was memory-limited. I judge it low risk because has_offset_axes is public on 1.12, but I did not run it.
• The PR's GC-segfault claim was not reproduced. On master with default bounds I saw only silently wrong reads.
• Not exercised: Julia pre, the 32-bit Core job, GPU arrays, and downstream packages.
• Linked or tracking issue: the PR references none (closingIssuesReferences is empty), so there was no sibling work to check against.
• New test dependency: I did not check the OffsetArrays license myself. It is widely known to be MIT, and it is only a test dependency with [compat] OffsetArrays = "1".

---

🤖 Posted by an AI agent — harness: Devin CLI 3000.11.3 (review), Claude Code 2.1.285 (posting) · model: fusion-claude-opus-5-5-high-sidekick-swe-2-medium
Conversation: local Claude session 3c311569-dda6-4687-bc9f-65febb212c33 (fleet master)