From 11df15986a7f0e5514b5429b95910ba37c468a5e Mon Sep 17 00:00:00 2001 From: Vasily Ilin Date: Fri, 22 Sep 2023 19:20:54 -0400 Subject: [PATCH 1/4] Choose an SSA if no SSA is passed in `JumpProblem`. --- src/aggregators/aggregators.jl | 26 ++++++++++++++++++++++++++ src/problem.jl | 7 +++++-- test/geneexpr_test.jl | 2 +- test/linearreaction_test.jl | 2 +- test/spatial/ABC.jl | 3 +++ 5 files changed, 36 insertions(+), 4 deletions(-) diff --git a/src/aggregators/aggregators.jl b/src/aggregators/aggregators.jl index c1553d03..a2cdacc4 100644 --- a/src/aggregators/aggregators.jl +++ b/src/aggregators/aggregators.jl @@ -184,6 +184,32 @@ needs_vartojumps_map(aggregator::RSSACR) = true supports_variablerates(aggregator::AbstractAggregatorAlgorithm) = false supports_variablerates(aggregator::Coevolve) = true +# true if aggregator supports hops, e.g. diffusion is_spatial(aggregator::AbstractAggregatorAlgorithm) = false is_spatial(aggregator::NSM) = true is_spatial(aggregator::DirectCRDirect) = true + +# return the fastest aggregator out of the available ones +function select_aggregator(jumps::JumpSet; vartojumps_map=nothing, jumptovars_map=nothing, dep_graph=nothing, spatial_system=nothing, hopping_constants=nothing) + + # detect if a spatial SSA should be used + !isnothing(spatial_system) && !isnothing(hopping_constants) && return DirectCRDirect + + # if variable rate jumps are present, return the only SSA that supports them + num_vrjs(jumps) != 0 && return Coevolve + + # if the number of jumps is small, return the Direct + num_majumps(jumps) + num_crjs(jumps) < 10 && return Direct + + # if there are only massaction jumps, we can any build dependency graph, so return the fastest aggregator + num_crjs(jumps) == 0 && num_vrjs(jumps) == 0 && return RSSACR + + # in the presence of constant rate jumps, we rely on the given dependency graphs + if !isnothing(vartojumps_map) && !isnothing(jumptovars_map) + return RSSACR + elseif !isnothing(dep_graph) + return DirectCR + else + return Direct + end +end \ No newline at end of file diff --git a/src/problem.jl b/src/problem.jl index 192c7b8a..84b07582 100644 --- a/src/problem.jl +++ b/src/problem.jl @@ -171,9 +171,12 @@ function JumpProblem(prob, aggregator::AbstractAggregatorAlgorithm, jumps::Abstr kwargs...) JumpProblem(prob, aggregator, JumpSet(jumps...); kwargs...) end -function JumpProblem(prob, jumps::JumpSet; kwargs...) - JumpProblem(prob, NullAggregator(), jumps; kwargs...) +function JumpProblem(prob, jumps::JumpSet; vartojumps_map=nothing, jumptovars_map=nothing, dep_graph=nothing, spatial_system=nothing, hopping_constants=nothing, kwargs...) + aggregator = select_aggregator(jumps::JumpSet; vartojumps_map=vartojumps_map, jumptovars_map=jumptovars_map, dep_graph=dep_graph, spatial_system=spatial_system, hopping_constants=hopping_constants) + return JumpProblem(prob, aggregator(), jumps; vartojumps_map=vartojumps_map, jumptovars_map=jumptovars_map, dep_graph=dep_graph, spatial_system=spatial_system, hopping_constants=hopping_constants, kwargs...) end +# this makes it easier to test the aggregator selection +JumpProblem(prob, aggregator::NullAggregator, jumps::JumpSet; kwargs...) = JumpProblem(prob, jumps; kwargs...) make_kwarg(; kwargs...) = kwargs diff --git a/test/geneexpr_test.jl b/test/geneexpr_test.jl index 24e1e414..d4e0cd9b 100644 --- a/test/geneexpr_test.jl +++ b/test/geneexpr_test.jl @@ -12,7 +12,7 @@ dotestmean = true doprintmeans = false # SSAs to test -SSAalgs = (RDirect(), RSSACR(), Direct(), DirectFW(), FRM(), FRMFW(), SortingDirect(), +SSAalgs = (JumpProcesses.NullAggregator(), RDirect(), RSSACR(), Direct(), DirectFW(), FRM(), FRMFW(), SortingDirect(), NRM(), RSSA(), DirectCR(), Coevolve()) # numerical parameters diff --git a/test/linearreaction_test.jl b/test/linearreaction_test.jl index d169b571..e0e2112c 100644 --- a/test/linearreaction_test.jl +++ b/test/linearreaction_test.jl @@ -16,7 +16,7 @@ tf = 0.1 baserate = 0.1 A0 = 100 exactmean = (t, ratevec) -> A0 * exp(-sum(ratevec) * t) -SSAalgs = [RSSACR(), Direct(), RSSA()] +SSAalgs = [RSSACR(), Direct(), RSSA(), JumpProcesses.NullAggregator()] spec_to_dep_jumps = [collect(1:Nrxs), []] jump_to_dep_specs = [[1, 2] for i in 1:Nrxs] diff --git a/test/spatial/ABC.jl b/test/spatial/ABC.jl index c44358e8..50934d1f 100644 --- a/test/spatial/ABC.jl +++ b/test/spatial/ABC.jl @@ -56,6 +56,9 @@ jump_problems = JumpProblem[JumpProblem(prob, NSM(), majumps, push!(jump_problems, JumpProblem(prob, DirectCRDirect(), majumps, hopping_constants = hopping_constants, spatial_system = grids[1], save_positions = (false, false), rng = rng)) +push!(jump_problems, +JumpProblem(prob, majumps, hopping_constants = hopping_constants, + spatial_system = grids[1], save_positions = (false, false), rng = rng)) # setup flattenned jump prob push!(jump_problems, JumpProblem(prob, NRM(), majumps, hopping_constants = hopping_constants, From ed966e9c9b77649868d1a523dbab3a269634a7ab Mon Sep 17 00:00:00 2001 From: Vasily Ilin Date: Fri, 12 Jan 2024 22:41:01 -0800 Subject: [PATCH 2/4] Address comments. --- src/aggregators/aggregators.jl | 20 +++++++++++++------- test/spatial/ABC.jl | 2 +- 2 files changed, 14 insertions(+), 8 deletions(-) diff --git a/src/aggregators/aggregators.jl b/src/aggregators/aggregators.jl index a2cdacc4..280a565f 100644 --- a/src/aggregators/aggregators.jl +++ b/src/aggregators/aggregators.jl @@ -195,18 +195,24 @@ function select_aggregator(jumps::JumpSet; vartojumps_map=nothing, jumptovars_ma # detect if a spatial SSA should be used !isnothing(spatial_system) && !isnothing(hopping_constants) && return DirectCRDirect - # if variable rate jumps are present, return the only SSA that supports them - num_vrjs(jumps) != 0 && return Coevolve + # if variable rate jumps are present, return one of the two SSAs that support them + if num_vrjs(jumps) != 0 + any(isbounded, vrjs) && return Coevolve + return Direct + end # if the number of jumps is small, return the Direct - num_majumps(jumps) + num_crjs(jumps) < 10 && return Direct + num_jumps(jumps) < 10 && return Direct - # if there are only massaction jumps, we can any build dependency graph, so return the fastest aggregator - num_crjs(jumps) == 0 && num_vrjs(jumps) == 0 && return RSSACR + # if there are only massaction jumps, we can any build dependency graph + can_build_dependency_graphs = num_crjs(jumps) == 0 && num_vrjs(jumps) == 0 + have_dependency_graphs = !isnothing(vartojumps_map) && !isnothing(jumptovars_map) - # in the presence of constant rate jumps, we rely on the given dependency graphs - if !isnothing(vartojumps_map) && !isnothing(jumptovars_map) + # if we have the species-jumps dependency graphs or can build them, use one of the Rejection-based methods + if can_build_dependency_graphs || have_dependency_graphs + num_jumps(jumps) < 100 && return RSSA return RSSACR + # if we have the jumps-jumps dependency graph, use the Composition-Rejection Direct method elseif !isnothing(dep_graph) return DirectCR else diff --git a/test/spatial/ABC.jl b/test/spatial/ABC.jl index 50934d1f..2b706d84 100644 --- a/test/spatial/ABC.jl +++ b/test/spatial/ABC.jl @@ -57,7 +57,7 @@ push!(jump_problems, JumpProblem(prob, DirectCRDirect(), majumps, hopping_constants = hopping_constants, spatial_system = grids[1], save_positions = (false, false), rng = rng)) push!(jump_problems, -JumpProblem(prob, majumps, hopping_constants = hopping_constants, + JumpProblem(prob, majumps, hopping_constants = hopping_constants, spatial_system = grids[1], save_positions = (false, false), rng = rng)) # setup flattenned jump prob push!(jump_problems, From 6762a72e7b3473117803dcd41b181e2eeec75b31 Mon Sep 17 00:00:00 2001 From: Vasily Ilin Date: Sat, 13 Jan 2024 11:33:14 -0800 Subject: [PATCH 3/4] Reword variables and add a no-aggregator test. --- src/aggregators/aggregators.jl | 6 +++--- test/geneexpr_test.jl | 9 +++++++++ 2 files changed, 12 insertions(+), 3 deletions(-) diff --git a/src/aggregators/aggregators.jl b/src/aggregators/aggregators.jl index 280a565f..36d71a65 100644 --- a/src/aggregators/aggregators.jl +++ b/src/aggregators/aggregators.jl @@ -204,12 +204,12 @@ function select_aggregator(jumps::JumpSet; vartojumps_map=nothing, jumptovars_ma # if the number of jumps is small, return the Direct num_jumps(jumps) < 10 && return Direct - # if there are only massaction jumps, we can any build dependency graph + # if there are only massaction jumps, we can any build the species-jumps dependency graphs can_build_dependency_graphs = num_crjs(jumps) == 0 && num_vrjs(jumps) == 0 - have_dependency_graphs = !isnothing(vartojumps_map) && !isnothing(jumptovars_map) + have_species_to_jumps_dependency_graphs = !isnothing(vartojumps_map) && !isnothing(jumptovars_map) # if we have the species-jumps dependency graphs or can build them, use one of the Rejection-based methods - if can_build_dependency_graphs || have_dependency_graphs + if can_build_dependency_graphs || have_species_to_jumps_dependency_graphs num_jumps(jumps) < 100 && return RSSA return RSSACR # if we have the jumps-jumps dependency graph, use the Composition-Rejection Direct method diff --git a/test/geneexpr_test.jl b/test/geneexpr_test.jl index d4e0cd9b..56240348 100644 --- a/test/geneexpr_test.jl +++ b/test/geneexpr_test.jl @@ -116,3 +116,12 @@ end # end # println() # end + +# no-aggregator tests +jump_prob = JumpProblem(prob, majumps, save_positions = (false, false), + vartojumps_map = spec_to_dep_jumps, + jumptovars_map = jump_to_dep_specs, rng = rng) +@test abs(runSSAs(jump_prob) - expected_avg) < reltol * expected_avg + +jump_prob = JumpProblem(prob, majumps, save_positions = (false, false), rng = rng) +@test abs(runSSAs(jump_prob) - expected_avg) < reltol * expected_avg From 19ff59df2b92f93f0ccae819c121544ee0f50a82 Mon Sep 17 00:00:00 2001 From: Vasily Ilin Date: Sat, 13 Jan 2024 11:34:27 -0800 Subject: [PATCH 4/4] Fix typo in comment. --- src/aggregators/aggregators.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/aggregators/aggregators.jl b/src/aggregators/aggregators.jl index 36d71a65..b577df61 100644 --- a/src/aggregators/aggregators.jl +++ b/src/aggregators/aggregators.jl @@ -204,7 +204,7 @@ function select_aggregator(jumps::JumpSet; vartojumps_map=nothing, jumptovars_ma # if the number of jumps is small, return the Direct num_jumps(jumps) < 10 && return Direct - # if there are only massaction jumps, we can any build the species-jumps dependency graphs + # if there are only massaction jumps, we can build the species-jumps dependency graphs can_build_dependency_graphs = num_crjs(jumps) == 0 && num_vrjs(jumps) == 0 have_species_to_jumps_dependency_graphs = !isnothing(vartojumps_map) && !isnothing(jumptovars_map)