Merge branch 'master' into ms/observed_lowering

Author: YingboMa

Project.toml

@@ -1,7 +1,7 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Chris Rackauckas <[email protected]>"]
-version = "5.0.1"
+version = "5.1.0"
 
 [deps]
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
@@ -36,7 +36,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 [compat]
 ArrayInterface = "2.8"
 DataStructures = "0.17, 0.18"
-DiffEqBase = "6.48.1"
+DiffEqBase = "6.54.0"
 DiffEqJump = "6.7.5"
 DiffRules = "0.1, 1.0"
 Distributions = "0.24"
@@ -64,6 +64,7 @@ julia = "1.2"
 Dagger = "d58978e5-989f-55fb-8d15-ea34adc7bf54"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 GalacticOptim = "a75be94c-b780-496d-a8a9-0878b188d577"
+NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -72,4 +73,4 @@ StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Dagger", "ForwardDiff", "GalacticOptim", "OrdinaryDiffEq", "Optim", "Random", "SteadyStateDiffEq", "Test", "StochasticDiffEq"]
+test = ["Dagger", "ForwardDiff", "GalacticOptim", "NonlinearSolve", "OrdinaryDiffEq", "Optim", "Random", "SteadyStateDiffEq", "Test", "StochasticDiffEq"]

docs/src/tutorials/modelingtoolkitize.md

@@ -1,4 +1,4 @@
-# Symbolic Extensions to ODEProblem via Modelingtoolkize
+# Symbolic Extensions to ODEProblem via Modelingtoolkitize
 
 For some `DEProblem` types, automatic tracing functionality is already included
 via the `modelingtoolkitize` function. Take, for example, the Robertson ODE

docs/src/tutorials/symbolic_functions.md

@@ -287,7 +287,7 @@ just like any other function. To build a differential operator, use
 
 ```julia
 @variables t
-Differential(t)
+D = Differential(t)
 ```
 
 This is the differential operator ``D = \frac{\partial}{\partial t}``: the number of

src/systems/abstractsystem.jl

@@ -299,3 +299,6 @@ function (f::AbstractSysToExpr)(O)
     end
     return build_expr(:call, Any[operation(O); f.(arguments(O))])
 end
+
+get_default_p(sys) = sys.default_p
+get_default_u0(sys) = sys.default_u0

src/systems/control/controlsystem.jl

