remove MassActionJump evals -> substitute

Author: ChrisRackauckas

src/systems/jumps/jumpsystem.jl

@@ -9,24 +9,24 @@ struct JumpSystem <: AbstractSystem
     systems::Vector{JumpSystem}
 end
 
-function JumpSystem(eqs, iv, states, ps; systems = JumpSystem[], 
+function JumpSystem(eqs, iv, states, ps; systems = JumpSystem[],
                                           name = gensym(:JumpSystem))
     JumpSystem(eqs, iv, convert.(Variable, states), convert.(Variable, ps), name, systems)
 end
 
 
 
-generate_rate_function(js, rate) = build_function(rate, states(js), parameters(js), 
-                                                  independent_variable(js), 
+generate_rate_function(js, rate) = build_function(rate, states(js), parameters(js),
+                                                  independent_variable(js),
                                                   expression=Val{false})
 
-generate_affect_function(js, affect, outputidxs) = build_function(affect, states(js), 
-                                                      parameters(js), 
+generate_affect_function(js, affect, outputidxs) = build_function(affect, states(js),
+                                                      parameters(js),
                                                       independent_variable(js),
                                                       expression=Val{false},
-                                                      headerfun=add_integrator_header, 
+                                                      headerfun=add_integrator_header,
                                                       outputidxs=outputidxs)[2]
-function assemble_vrj(js, vrj, statetoid) 
+function assemble_vrj(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)...,)
@@ -42,18 +42,20 @@ function assemble_crj(js, crj, statetoid)
     ConstantRateJump(rate, affect)
 end
 
-function assemble_maj(js, maj::MassActionJump{U,Vector{Pair{V,W}},Vector{Pair{V2,W2}}}, 
-                      statetoid, ptoid, p, pcontext) where {U,V,W,V2,W2}
+function assemble_maj(js, maj::MassActionJump{U,Vector{Pair{V,W}},Vector{Pair{V2,W2}}},
+                      statetoid, ptoid, parammap) where {U,V,W,V2,W2}
     sr = maj.scaled_rates
-    if sr isa Operation || sr isa Variable
-        pval = Base.eval(pcontext, Expr(maj.scaled_rates))
-    else   
+    if sr isa Operation
+        pval = substitute(sr,parammap)
+    elseif sr isa Variable
+        pval = Dict(parammap)[sr()]
+    else
         pval = maj.scaled_rates
     end
-    
+
     rs = Vector{Pair{Int,W}}()
     for (spec,stoich) in maj.reactant_stoch
-        if iszero(spec) 
+        if iszero(spec)
             push!(rs, 0 => stoich)
         else
             push!(rs, statetoid[convert(Variable,spec)] => stoich)
@@ -73,13 +75,13 @@ end
 
 """
 ```julia
-function DiffEqBase.DiscreteProblem(sys::AbstractSystem, u0map, tspan, 
+function DiffEqBase.DiscreteProblem(sys::AbstractSystem, u0map, tspan,
                                     parammap=DiffEqBase.NullParameters; kwargs...)
 ```
 
 Generates a DiscreteProblem from an AbstractSystem
 """
-function DiffEqBase.DiscreteProblem(sys::AbstractSystem, u0map, tspan::Tuple, 
+function DiffEqBase.DiscreteProblem(sys::AbstractSystem, u0map, tspan::Tuple,
                                     parammap=DiffEqBase.NullParameters(); kwargs...)
     u0 = varmap_to_vars(u0map, states(sys))
     p = varmap_to_vars(parammap, parameters(sys))
@@ -99,29 +101,20 @@ function DiffEqJump.JumpProblem(js::JumpSystem, prob, aggregator; kwargs...)
     majs = Vector{MassActionJump}()
     pvars = parameters(js)
     statetoid = Dict(convert(Variable,state) => i for (i,state) in enumerate(states(js)))
-    ptoid     = Dict(convert(Variable,par) => i for (i,par) in enumerate(parameters(js)))
-    
-    # for mass action jumps might need to evaluate parameter expressions
-    # populate dummy module with params as local variables
-    # (for eval-ing parameter expressions)
-    param_context = Module()
-    for (i, pval) in enumerate(prob.p)        
-        psym = Symbol(pvars[i])
-        Base.eval(param_context, :($psym = $pval))
-    end
+    parammap = map(Pair,pvars,prob.p)
 
     for j in equations(js)
         if j isa ConstantRateJump
             push!(crjs, assemble_crj(js, j, statetoid))
         elseif j isa VariableRateJump
             push!(vrjs, assemble_vrj(js, j, statetoid))
         elseif j isa MassActionJump
-            push!(majs, assemble_maj(js, j, statetoid, ptoid, prob.p, param_context))
+            push!(majs, assemble_maj(js, j, statetoid, parammap))
         else
             error("JumpSystems should only contain Constant, Variable or Mass Action Jumps.")
         end
     end
-    ((prob isa DiscreteProblem) && !isempty(vrjs)) && error("Use continuous problems such as an ODEProblem or a SDEProblem with VariableRateJumps")    
+    ((prob isa DiscreteProblem) && !isempty(vrjs)) && error("Use continuous problems such as an ODEProblem or a SDEProblem with VariableRateJumps")
     jset = JumpSet(Tuple(vrjs), Tuple(crjs), nothing, isempty(majs) ? nothing : majs)
     JumpProblem(prob, aggregator, jset)
-end
+end