From 8bf9140a2244a7b69be039dd8b6a19a327593a53 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Sun, 1 Sep 2024 20:27:41 -0400 Subject: [PATCH 1/8] fix vrj reinitialization bug --- src/solve.jl | 4 ++-- test/variable_rate.jl | 42 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 44 insertions(+), 2 deletions(-) diff --git a/src/solve.jl b/src/solve.jl index bc9ed6c8..bdcb03d3 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -58,7 +58,7 @@ function resetted_jump_problem(_jump_prob, seed) if !isempty(jump_prob.variable_jumps) @assert jump_prob.prob.u0 isa ExtendedJumpArray - jump_prob.prob.u0.jump_u .= -randexp.(_jump_prob.rng, eltype(_jump_prob.prob.tspan)) + @. jump_prob.prob.u0.jump_u = -randexp(_jump_prob.rng, eltype(_jump_prob.prob.tspan)) end jump_prob end @@ -70,6 +70,6 @@ function reset_jump_problem!(jump_prob, seed) if !isempty(jump_prob.variable_jumps) @assert jump_prob.prob.u0 isa ExtendedJumpArray - jump_prob.prob.u0.jump_u .= -randexp.(jump_prob.rng, eltype(jump_prob.prob.tspan)) + @. jump_prob.prob.u0.jump_u = -randexp(jump_prob.rng, eltype(jump_prob.prob.tspan)) end end diff --git a/test/variable_rate.jl b/test/variable_rate.jl index 3f2b3853..ad6c24d0 100644 --- a/test/variable_rate.jl +++ b/test/variable_rate.jl @@ -229,3 +229,45 @@ let solve(jprob, SSAStepper()) end end + +# test u0 resets correctly +let + b = 2.0 + d = 1.0 + n0 = 1 + tspan = (0.0, 4.0) + Nsims = 10 + u0 = [n0] + p = [b,d] + + function ode_fxn(du, u, p, t) + du .= 0 + nothing + end + b_rate(u, p, t) = (u[1] * p[1]) + function birth!(integrator) + integrator.u[1] += 1 + nothing + end + b_jump = VariableRateJump(b_rate, birth!) + + d_rate(u, p, t) = (u[1] * p[2]) + function death!(integrator) + integrator.u[1] -= 1 + nothing + end + d_jump = VariableRateJump(d_rate, death!) + + ode_prob = ODEProblem(ode_fxn, u0, tspan, p) + sjm_prob = JumpProblem(ode_prob, b_jump, d_jump; rng) + u0old = copy(sjm_prob.prob.u0.jump_u) + dt = .1 + tsave = range(tspan[1], tspan[2]; step = dt) + umean = zeros(length(tsave)) + for i in 1:Nsims + sol = solve(sjm_prob, Tsit5(); saveat = dt) + @test allunique(sjm_prob.prob.u0.jump_u) + @test all(u0old != sjm_prob.prob.u0.jump_u) + u0old .= sjm_prob.prob.u0.jump_u + end +end \ No newline at end of file From 4a6cd06a261ed591334cb3427e28ec5ca236d5b9 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Sun, 1 Sep 2024 20:30:18 -0400 Subject: [PATCH 2/8] format and drop unneeded parts of test --- src/solve.jl | 3 ++- test/variable_rate.jl | 21 +++++++++------------ 2 files changed, 11 insertions(+), 13 deletions(-) diff --git a/src/solve.jl b/src/solve.jl index bdcb03d3..d23775c1 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -58,7 +58,8 @@ function resetted_jump_problem(_jump_prob, seed) if !isempty(jump_prob.variable_jumps) @assert jump_prob.prob.u0 isa ExtendedJumpArray - @. jump_prob.prob.u0.jump_u = -randexp(_jump_prob.rng, eltype(_jump_prob.prob.tspan)) + @. jump_prob.prob.u0.jump_u = -randexp( + _jump_prob.rng, eltype(_jump_prob.prob.tspan)) end jump_prob end diff --git a/test/variable_rate.jl b/test/variable_rate.jl index ad6c24d0..81351b88 100644 --- a/test/variable_rate.jl +++ b/test/variable_rate.jl @@ -234,11 +234,11 @@ end let b = 2.0 d = 1.0 - n0 = 1 + n0 = 1 tspan = (0.0, 4.0) - Nsims = 10 + Nsims = 10 u0 = [n0] - p = [b,d] + p = [b, d] function ode_fxn(du, u, p, t) du .= 0 @@ -246,28 +246,25 @@ let end b_rate(u, p, t) = (u[1] * p[1]) function birth!(integrator) - integrator.u[1] += 1 + integrator.u[1] += 1 nothing end b_jump = VariableRateJump(b_rate, birth!) - + d_rate(u, p, t) = (u[1] * p[2]) function death!(integrator) integrator.u[1] -= 1 nothing end d_jump = VariableRateJump(d_rate, death!) - + ode_prob = ODEProblem(ode_fxn, u0, tspan, p) sjm_prob = JumpProblem(ode_prob, b_jump, d_jump; rng) u0old = copy(sjm_prob.prob.u0.jump_u) - dt = .1 - tsave = range(tspan[1], tspan[2]; step = dt) - umean = zeros(length(tsave)) for i in 1:Nsims - sol = solve(sjm_prob, Tsit5(); saveat = dt) + sol = solve(sjm_prob, Tsit5(); saveat = tspan[2]) @test allunique(sjm_prob.prob.u0.jump_u) @test all(u0old != sjm_prob.prob.u0.jump_u) u0old .= sjm_prob.prob.u0.jump_u - end -end \ No newline at end of file + end +end From 3bd2cac368007bbb1cbd37d2728bb81f0bab8544 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Sun, 1 Sep 2024 20:31:02 -0400 Subject: [PATCH 3/8] drop bad format --- src/solve.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/solve.jl b/src/solve.jl index d23775c1..b9f552a4 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -58,8 +58,8 @@ function resetted_jump_problem(_jump_prob, seed) if !isempty(jump_prob.variable_jumps) @assert jump_prob.prob.u0 isa ExtendedJumpArray - @. jump_prob.prob.u0.jump_u = -randexp( - _jump_prob.rng, eltype(_jump_prob.prob.tspan)) + @. jump_prob.prob.u0.jump_u = -randexp(_jump_prob.rng, + eltype(_jump_prob.prob.tspan)) end jump_prob end From 94e4558103bc4f4dbe7fa66786607655aa78eb3c Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Sun, 1 Sep 2024 20:31:40 -0400 Subject: [PATCH 4/8] format --- src/solve.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solve.jl b/src/solve.jl index b9f552a4..c64046c7 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -58,7 +58,7 @@ function resetted_jump_problem(_jump_prob, seed) if !isempty(jump_prob.variable_jumps) @assert jump_prob.prob.u0 isa ExtendedJumpArray - @. jump_prob.prob.u0.jump_u = -randexp(_jump_prob.rng, + @. jump_prob.prob.u0.jump_u = -randexp(_jump_prob.rng, eltype(_jump_prob.prob.tspan)) end jump_prob From 103c487c31a6cdea178e0912bf3597f2256607b9 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Sun, 1 Sep 2024 20:43:41 -0400 Subject: [PATCH 5/8] accuracy test --- test/variable_rate.jl | 47 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 47 insertions(+) diff --git a/test/variable_rate.jl b/test/variable_rate.jl index 81351b88..64b3873f 100644 --- a/test/variable_rate.jl +++ b/test/variable_rate.jl @@ -1,4 +1,5 @@ using DiffEqBase, JumpProcesses, OrdinaryDiffEq, StochasticDiffEq, Test +using Random using StableRNGs rng = StableRNG(12345) @@ -268,3 +269,49 @@ let u0old .= sjm_prob.prob.u0.jump_u end end + +# accuracy test based on +# https://github.com/SciML/JumpProcesses.jl/issues/320 +# note that testing that precisely is not trivial +let + rng = StableRNG(12345) + b = 2.0 + d = 1.0 + n0 = 1 + tspan = (0.0, 4.0) + Nsims = 10000 + n(t) = n0 * exp((b - d)*t) + u0 = [n0] + p = [b,d] + + function ode_fxn(du, u, p, t) + du .= 0 + nothing + end + + b_rate(u, p, t) = (u[1] * p[1]) + function birth!(integrator) + integrator.u[1] += 1 + nothing + end + b_jump = VariableRateJump(b_rate, birth!) + + d_rate(u, p, t) = (u[1] * p[2]) + function death!(integrator) + integrator.u[1] -= 1 + nothing + end + d_jump = VariableRateJump(d_rate, death!) + + ode_prob = ODEProblem(ode_fxn, u0, tspan, p) + sjm_prob = JumpProblem(ode_prob, b_jump, d_jump) + dt = .1 + tsave = range(tspan[1], tspan[2]; step = dt) + umean = zeros(length(tsave)) + for i in 1:Nsims + sol = solve(sjm_prob, Tsit5(); saveat = dt) + umean .+= Array(sol(tsave; idxs = 1)) + end + umean ./= Nsims + @test all(abs.(umean .- n.(tsave)) .< .05 * n.(tsave)) +end \ No newline at end of file From 3f74888adf0519e54a31ad41999535ba1188b9e4 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Sun, 1 Sep 2024 21:21:25 -0400 Subject: [PATCH 6/8] fix broadcast issue --- src/solve.jl | 7 ++++--- test/variable_rate.jl | 8 +++++--- 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/src/solve.jl b/src/solve.jl index c64046c7..6ba2fe23 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -58,8 +58,8 @@ function resetted_jump_problem(_jump_prob, seed) if !isempty(jump_prob.variable_jumps) @assert jump_prob.prob.u0 isa ExtendedJumpArray - @. jump_prob.prob.u0.jump_u = -randexp(_jump_prob.rng, - eltype(_jump_prob.prob.tspan)) + ttype = eltype(_jump_prob.prob.tspan) + @. jump_prob.prob.u0.jump_u = -randexp(_jump_prob.rng, ttype) end jump_prob end @@ -71,6 +71,7 @@ function reset_jump_problem!(jump_prob, seed) if !isempty(jump_prob.variable_jumps) @assert jump_prob.prob.u0 isa ExtendedJumpArray - @. jump_prob.prob.u0.jump_u = -randexp(jump_prob.rng, eltype(jump_prob.prob.tspan)) + ttype = eltype(jump_prob.prob.tspan) + @. jump_prob.prob.u0.jump_u = -randexp(jump_prob.rng, ttype) end end diff --git a/test/variable_rate.jl b/test/variable_rate.jl index 64b3873f..d7bf2d9c 100644 --- a/test/variable_rate.jl +++ b/test/variable_rate.jl @@ -261,9 +261,10 @@ let ode_prob = ODEProblem(ode_fxn, u0, tspan, p) sjm_prob = JumpProblem(ode_prob, b_jump, d_jump; rng) + @test allunique(sjm_prob.prob.u0.jump_u) u0old = copy(sjm_prob.prob.u0.jump_u) for i in 1:Nsims - sol = solve(sjm_prob, Tsit5(); saveat = tspan[2]) + sol = solve(sjm_prob, Tsit5(); saveat = tspan[2]) @test allunique(sjm_prob.prob.u0.jump_u) @test all(u0old != sjm_prob.prob.u0.jump_u) u0old .= sjm_prob.prob.u0.jump_u @@ -272,7 +273,8 @@ end # accuracy test based on # https://github.com/SciML/JumpProcesses.jl/issues/320 -# note that testing that precisely is not trivial +# note that even with the seeded StableRNG this test is not +# deterministic for some reason. let rng = StableRNG(12345) b = 2.0 @@ -304,7 +306,7 @@ let d_jump = VariableRateJump(d_rate, death!) ode_prob = ODEProblem(ode_fxn, u0, tspan, p) - sjm_prob = JumpProblem(ode_prob, b_jump, d_jump) + sjm_prob = JumpProblem(ode_prob, b_jump, d_jump; rng) dt = .1 tsave = range(tspan[1], tspan[2]; step = dt) umean = zeros(length(tsave)) From 7be0f0177d9d53d32fe459c45540a9c80d769be9 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Sun, 1 Sep 2024 21:22:16 -0400 Subject: [PATCH 7/8] format --- test/variable_rate.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/test/variable_rate.jl b/test/variable_rate.jl index d7bf2d9c..17800adb 100644 --- a/test/variable_rate.jl +++ b/test/variable_rate.jl @@ -264,7 +264,7 @@ let @test allunique(sjm_prob.prob.u0.jump_u) u0old = copy(sjm_prob.prob.u0.jump_u) for i in 1:Nsims - sol = solve(sjm_prob, Tsit5(); saveat = tspan[2]) + sol = solve(sjm_prob, Tsit5(); saveat = tspan[2]) @test allunique(sjm_prob.prob.u0.jump_u) @test all(u0old != sjm_prob.prob.u0.jump_u) u0old .= sjm_prob.prob.u0.jump_u @@ -279,12 +279,12 @@ let rng = StableRNG(12345) b = 2.0 d = 1.0 - n0 = 1 + n0 = 1 tspan = (0.0, 4.0) Nsims = 10000 - n(t) = n0 * exp((b - d)*t) + n(t) = n0 * exp((b - d) * t) u0 = [n0] - p = [b,d] + p = [b, d] function ode_fxn(du, u, p, t) du .= 0 @@ -293,7 +293,7 @@ let b_rate(u, p, t) = (u[1] * p[1]) function birth!(integrator) - integrator.u[1] += 1 + integrator.u[1] += 1 nothing end b_jump = VariableRateJump(b_rate, birth!) @@ -307,7 +307,7 @@ let ode_prob = ODEProblem(ode_fxn, u0, tspan, p) sjm_prob = JumpProblem(ode_prob, b_jump, d_jump; rng) - dt = .1 + dt = 0.1 tsave = range(tspan[1], tspan[2]; step = dt) umean = zeros(length(tsave)) for i in 1:Nsims @@ -315,5 +315,5 @@ let umean .+= Array(sol(tsave; idxs = 1)) end umean ./= Nsims - @test all(abs.(umean .- n.(tsave)) .< .05 * n.(tsave)) -end \ No newline at end of file + @test all(abs.(umean .- n.(tsave)) .< 0.05 * n.(tsave)) +end From e6ff52d3d32ff8bb0eab986f7223c0c2e0e16f17 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Sun, 1 Sep 2024 21:40:03 -0400 Subject: [PATCH 8/8] use inplace randexp --- src/solve.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/solve.jl b/src/solve.jl index 6ba2fe23..ee51ce98 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -58,8 +58,8 @@ function resetted_jump_problem(_jump_prob, seed) if !isempty(jump_prob.variable_jumps) @assert jump_prob.prob.u0 isa ExtendedJumpArray - ttype = eltype(_jump_prob.prob.tspan) - @. jump_prob.prob.u0.jump_u = -randexp(_jump_prob.rng, ttype) + randexp!(_jump_prob.rng, jump_prob.prob.u0.jump_u) + jump_prob.prob.u0.jump_u .*= -1 end jump_prob end @@ -71,7 +71,7 @@ function reset_jump_problem!(jump_prob, seed) if !isempty(jump_prob.variable_jumps) @assert jump_prob.prob.u0 isa ExtendedJumpArray - ttype = eltype(jump_prob.prob.tspan) - @. jump_prob.prob.u0.jump_u = -randexp(jump_prob.rng, ttype) + randexp!(jump_prob.rng, jump_prob.prob.u0.jump_u) + jump_prob.prob.u0.jump_u .*= -1 end end