`MTKParameters` type mismatch in `init_dae!` when tunable buffer is a `ReinterpretArray` (DiffCache + setsym parameter-estimation pattern)

This report was drafted by Claude based on a sandboxed investigation, then reviewed by me before posting.

Summary

When a user follows a parameter-estimation pattern that puts a PreallocationTools.DiffCache buffer into the tunable slot of MTKParameters (e.g. via SciMLStructures.replace(Tunable(), p, get_tmp(cache, x))), the first type parameter of MTKParameters becomes Base.ReinterpretArray{...} instead of Vector{...}. The primal solve works, but ForwardDiff.gradient on the loss then crashes inside init_dae! with a MethodError: Cannot convert between two MTKParameters types that differ only in the array-type of the tunable field.

The mismatch is between:

• the integrator slot type — MTKParameters{ReinterpretArray{Dual,1,Float64,SubArray{...},false}, ...} (set when the integrator was built from prob.p whose tunable is the ReinterpretArray), and
• the p reconstructed by OverrideInit — MTKParameters{Vector{Dual}, ...} (a plain Vector because the MTK reconstructor’s PConstructorApplicator calls collect(::SubArray)).

MWE

using ModelingToolkit
using ModelingToolkit: t_nounits as t, subset_tunables
using ModelingToolkitStandardLibrary.Electrical
using ModelingToolkitStandardLibrary.Blocks: Sine
using OrdinaryDiffEqRosenbrock
using OrdinaryDiffEqNonlinearSolve
using PreallocationTools: DiffCache, get_tmp
using SymbolicIndexingInterface: setsym
using SciMLStructures: Tunable, replace as smlreplace
using SciMLBase: remake
using ForwardDiff

function create_model(; C₁ = 3e-5, C₂ = 1e-6)
    @named resistor1   = Resistor(R = 5.0)
    @named resistor2   = Resistor(R = 2.0)
    @named capacitor1  = Capacitor(C = C₁)
    @named capacitor2  = Capacitor(C = C₂)
    @named source       = Voltage()
    @named input_signal = Sine(frequency = 100.0)
    @named ground       = Ground()
    @named ampermeter   = CurrentSensor()
    eqs = [connect(input_signal.output, source.V)
        connect(source.p, capacitor1.n, capacitor2.n)
        connect(source.n, resistor1.p, resistor2.p, ground.g)
        connect(resistor1.n, capacitor1.p, ampermeter.n)
        connect(resistor2.n, capacitor2.p, ampermeter.p)]
    @named circuit_model = System(eqs, t,
        systems = [resistor1, resistor2, capacitor1, capacitor2,
                   source, input_signal, ground, ampermeter])
end

sys = mtkcompile(create_model())

guesses = Dict(
    sys.resistor1.p.i  => 0.0, sys.resistor2.p.i  => 0.0,
    sys.capacitor1.p.i => 0.0, sys.capacitor2.p.i => 0.0,
    sys.ampermeter.p.i => 0.0,
)

tspan, saveat = (0.0, 0.05), 0.0005
true_prob = ODEProblem(sys, Dict(), tspan; guesses)
data      = solve(true_prob, Rodas5P(); saveat)[sys.ampermeter.i]

sys_sub = subset_tunables(sys, [sys.capacitor2.C])
prob    = ODEProblem(sys_sub, Dict(), tspan; guesses = (
    sys_sub.resistor1.p.i  => 0.0, sys_sub.resistor2.p.i  => 0.0,
    sys_sub.capacitor1.p.i => 0.0, sys_sub.capacitor2.p.i => 0.0,
    sys_sub.ampermeter.p.i => 0.0) |> Dict)

cache   = DiffCache(copy(prob.p.tunable))
set_C2! = setsym(sys_sub, sys_sub.capacitor2.C)

function loss(x)
    tunables = get_tmp(cache, x)             # ReinterpretArray when x::Vector{Dual}
    copyto!(tunables, prob.p.tunable)
    new_p = smlreplace(Tunable(), prob.p, tunables)
    set_C2!(new_p, x[1])                     # in-place setsym update
    sol = solve(remake(prob; p = new_p), Rodas5P(); saveat)
    return sum(abs2, sol[sys_sub.ampermeter.i] .- data)
