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