random-polytope generation

docs/src/lib/interfaces.md

@@ -132,6 +132,7 @@ This interface defines the following functions:
 ```@docs
 an_element(::AbstractHPolygon{N}) where {N<:Real}
 ∈(::AbstractVector{Real}, ::AbstractHPolygon{Real})
+rand(::Type{HPOLYGON}) where {HPOLYGON<:AbstractHPolygon}
 vertices_list(::AbstractHPolygon{Real})
 tohrep(::AbstractHPolygon{Real})
 tovrep(::AbstractHPolygon{Real})

docs/src/lib/representations.md

@@ -351,6 +351,7 @@ VPolygon
 σ(::AbstractVector{Real}, ::VPolygon{Real})
 ∈(::AbstractVector{Real}, ::VPolygon{Real})
 an_element(::VPolygon{N}) where {N<:Real}
+rand(::Type{VPolygon})
 vertices_list(::VPolygon)
 tohrep(::VPolygon)
 tovrep(::VPolygon)
@@ -391,7 +392,7 @@ HPolytope
 HPolyhedron
 ```
 
-The following methods are shared between polyhedra and polytopes. 
+The following methods are shared between `HPolytope` and `HPolyhedron`.
 
 ```@docs
 dim(::HPoly{Real})
@@ -406,8 +407,6 @@ isempty(::HPoly{N}) where {N<:Real}
 convex_hull(::HPoly{Real}, ::HPoly{Real})
 cartesian_product(::HPoly{N}, ::HPoly{N}) where {N<:Real}
 tovrep(::HPoly{Real})
-vertices_list(::HPolyhedron{Real})
-singleton_list(::HPolyhedron{N}) where {N<:Real}
 polyhedron(::HPoly)
 ```
 
@@ -417,15 +416,37 @@ Inherited from [`LazySet`](@ref):
 * [`diameter`](@ref diameter(::LazySet, ::Real))
 
 Inherited from [`AbstractPolytope`](@ref):
-* [`singleton_list`](@ref singleton_list(::AbstractPolytope))
 * [`linear_map`](@ref linear_map(::AbstractMatrix, ::AbstractPolytope))
 
