Skip to content

Default MTK problems to AutoDespecialize - #4919

Draft
ChrisRackauckas-Claude wants to merge 12 commits into
SciML:masterfrom
ChrisRackauckas-Claude:agent/opaque-mtk-parameters
Draft

Default MTK problems to AutoDespecialize#4919
ChrisRackauckas-Claude wants to merge 12 commits into
SciML:masterfrom
ChrisRackauckas-Claude:agent/opaque-mtk-parameters

Conversation

@ChrisRackauckas-Claude

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

Copy link
Copy Markdown
Member

Ignore this draft until it has been reviewed by @ChrisRackauckas.

What changed and why

ModelingToolkit-generated problems now default to the new dynamic parameter policy, AutoDespecialize. AutoSpecialize is unchanged. DiffEqBase performs the generic solve-time promotion to DespecializedParameters; MTK constructs ordinary concrete MTKParameters and only supplies the unwrapping needed at generated-function, symbolic-indexing, reconstruction, initialization, parameter time-series, jump, linearization, and cache boundaries.

The distinction between the policies is:

  • AutoSpecialize: existing function-specialization behavior, unchanged by this work.
  • AutoDespecialize: dynamic p; one solver-facing outer type can hold changing concrete parameter layouts, with a function barrier restoring the concrete value immediately before the SciMLFunction call.
  • AutoRespecialize: the constrained, non-dynamic opaque-closure policy formerly named AutoDePSpecialize; the old name remains a deprecated alias.
  • FullSpecialize: fully concrete opt-out.

SciMLBase owns the public marker, DespecializedParameters, generic unwrapping, and barriers for the primary call of every built-in SciMLFunction family. DiffEqBase owns generic solve-time promotion and its precompile entry point; OrdinaryDiffEq packages precompile representative solves. MTK contains no solver wrapping or solve precompile implementation.

The prerequisite implementations are merged and registered:

Failing before / passing after

The same focused test was run against the unfixed and fixed implementations. With registered SciMLBase 3.46 but without the generic DiffEqBase implementation, 3 of 15 assertions fail because the second concrete parameter layout remains solver-facing and changes problem/function types:

Test Summary:                          | Pass  Fail  Total
AutoDespecialize generic call barrier |   12     3     15
ERROR: Some tests did not pass: 12 passed, 3 failed, 0 errored, 0 broken.

With the generic DiffEqBase implementation, the identical test passes:

Test Summary:                          | Pass  Total
AutoDespecialize generic call barrier |   15     15

The complete MTK discriminator passes 53/53. It covers default and opt-out policies, stable solver-facing problem types across different concrete MTK parameter-buffer layouts, generated RHS/Jacobian calls, symbolic getp/setp, SciMLStructures canonicalize/replace/replace!, reconstruction, initialization, remake, solving, observed values, jumps, ForwardDiff, Zygote, and SciMLSensitivity.

The SCC initialization path exposed a separate constructor-parameter mismatch on the
first registered-stack run. Before the fix, constructing the existing SCC test problem
errored with MethodError: no method matching SCCNonlinearProblem{true, AutoDespecialize}(...). The identical construction against registered SciMLBase 3.46.0
and DiffEqBase 7.14.0 passes after the fix and prints
REGISTERED_SCC_AUTODESPECIALIZE_PASS.

The first full root InterfaceII run then exposed the custom semilinear constructor
marker as the only remaining unsupported specialization query. The exact added assertion
failed before with MethodError: no method matching specialization(::Type{SemilinearODEFunction{true, AutoDespecialize}}) and passes after
adding the local interface method, printing
REGISTERED_SEMILINEAR_SPECIALIZATION_PASS.

The released-stack InterfaceI run exposed one non-split initialization path that rebuilt
p as a raw vector after DiffEqBase had promoted it. Before the fix, the symbolic-event
problem had a DespecializedParameters function wrapper but a Vector{Float64} concrete
parameter value and failed with No matching function wrapper was found!. The same
problem after the fix solves successfully, keeps DespecializedParameters in the concrete
problem and solution, and locates the event at t = 1 to floating-point precision.

The Julia 1.10 initialization group also showed that two existing @test_broken @inferred
checks now succeed. On the pushed head they therefore error with Expression evaluated to non-Boolean (660 pass / 2 error / 12 broken). Testing the inference directly instead gives
660 pass / 12 broken, followed by 65/65 initial-value tests.

Verification completed

# focused MTK AutoDespecialize test
53 passed, 53 total; exit 0

# ModelingToolkitBase SymbolicIndexingInterface group
Symbolic indexing: 64/64
SciML Problem Input: 1824/1824
MTKParameters: 118/118
Despecialized parameters: 53/53
package tests passed; exit 0

# registered release InterfaceII group
1040 passed, 7 pre-existing broken; exit 0

# registered release Initialization group, Julia 1.10
InitializationSystem: 660 passed, 12 pre-existing broken
Initial Values: 65/65
package tests passed; exit 0

# registered release discriminator
SciMLBase 3.46.0; DiffEqBase 7.14.0
REGISTERED_SCC_AUTODESPECIALIZE_PASS; exit 0
REGISTERED_SEMILINEAR_SPECIALIZATION_PASS; exit 0

