Skip to content

Fix sampled evalrule silent wrong results on OffsetArrays - #379

Merged
ChrisRackauckas merged 1 commit into
SciML:masterfrom
ChrisRackauckas-Claude:int-sampled-offset-axes
Oct 16, 2026
Merged

ChrisRackauckas merged 1 commit into
SciML:masterfrom
ChrisRackauckas-Claude:int-sampled-offset-axes

Conversation

@ChrisRackauckas-Claude

@ChrisRackauckas-Claude ChrisRackauckas-Claude commented Oct 9, 2026 •

Copy link
Copy Markdown
Member

Summary

The @inbounds fast paths in evalrule for AbstractVector/AbstractMatrix loop over 1:m / 1:n, which is incorrect when data has offset axes. With an OffsetArray whose non-integrated axis is offset (e.g. OffsetArray(ones(2,5), 0:1, 1:5) with dim=2), this silently returns wrong values such as [1.0, 0.875] instead of [1, 1] under default bounds checking (and BoundsError with --check-bounds=yes). This PR adds && !Base.has_offset_axes(data) to both fast-path gates so OffsetArrays take _evalrule_general, while keeping @inbounds on the dense 1-based path.

Evidence

Regression (default bounds checking):

  • Before (src fix stashed): collect(sol.u) == [1.0, 1.0] failed with [1.0, 0.875]; OffsetVector path also BoundsError at index 6.
  • After: OffsetArrays testset 2/2 pass; full sampled_tests.jl green on Julia 1.12 and julia +1.10.

Matrix fast-path timing (JULIA_NUM_THREADS=1 taskset -c 47, @benchmark minima, 4×10001 trapz dim=2):

  • Before: 44189 ns (5 allocs)
  • After (warm): 44210 ns / 44209 ns (5 allocs) — unchanged within noise.

Test groups

GROUP=Core: all passed, including Core/sampled_tests.jl | 32 32 10.4s. Final: Testing Integrals tests passed.

GROUP=QA: QA/qa.jl | 27 27 2m34.2s. Final: Testing Integrals tests passed.

Runic (inplace) and typos clean on changed files.

Not verified

  • Full GPU / non-CPU array paths beyond the existing fast_scalar_indexing gate.
  • Whether _evalrule_general preserving OffsetArray axes on the result (vs always collecting) is the preferred API; the regression uses collect(sol.u).
  • CI matrix on all declared Julia versions beyond local 1.12 Core/QA and local 1.10 sampled regression.

Reviewer pushback

  • Prefer fixing the fast path to use axes/eachindex instead of bailing to the general path (may keep more performance for OffsetArrays).
  • Whether OffsetArrays should be a documented supported input for SampledIntegralProblem (SciMLBase already checks integrated-axis index agreement).
  • Whether a patch version bump is warranted for a silent correctness fix (no new public names; left at 5.5.2).

Please ignore until reviewed by @ChrisRackauckas.

Risk assessment

Independent review: pending

🤖 Generated with Cursor Agent 2026.10.01-14929f9 (model: unknown (Cursor auto)), transcript /home/crackauc/sandbox/goals/performance/jobs/integrals/int-offset/log.txt on amdci2; orchestrated by Claude Code (claude-opus-5-5[1m]) https://claude.ai/code/session_01LPHREnnonfLg1VcE1EJovv

Made with Cursor

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

Risk assessment

  • Risk: low
  • Blast radius: Offset-indexed sampled vectors and matrices now use the existing generic integration path. Ordinary one-based arrays retain their fast path. Offset matrix results can preserve the remaining offset axis. OffsetArrays is an MIT-licensed test-only dependency.
  • Evidence: All 15 reported PR checks pass. Locally, the head sampled test file passes 32/32 assertions. Running the same new test file on the PR base with --check-bounds=yes raises BoundsError at matrix index [2, 1] in the offset testset, confirming that it catches the original unsafe indexing. Five additional kernel assertions with nonconstant offset data pass for vectors and matrices integrated along either dimension, comparing to equivalent one-based data.
  • Independent review: codex-cli gpt-6.1-sol rated it low, high confidence: the guard removes unsafe one-based indexing, reuses the established iterator path, and passes package and nonconstant-value probes.
  • Merge: auto-merge candidate — narrow internal correction with targeted regression tests and all reported checks passing; no blocking findings.

Independent risk assessment by codex-cli (gpt-6.1-sol), reviewed head 5bc85af7e1; dispatched by the fleet master, Claude Code (claude-opus-5-5), https://claude.ai/code/session_01VZjwf7Ux9YtTzwHGHMqnsq

…xes.

Fast paths looped 1:m/1:n under @inbounds; gate them out when data has offset axes so OffsetArrays take the general path.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Cursor Agent <noreply@cursor.com>
Agent-Harness: Cursor Agent 2026.10.01-14929f9
Agent-Model: unknown (Cursor auto)
Agent-Session: local session, transcript /home/crackauc/sandbox/goals/performance/jobs/integrals/int-offset/log.txt on amdci2
@ChrisRackauckas

Copy link
Copy Markdown
Member

🤖 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

Two-line source change: && !Base.has_offset_axes(data) on the AbstractVector and AbstractMatrix fast-path gates in src/sampled.jl:53,70. Arrays with offset axes now go through _evalrule_general. That path already iterates eachslice zipped with the weights' own iterate, so it respects axes. OffsetArrays (MIT) is added only as a test dependency, in [extras], [compat] and [targets]. The PR links no issue.

Blocking findings

None.

