[RFC/WIP] MOI wrapper

Project.toml

@@ -4,6 +4,7 @@ authors = ["Mathieu Tanneau and contributors"]
 version = "0.1.0"
 
 [deps]
+MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]

src/MathOptPresolve.jl

@@ -1,6 +1,8 @@
 module MathOptPresolve
 
 using SparseArrays
+using MathOptInterface
+const MOI = MathOptInterface
 
 include("util.jl")
 include("status.jl")
@@ -9,4 +11,6 @@ include("problem_data.jl")
 include("solution.jl")
 include("presolve.jl")
 
+include("moi.jl")
+
 end

src/lp/empty_column.jl

@@ -5,7 +5,7 @@ struct EmptyColumn{T} <: AbstractReduction{T}
 end
 
 function remove_empty_column!(ps::PresolveData{T}, j::Int) where {T}
-    ps.pb0.is_continuous || error("Empty column routine currently only supported for LPs.")
+    # ps.pb0.is_continuous || error("Empty column routine currently only supported for LPs.")
 
     # Sanity check
     (ps.colflag[j] && (ps.nzcol[j] == 0)) || return nothing

src/lp/free_column_singleton.jl

@@ -9,7 +9,7 @@ struct FreeColumnSingleton{T} <: AbstractReduction{T}
 end
 
 function remove_free_column_singleton!(ps::PresolveData{T}, j::Int) where {T}
-    ps.pb0.is_continuous || error("Free column routine currently only supported for LPs.")
+    # ps.pb0.is_continuous || error("Free column routine currently only supported for LPs.")
 
     ps.colflag[j] && ps.nzcol[j] == 1 || return nothing  # only column singletons
 

src/lp/row_singleton.jl

@@ -8,7 +8,7 @@ end
 
 
 function remove_row_singleton!(ps::PresolveData{T}, i::Int) where {T}
-    ps.pb0.is_continuous || error("Row singleton routine currently only supported for LPs.")
+    # ps.pb0.is_continuous || error("Row singleton routine currently only supported for LPs.")
 
     # Sanity checks
     (ps.rowflag[i] && ps.nzrow[i] == 1) || return nothing

src/moi.jl

