Merge #281

281: pass RNG instead of seeds r=costachris a=costachris

Closes #260 

Pass RNG  objects instead of seeds to make randomness more traceable. 

Co-authored-by: costachris <[email protected]>

examples/AerosolActivation/aerosol_activation.jl

@@ -39,7 +39,9 @@ using Plots
 using Distributions
 using LinearAlgebra
 using Random
+
 rng_seed = 44
+rng = Random.seed!(Random.GLOBAL_RNG, rng_seed)
 
 using EnsembleKalmanProcesses
 using EnsembleKalmanProcesses.ParameterDistributions
@@ -169,7 +171,7 @@ nothing # hide
 N_ens = 50
 N_iter = 10
 
-initial_par = EKP.construct_initial_ensemble(priors, N_ens; rng_seed)
+initial_par = EKP.construct_initial_ensemble(rng, priors, N_ens)
 ekiobj = EKP.EnsembleKalmanProcess(initial_par, truth_sample, truth_array.obs_noise_cov, EKP.Inversion())
 nothing # hide
 

examples/Cloudy/Cloudy_example_eki.jl

@@ -27,6 +27,7 @@ using .DynamicalModel
 
 rng_seed = 41
 Random.seed!(rng_seed)
+rng = Random.seed!(Random.GLOBAL_RNG, rng_seed)
 
 homedir = pwd()
 figure_save_directory = homedir * "/output/"
@@ -121,7 +122,7 @@ truth_sample = truth.mean
 N_ens = 50 # number of ensemble members
 N_iter = 8 # number of EKI iterations
 # initial parameters: N_par x N_ens
-initial_par = construct_initial_ensemble(priors, N_ens; rng_seed)
+initial_par = construct_initial_ensemble(rng, priors, N_ens)
 ekiobj = EnsembleKalmanProcess(initial_par, truth_sample, truth.obs_noise_cov, Inversion(), Δt = 0.1)
 
 # Initialize a ParticleDistribution with dummy parameters. The parameters 

examples/Darcy/calibrate.jl

@@ -95,7 +95,7 @@ function main()
     N_ens = dofs + 2 # number of ensemble members
     N_iter = 10 # number of EKI iterations
     # initial parameters: N_params x N_ens
-    initial_params = construct_initial_ensemble(prior, N_ens; rng_seed = seed)
+    initial_params = construct_initial_ensemble(rng, prior, N_ens)
     ekiobj = EKP.EnsembleKalmanProcess(initial_params, truth_sample, obs_noise_cov, Inversion())
 
     println("Begin inversion")

examples/Lorenz/Lorenz_example.jl

@@ -17,7 +17,7 @@ using EnsembleKalmanProcesses.ParameterDistributions
 const EKP = EnsembleKalmanProcesses
 
 rng_seed = 4137
-Random.seed!(rng_seed)
+rng = Random.seed!(Random.GLOBAL_RNG, rng_seed)
 
 # Output figure save directory
 homedir = pwd()
@@ -192,7 +192,7 @@ lorenz_settings_G = lorenz_settings; # initialize to truth settings
 N_ens = 20 # number of ensemble members
 N_iter = 5 # number of EKI iterations
 # initial parameters: N_params x N_ens
-initial_params = construct_initial_ensemble(priors, N_ens; rng_seed = rng_seed)
+initial_params = construct_initial_ensemble(rng, priors, N_ens)
 
 ekiobj = EKP.EnsembleKalmanProcess(initial_params, truth_sample, truth.obs_noise_cov, Inversion())
 

examples/Lorenz/distributed_Lorenz_example.jl

@@ -41,7 +41,7 @@ using EnsembleKalmanProcesses.ParameterDistributions
 const EKP = EnsembleKalmanProcesses
 
 rng_seed = 4137
-Random.seed!(rng_seed)
+rng = Random.seed!(Random.GLOBAL_RNG, rng_seed)
 
 # Output figure save directory
 homedir = pwd()
@@ -216,7 +216,7 @@ lorenz_settings_G = lorenz_settings; # initialize to truth settings
 N_ens = 20 # number of ensemble members
 N_iter = 5 # number of EKI iterations
 # initial parameters: N_params x N_ens
-initial_params = construct_initial_ensemble(priors, N_ens; rng_seed = rng_seed)
+initial_params = construct_initial_ensemble(rng, priors, N_ens)
 
 ekiobj = EKP.EnsembleKalmanProcess(initial_params, truth_sample, truth.obs_noise_cov, Inversion())
 

examples/LossMinimization/loss_minimization.jl

@@ -22,7 +22,7 @@ nothing # hide
 
 # We set the seed for pseudo-random number generator for reproducibility.
 rng_seed = 41
-Random.seed!(rng_seed)
+rng = Random.seed!(Random.GLOBAL_RNG, rng_seed)
 nothing # hide
 
 # We set a stabilization level, which can aid the algorithm convergence
@@ -54,7 +54,7 @@ N_iterations = 10
 nothing # hide
 
 # The initial ensemble is constructed by sampling the prior
