removevredundancy! ignores the specified solver

[schillic]

I have a simple polyhedron P (below). The function removevredundancy! does not use the solver associated with P (as can be seen in the stack trace) and then crashes. The reason is that solver is initialized to nothing and the line solver = default_solver(p) is only executed under some conditions. There is no other way to pass the solver to this function.

julia> P
Polyhedron DefaultPolyhedron{Float64,Polyhedra.Intersection{Float64,Array{Float64,1},Int64},Polyhedra.Hull{Float64,Array{Float64,1},Int64}}:
2-element iterator of HyperPlane{Float64,Array{Float64,1}}:
 HyperPlane([-0.0, 0.3333333333333333], 0.3333333333333333)
 HyperPlane([0.18181818181818182, -0.0], 0.18181818181818182),
4-element iterator of Polyhedra.HalfSpace{Float64,Array{Float64,1}}:
 HalfSpace([-0.0, -0.15], 0.15)
 HalfSpace([-0.13953488372093026, -0.0], 0.13953488372093026)
 HalfSpace([1.460819769243627e-17, 0.13157894736842105], 0.39473684210526316)
 HalfSpace([0.07809788267962511, -0.0], 0.23429364803887537)

julia> @which removevredundancy!(P)
removevredundancy!(p::Polyhedron; strongly) in Polyhedra at .julia/packages/Polyhedra/Wu1SI/src/redundancy.jl:59

julia> removevredundancy!(P)
Error: primal simplex failed
glp_simplex: unable to recover undefined or non-optimal solution
ERROR: Solver returned NUMERICAL_ERROR when detecting new linearity
 (you may need to activate presolve for some solvers, e.g. by replacing `GLPK.Optimizer` with
 `optimizer_with_attributes(GLPK.Optimizer, "presolve" => GLPK.ON)`). Solver specific status:
 The search was prematurely terminated due to the solver failure.
Stacktrace:
 [1] error(::String, ::MathOptInterface.TerminationStatusCode, ::String, ::String, ::String, ::String) at ./error.jl:42
 [2] _unknown_status(::MathOptInterface.Bridges.LazyBridgeOptimizer{GLPK.Optimizer}, ::MathOptInterface.TerminationStatusCode, ::String) at .julia/packages/Polyhedra/Wu1SI/src/opt.jl:130
 [3] is_feasible(::MathOptInterface.Bridges.LazyBridgeOptimizer{GLPK.Optimizer}, ::String) at .julia/packages/Polyhedra/Wu1SI/src/opt.jl:144
 [4] detect_new_linearities(::Polyhedra.Intersection{Float64,Array{Float64,1},Int64}, ::MathOptInterface.OptimizerWithAttributes; verbose::Int64) at .julia/packages/Polyhedra/Wu1SI/src/linearity.jl:163
 [5] _detect_linearity(::Polyhedra.Intersection{Float64,Array{Float64,1},Int64}, ::MathOptInterface.OptimizerWithAttributes; verbose::Int64) at .julia/packages/Polyhedra/Wu1SI/src/linearity.jl:219
 [6] detecthlinearity(::Polyhedra.Intersection{Float64,Array{Float64,1},Int64}, ::MathOptInterface.OptimizerWithAttributes; verbose::Int64) at .julia/packages/Polyhedra/Wu1SI/src/linearity.jl:239
 [7] detecthlinearity!(::DefaultPolyhedron{Float64,Polyhedra.Intersection{Float64,Array{Float64,1},Int64},Polyhedra.Hull{Float64,Array{Float64,1},Int64}}, ::MathOptInterface.OptimizerWithAttributes; verbose::Int64) at .julia/packages/Polyhedra/Wu1SI/src/linearity.jl:35
 [8] detecthlinearity!(::DefaultPolyhedron{Float64,Polyhedra.Intersection{Float64,Array{Float64,1},Int64},Polyhedra.Hull{Float64,Array{Float64,1},Int64}}, ::MathOptInterface.OptimizerWithAttributes) at .julia/packages/Polyhedra/Wu1SI/src/linearity.jl:35 (repeats 2 times)
 [9] removevredundancy!(::DefaultPolyhedron{Float64,Polyhedra.Intersection{Float64,Array{Float64,1},Int64},Polyhedra.Hull{Float64,Array{Float64,1},Int64}}; strongly::Bool) at .julia/packages/Polyhedra/Wu1SI/src/redundancy.jl:72
 [10] removevredundancy!(::DefaultPolyhedron{Float64,Polyhedra.Intersection{Float64,Array{Float64,1},Int64},Polyhedra.Hull{Float64,Array{Float64,1},Int64}}) at .julia/packages/Polyhedra/Wu1SI/src/redundancy.jl:59
 [11] top-level scope at REPL[38]:1
 [12] eval(::Module, ::Any) at ./boot.jl:331
 [13] eval_user_input(::Any, ::REPL.REPLBackend) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.4/REPL/src/REPL.jl:86
 [14] run_backend(::REPL.REPLBackend) at .julia/packages/Revise/BqeJF/src/Revise.jl:1184
 [15] top-level scope at none:0

