Merge branch 'jump-construction-fix'

Author: ChrisRackauckas

src/systems/jumps/jumpsystem.jl

@@ -91,13 +91,18 @@ function assemble_crj(js, crj, statetoid)
 end
 
 function assemble_maj(js, maj::MassActionJump{U,Vector{Pair{V,W}},Vector{Pair{V2,W2}}},
-                      statetoid, parammap) where {U,V,W,V2,W2}
+                      statetoid, parammap, pcontext) where {U,V,W,V2,W2}
+    
     sr = maj.scaled_rates
-    if sr isa Operation
-        pval = simplify(substitute(sr,parammap)).value
+    if sr isa Operation 
+        if isempty(sr.args)
+            pval = parammap[sr.op]
+        else
+            pval = Base.eval(pcontext, Expr(maj.scaled_rates))
+        end
     elseif sr isa Variable
-        pval = Dict(parammap)[sr()]
-    else
+        pval = parammap[sr]
+    else   
         pval = maj.scaled_rates
     end
 
@@ -164,10 +169,20 @@ 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  = map((x,y)->Pair(x(),y), parameters(js), prob.p)
+    parammap  = Dict(convert(Variable,param) => prob.p[i] for (i,param) in enumerate(parameters(js)))
     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) for j in eqs.x[1]]
+    majs = MassActionJump[assemble_maj(js, j, statetoid, parammap, param_context) 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") 

src/systems/reaction/reactionsystem.jl

@@ -244,8 +244,10 @@ function assemble_jumps(rs)
         haveivdep = any(var -> isequal(rs.iv,convert(Variable,var)), rxvars)
         if ismassaction(rx, rs; rxvars=rxvars, haveivdep=haveivdep)
             reactant_stoch = isempty(rx.substoich) ? [0 => 1] : [var2op(sub.op) => stoich for (sub,stoich) in zip(rx.substrates,rx.substoich)]
+            coef           = isempty(rx.substoich) ? one(eltype(rx.substoich)) : prod(stoich -> factorial(stoich), rx.substoich)
+            rate           = isone(coef) ? rx.rate : rx.rate/coef
             net_stoch      = [Pair(var2op(p[1]),p[2]) for p in rx.netstoich]
-            push!(eqs, MassActionJump(rx.rate, reactant_stoch, net_stoch))
+            push!(eqs, MassActionJump(rate, reactant_stoch, net_stoch, scale_rates=false))
         else
             rl     = jumpratelaw(rx, rxvars=rxvars)
             affect = Vector{Equation}()

test/reactionsystem.jl

@@ -115,9 +115,15 @@ 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 = map((x,y)->Pair(x(),y),parameters(js),pars)
+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
 for i = 1:14
-  maj = MT.assemble_maj(js, js.eqs[i], statetoid,parammap)
+  maj = MT.assemble_maj(js, js.eqs[i], statetoid, parammap, param_context)
   @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