Removing constraints of empty polyhedron results in non-empty set

[mforets]

The following example shows what i think is a bug in removehredundancy!:

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

(reported here @schillic)

---

(EDITED)

I had a quick look at a workaround using a Polyhedra.Ray:

julia> v = polyhedron(vrep([Polyhedra.Ray([-0.714286, 1.0])]), CDDLib.Library())
Polyhedron CDDLib.Polyhedron{Float64}:
1-element iterator of Array{Float64,1}:
 [0.0, 0.0],
1-element iterator of Ray{Float64,Array{Float64,1}}:
 Ray([-0.714286, 1.0])

julia> P = Polyhedra.intersect(p, v)
Polyhedron CDDLib.Polyhedron{Float64}:
1-element iterator of HyperPlane{Float64,Array{Float64,1}}:
 HyperPlane([-1.4, -1.0], 0.0),
5-element iterator of Polyhedra.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)

julia> Polyhedra.removehredundancy!(P)

julia> P 
Polyhedron CDDLib.Polyhedron{Float64}:
3-element iterator of HyperPlane{Float64,Array{Float64,1}}:
 HyperPlane([-1.4, -1.0], 0.0)
 HyperPlane([1.0, 0.0], 0.308261)
 HyperPlane([-1.0, 0.0], -0.2)

i just plotted this in geogebra; the intersection is actually empty :)

indeed:

julia> Polyhedra.isempty(Polyhedra.intersect(p, q))
true

and

julia> Polyhedra.isempty(Polyhedra.intersect(p, v))
true

---

given this, i take back my assumption that removehredundancy! is buggy.

[schillic]

Is it not even worse if the removal of redundant constraints from an empty set results in a non-empty set?

[mforets]

@schillic true, i've reopened this issue with a different title to suggest that some action should be taken, either make it clear in removehredundancy! docs that the input polyhedron should not be empty, or revise the algorithm.

[blegat]

Isn't the bug in CDD ?

[schillic]

Isn't the bug in CDD ?

The problem also occurs with the "default library". (EDIT: simplified the numbers)

for backend in [CDDLib.Library(), default_library(2, Float64)]
    p = polyhedron(hrep([HalfSpace([1.0, 0.0], 0.3),
                         HalfSpace([-1.0, 0.0], -0.2),
                         HalfSpace([0.0, 1.0], 0.1),
                         HalfSpace([0.0, -1.0], 0.1)]), backend)
    q = polyhedron(hrep([HalfSpace([-1.0, 0.0], 0.0),
                         HalfSpace([-0.7, -1.0], 0.0),
                         HalfSpace([0.7, 1.0], 0.0)]), backend)
    P = Polyhedra.intersect(p, q)
    Polyhedra.removehredundancy!(P)
    println(P)
end

HyperPlane([1.0, 0.0], 0.3)
 ∩ HyperPlane([-1.0, 0.0], -0.2)
 ∩ HyperPlane([0.0, 1.0], 0.1)

HyperPlane([0.7, 1.0], 0.0) :