end

loss([1e-6])                                  # 0.0  — primal works
loss([2e-6])                                  # ~1.6e-6 — primal works
ForwardDiff.gradient(loss, [2e-6])            # 💥 see error below

Error

ERROR: MethodError: Cannot `convert` an object of type
  MTKParameters{Vector{ForwardDiff.Dual{Tag{typeof(loss),Float64},Float64,1}}, ...}
to an object of type
  MTKParameters{Base.ReinterpretArray{ForwardDiff.Dual{Tag{typeof(loss),Float64},Float64,1},
                                      1, Float64,
                                      SubArray{Float64,1,Vector{Float64},
                                               Tuple{UnitRange{Int64}}, true},
                                      false}, ...}
Stacktrace:
 [1] setproperty!(x::OrdinaryDiffEqCore.ODEIntegrator{…, MTKParameters{ReinterpretArray,…}, …}, …)
 [3] @ OrdinaryDiffEqCore  OrdinaryDiffEq.jl/lib/OrdinaryDiffEqCore/src/initialize_dae.jl:148
 [4] @ initialize_dae.jl:35 [inlined]   # _initialize_dae!  for OverrideInit
 [6] @ OrdinaryDiffEqCore  solve.jl:670
 ...

Diagnosis

1. SciMLStructures.replace(::Tunable, p, newvals) at lib/ModelingToolkitBase/src/systems/parameter_buffer.jl:373 is literally @set! p.tunable = newvals — it stores the ReinterpretArray as-is, so the returned MTKParameters’ first type parameter is ReinterpretArray{...}.
2. The integrator gets typed against that. So far so good — primal solve works.
3. During AD-typed solve, _initialize_dae!’s OverrideInit branch (lib/OrdinaryDiffEqCore/src/initialize_dae.jl:115-159) rebuilds parameters via initdata.initializeprobpmap, which goes through get_mtkparameters_reconstructor (lib/ModelingToolkitBase/src/systems/problem_utils.jl:735).
4. The reconstructor uses PConstructorApplicator (problem_utils.jl:693-719); for the SubArray branch (lines 705-707) with the default p_constructor = identity, it calls collect(x) and returns a plain Vector. That makes the reconstructed MTKParameters carry Vector{Dual} as its tunable type.
5. Back in initialize_dae.jl:148, integrator.p = p triggers convert(MTKParameters{ReinterpretArray,…}, MTKParameters{Vector{Dual},…}) — no such method exists.

Possible fixes

In rough order of how surgical they are:

1. MTK — make get_mtkparameters_reconstructor / PConstructorApplicator preserve the array-type of the original p.tunable instead of collapsing SubArrays with collect to Vector. Most direct fix at the root cause.
2. OrdinaryDiffEqCore — in _initialize_dae!’s OverrideInit branch, coerce the returned p.tunable (and friends) into the same array type as integrator.p.tunable before the assignment at initialize_dae.jl:148. More defensive — also protects against any other reconstructor losing the array type.
3. MTK — make SciMLStructures.replace(::Tunable, p, newvals) collect newvals to a Vector so a ReinterpretArray can never become the T parameter of MTKParameters. Removes the in-place AD buffer benefit users are after, so probably not desirable.

Environment

• ModelingToolkit master (commit 8b82464)
• OrdinaryDiffEq master (OrdinaryDiffEqCore / OrdinaryDiffEqRosenbrock / OrdinaryDiffEqNonlinearSolve dev'd from local checkout)
• SciMLBase master (dev'd locally)
• ModelingToolkitStandardLibrary 2.28.0, PreallocationTools 1.2.0, SciMLStructures 1.10.0, SymbolicIndexingInterface 0.3.47, ForwardDiff 1.3.3
• Julia 1.12.6

Reproduce the MWE first, then inspect parameter_buffer.jl:373 and the PConstructorApplicator and get_mtkparameters_reconstructor paths in problem_utils.jl:693-735. Compare them with the OverrideInit assignment in OrdinaryDiffEqCore's initialize_dae.jl:115-159. Done means the ForwardDiff.gradient call completes without the MTKParameters conversion error, with regression coverage for a ReinterpretArray tunable buffer.