Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

more methods working with a non-global RNG #8854

Merged
merged 3 commits into from
Nov 1, 2014
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions base/bitarray.jl
Original file line number Diff line number Diff line change
Expand Up @@ -331,10 +331,10 @@ bitpack{T,N}(A::AbstractArray{T,N}) = convert(BitArray{N}, A)

## Random ##

function bitarray_rand_fill!(B::BitArray)
function bitarray_rand_fill!(rng, B::BitArray) # rng is an AbstractRNG
length(B) == 0 && return B
Bc = B.chunks
rand!(Bc)
rand!(rng, Bc)
Bc[end] &= @_msk_end length(B)
return B
end
Expand Down
144 changes: 72 additions & 72 deletions base/random.jl
Original file line number Diff line number Diff line change
Expand Up @@ -13,17 +13,13 @@ abstract AbstractRNG

type MersenneTwister <: AbstractRNG
state::DSFMT_state
seed::Union(Uint32,Vector{Uint32})
vals::Vector{Float64}
idx::Int
seed::Vector{Uint32}

function MersenneTwister(seed::Vector{Uint32})
state = DSFMT_state()
dsfmt_init_by_array(state, seed)
return new(state, seed, Array(Float64, dsfmt_get_min_array_size()), dsfmt_get_min_array_size())
end

MersenneTwister(seed=0) = MersenneTwister(make_seed(seed))
MersenneTwister(seed) = srand(new(DSFMT_state(), Array(Float64, dsfmt_get_min_array_size())),
seed)
MersenneTwister() = MersenneTwister(0)
end

## Low level API for MersenneTwister
Expand All @@ -47,50 +43,43 @@ end
@inline rand_ui32(r::MersenneTwister) = reinterpret(Uint64, rand_close1_open2(r)) % Uint32


function srand(r::MersenneTwister, seed)
function srand(r::MersenneTwister, seed::Vector{Uint32})
r.seed = seed
dsfmt_init_gen_rand(r.state, seed)
dsfmt_init_by_array(r.state, r.seed)
r.idx = length(r.vals)
return r
end

## initialization

function srand()
__init__() = srand()

## make_seed()
# make_seed methods produce values of type Array{Uint32}, suitable for MersenneTwister seeding

function make_seed()

@unix_only begin
try
srand("/dev/urandom")
return make_seed("/dev/urandom", 4)
catch
println(STDERR, "Entropy pool not available to seed RNG; using ad-hoc entropy sources.")
seed = reinterpret(Uint64, time())
seed = hash(seed, uint64(getpid()))
try
seed = hash(seed, parseint(Uint64, readall(`ifconfig` |> `sha1sum`)[1:40], 16))
end
srand(seed)
return make_seed(seed)
end
end

@windows_only begin
a = zeros(Uint32, 2)
win32_SystemFunction036!(a)
srand(a)
return a
end
end

__init__() = srand()

## srand()

function srand(seed::Vector{Uint32})
GLOBAL_RNG.seed = seed
dsfmt_init_by_array(GLOBAL_RNG.state, seed)
GLOBAL_RNG.idx = length(GLOBAL_RNG.vals)
return GLOBAL_RNG
end
srand(n::Integer) = srand(make_seed(n))

function make_seed(n::Integer)
n < 0 && throw(DomainError())
seed = Uint32[]
Expand All @@ -103,63 +92,76 @@ function make_seed(n::Integer)
end
end

function srand(filename::String, n::Integer)
function make_seed(filename::String, n::Integer)
open(filename) do io
a = Array(Uint32, int(n))
read!(io, a)
srand(a)
a
end
end
srand(filename::String) = srand(filename, 4)

## srand()

srand(r::MersenneTwister) = srand(r, make_seed())
srand(r::MersenneTwister, n::Integer) = srand(r, make_seed(n))
srand(r::MersenneTwister, filename::String, n::Integer=4) = srand(r, make_seed(filename, n))

