Merge branch 'master' into id-tutorial

Author: iliailmer

Project.toml

@@ -1,7 +1,7 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Chris Rackauckas <[email protected]>"]
-version = "6.4.8"
+version = "6.6.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, 0.15"
+JuliaFormatter = "0.12, 0.13, 0.14, 0.15, 0.16, 0.17"
 LabelledArrays = "1.3"
 Latexify = "0.11, 0.12, 0.13, 0.14, 0.15"
 LightGraphs = "1.3"
@@ -69,11 +69,11 @@ Requires = "1.0"
 RuntimeGeneratedFunctions = "0.4.3, 0.5"
 SafeTestsets = "0.0.1"
 SciMLBase = "1.3"
-Setfield = "0.7"
+Setfield = "0.7, 0.8"
 SpecialFunctions = "0.7, 0.8, 0.9, 0.10, 1.0"
 StaticArrays = "0.10, 0.11, 0.12, 1.0"
-SymbolicUtils = "0.13.4"
-Symbolics = "3.1"
+SymbolicUtils = "0.16"
+Symbolics = "3.3.0"
 UnPack = "0.1, 1.0"
 Unitful = "1.1"
 julia = "1.2"

src/systems/abstractsystem.jl

@@ -397,7 +397,7 @@ end
 function parameters(sys::AbstractSystem)
     ps = get_ps(sys)
     systems = get_systems(sys)
-    isempty(systems) ? ps : [ps; reduce(vcat,namespace_parameters.(systems))]
+    unique(isempty(systems) ? ps : [ps; reduce(vcat,namespace_parameters.(systems))])
 end
 
 function controls(sys::AbstractSystem)
@@ -420,7 +420,13 @@ Base.@deprecate default_p(x) defaults(x) false
 function defaults(sys::AbstractSystem)
     systems = get_systems(sys)
     defs = get_defaults(sys)
-    isempty(systems) ? defs : mapreduce(namespace_defaults, merge, systems; init=defs)
+    # `mapfoldr` is really important!!! We should prefer the base model for
+    # defaults, because people write:
+    #
+    # `compose(ODESystem(...; defaults=defs), ...)`
+    #
+    # Thus, right associativity is required and crucial for correctness.
+    isempty(systems) ? defs : mapfoldr(namespace_defaults, merge, systems; init=defs)
 end
 
 states(sys::AbstractSystem, v) = renamespace(sys, v)
@@ -760,7 +766,7 @@ end
 function _config(expr, namespace)
     cn = Base.Fix2(_config, namespace)
     if Meta.isexpr(expr, :.)
-        return :($getvar($(map(cn, expr.args)...); namespace=$namespace))
+        return :($getproperty($(map(cn, expr.args)...); namespace=$namespace))
     elseif Meta.isexpr(expr, :function)
         def = splitdef(expr)
         def[:args] = map(cn, def[:args])

src/systems/alias_elimination.jl

@@ -3,7 +3,7 @@ using SymbolicUtils: Rewriters
 const KEEP = typemin(Int)
 
 function alias_elimination(sys)
-    sys = initialize_system_structure(sys)
+    sys = initialize_system_structure(sys; quick_cancel=true)
     s = structure(sys)
     is_linear_equations, eadj, cadj = find_linear_equations(sys)
 

src/systems/diffeqs/abstractodesystem.jl

