Блог

Я хочу генерировать случайные числа в Julia, используя многопоточность. Я использую Threads.@threads макрос для этого. Тем не менее, я изо всех сил пытаюсь зафиксировать количество семян, чтобы получать один и тот же результат каждый раз, когда я запускаю код. Вот мое испытание:

 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))
 

Приведенный выше сценарий дает разные результаты каждый раз, когда я его запускаю. Напротив, я получаю те же результаты, если использую простой цикл for:

 Random.seed!(12445)
b = []

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

sum(b)
 

У меня сложилось впечатление, что решение этого вопроса должно быть легким. И все же я не мог его найти. Мы очень ценим любую помощь.

Спасибо.

Вам нужно сгенерировать отдельный случайный поток для каждого потока. Самый простой способ-иметь генератор случайных чисел с другим начальным значением:

 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
 

Если вы не хотите рисковать корреляцией для разных семян случайных чисел, вы можете фактически обойти один генератор чисел:

 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))
 

Ciao Fabrizio. В BetaML я решил эту проблему с помощью:

 """
    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
 

Приведенная выше функция дает те же результаты, также независимо от количества потоков, используемых в Julia, и затем может быть использована, например, как здесь:

 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 amp;amp; a == A amp;amp; b == B


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

Предполагая, что вы находитесь на Julia 1.6, вы можете сделать, например, следующее:

 julia> using Random

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

Дело в том, что в настоящее время у Джулии уже есть отдельный генератор случайных чисел для каждого потока, поэтому вам не нужно создавать свой собственный (конечно, вы могли бы сделать это, как в других ответах, но вам это не нужно).

Также обратите внимание, что в будущих версиях Джулии:

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

часть не гарантирует получения воспроизводимых результатов. В Julia 1.6 Threads.@threads используется статическое планирование, но, как вы можете прочитать в его строке документа, оно может быть изменено.