# Runic on changed Julia files
# typos over the diff
# git diff --check
all exit 0

Root QA was also run on the exact stack. It reported 19 pass / 4 fail / 2 error:
the current-master JET, ExplicitImports, public-docstring, and reexport failures tracked at
#4670, plus one Aqua persistent-task
deadline while the cold wrapper cache was still compiling. Aqua's identical persistent-task
probe passed in a warm standalone released environment. Piracy and all public-owner checks
pass. No test was skipped or silenced.

The registered-release InterfaceI group reached 1500 pass / 5 pre-existing broken with
all 366 symbolic-event assertions passing. Its sole failure was the independently known
load-sensitive full_equations subprocess exceeding its fixed 180-second deadline during
concurrent compilation. The identical subprocess passed in isolation in 10.74 seconds.

The exact full documentation command completed all examples. It contains zero SCC
AutoDespecialize constructor errors and zero unresolved AutoDespecialize or
DespecializedParameters references. It exits 1 on the independently reproduced
current-master docs_block, missing_docs, cross_references, and linkcheck failures;
the clean-master comparison is
https://github.com/SciML/ModelingToolkit.jl/actions/runs/31478304896/job/93752799681.

Performance

A fresh-process, one-thread OrdinaryDiffEq benchmark used a second concrete parameter layout and a fixed 1,000-step solve:

Policy Same solver-facing type Second-layout compile Warm median
AutoDespecialize yes 0.142 s 0.555 ms
AutoSpecialize no 2.483 s 0.296 ms
FullSpecialize no 2.273 s 0.184 ms

AutoDespecialize reduced second-layout compile time by 94.3% versus AutoSpecialize. The deliberately tiny RHS exposes the fixed barrier overhead: warm execution was 1.88× AutoSpecialize, an additional 0.259 ms per 1,000-step solve.

Generated MTK systems show the intended larger-workload behavior:

States Policy Second-layout compile Warm median
64 AutoDespecialize 0.572 s 1.203 ms
64 AutoSpecialize 0.550 s 0.827 ms
64 FullSpecialize 1.457 s 0.890 ms
512 AutoDespecialize 0.575 s 3.880 ms
512 AutoSpecialize 0.580 s 5.273 ms
512 FullSpecialize 1.400 s 3.331 ms

For generated MTK functions, compile time is substantially lower than FullSpecialize and approximately tied with AutoSpecialize, because generated-function compilation dominates this benchmark. Warm performance is workload-dependent; at 512 states the wrapper costs about 0.55 ms versus FullSpecialize, while the measured AutoSpecialize median was noisier and slower. These results support a compile-time default, not a claim that the dynamic path is runtime-free.

Review points

Prior art

Copy link
Copy Markdown
Member Author

Follow-up commit dc26a82 addresses the Julia 1.10 allocation-sensitive test exposed by the first CI run.

  • Prior default opaque AutoSpecialize result: the unchanged allocation assertion evaluated as 16 == 0.
  • Focused Julia 1.10 reproduction: 16 allocated bytes.
  • Changed only that performance-specific problem construction to ODEProblem{true, SciMLBase.FullSpecialize}.
  • Focused result after: 0 allocated bytes. The == 0 assertion remains unchanged.
  • Full local Julia 1.10 GROUP=InterfaceI: 1,487 passed, 3 pre-existing broken, 1,490 total in 27m15.4s; exit 0.
  • Runic, typos, and git diff --check: exit 0.

Prior failing CI job: https://github.com/SciML/ModelingToolkit.jl/actions/runs/31363657429/job/93377585293

Copy link
Copy Markdown
Member Author

Architecture update: this draft will be refactored before review. The generic parameter container and function barrier are moving to SciMLBase as a reusable public API; ModelingToolkit will retain only the AutoSpecialize selection and MTK-specific interface hooks. The current MTK-owned OpaqueMTKParameters implementation and its CI/benchmark evidence should be treated as an implementation prototype until the prerequisite SciMLBase PR is validated and stacked here.

@ChrisRackauckas-Claude

Copy link
Copy Markdown
Member Author

The SciMLBase naming prerequisite is now visible as draft SciML/SciMLBase.jl#1510. The next stacked SciMLBase change will add the general DespecializedParameters dynamic barrier; this ModelingToolkit prototype will then be refactored to use that public API, removing most of its local OpaqueMTKParameters implementation. Local QA/docs for the naming prerequisite are still running, as recorded in that draft.

@ChrisRackauckas-Claude

Copy link
Copy Markdown
Member Author

Pushed the SciMLBase-backed refactor in commit 2ed51bb1db. The MTK-local OpaqueMTKParameters implementation is removed; default AutoSpecialize parameter wrapping now happens once in central MTK problem construction and uses SciMLBase.DespecializedParameters. This covers specialization-bearing MTK SciMLFunction families rather than only ODE problems; the focused final test passes 32/32, including ODE, generated expressions, nonlinear problems, symbolic indexing, SciMLStructures, remake, solve, jumps, and FullSpecialize. Runic, typos, and git diff --check pass. Full group/docs/AD/benchmark reruns on the refactor are still in progress. Dependencies: SciML/SciMLBase.jl#1510 and SciML/SciMLBase.jl#1511.

