Merge branch 'SciML:master' into func-affect

Author: dodoplus

Project.toml

@@ -1,7 +1,7 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Chris Rackauckas <[email protected]>"]
-version = "8.15.2"
+version = "8.16.0"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"

docs/src/basics/Variable_metadata.md

@@ -1,6 +1,38 @@
 # Symbolic metadata
-It is possible to add metadata to symbolic variables. The following
-information can be added (note, it's possible to extend this to user-defined metadata as well)
+It is possible to add metadata to symbolic variables, the metadata will be displayed when calling help on a variable.
+
+The following information can be added (note, it's possible to extend this to user-defined metadata as well)
+
+## Variable descriptions
+Descriptive strings can be attached to variables using the `[description = "descriptive string"]` syntax:
+```@example metadata
+using ModelingToolkit
+@variables u [description = "This is my input"]
+getdescription(u)
+```
+
+When variables with descriptions are present in systems, they will be printed when the system is shown in the terminal:
+```@example metadata
+@parameters t
+@variables u(t) [description = "A short description of u"]
+@parameters p   [description = "A description of p"]
+@named sys = ODESystem([u ~ p], t)
+show(stdout, "text/plain", sys) # hide
+```
+
+Calling help on the variable `u` displays the description, alongside other metadata:
+```julia
+help?> u
+
+  A variable of type Symbolics.Num (Num wraps anything in a type that is a subtype of Real)
+
+  Metadata
+  ≡≡≡≡≡≡≡≡≡≡
+
+  ModelingToolkit.VariableDescription: This is my input
+
+  Symbolics.VariableSource: (:variables, :u)
+```
 
 ## Input or output
 Designate a variable as either an input or an output using the following
@@ -34,7 +66,7 @@ isdisturbance(u)
 ```
 
 ## Mark parameter as tunable
-Indicate that a parameter can be automatically tuned by automatic control tuning apps.
+Indicate that a parameter can be automatically tuned by parameter optimization or automatic control tuning apps.
 
 ```@example metadata
 @parameters Kp [tunable=true]

docs/src/tutorials/ode_modeling.md

@@ -125,7 +125,7 @@ algebraic variables as "observables" (see
 That means, MTK still knows how to calculate them out of the information available
 in a simulation result. The intermediate variable `RHS` therefore can be plotted
 along with the state variable. Note that this has to be requested explicitly,
-though:
+through:
 
 ```julia
 prob = ODEProblem(fol_simplified, [x => 0.0], (0.0,10.0), [τ => 3.0])

src/ModelingToolkit.jl

@@ -41,6 +41,7 @@ import SymbolicUtils.Rewriters: Chain, Postwalk, Prewalk, Fixpoint
 import JuliaFormatter
 
 using Reexport
+using Symbolics: degree
 @reexport using Symbolics
 export @derivatives
 using Symbolics: _parse_vars, value, @derivatives, get_variables,
@@ -172,7 +173,7 @@ export NonlinearSystem, OptimizationSystem
 export alias_elimination, flatten
 export connect, @connector, Connection, Flow, Stream, instream
 export isinput, isoutput, getbounds, hasbounds, isdisturbance, istunable, getdist, hasdist,
-       tunable_parameters, isirreducible
+       tunable_parameters, isirreducible, getdescription, hasdescription
 export ode_order_lowering, dae_order_lowering, liouville_transform
 export PDESystem
 export Differential, expand_derivatives, @derivatives

src/inputoutput.jl

@@ -170,38 +170,28 @@ The return values also include the remaining states and parameters, in the order
 # Example
 ```
 using ModelingToolkit: generate_control_function, varmap_to_vars, defaults
-f, dvs, ps = generate_control_function(sys, expression=Val{false}, simplify=true)
+f, dvs, ps = generate_control_function(sys, expression=Val{false}, simplify=false)
 p = varmap_to_vars(defaults(sys), ps)
 x = varmap_to_vars(defaults(sys), dvs)
 t = 0
 f[1](x, inputs, p, t)
 ```
 """
-function generate_control_function(sys::AbstractODESystem;
+function generate_control_function(sys::AbstractODESystem, inputs = unbound_inputs(sys);
                                    implicit_dae = false,
-                                   has_difference = false,
-                                   simplify = true,
+                                   simplify = false,
                                    kwargs...)
-    ctrls = unbound_inputs(sys)
-    if isempty(ctrls)
+    if isempty(inputs)
         error("No unbound inputs were found in system.")
     end
 
-    # One can either connect unbound inputs to new parameters and allow structural_simplify, but then the unbound inputs appear as states :( .
-    # One can also just remove them from the states and parameters for the purposes of code generation, but then structural_simplify fails :(
-    # To have the best of both worlds, all unbound inputs must be converted to `@parameters` in which case structural_simplify handles them correctly :)
-    sys = toparam(sys, ctrls)
-
-    if simplify
-        sys = structural_simplify(sys)
-    end
+    sys, diff_idxs, alge_idxs = io_preprocessing(sys, inputs, []; simplify,
+                                                 check_bound = false, kwargs...)
 
     dvs = states(sys)
     ps = parameters(sys)
-
-    dvs = setdiff(dvs, ctrls)
-    ps = setdiff(ps, ctrls)
-    inputs = map(x -> time_varying_as_func(value(x), sys), ctrls)
+    ps = setdiff(ps, inputs)
+    inputs = map(x -> time_varying_as_func(value(x), sys), inputs)
 
     eqs = [eq for eq in equations(sys) if !isdifferenceeq(eq)]
     check_operator_variables(eqs, Differential)
@@ -223,24 +213,10 @@ function generate_control_function(sys::AbstractODESystem;
     end
     pre, sol_states = get_substitutions_and_solved_states(sys)
     f = build_function(rhss, args...; postprocess_fbody = pre, states = sol_states,
-                       kwargs...)
+                       expression = Val{false}, kwargs...)
     f, dvs, ps
 end
 
-"""
-    toparam(sys, ctrls::AbstractVector)
-
-Transform all instances of `@varibales` in `ctrls` appearing as states and in equations of `sys` with similarly named `@parameters`. This allows [`structural_simplify`](@ref)(sys) in the presence unbound inputs.
-"""
-function toparam(sys, ctrls::AbstractVector)
-    eqs = equations(sys)
-    subs = Dict(ctrls .=> toparam.(ctrls))
-    eqs = map(eqs) do eq
-        substitute(eq.lhs, subs) ~ substitute(eq.rhs, subs)
-    end
-    ODESystem(eqs, name = nameof(sys))
-end
-
 function inputs_to_parameters!(state::TransformationState, check_bound = true)
     @unpack structure, fullvars, sys = state
     @unpack var_to_diff, graph, solvable_graph = structure

src/parameters.jl

@@ -3,9 +3,16 @@ struct MTKParameterCtx end
 
 function isparameter(x)
     x = unwrap(x)
-    if istree(x) && operation(x) isa Symbolic
+
+    #TODO: Delete this branch
+    if x isa Symbolic && Symbolics.getparent(x, false) !== false
+        p = Symbolics.getparent(x)
+        isparameter(p) ||
+        (hasmetadata(p, Symbolics.VariableSource) &&
+         getmetadata(p, Symbolics.VariableSource)[1] == :parameters)
+    elseif istree(x) && operation(x) isa Symbolic
         getmetadata(x, MTKParameterCtx, false) ||
-        isparameter(operation(x))
+            isparameter(operation(x))
     elseif istree(x) && operation(x) == (getindex)
         isparameter(arguments(x)[1])
     elseif x isa Symbolic

src/systems/abstractsystem.jl

@@ -752,6 +752,10 @@ function Base.show(io::IO, mime::MIME"text/plain", sys::AbstractSystem)
                                :displaysize => (1, displaysize(io)[2])), val)
                 print(io, "]")
             end
+            description = getdescription(s)
+            if description !== nothing && description != ""
+                print(io, ": ", description)
+            end
         end
     end
     limited && print(io, "\n⋮")
@@ -774,6 +778,10 @@ function Base.show(io::IO, mime::MIME"text/plain", sys::AbstractSystem)
                                :displaysize => (1, displaysize(io)[2])), val)
                 print(io, "]")
             end
+            description = getdescription(s)
+            if description !== nothing && description != ""
+                print(io, ": ", description)
+            end
         end
     end
     limited && print(io, "\n⋮")
@@ -948,10 +956,17 @@ function will be applied during the tearing process. It also takes kwargs
 `allow_symbolic=false` and `allow_parameter=true` which limits the coefficient
 types during tearing.
 """
-function structural_simplify(sys::AbstractSystem; simplify = false, kwargs...)
+function structural_simplify(sys::AbstractSystem, args...; kwargs...)
     sys = expand_connections(sys)
     state = TearingState(sys)
-    state, = inputs_to_parameters!(state)
+    sys, input_idxs = _structural_simplify(sys, state, args...; kwargs...)
+    sys
+end
+
+function _structural_simplify(sys::AbstractSystem, state; simplify = false,
+                              check_bound = true,
+                              kwargs...)
+    state, input_idxs = inputs_to_parameters!(state, check_bound)
     sys = alias_elimination!(state)
     state = TearingState(sys)
     check_consistency(state)
@@ -964,7 +979,31 @@ function structural_simplify(sys::AbstractSystem; simplify = false, kwargs...)
     fullstates = [map(eq -> eq.lhs, observed(sys)); states(sys)]
     @set! sys.observed = topsort_equations(observed(sys), fullstates)
     invalidate_cache!(sys)
-    return sys
+    return sys, input_idxs
+end
+
+function io_preprocessing(sys::AbstractSystem, inputs,
+                          outputs; simplify = false, kwargs...)
+    sys = expand_connections(sys)
+    state = TearingState(sys)
+    markio!(state, inputs, outputs)
+    sys, input_idxs = _structural_simplify(sys, state; simplify, check_bound = false,
+                                           kwargs...)
+
+    eqs = equations(sys)
+    check_operator_variables(eqs, Differential)
+    # Sort equations and states such that diff.eqs. match differential states and the rest are algebraic
+    diffstates = collect_operator_variables(sys, Differential)
+    eqs = sort(eqs, by = e -> !isoperator(e.lhs, Differential),
+               alg = Base.Sort.DEFAULT_STABLE)
+    @set! sys.eqs = eqs
+    diffstates = [arguments(e.lhs)[1] for e in eqs[1:length(diffstates)]]
+    sts = [diffstates; setdiff(states(sys), diffstates)]
+    @set! sys.states = sts
+    diff_idxs = 1:length(diffstates)
+    alge_idxs = (length(diffstates) + 1):length(sts)
+
+    sys, diff_idxs, alge_idxs, input_idxs
 end
 
 """
@@ -994,36 +1033,9 @@ See also [`linearize`](@ref) which provides a higher-level interface.
 function linearization_function(sys::AbstractSystem, inputs,
                                 outputs; simplify = false,
                                 kwargs...)
-    sys = expand_connections(sys)
-    state = TearingState(sys)
-    markio!(state, inputs, outputs)
-    state, input_idxs = inputs_to_parameters!(state, false)
-    sys = alias_elimination!(state)
-    state = TearingState(sys)
-    check_consistency(state)
-    if sys isa ODESystem
-        sys = dae_order_lowering(dummy_derivative(sys, state))
-    end
-    state = TearingState(sys)
-    find_solvables!(state; kwargs...)
-    sys = tearing_reassemble(state, tearing(state), simplify = simplify)
-    fullstates = [map(eq -> eq.lhs, observed(sys)); states(sys)]
-    @set! sys.observed = topsort_equations(observed(sys), fullstates)
-    invalidate_cache!(sys)
-
-    eqs = equations(sys)
-    check_operator_variables(eqs, Differential)
-    # Sort equations and states such that diff.eqs. match differential states and the rest are algebraic
-    diffstates = collect_operator_variables(sys, Differential)
-    eqs = sort(eqs, by = e -> !isoperator(e.lhs, Differential),
-               alg = Base.Sort.DEFAULT_STABLE)
-    @set! sys.eqs = eqs
-    diffstates = [arguments(e.lhs)[1] for e in eqs[1:length(diffstates)]]
-    sts = [diffstates; setdiff(states(sys), diffstates)]
-    @set! sys.states = sts
-
-    diff_idxs = 1:length(diffstates)
-    alge_idxs = (length(diffstates) + 1):length(sts)
+    sys, diff_idxs, alge_idxs, input_idxs = io_preprocessing(sys, inputs, outputs; simplify,
+                                                             kwargs...)
+    sts = states(sys)
     fun = ODEFunction(sys)
     lin_fun = let fun = fun,
         h = ModelingToolkit.build_explicit_observed_function(sys, outputs)
@@ -1210,11 +1222,11 @@ function linearize(sys, lin_fun; t = 0.0, op = Dict(), allow_input_derivatives =
         C = [
         h_x h_z
 ]
-        Bs = -(gz \ (f_x * f_u + g_u))
+        Bs = -(gz \ g_u) # This equation differ from the cited paper, the paper is likely wrong since their equaiton leads to a dimension mismatch.
         if !iszero(Bs)
             if !allow_input_derivatives
                 der_inds = findall(vec(any(!=(0), Bs, dims = 1)))
-                error("Input derivatives appeared in expressions (-g_z\\(f_x*f_u + g_u) != 0), the following inputs appeared differentiated: $(inputs(sys)[der_inds]). Call `linear_staespace` with keyword argument `allow_input_derivatives = true` to allow this and have the returned `B` matrix be of double width ($(2nu)), where the last $nu inputs are the derivatives of the first $nu inputs.")
+                error("Input derivatives appeared in expressions (-g_z\\g_u != 0), the following inputs appeared differentiated: $(inputs(sys)[der_inds]). Call `linear_staespace` with keyword argument `allow_input_derivatives = true` to allow this and have the returned `B` matrix be of double width ($(2nu)), where the last $nu inputs are the derivatives of the first $nu inputs.")
             end
             B = [B Bs]
         end

src/utils.jl

@@ -223,7 +223,8 @@ function iv_from_nested_derivative(x, op = Differential)
     if istree(x) && operation(x) == getindex
         iv_from_nested_derivative(arguments(x)[1], op)
     elseif istree(x)
-        operation(x) isa op ? iv_from_nested_derivative(arguments(x)[1], op) : arguments(x)[1]
+        operation(x) isa op ? iv_from_nested_derivative(arguments(x)[1], op) :
+        arguments(x)[1]
     elseif issym(x)
         x
     else

src/variables.jl

@@ -1,14 +1,12 @@
 struct VariableUnit end
 struct VariableConnectType end
 struct VariableNoiseType end
-struct VariableDescriptionType end
 struct VariableInput end
 struct VariableOutput end
 struct VariableIrreducible end
 Symbolics.option_to_metadata_type(::Val{:unit}) = VariableUnit
 Symbolics.option_to_metadata_type(::Val{:connect}) = VariableConnectType
 Symbolics.option_to_metadata_type(::Val{:noise}) = VariableNoiseType
-Symbolics.option_to_metadata_type(::Val{:description}) = VariableDescriptionType
 Symbolics.option_to_metadata_type(::Val{:input}) = VariableInput
 Symbolics.option_to_metadata_type(::Val{:output}) = VariableOutput
 Symbolics.option_to_metadata_type(::Val{:irreducible}) = VariableIrreducible
@@ -267,3 +265,24 @@ function getbounds(p::AbstractVector)
     ub = last.(bounds)
     (; lb, ub)
 end
+
+## Description =================================================================
+struct VariableDescription end
+Symbolics.option_to_metadata_type(::Val{:description}) = VariableDescription
+
+getdescription(x::Num) = getdescription(Symbolics.unwrap(x))
+
+"""
+    getdescription(x)
+
+Return any description attached to variables `x`. If no description is attached, an empty string is returned.
+"""
+function getdescription(x)
+    p = Symbolics.getparent(x, nothing)
+    p === nothing || (x = p)
+    Symbolics.getmetadata(x, VariableDescription, "")
+end
+
+function hasdescription(x)
+    getdescription(x) != ""
+end

test/input_output_handling.jl

@@ -108,8 +108,7 @@ eqs = [
 ]
 
 @named sys = ODESystem(eqs)
-f, dvs, ps = ModelingToolkit.generate_control_function(sys, expression = Val{false},
-                                                       simplify = true)
+f, dvs, ps = ModelingToolkit.generate_control_function(sys, simplify = true)
 
 @test isequal(dvs[], x)
 @test isempty(ps)
@@ -170,8 +169,7 @@ eqs = [connect_sd(sd, mass1, mass2)
 @named _model = ODESystem(eqs, t)
 @named model = compose(_model, mass1, mass2, sd);
 
-f, dvs, ps = ModelingToolkit.generate_control_function(model, expression = Val{false},
-                                                       simplify = true)
+f, dvs, ps = ModelingToolkit.generate_control_function(model, simplify = true)
 @test length(dvs) == 4
 @test length(ps) == length(parameters(model))
 p = ModelingToolkit.varmap_to_vars(ModelingToolkit.defaults(model), ps)