#235 - Random set generation functions

[schillic]

Closes #235.

• The only basic set type missing is HPolyhedron. I am not sure how to randomly generate that. One idea is to create an HPolytope and randomly remove constraints.

[schillic]

Side comment: When we merge this, we will locally get the following warning whenever LazySets is loaded in versions >= v0.7:

┌Warning: Package LazySets does not have Random in its dependencies:
│ - If you have LazySets checked out for development and have
│   added Random as a dependency but haven't updated your primary
│   environment's manifest file, try `Pkg.resolve()`.
│ - Otherwise you may need to report an issue with LazySets
└ Loading Random into LazySets from project dependency, future warnings for LazySets are suppressed.

This is resolved by

pkg> resolve

[mforets]

for debugging i find useful if the error message is slightly more informative, such as: "cannot create a random HPolygon of dimension $dim"

[mforets]

i suggest to take at least rand(3:10), to avoid some corner cases

[schillic]

I was also thinking about that, but sometimes it is helpful to have these corner cases for testing (my main application for random sets is random testing). Otherwise you will systematically exclude these corner cases. You can control this number with num_vertices if you want to, so I think we can keep it like that.

[schillic]

Thinking about it, my argument also holds in the other direction. So I will change it.

[mforets]

the empty constructor (case num_vertices=0) does not work with the kwarg apply_convex_hull

julia> VPolygon()
VPolygon{Float64}(Array{Float64,1}[])

julia> VPolygon([]; apply_convex_hull=false)
ERROR: MethodError: no method matching VPolygon(::Array{Any,1}; apply_convex_hull=false)
Closest candidates are:
  VPolygon() at /Users/forets/.julia/dev/LazySets/src/VPolygon.jl:53 got unsupported keyword argument "apply_convex_hull"
  VPolygon(::Array{Array{N<:Real,1},1}; apply_convex_hull, algorithm) where N<:Real at /Users/forets/.julia/dev/LazySets/src/VPolygon.jl:44
Stacktrace:
 [1] top-level scope at none:0

don't know if this is a bug; but we can do here

if num_vertices == 0
    return VPolygon()    # convex hull is not applied
elseif num_vertices == 1
    return VPolygon([randn(rng, N, 2)])
end

[schillic]

You have to add a type to the array (as the message suggested):

julia> VPolygon(Vector{N}[]; apply_convex_hull=false)
VPolygon{Float64}(Array{Float64,1}[])

[mforets]

can one substitute random_axis_aligned_vectors to picking n-elements randomly from a normal distribution? or what do these have in common?

julia> random_axis_aligned_vectors(rng, N, 5)
5-element Array{Float64,1}:
  0.12830892377309455
  2.387169565427545
  0.3232363202950286
 -0.29623066338265214
 -2.542484146113016

julia> random_axis_aligned_vectors(rng, N, 5)
5-element Array{Float64,1}:
  0.9638834005050456
  1.0636721009415822
  0.8336680989665766
 -1.5437334714482494
 -1.317490128964955

julia> random_axis_aligned_vectors(rng, N, 5)
5-element Array{Float64,1}:
  1.8457700493139375
  0.9765805620918393
 -0.40900877044843353
 -0.5103436430987645
 -1.9029981978585788

[schillic]

The sum is zero. Have a look at the video here.

[mforets]

Very nice video.

I would suggest to attach a docstring and rename random_axis_aligned_vectors, maybe to random_zero_sum_vector.

[mforets]

would this work?

function random_zero_sum_vector(rng, N, n)
    list = randn(rng, N, n)
    return [list[1:1:end-1] - list[2:1:end];
            list[end] - list[1]]
end

[schillic]

Yes, but the numbers are not uniformly distributed. The article I linked to claims that the algorithm guarantees that, and I think it is good to have this property (if you can have it).

[mforets]

alright, i think this function is interesting in its own right and i think we can put it outside (eventually call it _random_zero_sum_vector, we had this discussion 😄 ) and add docstring with reference to the article.

[mforets]

A general comment (that could be addressed in another PR): rand(BallInf, 5) has the same behavior as rand(Int, 5), ie. to create 5 BallInf's in an array. This is often useful because one can also create arrays of any order easily. However, one is subject to the default options, in particular rand(BallInf, 5, dim=2) doesn't work.

[schillic]

Yes, this feature is thanks to

# hook into random API
import Random.rand
function rand(rng::AbstractRNG, ::SamplerType{T}) where T<:LazySet
    rand(T, rng=rng)
end

I am not sure if this API is made for having keyword arguments.
You can always do [rand(BallInf; dim=mydim) for i in 1:5]. Not too bad I think.

[schillic]

The only basic set type missing is HPolyhedron. I am not sure how to randomly generate that. One idea is to create an HPolytope and randomly remove constraints.

What do you think about that?

[mforets]

Since HPolyhedron are HPolytope, one can just run the rand(HPolytope) and do a convert(HPolyhedron, ..). Doing this and then removing constraints sounds good as well. In general we don't give guarantees on the distribution of the constructed set (that's another story).

[schillic]

I added all the fixes you proposed.