Re `seed!` ing the RNG at each step of a simulation?

I am running long stochastic simulations that look like:

sim_name = "example/1"
Random.seed!(sim_name)
ctx = setup(sim_name)
for frame in 1:1000
  ctx = loop(ctx)
end

I need to be able to restart simulations that crash (for example, due to hardware issues). I have defined serialization functions for ctx, but from what I can tell, there is no stable way to serialize the state of the default RNG because its implementation can change between minor Julia versions.

Is reseeding the default RNG on each step a valid alternative?

sim_name = "example/1"
Random.seed!("$(sim_name)/$(0)")
ctx = setup(sim_name)
for frame in 1:1000
  Random.seed!("$(sim_name)/$(frame)")
  ctx = loop(ctx)
end

No, since this is also not reproducible across Julia versions (and may be less random).

The best approach depends a bit on what you are doing. How is randomness used in your simulation?

e.g. you could use StableRNGs.jl if you restrict yourself to random sampling functions in its stable API.

3 Likes

I’m okay with it only being reproducible in a single Julia version, but what do you mean by it being less random?

I am simulating a jump process; the RNG is used to randomly select the next event, and when that event happens, the event callback then also uses the RNG to determine what happened or reject the event.

I could store the RNG in ctx and then rewrite all the event callbacks to use the RNG in ctx instead of the default one.

Then what’s the problem with serializing the RNG state?

If you re-seed too often it changes the statistics. (Imagine the extreme case where you re-seed after each random number.)

2 Likes

Yes, that works:

julia> using Base64, Serialization, Random

julia> function rng2str()
         base64encode(sprint(serialize, copy(Random.default_rng())))
       end
rng2str (generic function with 1 method)

julia> function str2rng(s)
         copy!(Random.default_rng(), deserialize(IOBuffer(base64decode(s))))
         nothing
       end
str2rng (generic function with 1 method)

julia> s = rng2str()
"N0pMGgQAAAA1EAEHWG9zaGlybx8LXP2hgFmEEZoCX8mihII/mgEGUmFuZG9tRAlpUUqIGas5+QnpupKUwWTyqAms62WLbMY8vAngTB2njJ+GtwmHsEKhRY8d/g=="

julia> rand()
0.29811661892484653

julia> str2rng(s)

julia> rand()
0.29811661892484653

Though using deserialize doesn’t seem ideal because its docstring says “Malformed data can result in process termination.”