Skip to content

Generate FullSpecialize initialization maps with static buffers - #5045

Merged
ChrisRackauckas merged 12 commits into
SciML:masterfrom
ChrisRackauckas-Claude:agent/full-specialize-initialization-maps
Sep 13, 2026
Merged

Generate FullSpecialize initialization maps with static buffers#5045
ChrisRackauckas merged 12 commits into
SciML:masterfrom
ChrisRackauckas-Claude:agent/full-specialize-initialization-maps

Conversation

@ChrisRackauckas-Claude

@ChrisRackauckas-Claude ChrisRackauckas-Claude commented Aug 30, 2026

Copy link
Copy Markdown
Member

Scope

Generate the state and parameter initialization maps as individual runtime-generated functions only for explicitly requested FullSpecialize. Preserve the existing default host-map path. Static buffers reduce warm map allocations; this is not an in-place preallocation API.

The implementation shares one buffer constructor and one reverse-mode rule. It preserves mutable parameter buffers, copies cache buffers instead of aliasing them, handles split=false, and retains array-valued callable parameters. Immutable static constructors produce isbits results in the tested case; mutable non-isbits buffers use SizedVector storage.

Addresses #5043. Related GPU workaround: SciML/DiffEqGPU.jl#516.

Review round (@AayushSabharwal, 2026-09-09)

Seven inline comments, all addressed in fed8cdb77b, and each answered on its own thread.

1. gensym — this was a real defect, not style. The five gensymed locals in the new codegen are now fixed sentinels (INITMAP_VALUES, INITMAP_ELTYPE, INITMAP_SOLUTION, INITMAP_PROBLEM, INITMAP_OUTER_PARAMETERS), following the rationale already recorded for HOMOTOPY_LAMBDA in homotopy_operator.jl: a gensym embeds a process-global counter, so the same system lowers to a different Expr every time and the RuntimeGeneratedFunctions Expr-hash cache never hits.

This is the cause of the generated-map type churn the earlier revision of this body reported as an open problem.

Failing on the unfixed branch (git stash of the fix only, same probe, same session):

generated map type reuse across identical systems: Test Failed
  Expression: typeof(d1.initializeprobpmap) === typeof(d2.initializeprobpmap)
   Evaluated: RuntimeGeneratedFunction{(Symbol("##problem#280"), Symbol("##initialization_solution#281")), …, (0x8298bbd3, 0xbcbc4d5b, 0x7a096718, 0x43533499, 0x57e4022b), Nothing} ===
              RuntimeGeneratedFunction{(Symbol("##problem#287"), Symbol("##initialization_solution#288")), …, (0x44632bf6, 0x42fe7ec6, 0x4d5618c1, 0xb6d83a49, 0x2afa037a), Nothing}
Test Summary:                                     | Pass  Fail  Total   Time
generated map type reuse across identical systems |    2     2      4  10.6s

Passing with the fix:

Test Summary:                                     | Pass  Total  Time
generated map type reuse across identical systems |    4      4  7.0s

Checked in as three assertions at the end of the FullSpecialize initialization maps are generated functions testset, which is why that group's count moved from 698 to 701.

2-6. Splats in Expr construction. _generated_map_expr, _parameter_buffer_expr, the discrete block sizes, the portion tuple, and the MTKParameters call all build their Exprs by push!/append! onto .args now. _static_initialization_buffer's T[values...] became collect(T, values). No splats remain in the added code other than kwargs... forwarding.

7. groups as a tuple of tuples. Now a Vector{Vector{Vector{SymbolicT}}} built with type asserts against ReorderedParametersT's Union element type, so both the flat_syms loop and the buffer loop are concretely typed instead of specializing on each system's shape. Nothing indexed groups heterogeneously — it is only iterated alongside the field names — so the tuple bought nothing.

Rebased onto 35818b8328. The CasADi compat bump and the Optimization runner pin that used to show in this diff have landed on master separately and are gone from it.

Second review round (@AayushSabharwal, 2026-09-11 — approved)

One inline comment: pass an options struct instead of kwargs.... Done in
16ef998c28 — both map builders take a
GeneratedFunctionOptions positionally and call the struct method of
build_explicit_observed_function, with _fullspecialize_map_options deriving it from
the problem's own opts.fn_opts.codegen.

Two behaviour fixes fall out of it. The old kwargs... path went through the keyword
compat shim, which built a fresh CodegenFunctionOptions from whatever was left in
the bag — so the maps were not inheriting the problem's checkbounds/cse. And
eval_expression is now honoured; previously the builders never received it, so the maps
were unconditionally RuntimeGeneratedFunctions even with eval_expression = true,
unlike the default map path. Both are called out on the thread in case the second is not
wanted.

Local verification

Julia 1.12.7, on 16ef998c28. JULIA_NUM_PRECOMPILE_TASKS=1, --heap-size-hint=8G.

julia +1.12 --project=<runicenv> -m Runic --check <changed Julia files>   # exit 0
git diff upstream/master --check                                          # exit 0
git diff upstream/master --unified=0 | sed -n '/^+++ /d; s/^+//p' | typos -  # exit 0
GROUP=Initialization
  Guess Propagation         |  11 pass |  11 total | 3m13.3s
  InitializationSystem Test | 706 pass | 12 broken | 718 total | 20m43.6s
  Initial Values Test       |  65 pass |  65 total | 6m00.6s     (0 errors)

