1
1
using ModelingToolkit, DiffEqBase, JumpProcesses, Test, LinearAlgebra, StableRNGs
2
+ using OrdinaryDiffEq
2
3
using ModelingToolkit: t_nounits as t, D_nounits as D
3
4
MT = ModelingToolkit
4
5
@@ -67,7 +68,7 @@ tspan = (0.0, 250.0);
67
68
u₀map = [S => 999 , I => 1 , R => 0 ]
68
69
parammap = [β => 0.1 / 1000 , γ => 0.01 ]
69
70
dprob = DiscreteProblem (js2, u₀map, tspan, parammap)
70
- jprob = JumpProblem (js2, dprob, Direct (), save_positions = (false , false ), rng)
71
+ jprob = JumpProblem (js2, dprob, Direct (); save_positions = (false , false ), rng)
71
72
Nsims = 30000
72
73
function getmean (jprob, Nsims; use_stepper = true )
73
74
m = 0.0
@@ -89,13 +90,13 @@ obs = [S2 ~ 2 * S]
89
90
@named js2b = JumpSystem ([j₁, j₃], t, [S, I, R], [β, γ], observed = obs)
90
91
js2b = complete (js2b)
91
92
dprob = DiscreteProblem (js2b, u₀map, tspan, parammap)
92
- jprob = JumpProblem (js2b, dprob, Direct (), save_positions = (false , false ), rng)
93
- sol = solve (jprob, SSAStepper (), saveat = tspan[2 ] / 10 )
93
+ jprob = JumpProblem (js2b, dprob, Direct (); save_positions = (false , false ), rng)
94
+ sol = solve (jprob, SSAStepper (); saveat = tspan[2 ] / 10 )
94
95
@test all (2 .* sol[S] .== sol[S2])
95
96
96
97
# test save_positions is working
97
- jprob = JumpProblem (js2, dprob, Direct (), save_positions = (false , false ), rng)
98
- sol = solve (jprob, SSAStepper (), saveat = 1.0 )
98
+ jprob = JumpProblem (js2, dprob, Direct (); save_positions = (false , false ), rng)
99
+ sol = solve (jprob, SSAStepper (); saveat = 1.0 )
99
100
@test all ((sol. t) .== collect (0.0 : tspan[2 ]))
100
101
101
102
# test the MT JumpProblem rates/affects are correct
@@ -129,7 +130,7 @@ function a2!(integrator)
129
130
end
130
131
j2 = ConstantRateJump (r2, a2!)
131
132
jset = JumpSet ((), (j1, j2), nothing , nothing )
132
- jprob = JumpProblem (prob, Direct (), jset, save_positions = (false , false ), rng)
133
+ jprob = JumpProblem (prob, Direct (), jset; save_positions = (false , false ), rng)
133
134
m2 = getmean (jprob, Nsims)
134
135
135
136
# test JumpSystem solution agrees with direct version
@@ -141,17 +142,17 @@ maj2 = MassActionJump(γ, [I => 1], [I => -1, R => 1])
141
142
@named js3 = JumpSystem ([maj1, maj2], t, [S, I, R], [β, γ])
142
143
js3 = complete (js3)
143
144
dprob = DiscreteProblem (js3, u₀map, tspan, parammap)
144
- jprob = JumpProblem (js3, dprob, Direct (), rng)
145
+ jprob = JumpProblem (js3, dprob, Direct (); rng)
145
146
m3 = getmean (jprob, Nsims)
146
147
@test abs (m - m3) / m < 0.01
147
148
148
149
# maj jump test with various dep graphs
149
150
@named js3b = JumpSystem ([maj1, maj2], t, [S, I, R], [β, γ])
150
151
js3b = complete (js3b)
151
- jprobb = JumpProblem (js3b, dprob, NRM (), rng)
152
+ jprobb = JumpProblem (js3b, dprob, NRM (); rng)
152
153
m4 = getmean (jprobb, Nsims)
153
154
@test abs (m - m4) / m < 0.01
154
- jprobc = JumpProblem (js3b, dprob, RSSA (), rng)
155
+ jprobc = JumpProblem (js3b, dprob, RSSA (); rng)
155
156
m4 = getmean (jprobc, Nsims)
156
157
@test abs (m - m4) / m < 0.01
157
158
@@ -161,7 +162,7 @@ maj2 = MassActionJump(γ, [S => 1], [S => -1])
161
162
@named js4 = JumpSystem ([maj1, maj2], t, [S], [β, γ])
162
163
js4 = complete (js4)
163
164
dprob = DiscreteProblem (js4, [S => 999 ], (0 , 1000.0 ), [β => 100.0 , γ => 0.01 ])
164
- jprob = JumpProblem (js4, dprob, Direct (), rng)
165
+ jprob = JumpProblem (js4, dprob, Direct (); rng)
165
166
m4 = getmean (jprob, Nsims)
166
167
@test abs (m4 - 2.0 / 0.01 ) * 0.01 / 2.0 < 0.01
167
168
@@ -171,7 +172,7 @@ maj2 = MassActionJump(γ, [S => 2], [S => -1])
171
172
@named js4 = JumpSystem ([maj1, maj2], t, [S], [β, γ])
172
173
js4 = complete (js4)
173
174
dprob = DiscreteProblem (js4, [S => 999 ], (0 , 1000.0 ), [β => 100.0 , γ => 0.01 ])
174
- jprob = JumpProblem (js4, dprob, Direct (), rng)
175
+ jprob = JumpProblem (js4, dprob, Direct (); rng)
175
176
sol = solve (jprob, SSAStepper ());
176
177
177
178
# issue #819
295
296
@test jprob. aggregator isa algtype
296
297
end
297
298
end
299
+
300
+ # basic VariableRateJump test
301
+ let
302
+ @variables A (t)
303
+ vrj = VariableRateJump (sin (t) + 1 , [A ~ A + 1 ])
304
+ js = complete (JumpSystem ([vrj], t, [A], []; name = :js ))
305
+ oprob = ODEProblem (js, [A => 0 ], (0.0 , 10.0 ))
306
+ jprob = JumpProblem (js, oprob, Direct ())
307
+ sol = solve (jprob, Tsit5 ())
308
+
309
+
310
+
311
+ end
0 commit comments