Conversion of rational vrep polytope to hrep seems buggy

[aquohn]

Julia throws an error when I try to get the constraints of a vrep polytope (using the LazySets.jl interface); looks like a bug in the Polyhedra.jl implementation?

julia> rpoly
VPolytope{Rational, Vector{Rational}}(Vector{Rational}[[1//1, 1//1, 1//1, 1//1, 1//1, 1//1, 1//1, 1//1, 1//1, 1//1, 1//1], [0//1, 1//1, 1//1, 1//1, 1//1, 0//1, 1//1, 0//1, 1//1, 0//1, 1//1], [1//1, 0//1, 1//1, 1//1, 1//1, 1//1, 0//1, 1//1, 0//1, 1//1, 0//1], [0//1, 0//1, 1//1, 1//1, 1//1, 0//1, 0//1, 0//1, 0//1, 0//1, 0//1], [1//1, 1//1, 0//1, 1//1, 1//1, 0//1, 0//1, 1//1, 1//1, 1//1, 1//1], [0//1, 1//1, 0//1, 1//1, 1//1, 0//1, 0//1, 0//1, 1//1, 0//1, 1//1], [1//1, 0//1, 0//1, 1//1, 1//1, 0//1, 0//1, 1//1, 0//1, 1//1, 0//1], [0//1, 0//1, 0//1, 1//1, 1//1, 0//1, 0//1, 0//1, 0//1, 0//1, 0//1], [1//1, 1//1, 1//1, 0//1, 1//1, 1//1, 1//1, 0//1, 0//1, 1//1, 1//1], [0//1, 1//1, 1//1, 0//1, 1//1, 0//1, 1//1, 0//1, 0//1, 0//1, 1//1]  …  [1//1, 0//1, 0//1, 1//1, 0//1, 0//1, 0//1, 1//1, 0//1, 0//1, 0//1], [0//1, 0//1, 0//1, 1//1, 0//1, 0//1, 0//1, 0//1, 0//1, 0//1, 0//1], [1//1, 1//1, 1//1, 0//1, 0//1, 1//1, 1//1, 0//1, 0//1, 0//1, 0//1], [0//1, 1//1, 1//1, 0//1, 0//1, 0//1, 1//1, 0//1, 0//1, 0//1, 0//1], [1//1, 0//1, 1//1, 0//1, 0//1, 1//1, 0//1, 0//1, 0//1, 0//1, 0//1], [0//1, 0//1, 1//1, 0//1, 0//1, 0//1, 0//1, 0//1, 0//1, 0//1, 0//1], [1//1, 1//1, 0//1, 0//1, 0//1, 0//1, 0//1, 0//1, 0//1, 0//1, 0//1], [0//1, 1//1, 0//1, 0//1, 0//1, 0//1, 0//1, 0//1, 0//1, 0//1, 0//1], [1//1, 0//1, 0//1, 0//1, 0//1, 0//1, 0//1, 0//1, 0//1, 0//1, 0//1], [0//1, 0//1, 0//1, 0//1, 0//1, 0//1, 0//1, 0//1, 0//1, 0//1, 0//1]])

julia> rpoly |> constraints_list
ERROR: MethodError: no method matching Ray{Rational, Vector{Rational}}(::Vector{Rational{Int64}})
Closest candidates are:
  Ray{T, AT}(::AT) where {T, AT<:AbstractVector{T}} at /home/aquohn/.julia/packages/Polyhedra/NGGvZ/src/elements.jl:118
  Ray{T, AT}(::Ray) where {T, AT} at /home/aquohn/.julia/packages/Polyhedra/NGGvZ/src/elements.jl:122
Stacktrace:
  [1] -(elem::Ray{Rational, Vector{Rational}})
    @ Polyhedra ~/.julia/packages/Polyhedra/NGGvZ/src/elements.jl:153
  [2] _λ_proj(r::Ray{Rational, Vector{Rational}}, line::Polyhedra.Line{Rational, Vector{Rational}}, h::Polyhedra.HalfSpace{Rational, Vector{Rational}})
    @ Polyhedra ~/.julia/packages/Polyhedra/NGGvZ/src/doubledescription.jl:422
  [3] line_project(el::Ray{Rational, Vector{Rational}}, line::Polyhedra.Line{Rational, Vector{Rational}}, h::Polyhedra.HalfSpace{Rational, Vector{Rational}})
    @ Polyhedra ~/.julia/packages/Polyhedra/NGGvZ/src/doubledescription.jl:431
  [4] add_element!(data::Polyhedra.DoubleDescriptionData{Vector{Rational}, Ray{Rational, Vector{Rational}}, Polyhedra.Line{Rational, Vector{Rational}}, Polyhedra.HalfSpace{Rational, Vector{Rational}}}, k::Int64, el::Ray{Rational, Vector{Rational}}, tight::BitSet)
    @ Polyhedra ~/.julia/packages/Polyhedra/NGGvZ/src/doubledescription.jl:311
  [5] doubledescription(hr::MixedMatHRep{Rational, Matrix{Rational}}, #unused#::MathOptInterface.OptimizerWithAttributes)
    @ Polyhedra ~/.julia/packages/Polyhedra/NGGvZ/src/doubledescription.jl:477
  [6] doubledescription(v::Polyhedra.Hull{Rational, Vector{Rational}, Int64}, solver::MathOptInterface.OptimizerWithAttributes; kws::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Polyhedra ~/.julia/packages/Polyhedra/NGGvZ/src/doubledescription.jl:538
  [7] doubledescription(v::Polyhedra.Hull{Rational, Vector{Rational}, Int64}, solver::MathOptInterface.OptimizerWithAttributes)
    @ Polyhedra ~/.julia/packages/Polyhedra/NGGvZ/src/doubledescription.jl:535
  [8] computehrep!(p::DefaultPolyhedron{Rational, Polyhedra.Intersection{Rational, Vector{Rational}, Int64}, Polyhedra.Hull{Rational, Vector{Rational}, Int64}})
    @ Polyhedra ~/.julia/packages/Polyhedra/NGGvZ/src/defaultlibrary.jl:91
  [9] hrep
    @ ~/.julia/packages/Polyhedra/NGGvZ/src/defaultlibrary.jl:96 [inlined]
 [10] hyperplanetype(p::DefaultPolyhedron{Rational, Polyhedra.Intersection{Rational, Vector{Rational}, Int64}, Polyhedra.Hull{Rational, Vector{Rational}, Int64}})
    @ Polyhedra ~/.julia/packages/Polyhedra/NGGvZ/src/iterators.jl:175
 [11] _broadcast_getindex_evalf
    @ ./broadcast.jl:648 [inlined]
 [12] _broadcast_getindex
    @ ./broadcast.jl:621 [inlined]
 [13] #19
    @ ./broadcast.jl:1098 [inlined]
 [14] ntuple
    @ ./ntuple.jl:48 [inlined]
 [15] copy
    @ ./broadcast.jl:1098 [inlined]
 [16] materialize
    @ ./broadcast.jl:883 [inlined]
 [17] hyperplanes
    @ ~/.julia/packages/Polyhedra/NGGvZ/src/iterators.jl:183 [inlined]
 [18] allhalfspaces
    @ ~/.julia/packages/Polyhedra/NGGvZ/src/iterators.jl:299 [inlined]
 [19] convert(#unused#::Type{HPolytope}, P::DefaultPolyhedron{Rational, Polyhedra.Intersection{Rational, Vector{Rational}, Int64}, Polyhedra.Hull{Rational, Vector{Rational}, Int64}})
    @ LazySets ~/.julia/packages/LazySets/qnBZn/src/Sets/HPolytope.jl:160
 [20] #tohrep#215
    @ ~/.julia/packages/LazySets/qnBZn/src/Sets/VPolytope.jl:490 [inlined]
 [21] tohrep(P::VPolytope{Rational, Vector{Rational}})
    @ LazySets ~/.julia/packages/LazySets/qnBZn/src/Sets/VPolytope.jl:485
 [22] constraints_list(P::VPolytope{Rational, Vector{Rational}})
    @ LazySets ~/.julia/packages/LazySets/qnBZn/src/Sets/VPolytope.jl:406
 [23] |>(x::VPolytope{Rational, Vector{Rational}}, f::typeof(constraints_list))
    @ Base ./operators.jl:858
 [24] top-level scope
    @ REPL[31]:1

[schillic]

The problem is that rpoly has insufficient type information: VPolytope{Rational, Vector{Rational}} does not fully reveal the numeric type (the {Int} in Rational{Int} is missing). This could be either a LazySets bug or a problem in how you create the VPolytope (not possible to say from the one line).

[aquohn]

Thanks for the prompt response! It was indeed the missing numeric type: I supplied the points for the vrep as Vector{Rational}; changing to Vector{Rational{Int64}} fixed it. The Int64 seems to be introduced at [19], which is given by

function convert(::Type{HPolytope}, P::HRep{N}) where {N}
  VT = Polyhedra.hvectortype(P)
  constraints = Vector{LinearConstraint{N, VT}}()
  for hi in Polyhedra.allhalfspaces(P)
    a, b = hi.a, hi.β
    if isapproxzero(norm(a))
      @assert b >= zero(N) "the half-space is inconsistent since it has a " *
      "zero normal direction but the constraint is negative"
      continue
    end
    push!(constraints, HalfSpace(hi.a, hi.β))
  end
  return HPolytope(constraints)
end

which in turn is called by

function tohrep(P::VPolytope{N};
    backend=default_polyhedra_backend(P)) where {N}
  vl = P.vertices
  if isempty(vl)
    return EmptySet{N}(dim(P))
  end
  require(:Polyhedra; fun_name="tohrep")
  return convert(HPolytope, polyhedron(P; backend=backend))
end

It seems the assumption regarding the numeric type is being made by the call to polyhedron(P) from Polyhedra.jl; is it supposed to assume the numeric type is Int64?

[schillic]

Sorry, I do not follow. Does it work or not? If you create the Vector with the full type, LazySets is supposed to preserve it through all calls and eventually to Polyhedra. If this is not the case, I suggest to open an issue with LazySets.

About guessing Int64: I am not sure this should be the job of the library.

[aquohn]

Yes, it works with the full type. I'm trying to trace where the type assumption Int64 was introduced, since that breaks the conversion. It seems to have been introduced at [19] in the stack trace.