@ChrisRackauckas-Claude ChrisRackauckas-Claude changed the title Reduce MTK AutoSpecialize parameter recompilation Despecialize MTK parameters through SciMLBase Aug 10, 2026
@ChrisRackauckas-Claude

Copy link
Copy Markdown
Member Author

Additional local validation:

# clean upstream/master, direct MTKBase group
GROUP=SymbolicIndexingInterface julia +release --project=lib/ModelingToolkitBase -e 'using Pkg; Pkg.test()'\nSymbolicIndexingInterface test | 64/64\nSciML Problem Input Test | 1824/1824\nMTKParameters Test | 118/118\nTesting ModelingToolkitBase tests passed\n```\n\n```text\n# feature stack, focused ForwardDiff + Zygote + SciMLSensitivity test\ndespecialized parameter sensitivities | 2/2 | 5m37.1s\nfocused despecialized AD: pass\n```\n\nThe AD test exercises explicit unwrapping/reconstruction through the public SciMLBase interface. The full feature group is still running.

@ChrisRackauckas-Claude

Copy link
Copy Markdown
Member Author

The first CI failures are the expected stacked-dependency boundary, not test failures: each inspected job exits during resolution because this PR requires SciMLBase = "3.46", which is the unreleased dependency in SciML/SciMLBase.jl#1511. Example run: https://github.com/SciML/ModelingToolkit.jl/actions/runs/31378418992.

Local validations explicitly develop that SciMLBase branch and therefore exercise the intended stack. I am leaving the compat floor intact rather than weakening it or adding a temporary non-release source to the shipping project.

@AayushSabharwal AayushSabharwal left a comment

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.

Apart from these, why does this need to live in MTK? Why can't SciMLBase/DiffEqBase/whatever automatically wrap prob.p in DespecializedParams the same way it adds the FWWs?

::AbstractSystem, ps::Union{MTKParameters, SciMLBase.DespecializedParameters},
args::Pair{<:Any, <:NestedGetIndex}...
)
unwrapped_ps = _unwrap_mtk_parameters(ps)

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.

This should recurse and dispatch on unwrapped_ps, otherwise the loop body is slow.

ic = get_index_cache(sys) # this exists because the parameters are `MTKParameters`
isempty(ps.discrete) && return nothing
num_discretes = only(blocksize(ps.discrete[1]))
unwrapped_ps = _unwrap_mtk_parameters(ps)

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.

Similarly, recurse on unwrapped_ps

src_ps = parameter_values(src)
dst_ps = parameter_values(dst)
oldcache = dst_ps.caches
unwrapped_dst_ps = _unwrap_mtk_parameters(dst_ps)

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.

Recurse on unwrapped_dst_ps

@ChrisRackauckas

Copy link
Copy Markdown
Member

It is pretty much just in SciMLBase now. It's defined there, and all of the solver interactions are defined there. What's defined here are just things that assume you have an MTKParameters: now it's possible for it to be a DespecializedParameters, so MTK needs to unwrap them in order to do what it's doing, like for observed functions and such. So the only things that are here are (1) wrapping it for the solver and (2) unwrapping it where it previously assumed it would have an MTKParameters.

So if you go through it's like https://github.com/SciML/ModelingToolkit.jl/pull/4919/changes#diff-03e65df124406336ab53d3df653e5f77efb7e6dc09c13462e4d4e66f8dd61cffR194 unwrapping, https://github.com/SciML/ModelingToolkit.jl/pull/4919/changes#diff-af7ea5152025e34e41b4af797ced76c1c311b16bc1714c659050ebb3f73c6a0fR720 unwrapping, etc.

Now yeah, one last thing is, is it possible for this to work with a DiffEqBase autowrap? Maybe. It would still need some unwrapping in routines here so that it's aware of DespecializedParameters. But that will naturally have to compile at least a little bit more (because it would despecialize later). I can figure out if that can work out though.

@ChrisRackauckas-Claude ChrisRackauckas-Claude changed the title Despecialize MTK parameters through SciMLBase Default MTK problems to AutoDespecialize Aug 11, 2026
@ChrisRackauckas-Claude

Copy link
Copy Markdown
Member Author

Pushed the corrected architecture as a normal fast-forward in commit 2849525c78.

  • SciMLBase v3.46.0 is now registered and owns AutoDespecialize, DespecializedParameters, generic unwrapping, and the dynamic call barrier for every built-in SciMLFunction family: Add AutoDespecialize parameter policy SciMLBase.jl#1514.
  • DiffEqBase owns solve-time parameter promotion and the solver precompile entry point. The OrdinaryDiffEq implementation is locally validated and will be linked here as a separate draft; it is not implemented in MTK.
  • MTK defaults its specialization-bearing generated functions/problems to AutoDespecialize, keeps explicit AutoSpecialize unchanged, constructs ordinary concrete MTKParameters, and adds only MTK-specific unwrapping/reconstruction/indexing/initialization hooks.
  • AutoRespecialize remains the constrained non-dynamic policy; AutoDePSpecialize is its compatibility alias.

