Merge branch 'master' of github.com:SciML/ModelingToolkit.jl into domain-sets

Author: valentinsulzer

Project.toml

@@ -1,7 +1,7 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Chris Rackauckas <[email protected]>"]
-version = "5.15.0"
+version = "5.17.1"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -16,6 +16,7 @@ Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
 IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
+JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
 LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800"
 Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
 Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
@@ -48,9 +49,10 @@ DataStructures = "0.17, 0.18"
 DiffEqBase = "6.54.0"
 DiffEqJump = "6.7.5"
 DiffRules = "0.1, 1.0"
-Distributions = "0.23, 0.24"
+Distributions = "0.23, 0.24, 0.25"
 DocStringExtensions = "0.7, 0.8"
 IfElse = "0.1"
+JuliaFormatter = "0.12, 0.13"
 LabelledArrays = "1.3"
 Latexify = "0.11, 0.12, 0.13, 0.14, 0.15"
 LightGraphs = "1.3"

docs/src/basics/Composition.md

@@ -125,6 +125,40 @@ done even if the variable `x` is eliminated from the system from
 transformations like `alias_elimination` or `tearing`: the variable
 will be lazily reconstructed on demand.
 
+### Variable scope and parameter expressions
+
+In some scenarios, it could be useful for model parameters to be expressed
+in terms of other parameters, or shared between common subsystems.
+To fascilitate this, ModelingToolkit supports sybmolic expressions
+in default values, and scoped variables.
+
+With symbolic parameters, it is possible to set the default value of a parameter or initial condition to an expression of other variables.
+
+```julia
+# ...
+sys = ODESystem(
+    # ...
+    # directly in the defauls argument
+    defaults=Pair{Num, Any}[
+    x => u,
+    y => σ,
+    z => u-0.1,
+])
+# by assigning to the parameter
+sys.y = u*1.1
+```
+
+In a hierarchical system, variables of the subsystem get namespaced by the name of the system they are in. This prevents naming clashes, but also enforces that every state and parameter is local to the subsystem it is used in. In some cases it might be desirable to have variables and parameters that are shared between subsystems, or even global. This can be accomplished as follows.
+
+```julia
+@variables a b c d
+
+# a is a local variable
+b = ParentScope(b) # b is a variable that belongs to one level up in the hierarchy
+c = ParentScope(ParentScope(c)) # ParentScope can be nested
+d = GlobalScope(d) # global variables will never be namespaced
+```
+
 ## Structural Simplify
 
 In many cases, the nicest way to build a model may leave a lot of

docs/src/tutorials/ode_modeling.md

