# Random number generation within Compile and Parallelization->True

According to this thread RandomChoice is a compilable function. I use it to generate random numbers in a compiled function with parallelization enabled.

I wrote another variant of my program where I generate these random numbers outside of the compiled function and then use the generated array as an input argument of the compiled function.

Now it seems that when computing moments of my random trajectories with the parallelized compiled function (where the random numbers are generated within Compile), the outcome of the averaged result for a fixed number of steps is always the same as if it always uses the same sequence of random numbers.

On the other hand, if I generate the random numbers outside of Compile the averaged result fluctuates when repeating the calculation for a fixed number of steps as expected.

Anybody has a clue of what could be going on in the parallelized version of the program? Is it maybe due to the fact that if I generate the random numbers via Compile in parallel, the sequence of random numbers per thread is the same on all threads? And if yes, how would I avoid this but still use the parallelization capabilities of Compile?

I am aware of this tutorial in the documentation, but cannot apply it. What I want is not to predefine a seed, but genuinely independent random number generation on each thread when calling the parallelized compiled function.

[Edit2]

Ok, now I found the culprit. Probably it is a trivial error which I don't see?

If I execute the following in a single code block

x0 = Table[{RandomVariate[pdf1], RandomVariate[pdf2]}, {nTraj}];
list = av[cf[x0, nSteps], rangeav];


then list remains almost the same after several runs of the code block.

If I execute the commands in two code blocks, i.e.

x0 = Table[{RandomVariate[pdf1], RandomVariate[pdf2]}, {nTraj}];


and subsequently

list = av[cf[x0, nSteps], rangeav];


then it behaves as supposed and changes drastically after each run of the code blocks. I am confused now whether this is normal behavior or not?

So sometimes it seems to take over the new argument x0 (which is updated correctly) and sometimes it seems to be stuck and then it uses the same random numbers within Compile. The random number generation, however, does not depend on x0 (x0 is only used to pre-allocate an array x).

Or in other words, it cannot update the argument x0 and create a different set of random numbers within Compile at the same time.

[Edit3]

As outlined in the comments, I think this behavior is a bug. To make the code lines now completely unrelated, I define a symbolic continous pdf by

pdf = ProbabilityDistribution[
Sqrt[2]/Pi/(1 + x^4), {x, -Infinity, Infinity}];


The specific appearence of the pdf doesn't matter. Now if I call

RandomVariate[pdf];
list = av[cf[x0, nSteps], rangeav];


in the same code block and have Parallelization->True in the compiled function cf then only the sequence of results for the different threads changes for different runs, so that the position of the sublists are interchanged but the overall list gives the same result for different runs.

However, if I call

RandomVariate[NormalDistribution[0.,1.]];
list = av[cf[x0, nSteps], rangeav];


Then the sublists in list are created anew every time I run the code block.

So it seems as if RandomVariate[pdf] of a pre-defined pdf triggers the seed of a compiled function (which itself generates random numbers) if it is called in the same code block. If one calls RandomVariate[pdf] in a separate code block everything is fine again.

Very, very strange.

Now the question remains whether this kind of behaviour is reproducable in a minimal working example or appears only in the specific case of my compiled function that I use.

[Edit4]

Here is a minimal working example that produces the same issue on my machine with MMA 11.3

(*Define pdf*)
pdf = ProbabilityDistribution[
Sqrt[2]/Pi/(1 + x^4), {x, -Infinity, Infinity}];

(*Define compiled funcion*)
cf := Compile[{{nSteps, _Integer, 0}},
Block[{rand},
(*Random number generation*)
rand = RandomChoice[{1./6, 1./6, 2./3} -> {Sqrt[3.], -Sqrt[3.],
0.}, {nSteps}];
rand], CompilationTarget -> "C", RuntimeAttributes -> {Listable},
Parallelization -> True, RuntimeOptions -> "Speed"];


Now investigate the two cases (execute the two lines in one code block respectively):

list=Table[RandomVariate[pdf];
cf[Table[10, {$ProcessorCount}]], {2}];  gives output {{{1.73205, -1.73205, 1.73205, 1.73205, 0., 1.73205, 0., 0., 0., 1.73205}, {0., -1.73205, 1.73205, 0., 0., -1.73205, 0., 0., 0., 0.}, {0., 0., -1.73205, 0., 0., 0., 0., 0., 0., 0.}, {0., 0., 0., 1.73205, 0., 0., 0., -1.73205, 0., 0.}}, {{1.73205, -1.73205, 1.73205, 1.73205, 0., 1.73205, 0., 0., 0., 1.73205}, {0., 0., 0., 1.73205, 0., 0., 0., -1.73205, 0., 0.}, {0., 0., -1.73205, 0., 0., 0., 0., 0., 0., 0.}, {0., -1.73205, 1.73205, 0., 0., -1.73205, 0., 0., 0., 0.}}}  and list=Table[RandomVariate[NormalDistribution[0., 1.]]; cf[Table[10, {$ProcessorCount}]], {2}]


gives output

{{{0., 0., -1.73205, 0., 0., 0., 0., 0., 0., 0.}, {1.73205, 0.,
1.73205, 1.73205, 0., 0., 0., 0., 0., 0.}, {0., -1.73205, -1.73205,
0., 0., -1.73205, 0., 0., 1.73205, 0.}, {1.73205, -1.73205, 0.,
0., 0., 0., -1.73205, 0., 0., 0.}}, {{0., -1.73205, 0., -1.73205,
0., 1.73205, 1.73205, -1.73205, 0., 1.73205}, {0., 0., 0., 0., 0.,
0., 0., 0., 0., 0.}, {0., 0., 0., 0., 0., 0., 1.73205, 1.73205, 0.,
0.}, {1.73205, 0., 0., 0., 0., 0., 0., 0., 0., 0.}}}


Now compare list[[1]] with list[[2]] for the two code blocks. With the second code block the random numbers are generated anew each time the command is executed. With the first code block only the sublists are interchanged, but the total list[[1]] and list[[2]] remains the same.

Or in shorter and possibly more convenient form:

(*Function for comparison of parts of list*)
same[{a_, b__}] := SameQ[a, b]


and then

m = 10;
list = Table[RandomVariate[pdf];
cf[Table[10, {$ProcessorCount}]], {m}]; same[Map[Sort, list]]  which gives True and m = 10; list = Table[RandomVariate[NormalDistribution[0., 1.]]; cf[Table[10, {$ProcessorCount}]], {m}];
same[Map[Sort, list]]


which gives False.

• My suggestion would be to write your RNG-using Compile function to take a second argument for a seed or something like that, then pass in a set of random seeds. On the other hand without some code to work with it'll be hard to debug fully for you. – b3m2a1 Nov 10 '18 at 22:31
• Thanks for the suggestion. It is sort of strange because I can just let MMA give the random numbers as an output and they differ for the different threads. I wonder what leads then to this strange behavior. Mhh – Display Name Nov 10 '18 at 22:38
• @b3m2a1: See my [Edit2] it has to do with some outside definitions which are then passed over to the compiled function as an argument. – Display Name Nov 10 '18 at 23:08
• I think we'll need a more complete sense of your code to understand what's going on. Can you make a minimum working example with all functions defined? (e.g. af, pdf1, etc.) – b3m2a1 Nov 10 '18 at 23:30
• The problem is that the whole computation is a little bit long and consists of different parts. I will try to come up with a simple example which is shorter tomorrow. – Display Name Nov 10 '18 at 23:31