Problem and Function expr construction

`probexpr = ODEProblemExpr(sys,u0,tspan,p,jac=true)` gives:

```julia
quote
    f = begin
            f = ((var"##MTIIPVar#285", var"##MTKArg#281", var"##MTKArg#282", var"##MTKArg#283")->begin
                        @inbounds begin
                                let (xˍt, y, z, x, σ, β, ρ, t) = (var"##MTKArg#281"[1], var"##MTKArg#281"[2], var"##MTKArg#281"[3], var"##MTKArg#281"[4], var"##MTKArg#282"[1], var"##MTKArg#282"[2], var"##MTKArg#282"[3], var"##MTKArg#283")
                                    var"##MTIIPVar#285"[1] = σ * (y - x)
                                    var"##MTIIPVar#285"[2] = x * (ρ - z) - y
                                    var"##MTIIPVar#285"[3] = x * y - β * z
                                    var"##MTIIPVar#285"[4] = xˍt
                                end
                            end
                        nothing
                    end)
            tgrad = nothing
            jac = ((var"##MTIIPVar#291", var"##MTKArg#287", var"##MTKArg#288", var"##MTKArg#289")->begin
                        @inbounds begin
                                let (xˍt, y, z, x, σ, β, ρ, t) = (var"##MTKArg#287"[1], var"##MTKArg#287"[2], var"##MTKArg#287"[3], var"##MTKArg#287"[4], var"##MTKArg#288"[1], var"##MTKArg#288"[2], var"##MTKArg#288"[3], var"##MTKArg#289")
                                    var"##MTIIPVar#291"[1] = 0
                                    var"##MTIIPVar#291"[2] = 0
                                    var"##MTIIPVar#291"[3] = 0
                                    var"##MTIIPVar#291"[4] = 1
                                    var"##MTIIPVar#291"[5] = σ
                                    var"##MTIIPVar#291"[6] = -1
                                    var"##MTIIPVar#291"[7] = x
                                    var"##MTIIPVar#291"[8] = 0
                                    var"##MTIIPVar#291"[9] = 0
                                    var"##MTIIPVar#291"[10] = -1x
                                    var"##MTIIPVar#291"[11] = -1β
                                    var"##MTIIPVar#291"[12] = 0
                                    var"##MTIIPVar#291"[13] = -1σ
                                    var"##MTIIPVar#291"[14] = -1z + ρ
                                    var"##MTIIPVar#291"[15] = y
                                    var"##MTIIPVar#291"[16] = 0
                                end
                            end
                        nothing
                    end)
            Wfact = nothing
            Wfact_t = nothing
            M = UniformScaling{Bool}(true)
            ODEFunction{iip}(f, jac = jac, tgrad = tgrad, Wfact = Wfact, Wfact_t = Wfact_t, mass_matrix = M, syms = [:xˍt, :y, :z, :x])
        end
    u0 = [2.0, 0.0, 0.0, 1.0]
    tspan = (0.0, 100.0)
    p = [28.0, 2.6666666666666665, 10.0]
    ODEProblem{iip}(f, u0, tspan, p; )
end
```

Author: ChrisRackauckas

src/ModelingToolkit.jl

@@ -113,7 +113,7 @@ include("systems/dependency_graphs.jl")
 include("latexify_recipes.jl")
 include("build_function.jl")
 
-export ODESystem, ODEFunction
+export ODESystem, ODEFunction, ODEFunctionExpr, ODEProblemExpr
 export SDESystem, SDEFunction
 export JumpSystem
 export ODEProblem, SDEProblem, NonlinearProblem, OptimizationProblem, SteadyStateProblem

src/systems/diffeqs/abstractodesystem.jl

@@ -190,6 +190,79 @@ function DiffEqBase.ODEFunction{iip}(sys::AbstractODESystem, dvs = states(sys),
                       syms = Symbol.(states(sys)))
 end
 
