Merge pull request #905 from SciML/myb/defaults

Use `defaults` and variable mapping clean up

Author: YingboMa

docs/src/basics/AbstractSystem.md

@@ -22,7 +22,7 @@ same keyword arguments, which are:
 
 - `system`: This is used for specifying subsystems for hierarchical modeling with
   reusable components. For more information, see the [components page](@ref components)
-- Defaults: Keyword arguments like `default_u0` are used for specifying default
+- Defaults: Keyword arguments like `defaults` are used for specifying default
   values which are used. If a value is not given at the `SciMLProblem` construction
   time, its numerical value will be the default.
 
@@ -49,9 +49,7 @@ Optionally, a system could have:
 
 - `observed(sys)`: All observed equations of the system and its subsystems.
 - `get_observed(sys)`: Observed equations of the current-level system.
-- `get_default_u0(sys)`: A `Dict` that maps states into their default initial
-  condition.
-- `get_default_p(sys)`: A `Dict` that maps parameters into their default value.
+- `get_defaults(sys)`: A `Dict` that maps variables into their default values.
 - `independent_variable(sys)`: The independent variable of a system.
 - `get_noiseeqs(sys)`: Noise equations of the current-level system.
 
@@ -132,7 +130,7 @@ u0 = [
 ## Default Value Handling
 
 The `AbstractSystem` types allow for specifying default values, for example
-`default_p` inside of them. At problem construction time, these values are merged
+`defaults` inside of them. At problem construction time, these values are merged
 into the value maps, where for any repeats the value maps override the default.
 In addition, defaults of a higher level in the system override the defaults of
-a lower level in the system. 
+a lower level in the system.

docs/src/tutorials/acausal_components.md

@@ -19,7 +19,7 @@ using ModelingToolkit, Plots, DifferentialEquations
 # Basic electric components
 function Pin(;name)
     @variables v(t) i(t)
-    ODESystem(Equation[], t, [v, i], [], name=name, default_u0=[v=>1.0, i=>1.0])
+    ODESystem(Equation[], t, [v, i], [], name=name, defaults=[v=>1.0, i=>1.0])
 end
 
 function Ground(;name)
@@ -39,7 +39,7 @@ function Resistor(;name, R = 1.0)
            0 ~ p.i + n.i
            v ~ p.i * R
           ]
-    ODESystem(eqs, t, [v], [R], systems=[p, n], default_p=Dict(R => val), name=name)
+    ODESystem(eqs, t, [v], [R], systems=[p, n], defaults=Dict(R => val), name=name)
 end
 
 function Capacitor(; name, C = 1.0)
@@ -54,7 +54,7 @@ function Capacitor(; name, C = 1.0)
            0 ~ p.i + n.i
            D(v) ~ p.i / C
           ]
-    ODESystem(eqs, t, [v], [C], systems=[p, n], default_p=Dict(C => val), name=name)
+    ODESystem(eqs, t, [v], [C], systems=[p, n], defaults=Dict(C => val), name=name)
 end
 
 function ConstantVoltage(;name, V = 1.0)
@@ -66,7 +66,7 @@ function ConstantVoltage(;name, V = 1.0)
            V ~ p.v - n.v
            0 ~ p.i + n.i
           ]
-    ODESystem(eqs, t, [], [V], systems=[p, n], default_p=Dict(V => val), name=name)
+    ODESystem(eqs, t, [], [V], systems=[p, n], defaults=Dict(V => val), name=name)
 end
 
 R = 1.0
@@ -122,7 +122,7 @@ component to simply be the values there:
 ```julia
 function Pin(;name)
     @variables v(t) i(t)
-    ODESystem(Equation[], t, [v, i], [], name=name, default_u0=[v=>1.0, i=>1.0])
+    ODESystem(Equation[], t, [v, i], [], name=name, defaults=[v=>1.0, i=>1.0])
 end
 ```
 
@@ -176,11 +176,11 @@ function Resistor(;name, R = 1.0)
            0 ~ p.i + n.i
            v ~ p.i * R
           ]
-    ODESystem(eqs, t, [v], [R], systems=[p, n], default_p=Dict(R => val), name=name)
+    ODESystem(eqs, t, [v], [R], systems=[p, n], defaults=Dict(R => val), name=name)
 end
 ```
 
-Notice that we have created this system with a `default_p` for the resistor's
+Notice that we have created this system with a `defaults` for the resistor's
 resistance. By doing so, if the resistance of this resistor is not overridden
 by a higher level default or overridden at `ODEProblem` construction time, this
 will be the value of the resistance.
@@ -200,7 +200,7 @@ function Capacitor(; name, C = 1.0)
            0 ~ p.i + n.i
            D(v) ~ p.i / C
           ]
-    ODESystem(eqs, t, [v], [C], systems=[p, n], default_p=Dict(C => val), name=name)
+    ODESystem(eqs, t, [v], [C], systems=[p, n], defaults=Dict(C => val), name=name)
 end
 ```
 
