ModelingToolkit extension completely freezing Julia

[isaacsas]

using ModelingToolkit
using StructuralIdentifiability
const SI = StructuralIdentifiability

function testSI()
    t = ModelingToolkit.t_nounits
    D = ModelingToolkit.D_nounits
    @variables X(t)
    @parameters p d
    eqs = [
        D(X) ~ p - d*X
    ]
    @mtkbuild osys = ODESystem(eqs, t)
    measured_quantities = [X]    
    funcs_to_check = [osys.p]

    SI.assess_identifiability(osys; measured_quantities, funcs_to_check)
end

causes Julia to completely hang (as reported on Slack by @TorkelE).

Note that changing to funcs_to_check = [p] fixes the problem.

Environment:

(StructuralIdentifiability) pkg> st
Project StructuralIdentifiability v0.5.9
Status `~/.julia/dev/StructuralIdentifiability/Project.toml`
⌅ [c3fe647b] AbstractAlgebra v0.41.11
  [861a8166] Combinatorics v1.0.2
  [864edb3b] DataStructures v0.18.20
⌅ [0b43b601] Groebner v0.7.5
  [c8e1da08] IterTools v1.10.0
  [1914dd2f] MacroTools v0.5.13
⌅ [2edaba10] Nemo v0.45.7
⌅ [3e851597] ParamPunPam v0.4.1
  [aea7be01] PrecompileTools v1.2.1
  [27ebfcd6] Primes v0.5.6
  [a759f4b9] TimerOutputs v0.5.24
  [ade2ca70] Dates
  [37e2e46d] LinearAlgebra
  [56ddb016] Logging
  [9a3f8284] Random
Info Packages marked with ⌅ have new versions available but compatibility constraints restrict them from upgrading. To see why use `status --outdated`

on Julia 1.10.5.

[pogudingleb]

Many thanks for reporting ! I am looking into this

[pogudingleb]

I just tried with Julia 1.10.0 and was not able to reproduce the problem: both versions (osys.p and p) work fine. I will try to upgrade my Julia to see if I can get the hanging behaviour.

Environment:

Status `~/.julia/environments/v1.10/Project.toml`
⌅ [c3fe647b] AbstractAlgebra v0.41.11
⌅ [66b61cbe] AlgebraicSolving v0.4.16
  [4c88cf16] Aqua v0.8.7
  [6e4b80f9] BenchmarkTools v1.5.0
  [0c46a032] DifferentialEquations v7.14.0
  [7c1d4256] DynamicPolynomials v0.6.0
  [60bf3e95] GLPK v1.2.1
⌅ [0b43b601] Groebner v0.7.5
  [7073ff75] IJulia v1.25.0
  [c8e1da08] IterTools v1.10.0
  [682c06a0] JSON v0.21.4
  [98e50ef6] JuliaFormatter v1.0.60
  [b964fa9f] LaTeXStrings v1.3.1
  [16fef848] LiveServer v1.3.1
⌃ [291d046c] MixedSubdivisions v1.1.4
⌃ [961ee093] ModelingToolkit v9.39.1
⌅ [2edaba10] Nemo v0.45.7
  [bac558e1] OrderedCollections v1.6.3
  [f1435218] Oscar v1.1.1
⌅ [3e851597] ParamPunPam v0.4.1
  [91a5bcdd] Plots v1.40.8
  [67491407] Polyhedra v0.7.8
  [731186ca] RecursiveArrayTools v3.27.0
⌃ [295af30f] Revise v3.5.18
  [276daf66] SpecialFunctions v2.4.0
  [2913bbd2] StatsBase v0.34.3
  [220ca800] StructuralIdentifiability v0.5.9 `~/.julia/dev/StructuralIdentifiability`
⌃ [d1185830] SymbolicUtils v3.6.0
  [98d24dd4] TestSetExtensions v3.0.0
  [9a3f8284] Random
  [8dfed614] Test

[pogudingleb]

Two observations so far:

• there is no problem in Julia 1.10.0 but there is on 1.10.5
• if the osys.p version, if I interrupt by Ctrl+C, then the structural identifiability computation runs smoothly

[vyudu]

It seems like it is hanging due to some kind of type inference. Cutting the compiler shows that the compiler is hanging in specialize_method, which does some type inference. When I snoop the compiling process it seems like it is compiling reduce_ode_mod_p when it hangs (I called this with a ReactionSystem but essentially the same trace shows up if you call it with an ODESystem).