+"""
+```julia
+function DiffEqBase.ODEFunctionExpr{iip}(sys::AbstractODESystem, dvs = states(sys),
+                                     ps = parameters(sys);
+                                     version = nothing, tgrad=false,
+                                     jac = false, Wfact = false,
+                                     sparse = false,
+                                     kwargs...) where {iip}
+```
+
+Create a Julia expression for an `ODEFunction` from the [`ODESystem`](@ref).
+The arguments `dvs` and `ps` are used to set the order of the dependent
+variable and parameter vectors, respectively.
+"""
+struct ODEFunctionExpr{iip} end
+
+function ODEFunctionExpr{iip}(sys::AbstractODESystem, dvs = states(sys),
+                                     ps = parameters(sys), u0 = nothing;
+                                     version = nothing, tgrad=false,
+                                     jac = false, Wfact = false,
+                                     sparse = false,
+                                     kwargs...) where {iip}
+
+    idx = iip ? 2 : 1
+    f = generate_function(sys, dvs, ps; expression=Val{true}, kwargs...)[idx]
+    if tgrad
+        _tgrad = generate_tgrad(sys, dvs, ps; expression=Val{true}, kwargs...)[idx]
+    else
+        _tgrad = :nothing
+    end
+
+    if jac
+        _jac = generate_jacobian(sys, dvs, ps; sparse = sparse, expression=Val{true}, kwargs...)[idx]
+    else
+        _jac = :nothing
+    end
+
+    if Wfact
+        tmp_Wfact,tmp_Wfact_t = generate_factorized_W(sys, dvs, ps; expression=Val{true}, kwargs...)
+        _Wfact =  tmp_Wfact[idx]
+        _Wfact_t =  tmp_Wfact_t[idx]
+    else
+        _Wfact,_Wfact_t = :nothing,:nothing
+    end
+
+    M = calculate_massmatrix(sys)
+
+    _M = (u0 === nothing || M == I) ? M : ArrayInterface.restructure(u0 .* u0',M)
+
+    quote
+        f = $f
+        tgrad = $_tgrad
+        jac = $_jac
+        Wfact = $_Wfact
+        Wfact_t = $_Wfact_t
+        M = $_M
+
+        ODEFunction{iip}(f,
+                         jac = jac,
+                         tgrad = tgrad,
+                         Wfact = Wfact,
+                         Wfact_t = Wfact_t,
+                         mass_matrix = M,
+                         syms = $(Symbol.(states(sys))))
+    end
+end
+
+
+function ODEFunctionExpr(sys::AbstractODESystem, args...; kwargs...)
+    ODEFunctionExpr{true}(sys, args...; kwargs...)
+end
+
+
 function DiffEqBase.ODEProblem(sys::AbstractODESystem, args...; kwargs...)
     ODEProblem{true}(sys, args...; kwargs...)
 end
@@ -225,6 +298,50 @@ function DiffEqBase.ODEProblem{iip}(sys::AbstractODESystem,u0map,tspan,
     ODEProblem{iip}(f,u0,tspan,p;kwargs...)
 end
 
+"""
+```julia
+function DiffEqBase.ODEProblemExpr{iip}(sys::AbstractODESystem,u0map,tspan,
+                                    parammap=DiffEqBase.NullParameters();
+                                    version = nothing, tgrad=false,
+                                    jac = false, Wfact = false,
+                                    checkbounds = false, sparse = false,
+                                    linenumbers = true, parallel=SerialForm(),
+                                    kwargs...) where iip
+```
+
+Generates a Julia expression for constructing an ODEProblem from an
+ODESystem and allows for automatically symbolically calculating
+numerical enhancements.
+"""
+struct ODEProblemExpr{iip} end
+
+function ODEProblemExpr{iip}(sys::AbstractODESystem,u0map,tspan,
+                                    parammap=DiffEqBase.NullParameters();
+                                    version = nothing, tgrad=false,
+                                    jac = false, Wfact = 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{iip}(sys,dvs,ps,u0;tgrad=tgrad,jac=jac,Wfact=Wfact,checkbounds=checkbounds,
+                        linenumbers=linenumbers,parallel=parallel,
+                        sparse=sparse)
+    quote
+        f = $f
+        u0 = $u0
+        tspan = $tspan
+        p = $p
+        ODEProblem{iip}(f,u0,tspan,p;$(kwargs...))
+    end
+end
+
+function ODEProblemExpr(sys::AbstractODESystem, args...; kwargs...)
+    ODEProblemExpr{true}(sys, args...; kwargs...)
+end
+
 
 ### Enables Steady State Problems ###
 function DiffEqBase.SteadyStateProblem(sys::AbstractODESystem, args...; kwargs...)

test/lowering_solving.jl

@@ -31,7 +31,10 @@ p  = [σ => 28.0,
 
 tspan = (0.0,100.0)
 prob = ODEProblem(sys,u0,tspan,p,jac=true)
+probexpr = ODEProblemExpr(sys,u0,tspan,p,jac=true)
 sol = solve(prob,Tsit5())
+solexpr = solve(eval(prob),Tsit5())
+@test all(x->x==0,Array(sol - solexpr))
 #using Plots; plot(sol,vars=(:x,:y))
 
 @parameters t σ ρ β