Skip to content

Commit 7396a2c

Browse files
committed
indexing works with JumpProcesses 9.13.2
1 parent 64fb99d commit 7396a2c

File tree

1 file changed

+37
-6
lines changed

1 file changed

+37
-6
lines changed

Diff for: test/jumpsystem.jl

+37-6
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,5 @@
1-
using ModelingToolkit, DiffEqBase, JumpProcesses, Test, LinearAlgebra, StableRNGs
1+
using ModelingToolkit, DiffEqBase, JumpProcesses, Test, LinearAlgebra
2+
using Random, StableRNGs
23
using OrdinaryDiffEq
34
using ModelingToolkit: t_nounits as t, D_nounits as D
45
MT = ModelingToolkit
@@ -299,13 +300,43 @@ end
299300

300301
# basic VariableRateJump test
301302
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())
303+
N = 1000 # number of simulations for testing solve accuracy
304+
Random.seed!(rng, 1111)
305+
@variables A(t) B(t) C(t)
306+
@parameters k
307+
vrj = VariableRateJump(k * (sin(t) + 1), [A ~ A + 1, C ~ C + 2])
308+
js = complete(JumpSystem([vrj], t, [A,C], [k]; name = :js, observed = [B ~ C*A]))
309+
oprob = ODEProblem(js, [A => 0, C => 0], (0.0, 10.0), [k => 1.0])
310+
jprob = JumpProblem(js, oprob, Direct(); rng)
307311
sol = solve(jprob, Tsit5())
308312

313+
# test observed and symbolic indexing work
314+
@test all(sol[:A] .* sol[:C] .== sol[:B])
315+
316+
dt = 1.0
317+
tv = range(0.0, 10.0; step = 1.0)
318+
cmean = zeros(11)
319+
for n in 1:N
320+
sol = solve(jprob, Tsit5(); save_everystep = false, saveat = dt)
321+
cmean += Array(sol(tv; idxs = :C))
322+
end
323+
cmean ./= N
309324

325+
vrjrate(u, p, t) = p[1] * (sin(t) + 1)
326+
function vrjaffect!(integ)
327+
integ.u[1] += 1
328+
integ.u[2] += 2
329+
nothing
330+
end
331+
vrj2 = VariableRateJump(vrjrate, vrjaffect!)
332+
oprob2 = ODEProblem((du,u,p,t) -> (du .= 0; nothing), [0, 0], (0.0, 10.0), (1.0,))
333+
jprob2 = JumpProblem(oprob2, Direct(), vrj2; rng)
334+
cmean2 = zeros(11)
335+
for n in 1:N
336+
sol2 = solve(jprob2, Tsit5(); saveat = dt)
337+
cmean2 += Array(sol2(tv; idxs = 2))
338+
end
339+
cmean2 ./= N
310340

341+
@test all( abs.(cmean .- cmean2) .<= .05 .* cmean)
311342
end

0 commit comments

Comments
 (0)