@@ -119,8 +119,8 @@ equations(fol_simplified) == equations(fol_model)
 ```
 
 You can extract the equations from a system using `equations` (and, in the same
-way, `variables` and `parameters`). The simplified equation is exactly the same
-as the original one, so the simulation performence will also be the same.
+way, `states` and `parameters`). The simplified equation is exactly the same
+as the original one, so the simulation performance will also be the same.
 However, there is one difference. MTK does keep track of the eliminated
 algebraic variables as "observables" (see
 [Observables and Variable Elimination](@ref)).

src/ModelingToolkit.jl

@@ -35,6 +35,7 @@ import SymbolicUtils: istree, arguments, operation, similarterm, promote_symtype
 using SymbolicUtils.Code
 import SymbolicUtils.Code: toexpr
 import SymbolicUtils.Rewriters: Chain, Postwalk, Prewalk, Fixpoint
+import JuliaFormatter
 
 using Reexport
 @reexport using Symbolics
@@ -129,14 +130,16 @@ include("systems/pde/pdesystem.jl")
 include("systems/reaction/reactionsystem.jl")
 include("systems/dependency_graphs.jl")
 
+include("systems/discrete_system/discrete_system.jl")
+
 include("systems/systemstructure.jl")
 using .SystemStructures
 
 include("systems/alias_elimination.jl")
 include("structural_transformation/StructuralTransformations.jl")
 @reexport using .StructuralTransformations
 
-export ODESystem, ODEFunction, ODEFunctionExpr, ODEProblemExpr
+export ODESystem, ODEFunction, ODEFunctionExpr, ODEProblemExpr, convert_system
 export DAEFunctionExpr, DAEProblemExpr
 export SDESystem, SDEFunction, SDEFunctionExpr, SDESystemExpr
 export SystemStructure
@@ -161,6 +164,7 @@ export Term, Sym
 export SymScope, LocalScope, ParentScope, GlobalScope
 export independent_variable, states, parameters, equations, controls, observed, structure
 export structural_simplify
+export DiscreteSystem, DiscreteProblem
 
 export calculate_jacobian, generate_jacobian, generate_function
 export calculate_tgrad, generate_tgrad

src/structural_transformation/codegen.jl

@@ -68,7 +68,8 @@ function torn_system_jacobian_sparsity(sys)
         end
     end
 
-    dvar2idx(idx) = idx # maps `dvar` to the index of the states
+    dvrange = diffvars_range(s)
+    dvar2idx = Dict(v=>i for (i, v) in enumerate(dvrange))
     I = Int[]; J = Int[]
     eqidx = 0
     for ieq in 𝑠vertices(graph)
@@ -77,11 +78,11 @@ function torn_system_jacobian_sparsity(sys)
         for ivar in 𝑠neighbors(graph, ieq)
             if isdiffvar(s, ivar)
                 push!(I, eqidx)
-                push!(J, dvar2idx(ivar))
+                push!(J, dvar2idx[ivar])
             elseif isalgvar(s, ivar)
                 for dvar in avars2dvars[ivar]
                     push!(I, eqidx)
-                    push!(J, dvar2idx(dvar))
+                    push!(J, dvar2idx[dvar])
                 end
             end
         end

src/systems/abstractsystem.jl

@@ -411,6 +411,103 @@ function (f::AbstractSysToExpr)(O)
     return build_expr(:call, Any[operation(O); f.(arguments(O))])
 end
 
+###
+### System utils
+###
+function push_vars!(stmt, name, typ, vars)
+    isempty(vars) && return
+    vars_expr = Expr(:macrocall, typ, nothing)
+    for s in vars
+        if istree(s)
+            f = nameof(operation(s))
+            args = arguments(s)
+            ex = :($f($(args...)))
+        else
+            ex = nameof(s)
+        end
+        push!(vars_expr.args, ex)
+    end
+    push!(stmt, :($name = $collect($vars_expr)))
+    return
+end
+
+function round_trip_expr(t, var2name)
+    name = get(var2name, t, nothing)
+    name !== nothing && return name
+    t isa Sym && return nameof(t)
+    istree(t) || return t
+    f = round_trip_expr(operation(t), var2name)
+    args = map(Base.Fix2(round_trip_expr, var2name), arguments(t))
+    return :($f($(args...)))
+end
+round_trip_eq(eq, var2name) = Expr(:call, :~, round_trip_expr(eq.lhs, var2name), round_trip_expr(eq.rhs, var2name))
+
+function push_eqs!(stmt, eqs, var2name)
+    eqs_name = gensym(:eqs)
+    eqs_expr = Expr(:vcat)
+    eqs_blk = Expr(:(=), eqs_name, eqs_expr)
+    for eq in eqs
+        push!(eqs_expr.args, round_trip_eq(eq, var2name))
+    end
+
+    push!(stmt, eqs_blk)
+    return eqs_name
+end
+
+function push_defaults!(stmt, defs, var2name)
+    defs_name = gensym(:defs)
+    defs_expr = Expr(:call, Dict)
+    defs_blk = Expr(:(=), defs_name, defs_expr)
+    for d in defs
+        n = round_trip_expr(d.first, var2name)
+        v = round_trip_expr(d.second, var2name)
+        push!(defs_expr.args, :($(=>)($n, $v)))
+    end
+
+    push!(stmt, defs_blk)
+    return defs_name
+end
+
+###
+### System I/O
+###
+function toexpr(sys::AbstractSystem)
+    sys = flatten(sys)
+    expr = Expr(:block)
+    stmt = expr.args
+
+    iv = independent_variable(sys)
+    ivname = gensym(:iv)
+    if iv !== nothing
+        push!(stmt, :($ivname = (@variables $(getname(iv)))[1]))
+    end
+
+    stsname = gensym(:sts)
+    sts = states(sys)
+    push_vars!(stmt, stsname, Symbol("@variables"), sts)
+    psname = gensym(:ps)
+    ps = parameters(sys)
+    push_vars!(stmt, psname, Symbol("@parameters"), ps)
+
+    var2name = Dict{Any,Symbol}()
+    for v in Iterators.flatten((sts, ps))
+        var2name[v] = getname(v)
+    end
+
+    eqs_name = push_eqs!(stmt, equations(sys), var2name)
+    defs_name = push_defaults!(stmt, defaults(sys), var2name)
+
+    if sys isa ODESystem
+        push!(stmt, :($ODESystem($eqs_name, $ivname, $stsname, $psname; defaults=$defs_name)))
+    elseif sys isa NonlinearSystem
+        push!(stmt, :($NonlinearSystem($eqs_name, $stsname, $psname; defaults=$defs_name)))
+    end
+
+    striplines(expr) # keeping the line numbers is never helpful
+end
+
+Base.write(io::IO, sys::AbstractSystem) = write(io, readable_code(toexpr(sys)))
+
 function Base.show(io::IO, ::MIME"text/plain", sys::AbstractSystem)
     eqs = equations(sys)
     if eqs isa AbstractArray
@@ -581,6 +678,10 @@ function check_eqs_u0(eqs, dvs, u0)
     return nothing
 end
 
+###
+### Connectors
+###
+
 function with_connection_type(expr)
     @assert expr isa Expr && (expr.head == :function || (expr.head == :(=) &&
                                        expr.args[1] isa Expr &&

src/systems/diffeqs/abstractodesystem.jl

@@ -221,7 +221,7 @@ function DiffEqBase.ODEFunction{iip}(sys::AbstractODESystem, dvs = states(sys),
                      jac = _jac === nothing ? nothing : _jac,
                      tgrad = _tgrad === nothing ? nothing : _tgrad,
                      mass_matrix = _M,
-                     jac_prototype = (!isnothing(u0) && sparse) ? (!jac ? similar(jacobian_sparsity(sys),uElType) : similar(get_jac(sys)[],uElType)) : nothing,
+                     jac_prototype = (!isnothing(u0) && sparse) ? (!jac ? similar(jacobian_sparsity(sys),uElType) : SparseArrays.sparse(similar(get_jac(sys)[],uElType))) : nothing,
                      syms = Symbol.(states(sys)),
                      indepsym = Symbol(independent_variable(sys)),
                      observed = observedfun,
@@ -397,16 +397,34 @@ function process_DEProblem(constructor, sys::AbstractODESystem,u0map,parammap;
     ps = parameters(sys)
     defs = defaults(sys)
     iv = independent_variable(sys)
+    if parammap isa Dict
+        u0defs = merge(parammap, defs)
+    elseif eltype(parammap) <: Pair
+        u0defs = merge(Dict(parammap), defs)
+    elseif eltype(parammap) <: Number
+        u0defs = merge(Dict(zip(ps, parammap)), defs)
+    else
+        u0defs = defs
+    end
+    if u0map isa Dict
+        pdefs = merge(u0map, defs)
+    elseif eltype(u0map) <: Pair
+        pdefs = merge(Dict(u0map), defs)
+    elseif eltype(u0map) <: Number
+        pdefs = merge(Dict(zip(dvs, u0map)), defs)
+    else
+        pdefs = defs
+    end
 
-    u0 = varmap_to_vars(u0map,dvs; defaults=defs)
+    u0 = varmap_to_vars(u0map,dvs; defaults=u0defs)
     if implicit_dae && du0map !== nothing
         ddvs = map(Differential(iv), dvs)
         du0 = varmap_to_vars(du0map, ddvs; defaults=defaults, toterm=identity)
     else
         du0 = nothing
         ddvs = nothing
     end
-    p = varmap_to_vars(parammap,ps; defaults=defs)
+    p = varmap_to_vars(parammap,ps; defaults=pdefs)
 
     check_eqs_u0(eqs, dvs, u0)
 

src/systems/diffeqs/modelingtoolkitize.jl

@@ -17,8 +17,15 @@ function modelingtoolkitize(prob::DiffEqBase.ODEProblem)
     has_p = !(p isa Union{DiffEqBase.NullParameters,Nothing})
 
     var(x, i) = Num(Sym{FnType{Tuple{symtype(t)}, Real}}(nameof(Variable(x, i))))
-    vars = ArrayInterface.restructure(prob.u0,[var(:x, i)(ModelingToolkit.value(t)) for i in eachindex(prob.u0)])
-    params = has_p ? reshape([Num(toparam(Sym{Real}(nameof(Variable(:α, i))))) for i in eachindex(p)],size(p)) : []
+    _vars = [var(:x, i)(ModelingToolkit.value(t)) for i in eachindex(prob.u0)]
+    vars = prob.u0 isa Number ? _vars : ArrayInterface.restructure(prob.u0,_vars)
+    params = if has_p
+        _params = [Num(toparam(Sym{Real}(nameof(Variable(:α, i))))) for i in eachindex(p)]
+        p isa Number ? _params[1] : reshape(_params,size(p))
+    else
+        []
+    end
+
     var_set = Set(vars)
 
     D = Differential(t)
@@ -64,8 +71,6 @@ function modelingtoolkitize(prob::DiffEqBase.ODEProblem)
     de
 end
 
-
-
 """
 $(TYPEDSIGNATURES)
 

src/systems/diffeqs/odesystem.jl

@@ -300,3 +300,42 @@ function collect_differential_variables(sys::ODESystem)
     end
     return diffvars
 end
+
+# We have a stand-alone function to convert a `NonlinearSystem` or `ODESystem`
+# to an `ODESystem` to connect systems, and we later can reply on
+# `structural_simplify` to convert `ODESystem`s to `NonlinearSystem`s.
+"""
+$(TYPEDSIGNATURES)
+
+Convert a `NonlinearSystem` to an `ODESystem` or converts an `ODESystem` to a
+new `ODESystem` with a different independent variable.
+"""
+function convert_system(::Type{<:ODESystem}, sys, t; name=nameof(sys))
+    isempty(observed(sys)) || throw(ArgumentError("`convert_system` cannot handle reduced model (i.e. observed(sys) is non-empty)."))
+    t = value(t)
+    varmap = Dict()
+    sts = states(sys)
+    newsts = similar(sts, Any)
+    for (i, s) in enumerate(sts)
+        if istree(s)
+            args = arguments(s)
+            length(args) == 1 || throw(InvalidSystemException("Illegal state: $s. The state can have at most one argument like `x(t)`."))
+            arg = args[1]
+            if isequal(arg, t)
+                newsts[i] = s
+                continue
+            end
+            ns = operation(s)(t)
+            newsts[i] = ns
+            varmap[s] = ns
+        else
+            ns = indepvar2depvar(s, t)
+            newsts[i] = ns
+            varmap[s] = ns
+        end
+    end
+    sub = Base.Fix2(substitute, varmap)
+    neweqs = map(sub, equations(sys))
+    defs = Dict(sub(k) => sub(v) for (k, v) in defaults(sys))
+    return ODESystem(neweqs, t, newsts, parameters(sys); defaults=defs, name=name)
+end