add a lot of expr construction pieces

Author: ChrisRackauckas

src/systems/jumps/jumpsystem.jl

@@ -63,33 +63,55 @@ function JumpSystem(eqs, iv, states, ps; systems = JumpSystem[],
 end
 
 
-generate_rate_function(js, rate) = ModelingToolkit.eval(
-                                        build_function(rate, states(js), parameters(js),
+generate_rate_function(js, rate) = build_function(rate, states(js), parameters(js),
                                         independent_variable(js),
-                                        expression=Val{true}))
+                                        expression=Val{true})
 
-generate_affect_function(js, affect, outputidxs) = ModelingToolkit.eval(
-                                                      build_function(affect, states(js),
+generate_affect_function(js, affect, outputidxs) = build_function(affect, states(js),
                                                       parameters(js),
                                                       independent_variable(js),
                                                       expression=Val{true},
                                                       headerfun=add_integrator_header,
-                                                      outputidxs=outputidxs)[2])
+                                                      outputidxs=outputidxs)[2]
 
 function assemble_vrj(js, vrj, statetoid)
+    rate   = eval(generate_rate_function(js, vrj.rate))
+    outputvars = (convert(Variable,affect.lhs) for affect in vrj.affect!)
+    outputidxs = ((statetoid[var] for var in outputvars)...,)
+    affect = eval(generate_affect_function(js, vrj.affect!, outputidxs))
+    VariableRateJump(rate, affect)
+end
+
+function assemble_vrj_expr(js, vrj, statetoid)
     rate   = generate_rate_function(js, vrj.rate)
     outputvars = (convert(Variable,affect.lhs) for affect in vrj.affect!)
     outputidxs = ((statetoid[var] for var in outputvars)...,)
     affect = generate_affect_function(js, vrj.affect!, outputidxs)
-    VariableRateJump(rate, affect)
+    quote
+        rate = $rate
+        affect = $affect
+        VariableRateJump(rate, affect)
+    end
 end
 
 function assemble_crj(js, crj, statetoid)
+    rate   = eval(generate_rate_function(js, crj.rate))
+    outputvars = (convert(Variable,affect.lhs) for affect in crj.affect!)
+    outputidxs = ((statetoid[var] for var in outputvars)...,)
+    affect = eval(generate_affect_function(js, crj.affect!, outputidxs))
+    ConstantRateJump(rate, affect)
+end
+
+function assemble_crj_expr(js, crj, statetoid)
     rate   = generate_rate_function(js, crj.rate)
     outputvars = (convert(Variable,affect.lhs) for affect in crj.affect!)
     outputidxs = ((statetoid[var] for var in outputvars)...,)
     affect = generate_affect_function(js, crj.affect!, outputidxs)
-    ConstantRateJump(rate, affect)
+    quote
+        rate = $rate
+        affect = $affect
+        ConstantRateJump(rate, affect)
+    end
 end
 
 function assemble_maj(js, maj::MassActionJump{U,Vector{Pair{V,W}},Vector{Pair{V2,W2}}},
@@ -124,6 +146,38 @@ function assemble_maj(js, maj::MassActionJump{U,Vector{Pair{V,W}},Vector{Pair{V2
     return maj
 end
 
+function assemble_maj_expr(js, maj::MassActionJump{U,Vector{Pair{V,W}},Vector{Pair{V2,W2}}},
+                      statetoid, subber, invttype) where {U,V,W,V2,W2}
+    sr = maj.scaled_rates
+    if sr isa Operation
+        pval = subber(sr).value
+    elseif sr isa Variable
+        pval = subber(sr()).value
+    else
+        pval = maj.scaled_rates
+    end
+
+    rs = Vector{Pair{Int,W}}()
+    for (spec,stoich) in maj.reactant_stoch
+        if iszero(spec)
+            push!(rs, 0 => stoich)
+        else
+            push!(rs, statetoid[convert(Variable,spec)] => stoich)
+        end
+    end
+    sort!(rs)
+
+    ns = Vector{Pair{Int,W2}}()
+    for (spec,stoich) in maj.net_stoch
+        iszero(spec) && error("Net stoichiometry can not have a species labelled 0.")
+        push!(ns, statetoid[convert(Variable,spec)] => stoich)
+    end
+    sort!(ns)
+
+    maj = MassActionJump(convert(invttype, pval), rs, ns, scale_rates = false)
+    return maj
+end
+
 """
 ```julia
 function DiffEqBase.DiscreteProblem(sys::JumpSystem, u0map, tspan,
@@ -154,6 +208,40 @@ function DiffEqBase.DiscreteProblem(sys::JumpSystem, u0map, tspan::Tuple,
     DiscreteProblem(df, u0, tspan, p; kwargs...)
 end
 
+"""
+```julia
+function DiffEqBase.DiscreteProblemExpr(sys::JumpSystem, u0map, tspan,
+                                    parammap=DiffEqBase.NullParameters; kwargs...)
+```
+
+Generates a black DiscreteProblem for a JumpSystem to utilize as its
+solving `prob.prob`. This is used in the case where there are no ODEs
+and no SDEs associated with the system.
+
+Continuing the example from the [`JumpSystem`](@ref) definition:
+```julia
+using DiffEqBase, DiffEqJump
+u₀map = [S => 999, I => 1, R => 0]
+parammap = [β => .1/1000, γ => .01]
+tspan = (0.0, 250.0)
+dprob = DiscreteProblem(js, u₀map, tspan, parammap)
+```
+"""
+function DiffEqBase.DiscreteProblemExpr(sys::JumpSystem, u0map, tspan::Tuple,
+                                    parammap=DiffEqBase.NullParameters(); kwargs...)
+    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)
+        u0 = $u0
+        p = $p
+        tspan = $tspan
+        df = DiscreteFunction(f, syms=$(Symbol.(states(sys))))
+        DiscreteProblem(df, u0, tspan, p; kwargs...)
+    end
+end
 
 """
 ```julia