Merge pull request #432 from SciML/s/substituter

`substituter` to avoid Dict re-construction

Author: ChrisRackauckas

src/direct.jl

@@ -43,6 +43,11 @@ function simplified_expr(O::Operation)
     isempty(O.args) && return O.op.name
     return Expr(:call, Symbol(O.op), simplified_expr.(O.args)...)
   end
+  if O.op === (^)
+      if length(O.args) > 1  && O.args[2] isa Constant && O.args[2].value < 0
+          return Expr(:call, :^, Expr(:call, :inv, simplified_expr(O.args[1])), -(O.args[2].value))
+      end
+  end
   return Expr(:call, Symbol(O.op), simplified_expr.(O.args)...)
 end
 

src/systems/jumps/jumpsystem.jl

@@ -91,18 +91,13 @@ function assemble_crj(js, crj, statetoid)
 end
 
 function assemble_maj(js, maj::MassActionJump{U,Vector{Pair{V,W}},Vector{Pair{V2,W2}}},
-                      statetoid, parammap, pcontext) where {U,V,W,V2,W2}
-    
+                      statetoid, subber) where {U,V,W,V2,W2}
     sr = maj.scaled_rates
-    if sr isa Operation 
-        if isempty(sr.args)
-            pval = parammap[sr.op]
-        else
-            pval = Base.eval(pcontext, Expr(maj.scaled_rates))
-        end
+    if sr isa Operation
+        pval = subber(sr).value
     elseif sr isa Variable
-        pval = parammap[sr]
-    else   
+        pval = subber(sr()).value
+    else
         pval = maj.scaled_rates
     end
 
@@ -169,20 +164,11 @@ sol = solve(jprob, SSAStepper())
 function DiffEqJump.JumpProblem(js::JumpSystem, prob, aggregator; kwargs...)
 
     statetoid = Dict(convert(Variable,state) => i for (i,state) in enumerate(states(js)))
-    parammap  = Dict(convert(Variable,param) => prob.p[i] for (i,param) in enumerate(parameters(js)))
+    parammap  = map((x,y)->Pair(x(),y), parameters(js), prob.p)
     eqs       = equations(js)
-    
-    # for mass action jumps might need to evaluate parameter expressions
-    # populate dummy module with params as local variables
-    # (for eval-ing parameter expressions)
-    pvars         = parameters(js)
-    param_context = Module()
-    for (i, pval) in enumerate(prob.p)        
-        psym = Symbol(pvars[i])
-        Base.eval(param_context, :($psym = $pval))
-    end
 
-    majs = MassActionJump[assemble_maj(js, j, statetoid, parammap, param_context) for j in eqs.x[1]]
+    subber = substituter(first.(parammap), last.(parammap))
+    majs = MassActionJump[assemble_maj(js, j, statetoid, subber) for j in eqs.x[1]]
     crjs = ConstantRateJump[assemble_crj(js, j, statetoid) for j in eqs.x[2]]
     vrjs = VariableRateJump[assemble_vrj(js, j, statetoid) for j in eqs.x[3]]
     ((prob isa DiscreteProblem) && !isempty(vrjs)) && error("Use continuous problems such as an ODEProblem or a SDEProblem with VariableRateJumps") 
@@ -225,4 +211,4 @@ function modified_states!(mstates, jump::MassActionJump, sts)
     for (state,stoich) in jump.net_stoch
         (state.op in sts) && push!(mstates, state)
     end
-end
+end

src/utils.jl

@@ -121,12 +121,16 @@ substitute(expr::Operation, s::Pair) = _substitute(expr, [s[1]], [s[2]])
 substitute(expr::Operation, dict::Dict) = _substitute(expr, keys(dict), values(dict))
 substitute(expr::Operation, s::Vector) = _substitute(expr, first.(s), last.(s))
 
-function _substitute(expr, ks, vs)
-    _substitute(expr, Dict(map(Pair, map(to_symbolic, ks), map(to_symbolic, vs))))
+function _substitute(ks, vs)
+    expr -> _substitute(expr, Dict(map(Pair, map(to_symbolic, ks), map(to_symbolic, vs))))
 end
 
-function _substitute(expr, dict::Dict)
-    simplify(SymbolicUtils.substitute(expr, dict))
+function substituter(ks, vs)
+    dict = Dict(map(Pair, map(to_symbolic, ks), map(to_symbolic, vs)))
+    expr -> to_mtk(SymbolicUtils.simplify(SymbolicUtils.substitute(expr, dict)))
 end
 
-@deprecate substitute_expr!(expr,s) substitute(expr,s)
+_substitute(expr, ks, vs) = substituter(ks, vs)(expr)
+
+@deprecate substitute_expr!(expr,s) substitute(expr,s)
+

test/reactionsystem.jl

@@ -115,15 +115,9 @@ jumps[19] = VariableRateJump((u,p,t) -> p[19]*u[1]*t, integrator -> (integrator.
 jumps[20] = VariableRateJump((u,p,t) -> p[20]*t*u[1]*binomial(u[2],2)*u[3], integrator -> (integrator.u[2] -= 2; integrator.u[3] -= 1; integrator.u[4] += 2))
 
 statetoid = Dict(convert(Variable,state) => i for (i,state) in enumerate(states(js)))
-parammap  = Dict(convert(Variable,param) => pars[i] for (i,param) in enumerate(parameters(js)))
-pvars = parameters(js)
-param_context = Module()
-for (i, pval) in enumerate(pars)        
-    psym = Symbol(pvars[i])
-    Base.eval(param_context, :($psym = $pval))
-end
+parammap = map((x,y)->Pair(x(),y),parameters(js),pars)
 for i = 1:14
-  maj = MT.assemble_maj(js, js.eqs[i], statetoid, parammap, param_context)
+  maj = MT.assemble_maj(js, js.eqs[i], statetoid, ModelingToolkit.substituter(first.(parammap), last.(parammap)))
   @test abs(jumps[i].scaled_rates - maj.scaled_rates) < 100*eps()
   @test jumps[i].reactant_stoch == maj.reactant_stoch
   @test jumps[i].net_stoch == maj.net_stoch