Merge pull request #1053 from JuliaReach/mforets/1051

#1051 - Add more list of constraints methods

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

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))

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}

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")

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}}

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
 # ==================================

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))

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
     # ================