Merge branch 'SciML:master' into dg/metamacro

Author: DhairyaLGandhi

Project.toml

@@ -1,7 +1,7 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Yingbo Ma <[email protected]>", "Chris Rackauckas <[email protected]> and contributors"]
-version = "10.19.0"
+version = "10.21.0"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"

docs/src/API/model_building.md

@@ -227,6 +227,7 @@ add_accumulations
 noise_to_brownians
 convert_system_indepvar
 subset_tunables
+respecialize
 ```
 
 ## Hybrid systems

src/ModelingToolkit.jl

@@ -268,7 +268,8 @@ export isinput, isoutput, getbounds, hasbounds, getguess, hasguess, isdisturbanc
        hasmisc, getmisc, state_priority,
        subset_tunables
 export liouville_transform, change_independent_variable, substitute_component,
-       add_accumulations, noise_to_brownians, Girsanov_transform, change_of_variables
+       add_accumulations, noise_to_brownians, Girsanov_transform, change_of_variables,
+       respecialize
 export PDESystem
 export Differential, expand_derivatives, @derivatives
 export Equation, ConstrainedEquation

src/inputoutput.jl

@@ -5,7 +5,7 @@ using Symbolics: get_variables
 Return all variables that mare marked as inputs. See also [`unbound_inputs`](@ref)
 See also [`bound_inputs`](@ref), [`unbound_inputs`](@ref)
 """
-inputs(sys) = [filter(isinput, unknowns(sys)); filter(isinput, parameters(sys))]
+inputs(sys) = collect(get_inputs(sys))
 
 """
     outputs(sys)
@@ -14,13 +14,7 @@ Return all variables that mare marked as outputs. See also [`unbound_outputs`](@
 See also [`bound_outputs`](@ref), [`unbound_outputs`](@ref)
 """
 function outputs(sys)
-    o = observed(sys)
-    rhss = [eq.rhs for eq in o]
-    lhss = [eq.lhs for eq in o]
-    unique([filter(isoutput, unknowns(sys))
-            filter(isoutput, parameters(sys))
-            filter(x -> iscall(x) && isoutput(x), rhss) # observed can return equations with complicated expressions, we are only looking for single Terms
-            filter(x -> iscall(x) && isoutput(x), lhss)])
+    return collect(get_outputs(sys))
 end
 
 """
@@ -288,7 +282,12 @@ function inputs_to_parameters!(state::TransformationState, inputsyms)
             push!(new_fullvars, v)
         end
     end
-    ninputs == 0 && return state
+    if ninputs == 0
+        @set! sys.inputs = OrderedSet{BasicSymbolic}()
+        @set! sys.outputs = OrderedSet{BasicSymbolic}(filter(isoutput, fullvars))
+        state.sys = sys
+        return state
+    end
 
     nvars = ndsts(graph) - ninputs
     new_graph = BipartiteGraph(nsrcs(graph), nvars, Val(false))
@@ -318,6 +317,8 @@ function inputs_to_parameters!(state::TransformationState, inputsyms)
     ps = parameters(sys)
 
     @set! sys.ps = [ps; new_parameters]
+    @set! sys.inputs = OrderedSet{BasicSymbolic}(new_parameters)
+    @set! sys.outputs = OrderedSet{BasicSymbolic}(filter(isoutput, fullvars))
     @set! state.sys = sys
     @set! state.fullvars = Vector{BasicSymbolic}(new_fullvars)
     @set! state.structure = structure

src/linearization.jl

@@ -572,8 +572,7 @@ struct IONotFoundError <: Exception
 end
 
 function Base.showerror(io::IO, err::IONotFoundError)
-    println(io,
-        "The following $(err.variant) provided to `mtkcompile` were not found in the system:")
+    println(io, "The following $(err.variant) provided to `mtkcompile` were not found in the system:")
     maybe_namespace_issue = false
     for var in err.not_found
         println(io, "  ", var)

src/problems/odeproblem.jl

@@ -3,7 +3,8 @@
         t = nothing, eval_expression = false, eval_module = @__MODULE__, sparse = false,
         steady_state = false, checkbounds = false, sparsity = false, analytic = nothing,
         simplify = false, cse = true, initialization_data = nothing, expression = Val{false},
-        check_compatibility = true, nlstep = false, nlstep_compile = true, kwargs...) where {iip, spec}
+        check_compatibility = true, nlstep = false, nlstep_compile = true, nlstep_scc = false,
+        kwargs...) where {iip, spec}
     check_complete(sys, ODEFunction)
     check_compatibility && check_compatible_system(ODEFunction, sys)
 
@@ -42,7 +43,7 @@
     _M = concrete_massmatrix(M; sparse, u0)
 
     if nlstep
-        ode_nlstep = generate_ODENLStepData(sys, u0, p, M, nlstep_compile)
+        ode_nlstep = generate_ODENLStepData(sys, u0, p, M, nlstep_compile, nlstep_scc)
     else
         ode_nlstep = nothing
     end

src/systems/abstractsystem.jl

@@ -784,6 +784,8 @@ const SYS_PROPS = [:eqs
                    :parent
                    :is_dde
                    :tstops
+                   :inputs
+                   :outputs
                    :index_cache
                    :isscheduled
                    :costs
@@ -1820,6 +1822,17 @@ function push_vars!(stmt, name, typ, vars)
             ex = nameof(s)
         end
         push!(vars_expr.args, ex)
+
+        meta_kvps = Expr[]
+        if isinput(s)
+            push!(meta_kvps, :(input = true))
+        end
+        if isoutput(s)
+            push!(meta_kvps, :(output = true))
+        end
+        if !isempty(meta_kvps)
+            push!(vars_expr.args, Expr(:vect, meta_kvps...))
+        end
     end
     push!(stmt, :($name = $collect($vars_expr)))
     return

src/systems/callbacks.jl

@@ -25,6 +25,12 @@ function SymbolicAffect(affect::SymbolicAffect; kwargs...)
 end
 SymbolicAffect(affect; kwargs...) = make_affect(affect; kwargs...)
 
+function Symbolics.fast_substitute(aff::SymbolicAffect, rules)
+    substituter = Base.Fix2(fast_substitute, rules)
+    SymbolicAffect(map(substituter, aff.affect), map(substituter, aff.alg_eqs),
+        map(substituter, aff.discrete_parameters))
+end
+
 struct AffectSystem
     """The internal implicit discrete system whose equations are solved to obtain values after the affect."""
     system::AbstractSystem
@@ -36,6 +42,19 @@ struct AffectSystem
     discretes::Vector
 end
 
+function Symbolics.fast_substitute(aff::AffectSystem, rules)
+    substituter = Base.Fix2(fast_substitute, rules)
+    sys = aff.system
+    @set! sys.eqs = map(substituter, get_eqs(sys))
+    @set! sys.parameter_dependencies = map(substituter, get_parameter_dependencies(sys))
+    @set! sys.defaults = Dict([k => substituter(v) for (k, v) in defaults(sys)])
+    @set! sys.guesses = Dict([k => substituter(v) for (k, v) in guesses(sys)])
+    @set! sys.unknowns = map(substituter, get_unknowns(sys))
+    @set! sys.ps = map(substituter, get_ps(sys))
+    AffectSystem(sys, map(substituter, aff.unknowns),
+        map(substituter, aff.parameters), map(substituter, aff.discretes))
+end
+
 function AffectSystem(spec::SymbolicAffect; iv = nothing, alg_eqs = Equation[], kwargs...)
     AffectSystem(spec.affect; alg_eqs = vcat(spec.alg_eqs, alg_eqs), iv,
         discrete_parameters = spec.discrete_parameters, kwargs...)
@@ -855,7 +874,7 @@ end
 function default_operating_point(affsys::AffectSystem)
     sys = system(affsys)
 
-    op = Dict(unknowns(sys) .=> 0.0)
+    op = AnyDict(unknowns(sys) .=> 0.0)
     for p in parameters(sys)
         T = symtype(p)
         if T <: Number

src/systems/connectors.jl

@@ -759,6 +759,7 @@ function generate_connection_equations_and_stream_connections(
         var = variable_from_vertex(sys, cvert)::BasicSymbolic
         vtype = cvert.type
         if vtype <: Union{InputVar, OutputVar}
+            length(cset) > 1 || continue
             inner_output = nothing
             outer_input = nothing
             for cvert in cset
@@ -780,11 +781,11 @@ function generate_connection_equations_and_stream_connections(
                     inner_output = cvert
                 end
             end
-            root, rest = Iterators.peel(cset)
-            root_var = variable_from_vertex(sys, root)
-            for cvert in rest
-                var = variable_from_vertex(sys, cvert)
-                push!(eqs, root_var ~ var)
+            root_vert = something(inner_output, outer_input)
+            root_var = variable_from_vertex(sys, root_vert)
+            for cvert in cset
+                isequal(cvert, root_vert) && continue
+                push!(eqs, variable_from_vertex(sys, cvert) ~ root_var)
             end
         elseif vtype === Stream
             push!(stream_connections, cset)
@@ -807,10 +808,37 @@ function generate_connection_equations_and_stream_connections(
                 push!(eqs, 0 ~ rhs)
             end
         else # Equality
-            base = variable_from_vertex(sys, cset[1])
-            for i in 2:length(cset)
-                v = variable_from_vertex(sys, cset[i])
-                push!(eqs, base ~ v)
+            vars = map(Base.Fix1(variable_from_vertex, sys), cset)
+            outer_input = inner_output = nothing
+            all_io = true
+            # attempt to interpret the equality as a causal connectionset if
+            # possible
+            for (cvert, vert) in zip(cset, vars)
+                is_i = isinput(vert)
+                is_o = isoutput(vert)
+                all_io &= is_i || is_o
+                all_io || break
+                if cvert.isouter && is_i && outer_input === nothing
+                    outer_input = cvert
+                elseif !cvert.isouter && is_o && inner_output === nothing
+                    inner_output = cvert
+                end
+            end
+            # this doesn't necessarily mean this is a well-structured causal connection,
+            # but it is sufficient and we're generating equalities anyway.
+            if all_io && xor(outer_input !== nothing, inner_output !== nothing)
+                root_vert = something(inner_output, outer_input)
+                root_var = variable_from_vertex(sys, root_vert)
+                for (cvert, var) in zip(cset, vars)
+                    isequal(cvert, root_vert) && continue
+                    push!(eqs, var ~ root_var)
+                end
+            else
+                base = variable_from_vertex(sys, cset[1])
+                for i in 2:length(cset)
+                    v = vars[i]
+                    push!(eqs, base ~ v)
+                end
             end
         end
     end

src/systems/diffeqs/basic_transformations.jl

@@ -706,3 +706,143 @@ function convert_system_indepvar(sys::System, t; name = nameof(sys))
     @set! sys.var_to_name = var_to_name
     return sys
 end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Shorthand for `respecialize(sys, []; all = true)`
+"""
+respecialize(sys::AbstractSystem) = respecialize(sys, []; all = true)
+
+"""
+    $(TYPEDSIGNATURES)
+
+Specialize nonnumeric parameters in `sys` by changing their symtype to a concrete type.
+`mapping` is an iterable, where each element can be a parameter or a pair mapping a parameter
+to a value. If the element is a parameter, it must have a default. Each specified parameter
+is updated to have the symtype of the value associated with it (either in `mapping` or in
+the defaults). This operation can only be performed on nonnumeric, non-array parameters. The
+defaults of respecialized parameters are set to the associated values.
+
+This operation can only be performed on `complete`d systems.
+
+# Keyword arguments
+
+- `all`: Specialize all nonnumeric parameters in the system. This will error if any such
+  parameter does not have a default.
+"""
+function respecialize(sys::AbstractSystem, mapping; all = false)
+    if !iscomplete(sys)
+        error("""
+        This operation can only be performed on completed systems. Use `complete(sys)` or
+        `mtkcompile(sys)`.
+        """)
+    end
+    if !is_split(sys)
+        error("""
+        This operation can only be performed on split systems. Use `complete(sys)` or
+        `mtkcompile(sys)` with the `split = true` keyword argument.
+        """)
+    end
+
+    new_ps = copy(get_ps(sys))
+    @set! sys.ps = new_ps
+
+    extras = []
+    if all
+        for x in filter(!is_variable_numeric, get_ps(sys))
+            if any(y -> isequal(x, y) || y isa Pair && isequal(x, y[1]), mapping) ||
+               symbolic_type(x) === ArraySymbolic() ||
+               iscall(x) && operation(x) === getindex
+                continue
+            end
+            push!(extras, x)
+        end
+    end
+    ps_to_specialize = Iterators.flatten((extras, mapping))
+
+    defs = copy(defaults(sys))
+    @set! sys.defaults = defs
+    final_defs = copy(defs)
+    evaluate_varmap!(final_defs, ps_to_specialize)
+
+    subrules = Dict()
+
+    for element in ps_to_specialize
+        if element isa Pair
+            k, v = element
+        else
+            k = element
+            v = get(final_defs, k, nothing)
+            @assert v !== nothing """
+            Parameter $k needs an associated value to be respecialized.
+            """
+            @assert symbolic_type(v) == NotSymbolic() && !is_array_of_symbolics(v) """
+            Parameter $k needs an associated value to be respecialized. Found symbolic \
+            default $v.
+            """
+        end
+
+        k = unwrap(k)
+        T = typeof(v)
+
+        @assert !is_variable_numeric(k) """
+        Numeric types cannot be respecialized - tried to respecialize $k.
+        """
+        @assert symbolic_type(k) !== ArraySymbolic() """
+        Cannot respecialize array symbolics - tried to respecialize $k.
+        """
+        @assert !iscall(k) || operation(k) !== getindex """
+        Cannot respecialized scalarized array variables - tried to respecialize $k.
+        """
+        idx = findfirst(isequal(k), get_ps(sys))
+        @assert idx !== nothing """
+        Parameter $k does not exist in the system.
+        """
+
+        if iscall(k)
+            op = operation(k)
+            args = arguments(k)
+            new_p = SymbolicUtils.term(op, args...; type = T)
+        else
+            new_p = SymbolicUtils.Sym{T}(getname(k))
+        end
+
+        get_ps(sys)[idx] = new_p
+        defaults(sys)[new_p] = v
+        subrules[unwrap(k)] = unwrap(new_p)
+    end
+
+    substituter = Base.Fix2(fast_substitute, subrules)
+    @set! sys.eqs = map(substituter, get_eqs(sys))
+    @set! sys.observed = map(substituter, get_observed(sys))
+    @set! sys.initialization_eqs = map(substituter, get_initialization_eqs(sys))
+    if get_noise_eqs(sys) !== nothing
+        @set! sys.noise_eqs = map(substituter, get_noise_eqs(sys))
+    end
+    @set! sys.assertions = Dict([substituter(k) => v for (k, v) in assertions(sys)])
+    @set! sys.parameter_dependencies = map(substituter, get_parameter_dependencies(sys))
+    @set! sys.defaults = Dict([substituter(k) => substituter(v) for (k, v) in defaults(sys)])
+    @set! sys.guesses = Dict([k => substituter(v) for (k, v) in guesses(sys)])
+    @set! sys.continuous_events = map(get_continuous_events(sys)) do cev
+        SymbolicContinuousCallback(
+            map(substituter, cev.conditions), substituter(cev.affect),
+            substituter(cev.affect_neg), substituter(cev.initialize),
+            substituter(cev.finalize), cev.rootfind,
+            cev.reinitializealg, cev.zero_crossing_id)
+    end
+    @set! sys.discrete_events = map(get_discrete_events(sys)) do dev
+        SymbolicDiscreteCallback(map(substituter, dev.conditions), substituter(dev.affect),
+            substituter(dev.initialize), substituter(dev.finalize), dev.reinitializealg)
+    end
+    if get_schedule(sys) !== nothing
+        sched = get_schedule(sys)
+        @set! sys.schedule = Schedule(
+            sched.var_sccs, AnyDict(k => substituter(v) for (k, v) in sched.dummy_sub))
+    end
+    @set! sys.constraints = map(substituter, get_constraints(sys))
+    @set! sys.tstops = map(substituter, get_tstops(sys))
+    @set! sys.costs = Vector{Union{Real, BasicSymbolic}}(map(substituter, get_costs(sys)))
+    sys = complete(sys; split = is_split(sys))
+    return sys
+end