@@ -219,7 +219,7 @@ function ConstantVoltage(;name, V = 1.0)
            V ~ p.v - n.v
            0 ~ p.i + n.i
           ]
-    ODESystem(eqs, t, [], [V], systems=[p, n], default_p=Dict(V => val), name=name)
+    ODESystem(eqs, t, [], [V], systems=[p, n], defaults=Dict(V => val), name=name)
 end
 ```
 

docs/src/tutorials/tearing_parallelism.md

@@ -41,7 +41,7 @@ const t = Sym{ModelingToolkit.Parameter{Real}}(:t)
 const D = Differential(t)
 function Pin(;name)
     @variables v(t) i(t)
-    ODESystem(Equation[], t, [v, i], [], name=name, default_u0=Dict([v=>1.0, i=>1.0]))
+    ODESystem(Equation[], t, [v, i], [], name=name, defaults=Dict([v=>1.0, i=>1.0]))
 end
 
 function Ground(;name)
@@ -59,12 +59,12 @@ function ConstantVoltage(;name, V = 1.0)
            V ~ p.v - n.v
            0 ~ p.i + n.i
           ]
-    ODESystem(eqs, t, [], [V], systems=[p, n], default_p=Dict(V => val), name=name)
+    ODESystem(eqs, t, [], [V], systems=[p, n], defaults=Dict(V => val), name=name)
 end
 
 function HeatPort(;name)
     @variables T(t) Q_flow(t)
-    return ODESystem(Equation[], t, [T, Q_flow], [], default_u0=Dict(T=>293.15, Q_flow=>0.0), name=name)
+    return ODESystem(Equation[], t, [T, Q_flow], [], defaults=Dict(T=>293.15, Q_flow=>0.0), name=name)
 end
 
 function HeatingResistor(;name, R=1.0, TAmbient=293.15, alpha=1.0)
@@ -83,8 +83,10 @@ function HeatingResistor(;name, R=1.0, TAmbient=293.15, alpha=1.0)
           ]
     ODESystem(
         eqs, t, [v, RTherm], [R, TAmbient, alpha], systems=[p, n, h],
-        default_p=Dict(R=>R_val, TAmbient=>TAmbient_val, alpha=>alpha_val),
-        default_u0=Dict(v=>0.0, RTherm=>R_val),
+        defaults=Dict(
+            R=>R_val, TAmbient=>TAmbient_val, alpha=>alpha_val,
+            v=>0.0, RTherm=>R_val
+        ),
         name=name,
     )
 end
@@ -99,7 +101,7 @@ function HeatCapacitor(;name, rho=8050, V=1, cp=460, TAmbient=293.15)
           ]
     ODESystem(
         eqs, t, [], [rho, V, cp], systems=[h],
-        default_p=Dict(rho=>rho_val, V=>V_val, cp=>cp_val),
+        defaults=Dict(rho=>rho_val, V=>V_val, cp=>cp_val),
         name=name,
     )
 end
@@ -117,8 +119,7 @@ function Capacitor(;name, C = 1.0)
           ]
     ODESystem(
         eqs, t, [v], [C], systems=[p, n],
-        default_u0=Dict(v => 0.0),
-        default_p=Dict(C => val),
+        defaults=Dict(v => 0.0, C => val),
         name=name
     )
 end
@@ -159,7 +160,7 @@ end
 eqs = [
        D(E) ~ sum(((i, sys),)->getproperty(sys, Symbol(:resistor, i)).h.Q_flow, enumerate(rc_systems))
       ]
-big_rc = ODESystem(eqs, t, [], [], systems=rc_systems, default_u0=Dict(E=>0.0))
+big_rc = ODESystem(eqs, t, [], [], systems=rc_systems, defaults=Dict(E=>0.0))
 ```
 
 Now let's say we want to expose a bit more parallelism via running tearing.

src/structural_transformation/StructuralTransformations.jl

@@ -16,7 +16,7 @@ using ModelingToolkit: ODESystem, var_from_nested_derivative, Differential,
                        states, equations, vars, Symbolic, diff2term, value,
                        operation, arguments, Sym, Term, simplify, solve_for,
                        isdiffeq, isdifferential,
