RNG object that provides random numbers from a pre-defined vector

General Usage
1

Hi all,

I could need some help to understand how I can implement (clean) new RNGs. I have the following in mind: Given I have a (large) set of random numbers (e.g. real random numbers from some experiment) or a sequence of random numbers that I want to be able to manipulate [yes, I really want to do this and I know what I am doing ;)], how could I construct an RNG object which simply goes through the array and returns the next random number, which I can pass to my julia code that is already running with MersenneTwister?

I have the following “working” solution derived from Random.jl, but it is not complete and possibly not the cleanest solution:

using Random
import Random: rand, SamplerTrivial, CloseOpen12_64, CloseOpen01_64, BitInteger, UInt52Raw, CloseOpen12, SamplerUnion, SamplerType

mutable struct CustomRandomNumbers <: AbstractRNG
    random_numbers::Vector{Float64}
    index_current::Int

    function CustomRandomNumbers(random_numbers::Vector{Float64})
        new(random_numbers, 0)
    end
end

function gen_rand01(r::CustomRandomNumbers)
    r.index_current += 1
    if checkbounds(Bool, r.random_numbers, r.index_current)
        return @inbounds r.random_numbers[r.index_current]
    else
        throw(DomainError(r.index_current, "CustomRandomNumbers object has no further random numbers left and is configured in mode `:static`"))
    end
end

### helper functions

rand_inbounds(r::CustomRandomNumbers, ::CloseOpen12_64) = gen_rand01(r) + 1
rand_inbounds(r::CustomRandomNumbers, ::CloseOpen01_64=CloseOpen01()) =
    rand_inbounds(r, CloseOpen12()) - 1.0

rand_inbounds(r::CustomRandomNumbers, ::UInt52Raw{T}) where {T<:BitInteger} =
    reinterpret(UInt64, rand_inbounds(r, CloseOpen12())) % T
    
#### generation

function rand(r::CustomRandomNumbers, x::SamplerTrivial{UInt52Raw{UInt64}})
    rand_inbounds(r, x[])
end

rand(r::CustomRandomNumbers, sp::SamplerTrivial{CloseOpen12_64}) = rand_inbounds(r, sp[])

function test_custom_rng()
    seed = 1000
    rng = CustomRandomNumbers(rand(MersenneTwister(seed),100))

    rng_MT = MersenneTwister(seed)

    println("uniform floats:")
    for (r1, r2) in zip(rand(rng, 10), rand(rng_MT,10))
        println(r1, " vs ", r2)
    end

    println("\nexponentially distributed floats:")
    for (r1, r2) in zip(randexp(rng, 10), randexp(rng_MT,10))
        println(r1, " vs ", r2)
    end

    println("\nrandom int from range:")
    rand(rng, 1:10)
end

The first two tests work, the last sampling from a range does not work. Obviously I failed in porting all generating functions from MersenneTwister to my example, in part because I did not understand what they were really doing. Can somebody help me?

Or am I approaching the problem from the wrong angle and there is a simple and straight-forward alternative?

Thanks for any help.

2

Did you look at the docs? See

https://docs.julialang.org/en/v1/stdlib/Random/#Creating-new-generators-1

The API is not 100% specified, but neither is your problem: do you just want to recover

as eg Float64, or use them as a source of entropy in some way?

3

Thanks for the hint, but of course I looked at the docs. I have to say though that the docs are not very comprehensible to me at this point. It seemed to me that the Sampler concept is more on the level of “drawing from some specific distribution”.

What I want is literally what I present in the example: I start with a set of random numbers (Float64) and I want to use these numbers one after another for some function that generates some stochastic dynamics. Let me be more specific maybe. I have a stochastic process that creates a trajectory using an RNG. It is defined as

function my_stochastic_process(rng::AbstractRNG, some_initial_condition)
...
end

The rng is used to generate random uniform Float64, random exponentially distributed Float64, random integers in the range 1:N. I could now rewrite the function to work with random numbers not from rng but from some vector, but I thought it would be much cooler if I could use an RNG object instead. Indeed, the code that I presented in my example above works! Up to the point where I draw rand(rng, 1:N), which for now I solved as ceil(Int, rand(rng)*N). But I would consider this not optimal.

4

I would recommend the following:

  1. return Float64s directly (ie implement rand(::YourRNG, ::Type{Float64}) as popping elements from the vector,
  2. use the generic randexp fallback by implementing a method for UInt52Raw, which is an internal type for basically the mantissa of floats (look at the source),
  3. to sample for ranges, you basically need to implement rand for integers. look at the code surrounding SamplerRangeNDL in generation.jl.

As I said, this part of the API is somewhat underdocumented, so you may need to read some source code. If you believe that it could be clarified, please do make a PR, fresh eyes are important.

1 Like