@@ -0,0 +1,259 @@
+function _postsolve_fn(ps::PresolveData)
+    return (orig_sol, transformed_sol) -> postsolve!(ps, orig_sol, transformed_sol)
+end
+
+mutable struct IntermediaryRepresentation{T<:Real}
+    I::Vector{Int}
+    J::Vector{Int}
+    V::Vector{T}
+    lcon::Vector{T}
+    ucon::Vector{T}
+    lvar::Vector{T}
+    uvar::Vector{T}
+    var_types::Vector{VariableType}
+
+    function IntermediaryRepresentation{T}(n::Int) where {T<:Real}
+        return new{T}(
+            Int[],
+            Int[],
+            T[],
+            T[],
+            T[],
+            fill(T(-Inf), n),
+            fill(T(Inf), n),
+            fill(CONTINUOUS, n),
+        )
+    end
+end
+
+num_rows(ir::IntermediaryRepresentation) = length(ir.lcon)
+
+function add_row!(ir::IntermediaryRepresentation)::Int
+    if isempty(ir.I)
+        @assert isempty(ir.J)
+        @assert isempty(ir.V)
+        @assert isempty(ir.lcon)
+        @assert isempty(ir.ucon)
+        return 1
+    else
+        row_num = I[end]
+        @assert row_num == length(lcon) == length(ucon)
+        return row_num + 1
+    end
+end
+
+function process_constraint!(
+    ir::IntermediaryRepresentation,
+    f::MOI.SingleVariable,
+    s::MOI.EqualTo{T},
+) where {T<:Real}
+    i = f.variable.value
+    ir.lvar[i] = s.value
+    ir.uvar[i] = s.value
+    return nothing
+end
+
+function process_constraint!(
+    ir::IntermediaryRepresentation,
+    f::MOI.SingleVariable,
+    s::MOI.LessThan{T},
+) where {T<:Real}
+    i = f.variable.value
+    ir.uvar[i] = s.upper
+    return nothing
+end
+
+function process_constraint!(
+    ir::IntermediaryRepresentation,
+    f::MOI.SingleVariable,
+    s::MOI.GreaterThan{T},
+) where {T<:Real}
+    i = f.variable.value
+    ir.lvar[i] = s.lower
+    return nothing
+end
+
+function process_constraint!(
+    ir::IntermediaryRepresentation,
+    f::MOI.SingleVariable,
+    s::MOI.Interval{T},
+) where {T<:Real}
+    i = f.variable.value
+    ir.lvar[i] = s.lower
+    ir.uvar[i] = s.upper
+    return nothing
+end
+
+function process_constraint!(
+    ir::IntermediaryRepresentation,
+    f::MOI.SingleVariable,
+    ::MOI.ZeroOne,
+)
+    i = f.variable.value
+    ir.var_types[i] = BINARY
+    return nothing
+end
+
+function process_constraint!(
+    ir::IntermediaryRepresentation,
+    f::MOI.SingleVariable,
+    ::MOI.Integer,
+)
+    i = f.variable.value
+    ir.var_types[i] = GENERAL_INTEGER
+    return nothing
+end
+
+_get_bounds(s::MOI.EqualTo) = (s.value, s.value)
+_get_bounds(s::MOI.LessThan{T}) where {T} = (T(-Inf), s.upper)
+_get_bounds(s::MOI.GreaterThan{T}) where {T} = (s.lower, T(Inf))
+_get_bounds(s::MOI.Interval) = (s.lower, s.upper)
+_get_bounds(s) = error("Unexpected set $s.")
+
+function process_constraint!(
+    ir::IntermediaryRepresentation,
+    f::MOI.ScalarAffineFunction{T},
+    s::MOI.AbstractScalarSet,
+) where {T<:Real}
+    row_id = add_row!(ir)
+    for term in f.terms
+        vi = term.variable_index
+        push!(ir.I, row_id)
+        push!(ir.J, term.variable_index.value)
+        push!(ir.V, term.coefficient)
+    end
+    lb, ub = _get_bounds(s)
+    push!(ir.lcon, lb - f.constant)
+    push!(ir.ucon, ub - f.constant)
+    return nothing
+end
+
+process_constraint!(ir::IntermediaryRepresentation, f, s) =
+    error("Unsupported constraint of type $(typeof(f))-in-$(typeof(s)).")
+
+function process_constraint!(
+    ir::IntermediaryRepresentation,
+    f::MOI.AbstractVectorFunction,
+    s::MOI.AbstractVectorSet,
+)
+    return process_constraint!(ir, MOIU.scalarize(f), MOIU.scalarize(s))
+end
+
+function presolve!(dest::MOI.ModelLike, src::MOI.ModelLike, T::Type{<:Real})
+    @assert MOI.is_empty(dest)
+
+    objsense = MOI.get(src, MOI.ObjectiveSense())
+    obj = MOI.get(src, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}())
+
+    n = MOI.get(src, MOI.NumberOfVariables())
+    ir = IntermediaryRepresentation{T}(n)
+    for (F, S) in MOI.get(src, MOI.ListOfConstraints())
+        for ci in MOI.get(src, MOI.ListOfConstraintIndices{F,S}())
+            f = MOI.get(src, MOI.ConstraintFunction(), ci)
+            s = MOI.get(src, MOI.ConstraintSet(), ci)
+            process_constraint!(ir, f, s)
+        end
+    end
+    obj_vec = zeros(T, n)
+    for term in obj.terms
+        obj_vec[term.variable_index.value] += term.coefficient
+    end
+    pd = ProblemData{T}()
+    load_problem!(
+        pd,
+        "model-from-MOI",
+        objsense == MOI.MIN_SENSE,
+        obj_vec,
+        MOI.constant(obj),
+        SparseArrays.dropzeros!(SparseArrays.sparse(ir.I, ir.J, ir.V, num_rows(ir), n)),
+        ir.lcon,
+        ir.ucon,
+        ir.lvar,
+        ir.uvar,
+        ir.var_types,
+    )
+    ps = PresolveData(pd)
+    presolve!(ps)
+    MOI.empty!(dest)
+    model_status = MOI.OPTIMIZE_NOT_CALLED
+    if ps.status == OPTIMAL
+        model_status = MOI.OPTIMAL
+    elseif ps.status == PRIMAL_INFEASIBLE
+        model_status = MOI.PRIMAL_INFEASIBLE
+    elseif ps.status == DUAL_INFEASIBLE
+        model_status = MOI.DUAL_INFEASIBLE
+    else
+        @assert ps.status == NOT_INFERRED
+        @assert moi_status == MOI.OPTIMIZE_NOT_CALLED
+        extract_reduced_problem!(ps)
+
+        pd = ps.pb_red
+        @assert pd.nvar > 0
+        MOI.empty!(dest)
+
+        vis = MOI.add_variables(dest, pd.nvar)
+        x = [MOI.SingleVariable(vis[j]) for j = 1:pd.nvar]
+        for j = 1:pd.nvar
+            MOI.add_constraint(dest, x[j], MOI.Interval{T}(pd.lvar[j], pd.lcon[j]))
+            if pd.var_types[j] == BINARY
+                MOI.add_constraint(dest, x[j], MOI.ZeroOne())
+            elseif pd.var_types[j] == GENERAL_INTEGER
+                MOI.add_constraint(dest, x[j], MOI.Integer())
+            end
+        end
+        MOI.set(dest, MOI.ObjectiveSense(), pd.objsense ? MOI.MIN_SENSE : MOI.MAX_SENSE)
+        MOI.set(
+            dest,
+            MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}},
+            sum(pd.obj[i] * x[i] for i = 1:pd.nvar) + pd.obj0,
+        )
+        for i = 1:pd.ncon
+            row = pd.arows[i]
+            lb, ub = pd.lcon[i], pd.ucon[i]
+            set = (
+                if lb == T(-Inf)
+                    MOI.LessThan{T}(ub)
+                elseif ub == T(Inf)
+                    MOI.GreaterThan{T}(lb)
+                elseif lb ≈ ub
+                    MOI.EqualTo{T}(lb)
+                else
+                    MOI.Interval{T}(lb, ub)
+                end
+            )
+            MOI.add_constraint(
+                dest,
+                sum(row.nzval[k] * x[row.nzind[k]] for k = 1:length(row.nzind)),
+                scalar_set,
+            )
+        end
+    end
+    return model_status, x -> _postsolve_fn(ps, x)
+end
+
+function _postsolve_fn(ps::PresolveData{T}, x::Vector{T}) where {T}
+    if ps.model_status == OPTIMAL
+        @assert ps.primal_status == FEASIBLE_POINT
+        if !isempty(x)
+            throw(ArgumentError("Presolve solved model to optimality; postsolve expects an empty input argument."))
+        end
+        return ps.solution.x
+    elseif ps.model_status == PRIMAL_INFEASIBLE
+        throw(ArgumentError("Presolve solved proven infeasible; cannot postsolve."))
+    elseif ps.model_status == DUAL_INFEASIBLE
+        @assert ps.is_primal_ray
+        # QUESTION: What do we know about ps.primal_status here?
+        return ps.solution.x
+    else
+        @assert ps.model_status == NOT_INFERRED
+        orig_sol = Solution(ps.pb0.m, ps.pb0.n)
+        trans_sol = Solution(ps.nrow, ps.ncol)
+        trans_sol.primal_status = FEASIBLE_POINT
+        trans_sol.x = x
+        postsolve!(orig_sol, trans_sol, ps)
+        @assert orig_sol.primal_status == FEASIBLE_POINT
+        @assert !ps.solution.is_primal_ray
+        @assert !ps.solution.is_dual_ray
+        return orig_sol.x
+    end
+end

