test rx edge cases

Author: isaacsas

Project.toml

@@ -27,7 +27,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 [compat]
 ArrayInterface = "2.8"
 DiffEqBase = "6.28"
-DiffEqJump = "6.6"
+DiffEqJump = "6.6.2"
 DiffRules = "0.1, 1.0"
 DocStringExtensions = "0.7, 0.8"
 GeneralizedGenerated = "0.1.4, 0.2"

src/systems/jumps/jumpsystem.jl

@@ -44,8 +44,13 @@ 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}
-    pval = Base.eval(pcontext, Expr(maj.scaled_rates))
-
+    sr = maj.scaled_rates
+    if sr isa Operation || sr isa Variable
+        pval = Base.eval(pcontext, Expr(maj.scaled_rates))
+    else   
+        pval = maj.scaled_rates
+    end
+    
     rs = Vector{Pair{Int,W}}()
     for (spec,stoich) in maj.reactant_stoch
         if iszero(spec) 

test/jumpsystem.jl

@@ -105,19 +105,37 @@ jprob = JumpProblem(prob,Direct(),jset, save_positions=(false,false))
 m2 = getmean(jprob,Nsims)
 
 # test JumpSystem solution agrees with direct version
-@test abs(m-m2) ./ m < .01
+@test abs(m-m2)/m < .01
 
 
 # mass action jump tests for SIR model
 maj1 = MassActionJump(2*β/2, [S => 1, I => 1], [S => -1, I => 1])
 maj2 = MassActionJump(γ, [I => 1], [I => -1, R => 1])
-js   = JumpSystem([maj1,maj2], t, [S,I,R], [β,γ])
+js3   = JumpSystem([maj1,maj2], t, [S,I,R], [β,γ])
 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)))
-dprob = DiscreteProblem(js, u₀map, tspan, parammap)
-jprob = JumpProblem(js, dprob, Direct())
+dprob = DiscreteProblem(js3, u₀map, tspan, parammap)
+jprob = JumpProblem(js3, dprob, Direct())
 m3 = getmean(jprob,Nsims)
-@test abs(m2-m3) < .01
+@test abs(m-m3)/m < .01
 
-# mass action jump tests for other reaction types (zero order, second order, decay)
-# TODO
+# mass action jump tests for other reaction types (zero order, decay)
+maj1 = MassActionJump(2.0, [0 => 1], [S => 1])
+maj2 = MassActionJump(γ, [S => 1], [S => -1])
+js4   = JumpSystem([maj1,maj2], t, [S], [β,γ])
+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)))
+dprob = DiscreteProblem(js4, [S => 999], (0,1000.), [β => 100.,γ => .01])
+jprob = JumpProblem(js4, dprob, Direct())
+m4 = getmean(jprob,Nsims)
+@test abs(m4 - 2.0/.01)*.01/2.0 < .01
+
+# test second order rx runs 
+maj1 = MassActionJump(2.0, [0 => 1], [S => 1])
+maj2 = MassActionJump(γ, [S => 2], [S => -1])
+js4   = JumpSystem([maj1,maj2], t, [S], [β,γ])
+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)))
+dprob = DiscreteProblem(js4, [S => 999], (0,1000.), [β => 100.,γ => .01])
+jprob = JumpProblem(js4, dprob, Direct())
+sol = solve(jprob, SSAStepper())