GROUP=Extensions
  HomotopyContinuation Extension Test |  76 pass |  76 total | 10m07.9s
  LabelledArrays Test                 |  18 pass |  18 total | 31.1s
  BifurcationKit Extension Test       |   7 pass |   7 total | 1m25.4s
  Despecialized MTKParameters AD Test |   2 pass |   2 total | 7m07.9s
  Initialization maps AD              |   7 pass |   7 total | 17.2s

GROUP=FMI
  FMI | 60 pass | 60 total | 15m35.4s

All three groups exited 0. The 12 broken tests in InitializationSystem Test are
pre-existing @test_brokens on master; none were added, removed, skipped or loosened.

The MTKStdlib reproducer below also runs clean against this commit (solveremake
solve), outside CI.

Allocation measurement and limits

A warmed, two-state initialization example without user parameters measured allocations for both map calls together, taking the minimum of ten calls behind a function barrier:

FullSpecialize map configuration Bytes/call pair
Existing host maps 784
Generated maps, mutable buffers 80
Generated maps, immutable static constructors 0

This is a map microbenchmark, not a solver benchmark or a general zero-allocation guarantee.

On TTFX. The generated map types are now reused across identical systems, which is what the earlier revision of this body reported as an unresolved code-generation/cache concern. That measurement was taken before the gensym fix and no longer describes the branch. I have not re-run the fresh-process TTFX comparison since the fix, so this PR still makes no TTFX claim in either direction — only that the specific type-churn defect behind it is gone and now has a regression guarding it. The default path remains unchanged.

What was not verified

  • GPU execution. This machine has no GPU. The static/isbits CPU assertions are not a substitute for a device execution test.
  • Fresh-process TTFX, as above.
  • GROUP=QA is red, from a pre-existing master failure. Rerun on the rebased branch: targeted JET Tests pass 54/54, Aqua Tests are 20 passed / 1 failed / 21, and the single failure is the package-wide JET typo-mode scan reporting 234 possible errors. Every one of the 234 is the known ##And# (182) / ##Call# (52) Moshi @match codegen category tracked in ModelingToolkitBase QA lane red since #4832: JET typo mode reports 263 errors, 234 of them from Moshi @match codegen #5063, which is still open after Fix JET typo-mode findings in ModelingToolkitBase #5065 cleared the non-Moshi ones. None of the 234 diagnostics point anywhere inside the code this PR adds (problem_utils.jl:1636-1810). No QA suppression or dependency override is included here.
  • The StructuralIdentifiability.jl/Core/1/ downstream job fails, and it is not this PR. It passed on the pre-rebase head 438703fbe4 and fails on clean master 35818b8328 — the exact base this branch was rebased onto. It regressed on master between f21201f889 (pass) and 9b458b059d (fail); a separate bisect is running and will get its own issue. Nothing in this diff touches it.
  • Not rerun since the rebase: LTS/prerelease Julia, FMI, Optimization, SymbolicIndexingInterface, the neural-network downstream suite, and the full downstream matrix. Their pre-rebase results are listed above.
  • Docs build — no docs, docstrings, or public API are changed.

An initial source-loading test run exposed an independently reproduced LinearSolve world-age failure, reported separately at SciML/LinearSolve.jl#1282; no LinearSolve workaround is included here.

Anything a reviewer should push back on

  • The sentinel names (__mtk_initialization_*) are module-level consts in ModelingToolkitBase. They are internal and unexported, but they are five new names in the module namespace where a single NamedTuple or a naming function would also work.
  • _static_initialization_buffer's collect(T, values) change was not one of the review comments; it is a splat in the same function and I removed it for consistency. It is on the runtime path (the non-isbits SizedVector branch), covered by the Any[Ref(1)] assertion in the Initialization group.
  • The groups rewrite carries the same behaviour as before for an empty parameters(sys; initial_parameters = true): grouped[1] would throw, because reorder_parameters returns early on empty ps. I did not add a guard, since that would be a behaviour change beyond the review, but it is reachable in principle.
  • @ChrisRackauckas asked for the gensym rule to be written down. Add AGENTS.md with code generation conventions #5105 adds a root AGENTS.md (the repo had none) with that rule plus the splat and container-typing ones. Kept separate so this diff stays behaviour-only.

Buffer types: only isbits buffers go static

FullSpecialize was turning every parameter buffer into a StaticArray, heap-element
ones included. That is what broke ModelingToolkitStandardLibrary.jl/Core/1/, and the
SizedVector branch added earlier for FMI was aimed at the wrong buffer. Fixed in
b8c603d04a.