srand() = srand(GLOBAL_RNG)
srand(seed::Union(Integer, Vector{Uint32})) = srand(GLOBAL_RNG, seed)
srand(filename::String, n::Integer=4) = srand(GLOBAL_RNG, filename, n)

## Global RNG

const GLOBAL_RNG = MersenneTwister()
globalRNG() = GLOBAL_RNG

## random floating point values
# rand: a non-specified RNG defaults to GLOBAL_RNG

rand(r::MersenneTwister=GLOBAL_RNG) = rand_close_open(r)
rand() = rand(GLOBAL_RNG)
rand(T::Type) = rand(GLOBAL_RNG, T)
rand(::()) = rand(GLOBAL_RNG, ()) # needed to resolve ambiguity
rand(dims::Dims) = rand(GLOBAL_RNG, dims)
rand(dims::Int...) = rand(dims)
rand(T::Type, dims::Dims) = rand(GLOBAL_RNG, T, dims)
rand(T::Type, d1::Int, dims::Int...) = rand(T, tuple(d1, dims...))
rand!(A::AbstractArray) = rand!(GLOBAL_RNG, A)

rand(::Type{Float64}) = rand()
## random floating point values

rand(::Type{Float32}) = float32(rand())
rand(::Type{Float16}) = float16(rand())
rand(r::AbstractRNG) = rand(r, Float64)

rand{T<:Real}(::Type{Complex{T}}) = complex(rand(T),rand(T))
# MersenneTwister
rand(r::MersenneTwister, ::Type{Float64}) = rand_close_open(r)
rand{T<:Union(Float16, Float32)}(r::MersenneTwister, ::Type{T}) = convert(T, rand(r, Float64))

## random integers
## random integers (MersenneTwister)

rand(::Type{Uint8}) = rand(Uint32) % Uint8
rand(::Type{Uint16}) = rand(Uint32) % Uint16
rand(::Type{Uint32}) = rand_ui32(GLOBAL_RNG)
rand(::Type{Uint64}) = uint64(rand(Uint32)) <<32 | rand(Uint32)
rand(::Type{Uint128}) = uint128(rand(Uint64))<<64 | rand(Uint64)
rand(r::MersenneTwister, ::Type{Uint8}) = rand(r, Uint32) % Uint8
rand(r::MersenneTwister, ::Type{Uint16}) = rand(r, Uint32) % Uint16
rand(r::MersenneTwister, ::Type{Uint32}) = rand_ui32(r)
rand(r::MersenneTwister, ::Type{Uint64}) = uint64(rand(r, Uint32)) <<32 | rand(r, Uint32)
rand(r::MersenneTwister, ::Type{Uint128}) = uint128(rand(r, Uint64))<<64 | rand(r, Uint64)

rand(::Type{Int8}) = rand(Uint32) % Int8
rand(::Type{Int16}) = rand(Uint32) % Int16
rand(::Type{Int32}) = reinterpret(Int32,rand(Uint32))
rand(::Type{Int64}) = reinterpret(Int64,rand(Uint64))
rand(::Type{Int128}) = reinterpret(Int128,rand(Uint128))
rand(r::MersenneTwister, ::Type{Int8}) = rand(r, Uint32) % Int8
rand(r::MersenneTwister, ::Type{Int16}) = rand(r, Uint32) % Int16
rand(r::MersenneTwister, ::Type{Int32}) = reinterpret(Int32, rand(r, Uint32))
rand(r::MersenneTwister, ::Type{Int64}) = reinterpret(Int64, rand(r, Uint64))
rand(r::MersenneTwister, ::Type{Int128}) = reinterpret(Int128, rand(r, Uint128))

# Arrays of random numbers
## random complex values

rand(::Type{Float64}, dims::Dims) = rand!(Array(Float64, dims))
rand(::Type{Float64}, dims::Int...) = rand(Float64, dims)
rand{T<:Real}(r::AbstractRNG, ::Type{Complex{T}}) = complex(rand(r, T), rand(r, T))