@@ -68,19 +68,33 @@ struct ControlSystem <: AbstractControlSystem
     systems: The internal systems
     """
     systems::Vector{ControlSystem}
+    """
+    default_u0: The default initial conditions to use when initial conditions
+    are not supplied in `ODEProblem`.
+    """
+    default_u0::Dict
+    """
+    default_p: The default parameters to use when parameters are not supplied
+    in `ODEProblem`.
+    """
+    default_p::Dict
 end
 
 function ControlSystem(loss, deqs::AbstractVector{<:Equation}, iv, dvs, controls, ps;
-                   pins = [],
-                   observed = [],
-                   systems = ODESystem[],
-                   name=gensym(:ControlSystem))
+                       pins = [],
+                       observed = [],
+                       systems = ODESystem[],
+                       default_u0=Dict(),
+                       default_p=Dict(),
+                       name=gensym(:ControlSystem))
     iv′ = value(iv)
     dvs′ = value.(dvs)
     controls′ = value.(controls)
     ps′ = value.(ps)
+    default_u0 isa Dict || (default_u0 = Dict(default_u0))
+    default_p isa Dict || (default_p = Dict(default_p))
     ControlSystem(value(loss), deqs, iv′, dvs′, controls′,
-                  ps′, pins, observed, name, systems)
+                  ps′, pins, observed, name, systems, default_u0, default_p)
 end
 
 struct ControlToExpr
@@ -102,7 +116,7 @@ end
 (f::ControlToExpr)(x::Sym) = x.name
 
 function constructRadauIIA5(T::Type = Float64)
-  sq6 = sqrt(6)
+  sq6 = sqrt(convert(T, 6))
   A = [11//45-7sq6/360 37//225-169sq6/1800 -2//225+sq6/75
        37//225+169sq6/1800 11//45+7sq6/360 -2//225-sq6/75
        4//9-sq6/36 4//9+sq6/36 1//9]
@@ -111,7 +125,7 @@ function constructRadauIIA5(T::Type = Float64)
   A = map(T,A)
   α = map(T,α)
   c = map(T,c)
-  return(DiffEqBase.ImplicitRKTableau(A,c,α,5))
+  return DiffEqBase.ImplicitRKTableau(A,c,α,5)
 end
 
 

src/systems/diffeqs/abstractodesystem.jl

@@ -232,6 +232,31 @@ function ODEFunctionExpr{iip}(sys::AbstractODESystem, dvs = states(sys),
     !linenumbers ? striplines(ex) : ex
 end
 
+function process_DEProblem(constructor, sys::AbstractODESystem,u0map,parammap;
+                           version = nothing, tgrad=false,
+                           jac = false,
+                           checkbounds = false, sparse = false,
+                           simplify=false,
+                           linenumbers = true, parallel=SerialForm(),
+                           eval_expression = true,
+                           kwargs...)
+    dvs = states(sys)
+    ps = parameters(sys)
+    u0map′ = lower_mapnames(u0map,sys.iv)
+    u0 = varmap_to_vars(u0map′,dvs; defaults=get_default_u0(sys))
+
+    if !(parammap isa DiffEqBase.NullParameters)
+        parammap′ = lower_mapnames(parammap)
+        p = varmap_to_vars(parammap′,ps; defaults=get_default_p(sys))
+    else
+        p = ps
+    end
+
+    f = constructor(sys,dvs,ps,u0;tgrad=tgrad,jac=jac,checkbounds=checkbounds,
+                    linenumbers=linenumbers,parallel=parallel,simplify=simplify,
+                    sparse=sparse,eval_expression=eval_expression,kwargs...)
+    return f, u0, p
+end
 
 function ODEFunctionExpr(sys::AbstractODESystem, args...; kwargs...)
     ODEFunctionExpr{true}(sys, args...; kwargs...)
@@ -258,29 +283,8 @@ Generates an ODEProblem from an ODESystem and allows for automatically
 symbolically calculating numerical enhancements.
 """
 function DiffEqBase.ODEProblem{iip}(sys::AbstractODESystem,u0map,tspan,
-                                    parammap=DiffEqBase.NullParameters();
-                                    version = nothing, tgrad=false,
-                                    jac = false,
-                                    checkbounds = false, sparse = false,
-                                    simplify=false,
-                                    linenumbers = true, parallel=SerialForm(),
-                                    eval_expression = true,
-                                    kwargs...) where iip
-    dvs = states(sys)
-    ps = parameters(sys)
-    u0map′ = lower_mapnames(u0map,sys.iv)
-    u0 = varmap_to_vars(u0map′,dvs)
-
-    if !(parammap isa DiffEqBase.NullParameters)
-        parammap′ = lower_mapnames(parammap)
-        p = varmap_to_vars(parammap′,ps)
-    else
-        p = ps
-    end
-
-    f = ODEFunction{iip}(sys,dvs,ps,u0;tgrad=tgrad,jac=jac,checkbounds=checkbounds,
-                        linenumbers=linenumbers,parallel=parallel,simplify=simplify,
-                        sparse=sparse,eval_expression=eval_expression,kwargs...)
+                                    parammap=DiffEqBase.NullParameters();kwargs...) where iip
+    f, u0, p = process_DEProblem(ODEFunction{iip}, sys, u0map, parammap; kwargs...)
     ODEProblem{iip}(f,u0,tspan,p;kwargs...)
 end
 
@@ -304,30 +308,12 @@ numerical enhancements.
 struct ODEProblemExpr{iip} end
 
 function ODEProblemExpr{iip}(sys::AbstractODESystem,u0map,tspan,
-                                    parammap=DiffEqBase.NullParameters();
-                                    version = nothing, tgrad=false,
-                                    jac = false,
-                                    checkbounds = false, sparse = false,
-                                    simplify=false,
-                                    linenumbers = false, parallel=SerialForm(),
-                                    kwargs...) where iip
-
-    dvs = states(sys)
-    ps = parameters(sys)
-    u0map′ = lower_mapnames(u0map,sys.iv)
-    u0 = varmap_to_vars(u0map′,dvs)
+                             parammap=DiffEqBase.NullParameters();
+                             kwargs...) where iip
 
-    if !(parammap isa DiffEqBase.NullParameters)
-        parammap′ = lower_mapnames(parammap)
-        p = varmap_to_vars(parammap′,ps)
-    else
-        p = ps
-    end
+    f, u0, p = process_DEProblem(ODEFunctionExpr{iip}, sys, u0map, parammap; kwargs...)
+    linenumbers = get(kwargs, :linenumbers, true)
 
