Merge pull request #1866 from SciML/myb/repeated

Repeated components optimization and make `@named` more intuitive

Author: YingboMa

src/systems/abstractsystem.jl

@@ -1,3 +1,5 @@
+const SYSTEM_COUNT = Threads.Atomic{UInt}(0)
+
 """
 ```julia
 calculate_tgrad(sys::AbstractTimeDependentSystem)
@@ -175,6 +177,7 @@ function complete(sys::AbstractSystem)
 end
 
 for prop in [:eqs
+             :tag
              :noiseeqs
              :iv
              :states
@@ -255,7 +258,7 @@ end
     end
 end
 
-rename(x::AbstractSystem, name) = @set x.name = name
+rename(x, name) = @set x.name = name
 
 function Base.propertynames(sys::AbstractSystem; private = false)
     if private
@@ -842,12 +845,26 @@ function _named(name, call, runtime = false)
         end
     end
 
+    is_sys_construction = Symbol("###__is_system_construction###")
     kws = call.args[2].args
+    for (i, kw) in enumerate(kws)
+        if Meta.isexpr(kw, (:(=), :kw))
+            kw.args[2] = :($is_sys_construction ? $(kw.args[2]) :
+                           $default_to_parentscope($(kw.args[2])))
+        elseif kw isa Symbol
+            rhs = :($is_sys_construction ? $(kw) : $default_to_parentscope($(kw)))
+            kws[i] = Expr(:kw, kw, rhs)
+        end
+    end
 
     if !any(kw -> (kw isa Symbol ? kw : kw.args[1]) == :name, kws) # don't overwrite `name` kwarg
         pushfirst!(kws, Expr(:kw, :name, runtime ? name : Meta.quot(name)))
     end
-    call
+    op = call.args[1]
+    quote
+        $is_sys_construction = ($op isa $DataType) && ($op <: $AbstractSystem)
+        $call
+    end
 end
 
 function _named_idxs(name::Symbol, idxs, call)
@@ -872,46 +889,44 @@ end
 """
     @named y = foo(x)
     @named y[1:10] = foo(x)
-    @named y 1:10 i -> foo(x*i)
-
-Rewrite `@named y = foo(x)` to `y = foo(x; name=:y)`.
-
-Rewrite `@named y[1:10] = foo(x)` to `y = map(i′->foo(x; name=Symbol(:y_, i′)), 1:10)`.
+    @named y 1:10 i -> foo(x*i)  # This is not recommended
 
-Rewrite `@named y 1:10 i -> foo(x*i)` to `y = map(i->foo(x*i; name=Symbol(:y_, i)), 1:10)`.
+Pass the LHS name to the model. When it's calling anything that's not an
+AbstractSystem, it wraps all keyword arguments in `default_to_parentscope` so
+that namespacing works intuitively when passing a symbolic default into a
+component.
 
 Examples:
-```julia
+```julia-repl
 julia> using ModelingToolkit
 
-julia> foo(i; name) = i, name
+julia> foo(i; name) = (; i, name)
 foo (generic function with 1 method)
 
 julia> x = 41
 41
 
 julia> @named y = foo(x)
-(41, :y)
+(i = 41, name = :y)
 
 julia> @named y[1:3] = foo(x)
-3-element Vector{Tuple{Int64, Symbol}}:
- (41, :y_1)
- (41, :y_2)
- (41, :y_3)
-
-julia> @named y 1:3 i -> foo(x*i)
-3-element Vector{Tuple{Int64, Symbol}}:
- (41, :y_1)
- (82, :y_2)
- (123, :y_3)
+3-element Vector{NamedTuple{(:i, :name), Tuple{Int64, Symbol}}}:
+ (i = 41, name = :y_1)
+ (i = 41, name = :y_2)
+ (i = 41, name = :y_3)
 ```
 """
 macro named(expr)
     name, call = split_assign(expr)
     if Meta.isexpr(name, :ref)
         name, idxs = name.args
         check_name(name)
-        esc(_named_idxs(name, idxs, :($(gensym()) -> $call)))
+        var = gensym(name)
+        ex = quote
+            $var = $(_named(name, call))
+            $name = map(i -> $rename($var, Symbol($(Meta.quot(name)), :_, i)), $idxs)
+        end
+        esc(ex)
     else
         check_name(name)
         esc(:($name = $(_named(name, call))))
@@ -922,6 +937,11 @@ macro named(name::Symbol, idxs, call)
     esc(_named_idxs(name, idxs, call))
 end
 
+function default_to_parentscope(v)
+    uv = unwrap(v)
+    uv isa Symbolic && !hasmetadata(uv, SymScope) ? ParentScope(v) : v
+end
+
 function _config(expr, namespace)
     cn = Base.Fix2(_config, namespace)
     if Meta.isexpr(expr, :.)

src/systems/diffeqs/odesystem.jl

@@ -23,6 +23,11 @@ eqs = [D(x) ~ σ*(y-x),
 ```
 """
 struct ODESystem <: AbstractODESystem
+    """
+    tag: a tag for the system. If two system have the same tag, then they are
+    structurally identical.
+    """
+    tag::UInt
     """The ODEs defining the system."""
     eqs::Vector{Equation}
     """Independent variable."""
@@ -120,7 +125,7 @@ struct ODESystem <: AbstractODESystem
     """
     complete::Bool
 
