Merge branch 'master' into speed

Author: SebastianM-C

src/systems/abstractsystem.jl

@@ -746,38 +746,64 @@ end
 abstract type SymScope end
 
 struct LocalScope <: SymScope end
-function LocalScope(sym::Union{Num, Symbolic})
+function LocalScope(sym::Union{Num, Symbolic, Symbolics.Arr{Num}})
     apply_to_variables(sym) do sym
-        setmetadata(sym, SymScope, LocalScope())
+        if istree(sym) && operation(sym) === getindex
+            args = arguments(sym)
+            a1 = setmetadata(args[1], SymScope, LocalScope())
+            similarterm(sym, operation(sym), [a1, args[2:end]...])
+        else
+            setmetadata(sym, SymScope, LocalScope())
+        end
     end
 end
 
 struct ParentScope <: SymScope
     parent::SymScope
 end
-function ParentScope(sym::Union{Num, Symbolic})
+function ParentScope(sym::Union{Num, Symbolic, Symbolics.Arr{Num}})
     apply_to_variables(sym) do sym
-        setmetadata(sym, SymScope,
-            ParentScope(getmetadata(value(sym), SymScope, LocalScope())))
+        if istree(sym) && operation(sym) === getindex
+            args = arguments(sym)
+            a1 = setmetadata(args[1], SymScope,
+                ParentScope(getmetadata(value(args[1]), SymScope, LocalScope())))
+            similarterm(sym, operation(sym), [a1, args[2:end]...])
+        else
+            setmetadata(sym, SymScope,
+                ParentScope(getmetadata(value(sym), SymScope, LocalScope())))
+        end
     end
 end
 
 struct DelayParentScope <: SymScope
     parent::SymScope
     N::Int
 end
-function DelayParentScope(sym::Union{Num, Symbolic}, N)
+function DelayParentScope(sym::Union{Num, Symbolic, Symbolics.Arr{Num}}, N)
     apply_to_variables(sym) do sym
-        setmetadata(sym, SymScope,
-            DelayParentScope(getmetadata(value(sym), SymScope, LocalScope()), N))
+        if istree(sym) && operation(sym) == getindex
+            args = arguments(sym)
+            a1 = setmetadata(args[1], SymScope,
+                DelayParentScope(getmetadata(value(args[1]), SymScope, LocalScope()), N))
+            similarterm(sym, operation(sym), [a1, args[2:end]...])
+        else
+            setmetadata(sym, SymScope,
+                DelayParentScope(getmetadata(value(sym), SymScope, LocalScope()), N))
+        end
     end
 end
-DelayParentScope(sym::Union{Num, Symbolic}) = DelayParentScope(sym, 1)
+DelayParentScope(sym::Union{Num, Symbolic, Symbolics.Arr{Num}}) = DelayParentScope(sym, 1)
 
 struct GlobalScope <: SymScope end
-function GlobalScope(sym::Union{Num, Symbolic})
+function GlobalScope(sym::Union{Num, Symbolic, Symbolics.Arr{Num}})
     apply_to_variables(sym) do sym
-        setmetadata(sym, SymScope, GlobalScope())
+        if istree(sym) && operation(sym) == getindex
+            args = arguments(sym)
+            a1 = setmetadata(args[1], SymScope, GlobalScope())
+            similarterm(sym, operation(sym), [a1, args[2:end]...])
+        else
+            setmetadata(sym, SymScope, GlobalScope())
+        end
     end
 end
 
@@ -793,6 +819,11 @@ function renamespace(sys, x)
             return similarterm(x, operation(x),
                 Any[renamespace(sys, only(arguments(x)))])::T
         end
+        if istree(x) && operation(x) === getindex
+            args = arguments(x)
+            return similarterm(
+                x, operation(x), vcat(renamespace(sys, args[1]), args[2:end]))::T
+        end
         let scope = getmetadata(x, SymScope, LocalScope())
             if scope isa LocalScope
                 rename(x, renamespace(getname(sys), getname(x)))::T
@@ -849,7 +880,8 @@ function namespace_assignment(eq::Assignment, sys)
     Assignment(_lhs, _rhs)
 end
 
-function namespace_expr(O, sys, n = nameof(sys); ivs = independent_variables(sys))
+function namespace_expr(
+        O, sys, n = nameof(sys); ivs = independent_variables(sys))
     O = unwrap(O)
     if any(isequal(O), ivs)
         return O