https://github.com/JuliaPolyhedra/Polyhedra.jl/blob/485f7dfb2abb93b9032d711bc017c355796ad02e/src/redundancy.jl#L58-L79

[blegat]

IIUC, you suggest adding a keyword argument solver = default_solver(p) and pass the solver to detectvlinearity! in line 76 and 77, is that correct ?

[schillic]

That or always use solver = default_solver(p) (but maybe there is a reason for not doing that).

[blegat]

Which solver is used in your example? It seems to be GLPK, isn't it the default solver?

[schillic]

Yes, I used GLPK but the default is nothing. https://github.com/JuliaPolyhedra/Polyhedra.jl/blob/485f7dfb2abb93b9032d711bc017c355796ad02e/src/default.jl#L45-L51 https://github.com/JuliaPolyhedra/Polyhedra.jl/blob/485f7dfb2abb93b9032d711bc017c355796ad02e/src/defaultlibrary.jl#L11-L13

[blegat]

Then I don't get your comment:

The function removevredundancy! does not use the solver associated with P (as can be seen in the stack trace) and then crashes.

In the stacktrace, we see that it uses GLPK which is the solver associated with P.

[schillic]

It goes to line 72 which is the solver === nothing branch. The error message say that presolve was not activated but I think the solver I passed had it activated. I am not 100% sure but I guess in the function in line 72 it uses a different solver.

[blegat]

It is using default_solver too https://github.com/JuliaPolyhedra/Polyhedra.jl/blob/485f7dfb2abb93b9032d711bc017c355796ad02e/src/linearity.jl#L35

[schillic]

But that one does not have presolve activated.

[blegat]

Then that means that p.solver does not have presolve. There is no default solver in Polyhedra if p.solver is nothing so if a solver is used, the only possibility is that you set it to the library.

[schillic]

~I tried to reproduce what I did but the behavior is not the same as in the OP. I guess I did not pass the solver to the correct polyhedron :sweat_smile:~ Now I could reproduce.

[schillic]

julia> using Polyhedra, JuMP, GLPK
julia> v1 = [[1.0, 1.0], [-1.0, 1.0], [1.0, -1.0], [-1.0, -1.0]];
julia> v2 = [[3.0, 3.0], [1.0, 3.0], [1.0, 1.0], [3.0, 1.0]];
julia> solver = Polyhedra.DefaultLibrary{Float64}(
    JuMP.optimizer_with_attributes(GLPK.Optimizer, "presolve" => GLPK.ON));
julia> p1 = polyhedron(Polyhedra.vrep(v1), solver);
julia> p2 = polyhedron(Polyhedra.vrep(v2), solver);
julia> p = Polyhedra.intersect(p1, p2);

julia> Polyhedra.default_solver(p1)
MathOptInterface.OptimizerWithAttributes(GLPK.Optimizer,
 Pair{MathOptInterface.AbstractOptimizerAttribute,Any}[MathOptInterface.RawParameter("presolve") => 1])

julia> Polyhedra.default_solver(p2)
MathOptInterface.OptimizerWithAttributes(GLPK.Optimizer,
 Pair{MathOptInterface.AbstractOptimizerAttribute,Any}[MathOptInterface.RawParameter("presolve") => 1])

julia> Polyhedra.default_solver(p)
MathOptInterface.OptimizerWithAttributes(GLPK.Optimizer,
 Pair{MathOptInterface.AbstractOptimizerAttribute,Any}[MathOptInterface.RawParameter("presolve") => 1])