src/presolve.jl

@@ -454,8 +454,10 @@ function presolve!(ps::PresolveData{T}) where {T}
         # TODO: remove column singleton with doubleton equation
 
         # Dual reductions
-        @_return_if_inferred remove_row_singletons!(ps)
-        @_return_if_inferred remove_dominated_columns!(ps)
+        if ps.pb0.is_continuous
+            @_return_if_inferred remove_row_singletons!(ps)
+            @_return_if_inferred remove_dominated_columns!(ps)
+        end
     end
 
     @_return_if_inferred remove_empty_columns!(ps)
@@ -681,7 +683,7 @@ function remove_dominated_columns!(ps::PresolveData{T}) where {T}
         iszero(aij) && continue  # empty column
 
         # Strengthen dual bounds
-        #= 
+        #=
 
         =#
         cj = ps.obj[j]

test/Project.toml

@@ -1,3 +1,4 @@
 [deps]
+MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

test/moi.jl

@@ -0,0 +1,35 @@
+function _build_model(T::Type)
+    # min  x[1]
+    # s.t. x[1] >= -1
+    #      -1 <= x[2] <= 3
+    #      x[3] <= 1.2
+    # .    x[3] in Z
+    #      2.0 x[2] == 2.5
+    model = MOIU.Model{T}()
+    n = 3
+    vis = MOI.add_variables(model, n)
+    x = [MOI.SingleVariable(vis[i]) for i = 1:n]
+    MOI.add_constraint(model, x[3], MOI.Integer())
+    MOI.add_constraint(model, x[1], MOI.GreaterThan{T}(-1.0))
+    MOI.add_constraint(model, x[2], MOI.Interval{T}(-1.0, 3.0))
+    MOI.add_constraint(model, x[3], MOI.LessThan{T}(1.2))
+    MOI.add_constraint(model, T(2.0) * x[2], MOI.EqualTo{T}(2.5))
+    # TODO: Test when objective is not set
+    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}(), x[1])
+    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
+    return model
+end
+
+@testset "presolve!" begin
+    for T in COEFF_TYPES
+        src = _build_model(T)
+        dest = MOIU.Model{T}()
+        MOP.presolve!(dest, src, T)
+        @test MOI.get(dest, MOI.NumberOfVariables()) == 0
+        @test MOI.get(dest, MOI.ListOfConstraints()) == []
+        @test MOI.get(dest, MOI.ObjectiveSense()) == MOI.FEASIBILITY_SENSE
+        obj = MOI.get(dest, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}())
+        @test obj.terms == []
+        @test obj.constant ≈ -1.0
+    end
+end

test/runtests.jl

@@ -5,6 +5,10 @@ using SparseArrays
 
 const MOP = MathOptPresolve
 
+using MathOptInterface
+const MOI = MathOptInterface
+const MOIU = MOI.Utilities
+
 const COEFF_TYPES = [Float64, BigFloat]
 
 @testset "Problem data" begin
@@ -28,3 +32,7 @@ end
         include("mip/round_integer_bounds.jl")
     end
 end
+
+@testset "MOI" begin
+    include("moi.jl")
+end