Empty polyhedra intersection is broken

[schillic]

julia> using LazySets, Polyhedra

julia> H = LazySets.HalfSpace([-1.0], -21.0); # x >= 21

julia> P = HPolytope([LazySets.HalfSpace([1.0], 18.1), # x <= 18.1
                      LazySets.HalfSpace([-1.0], -18.0)]); # x >= 18

In v0.6 this results in:

julia> cap = intersection(H, P)
LazySets.HPolytope{Float64}(LazySets.HalfSpace{Float64}[
 LazySets.HalfSpace{Float64}([1.0], 18.1),
 LazySets.HalfSpace{Float64}([-1.0], -18.0),
 LazySets.HalfSpace{Float64}([-1.0], -18.1),
 LazySets.HalfSpace{Float64}([1.0], 18.0)])

julia> isempty(cap)
true

In v0.7 this results in:

julia> cap = intersection(H, P)
HPolytope{Float64}(LazySets.HalfSpace{Float64}[])

julia> isempty(cap)
false

Problems:

1. Why does v0.6 not simplify?
   It might be because the simplification with CDDLib is hard to interpret (contains two hyperplanes):

julia> Q = polyhedron(cap, CDDLib.CDDLibrary())
Polyhedron CDDLib.CDDPolyhedron{1,Float64}:
4-element iterator of Polyhedra.HalfSpace{1,Float64,Array{Float64,1}}:
 HalfSpace([1.0], 18.1)
 HalfSpace([-1.0], -18.0)
 HalfSpace([-1.0], -18.1)
 HalfSpace([1.0], 18.0)

julia> removehredundancy!(Q)

julia> Q
Polyhedron CDDLib.CDDPolyhedron{1,Float64}:
2-element iterator of Polyhedra.HyperPlane{1,Float64,Array{Float64,1}}:
 HyperPlane([1.0], 18.1)
 HyperPlane([-1.0], -18.0)

2. The v0.7 result (polyhedron without constraints) is due to Dimension error for default Polyhedra backend in >= v0.7 #842. By hard-coding a 1 in the code, the result is as follows:

julia> cap = intersection(H, P)
HPolytope{Float64}(LazySets.HalfSpace{Float64}[HalfSpace{Float64}([1.0], 0.0), HalfSpace{Float64}([0.0], 1.0), HalfSpace{Float64}([-1.0], -0.0), HalfSpace{Float64}([-0.0], -1.0)])

This is still a strange result.

3. Assuming that the v0.7 result (polyhedron without constraints) can occur after fixing Dimension error for default Polyhedra backend in >= v0.7 #842: This polytope is not empty by definition. But it seems that Polyhedra interprets it like that. We should catch this case in the intersection methods and return an EmptySet in this case. Maybe we should check if Polyhedra returns a specific status flag or something to indicate empty intersections.

[mforets]

in #master and v0.7 one now gets:

julia> using LazySets, Polyhedra

julia> H = LazySets.HalfSpace([-1.0], -21.0); # x >= 21

julia> P = HPolytope([LazySets.HalfSpace([1.0], 18.1), # x <= 18.1
                      LazySets.HalfSpace([-1.0], -18.0)]); # x >= 18

julia> cap = intersection(H, P);

julia> isempty(cap)
true # correct

julia> cap.constraints
4-element Array{LazySets.HalfSpace{Float64},1}:
 LazySets.HalfSpace{Float64}([1.0], 0.0)  
 LazySets.HalfSpace{Float64}([0.0], 1.0)  
 LazySets.HalfSpace{Float64}([-1.0], -0.0)
 LazySets.HalfSpace{Float64}([-0.0], -1.0)

julia> remove_redundant_constraints(cap)
LP is not optimal

this is a case where "LP is not optimal" is indicating that it is empty.

[mforets]

After PR #1050

julia> using Revise, LazySets
[ Info: Recompiling stale cache file /Users/forets/.julia/compiled/v1.0/LazySets/NjrGc.ji for LazySets [b4f0291d-fe17-52bc-9479-3d1a343d9043]

julia> H = LazySets.HalfSpace([-1.0], -21.0); # x >= 21

julia> P = HPolytope([LazySets.HalfSpace([1.0], 18.1), # x <= 18.1

                             LazySets.HalfSpace([-1.0], -18.0)]); # x >= 18

julia> cap = intersection(H, P)
EmptySet{Float64}()

julia> using Polyhedra

julia> isempty(H∩P)
true

julia> isdisjoint(H, P)
true

julia> isdisjoint(P, H)
true

julia> using BenchmarkTools

julia> @btime intersection($H, $P)
  18.596 μs (55 allocations: 4.00 KiB)
EmptySet{Float64}()

julia> @btime isdisjoint($H, $P)
  17.676 μs (56 allocations: 3.73 KiB)
true