Fix StructuralIdentifiability extension typing, DAE support, and test coverage

[devmotion]

Motivation

While debugging the StructuralIdentifiability extension on an older Catalyst version, I ran into a Vector{Vector{Pair{Any, Any}}} type error for models without conservation laws — caused by a typo in make_osys where Vector{Pair{Any, Any}}[] produces a nested vector instead of the intended Vector{Pair{Any, Any}}. Investigating further, I discovered that the latest release of StructuralIdentifiability (v0.5.20) breaks the extension tests on master with a TypeError: in typeassert, expected DataType, got Type{Union{Missing, Real}}.

The root cause of the TypeError turned out to be on the Catalyst side: cache_conservationlaw_eqs! declared the conservation law constant Γ with default = missing, which caused its binding to carry symtype = Missing. When SI.jl v0.5.20's __mtk_to_si substitutes parameter bindings into observable expressions, this produced Union{Missing, Real} terms that SymbolicUtils' promote_symtype rejected.

Additionally, the test file used raw let..end blocks, so the first failure silently hid all subsequent tests — including the nested-vector typo which was latent but real.

Changes

src/network_analysis.jl

• Remove = missing default from conservation law constant Γ. The guess metadata already handles initialization; the missing default was redundant and actively harmful.

src/reactionsystem_conversions.jl

• Γ initialization via System-level binding: Instead of = missing on the variable metadata (which pollutes the symtype), the Γ => missing binding and initialization equations are now placed at the System level following MTK's official parameter initialization pattern. This keeps the symtype clean for SI.jl while correctly telling MTK to solve for Γ during numerical initialization.
• addconstraints! now always computes the binding and initialization equations (removing the compute_cl_initeqs kwarg), returning them for callers to pass to the System constructor.
• hybrid_model gains an add_cl_bindings kwarg (default true) to control whether the Γ binding and initialization equations are included in the constructed System. Callers that need purely symbolic analysis (SI) can set this to false.
• ss_ode_model and sde_model (legacy path) updated to pass the binding and initialization equations through to the System constructor.

ext/.../structural_identifiability_extension.jl

• Type fix: Replace Vector{Pair{Any, Any}}[] (which produces Vector{Vector{Pair{Any, Any}}}) with type-annotated comprehensions Pair{Any, Any}[... for ...] that guarantee the correct element type even when the iterable is empty. The four isempty(...) && (...) fixup lines collapse away.
• DAE support: Add mtkcompile::Bool = false keyword argument to all four public entry points (make_si_ode, assess_identifiability, assess_local_identifiability, find_identifiable_functions) and to the internal make_osys. Follows the same conditional pattern as ODEProblem: calls mtkcompile when true, errors informatively when algebraic equations are present without it, otherwise calls complete.
• Skip conservation law bindings for SI: make_osys passes add_cl_bindings = false to ode_model since SI treats Γ as a free symbolic parameter and does not need the numerical initialization binding.

test/extensions/structural_identifiability.jl

• Convert all 8 let..end blocks to @testset with descriptive names so individual failures no longer cascade-hide subsequent tests.
• Update the coupled CRN/DAE test to use mtkcompile = true with explicit funcs_to_check (SI.jl currently cannot handle observed variables in funcs_to_check after mtkcompile; this is an upstream limitation).
• New tests: make_osys return type invariants (regression for the nested-vector fix), funcs_to_check kwarg, ignore_no_measured_warn kwarg.

Test plan

• Extension tests pass locally (11/11 testsets, 57 tests)
• NonlinearProblem with remove_conserved = true (conservation_laws.jl:131 path)
• hc_steady_states with conservation laws (homotopy_continuation.jl:28 path)
• CI TestExtensions.yml passes on Julia 1, lts, pre
• CI Tests passes

🤖 Generated with Claude Code

[isaacsas]

Looks like there are now a lot of initialization related test failures. Not sure if they are due to this PR or are new from it.