-    function ODESystem(deqs, iv, dvs, ps, var_to_name, ctrls, observed, tgrad,
+    function ODESystem(tag, deqs, iv, dvs, ps, var_to_name, ctrls, observed, tgrad,
                        jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults,
                        torn_matching, connector_type, preface, cevents,
                        devents, metadata = nothing, tearing_state = nothing,
@@ -135,7 +140,7 @@ struct ODESystem <: AbstractODESystem
         if checks == true || (checks & CheckUnits) > 0
             all_dimensionless([dvs; ps; iv]) || check_units(deqs)
         end
-        new(deqs, iv, dvs, ps, var_to_name, ctrls, observed, tgrad, jac,
+        new(tag, deqs, iv, dvs, ps, var_to_name, ctrls, observed, tgrad, jac,
             ctrl_jac, Wfact, Wfact_t, name, systems, defaults, torn_matching,
             connector_type, preface, cevents, devents, metadata, tearing_state,
             substitutions, complete)
@@ -189,7 +194,8 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps;
     end
     cont_callbacks = SymbolicContinuousCallbacks(continuous_events)
     disc_callbacks = SymbolicDiscreteCallbacks(discrete_events)
-    ODESystem(deqs, iv′, dvs′, ps′, var_to_name, ctrl′, observed, tgrad, jac,
+    ODESystem(Threads.atomic_add!(SYSTEM_COUNT, UInt(1)),
+              deqs, iv′, dvs′, ps′, var_to_name, ctrl′, observed, tgrad, jac,
               ctrl_jac, Wfact, Wfact_t, name, systems, defaults, nothing,
               connector_type, preface, cont_callbacks, disc_callbacks,
               metadata, checks = checks)

src/systems/diffeqs/sdesystem.jl

@@ -27,6 +27,11 @@ noiseeqs = [0.1*x,
 ```
 """
 struct SDESystem <: AbstractODESystem
+    """
+    tag: a tag for the system. If two system have the same tag, then they are
+    structurally identical.
+    """
+    tag::UInt
     """The expressions defining the drift term."""
     eqs::Vector{Equation}
     """The expressions defining the diffusion term."""
@@ -105,7 +110,8 @@ struct SDESystem <: AbstractODESystem
     """
     complete::Bool
 
-    function SDESystem(deqs, neqs, iv, dvs, ps, var_to_name, ctrls, observed, tgrad, jac,
+    function SDESystem(tag, deqs, neqs, iv, dvs, ps, var_to_name, ctrls, observed, tgrad,
+                       jac,
                        ctrl_jac, Wfact, Wfact_t, name, systems, defaults, connector_type,
                        cevents, devents, metadata = nothing, complete = false;
                        checks::Union{Bool, Int} = true)
@@ -118,7 +124,8 @@ struct SDESystem <: AbstractODESystem
         if checks == true || (checks & CheckUnits) > 0
             all_dimensionless([dvs; ps; iv]) || check_units(deqs, neqs)
         end
-        new(deqs, neqs, iv, dvs, ps, var_to_name, ctrls, observed, tgrad, jac, ctrl_jac,
+        new(tag, deqs, neqs, iv, dvs, ps, var_to_name, ctrls, observed, tgrad, jac,
+            ctrl_jac,
             Wfact, Wfact_t, name, systems, defaults, connector_type, cevents, devents,
             metadata, complete)
     end
@@ -169,7 +176,8 @@ function SDESystem(deqs::AbstractVector{<:Equation}, neqs, iv, dvs, ps;
     cont_callbacks = SymbolicContinuousCallbacks(continuous_events)
     disc_callbacks = SymbolicDiscreteCallbacks(discrete_events)
 
-    SDESystem(deqs, neqs, iv′, dvs′, ps′, var_to_name, ctrl′, observed, tgrad, jac,
+    SDESystem(Threads.atomic_add!(SYSTEM_COUNT, UInt(1)),
+              deqs, neqs, iv′, dvs′, ps′, var_to_name, ctrl′, observed, tgrad, jac,
               ctrl_jac, Wfact, Wfact_t, name, systems, defaults, connector_type,
               cont_callbacks, disc_callbacks, metadata; checks = checks)
 end

src/systems/discrete_system/discrete_system.jl

@@ -24,6 +24,11 @@ eqs = [D(x) ~ σ*(y-x),
 ```
 """
 struct DiscreteSystem <: AbstractTimeDependentSystem
+    """
+    tag: a tag for the system. If two system have the same tag, then they are
+    structurally identical.
+    """
+    tag::UInt
     """The differential equations defining the discrete system."""
     eqs::Vector{Equation}
     """Independent variable."""
@@ -76,7 +81,8 @@ struct DiscreteSystem <: AbstractTimeDependentSystem
     """
     complete::Bool
 
-    function DiscreteSystem(discreteEqs, iv, dvs, ps, var_to_name, ctrls, observed, name,
+    function DiscreteSystem(tag, discreteEqs, iv, dvs, ps, var_to_name, ctrls, observed,
+                            name,
                             systems, defaults, preface, connector_type,
                             metadata = nothing,
                             tearing_state = nothing, substitutions = nothing,
@@ -88,7 +94,8 @@ struct DiscreteSystem <: AbstractTimeDependentSystem
         if checks == true || (checks & CheckUnits) > 0
             all_dimensionless([dvs; ps; iv; ctrls]) || check_units(discreteEqs)
         end
-        new(discreteEqs, iv, dvs, ps, var_to_name, ctrls, observed, name, systems, defaults,
+        new(tag, discreteEqs, iv, dvs, ps, var_to_name, ctrls, observed, name, systems,
+            defaults,
             preface, connector_type, metadata, tearing_state, substitutions, complete)
     end
 end
@@ -134,7 +141,8 @@ function DiscreteSystem(eqs::AbstractVector{<:Equation}, iv, dvs, ps;
     if length(unique(sysnames)) != length(sysnames)
         throw(ArgumentError("System names must be unique."))
     end
-    DiscreteSystem(eqs, iv′, dvs′, ps′, var_to_name, ctrl′, observed, name, systems,
+    DiscreteSystem(Threads.atomic_add!(SYSTEM_COUNT, UInt(1)),
+                   eqs, iv′, dvs′, ps′, var_to_name, ctrl′, observed, name, systems,
                    defaults, preface, connector_type, metadata, kwargs...)
 end
 

src/systems/jumps/jumpsystem.jl

@@ -48,6 +48,11 @@ j₃      = MassActionJump(2*β+γ, [R => 1], [S => 1, R => -1])
 ```
 """
 struct JumpSystem{U <: ArrayPartition} <: AbstractTimeDependentSystem
+    """
+    tag: a tag for the system. If two system have the same tag, then they are
+    structurally identical.
+    """
+    tag::UInt
     """
     The jumps of the system. Allowable types are `ConstantRateJump`,
     `VariableRateJump`, `MassActionJump`.
@@ -92,7 +97,7 @@ struct JumpSystem{U <: ArrayPartition} <: AbstractTimeDependentSystem
     """
     complete::Bool
 
-    function JumpSystem{U}(ap::U, iv, states, ps, var_to_name, observed, name, systems,
+    function JumpSystem{U}(tag, ap::U, iv, states, ps, var_to_name, observed, name, systems,
                            defaults, connector_type, devents,
                            metadata = nothing, complete = false;
                            checks::Union{Bool, Int} = true) where {U <: ArrayPartition}
@@ -103,7 +108,7 @@ struct JumpSystem{U <: ArrayPartition} <: AbstractTimeDependentSystem
         if checks == true || (checks & CheckUnits) > 0
             all_dimensionless([states; ps; iv]) || check_units(ap, iv)
         end
-        new{U}(ap, iv, states, ps, var_to_name, observed, name, systems, defaults,
+        new{U}(tag, ap, iv, states, ps, var_to_name, observed, name, systems, defaults,
                connector_type, devents, metadata, complete)
     end
 end
@@ -156,7 +161,8 @@ function JumpSystem(eqs, iv, states, ps;
         error("JumpSystems currently only support discrete events.")
     disc_callbacks = SymbolicDiscreteCallbacks(discrete_events)
 
-    JumpSystem{typeof(ap)}(ap, value(iv), states, ps, var_to_name, observed, name, systems,
+    JumpSystem{typeof(ap)}(Threads.atomic_add!(SYSTEM_COUNT, UInt(1)),
+                           ap, value(iv), states, ps, var_to_name, observed, name, systems,
                            defaults, connector_type, disc_callbacks, metadata,
                            checks = checks)
 end

src/systems/nonlinear/nonlinearsystem.jl

@@ -19,6 +19,11 @@ eqs = [0 ~ σ*(y-x),
 ```
 """
 struct NonlinearSystem <: AbstractTimeIndependentSystem
+    """
+    tag: a tag for the system. If two system have the same tag, then they are
+    structurally identical.
+    """
+    tag::UInt
     """Vector of equations defining the system."""
     eqs::Vector{Equation}
     """Unknown variables."""
@@ -67,14 +72,15 @@ struct NonlinearSystem <: AbstractTimeIndependentSystem
     """
     complete::Bool
 
-    function NonlinearSystem(eqs, states, ps, var_to_name, observed, jac, name, systems,
+    function NonlinearSystem(tag, eqs, states, ps, var_to_name, observed, jac, name,
+                             systems,
                              defaults, connector_type, metadata = nothing,
                              tearing_state = nothing, substitutions = nothing,
                              complete = false; checks::Union{Bool, Int} = true)
         if checks == true || (checks & CheckUnits) > 0
             all_dimensionless([states; ps]) || check_units(eqs)
         end
-        new(eqs, states, ps, var_to_name, observed, jac, name, systems, defaults,
+        new(tag, eqs, states, ps, var_to_name, observed, jac, name, systems, defaults,
             connector_type, metadata, tearing_state, substitutions, complete)
     end
 end
@@ -124,7 +130,8 @@ function NonlinearSystem(eqs, states, ps;
     process_variables!(var_to_name, defaults, ps)
     isempty(observed) || collect_var_to_name!(var_to_name, (eq.lhs for eq in observed))
 
-    NonlinearSystem(eqs, states, ps, var_to_name, observed, jac, name, systems, defaults,
+    NonlinearSystem(Threads.atomic_add!(SYSTEM_COUNT, UInt(1)),
+                    eqs, states, ps, var_to_name, observed, jac, name, systems, defaults,
                     connector_type, metadata, checks = checks)
 end
 

src/systems/optimization/optimizationsystem.jl

@@ -17,6 +17,11 @@ op = a*(y-x) + x*(b-z)-y + x*y - c*z
 ```
 """
 struct OptimizationSystem <: AbstractTimeIndependentSystem
+    """
+    tag: a tag for the system. If two system have the same tag, then they are
+    structurally identical.
+    """
+    tag::UInt
     """Objective function of the system."""
     op::Any
     """Unknown variables."""
@@ -46,15 +51,15 @@ struct OptimizationSystem <: AbstractTimeIndependentSystem
     """
     complete::Bool
 
-    function OptimizationSystem(op, states, ps, var_to_name, observed,
+    function OptimizationSystem(tag, op, states, ps, var_to_name, observed,
                                 constraints, name, systems, defaults, metadata = nothing,
                                 complete = false; checks::Union{Bool, Int} = true)
         if checks == true || (checks & CheckUnits) > 0
             unwrap(op) isa Symbolic && check_units(op)
             check_units(observed)
             all_dimensionless([states; ps]) || check_units(constraints)
         end
-        new(op, states, ps, var_to_name, observed,
+        new(tag, op, states, ps, var_to_name, observed,
             constraints, name, systems, defaults, metadata, complete)
     end
 end
@@ -92,7 +97,8 @@ function OptimizationSystem(op, states, ps;
     process_variables!(var_to_name, defaults, ps′)
     isempty(observed) || collect_var_to_name!(var_to_name, (eq.lhs for eq in observed))
 
-    OptimizationSystem(value(op), states′, ps′, var_to_name,
+    OptimizationSystem(Threads.atomic_add!(SYSTEM_COUNT, UInt(1)),
+                       value(op), states′, ps′, var_to_name,
                        observed,
                        constraints,
                        name, systems, defaults, metadata; checks = checks)