Skip to content

fill!(::AbstractVectorOfArray, x): fill inner arrays directly, fixing a silent ragged bug - #669

Merged
ChrisRackauckas merged 1 commit into
SciML:masterfrom
ChrisRackauckas-Claude:voa-fill-fastpath
Oct 10, 2026
Merged

ChrisRackauckas merged 1 commit into
SciML:masterfrom
ChrisRackauckas-Claude:voa-fill-fastpath

Conversation

@ChrisRackauckas-Claude

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

Copy link
Copy Markdown
Member

What changed and why

fill!(VA::AbstractVectorOfArray, x) looped for i in 1:length(VA.u) and read/wrote through VA[:, i]. For a ragged inner array, VA[:, i] (_getindex(A, ::NotSymbolic, ::Colon, I::Int)) builds and returns a zero-padded copy when the inner array's size doesn't match the ragged maximum — it is not the real storage in that case. fill!(VA[:, i], x) then fills that throwaway copy and discards it; the real VA.u[i] is never touched. In the common non-ragged case, VA[:, i] does return the real storage directly, so this only costs an unnecessary allocating round-trip through _getindex/checkbounds/size(A) on every iteration (no correctness bug there).

Fixed by iterating VA.u[i] directly, which is always the actual storage regardless of raggedness. This is both faster (no reconstruction) and fixes the silent correctness bug for ragged arrays described below.

Before/after

n = 1000 inner (non-ragged) Vector{Float64} of length 2:

before  fill! n=1000:  min  2430186.0 ns (2.43 ms)  allocs 0
after   fill! n=1000:  min     4848.4 ns (4.8 µs)   allocs 0

