search

JuliaPolyhedra / Polyhedra.jl

Polyhedral Computation Interface
Other
172 stars 27 forks source link

AssertionError in intersect #251

Open schillic opened 3 years ago

schillic commented 3 years ago

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 commented 1 year ago

@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 commented 1 year ago

Nope, sorry.

blegat commented 1 year ago

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 commented 1 year ago

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:

https://github.com/JuliaPolyhedra/Polyhedra.jl/blob/e56953d40dc0a065da634777a2ded853880d0a2b/src/doubledescription.jl#L358-L361

which is called from here:

https://github.com/JuliaPolyhedra/Polyhedra.jl/blob/e56953d40dc0a065da634777a2ded853880d0a2b/src/doubledescription.jl#L372

(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).