[ChrisRackauckas]

@AayushSabharwal

[AayushSabharwal]

That's weird. I don't remember seeing Catalyst tests fail on MTK. Will try to find what caused this.

[ChrisRackauckas]

They might've been blocked from updated due to some upper bounding.

[devmotion]

It's due to the removal of missing default values I assume, which however caused problems for SI. IMO it seems also a bit strange to specify any default values of missing - why specify any at all if they are unknown/missing?

[isaacsas]

Catalyst doesn't upper bound anything currently. We removed that prior to the last release a week ago.

The use of missing is part of what @AayushSabharwal told us we need to do to have initialization and remake work with conservation laws and their conserved constants.

[devmotion]

Is this special missing behaviour documented MTK API? The resulting symbolic union type of Real and Missing is not supported by the latest SI release, and IMO the use of missing seems generally a surprising and unintuitve API (if it is public API). Couldn't the desired behaviour here be achieved with eg a binding?

[TorkelE]

I'm currently traveling, but will try to have a proper look at this while doing so or when I get back next week. I remember there were lots of back and forth regarding the conservation laws to ensure we set them up right, so need to remember all that discussion.

Also, checking, do we have a MWE of what actually goes wrong the help me understand?

Thanks a lot for the help, sorry you caught be at an inconvenient timing!

[AayushSabharwal]

I just read the description of the first commit

Remove = missing default from conservation law constant Γ in
cache_conservationlaw_eqs!. The missing default caused Γ's binding
to carry symtype = Missing, which poisoned SI.jl's __mtk_to_si
substitution with Union{Missing, Real} → TypeError in
SymbolicUtils.promote_symtype. The guess metadata already handles
initialization.

Yeah this is wrong. Defaults, including missing, do not affect the symtype of the created variable.

Is this special missing behaviour documented MTK API?

Yes, in this section (including the table at the end) https://docs.sciml.ai/ModelingToolkit/stable/tutorials/initialization/#bindings_and_ics and in https://docs.sciml.ai/ModelingToolkit/stable/tutorials/initialization/#Initialization-of-parameters

[devmotion]

Yes, that was new to me but I adjusted it according to the documentation in the second commit (or at least according to how I understood the docs).

[devmotion]

Also, checking, do we have a MWE of what actually goes wrong the help me understand?

There's actually different issues and MWE involved here.

---

Originally, I hit the following issue:

I was working in an environment in which Symbolics was pinned to 6.55.0 and SymbolicUtils to 3.32.0 (the resolver picked StructuralIdentifiability 0.5.14, the latest compatible with those Symbolics pins). I ran into an ambiguity error with Symbolics.substitute, a simple MWE is:

using Pkg
Pkg.activate(; temp = true)
Pkg.add([
    PackageSpec(name = "Catalyst",        version = "15.0.11"),
    PackageSpec(name = "ModelingToolkit", version = "9.84.0"),
    PackageSpec(name = "Symbolics",       version = "6.55.0"),
    PackageSpec(name = "SymbolicUtils",   version = "3.32.0"),
])
Pkg.add("StructuralIdentifiability")   # resolver picks 0.5.14

using Catalyst, StructuralIdentifiability

rs = @reaction_network begin
    p, 0 --> X
    d, X --> 0
end
rs = complete(rs)
@unpack X, p = rs
@variables Y(Catalyst.get_iv(rs))

assess_identifiability(rs; measured_quantities = [Y ~ X*p])

This results in:

MethodError: substitute(::Equation, ::Vector{Vector{Pair{Any, Any}}}) is ambiguous.

Candidates:
  substitute(x::Equation, rules; kw...)        @ Symbolics .../equations.jl:134
  substitute(expr, s::Vector; kw...)            @ Symbolics .../num.jl:90

This made me look into the source code of the SI extension in Catalyst, which then made me realize that types are inconsistent and the nested vector was a (seemingly) unintentional bug.