-initial_ensemble = EKP.construct_initial_ensemble(prior, N_ensemble; rng_seed = rng_seed)
+initial_ensemble = EKP.construct_initial_ensemble(rng, prior, N_ensemble)
 
 # We then initialize the Ensemble Kalman Process algorithm, with the initial ensemble, the
 # target, the stabilization and the process type (for EKI this is `Inversion`, initialized 
@@ -127,7 +127,7 @@ nothing # hide
 
 # We set the seed for pseudo-random number generator for reproducibility,
 rng_seed = 10
-Random.seed!(rng_seed)
+rng = Random.seed!(Random.GLOBAL_RNG, rng_seed)
 nothing # hide
 
 # A positive function can be minimized with a target of 0,
@@ -152,7 +152,7 @@ prior = combine_distributions([prior_u1, prior_u2])
 N_ensemble = 20
 N_iterations = 20
 
-initial_ensemble = EKP.construct_initial_ensemble(prior, N_ensemble; rng_seed = rng_seed)
+initial_ensemble = EKP.construct_initial_ensemble(rng, prior, N_ensemble)
 
 ensemble_kalman_process = EKP.EnsembleKalmanProcess(initial_ensemble, G_target, Γ_stabilization, Inversion())
 

examples/SparseLossMinimization/loss_minimization_sparse_eki.jl

@@ -20,7 +20,7 @@ nothing # hide
 
 # We set the seed for pseudo-random number generator for reproducibility.
 rng_seed = 41
-Random.seed!(rng_seed)
+rng = Random.seed!(Random.GLOBAL_RNG, rng_seed)
 nothing # hide
 
 # We set a stabilization level, which can aid the algorithm convergence
@@ -52,7 +52,7 @@ N_iterations = 10
 nothing # hide
 
 # The initial ensemble is constructed by sampling the prior
-initial_ensemble = EKP.construct_initial_ensemble(prior, N_ensemble; rng_seed = rng_seed)
+initial_ensemble = EKP.construct_initial_ensemble(rng, prior, N_ensemble)
 
 # Sparse EKI parameters
 γ = 1.0

examples/example_template.jl

@@ -44,7 +44,7 @@ priors = ParameterDistribution(prior_dict)
 ##
 
 # initial ensemble of parameters
-params = construct_initial_ensemble(priors, N_ens; rng_seed = rng_seed)
+params = construct_initial_ensemble(priors, N_ens)
 
 # constructor for the EKI object
 ekiobj =

src/EnsembleKalmanProcess.jl

@@ -372,30 +372,19 @@ end
     construct_initial_ensemble(
         rng::AbstractRNG,
         prior::ParameterDistribution,
-        N_ens::IT;
-        rng_seed::Union{IT, Nothing} = nothing,
+        N_ens::IT
     ) where {IT <: Int}
-    construct_initial_ensemble(prior::ParameterDistribution, N_ens::IT; kwargs...) where {IT <: Int}
+    construct_initial_ensemble(prior::ParameterDistribution, N_ens::IT) where {IT <: Int}
 
 Construct the initial parameters, by sampling `N_ens` samples from specified
 prior distribution. Returned with parameters as columns.
 """
-function construct_initial_ensemble(
-    rng::AbstractRNG,
-    prior::ParameterDistribution,
-    N_ens::IT;
-    rng_seed::Union{IT, Nothing} = nothing,
-) where {IT <: Int}
-    # Ensuring reproducibility of the sampled parameter values: 
-    # re-seed the rng *only* if we're given a seed
-    if rng_seed !== nothing
-        rng = Random.seed!(rng, rng_seed)
-    end
+function construct_initial_ensemble(rng::AbstractRNG, prior::ParameterDistribution, N_ens::IT) where {IT <: Int}
     return sample(rng, prior, N_ens) #of size [dim(param space) N_ens]
 end
 # first arg optional; defaults to GLOBAL_RNG (as in Random, StatsBase)
-construct_initial_ensemble(prior::ParameterDistribution, N_ens::IT; kwargs...) where {IT <: Int} =
-    construct_initial_ensemble(Random.GLOBAL_RNG, prior, N_ens; kwargs...)
+construct_initial_ensemble(prior::ParameterDistribution, N_ens::IT) where {IT <: Int} =
+    construct_initial_ensemble(Random.GLOBAL_RNG, prior, N_ens)
 
 """
     compute_error!(ekp::EnsembleKalmanProcess)

test/EnsembleKalmanProcess/runtests.jl

@@ -456,8 +456,5 @@ end
     u = ParameterDistribution(Dict("distribution" => d, "constraint" => no_constraint(), "name" => "test"))
     draw_1 = construct_initial_ensemble(rng, u, 1)
     draw_2 = construct_initial_ensemble(u, 1)
-    # Re-seeded draw should be the same as first draw
-    draw_3 = construct_initial_ensemble(rng, u, 1; rng_seed = rng_seed)
     @test !isapprox(draw_1, draw_2)
-    @test isapprox(draw_1, draw_3)
 end