Merge branch 'master' into myb/strong_alias

Author: YingboMa

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"
@@ -91,11 +91,12 @@ OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 ReferenceTests = "324d217c-45ce-50fc-942e-d289b448e8cf"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
 StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["AmplNLWriter", "BenchmarkTools", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "OrdinaryDiffEq", "Random", "ReferenceTests", "SafeTestsets", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials"]
+test = ["AmplNLWriter", "BenchmarkTools", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "OrdinaryDiffEq", "Random", "ReferenceTests", "SafeTestsets", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials"]

docs/src/basics/FAQ.md

@@ -41,5 +41,58 @@ ERROR: TypeError: non-boolean (Num) used in boolean context
 
 then it's likely you are trying to trace through a function which cannot be
 directly represented in Julia symbols. The techniques to handle this problem,
-such as `@register_symbolic`, are described in detail 
+such as `@register_symbolic`, are described in detail
 [in the Symbolics.jl documentation](https://symbolics.juliasymbolics.org/dev/manual/faq/#Transforming-my-function-to-a-symbolic-equation-has-failed.-What-do-I-do?-1).
+
+## Using ModelingToolkit with Optimization / Automatic Differentiation
+
+If you are using ModelingToolkit inside of a loss function and are having issues with
+mixing MTK with automatic differentiation, getting performance, etc... don't! Instead, use
+MTK outside of the loss function to generate the code, and then use the generated code
+inside of the loss function.
+
+For example, let's say you were building ODEProblems in the loss function like:
+
+```julia
+function loss(p)
+    prob = ODEProblem(sys, [], [p1 => p[1], p2 => p[2]])
+    sol = solve(prob, Tsit5())
+    sum(abs2,sol)
+end
+```
+
+Since `ODEProblem` on a MTK `sys` will have to generate code, this will be slower than
+caching the generated code, and will required automatic differentiation to go through the
+code generation process itself. All of this is unnecessary. Instead, generate the problem
+once outside of the loss function, and remake the prob inside of the loss function:
+
+```julia
+prob = ODEProblem(sys, [], [p1 => p[1], p2 => p[2]])
+function loss(p)
+    remake(prob,p = ...)
+    sol = solve(prob, Tsit5())
+    sum(abs2,sol)
+end
+```
+
+Now, one has to be careful with `remake` to ensure that the parameters are in the right
+order. One can use the previously mentioned indexing functionality to generate index
+maps for reordering `p` like:
+
+```julia
+p = @parameters x y z
+idxs = ModelingToolkit.varmap_to_vars([p[1] => 1, p[2] => 2, p[3] => 3], p)
+p[idxs]
+```
+
+Using this, the fixed index map can be used in the loss function. This would look like:
+
+```julia
+prob = ODEProblem(sys, [], [p1 => p[1], p2 => p[2]])
+idxs = Int.(ModelingToolkit.varmap_to_vars([p1 => 1, p2 => 2], p))
+function loss(p)
+    remake(prob,p = p[idxs])
+    sol = solve(prob, Tsit5())
+    sum(abs2,sol)
+end
+```

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/systems/SDESystem.md

@@ -24,6 +24,7 @@ sde = SDESystem(ode, noiseeqs)
 ```@docs
 structural_simplify
 alias_elimination
+Girsanov_transform
 ```
 
 ## Analyses

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

@@ -173,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,7 +3,14 @@ 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))
     elseif istree(x) && operation(x) == (getindex)

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⋮")
@@ -965,6 +973,21 @@ function structural_simplify(sys::AbstractSystem, io = nothing; simplify = false
     return has_io ? (sys, input_idxs) : sys
 end
 
+function io_preprocessing(sys::AbstractSystem, inputs,
+                          outputs; simplify = false, kwargs...)
+    sys, input_idxs = structural_simplify(sys, (inputs, outputs); simplify, kwargs...)
+
+    eqs = equations(sys)
+    alg_start_idx = findfirst(!isdiffeq, eqs)
+    if alg_start_idx === nothing
+        alg_start_idx = length(eqs) + 1
+    end
+    diff_idxs = 1:alg_start_idx - 1
+    alge_idxs = alg_start_idx:length(eqs)
+
+    sys, diff_idxs, alge_idxs, input_idxs
+end
+
 """
     lin_fun, simplified_sys = linearization_function(sys::AbstractSystem, inputs, outputs; simplify = false, kwargs...)
 
@@ -990,21 +1013,11 @@ The `simplified_sys` has undergone [`structural_simplify`](@ref) and had any occ
 See also [`linearize`](@ref) which provides a higher-level interface.
 """
 function linearization_function(sys::AbstractSystem, inputs,
-                                outputs; kwargs...)
-    sys, input_idxs = structural_simplify(sys, (inputs, outputs); 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)
+                                outputs; simplify = false,
+                                kwargs...)
+    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)