Later on, I actually also noticed that with Symbolics 6.55, substitute(::Equation, ::Vector) is ambiguous regardless of the inner element type, so the PR's Catalyst-side typo fix (flat Vector{Pair{Any, Any}} instead of nested Vector{Vector{Pair{Any, Any}}}) doesn't by itself unblock that older Symbolics — the full resolution also needs Symbolics 7+ . What the typo fix does clean up is the root type-correctness of make_osys, which can be observed independent of the Symbolics version:

using Pkg
Pkg.activate(; temp = true)
Pkg.add(PackageSpec(name = "Catalyst", rev = "master"))
Pkg.add(PackageSpec(name = "StructuralIdentifiability", version = "0.5.19"))

using Catalyst, StructuralIdentifiability

# `const` lets inference see through the extension-module lookup.
const Ext = Base.get_extension(Catalyst, :CatalystStructuralIdentifiabilityExtension)
conseqs_of(rs; rc = true) = Ext.make_osys(rs; remove_conserved = rc)[2:3]

rs = @reaction_network begin
    p, 0 --> X
    d, X --> 0
end
rs = complete(rs)

using Test
@inferred conseqs_of(rs; rc = true)
@inferred conseqs_of(rs; rc = false)

Inference fails on master due to the typo:

return type Tuple{Vector{Vector{Pair{Any, Any}}}, Vector{Vector{Pair{Any, Any}}}} does not match
inferred return type Tuple{Union{Vector{Vector{Pair{Any, Any}}}, Vector{Pair{…BasicSymbolic{SymReal}, …}}},
                           Union{Vector{Vector{Pair{Any, Any}}}, Vector{Pair{…BasicSymbolic{SymReal}, …}}}}

On this PR, both inference checks pass.

---

When preparing the PR and running tests, I then noticed that SI 0.5.20 (released 3 days ago) causes additional new problems. SI 0.5.20's __mtk_to_si substitutes System-level parameter bindings into observable expressions before handing them to SymbolicUtils. cache_conservationlaw_eqs! on master declares the conservation-law constant Γ with = missing as a variable-metadata default, and the new substitution in SI then produces Union{Missing, Real} terms that SymbolicUtils.promote_symtype rejects. Under 0.5.19 the substitution never happened, so this was not a problem (and hence CI also passed a few days ago, prior to the release of SI 0.5.20).

using Pkg
Pkg.activate(; temp = true)
Pkg.add(PackageSpec(name = "Catalyst", rev = "master"))
Pkg.add(PackageSpec(name = "StructuralIdentifiability", version = "0.5.20"))

using Catalyst, StructuralIdentifiability

rs = @reaction_network begin
    k1, x1 --> x2
end
rs = complete(rs)

assess_identifiability(rs; measured_quantities = [:x1])

This errors with:

TypeError: in typeassert, expected DataType, got Type{Union{Missing, Real}}
Stacktrace:
  [1] promote_symtype(::typeof(+), T::DataType, S::DataType)
      @ SymbolicUtils .../methods.jl:208
  [2] _promote_symtype(f::Function, args::SmallVec{…})
      @ SymbolicUtils .../types.jl:1077
  [3] maketerm(::Type{…BasicSymbolic{SymReal}}, f::Any, args::…, metadata::Nothing)
      @ SymbolicUtils .../terminterface.jl:259
  [4] DefaultSubstituter(ex::…BasicSymbolic{SymReal})
      @ SymbolicUtils .../substitute.jl:198
  [5] substitute(expr, dict; fold=Val{false}, filterer)
      @ SymbolicUtils .../substitute.jl:263
  ⋮
 [12] __mtk_to_si(de::System, measured_quantities::Vector{…})
      @ ModelingToolkitSIExt .../ModelingToolkitSIExt.jl:283
 [21] assess_identifiability(::ReactionSystem; …)
      @ CatalystStructuralIdentifiabilityExtension .../structural_identifiability_extension.jl:118