The FMI setindex! failure is not about size — MArray's setindex! bounds-checks first,
then hard-error()s on any non-isbits eltype, because it mutates an inline NTuple
through a raw pointer (StaticArrays #27). The buffer was an array-valued discrete,
MVector{1, Vector{Float64}}, written by an FMI imperative_affect:

[3] setindex!  @ BlockArrays/blockedarray.jl:249 [inlined]
[4] set_parameter!(p::MTKParameters{…, Tuple{BlockedVector{Vector{Float64},
      MVector{1, Vector{Float64}}, …}}, …},
      val::Vector{Float64}, pidx::ParameterIndex{Discrete, Tuple{Int64, Int64}})

Making such a buffer static buys nothing: the elements are still pointers, so an
MVector{1, Vector{Float64}} is no more GPU-resident than a Vector{Vector{Float64}}.
Gating on isbitstype therefore costs nothing that matters, and every buffer that could
be GPU-resident stays an MVector. Same FMI model, only spec changed:

buffer AutoSpecialize FullSpecialize after this commit
tunable Vector{Float64} MVector{7, Float64}
initials Vector{Float64} MVector{8, Float64}
discrete block sizes Vector{Int64} SVector{1, Int64}
discrete storage Vector{Vector{Float64}} Vector{Vector{Float64}}
nonnumeric Vector{FMI2CSFunctor} Vector{FMI2CSFunctor}

The gate has to go in two places or the generated map and MTKParametersReconstructor
disagree on the buffer type: _static_initialization_buffer, and the get_p_constructor
closure the reconstructor routes through. get_p_constructor returns identity unless
pType <: StaticArray, and this PR is what makes pType static — so it only started
applying create_array to non-isbits buffers because of this PR.

Each variant measured separately, one change at a time:

non-isbits buffer MTKStdlib remake FMI
MVector passes setindex! error
SizedVector (previous revision) StackOverflowError 60/60
Vector, gate in one place typeassert, Vector vs MVector 60/60
Vector, gate in both passes 60/60

Failing before / passing after, same testset, fix stashed and restored in one session:

before:  FullSpecialize initialization maps are generated functions | 38 pass | 2 fail | 40 total
after:   FullSpecialize initialization maps are generated functions | 40 pass |  40 total | 26.8s

The two new assertions pin both halves — non-isbits buffers stay on the heap, and
tunable/initials stay StaticArray so the gate cannot silently over-apply. SizedVector
is now unused and its import is dropped.

Ignore this draft until reviewed by @ChrisRackauckas.


Commits through 9a652005c8 were produced with Codex CLI 0.153.4 (model gpt-6-astra), local session 01a070db-f817-7f93-9295-aae2232fc2c7 (no shareable conversation URL exposed).

The review-response commit fed8cdb77b was produced with 🤖 Claude Code 2.0.14 (model: claude-opus-5[1m])

https://claude.ai/code/session_01GdSpCLd7NBZuuePJmcDzU7

@ChrisRackauckas-Claude

Copy link
Copy Markdown
Member Author

Plan/status: generate single dropped-body RuntimeGeneratedFunctions only for explicit FullSpecialize; preserve AutoSpecialize and AutoDespecialize behavior; verify map structure and output against existing maps; run Initialization and local hygiene checks; track the independently reproduced upstream/master QA failures in the required parallel investigation before final handoff.

@ChrisRackauckas-Claude

Copy link
Copy Markdown
Member Author

CI Runic failure is pre-existing formatter-version drift in three files untouched by this PR. The mechanical cleanup is isolated in draft PR #5046. This feature PR should be rebased after that cleanup lands; no formatting-only changes will be mixed into this behavior diff.

@ChrisRackauckas-Claude

Copy link
Copy Markdown
Member Author

Follow-up: the isolated formatting PR #5046 now passes CI Runic, Runic Suggestions, and spelling. This confirms the formatter failure here is resolved by that base cleanup rather than a change to this feature diff.

@ChrisRackauckas-Claude

Copy link
Copy Markdown
Member Author

Follow-up commit b3256600f9d2f0dae9ebdbb25f0841ff375a03c8 addresses the CI-discovered DiffCache path. The original run failed with the explicit unsupported-parameter ArgumentError; removing that guard exposed an existing immutable nonnumeric-buffer typo (fieldtype(C, i) while reconstructing N) and produced 66 passes / 2 errors in IndexCache. Correcting it to fieldtype(N, i) produced IndexCache Test | 68 pass | 68 total locally. The broader InterfaceII run later hit two unrelated stochastic jump assertions after 4,202 passes; clean-base reproduction/bisect is running separately.

@ChrisRackauckas-Claude
ChrisRackauckas-Claude force-pushed the agent/full-specialize-initialization-maps branch from b325660 to d3c8a49 Compare August 31, 2026 04:02
@ChrisRackauckas-Claude

Copy link
Copy Markdown
Member Author

Added array-valued callable-parameter support in d3c8a49.

Failure on the preceding implementation (b3256600f9):

ArgumentError: FullSpecialize initialization maps do not yet support array-valued parameter storage slots.
Test Summary: FullSpecialize array-valued storage slot | 1 Error

The same namespaced ArraySymbolic reproduction now passes:

Test Summary: array callable | 3 pass | 3 total | 14.9s

The complete ModelingToolkitNeuralNets suite also passes against this branch (10 + 1 + 37 + 55 + 7 passing assertions, with its 2 pre-existing broken tests). This directly covers the lux_apply parameter shape that failed in https://github.com/SciML/ModelingToolkit.jl/actions/runs/33300901287/job/99228760344.

@ChrisRackauckas-Claude

Copy link
Copy Markdown
Member Author

CI exposed a mutable-buffer regression in the first generated-parameter implementation: https://github.com/SciML/ModelingToolkit.jl/actions/runs/33355753893/job/99377920719 failed six MTKParameters cases because SciMLStructures.replace! attempted setindex! on generated SVector tunable/constant buffers.

Fixed in 93854e0 by selecting MVector when the source buffer is mutable while retaining SVector for immutable prototypes. The generated map itself remains an isbits, body-dropped RuntimeGeneratedFunction.

Local verification after the fix:

Focused Tunable/Constants replace! reproducer: exit 0
SymbolicIndexingInterface | 1984 pass | 1984 total | 8m06.7s
Testing ModelingToolkit tests passed
Runic.jl v1.10, typos, git diff --check: exit 0

@ChrisRackauckas-Claude

Copy link
Copy Markdown
Member Author

Commit 613bdd1 fixes the non-isbits mutable-buffer regression found by the FMI LTS job: https://github.com/SciML/ModelingToolkit.jl/actions/runs/33358934228/job/99386604733.

Failing before (CI):

FMI | 56 pass | 2 errors | 58 total
setindex!() with non-isbitstype eltype is not supported by StaticArrays. Consider using SizedArray.

Passing after (local):

GROUP=FMI ... Pkg.test()
FMI | 60 pass | 60 total | 15m53.9s
Testing ModelingToolkit tests passed

GROUP=Initialization ... Pkg.test()
Initialization | 870 pass | 12 broken | 882 total | 36m22.1s
Testing ModelingToolkit tests passed

Runic 1.10, typos over the diff, and git diff --check also exited successfully.

@ChrisRackauckas-Claude

Copy link
Copy Markdown
Member Author

Commit 010bc30 addresses the downgraded Initialization failure: https://github.com/SciML/ModelingToolkit.jl/actions/runs/33366256108/job/99407718377.

The failed job used SymbolicIndexingInterface 0.3.46, where the callable parameter is classified as ScalarSymbolic despite retaining its declared [1:2] shape. The generated map already returned the exact callable; the failing enum assertion encoded a newer dependency detail.

Local Julia 1.10 reproduction pinned the relevant CI floors (ArrayInterface 7.28.0, SciMLBase 3.48.0, StaticArraysCore 1.4.3, SymbolicUtils 4.37.0, SymbolicIndexingInterface 0.3.46, Symbolics 7.37.0):

classifier=ScalarSymbolic
shape/generated-map/callable-identity assertions: 3/3 pass

The same semantic assertions also pass 3/3 on the current dependency stack. Runic 1.10, typos, and git diff --check pass for the follow-up diff.

@ChrisRackauckas-Claude

Copy link
Copy Markdown
Member Author

Follow-up commit 258f72c fixes the remaining downgrade error from https://github.com/SciML/ModelingToolkit.jl/actions/runs/33387712747/job/99474377835. SymbolicUtils 4.37 provides shape at its owning module but does not re-export it through ModelingToolkitBase.

Exact CI-floor local result after using the owner API:

classifier=ScalarSymbolic
SU.shape/generated-map/callable-identity assertions: 3/3 pass

The current-stack owner-API shape assertions pass 2/2. Runic 1.10, typos, and git diff --check pass.

@ChrisRackauckas-Claude

Copy link
Copy Markdown
Member Author

CI confirmation on head 258f72c:

At this update, 38 checks pass, 66 remain in progress, 2 known base checks fail, and 1 is skipped. No completed branch-relevant check fails.

@ChrisRackauckas-Claude ChrisRackauckas-Claude changed the title Generate device-compatible FullSpecialize initialization maps Generate FullSpecialize initialization maps with static buffers Sep 5, 2026
@ChrisRackauckas-Claude
ChrisRackauckas-Claude force-pushed the agent/full-specialize-initialization-maps branch from 258f72c to 66a175d Compare September 5, 2026 13:38
@ChrisRackauckas-Claude

Copy link
Copy Markdown
Member Author

Updated the draft with a shared static-buffer constructor and reverse-mode rule, independent cache ownership, and split=false support. Added discriminating regressions for cache aliasing, unsplit parameters, and state/parameter gradients. Final local results: root QA 40/40; Initialization 882 passed plus 12 existing broken; SymbolicIndexingInterface 1984/1984; focused map/AD tests 39/39. Full Extensions and neural-network Core also passed as detailed in the body.

The allocation microbenchmark improved from 784 to 80 bytes per map-call pair, or zero with immutable static constructors. This is not a TTFX fix: repeated FullSpecialize construction still creates distinct generated map types. The title and body now state the narrower scope.

Registered-dependency Base QA remains blocked by the independently reproduced master Moshi/JET issue: #5063 . No QA suppression or dependency workaround is included. Ignore this draft until reviewed by @ChrisRackauckas.

AI update: Codex CLI 0.153.4; model gpt-6-astra; local session 01a070db-f817-7f93-9295-aae2232fc2c7 (no shareable conversation URL exposed).

@ChrisRackauckas-Claude
ChrisRackauckas-Claude force-pushed the agent/full-specialize-initialization-maps branch from 66a175d to 438703f Compare September 6, 2026 06:53
for (arg, source) in zip(generated_args, sources)
push!(body.args, :(local $arg = $source))
end
raw = gensym(:initialization_map_values)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I thought we want to avoid gensym?

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah we do want to avoid gensym, good catch, should probably agents.md that.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fixed in fed8cdb77b. All five gensyms in the new codegen are now fixed sentinels (INITMAP_VALUES, INITMAP_ELTYPE, INITMAP_SOLUTION, INITMAP_PROBLEM, INITMAP_OUTER_PARAMETERS), following the rationale already written down for HOMOTOPY_LAMBDA in homotopy_operator.jl.

This was not just style — it was a real defect in the PR, and it explains the TTFX non-improvement the body reported. Building the same system twice produced two distinct RuntimeGeneratedFunction types, so the RGF Expr-hash cache never hit:

generated map type reuse across identical systems: Test Failed
  Expression: typeof(d1.initializeprobpmap) === typeof(d2.initializeprobpmap)
   Evaluated: RuntimeGeneratedFunction{(Symbol("##problem#280"), Symbol("##initialization_solution#281")), …, (0x8298bbd3, 0xbcbc4d5b, 0x7a096718, 0x43533499, 0x57e4022b), Nothing} ===
              RuntimeGeneratedFunction{(Symbol("##problem#287"), Symbol("##initialization_solution#288")), …, (0x44632bf6, 0x42fe7ec6, 0x4d5618c1, 0xb6d83a49, 0x2afa037a), Nothing}
Test Summary:                                     | Pass  Fail  Total   Time
generated map type reuse across identical systems |    2     2      4  10.6s

With the sentinels:

Test Summary:                                     | Pass  Total  Time
generated map type reuse across identical systems |    4      4  7.0s

Checked in as three assertions at the end of the FullSpecialize initialization maps are generated functions testset, so the two lowerings are required to share a type.

@ChrisRackauckas — written up in #5105, which adds a root AGENTS.md (the repo had none) with this rule plus the two below.

raw = gensym(:initialization_map_values)
push!(body.args, :(local $raw = $(expr.args[2])))
push!(body.args, finish(raw))
return Expr(:function, Expr(:tuple, map_args...), body)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Shouldn't splat. Prefer

Suggested change
return Expr(:function, Expr(:tuple, map_args...), body)
fn_args = Expr(:tuple)
append!(fn_args.args, map_args)
return Expr(:function, fn_args, body)

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Applied your suggestion verbatim in fed8cdb77b.

end

function _parameter_buffer_expr(prototype, raw::Symbol, idxs, p_constructor)
values = Expr(:tuple, [Expr(:ref, raw, i) for i in idxs]...)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Shouldn't splat

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fixed. _parameter_buffer_expr now builds the tuple by pushing:

values = Expr(:tuple)
for i in idxs
    push!(values.args, Expr(:ref, raw, i))
end

The one other splat in that function, SizedVector{length(values), T}(T[values...]) in _static_initialization_buffer, is now collect(T, values). There are no splats left in the added code other than kwargs... forwarding.

compiler_options::CompilerOptions = CompilerOptions(), kwargs...
)
ps = parameters(sys; initial_parameters = true)
groups = if is_split(sys)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why does groups need to be a tuple? It can just be an appropriately typed vector? This is catastrophic for inference.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Changed in fed8cdb77b. It does not need to be a tuple — nothing indexes groups heterogeneously, it is only iterated alongside the field names. It is now a Vector{Vector{Vector{SymbolicT}}}:

