add test for generated jump system

Author: isaacsas

test/jumpsystem.jl

@@ -57,42 +57,50 @@ j₃ = ConstantRateJump(rate₃,affect₃)
 js2 = JumpSystem([j₁,j₃], t, [S,I,R], [β,γ])
 u₀ = [999,1,0]; p = (0.1/1000,0.01); tspan = (0.,250.)
 dprob = DiscreteProblem(u₀,tspan,p)
-jprob = JumpProblem(js2, dprob, Direct())
-sol = solve(jprob, SSAStepper())
-plot(sol)
+jprob = JumpProblem(js2, dprob, Direct(), save_positions=(false,false))
+Nsims = 10000
+function getmean(jprob,Nsims)
+  m = 0.0
+  for i = 1:Nsims
+    sol = solve(jprob, SSAStepper())
+    m += sol[end,end]
+  end
+  m/Nsims
+end
+m = getmean(jprob,Nsims)
 
-# test the MT JumpProblem rates/affects are correct
-# rate2(u,p,t) = 0.01u[2]
-# jump2 = ConstantRateJump(rate2,affect2!)
-# mtjumps = jprob.discrete_jump_aggregation
-# @test abs(mtjumps.rates[1](u,p,tf) - jump1.rate(u,p,tf)) < 10*eps()
-# @test abs(mtjumps.rates[2](u,p,tf) - jump2.rate(u,p,tf)) < 10*eps()
-# mtjumps.affects![1](mtintegrator)
-# jump1.affect!(integrator)
-# @test all(integrator.u .== mtintegrator.u) 
-# mtintegrator.u .= u; integrator.u .= u
-# mtjumps.affects![2](mtintegrator)
-# jump2.affect!(integrator)
-# @test all(integrator.u .== mtintegrator.u)
+#test the MT JumpProblem rates/affects are correct
+rate2(u,p,t) = 0.01u[2]
+jump2 = ConstantRateJump(rate2,affect2!)
+mtjumps = jprob.discrete_jump_aggregation
+@test abs(mtjumps.rates[1](u,p,tf) - jump1.rate(u,p,tf)) < 10*eps()
+@test abs(mtjumps.rates[2](u,p,tf) - jump2.rate(u,p,tf)) < 10*eps()
+mtjumps.affects![1](mtintegrator)
+jump1.affect!(integrator)
+@test all(integrator.u .== mtintegrator.u) 
+mtintegrator.u .= u; integrator.u .= u
+mtjumps.affects![2](mtintegrator)
+jump2.affect!(integrator)
+@test all(integrator.u .== mtintegrator.u)
 
-# # # direct vers
-# p = (0.1/1000,0.01)
-# prob = DiscreteProblem([999,1,0],(0.0,250.0),p)
-# r1(u,p,t) = (0.1/1000.0)*u[1]*u[2]
-# function a1!(integrator)
-#   integrator.u[1] -= 1
-#   integrator.u[2] += 1
-# end
-# j1 = ConstantRateJump(r1,a1!)
-# r2(u,p,t) = 0.01u[2]
-# function a2!(integrator)
-#   integrator.u[2] -= 1
-#   integrator.u[3] += 1
-# end
-# j2 = ConstantRateJump(r2,a2!)
-# jset = JumpSet((),(j1,j2),nothing,nothing)
-# jprob = JumpProblem(prob,Direct(),jset)
-# sol2 = solve(jprob, SSAStepper())
+# direct vers
+p = (0.1/1000,0.01)
+prob = DiscreteProblem([999,1,0],(0.0,250.0),p)
+r1(u,p,t) = (0.1/1000.0)*u[1]*u[2]
+function a1!(integrator)
+  integrator.u[1] -= 1
+  integrator.u[2] += 1
+end
+j1 = ConstantRateJump(r1,a1!)
+r2(u,p,t) = 0.01u[2]
+function a2!(integrator)
+  integrator.u[2] -= 1
+  integrator.u[3] += 1
+end
+j2 = ConstantRateJump(r2,a2!)
+jset = JumpSet((),(j1,j2),nothing,nothing)
+jprob = JumpProblem(prob,Direct(),jset, save_positions=(false,false))
+m2 = getmean(jprob,Nsims)
 
-# using Plots
-# plot(sol)
+# test JumpSystem solution agrees with direct version
+@test abs(m-m2) ./ m < .01