Spatial SSAs progress and TODOs

This summer I am working on implementing spatial SSAs. There are currently two issues open about it: https://github.com/SciML/DiffEqJump.jl/issues/107 and https://github.com/SciML/DiffEqJump.jl/issues/130. The first PR https://github.com/SciML/DiffEqJump.jl/pull/183 implements the NSM, adds two tests and some utility functions, essentially setting up the interface for spatial SSAs. LightGraphs are supported as well. There are a few TODOs that should not be forgotten:

• The citation must be added to the just-implemented NSM.
• Store numspecies in aggregation?
• More comprehensive tests for spatial SSAs are needed. The two current tests only test spatial systems on 1D grids.
• The CartesianGrid struct needs to be optimized -- probably by storing the number of neighbors for each site so num_neighbors is O(1) and making neighbors return an iterator instead of allocating an array.
• Make CartesianGrid use CartesianIndices instead of doing all the calculations redundantly.
• An interface for sampling a neighbor to hop to can be a nice but not necessary perk. For example, using https://docs.julialang.org/en/v1/stdlib/Random/#Random.Sampler.
• Various modes of hopping constants need to be supported. In particular, we should support hopping rates of forms D_{s,i}, D_s * L_{i,j}, L_{s,i,j}, D_{s,i}*L_{i,j}, where D_s is the diffusivity constant of species s and L is the discrete Laplacian matrix. Only the first of these is implemented in the PR. This change will require tweaking the code in nsm.jl.
• The escaping boundary condition for CartesianGrid needs to be supported. This can be done with a special sink site: it neighbors all boundary sites and does not get updated if it is the target site of a hop. This change will require tweaking the code in nsm.jl.
• Constant-rate jumps should be supported. Currently, only mass-action jumps are allowed. This change will require tweaking the code in nsm.jl.
• Implement the flattening SSA: turn all spatial hops into reactions and simulate using a non-spatial SSA. Can modify https://github.com/SciML/DiffEqJump.jl/pull/131 to do that.
• Accept variable rx rates in flattened SSA.
• Optimize allocations in HopRatesGeneral: switch to using cumulative sums to sample neighbor.
• implement spatially varying massaction jump rates, including flattening.
• add docstrings on all structs and functions explaining what things are.
• make a tutorial for setting up a spatial problem in different ways (cartesian grid, arbitrary graph) and solving with different solvers (NSM, DirectCRonDirect, flattening). For example, use A+B <--> C. Make a PR to SciMLTutorials.jl.
• setup a comprehensive benchmarking suite for spatial solvers and make a PR to SciMLBenchmarks.jl
• add jump counter to SSAStepper to count jumps. This can be easily done with a callback as in the benchmark.
• Choose a default SSA (probably DirectCRRSSA). Can also do it for non-spatial solvers.
• compare with other software like urdme to see how the spatial solvers in DiffEqJump stack up.
• Optimize CartesianGrid. It should be faster than a Graph.
• Compute the hopping rates out of the site for each site at the aggregation step of the simulation (instead of essentially recomputing it every time we evaluate the hop rate, e.g. here https://github.com/SciML/DiffEqJump.jl/blob/master/src/spatial/hop_rates.jl#L115). Take, for example, hopping rates of form Dsi. The quantity hop_rates.hopping_constants[species,site]*outdegree(spatial_system, site) is always the same for any given site. If we precompute this quantity we won't need to store the number of neighbors for each site in the CartesianGrid struct resulting in less memory used. And we would also avoid one multiplication and one pointer access for each update of hopping rates. [EDIT]: maybe store the outdegrees as a dict, where all the interior sites do not appear. This would make the nums_neighbors array O(n^2) instead of O(n^3) where n is linear dimension of the grid.
• Add more tests to thoroughly test the spatial functionality.
• Add more spatial solvers. DirectCR-RSSA.
• (Longer-term) Update the tutorial https://tutorials.sciml.ai/html/jumps/spatial.html. Add more known simple examples with pretty pictures. More visual examples, comparing stochastic solution to deterministic solution. Maybe add a spatial SIR example. Maybe use Catalyst for ease of understanding.
• (Longer-term) Updating the DifferentialEquations.jl https://diffeq.sciml.ai/stable/tutorials/discrete_stochastic_example/ with some spatial examples.
• (Longer-term) Add benchmarks of models from other papers like https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1005865.

@isaacsas and @ChrisRackauckas might have more TODOs to add to this list.

Start by selecting one unchecked TODO rather than treating this issue as a single task. The payload points to src/spatial/nsm.jl and src/spatial/hop_rates.jl, along with spatial tests and tutorials; read the relevant entry point and existing tests first. Done should mean one scoped TODO is implemented, covered by tests, and its documentation or benchmark updates are made where applicable.