From 8c131a6cc883c541922a8b8efe835932e8b2593f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Fri, 5 Apr 2024 16:26:29 +0200 Subject: [PATCH] Fix order of dimensions in LPHRep (#339) * Fix order of dimensions in LPHRep * Fix dimension_names * Fix * Fix * Fix * Fix --- docs/src/optimization.md | 28 ++++++++++++++++++---- src/lphrep.jl | 47 +++++++++++++++++++++++++++++++++---- test/lphrep.jl | 15 ++++++++++++ test/model_to_polyhedron.jl | 8 ++++++- 4 files changed, 89 insertions(+), 9 deletions(-) diff --git a/docs/src/optimization.md b/docs/src/optimization.md index 60e9a024..a4befbeb 100644 --- a/docs/src/optimization.md +++ b/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))) ``` diff --git a/src/lphrep.jl b/src/lphrep.jl index e6ce99ec..1d100068 100644 --- a/src/lphrep.jl +++ b/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) diff --git a/test/lphrep.jl b/test/lphrep.jl index 93f6697e..27572e4d 100644 --- a/test/lphrep.jl +++ b/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 diff --git a/test/model_to_polyhedron.jl b/test/model_to_polyhedron.jl index 08efb49f..03e7271f 100644 --- a/test/model_to_polyhedron.jl +++ b/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