There are several instances in the codebase where we have to create N pseudo-random number generators from one, for example:
|
""" |
|
update_keys!(pc::ParticleContainer) |
|
|
|
Create new unique keys for the particles in the ParticleContainer |
|
""" |
|
function update_keys!(pc::ParticleContainer, ref::Union{Particle,Nothing}=nothing) |
|
# Update keys to new particle ids |
|
nparticles = length(pc) |
|
n = ref === nothing ? nparticles : nparticles - 1 |
|
for i in 1:n |
|
pi = pc.vals[i] |
|
k = split(state(pi.rng.rng)) |
|
Random.seed!(pi.rng, k[1]) |
|
end |
|
return nothing |
|
end |
Right now, this is handled by split, which generates new seeds to be used with new PRNGs.
|
""" |
|
split(key::Integer, n::Integer=1) |
|
|
|
Split `key` into `n` new keys |
|
""" |
|
function split(key::Integer, n::Integer=1) |
|
T = typeof(key) |
|
inner_rng = Random.MersenneTwister(key) |
|
return rand(inner_rng, T, n) |
|
end |
This is not always a reliable method as the new seeds may end up accidentally generating highly correlated sequences of numbers. See, e.g., https://gee.cs.oswego.edu/dl/papers/oopsla14.pdf
We should therefore switch to an underlying RNG type that supports splitting. (Note that right now, we are using Philox2x [paper, library] but the library does not implement a split function, even though it should theoretically be possible. Compare e.g. with Numpy's implementation, which provides a jumped method for this purpose.)
There are several instances in the codebase where we have to create N pseudo-random number generators from one, for example:
AdvancedPS.jl/src/container.jl
Lines 121 to 136 in d0d180f
Right now, this is handled by
split, which generates new seeds to be used with new PRNGs.AdvancedPS.jl/src/rng.jl
Lines 33 to 42 in d0d180f
This is not always a reliable method as the new seeds may end up accidentally generating highly correlated sequences of numbers. See, e.g., https://gee.cs.oswego.edu/dl/papers/oopsla14.pdf
We should therefore switch to an underlying RNG type that supports splitting. (Note that right now, we are using Philox2x [paper, library] but the library does not implement a split function, even though it should theoretically be possible. Compare e.g. with Numpy's implementation, which provides a
jumpedmethod for this purpose.)