Merge pull request #2798 from SciML/lowerplus

Bye bye ₊, use var symbols with .

Author: ChrisRackauckas

Project.toml

@@ -8,6 +8,7 @@ AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
+ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"
 ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
@@ -67,6 +68,7 @@ ArrayInterface = "6, 7"
 BifurcationKit = "0.3"
 Combinatorics = "1"
 Compat = "3.42, 4"
+ComponentArrays = "0.15"
 ConstructionBase = "1"
 DataStructures = "0.17, 0.18"
 DeepDiffs = "1"

docs/src/basics/Composition.md

@@ -37,10 +37,10 @@ connected = compose(
 equations(connected)
 
 #4-element Vector{Equation}:
-# Differential(t)(decay1₊f(t)) ~ 0
-# decay2₊f(t) ~ decay1₊x(t)
-# Differential(t)(decay1₊x(t)) ~ decay1₊f(t) - (decay1₊a*(decay1₊x(t)))
-# Differential(t)(decay2₊x(t)) ~ decay2₊f(t) - (decay2₊a*(decay2₊x(t)))
+# Differential(t)(decay1.f(t)) ~ 0
+# decay2.f(t) ~ decay1.x(t)
+# Differential(t)(decay1.x(t)) ~ decay1.f(t) - (decay1.a*(decay1.x(t)))
+# Differential(t)(decay2.x(t)) ~ decay2.f(t) - (decay2.a*(decay2.x(t)))
 
 simplified_sys = structural_simplify(connected)
 
@@ -149,27 +149,27 @@ p = [a, b, c, d, e, f]
 level0 = ODESystem(Equation[], t, [], p; name = :level0)
 level1 = ODESystem(Equation[], t, [], []; name = :level1) ∘ level0
 parameters(level1)
-#level0₊a
+#level0.a
 #b
 #c
-#level0₊d
-#level0₊e
+#level0.d
+#level0.e
 #f
 level2 = ODESystem(Equation[], t, [], []; name = :level2) ∘ level1
 parameters(level2)
-#level1₊level0₊a
-#level1₊b
+#level1.level0.a
+#level1.b
 #c
-#level0₊d
-#level1₊level0₊e
+#level0.d
+#level1.level0.e
 #f
 level3 = ODESystem(Equation[], t, [], []; name = :level3) ∘ level2
 parameters(level3)
-#level2₊level1₊level0₊a
-#level2₊level1₊b
-#level2₊c
-#level2₊level0₊d
-#level1₊level0₊e
+#level2.level1.level0.a
+#level2.level1.b
+#level2.c
+#level2.level0.d
+#level1.level0.e
 #f
 ```
 

docs/src/basics/Events.md

@@ -172,10 +172,10 @@ When accessing variables of a sub-system, it can be useful to rename them
 (alternatively, an affect function may be reused in different contexts):
 
 ```julia
-[x ~ 0] => (affect!, [resistor₊v => :v, x], [p, q => :p2], [], ctx)
+[x ~ 0] => (affect!, [resistor.v => :v, x], [p, q => :p2], [], ctx)
 ```
 
-Here, the symbolic variable `resistor₊v` is passed as `v` while the symbolic
+Here, the symbolic variable `resistor.v` is passed as `v` while the symbolic
 parameter `q` has been renamed `p2`.
 
 As an example, here is the bouncing ball example from above using the functional

docs/src/tutorials/domain_connections.md

@@ -115,7 +115,7 @@ end
 nothing #hide
 ```
 
-To see how the domain works, we can examine the set parameter values for each of the ports `src.port` and `vol.port`.  First we assemble the system using `structural_simplify()` and then check the default value of `vol.port.ρ`, whichs points to the setter value `fluid₊ρ`.  Likewise, `src.port.ρ`, will also point to the setter value `fluid₊ρ`.  Therefore, there is now only 1 defined density value `fluid₊ρ` which sets the density for the connected network.
+To see how the domain works, we can examine the set parameter values for each of the ports `src.port` and `vol.port`.  First we assemble the system using `structural_simplify()` and then check the default value of `vol.port.ρ`, whichs points to the setter value `fluid.ρ`.  Likewise, `src.port.ρ`, will also point to the setter value `fluid.ρ`.  Therefore, there is now only 1 defined density value `fluid.ρ` which sets the density for the connected network.
 
 ```@repl domain
 sys = structural_simplify(odesys)
@@ -181,7 +181,7 @@ end
 nothing #hide
 ```
 
-After running `structural_simplify()` on `actsys2`, the defaults will show that `act.port_a.ρ` points to `fluid_a₊ρ` and `act.port_b.ρ` points to `fluid_b₊ρ`.  This is a special case, in most cases a hydraulic system will have only 1 fluid, however this simple system has 2 separate domain networks.  Therefore, we can connect a single fluid to both networks.  This does not interfere with the mathematical equations of the system, since no unknown variables are connected.
+After running `structural_simplify()` on `actsys2`, the defaults will show that `act.port_a.ρ` points to `fluid_a.ρ` and `act.port_b.ρ` points to `fluid_b.ρ`.  This is a special case, in most cases a hydraulic system will have only 1 fluid, however this simple system has 2 separate domain networks.  Therefore, we can connect a single fluid to both networks.  This does not interfere with the mathematical equations of the system, since no unknown variables are connected.
 
 ```@example domain
 @component function ActuatorSystem1(; name)

src/ModelingToolkit.jl

@@ -22,6 +22,7 @@ using DiffEqCallbacks
 using Graphs
 import ExprTools: splitdef, combinedef
 import OrderedCollections
+import ComponentArrays
 
 using SymbolicIndexingInterface
 using LinearAlgebra, SparseArrays, LabelledArrays

src/inputoutput.jl

@@ -71,11 +71,11 @@ function is_bound(sys, u, stack = [])
     In the following scenario
     julia> observed(syss)
         2-element Vector{Equation}:
-        sys₊y(tv) ~ sys₊x(tv)
-        y(tv) ~ sys₊x(tv)
-    sys₊y(t) is bound to the outer y(t) through the variable sys₊x(t) and should thus return is_bound(sys₊y(t)) = true.
-    When asking is_bound(sys₊y(t)), we know that we are looking through observed equations and can thus ask
-    if var is bound, if it is, then sys₊y(t) is also bound. This can lead to an infinite recursion, so we maintain a stack of variables we have previously asked about to be able to break cycles
+        sys.y(tv) ~ sys.x(tv)
+        y(tv) ~ sys.x(tv)
+    sys.y(t) is bound to the outer y(t) through the variable sys.x(t) and should thus return is_bound(sys.y(t)) = true.
+    When asking is_bound(sys.y(t)), we know that we are looking through observed equations and can thus ask
+    if var is bound, if it is, then sys.y(t) is also bound. This can lead to an infinite recursion, so we maintain a stack of variables we have previously asked about to be able to break cycles
     =#
     u ∈ Set(stack) && return false # Cycle detected
     eqs = equations(sys)
@@ -119,17 +119,17 @@ function same_or_inner_namespace(u, var)
     nv = get_namespace(var)
     nu == nv ||           # namespaces are the same
         startswith(nv, nu) || # or nv starts with nu, i.e., nv is an inner namespace to nu
-        occursin('₊', string(getname(var))) &&
-            !occursin('₊', string(getname(u))) # or u is top level but var is internal
+        occursin('.', string(getname(var))) &&
+            !occursin('.', string(getname(u))) # or u is top level but var is internal
 end
 
 function inner_namespace(u, var)
     nu = get_namespace(u)
     nv = get_namespace(var)
     nu == nv && return false
     startswith(nv, nu) || # or nv starts with nu, i.e., nv is an inner namespace to nu
-        occursin('₊', string(getname(var))) &&
-            !occursin('₊', string(getname(u))) # or u is top level but var is internal
+        occursin('.', string(getname(var))) &&
+            !occursin('.', string(getname(u))) # or u is top level but var is internal
 end
 
 """
@@ -139,11 +139,11 @@ Return the namespace of a variable as a string. If the variable is not namespace
 """
 function get_namespace(x)
     sname = string(getname(x))
-    parts = split(sname, '₊')
+    parts = split(sname, '.')
     if length(parts) == 1
         return ""
     end
-    join(parts[1:(end - 1)], '₊')
+    join(parts[1:(end - 1)], '.')
 end
 
 """

src/systems/abstractsystem.jl

@@ -365,8 +365,8 @@ function SymbolicIndexingInterface.is_variable(sys::AbstractSystem, sym::Symbol)
         return is_variable(ic, sym)
     end
     return any(isequal(sym), getname.(variable_symbols(sys))) ||
-           count('₊', string(sym)) == 1 &&
-           count(isequal(sym), Symbol.(nameof(sys), :₊, getname.(variable_symbols(sys)))) ==
+           count('.', string(sym)) == 1 &&
+           count(isequal(sym), Symbol.(nameof(sys), :., getname.(variable_symbols(sys)))) ==
            1
 end
 
@@ -399,9 +399,9 @@ function SymbolicIndexingInterface.variable_index(sys::AbstractSystem, sym::Symb
     idx = findfirst(isequal(sym), getname.(variable_symbols(sys)))
     if idx !== nothing
         return idx
-    elseif count('₊', string(sym)) == 1
+    elseif count('.', string(sym)) == 1
         return findfirst(isequal(sym),
-            Symbol.(nameof(sys), :₊, getname.(variable_symbols(sys))))
+            Symbol.(nameof(sys), :., getname.(variable_symbols(sys))))
     end
     return nothing
 end
@@ -431,9 +431,9 @@ function SymbolicIndexingInterface.is_parameter(sys::AbstractSystem, sym::Symbol
         return is_parameter(ic, sym)
     end
     return any(isequal(sym), getname.(parameter_symbols(sys))) ||
-           count('₊', string(sym)) == 1 &&
+           count('.', string(sym)) == 1 &&
            count(isequal(sym),
-        Symbol.(nameof(sys), :₊, getname.(parameter_symbols(sys)))) == 1
+        Symbol.(nameof(sys), :., getname.(parameter_symbols(sys)))) == 1
 end
 
 function SymbolicIndexingInterface.parameter_index(sys::AbstractSystem, sym)
@@ -466,9 +466,9 @@ function SymbolicIndexingInterface.parameter_index(sys::AbstractSystem, sym::Sym
     idx = findfirst(isequal(sym), getname.(parameter_symbols(sys)))
     if idx !== nothing
         return idx
-    elseif count('₊', string(sym)) == 1
+    elseif count('.', string(sym)) == 1
         return findfirst(isequal(sym),
-            Symbol.(nameof(sys), :₊, getname.(parameter_symbols(sys))))
+            Symbol.(nameof(sys), :., getname.(parameter_symbols(sys))))
     end
     return nothing
 end
@@ -889,7 +889,7 @@ function renamespace(sys, x)
     elseif x isa AbstractSystem
         rename(x, renamespace(sys, nameof(x)))
     else
-        Symbol(getname(sys), :₊, x)
+        Symbol(getname(sys), :., x)
     end
 end
 
@@ -1248,7 +1248,7 @@ function round_trip_eq(eq::Equation, var2name)
         syss = get_systems(eq.rhs)
         call = Expr(:call, connect)
         for sys in syss
-            strs = split(string(nameof(sys)), "₊")
+            strs = split(string(nameof(sys)), ".")
             s = Symbol(strs[1])
             for st in strs[2:end]
                 s = Expr(:., s, Meta.quot(Symbol(st)))

src/systems/callbacks.jl

@@ -512,6 +512,24 @@ function generate_rootfinding_callback(cbs, sys::AbstractODESystem, dvs = unknow
     end
 end
 
+# Put a wrapper on NamedTuple so that u.resistor.v indexes like u.var"resistor.v"
+# Required for hierarchical, but a hack that should be fixed in the future
+struct NamedTupleSymbolFix{T}
+    x::T
+    sym::Symbol
+end
+NamedTupleSymbolFix(x) = NamedTupleSymbolFix(x, Symbol(""))
+function Base.getproperty(u::NamedTupleSymbolFix, s::Symbol)
+    newsym = getfield(u, :sym) == Symbol("") ? s : Symbol(getfield(u, :sym), ".", s)
+    x = getfield(u, :x)
+    if newsym in keys(x)
+        getproperty(x, newsym)
+    else
+        NamedTupleSymbolFix(x, newsym)
+    end
+end
+Base.getindex(u::NamedTupleSymbolFix, idxs::Int...) = getfield(u, :x)[idxs...]
+
 function compile_user_affect(affect::FunctionalAffect, sys, dvs, ps; kwargs...)
     dvs_ind = Dict(reverse(en) for en in enumerate(dvs))
     v_inds = map(sym -> dvs_ind[sym], unknowns(affect))
@@ -526,9 +544,10 @@ function compile_user_affect(affect::FunctionalAffect, sys, dvs, ps; kwargs...)
     # HACK: filter out eliminated symbols. Not clear this is the right thing to do
     # (MTK should keep these symbols)
     u = filter(x -> !isnothing(x[2]), collect(zip(unknowns_syms(affect), v_inds))) |>
-        NamedTuple
+        NamedTuple |> NamedTupleSymbolFix
+
     p = filter(x -> !isnothing(x[2]), collect(zip(parameters_syms(affect), p_inds))) |>
-        NamedTuple
+        NamedTuple |> NamedTupleSymbolFix
 
     let u = u, p = p, user_affect = func(affect), ctx = context(affect)
         function (integ)

src/systems/connectors.jl

@@ -129,7 +129,7 @@ function generate_isouter(sys::AbstractSystem)
         function isouter(sys)::Bool
             s = string(nameof(sys))
             isconnector(sys) || error("$s is not a connector!")
-            idx = findfirst(isequal('₊'), s)
+            idx = findfirst(isequal('.'), s)
             parent_name = Symbol(idx === nothing ? s : s[1:prevind(s, idx)])
             parent_name in outer_connectors
         end

test/funcaffect.jl

@@ -123,7 +123,7 @@ i8 = findfirst(==(8.0), sol[:t])
 ctx = [0]
 function affect4!(integ, u, p, ctx)
     ctx[1] += 1
-    @test u.resistor₊v == 1
+    @test u.resistor.v == 1
 end
 s1 = compose(
     ODESystem(Equation[], t, [], [], name = :s1,
@@ -137,7 +137,7 @@ sol = solve(prob, Tsit5())
 include("../examples/rc_model.jl")
 
 function affect5!(integ, u, p, ctx)
-    @test integ.u[u.capacitor₊v] ≈ 0.3
+    @test integ.u[u.capacitor.v] ≈ 0.3
     integ.ps[p.C] *= 200
 end