Non-blocking findings

  1. test/sampled_tests.jl:45-55: the regression test uses constant data (ones). Confirmed by running: it fails on master ([1.0, 0.875], plus a BoundsError at index 6) and passes on the PR. But constant data cannot catch a weight/data misalignment along the integrated axis, such as a shifted or reversed weight index. My probes (below) cover that with random data. Suggested: use rand data and compare against the same problem built from 1-based arrays (≈, or == where the summation order is identical). I would not block on this.
  2. The return type changes for offset input (confirmed by running bench.jl). On master, a matrix with an offset non-integrated axis returned a 1-based Vector, though its values were wrong (memory read out of bounds under @inbounds). On the PR it returns an OffsetVector with axes 0:3, which matches what the existing N-d general path already did. I think that is the right behaviour, but the test's collect(sol.u) hides it. Asserting axes(sol.u) == axes(y, 1) would pin the behaviour.
  3. Sibling PRs open on the same file: SIMD the 1-D sampled-data weight×data reduction #381 (adds @simd inside the vector fast path) and Stabilize sampled eachslice integration with a function barrier #380 (adds a function barrier in _evalrule_general). Suspected from reading their diffs: neither touches the gate lines this PR edits, and the semantics compose. Stabilize sampled eachslice integration with a function barrier #380 speeds up the path that OffsetArrays now take. I did not test a merge.
  4. No version bump. The PR body leaves that to the maintainer. A patch bump is reasonable for a silent-wrong-result fix.

What I ran

All runs used Julia 1.12.7 unless noted, on macOS aarch64 (M2 Max). Environments were in scratch/: env develops the PR checkout, envm develops a copy of the repo with src/sampled.jl reverted to origin/master.

  • Fails before the fix (confirmed): julia --project=envm repo/test/sampled_tests.jl.

    • Sampled integration with OffsetArrays: Test Failed ... Evaluated: [1.0, 0.875] == [1.0, 1.0].
    • BoundsError: attempt to access 5-element OffsetArray(::Vector{Float64}, 0:4) ... at index [6].
  • Passes after the fix (confirmed): julia --project=env repo/test/sampled_tests.jl passes all testsets: 24, 1, 2 (OffsetArrays) and 5. The same file passes with --check-bounds=yes, and under julia +1.10 (1.10.12).

  • Adversarial probes (scratch/probe.jl, all pass on the PR):

    • Rules and grids: random data, TrapezoidalRule and SimpsonsRule, on both a uniform range and a nonuniform Vector grid.
    • Array cases:
      • matrix with an offset non-integrated axis, dim=2;
      • matrix offset on both axes with an offset x;
      • offset vector with an offset x;
      • offset matrix with dim=1;
      • 3-D array with offset non-integrated axes;
      • the init/solve! caching path, including replacing cache.y.
    • Comparison: every result was compared (≈) against the same data on 1-based arrays with the same weight type.
    • One false alarm, worth recording: my first version compared Simpson on an offset collect(range) against Simpson on the range. Those produce different weight types, nonuniform composite 1/3 versus uniform 1/3–3/8, which are different rules. That was my error, not the PR's. With matching weight types, every case agrees.
  • Axis-mismatch inputs never reach the fast path (confirmed): 1-based data with an offset x, and offset data with a 1-based x, are both rejected by SciMLBase 3.57.0 with AssertionError: The integrand y must obey the same indexing as the sampling points x along the integrated dimension. So the 1-based-weights fast path cannot be reached with mismatched axes.

  • IdOffsetRange grid (confirmed): x = OffsetArrays.IdOffsetRange(1:9, -1) is an AbstractRange and so gets 1-based UniformWeights. With offset data, Trapezoid and Simpson match the 1-based result for both vector and matrix input.

  • Performance (scratch/bench.jl, evalrule minima, master vs PR):

    Case master PR
    4×10001 matrix 57.458 µs 57.500 µs
    vector 12.375 µs 12.416 µs
    nonuniform-weight matrix 57.000 µs 56.917 µs

    No regression, consistent with the PR body's claim. My absolute numbers differ from its figures because the hardware differs.

  • has_offset_axes is public: Base.ispublic(Base, :has_offset_axes) returns true on 1.12.

  • QA (confirmed): GROUP=QA julia +1.12 --project -e 'using Pkg; Pkg.test()', run on a copy of the PR head, gives QA/qa.jl | 27 27, Testing Integrals tests passed.

  • Formatting and spelling (confirmed): julia +1.12 --project=@runic -m Runic --check --diff src/sampled.jl test/sampled_tests.jl exits 0. git diff origin/master...HEAD | typos - exits 0.

What I did not verify

  • The full GROUP=Core suite. I ran only sampled_tests.jl, which is the only test file touched. The source change only affects arrays with offset axes, because has_offset_axes is constant-false for Array.
  • QA on Julia 1.10. ExplicitImports' public-access check does not apply there.
  • The Core 32-bit CI job.
  • GPU arrays: they still take the fast_scalar_indexing == false branch, which is unchanged.
  • A trial merge with Stabilize sampled eachslice integration with a function barrier #380 or SIMD the 1-D sampled-data weight×data reduction #381.
  • That master gives a BoundsError under --check-bounds=yes (a claim in the PR body). I only ran the PR side with check-bounds.

🤖 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)

@ChrisRackauckas
ChrisRackauckas marked this pull request as ready for review October 16, 2026 04:24
@ChrisRackauckas
ChrisRackauckas merged commit 49f8b10 into SciML:master Oct 16, 2026
15 checks passed
ChrisRackauckas referenced this pull request Oct 17, 2026
Bump version to 5.5.4
Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Agent-Harness: Claude Code 2.1.301
Agent-Model: claude-opus-5-5
Claude-Session: https://claude.ai/code/session_01VZjwf7Ux9YtTzwHGHMqnsq
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants