Merge remote-tracking branch 'origin/master' into sm/symutils3

Author: shashi

Project.toml

@@ -1,7 +1,7 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Chris Rackauckas <[email protected]>"]
-version = "3.20.1"
+version = "3.21.0"
 
 [deps]
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
@@ -11,7 +11,6 @@ DiffEqJump = "c894b116-72e5-5b58-be3c-e6d8d4ac2b12"
 DiffRules = "b552c78f-8df3-52c6-915a-8e097449b14b"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
-GeneralizedGenerated = "6b9d7cbe-bcb9-11e9-073f-15a7a543e2eb"
 IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
 LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800"
 Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
@@ -22,6 +21,7 @@ MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
+RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
@@ -34,11 +34,10 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 [compat]
 ArrayInterface = "2.8"
 DataStructures = "0.17, 0.18"
-DiffEqBase = "6.38"
+DiffEqBase = "6.48"
 DiffEqJump = "6.7.5"
 DiffRules = "0.1, 1.0"
 DocStringExtensions = "0.7, 0.8"
-GeneralizedGenerated = "0.1.4, 0.2"
 IfElse = "0.1"
 LabelledArrays = "1.3"
 Latexify = "0.11, 0.12, 0.13, 0.14"
@@ -47,6 +46,7 @@ MacroTools = "0.5"
 NaNMath = "0.3"
 RecursiveArrayTools = "2.3"
 Requires = "1.0"
+RuntimeGeneratedFunctions = "0.4"
 SafeTestsets = "0.0.1"
 SpecialFunctions = "0.7, 0.8, 0.9, 0.10"
 StaticArrays = "0.10, 0.11, 0.12"
@@ -59,11 +59,13 @@ julia = "1.2"
 [extras]
 Dagger = "d58978e5-989f-55fb-8d15-ea34adc7bf54"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+GalacticOptim = "a75be94c-b780-496d-a8a9-0878b188d577"
+Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
 StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Dagger", "ForwardDiff", "OrdinaryDiffEq", "Random", "SteadyStateDiffEq", "Test", "StochasticDiffEq"]
+test = ["Dagger", "ForwardDiff", "GalacticOptim", "OrdinaryDiffEq", "Optim", "Random", "SteadyStateDiffEq", "Test", "StochasticDiffEq"]

src/ModelingToolkit.jl

@@ -8,14 +8,16 @@ using UnPack: @unpack
 using DiffEqJump
 using DataStructures: OrderedDict, OrderedSet
 using SpecialFunctions, NaNMath
+using RuntimeGeneratedFunctions
 using Base.Threads
 import MacroTools: splitdef, combinedef, postwalk, striplines
-import GeneralizedGenerated
 import Libdl
 using DocStringExtensions
 using Base: RefValue
 import IfElse
 
+RuntimeGeneratedFunctions.init(@__MODULE__)
+
 using RecursiveArrayTools
 
 import SymbolicUtils

src/build_function.jl

