Remove Random.seed! from test suite, use Xoshiro/StableRNG instead (#2779)

Closes #2726

Global state is bad, so I went through the test suite and replaced all 
`Random.seed!` calls with local alternatives.

The replacement depends on what the test actually needs:

- If the test only checks that results are roughly correct (using
`atol=`
or helpers like `check_gdemo`), I used `Xoshiro` passed directly to
`sample`.
- If the test checks exact reproducibility (using `==` or `isequal`), 
  I used `StableRNG` instead.
- If the seed was before optimization code like MLE/MAP (which is 
  deterministic anyway), I just deleted it.
- For `MCMCThreads` tests that reset the same rng between runs, I kept 
  the `Random.seed!(rng, seed)` pattern since that's already using an 
  explicit rng and not global state.

I left the single `Random.seed!` in `runtests.jl` untouched since 
that one is intentional (as discussed in the issue).

---------

Co-authored-by: Penelope Yong <penelopeysm@gmail.com>

Author: anurag-mds

test/ext/dynamichmc.jl

@@ -7,12 +7,13 @@ using Distributions: sample
 using DynamicHMC: DynamicHMC
 using DynamicPPL: DynamicPPL
 using Random: Random
+using StableRNGs: StableRNG
 using Turing
 
 @testset "TuringDynamicHMCExt" begin
-    Random.seed!(100)
+    rng = StableRNG(100)
     spl = externalsampler(DynamicHMC.NUTS())
-    chn = sample(gdemo_default, spl, 10_000)
+    chn = sample(rng, gdemo_default, spl, 10_000)
     check_gdemo(chn)
 end
 

test/mcmc/Inference.jl

@@ -10,6 +10,7 @@ import ForwardDiff
 using LinearAlgebra: I
 import MCMCChains
 import Random
+using Random: Xoshiro
 import ReverseDiff
 using StableRNGs: StableRNG
 using StatsFuns: logsumexp
@@ -35,11 +36,11 @@ using Turing
                 Gibbs(:s => HMC(0.1, 5), :m => ESS()),
             )
             for sampler in samplers
-                Random.seed!(5)
-                chain1 = sample(model, sampler, MCMCThreads(), 10, 4)
+                rng1 = Xoshiro(5)
+                chain1 = sample(rng1, model, sampler, MCMCThreads(), 10, 4)
 
-                Random.seed!(5)
-                chain2 = sample(model, sampler, MCMCThreads(), 10, 4)
+                rng2 = Xoshiro(5)
+                chain2 = sample(rng2, model, sampler, MCMCThreads(), 10, 4)
 
                 # For HMC, the first step does not have stats, so we need to use isequal to
                 # avoid comparing `missing`s

test/mcmc/emcee.jl

@@ -4,19 +4,20 @@ using ..Models: gdemo_default
 using ..NumericalTests: check_gdemo
 using Distributions: sample
 using DynamicPPL: DynamicPPL
-using Random: Random
+using Random: Random, Xoshiro
+using StableRNGs: StableRNG
 using Test: @test, @test_throws, @testset
 using Turing
 
 @testset "emcee.jl" begin
     @testset "gdemo" begin
-        Random.seed!(9876)
+        rng = StableRNG(9876)
 
         n_samples = 1000
         n_walkers = 250
 
         spl = Emcee(n_walkers, 2.0)
-        chain = sample(gdemo_default, spl, n_samples)
+        chain = sample(rng, gdemo_default, spl, n_samples)
         check_gdemo(chain)
     end
 
@@ -25,18 +26,18 @@ using Turing
         @info "Testing emcee with large number of iterations"
         spl = Emcee(10, 2.0)
         n_samples = 10_000
-        chain = sample(gdemo_default, spl, n_samples)
+        chain = sample(StableRNG(5), gdemo_default, spl, n_samples)
         check_gdemo(chain)
     end
 
     @testset "initial parameters" begin
         nwalkers = 250
         spl = Emcee(nwalkers, 2.0)
 
-        Random.seed!(1234)
-        chain1 = sample(gdemo_default, spl, 1)
-        Random.seed!(1234)
-        chain2 = sample(gdemo_default, spl, 1)
+        rng1 = Xoshiro(1234)
+        chain1 = sample(rng1, gdemo_default, spl, 1)
+        rng2 = Xoshiro(1234)
+        chain2 = sample(rng2, gdemo_default, spl, 1)
         @test Array(chain1) == Array(chain2)
 
         initial_nt = DynamicPPL.InitFromParams((s=2.0, m=1.0))

test/mcmc/gibbs.jl

@@ -13,7 +13,7 @@ using Distributions: InverseGamma, Normal
 using Distributions: sample
 using DynamicPPL: DynamicPPL
 using ForwardDiff: ForwardDiff
-using Random: Random
+using Random: Random, Xoshiro
 using ReverseDiff: ReverseDiff
 using StableRNGs: StableRNG
 using Test: @inferred, @test, @test_broken, @test_throws, @testset
@@ -262,10 +262,8 @@ end
     sampler2 = Gibbs(
         @varname(s) => MH(), @varname(s) => MH(), @varname(s) => MH(), @varname(m) => ESS()
     )
