From 0e829d9a34bfb972c4b125a9501ae936a7b63a11 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Thu, 29 May 2025 18:17:55 +0530 Subject: [PATCH 01/42] fix: fix infinite recursion in `full_equations` --- .../symbolics_tearing.jl | 25 +++++----------- test/odesystem.jl | 29 +++++++++++++++++++ 2 files changed, 36 insertions(+), 18 deletions(-) diff --git a/src/structural_transformation/symbolics_tearing.jl b/src/structural_transformation/symbolics_tearing.jl index d41cbf12d2..3f46b594df 100644 --- a/src/structural_transformation/symbolics_tearing.jl +++ b/src/structural_transformation/symbolics_tearing.jl @@ -748,12 +748,14 @@ function update_simplified_system!( unknowns = [unknowns; extra_unknowns] @set! sys.unknowns = unknowns - obs, subeqs, deps = cse_and_array_hacks( - sys, obs, solved_eqs, unknowns, neweqs; cse = cse_hack, array = array_hack) + obs = cse_and_array_hacks( + sys, obs, unknowns, neweqs; cse = cse_hack, array = array_hack) + deps = Vector{Int}[i == 1 ? Int[] : collect(1:(i - 1)) + for i in 1:length(solved_eqs)] @set! sys.eqs = neweqs @set! sys.observed = obs - @set! sys.substitutions = Substitutions(subeqs, deps) + @set! sys.substitutions = Substitutions(solved_eqs, deps) # Only makes sense for time-dependent # TODO: generalize to SDE @@ -850,7 +852,7 @@ if all `p[i]` are present and the unscalarized form is used in any equation (obs not) we first count the number of times the scalarized form of each observed variable occurs in observed equations (and unknowns if it's split). """ -function cse_and_array_hacks(sys, obs, subeqs, unknowns, neweqs; cse = true, array = true) +function cse_and_array_hacks(sys, obs, unknowns, neweqs; cse = true, array = true) # HACK 1 # mapping of rhs to temporary CSE variable # `f(...) => tmpvar` in above example @@ -878,7 +880,6 @@ function cse_and_array_hacks(sys, obs, subeqs, unknowns, neweqs; cse = true, arr tempeq = tempvar ~ rhs_arr rhs_to_tempvar[rhs_arr] = tempvar push!(obs, tempeq) - push!(subeqs, tempeq) end # getindex_wrapper is used because `observed2graph` treats `x` and `x[i]` as different, @@ -887,10 +888,6 @@ function cse_and_array_hacks(sys, obs, subeqs, unknowns, neweqs; cse = true, arr neweq = lhs ~ getindex_wrapper( rhs_to_tempvar[rhs_arr], Tuple(arguments(rhs)[2:end])) obs[i] = neweq - subeqi = findfirst(isequal(eq), subeqs) - if subeqi !== nothing - subeqs[subeqi] = neweq - end end # end HACK 1 @@ -920,7 +917,6 @@ function cse_and_array_hacks(sys, obs, subeqs, unknowns, neweqs; cse = true, arr tempeq = tempvar ~ rhs_arr rhs_to_tempvar[rhs_arr] = tempvar push!(obs, tempeq) - push!(subeqs, tempeq) end # don't need getindex_wrapper, but do it anyway to know that this # hack took place @@ -960,15 +956,8 @@ function cse_and_array_hacks(sys, obs, subeqs, unknowns, neweqs; cse = true, arr push!(obs_arr_eqs, arrvar ~ rhs) end append!(obs, obs_arr_eqs) - append!(subeqs, obs_arr_eqs) - - # need to re-sort subeqs - subeqs = ModelingToolkit.topsort_equations(subeqs, [eq.lhs for eq in subeqs]) - - deps = Vector{Int}[i == 1 ? Int[] : collect(1:(i - 1)) - for i in 1:length(subeqs)] - return obs, subeqs, deps + return obs end function is_getindexed_array(rhs) diff --git a/test/odesystem.jl b/test/odesystem.jl index 3220555e62..83bb9b3621 100644 --- a/test/odesystem.jl +++ b/test/odesystem.jl @@ -1755,3 +1755,32 @@ end sol = solve(prob, Tsit5()) @test SciMLBase.successful_retcode(sol) end + +@testset "`full_equations` doesn't recurse infinitely" begin + code = """ + using ModelingToolkit + using ModelingToolkit: t_nounits as t, D_nounits as D + @variables x(t)[1:3]=[0,0,1] + @variables u1(t)=0 u2(t)=0 + y₁, y₂, y₃ = x + k₁, k₂, k₃ = 1,1,1 + eqs = [ + D(y₁) ~ -k₁*y₁ + k₃*y₂*y₃ + u1 + D(y₂) ~ k₁*y₁ - k₃*y₂*y₃ - k₂*y₂^2 + u2 + y₁ + y₂ + y₃ ~ 1 + ] + + @named sys = ODESystem(eqs, t) + + inputs = [u1, u2] + outputs = [y₁, y₂, y₃] + ss, = structural_simplify(sys, (inputs, [])) + full_equations(ss) + """ + + cmd = `$(Base.julia_cmd()) --project=$(@__DIR__) -e $code` + proc = run(cmd, stdin, stdout, stderr; wait = false) + sleep(120) + @test !process_running(proc) + kill(proc, Base.SIGKILL) +end From 338191663bcf19f16765a02bc39706086a341f0f Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Thu, 29 May 2025 18:49:46 +0530 Subject: [PATCH 02/42] fix: fix type promotion in `late_binding_update_u0_p` with non-dual types --- src/systems/nonlinear/initializesystem.jl | 44 ++++++++++++++++------- test/initial_values.jl | 10 ++++++ 2 files changed, 41 insertions(+), 13 deletions(-) diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl index d2a988dc07..8812a060ad 100644 --- a/src/systems/nonlinear/initializesystem.jl +++ b/src/systems/nonlinear/initializesystem.jl @@ -598,24 +598,42 @@ function SciMLBase.remake_initialization_data( return SciMLBase.remake_initialization_data(sys, odefn, newu0, t0, newp, newu0, newp) end -function promote_u0_p(u0, p::MTKParameters, t0) - u0 = DiffEqBase.promote_u0(u0, p.tunable, t0) - u0 = DiffEqBase.promote_u0(u0, p.initials, t0) +promote_type_with_nothing(::Type{T}, ::Nothing) where {T} = T +promote_type_with_nothing(::Type{T}, ::SizedVector{0}) where {T} = T +function promote_type_with_nothing(::Type{T}, ::AbstractArray{T2}) where {T, T2} + promote_type(T, T2) +end - if !isempty(p.tunable) - tunables = DiffEqBase.promote_u0(p.tunable, u0, t0) - p = SciMLStructures.replace(SciMLStructures.Tunable(), p, tunables) - end - if !isempty(p.initials) - initials = DiffEqBase.promote_u0(p.initials, u0, t0) - p = SciMLStructures.replace(SciMLStructures.Initials(), p, initials) +promote_with_nothing(::Type, ::Nothing) = nothing +promote_with_nothing(::Type, x::SizedVector{0}) = x +promote_with_nothing(::Type{T}, x::AbstractArray{T}) where {T} = x +function promote_with_nothing(::Type{T}, x::AbstractArray{T2}) where {T, T2} + if ArrayInterface.ismutable(x) + y = similar(x, T) + copyto!(y, x) + return y + else + yT = similar_type(x, T) + return yT(x) end - - return u0, p +end +function promote_with_nothing(::Type{T}, p::MTKParameters) where {T} + tunables = promote_with_nothing(T, p.tunable) + p = SciMLStructures.replace(SciMLStructures.Tunable(), p, tunables) + initials = promote_with_nothing(T, p.initials) + p = SciMLStructures.replace(SciMLStructures.Initials(), p, initials) + return p end function promote_u0_p(u0, p, t0) - return DiffEqBase.promote_u0(u0, p, t0), DiffEqBase.promote_u0(p, u0, t0) + T = Union{} + T = promote_type_with_nothing(T, u0) + T = promote_type_with_nothing(T, p.tunable) + T = promote_type_with_nothing(T, p.initials) + + u0 = promote_with_nothing(T, u0) + p = promote_with_nothing(T, p) + return u0, p end function SciMLBase.late_binding_update_u0_p( diff --git a/test/initial_values.jl b/test/initial_values.jl index 0ed8f7bffe..f66d103cc9 100644 --- a/test/initial_values.jl +++ b/test/initial_values.jl @@ -362,3 +362,13 @@ end @test state_values(initdata.initializeprob) isa SVector @test parameter_values(initdata.initializeprob) isa SVector end + +@testset "Type promotion of `p` works with non-dual types" begin + @variables x(t) y(t) + @mtkbuild sys = ODESystem([D(x) ~ x + y, x^3 + y^3 ~ 5], t; guesses = [y => 1.0]) + prob = ODEProblem(sys, [x => 1.0], (0.0, 1.0)) + prob2 = remake(prob; u0 = BigFloat.(prob.u0)) + @test prob2.p.initials isa Vector{BigFloat} + sol = solve(prob2) + @test SciMLBase.successful_retcode(sol) +end From 62b0b666e13a84b233a28aaeaefd85e1b9e4e06a Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Fri, 30 May 2025 13:14:00 +0530 Subject: [PATCH 03/42] build: bump patch version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 305f29c691..ad9a2e4f74 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ModelingToolkit" uuid = "961ee093-0014-501f-94e3-6117800e7a78" authors = ["Yingbo Ma ", "Chris Rackauckas and contributors"] -version = "9.80.1" +version = "9.80.2" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" From 1ac273fb06c978bc214a6db8913effcc371c16a4 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Wed, 4 Jun 2025 13:55:47 +0530 Subject: [PATCH 04/42] fix: handle non-split systems in initial value promotion --- src/systems/nonlinear/initializesystem.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl index 8812a060ad..4b0086840a 100644 --- a/src/systems/nonlinear/initializesystem.jl +++ b/src/systems/nonlinear/initializesystem.jl @@ -603,6 +603,9 @@ promote_type_with_nothing(::Type{T}, ::SizedVector{0}) where {T} = T function promote_type_with_nothing(::Type{T}, ::AbstractArray{T2}) where {T, T2} promote_type(T, T2) end +function promote_type_with_nothing(::Type{T}, p::MTKParameters) where {T} + promote_type_with_nothing(promote_type_with_nothing(T, p.tunable), p.initials) +end promote_with_nothing(::Type, ::Nothing) = nothing promote_with_nothing(::Type, x::SizedVector{0}) = x @@ -628,8 +631,7 @@ end function promote_u0_p(u0, p, t0) T = Union{} T = promote_type_with_nothing(T, u0) - T = promote_type_with_nothing(T, p.tunable) - T = promote_type_with_nothing(T, p.initials) + T = promote_type_with_nothing(T, p) u0 = promote_with_nothing(T, u0) p = promote_with_nothing(T, p) From c595416864bded32a4aab6ad3a73222c41e432dc Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Wed, 4 Jun 2025 18:48:33 +0530 Subject: [PATCH 05/42] fix: fix major compile time regression due to `concrete_getu` --- src/systems/problem_utils.jl | 50 +++++++++++++++++++++--------------- 1 file changed, 30 insertions(+), 20 deletions(-) diff --git a/src/systems/problem_utils.jl b/src/systems/problem_utils.jl index 5ff00b4845..f79e8f6286 100644 --- a/src/systems/problem_utils.jl +++ b/src/systems/problem_utils.jl @@ -617,30 +617,40 @@ struct ReconstructInitializeprob{GP, GU} ugetter::GU end +""" + $(TYPEDEF) + +A wrapper over an observed function which allows calling it on a problem-like object. +`TD` determines whether the getter function is `(u, p, t)` (if `true`) or `(u, p)` (if +`false`). +""" +struct ObservedWrapper{TD, F} + f::F +end + +ObservedWrapper{TD}(f::F) where {TD, F} = ObservedWrapper{TD, F}(f) + +function (ow::ObservedWrapper{true})(prob) + ow.f(state_values(prob), parameter_values(prob), current_time(prob)) +end + +function (ow::ObservedWrapper{false})(prob) + ow.f(state_values(prob), parameter_values(prob)) +end + """ $(TYPEDSIGNATURES) Given an index provider `indp` and a vector of symbols `syms` return a type-stable getter -function by splitting `syms` into contiguous buffers where the getter of each buffer -is type-stable and constructing a function that calls and concatenates the results. -""" -function concrete_getu(indp, syms::AbstractVector) - # a list of contiguous buffer - split_syms = [Any[syms[1]]] - # the type of the getter of the last buffer - current = typeof(getu(indp, syms[1])) - for sym in syms[2:end] - getter = getu(indp, sym) - if typeof(getter) != current - # if types don't match, build a new buffer - push!(split_syms, []) - current = typeof(getter) - end - push!(split_syms[end], sym) - end - split_syms = Tuple(split_syms) - # the getter is now type-stable, and we can vcat it to get the full buffer - return Base.Fix1(reduce, vcat) ∘ getu(indp, split_syms) +function. + +Note that the getter ONLY works for problem-like objects, since it generates an observed +function. It does NOT work for solutions. +""" +Base.@nospecializeinfer function concrete_getu(indp, syms::AbstractVector) + @nospecialize + obsfn = SymbolicIndexingInterface.observed(indp, syms) + return ObservedWrapper{is_time_dependent(indp)}(obsfn) end """ From 7fc7a41659ac3709606262dc0bbd4239d27fdd94 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Thu, 5 Jun 2025 11:44:30 +0530 Subject: [PATCH 06/42] build: bump patch version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index ad9a2e4f74..169f7abb5e 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ModelingToolkit" uuid = "961ee093-0014-501f-94e3-6117800e7a78" authors = ["Yingbo Ma ", "Chris Rackauckas and contributors"] -version = "9.80.2" +version = "9.80.3" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" From ae95c7decbbe263505b99a6744735d44b36a0413 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Thu, 5 Jun 2025 13:04:31 +0530 Subject: [PATCH 07/42] feat: allow passing `wrap_delays` to `build_explicit_observed_function` --- src/systems/diffeqs/odesystem.jl | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl index 4a13d7ccf1..334835af58 100644 --- a/src/systems/diffeqs/odesystem.jl +++ b/src/systems/diffeqs/odesystem.jl @@ -476,6 +476,8 @@ Generates a function that computes the observed value(s) `ts` in the system `sys - `throw = true` if true, throw an error when generating a function for `ts` that reference variables that do not exist. - `mkarray`: only used if the output is an array (that is, `!isscalar(ts)` and `ts` is not a tuple, in which case the result will always be a tuple). Called as `mkarray(ts, output_type)` where `ts` are the expressions to put in the array and `output_type` is the argument of the same name passed to build_explicit_observed_function. - `cse = true`: Whether to use Common Subexpression Elimination (CSE) to generate a more efficient function. +- `wrap_delays = is_dde(sys)`: Whether to add an argument for the history function and use + it to calculate all delayed variables. ## Returns @@ -514,7 +516,8 @@ function build_explicit_observed_function(sys, ts; op = Operator, throw = true, cse = true, - mkarray = nothing) + mkarray = nothing, + wrap_delays = is_dde(sys)) is_tuple = ts isa Tuple if is_tuple ts = collect(ts) @@ -600,14 +603,15 @@ function build_explicit_observed_function(sys, ts; p_end = length(dvs) + length(inputs) + length(ps) fns = build_function_wrapper( sys, ts, args...; p_start, p_end, filter_observed = obsfilter, - output_type, mkarray, try_namespaced = true, expression = Val{true}, cse) + output_type, mkarray, try_namespaced = true, expression = Val{true}, cse, + wrap_delays) if fns isa Tuple if expression return return_inplace ? fns : fns[1] end oop, iip = eval_or_rgf.(fns; eval_expression, eval_module) f = GeneratedFunctionWrapper{( - p_start + is_dde(sys), length(args) - length(ps) + 1 + is_dde(sys), is_split(sys))}( + p_start + wrap_delays, length(args) - length(ps) + 1 + wrap_delays, is_split(sys))}( oop, iip) return return_inplace ? (f, f) : f else @@ -616,7 +620,7 @@ function build_explicit_observed_function(sys, ts; end f = eval_or_rgf(fns; eval_expression, eval_module) f = GeneratedFunctionWrapper{( - p_start + is_dde(sys), length(args) - length(ps) + 1 + is_dde(sys), is_split(sys))}( + p_start + wrap_delays, length(args) - length(ps) + 1 + wrap_delays, is_split(sys))}( f, nothing) return f end From c22f2db69029a76fadda793d5388d4949671001e Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Wed, 4 Jun 2025 18:48:33 +0530 Subject: [PATCH 08/42] fix: fix `concrete_getu` for DDEs --- src/systems/problem_utils.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/systems/problem_utils.jl b/src/systems/problem_utils.jl index f79e8f6286..2f35ec69ad 100644 --- a/src/systems/problem_utils.jl +++ b/src/systems/problem_utils.jl @@ -649,7 +649,7 @@ function. It does NOT work for solutions. """ Base.@nospecializeinfer function concrete_getu(indp, syms::AbstractVector) @nospecialize - obsfn = SymbolicIndexingInterface.observed(indp, syms) + obsfn = build_explicit_observed_function(indp, syms; wrap_delays = false) return ObservedWrapper{is_time_dependent(indp)}(obsfn) end From ca8af6a646b1d2ccb7db0f04a11a5f895bfa7f49 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Thu, 5 Jun 2025 12:57:11 +0530 Subject: [PATCH 09/42] test: fix preface test --- test/odesystem.jl | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/test/odesystem.jl b/test/odesystem.jl index 83bb9b3621..6065e9e70f 100644 --- a/test/odesystem.jl +++ b/test/odesystem.jl @@ -5,10 +5,14 @@ using OrdinaryDiffEq, Sundials using DiffEqBase, SparseArrays using StaticArrays using Test -using SymbolicUtils: issym +using SymbolicUtils.Code +using SymbolicUtils: Sym, issym using ForwardDiff using ModelingToolkit: value using ModelingToolkit: t_nounits as t, D_nounits as D +using Symbolics +using Symbolics: unwrap +using DiffEqBase: isinplace # Define some variables @parameters σ ρ β @@ -607,13 +611,6 @@ sys = complete(sys) @test_throws Any ODEFunction(sys) @testset "Preface tests" begin - using OrdinaryDiffEq - using Symbolics - using DiffEqBase: isinplace - using ModelingToolkit - using SymbolicUtils.Code - using SymbolicUtils: Sym - c = [0] function f(c, du::AbstractVector{Float64}, u::AbstractVector{Float64}, p, t::Float64) c .= [c[1] + 1] @@ -656,7 +653,9 @@ sys = complete(sys) @named sys = ODESystem(eqs, t, us, ps; defaults = defs, preface = preface) sys = complete(sys) - prob = ODEProblem(sys, [], (0.0, 1.0)) + # don't build initializeprob because it will use preface in other functions and + # affect `c` + prob = ODEProblem(sys, [], (0.0, 1.0); build_initializeprob = false) sol = solve(prob, Euler(); dt = 0.1) @test c[1] == length(sol) From d77931556bece7fcfb3b400415f1b928a79c2655 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Thu, 5 Jun 2025 13:14:59 +0530 Subject: [PATCH 10/42] ci: run CI on PRs to backport-v9 --- .github/workflows/Downstream.yml | 2 +- .github/workflows/ReleaseTest.yml | 2 +- .github/workflows/Tests.yml | 2 ++ 3 files changed, 4 insertions(+), 2 deletions(-) diff --git a/.github/workflows/Downstream.yml b/.github/workflows/Downstream.yml index 6d79ee746a..0c53a7779e 100644 --- a/.github/workflows/Downstream.yml +++ b/.github/workflows/Downstream.yml @@ -1,7 +1,7 @@ name: IntegrationTest on: push: - branches: [master] + branches: [master, 'backport-v9'] tags: [v*] pull_request: paths-ignore: diff --git a/.github/workflows/ReleaseTest.yml b/.github/workflows/ReleaseTest.yml index f8b592e9d0..96b84f4d79 100644 --- a/.github/workflows/ReleaseTest.yml +++ b/.github/workflows/ReleaseTest.yml @@ -1,7 +1,7 @@ name: ReleaseTest on: push: - branches: [master] + branches: [master, 'backport-v9'] tags: [v*] pull_request: paths-ignore: diff --git a/.github/workflows/Tests.yml b/.github/workflows/Tests.yml index 52c5482970..835ee358a6 100644 --- a/.github/workflows/Tests.yml +++ b/.github/workflows/Tests.yml @@ -5,11 +5,13 @@ on: branches: - master - 'release-' + - 'backport-v9' paths-ignore: - 'docs/**' push: branches: - master + - 'backport-v9' paths-ignore: - 'docs/**' From 5bfb166f7cadc95a39539723663ddd9cae4da18a Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Thu, 5 Jun 2025 14:57:57 +0530 Subject: [PATCH 11/42] fix: fix `ObservedWrapper` for steady state systems --- src/systems/problem_utils.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/systems/problem_utils.jl b/src/systems/problem_utils.jl index 2f35ec69ad..e099ffafbe 100644 --- a/src/systems/problem_utils.jl +++ b/src/systems/problem_utils.jl @@ -631,7 +631,9 @@ end ObservedWrapper{TD}(f::F) where {TD, F} = ObservedWrapper{TD, F}(f) function (ow::ObservedWrapper{true})(prob) - ow.f(state_values(prob), parameter_values(prob), current_time(prob)) + # Edge case for steady state problems + t = applicable(current_time, prob) ? current_time(prob) : Inf + ow.f(state_values(prob), parameter_values(prob), t) end function (ow::ObservedWrapper{false})(prob) From e9b27f9e3a6ef98d15d78cf480caf2459f67a52d Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Thu, 5 Jun 2025 14:58:25 +0530 Subject: [PATCH 12/42] fix: fix `observed_equations_used_by` for `Initial`s --- src/utils.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils.jl b/src/utils.jl index 90431a7749..6fa574ba35 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1083,7 +1083,7 @@ Keyword arguments: `available_vars` will not be searched for in the observed equations. """ function observed_equations_used_by(sys::AbstractSystem, exprs; - involved_vars = vars(exprs; op = Union{Shift, Differential}), obs = observed(sys), available_vars = []) + involved_vars = vars(exprs; op = Union{Shift, Differential, Initial}), obs = observed(sys), available_vars = []) obsvars = getproperty.(obs, :lhs) graph = observed_dependency_graph(obs) if !(available_vars isa Set) From ce655ae279f14c206df589e6ca26b82302b0aea3 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Thu, 5 Jun 2025 12:27:37 +0530 Subject: [PATCH 13/42] refactor: remove CSE hack --- .../symbolics_tearing.jl | 97 +++---------------- src/systems/nonlinear/initializesystem.jl | 14 --- test/code_generation.jl | 31 ++++++ test/structural_transformation/utils.jl | 73 +------------- 4 files changed, 47 insertions(+), 168 deletions(-) diff --git a/src/structural_transformation/symbolics_tearing.jl b/src/structural_transformation/symbolics_tearing.jl index 3f46b594df..cc582cd473 100644 --- a/src/structural_transformation/symbolics_tearing.jl +++ b/src/structural_transformation/symbolics_tearing.jl @@ -724,7 +724,7 @@ Update the system equations, unknowns, and observables after simplification. """ function update_simplified_system!( state::TearingState, neweqs, solved_eqs, dummy_sub, var_eq_matching, extra_unknowns; - cse_hack = true, array_hack = true) + array_hack = true) @unpack solvable_graph, var_to_diff, eq_to_diff, graph = state.structure diff_to_var = invview(var_to_diff) @@ -748,8 +748,7 @@ function update_simplified_system!( unknowns = [unknowns; extra_unknowns] @set! sys.unknowns = unknowns - obs = cse_and_array_hacks( - sys, obs, unknowns, neweqs; cse = cse_hack, array = array_hack) + obs = tearing_hacks(sys, obs, unknowns, neweqs; array = array_hack) deps = Vector{Int}[i == 1 ? Int[] : collect(1:(i - 1)) for i in 1:length(solved_eqs)] @@ -793,7 +792,7 @@ appear in the system. Algebraic variables are variables that are not differential variables. """ function tearing_reassemble(state::TearingState, var_eq_matching, - full_var_eq_matching = nothing; simplify = false, mm = nothing, cse_hack = true, array_hack = true) + full_var_eq_matching = nothing; simplify = false, mm = nothing, array_hack = true) extra_vars = Int[] if full_var_eq_matching !== nothing for v in 𝑑vertices(state.structure.graph) @@ -829,7 +828,7 @@ function tearing_reassemble(state::TearingState, var_eq_matching, state, var_eq_matching, eq_ordering, var_ordering, nelim_eq, nelim_var) sys = update_simplified_system!(state, neweqs, solved_eqs, dummy_sub, var_eq_matching, - extra_unknowns; cse_hack, array_hack) + extra_unknowns; array_hack) @set! state.sys = sys @set! sys.tearing_state = state @@ -837,14 +836,7 @@ function tearing_reassemble(state::TearingState, var_eq_matching, end """ -# HACK 1 - -Since we don't support array equations, any equation of the sort `x[1:n] ~ f(...)[1:n]` -gets turned into `x[1] ~ f(...)[1], x[2] ~ f(...)[2]`. Repeatedly calling `f` gets -_very_ expensive. this hack performs a limited form of CSE specifically for this case to -avoid the unnecessary cost. This and the below hack are implemented simultaneously - -# HACK 2 +# HACK Add equations for array observed variables. If `p[i] ~ (...)` are equations, add an equation `p ~ [p[1], p[2], ...]` allow topsort to reorder them only add the new equation @@ -852,13 +844,7 @@ if all `p[i]` are present and the unscalarized form is used in any equation (obs not) we first count the number of times the scalarized form of each observed variable occurs in observed equations (and unknowns if it's split). """ -function cse_and_array_hacks(sys, obs, unknowns, neweqs; cse = true, array = true) - # HACK 1 - # mapping of rhs to temporary CSE variable - # `f(...) => tmpvar` in above example - rhs_to_tempvar = Dict() - - # HACK 2 +function tearing_hacks(sys, obs, unknowns, neweqs; array = true) # map of array observed variable (unscalarized) to number of its # scalarized terms that appear in observed equations arr_obs_occurrences = Dict() @@ -866,31 +852,6 @@ function cse_and_array_hacks(sys, obs, unknowns, neweqs; cse = true, array = tru lhs = eq.lhs rhs = eq.rhs - # HACK 1 - if cse && is_getindexed_array(rhs) - rhs_arr = arguments(rhs)[1] - iscall(rhs_arr) && operation(rhs_arr) isa Symbolics.Operator && continue - if !haskey(rhs_to_tempvar, rhs_arr) - tempvar = gensym(Symbol(lhs)) - N = length(rhs_arr) - tempvar = unwrap(Symbolics.variable( - tempvar; T = Symbolics.symtype(rhs_arr))) - tempvar = setmetadata( - tempvar, Symbolics.ArrayShapeCtx, Symbolics.shape(rhs_arr)) - tempeq = tempvar ~ rhs_arr - rhs_to_tempvar[rhs_arr] = tempvar - push!(obs, tempeq) - end - - # getindex_wrapper is used because `observed2graph` treats `x` and `x[i]` as different, - # so it doesn't find a dependency between this equation and `tempvar ~ rhs_arr` - # which fails the topological sort - neweq = lhs ~ getindex_wrapper( - rhs_to_tempvar[rhs_arr], Tuple(arguments(rhs)[2:end])) - obs[i] = neweq - end - # end HACK 1 - array || continue iscall(lhs) || continue operation(lhs) === getindex || continue @@ -901,31 +862,6 @@ function cse_and_array_hacks(sys, obs, unknowns, neweqs; cse = true, array = tru continue end - # Also do CSE for `equations(sys)` - if cse - for (i, eq) in enumerate(neweqs) - (; lhs, rhs) = eq - is_getindexed_array(rhs) || continue - rhs_arr = arguments(rhs)[1] - if !haskey(rhs_to_tempvar, rhs_arr) - tempvar = gensym(Symbol(lhs)) - N = length(rhs_arr) - tempvar = unwrap(Symbolics.variable( - tempvar; T = Symbolics.symtype(rhs_arr))) - tempvar = setmetadata( - tempvar, Symbolics.ArrayShapeCtx, Symbolics.shape(rhs_arr)) - tempeq = tempvar ~ rhs_arr - rhs_to_tempvar[rhs_arr] = tempvar - push!(obs, tempeq) - end - # don't need getindex_wrapper, but do it anyway to know that this - # hack took place - neweq = lhs ~ getindex_wrapper( - rhs_to_tempvar[rhs_arr], Tuple(arguments(rhs)[2:end])) - neweqs[i] = neweq - end - end - # count variables in unknowns if they are scalarized forms of variables # also present as observed. e.g. if `x[1]` is an unknown and `x[2] ~ (..)` # is an observed equation. @@ -960,18 +896,7 @@ function cse_and_array_hacks(sys, obs, unknowns, neweqs; cse = true, array = tru return obs end -function is_getindexed_array(rhs) - (!ModelingToolkit.isvariable(rhs) || ModelingToolkit.iscalledparameter(rhs)) && - iscall(rhs) && operation(rhs) === getindex && - Symbolics.shape(rhs) != Symbolics.Unknown() -end - -# PART OF HACK 1 -getindex_wrapper(x, i) = x[i...] - -@register_symbolic getindex_wrapper(x::AbstractArray, i::Tuple{Vararg{Int}}) - -# PART OF HACK 2 +# PART OF HACK function change_origin(origin, arr) if all(isone, Tuple(origin)) return arr @@ -999,10 +924,10 @@ new residual equations after tearing. End users are encouraged to call [`structu instead, which calls this function internally. """ function tearing(sys::AbstractSystem, state = TearingState(sys); mm = nothing, - simplify = false, cse_hack = true, array_hack = true, kwargs...) + simplify = false, array_hack = true, kwargs...) var_eq_matching, full_var_eq_matching = tearing(state) invalidate_cache!(tearing_reassemble( - state, var_eq_matching, full_var_eq_matching; mm, simplify, cse_hack, array_hack)) + state, var_eq_matching, full_var_eq_matching; mm, simplify, array_hack)) end """ @@ -1024,7 +949,7 @@ Perform index reduction and use the dummy derivative technique to ensure that the system is balanced. """ function dummy_derivative(sys, state = TearingState(sys); simplify = false, - mm = nothing, cse_hack = true, array_hack = true, kwargs...) + mm = nothing, array_hack = true, kwargs...) jac = let state = state (eqs, vars) -> begin symeqs = EquationsView(state)[eqs] @@ -1048,5 +973,5 @@ function dummy_derivative(sys, state = TearingState(sys); simplify = false, end var_eq_matching = dummy_derivative_graph!(state, jac; state_priority, kwargs...) - tearing_reassemble(state, var_eq_matching; simplify, mm, cse_hack, array_hack) + tearing_reassemble(state, var_eq_matching; simplify, mm, array_hack) end diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl index 4b0086840a..cf38b94687 100644 --- a/src/systems/nonlinear/initializesystem.jl +++ b/src/systems/nonlinear/initializesystem.jl @@ -739,20 +739,6 @@ function unhack_observed(obseqs::Vector{Equation}, eqs::Vector{Equation}) push!(rm_idxs, i) continue end - if operation(eq.rhs) == StructuralTransformations.getindex_wrapper - var, idxs = arguments(eq.rhs) - subs[eq.rhs] = var[idxs...] - push!(tempvars, var) - end - end - - for (i, eq) in enumerate(eqs) - iscall(eq.rhs) || continue - if operation(eq.rhs) == StructuralTransformations.getindex_wrapper - var, idxs = arguments(eq.rhs) - subs[eq.rhs] = var[idxs...] - push!(tempvars, var) - end end for (i, eq) in enumerate(obseqs) diff --git a/test/code_generation.jl b/test/code_generation.jl index cf3d660b81..05b945f503 100644 --- a/test/code_generation.jl +++ b/test/code_generation.jl @@ -78,3 +78,34 @@ end @test SciMLBase.successful_retcode(sol) end end + +@testset "scalarized array observed calling same function multiple times" begin + @variables x(t) y(t)[1:2] + @parameters foo(::Real)[1:2] + val = Ref(0) + function _tmp_fn2(x) + val[] += 1 + return [x, 2x] + end + @mtkbuild sys = ODESystem([D(x) ~ y[1] + y[2], y ~ foo(x)], t) + @test length(equations(sys)) == 1 + @test length(observed(sys)) == 3 + prob = ODEProblem(sys, [x => 1.0], (0.0, 1.0), [foo => _tmp_fn2]) + val[] = 0 + @test_nowarn prob.f(prob.u0, prob.p, 0.0) + @test val[] == 1 + + @testset "CSE in equations(sys)" begin + val[] = 0 + @variables z(t)[1:2] + @mtkbuild sys = ODESystem( + [D(y) ~ foo(x), D(x) ~ sum(y), zeros(2) ~ foo(prod(z))], t) + @test length(equations(sys)) == 5 + @test length(observed(sys)) == 0 + prob = ODEProblem( + sys, [y => ones(2), z => 2ones(2), x => 3.0], (0.0, 1.0), [foo => _tmp_fn2]) + val[] = 0 + @test_nowarn prob.f(prob.u0, prob.p, 0.0) + @test val[] == 2 + end +end diff --git a/test/structural_transformation/utils.jl b/test/structural_transformation/utils.jl index 67cf2f72a0..24fb166755 100644 --- a/test/structural_transformation/utils.jl +++ b/test/structural_transformation/utils.jl @@ -51,7 +51,7 @@ end @mtkbuild sys = ODESystem( [D(x) ~ z[1] + z[2] + foo(z)[1], y[1] ~ 2t, y[2] ~ 3t, z ~ foo(y)], t) @test length(equations(sys)) == 1 - @test length(observed(sys)) == 7 + @test length(observed(sys)) == 6 @test any(obs -> isequal(obs, y), observables(sys)) @test any(obs -> isequal(obs, z), observables(sys)) prob = ODEProblem(sys, [x => 1.0], (0.0, 1.0), [foo => _tmp_fn]) @@ -61,60 +61,11 @@ end @test length(unknowns(isys)) == 5 @test length(equations(isys)) == 4 @test !any(equations(isys)) do eq - iscall(eq.rhs) && operation(eq.rhs) in [StructuralTransformations.getindex_wrapper, - StructuralTransformations.change_origin] + iscall(eq.rhs) && operation(eq.rhs) in [StructuralTransformations.change_origin] end end -@testset "scalarized array observed calling same function multiple times" begin - @variables x(t) y(t)[1:2] - @parameters foo(::Real)[1:2] - val = Ref(0) - function _tmp_fn2(x) - val[] += 1 - return [x, 2x] - end - @mtkbuild sys = ODESystem([D(x) ~ y[1] + y[2], y ~ foo(x)], t) - @test length(equations(sys)) == 1 - @test length(observed(sys)) == 4 - prob = ODEProblem(sys, [x => 1.0], (0.0, 1.0), [foo => _tmp_fn2]) - val[] = 0 - @test_nowarn prob.f(prob.u0, prob.p, 0.0) - @test val[] == 1 - - isys = ModelingToolkit.generate_initializesystem(sys) - @test length(unknowns(isys)) == 3 - @test length(equations(isys)) == 2 - @test !any(equations(isys)) do eq - iscall(eq.rhs) && operation(eq.rhs) in [StructuralTransformations.getindex_wrapper, - StructuralTransformations.change_origin] - end - - @testset "CSE hack in equations(sys)" begin - val[] = 0 - @variables z(t)[1:2] - @mtkbuild sys = ODESystem( - [D(y) ~ foo(x), D(x) ~ sum(y), zeros(2) ~ foo(prod(z))], t) - @test length(equations(sys)) == 5 - @test length(observed(sys)) == 2 - prob = ODEProblem( - sys, [y => ones(2), z => 2ones(2), x => 3.0], (0.0, 1.0), [foo => _tmp_fn2]) - val[] = 0 - @test_nowarn prob.f(prob.u0, prob.p, 0.0) - @test val[] == 2 - - isys = ModelingToolkit.generate_initializesystem(sys) - @test length(unknowns(isys)) == 5 - @test length(equations(isys)) == 2 - @test !any(equations(isys)) do eq - iscall(eq.rhs) && - operation(eq.rhs) in [StructuralTransformations.getindex_wrapper, - StructuralTransformations.change_origin] - end - end -end - -@testset "array and cse hacks can be disabled" begin +@testset "array hack can be disabled" begin @testset "fully_determined = true" begin @variables x(t) y(t)[1:2] z(t)[1:2] @parameters foo(::AbstractVector)[1:2] @@ -122,15 +73,8 @@ end @named sys = ODESystem( [D(x) ~ z[1] + z[2] + foo(z)[1], y[1] ~ 2t, y[2] ~ 3t, z ~ foo(y)], t) - sys1 = structural_simplify(sys; cse_hack = false) - @test length(observed(sys1)) == 6 - @test !any(observed(sys1)) do eq - iscall(eq.rhs) && - operation(eq.rhs) == StructuralTransformations.getindex_wrapper - end - sys2 = structural_simplify(sys; array_hack = false) - @test length(observed(sys2)) == 5 + @test length(observed(sys2)) == 4 @test !any(observed(sys2)) do eq iscall(eq.rhs) && operation(eq.rhs) == StructuralTransformations.change_origin end @@ -143,15 +87,8 @@ end @named sys = ODESystem( [D(x) ~ z[1] + z[2] + foo(z)[1] + w, y[1] ~ 2t, y[2] ~ 3t, z ~ foo(y)], t) - sys1 = structural_simplify(sys; cse_hack = false, fully_determined = false) - @test length(observed(sys1)) == 6 - @test !any(observed(sys1)) do eq - iscall(eq.rhs) && - operation(eq.rhs) == StructuralTransformations.getindex_wrapper - end - sys2 = structural_simplify(sys; array_hack = false, fully_determined = false) - @test length(observed(sys2)) == 5 + @test length(observed(sys2)) == 4 @test !any(observed(sys2)) do eq iscall(eq.rhs) && operation(eq.rhs) == StructuralTransformations.change_origin end From 74c38c5ae0a4990a26a29c58d6c6257727010336 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Thu, 5 Jun 2025 12:32:06 +0530 Subject: [PATCH 14/42] test: test array unknowns occurring unscalarized in initializeprobpmap --- test/initializationsystem.jl | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl index 7c512d37af..b1ba21b254 100644 --- a/test/initializationsystem.jl +++ b/test/initializationsystem.jl @@ -1650,3 +1650,24 @@ end @test !SciMLBase.isinplace(prob) @test !SciMLBase.isinplace(prob.f.initialization_data.initializeprob) end + +@testset "Array unknowns occurring unscalarized in initializeprobpmap" begin + @variables begin + u(t)[1:2] = 0.9ones(2) + x(t)[1:2], [guess = 0.01ones(2)] + o(t)[1:2] + end + @parameters p[1:4] = [2.0, 1.875, 2.0, 1.875] + + eqs = [D(u[1]) ~ p[1] * u[1] - p[2] * u[1] * u[2] + x[1] + 0.1 + D(u[2]) ~ p[4] * u[1] * u[2] - p[3] * u[2] - x[2] + o[1] ~ sum(p) * sum(u) + o[2] ~ sum(p) * sum(x) + x[1] ~ 0.01exp(-1) + x[2] ~ 0.01cos(t)] + + @mtkbuild sys = ODESystem(eqs, t) + prob = ODEProblem(sys, [], (0.0, 1.0)) + sol = solve(prob, Tsit5()) + @test SciMLBase.successful_retcode(sol) +end From 263c8701cdbab1c385847bf32e233474086ce997 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Thu, 5 Jun 2025 18:41:20 +0530 Subject: [PATCH 15/42] build: bump patch version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 169f7abb5e..4e9999361b 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ModelingToolkit" uuid = "961ee093-0014-501f-94e3-6117800e7a78" authors = ["Yingbo Ma ", "Chris Rackauckas and contributors"] -version = "9.80.3" +version = "9.80.4" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" From 6e183b8e16ab86b944bc711505d56b073ea6d017 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Fri, 6 Jun 2025 15:22:42 +0530 Subject: [PATCH 16/42] fix: fix `get_mtkparameters_reconstructor` handling of nonnumerics --- src/systems/problem_utils.jl | 40 ++++++++++++++++++++++++++++-------- 1 file changed, 32 insertions(+), 8 deletions(-) diff --git a/src/systems/problem_utils.jl b/src/systems/problem_utils.jl index e099ffafbe..b893d9ffc6 100644 --- a/src/systems/problem_utils.jl +++ b/src/systems/problem_utils.jl @@ -659,7 +659,8 @@ end $(TYPEDEF) A callable struct which applies `p_constructor` to possibly nested arrays. It also -ensures that views (including nested ones) are concretized. +ensures that views (including nested ones) are concretized. This is implemented manually +of using `narrow_buffer_type` to preserve type-stability. """ struct PConstructorApplicator{F} p_constructor::F @@ -669,10 +670,18 @@ function (pca::PConstructorApplicator)(x::AbstractArray) pca.p_constructor(x) end +function (pca::PConstructorApplicator)(x::AbstractArray{Bool}) + pca.p_constructor(BitArray(x)) +end + function (pca::PConstructorApplicator{typeof(identity)})(x::SubArray) collect(x) end +function (pca::PConstructorApplicator{typeof(identity)})(x::SubArray{Bool}) + BitArray(x) +end + function (pca::PConstructorApplicator{typeof(identity)})(x::SubArray{<:AbstractArray}) collect(pca.(x)) end @@ -695,6 +704,7 @@ takes a value provider of `srcsys` and a value provider of `dstsys` and returns """ function get_mtkparameters_reconstructor(srcsys::AbstractSystem, dstsys::AbstractSystem; initials = false, unwrap_initials = false, p_constructor = identity) + _p_constructor = p_constructor p_constructor = PConstructorApplicator(p_constructor) # if we call `getu` on this (and it were able to handle empty tuples) we get the # fields of `MTKParameters` except caches. @@ -748,14 +758,24 @@ function get_mtkparameters_reconstructor(srcsys::AbstractSystem, dstsys::Abstrac Base.Fix1(broadcast, p_constructor) ∘ getu(srcsys, syms[3]) end - rest_getters = map(Base.tail(Base.tail(Base.tail(syms)))) do buf - if buf == () - return Returns(()) - else - return Base.Fix1(broadcast, p_constructor) ∘ getu(srcsys, buf) - end + const_getter = if syms[4] == () + Returns(()) + else + Base.Fix1(broadcast, p_constructor) ∘ getu(srcsys, syms[4]) end - getters = (tunable_getter, initials_getter, discs_getter, rest_getters...) + nonnumeric_getter = if syms[5] == () + Returns(()) + else + ic = get_index_cache(dstsys) + buftypes = Tuple(map(ic.nonnumeric_buffer_sizes) do bufsize + Vector{bufsize.type} + end) + # nonnumerics retain the assigned buffer type without narrowing + Base.Fix1(broadcast, _p_constructor) ∘ + Base.Fix1(Broadcast.BroadcastFunction(call), buftypes) ∘ getu(srcsys, syms[5]) + end + getters = ( + tunable_getter, initials_getter, discs_getter, const_getter, nonnumeric_getter) getter = let getters = getters function _getter(valp, initprob) oldcache = parameter_values(initprob).caches @@ -768,6 +788,10 @@ function get_mtkparameters_reconstructor(srcsys::AbstractSystem, dstsys::Abstrac return getter end +function call(f, args...) + f(args...) +end + """ $(TYPEDSIGNATURES) From 713a0928b6e6af662da654ea5ecf367b7f7f298b Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Fri, 6 Jun 2025 14:23:30 +0530 Subject: [PATCH 17/42] test: test nonnumerics aren't narrowed in `ODEProblem` and `init` --- test/initializationsystem.jl | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl index b1ba21b254..fd84c346f9 100644 --- a/test/initializationsystem.jl +++ b/test/initializationsystem.jl @@ -1671,3 +1671,23 @@ end sol = solve(prob, Tsit5()) @test SciMLBase.successful_retcode(sol) end + +@testset "Nonnumerics aren't narrowed" begin + @mtkmodel Foo begin + @variables begin + x(t) = 1.0 + end + @parameters begin + p::AbstractString + r = 1.0 + end + @equations begin + D(x) ~ r * x + end + end + @mtkbuild sys = Foo(p = "a") + prob = ODEProblem(sys, [], (0.0, 1.0)) + @test prob.p.nonnumeric[1] isa Vector{AbstractString} + integ = init(prob) + @test integ.p.nonnumeric[1] isa Vector{AbstractString} +end From 81596b4528bb2bcff04c84f4dae21c5187051459 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Fri, 6 Jun 2025 15:23:05 +0530 Subject: [PATCH 18/42] test: fix codegen test --- test/code_generation.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/code_generation.jl b/test/code_generation.jl index 05b945f503..b48e4acf64 100644 --- a/test/code_generation.jl +++ b/test/code_generation.jl @@ -89,7 +89,7 @@ end end @mtkbuild sys = ODESystem([D(x) ~ y[1] + y[2], y ~ foo(x)], t) @test length(equations(sys)) == 1 - @test length(observed(sys)) == 3 + @test length(ModelingToolkit.observed(sys)) == 3 prob = ODEProblem(sys, [x => 1.0], (0.0, 1.0), [foo => _tmp_fn2]) val[] = 0 @test_nowarn prob.f(prob.u0, prob.p, 0.0) @@ -101,7 +101,7 @@ end @mtkbuild sys = ODESystem( [D(y) ~ foo(x), D(x) ~ sum(y), zeros(2) ~ foo(prod(z))], t) @test length(equations(sys)) == 5 - @test length(observed(sys)) == 0 + @test length(ModelingToolkit.observed(sys)) == 0 prob = ODEProblem( sys, [y => ones(2), z => 2ones(2), x => 3.0], (0.0, 1.0), [foo => _tmp_fn2]) val[] = 0 From 897524e372ae2b324e6bf22e38e2e8f05dfcf664 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Fri, 6 Jun 2025 11:06:26 +0000 Subject: [PATCH 19/42] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 4e9999361b..8b27d2f2df 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ModelingToolkit" uuid = "961ee093-0014-501f-94e3-6117800e7a78" authors = ["Yingbo Ma ", "Chris Rackauckas and contributors"] -version = "9.80.4" +version = "9.80.5" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" From a91cb2272348cf30662cbb3e9281f39681636f73 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Fri, 20 Jun 2025 16:07:55 +0530 Subject: [PATCH 20/42] fix: fix `get_updated_symbolic_problem` Respect `u0` and `p` passed as kwargs --- src/systems/nonlinear/initializesystem.jl | 10 +++--- test/extensions/Project.toml | 1 + test/extensions/ad.jl | 44 +++++++++++++++++++++++ 3 files changed, 50 insertions(+), 5 deletions(-) diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl index cf38b94687..5a8687984f 100644 --- a/src/systems/nonlinear/initializesystem.jl +++ b/src/systems/nonlinear/initializesystem.jl @@ -685,7 +685,9 @@ function SciMLBase.late_binding_update_u0_p( return newu0, newp end -function DiffEqBase.get_updated_symbolic_problem(sys::AbstractSystem, prob; kw...) +function DiffEqBase.get_updated_symbolic_problem( + sys::AbstractSystem, prob; u0 = state_values(prob), + p = parameter_values(prob), kw...) supports_initialization(sys) || return prob initdata = prob.f.initialization_data initdata isa SciMLBase.OverrideInitData || return prob @@ -693,10 +695,8 @@ function DiffEqBase.get_updated_symbolic_problem(sys::AbstractSystem, prob; kw.. meta isa InitializationMetadata || return prob meta.get_updated_u0 === nothing && return prob - u0 = state_values(prob) - u0 === nothing && return prob + u0 === nothing && return remake(prob; p) - p = parameter_values(prob) t0 = is_time_dependent(prob) ? current_time(prob) : nothing if p isa MTKParameters @@ -713,7 +713,7 @@ function DiffEqBase.get_updated_symbolic_problem(sys::AbstractSystem, prob; kw.. T = StaticArrays.similar_type(u0) end - return remake(prob; u0 = T(meta.get_updated_u0(prob, initdata.initializeprob))) + return remake(prob; u0 = T(meta.get_updated_u0(prob, initdata.initializeprob)), p) end """ diff --git a/test/extensions/Project.toml b/test/extensions/Project.toml index 9f43e6f4a4..f1647aef1a 100644 --- a/test/extensions/Project.toml +++ b/test/extensions/Project.toml @@ -23,6 +23,7 @@ OrdinaryDiffEqVerner = "79d7bb75-1356-48c1-b8c0-6832512096c2" SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1" SciMLStructures = "53ae85a6-f571-4167-b2af-e1d143709226" SimpleDiffEq = "05bca326-078c-5bf0-a5bf-ce7c7982d7fd" +StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" diff --git a/test/extensions/ad.jl b/test/extensions/ad.jl index 14649b6bb6..b5e71bf0af 100644 --- a/test/extensions/ad.jl +++ b/test/extensions/ad.jl @@ -8,6 +8,7 @@ using OrdinaryDiffEqNonlinearSolve using NonlinearSolve using SciMLSensitivity using ForwardDiff +using StableRNGs using ChainRulesCore using ChainRulesCore: NoTangent using ChainRulesTestUtils: test_rrule, rand_tangent @@ -135,3 +136,46 @@ end prob[sys.x] end end + +@testset "`p` provided to `solve` is respected" begin + @mtkmodel Linear begin + @variables begin + x(t) = 1.0, [description = "Prey"] + end + @parameters begin + α = 1.5 + end + @equations begin + D(x) ~ -α * x + end + end + + @mtkbuild linear = Linear() + problem = ODEProblem(linear, [], (0.0, 1.0)) + solution = solve(problem, Tsit5(), saveat = 0.1) + rng = StableRNG(42) + data = (; + t = solution.t, + # [[y, x], :] + measurements = Array(solution) + ) + data.measurements .+= 0.05 * randn(rng, size(data.measurements)) + + p0, repack, _ = SciMLStructures.canonicalize(SciMLStructures.Tunable(), problem.p) + + objective = let repack = repack, problem = problem + (p, data) -> begin + pnew = repack(p) + sol = solve(problem, Tsit5(), p = pnew, saveat = data.t) + sum(abs2, sol .- data.measurements) / size(data.t, 1) + end + end + + # Check 0.0031677344878386607 + @test_nowarn objective(p0, data) + + fd = ForwardDiff.gradient(Base.Fix2(objective, data), p0) + zg = Zygote.gradient(Base.Fix2(objective, data), p0) + + @test fd≈zg[1] atol=1e-6 +end From 4578b669dc9b7270d787b8ad0dd349564bc58253 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 23 Jun 2025 12:20:07 +0530 Subject: [PATCH 21/42] build: bump patch version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 8b27d2f2df..130e0541f3 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ModelingToolkit" uuid = "961ee093-0014-501f-94e3-6117800e7a78" authors = ["Yingbo Ma ", "Chris Rackauckas and contributors"] -version = "9.80.5" +version = "9.80.6" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" From 278fe05d982ed890ab7a6d3f68bf11443dd8f1ff Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Tue, 24 Jun 2025 22:08:54 +0200 Subject: [PATCH 22/42] fix: enable support for complex ODEProblem again --- src/systems/index_cache.jl | 4 ++-- src/systems/problem_utils.jl | 6 +++--- test/complex.jl | 28 ++++++++++++++++++++++++++++ 3 files changed, 33 insertions(+), 5 deletions(-) diff --git a/src/systems/index_cache.jl b/src/systems/index_cache.jl index d0b687c212..ae4feec62b 100644 --- a/src/systems/index_cache.jl +++ b/src/systems/index_cache.jl @@ -388,8 +388,8 @@ function IndexCache(sys::AbstractSystem) observed_syms_to_timeseries, dependent_pars_to_timeseries, disc_buffer_templates, - BufferTemplate(Real, tunable_buffer_size), - BufferTemplate(Real, initials_buffer_size), + BufferTemplate(Number, tunable_buffer_size), + BufferTemplate(Number, initials_buffer_size), const_buffer_sizes, nonnumeric_buffer_sizes, symbol_to_variable diff --git a/src/systems/problem_utils.jl b/src/systems/problem_utils.jl index b893d9ffc6..ec72e41fe7 100644 --- a/src/systems/problem_utils.jl +++ b/src/systems/problem_utils.jl @@ -962,7 +962,7 @@ end $(TYPEDEF) A callable struct to use as the `get_updated_u0` field of `InitializationMetadata`. -Returns the value to use for the `u0` of the problem. +Returns the value to use for the `u0` of the problem. # Fields @@ -1160,7 +1160,7 @@ function float_type_from_varmap(varmap, floatT = Bool) if v isa AbstractArray floatT = promote_type(floatT, eltype(v)) - elseif v isa Real + elseif v isa Number floatT = promote_type(floatT, typeof(v)) end end @@ -1432,7 +1432,7 @@ function check_inputmap_keys(sys, u0map, pmap) end const BAD_KEY_MESSAGE = """ - Undefined keys found in the parameter or initial condition maps. Check if symbolic variable names have been reassigned. + Undefined keys found in the parameter or initial condition maps. Check if symbolic variable names have been reassigned. The following keys are invalid: """ diff --git a/test/complex.jl b/test/complex.jl index 69cc22c985..04be8e4dac 100644 --- a/test/complex.jl +++ b/test/complex.jl @@ -1,4 +1,5 @@ using ModelingToolkit +using OrdinaryDiffEq using ModelingToolkit: t_nounits as t using Test @@ -14,3 +15,30 @@ using Test end @named mixed = ComplexModel() @test length(equations(mixed)) == 2 + +@testset "Complex ODEProblem" begin + using ModelingToolkit: t_nounits as t, D_nounits as D + + vars = @variables x(t) y(t) z(t) + pars = @parameters a b + + eqs = [ + D(x) ~ y - x, + D(y) ~ -x * z + b * abs(z), + D(z) ~ x * y - a + ] + @named modlorenz = System(eqs, t) + sys = structural_simplify(modlorenz) + + ic = ModelingToolkit.get_index_cache(sys) + @test ic.tunable_buffer_size.type == Number + + u0 = ComplexF64[-4.0, 5.0, 0.0] .+ randn(ComplexF64, 3) + p = ComplexF64[5.0, 0.1] + dict = merge(Dict(unknowns(sys) .=> u0), Dict(parameters(sys) .=> p)) + prob = ODEProblem(sys, dict, (0.0, 1.0)) + + sol = solve(prob, Tsit5(), saveat = 0.1) + + @test sol.u[1] isa Vector{ComplexF64} +end From ff98434602b92993880e47e2c13fa43c6e266137 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Tue, 17 Jun 2025 16:57:17 +0530 Subject: [PATCH 23/42] feat: add fast path for constructing `MTKParameters` in `process_SciMLProblem` --- src/systems/parameter_buffer.jl | 14 ++++++++++---- src/systems/problem_utils.jl | 2 +- 2 files changed, 11 insertions(+), 5 deletions(-) diff --git a/src/systems/parameter_buffer.jl b/src/systems/parameter_buffer.jl index c3d2a0e831..00742f028c 100644 --- a/src/systems/parameter_buffer.jl +++ b/src/systems/parameter_buffer.jl @@ -29,7 +29,7 @@ the default behavior). function MTKParameters( sys::AbstractSystem, p, u0 = Dict(); tofloat = false, t0 = nothing, substitution_limit = 1000, floatT = nothing, - p_constructor = identity) + p_constructor = identity, fast_path = false) ic = if has_index_cache(sys) && get_index_cache(sys) !== nothing get_index_cache(sys) else @@ -50,9 +50,15 @@ function MTKParameters( is_time_dependent(sys) && add_observed!(sys, u0) add_parameter_dependencies!(sys, p) - op, missing_unknowns, missing_pars = build_operating_point!(sys, - u0, p, defs, cmap, dvs, ps) - + u0map = anydict() + pmap = anydict() + if fast_path + missing_pars = missingvars(p, ps) + op = p + else + op, _, missing_pars = build_operating_point!(sys, + u0, p, defs, cmap, dvs, ps) + end if t0 !== nothing op[get_iv(sys)] = t0 end diff --git a/src/systems/problem_utils.jl b/src/systems/problem_utils.jl index b893d9ffc6..e1d2059657 100644 --- a/src/systems/problem_utils.jl +++ b/src/systems/problem_utils.jl @@ -1371,7 +1371,7 @@ function process_SciMLProblem( if !(pType <: AbstractArray) pType = Array end - p = MTKParameters(sys, op; floatT = floatT, p_constructor) + p = MTKParameters(sys, op; floatT = floatT, p_constructor, fast_path = true) else p = p_constructor(better_varmap_to_vars(op, ps; tofloat, container_type = pType)) end From 95aa8801ce75e4673b443c0527611cde480387c4 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Tue, 17 Jun 2025 16:57:38 +0530 Subject: [PATCH 24/42] refactor: avoid unnecessary computation in `evaluate_varmap!` --- src/systems/problem_utils.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/systems/problem_utils.jl b/src/systems/problem_utils.jl index e1d2059657..5fdd61d903 100644 --- a/src/systems/problem_utils.jl +++ b/src/systems/problem_utils.jl @@ -450,8 +450,11 @@ in `varmap`, it is ignored. """ function evaluate_varmap!(varmap::AbstractDict, vars; limit = 100) for k in vars + v = get(varmap, k, nothing) + v === nothing && continue + symbolic_type(v) == NotSymbolic() && !is_array_of_symbolics(v) && continue haskey(varmap, k) || continue - varmap[k] = fixpoint_sub(varmap[k], varmap; maxiters = limit) + varmap[k] = fixpoint_sub(v, varmap; maxiters = limit) end end From 3b48e63c198af2c4bc674a9c24109ca8a7bd1e96 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Tue, 17 Jun 2025 16:58:00 +0530 Subject: [PATCH 25/42] refactor: avoid unnecesary computation in `build_operating_point!` --- src/systems/problem_utils.jl | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/systems/problem_utils.jl b/src/systems/problem_utils.jl index 5fdd61d903..2482ff3b9c 100644 --- a/src/systems/problem_utils.jl +++ b/src/systems/problem_utils.jl @@ -583,13 +583,15 @@ function build_operating_point!(sys::AbstractSystem, end end - for k in keys(u0map) - v = fixpoint_sub(u0map[k], neithermap; operator = Symbolics.Operator) + for (k, v) in u0map + symbolic_type(v) == NotSymbolic() && !is_array_of_symbolics(v) && continue + v = fixpoint_sub(v, neithermap; operator = Symbolics.Operator) isequal(k, v) && continue u0map[k] = v end - for k in keys(pmap) - v = fixpoint_sub(pmap[k], neithermap; operator = Symbolics.Operator) + for (k, v) in pmap + symbolic_type(v) == NotSymbolic() && !is_array_of_symbolics(v) && continue + v = fixpoint_sub(v, neithermap; operator = Symbolics.Operator) isequal(k, v) && continue pmap[k] = v end From 40f0587c4a5e39529eb6477d7b337b16af83d419 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Tue, 17 Jun 2025 16:58:29 +0530 Subject: [PATCH 26/42] refactor: batch computation of temporary values in `maybe_build_initialization_problem` --- src/systems/problem_utils.jl | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/src/systems/problem_utils.jl b/src/systems/problem_utils.jl index 2482ff3b9c..377b3bdc7c 100644 --- a/src/systems/problem_utils.jl +++ b/src/systems/problem_utils.jl @@ -1128,20 +1128,24 @@ function maybe_build_initialization_problem( update_initializeprob! = ModelingToolkit.update_initializeprob! end - for p in punknowns - is_parameter_solvable(p, pmap, defs, guesses) || continue - get(op, p, missing) === missing || continue + filter!(punknowns) do p + is_parameter_solvable(p, op, defs, guesses) && get(op, p, missing) === missing + end + pvals = getu(initializeprob, punknowns)(initializeprob) + for (p, pval) in zip(punknowns, pvals) p = unwrap(p) - op[p] = getu(initializeprob, p)(initializeprob) + op[p] = pval if iscall(p) && operation(p) === getindex arrp = arguments(p)[1] + get(op, arrp, nothing) !== missing && continue op[arrp] = collect(arrp) end end if is_time_dependent(sys) - for v in missing_unknowns - op[v] = getu(initializeprob, v)(initializeprob) + uvals = getu(initializeprob, collect(missing_unknowns))(initializeprob) + for (v, val) in zip(missing_unknowns, uvals) + op[v] = val end empty!(missing_unknowns) end From 592f3f3ac6d7d8ceeacee1240cf2467672af0d03 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Thu, 19 Jun 2025 15:40:49 +0530 Subject: [PATCH 27/42] refactor: add fast path in `build_operating_point!` --- src/systems/problem_utils.jl | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/src/systems/problem_utils.jl b/src/systems/problem_utils.jl index 377b3bdc7c..38b38d89af 100644 --- a/src/systems/problem_utils.jl +++ b/src/systems/problem_utils.jl @@ -583,17 +583,19 @@ function build_operating_point!(sys::AbstractSystem, end end - for (k, v) in u0map - symbolic_type(v) == NotSymbolic() && !is_array_of_symbolics(v) && continue - v = fixpoint_sub(v, neithermap; operator = Symbolics.Operator) - isequal(k, v) && continue - u0map[k] = v - end - for (k, v) in pmap - symbolic_type(v) == NotSymbolic() && !is_array_of_symbolics(v) && continue - v = fixpoint_sub(v, neithermap; operator = Symbolics.Operator) - isequal(k, v) && continue - pmap[k] = v + if !isempty(neithermap) + for (k, v) in u0map + symbolic_type(v) == NotSymbolic() && !is_array_of_symbolics(v) && continue + v = fixpoint_sub(v, neithermap; operator = Symbolics.Operator) + isequal(k, v) && continue + u0map[k] = v + end + for (k, v) in pmap + symbolic_type(v) == NotSymbolic() && !is_array_of_symbolics(v) && continue + v = fixpoint_sub(v, neithermap; operator = Symbolics.Operator) + isequal(k, v) && continue + pmap[k] = v + end end return op, missing_unknowns, missing_pars From 2922aa7d5434835344adabdd8f0a0122123d968e Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Fri, 20 Jun 2025 14:45:14 +0530 Subject: [PATCH 28/42] fix: properly handle values given to parameter dependencies in `late_binding_update_u0_p` --- src/systems/nonlinear/initializesystem.jl | 60 ++++++++++++++++------- test/initializationsystem.jl | 4 +- 2 files changed, 43 insertions(+), 21 deletions(-) diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl index 5a8687984f..bd2c7de753 100644 --- a/src/systems/nonlinear/initializesystem.jl +++ b/src/systems/nonlinear/initializesystem.jl @@ -648,7 +648,26 @@ function SciMLBase.late_binding_update_u0_p( newu0, newp = promote_u0_p(newu0, newp, t0) # non-symbolic u0 updates initials... - if !(eltype(u0) <: Pair) + if eltype(u0) <: Pair + syms = [] + vals = [] + allsyms = all_symbols(sys) + for (k, v) in u0 + v === nothing && continue + (symbolic_type(v) == NotSymbolic() && !is_array_of_symbolics(v)) || continue + if k isa Symbol + k2 = symbol_to_symbolic(sys, k; allsyms) + # if it is returned as-is, there is no match so skip it + k2 === k && continue + k = k2 + end + is_parameter(sys, Initial(k)) || continue + push!(syms, Initial(k)) + push!(vals, v) + end + newp = setp_oop(sys, syms)(newp, vals) + else + allsyms = nothing # if `p` is not provided or is symbolic p === missing || eltype(p) <: Pair || return newu0, newp (newu0 === nothing || isempty(newu0)) && return newu0, newp @@ -661,27 +680,30 @@ function SciMLBase.late_binding_update_u0_p( throw(ArgumentError("Expected `newu0` to be of same length as unknowns ($(length(prob.u0))). Got $(typeof(newu0)) of length $(length(newu0))")) end newp = meta.set_initial_unknowns!(newp, newu0) - return newu0, newp - end - - syms = [] - vals = [] - allsyms = all_symbols(sys) - for (k, v) in u0 - v === nothing && continue - (symbolic_type(v) == NotSymbolic() && !is_array_of_symbolics(v)) || continue - if k isa Symbol - k2 = symbol_to_symbolic(sys, k; allsyms) - # if it is returned as-is, there is no match so skip it - k2 === k && continue - k = k2 + end + + if eltype(p) <: Pair + syms = [] + vals = [] + if allsyms === nothing + allsyms = all_symbols(sys) + end + for (k, v) in p + v === nothing && continue + (symbolic_type(v) == NotSymbolic() && !is_array_of_symbolics(v)) || continue + if k isa Symbol + k2 = symbol_to_symbolic(sys, k; allsyms) + # if it is returned as-is, there is no match so skip it + k2 === k && continue + k = k2 + end + is_parameter(sys, Initial(k)) || continue + push!(syms, Initial(k)) + push!(vals, v) end - is_parameter(sys, Initial(k)) || continue - push!(syms, Initial(k)) - push!(vals, v) + newp = setp_oop(sys, syms)(newp, vals) end - newp = setp_oop(sys, syms)(newp, vals) return newu0, newp end diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl index fd84c346f9..a7c6f4a27e 100644 --- a/test/initializationsystem.jl +++ b/test/initializationsystem.jl @@ -1253,9 +1253,9 @@ end @test init(prob3)[x] ≈ 1.0 prob4 = remake(prob; p = [p => 1.0]) test_dummy_initialization_equation(prob4, x) - prob5 = remake(prob; p = [p => missing, q => 2.0]) + prob5 = remake(prob; p = [p => missing, q => 4.0]) @test prob5.f.initialization_data !== nothing - @test init(prob5).ps[p] ≈ 1.0 + @test init(prob5).ps[p] ≈ 2.0 end @testset "Variables provided as symbols" begin From d37b2a8eca64883bc0f138314a16297303f84132 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Thu, 26 Jun 2025 11:41:01 +0530 Subject: [PATCH 29/42] fix: remove source of allocations in `InitializationProblem` --- src/systems/diffeqs/abstractodesystem.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index d802f49fee..f87c4d0071 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -1495,7 +1495,7 @@ function InitializationProblem{iip, specialize}(sys::AbstractSystem, @warn errmsg end - uninit = setdiff(unknowns(sys), [unknowns(isys); observables(isys)]) + uninit = setdiff(unknowns(sys), unknowns(isys), observables(isys)) # TODO: throw on uninitialized arrays filter!(x -> !(x isa Symbolics.Arr), uninit) From 39c62deb1e06d3d1a4478cdbd2f65a83dceb49c5 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Tue, 24 Jun 2025 23:26:37 +0530 Subject: [PATCH 30/42] feat: add the ability to completely remove vertices from `BipartiteGraph` --- src/bipartite_graph.jl | 30 ++++++++++++++++++++++++++++-- 1 file changed, 28 insertions(+), 2 deletions(-) diff --git a/src/bipartite_graph.jl b/src/bipartite_graph.jl index b6665646c9..8cdb76cca0 100644 --- a/src/bipartite_graph.jl +++ b/src/bipartite_graph.jl @@ -535,13 +535,39 @@ function set_neighbors!(g::BipartiteGraph, i::Integer, new_neighbors) end end -function delete_srcs!(g::BipartiteGraph, srcs) +function delete_srcs!(g::BipartiteGraph{I}, srcs; rm_verts = false) where {I} for s in srcs set_neighbors!(g, s, ()) end + if rm_verts + old_to_new_idxs = collect(one(I):I(nsrcs(g))) + for s in srcs + old_to_new_idxs[s] = zero(I) + end + offset = zero(I) + for i in eachindex(old_to_new_idxs) + if iszero(old_to_new_idxs[i]) + offset += one(I) + continue + end + old_to_new_idxs[i] -= offset + end + + if g.badjlist isa AbstractVector + for i in 1:ndsts(g) + for j in eachindex(g.badjlist[i]) + g.badjlist[i][j] = old_to_new_idxs[g.badjlist[i][j]] + end + filter!(!iszero, g.badjlist[i]) + end + end + deleteat!(g.fadjlist, srcs) + end g end -delete_dsts!(g::BipartiteGraph, srcs) = delete_srcs!(invview(g), srcs) +function delete_dsts!(g::BipartiteGraph, srcs; rm_verts = false) + delete_srcs!(invview(g), srcs; rm_verts) +end ### ### Edges iteration From 8216332cb80007192aba68abccacebb93f677bd5 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Thu, 19 Jun 2025 11:51:25 +0530 Subject: [PATCH 31/42] feat: preemptively tear some trivial equations in `mtkcompile` --- .../symbolics_tearing.jl | 3 +- src/systems/systemstructure.jl | 116 +++++++++++++++++- 2 files changed, 117 insertions(+), 2 deletions(-) diff --git a/src/structural_transformation/symbolics_tearing.jl b/src/structural_transformation/symbolics_tearing.jl index cc582cd473..244a737705 100644 --- a/src/structural_transformation/symbolics_tearing.jl +++ b/src/structural_transformation/symbolics_tearing.jl @@ -740,7 +740,8 @@ function update_simplified_system!( obs_sub[eq.lhs] = eq.rhs end # TODO: compute the dependency correctly so that we don't have to do this - obs = [fast_substitute(observed(sys), obs_sub); solved_eqs] + obs = [fast_substitute(observed(sys), obs_sub); solved_eqs; + fast_substitute(state.additional_observed, obs_sub)] unknowns = Any[v for (i, v) in enumerate(state.fullvars) diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl index e0feb0d34d..48214e01a4 100644 --- a/src/systems/systemstructure.jl +++ b/src/systems/systemstructure.jl @@ -208,12 +208,19 @@ mutable struct TearingState{T <: AbstractSystem} <: AbstractTearingState{T} structure::SystemStructure extra_eqs::Vector param_derivative_map::Dict{BasicSymbolic, Any} + original_eqs::Vector{Equation} + """ + Additional user-provided observed equations. The variables calculated here + are not used in the rest of the system. + """ + additional_observed::Vector{Equation} end TransformationState(sys::AbstractSystem) = TearingState(sys) function system_subset(ts::TearingState, ieqs::Vector{Int}) eqs = equations(ts) @set! ts.sys.eqs = eqs[ieqs] + @set! ts.original_eqs = ts.original_eqs[ieqs] @set! ts.structure = system_subset(ts.structure, ieqs) ts end @@ -266,6 +273,7 @@ function TearingState(sys; quick_cancel = false, check = true, sort_eqs = true) iv = length(ivs) == 1 ? ivs[1] : nothing # scalarize array equations, without scalarizing arguments to registered functions eqs = flatten_equations(copy(equations(sys))) + original_eqs = copy(eqs) neqs = length(eqs) dervaridxs = OrderedSet{Int}() var2idx = Dict{Any, Int}() @@ -378,6 +386,7 @@ function TearingState(sys; quick_cancel = false, check = true, sort_eqs = true) end end eqs = eqs[eqs_to_retain] + original_eqs = original_eqs[eqs_to_retain] neqs = length(eqs) symbolic_incidence = symbolic_incidence[eqs_to_retain] @@ -386,6 +395,7 @@ function TearingState(sys; quick_cancel = false, check = true, sort_eqs = true) # depending on order due to NP-completeness of tearing. sortidxs = Base.sortperm(eqs, by = string) eqs = eqs[sortidxs] + original_eqs = original_eqs[sortidxs] symbolic_incidence = symbolic_incidence[sortidxs] end @@ -475,13 +485,116 @@ function TearingState(sys; quick_cancel = false, check = true, sort_eqs = true) ts = TearingState(sys, fullvars, SystemStructure(complete(var_to_diff), complete(eq_to_diff), complete(graph), nothing, var_types, sys isa AbstractDiscreteSystem), - Any[], param_derivative_map) + Any[], param_derivative_map, original_eqs, Equation[]) if sys isa DiscreteSystem ts = shift_discrete_system(ts) end return ts end +""" + $(TYPEDSIGNATURES) + +Preemptively identify observed equations in the system and tear them. This happens before +any simplification. The equations torn by this process are ones that are already given in +an explicit form in the system and where the LHS is not present in any other equation of +the system except for other such preempitvely torn equations. +""" +function trivial_tearing!(ts::TearingState) + @assert length(ts.original_eqs) == length(equations(ts)) + # equations that can be trivially torn an observed equations + trivial_idxs = BitSet() + # equations to never check + blacklist = BitSet() + torn_eqs = Equation[] + # variables that have been matched to trivially torn equations + matched_vars = BitSet() + # variable to index in fullvars + var_to_idx = Dict{Any, Int}(ts.fullvars .=> eachindex(ts.fullvars)) + + complete!(ts.structure) + var_to_diff = ts.structure.var_to_diff + graph = ts.structure.graph + while true + # track whether we added an equation to the trivial list this iteration + added_equation = false + for (i, eq) in enumerate(ts.original_eqs) + # don't check already torn equations + i in trivial_idxs && continue + i in blacklist && continue + # ensure it is an observed equation matched to a variable in fullvars + vari = get(var_to_idx, eq.lhs, 0) + iszero(vari) && continue + # don't tear irreducible variables + if isirreducible(eq.lhs) + push!(blacklist, i) + continue + end + # if a variable was the LHS of two trivial observed equations, we wouldn't have + # included it in the list. Error if somehow it made it through. + @assert !(vari in matched_vars) + # don't tear differential/shift equations (or differentiated/shifted variables) + var_to_diff[vari] === nothing || continue + invview(var_to_diff)[vari] === nothing || continue + # get the equations that the candidate matched variable is present in, except + # those equations which have already been torn as observed + eqidxs = setdiff(𝑑neighbors(graph, vari), trivial_idxs) + # it should only be present in this equation + length(eqidxs) == 1 || continue + eqi = only(eqidxs) + @assert eqi == i + + # for every variable present in this equation, make sure it isn't _only_ + # present in trivial equations + isvalid = true + for v in 𝑠neighbors(graph, eqi) + v == vari && continue + v in matched_vars && continue + # `> 1` and not `0` because one entry will be this equation (`eqi`) + isvalid &= count(!in(trivial_idxs), 𝑑neighbors(graph, v)) > 1 + isvalid || break + end + isvalid || continue + # skip if the LHS is present in the RHS, since then this isn't explicit + if occursin(eq.lhs, eq.rhs) + push!(blacklist, i) + continue + end + + added_equation = true + push!(trivial_idxs, eqi) + push!(torn_eqs, eq) + push!(matched_vars, vari) + end + + # if we didn't add an equation this iteration, we won't add one next iteration + added_equation || break + end + + deleteat!(var_to_diff.primal_to_diff, matched_vars) + deleteat!(var_to_diff.diff_to_primal, matched_vars) + deleteat!(ts.structure.eq_to_diff.primal_to_diff, trivial_idxs) + deleteat!(ts.structure.eq_to_diff.diff_to_primal, trivial_idxs) + delete_srcs!(ts.structure.graph, trivial_idxs; rm_verts = true) + delete_dsts!(ts.structure.graph, matched_vars; rm_verts = true) + if ts.structure.solvable_graph !== nothing + delete_srcs!(ts.structure.solvable_graph, trivial_idxs; rm_verts = true) + delete_dsts!(ts.structure.solvable_graph, matched_vars; rm_verts = true) + end + if ts.structure.var_types !== nothing + deleteat!(ts.structure.var_types, matched_vars) + end + deleteat!(ts.fullvars, matched_vars) + deleteat!(ts.original_eqs, trivial_idxs) + ts.additional_observed = torn_eqs + sys = ts.sys + eqs = copy(get_eqs(sys)) + deleteat!(eqs, trivial_idxs) + @set! sys.eqs = eqs + ts.sys = sys + return ts +end + function lower_order_var(dervar, t) if isdifferential(dervar) diffvar = arguments(dervar)[1] @@ -739,6 +852,7 @@ function _structural_simplify!(state::TearingState, io; simplify = false, else input_idxs = 0:-1 # Empty range end + trivial_tearing!(state) sys, mm = ModelingToolkit.alias_elimination!(state; kwargs...) if check_consistency fully_determined = ModelingToolkit.check_consistency( From 396d9868ce451b74429629c695ba633ef116b9ce Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Wed, 25 Jun 2025 11:58:07 +0530 Subject: [PATCH 32/42] fix: ensure `initializeprobpmap` returns floats --- src/systems/problem_utils.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/systems/problem_utils.jl b/src/systems/problem_utils.jl index 38b38d89af..e4d043d334 100644 --- a/src/systems/problem_utils.jl +++ b/src/systems/problem_utils.jl @@ -1043,6 +1043,9 @@ function (siu::SetInitialUnknowns)(p::AbstractVector, u0) return p end +safe_float(x) = x +safe_float(x::AbstractArray) = isempty(x) ? x : float(x) + """ $(TYPEDSIGNATURES) @@ -1107,7 +1110,8 @@ function maybe_build_initialization_problem( if is_time_dependent(sys) all_init_syms = Set(all_symbols(initializeprob)) solved_unknowns = filter(var -> var in all_init_syms, unknowns(sys)) - initializeprobmap = u0_constructor ∘ getu(initializeprob, solved_unknowns) + initializeprobmap = u0_constructor ∘ safe_float ∘ + getu(initializeprob, solved_unknowns) else initializeprobmap = nothing end From a2deaa4c12d7842518174f5b6c408e6804a91a28 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Thu, 26 Jun 2025 11:44:51 +0530 Subject: [PATCH 33/42] fix: fix potential infinite recursion in simplification of `OptimizationSystem` --- src/systems/optimization/optimizationsystem.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/systems/optimization/optimizationsystem.jl b/src/systems/optimization/optimizationsystem.jl index bfe15b62d7..fcb502efdf 100644 --- a/src/systems/optimization/optimizationsystem.jl +++ b/src/systems/optimization/optimizationsystem.jl @@ -741,8 +741,9 @@ function structural_simplify(sys::OptimizationSystem; split = true, kwargs...) nlsys = NonlinearSystem(econs, unknowns(sys), parameters(sys); name = :___tmp_nlsystem) snlsys = structural_simplify(nlsys; fully_determined = false, kwargs...) obs = observed(snlsys) - subs = Dict(eq.lhs => eq.rhs for eq in observed(snlsys)) seqs = equations(snlsys) + trueobs, _ = unhack_observed(obs, seqs) + subs = Dict(eq.lhs => eq.rhs for eq in trueobs) cons_simplified = similar(cons, length(icons) + length(seqs)) for (i, eq) in enumerate(Iterators.flatten((seqs, icons))) cons_simplified[i] = fixpoint_sub(eq, subs) From f2ddf79451087605b33747a073c2ebacb1585b93 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Wed, 25 Jun 2025 23:22:01 +0530 Subject: [PATCH 34/42] fix: retain events and metadata in `substitute` --- src/systems/abstractsystem.jl | 19 ++++++++++++------- test/odesystem.jl | 23 +++++++++++++++++++++++ 2 files changed, 35 insertions(+), 7 deletions(-) diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl index 6c00806b7d..d78c48a857 100644 --- a/src/systems/abstractsystem.jl +++ b/src/systems/abstractsystem.jl @@ -2835,15 +2835,20 @@ function Symbolics.substitute(sys::AbstractSystem, rules::Union{Vector{<:Pair}, elseif sys isa ODESystem rules = todict(map(r -> Symbolics.unwrap(r[1]) => Symbolics.unwrap(r[2]), collect(rules))) - eqs = fast_substitute(get_eqs(sys), rules) - pdeps = fast_substitute(get_parameter_dependencies(sys), rules) - defs = Dict(fast_substitute(k, rules) => fast_substitute(v, rules) + newsys = @set sys.eqs = fast_substitute(get_eqs(sys), rules) + @set! newsys.unknowns = map(get_unknowns(sys)) do var + get(rules, var, var) + end + @set! newsys.ps = map(get_ps(sys)) do var + get(rules, var, var) + end + @set! newsys.parameter_dependencies = fast_substitute( + get_parameter_dependencies(sys), rules) + @set! newsys.defaults = Dict(fast_substitute(k, rules) => fast_substitute(v, rules) for (k, v) in get_defaults(sys)) - guess = Dict(fast_substitute(k, rules) => fast_substitute(v, rules) + @set! newsys.guesses = Dict(fast_substitute(k, rules) => fast_substitute(v, rules) for (k, v) in get_guesses(sys)) - subsys = map(s -> substitute(s, rules), get_systems(sys)) - ODESystem(eqs, get_iv(sys); name = nameof(sys), defaults = defs, - guesses = guess, parameter_dependencies = pdeps, systems = subsys) + @set! newsys.systems = map(s -> substitute(s, rules), get_systems(sys)) else error("substituting symbols is not supported for $(typeof(sys))") end diff --git a/test/odesystem.jl b/test/odesystem.jl index 6065e9e70f..5cdce4f2eb 100644 --- a/test/odesystem.jl +++ b/test/odesystem.jl @@ -1783,3 +1783,26 @@ end @test !process_running(proc) kill(proc, Base.SIGKILL) end + +@testset "`substitute` retains events and metadata" begin + @parameters p(t) = 1.0 + @variables x(t) = 0.0 + event = [0.5] => [p ~ t] + event2 = [x ~ 0.75] => [p ~ 2 * t] + + eq = [ + D(x) ~ p + ] + @named sys = ODESystem(eq, t, [x], [p], discrete_events = [event], + continuous_events = [event2], metadata = "TEST") + + @variables x2(t) = 0.0 + sys2 = substitute(sys, [x => x2]) + + @test length(ModelingToolkit.get_discrete_events(sys)) == 1 + @test length(ModelingToolkit.get_discrete_events(sys2)) == 1 + @test length(ModelingToolkit.get_continuous_events(sys)) == 1 + @test length(ModelingToolkit.get_continuous_events(sys2)) == 1 + @test ModelingToolkit.get_metadata(sys) == "TEST" + @test ModelingToolkit.get_metadata(sys2) == "TEST" +end From 5b6aff223e9bf9ab2589d7f820b29d45891348a2 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Thu, 3 Jul 2025 20:42:06 +0530 Subject: [PATCH 35/42] build: bump minor version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 130e0541f3..e196aa7ab0 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ModelingToolkit" uuid = "961ee093-0014-501f-94e3-6117800e7a78" authors = ["Yingbo Ma ", "Chris Rackauckas and contributors"] -version = "9.80.6" +version = "9.81.0" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" From 9f68ac7a1f7462379dbb4d658fb8c9fe92b00c17 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Fri, 4 Jul 2025 17:13:54 +0530 Subject: [PATCH 36/42] build: bump patch version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index e196aa7ab0..b5c874fd54 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ModelingToolkit" uuid = "961ee093-0014-501f-94e3-6117800e7a78" authors = ["Yingbo Ma ", "Chris Rackauckas and contributors"] -version = "9.81.0" +version = "9.81.1" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" From f9ccc77eed9a444092356ac891558514cf4d7abc Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Wed, 9 Jul 2025 11:25:53 +0530 Subject: [PATCH 37/42] feat: allow using `System` as a type --- src/systems/systems.jl | 2 ++ test/odesystem.jl | 6 ++++++ 2 files changed, 8 insertions(+) diff --git a/src/systems/systems.jl b/src/systems/systems.jl index 52f93afb9b..aa7ccce917 100644 --- a/src/systems/systems.jl +++ b/src/systems/systems.jl @@ -1,3 +1,5 @@ +const System = AbstractODESystem + function System(eqs::AbstractVector{<:Equation}, iv, args...; name = nothing, kw...) ODESystem(eqs, iv, args...; name, kw..., checks = false) diff --git a/test/odesystem.jl b/test/odesystem.jl index 5cdce4f2eb..dcda766177 100644 --- a/test/odesystem.jl +++ b/test/odesystem.jl @@ -1806,3 +1806,9 @@ end @test ModelingToolkit.get_metadata(sys) == "TEST" @test ModelingToolkit.get_metadata(sys2) == "TEST" end + +@testset "`System` works as a type" begin + @variables x(t) + @named sys = System([D(x) ~ 2x], t) + @test sys isa System +end From 206db0b783cab9209d96c23f18e75d064a3b1014 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Wed, 9 Jul 2025 19:01:01 +0530 Subject: [PATCH 38/42] build: bump minor version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index b5c874fd54..c04793d1e2 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ModelingToolkit" uuid = "961ee093-0014-501f-94e3-6117800e7a78" authors = ["Yingbo Ma ", "Chris Rackauckas and contributors"] -version = "9.81.1" +version = "9.82.0" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" From c46322df5fe771c1f9dc1f6f7c7f03cb8146575b Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 26 Jul 2025 21:55:11 -0400 Subject: [PATCH 39/42] Merge pull request #3843 from oscardssmith/os/ForwardDiff@1 update ForwardDiff to 1 (cherry picked from commit c2281f8c359503f54bef2d2ab5694c1a76a49894) --- Project.toml | 4 ++-- test/extensions/Project.toml | 3 +++ 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index c04793d1e2..3e84228112 100644 --- a/Project.toml +++ b/Project.toml @@ -111,7 +111,7 @@ EnumX = "1.0.4" ExprTools = "0.1.10" FMI = "0.14" FindFirstFunctions = "1" -ForwardDiff = "0.10.3" +ForwardDiff = "0.10.3, 1" FunctionWrappers = "1.1" FunctionWrappersWrappers = "0.1" Graphs = "1.5.2" @@ -147,7 +147,7 @@ Serialization = "1" Setfield = "0.7, 0.8, 1" SimpleNonlinearSolve = "0.1.0, 1, 2" SparseArrays = "1" -SpecialFunctions = "0.7, 0.8, 0.9, 0.10, 1.0, 2" +SpecialFunctions = "1, 2" StaticArrays = "0.10, 0.11, 0.12, 1.0" StochasticDelayDiffEq = "1.10" StochasticDiffEq = "6.72.1" diff --git a/test/extensions/Project.toml b/test/extensions/Project.toml index f1647aef1a..f91231fe05 100644 --- a/test/extensions/Project.toml +++ b/test/extensions/Project.toml @@ -1,3 +1,6 @@ +[sources] +ModelingToolkit = {path = "../.."} + [deps] BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665" CasADi = "c49709b8-5c63-11e9-2fb2-69db5844192f" From 591e486f83ff72c4595b3b7479a244c2eb9dcd84 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Tue, 7 Oct 2025 11:28:05 +0200 Subject: [PATCH 40/42] Update version number to 9.83.0 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 3e84228112..6409dca9dc 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ModelingToolkit" uuid = "961ee093-0014-501f-94e3-6117800e7a78" authors = ["Yingbo Ma ", "Chris Rackauckas and contributors"] -version = "9.82.0" +version = "9.83.0" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" From 5b9e7abb210aba9a995cfb13475de97a0d8bc537 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 2 Jul 2025 12:53:28 +0000 Subject: [PATCH 41/42] Merge pull request #3733 from oscardssmith/os/bump-di fix linearization t0 and bump di (cherry picked from commit d52b73540c1dab6550f3958b489880b860e389e4) --- Project.toml | 2 +- src/linearization.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index 6409dca9dc..981226aa15 100644 --- a/Project.toml +++ b/Project.toml @@ -101,7 +101,7 @@ DiffEqBase = "6.170.1" DiffEqCallbacks = "2.16, 3, 4" DiffEqNoiseProcess = "5" DiffRules = "0.1, 1.0" -DifferentiationInterface = "0.6.47" +DifferentiationInterface = "0.6.47, 0.7" Distributed = "1" Distributions = "0.23, 0.24, 0.25" DocStringExtensions = "0.7, 0.8, 0.9" diff --git a/src/linearization.jl b/src/linearization.jl index b30d275818..fec5711dde 100644 --- a/src/linearization.jl +++ b/src/linearization.jl @@ -163,13 +163,13 @@ struct PreparedJacobian{iip, P, F, B, A} end function PreparedJacobian{true}(f, buf, autodiff, args...) - prep = DI.prepare_jacobian(f, buf, autodiff, args...) + prep = DI.prepare_jacobian(f, buf, autodiff, args...; strict = Val(false)) return PreparedJacobian{true, typeof(prep), typeof(f), typeof(buf), typeof(autodiff)}( prep, f, buf, autodiff) end function PreparedJacobian{false}(f, autodiff, args...) - prep = DI.prepare_jacobian(f, autodiff, args...) + prep = DI.prepare_jacobian(f, autodiff, args...; strict = Val(false)) return PreparedJacobian{true, typeof(prep), typeof(f), Nothing, typeof(autodiff)}( prep, f, nothing) end From 62dfa1f1254121bcaf66e99d448ecb198de5b905 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Thu, 6 Nov 2025 09:10:02 +0100 Subject: [PATCH 42/42] Bump version to 9.84.0 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 981226aa15..d68990f6e4 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ModelingToolkit" uuid = "961ee093-0014-501f-94e3-6117800e7a78" authors = ["Yingbo Ma ", "Chris Rackauckas and contributors"] -version = "9.83.0" +version = "9.84.0" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"