From 9b4a18c902efa382d095512d93051e8f57b5d438 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Wed, 2 Apr 2025 16:16:50 +0530
Subject: [PATCH 1/6] docs: document `SymScope` variants
---
src/systems/abstractsystem.jl | 73 +++++++++++++++++++++++++++++++++++
1 file changed, 73 insertions(+)
diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index 1987cb9a60..9bbe05dbeb 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -1114,9 +1114,25 @@ function _apply_to_variables(f::F, ex) where {F}
metadata(ex))
end
+"""
+Variable metadata key which contains information about scoping/namespacing of the
+variable in a hierarchical system.
+"""
abstract type SymScope end
+"""
+ $(TYPEDEF)
+
+The default scope of a variable. It belongs to the system whose equations it is involved
+in and is namespaced by every level of the hierarchy.
+"""
struct LocalScope <: SymScope end
+
+"""
+ $(TYPEDSIGNATURES)
+
+Apply `LocalScope` to `sym`.
+"""
function LocalScope(sym::Union{Num, Symbolic, Symbolics.Arr{Num}})
apply_to_variables(sym) do sym
if iscall(sym) && operation(sym) === getindex
@@ -1130,9 +1146,25 @@ function LocalScope(sym::Union{Num, Symbolic, Symbolics.Arr{Num}})
end
end
+"""
+ $(TYPEDEF)
+
+Denotes that the variable does not belong to the system whose equations it is involved
+in. It is not namespaced by this system. In the immediate parent of this system, the
+scope of this variable is given by `parent`.
+
+# Fields
+
+$(TYPEDFIELDS)
+"""
struct ParentScope <: SymScope
parent::SymScope
end
+"""
+ $(TYPEDSIGNATURES)
+
+Apply `ParentScope` to `sym`, with `parent` being `LocalScope`.
+"""
function ParentScope(sym::Union{Num, Symbolic, Symbolics.Arr{Num}})
apply_to_variables(sym) do sym
if iscall(sym) && operation(sym) === getindex
@@ -1148,10 +1180,31 @@ function ParentScope(sym::Union{Num, Symbolic, Symbolics.Arr{Num}})
end
end
+"""
+ $(TYPEDEF)
+
+Denotes that a variable belongs to a system that is at least `N + 1` levels up in the
+hierarchy from the system whose equations it is involved in. It is namespaced by the
+first `N` parents and not namespaced by the `N+1`th parent in the hierarchy. The scope
+of the variable after this point is given by `parent`.
+
+In other words, this scope delays applying `ParentScope` by `N` levels, and applies
+`LocalScope` in the meantime.
+
+# Fields
+
+$(TYPEDFIELDS)
+"""
struct DelayParentScope <: SymScope
parent::SymScope
N::Int
end
+
+"""
+ $(TYPEDSIGNATURES)
+
+Apply `DelayParentScope` to `sym`, with a delay of `N` and `parent` being `LocalScope`.
+"""
function DelayParentScope(sym::Union{Num, Symbolic, Symbolics.Arr{Num}}, N)
apply_to_variables(sym) do sym
if iscall(sym) && operation(sym) == getindex
@@ -1166,9 +1219,29 @@ function DelayParentScope(sym::Union{Num, Symbolic, Symbolics.Arr{Num}}, N)
end
end
end
+
+"""
+ $(TYPEDSIGNATURES)
+
+Apply `DelayParentScope` to `sym`, with a delay of `1` and `parent` being `LocalScope`.
+"""
DelayParentScope(sym::Union{Num, Symbolic, Symbolics.Arr{Num}}) = DelayParentScope(sym, 1)
+"""
+ $(TYPEDEF)
+
+Denotes that a variable belongs to the root system in the hierarchy, regardless of which
+equations of subsystems in the hierarchy it is involved in. Variables with this scope
+are never namespaced and only added to the unknowns/parameters of a system when calling
+`complete` or `structural_simplify`.
+"""
struct GlobalScope <: SymScope end
+
+"""
+ $(TYPEDSIGNATURES)
+
+Apply `GlobalScope` to `sym`.
+"""
function GlobalScope(sym::Union{Num, Symbolic, Symbolics.Arr{Num}})
apply_to_variables(sym) do sym
if iscall(sym) && operation(sym) == getindex
From 0166a34bdaff285a6a00d1f69efa8d2ce6c65797 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Wed, 2 Apr 2025 16:17:29 +0530
Subject: [PATCH 2/6] fix: fix incorrect namespacing of `DelayParentScope`
variables in `collect_scoped_vars!`
---
src/utils.jl | 15 ++++++---------
1 file changed, 6 insertions(+), 9 deletions(-)
diff --git a/src/utils.jl b/src/utils.jl
index 2a0009b644..1884a91c19 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -535,7 +535,7 @@ recursively searches through all subsystems of `sys`, increasing the depth if it
"""
function collect_scoped_vars!(unknowns, parameters, sys, iv; depth = 1, op = Differential)
if has_eqs(sys)
- for eq in get_eqs(sys)
+ for eq in equations(sys)
eqtype_supports_collect_vars(eq) || continue
if eq isa Equation
eq.lhs isa Union{Symbolic, Number} || continue
@@ -544,7 +544,7 @@ function collect_scoped_vars!(unknowns, parameters, sys, iv; depth = 1, op = Dif
end
end
if has_parameter_dependencies(sys)
- for eq in get_parameter_dependencies(sys)
+ for eq in parameter_dependencies(sys)
if eq isa Pair
collect_vars!(unknowns, parameters, eq, iv; depth, op)
else
@@ -553,17 +553,13 @@ function collect_scoped_vars!(unknowns, parameters, sys, iv; depth = 1, op = Dif
end
end
if has_constraints(sys)
- for eq in get_constraints(sys)
+ for eq in constraints(sys)
eqtype_supports_collect_vars(eq) || continue
collect_vars!(unknowns, parameters, eq, iv; depth, op)
end
end
if has_op(sys)
- collect_vars!(unknowns, parameters, get_op(sys), iv; depth, op)
- end
- newdepth = depth == -1 ? depth : depth + 1
- for ssys in get_systems(sys)
- collect_scoped_vars!(unknowns, parameters, ssys, iv; depth = newdepth, op)
+ collect_vars!(unknowns, parameters, objective(sys), iv; depth, op)
end
end
@@ -608,6 +604,7 @@ end
function collect_var!(unknowns, parameters, var, iv; depth = 0)
isequal(var, iv) && return nothing
check_scope_depth(getmetadata(var, SymScope, LocalScope()), depth) || return nothing
+ var = setmetadata(var, SymScope, LocalScope())
if iscalledparameter(var)
callable = getcalledparameter(var)
push!(parameters, callable)
@@ -636,7 +633,7 @@ function check_scope_depth(scope, depth)
elseif scope isa ParentScope
return depth > 0 && check_scope_depth(scope.parent, depth - 1)
elseif scope isa DelayParentScope
- return depth >= scope.N && check_scope_depth(scope.parent, depth - scope.N)
+ return depth >= scope.N && check_scope_depth(scope.parent, depth - scope.N - 1)
elseif scope isa GlobalScope
return depth == -1
end
From 439e8623fe1f43aa96fbb28726cfefe0baf36134 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Wed, 2 Apr 2025 16:17:46 +0530
Subject: [PATCH 3/6] test: test proper scoping of `DelayParentScope` during
variable discovery
---
test/variable_scope.jl | 8 ++++----
1 file changed, 4 insertions(+), 4 deletions(-)
diff --git a/test/variable_scope.jl b/test/variable_scope.jl
index e90f7e6fbf..bd1d3cb0cf 100644
--- a/test/variable_scope.jl
+++ b/test/variable_scope.jl
@@ -106,12 +106,12 @@ defs = ModelingToolkit.defaults(bar)
@variables x1(t) x2(t) x3(t) x4(t) x5(t)
x2 = ParentScope(x2)
x3 = ParentScope(ParentScope(x3))
-x4 = DelayParentScope(x4, 2)
+x4 = DelayParentScope(x4)
x5 = GlobalScope(x5)
@parameters p1 p2 p3 p4 p5
p2 = ParentScope(p2)
p3 = ParentScope(ParentScope(p3))
-p4 = DelayParentScope(p4, 2)
+p4 = DelayParentScope(p4)
p5 = GlobalScope(p5)
@named sys1 = ODESystem([D(x1) ~ p1, D(x2) ~ p2, D(x3) ~ p3, D(x4) ~ p4, D(x5) ~ p5], t)
@@ -126,10 +126,10 @@ p5 = GlobalScope(p5)
sys3 = sys3 ∘ sys2
@test length(unknowns(sys3)) == 4
@test any(isequal(x3), unknowns(sys3))
-@test any(isequal(x4), unknowns(sys3))
+@test any(isequal(ModelingToolkit.renamespace(sys1, x4)), unknowns(sys3))
@test length(parameters(sys3)) == 4
@test any(isequal(p3), parameters(sys3))
-@test any(isequal(p4), parameters(sys3))
+@test any(isequal(ModelingToolkit.renamespace(sys1, p4)), parameters(sys3))
sys4 = complete(sys3)
@test length(unknowns(sys3)) == 4
@test length(parameters(sys4)) == 5
From 441ab0205a592c5d7f1e4c83494f67cc85ddb729 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Thu, 3 Apr 2025 16:17:44 +0530
Subject: [PATCH 4/6] feat: add `expected_scope_depth`
---
src/utils.jl | 17 +++++++++++++++++
1 file changed, 17 insertions(+)
diff --git a/src/utils.jl b/src/utils.jl
index 1884a91c19..ca62316c2e 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -639,6 +639,23 @@ function check_scope_depth(scope, depth)
end
end
+"""
+ $(TYPEDSIGNATURES)
+
+Return the expected depth of the given `SymScope` from the root system.
+"""
+function expected_scope_depth(scope)
+ if scope isa LocalScope
+ return 0
+ elseif scope isa ParentScope
+ return expected_scope_depth(scope.parent) + 1
+ elseif scope isa DelayParentScope
+ return expected_scope_depth(scope.parent) + scope.N + 1
+ elseif scope isa GlobalScope
+ return -1
+ end
+end
+
"""
Find all the symbolic constants of some equations or terms and return them as a vector.
"""
From 25795584d7fce0296446439eb5eb9c7197dd6500 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Fri, 4 Apr 2025 13:58:01 +0530
Subject: [PATCH 5/6] fix: don't `toparam` inside `Initial`
---
src/systems/abstractsystem.jl | 6 +++---
1 file changed, 3 insertions(+), 3 deletions(-)
diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index 9bbe05dbeb..6e1304260d 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -668,14 +668,14 @@ function (f::Initial)(x)
iscall(x) && operation(x) isa Initial && return x
result = if symbolic_type(x) == ArraySymbolic()
# create an array for `Initial(array)`
- Symbolics.array_term(f, toparam(x))
+ Symbolics.array_term(f, x)
elseif iscall(x) && operation(x) == getindex
# instead of `Initial(x[1])` create `Initial(x)[1]`
# which allows parameter indexing to handle this case automatically.
arr = arguments(x)[1]
- term(getindex, f(toparam(arr)), arguments(x)[2:end]...)
+ term(getindex, f(arr), arguments(x)[2:end]...)
else
- term(f, toparam(x))
+ term(f, x)
end
# the result should be a parameter
result = toparam(result)
From 2a2b62070d8946e8ff5d6c49ec3143fd7210e6e3 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Fri, 4 Apr 2025 13:58:22 +0530
Subject: [PATCH 6/6] feat: reduce reliance on metadata in
`structural_simplify`
---
src/systems/systemstructure.jl | 183 +++++++++++++++++++--------------
1 file changed, 104 insertions(+), 79 deletions(-)
diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl
index d27e5c93a1..d91f387ff7 100644
--- a/src/systems/systemstructure.jl
+++ b/src/systems/systemstructure.jl
@@ -253,105 +253,106 @@ function Base.push!(ev::EquationsView, eq)
push!(ev.ts.extra_eqs, eq)
end
+function symbolic_contains(var, set)
+ var in set || symbolic_type(var) == ArraySymbolic() && Symbolics.shape(var) != Symbolics.Unknown() && all(i -> var[i] in set, eachindex(var))
+end
+
function TearingState(sys; quick_cancel = false, check = true)
+ # flatten system
sys = flatten(sys)
ivs = independent_variables(sys)
iv = length(ivs) == 1 ? ivs[1] : nothing
- # scalarize array equations, without scalarizing arguments to registered functions
- eqs = flatten_equations(copy(equations(sys)))
+ # flatten array equations
+ eqs = flatten_equations(equations(sys))
neqs = length(eqs)
- dervaridxs = OrderedSet{Int}()
- var2idx = Dict{Any, Int}()
- symbolic_incidence = []
- fullvars = []
- var_counter = Ref(0)
- var_types = VariableType[]
- addvar! = let fullvars = fullvars, var_counter = var_counter, var_types = var_types
+ # * Scalarize unknowns
+ dvs = Set{BasicSymbolic}()
+ fullvars = BasicSymbolic[]
+ for x in unknowns(sys)
+ push!(dvs, x)
+ xx = Symbolics.scalarize(x)
+ if xx isa AbstractArray
+ union!(dvs, xx)
+ append!(fullvars, xx)
+ else
+ push!(fullvars, xx)
+ end
+ end
+ var2idx = Dict{BasicSymbolic, Int}(v => k for (k, v) in enumerate(fullvars))
+ addvar! = let fullvars = fullvars, dvs = dvs, var2idx = var2idx
var -> get!(var2idx, var) do
+ push!(dvs, var)
push!(fullvars, var)
- push!(var_types, getvariabletype(var))
- var_counter[] += 1
+ return length(fullvars)
end
end
- vars = OrderedSet()
- varsvec = []
- for (i, eq′) in enumerate(eqs)
- if eq′.lhs isa Connection
- check ? error("$(nameof(sys)) has unexpanded `connect` statements") :
- return nothing
- end
- if _iszero(eq′.lhs)
- rhs = quick_cancel ? quick_cancel_expr(eq′.rhs) : eq′.rhs
- eq = eq′
- else
- lhs = quick_cancel ? quick_cancel_expr(eq′.lhs) : eq′.lhs
- rhs = quick_cancel ? quick_cancel_expr(eq′.rhs) : eq′.rhs
- eq = 0 ~ rhs - lhs
+ # build symbolic incidence
+ symbolic_incidence = Vector{BasicSymbolic}[]
+ varsbuf = Set()
+ for (i, eq) in enumerate(eqs)
+ rhs = quick_cancel ? quick_cancel_expr(eq.rhs) : eq.rhs
+ if !_iszero(eq.lhs)
+ lhs = quick_cancel ? quick_cancel_expr(eq.lhs) : eq.lhs
+ eq = eqs[i] = 0 ~ rhs - lhs
end
- vars!(vars, eq.rhs, op = Symbolics.Operator)
- for v in vars
- _var, _ = var_from_nested_derivative(v)
- any(isequal(_var), ivs) && continue
- if isparameter(_var) ||
- (iscall(_var) && isparameter(operation(_var)) || isconstant(_var))
- continue
+ empty!(varsbuf)
+ vars!(varsbuf, eq; op = Symbolics.Operator)
+ incidence = Set{BasicSymbolic}()
+ for v in varsbuf
+ # FIXME: This check still needs to rely on metadata
+ isconstant(v) && continue
+ vtype = getvariabletype(v)
+ # additionally track brownians in fullvars
+ # TODO: When uniting system types, track brownians in their own field
+ if vtype == BROWNIAN
+ i = addvar!(v)
+ push!(incidence, v)
end
- v = scalarize(v)
- if v isa AbstractArray
- append!(varsvec, v)
- else
- push!(varsvec, v)
- end
- end
- isalgeq = true
- unknownvars = []
- for var in varsvec
- ModelingToolkit.isdelay(var, iv) && continue
- set_incidence = true
- @label ANOTHER_VAR
- _var, _ = var_from_nested_derivative(var)
- any(isequal(_var), ivs) && continue
- if isparameter(_var) ||
- (iscall(_var) && isparameter(operation(_var)) || isconstant(_var))
- continue
- end
- varidx = addvar!(var)
- set_incidence && push!(unknownvars, var)
-
- dvar = var
- idx = varidx
- while isdifferential(dvar)
- if !(idx in dervaridxs)
- push!(dervaridxs, idx)
+
+ vtype == VARIABLE || continue
+
+ if !symbolic_contains(v, dvs)
+ isvalid = iscall(v) && operation(v) isa Union{Shift, Sample, Hold}
+ v′ = v
+ while !isvalid && iscall(v′) && operation(v′) isa Union{Differential, Shift}
+ v′ = arguments(v)[1]
+ if v′ in dvs || getmetadata(v′, SymScope, LocalScope()) isa GlobalScope
+ isvalid = true
+ break
+ end
+ end
+ if !isvalid
+ throw(ArgumentError("$v is present in the system but $v′ is not an unknown."))
end
- isalgeq = false
- dvar = arguments(dvar)[1]
- idx = addvar!(dvar)
- end
- dvar = var
- idx = varidx
+ addvar!(v)
+ if iscall(v) && operation(v) isa Symbolics.Operator && !isdifferential(v) && (it = input_timedomain(v)) !== nothing
+ v′ = only(arguments(v))
+ addvar!(setmetadata(v′, VariableTimeDomain, it))
+ end
+ end
- if iscall(var) && operation(var) isa Symbolics.Operator &&
- !isdifferential(var) && (it = input_timedomain(var)) !== nothing
- set_incidence = false
- var = only(arguments(var))
- var = setmetadata(var, VariableTimeDomain, it)
- @goto ANOTHER_VAR
+ if symbolic_type(v) == ArraySymbolic()
+ union!(incidence, collect(v))
+ else
+ push!(incidence, v)
end
end
- push!(symbolic_incidence, copy(unknownvars))
- empty!(unknownvars)
- empty!(vars)
- empty!(varsvec)
- if isalgeq
- eqs[i] = eq
- else
- eqs[i] = eqs[i].lhs ~ rhs
+
+ push!(symbolic_incidence, collect(incidence))
+ end
+
+ dervaridxs = Int[]
+ for (i, v) in enumerate(fullvars)
+ while isdifferential(v)
+ push!(dervaridxs, i)
+ v = arguments(v)[1]
+ i = addvar!(v)
end
end
+ # Handle shifts - find lowest shift and add intermediates with derivative edges
### Handle discrete variables
lowest_shift = Dict()
for var in fullvars
@@ -391,6 +392,9 @@ function TearingState(sys; quick_cancel = false, check = true)
end
end
end
+
+ var_types = Vector{VariableType}(getvariabletype.(fullvars))
+
# sort `fullvars` such that the mass matrix is as diagonal as possible.
dervaridxs = collect(dervaridxs)
sorted_fullvars = OrderedSet(fullvars[dervaridxs])
@@ -414,6 +418,7 @@ function TearingState(sys; quick_cancel = false, check = true)
var2idx = Dict(fullvars .=> eachindex(fullvars))
dervaridxs = 1:length(dervaridxs)
+ # build `var_to_diff`
nvars = length(fullvars)
diffvars = []
var_to_diff = DiffGraph(nvars, true)
@@ -425,6 +430,7 @@ function TearingState(sys; quick_cancel = false, check = true)
var_to_diff[diffvaridx] = dervaridx
end
+ # build incidence graph
graph = BipartiteGraph(neqs, nvars, Val(false))
for (ie, vars) in enumerate(symbolic_incidence), v in vars
jv = var2idx[v]
@@ -432,6 +438,7 @@ function TearingState(sys; quick_cancel = false, check = true)
end
@set! sys.eqs = eqs
+ @set! sys.unknowns = [v for (i, v) in enumerate(fullvars) if var_types[i] != BROWNIAN]
eq_to_diff = DiffGraph(nsrcs(graph))
@@ -439,6 +446,8 @@ function TearingState(sys; quick_cancel = false, check = true)
SystemStructure(complete(var_to_diff), complete(eq_to_diff),
complete(graph), nothing, var_types, sys isa AbstractDiscreteSystem),
Any[])
+
+ # `shift_discrete_system`
if sys isa DiscreteSystem
ts = shift_discrete_system(ts)
end
@@ -726,3 +735,19 @@ function _structural_simplify!(state::TearingState, io; simplify = false,
ModelingToolkit.invalidate_cache!(sys), input_idxs
end
+
+struct DifferentiatedVariableNotUnknownError <: Exception
+ differentiated
+ undifferentiated
+end
+
+function Base.showerror(io::IO, err::DifferentiatedVariableNotUnknownError)
+ undiff = err.undifferentiated
+ diff = err.differentiated
+ print(io, "Variable $undiff occurs differentiated as $diff but is not an unknown of the system.")
+ scope = getmetadata(undiff, SymScope, LocalScope())
+ depth = expected_scope_depth(scope)
+ if depth > 0
+ print(io, "\nVariable $undiff expects $depth more levels in the hierarchy to be an unknown.")
+ end
+end