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,