Convex hull algorithm from Polyhedra.jl produces invalid constraints

[schillic]

julia> X1 = HPolytope([HalfSpace([0.5, 1.5], 4.0), HalfSpace([0.0, -1.0], -1.0),
                       HalfSpace([0.5, -0.5], 0.0), HalfSpace([-1.0, 0.0], 0.0)])
julia> X2 = HPolytope([HalfSpace([-1.0, 0.0], 1.0), HalfSpace([0.5, 1.5], 4.0),
                       HalfSpace([0.0, -1.0], -1.0), HalfSpace([1.0, 0.0], 0.0)]);
julia> convex_hull(X1, X2)

ERROR: AssertionError: the half-space is inconsistent since it has a zero normal direction but the constraint is negative
Stacktrace:
 [1] macro expansion
   @ ~/.julia/packages/ReachabilityBase/fGisH/src/Assertions/Assertions.jl:21 [inlined]
 [2] convert(::Type{…}, P::Polyhedra.DefaultPolyhedron{…})
   @ LazySets ~/.julia/dev/LazySets/src/Sets/HPolytope.jl:169
 [3] #convex_hull#453
   @ ~/.julia/dev/LazySets/src/ConcreteOperations/convex_hull.jl:561 [inlined]
 [4] convex_hull(P1::HPolytope{Float64, Vector{Float64}}, P2::HPolytope{Float64, Vector{Float64}})
   @ LazySets ~/.julia/dev/LazySets/src/ConcreteOperations/convex_hull.jl:555
 [5] top-level scope
   @ REPL[10]:1
Some type information was truncated. Use `show(err)` to see complete types.

The half-space in question is ⟨a, x⟩ <= b for a = [-0.0, -0.0], b = -1.0.

[schillic]

I debugged this and the problem is rooted in a combination with Polyhedra and GLPK's presolver.

We have two polygons (with shared facet, which may be relevant), for which we can compute the convex hull. Note that I use the solver option GLPK.GLP_ON.

julia> using Polyhedra, JuMP, GLPK
julia> solver = JuMP.optimizer_with_attributes(GLPK.Optimizer, "presolve" => GLPK.GLP_ON);
julia> backend = DefaultLibrary{Float64}(solver);
julia> A, b = ([0.5 1.5; 0.0 -1.0; 0.5 -0.5; -1.0 0.0], [4.0, -1.0, 0.0, 0.0]);
julia> P = polyhedron(hrep(A, b), backend);
julia> A, b = ([-1.0 0.0; 0.5 1.5; 0.0 -1.0; 1.0 0.0], [1.0, 4.0, -1.0, 0.0]);
julia> Q = polyhedron(hrep(A, b), backend);
julia> R = convexhull(P, Q)
Polyhedron DefaultPolyhedron{Float64, MixedMatHRep{Float64, Matrix{Float64}}, MixedMatVRep{Float64, Matrix{Float64}}}:
8-element iterator of Vector{Float64}:
 [0.0, 1.0]
 [0.0, 2.666666666666667]
 [1.0, 1.0]
 [2.0, 2.0]
 [0.0, 1.0]
 [0.0, 2.666666666666667]
 [-1.0, 1.0]
 [-1.0, 3.0000000000000004]

First, if we remove redundant vertices, we can then also remove redundant constraints without problems - up to some warning from GLPK.

julia> removevredundancy!(R); R
Polyhedron DefaultPolyhedron{Float64, MixedMatHRep{Float64, Matrix{Float64}}, MixedMatVRep{Float64, Matrix{Float64}}}:
5-element iterator of Vector{Float64}:
 [0.0, 1.0]
 [-1.0, 1.0]
 [-1.0, 3.0000000000000004]
 [2.0, 2.0]
 [1.0, 1.0]

julia> removehredundancy!(R); R
glp_simplex: unable to recover undefined or non-optimal solution
Polyhedron DefaultPolyhedron{Float64, MixedMatHRep{Float64, Matrix{Float64}}, MixedMatVRep{Float64, Matrix{Float64}}}:
4-element iterator of Polyhedra.HalfSpace{Float64, Vector{Float64}}:
 HalfSpace([1.0, -1.0], 0.0)
 HalfSpace([-0.0, -0.6666666666666666], -0.6666666666666666)
 HalfSpace([0.25, 0.75], 2.0)
 HalfSpace([-0.2, -0.0], 0.2):
5-element iterator of Vector{Float64}:
 [0.0, 1.0]
 [-1.0, 1.0]
 [-1.0, 3.0000000000000004]
 [2.0, 2.0]
 [1.0, 1.0]

But if we remove redundant constraints before removing redundant vertices, we obtain an empty polyhedron. This is obviously wrong because the convex hull of points can never be empty. Also note that the set itself is not small, so numerical issues should not be to blame. But what may be relevant is that some vertices at the shared facets occur twice.

julia> R = convexhull(P, Q);  # redefine R
julia> removehredundancy!(R); R
Polyhedron DefaultPolyhedron{Float64, MixedMatHRep{Float64, Matrix{Float64}}, MixedMatVRep{Float64, Matrix{Float64}}}:
3-element iterator of HyperPlane{Float64, Vector{Float64}}:
 HyperPlane([1.0, 0.0], 0.0)
 HyperPlane([0.0, 1.0], 0.0)
 HyperPlane([0.0, 0.0], 1.0):
8-element iterator of Vector{Float64}:
 [0.0, 1.0]
 [0.0, 2.666666666666667]
 [1.0, 1.0]
 [2.0, 2.0]
 [0.0, 1.0]
 [0.0, 2.666666666666667]
 [-1.0, 1.0]
 [-1.0, 3.0000000000000004]

Finally, if I do not use the presolver in GLPK, then we also do not see this problem.

julia> solver = JuMP.optimizer_with_attributes(GLPK.Optimizer);
julia> backend = DefaultLibrary{Float64}(solver);
julia> A, b = ([0.5 1.5; 0.0 -1.0; 0.5 -0.5; -1.0 0.0], [4.0, -1.0, 0.0, 0.0]);
julia> P = polyhedron(hrep(A, b), backend);
julia> A, b = ([-1.0 0.0; 0.5 1.5; 0.0 -1.0; 1.0 0.0], [1.0, 4.0, -1.0, 0.0]);
julia> Q = polyhedron(hrep(A, b), backend);
julia> R = convexhull(P, Q);

julia> removehredundancy!(R); R
Polyhedron DefaultPolyhedron{Float64, MixedMatHRep{Float64, Matrix{Float64}}, MixedMatVRep{Float64, Matrix{Float64}}}:
4-element iterator of HalfSpace{Float64, Vector{Float64}}:
 HalfSpace([-1.0, -0.0], 1.0)
 HalfSpace([0.3333333333333335, 1.0], 2.666666666666667)
 HalfSpace([6.055761952500853e-17, -0.27272727272727265], -0.27272727272727265)
 HalfSpace([0.13043478260869565, -0.13043478260869568], -7.675020039799994e-17):
8-element iterator of Vector{Float64}:
 [0.0, 1.0]
 [0.0, 2.666666666666667]
 [1.0, 1.0]
 [2.0, 2.0]
 [0.0, 1.0]
 [0.0, 2.666666666666667]
 [-1.0, 1.0]
 [-1.0, 3.0000000000000004]

[schillic]

This turns out to be another instance of #3039, so I close this as a duplicate. I reported the example upstream.