diff --git a/docs/src/lib/interfaces.md b/docs/src/lib/interfaces.md index 92f14364f8..2138fd7c68 100644 --- a/docs/src/lib/interfaces.md +++ b/docs/src/lib/interfaces.md @@ -63,6 +63,7 @@ RecipesBase.apply_recipe(::Dict{Symbol,Any}, ::LazySet) RecipesBase.apply_recipe(::Dict{Symbol,Any}, ::Vector{S}) where {S<:LazySet} RecipesBase.apply_recipe(::Dict{Symbol,Any}, ::LazySet, ::Float64) RecipesBase.apply_recipe(::Dict{Symbol,Any}, ::Vector{S}, ::Float64) where {S<:LazySet} +tosimplehrep(::LazySet) ``` ### Set functions that override Base functions diff --git a/docs/src/lib/operations.md b/docs/src/lib/operations.md index e2b94b6035..cbad8f9823 100644 --- a/docs/src/lib/operations.md +++ b/docs/src/lib/operations.md @@ -120,6 +120,7 @@ dim(::Intersection) isbounded(::Intersection) isempty(::Intersection) ∈(::AbstractVector{N}, ::Intersection{N}) where {N<:Real} +constraints_list(::Intersection{N}) where {N<:Real} isempty_known(::Intersection) set_isempty!(::Intersection, ::Bool) swap(::Intersection) @@ -149,6 +150,7 @@ dim(::IntersectionArray) isbounded(::IntersectionArray) ∈(::AbstractVector{N}, ::IntersectionArray{N}) where {N<:Real} array(::IntersectionArray{N, S}) where {N<:Real, S<:LazySet{N}} +constraints_list(::IntersectionArray{N}) where {N<:Real} ``` Inherited from [`LazySet`](@ref): * [`norm`](@ref norm(::LazySet, ::Real)) diff --git a/docs/src/lib/representations.md b/docs/src/lib/representations.md index 415699181e..0bab47de04 100644 --- a/docs/src/lib/representations.md +++ b/docs/src/lib/representations.md @@ -173,6 +173,7 @@ constraints_list(::HalfSpace{N}) where {N<:Real} constrained_dimensions(::HalfSpace{N}) where {N<:Real} halfspace_left(::AbstractVector{N}, ::AbstractVector{N}) where {N<:Real} halfspace_right(::AbstractVector{N}, ::AbstractVector{N}) where {N<:Real} +tosimplehrep(::AbstractVector{HalfSpace{N}}) where {N<:Real} ``` Inherited from [`LazySet`](@ref): * [`norm`](@ref norm(::LazySet, ::Real)) @@ -437,7 +438,6 @@ dim(::HPoly{N}) where {N<:Real} ∈(::AbstractVector{N}, ::HPoly{N}) where {N<:Real} addconstraint!(::HPoly{N}, ::LinearConstraint{N}) where {N<:Real} constraints_list(::HPoly{N}) where {N<:Real} -tosimplehrep(::HPoly{N}) where {N<:Real} tohrep(::HPoly{N}) where {N<:Real} isempty(::HPoly{N}, ::Bool=false) where {N<:Real} cartesian_product(::HPoly{N}, ::HPoly{N}) where {N<:Real} diff --git a/src/AbstractPolytope.jl b/src/AbstractPolytope.jl index ddb29eb735..bbc074597c 100644 --- a/src/AbstractPolytope.jl +++ b/src/AbstractPolytope.jl @@ -124,11 +124,6 @@ function isempty(P::AbstractPolytope)::Bool return isempty(vertices_list(P)) end -# TODO simplify this -function tosimplehrep(P::LazySet) - return tosimplehrep(HPolyhedron(constraints_list(P))) -end - function default_polyhedra_backend(P, N) @assert isdefined(@__MODULE__, :Polyhedra) "this function needs the package 'Polyhedra' to be loaded" error("no default backend for numeric type $N") diff --git a/src/HPolyhedron.jl b/src/HPolyhedron.jl index b375dba3c2..e583ed2e3a 100644 --- a/src/HPolyhedron.jl +++ b/src/HPolyhedron.jl @@ -8,7 +8,6 @@ export HPolyhedron, dim, σ, ∈, addconstraint!, constraints_list, - tosimplehrep, tohrep, tovrep, convex_hull, cartesian_product, @@ -369,31 +368,33 @@ function constraints_list(P::HPoly{N} end """ - tosimplehrep(P::HPoly{N}) where {N} + tosimplehrep(constraints::AbstractVector{LinearConstraint{N}}) where {N<:Real} -Return the simple H-representation ``Ax ≤ b`` of a polyhedron. +Return the simple H-representation ``Ax ≤ b`` of a list of constraints. ### Input -- `P` -- polyhedron +- `constraints` -- a list of constraints ### Output The tuple `(A, b)` where `A` is the matrix of normal directions and `b` are the offsets. """ -function tosimplehrep(P::HPoly{N}) where {N<:Real} - n = length(constraints_list(P)) +function tosimplehrep(constraints::AbstractVector{LinearConstraint{N}}) where {N<:Real} + n = length(constraints) if n == 0 A = Matrix{N}(undef, 0, 0) b = Vector{N}(undef, 0) return (A, b) end - A = zeros(N, n, dim(P)) + A = zeros(N, n, dim(first(constraints))) b = zeros(N, n) - for (i, Pi) in enumerate(constraints_list(P)) - A[i, :] = Pi.a - b[i] = Pi.b + @inbounds begin + for (i, Pi) in enumerate(constraints) + A[i, :] = Pi.a + b[i] = Pi.b + end end return (A, b) end @@ -430,11 +431,13 @@ Remove the redundant constraints in a polyhedron in H-representation. ### Output A polyhedron equivalent to `P` but with no redundant constraints, or an empty set -if `P` is detected to be empty (which happens if the constraints are infeasible). +if `P` is detected to be empty, which may happen if the constraints are infeasible. ### Algorithm -See [`remove_redundant_constraints!`](@ref) for details. +See +[`remove_redundant_constraints!(::Vector{LinearConstraint{N}}) where {N<:Real}`](@ref) +for details. """ function remove_redundant_constraints(P::PT; backend=GLPKSolverLP())::Union{PT, EmptySet{N}} where {N, PT<:HPoly{N}} @@ -461,23 +464,59 @@ is updated in-place. ### 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 happens if the constraints are infeasible). +removing its redundant constraints, and `false` if `P` is detected to be empty, +which may happen if the constraints are infeasible. + +### Algorithm + +See +[`remove_redundant_constraints!(::Vector{LinearConstraint{N}}) where {N<:Real}`](@ref) +for details. +""" +function remove_redundant_constraints!(P::PT; + backend=GLPKSolverLP())::Bool where {N, PT<:HPoly{N}} + remove_redundant_constraints!(P.constraints, backend=backend) +end + +""" + remove_redundant_constraints!(constraints::Vector{LinearConstraint{N}}; + backend=GLPKSolverLP())::Bool where {N} + +Remove the redundant constraints of a given list of linear constraints; the list +is updated in-place. + +### Input + +- `constraints` -- list of constraints +- `backend` -- (optional, default: `GLPKSolverLP`) the numeric LP solver backend + +### Output + +`true` if the method was successful and the list of constraints `constraints` is +modified by removing the redundant constraints, and `false` if the constraints +are infeasible. ### Algorithm -If the polyhedron `P` has `m` constraints and its dimension is `n`, -this function checks one by one if each of the `m` constraints is -implied by the remaining ones. To check if the `k`-th constraint -is redundant, an LP is formulated. +If there are `m` constraints in `n` dimensions, this function checks one by one +if each of the `m` constraints is implied by the remaining ones. + +To check if the `k`-th constraint is redundant, an LP is formulated using the +constraints that have not yet being removed. If, at an intermediate step, it is +detected that a subgruop of the constraints is infeasible, this function returns +`false`; if the calculation finished successfully it returns `true`. + +Note that the constraints being infeasible does not imply that `false` is returned. +For example, `x <= 0 && x >= 1` will return `true` without removing any constraint. +To check if the constraints are infeasible use `isempty(HPolyhedron(constraints)`. For details, see [Fukuda's Polyhedra FAQ](https://www.cs.mcgill.ca/~fukuda/soft/polyfaq/node24.html). """ -function remove_redundant_constraints!(P::PT; - backend=GLPKSolverLP())::Bool where {N, PT<:HPoly{N}} +function remove_redundant_constraints!(constraints::AbstractVector{LinearConstraint{N}}; + backend=GLPKSolverLP())::Bool where {N} - A, b = tosimplehrep(P) + A, b = tosimplehrep(constraints) m, n = size(A) non_redundant_indices = 1:m @@ -506,10 +545,42 @@ function remove_redundant_constraints!(P::PT; end end - deleteat!(P.constraints, setdiff(1:m, non_redundant_indices)) + deleteat!(constraints, setdiff(1:m, non_redundant_indices)) return true end +""" + remove_redundant_constraints(constraints::AbstractVector{LinearConstraint{N}}; + backend=GLPKSolverLP())::Union{AbstractVector{LinearConstraint{N}}, EmptySet{N}} where {N} + +Remove the redundant constraints of a given list of linear constraints. + +### Input + +- `constraints` -- list of constraints +- `backend` -- (optional, default: `GLPKSolverLP`) the numeric LP solver backend + +### Output + +The list of constraints with the redundant ones removed, or an empty set +if the given constraints are infeasible. + +### Algorithm + +See +[`remove_redundant_constraints!(::AbstractVector{LinearConstraint{N}}) where {N<:Real}`](@ref) +for details. +""" +function remove_redundant_constraints(constraints::AbstractVector{LinearConstraint{N}}; + backend=GLPKSolverLP())::Union{AbstractVector{LinearConstraint{N}}, EmptySet{N}} where {N} + constraints_copy = copy(constraints) + if remove_redundant_constraints!(constraints_copy, backend=backend) + return constraints_copy + else # the constraints are infeasible + return EmptySet{N}() + end +end + """ linear_map(M::AbstractMatrix{N}, P::PT; [cond_tol=DEFAULT_COND_TOL]::Number) where {N<:Real, PT<:HPoly{N}} diff --git a/src/Intersection.jl b/src/Intersection.jl index 794c256b1b..675c666030 100644 --- a/src/Intersection.jl +++ b/src/Intersection.jl @@ -6,7 +6,8 @@ export Intersection, swap, use_precise_ρ, IntersectionArray, - array + array, + constraints_list """ IntersectionCache @@ -510,6 +511,39 @@ function ∈(x::AbstractVector{N}, cap::Intersection{N})::Bool where {N<:Real} return (x ∈ cap.X) && (x ∈ cap.Y) end +""" + constraints_list(cap::Intersection{N}) where {N<:Real} + +Return the list of constraints of an intersection of two (polyhedral) sets. + +### Input + +- `cap` -- intersection of two (polyhedral) sets + +### Output + +The list of constraints of the intersection. + +### Notes + +We assume that the underlying sets are polyhedral, i.e., offer a method +`constraints_list`. + +### Algorithm + +We create the polyhedron by taking the intersection of the `constraints_list`s of +the sets and remove redundant constraints. + +This function ignores the boolean output from the in-place `remove_redundant_constraints!`, +which may inform the user that the constraints are infeasible. In that case, the +list of constraints at the moment when the infeasibility was detected is returned. +""" +function constraints_list(cap::Intersection{N}) where {N<:Real} + constraints = [constraints_list(cap.X); constraints_list(cap.Y)] + remove_redundant_constraints!(constraints) + return constraints +end + # --- Intersection functions --- @@ -704,6 +738,39 @@ function ∈(x::AbstractVector{N}, ia::IntersectionArray{N})::Bool where {N<:Rea return true end +""" + constraints_list(ia::IntersectionArray{N}) where {N<:Real} + +Return the list of constraints of an intersection of a finite number of +(polyhedral) sets. + +### Input + +- `ia` -- intersection of a finite number of (polyhedral) sets + +### Output + +The list of constraints of the intersection. + +### Notes + +We assume that the underlying sets are polyhedral, i.e., offer a method +`constraints_list`. + +### Algorithm + +We create the polyhedron from the `constraints_list`s of the sets and remove +redundant constraints. +""" +function constraints_list(ia::IntersectionArray{N}) where {N<:Real} + constraints = Vector{LinearConstraint{N}}() + for X in array(ia) + append!(constraints, constraints_list(X)) + end + remove_redundant_constraints!(constraints) + return constraints +end + # ================================== # Algorithms for lazy intersection # ================================== diff --git a/src/LazySet.jl b/src/LazySet.jl index 061cd4b6c9..cb0ba26b72 100644 --- a/src/LazySet.jl +++ b/src/LazySet.jl @@ -10,7 +10,8 @@ export LazySet, an_element, isbounded, isbounded_unit_dimensions, neutral, - absorbing + absorbing, + tosimplehrep """ LazySet{N} @@ -320,3 +321,27 @@ completely independent object. See the documentation of `?deepcopy` for further details. """ copy(S::LazySet) = deepcopy(S) + +""" + tosimplehrep(S::LazySet) + +Return the simple H-representation ``Ax ≤ b`` of a set from its list of +constraints. + +### Input + +- `S` -- set + +### Output + +The tuple `(A, b)` where `A` is the matrix of normal directions and `b` are the +offsets. + +### Notes + +This function uses `constraints_list(S)`. It is a fallback implementation that +works only for those sets that can be represented exactly by a list of linear +constraints, which is available through the `constraints_list(S)` +function. +""" +tosimplehrep(S::LazySet) = tosimplehrep(constraints_list(S)) diff --git a/test/unit_Intersection.jl b/test/unit_Intersection.jl index ba9528f6c8..efc314b1da 100644 --- a/test/unit_Intersection.jl +++ b/test/unit_Intersection.jl @@ -27,6 +27,13 @@ for N in [Float64, Rational{Int}, Float32] @test isempty_known(I) @test !isempty(I) + # constraints_list for polytopic intersection + @test ispermutation(constraints_list(I), + [HalfSpace(N[1, 0], N(2)), + HalfSpace(N[0, 1], N(2)), + HalfSpace(N[-1, 0], N(0)), + HalfSpace(N[0, -1], N(0))]) + # ================= # IntersectionArray # ================= @@ -64,6 +71,13 @@ for N in [Float64, Rational{Int}, Float32] # constructor with size hint and type IntersectionArray(10, N) + # constraints_list for polytopic intersection + @test ispermutation(constraints_list(IA), + [HalfSpace(N[1, 0], N(2)), + HalfSpace(N[0, 1], N(2)), + HalfSpace(N[-1, 0], N(0)), + HalfSpace(N[0, -1], N(0))]) + # ================ # common functions # ================ diff --git a/test/unit_Polyhedron.jl b/test/unit_Polyhedron.jl index ef7b221fc6..68f2de0101 100644 --- a/test/unit_Polyhedron.jl +++ b/test/unit_Polyhedron.jl @@ -28,6 +28,13 @@ for N in [Float64, Rational{Int}, Float32] addconstraint!(p, c4) addconstraint!(p, c2) + # tosimplehrep for other polytopic types + A, b = tosimplehrep(BallInf(N[0], N(2))) + @test (A == hcat(N[1; -1]) || A == hcat(N[-1; 1])) && b == N[2, 2] + # tosimplehrep from list of constraints + A, b = tosimplehrep([HalfSpace(N[1], N(2)), HalfSpace(N[-1], N(2))]) + @test (A == hcat(N[1; -1]) || A == hcat(N[-1; 1])) && b == N[2, 2] + # support vector d = N[1, 0] @test σ(d, p) == N[4, 2] @@ -154,6 +161,16 @@ if test_suite_polyhedra Ar, br = tosimplehrep(p1) @test Ar == A[1:2, :] && br == b[1:2] + # removing redundant constraints from a list of constraints + constraints = [HalfSpace(N[1], N(1)), HalfSpace(N[1], N(0))] + constraints2 = remove_redundant_constraints(constraints) + result = remove_redundant_constraints!(constraints) + @test result + @test ispermutation(constraints, constraints2) + constraints = [HalfSpace(N[1], N(0)), HalfSpace(N[-1], N(-1)), HalfSpace(N[-1], N(-1))] + result = remove_redundant_constraints!(constraints) + @test !result + # checking for empty intersection (also test symmetric calls) P = convert(HPolytope, BallInf(zeros(N, 2), N(1))) Q = convert(HPolytope, BallInf(ones(N, 2), N(1)))