groups = Vector{Vector{SymbolicT}}[]
if is_split(sys)
    grouped = reorder_parameters(sys, ps; flatten = false)
    push!(groups, Vector{SymbolicT}[grouped[1]::Vector{SymbolicT}])
    initial_syms = _unwrap_initial_symbols!(copy(grouped[2]::Vector{SymbolicT}), initsys)
    push!(groups, Vector{SymbolicT}[initial_syms])
    for i in 3:5
        push!(groups, grouped[i]::Vector{Vector{SymbolicT}})
    end
else
    push!(groups, Vector{SymbolicT}[ps])
end

The type asserts pin the Union element type coming out of ReorderedParametersT, so both the flat_syms loop and the buffer loop below it are now concretely typed instead of specializing on each system's shape.

sizes = get_index_cache(sys).discrete_buffer_sizes[buffer_idx]
block_sizes = Expr(
:call, Expr(:curly, SVector, length(sizes), Int),
map(size -> size.length, sizes)...

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Avoid the splat

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fixed — block_sizes starts as Expr(:call, Expr(:curly, SVector, length(sizes), Int)) and the lengths are pushed onto .args in a loop.

push!(buffers, buffer)
end
portion = field === :tunable || field === :initials ? only(buffers) :
Expr(:tuple, buffers...)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Avoid the splat

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fixed — the non-tunable/initials portion is now built with Expr(:tuple) + append!(tup.args, buffers).

push!(portions, portion)
end
result = Expr(
:call, MTKParameters, portions...,

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Avoid the splat

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fixed — result = Expr(:call, MTKParameters), then append!(result.args, portions) and push! for the caches argument.

ChrisRackauckas and others added 10 commits September 10, 2026 08:26
Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: OpenAI Codex <noreply@openai.com>
Agent-Harness: Codex (version unknown)
Agent-Model: unknown
Agent-Session: local session ID 01a04feb-43da-7ce1-99e4-a68438c69833
Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: OpenAI Codex <noreply@openai.com>
Agent-Harness: Codex (version unknown)
Agent-Model: unknown
Agent-Session: local session ID 01a04feb-43da-7ce1-99e4-a68438c69833
Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: OpenAI Codex <noreply@openai.com>
Agent-Harness: Codex (version unknown)
Agent-Model: unknown
Agent-Session: local session ID 01a04feb-43da-7ce1-99e4-a68438c69833
Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: OpenAI Codex <noreply@openai.com>
Agent-Harness: Codex (version unknown)
Agent-Model: unknown
Agent-Session: local session ID 01a04feb-43da-7ce1-99e4-a68438c69833
Use SizedVector-backed storage when a mutable FullSpecialize parameter buffer contains non-isbits values, while retaining MVector for isbits element types.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: OpenAI Codex <noreply@openai.com>
Agent-Harness: Codex (version unknown)
Agent-Model: unknown
Agent-Session: local session ID 01a04feb-43da-7ce1-99e4-a68438c69833
Assert the declared array shape instead of the SymbolicIndexingInterface classification enum, which differs across supported dependency versions while the generated map behavior remains identical.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: OpenAI Codex <noreply@openai.com>
Agent-Harness: Codex (version unknown)
Agent-Model: unknown
Agent-Session: local session ID 01a04feb-43da-7ce1-99e4-a68438c69833
Read the callable parameter shape from SymbolicUtils directly so the regression works across supported versions that do not re-export shape through ModelingToolkitBase.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: OpenAI Codex <noreply@openai.com>
Agent-Harness: Codex (version unknown)
Agent-Model: unknown
Agent-Session: local session ID 01a04feb-43da-7ce1-99e4-a68438c69833
Share the state and parameter buffer constructor and give it a shape-preserving
pullback. Handle unsplit parameter storage before indexing split groups, copy
mutable cache buffers, and consolidate parameter reconstruction.

Cover cache ownership, nonempty mutable parameter portions, split and unsplit
static construction, and reverse-mode differentiation in the test suites.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: OpenAI Codex <noreply@openai.com>
Agent-Harness: Codex CLI 0.153.4
Agent-Model: gpt-6-astra
Agent-Session: local transcript /home/crackauc/.codex/sessions/2026/09/05/rollout-2026-09-05T05-17-39-01a070db-f817-7f93-9295-aae2232fc2c7.jsonl
Compare the generated parameter-map pullback with ForwardDiff and the expected
nonzero derivative, alongside the generated state-map gradient check.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: OpenAI Codex <noreply@openai.com>
Agent-Harness: Codex CLI 0.153.4
Agent-Model: gpt-6-astra
Agent-Session: local transcript /home/crackauc/.codex/sessions/2026/09/05/rollout-2026-09-05T05-17-39-01a070db-f817-7f93-9295-aae2232fc2c7.jsonl
Names embedded in the generated initialization maps are now fixed sentinels
instead of `gensym`s. A `gensym` embeds a process-global counter, so lowering
the same system twice produced two distinct `RuntimeGeneratedFunction` types
and never hit the RGF Expr-hash cache. Added a regression asserting the two
lowerings share a type.

`Expr`s are built by pushing onto `args` rather than splatting, and the
parameter groups are a concretely typed `Vector{Vector{Vector{SymbolicT}}}`
rather than a tuple of tuples, which forced the buffer loop to specialize on
each system's shape.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Agent-Harness: Claude Code 2.0.14
Agent-Model: claude-opus-5[1m]
Agent-Session: https://claude.ai/code/session_01GdSpCLd7NBZuuePJmcDzU7
Claude-Session: https://claude.ai/code/session_01GdSpCLd7NBZuuePJmcDzU7
@ChrisRackauckas-Claude
ChrisRackauckas-Claude force-pushed the agent/full-specialize-initialization-maps branch from 438703f to fed8cdb Compare September 10, 2026 13:12
@ChrisRackauckas-Claude

Copy link
Copy Markdown
Member Author

The rebase onto master picked up a new downstream failure that is this PR's bug, not a master one. Flagging it rather than fixing it, because the fix is a design call for this PR.

ModelingToolkitStandardLibrary.jl/Core/1/ fails: https://github.com/SciML/ModelingToolkit.jl/actions/runs/34481227325/job/102884238601

LinearInterpolation: Error During Test at downstream/test/sources.jl:652
  Got exception outside of a @test
  StackOverflowError:
   [1] StaticArraysCore.SizedVector{1, UnionAll, Vector{UnionAll}}(a::Vector{Any})
       @ StaticArrays ~/.julia/packages/StaticArrays/k7bbZ/src/SizedArray.jl:15
   repeated 79984 times

It is not pre-existing on master: the same job passes on #5105, which is master plus a single markdown file. It passed on this branch's pre-rebase head 438703fbe4 only because the test did not exist yet — the earlier job log has zero occurrences of ParametrizedInterpolation. The new MTKStdlib test builds ODEProblem{true, SciMLBase.FullSpecialize}, which is exactly the path this PR adds, and then remakes it.

Reproduced locally against this branch, outside CI:

built prob
nonnumeric buffer type: Tuple{StaticArraysCore.SizedVector{1, UnionAll, Vector{UnionAll}}}
solved
ERROR: LoadError: StackOverflowError:
Stacktrace:
 [1] StaticArraysCore.SizedVector{1, UnionAll, Vector{UnionAll}}(a::Vector{Any}) (repeats 79984 times)
   @ StaticArrays ~/.julia/packages/StaticArrays/k7bbZ/src/SizedArray.jl:15

Reduced to StaticArrays alone, no MTK involved:

using StaticArrays
struct Foo{T} end
SizedVector{1, UnionAll, Vector{UnionAll}}(Any[Foo])   # StackOverflowError

Mechanism. _static_initialization_buffer's third branch wraps a mutable non-isbits buffer in a SizedVector. Nonnumeric parameter buffers are always non-isbits, so a nonnumeric buffer holding a UnionAll (here DataInterpolations.LinearInterpolation) becomes SizedVector{1, UnionAll, Vector{UnionAll}}. remake then reconstructs it through similar_type(fieldtype(N, i), nonnumericT[i])(nonnumerics[i]) in parameter_buffer.jl. similar_type is a no-op on that fully-parameterized SizedVector, so the constructor is re-entered with the same type and a Vector argument forever.

That SizedVector branch came from 613bdd13 to fix the FMI setindex!-on-non-isbits failure, so it is load-bearing — this is not a one-line revert.

Why I did not fix it here. I tried the obvious fix (return a plain mutable Vector{T} for the non-isbits branch, and drop the now-unused SizedVector import). That does clear the StackOverflowError — the buffer becomes Tuple{Vector{UnionAll}} and the solve succeeds — but remake then fails one layer down:

ERROR: LoadError: TypeError: in typeassert, expected Tuple{Vector{UnionAll}}, got a value of type Tuple{StaticArraysCore.MVector{1, UnionAll}}
  [1] (::ModelingToolkitBase.MTKParametersReconstructor{...})

The generated map picks the nonnumeric buffer type from the outer problem's prototype, while MTKParametersReconstructor rebuilds it from buftypes captured off the initialization problem, and the two only agreed by accident before. Making them agree is a decision about which prototype the generated map should follow, so I have left the branch at the reviewed state rather than shipping a change that swaps one error for another.

Three buffer choices, three failures, for whoever picks this up:

non-isbits nonnumeric buffer fails as
MVector FMI: setindex! unsupported for non-isbits eltype
SizedVector (current) MTKStdlib: remake StackOverflowError
Vector MTKStdlib: remake typeassert, Vector vs MVector

Standalone reproducer (needs MTKStdlib + DataInterpolations dev-ed against this branch):

using ModelingToolkitStandardLibrary.Blocks
using ModelingToolkit, DataInterpolations, OrdinaryDiffEq, SciMLBase
using ModelingToolkit: t_nounits as t, D_nounits as D

@variables y(t) = 0
u = rand(15); x = 0:14.0
@named i = ParametrizedInterpolation(LinearInterpolation, u, x)
@named model = System([i.input.u ~ t, D(y) ~ i.output.u], t, systems = [i])
sys = mtkcompile(model)

prob = ODEProblem{true, SciMLBase.FullSpecialize}(sys, [], (0.0, 4))
solve(prob, Tsit5())                              # fine
remake(prob, p = [i.data => ones(15)])            # StackOverflowError

The other three red checks on this PR (NeuralPDE.jl/NNPDE1, SciMLSensitivity.jl/Core8, StructuralIdentifiability.jl/Core) all fail identically on the master-only #5105, so those are not this PR. A separate bisect for the StructuralIdentifiability master regression is running.

ChrisRackauckas added a commit that referenced this pull request Sep 10, 2026
Records three conventions that came up in review of #5045: no `gensym` in
names that reach generated code (it defeats the RuntimeGeneratedFunctions
Expr-hash cache), build `Expr`s by pushing onto `args` instead of splatting,
and keep codegen-side containers concretely typed vectors rather than tuples.



Agent-Harness: Claude Code 2.0.14
Agent-Model: claude-opus-5[1m]
Agent-Session: https://claude.ai/code/session_01GdSpCLd7NBZuuePJmcDzU7
Claude-Session: https://claude.ai/code/session_01GdSpCLd7NBZuuePJmcDzU7

Co-authored-by: ChrisRackauckas-Claude <accounts@chrisrackauckas.com>
Co-authored-by: Claude Opus 5 (1M context) <noreply@anthropic.com>
`FullSpecialize` was turning every parameter buffer into a `StaticArray`,
including ones whose elements are heap objects. That buys nothing — an
`MVector{1, Vector{Float64}}` still holds a pointer, so it is no more
GPU-resident than a `Vector{Vector{Float64}}` — and it breaks two ways:
`MArray` refuses `setindex!` on a non-isbits eltype (the FMI array-valued
discrete write), and the `SizedVector` workaround for that makes `remake`
recurse forever, because `similar_type` is a no-op on a fully parameterized
`SizedVector{N, T, Vector{T}}`.

Gate on `isbitstype` in both places that pick a buffer type: the generated
map's `_static_initialization_buffer`, and the `get_p_constructor` closure the
reconstructor routes through. Both sides have to use the same rule or the map
and `MTKParametersReconstructor` disagree on the buffer type. `get_p_constructor`
only started seeing non-isbits buffers because this PR is what makes `pType`
static in the first place.

Tunable and initials buffers are isbits and stay `MVector`, which is the
allocation and GPU win this PR is for.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Agent-Harness: Claude Code 2.0.14
Agent-Model: claude-opus-5[1m]
Agent-Session: https://claude.ai/code/session_01GdSpCLd7NBZuuePJmcDzU7
Claude-Session: https://claude.ai/code/session_01GdSpCLd7NBZuuePJmcDzU7
end
expr = build_explicit_observed_function(
initsys, Tuple(flat_syms);
expression = Val(true), compiler_options, kwargs...

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Relatively minor, but I think we should be able to just pass an options struct here and avoid kwargs....

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done in 16ef998c28. Both map builders now take a GeneratedFunctionOptions positionally and call the struct method of build_explicit_observed_function directly, so there is no kwargs... on either of them:

function _construct_fullspecialize_initializeprobpmap(
        sys::AbstractSystem, initsys::AbstractSystem, gen_opts::GeneratedFunctionOptions;
        p_constructor
    )
    ...
    expr = build_explicit_observed_function(initsys, Tuple(flat_syms), gen_opts)

The options come from the problem's own opts.fn_opts.codegen rather than being rebuilt:

function _fullspecialize_map_options(opts::SciMLProblemOptions)
    codegen = opts.fn_opts.codegen
    return GeneratedFunctionOptions(;
        expression = Val{true}, codegen.eval_expression, codegen.eval_module,
        compiler_options = opts.init_compiler_options,
        codegen_function_options = codegen.codegen
    )
end

This is a small behaviour fix as well as a cleanup, in two ways worth flagging:

  • The old kwargs... path went through the keyword compat shim, which funnelled whatever was left in the bag into a fresh Symbolics.CodegenFunctionOptions(; checkbounds, kwargs...). So the generated maps were not seeing the checkbounds/cse/... the rest of the problem was built with — only whatever happened to survive as loose keywords. Reusing codegen.codegen fixes that.
  • eval_expression is now honoured. Previously the builders took eval_module but never received eval_expression, so eval_or_rgf always defaulted to false and the maps were unconditionally RuntimeGeneratedFunctions even with eval_expression = true — inconsistent with the default map path, which forwards it. Say the word if you would rather keep the maps always-RGF and I will pin it back.

init_compiler_options is still what the maps compile under; it is just read off opts in the helper now, so it is dropped from the destructuring in maybe_build_initialization_problem where it had become unused.

Local on 16ef998c28, Julia 1.12.7:

GROUP=Initialization   InitializationSystem Test | 706 pass | 12 broken | 718 total | 20m43.6s
                       Guess Propagation         |  11 pass |  11 total
                       Initial Values Test       |  65 pass |  65 total     (0 errors)
GROUP=Extensions       HomotopyContinuation 76/76, LabelledArrays 18/18, BifurcationKit 7/7,
                       Despecialized MTKParameters AD 2/2, Initialization maps AD 7/7
GROUP=FMI              FMI | 60 pass | 60 total | 15m35.4s

All three exited 0. Runic, typos over the diff, and git diff --check are clean.

The two `FullSpecialize` map builders forwarded `kwargs...` into
`build_explicit_observed_function`, which routed through the keyword compat
shim and funnelled whatever was left in the bag into a fresh
`CodegenFunctionOptions`. Take a `GeneratedFunctionOptions` instead and call
the struct method directly.

`_fullspecialize_map_options` derives it from the problem's own
`opts.fn_opts.codegen`, so the maps now inherit the same `checkbounds`/`cse`
settings the rest of the problem was built with rather than whatever happened
to survive in `kwargs`, while still emitting an `Expr` and compiling under
`init_compiler_options`. `eval_expression` is now honoured too; previously the
generated maps ignored it and were always `RuntimeGeneratedFunction`s, unlike
the default map path.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Agent-Harness: Claude Code 2.0.14
Agent-Model: claude-opus-5[1m]
Agent-Session: https://claude.ai/code/session_01GdSpCLd7NBZuuePJmcDzU7
Claude-Session: https://claude.ai/code/session_01GdSpCLd7NBZuuePJmcDzU7
@ChrisRackauckas
ChrisRackauckas marked this pull request as ready for review September 13, 2026 10:18
@ChrisRackauckas
ChrisRackauckas merged commit 4764cc4 into SciML:master Sep 13, 2026
41 of 45 checks passed
ChrisRackauckas added a commit that referenced this pull request Sep 13, 2026
- Generate FullSpecialize initialization maps with static buffers (#5045)



Agent-Harness: Claude Code
Agent-Model: claude-opus-5[1m]
Claude-Session: https://claude.ai/code/session_014FEzNTLFutCmTEAZ3zBg5R

Co-authored-by: Claude Opus 5 (1M context) <noreply@anthropic.com>
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.

3 participants