-    f = ODEFunctionExpr{iip}(sys,dvs,ps,u0;tgrad=tgrad,jac=jac,checkbounds=checkbounds,
-                        linenumbers=linenumbers,parallel=parallel,
-                        simplify=simplify,
-                        sparse=sparse,kwargs...)
     ex = quote
         f = $f
         u0 = $u0
@@ -362,19 +348,9 @@ Generates an SteadyStateProblem from an ODESystem and allows for automatically
 symbolically calculating numerical enhancements.
 """
 function DiffEqBase.SteadyStateProblem{iip}(sys::AbstractODESystem,u0map,
-                                    parammap=DiffEqBase.NullParameters();
-                                    version = nothing, tgrad=false,
-                                    jac = false,
-                                    checkbounds = false, sparse = false,
-                                    linenumbers = true, parallel=SerialForm(),
-                                    kwargs...) where iip
-    dvs = states(sys)
-    ps = parameters(sys)
-    u0 = varmap_to_vars(u0map,dvs)
-    p = varmap_to_vars(parammap,ps)
-    f = ODEFunction(sys,dvs,ps,u0;tgrad=tgrad,jac=jac,checkbounds=checkbounds,
-                        linenumbers=linenumbers,parallel=parallel,
-                        sparse=sparse,kwargs...)
+                                            parammap=DiffEqBase.NullParameters();
+                                            kwargs...) where iip
+    f, u0, p = process_DEProblem(ODEFunction{iip}, sys, u0map, parammap; kwargs...)
     SteadyStateProblem(f,u0,p;kwargs...)
 end
 
@@ -397,18 +373,9 @@ struct SteadyStateProblemExpr{iip} end
 
 function SteadyStateProblemExpr{iip}(sys::AbstractODESystem,u0map,
                                     parammap=DiffEqBase.NullParameters();
-                                    version = nothing, tgrad=false,
-                                    jac = false,
-                                    checkbounds = false, sparse = false,
-                                    linenumbers = true, parallel=SerialForm(),
                                     kwargs...) where iip
-    dvs = states(sys)
-    ps = parameters(sys)
-    u0 = varmap_to_vars(u0map,dvs)
-    p = varmap_to_vars(parammap,ps)
-    f = ODEFunctionExpr(sys,dvs,ps,u0;tgrad=tgrad,jac=jac,checkbounds=checkbounds,
-                        linenumbers=linenumbers,parallel=parallel,
-                        sparse=sparse,kwargs...)
+    f, u0, p = process_DEProblem(ODEFunctionExpr{iip}, sys, u0map, parammap; kwargs...)
+    linenumbers = get(kwargs, :linenumbers, true)
     ex = quote
         f = $f
         u0 = $u0

src/systems/diffeqs/odesystem.jl

@@ -61,22 +61,39 @@ struct ODESystem <: AbstractODESystem
     systems: The internal systems
     """
     systems::Vector{ODESystem}
+    """
+    default_u0: The default initial conditions to use when initial conditions
+    are not supplied in `ODEProblem`.
+    """
+    default_u0::Dict
+    """
+    default_p: The default parameters to use when parameters are not supplied
+    in `ODEProblem`.
+    """
+    default_p::Dict
 end
 
-function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps;
+function ODESystem(
+                   deqs::AbstractVector{<:Equation}, iv, dvs, ps;
                    pins = Num[],
                    observed = Num[],
                    systems = ODESystem[],
-                   name=gensym(:ODESystem))
+                   name=gensym(:ODESystem),
+                   default_u0=Dict(),
+                   default_p=Dict(),
+                  )
     iv′ = value(iv)
     dvs′ = value.(dvs)
     ps′ = value.(ps)
-    pins′ = value.(pins)
+
+    default_u0 isa Dict || (default_u0 = Dict(default_u0))
+    default_p isa Dict || (default_p = Dict(default_p))
+
     tgrad = RefValue(Vector{Num}(undef, 0))
     jac = RefValue{Any}(Matrix{Num}(undef, 0, 0))
     Wfact   = RefValue(Matrix{Num}(undef, 0, 0))
     Wfact_t = RefValue(Matrix{Num}(undef, 0, 0))
-    ODESystem(deqs, iv′, dvs′, ps′, pins′, observed, tgrad, jac, Wfact, Wfact_t, name, systems)
+    ODESystem(deqs, iv′, dvs′, ps′, pins, observed, tgrad, jac, Wfact, Wfact_t, name, systems, default_u0, default_p)
 end
 
 var_from_nested_derivative(x, i=0) = (missing, missing)