@@ -261,7 +261,7 @@ function DiffEqBase.ODEFunction{iip}(sys::AbstractODESystem, dvs = states(sys),
 
     obs = observed(sys)
     observedfun = if steady_state
-        isempty(obs) ? SciMLBase.DEFAULT_OBSERVED_NO_TIME : let sys = sys, dict = Dict()
+        let sys = sys, dict = Dict()
             function generated_observed(obsvar, u, p, t=Inf)
                 obs = get!(dict, value(obsvar)) do
                     build_explicit_observed_function(sys, obsvar)
@@ -270,7 +270,7 @@ function DiffEqBase.ODEFunction{iip}(sys::AbstractODESystem, dvs = states(sys),
             end
         end
     else
-        isempty(obs) ? SciMLBase.DEFAULT_OBSERVED : let sys = sys, dict = Dict()
+        let sys = sys, dict = Dict()
             function generated_observed(obsvar, u, p, t)
                 obs = get!(dict, value(obsvar)) do
                     build_explicit_observed_function(sys, obsvar; checkbounds=checkbounds)
@@ -338,16 +338,7 @@ function DiffEqBase.DAEFunction{iip}(sys::AbstractODESystem, dvs = states(sys),
     # TODO: Jacobian sparsity / sparse Jacobian / dense Jacobian
 
     #=
-    observedfun = let sys = sys, dict = Dict()
         # TODO: We don't have enought information to reconstruct arbitrary state
-        # in general from `(u, p, t)`, e.g. `a ~ D(x)`.
-        function generated_observed(obsvar, u, p, t)
-            obs = get!(dict, value(obsvar)) do
-                build_explicit_observed_function(sys, obsvar)
-            end
-            obs(u, p, t)
-        end
-    end
     =#
 
     DAEFunction{iip}(
@@ -394,23 +385,6 @@ function ODEFunctionExpr{iip}(sys::AbstractODESystem, dvs = states(sys),
     f_oop, f_iip = generate_function(sys, dvs, ps; expression=Val{true}, kwargs...)
 
     dict = Dict()
-    #=
-    observedfun = if steady_state
-        :(function generated_observed(obsvar, u, p, t=Inf)
-              obs = get!($dict, value(obsvar)) do
-                  build_explicit_observed_function($sys, obsvar)
-              end
-              obs(u, p, t)
-          end)
-    else
-        :(function generated_observed(obsvar, u, p, t)
-              obs = get!($dict, value(obsvar)) do
-                  build_explicit_observed_function($sys, obsvar)
-              end
-              obs(u, p, t)
-          end)
-    end
-    =#
 
     fsym = gensym(:f)
     _f = :($fsym = $ODEFunctionClosure($f_oop, $f_iip))

src/systems/diffeqs/modelingtoolkitize.jl

@@ -51,7 +51,7 @@ function modelingtoolkitize(prob::DiffEqBase.ODEProblem; kwargs...)
 
     sts = vec(collect(vars))
 
-    params = if params isa Array && ndims(params) == 0
+    params = if params isa Number || (params isa Array && ndims(params) == 0)
         [params[1]]
     else
         vec(collect(params))

src/systems/diffeqs/odesystem.jl

@@ -19,7 +19,7 @@ eqs = [D(x) ~ σ*(y-x),
        D(y) ~ x*(ρ-z)-y,
        D(z) ~ x*y - β*z]
 
-de = ODESystem(eqs,t,[x,y,z],[σ,ρ,β])
+@named de = ODESystem(eqs,t,[x,y,z],[σ,ρ,β])
 ```
 """
 struct ODESystem <: AbstractODESystem
@@ -225,42 +225,52 @@ Build the observed function assuming the observed equations are all explicit,
 i.e. there are no cycles.
 """
 function build_explicit_observed_function(
-        sys, syms;
+        sys, ts;
         expression=false,
         output_type=Array,
         checkbounds=true)
 
-    if (isscalar = !(syms isa Vector))
-        syms = [syms]
+    if (isscalar = !(ts isa AbstractVector))
+        ts = [ts]
     end
-    syms = value.(syms)
+    ts = Symbolics.scalarize.(value.(ts))
+
+    vars = Set()
+    syms = foreach(Base.Fix1(vars!, vars), ts)
+    ivs = independent_variables(sys)
+    dep_vars = collect(setdiff(vars, ivs))
 
     obs = observed(sys)
+    sts = Set(states(sys))
     observed_idx = Dict(map(x->x.lhs, obs) .=> 1:length(obs))
-    output = similar(syms, Any)
-    # FIXME: this is a rather rough estimate of dependencies.
+
+    # FIXME: This is a rather rough estimate of dependencies. We assume
+    # the expression depends on everything before the `maxidx`.
     maxidx = 0
-    for (i, s) in enumerate(syms)
+    for (i, s) in enumerate(dep_vars)
         idx = get(observed_idx, s, nothing)
-        idx === nothing && throw(ArgumentError("$s is not an observed variable."))
+        if idx === nothing
+            if !(s in sts)
+                throw(ArgumentError("$s is either an observed nor a state variable."))
+            end
+            continue
+        end
         idx > maxidx && (maxidx = idx)
-        output[i] = obs[idx].rhs
     end
+    obsexprs = map(eq -> eq.lhs←eq.rhs, obs[1:maxidx])
 
     dvs = DestructuredArgs(states(sys), inbounds=!checkbounds)
     ps = DestructuredArgs(parameters(sys), inbounds=!checkbounds)
-    ivs = independent_variables(sys)
     args = [dvs, ps, ivs...]
     pre = get_postprocess_fbody(sys)
 
     ex = Func(
         args, [],
         pre(Let(
-            map(eq -> eq.lhs←eq.rhs, obs[1:maxidx]),
-            isscalar ? output[1] : MakeArray(output, output_type)
+            obsexprs,
+            isscalar ? ts[1] : MakeArray(ts, output_type)
            ))
     ) |> toexpr
-
     expression ? ex : @RuntimeGeneratedFunction(ex)
 end
 

src/systems/diffeqs/sdesystem.jl

@@ -23,7 +23,7 @@ noiseeqs = [0.1*x,
             0.1*y,
             0.1*z]
 
-de = SDESystem(eqs,noiseeqs,t,[x,y,z],[σ,ρ,β])
+@named de = SDESystem(eqs,noiseeqs,t,[x,y,z],[σ,ρ,β])
 ```
 """
 struct SDESystem <: AbstractODESystem

src/systems/discrete_system/discrete_system.jl

@@ -11,14 +11,16 @@ $(FIELDS)
 ```
 using ModelingToolkit
 
-@parameters σ ρ β
-@variables t x(t) y(t) z(t) next_x(t) next_y(t) next_z(t)
+@parameters σ=28.0 ρ=10.0 β=8/3 δt=0.1
+@variables t x(t)=1.0 y(t)=0.0 z(t)=0.0
+D = Difference(t; dt=δt)
 
-eqs = [next_x ~ σ*(y-x),
-       next_y ~ x*(ρ-z)-y,
-       next_z ~ x*y - β*z]
+eqs = [D(x) ~ σ*(y-x),
+       D(y) ~ x*(ρ-z)-y,
+       D(z) ~ x*y - β*z]
 
-de = DiscreteSystem(eqs,t,[x,y,z],[σ,ρ,β])
+@named de = DiscreteSystem(eqs,t,[x,y,z],[σ,ρ,β]) # or 
+@named de = DiscreteSystem(eqs)
 ```
 """
 struct DiscreteSystem <: AbstractTimeDependentSystem
@@ -45,22 +47,21 @@ struct DiscreteSystem <: AbstractTimeDependentSystem
     """
     systems::Vector{DiscreteSystem}
     """
-    default_u0: The default initial conditions to use when initial conditions
-    are not supplied in `DiscreteSystem`.
+    defaults: The default values to use when initial conditions and/or
+    parameters are not supplied in `DiscreteProblem`.
     """
-    default_u0::Dict
+    defaults::Dict
     """
-    default_p: The default parameters to use when parameters are not supplied
-    in `DiscreteSystem`.
+    type: type of the system
     """
-    default_p::Dict
-    function DiscreteSystem(discreteEqs, iv, dvs, ps, var_to_name, ctrls, observed, name, systems, default_u0, default_p; checks::Bool = true)
+    connection_type::Any
+    function DiscreteSystem(discreteEqs, iv, dvs, ps, var_to_name, ctrls, observed, name, systems, defaults, connection_type; checks::Bool = true)
         if checks
             check_variables(dvs, iv)
             check_parameters(ps, iv)
             all_dimensionless([dvs;ps;iv;ctrls]) ||check_units(discreteEqs)
         end
-        new(discreteEqs, iv, dvs, ps, var_to_name, ctrls, observed, name, systems, default_u0, default_p)
+        new(discreteEqs, iv, dvs, ps, var_to_name, ctrls, observed, name, systems, defaults, connection_type)
     end
 end
 
@@ -78,6 +79,7 @@ function DiscreteSystem(
                    default_u0=Dict(),
                    default_p=Dict(),
                    defaults=_merge(Dict(default_u0), Dict(default_p)),
+                   connection_type=nothing,
                    kwargs...,
                   )
     name === nothing && throw(ArgumentError("The `name` keyword must be provided. Please consider using the `@named` macro"))
@@ -88,7 +90,7 @@ function DiscreteSystem(
     ctrl′ = value.(controls)
 
     if !(isempty(default_u0) && isempty(default_p))
-        Base.depwarn("`default_u0` and `default_p` are deprecated. Use `defaults` instead.", :ODESystem, force=true)
+        Base.depwarn("`default_u0` and `default_p` are deprecated. Use `defaults` instead.", :DiscreteSystem, force=true)
     end
     defaults = todict(defaults)
     defaults = Dict(value(k) => value(v) for (k, v) in pairs(defaults))
@@ -101,7 +103,46 @@ function DiscreteSystem(
     if length(unique(sysnames)) != length(sysnames)
         throw(ArgumentError("System names must be unique."))
     end
-    DiscreteSystem(eqs, iv′, dvs′, ps′, var_to_name, ctrl′, observed, name, systems, default_u0, default_p, kwargs...)
+    DiscreteSystem(eqs, iv′, dvs′, ps′, var_to_name, ctrl′, observed, name, systems, defaults, connection_type, kwargs...)
+end
+
+
+function DiscreteSystem(eqs, iv=nothing; kwargs...)
+    eqs = collect(eqs)
+    # NOTE: this assumes that the order of algebric equations doesn't matter
+    diffvars = OrderedSet()
+    allstates = OrderedSet()
+    ps = OrderedSet()
+    # reorder equations such that it is in the form of `diffeq, algeeq`
+    diffeq = Equation[]
+    algeeq = Equation[]
+    # initial loop for finding `iv`
+    if iv === nothing
+        for eq in eqs
+            if !(eq.lhs isa Number) # assume eq.lhs is either Differential or Number
+                iv = iv_from_nested_difference(eq.lhs)
+                break
+            end
+        end
+    end
+    iv = value(iv)
+    iv === nothing && throw(ArgumentError("Please pass in independent variables."))
+    for eq in eqs
+        collect_vars_difference!(allstates, ps, eq.lhs, iv)
+        collect_vars_difference!(allstates, ps, eq.rhs, iv)
+        if isdifferenceeq(eq)
+            diffvar, _ = var_from_nested_difference(eq.lhs)
+            isequal(iv, iv_from_nested_difference(eq.lhs)) || throw(ArgumentError("A DiscreteSystem can only have one independent variable."))
+            diffvar in diffvars && throw(ArgumentError("The difference variable $diffvar is not unique in the system of equations."))
+            push!(diffvars, diffvar)
+            push!(diffeq, eq)
+        else
+            push!(algeeq, eq)
+        end
+    end
+    algevars = setdiff(allstates, diffvars)
+    # the orders here are very important!
+    return DiscreteSystem(append!(diffeq, algeeq), iv, vcat(collect(diffvars), collect(algevars)), ps; kwargs...)
 end
 
 """
@@ -118,16 +159,37 @@ function DiffEqBase.DiscreteProblem(sys::DiscreteSystem,u0map,tspan,
     ps = parameters(sys)
     eqs = equations(sys)
     eqs = linearize_eqs(sys, eqs)
-    # defs = defaults(sys)
-    t = get_iv(sys)
-    u0 = varmap_to_vars(u0map,dvs)
+    defs = defaults(sys)
+    iv = get_iv(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=u0defs)
+    
     rhss = [eq.rhs for eq in eqs]
     u = dvs
-    p = varmap_to_vars(parammap,ps)
+    p = varmap_to_vars(parammap,ps; defaults=pdefs)
 
     f_gen = generate_function(sys; expression=Val{eval_expression}, expression_module=eval_module)
     f_oop, _ = (@RuntimeGeneratedFunction(eval_module, ex) for ex in f_gen)
-    f(u,p,t) = f_oop(u,p,t)
+    f(u,p,iv) = f_oop(u,p,iv)
     DiscreteProblem(f,u0,tspan,p;kwargs...)
 end