-                       get_structure, get_reduced_states, default_u0, default_p
+                       get_structure, get_reduced_states, defaults
 
 using ModelingToolkit.BipartiteGraphs
 using ModelingToolkit.SystemStructures

src/structural_transformation/codegen.jl

@@ -131,7 +131,7 @@ function gen_nlsolve(sys, eqs, vars)
     allvars = unique(collect(Iterators.flatten(map(ModelingToolkit.vars, rhss))))
     params = setdiff(allvars, vars)
 
-    u0map = default_u0(sys)
+    u0map = defaults(sys)
     # splatting to tighten the type
     u0 = [map(var->get(u0map, var, 1e-3), vars)...]
     # specialize on the scalar case
@@ -338,22 +338,11 @@ function ODAEProblem{iip}(
     s = structure(sys)
     @unpack fullvars = s
     dvs = fullvars[diffvars_range(s)]
-    defaults = merge(default_p(sys), default_u0(sys))
-    u0map′ = ModelingToolkit.lower_mapnames(u0map, independent_variable(sys))
-    u0 = ModelingToolkit.varmap_to_vars(u0map′, dvs; defaults=defaults)
-
     ps = parameters(sys)
-    if parammap isa DiffEqBase.NullParameters && isempty(default_p(sys))
-        isempty(ps) || throw(ArgumentError("The model has non-empty parameters but no parameters are specified in the problem."))
-        p = parammap
-    else
-        if parammap isa DiffEqBase.NullParameters
-            pp = Pair[]
-        else
-            pp = ModelingToolkit.lower_mapnames(parammap)
-        end
-        p = ModelingToolkit.varmap_to_vars(pp, ps; defaults=defaults)
-    end
+    defs = defaults(sys)
+
+    u0 = ModelingToolkit.varmap_to_vars(u0map, dvs; defaults=defs)
+    p = ModelingToolkit.varmap_to_vars(parammap, ps; defaults=defs)
 
     ODEProblem{iip}(build_torn_function(sys; kw...), u0, tspan, p; kw...)
 end

src/systems/abstractsystem.jl

@@ -140,8 +140,7 @@ for prop in [
              :iv
              :states
              :ps
-             :default_p
-             :default_u0
+             :defaults
              :observed
              :tgrad
              :jac
@@ -155,6 +154,10 @@ for prop in [
              :controls
              :loss
              :reduced_states
+             :bcs
+             :domain
+             :depvars
+             :indvars
             ]
     fname1 = Symbol(:get_, prop)
     fname2 = Symbol(:has_, prop)
@@ -253,13 +256,13 @@ function Base.setproperty!(sys::AbstractSystem, prop::Symbol, val)
         idx = findfirst(s->getname(s) == prop, params);
         idx !== nothing;
        )
-        get_default_p(sys)[params[idx]] = value(val)
+        get_defaults(sys)[params[idx]] = value(val)
     elseif (
             sts = states(sys);
             idx = findfirst(s->getname(s) == prop, sts);
             idx !== nothing;
            )
-        get_default_u0(sys)[sts[idx]] = value(val)
+        get_defaults(sys)[sts[idx]] = value(val)
     else
         setfield!(sys, prop, val)
     end
@@ -280,14 +283,9 @@ end
 namespace_variables(sys::AbstractSystem) = states(sys, states(sys))
 namespace_parameters(sys::AbstractSystem) = parameters(sys, parameters(sys))
 
-function namespace_default_u0(sys)
-    d_u0 = default_u0(sys)
-    Dict(states(sys, k) => namespace_expr(d_u0[k], nameof(sys), independent_variable(sys)) for k in keys(d_u0))
-end
-
-function namespace_default_p(sys)
-    d_p = default_p(sys)
-    Dict(parameters(sys, k) => namespace_expr(d_p[k], nameof(sys), independent_variable(sys)) for k in keys(d_p))
+function namespace_defaults(sys)
+    defs = defaults(sys)
+    Dict((isparameter(k) ? parameters(sys, k) : states(sys, k)) => namespace_expr(defs[k], nameof(sys), independent_variable(sys)) for k in keys(defs))
 end
 
 function namespace_equations(sys::AbstractSystem)
@@ -344,16 +342,12 @@ function observed(sys::AbstractSystem)
             init=Equation[])]
 end
 
-function default_u0(sys::AbstractSystem)
+Base.@deprecate default_u0(x) defaults(x) false
+Base.@deprecate default_p(x) defaults(x) false
+function defaults(sys::AbstractSystem)
     systems = get_systems(sys)