On this PR, the example does not error anymore. I fixed it by moving Γ => missing from variable metadata to a System-level bindings entry (plus initialization_eqs), so in the non-SI use cases MTK still initialises Γ numerically while in the SI extension the symbolic path SI walks no longer sees missing default values since the bindings are not propagated to SI.

[devmotion]

Bump 🙂

[TorkelE]

I just got back last night, will look over the weekend!

[AayushSabharwal]

https://github.com/SciML/StructuralIdentifiability.jl/blob/master/ext/ModelingToolkitSIExt.jl#L274

-            if ModelingToolkitBase.isparameter(k)
+            if ModelingToolkitBase.isparameter(k) && Symbolics.value(v) !== missing

This fixes the issue. bindings(sys) can have entries for which the value is missing. SI needs to handle that case.

[pogudingleb]

https://github.com/SciML/StructuralIdentifiability.jl/blob/master/ext/ModelingToolkitSIExt.jl#L274

-            if ModelingToolkitBase.isparameter(k)
+            if ModelingToolkitBase.isparameter(k) && Symbolics.value(v) !== missing

This fixes the issue. bindings(sys) can have entries for which the value is missing. SI needs to handle that case.

@AayushSabharwal
Thanks a lot for looking into this! I will double check and make a PR to SI.jl.

[devmotion]

I have to recheck later (away from my computer right now) but I'm pretty sure the missings weren't the only issue I ran into, and definitely the bindings are not the only change in this PR. Additionally, while I agree that it would be good to improve the MTK support in SI if possible, IMO Catalyst shouldn't forward information to SI in the first place that is known to be irrelevant and has to be filtered out in SI.

[isaacsas]

The current Catalyst approach for handling conservation laws was developed and refined with an enormous amount of time and effort on @TorkelE, @AayushSabharwal, and my parts. Torkel and I have checked with Aayush and @ChrisRackauckas and the conclusion is we are handling them correctly for the current MTK versions, and this is an issue on the SI side.

So while we are happy to discuss improving the SI extension, we are not comfortable making changes to the core ReactionSystem conversion code that changes the approach for handling conservation laws / constants unless we are told that we are not using the correct approach by MTK devs (and again, in our discussions with Chris and Aayush the conclusion is that we are handling them correctly with how we forward them to MTK).

[devmotion]

As I said, how the bindings are constructed and whether and how they are forwarded to SI is only one part of the PR, and not even the issue I ran into initially.

IMO it's still conceptually wrong to forward things to SI we know have to be filtered out/ignored in SI - it just makes this path more brittle and complex, as the current bugs show, without any benefit. AFICT the PR also doesn't change the system in any way, it just arguably more cleanly separates the missing bindings from the definition of the onserved variables which makes it easy to include them or not (as in the SI case).

[pogudingleb]

I have fixed the missing's issue (the handling of this case was indeed suboptimal in SI), will make a new release right away. Let me know if there are any more issues which appear to be more on SI.jl side.

[isaacsas]

@pogudingleb thanks!

[AayushSabharwal]

how the bindings are constructed and whether and how they are forwarded to SI is only one part of the PR, and not even the issue I ran into initially.

Just to make sure we're on the same page, the initial issue was substitute not working on Symbolics@6, right? Is there a way to reproduce that on the latest versions now? And then the subsequent error was type-stability?

IMO it's still conceptually wrong to forward things to SI we know have to be filtered out/ignored in SI - it just makes this path more brittle and complex, as the current bugs show, without any benefit.

Respectfully, I believe the contrary. SI has an extension for MTK, not Catalyst. The MTK extension should be able to handle valid MTK systems, which include ones where the binding is missing. Please correct me if I'm misunderstanding something here, I'm not very familiar with SI. As I understand it, the main issue described here that may require changes to Catalyst is the type-instability one in case Catalyst is passing a type-unstable data structure to SI.