+#### Polytopes in constraint representation
+
+The following methods are specific for `HPolytope`.
+
+```@docs
+rand(::Type{HPolytope})
+vertices_list(::HPolytope{Real})
+```
+
+Inherited from [`AbstractPolytope`](@ref):
+* [`singleton_list`](@ref singleton_list(::AbstractPolytope))
+The following methods are specific for polytopes.
+
+#### Polyhedra
+
+The following methods are specific for `HPolyhedron`.
+
+```@docs
+vertices_list(::HPolyhedron{Real})
+singleton_list(::HPolyhedron{N}) where {N<:Real}
+```
+
 ### Vertex representation
 
 ```@docs
 VPolytope
 dim(::VPolytope)
 σ(::AbstractVector{Real}, ::VPolytope{Real})
+rand(::Type{VPolytope})
 vertices_list(::VPolytope)
 cartesian_product(::VPolytope{N}, ::VPolytope{N}) where N
 polyhedron(::VPolytope)

src/AbstractHPolygon.jl

@@ -1,4 +1,5 @@
-import Base.∈
+import Base: rand,
+             ∈
 
 export AbstractHPolygon,
        an_element,
@@ -128,6 +129,7 @@ function constraints_list(P::AbstractHPolygon{N})::Vector{LinearConstraint{N}} w
     return P.constraints
 end
 
+
 # --- LazySet interface functions ---
 
 
@@ -181,6 +183,47 @@ function ∈(x::AbstractVector{N}, P::AbstractHPolygon{N})::Bool where {N<:Real}
     return true
 end
 
+"""
+    rand(::Type{HPOLYGON}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
+         [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing
+        )::HPOLYGON{N} where {HPOLYGON<:AbstractHPolygon}
+
+Create a random polygon in constraint representation.
+
+### Input
+
+- `HPOLYGON`        -- type for dispatch
+- `N`               -- (optional, default: `Float64`) numeric type
+- `dim`             -- (optional, default: 2) dimension
+- `rng`             -- (optional, default: `GLOBAL_RNG`) random number generator
+- `seed`            -- (optional, default: `nothing`) seed for reseeding
+- `num_constraints` -- (optional, default: `-1`) number of constraints of the
+                       polygon (see comment below)
+
+### Output
+
+A random polygon in constraint representation.
+
+### Algorithm
+
+We create a random polygon in vertex representation and convert it to constraint
+representation.
+See [`rand(::Type{VPolygon})`](@ref).
+"""
+function rand(::Type{HPOLYGON};
+              N::Type{<:Real}=Float64,
+              dim::Int=2,
+              rng::AbstractRNG=GLOBAL_RNG,
+              seed::Union{Int, Nothing}=nothing,
+              num_constraints::Int=-1
+             )::HPOLYGON{N} where {HPOLYGON<:AbstractHPolygon}
+    @assert dim == 2 "a polygon must have dimension 2"
+    rng = reseed(rng, seed)
+    vpolygon = rand(VPolygon; N=N, dim=dim, rng=rng, seed=seed,
+                    num_vertices=num_constraints)
+    return convert(HPOLYGON, vpolygon)
+end
+
 
 # --- common AbstractHPolygon functions ---
 

src/HPolytope.jl

@@ -1,3 +1,5 @@
+import Base.rand
+
 export HPolytope,
        vertices_list
 
@@ -38,9 +40,55 @@ function HPolytope(A::AbstractMatrix{N}, b::AbstractVector{N}) where {N<:Real}
     return HPolytope(constraints)
 end
 
-# ==========================================
-# Lower level methods that use Polyhedra.jl
-# ==========================================
+
+# --- LazySet interface functions ---
+
+
+"""
+    rand(::Type{HPolytope}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
+         [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing
+        )::HPolytope{N}
+
+Create a random polytope in constraint representation.
+
+### Input
+
+- `HPolytope`    -- type for dispatch
+- `N`            -- (optional, default: `Float64`) numeric type
+- `dim`          -- (optional, default: 2) dimension
+- `rng`          -- (optional, default: `GLOBAL_RNG`) random number generator
+- `seed`         -- (optional, default: `nothing`) seed for reseeding
+- `num_vertices` -- (optional, default: `-1`) upper bound on the number of
+                    vertices of the polytope (see comment below)
+
+### Output
+
+A random polytope in constraint representation.
+
+### Algorithm
+
+We create a random polytope in vertex representation and convert it to
+constraint representation (hence the argument `num_vertices`).
+See [`rand(::Type{VPolytope})`](@ref).
+"""
+function rand(::Type{HPolytope};
+              N::Type{<:Real}=Float64,
+              dim::Int=2,
+              rng::AbstractRNG=GLOBAL_RNG,
+              seed::Union{Int, Nothing}=nothing,
+              num_vertices::Int=-1
+             )::HPolytope{N}
+    @assert isdefined(Main, :Polyhedra) "the function `rand` needs the " *
+                                        "package 'Polyhedra' loaded"
+    rng = reseed(rng, seed)
+    vpolytope = rand(VPolytope; N=N, dim=dim, rng=rng, seed=seed,
+                    num_vertices=num_vertices)
+    return convert(HPolytope, vpolytope)
+end
+
+
+# --- functions that use Polyhedra.jl ---
+
 
 function load_polyhedra_hpolytope() # function to be loaded by Requires
 return quote

src/VPolygon.jl

@@ -1,4 +1,5 @@
-import Base: ∈
+import Base: rand,
+             ∈
 
 export VPolygon
 
@@ -276,6 +277,123 @@ function ∈(x::AbstractVector{N}, P::VPolygon{N})::Bool where {N<:Real}
     return true
 end
 
+"""
+    rand(::Type{VPolygon}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
+         [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing
+        )::VPolygon{N}
+
+Create a random polygon in vertex representation.
+
+### Input
+
+- `VPolygon`     -- type for dispatch
+- `N`            -- (optional, default: `Float64`) numeric type
+- `dim`          -- (optional, default: 2) dimension
+- `rng`          -- (optional, default: `GLOBAL_RNG`) random number generator
+- `seed`         -- (optional, default: `nothing`) seed for reseeding
+- `num_vertices` -- (optional, default: `-1`) number of vertices of the
+                    polygon (see comment below)
+
+### Output
+
+A random polygon in vertex representation.
+
+### Algorithm
+
+We follow the idea [here](https://stackoverflow.com/a/47358689) based on
+[P. Valtr. Probability that n random points are in convex
+position](https://link.springer.com/article/10.1007%2FBF01271274).
+
+The number of vertices can be controlled with the argument `num_vertices`.
+For a negative value we choose a random number in the range `1:10`.
+"""
+function rand(::Type{VPolygon};
+              N::Type{<:Real}=Float64,
+              dim::Int=2,
+              rng::AbstractRNG=GLOBAL_RNG,
+              seed::Union{Int, Nothing}=nothing,
+              num_vertices::Int=-1
+             )::VPolygon{N}
+    @assert dim == 2 "a polygon must have dimension 2"
+    rng = reseed(rng, seed)
+    if num_vertices < 0
+        num_vertices = rand(1:10)
+    end
+
+    # special cases, 0 or 1 vertex
+    if num_vertices <= 1
+        vertices = [randn(rng, N, 2) for i in 1:num_vertices]
+        return VPolygon(vertices; apply_convex_hull=false)
+    end
+
+    # general case, >= 2 vertices
+
+    function remove_duplicates_sorted!(v::AbstractVector)
+        for i in length(v)-1:-1:1
+            if v[i] == v[i+1]
+                splice!(v, i+1)
+            end
+        end
+        return v
+    end
+
+    function random_axis_aligned_vectors(rng, N, n)
+        # generate a sorted list of random x and y coordinates
+        list = sort!(randn(rng, N, n))
+        while (length(remove_duplicates_sorted!(list)) < n)
+            # make sure that no duplicates exist
+            list = sort!(append!(list, randn(rng, N, length(list) - n)))
+        end
+        # lists of consecutive points
+        in_l1 = rand(rng, Bool, n-2)
+        l1 = Vector{N}() # normal
+        l2 = Vector{N}() # inverted
+        push!(l1, list[1])
+        push!(l2, list[1])
+        for i in 2:n-1
+            push!(in_l1[i-1] ? l1 : l2, list[i])
+        end
+        push!(l1, list[end])
+        push!(l2, list[end])
+        # convert to vectors representing the distance (order does not matter)
+        dist = Vector{N}()
+        sizehint!(dist, n)
+        for i in 1:length(l1)-1
+            push!(dist, l1[i+1] - l1[i])
+        end
+        for i in 1:length(l2)-1
+            push!(dist, l2[i] - l2[i+1])
+        end
+        @assert isapprox(sum(dist), zero(N), atol=1e-6)
+        return dist
+    end
+
+    # get random horizontal and vertical vectors
+    horiz = random_axis_aligned_vectors(rng, N, num_vertices)
+    vert = random_axis_aligned_vectors(rng, N, num_vertices)
+
+    # randomly combine horizontal and vertical vectors
+    m = num_vertices
+    directions = Vector{Vector{N}}(undef, num_vertices)
+    shuffle(rng, vert)
+    for (i, x) in enumerate(horiz)
+        y = splice!(vert, rand(rng, 1:m))
+        directions[i] = [x, y]
+        m -= 1
+    end
+    sort!(directions, lt=<=) # sort by angle
+
+    # connect directions
+    vertices = Vector{Vector{N}}(undef, num_vertices)
+    # random starting point
+    vertices[1] = randn(rng, N, 2)
+    for i in 1:length(directions)-1
+        vertices[i+1] = vertices[i] + directions[i]
+    end
+    @assert isapprox(vertices[end] + directions[end], vertices[1], atol=1e-6)
+    return VPolygon(vertices; apply_convex_hull=true)
+end
+
 """
     constraints_list(P::VPolygon{N})::Vector{LinearConstraint{N}} where {N<:Real}