Create a custom Random.Sampler or Matrix sampler via Distribution.Sampleable

General Usage
1

I am trying to define rand methods to “randomize” points x represented in a matrix (n, s) (n is the number of point and s is the dimension).
I could directly define

function rand(RM::RandomizationMethod)
     x = RM.x
     return shift(x)
end

where shift(x) can be a random translation e.g. shift(x) = x+rand(size(x)) and RM is a structure containing information on the original set of point to randomize (here directly the points but it can be their bit representation in some base and other properties).
All my different randomization methods are gathered in a type RandomizationMethod.

My questions:

  • I have read about Random.Sampler but I am not sure to understand what they actually do (especially for nontrivial cases) and why one need them.
    In my case are they necessary?
    From what I understood, if I manage to define the appropriate Sampler then rand and all associate functions rand! or rand(RM, 5) (for 5 randomization so a (n,s,5) array) or rand(RM, 5,2) (for a (n,s,5,2) array) should work almost directly.
    In particular, I could have a Sampler for one randomization and a slightly different if I want to generate multiple ones more efficiently.

  • I saw Matrix Sampler from the Distribution.jl package. It seems to be doing sort of what I want with the Sampleable type.
    Is it more appropriate for my problem?

  • I am a bit lost on how to properly use that for matrix sampling: especially the part to specialize my code for multiple realization.
    Here is a sort of MWE of my issue

using Random
using Distributions: Matrixvariate, Continuous
import Distributions: Sampleable
import Random: _rand!, rand!

abstract type SamplerShifting{F<:Matrixvariate,S<:Continuous} <: Sampleable{F,S} end
Base.size(sampler::SamplerShifting) = size(sampler.points) # the size of each matrix sample
struct Shifter{F<:AbstractMatrix{<:Real}} <: SamplerShifting{Matrixvariate,Continuous}
    points::F
    #... potentially other things
end

function shift!(rng, sampler, x)
    for i in eachindex(x)
        x[i] = sampler.points[i] + rand(rng)
    end
end

function _rand!(rng::AbstractRNG, sampler::SamplerShifting, A::DenseMatrix{T}) where {T<:Real}
    # ... generate a single matrix sample to A
    #... potentially some initialization
    shift!(rng, sampler, A)
    return A #? is it not bad to return the value here for performance ? 
end

function rand!(rng::AbstractRNG, s::Sampleable{Matrixvariate}, A::AbstractMatrix)
    size(A) == size(s) ||
        throw(DimensionMismatch("Output size inconsistent with sample length."))
    _rand!(rng, s, A)
end

m = 5
n = 2^m
s = 10
x = zeros(n, s)
sp = Shifter(x)

rand(sp, 5, 2) #! this work! Returns a Matrix of my matrix samples

#* Now I want to generate mutliple samples to do all the initialization only once
#* so I try the following (which do not work)
function _rand!(rng::AbstractRNG, sampler::SamplerShifting, A::Array{DenseMatrix{T}}) where {T<:Real}
    # ... generate a multiple matrix sample to A
    #... potentially some initialization
    for i in eachindex(A)
        shift!(rng, sampler, view(A, i)) #? this in place replacment do not work proprely
    end
    return A
end

function rand!(rng::AbstractRNG, s::Sampleable{Matrixvariate}, A::Array{DenseMatrix{T}}) where {T<:Real}
    size(A[firstindex(A)]) == size(s) ||
        throw(DimensionMismatch("Output size inconsistent with sample length."))
    _rand!(rng, s, A)
end

rand(sp, 5, 2) #!! Run but this does NOT use the new method for multiple call 
A = [x, x]
rand!(sp, A) #! Run but produce the same realization