From 16148e648322a88fef1fe252c88241679aa5797a Mon Sep 17 00:00:00 2001 From: schillic Date: Tue, 16 Jul 2024 12:08:51 +0200 Subject: [PATCH 1/6] move HPolyhedron.jl to subfolder --- src/LazySets.jl | 2 +- src/Sets/{HPolyhedron.jl => HPolyhedron/HPolyhedronModule.jl} | 0 2 files changed, 1 insertion(+), 1 deletion(-) rename src/Sets/{HPolyhedron.jl => HPolyhedron/HPolyhedronModule.jl} (100%) diff --git a/src/LazySets.jl b/src/LazySets.jl index 35cb0bfbd2..c5c4a80237 100644 --- a/src/LazySets.jl +++ b/src/LazySets.jl @@ -126,7 +126,7 @@ include("Sets/HPolygonOpt.jl") include("Sets/HPolytope/HPolytopeModule.jl") @reexport using ..HPolytopeModule: HPolytope -include("Sets/HPolyhedron.jl") +include("Sets/HPolyhedron/HPolyhedronModule.jl") include("Sets/HParallelotope/HParallelotopeModule.jl") @reexport using ..HParallelotopeModule: HParallelotope, diff --git a/src/Sets/HPolyhedron.jl b/src/Sets/HPolyhedron/HPolyhedronModule.jl similarity index 100% rename from src/Sets/HPolyhedron.jl rename to src/Sets/HPolyhedron/HPolyhedronModule.jl From 6923d80c9e47c81d42f7d79f64c4a117c0f9ebdd Mon Sep 17 00:00:00 2001 From: schillic Date: Tue, 16 Jul 2024 12:40:14 +0200 Subject: [PATCH 2/6] define stubs (addconstraint, is_hyperplanar, _is_halfspace, _is_hyperplane) --- docs/src/lib/interfaces/AbstractPolyhedron.md | 7 +++ .../AbstractPolyhedron_functions.jl | 43 ++++++++++++++++++- src/Sets/HPolyhedron/HPolyhedronModule.jl | 8 +--- 3 files changed, 50 insertions(+), 8 deletions(-) diff --git a/docs/src/lib/interfaces/AbstractPolyhedron.md b/docs/src/lib/interfaces/AbstractPolyhedron.md index 366c1a6fd8..54980ae9e3 100644 --- a/docs/src/lib/interfaces/AbstractPolyhedron.md +++ b/docs/src/lib/interfaces/AbstractPolyhedron.md @@ -42,6 +42,13 @@ LazySets._isbounded_stiemke LazySets._linear_map_polyhedron ``` +Some common functions implemented by several subtypes: + +```@docs +addconstraint!(::AbstractPolyhedron, ::HalfSpace) +is_hyperplanar(::AbstractPolyhedron) +``` + Some common functions to work with linear constraints: ```@docs diff --git a/src/Interfaces/AbstractPolyhedron_functions.jl b/src/Interfaces/AbstractPolyhedron_functions.jl index 059ac608d9..0fb5f11bf7 100644 --- a/src/Interfaces/AbstractPolyhedron_functions.jl +++ b/src/Interfaces/AbstractPolyhedron_functions.jl @@ -2,7 +2,9 @@ export constrained_dimensions, tosimplehrep, remove_redundant_constraints, remove_redundant_constraints!, - isfeasible + isfeasible, + addconstraint!, + is_hyperplanar isconvextype(::Type{<:AbstractPolyhedron}) = true @@ -1254,3 +1256,42 @@ function isfeasible(constraints::AbstractVector{<:HalfSpace}, return !_isempty_polyhedron_lp(constraints, witness; solver=solver) end end + +""" + addconstraint!(P::AbstractPolyhedron, constraint::HalfSpace) + +Add a linear constraint to a set in constraint representation in-place. + +### Input + +- `P` -- set in constraint representation +- `constraint` -- linear constraint to add + +### Notes + +It is left to the user to guarantee that the dimension of all linear constraints +is the same. +""" +function addconstraint!(::AbstractPolyhedron, ::HalfSpace) end + +""" + is_hyperplanar(P::AbstractPolyhedron) + +Determine whether a polyhedron is equivalent to a hyperplane. + +### Input + +- `P` -- polyhedron + +### Output + +`true` iff `P` is hyperplanar, i.e., consists of two linear constraints +``a·x ≤ b`` and ``-a·x ≤ -b``. +""" +function is_hyperplanar(::AbstractPolyhedron) end + +# internal function; defined here due to dependency SymEngine and submodules +function _is_halfspace() end + +# internal function; defined here due to dependency SymEngine and submodules +function _is_hyperplane() end diff --git a/src/Sets/HPolyhedron/HPolyhedronModule.jl b/src/Sets/HPolyhedron/HPolyhedronModule.jl index 75a1911b05..c83609f013 100644 --- a/src/Sets/HPolyhedron/HPolyhedronModule.jl +++ b/src/Sets/HPolyhedron/HPolyhedronModule.jl @@ -1,10 +1,4 @@ -export HPolyhedron, - addconstraint!, - tohrep, tovrep, - remove_redundant_constraints, - remove_redundant_constraints!, - constrained_dimensions, - is_hyperplanar +export HPolyhedron """ HPolyhedron{N, VN<:AbstractVector{N}} <: AbstractPolyhedron{N} From b493ac2b6d8c2f3fbf06e06f5e5e03c10371fc37 Mon Sep 17 00:00:00 2001 From: schillic Date: Tue, 16 Jul 2024 13:23:01 +0200 Subject: [PATCH 3/6] add HPolyhedronModule --- docs/src/lib/sets/HPolyhedron.md | 45 +++++++++++++---------- src/Initialization/init_Polyhedra.jl | 1 - src/Initialization/init_Symbolics.jl | 1 - src/LazySets.jl | 2 + src/Sets/HPolyhedron/HPolyhedronModule.jl | 40 +++++++++++++++++++- src/Sets/HPolyhedron/init.jl | 5 +++ src/Sets/HPolyhedron/init_LazySets.jl | 2 + 7 files changed, 73 insertions(+), 23 deletions(-) create mode 100644 src/Sets/HPolyhedron/init.jl create mode 100644 src/Sets/HPolyhedron/init_LazySets.jl diff --git a/docs/src/lib/sets/HPolyhedron.md b/docs/src/lib/sets/HPolyhedron.md index 585694b101..b2b725f494 100644 --- a/docs/src/lib/sets/HPolyhedron.md +++ b/docs/src/lib/sets/HPolyhedron.md @@ -1,5 +1,5 @@ ```@meta -CurrentModule = LazySets +CurrentModule = LazySets.HPolyhedronModule ``` # [Polyhedron in constraint representation (HPolyhedron)](@id def_HPolyhedron) @@ -8,44 +8,49 @@ CurrentModule = LazySets HPolyhedron ``` +## Operations + The following methods are shared between `HPolytope` and `HPolyhedron`. ```@docs -dim(::HPoly) -ρ(::AbstractVector{M}, ::HPoly{N}) where {M, N} -σ(::AbstractVector{M}, ::HPoly{N}) where {M, N} -addconstraint!(::HPoly, ::HalfSpace) constraints_list(::HPoly) -tohrep(::HPoly) -tovrep(::HPoly) +dim(::HPoly) normalize(::HPoly{N}, p::Real=N(2)) where {N} -translate(::HPoly, ::AbstractVector) remove_redundant_constraints(::HPoly) remove_redundant_constraints!(::HPoly) +tohrep(::HPoly) +tovrep(::HPoly) +addconstraint!(::HPoly, ::HalfSpace) +ρ(::AbstractVector{M}, ::HPoly{N}) where {M, N} +σ(::AbstractVector{M}, ::HPoly{N}) where {M, N} +translate(::HPoly, ::AbstractVector) +``` + +The following methods are specific to `HPolyhedron`. + +```@docs +rand(::Type{HPolyhedron}) ``` + +```@meta +CurrentModule = LazySets +``` + Inherited from [`LazySet`](@ref): +* [`diameter`](@ref diameter(::LazySet, ::Real)) * [`high`](@ref high(::LazySet)) +* [`isempty`](@ref isempty(::LazySet{N}, ::Bool=false) where {N}) * [`low`](@ref low(::LazySet)) * [`norm`](@ref norm(::LazySet, ::Real)) * [`radius`](@ref radius(::LazySet, ::Real)) -* [`diameter`](@ref diameter(::LazySet, ::Real)) +* [`reflect`](@ref reflect(::LazySet)) * [`singleton_list`](@ref singleton_list(::LazySet)) -* [`isempty`](@ref isempty(::LazySet{N}, ::Bool=false) where {N}) * [`linear_map`](@ref linear_map(::AbstractMatrix, ::LazySet) -* [`reflect`](@ref reflect(::LazySet)) Inherited from [`AbstractPolyhedron`](@ref): -* [`∈`](@ref ∈(::AbstractVector, ::AbstractPolyhedron)) * [`an_element`](@ref an_element(::AbstractPolyhedron)) * [`constrained_dimensions`](@ref constrained_dimensions(::AbstractPolyhedron)) - -The following methods are specific to `HPolyhedron`. - -```@docs -rand(::Type{HPolyhedron}) -``` - -Inherited from [`AbstractPolyhedron`](@ref): * [`isbounded`](@ref isbounded(::AbstractPolyhedron{N}) where {N}) * [`isuniversal`](@ref isuniversal(::AbstractPolyhedron{N}, ::Bool=false) where {N}) * [`vertices_list`](@ref vertices_list(::AbstractPolyhedron, ::Bool=false)) +* [`∈`](@ref ∈(::AbstractVector, ::AbstractPolyhedron)) diff --git a/src/Initialization/init_Polyhedra.jl b/src/Initialization/init_Polyhedra.jl index ab14e09fbc..edd6dfeb92 100644 --- a/src/Initialization/init_Polyhedra.jl +++ b/src/Initialization/init_Polyhedra.jl @@ -37,6 +37,5 @@ function _is_polyhedra_backend(backend::Polyhedra.Library) return true end -eval(load_polyhedra_hpolyhedron()) eval(load_polyhedra_mesh()) eval(load_polyhedra_lazyset()) diff --git a/src/Initialization/init_Symbolics.jl b/src/Initialization/init_Symbolics.jl index 4c1c83ee56..7daf003e7b 100644 --- a/src/Initialization/init_Symbolics.jl +++ b/src/Initialization/init_Symbolics.jl @@ -20,4 +20,3 @@ end eval(load_symbolics_hyperplane()) eval(load_symbolics_halfspace()) -eval(load_symbolics_hpolyhedron()) diff --git a/src/LazySets.jl b/src/LazySets.jl index c5c4a80237..089f02d963 100644 --- a/src/LazySets.jl +++ b/src/LazySets.jl @@ -127,6 +127,8 @@ include("Sets/HPolytope/HPolytopeModule.jl") @reexport using ..HPolytopeModule: HPolytope include("Sets/HPolyhedron/HPolyhedronModule.jl") +@reexport using ..HPolyhedronModule: HPolyhedron +using ..HPolyhedronModule: HPoly include("Sets/HParallelotope/HParallelotopeModule.jl") @reexport using ..HParallelotopeModule: HParallelotope, diff --git a/src/Sets/HPolyhedron/HPolyhedronModule.jl b/src/Sets/HPolyhedron/HPolyhedronModule.jl index c83609f013..44f38b06b0 100644 --- a/src/Sets/HPolyhedron/HPolyhedronModule.jl +++ b/src/Sets/HPolyhedron/HPolyhedronModule.jl @@ -1,3 +1,29 @@ +module HPolyhedronModule + +using Reexport, Requires + +using ..LazySets: AbstractPolyhedron, HalfSpace, default_lp_solver, + default_polyhedra_backend, iscomplement, is_lp_infeasible, + is_lp_optimal, is_lp_unbounded, has_lp_infeasibility_ray, + linprog, tosimplehrep, _isempty_polyhedron, _normal_Vector +using ..HPolytopeModule: HPolytope +using Random: AbstractRNG, GLOBAL_RNG +using ReachabilityBase.Arrays: to_negative_vector +using ReachabilityBase.Basetype: basetype +using ReachabilityBase.Comparison: isapproxzero +using ReachabilityBase.Distribution: reseed! +using ReachabilityBase.Require: require +using LinearAlgebra: dot + +@reexport import ..API: constraints_list, dim, isempty, isoperationtype, rand, + permute, ρ, σ, translate +@reexport import ..LazySets: is_hyperplanar, normalize, + remove_redundant_constraints, + remove_redundant_constraints!, tohrep, tovrep, + addconstraint! +@reexport import ..Base: convert +@reexport using ..API + export HPolyhedron """ @@ -340,6 +366,8 @@ function remove_redundant_constraints(P::HPoly; backend=nothing) if remove_redundant_constraints!(Pred; backend=backend) return Pred else # the polyhedron P is empty + require(@__MODULE__, :LazySets; fun_name="remove_redundant_constraints") + N = eltype(P) return EmptySet{N}(dim(P)) end @@ -434,7 +462,9 @@ For further information on the supported backends see [Polyhedra's documentation](https://juliapolyhedra.github.io/). """ function tovrep(P::HPoly; backend=default_polyhedra_backend(P)) + require(@__MODULE__, :LazySets; fun_name="tovrep") require(@__MODULE__, :Polyhedra; fun_name="tovrep") + P = polyhedron(P; backend=backend) return VPolytope(P) end @@ -458,7 +488,8 @@ end function load_polyhedra_hpolyhedron() # function to be loaded by Requires return quote - # see the interface file init_Polyhedra.jl for the imports + using .Polyhedra: HRep, + polyhedron """ convert(::Type{HPolyhedron}, P::HRep{N}) where {N} @@ -515,6 +546,9 @@ end function load_symbolics_hpolyhedron() return quote + using .Symbolics: Num + using ..LazySets: _get_variables, _vec, _is_halfspace, _is_hyperplane + """ HPolyhedron(expr::Vector{<:Num}, vars=_get_variables(expr); N::Type{<:Real}=Float64) @@ -597,3 +631,7 @@ function permute(P::HPoly, p::AbstractVector{Int}) T = basetype(P) return T([permute(H, p) for H in P.constraints]) end + +include("init.jl") + +end # module diff --git a/src/Sets/HPolyhedron/init.jl b/src/Sets/HPolyhedron/init.jl new file mode 100644 index 0000000000..2893f98389 --- /dev/null +++ b/src/Sets/HPolyhedron/init.jl @@ -0,0 +1,5 @@ +function __init__() + @require LazySets = "b4f0291d-fe17-52bc-9479-3d1a343d9043" include("init_LazySets.jl") + @require Polyhedra = "67491407-f73d-577b-9b50-8179a7c68029" eval(load_polyhedra_hpolyhedron()) + @require Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" eval(load_symbolics_hpolyhedron()) +end diff --git a/src/Sets/HPolyhedron/init_LazySets.jl b/src/Sets/HPolyhedron/init_LazySets.jl new file mode 100644 index 0000000000..11d965d352 --- /dev/null +++ b/src/Sets/HPolyhedron/init_LazySets.jl @@ -0,0 +1,2 @@ +using .LazySets.EmptySetModule: EmptySet +using .LazySets.VPolytopeModule: VPolytope From 442bdce8ba8e60d5762db259634c2a55fedfa1a2 Mon Sep 17 00:00:00 2001 From: schillic Date: Tue, 16 Jul 2024 13:34:17 +0200 Subject: [PATCH 4/6] outsource struct and operations to files --- src/Sets/HPolyhedron/HPolyhedron.jl | 132 ++++ src/Sets/HPolyhedron/HPolyhedronModule.jl | 578 +----------------- src/Sets/HPolyhedron/addconstraint.jl | 19 + src/Sets/HPolyhedron/constraints_list.jl | 17 + src/Sets/HPolyhedron/dim.jl | 17 + src/Sets/HPolyhedron/is_hyperplanar.jl | 17 + src/Sets/HPolyhedron/isempty.jl | 16 + src/Sets/HPolyhedron/isoperationtype.jl | 1 + src/Sets/HPolyhedron/normalize.jl | 20 + src/Sets/HPolyhedron/permute.jl | 4 + src/Sets/HPolyhedron/rand.jl | 39 ++ .../remove_redundant_constraints.jl | 66 ++ src/Sets/HPolyhedron/support_function.jl | 32 + src/Sets/HPolyhedron/support_vector.jl | 86 +++ src/Sets/HPolyhedron/tohrep.jl | 17 + src/Sets/HPolyhedron/tovrep.jl | 31 + src/Sets/HPolyhedron/translate.jl | 33 + 17 files changed, 564 insertions(+), 561 deletions(-) create mode 100644 src/Sets/HPolyhedron/HPolyhedron.jl create mode 100644 src/Sets/HPolyhedron/addconstraint.jl create mode 100644 src/Sets/HPolyhedron/constraints_list.jl create mode 100644 src/Sets/HPolyhedron/dim.jl create mode 100644 src/Sets/HPolyhedron/is_hyperplanar.jl create mode 100644 src/Sets/HPolyhedron/isempty.jl create mode 100644 src/Sets/HPolyhedron/isoperationtype.jl create mode 100644 src/Sets/HPolyhedron/normalize.jl create mode 100644 src/Sets/HPolyhedron/permute.jl create mode 100644 src/Sets/HPolyhedron/rand.jl create mode 100644 src/Sets/HPolyhedron/remove_redundant_constraints.jl create mode 100644 src/Sets/HPolyhedron/support_function.jl create mode 100644 src/Sets/HPolyhedron/support_vector.jl create mode 100644 src/Sets/HPolyhedron/tohrep.jl create mode 100644 src/Sets/HPolyhedron/tovrep.jl create mode 100644 src/Sets/HPolyhedron/translate.jl diff --git a/src/Sets/HPolyhedron/HPolyhedron.jl b/src/Sets/HPolyhedron/HPolyhedron.jl new file mode 100644 index 0000000000..f716637a0e --- /dev/null +++ b/src/Sets/HPolyhedron/HPolyhedron.jl @@ -0,0 +1,132 @@ +""" + HPolyhedron{N, VN<:AbstractVector{N}} <: AbstractPolyhedron{N} + +Type that represents a convex polyhedron in constraint representation, that is, +a finite intersection of half-spaces, +```math +P = ⋂_{i = 1}^m H_i, +``` +where each ``H_i = \\{x ∈ ℝ^n : a_i^T x ≤ b_i \\}`` is a +half-space, ``a_i ∈ ℝ^n`` is the normal vector of the ``i``-th +half-space and ``b_i`` is the displacement. The set ``P`` may or may not be +bounded (see also [`HPolytope`](@ref), which assumes boundedness). + +### Fields + +- `constraints` -- vector of linear constraints +""" +struct HPolyhedron{N,VN<:AbstractVector{N}} <: AbstractPolyhedron{N} + constraints::Vector{HalfSpace{N,VN}} + + function HPolyhedron(constraints::Vector{HalfSpace{N,VN}}) where {N,VN<:AbstractVector{N}} + return new{N,VN}(constraints) + end +end + +# constructor with no constraints +function HPolyhedron{N,VN}() where {N,VN<:AbstractVector{N}} + return HPolyhedron(Vector{HalfSpace{N,VN}}()) +end + +# constructor with no constraints, given only the numeric type +function HPolyhedron{N}() where {N} + return HPolyhedron(Vector{HalfSpace{N,Vector{N}}}()) +end + +# constructor without explicit numeric type, defaults to Float64 +function HPolyhedron() + return HPolyhedron{Float64}() +end + +# constructor with constraints of mixed type +function HPolyhedron(constraints::Vector{<:HalfSpace}) + return HPolyhedron(_normal_Vector(constraints)) +end + +# constructor from a simple constraint representation +function HPolyhedron(A::AbstractMatrix, b::AbstractVector) + return HPolyhedron(constraints_list(A, b)) +end + +function load_symbolics_hpolyhedron() + return quote + using .Symbolics: Num + using ..LazySets: _get_variables, _vec, _is_halfspace, _is_hyperplane + + """ + HPolyhedron(expr::Vector{<:Num}, vars=_get_variables(expr); + N::Type{<:Real}=Float64) + + Return a polyhedron in constraint representation given by a list of symbolic + expressions. + + ### Input + + - `expr` -- vector of symbolic expressions that describes each half-space + - `vars` -- (optional, default: `_get_variables(expr)`), if an array of + variables is given, use those as the ambient variables in the set + with respect to which derivations take place; otherwise, use only + the variables that appear in the given expression (but be careful + because the order may be incorrect; it is advised to always pass + `vars` explicitly) + - `N` -- (optional, default: `Float64`) the numeric type of the returned set + + ### Output + + An `HPolyhedron`. + + ### Examples + + ```jldoctest + julia> using Symbolics + + julia> vars = @variables x y + 2-element Vector{Num}: + x + y + + julia> HPolyhedron([x + y <= 1, x + y >= -1], vars) + HPolyhedron{Float64, Vector{Float64}}(HalfSpace{Float64, Vector{Float64}}[HalfSpace{Float64, Vector{Float64}}([1.0, 1.0], 1.0), HalfSpace{Float64, Vector{Float64}}([-1.0, -1.0], 1.0)]) + + julia> X = HPolyhedron([x == 0, y <= 0], vars) + HPolyhedron{Float64, Vector{Float64}}(HalfSpace{Float64, Vector{Float64}}[HalfSpace{Float64, Vector{Float64}}([1.0, 0.0], -0.0), HalfSpace{Float64, Vector{Float64}}([-1.0, -0.0], 0.0), HalfSpace{Float64, Vector{Float64}}([0.0, 1.0], -0.0)]) + ``` + """ + function HPolyhedron(expr::Vector{<:Num}, vars::AbstractVector{Num}; + N::Type{<:Real}=Float64) + clist = Vector{HalfSpace{N,Vector{N}}}() + sizehint!(clist, length(expr)) + got_hyperplane = false + got_halfspace = false + zeroed_vars = Dict(v => zero(N) for v in vars) + vars_list = collect(vars) + for ex in expr + exval = Symbolics.value(ex) + got_hyperplane, sexpr = _is_hyperplane(exval) + if !got_hyperplane + got_halfspace, sexpr = _is_halfspace(exval) + if !got_halfspace + throw(ArgumentError("expected an expression describing either " * + "a half-space of a hyperplane, got $expr")) + end + end + + coeffs = [N(α.val) for α in gradient(sexpr, vars_list)] + β = -N(Symbolics.substitute(sexpr, zeroed_vars)) + + push!(clist, HalfSpace(coeffs, β)) + if got_hyperplane + push!(clist, HalfSpace(-coeffs, -β)) + end + end + return HPolyhedron(clist) + end + + function HPolyhedron(expr::Vector{<:Num}; N::Type{<:Real}=Float64) + return HPolyhedron(expr, _get_variables(expr); N=N) + end + function HPolyhedron(expr::Vector{<:Num}, vars; N::Type{<:Real}=Float64) + return HPolyhedron(expr, _vec(vars); N=N) + end + end +end # quote / load_symbolics_hpolyhedron() diff --git a/src/Sets/HPolyhedron/HPolyhedronModule.jl b/src/Sets/HPolyhedron/HPolyhedronModule.jl index 44f38b06b0..826882fac6 100644 --- a/src/Sets/HPolyhedron/HPolyhedronModule.jl +++ b/src/Sets/HPolyhedron/HPolyhedronModule.jl @@ -26,465 +26,27 @@ using LinearAlgebra: dot export HPolyhedron -""" - HPolyhedron{N, VN<:AbstractVector{N}} <: AbstractPolyhedron{N} - -Type that represents a convex polyhedron in constraint representation, that is, -a finite intersection of half-spaces, -```math -P = ⋂_{i = 1}^m H_i, -``` -where each ``H_i = \\{x ∈ ℝ^n : a_i^T x ≤ b_i \\}`` is a -half-space, ``a_i ∈ ℝ^n`` is the normal vector of the ``i``-th -half-space and ``b_i`` is the displacement. The set ``P`` may or may not be -bounded (see also [`HPolytope`](@ref), which assumes boundedness). - -### Fields - -- `constraints` -- vector of linear constraints -""" -struct HPolyhedron{N,VN<:AbstractVector{N}} <: AbstractPolyhedron{N} - constraints::Vector{HalfSpace{N,VN}} - - function HPolyhedron(constraints::Vector{HalfSpace{N,VN}}) where {N,VN<:AbstractVector{N}} - return new{N,VN}(constraints) - end -end - -isoperationtype(::Type{<:HPolyhedron}) = false - -# constructor with no constraints -function HPolyhedron{N,VN}() where {N,VN<:AbstractVector{N}} - return HPolyhedron(Vector{HalfSpace{N,VN}}()) -end - -# constructor with no constraints, given only the numeric type -function HPolyhedron{N}() where {N} - return HPolyhedron(Vector{HalfSpace{N,Vector{N}}}()) -end - -# constructor without explicit numeric type, defaults to Float64 -function HPolyhedron() - return HPolyhedron{Float64}() -end - -# constructor with constraints of mixed type -function HPolyhedron(constraints::Vector{<:HalfSpace}) - return HPolyhedron(_normal_Vector(constraints)) -end - -# constructor from a simple constraint representation -function HPolyhedron(A::AbstractMatrix, b::AbstractVector) - return HPolyhedron(constraints_list(A, b)) -end +include("HPolyhedron.jl") # convenience union type const HPoly{N} = Union{HPolytope{N},HPolyhedron{N}} -""" - dim(P::HPoly) - -Return the dimension of a polyhedron in constraint representation. - -### Input - -- `P` -- polyhedron in constraint representation - -### Output - -The ambient dimension of the polyhedron in constraint representation. -If it has no constraints, the result is ``-1``. -""" -function dim(P::HPoly) - return length(P.constraints) == 0 ? -1 : length(P.constraints[1].a) -end - -""" - ρ(d::AbstractVector{M}, P::HPoly{N}; - solver=default_lp_solver(M, N)) where {M, N} - -Evaluate the support function of a polyhedron in constraint representation in a -given direction. - -### Input - -- `d` -- direction -- `P` -- polyhedron in constraint representation -- `solver` -- (optional, default: `default_lp_solver(M, N)`) the backend used to - solve the linear program - -### Output - -The evaluation of the support function for the polyhedron. -If a polytope is unbounded in the given direction, we throw an error. -If a polyhedron is unbounded in the given direction, the result is `Inf`. -""" -function ρ(d::AbstractVector{M}, P::HPoly{N}; - solver=default_lp_solver(M, N)) where {M,N} - lp, unbounded = σ_helper(d, P, solver) - if unbounded - if P isa HPolytope - error("the support function in direction $(d) is undefined " * - "because the polytope is unbounded") - end - return N(Inf) - end - return dot(d, lp.sol) -end - -""" - σ(d::AbstractVector{M}, P::HPoly{N}; - solver=default_lp_solver(M, N) where {M, N} - -Return a support vector of a polyhedron in constraint representation in a given -direction. - -### Input - -- `d` -- direction -- `P` -- polyhedron in constraint representation -- `solver` -- (optional, default: `default_lp_solver(M, N)`) the backend used to - solve the linear program - -### Output - -The support vector in the given direction. -If a polytope is unbounded in the given direction, we throw an error. -If a polyhedron is unbounded in the given direction, the result contains `±Inf` -entries. -""" -function σ(d::AbstractVector{M}, P::HPoly{N}; - solver=default_lp_solver(M, N)) where {M,N} - lp, unbounded = σ_helper(d, P, solver) - if unbounded - if P isa HPolytope - error("the support vector in direction $(d) is undefined because " * - "the polytope is unbounded") - end - return _σ_unbounded_lp(d, P, lp) - else - return lp.sol - end -end - -# construct the solution from the solver's ray result -function _σ_unbounded_lp(d, P::HPoly{N}, lp) where {N} - if isnothing(lp) - ray = d - elseif has_lp_infeasibility_ray(lp.model) - ray = lp.sol # infeasibility ray is stored as the solution - else - error("LP solver did not return an infeasibility ray") - end - - res = Vector{N}(undef, length(ray)) - e = isempty(P.constraints) ? zeros(N, length(ray)) : an_element(P) - @inbounds for i in eachindex(ray) - if isapproxzero(ray[i]) - res[i] = e[i] - elseif ray[i] > zero(N) - res[i] = N(Inf) - else - res[i] = N(-Inf) - end - end - return res -end - -function σ_helper(d::AbstractVector, P::HPoly, solver) - # represent c = -d as a Vector since GLPK does not accept sparse vectors - # (see #1011) - c = to_negative_vector(d) - - (A, b) = tosimplehrep(P) - if length(b) == 0 - unbounded = true - lp = nothing - else - sense = '<' - l = -Inf - u = Inf - lp = linprog(c, A, sense, b, l, u, solver) - if is_lp_infeasible(lp.status; strict=true) - throw(ArgumentError("the support vector is undefined because " * - "the polyhedron is empty")) - elseif is_lp_unbounded(lp.status) - unbounded = true - elseif is_lp_optimal(lp.status) - unbounded = false - else - error("got unknown LP status $(lp.status)") - end - end - return (lp, unbounded) -end - -""" - rand(::Type{HPolyhedron}; [N]::Type{<:Real}=Float64, [dim]::Int=2, - [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing) - -Create a random polyhedron. - -### Input - -- `HPolyhedron` -- type for dispatch -- `N` -- (optional, default: `Float64`) numeric type -- `dim` -- (optional, default: 2) dimension (is ignored) -- `rng` -- (optional, default: `GLOBAL_RNG`) random number generator -- `seed` -- (optional, default: `nothing`) seed for reseeding - -### Output - -A random polyhedron. - -### Algorithm - -We first create a random polytope and then for each constraint randomly (50%) -decide whether to include it. -""" -function rand(::Type{HPolyhedron}; - N::Type{<:Real}=Float64, - dim::Int=2, - rng::AbstractRNG=GLOBAL_RNG, - seed::Union{Int,Nothing}=nothing) - rng = reseed!(rng, seed) - P = rand(HPolytope; N=N, dim=dim, rng=rng) - constraints_P = constraints_list(P) - constraints_Q = Vector{eltype(constraints_P)}() - for i in eachindex(constraints_P) - if rand(rng, Bool) - push!(constraints_Q, constraints_P[i]) - end - end - return HPolyhedron(constraints_Q) -end - -""" - addconstraint!(P::HPoly, constraint::HalfSpace) - -Add a linear constraint to a polyhedron in constraint representation. - -### Input - -- `P` -- polyhedron in constraint representation -- `constraint` -- linear constraint to add - -### Notes - -It is left to the user to guarantee that the dimension of all linear constraints -is the same. -""" -function addconstraint!(P::HPoly, constraint::HalfSpace) - push!(P.constraints, constraint) - return nothing -end - -""" - constraints_list(P::HPoly) - -Return the list of constraints defining a polyhedron in constraint -representation. - -### Input - -- `P` -- polyhedron in constraint representation - -### Output - -The list of constraints of the polyhedron. -""" -function constraints_list(P::HPoly) - return P.constraints -end - -""" - tohrep(P::HPoly) - -Return a constraint representation of the given polyhedron in constraint -representation (no-op). - -### Input - -- `P` -- polyhedron in constraint representation - -### Output - -The same polyhedron instance. -""" -function tohrep(P::HPoly) - return P -end - -""" - normalize(P::HPoly{N}, p::Real=N(2)) where {N} - -Normalize a polyhedron in constraint representation. - -### Input - -- `P` -- polyhedron in constraint representation -- `p` -- (optional, default: `2`) norm - -### Output - -A new polyhedron in constraint representation whose normal directions ``a_i`` -are normalized, i.e., such that ``‖a_i‖_p = 1`` holds. -""" -function normalize(P::HPoly{N}, p::Real=N(2)) where {N} - constraints = [normalize(hs, p) for hs in constraints_list(P)] - T = basetype(P) - return T(constraints) -end - -""" - remove_redundant_constraints(P::HPoly{N}; [backend]=nothing) where {N} - -Remove the redundant constraints in a polyhedron in constraint representation. - -### Input - -- `P` -- polyhedron -- `backend` -- (optional, default: `nothing`) the backend used to solve the - linear program - -### Output - -A polyhedron equivalent to `P` but with no redundant constraints, or an empty -set if `P` is detected to be empty, which may happen if the constraints are -infeasible. - -### Notes - -If `backend` is `nothing`, it defaults to `default_lp_solver(N)`. - -### Algorithm - -See `remove_redundant_constraints!(::AbstractVector{<:HalfSpace})` for details. -""" -function remove_redundant_constraints(P::HPoly; backend=nothing) - Pred = copy(P) - if remove_redundant_constraints!(Pred; backend=backend) - return Pred - else # the polyhedron P is empty - require(@__MODULE__, :LazySets; fun_name="remove_redundant_constraints") - - N = eltype(P) - return EmptySet{N}(dim(P)) - end -end - -""" - remove_redundant_constraints!(P::HPoly{N}; [backend]=nothing) where {N} - -Remove the redundant constraints of a polyhedron in constraint representation; -the polyhedron is updated in-place. - -### Input - -- `P` -- polyhedron -- `backend` -- (optional, default: `nothing`) the backend used to solve the - linear program - -### Output - -`true` if the method was successful and the polyhedron `P` is modified by -removing its redundant constraints, and `false` if `P` is detected to be empty, -which may happen if the constraints are infeasible. - -### Notes - -If `backend` is `nothing`, it defaults to `default_lp_solver(N)`. - -### Algorithm - -See `remove_redundant_constraints!(::AbstractVector{<:HalfSpace})` for details. -""" -function remove_redundant_constraints!(P::HPoly; backend=nothing) - return remove_redundant_constraints!(P.constraints; backend=backend) -end - -""" - translate(P::HPoly, v::AbstractVector; [share]::Bool=false) - -Translate (i.e., shift) a polyhedron in constraint representation by a given -vector. - -### Input - -- `P` -- polyhedron in constraint representation -- `v` -- translation vector -- `share` -- (optional, default: `false`) flag for sharing unmodified parts of - the original set representation - -### Output - -A translated polyhedron in constraint representation. - -### Notes - -The normal vectors of the constraints (vector `a` in `a⋅x ≤ b`) are shared with -the original constraints if `share == true`. - -### Algorithm - -We translate every constraint. -""" -function translate(P::HPoly, v::AbstractVector; share::Bool=false) - @assert length(v) == dim(P) "cannot translate a $(dim(P))-dimensional " * - "set by a $(length(v))-dimensional vector" - constraints = [translate(c, v; share=share) for c in constraints_list(P)] - T = basetype(P) - return T(constraints) -end - -""" - tovrep(P::HPoly; [backend]=default_polyhedra_backend(P)) - -Transform a polytope in constraint representation to a polytope in vertex -representation. - -### Input - -- `P` -- polytope in constraint representation -- `backend` -- (optional, default: `default_polyhedra_backend(P)`) the - backend for polyhedral computations - -### Output - -A `VPolytope` which is a vertex representation of the given polytope in -constraint representation. - -### Notes - -The conversion may not preserve the numeric type (e.g., with `N == Float32`) -depending on the backend. -For further information on the supported backends see [Polyhedra's -documentation](https://juliapolyhedra.github.io/). -""" -function tovrep(P::HPoly; backend=default_polyhedra_backend(P)) - require(@__MODULE__, :LazySets; fun_name="tovrep") - require(@__MODULE__, :Polyhedra; fun_name="tovrep") - - P = polyhedron(P; backend=backend) - return VPolytope(P) -end - -# this method is required mainly for HPolytope (because the fallback for -# AbstractPolytope is incorrect with no constraints) -# -# the method also treats a corner case for problems with Rationals in LP solver -function isempty(P::HPoly, - witness::Bool=false; - use_polyhedra_interface::Bool=false, - solver=nothing, - backend=nothing) - if length(constraints_list(P)) < 2 - return witness ? (false, an_element(P)) : false - end - return _isempty_polyhedron(P, witness; - use_polyhedra_interface=use_polyhedra_interface, - solver=solver, backend=backend) -end +include("constraints_list.jl") +include("dim.jl") +include("isempty.jl") +include("isoperationtype.jl") +include("rand.jl") +include("permute.jl") +include("support_function.jl") +include("support_vector.jl") +include("translate.jl") + +include("is_hyperplanar.jl") +include("normalize.jl") +include("remove_redundant_constraints.jl") +include("tohrep.jl") +include("tovrep.jl") +include("addconstraint.jl") function load_polyhedra_hpolyhedron() # function to be loaded by Requires return quote @@ -526,112 +88,6 @@ function load_polyhedra_hpolyhedron() # function to be loaded by Requires end end # quote / load_polyhedra_hpolyhedron() -function is_hyperplanar(P::HPolyhedron) - clist = P.constraints - m = length(clist) - - # check that the number of constraints is fine - if m > 2 - # try to remove redundant constraints - clist = remove_redundant_constraints(clist) - m = length(clist) - end - if m != 2 - return false - end - - # check that the two half-spaces are complementary - return @inbounds iscomplement(clist[1], clist[2]) -end - -function load_symbolics_hpolyhedron() - return quote - using .Symbolics: Num - using ..LazySets: _get_variables, _vec, _is_halfspace, _is_hyperplane - - """ - HPolyhedron(expr::Vector{<:Num}, vars=_get_variables(expr); - N::Type{<:Real}=Float64) - - Return a polyhedron in constraint representation given by a list of symbolic - expressions. - - ### Input - - - `expr` -- vector of symbolic expressions that describes each half-space - - `vars` -- (optional, default: `_get_variables(expr)`), if an array of - variables is given, use those as the ambient variables in the set - with respect to which derivations take place; otherwise, use only - the variables that appear in the given expression (but be careful - because the order may be incorrect; it is advised to always pass - `vars` explicitly) - - `N` -- (optional, default: `Float64`) the numeric type of the returned set - - ### Output - - An `HPolyhedron`. - - ### Examples - - ```jldoctest - julia> using Symbolics - - julia> vars = @variables x y - 2-element Vector{Num}: - x - y - - julia> HPolyhedron([x + y <= 1, x + y >= -1], vars) - HPolyhedron{Float64, Vector{Float64}}(HalfSpace{Float64, Vector{Float64}}[HalfSpace{Float64, Vector{Float64}}([1.0, 1.0], 1.0), HalfSpace{Float64, Vector{Float64}}([-1.0, -1.0], 1.0)]) - - julia> X = HPolyhedron([x == 0, y <= 0], vars) - HPolyhedron{Float64, Vector{Float64}}(HalfSpace{Float64, Vector{Float64}}[HalfSpace{Float64, Vector{Float64}}([1.0, 0.0], -0.0), HalfSpace{Float64, Vector{Float64}}([-1.0, -0.0], 0.0), HalfSpace{Float64, Vector{Float64}}([0.0, 1.0], -0.0)]) - ``` - """ - function HPolyhedron(expr::Vector{<:Num}, vars::AbstractVector{Num}; - N::Type{<:Real}=Float64) - clist = Vector{HalfSpace{N,Vector{N}}}() - sizehint!(clist, length(expr)) - got_hyperplane = false - got_halfspace = false - zeroed_vars = Dict(v => zero(N) for v in vars) - vars_list = collect(vars) - for ex in expr - exval = Symbolics.value(ex) - got_hyperplane, sexpr = _is_hyperplane(exval) - if !got_hyperplane - got_halfspace, sexpr = _is_halfspace(exval) - if !got_halfspace - throw(ArgumentError("expected an expression describing either " * - "a half-space of a hyperplane, got $expr")) - end - end - - coeffs = [N(α.val) for α in gradient(sexpr, vars_list)] - β = -N(Symbolics.substitute(sexpr, zeroed_vars)) - - push!(clist, HalfSpace(coeffs, β)) - if got_hyperplane - push!(clist, HalfSpace(-coeffs, -β)) - end - end - return HPolyhedron(clist) - end - - function HPolyhedron(expr::Vector{<:Num}; N::Type{<:Real}=Float64) - return HPolyhedron(expr, _get_variables(expr); N=N) - end - function HPolyhedron(expr::Vector{<:Num}, vars; N::Type{<:Real}=Float64) - return HPolyhedron(expr, _vec(vars); N=N) - end - end -end # quote / load_symbolics_hpolyhedron() - -function permute(P::HPoly, p::AbstractVector{Int}) - T = basetype(P) - return T([permute(H, p) for H in P.constraints]) -end - include("init.jl") end # module diff --git a/src/Sets/HPolyhedron/addconstraint.jl b/src/Sets/HPolyhedron/addconstraint.jl new file mode 100644 index 0000000000..df4889b5e8 --- /dev/null +++ b/src/Sets/HPolyhedron/addconstraint.jl @@ -0,0 +1,19 @@ +""" + addconstraint!(P::HPoly, constraint::HalfSpace) + +Add a linear constraint to a polyhedron in constraint representation. + +### Input + +- `P` -- polyhedron in constraint representation +- `constraint` -- linear constraint to add + +### Notes + +It is left to the user to guarantee that the dimension of all linear constraints +is the same. +""" +function addconstraint!(P::HPoly, constraint::HalfSpace) + push!(P.constraints, constraint) + return nothing +end diff --git a/src/Sets/HPolyhedron/constraints_list.jl b/src/Sets/HPolyhedron/constraints_list.jl new file mode 100644 index 0000000000..cd7bda9c5a --- /dev/null +++ b/src/Sets/HPolyhedron/constraints_list.jl @@ -0,0 +1,17 @@ +""" + constraints_list(P::HPoly) + +Return the list of constraints defining a polyhedron in constraint +representation. + +### Input + +- `P` -- polyhedron in constraint representation + +### Output + +The list of constraints of the polyhedron. +""" +function constraints_list(P::HPoly) + return P.constraints +end diff --git a/src/Sets/HPolyhedron/dim.jl b/src/Sets/HPolyhedron/dim.jl new file mode 100644 index 0000000000..7aba4be608 --- /dev/null +++ b/src/Sets/HPolyhedron/dim.jl @@ -0,0 +1,17 @@ +""" + dim(P::HPoly) + +Return the dimension of a polyhedron in constraint representation. + +### Input + +- `P` -- polyhedron in constraint representation + +### Output + +The ambient dimension of the polyhedron in constraint representation. +If it has no constraints, the result is ``-1``. +""" +function dim(P::HPoly) + return length(P.constraints) == 0 ? -1 : length(P.constraints[1].a) +end diff --git a/src/Sets/HPolyhedron/is_hyperplanar.jl b/src/Sets/HPolyhedron/is_hyperplanar.jl new file mode 100644 index 0000000000..99607248df --- /dev/null +++ b/src/Sets/HPolyhedron/is_hyperplanar.jl @@ -0,0 +1,17 @@ +function is_hyperplanar(P::HPolyhedron) + clist = P.constraints + m = length(clist) + + # check that the number of constraints is fine + if m > 2 + # try to remove redundant constraints + clist = remove_redundant_constraints(clist) + m = length(clist) + end + if m != 2 + return false + end + + # check that the two half-spaces are complementary + return @inbounds iscomplement(clist[1], clist[2]) +end diff --git a/src/Sets/HPolyhedron/isempty.jl b/src/Sets/HPolyhedron/isempty.jl new file mode 100644 index 0000000000..a02e2595b1 --- /dev/null +++ b/src/Sets/HPolyhedron/isempty.jl @@ -0,0 +1,16 @@ +# this method is required mainly for HPolytope (because the fallback for +# AbstractPolytope is incorrect with no constraints) +# +# the method also treats a corner case for problems with Rationals in LP solver +function isempty(P::HPoly, + witness::Bool=false; + use_polyhedra_interface::Bool=false, + solver=nothing, + backend=nothing) + if length(constraints_list(P)) < 2 + return witness ? (false, an_element(P)) : false + end + return _isempty_polyhedron(P, witness; + use_polyhedra_interface=use_polyhedra_interface, + solver=solver, backend=backend) +end diff --git a/src/Sets/HPolyhedron/isoperationtype.jl b/src/Sets/HPolyhedron/isoperationtype.jl new file mode 100644 index 0000000000..c37399d223 --- /dev/null +++ b/src/Sets/HPolyhedron/isoperationtype.jl @@ -0,0 +1 @@ +isoperationtype(::Type{<:HPolyhedron}) = false diff --git a/src/Sets/HPolyhedron/normalize.jl b/src/Sets/HPolyhedron/normalize.jl new file mode 100644 index 0000000000..2d2c1781ed --- /dev/null +++ b/src/Sets/HPolyhedron/normalize.jl @@ -0,0 +1,20 @@ +""" + normalize(P::HPoly{N}, p::Real=N(2)) where {N} + +Normalize a polyhedron in constraint representation. + +### Input + +- `P` -- polyhedron in constraint representation +- `p` -- (optional, default: `2`) norm + +### Output + +A new polyhedron in constraint representation whose normal directions ``a_i`` +are normalized, i.e., such that ``‖a_i‖_p = 1`` holds. +""" +function normalize(P::HPoly{N}, p::Real=N(2)) where {N} + constraints = [normalize(hs, p) for hs in constraints_list(P)] + T = basetype(P) + return T(constraints) +end diff --git a/src/Sets/HPolyhedron/permute.jl b/src/Sets/HPolyhedron/permute.jl new file mode 100644 index 0000000000..e5fce92fb6 --- /dev/null +++ b/src/Sets/HPolyhedron/permute.jl @@ -0,0 +1,4 @@ +function permute(P::HPoly, p::AbstractVector{Int}) + T = basetype(P) + return T([permute(H, p) for H in P.constraints]) +end diff --git a/src/Sets/HPolyhedron/rand.jl b/src/Sets/HPolyhedron/rand.jl new file mode 100644 index 0000000000..73d12e7088 --- /dev/null +++ b/src/Sets/HPolyhedron/rand.jl @@ -0,0 +1,39 @@ +""" + rand(::Type{HPolyhedron}; [N]::Type{<:Real}=Float64, [dim]::Int=2, + [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing) + +Create a random polyhedron. + +### Input + +- `HPolyhedron` -- type for dispatch +- `N` -- (optional, default: `Float64`) numeric type +- `dim` -- (optional, default: 2) dimension (is ignored) +- `rng` -- (optional, default: `GLOBAL_RNG`) random number generator +- `seed` -- (optional, default: `nothing`) seed for reseeding + +### Output + +A random polyhedron. + +### Algorithm + +We first create a random polytope and then for each constraint randomly (50%) +decide whether to include it. +""" +function rand(::Type{HPolyhedron}; + N::Type{<:Real}=Float64, + dim::Int=2, + rng::AbstractRNG=GLOBAL_RNG, + seed::Union{Int,Nothing}=nothing) + rng = reseed!(rng, seed) + P = rand(HPolytope; N=N, dim=dim, rng=rng) + constraints_P = constraints_list(P) + constraints_Q = Vector{eltype(constraints_P)}() + for i in eachindex(constraints_P) + if rand(rng, Bool) + push!(constraints_Q, constraints_P[i]) + end + end + return HPolyhedron(constraints_Q) +end diff --git a/src/Sets/HPolyhedron/remove_redundant_constraints.jl b/src/Sets/HPolyhedron/remove_redundant_constraints.jl new file mode 100644 index 0000000000..35dd8da0cd --- /dev/null +++ b/src/Sets/HPolyhedron/remove_redundant_constraints.jl @@ -0,0 +1,66 @@ +""" + remove_redundant_constraints(P::HPoly{N}; [backend]=nothing) where {N} + +Remove the redundant constraints in a polyhedron in constraint representation. + +### Input + +- `P` -- polyhedron +- `backend` -- (optional, default: `nothing`) the backend used to solve the + linear program + +### Output + +A polyhedron equivalent to `P` but with no redundant constraints, or an empty +set if `P` is detected to be empty, which may happen if the constraints are +infeasible. + +### Notes + +If `backend` is `nothing`, it defaults to `default_lp_solver(N)`. + +### Algorithm + +See `remove_redundant_constraints!(::AbstractVector{<:HalfSpace})` for details. +""" +function remove_redundant_constraints(P::HPoly; backend=nothing) + Pred = copy(P) + if remove_redundant_constraints!(Pred; backend=backend) + return Pred + else # the polyhedron P is empty + require(@__MODULE__, :LazySets; fun_name="remove_redundant_constraints") + + N = eltype(P) + return EmptySet{N}(dim(P)) + end +end + +""" + remove_redundant_constraints!(P::HPoly{N}; [backend]=nothing) where {N} + +Remove the redundant constraints of a polyhedron in constraint representation; +the polyhedron is updated in-place. + +### Input + +- `P` -- polyhedron +- `backend` -- (optional, default: `nothing`) the backend used to solve the + linear program + +### Output + +`true` if the method was successful and the polyhedron `P` is modified by +removing its redundant constraints, and `false` if `P` is detected to be empty, +which may happen if the constraints are infeasible. + +### Notes + +If `backend` is `nothing`, it defaults to `default_lp_solver(N)`. + +### Algorithm + +See `remove_redundant_constraints!(::AbstractVector{<:HalfSpace})` for details. +""" +function remove_redundant_constraints!(P::HPoly; backend=nothing) + return remove_redundant_constraints!(P.constraints; backend=backend) +end diff --git a/src/Sets/HPolyhedron/support_function.jl b/src/Sets/HPolyhedron/support_function.jl new file mode 100644 index 0000000000..914df40189 --- /dev/null +++ b/src/Sets/HPolyhedron/support_function.jl @@ -0,0 +1,32 @@ +""" + ρ(d::AbstractVector{M}, P::HPoly{N}; + solver=default_lp_solver(M, N)) where {M, N} + +Evaluate the support function of a polyhedron in constraint representation in a +given direction. + +### Input + +- `d` -- direction +- `P` -- polyhedron in constraint representation +- `solver` -- (optional, default: `default_lp_solver(M, N)`) the backend used to + solve the linear program + +### Output + +The evaluation of the support function for the polyhedron. +If a polytope is unbounded in the given direction, we throw an error. +If a polyhedron is unbounded in the given direction, the result is `Inf`. +""" +function ρ(d::AbstractVector{M}, P::HPoly{N}; + solver=default_lp_solver(M, N)) where {M,N} + lp, unbounded = σ_helper(d, P, solver) + if unbounded + if P isa HPolytope + error("the support function in direction $(d) is undefined " * + "because the polytope is unbounded") + end + return N(Inf) + end + return dot(d, lp.sol) +end diff --git a/src/Sets/HPolyhedron/support_vector.jl b/src/Sets/HPolyhedron/support_vector.jl new file mode 100644 index 0000000000..7fcb5c75dd --- /dev/null +++ b/src/Sets/HPolyhedron/support_vector.jl @@ -0,0 +1,86 @@ +""" + σ(d::AbstractVector{M}, P::HPoly{N}; + solver=default_lp_solver(M, N) where {M, N} + +Return a support vector of a polyhedron in constraint representation in a given +direction. + +### Input + +- `d` -- direction +- `P` -- polyhedron in constraint representation +- `solver` -- (optional, default: `default_lp_solver(M, N)`) the backend used to + solve the linear program + +### Output + +The support vector in the given direction. +If a polytope is unbounded in the given direction, we throw an error. +If a polyhedron is unbounded in the given direction, the result contains `±Inf` +entries. +""" +function σ(d::AbstractVector{M}, P::HPoly{N}; + solver=default_lp_solver(M, N)) where {M,N} + lp, unbounded = σ_helper(d, P, solver) + if unbounded + if P isa HPolytope + error("the support vector in direction $(d) is undefined because " * + "the polytope is unbounded") + end + return _σ_unbounded_lp(d, P, lp) + else + return lp.sol + end +end + +# construct the solution from the solver's ray result +function _σ_unbounded_lp(d, P::HPoly{N}, lp) where {N} + if isnothing(lp) + ray = d + elseif has_lp_infeasibility_ray(lp.model) + ray = lp.sol # infeasibility ray is stored as the solution + else + error("LP solver did not return an infeasibility ray") + end + + res = Vector{N}(undef, length(ray)) + e = isempty(P.constraints) ? zeros(N, length(ray)) : an_element(P) + @inbounds for i in eachindex(ray) + if isapproxzero(ray[i]) + res[i] = e[i] + elseif ray[i] > zero(N) + res[i] = N(Inf) + else + res[i] = N(-Inf) + end + end + return res +end + +function σ_helper(d::AbstractVector, P::HPoly, solver) + # represent c = -d as a Vector since GLPK does not accept sparse vectors + # (see #1011) + c = to_negative_vector(d) + + (A, b) = tosimplehrep(P) + if length(b) == 0 + unbounded = true + lp = nothing + else + sense = '<' + l = -Inf + u = Inf + lp = linprog(c, A, sense, b, l, u, solver) + if is_lp_infeasible(lp.status; strict=true) + throw(ArgumentError("the support vector is undefined because " * + "the polyhedron is empty")) + elseif is_lp_unbounded(lp.status) + unbounded = true + elseif is_lp_optimal(lp.status) + unbounded = false + else + error("got unknown LP status $(lp.status)") + end + end + return (lp, unbounded) +end diff --git a/src/Sets/HPolyhedron/tohrep.jl b/src/Sets/HPolyhedron/tohrep.jl new file mode 100644 index 0000000000..f86444d269 --- /dev/null +++ b/src/Sets/HPolyhedron/tohrep.jl @@ -0,0 +1,17 @@ +""" + tohrep(P::HPoly) + +Return a constraint representation of the given polyhedron in constraint +representation (no-op). + +### Input + +- `P` -- polyhedron in constraint representation + +### Output + +The same polyhedron instance. +""" +function tohrep(P::HPoly) + return P +end diff --git a/src/Sets/HPolyhedron/tovrep.jl b/src/Sets/HPolyhedron/tovrep.jl new file mode 100644 index 0000000000..75bcdd41f7 --- /dev/null +++ b/src/Sets/HPolyhedron/tovrep.jl @@ -0,0 +1,31 @@ +""" + tovrep(P::HPoly; [backend]=default_polyhedra_backend(P)) + +Transform a polytope in constraint representation to a polytope in vertex +representation. + +### Input + +- `P` -- polytope in constraint representation +- `backend` -- (optional, default: `default_polyhedra_backend(P)`) the + backend for polyhedral computations + +### Output + +A `VPolytope` which is a vertex representation of the given polytope in +constraint representation. + +### Notes + +The conversion may not preserve the numeric type (e.g., with `N == Float32`) +depending on the backend. +For further information on the supported backends see [Polyhedra's +documentation](https://juliapolyhedra.github.io/). +""" +function tovrep(P::HPoly; backend=default_polyhedra_backend(P)) + require(@__MODULE__, :LazySets; fun_name="tovrep") + require(@__MODULE__, :Polyhedra; fun_name="tovrep") + + P = polyhedron(P; backend=backend) + return VPolytope(P) +end diff --git a/src/Sets/HPolyhedron/translate.jl b/src/Sets/HPolyhedron/translate.jl new file mode 100644 index 0000000000..50543c0e5c --- /dev/null +++ b/src/Sets/HPolyhedron/translate.jl @@ -0,0 +1,33 @@ +""" + translate(P::HPoly, v::AbstractVector; [share]::Bool=false) + +Translate (i.e., shift) a polyhedron in constraint representation by a given +vector. + +### Input + +- `P` -- polyhedron in constraint representation +- `v` -- translation vector +- `share` -- (optional, default: `false`) flag for sharing unmodified parts of + the original set representation + +### Output + +A translated polyhedron in constraint representation. + +### Notes + +The normal vectors of the constraints (vector `a` in `a⋅x ≤ b`) are shared with +the original constraints if `share == true`. + +### Algorithm + +We translate every constraint. +""" +function translate(P::HPoly, v::AbstractVector; share::Bool=false) + @assert length(v) == dim(P) "cannot translate a $(dim(P))-dimensional " * + "set by a $(length(v))-dimensional vector" + constraints = [translate(c, v; share=share) for c in constraints_list(P)] + T = basetype(P) + return T(constraints) +end From eaca6d189251658fa7777ec0fa3134a7f611c9ca Mon Sep 17 00:00:00 2001 From: schillic Date: Tue, 16 Jul 2024 13:58:15 +0200 Subject: [PATCH 5/6] add tests --- test/Sets/Polyhedron.jl | 29 ++++++++++++++++++++++++++++- 1 file changed, 28 insertions(+), 1 deletion(-) diff --git a/test/Sets/Polyhedron.jl b/test/Sets/Polyhedron.jl index 166ba7b4d6..c06ec40c0d 100644 --- a/test/Sets/Polyhedron.jl +++ b/test/Sets/Polyhedron.jl @@ -6,6 +6,8 @@ for N in [Float64, Rational{Int}, Float32] # random polyhedron rand(HPolyhedron) + @test HPolyhedron{N}() isa HPolyhedron{N} + @test HPolyhedron{N,Vector{N}}() isa HPolyhedron{N,Vector{N}} p_univ = HPolyhedron{N}() # constructor from matrix and vector @@ -45,6 +47,14 @@ for N in [Float64, Rational{Int}, Float32] # is_polyhedral @test is_polyhedral(p) + # is_hyperplanar + P = HPolyhedron([HalfSpace(N[1, 0], N(1))]) + @test !is_hyperplanar(P) + addconstraint!(P, HalfSpace(N[-1, 0], N(-1))) + @test is_hyperplanar(P) + addconstraint!(P, HalfSpace(N[0, 1], N(1))) + @test !is_hyperplanar(P) + # universality @test !isuniversal(p) res, w = isuniversal(p, true) @@ -80,9 +90,12 @@ for N in [Float64, Rational{Int}, Float32] @test isequivalent(p1, p1) && !isequivalent(p1, p2) if test_suite_polyhedra - # conversion to and from Polyhedra's VRep data structure + # conversion to and from Polyhedra's HRep data structure cl = constraints_list(HPolyhedron(polyhedron(p))) @test length(p.constraints) == length(cl) + # filtering of trivial constraint + P = Polyhedra.HalfSpace(N[0], N(1)) ∩ Polyhedra.HalfSpace(N[1], N(1)) + @test HPolyhedron(P).constraints == [HalfSpace(N[1], N(1))] # convert hyperrectangle to a HPolyhedron H = Hyperrectangle(N[1, 1], N[2, 2]) @@ -214,6 +227,10 @@ for N in [Float64] @test σ(d, p) == N[-1, 1] d = N[0, -1] @test σ(d, p) == N[0, 0] + # unbounded set + P = HPolyhedron([c1, c2, c3]) + d = N[1, 0] + @test σ(d, P) == N[Inf, -Inf] # membership @test [Inf, Inf] ∉ p @@ -332,6 +349,10 @@ for N in [Float64] P = HPolyhedron([HalfSpace(N[1], N(0)), HalfSpace(N[1], N(1))]) Q = remove_redundant_constraints(P) @test length(P.constraints) == 2 && length(Q.constraints) == 1 + # removing redundant constraints from an empty polyhedron + P = HPolyhedron([HalfSpace(N[1], N(0)), HalfSpace(N[-1], N(-1)), HalfSpace(N[-1], N(-1))]) + Q = remove_redundant_constraints(P) + @test Q isa EmptySet{N} && dim(Q) == 1 # checking for empty intersection (also test symmetric calls) P = convert(HPolytope, BallInf(zeros(N, 2), N(1))) @@ -433,5 +454,11 @@ for N in [Float64] h2 = HPolyhedron([HalfSpace([1.0, 0.0], 0.0), HalfSpace([-1.0, 0.0], 0.0), HalfSpace([0.0, 1.0], 0.0)]) @test p2 ⊆ h2 && h2 ⊆ p2 # isequivalent(p2, h2) see #2370 + + # invalid expression + @test_throws ArgumentError HPolyhedron([x + y != 1], vars) end end + +# isoperationtype +@test !isoperationtype(HPolyhedron) From a195a6ad9d6dd7b584f4f6404ba12f647ed7b9c1 Mon Sep 17 00:00:00 2001 From: schillic Date: Tue, 16 Jul 2024 21:11:58 +0200 Subject: [PATCH 6/6] prepend 'gradient' with 'Symbolics.' --- src/Initialization/init_Symbolics.jl | 3 +-- src/Sets/HPolyhedron/HPolyhedron.jl | 2 +- src/Sets/HalfSpace.jl | 2 +- src/Sets/Hyperplane.jl | 2 +- 4 files changed, 4 insertions(+), 5 deletions(-) diff --git a/src/Initialization/init_Symbolics.jl b/src/Initialization/init_Symbolics.jl index 7daf003e7b..bfd23ea741 100644 --- a/src/Initialization/init_Symbolics.jl +++ b/src/Initialization/init_Symbolics.jl @@ -1,5 +1,4 @@ -using .Symbolics: gradient, - simplify, +using .Symbolics: simplify, Num, # variable like, e.g. x[1] Term, # term like, eg. x[1] + x[2] == 1 Symbolic, diff --git a/src/Sets/HPolyhedron/HPolyhedron.jl b/src/Sets/HPolyhedron/HPolyhedron.jl index f716637a0e..2f9e565e92 100644 --- a/src/Sets/HPolyhedron/HPolyhedron.jl +++ b/src/Sets/HPolyhedron/HPolyhedron.jl @@ -111,7 +111,7 @@ function load_symbolics_hpolyhedron() end end - coeffs = [N(α.val) for α in gradient(sexpr, vars_list)] + coeffs = [N(α.val) for α in Symbolics.gradient(sexpr, vars_list)] β = -N(Symbolics.substitute(sexpr, zeroed_vars)) push!(clist, HalfSpace(coeffs, β)) diff --git a/src/Sets/HalfSpace.jl b/src/Sets/HalfSpace.jl index b0d6b6c614..404dec5559 100644 --- a/src/Sets/HalfSpace.jl +++ b/src/Sets/HalfSpace.jl @@ -627,7 +627,7 @@ function load_symbolics_halfspace() end # compute the linear coefficients by taking first-order derivatives - coeffs = [N(α.val) for α in gradient(sexpr, collect(vars))] + coeffs = [N(α.val) for α in Symbolics.gradient(sexpr, collect(vars))] # get the constant term by expression substitution zeroed_vars = Dict(v => zero(N) for v in vars) diff --git a/src/Sets/Hyperplane.jl b/src/Sets/Hyperplane.jl index cd1b713e02..e87285d782 100644 --- a/src/Sets/Hyperplane.jl +++ b/src/Sets/Hyperplane.jl @@ -596,7 +596,7 @@ function load_symbolics_hyperplane() end # compute the linear coefficients by taking first order derivatives - coeffs = [N(α.val) for α in gradient(sexpr, collect(vars))] + coeffs = [N(α.val) for α in Symbolics.gradient(sexpr, collect(vars))] # get the constant term by expression substitution zeroed_vars = Dict(v => zero(N) for v in vars)