#1488 - Pass solver option to prune polytopes (#1490)

* pass solver option in remove vertices

* update branch

* use import

* improve docstring and solver case

* add GLPK in extras:

* add rational default solver and cleanup functions

* refactor nothing backend

Project.toml

@@ -4,9 +4,11 @@ version = "1.13.0"
 
 [deps]
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
+GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
 GLPKMathProgInterface = "3c7084bd-78ad-589a-b5bb-dbd673274bea"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
+JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MathProgBase = "fdba3010-5040-5b88-9595-932c9decdf73"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -19,8 +21,10 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 julia = "≥ 1.0.0"
 CDDLib = "≥ 0.5.2"
 Expokit = "≥ 0.1.0"
+GLPK = "≥ 0.10.0"
 GLPKMathProgInterface = "≥ 0.4.0"
 IntervalArithmetic = "≥ 0.14.0"
+JuMP = "≥ 0.19.2"
 MathProgBase = "≥ 0.7.0"
 Optim = "≥ 0.14.1"
 Polyhedra = "≥ 0.3.4"
@@ -32,13 +36,15 @@ TaylorModels = "≥ 0.1.1"
 CDDLib = "3391f64e-dcde-5f30-b752-e11513730f60"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Expokit = "a1e7a1ef-7a5d-5822-a38c-be74e1bb89f4"
+GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
 GR = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
+JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Polyhedra = "67491407-f73d-577b-9b50-8179a7c68029"
 TaylorModels = "314ce334-5f6e-57ae-acf6-00b6e903104a"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["CDDLib", "Documenter", "Expokit", "GR", "Optim", "Plots", "Polyhedra",
+test = ["CDDLib", "Documenter", "Expokit", "GLPK", "GR", "JuMP", "Optim", "Plots", "Polyhedra",
         "TaylorModels", "Test"]

src/AbstractPolytope.jl

@@ -132,7 +132,8 @@ end
 """
     minkowski_sum(P1::AbstractPolytope{N}, P2::AbstractPolytope{N};
                   [apply_convex_hull]=true,
-                  [backend]=default_polyhedra_backend(P1, N)) where {N<:Real}
+                  [backend]=default_polyhedra_backend(P1, N),
+                  [solver]=default_lp_solver(N)) where {N<:Real}
 
 Compute the Minkowski sum between two polytopes using their vertex representation.
 
@@ -145,6 +146,8 @@ Compute the Minkowski sum between two polytopes using their vertex representatio
 - `backend`           -- (optional, default: `default_polyhedra_backend(P1, N)`)
                          the backend for polyhedral computations used to
                          post-process with a convex hull
+- `solver`            -- (optional, default: `default_lp_solver(N)`) the linear programming
+                         solver used in the backend
 
 ### Output
 
@@ -153,9 +156,12 @@ the sum of all possible sums of vertices of `P1` and `P2`.
 """
 function minkowski_sum(P1::AbstractPolytope{N}, P2::AbstractPolytope{N};
                        apply_convex_hull::Bool=true,
-                       backend=nothing) where {N<:Real}
+                       backend=nothing,
+                       solver=nothing) where {N<:Real}
+
     @assert dim(P1) == dim(P2) "cannot compute the Minkowski sum between a polyotope " *
         "of dimension $(dim(P1)) and a polytope of dimension $((dim(P2)))"
+
     vlist1, vlist2 = vertices_list(P1), vertices_list(P2)
     n, m = length(vlist1), length(vlist2)
     Vout = Vector{Vector{N}}()
@@ -168,8 +174,9 @@ function minkowski_sum(P1::AbstractPolytope{N}, P2::AbstractPolytope{N};
     if apply_convex_hull
         if backend == nothing
             backend = default_polyhedra_backend(P1, N)
+            solver = default_lp_solver(N)
         end
-        convex_hull!(Vout, backend=backend)
+        convex_hull!(Vout, backend=backend, solver=solver)
     end
     return VPolytope(Vout)
 end
@@ -190,6 +197,7 @@ using .Polyhedra: HRep, VRep,
                   hcartesianproduct, vcartesianproduct,
                   points,
                   default_library
+import JuMP, GLPK
 
 function default_polyhedra_backend(P, N::Type{<:AbstractFloat})
     return default_library(LazySets.dim(P), Float64)
@@ -198,5 +206,14 @@ end
 function default_polyhedra_backend(P, N::Type{<:Rational})
     return default_library(LazySets.dim(P), Rational{Int})
 end
+
+function default_lp_solver(N::Type{<:AbstractFloat})
+    return JuMP.with_optimizer(GLPK.Optimizer)
+end
+
+function default_lp_solver(N::Type{<:Rational})
+    return JuMP.with_optimizer(GLPK.Optimizer, method=:Exact)
+end
+
 end # quote
 end # function load_polyhedra_hpolytope()

src/VPolytope.jl

@@ -374,27 +374,49 @@ end
 
 """
     remove_redundant_vertices(P::VPolytope{N};
-                              [backend]=default_polyhedra_backend(P, N))::VPolytope{N} where {N<:Real}
+                              [backend]=nothing,
+                              [solver]=nothing)::VPolytope{N} where {N<:Real}
 
 Return the polytope obtained by removing the redundant vertices of the given polytope.
 
 ### Input
 
 - `P`       -- polytope in vertex representation
-- `backend` -- (optional, default: `default_polyhedra_backend(P1, N)`) the polyhedral
-               computations backend, see
+- `backend` -- (optional, default: `nothing`) the polyhedral
+               computations backend, see `default_polyhedra_backend(P1, N)` or
                [Polyhedra's documentation](https://juliapolyhedra.github.io/Polyhedra.jl/latest/installation.html#Getting-Libraries-1)
-               for further information
+               for further information on the available backends
+- `solver`  -- (optional, default: `nothing`) the linear programming
+               solver used in the backend, if needed; see `default_lp_solver(N)`
 
 ### Output
 
 A new polytope such that its vertices are the convex hull of the given polytope.
+
+### Notes
+
+The optimization problem associated to removing redundant vertices is handled
+by `Polyhedra`. If the polyhedral computations backend requires an LP solver but
+it has not been set, we use `default_lp_solver(N)` to define such solver.
+Otherwise, the redundancy removal function of the polyhedral backend is used.
 """
 function remove_redundant_vertices(P::VPolytope{N};
-                                   backend=default_polyhedra_backend(P, N))::VPolytope{N} where {N<:Real}
+                                   backend=nothing,
+                                   solver=nothing)::VPolytope{N} where {N<:Real}
     require(:Polyhedra; fun_name="remove_redundant_vertices")
