CODEWITHSUNDEEP
multithreadingfor-looprandomjulia
Fabrizio Leone
Fabrizio Leone

Reputation: 132

Setting seeds in multi-threading loop in Julia

I want to generate random numbers in Julia using multi-threading. I am using the Threads.@threads macro to accomplish it. However, I struggle fixing the number of seeds to obtain the same result every time I run the code. Here is my trial:

Random.seed!(1234)
a = [Float64[] for _ in 1:10]

Threads.@threads for i = 1:10
    push!(a[Threads.threadid()],rand())
end

sum(reduce(vcat, a))

The script above delivers different results every time I run it. By contrast, I get the same results if I use a plain for loop:

Random.seed!(12445)
b = []

for i = 1:10
    push!(b,rand())
end

sum(b)

I have the impression that the solution to this issue must be easy. Still, I couldn't find it. Any help is much appreciated.

Thank you.

Upvotes: 5

Views: 889

Answers (3)

Bogumił Kamiński
Bogumił Kamiński

Reputation: 69949

Assuming you are on Julia 1.6 you can do e.g. the following:

julia> using Random

julia> foreach(i -> Random.seed!(Random.default_rng(i), i), 1:Threads.nthreads())

The point is that currently Julia already has a separate random number generator per thread so you do not need to generate your own (of course you could do it as in the other answers, but you do not have to).

Also note that in the future versions of Julia the:

Threads.@threads for i = 1:10
    push!(a[Threads.threadid()],rand())
end

part is not guaranteed to produce reproducible results. In Julia 1.6 Threads.@threads uses static scheduling, but as you can read in its docstring it is subject to change.

Upvotes: 2

Antonello
Antonello

Reputation: 6423

Ciao Fabrizio. In BetaML I solved this problem with:

"""
    generateParallelRngs(rng::AbstractRNG, n::Integer;reSeed=false)
For multi-threaded models, return n independent random number generators (one per thread) to be used in threaded computations.
Note that each ring is a _copy_ of the original random ring. This means that code that _use_ these RNGs will not change the original RNG state.
Use it with `rngs = generateParallelRngs(rng,Threads.nthreads())` to have a separate rng per thread.
By default the function doesn't re-seed the RNG, as you may want to have a loop index based re-seeding strategy rather than a threadid-based one (to guarantee the same result independently of the number of threads).
If you prefer, you can instead re-seed the RNG here (using the parameter `reSeed=true`), such that each thread has a different seed. Be aware however that the stream  of number generated will depend from the number of threads at run time.
"""
function generateParallelRngs(rng::AbstractRNG, n::Integer;reSeed=false)
    if reSeed
        seeds = [rand(rng,100:18446744073709551615) for i in 1:n] # some RNGs have issues with too small seed
        rngs  = [deepcopy(rng) for i in 1:n]
        return Random.seed!.(rngs,seeds)
    else
        return [deepcopy(rng) for i in 1:n]
    end
end

The function above deliver the same results also independently of the number of threads used in Julia and can then be used for example like here:

using Test

TESTRNG = MersenneTwister(123)

println("** Testing generateParallelRngs()...")
x = rand(copy(TESTRNG),100)

function innerFunction(bootstrappedx; rng=Random.GLOBAL_RNG)
     sum(bootstrappedx .* rand(rng) ./ 0.5)
end
function outerFunction(x;rng = Random.GLOBAL_RNG)
    masterSeed = rand(rng,100:9999999999999) # important: with some RNG it is important to do this before the generateParallelRngs to guarantee independance from number of threads
    rngs       = generateParallelRngs(rng,Threads.nthreads()) # make new copy instances
    results    = Array{Float64,1}(undef,30)
    Threads.@threads for i in 1:30
        tsrng = rngs[Threads.threadid()]    # Thread safe random number generator: one RNG per thread
        Random.seed!(tsrng,masterSeed+i*10) # But the seeding depends on the i of the loop not the thread: we get same results indipendently of the number of threads
        toSample = rand(tsrng, 1:100,100)
        bootstrappedx = x[toSample]
        innerResult = innerFunction(bootstrappedx, rng=tsrng)
        results[i] = innerResult
    end
    overallResult = mean(results)
    return overallResult
end


# Different sequences..
@test outerFunction(x) != outerFunction(x)

# Different values, but same sequence
mainRng = copy(TESTRNG)
a = outerFunction(x, rng=mainRng)
b = outerFunction(x, rng=mainRng)

mainRng = copy(TESTRNG)
A = outerFunction(x, rng=mainRng)
B = outerFunction(x, rng=mainRng)

@test a != b && a == A && b == B


# Same value at each call
a = outerFunction(x,rng=copy(TESTRNG))
b = outerFunction(x,rng=copy(TESTRNG))
@test a == b

Upvotes: 2

Przemyslaw Szufel
Przemyslaw Szufel

Reputation: 42234

You need to generate a separate random stream for each thread. The simplest way is to have a random number generator with a different seed:

using Random

rngs = [MersenneTwister(i) for i in 1: Threads.nthreads()];

Threads.@threads for i = 1:10
     val = rand(rngs[Threads.threadid()])
     # do something with val
end

If you do not want to risk correlation for different random number seeds you could actually jump around a single number generator:

julia> rngs2 = Future.randjump.(Ref(MersenneTwister(0)), big(10)^20 .* (1:Threads.nthreads()))
4-element Vector{MersenneTwister}:
 MersenneTwister(0, (200000000000000000000, 0))
 MersenneTwister(0, (400000000000000000000, 0))
 MersenneTwister(0, (600000000000000000000, 0))
 MersenneTwister(0, (800000000000000000000, 0))

Upvotes: 3

Related Questions