diff --git a/src/spatial/hop_rates.jl b/src/spatial/hop_rates.jl
index 26471d38..85d96395 100644
--- a/src/spatial/hop_rates.jl
+++ b/src/spatial/hop_rates.jl
@@ -248,7 +248,7 @@ end
 Analogue of HopRatesGraphDsij, optimized for CartesianGrid
 """
 struct HopRatesGridDsij{F} <: AbstractHopRates
-    "hop_const_cumulative_sums[:,s,i] is the vector of cumulative sums of hopping constants of species s at site i"
+    "hop_const_cumulative_sums[:,s,i] is the vector of cumulative sums of hopping constants of species s at site i. Out-of-bounds neighbors are treated as having zero propensity."
     hop_const_cumulative_sums::Array{F, 3}
 
     "rates[s,i] is the total hopping rate of species s at site i"
diff --git a/test/runtests.jl b/test/runtests.jl
index 067b3525..1c89f66f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -21,7 +21,9 @@ using JumpProcesses, DiffEqBase, SafeTestsets
     @time @safetestset "A + B <--> C" begin include("reversible_binding.jl") end
     @time @safetestset "Remake tests" begin include("remake_test.jl") end
     @time @safetestset "Long time accuracy test" begin include("longtimes_test.jl") end
-    @time @safetestset "Spatial utilities" begin include("spatial/utils_test.jl") end
+    @time @safetestset "Reaction rates" begin include("spatial/reaction_rates.jl") end
+    @time @safetestset "Hop rates" begin include("spatial/hop_rates.jl") end
+    @time @safetestset "Topology" begin include("spatial/topology.jl") end
     @time @safetestset "Spatial A + B <--> C" begin include("spatial/ABC.jl") end
     @time @safetestset "Spatially Varying Reaction Rates" begin include("spatial/spatial_majump.jl") end
     @time @safetestset "Pure diffusion" begin include("spatial/diffusion.jl") end
diff --git a/test/spatial/hop_rates.jl b/test/spatial/hop_rates.jl
new file mode 100644
index 00000000..be2f2825
--- /dev/null
+++ b/test/spatial/hop_rates.jl
@@ -0,0 +1,194 @@
+using JumpProcesses, Graphs, Test, Random, StableRNGs
+const JP = JumpProcesses
+
+# Functions to test
+# HopRates(hop_constants, spatial_system)
+# update_hop_rates!(hop_rates::AbstractHopRates, species::AbstractArray, u, site, spatial_system)
+# total_site_hop_rate(hop_rates::AbstractHopRates, site)
+# reset!(hop_rates::AbstractHopRates)
+# sample_hop_at_site(hop_rates::AbstractHopRates, site, rng, spatial_system)
+
+function test_reset(hop_rates, num_nodes)
+    JP.reset!(hop_rates)
+    for site in 1:num_nodes
+        @test JP.total_site_hop_rate(hop_rates, site) == 0.0
+    end
+end
+
+function normalized(distribution::Dict)
+    d = copy(distribution)
+    s = sum(values(d))
+    for key in keys(d)
+        d[key] /= s
+    end
+    return d
+end
+
+normalized(distribution) = distribution / sum(distribution)
+
+function statistical_test(hop_rates, spec_propensities, target_propensities::Dict,
+                          num_species, u, site, g, rng, rel_tol)
+    spec_probs = normalized(spec_propensities)
+    target_probs = normalized(target_propensities)
+    JP.update_hop_rates!(hop_rates, 1:num_species, u, site, g)
+    @test JP.total_site_hop_rate(hop_rates, site) == sum(spec_propensities)
+    spec_dict = Dict{Int, Int}()
+    site_dict = Dict{Int, Int}()
+    for _ in 1:num_samples
+        spec, target = JP.sample_hop_at_site(hop_rates, site, rng, g)
+        spec_dict[spec] = get(spec_dict, spec, 0) + 1
+        site_dict[target] = get(site_dict, target, 0) + 1
+    end
+    for species in 1:num_species
+        @test abs(spec_dict[species]) / num_samples - spec_probs[species] < rel_tol
+    end
+    for target in JP.neighbors(g, site)
+        @test abs(site_dict[target]) / num_samples - target_probs[target] < rel_tol
+    end
+end
+
+io = IOBuffer()
+rng = StableRNG(12345)
+rel_tol = 0.05
+num_samples = 10^4
+dims = (3, 2)
+g = CartesianGridRej(dims)
+num_nodes = JP.num_sites(g)
+num_species = 2
+u = ones(Int, num_species, num_nodes)
+
+# Tests for HopRatesGraphDs
+hop_constants = [1.0, 3.0] # species
+hop_rates = JP.HopRates(hop_constants, g)
+@test hop_rates isa JP.HopRatesGraphDs
+show(io, "text/plain", hop_rates)
+for site in 1:num_nodes
+    spec_propensities = hop_constants * JP.outdegree(g, site)
+    target_propensities = Dict{Int, Float64}()
+    for (i, target) in enumerate(JP.neighbors(g, site))
+        target_propensities[target] = 1.0
+    end
+    statistical_test(hop_rates, spec_propensities, target_propensities, num_species, u,
+                     site, g, rng, rel_tol)
+end
+test_reset(hop_rates, num_nodes)
+
+# Tests for HopRatesGraphDsi
+hop_constants = [2.0 1.0 1.0 1.0 1.0 1.0; 6.0 3.0 3.0 3.0 3.0 3.0] # [species, site]
+hop_rates = JP.HopRates(hop_constants, g)
+@test hop_rates isa JP.HopRatesGraphDsi
+show(io, "text/plain", hop_rates)
+for site in 1:num_nodes
+    spec_propensities = hop_constants[:, site] * JP.outdegree(g, site)
+    target_propensities = Dict{Int, Float64}()
+    for (i, target) in enumerate(JP.neighbors(g, site))
+        target_propensities[target] = 1.0
+    end
+    statistical_test(hop_rates, spec_propensities, target_propensities, num_species, u,
+                     site, g, rng, rel_tol)
+end
+test_reset(hop_rates, num_nodes)
+
+# Tests for HopRatesGraphDsij
+hop_constants = Matrix{Vector{Float64}}(undef, num_species, num_nodes) # [species, site][target_site]
+hop_constants[1, :] = [
+    [2.0, 4.0],
+    [1.0, 2.0, 4.0],
+    [1.0, 2.0],
+    [1.0, 2.0],
+    [1.0, 2.0, 4.0],
+    [1.0, 2.0],
+]
+hop_constants[2, :] = [
+    [3.0, 12.0],
+    [3.0, 6.0, 12.0],
+    [3.0, 6.0],
+    [3.0, 6.0],
+    [3.0, 6.0, 12.0],
+    [3.0, 6.0],
+]
+hop_rates_structs = [
+    JP.HopRatesGraphDsij(hop_constants),
+    JP.HopRates(hop_constants, g),
+]
+@test hop_rates_structs[2] isa JP.HopRatesGridDsij
+for hop_rates in hop_rates_structs
+    show(io, "text/plain", hop_rates)
+    for site in 1:num_nodes
+        spec_propensities = [sum(hop_constants[species, site]) for species in 1:num_species]
+        target_propensities = Dict{Int, Float64}()
+        for (i, target) in enumerate(JP.neighbors(g, site))
+            target_propensities[target] = sum([hop_constants[species, site][i]
+                                               for species in 1:num_species])
+        end
+        statistical_test(hop_rates, spec_propensities, target_propensities, num_species, u,
+                         site, g, rng, rel_tol)
+    end
+end
+test_reset(hop_rates, num_nodes)
+
+# Tests for HopRatesGraphDsLij
+species_hop_constants = [1.0, 3.0] #species
+site_hop_constants = [
+    [2.0, 4.0],
+    [1.0, 2.0, 4.0],
+    [1.0, 2.0],
+    [1.0, 2.0],
+    [1.0, 2.0, 4.0],
+    [1.0, 2.0],
+] #[site][target_site]
+hop_rates_structs = [
+    JP.HopRatesGraphDsLij(species_hop_constants, site_hop_constants),
+    JP.HopRates((species_hop_constants => site_hop_constants), g),
+]
+@test hop_rates_structs[2] isa JP.HopRatesGridDsLij
+for hop_rates in hop_rates_structs
+    show(io, "text/plain", hop_rates)
+    for site in 1:num_nodes
+        spec_propensities = [species_hop_constants[species] * sum(site_hop_constants[site])
+                             for species in 1:num_species]
+        target_propensities = Dict{Int, Float64}()
+        for (i, target) in enumerate(JP.neighbors(g, site))
+            target_propensities[target] = sum([species_hop_constants[species] *
+                                               site_hop_constants[site][i]
+                                               for species in 1:num_species])
+        end
+        statistical_test(hop_rates, spec_propensities, target_propensities, num_species, u,
+                         site, g, rng, rel_tol)
+    end
+end
+test_reset(hop_rates, num_nodes)
+
+# Tests for HopRatesGraphDsiLij
+species_hop_constants = [2.0 1.0 1.0 1.0 1.0 1.0; 6.0 3.0 3.0 3.0 3.0 3.0] #[species, site]
+site_hop_constants = [
+    [2.0, 4.0],
+    [1.0, 2.0, 4.0],
+    [1.0, 2.0],
+    [1.0, 2.0],
+    [1.0, 2.0, 4.0],
+    [1.0, 2.0],
+] #[site][target_site]
+hop_rates_structs = [
+    JP.HopRatesGraphDsiLij(species_hop_constants, site_hop_constants),
+    JP.HopRates((species_hop_constants => site_hop_constants), g),
+]
+@test hop_rates_structs[2] isa JP.HopRatesGridDsiLij
+for hop_rates in hop_rates_structs
+    show(io, "text/plain", hop_rates)
+    for site in 1:num_nodes
+        spec_propensities = [species_hop_constants[species, site] *
+                             sum(site_hop_constants[site]) for species in 1:num_species]
+        target_propensities = Dict{Int, Float64}()
+        for (i, target) in enumerate(JP.neighbors(g, site))
+            target_propensities[target] = sum([species_hop_constants[species, site] *
+                                               site_hop_constants[site][i]
+                                               for species in 1:num_species])
+        end
+        statistical_test(hop_rates, spec_propensities, target_propensities, num_species, u,
+                         site, g, rng, rel_tol)
+    end
+end
+test_reset(hop_rates, num_nodes)
+
+@test String(take!(io)) !== nothing
diff --git a/test/spatial/reaction_rates.jl b/test/spatial/reaction_rates.jl
new file mode 100644
index 00000000..1f62e8eb
--- /dev/null
+++ b/test/spatial/reaction_rates.jl
@@ -0,0 +1,50 @@
+using JumpProcesses, Test, Random, StableRNGs
+const JP = JumpProcesses
+
+# Functions to test:
+# num_rxs
+# reset!
+# total_site_rx_rate
+# update_rx_rates!
+# update_rx_rates!
+# sample_rx_at_site
+
+io = IOBuffer()
+# setup of A + B <--> C
+rel_tol = 0.05
+num_samples = 10^4
+num_nodes = 27
+num_species = 3
+reactstoch = [[1 => 1, 2 => 1], [3 => 1]]
+netstoch = [[1 => -1, 2 => -1, 3 => 1], [1 => 1, 2 => 1, 3 => -1]]
+rates = [0.1, 1.0]
+num_rxs = length(rates)
+ma_jumps = MassActionJump(rates, reactstoch, netstoch)
+spatial_ma_jumps = SpatialMassActionJump(rates, reactstoch, netstoch)
+u = ones(Int, num_species, num_nodes)
+rng = StableRNG(12345)
+
+# Tests for RxRates
+rx_rates_list = [JP.RxRates(num_nodes, ma_jumps), JP.RxRates(num_nodes, spatial_ma_jumps)]
+for rx_rates in rx_rates_list
+    @test JP.num_rxs(rx_rates) == length(rates)
+    show(io, "text/plain", rx_rates)
+    for site in 1:num_nodes
+        JP.update_rx_rates!(rx_rates, 1:num_rxs, u, site)
+        @test JP.total_site_rx_rate(rx_rates, site) == 1.1
+        rx_props = [JP.evalrxrate(u[:, site], rx, ma_jumps) for rx in 1:num_rxs]
+        rx_probs = rx_props / sum(rx_props)
+        d = Dict{Int, Int}()
+        for i in 1:num_samples
+            rx = JP.sample_rx_at_site(rx_rates, site, rng)
+            rx in keys(d) ? d[rx] += 1 : d[rx] = 1
+        end
+        for (k, v) in d
+            @test abs(v / num_samples - rx_probs[k]) < rel_tol
+        end
+    end
+    JP.reset!(rx_rates)
+    for site in 1:num_nodes
+        @test JP.total_site_rx_rate(rx_rates, site) == 0.0
+    end
+end
diff --git a/test/spatial/run_spatial_tests.jl b/test/spatial/run_spatial_tests.jl
index 5d74ec8c..75f43161 100644
--- a/test/spatial/run_spatial_tests.jl
+++ b/test/spatial/run_spatial_tests.jl
@@ -1,7 +1,9 @@
 using JumpProcesses, DiffEqBase, SafeTestsets
 
 @time begin
-    @time @safetestset "Spatial utilities" begin include("utils_test.jl") end
+    @time @safetestset "Reaction rates" begin include("reaction_rates.jl") end
+    @time @safetestset "Hop rates" begin include("hop_rates.jl") end
+    @time @safetestset "Topology" begin include("topology.jl") end
     @time @safetestset "Spatial A + B <--> C" begin include("ABC.jl") end
     @time @safetestset "Pure diffusion" begin include("diffusion.jl") end
     @time @safetestset "Spatially Varying Reaction Rates" begin include("spatial_majump.jl") end
diff --git a/test/spatial/topology.jl b/test/spatial/topology.jl
new file mode 100644
index 00000000..7d77d325
--- /dev/null
+++ b/test/spatial/topology.jl
@@ -0,0 +1,44 @@
+# Tests for topology.jl
+using JumpProcesses, Graphs, Test, Random, StableRNGs
+const JP = JumpProcesses
+
+# Functions to test:
+# num_sites(grid)
+# neighbors(grid, site)
+# outdegree(grid, site)
+# rand_nbr(rng, grid, site)
+# nth_nbr(grid, site, n)
+
+io = IOBuffer()
+rng = StableRNG(12345)
+dims = (4, 3, 2)
+sites = 1:prod(dims)
+num_samples = 10^5
+rel_tol = 0.01
+grids = [
+    JP.CartesianGridRej(dims),
+    JP.CartesianGridIter(dims),
+    Graphs.grid(dims),
+]
+for grid in grids
+    show(io, "text/plain", grid)
+    @test String(take!(io)) !== nothing
+    @test JP.num_sites(grid) == prod(dims)
+    @test JP.outdegree(grid, 1) == 3
+    @test JP.outdegree(grid, 4) == 3
+    @test JP.outdegree(grid, 17) == 4
+    @test JP.outdegree(grid, 21) == 3
+    @test JP.outdegree(grid, 6) == 5
+    for site in sites
+        @test [JP.nth_nbr(grid, site, n) for n in 1:outdegree(grid, site)] ==
+              collect(neighbors(grid, site))
+        d = Dict{Int, Int}()
+        for i in 1:num_samples
+            nb = JP.rand_nbr(rng, grid, site)
+            nb in keys(d) ? d[nb] += 1 : d[nb] = 1
+        end
+        for val in values(d)
+            @test abs(val / num_samples - 1 / JP.outdegree(grid, site)) < rel_tol
+        end
+    end
+end
diff --git a/test/spatial/utils_test.jl b/test/spatial/utils_test.jl
deleted file mode 100644
index 2e38d626..00000000
--- a/test/spatial/utils_test.jl
+++ /dev/null
@@ -1,201 +0,0 @@
-# Tests for CartesianGrid
-using JumpProcesses, Graphs
-using Test, Random
-
-io = IOBuffer()
-rng = MersenneTwister()
-dims = (4, 3, 2)
-sites = rand(1:prod(dims), 10)
-num_samples = 10^5
-rel_tol = 0.01
-grids = [
-    JumpProcesses.CartesianGridRej(dims),
-    JumpProcesses.CartesianGridIter(dims),
-    Graphs.grid(dims),
-]
-for grid in grids
-    show(io, "text/plain", grid)
-    @test String(take!(io)) !== nothing
-    @test JumpProcesses.num_sites(grid) == prod(dims)
-    @test JumpProcesses.outdegree(grid, 1) == 3
-    @test JumpProcesses.outdegree(grid, 4) == 3
-    @test JumpProcesses.outdegree(grid, 17) == 4
-    @test JumpProcesses.outdegree(grid, 21) == 3
-    @test JumpProcesses.outdegree(grid, 6) == 5
-    for site in sites
-        d = Dict{Int, Int}()
-        for i in 1:num_samples
-            nb = JumpProcesses.rand_nbr(rng, grid, site)
-            nb in keys(d) ? d[nb] += 1 : d[nb] = 1
-        end
-        for val in values(d)
-            @test abs(val / num_samples - 1 / JumpProcesses.outdegree(grid, site)) < rel_tol
-        end
-    end
-end
-
-# setup
-rel_tol = 0.05
-num_samples = 10^4
-dims = (3, 3, 3)
-g = CartesianGridRej(dims)
-num_nodes = JumpProcesses.num_sites(g)
-num_species = 3
-reactstoch = [[1 => 1, 2 => 1], [3 => 1]]
-netstoch = [[1 => -1, 2 => -1, 3 => 1], [1 => 1, 2 => 1, 3 => -1]]
-rates = [0.1, 1.0]
-num_rxs = length(rates)
-ma_jumps = MassActionJump(rates, reactstoch, netstoch)
-u = ones(Int, num_species, num_nodes)
-rng = MersenneTwister()
-
-# Tests for RxRates
-rx_rates = JumpProcesses.RxRates(num_nodes, ma_jumps)
-show(io, "text/plain", rx_rates)
-for site in 1:num_nodes
-    JumpProcesses.update_rx_rates!(rx_rates, 1:num_rxs, u, site)
-    rx_props = [JumpProcesses.evalrxrate(u[:, site], rx, ma_jumps) for rx in 1:num_rxs]
-    rx_probs = rx_props / sum(rx_props)
-    d = Dict{Int, Int}()
-    for i in 1:num_samples
-        rx = JumpProcesses.sample_rx_at_site(rx_rates, site, rng)
-        rx in keys(d) ? d[rx] += 1 : d[rx] = 1
-    end
-    for (k, v) in d
-        @test abs(v / num_samples - rx_probs[k]) < rel_tol
-    end
-end
-
-# Tests for HopRatesGraphDs
-hopping_constants = ones(num_species)
-hop_rates = JumpProcesses.HopRatesGraphDs(hopping_constants, num_nodes)
-show(io, "text/plain", hop_rates)
-spec_probs = ones(num_species) / num_species
-
-for site in 1:num_nodes
-    JumpProcesses.update_hop_rates!(hop_rates, 1:num_species, u, site, g)
-    num_nbs = JumpProcesses.outdegree(g, site)
-    target_probs = ones(num_nbs) / num_nbs
-    d1 = Dict{Int, Int}()
-    d2 = Dict{Int, Int}()
-    for i in 1:num_samples
-        spec, target = JumpProcesses.sample_hop_at_site(hop_rates, site, rng, g)
-        d1[spec] = get(d1, spec, 0) + 1
-        d2[target] = get(d2, target, 0) + 1
-    end
-    @test maximum(abs.(collect(values(d1)) / num_samples - spec_probs)) < rel_tol
-    @test maximum(abs.(collect(values(d2)) / num_samples - target_probs)) < rel_tol
-end
-
-# Tests for HopRatesGraphDsi
-hopping_constants = ones(num_species, num_nodes)
-hop_rates = JumpProcesses.HopRatesGraphDsi(hopping_constants)
-show(io, "text/plain", hop_rates)
-spec_probs = ones(num_species) / num_species
-
-for site in 1:num_nodes
-    JumpProcesses.update_hop_rates!(hop_rates, 1:num_species, u, site, g)
-    num_nbs = JumpProcesses.outdegree(g, site)
-    target_probs = ones(num_nbs) / num_nbs
-    d1 = Dict{Int, Int}()
-    d2 = Dict{Int, Int}()
-    for i in 1:num_samples
-        spec, target = JumpProcesses.sample_hop_at_site(hop_rates, site, rng, g)
-        d1[spec] = get(d1, spec, 0) + 1
-        d2[target] = get(d2, target, 0) + 1
-    end
-    @test maximum(abs.(collect(values(d1)) / num_samples - spec_probs)) < rel_tol
-    @test maximum(abs.(collect(values(d2)) / num_samples - target_probs)) < rel_tol
-end
-
-# Tests for HopRatesGraphDsij
-hop_constants = Matrix{Vector{Float64}}(undef, num_species, num_nodes)
-for ci in CartesianIndices(hop_constants)
-    (species, site) = Tuple(ci)
-    hop_constants[ci] = repeat([1.0], JumpProcesses.outdegree(g, site))
-end
-spec_probs = ones(num_species) / num_species
-hop_rates_structs = [
-    JumpProcesses.HopRatesGraphDsij(hop_constants),
-    JumpProcesses.HopRatesGridDsij(hop_constants, g),
-]
-
-for hop_rates in hop_rates_structs
-    show(io, "text/plain", hop_rates)
-    for site in 1:num_nodes
-        JumpProcesses.update_hop_rates!(hop_rates, 1:num_species, u, site, g)
-        num_nbs = JumpProcesses.outdegree(g, site)
-        target_probs = ones(num_nbs) / num_nbs
-        d1 = Dict{Int, Int}()
-        d2 = Dict{Int, Int}()
-        for i in 1:num_samples
-            spec, target = JumpProcesses.sample_hop_at_site(hop_rates, site, rng, g)
-            d1[spec] = get(d1, spec, 0) + 1
-            d2[target] = get(d2, target, 0) + 1
-        end
-        @test maximum(abs.(collect(values(d1)) / num_samples - spec_probs)) < rel_tol
-        @test maximum(abs.(collect(values(d2)) / num_samples - target_probs)) < rel_tol
-    end
-end
-
-# Tests for HopRatesGraphDsLij
-species_hop_constants = ones(num_species)
-site_hop_constants = Vector{Vector{Float64}}(undef, num_nodes)
-for site in 1:num_nodes
-    site_hop_constants[site] = repeat([1.0], JumpProcesses.outdegree(g, site))
-end
-spec_probs = ones(num_species) / num_species
-hop_rates_structs = [
-    JumpProcesses.HopRatesGraphDsLij(species_hop_constants, site_hop_constants),
-    JumpProcesses.HopRatesGridDsLij(species_hop_constants, site_hop_constants, g),
-]
-
-for hop_rates in hop_rates_structs
-    show(io, "text/plain", hop_rates)
-    for site in 1:num_nodes
-        JumpProcesses.update_hop_rates!(hop_rates, 1:num_species, u, site, g)
-        num_nbs = JumpProcesses.outdegree(g, site)
-        target_probs = ones(num_nbs) / num_nbs
-        d1 = Dict{Int, Int}()
-        d2 = Dict{Int, Int}()
-        for i in 1:num_samples
-            spec, target = JumpProcesses.sample_hop_at_site(hop_rates, site, rng, g)
-            d1[spec] = get(d1, spec, 0) + 1
-            d2[target] = get(d2, target, 0) + 1
-        end
-        @test maximum(abs.(collect(values(d1)) / num_samples - spec_probs)) < rel_tol
-        @test maximum(abs.(collect(values(d2)) / num_samples - target_probs)) < rel_tol
-    end
-end
-
-# Tests for HopRatesGraphDsiLij
-species_hop_constants = ones(num_species, num_nodes)
-site_hop_constants = Vector{Vector{Float64}}(undef, num_nodes)
-for site in 1:num_nodes
-    site_hop_constants[site] = repeat([1.0], JumpProcesses.outdegree(g, site))
-end
-spec_probs = ones(num_species) / num_species
-hop_rates_structs = [
-    JumpProcesses.HopRatesGraphDsiLij(species_hop_constants, site_hop_constants),
-    JumpProcesses.HopRatesGridDsiLij(species_hop_constants, site_hop_constants, g),
-]
-
-for hop_rates in hop_rates_structs
-    show(io, "text/plain", hop_rates)
-    for site in 1:num_nodes
-        JumpProcesses.update_hop_rates!(hop_rates, 1:num_species, u, site, g)
-        num_nbs = JumpProcesses.outdegree(g, site)
-        target_probs = ones(num_nbs) / num_nbs
-        d1 = Dict{Int, Int}()
-        d2 = Dict{Int, Int}()
-        for i in 1:num_samples
-            spec, target = JumpProcesses.sample_hop_at_site(hop_rates, site, rng, g)
-            d1[spec] = get(d1, spec, 0) + 1
-            d2[target] = get(d2, target, 0) + 1
-        end
-        @test maximum(abs.(collect(values(d1)) / num_samples - spec_probs)) < rel_tol
-        @test maximum(abs.(collect(values(d2)) / num_samples - target_probs)) < rel_tol
-    end
-end
-
-@test String(take!(io)) !== nothing