-    Random.seed!(23)
-    chain1 = sample(gdemo_default, sampler1, 10)
-    Random.seed!(23)
-    chain2 = sample(gdemo_default, sampler1, 10)
+    chain1 = sample(Xoshiro(23), gdemo_default, sampler1, 10)
+    chain2 = sample(Xoshiro(23), gdemo_default, sampler1, 10)
     @test chain1.value == chain2.value
 end
 
@@ -681,8 +679,9 @@ end
                 # Sampler to use for Gibbs components.
                 hmc = HMC(0.1, 32)
                 sampler = Gibbs(@varname(s) => hmc, @varname(m) => hmc)
-                Random.seed!(42)
+                rng = StableRNG(42)
                 chain = sample(
+                    rng,
                     model,
                     sampler,
                     MCMCThreads(),
@@ -696,8 +695,9 @@ end
 
                 # "Ground truth" samples.
                 # TODO: Replace with closed-form sampling once that is implemented in DynamicPPL.
-                Random.seed!(42)
+
                 chain_true = sample(
+                    StableRNG(42),
                     model,
                     NUTS(),
                     MCMCThreads(),
@@ -742,8 +742,8 @@ end
             end
 
             # `sample`
-            Random.seed!(42)
-            chain = sample(model, spl, 1_000; progress=false)
+            rng = StableRNG(42)
+            chain = sample(rng, model, spl, 1_000; progress=false)
             check_numerical(chain, [:s, :m], [49 / 24, 7 / 6]; atol=0.4)
         end
 
@@ -820,8 +820,8 @@ end
         end
 
         # Sample!
-        Random.seed!(42)
-        chain = sample(MoGtest_default, spl, 1000; progress=false)
+        rng = StableRNG(42)
+        chain = sample(rng, MoGtest_default, spl, 1000; progress=false)
         check_MoGtest_default(chain; atol=0.2)
     end
 
@@ -846,8 +846,8 @@ end
         end
 
         # Sample!
-        Random.seed!(42)
-        chain = sample(model, spl, 1000; progress=false)
+        rng = StableRNG(42)
+        chain = sample(rng, model, spl, 1000; progress=false)
         check_MoGtest_default_z_vector(chain; atol=0.2)
     end
 
@@ -883,9 +883,9 @@ end
         ]
         @testset "$(sampler_inner)" for sampler_inner in samplers_inner
             sampler = Gibbs(@varname(m1) => sampler_inner, @varname(m2) => sampler_inner)
-            Random.seed!(42)
+            rng = StableRNG(42)
             chain = sample(
-                model, sampler, 1000; discard_initial=1000, thinning=10, n_adapts=0
+                rng, model, sampler, 1000; discard_initial=1000, thinning=10, n_adapts=0
             )
             check_numerical(chain, [:m1, :m2], [-0.2, 0.6]; atol=0.1)
             check_logp_correct(sampler_inner)

test/mcmc/hmc.jl

@@ -131,7 +131,6 @@ using Turing
     # easily make it fail, despite many more samples than taken by most other tests. Hence
     # explicitly specifying the seeds here.
     @testset "hmcda+gibbs inference" begin
-        Random.seed!(12345)
         alg = Gibbs(:s => PG(20), :m => HMCDA(500, 0.8, 0.25; init_ϵ=0.05))
         res = sample(StableRNG(123), gdemo_default, alg, 3000; discard_initial=1000)
         check_gdemo(res)

test/mcmc/particle_mcmc.jl

@@ -55,8 +55,6 @@ using Turing
     end
 
     @testset "logevidence" begin
-        Random.seed!(100)
-
         @model function test()
             a ~ Normal(0, 1)
             x ~ Bernoulli(1)
@@ -67,7 +65,7 @@ using Turing
             return x
         end
 
-        chains_smc = sample(test(), SMC(), 100)
+        chains_smc = sample(StableRNG(100), test(), SMC(), 100)
 
         @test all(isone, chains_smc[:x])
         # For SMC, the chain stores the collective logevidence of the sampled trajectories

test/stdlib/RandomMeasures.jl

@@ -2,6 +2,7 @@ module RandomMeasuresTests
 
 using Distributions: Normal, sample
 using Random: Random
+using StableRNGs: StableRNG
 using Test: @test, @testset
 using Turing
 using Turing.RandomMeasures: ChineseRestaurantProcess, DirichletProcess
@@ -49,13 +50,12 @@ using Turing.RandomMeasures: ChineseRestaurantProcess, DirichletProcess
         end
 
         # Generate some test data.
-        Random.seed!(1)
-        data = vcat(randn(10), randn(10) .- 5, randn(10) .+ 10)
+        rng = StableRNG(1)
+        data = vcat(randn(rng, 10), randn(rng, 10) .- 5, randn(rng, 10) .+ 10)
         data .-= mean(data)
         data /= std(data)
 
         # MCMC sampling
-        Random.seed!(2)
         iterations = 500
         model_fun = infiniteGMM(data)
         chain = sample(model_fun, SMC(), iterations)