(matches the earlier perf audit's 2.43ms -> 5.0µs on 1.12 / 2.45ms -> 2.1µs on 1.10)

Ragged correctness — this is the important part. fill! on a ragged VectorOfArray is silently broken on current master:

R = VectorOfArray([[1.0, 2.0], [3.0], [4.0, 5.0, 6.0]])
fill!(R, 9.0)

Before (unmodified master, 02c477c):

julia> R.u
3-element Vector{Vector{Float64}}:
 [1.0, 2.0]         # NOT filled — write silently dropped
 [3.0]              # NOT filled — write silently dropped
 [9.0, 9.0, 9.0]     # the only inner array that actually matches the ragged max size, so it alone gets filled

After (this branch):

julia> R.u
3-element Vector{Vector{Float64}}:
 [9.0, 9.0]
 [9.0]
 [9.0, 9.0, 9.0]

Every inner array is filled, matching what fill! obviously should do and what every non-ragged caller already (correctly) observed.

Failing-before / passing-after

Added a test to test/Core/interface_tests.jl (next to the existing fill!/testva2 test):

ragged_fill = VectorOfArray([[1.0, 2.0], [3.0], [4.0, 5.0, 6.0]])
fill!(ragged_fill, 9.0)
@test ragged_fill.u == [[9.0, 9.0], [9.0], [9.0, 9.0, 9.0]]

Before (reverted src/vector_of_array.jl to unmodified master, kept the new test):

Test Failed at .../test/Core/interface_tests.jl:257
  Expression: ragged_fill.u == [[9.0, 9.0], [9.0], [9.0, 9.0, 9.0]]
   Evaluated: [[1.0, 2.0], [3.0], [9.0, 9.0, 9.0]] == [[9.0, 9.0], [9.0], [9.0, 9.0, 9.0]]
ERROR: LoadError: There was an error during testing

After (this branch): test/Core/interface_tests.jl runs to completion with no failures, on both Julia 1.10.12 and 1.12.7.

Full GROUP=Core run (both Julia versions): all testsets green, including the existing non-ragged fill! coverage in interface_tests.jl and utils_test.jl (immutable/mixed SVector/MVector fill!). GROUP=QA on 1.12.7: 20/20 pass. GROUP=QA on 1.10.12 hits the same pre-existing, unrelated Aqua ambiguity failure documented in #667 (confirmed present on unmodified master).

What was not verified

  • Downstream, GPU, NoPre and AD groups were not run; fill! on a DiffEqArray/AbstractVectorOfArray wrapping GPU arrays was not specifically exercised (the ismutable/AbstractVectorOfArray branch is unchanged logic, just reading ui once instead of VA[:, i] repeatedly, so this shouldn't interact with GPU scalar-indexing restrictions differently than before — but it wasn't measured).

Anything a reviewer should push back on

  • This is a behavior change for ragged arrays, not just a speedup: before, fill! on a ragged VectorOfArray silently filled only the longest inner array(s); after, it fills every inner array. I believe the old behavior was simply a bug (a fill! that silently no-ops on part of its target with no error and no indication), and the fix is what any caller would expect, but flagging explicitly since it changes observable output for an existing (if undocumented) code path.

Please ignore until reviewed by @ChrisRackauckas.

Risk assessment

  • Risk: low
  • Blast radius: src/vector_of_array.jl, one fill! method. No public API signature change. The behavior change is a correctness fix (ragged arrays now get fully filled instead of silently partially filled) rather than an API change — anything that depended on the old silent-no-op behavior for part of a ragged array was already relying on a bug.
  • Evidence: failing-before/passing-after shown above for the ragged case; before/after timing for the non-ragged case; GROUP=Core full pass on Julia 1.10 and 1.12, including all pre-existing fill! tests.
  • Independent review: pending
  • Merge: needs human review — the ragged-array behavior change (even though it looks like a straightforward bug fix) is worth a second opinion before merge rather than auto-merging.

🤖 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: #669 (comment)

fill! looped over VA[:, i], which for a ragged inner array builds and
returns a zero-padded *copy* (src/vector_of_array.jl's _getindex for
(::Colon, ::Int)). Filling that copy touches nothing: the write is
silently dropped for every inner array shorter than the ragged
maximum. It also costs an unnecessary allocating reconstruction on
every iteration even in the common non-ragged case.

Iterate VA.u[i] directly instead: it is always the real storage, so
fill! correctly fills every inner array regardless of its size, and
the non-ragged case no longer reconstructs anything.

n=1000 inner arrays: fill! 2.43ms -> 4.8us (1.12), 2.45ms -> 2.1us
(1.10, per the earlier audit). Ragged correctness (previously silently
wrong) verified: filling a ragged VectorOfArray([[1.,2.],[3.],[4.,5.,6.]])
with 9.0 now gives [[9.,9.],[9.],[9.,9.,9.]] (every inner array filled)
instead of the pre-existing [[1.,2.],[3.],[9.,9.,9.]] (only the longest
inner array actually gets filled).

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Claude <noreply@anthropic.com>
Agent-Harness: Claude Code 2.1.294
Agent-Model: claude-sonnet-5
Agent-Session: https://claude.ai/code/session_01LPHREnnonfLg1VcE1EJovv
@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

PR: #669 (head 61250c9). It has no linked issue and no tracking issue, so there was no sibling interface to check beyond the other open perf PRs (#666–#678). None of them touch fill!.

Blocking findings

None.

The change at src/vector_of_array.jl:1169-1187 reads VA.u[i] where the old code went through VA[:, i]. I read _getindex(::NotSymbolic, ::Colon, ::Int) at src/vector_of_array.jl:622-661. For non-ragged inner arrays, and for inner arrays whose ndims differs, VA[:, i] already returned A.u[I] itself. For a ragged inner array with the same ndims, it returned a freshly allocated zero-padded copy, and that copy is what got filled. So the new code does the same thing as the old code in the non-ragged case, and fixes the ragged case. In the scalar-element branch, VA[:, i] = x dispatched to setindex!(VA, v, ::Colon, I::Int), which is just VA.u[I] = v (src/vector_of_array.jl:870-875), so swapping in VA.u[i] = x changes nothing. All of this was confirmed by running, below.

Non-blocking findings

  1. Comment volume and phrasing (read from the diff). 8 of the 19 added lines are comments (about 42%). CLAUDE.md asks for roughly 10%. The comment in src (src/vector_of_array.jl:1171-1174) explains the code by contrasting it with the alternative ("instead of through VA[:, i]"). The test comment (test/Core/interface_tests.jl:251-254) says the same thing again and describes how the old code would fail. Suggested fix: cut the src comment to one line in the present tense, e.g. # VA.u[i] is the real storage; VA[:, i] zero-pads ragged columns into a copy, and cut the test comment to one line or delete it.

  2. Test coverage is narrower than the fix (confirmed by running). The new test only covers a ragged Vector{Vector{Float64}}. My probe showed master is also wrong for the cases below, and all of them pass on the PR head:

    • ragged DiffEqArray, which goes through its own _getindex method at line 644 (master: [[1,2],[3],[7,7,7]]; head: all 7s, t untouched);
    • ragged matrices [2×2, 1×3, 3×1], where master filled none of the inner arrays;
    • nested ragged VectorOfArray;
    • ragged MVector.

    It would be cheap to loop the test over these cases and assert that D.t is untouched. Not required for merge.

  3. The PR body misdescribes the non-ragged cost (confirmed by running). It calls the old non-ragged path an "unnecessary allocating round-trip". I measured 0 bytes allocated before and after. The slowdown comes from size(A) inside _getindex, which scans every inner array, so the old loop was O(n²). The speedup itself reproduces: n=1000 non-ragged, Julia 1.12, min of 20 runs: master 1.67 ms → head 1.67 µs, 0 allocations both times.

  4. Unlinked reference in the PR body. The body cites an "earlier perf audit" and gives no link to it.

What I ran

All runs used Julia 1.12, with TMPDIR and scratch envs under .../SciML_RecursiveArrayTools.jl-669/{tmp,scratch}. Master was extracted with git archive origin/master into scratch/master.

  • scratch/probe.jl against master and head (confirmed by running):

    Case master head
    ragged vec [[1,2],[3],[4,5,6]], fill!(·, 9.0) [[1,2],[3],[9,9,9]] all 9s
    ragged DiffEqArray only the longest filled all filled, t unchanged
    ragged matrices nothing filled all filled
    nested ragged VoA partly filled all filled
    ragged MVector partly filled all filled

    These cases came out the same on master and head:

    • SVector (immutable branch);
    • Any-eltype ragged SVector;
    • scalar-element VoA;
    • non-ragged Array(NR) == fill(...);
    • fill! returns VA.

    On head, R[:,2] after the fill gives [9.0, 0.0, 0.0]: padding zeros appear only outside the stored data. Aliasing ([a, a, [0.0]]) fills everything on head; on master the [0.0] entry is left as is.

  • GROUP=Core Pkg.test("RecursiveArrayTools") with head dev'd: passed. Interface Tests 154/154, Indexing 180 pass + 2 pre-existing @test_broken (basic_indexing.jl is not touched by this PR). Log: scratch/core_head.log.

  • Fail-before check: master src plus the PR's test/Core/interface_tests.jl, GROUP=Core. Failed at interface_tests.jl:257 with Evaluated: [[1.0, 2.0], [3.0], [9.0, 9.0, 9.0]] == [[9.0, 9.0], [9.0], [9.0, 9.0, 9.0]], Interface Tests 153 pass / 1 fail. This is the output the PR body claims. Log: scratch/core_master_newtest.log.

  • julia +1.12 --project=@runic -m Runic --check --diff on both changed files: clean. typos on the diff: clean.

  • gh pr checks 669: everything green except SciMLSensitivity Downstream Core1. That job fails with a Mooncake MethodError: no constructors have been defined for Any in _build_output_tangent_cartesian (concrete_solve_derivatives.jl:451). PRs Fix two silent memory-safety bugs in recursivecopy!/copyat_or_push! #667 and Add O(1) fast path to AbstractVectorOfArray getindex #668 fail Core1 the same way, so it is not caused by this PR.

  • No new dependencies, exports or public-API signatures; license is unchanged.

What I did not verify

  • I did not run Julia 1.10/LTS locally. CI's Core (lts) job passed.
  • I did not run the QA, GPU, AD, NoPre or Downstream groups locally. All of them passed on CI except the unrelated Core1 failure above.
  • I did not test GPU-backed inner arrays myself. The new code does less indexing than the old (no VA[:, i], no size(A)), so I expect no new scalar-indexing risk, but that comes from reading the code, not from running it.
  • I did not separately recheck the 1.10 timing claim (2.45 ms → 2.1 µs).

🤖 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 10, 2026 18:44
@ChrisRackauckas
ChrisRackauckas merged commit 2234b23 into SciML:master Oct 10, 2026
45 of 46 checks passed
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