-    d_u0 = get_default_u0(sys)
-    isempty(systems) ? d_u0 : mapreduce(namespace_default_u0, merge, systems; init=d_u0)
-end
-
-function default_p(sys::AbstractSystem)
-    systems = get_systems(sys)
-    d_p = get_default_p(sys)
-    isempty(systems) ? d_p : mapreduce(namespace_default_p, merge, systems; init=d_p)
+    defs = get_defaults(sys)
+    isempty(systems) ? defs : mapreduce(namespace_defaults, merge, systems; init=defs)
 end
 
 states(sys::AbstractSystem, v) = renamespace(nameof(sys), v)
@@ -398,9 +392,13 @@ function (f::AbstractSysToExpr)(O)
     return build_expr(:call, Any[operation(O); f.(arguments(O))])
 end
 
-function Base.show(io::IO, sys::AbstractSystem)
+function Base.show(io::IO, ::MIME"text/plain", sys::AbstractSystem)
     eqs = equations(sys)
-    Base.printstyled(io, "Model $(nameof(sys)) with $(length(eqs)) equations\n"; bold=true)
+    if eqs isa AbstractArray
+        Base.printstyled(io, "Model $(nameof(sys)) with $(length(eqs)) equations\n"; bold=true)
+    else
+        Base.printstyled(io, "Model $(nameof(sys))\n"; bold=true)
+    end
     # The reduced equations are usually very long. It's not that useful to print
     # them.
     #Base.print_matrix(io, eqs)
@@ -413,13 +411,13 @@ function Base.show(io::IO, sys::AbstractSystem)
     Base.printstyled(io, "States ($nvars):"; bold=true)
     nrows = min(nvars, limit ? rows : nvars)
     limited = nrows < length(vars)
-    d_u0 = has_default_u0(sys) ? default_u0(sys) : nothing
+    defs = has_defaults(sys) ? defaults(sys) : nothing
     for i in 1:nrows
         s = vars[i]
         print(io, "\n  ", s)
 
-        if d_u0 !== nothing
-            val = get(d_u0, s, nothing)
+        if defs !== nothing
+            val = get(defs, s, nothing)
             if val !== nothing
                 print(io, " [defaults to $val]")
             end
@@ -432,13 +430,12 @@ function Base.show(io::IO, sys::AbstractSystem)
     Base.printstyled(io, "Parameters ($nvars):"; bold=true)
     nrows = min(nvars, limit ? rows : nvars)
     limited = nrows < length(vars)
-    d_p = has_default_p(sys) ? default_p(sys) : nothing
     for i in 1:nrows
         s = vars[i]
         print(io, "\n  ", s)
 
-        if d_p !== nothing
-            val = get(d_p, s, nothing)
+        if defs !== nothing
+            val = get(defs, s, nothing)
             if val !== nothing
                 print(io, " [defaults to $val]")
             end

src/systems/control/controlsystem.jl

@@ -68,33 +68,30 @@ struct ControlSystem <: AbstractControlSystem
     """
     systems::Vector{ControlSystem}
     """
-    default_u0: The default initial conditions to use when initial conditions
-    are not supplied in `ODEProblem`.
+    defaults: The default values to use when initial conditions and/or
+    parameters are not supplied in `ODEProblem`.
     """
-    default_u0::Dict
-    """
-    default_p: The default parameters to use when parameters are not supplied
-    in `ODEProblem`.
-    """
-    default_p::Dict
+    defaults::Dict
 end
 
 function ControlSystem(loss, deqs::AbstractVector{<:Equation}, iv, dvs, controls, ps;
                        observed = [],
                        systems = ODESystem[],
                        default_u0=Dict(),
                        default_p=Dict(),
+                       defaults=_merge(Dict(default_u0), Dict(default_p)),
                        name=gensym(:ControlSystem))
+    if !(isempty(default_u0) && isempty(default_p))
+        Base.depwarn("`default_u0` and `default_p` are deprecated. Use `defaults` instead.", :ControlSystem, force=true)
+    end
     iv′ = value(iv)
     dvs′ = value.(dvs)
     controls′ = value.(controls)
     ps′ = value.(ps)
-    default_u0 isa Dict || (default_u0 = Dict(default_u0))
-    default_p isa Dict || (default_p = Dict(default_p))
-    default_u0 = Dict(value(k) => value(default_u0[k]) for k in keys(default_u0))
-    default_p = Dict(value(k) => value(default_p[k]) for k in keys(default_p))
+    defaults = todict(defaults)
+    defaults = Dict(value(k) => value(v) for (k, v) in pairs(defaults))
     ControlSystem(value(loss), deqs, iv′, dvs′, controls′,
-                  ps′, observed, name, systems, default_u0, default_p)
+                  ps′, observed, name, systems, defaults)
 end
 
 struct ControlToExpr