From 11df15986a7f0e5514b5429b95910ba37c468a5e Mon Sep 17 00:00:00 2001
From: Vasily Ilin <vasilin97@gmail.com>
Date: Fri, 22 Sep 2023 19:20:54 -0400
Subject: [PATCH 01/21] 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 <vasilin97@gmail.com>
Date: Fri, 12 Jan 2024 22:41:01 -0800
Subject: [PATCH 02/21] 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 <vasilin97@gmail.com>
Date: Sat, 13 Jan 2024 11:33:14 -0800
Subject: [PATCH 03/21] 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 <vasilin97@gmail.com>
Date: Sat, 13 Jan 2024 11:34:27 -0800
Subject: [PATCH 04/21] 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)
 

From 8d4a5788b2176f1e9f094ea6143fc9723aa37f07 Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Tue, 6 Aug 2024 09:44:40 -0400
Subject: [PATCH 05/21] format and tweak

---
 src/aggregators/aggregators.jl | 23 +++++++++++++----------
 src/problem.jl                 | 14 ++++++++++----
 test/geneexpr_test.jl          |  9 ++++-----
 3 files changed, 27 insertions(+), 19 deletions(-)

diff --git a/src/aggregators/aggregators.jl b/src/aggregators/aggregators.jl
index 0a62dc40..8f3ab9b2 100644
--- a/src/aggregators/aggregators.jl
+++ b/src/aggregators/aggregators.jl
@@ -194,8 +194,10 @@ 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)
-    
+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
 
@@ -204,22 +206,23 @@ function select_aggregator(jumps::JumpSet; vartojumps_map=nothing, jumptovars_ma
         any(isbounded, vrjs) && return Coevolve
         return Direct
     end
-    
+
     # if the number of jumps is small, return the Direct
-    num_jumps(jumps) < 10 && return Direct
+    num_jumps(jumps) < 20 && return Direct
 
     # 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)
+    can_build_dgs = num_crjs(jumps) == 0 && num_vrjs(jumps) == 0
+    have_species_to_jumps_dgs = !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_species_to_jumps_dependency_graphs
-        num_jumps(jumps) < 100 && return RSSA
+    if can_build_dgs || have_species_to_jumps_dgs
+        (num_jumps(jumps) < 100) && return RSSA
         return RSSACR
-    # if we have the jumps-jumps dependency graph, use the Composition-Rejection Direct method
+        # if we have the jumps-jumps dependency graph, use the Composition-Rejection Direct method
     elseif !isnothing(dep_graph)
+        (num_jumps(jumps) < 200) && return SortingDirect
         return DirectCR
     else
         return Direct
     end
-end
\ No newline at end of file
+end
diff --git a/src/problem.jl b/src/problem.jl
index 9d83e7ef..6af1f1e3 100644
--- a/src/problem.jl
+++ b/src/problem.jl
@@ -176,12 +176,18 @@ function JumpProblem(prob, aggregator::AbstractAggregatorAlgorithm, jumps::Abstr
         kwargs...)
     JumpProblem(prob, aggregator, JumpSet(jumps...); kwargs...)
 end
-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...)
+function JumpProblem(prob, jumps::JumpSet; vartojumps_map = nothing,
+        jumptovars_map = nothing, dep_graph = nothing,
+        spatial_system = nothing, hopping_constants = nothing, kwargs...)
+    ps = (; vartojumps_map, jumptovars_map, dep_graph, spatial_system, hopping_constants)
+    aggtype = select_aggregator(jumps; ps...)
+    return JumpProblem(prob, aggtype(), jumps; ps..., kwargs...)
 end
+
 # this makes it easier to test the aggregator selection
-JumpProblem(prob, aggregator::NullAggregator, jumps::JumpSet; kwargs...) = JumpProblem(prob, jumps; kwargs...)
+function JumpProblem(prob, aggregator::NullAggregator, jumps::JumpSet; kwargs...)
+    JumpProblem(prob, jumps; kwargs...)
+end
 
 make_kwarg(; kwargs...) = kwargs
 
diff --git a/test/geneexpr_test.jl b/test/geneexpr_test.jl
index 4e50bd16..8e77fd11 100644
--- a/test/geneexpr_test.jl
+++ b/test/geneexpr_test.jl
@@ -12,8 +12,8 @@ dotestmean = true
 doprintmeans = false
 
 # SSAs to test
-SSAalgs = (JumpProcesses.NullAggregator(), RDirect(), RSSACR(), Direct(), DirectFW(), FRM(), FRMFW(), SortingDirect(),
-           NRM(), RSSA(), DirectCR(), Coevolve())
+SSAalgs = (JumpProcesses.NullAggregator(), RDirect(), RSSACR(), Direct(),
+    DirectFW(), FRM(), FRMFW(), SortingDirect(), NRM(), RSSA(), DirectCR(), Coevolve())
 
 # numerical parameters
 Nsims = 8000
@@ -118,9 +118,8 @@ end
 # 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)
+jump_prob = JumpProblem(prob, majumps; save_positions = (false, false),
+    vartojumps_map = spec_to_dep_jumps, jumptovars_map = jump_to_dep_specs, rng)
 @test abs(runSSAs(jump_prob) - expected_avg) < reltol * expected_avg
 
 jump_prob = JumpProblem(prob, majumps, save_positions = (false, false), rng = rng)

From e3256b01d6d4765371210beabd7e0336848d7c96 Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Tue, 6 Aug 2024 10:00:34 -0400
Subject: [PATCH 06/21] add default solver choice too

---
 src/aggregators/aggregators.jl | 9 ++++-----
 src/solve.jl                   | 7 +++++++
 2 files changed, 11 insertions(+), 5 deletions(-)

