Open schillic opened 1 year ago
Sorry for the delay on this. The output seems weird, as the constraint is LessThan
:
https://github.com/JuliaPolyhedra/Polyhedra.jl/blob/07567d48fa60e6acc0d0ed1af012eb82c82ff78c/src/linearity.jl#L170
The ConstraintPrimal
should be nonpositive for a feasible solution so if it is not zero, it should be negative. A MWE would definitely help.
Below is a simplified version. It definitely has to do with presolve in GLPK, so this might be the wrong place to report. If you look at this plot,
you see that the set is a normal polygon with no numerical issues. The only thing to note is that there are two almost-orthogonal diagonal constraints, and the difference in orthogonality is probably below numerical precision. But I would say that this should not matter because they are not parallel but bound in different directions.
Some explanation: This comes from an elimination algorithm that first extends the set with a dummy dimension and then projects some dimensions away.
julia> using Polyhedra, CDDLib, MathOptInterface, GLPK
julia> Ab = [
([1.0, 0.0], 10.0)
([0.0, 1.0], 1.5)
([0.0, -1.0], 1.5)
([1.0, 1.5], 11.125)
([-1.0, -1.5000000000000002], -6.875)
([-4.0, -1.0], -24.5)
];
julia> M = [1. 0];
julia> m, n = size(M);
julia> ₋Id_m = hcat(-1.0);
julia> Ax_leq_b = [Polyhedra.HalfSpace(vcat(zeros(m), Abi[1]), Abi[2]) for Abi in Ab];
julia> y_eq_Mx = [Polyhedra.HyperPlane(vcat(₋Id_m[i, :], Vector(M[i, :])), 0.0) for i in 1:m];
julia> Phrep = Polyhedra.hrep(y_eq_Mx, Ax_leq_b);
julia> backend = CDDLib.Library(:float);
julia> Phrep = polyhedron(Phrep, backend);
julia> method = Polyhedra.BlockElimination();
julia> Peli_block = Polyhedra.eliminate(Phrep, (m+1):(m+n), method);
### With presolve, the result is not correct
julia> solver = MathOptInterface.OptimizerWithAttributes(GLPK.Optimizer, Pair{MathOptInterface.AbstractOptimizerAttribute, Any}[MathOptInterface.RawOptimizerAttribute("presolve") => 1]);
julia> Polyhedra.removeduplicates(Polyhedra.hrep(Peli_block), solver)
H-representation CDDInequalityMatrix{Float64, Float64}:
2-element iterator of HyperPlane{Float64, Vector{Float64}}:
HyperPlane([1.0], 0.0)
HyperPlane([0.0], 1.0)
### Without presolve, the result is correct (just contains redundant constraints)
julia> solver2 = MathOptInterface.OptimizerWithAttributes(GLPK.Optimizer, Pair{MathOptInterface.AbstractOptimizerAttribute, Any}[]);
julia> Polyhedra.removeduplicates(Polyhedra.hrep(Peli_block), solver2)
H-representation CDDInequalityMatrix{Float64, Float64}:
7-element iterator of HalfSpace{Float64, Vector{Float64}}:
HalfSpace([-4.0], -23.0)
HalfSpace([-1.0], -4.625)
HalfSpace([-0.0], 3.0)
HalfSpace([1.0], 13.375)
HalfSpace([1.480297366166875e-16], 4.250000000000001)
HalfSpace([-5.0], -25.625)
HalfSpace([1.0], 10.0)
### The following part leads to the funny error message
julia> rep = Polyhedra.hrep(Peli_block);
julia> aff = Polyhedra._linearity_space(rep, true);
julia> Polyhedra._hasnonlinearity(rep)
true
julia> isnothing(Polyhedra.detect_new_linearities(rep, solver))
true
julia> Polyhedra.detect_new_linearities(rep, solver; verbose=1)
[ Info: The polyhedron is empty as 3.0 is negative.
Thanks for the detailed MWE, that should be easier to investigate
Here is another example. I first compute the convex hull of the two polygons below and then call detecthlinearity
. Again this is problem does not occur without the presolve
option.
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]
julia> detecthlinearity(hrep(R), solver)
H-representation MixedMatHRep{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)
I was debugging a bigger example (here) that leads to a (probably) wrong result when using
GLKP
with presolve activated. (I cannot say for sure that this is wrong because the input set is 4D, but when using the vertex representation, the result is a proper set and all 2D projections are proper too.)I wonder if the following check for
!isapproxzero(β_primal)
is correct because the error message indicates thatβ_primal
should just be nonnegative.https://github.com/JuliaPolyhedra/Polyhedra.jl/blob/07567d48fa60e6acc0d0ed1af012eb82c82ff78c/src/linearity.jl#L180-L181
In my example I get:
(I think
3.0
is not negative :sweat_smile:)So it would be useful to get a confirmation that the check is correct. If so, I guess I have to reduce the example to get a more helpful input.