allow u0 and p defaults to refer to each other

Author: pepijndevos

src/structural_transformation/codegen.jl

@@ -338,12 +338,12 @@ function ODAEProblem{iip}(
     s = structure(sys)
     @unpack fullvars = s
     dvs = fullvars[diffvars_range(s)]
+    defaults = merge(default_p(sys), default_u0(sys))
     u0map′ = ModelingToolkit.lower_mapnames(u0map, independent_variable(sys))
-    u0 = ModelingToolkit.varmap_to_vars(u0map′, dvs; defaults=default_u0(sys))
+    u0 = ModelingToolkit.varmap_to_vars(u0map′, dvs; defaults=defaults)
 
     ps = parameters(sys)
-    d_p = default_p(sys)
-    if parammap isa DiffEqBase.NullParameters && isempty(d_p)
+    if parammap isa DiffEqBase.NullParameters && isempty(default_p(sys))
         isempty(ps) || throw(ArgumentError("The model has non-empty parameters but no parameters are specified in the problem."))
         p = parammap
     else
@@ -352,7 +352,7 @@ function ODAEProblem{iip}(
         else
             pp = ModelingToolkit.lower_mapnames(parammap)
         end
-        p = ModelingToolkit.varmap_to_vars(pp, ps; defaults=d_p)
+        p = ModelingToolkit.varmap_to_vars(pp, ps; defaults=defaults)
     end
 
     ODEProblem{iip}(build_torn_function(sys; kw...), u0, tspan, p; kw...)

src/systems/diffeqs/abstractodesystem.jl

@@ -273,20 +273,20 @@ function process_DEProblem(constructor, sys::AbstractODESystem,u0map,parammap;
                            kwargs...)
     dvs = states(sys)
     ps = parameters(sys)
+    defaults = merge(default_p(sys), default_u0(sys))
 
     if u0map !== nothing
         u0map′ = lower_mapnames(u0map,get_iv(sys))
-        u0 = varmap_to_vars(u0map′,dvs; defaults=default_u0(sys))
+        u0 = varmap_to_vars(u0map′,dvs; defaults=defaults)
     else
         u0 = nothing
     end
 
-    defp = default_p(sys)
     if !(parammap isa DiffEqBase.NullParameters)
         parammap′ = lower_mapnames(parammap)
-        p = varmap_to_vars(parammap′,ps; defaults=defp)
-    elseif !isempty(defp)
-        p = varmap_to_vars(Dict(),ps; defaults=defp)
+        p = varmap_to_vars(parammap′,ps; defaults=defaults)
+    elseif !isempty(defaults)
+        p = varmap_to_vars(Dict(),ps; defaults=defaults)
     else
         p = ps
     end

src/systems/jumps/jumpsystem.jl

@@ -204,8 +204,9 @@ dprob = DiscreteProblem(js, u₀map, tspan, parammap)
 """
 function DiffEqBase.DiscreteProblem(sys::JumpSystem, u0map, tspan::Union{Tuple,Nothing},
                                     parammap=DiffEqBase.NullParameters(); kwargs...)
-    u0 = varmap_to_vars(u0map, states(sys); defaults=default_u0(sys))
-    p  = varmap_to_vars(parammap, parameters(sys); defaults=default_p(sys))
+    defaults = merge(default_p(sys), default_u0(sys))
+    u0 = varmap_to_vars(u0map, states(sys); defaults=defaults)
+    p  = varmap_to_vars(parammap, parameters(sys); defaults=defaults)
     f  = DiffEqBase.DISCRETE_INPLACE_DEFAULT
     df = DiscreteFunction{true,true}(f, syms=Symbol.(states(sys)))
     DiscreteProblem(df, u0, tspan, p; kwargs...)
@@ -232,8 +233,9 @@ dprob = DiscreteProblem(js, u₀map, tspan, parammap)
 """
 function DiscreteProblemExpr(sys::JumpSystem, u0map, tspan::Union{Tuple,Nothing},
                                     parammap=DiffEqBase.NullParameters(); kwargs...)
-    u0 = varmap_to_vars(u0map, states(sys); defaults=default_u0(sys))
-    p  = varmap_to_vars(parammap, parameters(sys); defaults=default_p(sys))
+    defaults = merge(default_p(sys), default_u0(sys))
+    u0 = varmap_to_vars(u0map, states(sys); defaults=defaults)
+    p  = varmap_to_vars(parammap, parameters(sys); defaults=defaults)
     # identity function to make syms works
     quote
         f  = DiffEqBase.DISCRETE_INPLACE_DEFAULT

src/systems/nonlinear/nonlinearsystem.jl

@@ -211,14 +211,14 @@ function process_NonlinearProblem(constructor, sys::NonlinearSystem,u0map,paramm
     dvs = states(sys)
     ps = parameters(sys)
     u0map′ = lower_mapnames(u0map)
-    u0 = varmap_to_vars(u0map′,dvs; defaults=default_u0(sys))
-    defp = default_p(sys)
+    defaults = merge(default_p(sys), default_u0(sys))
+    u0 = varmap_to_vars(u0map′,dvs; defaults=defaults)
 
     if !(parammap isa DiffEqBase.NullParameters)
         parammap′ = lower_mapnames(parammap)
-        p = varmap_to_vars(parammap′,ps; defaults=defp)
-    elseif !isempty(defp)
-        p = varmap_to_vars(Dict(),ps; defaults=defp)
+        p = varmap_to_vars(parammap′,ps; defaults=defaults)
+    elseif !isempty(default_p(sys))
+        p = varmap_to_vars(Dict(),ps; defaults=defaults)
     else
         p = ps
     end

src/systems/optimization/optimizationsystem.jl

@@ -157,8 +157,9 @@ function DiffEqBase.OptimizationProblem{iip}(sys::OptimizationSystem, u0,
 
     _f = DiffEqBase.OptimizationFunction{iip,AutoModelingToolkit,typeof(f),typeof(_grad),typeof(_hess),Nothing,Nothing,Nothing,Nothing}(f,AutoModelingToolkit(),_grad,_hess,nothing,nothing,nothing,nothing)
 
-    u0 = varmap_to_vars(u0,dvs; defaults=default_u0(sys))
-    p = varmap_to_vars(parammap,ps; defaults=default_p(sys))
+    defaults = merge(default_p(sys), default_u0(sys))
+    u0 = varmap_to_vars(u0,dvs; defaults=defaults)
+    p = varmap_to_vars(parammap,ps; defaults=defaults)
     lb = varmap_to_vars(lb,dvs)
     ub = varmap_to_vars(ub,dvs)
     OptimizationProblem{iip}(_f,u0,p;lb=lb,ub=ub,kwargs...)
@@ -212,8 +213,9 @@ function OptimizationProblemExpr{iip}(sys::OptimizationSystem, u0,
         _hess = :nothing
     end
 
-    u0 = varmap_to_vars(u0,dvs; defaults=default_u0(sys))
-    p = varmap_to_vars(parammap,ps; defaults=default_p(sys))
+    defaults = merge(default_p(sys), default_u0(sys))
+    u0 = varmap_to_vars(u0,dvs; defaults=defaults)
+    p = varmap_to_vars(parammap,ps; defaults=defaults)
     lb = varmap_to_vars(lb,dvs)
     ub = varmap_to_vars(ub,dvs)
     quote

test/symbolic_parameters.jl

@@ -16,7 +16,7 @@ par = [
 ]
 u0 = Pair{Num, Any}[
     x => u,
-    y => u,
+    y => σ, # default u0 from default p
     z => u-0.1,
 ]
 ns = NonlinearSystem(eqs, [x,y,z],[σ,ρ,β], name=:ns, default_p=par, default_u0=u0)