Merge branch 'master' of https://github.com/SciML/ModelingToolkit.jl into sy/difference_ODESys

Author: sharanry

.github/workflows/Downstream.yml

@@ -21,7 +21,8 @@ jobs:
           - {user: SciML, repo: CellMLToolkit.jl, group: All}
           - {user: SciML, repo: NeuralPDE.jl, group: NNPDE}
           - {user: SciML, repo: DataDrivenDiffEq.jl, group: Standard}
-
+          - {user: SciML, repo: StructuralIdentifiability.jl, group: All}
+          
     steps:
       - uses: actions/checkout@v2
       - uses: julia-actions/setup-julia@v1

Project.toml

@@ -1,7 +1,7 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Chris Rackauckas <[email protected]>"]
-version = "5.22.0"
+version = "5.23.0"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -56,7 +56,7 @@ Distributions = "0.23, 0.24, 0.25"
 DocStringExtensions = "0.7, 0.8"
 DomainSets = "0.5"
 IfElse = "0.1"
-JuliaFormatter = "0.12, 0.13, 0.14"
+JuliaFormatter = "0.12, 0.13, 0.14, 0.15"
 LabelledArrays = "1.3"
 Latexify = "0.11, 0.12, 0.13, 0.14, 0.15"
 LightGraphs = "1.3"

docs/src/basics/ContextualVariables.md

@@ -20,24 +20,36 @@ All modeling projects have some form of parameters. `@parameters` marks a variab
 as being the parameter of some system, which allows automatic detection algorithms
 to ignore such variables when attempting to find the states of a system.
 
-## Flow Variables (TODO)
+## Variable metadata [Experimental/TODO]
 
 In many engineering systems some variables act like "flows" while others do not.
 For example, in circuit models you have current which flows, and the related
 voltage which does not. Or in thermal models you have heat flows. In these cases,
 the `connect` statement enforces conservation of flow between all of the connected
 components.
 
-For example, the following specifies that `x` is a 2x2 matrix of flow variables:
+For example, the following specifies that `x` is a 2x2 matrix of flow variables
+with the unit m^3/s:
 
 ```julia
-@variables x[1:2,1:2] type=flow
+@variables x[1:2,1:2] [connect = Flow; unit = u"m^3/s"]
 ```
 
-## Stream Variables
+ModelingToolkit defines `connect`, `unit`, `noise`, and `description` keys for
+the metadata. One can get and set metadata by
 
-TODO
+```julia
+julia> @variables x [unit = u"m^3/s"];
+
+julia> hasmetadata(x, Symbolics.option_to_metadata_type(Val(:unit)))
+true
 
-## Brownian Variables
+julia> getmetadata(x, Symbolics.option_to_metadata_type(Val(:unit)))
+m³ s⁻¹
 
-TODO
+julia> x = setmetadata(x, Symbolics.option_to_metadata_type(Val(:unit)), u"m/s")
+x
+
+julia> getmetadata(x, Symbolics.option_to_metadata_type(Val(:unit)))
+m s⁻¹
+```

docs/src/tutorials/acausal_components.md

@@ -95,12 +95,11 @@ rc_eqs = [
           connect(capacitor.n, source.n, ground.g)
          ]
 
-@named rc_model = ODESystem(rc_eqs, t,
-                            systems=[resistor, capacitor, source, ground])
+@named rc_model = compose(ODESystem(rc_eqs, t),
+                          [resistor, capacitor, source, ground])
 sys = structural_simplify(rc_model)
 u0 = [
       capacitor.v => 0.0
-      capacitor.p.i => 0.0
      ]
 prob = ODAEProblem(sys, u0, (0, 10.0))
 sol = solve(prob, Tsit5())
