Merge pull request #748 from SciML/myb/default_params_states

Introduce default_ps and default_u0

Author: YingboMa

src/systems/abstractsystem.jl

@@ -297,3 +297,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

@@ -228,6 +228,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...)
@@ -254,29 +279,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
 
@@ -300,30 +304,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
@@ -358,19 +344,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
 
@@ -393,18 +369,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,21 +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)
+
+    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)

src/systems/diffeqs/sdesystem.jl

@@ -67,21 +67,37 @@ struct SDESystem <: AbstractODESystem
     Systems: the internal systems
     """
     systems::Vector{SDESystem}
+    """
+    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 SDESystem(deqs::AbstractVector{<:Equation}, neqs, iv, dvs, ps;
                    pins = [],
                    observed = [],
                    systems = SDESystem[],
+                   default_u0=Dict(),
+                   default_p=Dict(),
                    name = gensym(:SDESystem))
     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))
+
     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))
-    SDESystem(deqs, neqs, iv′, dvs′, ps′, pins, observed, tgrad, jac, Wfact, Wfact_t, name, systems)
+    SDESystem(deqs, neqs, iv′, dvs′, ps′, pins, observed, tgrad, jac, Wfact, Wfact_t, name, systems, default_u0, default_p)
 end
 
 function generate_diffusion_function(sys::SDESystem, dvs = sys.states, ps = sys.ps; kwargs...)
@@ -299,31 +315,11 @@ Generates an SDEProblem from an SDESystem and allows for automatically
 symbolically calculating numerical enhancements.
 """
 function DiffEqBase.SDEProblem{iip}(sys::SDESystem,u0map,tspan,parammap=DiffEqBase.NullParameters();
-                                    version = nothing, tgrad=false,
-                                    jac = false, Wfact = false,
-                                    checkbounds = false, sparse = false,
-                                    sparsenoise = sparse,
-                                    linenumbers = true, parallel=SerialForm(),
-                                    eval_expression = true,
+                                    sparsenoise = nothing,
                                     kwargs...) where iip
+    f, u0, p = process_DEProblem(SDEFunction{iip}, sys, u0map, parammap; kwargs...)
+    sparsenoise === nothing && (sparsenoise = get(kwargs, :sparse, false))
 
-    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 = SDEFunction{iip}(sys,dvs,ps,u0;tgrad=tgrad,jac=jac,Wfact=Wfact,
-                         checkbounds=checkbounds,
-                         linenumbers=linenumbers,parallel=parallel,
-                         sparse=sparse, eval_expression=eval_expression,kwargs...)
     if typeof(sys.noiseeqs) <: AbstractVector
         noise_rate_prototype = nothing
     elseif sparsenoise
@@ -358,29 +354,13 @@ numerical enhancements.
 struct SDEProblemExpr{iip} end
 
 function SDEProblemExpr{iip}(sys::SDESystem,u0map,tspan,
-                                    parammap=DiffEqBase.NullParameters();
-                                    version = nothing, tgrad=false,
-                                    jac = false, Wfact = false,
-                                    checkbounds = false, sparse = false,
-                                    linenumbers = false, parallel=SerialForm(),
-                                    kwargs...) where iip
-    dvs = states(sys)
-    ps = parameters(sys)
+                             parammap=DiffEqBase.NullParameters();
+                             sparsenoise = nothing,
+                             kwargs...) where iip
+    f, u0, p = process_DEProblem(SDEFunctionExpr{iip}, sys, u0map, parammap; kwargs...)
+    linenumbers = get(kwargs, :linenumbers, true)
+    sparsenoise === nothing && (sparsenoise = get(kwargs, :sparse, false))
 
-    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 = SDEFunctionExpr{iip}(sys,dvs,ps,u0;tgrad=tgrad,jac=jac,
-                        Wfact=Wfact,checkbounds=checkbounds,
-                        linenumbers=linenumbers,parallel=parallel,
-                        sparse=sparse,kwargs...)
     if typeof(sys.noiseeqs) <: AbstractVector
         noise_rate_prototype = nothing
     elseif sparsenoise