AssertionError in intersect

[schillic]

The following example is a reduced version (I removed dimensions and vertices until the error vanishes) of a more complex input.

using Polyhedra, JuMP, GLPK

vP = [[-0.20301280791607568, 0.007131220483095747, -0.02988255782931391, -0.17050077196699046, -0.10459728765273806, -0.3940570183286591, -0.1060778630887406, -0.2083884887166534]];

vQ = [[0.26944343620620703, 0.5820831768940575, 0.043944752271130684, 0.0, 0.0, 0.0, 0.0, 0.002290710865434633], [0.269764654943868, 0.5826024342879536, 0.04424939599453754, -0.3243733749380128, -0.27971766698656575, -0.040712222052006995, -0.4230541656638591, 0.0020389818225279343], [0.269764654943868, 0.5826024342879536, 0.04424939599453754, 0.0, 0.0, 0.0, 0.0, 0.0020389818225279343], [0.3608261278169685, 0.6795241675058606, 0.13674690705026338, -0.4026336283204254, -0.3521265118897102, -0.0985458649168468, -0.4867600723992247, 0.022141261211450644], [0.3608261278169685, 0.6795241675058606, 0.13674690705026338, 0.0, 0.0, 0.0, 0.0, 0.022141261211450644], [0.37879544708728713, 0.7586303128811639, 0.11059986791879892, -0.4390178801299712, -0.37783324554721326, -0.08204843892446032, -0.5496916242909842, 0.017745282399702605], [0.37879544708728713, 0.7586303128811639, 0.11059986791879892, 0.0, 0.0, 0.0, 0.0, 0.017745282399702605], [0.35699186847910663, 0.7349374913159256, 0.08851175424220448, -0.42014699571976977, -0.36009900240651904, -0.06804164876801723, -0.534431302171298, 0.013816230720892301], [0.35699186847910663, 0.7349374913159256, 0.08851175424220448, 0.0, 0.0, 0.0, 0.0, 0.013816230720892301], [0.3566093175846898, 0.734546246932556, 0.08812122743510388, -0.41982256106796995, -0.3598078125821968, -0.06780390511187857, -0.5341638922672701, 0.013702811188085457]];

backend = Polyhedra.DefaultLibrary{Float64}(JuMP.optimizer_with_attributes(GLPK.Optimizer));
P = polyhedron(Polyhedra.vrep(vP), backend);
Q = polyhedron(Polyhedra.vrep(vQ), backend);
Polyhedra.intersect(P, Q)

