more methods working with a non-global RNG

[rfourquet]

This is #8380 rebased and on a branch in the main repo.
However I simplified it a bit: the rand() methods working with the global RNG were constrained by the types with which rand(::MersenneTwister, ...) actually worked. But I think this did not add much value (the "no method" error is simply forwarded), and could easily fall out of sync, so now there is a fall-back rand(T::Type) = rand(GLOBAL_RNG, T).

I need to review this once more, hence the "RFC".

[ViralBShah]

Thanks @rfourquet. How close are we to the original performance with this PR merged?

[ViralBShah]

In julia/deps/Rmath/src/runif.c, the RNG was synced with the one in Julia. In that, Rmath would end up calling the same global RNG as in Julia. That may be tough to achieve now. The normal RNG is anyways already different given that we have a native implementation.

Perhaps we should at this point give up on trying to sync the Rmath RNGs with Julia, and just focus on reducing/removing dependence on Rmath in Distributions.jl.

Cc: @simonbyrne @johnmyleswhite @dmbates

[ViralBShah]

Cc: @JuliaBackports

[ViralBShah]

Single rand() is now 2x slower, whereas the array version is now 2x faster on my mac. I wonder what tricks we can use to improve the basic rand(). The code looks comparable to the C version in dSFMT, with perhaps GLOBAL_RNG not being passed around to all the routines. However, that may be getting inlined anyways.

Julia 0.3.2:

julia> @time for i=1:10^8; rand(); end
elapsed time: 0.315577166 seconds (0 bytes allocated)

julia> @time rand(10^8);
elapsed time: 0.734464713 seconds (800000168 bytes allocated)

master:

julia> @time for i=1:10^8; rand(); end
elapsed time: 0.734230661 seconds (0 bytes allocated)

julia> @time rand(10^8);
elapsed time: 0.408209218 seconds (800000168 bytes allocated)

@vtjnash Is there a way to check that rand() is getting properly inlined?

[vtjnash]

Can you repost with the scalar case being inside a function?

edit: scratch that, i see it doesn't affect timing

[johnmyleswhite]

My understanding from @simonbyrne was that we might not need R's RNG's for Distributions, only deterministic functions found in libRmath.

[vtjnash]

FWIW, the timing on my machine before and after this merged appears to be nearly identical for the scalar case, with this PR maybe being slightly faster (I didn't collect enough data to be statistically significant however)

[rfourquet]

I was expecting this PR to improve array versions for non Float64 types with the global RNG, but it seems that global and non global versions have in fact the same performances (so the "globalness" of GLOBAL_RNG have no speed effects on array versions). So I think the speed differences only come from the fact that the global RNG handled by libdSFMT was faster than a MersenneTwister instance (even a non-global one).
For the scalar version of rand() with the global RNG, I have no idea what to do (maybe a word in the docs to warn the user). For the non-Float64 array versions, I plan to write optimized algotithm.

[simonbyrne]

Yes, this seems to have broken Distributions:
JuliaStats/Distributions.jl#293

[lindahua]

Random number generation in Distributions.jl is seriously broken. It seems that in addition to JuliaStats/Distributions.jl#293, Binomial distribution also produces wrong results (on Julia 0.4 latest master):

julia> using Distributions

julia> p = Binomial(1, 0.5)
Binomial(size=1, prob=0.5)

julia> x = rand(p, 10^6);

julia> all(x .== 0)  # if it works well, half of x should be 1, now is produces all zeros!
true

[lindahua]

This should not be backported to 0.3 until we make sure that it works with Distributions.jl.

[staticfloat]

@lindahua We're closing the backport window on 0.3.3 soon. Unless you say otherwise, I'm assuming this is not a candidate for backporting at this time.

[ViralBShah]

Taking off the backport list, since a lot of work has happened on RNGs and much of it is API changes and still stabilizing.