Post-rebase focused validation is AutoDespecialize parameters | 50/50; Runic, typos, and git diff --check exit 0. Immediately before the final upstream QA/version sync, the complete MTKBase group passed 64/64 + 1820/1820 + 117/117 + 50/50. The identical group on the final pushed commit, root QA, docs, and serial benchmarks are running now; I will report their actual results rather than infer them.

@ChrisRackauckas-Claude

Copy link
Copy Markdown
Member Author

The first completed CI failures are confirmed registry-propagation failures, not test assertions. Representative job https://github.com/SciML/ModelingToolkit.jl/actions/runs/31481133357/job/93746089848 stops in julia-buildpkg before downstream checkout/tests:

SciMLBase log:
├─possible versions are: 1.0.0 - 3.45.0 or uninstalled
└─restricted to versions 3.46.0 - 3 by ModelingToolkitBase — no versions left

SciMLBase v3.46.0 is registered in General at JuliaRegistries/General#164118, but the fresh runner's package-server registry snapshot still ends at 3.45.0. Catalyst, NeuralPDE, and SciMLBase downstream jobs fail at this same build step. I am not weakening the required 3.46 floor. DiffEqBase 7.14 is the next stacked release prerequisite once the OrdinaryDiffEq draft is reviewed and merged.

@ChrisRackauckas-Claude

Copy link
Copy Markdown
Member Author

The generic solver path requested in review is now published as SciML/OrdinaryDiffEq.jl#4216.

The ownership boundary on the current stack is:

  • SciMLBase owns AutoDespecialize, DespecializedParameters, unwrapped_f, and call barriers for every built-in SciMLFunction family.
  • DiffEqBase promotes parameters generically at solve time for ODE/DAE/DDE/SDE/SDDE problems and exposes the precompile entry point.
  • OrdinaryDiffEq precompiles representative generic solves.
  • This MTK PR only selects AutoDespecialize by default and unwraps parameters at MTK indexing/reconstruction/generated-function boundaries.

The identical generic barrier test fails 3/15 assertions before the DiffEqBase implementation and passes 15/15 after it. The complete MTK discriminator passes 50/50, and the full ModelingToolkitBase SymbolicIndexingInterface group passes 64 symbolic-indexing, 1,820 SciML Problem Input, 117 MTKParameters, and 50 despecialized-parameter assertions. I have updated the PR body with these results and the final compile/runtime benchmarks. I have not resolved the existing review threads so their disposition remains explicit for the reviewer.

@ChrisRackauckas-Claude

Copy link
Copy Markdown
Member Author

Package-server propagation update (2026-08-11):

A fresh local probe using the package server eager registry preference now resolves and installs SciMLBase v3.46.0 at registered tree d7a1a327339865b422b89690d77bb135db810778.

The same probe with the default registry preference, after Pkg.Registry.update(), still fails:

SciMLBase log:
├─possible versions are: 1.0.0 - 3.44.0 or uninstalled
└─restricted to versions 3.46.0 by an explicit requirement — no versions left

The default registry tarball contains versions only through 3.44, whereas the eager tarball contains 3.46. I am therefore not retrying the failed standard CI jobs yet; they would still resolve against the stale default snapshot. Registration: JuliaRegistries/General#164118

@ChrisRackauckas-Claude
ChrisRackauckas-Claude force-pushed the agent/opaque-mtk-parameters branch from 2849525 to 0596aa9 Compare August 11, 2026 17:48
@ChrisRackauckas-Claude

Copy link
Copy Markdown
Member Author

Rebase and registered-release update:

  • A literal rebase of the branch's historical merge topology attempted to replay the obsolete OpaqueMTKParameters prototype and conflicted. I rebuilt the branch from the final six first-parent implementation commits instead. git range-diff reports all six patches as identical, and they rebased onto current master without conflicts.
  • The branch now targets 3061ef5 and the clean head is ChrisRackauckas-Claude@0596aa9.
  • The local registered graph resolves SciMLBase 3.46.0 and DiffEqBase 7.14.0 from General, with no source overrides.
  • The SCC initialization discriminator failed before with SCCNonlinearProblem{true, AutoDespecialize} and passes after (REGISTERED_SCC_AUTODESPECIALIZE_PASS). The existing SCC group passed 36/36.
  • The first full InterfaceII run found one additional feature bug: the custom SemilinearODEFunction constructor marker lacked the specialization query. The exact added assertion failed before with a MethodError and passes after (REGISTERED_SEMILINEAR_SPECIALIZATION_PASS). That run reached 995 pass / 2 fail / 12 error / 7 broken, with every failure in that semilinear path.
  • Runic, typos over the diff, and git diff --check pass. The exact registered-release InterfaceII rerun and full docs build are still active; I am not claiming those results before they finish.

No review thread was resolved by this update.

@ChrisRackauckas-Claude

Copy link
Copy Markdown
Member Author

Full local docs result: the exact command completed all examples and exited 1 on the independently reproduced current-master categories [:docs_block, :missing_docs, :cross_references, :linkcheck].