rand(dims::Dims) = rand(Float64, dims)
rand(dims::Int...) = rand(Float64, dims)
## Arrays of random numbers

rand(r::AbstractRNG) = rand(r, Float64)
rand(r::AbstractRNG, dims::Dims) = rand!(r, Array(Float64, dims))
rand(r::AbstractRNG, dims::Dims) = rand(r, Float64, dims)
rand(r::AbstractRNG, dims::Int...) = rand(r, dims)

function rand!{T}(A::Array{T})
for i = 1:length(A)
A[i] = rand(T)
end
A
end
rand(r::AbstractRNG, T::Type, dims::Dims) = rand!(r, Array(T, dims))
rand(r::AbstractRNG, T::Type, d1::Int, dims::Int...) = rand(r, T, tuple(d1, dims...))
# note: the above method would trigger an ambiguity warning if d1 was not separated out:
# rand(r, ()) would match both this method and rand(r, dims::Dims)
# moreover, a call like rand(r, NotImplementedType()) would be an infinite loop

function rand!{T}(r::AbstractRNG, A::AbstractArray{T})
for i = 1:length(A)
Expand All @@ -168,6 +170,8 @@ function rand!{T}(r::AbstractRNG, A::AbstractArray{T})
A
end

# MersenneTwister

function rand_AbstractArray_Float64!(r::MersenneTwister, A::AbstractArray{Float64})
n = length(A)
# what follows is equivalent to this simple loop but more efficient:
Expand Down Expand Up @@ -203,14 +207,6 @@ function rand!(r::MersenneTwister, A::Array{Float64})
A
end

rand!(A::AbstractArray{Float64}) = rand!(GLOBAL_RNG, A)
rand!(A::Array{Float64}) = rand!(GLOBAL_RNG, A)

rand(T::Type, dims::Dims) = rand!(Array(T, dims))
rand{T<:Number}(::Type{T}) = error("no random number generator for type $T; try a more specific type")
rand{T<:Number}(::Type{T}, dims::Int...) = rand(T, dims)


## Generate random integer within a range

# remainder function according to Knuth, where rem_knuth(a, 0) = a
Expand Down Expand Up @@ -297,16 +293,20 @@ end
rand{T}(r::Range{T}, dims::Dims) = rand!(r, Array(T, dims))
rand(r::Range, dims::Int...) = rand(r, dims)

## random BitArrays (AbstractRNG)

## random Bools
rand!(r::AbstractRNG, B::BitArray) = Base.bitarray_rand_fill!(r, B)

rand!(B::BitArray) = Base.bitarray_rand_fill!(B)
randbool(r::AbstractRNG, dims::Dims) = rand!(r, BitArray(dims))
randbool(r::AbstractRNG, dims::Int...) = rand!(r, BitArray(dims))

randbool(dims::Dims) = rand!(BitArray(dims))
randbool(dims::Dims) = rand!(BitArray(dims))
randbool(dims::Int...) = rand!(BitArray(dims))

randbool() = ((rand(Uint32) & 1) == 1)
rand(::Type{Bool}) = randbool()
randbool(r::MersenneTwister=GLOBAL_RNG) = ((rand(r, Uint32) & 1) == 1)

rand(r::MersenneTwister, ::Type{Bool}) = randbool(r)


## randn() - Normally distributed random numbers using Ziggurat algorithm

Expand Down
30 changes: 11 additions & 19 deletions doc/stdlib/base.rst
Original file line number Diff line number Diff line change
Expand Up @@ -3931,43 +3931,35 @@ Random Numbers

Random number generation in Julia uses the `Mersenne Twister library <http:https://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/#dSFMT>`_. Julia has a global RNG, which is used by default. Multiple RNGs can be plugged in using the ``AbstractRNG`` object, which can then be used to have multiple streams of random numbers. Currently, only ``MersenneTwister`` is supported.

