Merge remote-tracking branch 'origin/master' into s/array-fixes

Author: shashi

Project.toml

@@ -1,7 +1,7 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Chris Rackauckas <[email protected]>"]
-version = "5.16.0"
+version = "5.17.3"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -48,10 +48,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.13"
+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

@@ -130,6 +130,8 @@ 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
 
@@ -162,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/systems/abstractsystem.jl

@@ -665,15 +665,17 @@ AbstractTrees.children(sys::ModelingToolkit.AbstractSystem) = ModelingToolkit.ge
 AbstractTrees.printnode(io::IO, sys::ModelingToolkit.AbstractSystem) = print(io, nameof(sys))
 AbstractTrees.nodetype(::ModelingToolkit.AbstractSystem) = ModelingToolkit.AbstractSystem
 
-function check_eqs_u0(eqs, dvs, u0)
+function check_eqs_u0(eqs, dvs, u0; check_length=true, kwargs...)
     if u0 !== nothing
-        if !(length(eqs) == length(dvs) == length(u0))
-            throw(ArgumentError("Equations ($(length(eqs))), states ($(length(dvs))), and initial conditions ($(length(u0))) are of different lengths."))
-        end
-    else
-        if !(length(eqs) == length(dvs))
-            throw(ArgumentError("Equations ($(length(eqs))), states ($(length(dvs))) are of different lengths."))
+        if check_length
+            if !(length(eqs) == length(dvs) == length(u0))
+                throw(ArgumentError("Equations ($(length(eqs))), states ($(length(dvs))), and initial conditions ($(length(u0))) are of different lengths. To allow a different number of equations than states use kwarg check_length=false."))
+            end
+        elseif length(dvs) != length(u0)
+            throw(ArgumentError("States ($(length(dvs))) and initial conditions ($(length(u0))) are of different lengths."))
         end
+    elseif check_length && (length(eqs) != length(dvs))
+        throw(ArgumentError("Equations ($(length(eqs))) and states ($(length(dvs))) are of different lengths. To allow these to differ use kwarg check_length=false."))
     end
     return nothing
 end
@@ -713,7 +715,7 @@ end
 
 promote_connect_rule(::Type{T}, ::Type{S}) where {T, S} = Union{}
 promote_connect_rule(::Type{T}, ::Type{T}) where {T} = T
-promote_connect_type(t1::Type, t2::Type, ts::Type...) = promote_connect_rule(promote_connect_rule(t1, t2), ts...)
+promote_connect_type(t1::Type, t2::Type, ts::Type...) = promote_connect_type(promote_connect_rule(t1, t2), ts...)
 @inline function promote_connect_type(::Type{T}, ::Type{S}) where {T,S}
     promote_connect_result(
         T,

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,18 +397,36 @@ 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)
+    check_eqs_u0(eqs, dvs, u0; kwargs...)
 
     f = constructor(sys,dvs,ps,u0;ddvs=ddvs,tgrad=tgrad,jac=jac,checkbounds=checkbounds,
                     linenumbers=linenumbers,parallel=parallel,simplify=simplify,

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/discrete_system/discrete_system.jl

@@ -0,0 +1,107 @@
+"""
+$(TYPEDEF)
+
+A system of difference equations.
+
+# Fields
+$(FIELDS)
+
+# Example
+
+```
+using ModelingToolkit
+
+@parameters t σ ρ β
+@variables x(t) y(t) z(t) next_x(t) next_y(t) next_z(t)
+
+eqs = [next_x ~ σ*(y-x),
+       next_y ~ x*(ρ-z)-y,
+       next_z ~ x*y - β*z]
+
+de = DiscreteSystem(eqs,t,[x,y,z],[σ,ρ,β])
+```
+"""
+struct DiscreteSystem <: AbstractSystem
+    """The differential equations defining the discrete system."""
+    eqs::Vector{Equation}
+    """Independent variable."""
+    iv::Sym
+    """Dependent (state) variables."""
+    states::Vector
+    """Parameter variables."""
+    ps::Vector
+    observed::Vector{Equation}
+    """
+    Name: the name of the system
+    """
+    name::Symbol
+    """
+    systems: The internal systems. These are required to have unique names.
+    """
+    systems::Vector{DiscreteSystem}
+    """
+    default_u0: The default initial conditions to use when initial conditions
+    are not supplied in `DiscreteSystem`.
+    """
+    default_u0::Dict
+    """
+    default_p: The default parameters to use when parameters are not supplied
+    in `DiscreteSystem`.
+    """
+    default_p::Dict
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Constructs a DiscreteSystem.
+"""
+function DiscreteSystem(
+                   discreteEqs::AbstractVector{<:Equation}, iv, dvs, ps;
+                   observed = Num[],
+                   systems = DiscreteSystem[],
+                   name=gensym(:DiscreteSystem),
+                   default_u0=Dict(),
+                   default_p=Dict(),
+                  )
+    iv′ = value(iv)
+    dvs′ = value.(dvs)
+    ps′ = value.(ps)
+
+    default_u0 isa Dict || (default_u0 = Dict(default_u0))
+    default_p isa Dict || (default_p = Dict(default_p))
+    default_u0 = Dict(value(k) => value(default_u0[k]) for k in keys(default_u0))
+    default_p = Dict(value(k) => value(default_p[k]) for k in keys(default_p))
+
+    sysnames = nameof.(systems)
+    if length(unique(sysnames)) != length(sysnames)
+        throw(ArgumentError("System names must be unique."))
+    end
+    DiscreteSystem(discreteEqs, iv′, dvs′, ps′, observed, name, systems, default_u0, default_p)
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Generates an DiscreteProblem from an DiscreteSystem.
+"""
+function DiffEqBase.DiscreteProblem(sys::DiscreteSystem,u0map,tspan,
+                                    parammap=DiffEqBase.NullParameters();
+                                    eval_module = @__MODULE__,
+                                    eval_expression = true,
+                                    kwargs...)
+    dvs = states(sys)
+    ps = parameters(sys)
+    eqs = equations(sys)
+    # defs = defaults(sys)
+    t = get_iv(sys)
+    u0 = varmap_to_vars(u0map,dvs)
+    rhss = [eq.rhs for eq in eqs]
+    u = dvs
+    p = varmap_to_vars(parammap,ps)
+
+    f_gen = build_function(rhss, dvs, ps, t; expression=Val{eval_expression}, expression_module=eval_module)
+    f_oop,f_iip = (@RuntimeGeneratedFunction(eval_module, ex) for ex in f_gen)
+    f(u,p,t) = f_oop(u,p,t)
+    DiscreteProblem(f,u0,tspan,p;kwargs...)
+end

src/systems/nonlinear/nonlinearsystem.jl

@@ -220,7 +220,7 @@ function process_NonlinearProblem(constructor, sys::NonlinearSystem,u0map,paramm
     u0 = varmap_to_vars(u0map,dvs; defaults=defs)
     p = varmap_to_vars(parammap,ps; defaults=defs)
 
-    check_eqs_u0(eqs, dvs, u0)
+    check_eqs_u0(eqs, dvs, u0; kwargs...)
 
     f = constructor(sys,dvs,ps,u0;jac=jac,checkbounds=checkbounds,
                     linenumbers=linenumbers,parallel=parallel,simplify=simplify,