Generalize random sampling for other distributions

[mforets]

Add a methods to sample from a hyperrectangle given a multivariate Gaussian distribution (or other distributins), may use Distributions.jl as an optional dependency (it's already setup).

cc @AnderGray

[AnderGray]

Yes it should be quite easy. Although only when the set is a hyperrectange. Not sure otherwise. Also I guess this only works for bounded distributions? Or perhaps you can truncate the distribution if it's not bounded, like the Gaussian. Although this will change the distribution... Maybe we can display a warning and it's ok.

Currently I believe that only the multivariate gaussian is available in Distributions.jl. However, using a gaussian copula sampler (which is easy to implement), you can use any set of marginal distributions. Although the correlations will always be linear.

@mforets I would need some help. Can you point me to where the sampling is done/ what is the current function?

[mforets]

The sampling methods are in the file samples.jl. I've sent a (just starting) PR, #2731, which implements some planned changes in the way the sampler methods are used internally.

Regarding the new method of the current issue, the implementation might be:

# the new function should be inside the following block of code because Distributions.jl is an optional dependency

function load_distributions_samples()
return quote

# add a new structure that represents the distribution
struct MultivariateGaussian <: AbstractSampler
     # ...
end

# add a new in-place sampling function to dispatch on the set type on which it is applicable, in this case,
# hyperrectangular sets
function sample!(D::AbstractVector{VT}, X::AbstractHyperrectangle, sampler::MultivariateGaussian; kwargs...) where {N, VN<:AbstractVector{N}}
    m = length(D) # number of points to sample
    lo = low(X) # vector with lower bounds of X
    hi = high(X) # vector with higher bounds of X

    # ...
    
   return D
end

The new functionality can be tested like this:

using LazySets, Plots

H = rand(Hyperrectangle) # pass dim=3 to change the dimension, default is 2

Suni = sample(H, 100); # default: uses uniform sampling (doesn't require Distributions.jl)

using Distributions   # it's an optional dependency

Sgauss = sample(H, 100, sampler=MultivariateGaussian(...));

# compare
plot(H)
plot!(Singleton.(Suni))
plot!(Singleton.(Sgauss))

To get more familiar with the type hierarchy you may have a look at https://juliareach.github.io/LazySets.jl/dev/man/tour/#Exploring-the-type-hierarchy In the current issue i suggested to dispatch on AbstractHyperrectangle as it'll cover Hyperrectangle as a special case.

[mforets]

After the revision in #2731, I actually don't know if there is anything to be added to handle distributions which are not uniform.

Below is an example, where I passed one Gaussian per coordinate. @AnderGray is that what you had in mind?

using LazySets, Plots, Distributions

using LazySets: RejectionSampler, center, sample

H = rand(Hyperrectangle)

N = 1000

S = sample(H, N)
c = LazySets.center(H)

# note that if you pass tight=true to RejectionSampler, it will assume that you know what you're doing, ie.
# it won't check whether the samples drawn from the distribution actually belong to the set; by default we have tight=false
sampler = RejectionSampler([Normal(c[1], 0.15), Normal(c[2], 0.15)])

plot(H)
plot!(Singleton.(S))

S2 = sample(H, N, sampler=sampler);
plot!(Singleton.(S2), c=:green, ratio=1.)

[mforets]

I think we can consider this closed. We can open a new issue if we need to generalize the sample method further to work with Gaussian copula samplers.