Most functions related to random generation accept an optional ``AbstractRNG`` as the first argument, ``rng`` , which defaults to the global one if not provided. Morever, some of them accept optionally dimension specifications ``dims...`` (which can be given as a tuple) to generate arrays of random values.

A ``MersenneTwister`` RNG can generate random numbers of the following types: ``Float16, Float32, Float64, Bool, Int16, Uint16, Int32, Uint32, Int64, Uint64, Int128, Uint128`` (or complex numbers or arrays of those types). Random floating point numbers are generated uniformly in [0,1).

.. function:: srand([rng], [seed])

Reseed the random number generator. If a ``seed`` is provided, the RNG will give a reproducible sequence of numbers, otherwise Julia will get entropy from the system. The ``seed`` may be an unsigned integer, a vector of unsigned integers or a filename, in which case the seed is read from a file. If the argument ``rng`` is not provided, the default global RNG is seeded.
Reseed the random number generator. If a ``seed`` is provided, the RNG will give a reproducible sequence of numbers, otherwise Julia will get entropy from the system. The ``seed`` may be a non-negative integer, a vector of ``Uint32`` integers or a filename, in which case the seed is read from a file.

.. function:: MersenneTwister([seed])

Create a ``MersenneTwister`` RNG object. Different RNG objects can have their own seeds, which may be useful for generating different streams of random numbers.

.. function:: rand() -> Float64
.. function:: rand([rng], [t::Type], [dims...])

Generate a ``Float64`` random number uniformly in [0,1)
Generate a random value or an array of random values of the given type, ``t``, which defaults to ``Float64``.

.. function:: rand!([rng], A)

Populate the array A with random number generated from the specified RNG.

.. function:: rand(rng::AbstractRNG, [dims...])

Generate a random ``Float64`` number or array of the size specified by dims, using the specified RNG object. Currently, ``MersenneTwister`` is the only available Random Number Generator (RNG), which may be seeded using srand.

.. function:: rand(dims or [dims...])

Generate a random ``Float64`` array of the size specified by dims

.. function:: rand(t::Type, [dims...])

Generate a random number or array of random numbes of the given type.
Populate the array A with random values.

.. function:: rand(r, [dims...])

Pick a random element or array of random elements from range ``r`` (for example, ``1:n`` or ``0:2:10``).

.. function:: randbool([dims...])
.. function:: randbool([rng], [dims...])

Generate a random boolean value. Optionally, generate an array of random boolean values.
Generate a random boolean value. Optionally, generate a ``BitArray`` of random boolean values.

.. function:: randn([rng], dims or [dims...])
.. function:: randn([rng], [dims...])

Generate a normally-distributed random number with mean 0 and standard deviation 1. Optionally generate an array of normally-distributed random numbers.

Expand Down
7 changes: 6 additions & 1 deletion test/random.jl
Original file line number Diff line number Diff line change
Expand Up @@ -19,12 +19,17 @@ srand(0); rand(); x = rand(384);
# Try a seed larger than 2^32
@test rand(MersenneTwister(5294967296)) == 0.3498809918210497

# Test array filling, Issue #7643
# Test array filling, Issues #7643, #8360
@test rand(MersenneTwister(0), 1) == [0.8236475079774124]
A = zeros(2, 2)
rand!(MersenneTwister(0), A)
@test A == [0.8236475079774124 0.16456579813368521;
0.9103565379264364 0.17732884646626457]
@test rand(MersenneTwister(0), Int64, 1) == [172014471070449746]
A = zeros(Int64, 2, 2)
rand!(MersenneTwister(0), A)
@test A == [ 172014471070449746 -193283627354378518;
-4679130500738884555 -9008350441255501549]

# randn
@test randn(MersenneTwister(42)) == -0.5560268761463861
Expand Down