Intersection of polytope and polyhedron results in unbounded polyhedron

[schillic]

julia> P = HPolytope([
 HalfSpace([0.0, 0.0, 0.0, 0.0, -1.0], -1.0),
 HalfSpace([-0.0, -0.0, -0.0, -0.0, 1.0], 1.0),
 HalfSpace([1.0, 0.0, 0.0, 0.0, 0.0], 0.308261),
 HalfSpace([-1.0, 0.0, 0.0, 0.0, 0.0], -0.2),
 HalfSpace([0.0, 1.0, 0.0, 0.0, 0.0], 0.1),
 HalfSpace([0.0, -1.0, 0.0, 0.0, 0.0], 0.106045),
 HalfSpace([0.0, 0.0, 1.0, 0.0, 0.0], 0.0153778),
 HalfSpace([0.0, 0.0, -1.0, 0.0, 0.0], 0.0),
 HalfSpace([0.0, 0.0, 0.0, 1.0, 0.0], 0.00075663),
 HalfSpace([0.0, 0.0, 0.0, -1.0, 0.0], 0.000156122)
]);

julia> Q = HPolyhedron([
 HalfSpace([-1.0, 0.0, 0.0, 0.0, 0.0], 0.0),
 HalfSpace([-0.714286, -1.0, 0.0, 0.0, 0.0], 0.0),
 HalfSpace([0.714286, 1.0, 0.0, 0.0, 0.0], 0.0)
]);

julia> isdisjoint(P, Q)
false

julia> isdisjoint(Q, P)
true

# the above is already worrying

# the intersection should be bounded, as P is already bounded
julia> R = intersection(P, Q)
HPolytope([
 HalfSpace([0.714286, 1.0, 0.0, 0.0, 0.0], 0.0),
 HalfSpace([-0.0, -0.0, -0.0, -0.0, 1.0], 1.0),
 HalfSpace([-0.714286, -1.0, -0.0, -0.0, -0.0], -0.0),
 HalfSpace([0.0, 0.0, 0.0, 0.0, -1.0], -1.0)
])

[schillic]

The isdisjoint calls the same method in both cases.

[schillic]

Here is a 2D projection with the same problems:

julia> p = HPolytope([HalfSpace([1.0, 0.0], 0.308261),
        HalfSpace([-1.0, 0.0], -0.2),
        HalfSpace([0.0, 1.0], 0.1),
        HalfSpace([0.0, -1.0], 0.106045)]);

julia> q = HPolyhedron([
        HalfSpace([-1.0, 0.0], 0.0),
        HalfSpace([-0.714286, -1.0], 0.0),
        HalfSpace([0.714286, 1.0], 0.0)]);

julia> isdisjoint(p, q)
false

julia> isdisjoint(q, p)
true

julia> intersection(p, q)
HPolytope([
 HalfSpace([0.714286, 1.0], 0.0),
 HalfSpace([-0.714286, -1.0], -0.0)
])

[schillic]

Note that Q/q is flat.

[schillic]

The unbounded intersection is an issue with the Polyhedra backend.
In v0.6, where we use CDDLib, everything is fine.

julia> R = intersection(P, Q)
LazySets.HPolytope{Float64}(LazySets.HalfSpace{Float64}[
LazySets.HalfSpace{Float64}([0.0, 0.0, 0.0, 0.0, -1.0], -1.0),
LazySets.HalfSpace{Float64}([1.0, 0.0, 0.0, 0.0, 0.0], 0.308261),
LazySets.HalfSpace{Float64}([-1.0, 0.0, 0.0, 0.0, 0.0], -0.2),
LazySets.HalfSpace{Float64}([0.0, 1.0, 0.0, 0.0, 0.0], 0.1),
LazySets.HalfSpace{Float64}([0.0, 0.0, 1.0, 0.0, 0.0], 0.0153778),
LazySets.HalfSpace{Float64}([0.0, 0.0, 0.0, 1.0, 0.0], 0.00075663),
LazySets.HalfSpace{Float64}([-0.0, -0.0, -0.0, -0.0, 1.0], 1.0),
LazySets.HalfSpace{Float64}([-1.0, -0.0, -0.0, -0.0, -0.0], -0.308261),
LazySets.HalfSpace{Float64}([1.0, -0.0, -0.0, -0.0, -0.0], 0.2),
LazySets.HalfSpace{Float64}([-0.0, -1.0, -0.0, -0.0, -0.0], -0.1),
LazySets.HalfSpace{Float64}([-0.0, -0.0, -1.0, -0.0, -0.0], -0.0153778),
LazySets.HalfSpace{Float64}([-0.0, -0.0, -0.0, -1.0, -0.0], -0.00075663)])

