CODEWITHSUNDEEP
randomjuliamontecarlo
dapias
dapias

Reputation: 2632

Parallel random numbers julia

Consider the basic iteration to generate N random numbers and save them in an array (assume either that we are not interested in array comprehensions and also that we don't know the calling rand(N))

function random_numbers(N::Int)
array = zeros(N)
for i in 1:N
    array[i] = rand()
end
array
end

I am interested in a similar function that takes advantage of the cores of my laptop to generate the same array. I have checked this nice blog where the macros @everywhere, @spawn and @parallel are introduced but there the calculation is carried out "on-the-fly" and an array is not needed to save the data.

I have the impression that this is very basic and can be done easily using perhaps the function pmap but I am unfamiliar with parallel computing.

My aim is to apply this method to a function that I have built to generate random numbers drawn from an unusual distribution.

Upvotes: 2

Views: 787

Answers (2)

Bogumił Kamiński
Bogumił Kamiński

Reputation: 69949

I would recommend to do a more careful initialization of random number generators in parallel processes, e.g:

# choose the seed you want
@everywhere srand(1)
# replace 10 below by maximum process id in your case
@everywhere const LOCAL_R = randjump(Base.GLOBAL_RNG, 10)[myid()]
# here is an example usage
@everywhere f() = rand(LOCAL_R)

In this way you:

  • make sure that your results are reproducible;
  • have control that there is no overlap between random sequences generated by different processes.

Upvotes: 3

Dan Getz
Dan Getz

Reputation: 18217

As suggested in the comment more clarification in the question is always welcome. However, it seems pmap will do what is required. The relevant documentation is here.

The following is a an example. Note, the time spent in the pmap method is half of the regular map. With 16 cores, the situation might be substantially better:

julia> addprocs(2)
2-element Array{Int64,1}:
 2
 3

julia> @everywhere long_rand() = foldl(+,0,(randn() for i=1:10_000_000))

julia> long_rand()
-1165.9596619177153

julia> @time map(x->long_rand(), zeros(10,10))
  8.455930 seconds (204.89 k allocations: 11.069 MiB)
10×10 Array{Float64,2}:
   ⋮
   ⋮ 

julia> @time pmap(x->long_rand(), zeros(10,10));
  6.125479 seconds (773.08 k allocations: 42.242 MiB, 0.25% gc time)

julia> @time pmap(x->long_rand(), zeros(10,10))
  4.609745 seconds (20.99 k allocations: 954.991 KiB)
10×10 Array{Float64,2}:
   ⋮
   ⋮

Upvotes: 2

Related Questions