diff --git a/src/aggregators/aggregators.jl b/src/aggregators/aggregators.jl
index 8f3ab9b2..be167cb6 100644
--- a/src/aggregators/aggregators.jl
+++ b/src/aggregators/aggregators.jl
@@ -202,8 +202,8 @@ function select_aggregator(jumps::JumpSet; vartojumps_map = nothing,
     !isnothing(spatial_system) && !isnothing(hopping_constants) && return DirectCRDirect
 
     # 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
+    if num_vrjs(jumps) > 0
+        (num_bndvrjs(jumps) > 0) && return Coevolve
         return Direct
     end
 
@@ -214,12 +214,11 @@ function select_aggregator(jumps::JumpSet; vartojumps_map = nothing,
     can_build_dgs = num_crjs(jumps) == 0 && num_vrjs(jumps) == 0
     have_species_to_jumps_dgs = !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 we have the species-jumps dgs or can build them, use a Rejection-based methods
     if can_build_dgs || have_species_to_jumps_dgs
         (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)
+    elseif !isnothing(dep_graph)   # if only have a normal dg
         (num_jumps(jumps) < 200) && return SortingDirect
         return DirectCR
     else
diff --git a/src/solve.jl b/src/solve.jl
index 32908152..8add6ea8 100644
--- a/src/solve.jl
+++ b/src/solve.jl
@@ -7,6 +7,13 @@ function DiffEqBase.__solve(jump_prob::DiffEqBase.AbstractJumpProblem{P},
     integrator.sol
 end
 
+# if passed a JumpProblem and no aggregator is selected use SSAStepper
+function DiffEqBase.__solve(jump_prob::DiffEqBase.AbstractJumpProblem{P},
+        timeseries = [], ts = [], ks = [], recompile::Type{Val{recompile_flag}} = Val{true};
+        kwargs...) where {P, recompile_flag}
+    DiffEqBase.__solve(jump_prob, SSAStepper(), timeseries, ts, ks, recompile)
+end
+
 function DiffEqBase.__init(_jump_prob::DiffEqBase.AbstractJumpProblem{P},
         alg::DiffEqBase.DEAlgorithm, timeseries = [], ts = [], ks = [],
         recompile::Type{Val{recompile_flag}} = Val{true};

From aedf61153eb17eaf92cc00fe41c55f0499a28c2c Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Tue, 6 Aug 2024 10:28:51 -0400
Subject: [PATCH 07/21] clean up tests

---
 test/bimolerx_test.jl       | 40 +++++++++++-----------------------
 test/degenerate_rx_cases.jl |  3 +--
 test/extinction_test.jl     |  2 +-
 test/geneexpr_test.jl       | 42 ++++++++++++------------------------
 test/linearreaction_test.jl | 43 ++++++++++++++++++-------------------
 test/reversible_binding.jl  |  2 +-
 6 files changed, 51 insertions(+), 81 deletions(-)

diff --git a/test/bimolerx_test.jl b/test/bimolerx_test.jl
index 330b6946..f9df0460 100644
--- a/test/bimolerx_test.jl
+++ b/test/bimolerx_test.jl
@@ -14,7 +14,7 @@ dotestmean = true
 doprintmeans = false
 
 # SSAs to test
-SSAalgs = JumpProcesses.JUMP_AGGREGATORS
+SSAalgs = (JumpProcesses.JUMP_AGGREGATORS..., JumpProcesses.NullAggregator)
 
 Nsims = 32000
 tf = 0.01
@@ -55,10 +55,10 @@ jump_to_dep_specs = [[1, 2], [1, 2], [1, 2, 3], [1, 2, 3], [1, 3]]
 majumps = MassActionJump(rates, reactstoch, netstoch)
 
 # average number of proteins in a simulation
-function runSSAs(jump_prob)
+function runSSAs(jump_prob; use_stepper = true)
     Psamp = zeros(Int, Nsims)
     for i in 1:Nsims
-        sol = solve(jump_prob, SSAStepper())
+        sol = use_stepper ? solve(jump_prob, SSAStepper()) : solve(jump_prob)
         Psamp[i] = sol[1, end]
     end
     mean(Psamp)
@@ -81,37 +81,23 @@ end
 
 # test the means
 if dotestmean
-    means = zeros(Float64, length(SSAalgs))
     for (i, alg) in enumerate(SSAalgs)
         local jump_prob = JumpProblem(prob, alg, majumps, save_positions = (false, false),
             vartojumps_map = spec_to_dep_jumps,
             jumptovars_map = jump_to_dep_specs, rng = rng)
-        means[i] = runSSAs(jump_prob)
-        relerr = abs(means[i] - expected_avg) / expected_avg
-        if doprintmeans
-            println("Mean from method: ", typeof(alg), " is = ", means[i], ", rel err = ",
-                relerr)
-        end
-
-        # if dobenchmark
-        #      @btime (runSSAs($jump_prob);)
-        # end
-
-        @test abs(means[i] - expected_avg) < reltol * expected_avg
+        means = runSSAs(jump_prob)
+        relerr = abs(means - expected_avg) / expected_avg
+        doprintmeans && println("Mean from method: ", typeof(alg), " is = ", means,
+            ", rel err = ", relerr)
+        @test abs(means - expected_avg) < reltol * expected_avg
+
+        # test not specifying SSAStepper
+        means = runSSAs(jump_prob; use_stepper = false)
+        relerr = abs(means - expected_avg) / expected_avg
+        @test abs(means - expected_avg) < reltol * expected_avg
     end
 end
 
-# benchmark performance
-# if dobenchmark
-#     # exact methods
-#     for alg in SSAalgs
-#         println("Solving with method: ", typeof(alg), ", using SSAStepper")
-#         jump_prob = JumpProblem(prob, alg, majumps, vartojumps_map=spec_to_dep_jumps, jumptovars_map=jump_to_dep_specs, rng=rng)
-#         @btime solve($jump_prob, SSAStepper())
-#     end
-#     println()
-# end
-
 # add a test for passing MassActionJumps individually (tests combining)
 if dotestmean
     majump_vec = Vector{MassActionJump{Float64, Vector{Pair{Int, Int}}}}()
diff --git a/test/degenerate_rx_cases.jl b/test/degenerate_rx_cases.jl
index a965ae98..a56b86ce 100644
--- a/test/degenerate_rx_cases.jl
+++ b/test/degenerate_rx_cases.jl
@@ -12,8 +12,7 @@ doprint = false
 #using Plots; plotlyjs()
 doplot = false
 
-methods = (RDirect(), RSSACR(), Direct(), DirectFW(), FRM(), FRMFW(), SortingDirect(),
-    NRM(), RSSA(), DirectCR(), Coevolve())
+methods = (JumpProcesses.JUMP_AGGREGATORS..., JumpProcesses.NullAggregator)
 
 # one reaction case, mass action jump, vector of data
 rate = [2.0]
diff --git a/test/extinction_test.jl b/test/extinction_test.jl
index 41e03498..04ae0bd7 100644
--- a/test/extinction_test.jl
+++ b/test/extinction_test.jl
@@ -17,7 +17,7 @@ dg = [[1]]
 majump = MassActionJump(rates, reactstoch, netstoch)
 u0 = [100000]
 dprob = DiscreteProblem(u0, (0.0, 1e5), rates)
-algs = JumpProcesses.JUMP_AGGREGATORS
+algs = (JumpProcesses.JUMP_AGGREGATORS..., JumpProcesses.NullAggregator)
 
 for n in 1:Nsims
     for ssa in algs
diff --git a/test/geneexpr_test.jl b/test/geneexpr_test.jl
index 8e77fd11..de1ea543 100644
--- a/test/geneexpr_test.jl
+++ b/test/geneexpr_test.jl
@@ -12,8 +12,7 @@ dotestmean = true
 doprintmeans = false
 
 # SSAs to test
-SSAalgs = (JumpProcesses.NullAggregator(), RDirect(), RSSACR(), Direct(),
-    DirectFW(), FRM(), FRMFW(), SortingDirect(), NRM(), RSSA(), DirectCR(), Coevolve())
+SSAalgs = (JumpProcesses.JUMP_AGGREGATORS..., JumpProcesses.NullAggregator)
 
 # numerical parameters
 Nsims = 8000
@@ -23,10 +22,10 @@ expected_avg = 5.926553750000000e+02
 reltol = 0.01
 
 # average number of proteins in a simulation
-function runSSAs(jump_prob)
+function runSSAs(jump_prob; use_stepper = true)
     Psamp = zeros(Int, Nsims)
     for i in 1:Nsims
-        sol = solve(jump_prob, SSAStepper())
+        sol = use_stepper ? solve(jump_prob, SSAStepper()) : solve(jump_prob)
         Psamp[i] = sol[3, end]
     end
     mean(Psamp)
@@ -86,41 +85,28 @@ end
 
 # test the means
 if dotestmean
-    means = zeros(Float64, length(SSAalgs))
     for (i, alg) in enumerate(SSAalgs)
         local jump_prob = JumpProblem(prob, alg, majumps, save_positions = (false, false),
             vartojumps_map = spec_to_dep_jumps,
             jumptovars_map = jump_to_dep_specs, rng = rng)
-        means[i] = runSSAs(jump_prob)
-        relerr = abs(means[i] - expected_avg) / expected_avg
-        if doprintmeans
-            println("Mean from method: ", typeof(alg), " is = ", means[i], ", rel err = ",
-                relerr)
-        end
-
-        # if dobenchmark
-        #     @btime (runSSAs($jump_prob);)
-        # end
-
-        @test abs(means[i] - expected_avg) < reltol * expected_avg
+        means = runSSAs(jump_prob)
+        relerr = abs(means - expected_avg) / expected_avg
+        doprintmeans && println("Mean from method: ", typeof(alg), " is = ", means,
+            ", rel err = ", relerr)
+        @test abs(means - expected_avg) < reltol * expected_avg
+
+        means = runSSAs(jump_prob; use_stepper = false)
+        relerr = abs(means - expected_avg) / expected_avg
+        @test abs(means - expected_avg) < reltol * expected_avg
     end
 end
 
-# benchmark performance
-# if dobenchmark
-#     # exact methods
-#     for alg in SSAalgs
-#         println("Solving with method: ", typeof(alg), ", using SSAStepper")
-#         jump_prob = JumpProblem(prob, alg, majumps, vartojumps_map=spec_to_dep_jumps, jumptovars_map=jump_to_dep_specs, rng=rng)
-#         @btime solve($jump_prob, SSAStepper())
-#     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)
 @test abs(runSSAs(jump_prob) - expected_avg) < reltol * expected_avg
+@test abs(runSSAs(jump_prob; use_stepper = false) - 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
+@test abs(runSSAs(jump_prob; use_stepper = false) - expected_avg) < reltol * expected_avg
diff --git a/test/linearreaction_test.jl b/test/linearreaction_test.jl
index 4a8c2cf1..5ae0cc67 100644
--- a/test/linearreaction_test.jl
+++ b/test/linearreaction_test.jl
@@ -16,11 +16,12 @@ tf = 0.1
 baserate = 0.1
 A0 = 100
 exactmean = (t, ratevec) -> A0 * exp(-sum(ratevec) * t)
-SSAalgs = [RSSACR(), Direct(), RSSA(), JumpProcesses.NullAggregator()]
+SSAalgs = (JumpProcesses.JUMP_AGGREGATORS..., JumpProcesses.NullAggregator)
+SSAs_to_exclude = (DirectFW(), FRM(), FRMFW())
 
-spec_to_dep_jumps = [collect(1:Nrxs), []]
-jump_to_dep_specs = [[1, 2] for i in 1:Nrxs]
-namedpars = (vartojumps_map = spec_to_dep_jumps, jumptovars_map = jump_to_dep_specs)
+vartojumps_map = [collect(1:Nrxs), []]
+jumptovars_map = [[1, 2] for i in 1:Nrxs]
+namedpars = (vartojumps_map, jumptovars_map)
 rates = ones(Float64, Nrxs) * baserate;
 cumsum!(rates, rates)
 exactmeanval = exactmean(tf, rates)
@@ -51,7 +52,7 @@ function A_to_B_tuple(N, method)
     jumps = ((jump for jump in jumpvec)...,)
     jset = JumpSet((), jumps, nothing, nothing)
     prob = DiscreteProblem([A0, 0], (0.0, tf))
-    jump_prob = JumpProblem(prob, method, jset; save_positions = (false, false), rng = rng,
+    jump_prob = JumpProblem(prob, method, jset; save_positions = (false, false), rng,
         namedpars...)
 
     jump_prob
@@ -73,7 +74,7 @@ function A_to_B_vec(N, method)
     # convert jumpvec to tuple to send to JumpProblem...
     jset = JumpSet((), jumps, nothing, nothing)
     prob = DiscreteProblem([A0, 0], (0.0, tf))
-    jump_prob = JumpProblem(prob, method, jset; save_positions = (false, false), rng = rng,
+    jump_prob = JumpProblem(prob, method, jset; save_positions = (false, false), rng,
         namedpars...)
 
     jump_prob
@@ -91,7 +92,7 @@ function A_to_B_ma(N, method)
     majumps = MassActionJump(rates, reactstoch, netstoch)
     jset = JumpSet((), (), nothing, majumps)
     prob = DiscreteProblem([A0, 0], (0.0, tf))
-    jump_prob = JumpProblem(prob, method, jset; save_positions = (false, false), rng = rng,
+    jump_prob = JumpProblem(prob, method, jset; save_positions = (false, false), rng,
         namedpars...)
 
     jump_prob
@@ -125,7 +126,7 @@ function A_to_B_hybrid(N, method)
     majumps = MassActionJump(rates[1:switchidx], reactstoch, netstoch)
     jset = JumpSet((), jumps, nothing, majumps)
     prob = DiscreteProblem([A0, 0], (0.0, tf))
-    jump_prob = JumpProblem(prob, method, jset; save_positions = (false, false), rng = rng,
+    jump_prob = JumpProblem(prob, method, jset; save_positions = (false, false), rng,
         namedpars...)
 
     jump_prob
@@ -160,7 +161,7 @@ function A_to_B_hybrid_nojset(N, method)
     jumps = (constjumps..., majumps)
     prob = DiscreteProblem([A0, 0], (0.0, tf))
     jump_prob = JumpProblem(prob, method, jumps...; save_positions = (false, false),
-        rng = rng, namedpars...)
+        rng, namedpars...)
 
     jump_prob
 end
@@ -189,7 +190,7 @@ function A_to_B_hybrid_vecs(N, method)
     end
     jset = JumpSet((), jumpvec, nothing, majumps)
     prob = DiscreteProblem([A0, 0], (0.0, tf))
-    jump_prob = JumpProblem(prob, method, jset; save_positions = (false, false), rng = rng,
+    jump_prob = JumpProblem(prob, method, jset; save_positions = (false, false), rng,
         namedpars...)
 
     jump_prob
@@ -219,7 +220,7 @@ function A_to_B_hybrid_vecs_scalars(N, method)
     end
     jset = JumpSet((), jumpvec, nothing, majumps)
     prob = DiscreteProblem([A0, 0], (0.0, tf))
-    jump_prob = JumpProblem(prob, method, jset; save_positions = (false, false), rng = rng,
+    jump_prob = JumpProblem(prob, method, jset; save_positions = (false, false), rng,
         namedpars...)
 
     jump_prob
@@ -251,7 +252,7 @@ function A_to_B_hybrid_tups_scalars(N, method)
     jumps = ((maj for maj in majumpsv)..., (jump for jump in jumpvec)...)
     prob = DiscreteProblem([A0, 0], (0.0, tf))
     jump_prob = JumpProblem(prob, method, jumps...; save_positions = (false, false),
-        rng = rng, namedpars...)
+        rng, namedpars...)
 
     jump_prob
 end
@@ -281,7 +282,7 @@ function A_to_B_hybrid_tups(N, method)
     jumps = ((jump for jump in jumpvec)...,)
     jset = JumpSet((), jumps, nothing, majumps)
     prob = DiscreteProblem([A0, 0], (0.0, tf))
-    jump_prob = JumpProblem(prob, method, jset; save_positions = (false, false), rng = rng,
+    jump_prob = JumpProblem(prob, method, jset; save_positions = (false, false), rng,
         namedpars...)
 
     jump_prob
@@ -293,6 +294,10 @@ jump_prob_gens = [A_to_B_tuple, A_to_B_vec, A_to_B_ma, A_to_B_hybrid, A_to_B_hyb
 #jump_prob_gens = [A_to_B_tuple, A_to_B_ma, A_to_B_hybrid, A_to_B_hybrid_vecs, A_to_B_hybrid_vecs_scalars,A_to_B_hybrid_tups_scalars]
 
 for method in SSAalgs
+
+    # skip methods that use regular dep graphs
+    (needs_depgraph(method) || (method in SSAs_to_exclude)) && continue
+
     for jump_prob_gen in jump_prob_gens
         local jump_prob = jump_prob_gen(Nrxs, method)
         meanval = runSSAs(jump_prob)
@@ -301,17 +306,15 @@ for method in SSAalgs
                 ", sample mean = ", meanval, ", actual mean = ", exactmeanval)
         end
         @test abs(meanval - exactmeanval) < 1.0
-
-        # if dobenchmark
-        #     @btime (runSSAs($jump_prob);)
-        # end
     end
 end
 
 # for dependency graph methods just test with mass action jumps
-SSAalgs = [RDirect(), NRM(), SortingDirect(), DirectCR(), Coevolve()]
 jump_prob_gens = [A_to_B_ma]
 for method in SSAalgs
+    # skip methods that don't use a normal dep graph
+    (!needs_depgraph(method) || method in SSAs_to_exclude) && continue
+
     for jump_prob_gen in jump_prob_gens
         local jump_prob = jump_prob_gen(Nrxs, method)
         meanval = runSSAs(jump_prob)
@@ -320,9 +323,5 @@ for method in SSAalgs
                 ", sample mean = ", meanval, ", actual mean = ", exactmeanval)
         end
         @test abs(meanval - exactmeanval) < 1.0
-
-        # if dobenchmark
-        #     @btime (runSSAs($jump_prob);)
-        # end
     end
 end
diff --git a/test/reversible_binding.jl b/test/reversible_binding.jl
index 82e7a3dc..d97cbd0c 100644
--- a/test/reversible_binding.jl
+++ b/test/reversible_binding.jl
@@ -45,7 +45,7 @@ function mastereqmean(u, rates)
 end
 mastereq_mean = mastereqmean(u0, rates)
 
-algs = JumpProcesses.JUMP_AGGREGATORS
+algs = (JumpProcesses.JUMP_AGGREGATORS..., JumpProcesses.NullAggregator)
 relative_tolerance = 0.01
 for alg in algs
     local jprob = JumpProblem(prob, alg, majumps, save_positions = (false, false),

From acde1e976382a9f9eba4b5d35d148c40a6af5255 Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Tue, 6 Aug 2024 10:35:26 -0400
Subject: [PATCH 08/21] fix test

---
 test/linearreaction_test.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/linearreaction_test.jl b/test/linearreaction_test.jl
index 5ae0cc67..2372c751 100644
--- a/test/linearreaction_test.jl
+++ b/test/linearreaction_test.jl
@@ -21,7 +21,7 @@ SSAs_to_exclude = (DirectFW(), FRM(), FRMFW())
 
 vartojumps_map = [collect(1:Nrxs), []]
 jumptovars_map = [[1, 2] for i in 1:Nrxs]
-namedpars = (vartojumps_map, jumptovars_map)
+namedpars = (; vartojumps_map, jumptovars_map)
 rates = ones(Float64, Nrxs) * baserate;
 cumsum!(rates, rates)
 exactmeanval = exactmean(tf, rates)

From 43f765cf6532b85a0bdb07fc043ad112141f3c4c Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Tue, 6 Aug 2024 10:36:56 -0400
Subject: [PATCH 09/21] drop version check

---
 test/runtests.jl | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index 3e4fa79a..01e06ecb 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -16,9 +16,7 @@ using JumpProcesses, DiffEqBase, SafeTestsets
     @time @safetestset "Mass Action Jump Tests; Gene Expr Model" begin include("geneexpr_test.jl") end
     @time @safetestset "Mass Action Jump Tests; Nonlinear Rx Model" begin include("bimolerx_test.jl") end
     @time @safetestset "Mass Action Jump Tests; Special Cases" begin include("degenerate_rx_cases.jl") end
-    @static if VERSION >= v"1.9.0"
-        @time @safetestset "Direct allocations test" begin include("allocations.jl") end
-    end
+    @time @safetestset "Direct allocations test" begin include("allocations.jl") end
     @time @safetestset "Bracketing Tests" begin include("bracketing.jl") end
     @time @safetestset "Composition-Rejection Table Tests" begin include("table_test.jl") end
     @time @safetestset "Extinction test" begin include("extinction_test.jl") end

From 10411f6ac9dab3ac9c701c299283d9763c35ce0d Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Tue, 6 Aug 2024 10:40:40 -0400
Subject: [PATCH 10/21] test fix again

---
 test/bimolerx_test.jl       | 2 +-
 test/degenerate_rx_cases.jl | 2 +-
 test/extinction_test.jl     | 2 +-
 test/geneexpr_test.jl       | 2 +-
 test/linearreaction_test.jl | 2 +-
 test/reversible_binding.jl  | 2 +-
 6 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/test/bimolerx_test.jl b/test/bimolerx_test.jl
index f9df0460..befeefab 100644
--- a/test/bimolerx_test.jl
+++ b/test/bimolerx_test.jl
@@ -14,7 +14,7 @@ dotestmean = true
 doprintmeans = false
 
 # SSAs to test
-SSAalgs = (JumpProcesses.JUMP_AGGREGATORS..., JumpProcesses.NullAggregator)
+SSAalgs = (JumpProcesses.JUMP_AGGREGATORS..., JumpProcesses.NullAggregator())
 
 Nsims = 32000
 tf = 0.01
diff --git a/test/degenerate_rx_cases.jl b/test/degenerate_rx_cases.jl
index a56b86ce..de6938c0 100644
--- a/test/degenerate_rx_cases.jl
+++ b/test/degenerate_rx_cases.jl
@@ -12,7 +12,7 @@ doprint = false
 #using Plots; plotlyjs()
 doplot = false
 
-methods = (JumpProcesses.JUMP_AGGREGATORS..., JumpProcesses.NullAggregator)
+methods = (JumpProcesses.JUMP_AGGREGATORS..., JumpProcesses.NullAggregator())
 
 # one reaction case, mass action jump, vector of data
 rate = [2.0]
diff --git a/test/extinction_test.jl b/test/extinction_test.jl
index 04ae0bd7..da6f7583 100644
--- a/test/extinction_test.jl
+++ b/test/extinction_test.jl
@@ -17,7 +17,7 @@ dg = [[1]]
 majump = MassActionJump(rates, reactstoch, netstoch)
 u0 = [100000]
 dprob = DiscreteProblem(u0, (0.0, 1e5), rates)
-algs = (JumpProcesses.JUMP_AGGREGATORS..., JumpProcesses.NullAggregator)
+algs = (JumpProcesses.JUMP_AGGREGATORS..., JumpProcesses.NullAggregator())
 
 for n in 1:Nsims
     for ssa in algs
diff --git a/test/geneexpr_test.jl b/test/geneexpr_test.jl
index de1ea543..97d97c22 100644
--- a/test/geneexpr_test.jl
+++ b/test/geneexpr_test.jl
@@ -12,7 +12,7 @@ dotestmean = true
 doprintmeans = false
 
 # SSAs to test
-SSAalgs = (JumpProcesses.JUMP_AGGREGATORS..., JumpProcesses.NullAggregator)
+SSAalgs = (JumpProcesses.JUMP_AGGREGATORS..., JumpProcesses.NullAggregator())
 
 # numerical parameters
 Nsims = 8000
diff --git a/test/linearreaction_test.jl b/test/linearreaction_test.jl
index 2372c751..657982a6 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 = (JumpProcesses.JUMP_AGGREGATORS..., JumpProcesses.NullAggregator)
+SSAalgs = (JumpProcesses.JUMP_AGGREGATORS..., JumpProcesses.NullAggregator())
 SSAs_to_exclude = (DirectFW(), FRM(), FRMFW())
 
 vartojumps_map = [collect(1:Nrxs), []]
diff --git a/test/reversible_binding.jl b/test/reversible_binding.jl
index d97cbd0c..f872d92a 100644
--- a/test/reversible_binding.jl
+++ b/test/reversible_binding.jl
@@ -45,7 +45,7 @@ function mastereqmean(u, rates)
 end
 mastereq_mean = mastereqmean(u0, rates)
 
-algs = (JumpProcesses.JUMP_AGGREGATORS..., JumpProcesses.NullAggregator)
+algs = (JumpProcesses.JUMP_AGGREGATORS..., JumpProcesses.NullAggregator())
 relative_tolerance = 0.01
 for alg in algs
     local jprob = JumpProblem(prob, alg, majumps, save_positions = (false, false),

From 5b3a3c386d565455dd2585fcc57a330d3d3c0cac Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Tue, 6 Aug 2024 10:48:39 -0400
Subject: [PATCH 11/21] fix solve dispatch

---
 src/SSA_stepper.jl | 4 +---
 src/solve.jl       | 6 ++----
 2 files changed, 3 insertions(+), 7 deletions(-)

diff --git a/src/SSA_stepper.jl b/src/SSA_stepper.jl
index f118d132..8c28a866 100644
--- a/src/SSA_stepper.jl
+++ b/src/SSA_stepper.jl
@@ -107,9 +107,7 @@ function DiffEqBase.u_modified!(integrator::SSAIntegrator, bool::Bool)
     integrator.u_modified = bool
 end
 
-function DiffEqBase.__solve(jump_prob::JumpProblem,
-        alg::SSAStepper;
-        kwargs...)
+function DiffEqBase.__solve(jump_prob::JumpProblem, alg::SSAStepper; kwargs...)
     integrator = init(jump_prob, alg; kwargs...)
     solve!(integrator)
     integrator.sol
diff --git a/src/solve.jl b/src/solve.jl
index 8add6ea8..1f68a513 100644
--- a/src/solve.jl
+++ b/src/solve.jl
@@ -8,10 +8,8 @@ function DiffEqBase.__solve(jump_prob::DiffEqBase.AbstractJumpProblem{P},
 end
 
 # if passed a JumpProblem and no aggregator is selected use SSAStepper
-function DiffEqBase.__solve(jump_prob::DiffEqBase.AbstractJumpProblem{P},
-        timeseries = [], ts = [], ks = [], recompile::Type{Val{recompile_flag}} = Val{true};
-        kwargs...) where {P, recompile_flag}
-    DiffEqBase.__solve(jump_prob, SSAStepper(), timeseries, ts, ks, recompile)
+function DiffEqBase.__solve(jump_prob::DiffEqBase.AbstractJumpProblem; kwargs...)
+    DiffEqBase.__solve(jump_prob, SSAStepper(); kwargs...)
 end
 
 function DiffEqBase.__init(_jump_prob::DiffEqBase.AbstractJumpProblem{P},

From 4ae4352abd436290ba58a2ab9b5d0abfb776ffea Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Tue, 6 Aug 2024 11:13:13 -0400
Subject: [PATCH 12/21] narrow __solve dispatch

---
 src/solve.jl | 10 ++++++++--
 1 file changed, 8 insertions(+), 2 deletions(-)

diff --git a/src/solve.jl b/src/solve.jl
index 1f68a513..bc9ed6c8 100644
--- a/src/solve.jl
+++ b/src/solve.jl
@@ -7,11 +7,17 @@ function DiffEqBase.__solve(jump_prob::DiffEqBase.AbstractJumpProblem{P},
     integrator.sol
 end
 
-# if passed a JumpProblem and no aggregator is selected use SSAStepper
-function DiffEqBase.__solve(jump_prob::DiffEqBase.AbstractJumpProblem; kwargs...)
+# if passed a JumpProblem over a DiscreteProblem, and no aggregator is selected use
+# SSAStepper
+function DiffEqBase.__solve(jump_prob::DiffEqBase.AbstractJumpProblem{P};
+        kwargs...) where {P <: DiscreteProblem}
     DiffEqBase.__solve(jump_prob, SSAStepper(); kwargs...)
 end
 
+function DiffEqBase.__solve(jump_prob::DiffEqBase.AbstractJumpProblem; kwargs...)
+    error("Auto-solver selection is currently only implemented for JumpProblems defined over DiscreteProblems. Please explicitly specify a solver algorithm in calling solve.")
+end
+
 function DiffEqBase.__init(_jump_prob::DiffEqBase.AbstractJumpProblem{P},
         alg::DiffEqBase.DEAlgorithm, timeseries = [], ts = [], ks = [],
         recompile::Type{Val{recompile_flag}} = Val{true};

From 9a4fa0ff9aea23886cace0f2a26b78ce1438a464 Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Tue, 6 Aug 2024 12:43:46 -0400
Subject: [PATCH 13/21] update HISTORY

---
 HISTORY.md | 46 +++++++++++++++++++++++++++++++++++++---------
 1 file changed, 37 insertions(+), 9 deletions(-)

diff --git a/HISTORY.md b/HISTORY.md
index 403fea38..37987c8f 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -2,17 +2,45 @@
 
 ## JumpProcesses unreleased (master branch)
 
+## 9.12
+- Added a default aggregator selection algorithm based on the number of passed
+  in jumps. i.e. the following now auto-selects an aggregator (`Direct` in this
+  case):
+  ```julia
+  using JumpProcesses
+  rate(u,p,t) = u[1]
+  affect(integrator) = (integrator.u[1] -= 1; nothing)
+  crj = ConstantRateJump(rate,affect)
+  dprob = DiscreteProblem([10],(0.0, 10.0))
+  jprob = JumpProblem(dprob, crj)
+  sol = solve(jprob, SSAStepper())
+  ```
+- For `JumpProblem`s over `DiscreteProblem`s that only have `MassActionJump`s
+  and `ConstantRateJump`s one no longer needs to specify `SSAStepper()` when
+  calling `solve`, i.e. the following now works for the previous example and is
+  equivalent to manually passing `SSAStepper()`:
+  ```julia
+  sol = solve(jprob)
+  ```
+- Plotting a solution generated with `save_positions = (false, false)` no uses
+  piecewise linear plots between any saved time points specified via `saveat`
+  instead (previously the plots appeared piecewise constant even though each
+  jump was not being shown). Note that solution objects still use piecewise
+  constant interpolation, see [the
+  docs](https://docs.sciml.ai/JumpProcesses/stable/tutorials/discrete_stochastic_example/#save_positions_docs)
+  for details.
+
 ## 9.7
 
-  - `Coevolve` was updated to support use with coupled ODEs/SDEs. See the updated
-    documentation for details, and note the comments there about one needing to ensure
-    rate bounds hold however the ODE/SDE stepper could modify dependent variables during a timestep.
+- `Coevolve` was updated to support use with coupled ODEs/SDEs. See the updated
+  documentation for details, and note the comments there about one needing to ensure
+  rate bounds hold however the ODE/SDE stepper could modify dependent variables during a timestep.
 
 ## 9.3
 
-  - Support for "bounded" `VariableRateJump`s that can be used with the `Coevolve`
-    aggregator for faster simulation of jump processes with time-dependent rates.
-    In particular, if all `VariableRateJump`s in a pure-jump system are bounded one
-    can use `Coevolve` with `SSAStepper` for better performance. See the
-    documentation, particularly the first and second tutorials, for details on
-    defining and using bounded `VariableRateJump`s.
+- Support for "bounded" `VariableRateJump`s that can be used with the `Coevolve`
+  aggregator for faster simulation of jump processes with time-dependent rates.
+  In particular, if all `VariableRateJump`s in a pure-jump system are bounded one
+  can use `Coevolve` with `SSAStepper` for better performance. See the
+  documentation, particularly the first and second tutorials, for details on
+  defining and using bounded `VariableRateJump`s.

From 60b2a99949cc2ae946d88d52e1746c62f90a8767 Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Tue, 6 Aug 2024 12:53:37 -0400
Subject: [PATCH 14/21] bug fix and update HISTORY plus README

---
 README.md    | 14 +++++++-------
 src/jumps.jl |  2 +-
 2 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/README.md b/README.md
index 3c34471a..54c1b31c 100644
--- a/README.md
+++ b/README.md
@@ -74,10 +74,10 @@ using JumpProcesses, Plots
 # here we order S = 1, I = 2, and R = 3
 # substrate stoichiometry:
 substoich = [[1 => 1, 2 => 1],    # 1*S + 1*I
-    [2 => 1]]            # 1*I
+    [2 => 1]]                     # 1*I
 # net change by each jump type
 netstoich = [[1 => -1, 2 => 1],   # S -> S-1, I -> I+1
-    [2 => -1, 3 => 1]]   # I -> I-1, R -> R+1
+    [2 => -1, 3 => 1]]            # I -> I-1, R -> R+1
 # rate constants for each jump
 p = (0.1 / 1000, 0.01)
 
@@ -91,10 +91,10 @@ tspan = (0.0, 250.0)
 dprob = DiscreteProblem(u₀, tspan, p)
 
 # use the Direct method to simulate
-jprob = JumpProblem(dprob, Direct(), maj)
+jprob = JumpProblem(dprob, maj)
 
 # solve as a pure jump process, i.e. using SSAStepper
-sol = solve(jprob, SSAStepper())
+sol = solve(jprob)
 plot(sol)
 ```
 
@@ -117,8 +117,8 @@ function affect2!(integrator)
     integrator.u[3] += 1        # R -> R + 1
 end
 jump2 = ConstantRateJump(rate2, affect2!)
-jprob = JumpProblem(dprob, Direct(), jump, jump2)
-sol = solve(jprob, SSAStepper())
+jprob = JumpProblem(dprob, jump, jump2)
+sol = solve(jprob)
 ```
 
 ### Jump-ODE Example
@@ -164,7 +164,7 @@ plot(sol)
 
   - See the [SciML Style Guide](https://github.com/SciML/SciMLStyle) for common coding practices and other style decisions.
   - There are a few community forums for getting help and asking questions:
-    
+
       + The #diffeq-bridged and #sciml-bridged channels in the
         [Julia Slack](https://julialang.org/slack/)
       + The #diffeq-bridged and #sciml-bridged channels in the
diff --git a/src/jumps.jl b/src/jumps.jl
index e09236ff..f3b22f10 100644
--- a/src/jumps.jl
+++ b/src/jumps.jl
@@ -383,7 +383,7 @@ end
 using_params(maj::MassActionJump{T, S, U, Nothing}) where {T, S, U} = false
 using_params(maj::MassActionJump) = true
 using_params(maj::Nothing) = false
-@inline get_num_majumps(maj::MassActionJump) = length(maj.scaled_rates)
+@inline get_num_majumps(maj::MassActionJump) = length(maj.net_stoch)
 @inline get_num_majumps(maj::Nothing) = 0
 
 struct MassActionJumpParamMapper{U}

From 503224b13b04200f7c9896a1340d198d60ae68c0 Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Tue, 6 Aug 2024 12:56:43 -0400
Subject: [PATCH 15/21] tweak wording

---
 HISTORY.md | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/HISTORY.md b/HISTORY.md
index 37987c8f..c015c513 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -15,10 +15,10 @@
   jprob = JumpProblem(dprob, crj)
   sol = solve(jprob, SSAStepper())
   ```
-- For `JumpProblem`s over `DiscreteProblem`s that only have `MassActionJump`s
-  and `ConstantRateJump`s one no longer needs to specify `SSAStepper()` when
-  calling `solve`, i.e. the following now works for the previous example and is
-  equivalent to manually passing `SSAStepper()`:
+- For `JumpProblem`s over `DiscreteProblem`s that only have `MassActionJump`s,
+  `ConstantRateJump`s, and bounded `VariableRateJump`s, one no longer needs to
+  specify `SSAStepper()` when calling `solve`, i.e. the following now works for
+  the previous example and is equivalent to manually passing `SSAStepper()`:
   ```julia
   sol = solve(jprob)
   ```

From 71f2f93881ae6fc7503289415ea131dc34731bca Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Tue, 6 Aug 2024 12:57:47 -0400
Subject: [PATCH 16/21] update SIR tests

---
 test/sir_model.jl | 28 ++++++++++++++++++++++++++++
 1 file changed, 28 insertions(+)

diff --git a/test/sir_model.jl b/test/sir_model.jl
index a64ce2ad..3cd821ec 100644
--- a/test/sir_model.jl
+++ b/test/sir_model.jl
@@ -29,3 +29,31 @@ end
 cb = DiscreteCallback(condition, purge_affect!, save_positions = (false, false))
 sol = solve(jump_prob, FunctionMap(), callback = cb, tstops = [100])
 sol = solve(jump_prob, SSAStepper(), callback = cb, tstops = [100])
+
+# test README example using the auto-solver selection runs
+let
+    # here we order S = 1, I = 2, and R = 3
+    # substrate stoichiometry:
+    substoich = [[1 => 1, 2 => 1],    # 1*S + 1*I
+    [2 => 1]]                     # 1*I
+    # net change by each jump type
+    netstoich = [[1 => -1, 2 => 1],   # S -> S-1, I -> I+1
+    [2 => -1, 3 => 1]]            # I -> I-1, R -> R+1
+    # rate constants for each jump
+    p = (0.1 / 1000, 0.01)
+
+    # p[1] is rate for S+I --> 2I, p[2] for I --> R
+    pidxs = [1, 2]
+
+    maj = MassActionJump(substoich, netstoich; param_idxs = pidxs)
+
+    u₀ = [999, 1, 0]       #[S(0),I(0),R(0)]
+    tspan = (0.0, 250.0)
+    dprob = DiscreteProblem(u₀, tspan, p)
+
+    # use the Direct method to simulate
+    jprob = JumpProblem(dprob, maj)
+
+    # solve as a pure jump process, i.e. using SSAStepper
+    sol = solve(jprob)
+end
\ No newline at end of file

From 42773269aebcf3c1789a41dd656d2ba2de61e04d Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Tue, 6 Aug 2024 12:58:22 -0400
Subject: [PATCH 17/21] Update HISTORY.md

---
 HISTORY.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/HISTORY.md b/HISTORY.md
index c015c513..bac3bef3 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -22,7 +22,7 @@
   ```julia
   sol = solve(jprob)
   ```
-- Plotting a solution generated with `save_positions = (false, false)` no uses
+- Plotting a solution generated with `save_positions = (false, false)` now uses
   piecewise linear plots between any saved time points specified via `saveat`
   instead (previously the plots appeared piecewise constant even though each
   jump was not being shown). Note that solution objects still use piecewise

From ea3b7278bbebca12ea28caeefd24f8c7c89583ed Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Tue, 6 Aug 2024 12:59:13 -0400
Subject: [PATCH 18/21] format

---
 HISTORY.md        | 73 ++++++++++++++++++++++++-----------------------
 README.md         |  2 +-
 test/sir_model.jl |  6 ++--
 3 files changed, 42 insertions(+), 39 deletions(-)

diff --git a/HISTORY.md b/HISTORY.md
index c015c513..954786bc 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -3,44 +3,47 @@
 ## JumpProcesses unreleased (master branch)
 
 ## 9.12
-- Added a default aggregator selection algorithm based on the number of passed
-  in jumps. i.e. the following now auto-selects an aggregator (`Direct` in this
-  case):
-  ```julia
-  using JumpProcesses
-  rate(u,p,t) = u[1]
-  affect(integrator) = (integrator.u[1] -= 1; nothing)
-  crj = ConstantRateJump(rate,affect)
-  dprob = DiscreteProblem([10],(0.0, 10.0))
-  jprob = JumpProblem(dprob, crj)
-  sol = solve(jprob, SSAStepper())
-  ```
-- For `JumpProblem`s over `DiscreteProblem`s that only have `MassActionJump`s,
-  `ConstantRateJump`s, and bounded `VariableRateJump`s, one no longer needs to
-  specify `SSAStepper()` when calling `solve`, i.e. the following now works for
-  the previous example and is equivalent to manually passing `SSAStepper()`:
-  ```julia
-  sol = solve(jprob)
-  ```
-- Plotting a solution generated with `save_positions = (false, false)` no uses
-  piecewise linear plots between any saved time points specified via `saveat`
-  instead (previously the plots appeared piecewise constant even though each
-  jump was not being shown). Note that solution objects still use piecewise
-  constant interpolation, see [the
-  docs](https://docs.sciml.ai/JumpProcesses/stable/tutorials/discrete_stochastic_example/#save_positions_docs)
-  for details.
+
+  - Added a default aggregator selection algorithm based on the number of passed
+    in jumps. i.e. the following now auto-selects an aggregator (`Direct` in this
+    case):
+    
+    ```julia
+    using JumpProcesses
+    rate(u, p, t) = u[1]
+    affect(integrator) = (integrator.u[1] -= 1; nothing)
+    crj = ConstantRateJump(rate, affect)
+    dprob = DiscreteProblem([10], (0.0, 10.0))
+    jprob = JumpProblem(dprob, crj)
+    sol = solve(jprob, SSAStepper())
+    ```
+  - For `JumpProblem`s over `DiscreteProblem`s that only have `MassActionJump`s,
+    `ConstantRateJump`s, and bounded `VariableRateJump`s, one no longer needs to
+    specify `SSAStepper()` when calling `solve`, i.e. the following now works for
+    the previous example and is equivalent to manually passing `SSAStepper()`:
+    
+    ```julia
+    sol = solve(jprob)
+    ```
+  - Plotting a solution generated with `save_positions = (false, false)` no uses
+    piecewise linear plots between any saved time points specified via `saveat`
+    instead (previously the plots appeared piecewise constant even though each
+    jump was not being shown). Note that solution objects still use piecewise
+    constant interpolation, see [the
+    docs](https://docs.sciml.ai/JumpProcesses/stable/tutorials/discrete_stochastic_example/#save_positions_docs)
+    for details.
 
 ## 9.7
 
-- `Coevolve` was updated to support use with coupled ODEs/SDEs. See the updated
-  documentation for details, and note the comments there about one needing to ensure
-  rate bounds hold however the ODE/SDE stepper could modify dependent variables during a timestep.
+  - `Coevolve` was updated to support use with coupled ODEs/SDEs. See the updated
+    documentation for details, and note the comments there about one needing to ensure
+    rate bounds hold however the ODE/SDE stepper could modify dependent variables during a timestep.
 
 ## 9.3
 
-- Support for "bounded" `VariableRateJump`s that can be used with the `Coevolve`
-  aggregator for faster simulation of jump processes with time-dependent rates.
-  In particular, if all `VariableRateJump`s in a pure-jump system are bounded one
-  can use `Coevolve` with `SSAStepper` for better performance. See the
-  documentation, particularly the first and second tutorials, for details on
-  defining and using bounded `VariableRateJump`s.
+  - Support for "bounded" `VariableRateJump`s that can be used with the `Coevolve`
+    aggregator for faster simulation of jump processes with time-dependent rates.
+    In particular, if all `VariableRateJump`s in a pure-jump system are bounded one
+    can use `Coevolve` with `SSAStepper` for better performance. See the
+    documentation, particularly the first and second tutorials, for details on
+    defining and using bounded `VariableRateJump`s.
diff --git a/README.md b/README.md
index 54c1b31c..07adef60 100644
--- a/README.md
+++ b/README.md
@@ -164,7 +164,7 @@ plot(sol)
 
   - See the [SciML Style Guide](https://github.com/SciML/SciMLStyle) for common coding practices and other style decisions.
   - There are a few community forums for getting help and asking questions:
-
+    
       + The #diffeq-bridged and #sciml-bridged channels in the
         [Julia Slack](https://julialang.org/slack/)
       + The #diffeq-bridged and #sciml-bridged channels in the
diff --git a/test/sir_model.jl b/test/sir_model.jl
index 3cd821ec..e8cea455 100644
--- a/test/sir_model.jl
+++ b/test/sir_model.jl
@@ -35,10 +35,10 @@ let
     # here we order S = 1, I = 2, and R = 3
     # substrate stoichiometry:
     substoich = [[1 => 1, 2 => 1],    # 1*S + 1*I
-    [2 => 1]]                     # 1*I
+        [2 => 1]]                     # 1*I
     # net change by each jump type
     netstoich = [[1 => -1, 2 => 1],   # S -> S-1, I -> I+1
-    [2 => -1, 3 => 1]]            # I -> I-1, R -> R+1
+        [2 => -1, 3 => 1]]            # I -> I-1, R -> R+1
     # rate constants for each jump
     p = (0.1 / 1000, 0.01)
 
@@ -56,4 +56,4 @@ let
 
     # solve as a pure jump process, i.e. using SSAStepper
     sol = solve(jprob)
-end
\ No newline at end of file
+end

From 4ebc16f4654faa391383e5ddf1b3855d0199e28a Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Tue, 6 Aug 2024 13:01:26 -0400
Subject: [PATCH 19/21] fix history

---
 HISTORY.md | 29 -----------------------------
 1 file changed, 29 deletions(-)

diff --git a/HISTORY.md b/HISTORY.md
index e7f9b039..c6e97ee9 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -30,35 +30,6 @@
   docs](https://docs.sciml.ai/JumpProcesses/stable/tutorials/discrete_stochastic_example/#save_positions_docs)
   for details.
 
-  - Added a default aggregator selection algorithm based on the number of passed
-    in jumps. i.e. the following now auto-selects an aggregator (`Direct` in this
-    case):
-    
-    ```julia
-    using JumpProcesses
-    rate(u, p, t) = u[1]
-    affect(integrator) = (integrator.u[1] -= 1; nothing)
-    crj = ConstantRateJump(rate, affect)
-    dprob = DiscreteProblem([10], (0.0, 10.0))
-    jprob = JumpProblem(dprob, crj)
-    sol = solve(jprob, SSAStepper())
-    ```
-  - For `JumpProblem`s over `DiscreteProblem`s that only have `MassActionJump`s,
-    `ConstantRateJump`s, and bounded `VariableRateJump`s, one no longer needs to
-    specify `SSAStepper()` when calling `solve`, i.e. the following now works for
-    the previous example and is equivalent to manually passing `SSAStepper()`:
-    
-    ```julia
-    sol = solve(jprob)
-    ```
-  - Plotting a solution generated with `save_positions = (false, false)` no uses
-    piecewise linear plots between any saved time points specified via `saveat`
-    instead (previously the plots appeared piecewise constant even though each
-    jump was not being shown). Note that solution objects still use piecewise
-    constant interpolation, see [the
-    docs](https://docs.sciml.ai/JumpProcesses/stable/tutorials/discrete_stochastic_example/#save_positions_docs)
-    for details.
-
 ## 9.7
 
   - `Coevolve` was updated to support use with coupled ODEs/SDEs. See the updated

From 91ddeafd2f38ae6641387a185d3678c1f4253e31 Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Tue, 6 Aug 2024 13:02:12 -0400
Subject: [PATCH 20/21] format

---
 HISTORY.md | 55 ++++++++++++++++++++++++++++--------------------------
 1 file changed, 29 insertions(+), 26 deletions(-)

diff --git a/HISTORY.md b/HISTORY.md
index c6e97ee9..2e2f25b3 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -3,32 +3,35 @@
 ## JumpProcesses unreleased (master branch)
 
 ## 9.12
-- Added a default aggregator selection algorithm based on the number of passed
-  in jumps. i.e. the following now auto-selects an aggregator (`Direct` in this
-  case):
-  ```julia
-  using JumpProcesses
-  rate(u,p,t) = u[1]
-  affect(integrator) = (integrator.u[1] -= 1; nothing)
-  crj = ConstantRateJump(rate,affect)
-  dprob = DiscreteProblem([10],(0.0, 10.0))
-  jprob = JumpProblem(dprob, crj)
-  sol = solve(jprob, SSAStepper())
-  ```
-- For `JumpProblem`s over `DiscreteProblem`s that only have `MassActionJump`s,
-  `ConstantRateJump`s, and bounded `VariableRateJump`s, one no longer needs to
-  specify `SSAStepper()` when calling `solve`, i.e. the following now works for
-  the previous example and is equivalent to manually passing `SSAStepper()`:
-  ```julia
-  sol = solve(jprob)
-  ```
-- Plotting a solution generated with `save_positions = (false, false)` now uses
-  piecewise linear plots between any saved time points specified via `saveat`
-  instead (previously the plots appeared piecewise constant even though each
-  jump was not being shown). Note that solution objects still use piecewise
-  constant interpolation, see [the
-  docs](https://docs.sciml.ai/JumpProcesses/stable/tutorials/discrete_stochastic_example/#save_positions_docs)
-  for details.
+
+  - Added a default aggregator selection algorithm based on the number of passed
+    in jumps. i.e. the following now auto-selects an aggregator (`Direct` in this
+    case):
+    
+    ```julia
+    using JumpProcesses
+    rate(u, p, t) = u[1]
+    affect(integrator) = (integrator.u[1] -= 1; nothing)
+    crj = ConstantRateJump(rate, affect)
+    dprob = DiscreteProblem([10], (0.0, 10.0))
+    jprob = JumpProblem(dprob, crj)
+    sol = solve(jprob, SSAStepper())
+    ```
+  - For `JumpProblem`s over `DiscreteProblem`s that only have `MassActionJump`s,
+    `ConstantRateJump`s, and bounded `VariableRateJump`s, one no longer needs to
+    specify `SSAStepper()` when calling `solve`, i.e. the following now works for
+    the previous example and is equivalent to manually passing `SSAStepper()`:
+    
+    ```julia
+    sol = solve(jprob)
+    ```
+  - Plotting a solution generated with `save_positions = (false, false)` now uses
+    piecewise linear plots between any saved time points specified via `saveat`
+    instead (previously the plots appeared piecewise constant even though each
+    jump was not being shown). Note that solution objects still use piecewise
+    constant interpolation, see [the
+    docs](https://docs.sciml.ai/JumpProcesses/stable/tutorials/discrete_stochastic_example/#save_positions_docs)
+    for details.
 
 ## 9.7
 

From 331ea57d1efb0748116116b5c6122822ddc668b1 Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Tue, 6 Aug 2024 13:07:32 -0400
Subject: [PATCH 21/21] I hate julia formatter

---
 HISTORY.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/HISTORY.md b/HISTORY.md
index 2e2f25b3..1d2b0f7d 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -17,6 +17,7 @@
     jprob = JumpProblem(dprob, crj)
     sol = solve(jprob, SSAStepper())
     ```
+
   - For `JumpProblem`s over `DiscreteProblem`s that only have `MassActionJump`s,
     `ConstantRateJump`s, and bounded `VariableRateJump`s, one no longer needs to
     specify `SSAStepper()` when calling `solve`, i.e. the following now works for