From 2dffd8b92e2f52f94034294685719f60d70ad5f8 Mon Sep 17 00:00:00 2001 From: schillic Date: Sat, 3 Nov 2018 16:15:38 +0100 Subject: [PATCH] random-polytope generation --- docs/src/lib/representations.md | 31 +++++++-- src/AbstractHPolygon.jl | 45 +++++++++++- src/HPolytope.jl | 54 +++++++++++++- src/VPolygon.jl | 120 +++++++++++++++++++++++++++++++- src/VPolytope.jl | 55 ++++++++++++++- test/unit_Polygon.jl | 5 ++ test/unit_Polytope.jl | 8 +++ 7 files changed, 308 insertions(+), 10 deletions(-) diff --git a/docs/src/lib/representations.md b/docs/src/lib/representations.md index bbfbd5cb20..31588d9854 100644 --- a/docs/src/lib/representations.md +++ b/docs/src/lib/representations.md @@ -293,6 +293,7 @@ Inherited from [`AbstractCentrallySymmetricPolytope`](@ref): ```@docs HPolygon σ(::AbstractVector{Real}, ::HPolygon{Real}) +rand(::Type{HPolygon}) ``` Inherited from [`LazySet`](@ref): * [`norm`](@ref norm(::LazySet, ::Real)) @@ -321,6 +322,7 @@ Inherited from [`AbstractHPolygon`](@ref): ```@docs HPolygonOpt σ(::AbstractVector{Real}, ::HPolygonOpt{Real}) +rand(::Type{HPolygonOpt}) ``` Inherited from [`LazySet`](@ref): * [`norm`](@ref norm(::LazySet, ::Real)) @@ -351,6 +353,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 +394,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 +409,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 +418,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) diff --git a/src/AbstractHPolygon.jl b/src/AbstractHPolygon.jl index 4337cef797..bae0304ab7 100644 --- a/src/AbstractHPolygon.jl +++ b/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 --- diff --git a/src/HPolytope.jl b/src/HPolytope.jl index 3f29f89a32..58a6751227 100644 --- a/src/HPolytope.jl +++ b/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 diff --git a/src/VPolygon.jl b/src/VPolygon.jl index ff7ee1103f..5fc633361e 100644 --- a/src/VPolygon.jl +++ b/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} diff --git a/src/VPolytope.jl b/src/VPolytope.jl index 14395f0e73..46e59d6691 100644 --- a/src/VPolytope.jl +++ b/src/VPolytope.jl @@ -1,5 +1,7 @@ using MathProgBase, GLPKMathProgInterface +import Base.rand + export VPolytope, vertices_list @@ -97,7 +99,54 @@ function σ(d::AbstractVector{N}, end end -# --- VPolytope functions --- +""" + rand(::Type{VPolytope}; [N]::Type{<:Real}=Float64, [dim]::Int=2, + [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing + )::VPolytope{N} + +Create a random polytope in vertex representation. + +### Input + +- `VPolytope` -- 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 vertex representation. + +### Algorithm + +All numbers are normally distributed with mean 0 and standard deviation 1. + +The number of vertices can be controlled with the argument `num_vertices`. +For a negative value we choose a random number in the range `dim:5*dim` (except +if `dim == 1`, in which case we choose in the range `1:2`). +Note that we do not guarantee that the vertices are not redundant. +""" +function rand(::Type{VPolytope}; + N::Type{<:Real}=Float64, + dim::Int=2, + rng::AbstractRNG=GLOBAL_RNG, + seed::Union{Int, Nothing}=nothing, + num_vertices::Int=-1 + )::VPolytope{N} + rng = reseed(rng, seed) + if num_vertices < 0 + num_vertices = (dim == 1) ? rand(1:2) : rand(dim:5*dim) + end + vertices = [randn(rng, N, dim) for i in 1:num_vertices] + return VPolytope(vertices) +end + + +# --- AbstractPolytope interface functions --- + """ vertices_list(P::VPolytope{N})::Vector{Vector{N}} where {N<:Real} @@ -138,6 +187,10 @@ function constraints_list(P::VPolytope{N})::Vector{LinearConstraint{N}} where {N return constraints_list(tohrep(P)) end + +# --- functions that use Polyhedra.jl --- + + function load_polyhedra_vpolytope() # function to be loaded by Requires return quote # see the interface file AbstractPolytope.jl for the imports diff --git a/test/unit_Polygon.jl b/test/unit_Polygon.jl index 0642521ca9..5256899c25 100644 --- a/test/unit_Polygon.jl +++ b/test/unit_Polygon.jl @@ -1,4 +1,9 @@ for N in [Float64, Float32, Rational{Int}] + # random polygons + rand(HPolygon) + rand(HPolygonOpt) + rand(VPolygon) + # Empty polygon p = HPolygon{N}() diff --git a/test/unit_Polytope.jl b/test/unit_Polytope.jl index e2ce64953c..5e3466b5b6 100644 --- a/test/unit_Polytope.jl +++ b/test/unit_Polytope.jl @@ -1,6 +1,14 @@ global test_suite_polyhedra for N in [Float64, Rational{Int}, Float32] + # random polytopes + if test_suite_polyhedra + rand(HPolytope) + else + @test_throws AssertionError rand(HPolytope) + end + rand(VPolytope) + # ----- # H-rep # -----