From c7e3866945bc21c09f45517a510f56642f0a0d2e Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Tue, 6 Aug 2024 13:36:56 -0400 Subject: [PATCH 1/6] update JumpSystem for auto-alg support --- src/systems/jumps/jumpsystem.jl | 3 ++- test/jumpsystem.jl | 10 +++++++--- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/src/systems/jumps/jumpsystem.jl b/src/systems/jumps/jumpsystem.jl index da75b7dfd6..29f2906c6e 100644 --- a/src/systems/jumps/jumpsystem.jl +++ b/src/systems/jumps/jumpsystem.jl @@ -426,7 +426,8 @@ jprob = JumpProblem(complete(js), dprob, Direct()) sol = solve(jprob, SSAStepper()) ``` """ -function JumpProcesses.JumpProblem(js::JumpSystem, prob, aggregator; callback = nothing, +function JumpProcesses.JumpProblem(js::JumpSystem, prob, + aggregator = JumpProcesses.NullAggregator(); callback = nothing, eval_expression = false, eval_module = @__MODULE__, kwargs...) if !iscomplete(js) error("A completed `JumpSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `JumpProblem`") diff --git a/test/jumpsystem.jl b/test/jumpsystem.jl index 11c9fc1cd9..2af7adfb86 100644 --- a/test/jumpsystem.jl +++ b/test/jumpsystem.jl @@ -69,16 +69,21 @@ parammap = [β => 0.1 / 1000, γ => 0.01] dprob = DiscreteProblem(js2, u₀map, tspan, parammap) jprob = JumpProblem(js2, dprob, Direct(), save_positions = (false, false), rng = rng) Nsims = 30000 -function getmean(jprob, Nsims) +function getmean(jprob, Nsims; use_stepper = true) m = 0.0 for i in 1:Nsims - sol = solve(jprob, SSAStepper()) + sol = use_stepper ? solve(jprob, SSAStepper()) : solve(jprob) m += sol[end, end] end m / Nsims end m = getmean(jprob, Nsims) +# test auto-alg selection works +jprobb = JumpProblem(js2, dprob; save_positions = (false, false), rng) +mb = getmean(jprobb, Nsims; use_stepper = false) +@test abs(m - mb) / m < 0.01 + @variables S2(t) obs = [S2 ~ 2 * S] @named js2b = JumpSystem([j₁, j₃], t, [S, I, R], [β, γ], observed = obs) @@ -89,7 +94,6 @@ sol = solve(jprob, SSAStepper(), saveat = tspan[2] / 10) @test all(2 .* sol[S] .== sol[S2]) # test save_positions is working - jprob = JumpProblem(js2, dprob, Direct(), save_positions = (false, false), rng = rng) sol = solve(jprob, SSAStepper(), saveat = 1.0) @test all((sol.t) .== collect(0.0:tspan[2])) From 9103a670f8468076cd5b2c966eddab5cee093b1c Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Tue, 6 Aug 2024 14:14:48 -0400 Subject: [PATCH 2/6] update Project for JumpProcesses --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 9b54b5a676..4368fc9e82 100644 --- a/Project.toml +++ b/Project.toml @@ -89,7 +89,7 @@ FunctionWrappersWrappers = "0.1" Graphs = "1.5.2" InteractiveUtils = "1" JuliaFormatter = "1.0.47" -JumpProcesses = "9.1" +JumpProcesses = "9.13" LabelledArrays = "1.3" Latexify = "0.11, 0.12, 0.13, 0.14, 0.15, 0.16" Libdl = "1" From 7c05c99e0ccfd1292a3279bc25d51f410a93be05 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Tue, 6 Aug 2024 15:44:17 -0400 Subject: [PATCH 3/6] generate dep graphs --- src/systems/jumps/jumpsystem.jl | 3 ++- test/jumpsystem.jl | 20 ++++++++++++++++++++ 2 files changed, 22 insertions(+), 1 deletion(-) diff --git a/src/systems/jumps/jumpsystem.jl b/src/systems/jumps/jumpsystem.jl index 29f2906c6e..1a925922d2 100644 --- a/src/systems/jumps/jumpsystem.jl +++ b/src/systems/jumps/jumpsystem.jl @@ -449,7 +449,8 @@ function JumpProcesses.JumpProblem(js::JumpSystem, prob, error("Use continuous problems such as an ODEProblem or a SDEProblem with VariableRateJumps") jset = JumpSet(Tuple(vrjs), Tuple(crjs), nothing, majs) - if needs_vartojumps_map(aggregator) || needs_depgraph(aggregator) + if needs_vartojumps_map(aggregator) || needs_depgraph(aggregator) || + (aggregator isa JumpProcesses.NullAggregator) jdeps = asgraph(js) vdeps = variable_dependencies(js) vtoj = jdeps.badjlist diff --git a/test/jumpsystem.jl b/test/jumpsystem.jl index 2af7adfb86..a730c87bfa 100644 --- a/test/jumpsystem.jl +++ b/test/jumpsystem.jl @@ -274,3 +274,23 @@ affect = [X ~ X - 1] j1 = ConstantRateJump(k, [X ~ X - 1]) @test_nowarn @mtkbuild js1 = JumpSystem([j1], t, [X], [k]) + + +# test correct autosolver is selected +let + @parameters k + @variables X(t) + rate = k + affect = [X ~ X - 1] + j1 = ConstantRateJump(k, [X ~ X - 1]) + + Nv = [1, JumpProcesses.USE_DIRECT_THRESHOLD + 1, JumpProcesses.USE_RSSA_THRESHOLD + 1] + algtypes = [Direct, RSSA, RSSACR] + for (N, algtype) in zip(Nv, algtypes) + @named jsys = JumpSystem([deepcopy(j1) for _ in 1:N], t, [X], [k]) + jsys = complete(jsys) + dprob = DiscreteProblem(jsys, [X => 10], (0.0, 10.0), [k => 1]) + jprob = JumpProblem(jsys, dprob) + @test jprob.aggregator isa algtype + end +end \ No newline at end of file From 55a1f940ad56388dd22bd0d49aca123ba0d59d1f Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Tue, 6 Aug 2024 15:44:41 -0400 Subject: [PATCH 4/6] comment --- test/jumpsystem.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/jumpsystem.jl b/test/jumpsystem.jl index a730c87bfa..21e660d16e 100644 --- a/test/jumpsystem.jl +++ b/test/jumpsystem.jl @@ -276,7 +276,7 @@ j1 = ConstantRateJump(k, [X ~ X - 1]) @test_nowarn @mtkbuild js1 = JumpSystem([j1], t, [X], [k]) -# test correct autosolver is selected +# test correct autosolver is selected, which implies appropriate dep graphs are available let @parameters k @variables X(t) From 10489f508fba5a59bfc26dc17a9edd8603334344 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Tue, 6 Aug 2024 15:45:34 -0400 Subject: [PATCH 5/6] format --- src/systems/jumps/jumpsystem.jl | 2 +- test/jumpsystem.jl | 3 +-- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/src/systems/jumps/jumpsystem.jl b/src/systems/jumps/jumpsystem.jl index 1a925922d2..0ad1a009eb 100644 --- a/src/systems/jumps/jumpsystem.jl +++ b/src/systems/jumps/jumpsystem.jl @@ -450,7 +450,7 @@ function JumpProcesses.JumpProblem(js::JumpSystem, prob, jset = JumpSet(Tuple(vrjs), Tuple(crjs), nothing, majs) if needs_vartojumps_map(aggregator) || needs_depgraph(aggregator) || - (aggregator isa JumpProcesses.NullAggregator) + (aggregator isa JumpProcesses.NullAggregator) jdeps = asgraph(js) vdeps = variable_dependencies(js) vtoj = jdeps.badjlist diff --git a/test/jumpsystem.jl b/test/jumpsystem.jl index 21e660d16e..d14fa8b545 100644 --- a/test/jumpsystem.jl +++ b/test/jumpsystem.jl @@ -275,7 +275,6 @@ affect = [X ~ X - 1] j1 = ConstantRateJump(k, [X ~ X - 1]) @test_nowarn @mtkbuild js1 = JumpSystem([j1], t, [X], [k]) - # test correct autosolver is selected, which implies appropriate dep graphs are available let @parameters k @@ -293,4 +292,4 @@ let jprob = JumpProblem(jsys, dprob) @test jprob.aggregator isa algtype end -end \ No newline at end of file +end From 4ef25cbbde3dc11a2a090388bee5a468450a7d46 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Tue, 6 Aug 2024 15:52:04 -0400 Subject: [PATCH 6/6] bump JumpProcesses --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 4368fc9e82..5245347838 100644 --- a/Project.toml +++ b/Project.toml @@ -89,7 +89,7 @@ FunctionWrappersWrappers = "0.1" Graphs = "1.5.2" InteractiveUtils = "1" JuliaFormatter = "1.0.47" -JumpProcesses = "9.13" +JumpProcesses = "9.13.1" LabelledArrays = "1.3" Latexify = "0.11, 0.12, 0.13, 0.14, 0.15, 0.16" Libdl = "1"