JuliaPolyhedra / Polyhedra.jl

Polyhedral Computation Interface
Other
172 stars 27 forks source link

removevredundancy! applied to H-rep removes all vertices #270

Closed schillic closed 2 years ago

schillic commented 3 years ago

Below is an example where removevredundancy! removes all vertices. The (eight-dimensional) polytope is very small and the bug is likely related to floating-point issues.

julia> using Polyhedra
julia> import JuMP, GLPK

julia> A = [-1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; -0.00106131548486069 0.016634209829441304 -0.0006666769554416961 -1.0 -0.0 -0.0 -0.0 -0.0; -0.0004116155738065077 0.01693983068952557 -2.6779479610930705e-5 -0.0 -1.0 -0.0 -0.0 -0.0; -0.0024283354284051484 0.08467041915859919 -0.0005000035065344641 -0.0 -0.0 -1.0 -0.0 -0.0; -0.002024394302735674 0.08367543404800237 -0.00013321128253368072 -0.0 -0.0 -0.0 -1.0 -0.0; -0.0009309923368121566 0.009059293773203692 -0.0006675849264721348 -0.0 -0.0 -0.0 -0.0 -1.0; 0.00106131548486069 -0.016634209829441304 0.0006666769554416961 1.0 0.0 0.0 0.0 0.0; 0.0004116155738065077 -0.01693983068952557 2.6779479610930705e-5 0.0 1.0 0.0 0.0 0.0; 0.0024283354284051484 -0.08467041915859919 0.0005000035065344641 0.0 0.0 1.0 0.0 0.0; 0.002024394302735674 -0.08367543404800237 0.00013321128253368072 0.0 0.0 0.0 1.0 0.0; 0.0009309923368121566 -0.009059293773203692 0.0006675849264721348 0.0 0.0 0.0 0.0 1.0; -0.08320472072273968 2.548630897981053 0.001134866706464998 0.0 0.0 0.0 0.0 0.0; 1210.55555992692 -37126.858723322744 -13.95377106842899 -0.0 -0.0 -0.0 -0.0 -0.0; -1860.676468633397 57155.38024488611 18.21165578673308 -0.0 -0.0 -0.0 -0.0 -0.0];

julia> b = [-20.995, -1.007487377909021, -0.9877340129555137, -0.9705208833525414, -0.9548402631790559, -0.025711479466835108, 1.007487377909021, 0.9877340129555137, 0.9705208833525414, 0.9548402631790559, 0.025711479466835108, 0.8046759919340509, -11751.33121007926, 18148.933864538412];

julia> solver = JuMP.optimizer_with_attributes(GLPK.Optimizer, "presolve" => GLPK.ON);

julia> lib = DefaultLibrary{Float64}(solver);

julia> P = polyhedron(hrep(A, b), lib);

julia> removevredundancy!(P)
NO_REDUNDANCY::Redundancy = 2

julia> collect(points(P))
Vector{Float64}[]

If you first call points, then removevredundancy! does not remove the vertices anymore.

julia> P = polyhedron(hrep(A, b), lib);

julia> collect(points(P))
4-element Vector{Vector{Float64}}:
 [20.998401873559185, 1.0008152770580858, 0.998722721562126, 1.001183394739038, 0.9960176397894045, 1.0037698042890895, 0.9959418296347409, 0.014562075610188819]
 [20.995, 1.0007055006170342, 0.9956767683957485, 1.0011872098225525, 0.9960172620282571, 1.0037702933491626, 0.9959399365321673, 0.014566281663794974]
 [20.995, 1.0007053185484005, 0.9962481717932737, 1.0011868258525074, 0.9960172436421597, 1.0037699922296326, 0.9959398451801159, 0.01456589855408661]
 [20.995, 1.0007051089510135, 0.9967188758624705, 1.0011865085584646, 0.9960172274864054, 1.0037697391292488, 0.9959397649388708, 0.014565582420340884]

julia> removevredundancy!(P)

julia> collect(points(P))
4-element Vector{Vector{Float64}}:
 [20.998401873559185, 1.0008152770580858, 0.998722721562126, 1.001183394739038, 0.9960176397894045, 1.0037698042890895, 0.9959418296347409, 0.014562075610188819]
 [20.995, 1.0007055006170342, 0.9956767683957485, 1.0011872098225525, 0.9960172620282571, 1.0037702933491626, 0.9959399365321673, 0.014566281663794974]
 [20.995, 1.0007053185484005, 0.9962481717932737, 1.0011868258525074, 0.9960172436421597, 1.0037699922296326, 0.9959398451801159, 0.01456589855408661]
 [20.995, 1.0007051089510135, 0.9967188758624705, 1.0011865085584646, 0.9960172274864054, 1.0037697391292488, 0.9959397649388708, 0.014565582420340884]

So maybe the same operation used for converting to V-rep should be used in removevredundancy! applied to an H-rep?

Note that CDDLib also finds four vertices:


julia> import CDDLib

julia> lib = CDDLib.Library();

julia> P = polyhedron(hrep(A, b), lib);

julia> removevredundancy!(P)

julia> collect(points(P))
4-element Vector{Vector{Float64}}:
 [20.994999999999997, 1.0007051089510148, 0.9967188758598511, 1.0011865085584664, 0.9960172274864056, 1.0037697391292506, 0.9959397649388715, 0.014565582420342806]
 [20.994999999999997, 1.0007053185484, 0.996248171794203, 1.0011868258525067, 0.9960172436421596, 1.0037699922296324, 0.9959398451801159, 0.014565898554086151]
 [20.99840187353354, 1.0008152770572565, 0.9987227215443998, 1.0011833947390634, 0.9960176397894016, 1.0037698042890908, 0.995941829634726, 0.01456207561021718]
 [20.994999999999997, 1.0007055006170327, 0.9956767684009857, 1.0011872098225492, 0.9960172620282569, 1.0037702933491603, 0.9959399365321667, 0.014566281663791626]