Feature discrimination on the full log:

  • SCC AutoDespecialize constructor MethodErrors: 0
  • unresolved AutoDespecialize references: 0
  • unresolved DespecializedParameters references: 0
  • semilinear specialization MethodErrors: 0

The clean-master comparison is https://github.com/SciML/ModelingToolkit.jl/actions/runs/31478304896/job/93752799681. A separate clean-master audit is reducing those failures; no baseline docs change is included here.

@ChrisRackauckas-Claude

Copy link
Copy Markdown
Member Author

Pushed two follow-up commits after validating against the fully registered release graph (SciMLBase 3.46.0, DiffEqBase 7.14.0, OrdinaryDiffEq 7.6.0):

  • ce2902a422036a968cccceab7e3b20d12b165cb6 preserves DespecializedParameters when non-split MTK initialization reconstructs p. Before the fix, the symbolic-event path paired a DespecializedParameters function wrapper with a raw Vector{Float64} and errored with No matching function wrapper was found!. After the fix, the concrete problem and solution both retain DespecializedParameters, and the event is located at t = 1.
  • 8e5c44a310b98741acf0b9407a50adba25ec3ff6 directly tests two Julia 1.10 inference checks that now pass. Before this commit, the successful inferred values were incorrectly wrapped in @test_broken, producing two Expression evaluated to non-Boolean errors.

Local verification on the pushed tree:

SymbolicIndexingInterface: 64/64
SciML Problem Input: 1824/1824
MTKParameters: 118/118
Despecialized parameters: 53/53
InterfaceII: 1040 pass / 7 pre-existing broken
Initialization (Julia 1.10): 660 pass / 12 pre-existing broken, then 65/65
InterfaceI: 1500 pass / 1 load-sensitive failure / 5 pre-existing broken

The sole InterfaceI failure was the existing subprocess missing its fixed 180-second deadline under concurrent cold compilation. The identical subprocess passed in isolation in 10.74 seconds.

Root QA reported 19 pass / 4 fail / 2 error: the current-master JET, ExplicitImports, public-docstring, and reexport failures tracked in #4670, plus Aqua's persistent-task deadline while its cold wrapper cache was still compiling. Aqua's identical persistent-task probe passed in a warm standalone released environment. Piracy and all public-owner checks pass. Runic, typos over added lines, and git diff --check all exit 0. No test was skipped or silenced.

The corrected same-system Float64/Float32 parameter-layout discriminator also found that JumpProcesses rebuilds extended ODE/SDE/DDE/DAE functions as AutoSpecialize; this is an honest downstream compile-sharing limitation, not a parameter-correctness failure. It is tracked separately at #4944.

Commit links:


if p isa MTKParameters
buffer = p.initials
unwrapped_p = _unwrap_mtk_parameters(p)

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.

This should recurse on unwrapped_p

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?

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.

Wrong place, my bad

@ChrisRackauckas-Claude

Copy link
Copy Markdown
Member Author

Compile timings: default / Rodas5P / FBDF (before vs after MTK-shaped precompile)

