Fix order of dimensions in LPHRep

docs/src/optimization.md

@@ -122,9 +122,29 @@ we are either in the following three situations:
 - An extreme point is optimal.
 
 A JuMP model is treated by `polyhedron` just like any H-representation. For example, the hypercube of dimension `n` can be created as follows
+```@example lphrep
+using JuMP, Polyhedra
+model = Model()
+@variable(model, 0 ≤ x[1:2] ≤ 1)
+h = hrep(model)
+```
+The name of the variables for each dimension can be recovered as follows
+```@example lphrep
+dimension_names(h)
+```
+Note that the names of the JuMP variables are lost in the conversion to a
+polyhedron that does not support names, e.g.,
 ```julia
-m = Model()
-@variable(m, 0 ≤ x[1:n] ≤ 1)
-
-poly = polyhedron(m, CDDLib.Library(:exact))
+poly = polyhedron(model, CDDLib.Library(:exact))
+```
+However, the ordering of the dimension of the polyhedron is guaranteed to
+correspond to the order of the JuMP variables as listed by `all_variables`:
+```@example lphrep
+all_variables(model)
+```
+So `all_variables(model)[i]` provides the JuMP variable corresponding to the `i`th
+dimension.
+The reverse mapping can be constructed as follows:
+```@example lphrep
+Dict(v => i for (i, v) in enumerate(all_variables(model)))
 ```

src/lphrep.jl

@@ -5,8 +5,8 @@ MOI.Utilities.@model(_MOIModel,
                      (), (MOI.ScalarAffineFunction,), (), ())
 # We need the `VariableIndex` constraints to be bridged so we should say that
 # they are not supported. We notably exclude `Integer` as we just ignore
-# integrality constraints. Binary constraint should be bridged to integrality
-# once https://github.com/jump-dev/MathOptInterface.jl/issues/704 is done.
+# integrality constraints. Binary constraint are bridged to integrality
+# with `ZeroOneBridge`.
 function MOI.supports_constraint(
     ::_MOIModel{T}, ::Type{MOI.VariableIndex},
     ::Type{<:Union{MOI.EqualTo{T}, MOI.GreaterThan{T}, MOI.LessThan{T},
@@ -15,6 +15,7 @@ function MOI.supports_constraint(
 end
 
 
+
 mutable struct LPHRep{T} <: HRepresentation{T}
     model::_MOIModel{T}
     hyperplane_indices::Union{Nothing,
@@ -27,6 +28,15 @@ end
 function LPHRep(model::_MOIModel{T}) where T
     return LPHRep(model, nothing, nothing)
 end
+
+function _pass_constraints(dest, src, index_map, cis_src)
+    for ci_src in cis_src
+        f = MOI.get(src, MOI.ConstraintFunction(), ci_src)
+        s = MOI.get(src, MOI.ConstraintSet(), ci_src)
+        MOI.add_constraint(dest, MOI.Utilities.map_indices(index_map, f), s)
+    end
+end
+
 function LPHRep(model::MOI.ModelLike, T::Type = Float64)
     _model = _MOIModel{T}()
     bridged = MOI.Bridges.LazyBridgeOptimizer(_model)
@@ -42,12 +52,41 @@ function LPHRep(model::MOI.ModelLike, T::Type = Float64)
     MOI.Bridges.add_bridge(bridged, MOI.Bridges.Constraint.ScalarFunctionizeBridge{T})
     MOI.Bridges.add_bridge(bridged, MOI.Bridges.Constraint.VectorFunctionizeBridge{T})
     MOI.Bridges.add_bridge(bridged, MOI.Bridges.Constraint.SplitIntervalBridge{T})
+    MOI.Bridges.add_bridge(bridged, MOI.Bridges.Constraint.ZeroOneBridge{T})
     MOI.Bridges.add_bridge(bridged, MOI.Bridges.Constraint.NormInfinityBridge{T})
     # This one creates variables so the user need to consider the polyhedra as the
     # feasible set of the extended formulation.
     MOI.Bridges.add_bridge(bridged, MOI.Bridges.Constraint.NormOneBridge{T})
     MOI.Bridges.add_bridge(bridged, PolyhedraToLPBridge{T})
-    MOI.copy_to(bridged, model)
+    # If we call `MOI.copy_to` here, it will redirect to
+    # `MOI.Utilities.default_copy_to` which will copy constrained
+    # variables first.
+    # This might create a mismatch between the order of the dimensions
+    # and `JuMP.all_variables`. To avoid this, we copy it ourself.
+    # Note that this also prevents issue with errors of attributes not
+    # supported that would be thrown by `MOI.copy_to`.
+    vis_src = MOI.get(model, MOI.ListOfVariableIndices())
+    index_map = MOI.Utilities.IndexMap()
+    vis_dst = MOI.add_variables(bridged, length(vis_src))
+    attr = MOI.VariableName()
+    has_names = attr in MOI.get(model, MOI.ListOfVariableAttributesSet())
+    for (vi_src, vi_dst) in zip(vis_src, vis_dst)
+        index_map[vi_src] = vi_dst
+        if has_names
+            name = MOI.get(model, attr, vi_src)
+            if !isnothing(name)
+                MOI.set(bridged, attr, vi_dst, name)
+            end
+        end
+    end
+    for (F, S) in MOI.get(model, MOI.ListOfConstraintTypesPresent())
+        _pass_constraints(
+            bridged,
+            model,
+            index_map,
+            MOI.get(model, MOI.ListOfConstraintIndices{F,S}()),
+        )
+    end
     return LPHRep(_model)
 end
 # returns `Int64` so need to convert for 32-bit system
@@ -197,4 +236,4 @@ function Base.isvalid(lp::LPHRep{T}, idx::HIndex{T}) where {T}
 end
 
 dualtype(::Type{LPHRep{T}}, ::Type{AT}) where {T, AT} = dualtype(Intersection{T, AT, Int}, AT)
-dualfullspace(h::LPHRep, d::FullDim, ::Type{T}, ::Type{AT}) where {T, AT} = dualfullspace(Intersection{T, AT, Int}, d, T, AT)
+dualfullspace(::LPHRep, d::FullDim, ::Type{T}, ::Type{AT}) where {T, AT} = dualfullspace(Intersection{T, AT, Int}, d, T, AT)

test/lphrep.jl

@@ -10,3 +10,18 @@ function test_lphrep()
     @test nhyperplanes(lph) == 1
     @test first(hyperplanes(lph)) == hp
 end
+
+# Test that the order of variables is kept, see
+# https://github.com/JuliaPolyhedra/Polyhedra.jl/issues/338
+# `MOI.Utilities.default_copy_to` will copy `y` first
+# as it is a constrained variable. This tests that we do
+# our custom `copy` that preserves order.
+function test_order_lphrep()
+    model = Model()
+    @variable(model, x)
+    @variable(model, y <= 1)
+    h = hrep(model)
+    @show h
+    @test nhalfspaces(h) == 1
+    @test first(halfspaces(h)) == HalfSpace([0, 1], 1)
+end

test/model_to_polyhedron.jl

@@ -19,5 +19,11 @@ using Polyhedra
         @test p isa DefaultPolyhedron{Float64}
         @test nhalfspaces(p) == 3
     end
-    # TODO add test with binary variables once https://github.com/jump-dev/MathOptInterface.jl/issues/704 is done.
+    @testset "Binary" begin
+        model = Model()
+        @variable(model, x, Bin)
+        p = polyhedron(model)
+        @test p isa Interval{Float64}
+        @test nhalfspaces(p) == 2
+    end
 end