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

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.

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?

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.

I would recommend the following:

1. return `Float64`s 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