@@ -1500,8 +1532,7 @@ function default_to_parentscope(v)
     uv isa Symbolic || return v
     apply_to_variables(v) do sym
         if !hasmetadata(uv, SymScope)
-            setmetadata(sym, SymScope,
-                ParentScope(getmetadata(value(sym), SymScope, LocalScope())))
+            ParentScope(sym)
         else
             sym
         end

src/systems/diffeqs/abstractodesystem.jl

@@ -876,31 +876,34 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap;
             parammap = Dict(unwrap.(parameters(sys)) .=> parammap)
         end
     end
-    clockedparammap = Dict()
-    defs = ModelingToolkit.get_defaults(sys)
-    for v in ps
-        v = unwrap(v)
-        is_discrete_domain(v) || continue
-        op = operation(v)
-        if !isa(op, Symbolics.Operator) && parammap != SciMLBase.NullParameters() &&
-           haskey(parammap, v)
-            error("Initial conditions for discrete variables must be for the past state of the unknown. Instead of providing the condition for $v, provide the condition for $(Shift(iv, -1)(v)).")
+
+    if has_discrete_subsystems(sys) && get_discrete_subsystems(sys) !== nothing
+        clockedparammap = Dict()
+        defs = ModelingToolkit.get_defaults(sys)
+        for v in ps
+            v = unwrap(v)
+            is_discrete_domain(v) || continue
+            op = operation(v)
+            if !isa(op, Symbolics.Operator) && parammap != SciMLBase.NullParameters() &&
+               haskey(parammap, v)
+                error("Initial conditions for discrete variables must be for the past state of the unknown. Instead of providing the condition for $v, provide the condition for $(Shift(iv, -1)(v)).")
+            end
+            shiftedv = StructuralTransformations.simplify_shifts(Shift(iv, -1)(v))
+            if parammap != SciMLBase.NullParameters() &&
+               (val = get(parammap, shiftedv, nothing)) !== nothing
+                clockedparammap[v] = val
+            elseif op isa Shift
+                root = arguments(v)[1]
+                haskey(defs, root) || error("Initial condition for $v not provided.")
+                clockedparammap[v] = defs[root]
+            end
         end
-        shiftedv = StructuralTransformations.simplify_shifts(Shift(iv, -1)(v))
-        if parammap != SciMLBase.NullParameters() &&
-           (val = get(parammap, shiftedv, nothing)) !== nothing
-            clockedparammap[v] = val
-        elseif op isa Shift
-            root = arguments(v)[1]
-            haskey(defs, root) || error("Initial condition for $v not provided.")
-            clockedparammap[v] = defs[root]
+        parammap = if parammap == SciMLBase.NullParameters()
+            clockedparammap
+        else
+            merge(parammap, clockedparammap)
         end
     end
-    parammap = if parammap == SciMLBase.NullParameters()
-        clockedparammap
-    else
-        merge(parammap, clockedparammap)
-    end
     # TODO: make it work with clocks
     # ModelingToolkit.get_tearing_state(sys) !== nothing => Requires structural_simplify first
     if sys isa ODESystem && build_initializeprob &&
@@ -931,7 +934,12 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap;
     if has_index_cache(sys) && get_index_cache(sys) !== nothing
         u0, defs = get_u0(sys, trueinit, parammap; symbolic_u0)
         check_eqs_u0(eqs, dvs, u0; kwargs...)