@@ -130,19 +130,7 @@ function _build_function(target::JuliaTarget, op, args...;
 end
 
 function _build_and_inject_function(mod::Module, ex)
-    # Generate the function, which will process the expression
-    runtimefn = GeneralizedGenerated.mk_function(mod, ex)
-
-    # Extract the processed expression of the function body
-    params = typeof(runtimefn).parameters
-    fn_expr = GeneralizedGenerated.NGG.from_type(params[3])
-
-    # Inject our externally registered module functions
-    new_expr = ModelingToolkit.inject_registered_module_functions(fn_expr)
-
-    # Reconstruct the RuntimeFn's Body
-    new_body = GeneralizedGenerated.NGG.to_type(new_expr)
-    return GeneralizedGenerated.RuntimeFn{params[1:2]..., new_body, params[4]}()
+	@RuntimeGeneratedFunction(ModelingToolkit.inject_registered_module_functions(ex))
 end
 
 # Detect heterogeneous element types of "arrays of matrices/sparce matrices"
@@ -186,7 +174,7 @@ function _build_function(target::JuliaTarget, rhss, args...;
 
 Generates a Julia function which can then be utilized for further evaluations.
 If expression=Val{false}, the return is a Julia function which utilizes
-GeneralizedGenerated.jl in order to be free of world-age issues.
+RuntimeGeneratedFunctions.jl in order to be free of world-age issues.
 
 If the `Expression` is an `Operation`, the generated function is a function
 with a scalar output, otherwise if it's an `AbstractArray{Operation}`, the output
@@ -570,7 +558,7 @@ function _build_function(target::CTarget, eqs::Array{<:Equation}, args...;
 		open(`gcc -fPIC -O3 -msse3 -xc -shared -o $(libpath * "." * Libdl.dlext) -`, "w") do f
 			print(f, ex)
 		end
-		eval(:((du::Array{Float64},u::Array{Float64},p::Array{Float64},t::Float64) -> ccall(("diffeqf", $libpath), Cvoid, (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Float64), du, u, p, t)))
+		@RuntimeGeneratedFunction(:((du::Array{Float64},u::Array{Float64},p::Array{Float64},t::Float64) -> ccall(("diffeqf", $libpath), Cvoid, (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Float64), du, u, p, t)))
 	end
 end
 

src/systems/control/controlsystem.jl

@@ -132,8 +132,8 @@ function runge_kutta_discretize(sys::ControlSystem,dt,tspan;
     n = length(tspan[1]:dt:tspan[2]) - 1
     m = length(tab.α)
 
-    f = eval(build_function([x.rhs for x in equations(sys)],sys.states,sys.controls,sys.ps,sys.iv,conv = ModelingToolkit.ControlToExpr(sys))[1])
-    L = eval(build_function(sys.loss,sys.states,sys.controls,sys.ps,sys.iv,conv = ModelingToolkit.ControlToExpr(sys)))
+    f = @RuntimeGeneratedFunction(build_function([x.rhs for x in equations(sys)],sys.states,sys.controls,sys.ps,sys.iv,conv = ModelingToolkit.ControlToExpr(sys))[1])
+    L = @RuntimeGeneratedFunction(build_function(sys.loss,sys.states,sys.controls,sys.ps,sys.iv,conv = ModelingToolkit.ControlToExpr(sys)))
 
     var(n, i...) = var(nameof(n), i...)
     var(n::Symbol, i...) = Sym{FnType{Tuple{symtype(sys.iv)}, Number}}(nameof(Variable(n, i...)))

src/systems/diffeqs/abstractodesystem.jl

@@ -133,15 +133,15 @@ function DiffEqBase.ODEFunction{iip}(sys::AbstractODESystem, dvs = states(sys),
                                      kwargs...) where {iip}
 
     f_gen = generate_function(sys, dvs, ps; expression=Val{eval_expression}, kwargs...)
-    f_oop,f_iip = eval_expression ? ModelingToolkit.eval.(f_gen) : f_gen
+    f_oop,f_iip = eval_expression ? (@RuntimeGeneratedFunction(ex) for ex in f_gen) : f_gen
     f(u,p,t) = f_oop(u,p,t)
     f(du,u,p,t) = f_iip(du,u,p,t)
 
     if tgrad
         tgrad_gen = generate_tgrad(sys, dvs, ps;
                                    simplify=simplify,
                                    expression=Val{eval_expression}, kwargs...)
-        tgrad_oop,tgrad_iip = eval_expression ? ModelingToolkit.eval.(tgrad_gen) : tgrad_gen
+        tgrad_oop,tgrad_iip = eval_expression ? (@RuntimeGeneratedFunction(ex) for ex in tgrad_gen) : tgrad_gen
         _tgrad(u,p,t) = tgrad_oop(u,p,t)
         _tgrad(J,u,p,t) = tgrad_iip(J,u,p,t)
     else
@@ -152,7 +152,7 @@ function DiffEqBase.ODEFunction{iip}(sys::AbstractODESystem, dvs = states(sys),
         jac_gen = generate_jacobian(sys, dvs, ps;
                                     simplify=simplify, sparse = sparse,
                                     expression=Val{eval_expression}, kwargs...)
-        jac_oop,jac_iip = eval_expression ? ModelingToolkit.eval.(jac_gen) : jac_gen
+        jac_oop,jac_iip = eval_expression ? (@RuntimeGeneratedFunction(ex) for ex in jac_gen) : jac_gen
         _jac(u,p,t) = jac_oop(u,p,t)
         _jac(J,u,p,t) = jac_iip(J,u,p,t)
     else
@@ -163,12 +163,12 @@ function DiffEqBase.ODEFunction{iip}(sys::AbstractODESystem, dvs = states(sys),
 
     _M = (u0 === nothing || M == I) ? M : ArrayInterface.restructure(u0 .* u0',M)
 
-    ODEFunction{iip}(DiffEqBase.EvalFunc(f),
-                      jac = _jac === nothing ? nothing : DiffEqBase.EvalFunc(_jac),
-                      tgrad = _tgrad === nothing ? nothing : DiffEqBase.EvalFunc(_tgrad),
-                      mass_matrix = _M,
-                      jac_prototype = sparse ? similar(sys.jac[],Float64) : nothing,
-                      syms = Symbol.(states(sys)))
+    ODEFunction{iip}(f,
+                     jac = _jac === nothing ? nothing : _jac,
+                     tgrad = _tgrad === nothing ? nothing : _tgrad,
+                     mass_matrix = _M,
+                     jac_prototype = sparse ? similar(sys.jac[],Float64) : nothing,
+                     syms = Symbol.(states(sys)))
 end
 
 """

src/systems/diffeqs/modelingtoolkitize.jl

@@ -87,3 +87,27 @@ function modelingtoolkitize(prob::DiffEqBase.SDEProblem)
 
     de
 end
+
+
+"""
+$(TYPEDSIGNATURES)
+
+Generate `OptimizationSystem`, dependent variables, and parameters from an `OptimizationProblem`.
+"""
+function modelingtoolkitize(prob::DiffEqBase.OptimizationProblem)
+
+    if prob.p isa Tuple || prob.p isa NamedTuple
+        p = [x for x in prob.p]
+    else
+        p = prob.p
+    end
+
+    vars = reshape([Variable(:x, i)() for i in eachindex(prob.u0)],size(prob.u0))
+    params = p isa DiffEqBase.NullParameters ? [] :
+             reshape([Variable(:α,i)() for i in eachindex(p)],size(Array(p)))
+
+
+    eqs = prob.f(vars, params)
+    de = OptimizationSystem(eqs,vec(vars),vec(params))
+    de
+end

src/systems/diffeqs/sdesystem.jl

@@ -148,9 +148,9 @@ function DiffEqBase.SDEFunction{iip}(sys::SDESystem, dvs = sys.states, ps = sys.
                                      version = nothing, tgrad=false, sparse = false,
                                      jac = false, Wfact = false, eval_expression = true, kwargs...) where {iip}
     f_gen = generate_function(sys, dvs, ps; expression=Val{eval_expression}, kwargs...)
-    f_oop,f_iip = eval_expression ? ModelingToolkit.eval.(f_gen) : f_gen
+    f_oop,f_iip = eval_expression ? (@RuntimeGeneratedFunction(ex) for ex in f_gen) : f_gen
     g_gen = generate_diffusion_function(sys, dvs, ps; expression=Val{eval_expression}, kwargs...)
-    g_oop,g_iip = eval_expression ? ModelingToolkit.eval.(g_gen) : g_gen
+    g_oop,g_iip = eval_expression ? (@RuntimeGeneratedFunction(ex) for ex in g_gen) : g_gen
 
     f(u,p,t) = f_oop(u,p,t)
     f(du,u,p,t) = f_iip(du,u,p,t)
@@ -159,7 +159,7 @@ function DiffEqBase.SDEFunction{iip}(sys::SDESystem, dvs = sys.states, ps = sys.
 
     if tgrad
         tgrad_gen = generate_tgrad(sys, dvs, ps; expression=Val{eval_expression}, kwargs...)
-        tgrad_oop,tgrad_iip = eval_expression ? ModelingToolkit.eval.(tgrad_gen) : tgrad_gen
+        tgrad_oop,tgrad_iip = eval_expression ? (@RuntimeGeneratedFunction(ex) for ex in tgrad_gen) : tgrad_gen
         _tgrad(u,p,t) = tgrad_oop(u,p,t)
         _tgrad(J,u,p,t) = tgrad_iip(J,u,p,t)
     else
@@ -168,7 +168,7 @@ function DiffEqBase.SDEFunction{iip}(sys::SDESystem, dvs = sys.states, ps = sys.
 
     if jac
         jac_gen = generate_jacobian(sys, dvs, ps; expression=Val{eval_expression}, sparse=sparse, kwargs...)
-        jac_oop,jac_iip = eval_expression ? ModelingToolkit.eval.(jac_gen) : jac_gen
+        jac_oop,jac_iip = eval_expression ? (@RuntimeGeneratedFunction(ex) for ex in jac_gen) : jac_gen
         _jac(u,p,t) = jac_oop(u,p,t)
         _jac(J,u,p,t) = jac_iip(J,u,p,t)
     else
@@ -177,8 +177,8 @@ function DiffEqBase.SDEFunction{iip}(sys::SDESystem, dvs = sys.states, ps = sys.
 
     if Wfact
         tmp_Wfact,tmp_Wfact_t = generate_factorized_W(sys, dvs, ps, true; expression=Val{true}, kwargs...)
-        Wfact_oop, Wfact_iip = eval_expression ? ModelingToolkit.eval.(tmp_Wfact) : tmp_Wfact
-        Wfact_oop_t, Wfact_iip_t = eval_expression ? ModelingToolkit.eval.(tmp_Wfact_t) : tmp_Wfact_t
+        Wfact_oop, Wfact_iip = eval_expression ? (@RuntimeGeneratedFunction(ex) for ex in tmp_Wfact) : tmp_Wfact
+        Wfact_oop_t, Wfact_iip_t = eval_expression ? (@RuntimeGeneratedFunction(ex) for ex in tmp_Wfact_t) : tmp_Wfact_t
         _Wfact(u,p,dtgamma,t) = Wfact_oop(u,p,dtgamma,t)
         _Wfact(W,u,p,dtgamma,t) = Wfact_iip(W,u,p,dtgamma,t)
         _Wfact_t(u,p,dtgamma,t) = Wfact_oop_t(u,p,dtgamma,t)
@@ -190,13 +190,13 @@ function DiffEqBase.SDEFunction{iip}(sys::SDESystem, dvs = sys.states, ps = sys.
     M = calculate_massmatrix(sys)
     _M = (u0 === nothing || M == I) ? M : ArrayInterface.restructure(u0 .* u0',M)
 
-    SDEFunction{iip}(DiffEqBase.EvalFunc(f),DiffEqBase.EvalFunc(g),
-                      jac = _jac === nothing ? nothing : DiffEqBase.EvalFunc(_jac),
-                      tgrad = _tgrad === nothing ? nothing : DiffEqBase.EvalFunc(_tgrad),
-                      Wfact = _Wfact === nothing ? nothing : DiffEqBase.EvalFunc(_Wfact),
-                      Wfact_t = _Wfact_t === nothing ? nothing : DiffEqBase.EvalFunc(_Wfact_t),
-                      mass_matrix = _M,
-                      syms = Symbol.(sys.states))
+    SDEFunction{iip}(f,g,
+                     jac = _jac === nothing ? nothing : _jac,
+                     tgrad = _tgrad === nothing ? nothing : _tgrad,
+                     Wfact = _Wfact === nothing ? nothing : _Wfact,
+                     Wfact_t = _Wfact_t === nothing ? nothing : _Wfact_t,
+                     mass_matrix = _M,
+                     syms = Symbol.(sys.states))
 end
 
 function DiffEqBase.SDEFunction(sys::SDESystem, args...; kwargs...)

src/systems/jumps/jumpsystem.jl

@@ -85,10 +85,10 @@ function generate_affect_function(js, affect, outputidxs)
 end
 
 function assemble_vrj(js, vrj, statetoid)
-    rate   = eval(generate_rate_function(js, vrj.rate))
+    rate   = @RuntimeGeneratedFunction(generate_rate_function(js, vrj.rate))
     outputvars = (value(affect.lhs) for affect in vrj.affect!)
     outputidxs = ((statetoid[var] for var in outputvars)...,)
-    affect = eval(generate_affect_function(js, vrj.affect!, outputidxs))
+    affect = @RuntimeGeneratedFunction(generate_affect_function(js, vrj.affect!, outputidxs))
     VariableRateJump(rate, affect)
 end
 
@@ -105,10 +105,10 @@ function assemble_vrj_expr(js, vrj, statetoid)
 end
 
 function assemble_crj(js, crj, statetoid)
-    rate   = eval(generate_rate_function(js, crj.rate))
+    rate   = @RuntimeGeneratedFunction(generate_rate_function(js, crj.rate))
     outputvars = (value(affect.lhs) for affect in crj.affect!)
     outputidxs = ((statetoid[var] for var in outputvars)...,)
-    affect = eval(generate_affect_function(js, crj.affect!, outputidxs))
+    affect = @RuntimeGeneratedFunction(generate_affect_function(js, crj.affect!, outputidxs))
     ConstantRateJump(rate, affect)
 end
 
@@ -196,8 +196,7 @@ function DiffEqBase.DiscreteProblem(sys::JumpSystem, u0map, tspan::Tuple,
     else
         p = parammap
     end
-    # EvalFunc because we know that the jump functions are generated via eval
-    f  = DiffEqBase.EvalFunc(DiffEqBase.DISCRETE_INPLACE_DEFAULT)
+    f  = DiffEqBase.DISCRETE_INPLACE_DEFAULT
     df = DiscreteFunction(f, syms=Symbol.(states(sys)))
     DiscreteProblem(df, u0, tspan, p; kwargs...)
 end
@@ -226,9 +225,8 @@ function DiscreteProblemExpr(sys::JumpSystem, u0map, tspan::Tuple,
     u0 = varmap_to_vars(u0map, states(sys))
     p  = varmap_to_vars(parammap, parameters(sys))
     # identity function to make syms works
-    # EvalFunc because we know that the jump functions are generated via eval
     quote
-        f  = DiffEqBase.EvalFunc(DiffEqBase.DISCRETE_INPLACE_DEFAULT)
+        f  = DiffEqBase.DISCRETE_INPLACE_DEFAULT
         u0 = $u0
         p = $p
         tspan = $tspan
@@ -282,7 +280,7 @@ end
 
 
 ### Functions to determine which states a jump depends on
-function get_variables!(dep, jump::Union{ConstantRateJump,VariableRateJump}, variables) 
+function get_variables!(dep, jump::Union{ConstantRateJump,VariableRateJump}, variables)
     (jump.rate isa Symbolic) && get_variables!(dep, jump.rate, variables)
     dep
 end