ERROR: LoadError: AssertionError: 0 <= value2 / (value2 - value1) <= 1
Stacktrace:
  [1] combine(β::Float64, r1::Ray{Float64, Vector{Float64}}, value1::Float64, r2::Ray{Float64, Vector{Float64}}, value2::Float64)
    @ Polyhedra ~/.julia/packages/Polyhedra/lR2w5/src/doubledescription.jl:361
  [2] combine
    @ ~/.julia/packages/Polyhedra/lR2w5/src/doubledescription.jl:372 [inlined]
  [3] add_intersection!(data::Polyhedra.DoubleDescriptionData{Vector{Float64}, Ray{Float64, Vector{Float64}}, Line{Float64, Vector{Float64}}, HalfSpace{Float64, Vector{Float64}}}, idx1::Polyhedra.CutoffRayIndex, idx2::Polyhedra.CutoffRayIndex, hp_idx::Int64, hs_idx::Int64)
    @ Polyhedra ~/.julia/packages/Polyhedra/lR2w5/src/doubledescription.jl:404
  [4] add_intersection!
    @ ~/.julia/packages/Polyhedra/lR2w5/src/doubledescription.jl:391 [inlined]
  [5] doubledescription(hr::MixedMatHRep{Float64, Matrix{Float64}}, #unused#::MathOptInterface.OptimizerWithAttributes)
    @ Polyhedra ~/.julia/packages/Polyhedra/lR2w5/src/doubledescription.jl:514
  [6] doubledescription(v::Polyhedra.Hull{Float64, Vector{Float64}, Int64}, solver::MathOptInterface.OptimizerWithAttributes; kws::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Polyhedra ~/.julia/packages/Polyhedra/lR2w5/src/doubledescription.jl:538
  [7] doubledescription(v::Polyhedra.Hull{Float64, Vector{Float64}, Int64}, solver::MathOptInterface.OptimizerWithAttributes)
    @ Polyhedra ~/.julia/packages/Polyhedra/lR2w5/src/doubledescription.jl:535
  [8] computehrep!(p::DefaultPolyhedron{Float64, Polyhedra.Intersection{Float64, Vector{Float64}, Int64}, Polyhedra.Hull{Float64, Vector{Float64}, Int64}})
    @ Polyhedra ~/.julia/packages/Polyhedra/lR2w5/src/defaultlibrary.jl:91
  [9] hrep
    @ ~/.julia/packages/Polyhedra/lR2w5/src/defaultlibrary.jl:96 [inlined]
 [10] hyperplanetype(p::DefaultPolyhedron{Float64, Polyhedra.Intersection{Float64, Vector{Float64}, Int64}, Polyhedra.Hull{Float64, Vector{Float64}, Int64}})
    @ Polyhedra ~/.julia/packages/Polyhedra/lR2w5/src/iterators.jl:175
 [11] _broadcast_getindex_evalf
    @ ./broadcast.jl:648 [inlined]
 [12] _broadcast_getindex
    @ ./broadcast.jl:621 [inlined]
 [13] _getindex
    @ ./broadcast.jl:644 [inlined]
 [14] _broadcast_getindex
    @ ./broadcast.jl:620 [inlined]
 [15] #19
    @ ./broadcast.jl:1098 [inlined]
 [16] ntuple
    @ ./ntuple.jl:49 [inlined]
 [17] copy
    @ ./broadcast.jl:1098 [inlined]
 [18] materialize
    @ ./broadcast.jl:883 [inlined]
 [19] maphyperplanes
    @ ~/.julia/packages/Polyhedra/lR2w5/src/iterators.jl:188 [inlined]
 [20] hmap
    @ ~/.julia/packages/Polyhedra/lR2w5/src/iterators.jl:357 [inlined]
 [21] intersect(::DefaultPolyhedron{Float64, Polyhedra.Intersection{Float64, Vector{Float64}, Int64}, Polyhedra.Hull{Float64, Vector{Float64}, Int64}}, ::DefaultPolyhedron{Float64, Polyhedra.Intersection{Float64, Vector{Float64}, Int64}, Polyhedra.Hull{Float64, Vector{Float64}, Int64}})
    @ Polyhedra ~/.julia/packages/Polyhedra/lR2w5/src/repop.jl:21

Tested with Julia v1.6.1 and

  [67491407] Polyhedra v0.6.14
  [4076af6c] JuMP v0.21.8
  [60bf3e95] GLPK v0.13.0

[forrestlaine]

@schillic did you ever identify what the issue was here? I am also running into this assertion error when running the double description algorithm.

[schillic]

Nope, sorry.

[blegat]

IIRC, I tried debugging it and it didn't seem to be a bug, the problem is just very ill-conditioned but I might be confusing with another issue

[schillic]

I reduced the input a bit (unfortunately still 8-dimensional).

vP = [zeros(8)];

vQ = [ zeros(8),
[0.2697646549, 0.582602434, 0.044249395, 0.3243733749, 0.279717666, 0.040712222, 0.423054165, 0.00203898],
[0.2697646549, 0.582602434, 0.044249395, 0, 0, 0, 0, 0.00203898],
[0.3608261278, 0.679524167, 0.136746907, 0.402633628, 0.35212651, 0.098545864, 0.486760072, 0.02214126],
[0.3608261278, 0.679524167, 0.136746907, 0, 0, 0, 0, 0.02214126],
[0.378795447, 0.7586303128, 0.110599867, 0.43901788, 0.377833245, 0.0820484389, 0.549691624, 0.017745282],
[0.378795447, 0.7586303128, 0.110599867, 0, 0, 0, 0, 0.017745282],
[0.35699186847, 0.7349374913, 0.08851175424, 0.42014699571, 0.3600990024, 0.06804164876, 0.53443130217, 0.0138162307],
[0.35699186847, 0.7349374913, 0.08851175424, 0, 0, 0, 0, 0.0138162307],
[0.35660931758, 0.7345462469, 0.08812122743, 0.419822561067, 0.35980781258, 0.06780390511, 0.534163892267, 0.0137028111] ];

Interestingly, if I change the order of the vertices, the problem does not occur. (For instance, move the zeros(8) further down in vQ.)

The error is in this method:

Polyhedra.jl/src/doubledescription.jl

Lines 358 to 361 in e56953d

	function combine(β, r1::Polyhedra.Ray, value1, r2::Polyhedra.Ray, value2)
	# should take
	# λ = value2 / (value2 - value1)
	@assert 0 <= value2 / (value2 - value1) <= 1

which is called from here:

Polyhedra.jl/src/doubledescription.jl

Line 372 in e56953d

combine(h, el1, el2) = combine(h.β, el1, h.a ⋅ el1, el2, h.a ⋅ el2)

(Side note: The parameter β is not used in that method.)

Looking at typical values of value1 and value2, value1 is always negative and value2 is always positive - except when this assertion fails, then value1 is small but positive. For the above example I get:

h.a = [-1.0, -0.2697646549, -0.582602434, -0.044249395, -0.3243733749, -0.279717666, -0.040712222, -0.423054165, -0.00203898]
el1 = r1 = Ray([0.0, 29.090432340937888, -12.242070824013831, -16.36861456149422, 78.48690222274764, 20.316388217194387, -57.27032575277269, -68.10076386587039, 4.409582784205131])
(value1, value2) = (1.4249591747880763e-8, 0.00424417718061143)

So the problem is that value1 = h.a · el1 = 1.4249591747880763e-8. h.a is the negative second vertex of vQ plus a leading -1.0. I do not know where el1 comes from. Both these vectors look pretty normal to me, so numerical errors should not be the cause (except maybe further down the stack when el1 is computed - this I do not know).