-        p = MTKParameters(sys, parammap, trueinit)
+        p = if parammap === nothing ||
+               parammap == SciMLBase.NullParameters() && isempty(defs)
+            nothing
+        else
+            MTKParameters(sys, parammap, trueinit)
+        end
     else
         u0, p, defs = get_u0_p(sys,
             trueinit,
@@ -1592,7 +1600,6 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem,
     if !iscomplete(sys)
         error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `ODEProblem`")
     end
-
     if isempty(u0map) && get_initializesystem(sys) !== nothing
         isys = get_initializesystem(sys)
     elseif isempty(u0map) && get_initializesystem(sys) === nothing
@@ -1620,9 +1627,9 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem,
         @warn "Initialization system is underdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares. To suppress this warning pass warn_initialize_determined = false."
     end
 
-    parammap isa DiffEqBase.NullParameters || isempty(parammap) ?
-    [get_iv(sys) => t] :
-    merge(todict(parammap), Dict(get_iv(sys) => t))
+    parammap = parammap isa DiffEqBase.NullParameters || isempty(parammap) ?
+               [get_iv(sys) => t] :
+               merge(todict(parammap), Dict(get_iv(sys) => t))
 
     if neqs == nunknown
         NonlinearProblem(isys, guesses, parammap)

src/systems/nonlinear/nonlinearsystem.jl

@@ -241,7 +241,6 @@ function generate_function(
     pre, sol_states = get_substitutions_and_solved_unknowns(sys)
 
     p = reorder_parameters(sys, value.(ps))
-    @show p ps
     return build_function(rhss, value.(dvs), p...; postprocess_fbody = pre,
         states = sol_states, kwargs...)
 end
@@ -395,7 +394,6 @@ function process_NonlinearProblem(constructor, sys::NonlinearSystem, u0map, para
     eqs = equations(sys)
     dvs = unknowns(sys)
     ps = full_parameters(sys)
-
     if has_index_cache(sys) && get_index_cache(sys) !== nothing
         u0, defs = get_u0(sys, u0map, parammap)
         check_eqs_u0(eqs, dvs, u0; kwargs...)

src/systems/parameter_buffer.jl

@@ -43,14 +43,32 @@ function MTKParameters(
     p = merge(defs, p)
     p = merge(Dict(unwrap(k) => v for (k, v) in p),
         Dict(default_toterm(unwrap(k)) => v for (k, v) in p))
-    p = Dict(k => fixpoint_sub(v, p) for (k, v) in p)
+    p = Dict(unwrap(k) => fixpoint_sub(v, p) for (k, v) in p)
     for (sym, _) in p
         if istree(sym) && operation(sym) === getindex &&
            first(arguments(sym)) in all_ps
             error("Scalarized parameter values ($sym) are not supported. Instead of `[p[1] => 1.0, p[2] => 2.0]` use `[p => [1.0, 2.0]]`")
         end
     end
 
+    missing_params = Set()
+    for idxmap in (ic.tunable_idx, ic.discrete_idx, ic.constant_idx, ic.nonnumeric_idx)
+        for sym in keys(idxmap)
+            sym isa Symbol && continue
+            haskey(p, sym) && continue
+            hasname(sym) && haskey(p, getname(sym)) && continue
+            ttsym = default_toterm(sym)
+            haskey(p, ttsym) && continue
+            hasname(ttsym) && haskey(p, getname(ttsym)) && continue
+
+            istree(sym) && operation(sym) === getindex && haskey(p, arguments(sym)[1]) &&
+                continue
+            push!(missing_params, sym)
+        end
+    end
+
+    isempty(missing_params) || throw(MissingVariablesError(collect(missing_params)))
+
     tunable_buffer = Tuple(Vector{temp.type}(undef, temp.length)
     for temp in ic.tunable_buffer_sizes)
     disc_buffer = Tuple(Vector{temp.type}(undef, temp.length)

test/input_output_handling.jl

@@ -20,9 +20,9 @@ end
 @named sys = ODESystem([D(x) ~ -x + u], t) # both u and x are unbound
 @named sys1 = ODESystem([D(x) ~ -x + v[1] + v[2]], t) # both v and x are unbound
 @named sys2 = ODESystem([D(x) ~ -sys.x], t, systems = [sys]) # this binds sys.x in the context of sys2, sys2.x is still unbound
-@named sys21 = ODESystem([D(x) ~ -sys.x], t, systems = [sys1]) # this binds sys.x in the context of sys2, sys2.x is still unbound
+@named sys21 = ODESystem([D(x) ~ -sys1.x], t, systems = [sys1]) # this binds sys.x in the context of sys2, sys2.x is still unbound
 @named sys3 = ODESystem([D(x) ~ -sys.x + sys.u], t, systems = [sys]) # This binds both sys.x and sys.u
-@named sys31 = ODESystem([D(x) ~ -sys.x + sys1.v[1]], t, systems = [sys1]) # This binds both sys.x and sys1.v[1]
+@named sys31 = ODESystem([D(x) ~ -sys1.x + sys1.v[1]], t, systems = [sys1]) # This binds both sys.x and sys1.v[1]
 
 @named sys4 = ODESystem([D(x) ~ -sys.x, u ~ sys.u], t, systems = [sys]) # This binds both sys.x and sys3.u, this system is one layer deeper than the previous. u is directly forwarded to sys.u, and in this case sys.u is bound while u is not
 
@@ -43,7 +43,7 @@ end
 @test is_bound(sys2, sys.x)
 @test !is_bound(sys2, sys.u)
 @test !is_bound(sys2, sys2.sys.u)
-@test is_bound(sys21, sys.x)
+@test is_bound(sys21, sys1.x)
 @test !is_bound(sys21, sys1.v[1])
 @test !is_bound(sys21, sys1.v[2])
 @test is_bound(sys31, sys1.v[1])

test/mtkparameters.jl

@@ -123,6 +123,20 @@ newps = remake_buffer(sys, ps, Dict(p => 1.0f0))
 @test newps.tunable[1] isa Vector{Float32}
 @test newps.tunable[1] == [1.0f0, 2.0f0, 3.0f0]
 
+# Issue#2624
+@parameters p d
+@variables X(t)
+eqs = [D(X) ~ p - d * X]
+@mtkbuild sys = ODESystem(eqs, t)
+
+u0 = [X => 1.0]
+tspan = (0.0, 100.0)
+ps = [p => 1.0] # Value for `d` is missing
+
+@test_throws ModelingToolkit.MissingVariablesError ODEProblem(sys, u0, tspan, ps)
+@test_nowarn ODEProblem(sys, u0, tspan, [ps..., d => 1.0])
+
+  
 # JET tests
 
 # scalar parameters only

test/odesystem.jl

@@ -1066,3 +1066,57 @@ prob = SteadyStateProblem(sys, u0, p)
 @test prob isa SteadyStateProblem
 prob = SteadyStateProblem(ODEProblem(sys, u0, (0.0, 10.0), p))
 @test prob isa SteadyStateProblem
+
+# Issue#2344
+using ModelingToolkitStandardLibrary.Blocks
+
+function FML2(; name)
+    @parameters begin
+        k2[1:1] = [1.0]
+    end
+    systems = @named begin
+        constant = Constant(k = k2[1])
+    end
+    @variables begin
+        x(t) = 0
+    end
+    eqs = [
+        D(x) ~ constant.output.u + k2[1]
+    ]
+    ODESystem(eqs, t; systems, name)
+end
+
+@mtkbuild model = FML2()
+
+@test isequal(ModelingToolkit.defaults(model)[model.constant.k], model.k2[1])
+@test_nowarn ODEProblem(model, [], (0.0, 10.0))
+
+# Issue#2477
+function RealExpression(; name, y)
+    vars = @variables begin
+        u(t)
+    end
+    eqns = [
+        u ~ y
+    ]
+    sys = ODESystem(eqns, t, vars, []; name)
+end
+
+function RealExpressionSystem(; name)
+    vars = @variables begin
+        x(t)
+        z(t)[1:1]
+    end # doing a collect on z doesn't work either. 
+    @named e1 = RealExpression(y = x) # This works perfectly. 
+    @named e2 = RealExpression(y = z[1]) # This bugs. However, `full_equations(e2)` works as expected. 
+    systems = [e1, e2]
+    ODESystem(Equation[], t, Iterators.flatten(vars), []; systems, name)
+end
+
+@named sys = RealExpressionSystem()
+sys = complete(sys)
+@test Set(equations(sys)) == Set([sys.e1.u ~ sys.x, sys.e2.u ~ sys.z[1]])
+tearing_state = TearingState(expand_connections(sys))
+ts_vars = tearing_state.fullvars
+orig_vars = unknowns(sys)
+@test isempty(setdiff(ts_vars, orig_vars))

test/variable_scope.jl

@@ -73,3 +73,13 @@ ps = ModelingToolkit.getname.(parameters(level3))
 @test isequal(ps[4], :level2₊level0₊d)
 @test isequal(ps[5], :level1₊level0₊e)
 @test isequal(ps[6], :f)
+
+# Issue@2252
+# Tests from PR#2354
+@parameters xx[1:2]
+arr_p = [ParentScope(xx[1]), xx[2]]
+arr0 = ODESystem(Equation[], t, [], arr_p; name = :arr0)
+arr1 = ODESystem(Equation[], t, [], []; name = :arr1) ∘ arr0
+arr_ps = ModelingToolkit.getname.(parameters(arr1))
+@test isequal(arr_ps[1], Symbol("xx"))
+@test isequal(arr_ps[2], Symbol("arr0₊xx"))