17-element Vector{Core.Compiler.Timings.Timing}:
 Core.Compiler.Timings.Timing(InferenceFrameInfo for Core.Compiler.Timings.ROOT()) with 0 children
 Core.Compiler.Timings.Timing(InferenceFrameInfo for StructuralIdentifiability.assess_identifiability(::ReactionSystem{Catalyst.NetworkProperties{Int64, SymbolicUtils.BasicSymbolic{Real}}})) with 1 children
 Core.Compiler.Timings.Timing(InferenceFrameInfo for CatalystStructuralIdentifiabilityExtension.var"#assess_identifiability#3"(::Vector{Any}, ::Vector{Any}, ::Vector{Any}, ::Bool, ::Bool, ::Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, ::typeof(assess_identifiability), ::ReactionSystem{Catalyst.NetworkProperties{Int64, SymbolicUtils.BasicSymbolic{Real}}})) with 7 children
 Core.Compiler.Timings.Timing(InferenceFrameInfo for Core.kwcall(::NamedTuple{(:measured_quantities, :funcs_to_check), <:Tuple{Vector{Equation}, Any}}, ::typeof(assess_identifiability), ::ReactionSystem)) with 1 children
 Core.Compiler.Timings.Timing(InferenceFrameInfo for CatalystStructuralIdentifiabilityExtension.var"#assess_identifiability#3"(::Vector{Equation}, ::Vector{Any}, ::Any, ::Bool, ::Bool, ::Base.Pairs{Symbol, V, Tuple{Vararg{Symbol, N}}, NamedTuple{names, T}} where {V, N, names, T<:Tuple{Vararg{Any, N}}}, ::typeof(assess_identifiability), ::ReactionSystem)) with 7 children
 Core.Compiler.Timings.Timing(InferenceFrameInfo for Core.kwcall(::NamedTuple, ::typeof(assess_identifiability), ::ReactionSystem)) with 0 children
 Core.Compiler.Timings.Timing(InferenceFrameInfo for CatalystStructuralIdentifiabilityExtension.var"#assess_identifiability#3"(::Any, ::Any, ::Any, ::Any, ::Any, ::Base.Pairs{Symbol, V, Tuple{Vararg{Symbol, N}}, NamedTuple{names, T}} where {V, N, names, T<:Tuple{Vararg{Any, N}}}, ::typeof(assess_identifiability), ::ReactionSystem)) with 7 children
 Core.Compiler.Timings.Timing(InferenceFrameInfo for Core.kwcall(::NamedTuple, assess_identifiability::typeof(assess_identifiability), ::ODE{P}) where P<:AbstractAlgebra.MPolyRingElem{Nemo.QQFieldElem}) with 0 children
 Core.Compiler.Timings.Timing(InferenceFrameInfo for StructuralIdentifiability.var"#assess_identifiability#694"(::Any, ::Vector{<:Union{AbstractAlgebra.Generic.FracFieldElem{P}, P}}, ::Float64, ::Any, assess_identifiability::typeof(assess_identifiability), ::ODE{P}) where P<:AbstractAlgebra.MPolyRingElem{Nemo.QQFieldElem}) with 1 children
 Core.Compiler.Timings.Timing(InferenceFrameInfo for Base.CoreLogging.with_logger(::StructuralIdentifiability.var"#695#697"{_A, Vector{var"#s1850"}, Float64, ODE{P}} where {_A, P<:AbstractAlgebra.MPolyRingElem{Nemo.QQFieldElem}, var"#s1850"<:Union{AbstractAlgebra.Generic.FracFieldElem{P}, P}, P<:AbstractAlgebra.MPolyRingElem{Nemo.QQFieldElem}}, ::Logging.ConsoleLogger)) with 0 children
 Core.Compiler.Timings.Timing(InferenceFrameInfo for Base.CoreLogging.with_logstate(::StructuralIdentifiability.var"#695#697"{_A, Vector{var"#s1850"}, Float64, ODE{P}} where {_A, P<:AbstractAlgebra.MPolyRingElem{Nemo.QQFieldElem}, var"#s1850"<:Union{AbstractAlgebra.Generic.FracFieldElem{P}, P}, P<:AbstractAlgebra.MPolyRingElem{Nemo.QQFieldElem}}, ::Base.CoreLogging.LogState)) with 2 children
 Core.Compiler.Timings.Timing(InferenceFrameInfo for (::StructuralIdentifiability.var"#695#697"{_A, Vector{var"#s1850"}, Float64, ODE{P}} where {_A, P<:AbstractAlgebra.MPolyRingElem{Nemo.QQFieldElem}, var"#s1850"<:Union{AbstractAlgebra.Generic.FracFieldElem{P}, P}, P<:AbstractAlgebra.MPolyRingElem{Nemo.QQFieldElem}})()) with 1 children
 Core.Compiler.Timings.Timing(InferenceFrameInfo for Core.kwcall(::NamedTuple{(:funcs_to_check, :prob_threshold), <:Tuple{Any, Float64}}, _assess_identifiability::typeof(StructuralIdentifiability._assess_identifiability), ::ODE{P} where P<:AbstractAlgebra.MPolyRingElem{Nemo.QQFieldElem})) with 0 children
 Core.Compiler.Timings.Timing(InferenceFrameInfo for StructuralIdentifiability.var"#_assess_identifiability#699"(::Any, ::Float64, _assess_identifiability::typeof(StructuralIdentifiability._assess_identifiability), ::ODE{P} where P<:AbstractAlgebra.MPolyRingElem{Nemo.QQFieldElem})) with 9 children
 Core.Compiler.Timings.Timing(InferenceFrameInfo for Core.kwcall(::NamedTuple{(:funcs_to_check, :prob_threshold, :type, :trbasis), <:Tuple{Any, Float64, Symbol, Vector{Nemo.QQMPolyRingElem}}}, _assess_local_identifiability::typeof(StructuralIdentifiability._assess_local_identifiability), ::ODE{P} where P<:AbstractAlgebra.MPolyRingElem{Nemo.QQFieldElem})) with 1 children
 Core.Compiler.Timings.Timing(InferenceFrameInfo for StructuralIdentifiability.var"#_assess_local_identifiability#198"(::Vector, ::Float64, ::Symbol, ::Vector{Nemo.QQMPolyRingElem}, ::Vector{Any}, _assess_local_identifiability::typeof(StructuralIdentifiability._assess_local_identifiability), ::ODE{P}) where P<:AbstractAlgebra.MPolyRingElem{Nemo.QQFieldElem}) with 42 children
 Core.Compiler.Timings.Timing(InferenceFrameInfo for StructuralIdentifiability.reduce_ode_mod_p(::ODE{<:AbstractAlgebra.MPolyRingElem{Nemo.QQFieldElem}}, ::Int64)) with 6 children