[mforets]

For the 2D example in this comment, i noticed that (below i'm working with CDDLib.Library()):

# p is a rectangle
julia> p = polyhedron(hrep([HalfSpace([1.0, 0.0], 0.308261),
                            HalfSpace([-1.0, 0.0], -0.2),
                            HalfSpace([0.0, 1.0], 0.1),
                            HalfSpace([0.0, -1.0], 0.106045)]), CDDLib.Library());

# q is an unbounded unidimensional set  (a "ray" from the origin towards `x->+oo`, `y->-oo` and slope `-0.714286`)
julia> q = polyhedron(hrep([HalfSpace([-1.0, 0.0], 0.0),
                            HalfSpace([-0.714286, -1.0], 0.0),
                            HalfSpace([0.714286, 1.0], 0.0)]), CDDLib.Library());

julia> P = Polyhedra.intersect(p, q)
Polyhedron CDDLib.Polyhedron{Float64}:
7-element iterator of HalfSpace{Float64,Array{Float64,1}}:
 HalfSpace([1.0, 0.0], 0.308261)
 HalfSpace([-1.0, 0.0], -0.2)
 HalfSpace([0.0, 1.0], 0.1)
 HalfSpace([0.0, -1.0], 0.106045)
 HalfSpace([-1.0, 0.0], 0.0)
 HalfSpace([-0.714286, -1.0], 0.0)
 HalfSpace([0.714286, 1.0], 0.0)

julia> Polyhedra.removehredundancy!(P)

# the polyhedron P is unbounded in the y coordinate: 0.2 <= x <= 0.308261, y <= 0.1
# but it shouldn't: y should be bounded as well since `x` is bounded 
julia> P
Polyhedron CDDLib.Polyhedron{Float64}:
3-element iterator of HyperPlane{Float64,Array{Float64,1}}:
 HyperPlane([1.0, 0.0], 0.308261)
 HyperPlane([-1.0, 0.0], -0.2)
 HyperPlane([0.0, 1.0], 0.1)

Also:

julia> LazySets.intersection(convert(HPolytope, p), convert(HPolyhedron, q), prunefunc=remove_redundant_constraints!)
ERROR: LP is not optimal
Stacktrace:
 [1] error(::String) at ./error.jl:33
 [2] #remove_redundant_constraints!#51(::GLPKMathProgInterface.GLPKInterfaceLP.GLPKSolverLP, ::Function, ::HPolytope{Float64}) at /Users/forets/.julia/dev/LazySets/src/HPolyhedron.jl:488
 [3] remove_redundant_constraints! at /Users/forets/.julia/dev/LazySets/src/HPolyhedron.jl:475 [inlined]
 [4] #intersection#159 at /Users/forets/.julia/dev/LazySets/src/concrete_intersection.jl:386 [inlined]
 [5] (::getfield(LazySets, Symbol("#kw##intersection")))(::NamedTuple{(:prunefunc,),Tuple{typeof(remove_redundant_constraints!)}}, ::typeof(intersection), ::HPolytope{Float64}, ::HPolyhedron{Float64}) at ./none:0
 [6] top-level scope at none:0

[mforets]

Now on #master, the example in the initial report gives isdisjoint(P, Q) = isdisjoint(Q, P) = true, and the concrete intersection returns the empty set as expected.

The same results (isdisjoint true and intersection returns empty set) for the 2D example.

Since the intersection is empty (see this figure) i think that we can consider this issue closed.