1- using ModelingToolkit, DiffEqBase, JumpProcesses, Test, LinearAlgebra
1+ using ModelingToolkit, DiffEqBase, JumpProcesses, Test, LinearAlgebra, StableRNGs
22MT = ModelingToolkit
33
4+ rng = StableRNG (12345 )
5+
46# basic MT SIR model with tweaks
57@parameters β γ t
68@constants h = 1
@@ -63,7 +65,7 @@ tspan = (0.0, 250.0);
6365u₀map = [S => 999 , I => 1 , R => 0 ]
6466parammap = [β => 0.1 / 1000 , γ => 0.01 ]
6567dprob = DiscreteProblem (js2, u₀map, tspan, parammap)
66- jprob = JumpProblem (js2, dprob, Direct (), save_positions = (false , false ))
68+ jprob = JumpProblem (js2, dprob, Direct (), save_positions = (false , false ), rng = rng )
6769Nsims = 30000
6870function getmean (jprob, Nsims)
6971 m = 0.0
@@ -79,13 +81,13 @@ m = getmean(jprob, Nsims)
7981obs = [S2 ~ 2 * S]
8082@named js2b = JumpSystem ([j₁, j₃], t, [S, I, R], [β, γ], observed = obs)
8183dprob = DiscreteProblem (js2b, u₀map, tspan, parammap)
82- jprob = JumpProblem (js2b, dprob, Direct (), save_positions = (false , false ))
84+ jprob = JumpProblem (js2b, dprob, Direct (), save_positions = (false , false ), rng = rng )
8385sol = solve (jprob, SSAStepper (), saveat = tspan[2 ] / 10 )
8486@test all (2 .* sol[S] .== sol[S2])
8587
8688# test save_positions is working
8789
88- jprob = JumpProblem (js2, dprob, Direct (), save_positions = (false , false ))
90+ jprob = JumpProblem (js2, dprob, Direct (), save_positions = (false , false ), rng = rng )
8991sol = solve (jprob, SSAStepper (), saveat = 1.0 )
9092@test all ((sol. t) .== collect (0.0 : tspan[2 ]))
9193
@@ -120,7 +122,7 @@ function a2!(integrator)
120122end
121123j2 = ConstantRateJump (r2, a2!)
122124jset = JumpSet ((), (j1, j2), nothing , nothing )
123- jprob = JumpProblem (prob, Direct (), jset, save_positions = (false , false ))
125+ jprob = JumpProblem (prob, Direct (), jset, save_positions = (false , false ), rng = rng )
124126m2 = getmean (jprob, Nsims)
125127
126128# test JumpSystem solution agrees with direct version
@@ -131,16 +133,16 @@ maj1 = MassActionJump(2 * β / 2, [S => 1, I => 1], [S => -1, I => 1])
131133maj2 = MassActionJump (γ, [I => 1 ], [I => - 1 , R => 1 ])
132134@named js3 = JumpSystem ([maj1, maj2], t, [S, I, R], [β, γ])
133135dprob = DiscreteProblem (js3, u₀map, tspan, parammap)
134- jprob = JumpProblem (js3, dprob, Direct ())
136+ jprob = JumpProblem (js3, dprob, Direct (), rng = rng )
135137m3 = getmean (jprob, Nsims)
136138@test abs (m - m3) / m < 0.01
137139
138140# maj jump test with various dep graphs
139141@named js3b = JumpSystem ([maj1, maj2], t, [S, I, R], [β, γ])
140- jprobb = JumpProblem (js3b, dprob, NRM ())
142+ jprobb = JumpProblem (js3b, dprob, NRM (), rng = rng )
141143m4 = getmean (jprobb, Nsims)
142144@test abs (m - m4) / m < 0.01
143- jprobc = JumpProblem (js3b, dprob, RSSA ())
145+ jprobc = JumpProblem (js3b, dprob, RSSA (), rng = rng )
144146m4 = getmean (jprobc, Nsims)
145147@test abs (m - m4) / m < 0.01
146148
@@ -149,7 +151,7 @@ maj1 = MassActionJump(2.0, [0 => 1], [S => 1])
149151maj2 = MassActionJump (γ, [S => 1 ], [S => - 1 ])
150152@named js4 = JumpSystem ([maj1, maj2], t, [S], [β, γ])
151153dprob = DiscreteProblem (js4, [S => 999 ], (0 , 1000.0 ), [β => 100.0 , γ => 0.01 ])
152- jprob = JumpProblem (js4, dprob, Direct ())
154+ jprob = JumpProblem (js4, dprob, Direct (), rng = rng )
153155m4 = getmean (jprob, Nsims)
154156@test abs (m4 - 2.0 / 0.01 ) * 0.01 / 2.0 < 0.01
155157
@@ -158,7 +160,7 @@ maj1 = MassActionJump(2.0, [0 => 1], [S => 1])
158160maj2 = MassActionJump (γ, [S => 2 ], [S => - 1 ])
159161@named js4 = JumpSystem ([maj1, maj2], t, [S], [β, γ])
160162dprob = DiscreteProblem (js4, [S => 999 ], (0 , 1000.0 ), [β => 100.0 , γ => 0.01 ])
161- jprob = JumpProblem (js4, dprob, Direct ())
163+ jprob = JumpProblem (js4, dprob, Direct (), rng = rng )
162164sol = solve (jprob, SSAStepper ());
163165
164166# issue #819
@@ -179,7 +181,7 @@ p = [k1 => 2.0, k2 => 0.0, k3 => 0.5]
179181u₀ = [A => 100 , B => 0 ]
180182tspan = (0.0 , 2000.0 )
181183dprob = DiscreteProblem (js5, u₀, tspan, p)
182- jprob = JumpProblem (js5, dprob, Direct (), save_positions = (false , false ))
184+ jprob = JumpProblem (js5, dprob, Direct (), save_positions = (false , false ), rng = rng )
183185@test all (jprob. massaction_jump. scaled_rates .== [1.0 , 0.0 ])
184186
185187pcondit (u, t, integrator) = t == 1000.0
0 commit comments