The timings before I cut it are here:

[Int(timing.time) for timing in timings]
17-element Vector{Int64}:
113583
305501
659876
248916
278582
243250
513582
759375
397000
220875
341249
357750
628750
1068958
1226791
2985504
680757720458

[pogudingleb]

Many thanks ! I will look closer into this function

[sumiya11]

Could be related : #362

[isaacsas]

We are now seeing it crashing / hanging CI when just loading the extension on 1.11:

https://github.com/SciML/Catalyst.jl/actions/runs/11239742228/job/31248152175?pr=1053#step:6:1197

[sumiya11]

CI crash on 1.11 seems to be a bug in Julia from JuliaLang/julia#56040.

The freeze from the current issue seems to be of different nature: it manifests in 1.10.5 and the offending commit from JuliaLang/julia#56040 is not present in 1.10.5.

@pogudingleb How do we debug the freeze in the current issue? The functions reduce_ode_mod_p and friends hinted by @vyudu seem as legal Julia code to me.

[pogudingleb]

Very elusive behaviour. Even if I take out the original code from the function:

using ModelingToolkit
using StructuralIdentifiability
const SI = StructuralIdentifiability

t = ModelingToolkit.t_nounits
D = ModelingToolkit.D_nounits
@variables X(t)
@parameters p d
eqs = [
    D(X) ~ p - d*X
]
@mtkbuild osys = ODESystem(eqs, t)
measured_quantities = [X]    
funcs_to_check = [osys.p]

SI.assess_identifiability(osys; measured_quantities, funcs_to_check)

It also works fine ! But being enclosed in a function, it hangs...

[isaacsas]

Hi, is there any update on this? Is it fixed on 1.10.6 and 1.11.2 (I know the latter isn't released yet, but given it should be soon that would close this issue I think)?

[isaacsas]

(I'm just trying to figure out if we can reenable the Catalyst SI extension for our next release or will have to keep it disabled.)

[pogudingleb]

Sorry, I was busy recently and the bug is tricky so the several hours I could spend were not enough to understand the problem. I will double check that it works on 1.10.6 and 1.11.2 and, if yes, close the issue. It seems that the SI code is completely legal and was not precompiling correctly for some extraterrestrial reason.

[pogudingleb]

Nope, it hangs as well on 1.10.6.