julia> removevredundancy!(p);
*crash from OP*

[blegat]

Maybe the solver encounters numerical errors even if the presolve is on ? The warning is just a hint that might be useful, it does not mean that the presolve is off.

[schillic]

Maybe the solver encounters numerical errors even if the presolve is on ? The warning is just a hint that might be useful, it does not mean that the presolve is off.

Indeed, you are right, so the original discussion is wrong. Still, this example should somehow work as it is not particularly difficult. Geometrically it is the intersection of two squares in a single vertex - clearly a corner case but doable.

After replacing line 59 in

https://github.com/JuliaPolyhedra/Polyhedra.jl/blob/485f7dfb2abb93b9032d711bc017c355796ad02e/src/redundancy.jl#L58-L79

by solver = default_solver(p), it does not crash anymore. That is why I thought that passing the solver would be the solution. But this is not really the solution as it just skips the call to detecthlinearity!(p), which is the real source.

The delicate aspect here is that I can call detecthlinearity!(p) and detectvlinearity!(p) individually without a problem. But after calling detecthlinearity!(p) once, the second time it crashes. If I do not call detectvlinearity!(p) first, then that call will also crash (which is what happens in the OP; if it had been called first, I guess the second call ignores the H-representation and works on the V-representation).

Conclusion

My understanding is that detecthlinearity!(p) can make a stable representation unstable. Not sure if we can do anything about that.

Complete example

julia> using Polyhedra, JuMP, GLPK
julia> v1 = [[1.0, 1.0], [-1.0, 1.0], [1.0, -1.0], [-1.0, -1.0]];
julia> v2 = [[3.0, 3.0], [1.0, 3.0], [1.0, 1.0], [3.0, 1.0]];
julia> solver = Polyhedra.DefaultLibrary{Float64}(
    JuMP.optimizer_with_attributes(GLPK.Optimizer, "presolve" => GLPK.ON));
julia> p1 = polyhedron(Polyhedra.vrep(v1), solver);
julia> p2 = polyhedron(Polyhedra.vrep(v2), solver);
julia> p = Polyhedra.intersect(p1, p2)
Polyhedron DefaultPolyhedron{Float64,Polyhedra.Intersection{Float64,Array{Float64,1},Int64},Polyhedra.Hull{Float64,Array{Float64,1},Int64}}:
8-element iterator of HalfSpace{Float64,Array{Float64,1}}:
 HalfSpace([-0.0, 0.3333333333333333], 0.3333333333333333)
 HalfSpace([0.18181818181818182, -0.0], 0.18181818181818182)
 HalfSpace([-0.0, -0.15], 0.15)
 HalfSpace([-0.13953488372093026, -0.0], 0.13953488372093026)
 HalfSpace([1.460819769243627e-17, 0.13157894736842105], 0.39473684210526316)
 HalfSpace([-0.14251425142514254, 3.560008573561765e-17], -0.1425142514251425)
 HalfSpace([0.07809788267962511, -0.0], 0.23429364803887537)
 HalfSpace([-1.0367150553151396e-17, -0.0818244458905945], -0.08182444589059451)

julia> detectvlinearity!(p)  # works fine
V-representation Polyhedra.Hull{Float64,Array{Float64,1},Int64}:
1-element iterator of Array{Float64,1}:
 [1.0, 1.0]

julia> detecthlinearity!(p)  # works fine but changes the constraints
H-representation Polyhedra.Intersection{Float64,Array{Float64,1},Int64}:
2-element iterator of HyperPlane{Float64,Array{Float64,1}}:
 HyperPlane([-0.0, 0.3333333333333333], 0.3333333333333333)
 HyperPlane([0.18181818181818182, -0.0], 0.18181818181818182),
4-element iterator of HalfSpace{Float64,Array{Float64,1}}:
 HalfSpace([-0.0, -0.15], 0.15)
 HalfSpace([-0.13953488372093026, -0.0], 0.13953488372093026)
 HalfSpace([1.460819769243627e-17, 0.13157894736842105], 0.39473684210526316)
 HalfSpace([0.07809788267962511, -0.0], 0.23429364803887537)

julia> detecthlinearity!(p)  # second call crashes