Merge pull request #1548 from SciML/myb/promote

Promote to concrete type by default if it's an array

Author: ChrisRackauckas

Project.toml

@@ -82,6 +82,7 @@ julia = "1.6"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 GalacticOptim = "a75be94c-b780-496d-a8a9-0878b188d577"
+GalacticOptimJL = "9d3c5eb1-403b-401b-8c0f-c11105342e6b"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -93,4 +94,4 @@ Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["BenchmarkTools", "ForwardDiff", "GalacticOptim", "OrdinaryDiffEq", "Optim", "Random", "ReferenceTests", "SafeTestsets", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials"]
+test = ["BenchmarkTools", "ForwardDiff", "GalacticOptim", "GalacticOptimJL", "OrdinaryDiffEq", "Optim", "Random", "ReferenceTests", "SafeTestsets", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials"]

src/structural_transformation/codegen.jl

@@ -521,8 +521,9 @@ function ODAEProblem{iip}(
                           sys,
                           u0map,
                           tspan,
-                          parammap=DiffEqBase.NullParameters();
+                          parammap = DiffEqBase.NullParameters();
                           callback = nothing,
+                          use_union = false,
                           kwargs...
                          ) where {iip}
     fun, dvs = build_torn_function(sys; kwargs...)
@@ -531,8 +532,8 @@ function ODAEProblem{iip}(
 
     defs = ModelingToolkit.mergedefaults(defs,parammap,ps)
     defs = ModelingToolkit.mergedefaults(defs,u0map,dvs)
-    u0 = ModelingToolkit.varmap_to_vars(u0map, dvs; defaults=defs)
-    p = ModelingToolkit.varmap_to_vars(parammap, ps; defaults=defs)
+    u0 = ModelingToolkit.varmap_to_vars(u0map, dvs; defaults=defs, tofloat=true)
+    p = ModelingToolkit.varmap_to_vars(parammap, ps; defaults=defs, tofloat=!use_union, use_union)
 
     has_difference = any(isdifferenceeq, equations(sys))
     if has_continuous_events(sys)

src/systems/diffeqs/abstractodesystem.jl

@@ -599,6 +599,7 @@ function process_DEProblem(constructor, sys::AbstractODESystem,u0map,parammap;
                            simplify=false,
                            linenumbers = true, parallel=SerialForm(),
                            eval_expression = true,
+                           use_union = false,
                            kwargs...)
     eqs = equations(sys)
     dvs = states(sys)
@@ -609,12 +610,12 @@ function process_DEProblem(constructor, sys::AbstractODESystem,u0map,parammap;
     defs = mergedefaults(defs,parammap,ps)
     defs = mergedefaults(defs,u0map,dvs)
 
-    u0 = varmap_to_vars(u0map,dvs; defaults=defs, promotetoconcrete=true)
-    p = varmap_to_vars(parammap,ps; defaults=defs)
+    u0 = varmap_to_vars(u0map, dvs; defaults=defs, tofloat=true)
+    p = varmap_to_vars(parammap, ps; defaults=defs, tofloat=!use_union, use_union)
     if implicit_dae && du0map !== nothing
         ddvs = map(Differential(iv), dvs)
         defs = mergedefaults(defs,du0map, ddvs)
-        du0 = varmap_to_vars(du0map,ddvs; defaults=defs, toterm=identity, promotetoconcrete=true)
+        du0 = varmap_to_vars(du0map,ddvs; defaults=defs, toterm=identity, tofloat=true)
     else
         du0 = nothing
         ddvs = nothing

src/systems/discrete_system/discrete_system.jl

@@ -169,19 +169,20 @@ function DiffEqBase.DiscreteProblem(sys::DiscreteSystem,u0map,tspan,
                                     parammap=DiffEqBase.NullParameters();
                                     eval_module = @__MODULE__,
                                     eval_expression = true,
+                                    use_union = false,
                                     kwargs...)
     dvs = states(sys)
     ps = parameters(sys)
     eqs = equations(sys)
     eqs = linearize_eqs(sys, eqs)
     iv = get_iv(sys)
-    
+
     defs = defaults(sys)
     defs = mergedefaults(defs,parammap,ps)
     defs = mergedefaults(defs,u0map,dvs)
-    
-    u0 = varmap_to_vars(u0map,dvs; defaults=defs, promotetoconcrete=true)
-    p = varmap_to_vars(parammap,ps; defaults=defs)
+
+    u0 = varmap_to_vars(u0map, dvs; defaults=defs, tofloat=false)
+    p = varmap_to_vars(parammap, ps; defaults=defs, tofloat=false, use_union)
 
     rhss = [eq.rhs for eq in eqs]
     u = dvs

src/systems/jumps/jumpsystem.jl

@@ -202,7 +202,9 @@ end
 """
 ```julia
 function DiffEqBase.DiscreteProblem(sys::JumpSystem, u0map, tspan,
-                                    parammap=DiffEqBase.NullParameters; kwargs...)
+                                    parammap=DiffEqBase.NullParameters;
+                                    use_union=false,
+                                    kwargs...)
 ```
 
 Generates a blank DiscreteProblem for a pure jump JumpSystem to utilize as
@@ -219,20 +221,22 @@ dprob = DiscreteProblem(js, u₀map, tspan, parammap)
 ```
 """
 function DiffEqBase.DiscreteProblem(sys::JumpSystem, u0map, tspan::Union{Tuple,Nothing},
-                                    parammap=DiffEqBase.NullParameters(); checkbounds=false, kwargs...)
-    
+                                    parammap=DiffEqBase.NullParameters(); checkbounds=false,
+                                    use_union=false,
+                                    kwargs...)
+
     dvs = states(sys)
     ps = parameters(sys)
-    
+
     defs = defaults(sys)
     defs = mergedefaults(defs,parammap,ps)
-    defs = mergedefaults(defs,u0map,dvs) 
-    
-    u0 = varmap_to_vars(u0map,dvs; defaults=defs, promotetoconcrete=true)
-    p = varmap_to_vars(parammap,ps; defaults=defs)
-        
+    defs = mergedefaults(defs,u0map,dvs)
+
+    u0 = varmap_to_vars(u0map, dvs; defaults=defs, tofloat=false)
+    p = varmap_to_vars(parammap, ps; defaults=defs, tofloat=false, use_union)
+
     f  = DiffEqBase.DISCRETE_INPLACE_DEFAULT
-    
+
     # just taken from abstractodesystem.jl for ODEFunction def
     obs = observed(sys)
     observedfun = let sys = sys, dict = Dict()
@@ -268,10 +272,15 @@ dprob = DiscreteProblem(js, u₀map, tspan, parammap)
 ```
 """
 function DiscreteProblemExpr(sys::JumpSystem, u0map, tspan::Union{Tuple,Nothing},
-                                    parammap=DiffEqBase.NullParameters(); kwargs...)
+                                    parammap=DiffEqBase.NullParameters();
+                                    use_union=false,
+                                    kwargs...)
+    dvs = states(sys)
+    ps = parameters(sys)
     defs = defaults(sys)
-    u0 = varmap_to_vars(u0map, states(sys); defaults=defs)
-    p  = varmap_to_vars(parammap, parameters(sys); defaults=defs)
+
+    u0 = varmap_to_vars(u0map, dvs; defaults=defs, tofloat=false)
+    p = varmap_to_vars(parammap, ps; defaults=defs, tofloat=false, use_union)
     # identity function to make syms works
     quote
         f  = DiffEqBase.DISCRETE_INPLACE_DEFAULT

src/systems/nonlinear/nonlinearsystem.jl

@@ -260,17 +260,18 @@ function process_NonlinearProblem(constructor, sys::NonlinearSystem,u0map,paramm
                            simplify=false,
                            linenumbers = true, parallel=SerialForm(),
                            eval_expression = true,
+                           use_union = false,
                            kwargs...)
     eqs = equations(sys)
     dvs = states(sys)
     ps = parameters(sys)
-    
+
     defs = defaults(sys)
     defs = mergedefaults(defs,parammap,ps)
     defs = mergedefaults(defs,u0map,dvs)
-    
-    u0 = varmap_to_vars(u0map,dvs; defaults=defs, promotetoconcrete=true)
-    p = varmap_to_vars(parammap,ps; defaults=defs)
+
+    u0 = varmap_to_vars(u0map, dvs; defaults=defs, tofloat=true)
+    p = varmap_to_vars(parammap, ps; defaults=defs, tofloat=!use_union, use_union)
 
     check_eqs_u0(eqs, dvs, u0; kwargs...)
 

src/systems/optimization/optimizationsystem.jl

@@ -149,6 +149,7 @@ function DiffEqBase.OptimizationProblem{iip}(sys::OptimizationSystem, u0map,
                                           hess = false, sparse = false,
                                           checkbounds = false,
                                           linenumbers = true, parallel=SerialForm(),
+                                          use_union = false,
                                           kwargs...) where iip
     dvs = states(sys)
     ps = parameters(sys)
@@ -180,10 +181,10 @@ function DiffEqBase.OptimizationProblem{iip}(sys::OptimizationSystem, u0map,
     defs = mergedefaults(defs,parammap,ps)
     defs = mergedefaults(defs,u0map,dvs)
 
-    u0 = varmap_to_vars(u0map,dvs; defaults=defs)
-    p = varmap_to_vars(parammap,ps; defaults=defs)
-    lb = varmap_to_vars(lb,dvs; check=false)
-    ub = varmap_to_vars(ub,dvs; check=false)
+    u0 = varmap_to_vars(u0map, dvs; defaults=defs, tofloat=false)
+    p = varmap_to_vars(parammap, ps; defaults=defs, tofloat=false, use_union)
+    lb = varmap_to_vars(lb, dvs; check=false, tofloat=false, use_union)
+    ub = varmap_to_vars(ub, dvs; check=false, tofloat=false, use_union)
     OptimizationProblem{iip}(_f,u0,p;lb=lb,ub=ub,kwargs...)
 end
 
@@ -215,6 +216,7 @@ function OptimizationProblemExpr{iip}(sys::OptimizationSystem, u0,
                                           hess = false, sparse = false,
                                           checkbounds = false,
                                           linenumbers = false, parallel=SerialForm(),
+                                          use_union = false,
                                           kwargs...) where iip
     dvs = states(sys)
     ps = parameters(sys)
@@ -239,10 +241,10 @@ function OptimizationProblemExpr{iip}(sys::OptimizationSystem, u0,
     defs = mergedefaults(defs,parammap,ps)
     defs = mergedefaults(defs,u0map,dvs)
 
-    u0 = varmap_to_vars(u0map,dvs; defaults=defs, promotetoconcrete=true)
-    p = varmap_to_vars(parammap,ps; defaults=defs)
-    lb = varmap_to_vars(lb,dvs)
-    ub = varmap_to_vars(ub,dvs)
+    u0 = varmap_to_vars(u0map, dvs; defaults=defs, tofloat=false)
+    p = varmap_to_vars(parammap, ps; defaults=defs, tofloat=false, use_union)
+    lb = varmap_to_vars(lb, dvs; check=false, tofloat=false, use_union)
+    ub = varmap_to_vars(ub, dvs; check=false, tofloat=false, use_union)
     quote
         f = $f
         p = $p

src/utils.jl

@@ -472,15 +472,31 @@ function mergedefaults(defaults, varmap, vars)
     end
 end
 
-function promote_to_concrete(vs)
-    if isempty(vs) 
+function promote_to_concrete(vs; tofloat=true, use_union=false)
+    if isempty(vs)
         return vs
     end
     T = eltype(vs)
-    if Base.isconcretetype(T) # nothing to do
+    if Base.isconcretetype(T) && (!tofloat || T === float(T)) # nothing to do
         vs
     else
-        C = foldl((t, elem)->promote_type(t, eltype(elem)), vs; init=typeof(first(vs)))
+        C = typeof(first(vs))
+        has_int = false
+        I = Int8
+        for v in vs
+            E = eltype(v)
+            C = promote_type(C, E)
+            if E <: Integer
+                has_int = true
+                I = promote_type(I, E)
+            end
+        end
+        if tofloat
+            C = float(C)
+        elseif use_union && has_int && C !== I
+            C = Union{C, I}
+            return copyto!(similar(vs, C), vs)
+        end
         convert.(C, vs)
     end
-end
+end

src/variables.jl

@@ -32,7 +32,7 @@ Takes a list of pairs of `variables=>values` and an ordered list of variables
 and creates the array of values in the correct order with default values when
 applicable.
 """
-function varmap_to_vars(varmap, varlist; defaults=Dict(), check=true, toterm=Symbolics.diff2term, promotetoconcrete=false)
+function varmap_to_vars(varmap, varlist; defaults=Dict(), check=true, toterm=Symbolics.diff2term, promotetoconcrete=nothing, tofloat=true, use_union=false)
     varlist = map(unwrap, varlist)
     # Edge cases where one of the arguments is effectively empty.
     is_incomplete_initialization = varmap isa DiffEqBase.NullParameters || varmap === nothing
@@ -58,16 +58,16 @@ function varmap_to_vars(varmap, varlist; defaults=Dict(), check=true, toterm=Sym
         varmap
     end
 
+    promotetoconcrete === nothing && (promotetoconcrete = container_type <: AbstractArray)
     if promotetoconcrete
-        vals = promote_to_concrete(vals)
+        vals = promote_to_concrete(vals; tofloat=tofloat, use_union=use_union)
     end
 
     if isempty(vals)
         return nothing
     elseif container_type <: Tuple
         (vals...,)
     else
-        vals = identity.(vals)
         SymbolicUtils.Code.create_array(container_type, eltype(vals), Val{1}(), Val(length(vals)), vals...)
     end
 end
@@ -79,7 +79,7 @@ function _varmap_to_vars(varmap::Dict, varlist; defaults=Dict(), check=false, to
     for (p, v) in pairs(varmap)
         varmap[p] = fixpoint_sub(v, varmap)
     end
-    
+
     missingvars = setdiff(varlist, keys(varmap))
     check && (isempty(missingvars) || throw_missingvars(missingvars))
 

test/controlsystem.jl

@@ -1,4 +1,4 @@
-using ModelingToolkit, GalacticOptim, Optim
+using ModelingToolkit, GalacticOptim, Optim, GalacticOptimJL
 
 @variables t x(t) v(t) u(t)
 @parameters p[1:2]