@@ -289,15 +288,15 @@ rc_eqs = [
 Finally we build our four component model with these connection rules:
 
 ```julia
-@named rc_model = ODESystem(rc_eqs, t,
-                            systems=[resistor, capacitor, source, ground])
+@named rc_model = compose(ODESystem(rc_eqs, t)
+                          [resistor, capacitor, source, ground])
 ```
 
-Notice that this model is acasual because we have not specified anything about
-the causality of the model. We have simply specified what is true about each
-of the variables. This forms a system of differential-algebraic equations
-(DAEs) which define the evolution of each state of the system. The
-equations are:
+Note that we can also specify the subsystems in a vector. This model is acasual
+because we have not specified anything about the causality of the model. We have
+simply specified what is true about each of the variables. This forms a system
+of differential-algebraic equations (DAEs) which define the evolution of each
+state of the system. The equations are:
 
 ```julia
 equations(rc_model)
@@ -369,9 +368,10 @@ parameters(rc_model)
 This system could be solved directly as a DAE using [one of the DAE solvers
 from DifferentialEquations.jl](https://diffeq.sciml.ai/stable/solvers/dae_solve/).
 However, let's take a second to symbolically simplify the system before doing the
-solve. The function `structural_simplify` looks for all of the equalities and
-eliminates unnecessary variables to build the leanest numerical representation
-of the system. Let's see what it does here:
+solve. Although we can use ODE solvers that handles mass matrices to solve the
+above system directly, we want to run the `structural_simplify` function first,
+as it eliminates many unnecessary variables to build the leanest numerical
+representation of the system. Let's see what it does here:
 
 ```julia
 sys = structural_simplify(rc_model)
@@ -410,15 +410,13 @@ plot(sol)
 
 ![](https://user-images.githubusercontent.com/1814174/109416295-55184100-798b-11eb-96d1-5bb7e40135ba.png)
 
-However, we can also choose to use the "torn nonlinear system" to remove all
-of the algebraic variables from the solution of the system. Note that this
-requires having done `structural_simplify`. This is done by using `ODAEProblem`
-like:
+Since we have run `structural_simplify`, MTK can numerically solve all the
+unreduced algebraic equations numerically using the `ODAEProblem` (note the
+letter `A`):
 
 ```julia
 u0 = [
       capacitor.v => 0.0
-      capacitor.p.i => 0.0
      ]
 prob = ODAEProblem(sys, u0, (0, 10.0))
 sol = solve(prob, Rodas4())

docs/src/tutorials/tearing_parallelism.md

@@ -126,7 +126,8 @@ function rc_model(i; name, source, ground, R, C)
               connect_heat(resistor.h, heat_capacitor.h)
              ]
 
-    rc_model = ODESystem(rc_eqs, t, systems=[resistor, capacitor, source, ground, heat_capacitor], name=Symbol(name, i))
+    rc_model = compose(ODESystem(rc_eqs, t, name=Symbol(name, i)),
+                       [resistor, capacitor, source, ground, heat_capacitor])
 end
 ```
 
@@ -150,7 +151,7 @@ end;
 eqs = [
        D(E) ~ sum(((i, sys),)->getproperty(sys, Symbol(:resistor, i)).h.Q_flow, enumerate(rc_systems))
       ]
-big_rc = compose([ODESystem(eqs, t, [E], []); rc_systems])
+big_rc = compose(ODESystem(eqs, t, [E], []), rc_systems)
 ```
 
 Now let's say we want to expose a bit more parallelism via running tearing.

examples/rc_model.jl

@@ -14,4 +14,4 @@ rc_eqs = [
           connect(capacitor.n, source.n, ground.g)
          ]
 
-@named rc_model = compose(ODESystem(rc_eqs, t), resistor, capacitor, source, ground)
+@named rc_model = compose(ODESystem(rc_eqs, t), [resistor, capacitor, source, ground])

src/ModelingToolkit.jl

@@ -55,7 +55,7 @@ import Symbolics: rename, get_variables!, _solve, hessian_sparsity,
 
 import DiffEqBase: @add_kwonly
 
-import LightGraphs: SimpleDiGraph, add_edge!
+import LightGraphs: SimpleDiGraph, add_edge!, incidence_matrix
 
 using Requires
 
@@ -146,6 +146,8 @@ for S in subtypes(ModelingToolkit.AbstractSystem)
     @eval convert_system(::Type{<:$S}, sys::$S) = sys
 end
 
+struct Flow end
+
 export ODESystem, ODEFunction, ODEFunctionExpr, ODEProblemExpr, convert_system
 export DAEFunctionExpr, DAEProblemExpr
 export SDESystem, SDEFunction, SDEFunctionExpr, SDESystemExpr
@@ -191,7 +193,7 @@ export toexpr, get_variables
 export simplify, substitute
 export build_function
 export modelingtoolkitize
-export @variables, @parameters
+export @variables, @parameters, Flow
 export @named, @nonamespace, @namespace, extend, compose
 
 end # module

src/bipartite_graph.jl

@@ -2,14 +2,13 @@ module BipartiteGraphs
 
 export BipartiteEdge, BipartiteGraph
 
-using Reexport
 export 𝑠vertices, 𝑑vertices, has_𝑠vertex, has_𝑑vertex, 𝑠neighbors, 𝑑neighbors,
        𝑠edges, 𝑑edges, nsrcs, ndsts, SRC, DST
 
 using DocStringExtensions
 using UnPack
 using SparseArrays
-@reexport using LightGraphs
+using LightGraphs
 using Setfield
 
 ###

src/structural_transformation/StructuralTransformations.jl

@@ -21,6 +21,7 @@ using ModelingToolkit: ODESystem, AbstractSystem, var_from_nested_derivative,
                        ExtraVariablesSystemException
 
 using ModelingToolkit.BipartiteGraphs
+using LightGraphs
 using ModelingToolkit.SystemStructures
 
 using ModelingToolkit.DiffEqBase

src/systems/abstractsystem.jl

@@ -597,13 +597,16 @@ function Base.show(io::IO, ::MIME"text/plain", sys::AbstractSystem)
     return nothing
 end
 
-function _named(expr)
+function split_assign(expr)
     if !(expr isa Expr && expr.head === :(=) && expr.args[2].head === :call)
         throw(ArgumentError("expression should be of the form `sys = foo(a, b)`"))
     end
     name, call = expr.args
+end
 
+function _named(name, call, runtime=false)
     has_kw = false
+    call isa Expr || throw(Meta.ParseError("The rhs must be an Expr. Got $call."))
     if length(call.args) >= 2 && call.args[2] isa Expr
         # canonicalize to use `:parameters`
         if call.args[2].head === :kw
@@ -626,18 +629,78 @@ function _named(expr)
     kws = call.args[2].args
 
     if !any(kw->(kw isa Symbol ? kw : kw.args[1]) == :name, kws) # don't overwrite `name` kwarg
-        pushfirst!(kws, Expr(:kw, :name, Meta.quot(name)))
+        pushfirst!(kws, Expr(:kw, :name, runtime ? name : Meta.quot(name)))
+    end
+    call
+end
+
+function _named_idxs(name::Symbol, idxs, call)
+    if call.head !== :->
+        throw(ArgumentError("Not an anonymous function"))
+    end
+    if !isa(call.args[1], Symbol)
+        throw(ArgumentError("not a single-argument anonymous function"))
     end
-    :($name = $call)
+    sym, ex = call.args
+    ex = Base.Cartesian.poplinenum(ex)
+    ex = _named(:(Symbol($(Meta.quot(name)), :_, $sym)), ex, true)
+    ex = Base.Cartesian.poplinenum(ex)
+    :($name = $map($sym->$ex, $idxs))
 end
 
+check_name(name) = name isa Symbol || throw(Meta.ParseError("The lhs must be a symbol (a) or a ref (a[1:10]). Got $name."))
+
 """
-$(SIGNATURES)
+    @named y = foo(x)
+    @named y[1:10] = foo(x)
+    @named y 1:10 i -> foo(x*i)
 
 Rewrite `@named y = foo(x)` to `y = foo(x; name=:y)`.
+
+Rewrite `@named y[1:10] = foo(x)` to `y = map(i′->foo(x; name=Symbol(:y_, i′)), 1:10)`.
+
+Rewrite `@named y 1:10 i -> foo(x*i)` to `y = map(i->foo(x*i; name=Symbol(:y_, i)), 1:10)`.
+
+Examples:
+```julia
+julia> using ModelingToolkit
+
+julia> foo(i; name) = i, name
+foo (generic function with 1 method)
+
+julia> x = 41
+41
+
+julia> @named y = foo(x)
+(41, :y)
+
+julia> @named y[1:3] = foo(x)
+3-element Vector{Tuple{Int64, Symbol}}:
+ (41, :y_1)
+ (41, :y_2)
+ (41, :y_3)
+
+julia> @named y 1:3 i -> foo(x*i)
+3-element Vector{Tuple{Int64, Symbol}}:
+ (41, :y_1)
+ (82, :y_2)
+ (123, :y_3)
+```
 """
 macro named(expr)
-    esc(_named(expr))
+    name, call = split_assign(expr)
+    if Meta.isexpr(name, :ref)
+        name, idxs = name.args
+        check_name(name)
+        esc(_named_idxs(name, idxs, :($(gensym()) -> $call)))
+    else
+        check_name(name)
+        esc(:($name = $(_named(name, call))))
+    end
+end
+
+macro named(name::Symbol, idxs, call)
+    esc(_named_idxs(name, idxs, call))
 end
 
 function _config(expr, namespace)
@@ -826,7 +889,12 @@ by default.
 function extend(sys::AbstractSystem, basesys::AbstractSystem; name::Symbol=nameof(sys))
     T = SciMLBase.parameterless_type(basesys)
     iv = independent_variable(basesys)
-    sys = convert_system(T, sys, iv)
+    if iv === nothing
+        sys = convert_system(T, sys)
+    else
+        sys = convert_system(T, sys, iv)
+    end
+    
     eqs = union(equations(basesys), equations(sys))
     sts = union(states(basesys), states(sys))
     ps = union(parameters(basesys), parameters(sys))
@@ -849,15 +917,14 @@ Base.:(&)(sys::AbstractSystem, basesys::AbstractSystem; name::Symbol=nameof(sys)
 compose multiple systems together. The resulting system would inherit the first
 system's name.
 """
-compose(syss::AbstractSystem...; name=nameof(first(syss))) = compose(collect(syss); name=name)
-function compose(syss::AbstractArray{<:AbstractSystem}; name=nameof(first(syss)))
-    nsys = length(syss)
-    nsys >= 2 || throw(ArgumentError("There must be at least 2 systems. Got $nsys systems."))
-    sys = first(syss)
+function compose(sys::AbstractSystem, systems::AbstractArray{<:AbstractSystem}; name=nameof(sys))
+    nsys = length(systems)
+    nsys >= 1 || throw(ArgumentError("There must be at least 1 subsystem. Got $nsys subsystems."))
     @set! sys.name = name
-    @set! sys.systems = syss[2:end]
+    @set! sys.systems = systems
     return sys
 end
+compose(syss::AbstractSystem...; name=nameof(first(syss))) = compose(first(syss), collect(syss[2:end]); name=name)
 Base.:(∘)(sys1::AbstractSystem, sys2::AbstractSystem) = compose(sys1, sys2)
 
 UnPack.unpack(sys::ModelingToolkit.AbstractSystem, ::Val{p}) where p = getproperty(sys, p; namespace=false)