+    if backend == nothing
+        backend = default_polyhedra_backend(P, N)
+    end
     Q = polyhedron(P; backend=backend)
-    removevredundancy!(Q)
+    if Polyhedra.supportssolver(typeof(Q))
+        if solver == nothing
+            solver = default_lp_solver(N)
+        end
+        vQ = Polyhedra.vrep(Q)
+        Polyhedra.setvrep!(Q, Polyhedra.removevredundancy(vQ, solver))
+    else
+        Polyhedra.removevredundancy!(Q)
+    end
     return VPolytope(Q)
 end
 

src/concrete_convex_hull.jl

@@ -9,7 +9,8 @@ end
 """
     convex_hull(points::Vector{VN};
                 [algorithm]=default_convex_hull_algorithm(points),
-                [backend]=nothing
+                [backend]=nothing,
+                [solver]=nothing
                 )::Vector{VN} where {N<:Real, VN<:AbstractVector{N}}
 
 Compute the convex hull of the given points.
@@ -21,6 +22,8 @@ Compute the convex hull of the given points.
                  algorithm, see valid options below
 - `backend`   -- (optional, default: `"nothing"`) polyhedral computation backend
                  for higher-dimensional point sets
+- `solver`    -- (optional, default: `"nothing"`) the linear programming solver
+                 used in the backend
 
 ### Output
 
@@ -78,14 +81,16 @@ julia> plot!(VPolygon(hull), alpha=0.2);
 """
 function convex_hull(points::Vector{VN};
                      algorithm=default_convex_hull_algorithm(points),
-                     backend=nothing
+                     backend=nothing,
+                     solver=nothing
                      )::Vector{VN} where {N<:Real, VN<:AbstractVector{N}}
-    return convex_hull!(copy(points), algorithm=algorithm, backend=backend)
+    return convex_hull!(copy(points), algorithm=algorithm, backend=backend, solver=solver)
 end
 
 function convex_hull!(points::Vector{VN};
                       algorithm=default_convex_hull_algorithm(points),
-                      backend=nothing)::Vector{VN} where {N<:Real, VN<:AbstractVector{N}}
+                      backend=nothing,
+                      solver=nothing)::Vector{VN} where {N<:Real, VN<:AbstractVector{N}}
 
     m = length(points)
 
@@ -124,7 +129,7 @@ function convex_hull!(points::Vector{VN};
         end
     else
         # general case in nd
-        return _convex_hull_nd!(points, backend=backend)
+        return _convex_hull_nd!(points, backend=backend, solver=solver)
     end
 end
 
@@ -304,13 +309,11 @@ function _convex_hull_1d!(points::Vector{VN})::Vector{VN} where {N<:Real, VN<:Ab
 end
 
 function _convex_hull_nd!(points::Vector{VN};
-                          backend=nothing
+                          backend=nothing,
+                          solver=nothing
                           )::Vector{VN} where {N<:Real, VN<:AbstractVector{N}}
     V = VPolytope(points)
-    if backend == nothing
-        backend = default_polyhedra_backend(V, N)
-    end
-    Vch = remove_redundant_vertices(V, backend=backend)
+    Vch = remove_redundant_vertices(V, backend=backend, solver=solver)
     m = length(Vch.vertices)
     points[1:m] = Vch.vertices
     return resize!(points, m)

test/unit_Polytope.jl

@@ -351,6 +351,11 @@ if test_suite_polyhedra
 
         # concrete minkowski sum
         B = convert(VPolytope, BallInf(N[0, 0, 0], N(1)))
+        X = minkowski_sum(B, B)
+        twoB = 2.0*B
+        @test X ⊆ twoB && twoB ⊆ X
+
+        # same but specifying a custom polyhedral computations backend (CDDLib)
         X = minkowski_sum(B, B, backend=CDDLib.Library())
         twoB = 2.0*B
         @test X ⊆ twoB && twoB ⊆ X