Merge pull request #69 from spcornelius/fix_ma_rateconsts

Allow mass action rate consts involving parameter expressions

Author: isaacsas

src/massaction_jump_utils.jl

@@ -52,7 +52,7 @@ function get_net_stoich(rs, specmap)
 end
 
 # given a ReactionStruct and a species map construct a MassActionJump
-function make_majump(rs, specmap, ratemap, params)
+function make_majump(rs, specmap, ratemap, params, param_context)
     reactant_stoich = get_substrate_stoich(rs, specmap)
     net_stoich      = get_net_stoich(rs, specmap)
     if isempty(net_stoich)
@@ -62,7 +62,8 @@ function make_majump(rs, specmap, ratemap, params)
     if typeof(rs.rate_org) == Symbol
         rateconst = params[ratemap[rs.rate_org]]
     elseif typeof(rs.rate_org) == Expr
-        rateconst = eval(rs.rate_org)
+        # eval in param_context, in case Expr depends on params
+        rateconst = Base.eval(param_context, rs.rate_org)
     elseif typeof(rs.rate_org) <: Number
         rateconst = rs.rate_org
     else
@@ -79,9 +80,16 @@ function network_to_jumpset(rn, specmap, ratemap, params)
     majumpvec    = Vector{typeof(empty_majump)}()
     cjumpvec     = Vector{ConstantRateJump}()
 
+    # populate dummy module with params as local variables
+    # (for eval-ing parameter expressions)
+    param_context = Module()
+    for (param, index) in ratemap
+        Base.eval(param_context, :($param = $(params[index])))
+    end
+
     for (i,rs) in enumerate(rn.reactions)
         if rs.is_pure_mass_action
-            push!(majumpvec, make_majump(rs, specmap, ratemap, params))
+            push!(majumpvec, make_majump(rs, specmap, ratemap, params, param_context))
         else
             push!(cjumpvec, rn.jumps[i])
         end

test/mass_act_jump_tests.jl

@@ -113,3 +113,22 @@ for method in algs
     jump_prob = JumpProblem(prob, method, network)
     sol = solve(jump_prob,SSAStepper());
 end
+
+# make sure problem instantiation works when mass action jumps
+# rate consts are a nontrivial expression on the params
+network = @reaction_network begin
+    1, X --> 2*X
+    1/K, 2X --> X
+end K
+p = [1000.]
+prob = DiscreteProblem([500.], (0.,100.), p)
+jump_prob = JumpProblem(prob, Direct(), network)
+
+
+# same as above
+network = @reaction_network begin
+    p1*p2, X --> Y
+end p1 p2
+p = [1.,2.]
+prob = DiscreteProblem([10.,10.],(0.,10.), p)
+jump_prob = JumpProblem(prob, Direct(), network)