Implement block elimination

NEWS.md

@@ -0,0 +1,14 @@
+# Polyhedra.jl v0.8 Release Notes
+
+## Load time improvements
+- Dependencies on JuMP.jl, RecipesBase.jl, and GeometryBasics.jl were moved to
+  weak dependencies on Julia versions supporting package extensions, i.e. v1.9
+  and above. On v1.10 this reduces installation time by 15% and load time by
+  18% (see [#328]).
+
+## Breaking changes
+- `JuMP.optimizer_with_attributes` is no longer exported. Call it from JuMP.jl instead.
+- The following change is only breaking on Julia v1.9 and above:
+  `Polyhedra.Mesh` is now implemented in a package extension requiring
+  GeometryBasics.jl. It is sufficient to load your plotting package, i.e.
+  Makie.jl or MeshCat.jl, **before** calling `Polyhedra.Mesh`

Project.toml

@@ -6,18 +6,32 @@ version = "0.7.7"
 [deps]
 GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
-JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
+MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
 MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
+[weakdeps]
+GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
+JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
+RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+
+[extensions]
+PolyhedraGeometryBasicsExt = "GeometryBasics"
+PolyhedraJuMPExt = "JuMP"
+PolyhedraRecipesBaseExt = "RecipesBase"
+
 [compat]
 GenericLinearAlgebra = "0.2, 0.3"
 GeometryBasics = "0.2, 0.3, 0.4"
 JuMP = "0.23, 1"
+LinearAlgebra = "1.6"
+MathOptInterface = "1"
 MutableArithmetics = "1"
 RecipesBase = "0.7, 0.8, 1.0"
+SparseArrays = "1.6"
 StaticArrays = "0.12, 1.0"
 julia = "1.6"

docs/src/plot.md

@@ -63,26 +63,25 @@ Polyhedron DefaultPolyhedron{Rational{BigInt}, Polyhedra.Intersection{Rational{B
  Ray(Rational{BigInt}[0, 0, 1])
 ```
 
-Then, we need to create a mesh from the polyhedron:
-```jldoctest plots3
-julia> m = Polyhedra.Mesh(p)
-Polyhedra.Mesh{3, Rational{BigInt}, DefaultPolyhedron{Rational{BigInt}, Polyhedra.Intersection{Rational{BigInt}, Vector{Rational{BigInt}}, Int64}, Polyhedra.Hull{Rational{BigInt}, Vector{Rational{BigInt}}, Int64}}}(convexhull([0, 0, 0]) + convexhull(Ray(Rational{BigInt}[1, 0, 0]), Ray(Rational{BigInt}[0, 1, 0]), Ray(Rational{BigInt}[0, 0, 1])), nothing, nothing, nothing)
-```
-
-```@docs
-Polyhedra.Mesh
-```
-
-The polyhedron can be plotted with [MeshCat](https://github.com/rdeits/MeshCat.jl) as follows
+The polyhedron can then be plotted with [MeshCat](https://github.com/rdeits/MeshCat.jl) as follows
 ```julia
 julia> using MeshCat
 
+julia> m = Polyhedra.Mesh(p)
+
 julia> vis = Visualizer()
 
 julia> setobject!(vis, m)
 
 julia> open(vis)
 ```
+
+Note that the `Mesh` object should be created **after** loading the plotting package:
+
+```@docs
+Polyhedra.Mesh
+```
+
 To plot it in a notebook, replace `open(vis)` with `IJuliaCell(vis)`.
 
 To plot it with [Makie](https://github.com/JuliaPlots/Makie.jl) instead, you can use for instance `mesh` or `wireframe`.

examples/3D Plotting a projection of the 4D permutahedron.ipynb

@@ -77,23 +77,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "To get a plottable object, we transform the polyhedron into a mesh as follows."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 5,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "m = Polyhedra.Mesh(p3);"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "We can now plot this mesh with [MeshCat](https://github.com/rdeits/MeshCat.jl) or [Makie](https://github.com/JuliaPlots/Makie.jl) as follows:"
+    "We can now plot this polyhedron with [MeshCat](https://github.com/rdeits/MeshCat.jl) or [Makie](https://github.com/JuliaPlots/Makie.jl) as follows:"
    ]
   },
   {
@@ -129,6 +113,7 @@
    ],
    "source": [
     "using MeshCat\n",
+    "m = Polyhedra.Mesh(p3)\n",
     "vis = Visualizer()\n",
     "setobject!(vis, m)\n",
     "IJuliaCell(vis)"

examples/Polyhedral Function.ipynb

@@ -840,22 +840,14 @@
     "The top of the shape will have no face as it is unbounded in this direction."
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "m = Polyhedra.Mesh(p)"
-   ]
-  },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
     "using MeshCat\n",
+    "m = Polyhedra.Mesh(p)\n",
     "vis = Visualizer()\n",
     "setobject!(vis, m)\n",
     "IJuliaCell(vis)"

src/decompose.jl → ext/PolyhedraGeometryBasicsExt.jl

@@ -1,4 +1,10 @@
+module PolyhedraGeometryBasicsExt
+
+using LinearAlgebra
 import GeometryBasics
+using Polyhedra
+using Polyhedra: FullDim, isapproxzero, _planar_hull, counterclockwise, rotate
+using StaticArrays
 
 """
     struct Mesh{N, T, PT <: Polyhedron{T}} <: GeometryBasics.GeometryPrimitive{N, T}
@@ -29,7 +35,7 @@ function Mesh(polyhedron::Polyhedron, N::Int)
     # use polyhedron built from StaticArrays vector to avoid that.
     return Mesh{N}(polyhedron)
 end
-function Mesh(polyhedron::Polyhedron)
+function Polyhedra.Mesh(polyhedron::Polyhedron)
     return Mesh(polyhedron, FullDim(polyhedron))
 end
 
@@ -223,3 +229,5 @@ GeometryBasics.coordinates(poly::Mesh) = (fulldecompose!(poly); poly.coordinates
 GeometryBasics.faces(poly::Mesh) = (fulldecompose!(poly); poly.faces)
 GeometryBasics.texturecoordinates(poly::Mesh) = nothing
 GeometryBasics.normals(poly::Mesh) = (fulldecompose!(poly); poly.normals)
+
+end

ext/PolyhedraJuMPExt.jl

@@ -0,0 +1,53 @@
+module PolyhedraJuMPExt
+
+import JuMP
+import Polyhedra: hrep, LPHRep, polyhedron, _optimize!
+using Polyhedra: Rep, Projection, _moi_set, fulldim, dimension_names, PolyhedraToLPBridge, ProjectionBridge
+
+"""
+    hrep(model::JuMP.Model)
+
+Builds an H-representation from the feasibility set of the JuMP model `model`.
+Note that if non-linear constraint are present in the model, they are ignored.
+"""
+hrep(model::JuMP.Model) = LPHRep(model)
+LPHRep(model::JuMP.Model) = LPHRep(JuMP.backend(model))
+polyhedron(model::JuMP.Model, args...) = polyhedron(hrep(model), args...)
+_optimize!(model::JuMP.Model) = JuMP.optimize!(model)
+
+struct VariableInSet{V <: JuMP.ScalarVariable, S <: Union{Rep, Projection}} <: JuMP.AbstractVariable
+    variables::Vector{V}
+    set::S
+end
+function JuMP.build_variable(error_fun::Function, variables::Vector{<:JuMP.ScalarVariable}, set::Union{Rep, Projection})
+    if length(variables) != fulldim(set)
+        error("Number of variables ($(length(variables))) does not match the full dimension of the polyhedron ($(fulldim(set))).")
+    end
+    return VariableInSet(variables, set)
+end
+function JuMP.add_variable(model::JuMP.AbstractModel, v::VariableInSet, names)
+    dim_names = dimension_names(v.set)
+    if dim_names !== nothing
+        names = copy(names)
+        for i in eachindex(names)
+            if isempty(names[i]) && !isempty(dim_names[i])
+                names[i] = dim_names[i]
+            end
+        end
+    end
+    JuMP.add_bridge(model, PolyhedraToLPBridge)
+    JuMP.add_bridge(model, ProjectionBridge)
+    return JuMP.add_variable(model, JuMP.VariablesConstrainedOnCreation(v.variables, _moi_set(v.set)), names)
+end
+function JuMP.build_constraint(error_fun::Function, func, set::Rep)
+    return JuMP.BridgeableConstraint(
+        JuMP.build_constraint(error_fun, func, _moi_set(set)),
+        PolyhedraToLPBridge)
+end
+function JuMP.build_constraint(error_fun::Function, func, set::Projection)
+    return JuMP.BridgeableConstraint(JuMP.BridgeableConstraint(
+        JuMP.build_constraint(error_fun, func, _moi_set(set)),
+        ProjectionBridge), PolyhedraToLPBridge)
+end
+
+end

src/recipe.jl → ext/PolyhedraRecipesBaseExt.jl

@@ -1,4 +1,9 @@
+module PolyhedraRecipesBaseExt
+
+using LinearAlgebra
 using RecipesBase
+using Polyhedra
+using Polyhedra: basis, _semi_hull
 
 function planar_contour(p::Polyhedron)
     if fulldim(p) != 2
@@ -17,7 +22,7 @@ function planar_contour(p::Polyhedron)
         error("Plotting empty polyhedron is not supported.")
     end
     sort!(ps, by = x -> x[1])
-    counterclockwise(p1, p2) = dot(cross([p1; 0], [p2; 0]), [0, 0, 1])
+    counterclockwise = (p1, p2) -> dot(cross([p1; 0], [p2; 0]), [0, 0, 1])
     sweep_norm = basis(eltype(ps), fulldim(p), 1)
     top = _semi_hull(ps,  1, counterclockwise, sweep_norm)
     bot = _semi_hull(ps, -1, counterclockwise, sweep_norm)
@@ -39,3 +44,5 @@ end
     legend --> false
     planar_contour(p)
 end
+
+end

src/Polyhedra.jl

@@ -12,14 +12,14 @@ import GenericLinearAlgebra
 import MutableArithmetics
 const MA = MutableArithmetics
 
+import MathOptInterface as MOI
+import MathOptInterface.Utilities as MOIU
+
 export Polyhedron
 
 abstract type Library end
 abstract type Polyhedron{T} end
 
-using JuMP
-export optimizer_with_attributes
-
 coefficient_type(::Union{AbstractVector{T}, Type{<:AbstractVector{T}}}) where T = T
 similar_type(::Type{<:Vector}, ::Int, ::Type{T}) where T = Vector{T}
 similar_type(::Type{SparseVector{S, IT}}, ::Int, ::Type{T}) where {S, IT, T} = SparseVector{T, IT}
@@ -81,7 +81,6 @@ include("extended.jl")
 include("vecrep.jl")
 include("mixedrep.jl")
 include("lphrep.jl")
-include("jump.jl")
 include("matrep.jl")
 include("liftedrep.jl")
 include("doubledescription.jl") # FIXME move it after projection.jl once it stops depending on LiftedRep
@@ -100,7 +99,23 @@ include("projection_opt.jl")
 
 # Visualization
 include("show.jl")
-include("recipe.jl")
-include("decompose.jl")
+
+"""
+    Mesh(p::Polyhedron)
+
+Returns wrapper of a polyhedron suitable for plotting with MeshCat.jl and Makie.jl.
+
+!!! note "Extension in Julia 1.9 and above"
+    Although we require `using GeometryBasics` to use this function in Julia 1.9 and above,
+    in most use cases this extension dependency is loaded by the plotting package and no
+    further action is required.
+"""
+Mesh(p) = p isa Polyhedron ? error("this method requires using GeometryBasics") : throw(MethodError(Mesh, p))
+
+if !isdefined(Base, :get_extension)
+    include("../ext/PolyhedraJuMPExt.jl")
+    include("../ext/PolyhedraRecipesBaseExt.jl")
+    include("../ext/PolyhedraGeometryBasicsExt.jl")
+end
 
 end # module

src/center.jl

@@ -1,5 +1,4 @@
 export maximum_radius_with_center, hchebyshevcenter, vchebyshevcenter, chebyshevcenter
-using JuMP
 
 """
     maximum_radius_with_center(h::HRep, center)

src/liftedrep.jl

@@ -18,7 +18,17 @@ end
 FullDim(rep::LiftedHRepresentation) = size(rep.A, 2) - 1
 
 similar_type(::Type{LiftedHRepresentation{S, MT}}, ::FullDim, ::Type{T}) where {S, T, MT} = LiftedHRepresentation{T, similar_type(MT, T)}
-hvectortype(p::Type{LiftedHRepresentation{T, MT}}) where {T, MT} = vectortype(MT)
+hvectortype(::Type{LiftedHRepresentation{T, MT}}) where {T, MT} = vectortype(MT)
+fulltype(::Type{LiftedHRepresentation{T,MT}}) where {T,MT} = LiftedHRepresentation{T,MT}
+dualtype(::Type{LiftedHRepresentation{T,MT}}) where {T,MT} = LiftedVRepresentation{T,MT}
+function dualfullspace(::LiftedHRepresentation, d::FullDim, ::Type{T}) where {T}
+    N = fulldim(d)
+    return LiftedVRepresentation{T}(
+        [SparseArrays.spzeros(T, N) one(T) * I],
+        BitSet(1:N),
+    )
+end
+
 
 LiftedHRepresentation{T}(A::AbstractMatrix{T}, linset::BitSet=BitSet()) where {T} = LiftedHRepresentation{T, typeof(A)}(A, linset)
 LiftedHRepresentation{T}(A::AbstractMatrix, linset::BitSet=BitSet()) where {T} = LiftedHRepresentation{T}(AbstractMatrix{T}(A), linset)

src/linearity.jl

@@ -151,7 +151,6 @@ function detect_new_linearities(rep::Representation, solver; verbose=0, kws...)
     nonlins = _nonlinearity(rep)
     isempty(nonlins) && return Int[]
     model, T = layered_optimizer(solver)
-    verbose >= 2 && (_model = JuMP.direct_model(model))
     is_lin = falses(length(nonlins))
     active = trues(length(nonlins))
     # We pass `true` as we break homogeneity of the problem with `sum(λ) == 1`.
@@ -173,7 +172,7 @@ function detect_new_linearities(rep::Representation, solver; verbose=0, kws...)
     end
     for i in 1:fulldim(rep) # safer than `while ...`
         (count(is_lin) == length(is_lin) || iszero(count(active))) && break
-        verbose >= 2 && println(_model)
+        verbose >= 2 && println(model)
         # FIXME stopping when we have enough hyperplanes to prove that it's empty
         #       should also resolve the issue with presolve.
         is_feasible(model, "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.GLP_ON)`).") || break