Isolated first-solve compile, one fresh Julia 1.12.6 process per cell, JULIA_NUM_THREADS=1, save_everystep=false.

  • Before: 8e5c44a3 (this PR as-is). Root MTK already precompiles a 2-state ODE solve for default / Rodas5P / FBDF.
  • After: c4700d74 plus the MTKBase helpers moved to solver_precompile.jl so using ModelingToolkit actually loads (root MTK includes MTKBase's precompile.jl, and defining those helpers there was treated as a method extension).

Second-solve compile is 0 in every cell.

Isolated first solve compile

Solve Problem Before After Ratio
solve(prob) 2-state ODE 1.25 s 0.76 s 1.7×
Rodas5P() 2-state ODE 1.14 s 0.81 s 1.4×
FBDF() 2-state ODE 1.04 s 0.77 s 1.3×
solve(prob) index-1 mass-matrix DAE (ODEProblem) 18.50 s 1.49 s 12.4×
Rodas5P() index-1 mass-matrix DAE 9.73 s 1.67 s 5.8×
FBDF() index-1 mass-matrix DAE 5.88 s 1.63 s 3.6×
solve(prob) DAEProblem 7.45 s 5.33 s 1.4×

ODEProblem / DAEProblem construction is unchanged (~15 s ODE / ~19 s DAE) and is generated f + initialization, not the solver.

Breakout: wrap / first f / init / solve!

Same systems, still isolated. wrap is DiffEqBase.get_concrete_problem. solve! is the stepper after init.

wrap first f init solve!
after solve(ode) 0.67 0.05 0.50 0.00
after Rodas5P(ode) 0.63 0.05 0.18 0.05
after FBDF(ode) 0.67 0.05 0.26 0.03
after solve(mmdae) 1.00 0.05 0.91 0.00
after Rodas5P(mmdae) 1.01 0.05 0.63 0.05
after FBDF(mmdae) 1.00 0.08 0.64 0.03
after solve(dae) 0.74 0.09 2.47 2.68
after solve(dae) + NoInit 0.72 0.08 1.88 2.83
after DFBDF()(dae) 0.71 0.08 2.63 2.98
after DFBDF(AutoFiniteDiff())(dae) 0.70 0.08 2.34 2.78
before solve(ode) 0.63 0.05 1.06 0.00
before Rodas5P(ode) 0.63 0.05 1.00 0.05
before FBDF(ode) 0.66 0.05 0.93 0.02
before solve(mmdae) 1.22 0.05 4.74 11.45
before solve(dae) 0.80 0.11 3.11 4.50
before solve(dae) + NoInit 0.80 0.10 2.61 3.93

ODE / mass-matrix solve! is ~0. The leftover first-solve time there is this system's generated-function wrap plus OverrideInit, not the stepper.

Why DAEProblem does not compile down

It is the solver. After wrap, dae.f.f is still GeneratedFunctionWrapper{…, RGF hash, …}. DiffEqBase.promote_f only wraps ODEFunction:

wrap_path = f isa ODEFunction && isinplace(f) && !(f.f isa AbstractSciMLOperator) && 

DAEFunction misses that gate, so AutoDespecialize never installs FunctionWrappersWrappers + DespecializedParameters on the residual. DFBDF's integrator / cache / solve! therefore specialize on this system's RGF hash and recompile for every new DAE.

That is also why the mass-matrix DAE improved 12× and the residual-form DAEProblem did not: mmdae is an ODEProblem, so it takes the wrap path. NoInit only drops ~0.6 s of DAE init (OverrideInit); solve! stays ~2.7 s.

solve(dae) is DFBDF(autodiff=AutoFiniteDiff()). The BDF precompile used DFBDF() (AutoForwardDiff). Both were timed; both keep the unwrapped GFW, so the autodiff mismatch is secondary.

The next compile win for residual DAEs is wrapping DAEFunction in promote_f (signature (out, du, u, p, t)), not more MTK precompile statements.

Ignore until reviewed by @ChrisRackauckas.

@ChrisRackauckas-Claude

Copy link
Copy Markdown
Member Author

Why DAE init still compiles (~1.8–2.3 s)

Follow-up to the before/after + wrap/init/solve! breakout. The leftover init time is not missing MTK precompile statements. Two things are true:

  1. DiffEqBase is not wrapping the DAE residual.
  2. DFBDF init is not precompiled anywhere that transfers to a new f.

DiffEqBase does not wrap DAEFunction

get_concrete_problem(::DAEProblem) does call promote_f, but the wrap gate is f isa ODEFunction:

wrap_path = f isa ODEFunction && isinplace(f) && !(f.f isa AbstractSciMLOperator) && 
if !wrap_path
    return (f, p_out)   # residual unchanged
end

AutoDespecialize still wraps p (_promote_parametersDespecializedParameters). _despecialize_auxiliary_functions only touches f.g (SDE). The residual f.f is returned as-is.

Isolated get_concrete_problem on the after tree (wrap_f_same=true in every case):

Problem spec f after wrap p after wrap
Handwritten DAEProblem AutoSpecialize same closure, not FWW still Vector{Float64}
MTK DAEProblem AutoDespecialize same GeneratedFunctionWrapper{hash} DespecializedParameters

So the Newton cache built in init is still specialized on this system's residual type.

That cache is what init is compiling

DFBDF alg_cache calls build_nlsolver, which builds NLNewton + DAEResidualJacobianWrapper + FiniteDiff Jacobian prep, all parameterized by f. OrdinaryDiffEqBDF 2.4.1's own @compile_workload is FBDF on Lorenz plus a mass-matrix ODEProblem. There is no DFBDF and no DAEProblem in that workload.

Isolated init compile, DFBDF(autodiff=AutoFiniteDiff()), Julia 1.12.6, JULIA_NUM_THREADS=1:

Problem initializealg init compile solve! compile
Handwritten DAE default (no initialization_data) 1.85 s 3.82 s
Handwritten DAE NoInit 1.76 s 3.68 s
MTK DAE default (OverrideInit) 2.33 s 2.69 s
MTK DAE NoInit 1.74 s 2.72 s

Handwritten and MTK NoInit are the same ~1.8 s. That number is DFBDF cache construction, not ModelingToolkit. The extra ~0.6 s on default MTK init is OverrideInit of the NonlinearLeastSquaresProblem (another unwrapped GFW).

Our MTK solve(dae) / DFBDF() precompile only caches one GFW hash, so the next system's alg_cache still misses. More solve(dae) statements in ModelingToolkit will not move this.

What would move it

  • Wrap DAEFunction in DiffEqBase.promote_f with signature (out, du, u, p, t) (same FWW + DespecializedParameters path ODE already takes).
  • Add a DFBDF / DAEProblem precompile in OrdinaryDiffEqBDF so the Newton/Jacobian cache exists for the wrapped type.

Ignore until reviewed by @ChrisRackauckas.

@ChrisRackauckas-Claude

Copy link
Copy Markdown
Member Author

Full compile-time benches (ODE / mmdae / DAE, default / Rodas5P / FBDF)

Collecting the isolated numbers and the wrap / init / solve! breakout in one place, not just the DAE init follow-up. Julia 1.12.6, JULIA_NUM_THREADS=1, save_everystep=false. One fresh process per cell.

  • Before: 8e5c44a3 (this PR as-is). Root MTK already precompiles a 2-state ODE solve for default / Rodas5P / FBDF.
  • After: c4700d74 (MTK-shaped Default / Rodas5P / FBDF / DAE precompiles, plus the helper-file split so using ModelingToolkit loads).

Second-solve compile is 0 in every isolated cell.

Problem construction is unchanged and is not in the solve numbers: ODEProblem ~15 s compile, DAEProblem / mmdae ~19 s. That is generated f + initialization codegen.

Isolated first solve compile

Solve Problem Before After Ratio
solve(prob) 2-state ODE 1.25 s 0.76 s 1.7×
Rodas5P() 2-state ODE 1.14 s 0.81 s 1.4×
FBDF() 2-state ODE 1.04 s 0.77 s 1.3×
solve(prob) index-1 mass-matrix DAE (ODEProblem) 18.50 s 1.49 s 12.4×
Rodas5P() index-1 mass-matrix DAE 9.73 s 1.67 s 5.8×
FBDF() index-1 mass-matrix DAE 5.88 s 1.63 s 3.6×
solve(prob) DAEProblem 7.45 s 5.33 s 1.4×

ODE first-solves only drop ~0.3 s because before already precompiled that stepper. The win is the mass-matrix DAE path, which before did not precompile.

Breakout: wrap / first f / init / solve!

Same systems, still isolated. wrap = DiffEqBase.get_concrete_problem. solve! is the stepper after init.

wrap first f init solve!
after solve(ode) 0.67 0.05 0.50 0.00
after Rodas5P(ode) 0.63 0.05 0.18 0.05
after FBDF(ode) 0.67 0.05 0.26 0.03
after solve(mmdae) 1.00 0.05 0.91 0.00
after Rodas5P(mmdae) 1.01 0.05 0.63 0.05
after FBDF(mmdae) 1.00 0.08 0.64 0.03
after solve(dae) 0.74 0.09 2.47 2.68
after solve(dae) + NoInit 0.72 0.08 1.88 2.83
after DFBDF()(dae) 0.71 0.08 2.63 2.98
after DFBDF(AutoFiniteDiff())(dae) 0.70 0.08 2.34 2.78
before solve(ode) 0.63 0.05 1.06 0.00
before Rodas5P(ode) 0.63 0.05 1.00 0.05
before FBDF(ode) 0.66 0.05 0.93 0.02
before solve(mmdae) 1.22 0.05 4.74 11.45
before solve(dae) 0.80 0.11 3.11 4.50
before solve(dae) + NoInit 0.80 0.10 2.61 3.93

ODE / mass-matrix: f is wrapped to FunctionWrappersWrapper + DespecializedParameters, so solve! is ~0. Leftover first-solve time is this system's generated-function wrap plus OverrideInit, not the stepper.

DAEProblem: f is not wrapped. solve! is 2.7–4.5 s because DFBDF specializes on this system's GeneratedFunctionWrapper hash.

Same-process after (warmed; not isolated)

One process, after loading ModelingToolkitBase + the solver extensions. Later rows benefit from earlier solves. Useful as a “second system in the same session” check, not as a first-solve number.

first compile second compile
Tsit5 ODE 0.97 s 0
Rodas5P ODE 0.30 s 0
Rodas5P mmdae 1.14 s 0
FBDF ODE 0.067 s 0
FBDF mmdae 0.094 s 0
DFBDF() DAE 4.79 s 0
solve(ode) 0.059 s 0
solve(mmdae) 0.087 s 0
solve(dae) 2.79 s 0

Split on a new system in that same process (wrap already paid): Tsit5 init 0.16 s / solve! 0; Rodas5P init 0.09 s / solve! 0; FBDF init 0.04 s / solve! 0; DFBDF() init 0.39 s / solve! 0.

Why DAE init is still ~1.8–2.3 s

Two facts, both measured:

  1. DiffEqBase does not wrap the DAE residual. promote_f’s wrap gate is f isa ODEFunction. AutoDespecialize still wraps p (DespecializedParameters). f.f is returned unchanged (wrap_f_same=true).
  2. DFBDF init is not precompiled in a way that transfers. OrdinaryDiffEqBDF 2.4.1’s workload is FBDF on Lorenz plus a mass-matrix ODEProblem. No DFBDF, no DAEProblem. Our MTK solve(dae) / DFBDF() precompile only caches one GFW hash.

Isolated init of DFBDF(autodiff=AutoFiniteDiff()):

Problem initializealg init solve!
Handwritten DAE default (no initialization_data) 1.85 s 3.82 s
Handwritten DAE NoInit 1.76 s 3.68 s
MTK DAE default (OverrideInit) 2.33 s 2.69 s
MTK DAE NoInit 1.74 s 2.72 s

Handwritten and MTK NoInit are the same ~1.8 s. That is DFBDFCachebuild_nlsolver → NLNewton + DAEResidualJacobianWrapper + FiniteDiff jac prep, parameterized by this f. The extra ~0.6 s on default MTK init is OverrideInit of the NonlinearLeastSquaresProblem (another unwrapped GFW).

solve(dae) is DFBDF(autodiff=AutoFiniteDiff()). The BDF precompile used DFBDF() (AutoForwardDiff). Both were timed; both keep the unwrapped GFW, so the autodiff mismatch is secondary.

What would move the remaining numbers

  • ODE / mmdae leftover (~0.75–1.7 s first solve): this system's wrap + OverrideInit. Not the stepper.
  • DAEProblem: wrap DAEFunction in DiffEqBase.promote_f with signature (out, du, u, p, t), and add a DFBDF / DAEProblem precompile in OrdinaryDiffEqBDF. More MTK solve(dae) statements will not move init or solve!.

Ignore until reviewed by @ChrisRackauckas.

Copy link
Copy Markdown
Member Author

DAE coverage is now split into two focused OrdinaryDiffEq PRs:

I ran a two-system ModelingToolkit discriminator against this PR head with both changes developed locally. Before metadata widening it failed 16 pass / 4 fail / 20 total: the two MTK systems still produced different problem and DFBDF cache types under both ForwardDiff and FiniteDiff. With the final generic path it passed:

MTK DAE AutoDespecialize | 20 pass / 20 total

The test verifies that the source MTKParameters and generated residual types differ, while the solver-facing problem and DFBDF cache types are equal after init; both systems solve successfully and match the analytic endpoint under AutoForwardDiff() and AutoFiniteDiff().

The separate BDF precompile benchmark reduced fresh DFBDF initialization compilation from 2.426 s to 0.394 s (about 84%).

@ChrisRackauckas-Claude
ChrisRackauckas-Claude force-pushed the agent/opaque-mtk-parameters branch from 8e5c44a to f16f92b Compare August 15, 2026 12:22
@ChrisRackauckas-Claude

Copy link
Copy Markdown
Member Author

Review update: f16f92b now recursively processes unwrapped_p, then restores the original parameter wrapper without re-entering symbolic late binding.

Fail-before discriminator, using the identical new nested-wrapper test on the unfixed source:

Despecialized MTKParameters Test | 21 pass | 1 error
MethodError: no method matching promote_type_with_nothing(::Type{Float64}, ::_OpaqueParameterWrapper)

Passing-after on the rebased head:

SymbolicIndexingInterface test   |   64/64
SciML Problem Input Test         | 1824/1824
MTKParameters Test               |  118/118
Despecialized MTKParameters Test |    55/55
Testing ModelingToolkitBase tests passed

The full root initialization group also passes on the final code:

Initialization | 840 pass | 12 pre-existing broken | 852 total | 29m26.9s
Testing ModelingToolkit tests passed

Runic over every changed Julia file, typos over the full diff, and git diff --check all exit 0. The branch is rebased onto current master at f324831.

merge!(defaults, Dict(v => p[k] for (k, v) in pairs(paramobj)))
elseif p isa MTKParameters
elseif p isa Union{MTKParameters, SciMLBase.DespecializedParameters}
p = _unwrap_mtk_parameters(p)

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.

Should recurse on the unwrapped p

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
@ChrisRackauckas-Claude
ChrisRackauckas-Claude force-pushed the agent/opaque-mtk-parameters branch from f16f92b to 0678b4c Compare August 17, 2026 12:28
@ChrisRackauckas-Claude

Copy link
Copy Markdown
Member Author

Updated the branch for the latest review in commit 0678b4c.

The earlier recursion in initializesystem.jl was in the wrong place and has been removed. defaults_from_u0_p now unwraps SciMLBase.DespecializedParameters and recurses on the concrete wrapped parameter object before the MTKParameters/generic dispatch. This supports arbitrary wrapped parameter objects instead of assuming every wrapper contains an object with a tunable field.

Failing-before / passing-after discriminator:

  • Before the correction, the new modelingtoolkitize regression with DespecializedParameters((a = 10, b = 20)) exited 1 at common.jl:401: FieldError: type NamedTuple has no field tunable.
  • With the correction, the identical regression exits 0 and preserves both unknown and parameter initial conditions. A standalone check on the pushed tree printed despecialized_namedtuple_modelingtoolkitize=pass.

Local verification on Julia 1.12:

GROUP=InterfaceI julia +1.12 --startup-file=no --project=. -e 'using Pkg; Pkg.test()'
InterfaceI | 1599 passed, 5 broken, 1604 total | 42m47.8s
Testing ModelingToolkitBase tests passed

GROUP=QA julia +1.12 --startup-file=no --project=. -e 'using Pkg; Pkg.test()'
QA | 52 passed, 52 total | 4m30.0s
Testing ModelingToolkit tests passed

ModelingToolkitBase's concrete JET test passed 54/54. Its package-wide QA reproduced the current-master result (13 passed, 2 failed, 6 errored: existing JET/ExplicitImports/reexport debt), documented at #4670 (comment) and #4958 (comment).

Runic, added-line typos, and git diff --check passed. A strict clean Julia 1.12 docs build reached Documenter and failed only on the existing 34 unresolved cross-references; that exact current-master cluster is addressed separately by #4980. The latest review commit changes no documentation or public API, and no docs/test checks were silenced.

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