From ec41a4727ee0ca5b78ad0ce9d848a56d8b5e0a2a Mon Sep 17 00:00:00 2001 From: odow Date: Tue, 3 Nov 2020 16:32:10 +1300 Subject: [PATCH 1/2] Add farkas duals for variable bounds and refactor how we compute primal/dual rays. --- Project.toml | 2 +- src/MOI_wrapper/MOI_wrapper.jl | 122 ++++++++------ src/MOI_wrapper/infeasibility_certificates.jl | 153 ++++++++---------- test/MOI_wrapper.jl | 146 +++++++++++++++-- 4 files changed, 277 insertions(+), 146 deletions(-) diff --git a/Project.toml b/Project.toml index 108627f..bc98aec 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "GLPK" uuid = "60bf3e95-4087-53dc-ae20-288a0d20c6a6" repo = "https://github.com/jump-dev/GLPK.jl.git" -version = "0.14.2" +version = "0.14.3" [deps] BinaryProvider = "b99e7846-7c00-51b0-8f62-c81ae34c0232" diff --git a/src/MOI_wrapper/MOI_wrapper.jl b/src/MOI_wrapper/MOI_wrapper.jl index ce8e729..a26f1c7 100644 --- a/src/MOI_wrapper/MOI_wrapper.jl +++ b/src/MOI_wrapper/MOI_wrapper.jl @@ -389,7 +389,7 @@ function _indices_and_coefficients( ) i = 1 for term in f.terms - indices[i] = Cint(_info(model, term.variable_index).column) + indices[i] = Cint(column(model, term.variable_index)) coefficients[i] = term.coefficient i += 1 end @@ -422,6 +422,8 @@ function _info(model::Optimizer, key::MOI.VariableIndex) throw(MOI.InvalidIndex(key)) end +column(model, x::MOI.VariableIndex) = _info(model, x).column + function MOI.add_variable(model::Optimizer) # Initialize `VariableInfo` with a dummy `VariableIndex` and a column, # because we need `add_item` to tell us what the `VariableIndex` is. @@ -569,8 +571,8 @@ function MOI.set( num_vars = length(model.variable_info) obj = zeros(Float64, num_vars) for term in f.terms - column = _info(model, term.variable_index).column - obj[column] += term.coefficient + col = column(model, term.variable_index) + obj[col] += term.coefficient end for (col, coef) in enumerate(obj) glp_set_obj_coef(model, col, coef) @@ -619,6 +621,10 @@ function _info( return throw(MOI.InvalidIndex(c)) end +function column(model, c::MOI.ConstraintIndex{MOI.SingleVariable, <:Any}) + return _info(model, c).column +end + function MOI.is_valid( model::Optimizer, c::MOI.ConstraintIndex{MOI.SingleVariable, MOI.LessThan{Float64}} @@ -862,7 +868,7 @@ function MOI.get( c::MOI.ConstraintIndex{MOI.SingleVariable, MOI.GreaterThan{Float64}}, ) MOI.throw_if_not_valid(model, c) - lower = glp_get_col_lb(model, _info(model, c).column) + lower = glp_get_col_lb(model, column(model, c)) return MOI.GreaterThan(lower) end @@ -872,7 +878,7 @@ function MOI.get( c::MOI.ConstraintIndex{MOI.SingleVariable, MOI.LessThan{Float64}}, ) MOI.throw_if_not_valid(model, c) - upper = glp_get_col_ub(model, _info(model, c).column) + upper = glp_get_col_ub(model, column(model, c)) return MOI.LessThan(upper) end @@ -882,7 +888,7 @@ function MOI.get( c::MOI.ConstraintIndex{MOI.SingleVariable, MOI.EqualTo{Float64}}, ) MOI.throw_if_not_valid(model, c) - lower = glp_get_col_lb(model, _info(model, c).column) + lower = glp_get_col_lb(model, column(model, c)) return MOI.EqualTo(lower) end @@ -892,9 +898,9 @@ function MOI.get( c::MOI.ConstraintIndex{MOI.SingleVariable, MOI.Interval{Float64}}, ) MOI.throw_if_not_valid(model, c) - info = _info(model, c) - lower = glp_get_col_lb(model, info.column) - upper = glp_get_col_ub(model, info.column) + col = column(model, c) + lower = glp_get_col_lb(model, col) + upper = glp_get_col_ub(model, col) return MOI.Interval(lower, upper) end @@ -906,8 +912,7 @@ function MOI.set( ) where {S<:_SCALAR_SETS} MOI.throw_if_not_valid(model, c) lower, upper = _bounds(s) - info = _info(model, c) - _set_variable_bound(model, info.column, lower, upper) + _set_variable_bound(model, column(model, c), lower, upper) return end @@ -1411,14 +1416,18 @@ function MOI.optimize!(model::Optimizer) else _solve_linear_problem(model) end - if MOI.get(model, MOI.ResultCount()) > 0 - if MOI.get(model, MOI.PrimalStatus()) == MOI.INFEASIBILITY_CERTIFICATE - model.unbounded_ray = fill(NaN, glp_get_num_cols(model)) - _get_unbounded_ray(model, model.unbounded_ray) - end - if MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE - model.infeasibility_cert = fill(NaN, glp_get_num_rows(model)) - _get_infeasibility_ray(model, model.infeasibility_cert) + if _certificates_potentially_available(model) + (status, _) = _get_status(model) + if status == MOI.DUAL_INFEASIBLE + ray = zeros(glp_get_num_cols(model)) + if _get_unbounded_ray(model, ray) + model.unbounded_ray = ray + end + elseif status == MOI.INFEASIBLE + ray = zeros(glp_get_num_rows(model)) + if _get_infeasibility_ray(model, ray) + model.infeasibility_cert = ray + end end end model.solve_time = time() - start_time @@ -1561,7 +1570,7 @@ function MOI.get(model::Optimizer, attr::MOI.PrimalStatus) elseif status == MOI.LOCALLY_INFEASIBLE return MOI.INFEASIBLE_POINT elseif status == MOI.DUAL_INFEASIBLE - if _certificates_potentially_available(model) + if model.unbounded_ray !== nothing return MOI.INFEASIBILITY_CERTIFICATE end else @@ -1579,7 +1588,7 @@ function MOI.get(model::Optimizer, attr::MOI.DualStatus) if status == MOI.OPTIMAL return MOI.FEASIBLE_POINT elseif status == MOI.INFEASIBLE - if _certificates_potentially_available(model) + if model.infeasibility_cert !== nothing return MOI.INFEASIBILITY_CERTIFICATE end end @@ -1622,9 +1631,9 @@ function MOI.get(model::Optimizer, attr::MOI.VariablePrimal, x::MOI.VariableInde _throw_if_optimize_in_progress(model, attr) MOI.check_result_index_bounds(model, attr) if model.unbounded_ray !== nothing - return model.unbounded_ray[_info(model, x).column] + return model.unbounded_ray[column(model, x)] else - return _get_col_primal(model, _info(model, x).column) + return _get_col_primal(model, column(model, x)) end end @@ -1650,18 +1659,33 @@ function _dual_multiplier(model::Optimizer) return MOI.get(model, MOI.ObjectiveSense()) == MOI.MIN_SENSE ? 1.0 : -1.0 end +function _farkas_variable_dual(model::Optimizer, col::Cint) + nnz = glp_get_mat_col(model, col, C_NULL, C_NULL) + vind = Vector{Cint}(undef, nnz) + vval = Vector{Cdouble}(undef, nnz) + nnz = glp_get_mat_col(model, col, offset(vind), offset(vval)) + return sum( + model.infeasibility_cert[row] * val for (row, val) in zip(vind, vval) + ) +end + function MOI.get( - model::Optimizer, attr::MOI.ConstraintDual, - c::MOI.ConstraintIndex{MOI.SingleVariable, MOI.LessThan{Float64}} + model::Optimizer, + attr::MOI.ConstraintDual, + c::MOI.ConstraintIndex{MOI.SingleVariable, MOI.LessThan{Float64}}, ) _throw_if_optimize_in_progress(model, attr) MOI.check_result_index_bounds(model, attr) - column = _info(model, c).column + if model.infeasibility_cert !== nothing + dual = _farkas_variable_dual(model, Cint(column(model, c))) + return min(dual, 0.0) + end + col = column(model, c) reduced_cost = if model.method == SIMPLEX || model.method == EXACT - glp_get_col_dual(model, column) + glp_get_col_dual(model, col) else @assert model.method == INTERIOR - glp_ipt_col_dual(model, column) + glp_ipt_col_dual(model, col) end sense = MOI.get(model, MOI.ObjectiveSense()) # The following is a heuristic for determining whether the reduced cost @@ -1683,17 +1707,22 @@ function MOI.get( end function MOI.get( - model::Optimizer, attr::MOI.ConstraintDual, - c::MOI.ConstraintIndex{MOI.SingleVariable, MOI.GreaterThan{Float64}} + model::Optimizer, + attr::MOI.ConstraintDual, + c::MOI.ConstraintIndex{MOI.SingleVariable, MOI.GreaterThan{Float64}}, ) _throw_if_optimize_in_progress(model, attr) MOI.check_result_index_bounds(model, attr) - column = _info(model, c).column + if model.infeasibility_cert !== nothing + dual = _farkas_variable_dual(model, Cint(column(model, c))) + return max(dual, 0.0) + end + col = column(model, c) reduced_cost = if model.method == SIMPLEX || model.method == EXACT - glp_get_col_dual(model, column) + glp_get_col_dual(model, col) else @assert model.method == INTERIOR - glp_ipt_col_dual(model, column) + glp_ipt_col_dual(model, col) end sense = MOI.get(model, MOI.ObjectiveSense()) # The following is a heuristic for determining whether the reduced cost @@ -1715,23 +1744,28 @@ function MOI.get( end function MOI.get( - model::Optimizer, attr::MOI.ConstraintDual, - c::MOI.ConstraintIndex{MOI.SingleVariable, S} + model::Optimizer, + attr::MOI.ConstraintDual, + c::MOI.ConstraintIndex{MOI.SingleVariable, S}, ) where {S <: Union{MOI.EqualTo, MOI.Interval}} _throw_if_optimize_in_progress(model, attr) MOI.check_result_index_bounds(model, attr) - return _get_col_dual(model, _info(model, c).column) + if model.infeasibility_cert !== nothing + return _farkas_variable_dual(model, Cint(column(model, c))) + end + return _get_col_dual(model, column(model, c)) end function MOI.get( - model::Optimizer, attr::MOI.ConstraintDual, - c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, <:Any} + model::Optimizer, + attr::MOI.ConstraintDual, + c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, <:Any}, ) _throw_if_optimize_in_progress(model, attr) MOI.check_result_index_bounds(model, attr) row = _info(model, c).row if model.infeasibility_cert !== nothing - return model.infeasibility_cert[row] + return -model.infeasibility_cert[row] else @assert !model.last_solved_by_mip if model.method == SIMPLEX || model.method == EXACT @@ -1911,7 +1945,7 @@ function MOI.modify( chg::MOI.ScalarCoefficientChange{Float64} ) row = Cint(_info(model, c).row) - col = _info(model, chg.variable).column + col = column(model, chg.variable) nnz = glp_get_mat_row(model, row, C_NULL, C_NULL) indices, coefficients = zeros(Cint, nnz), zeros(Cdouble, nnz) glp_get_mat_row(model, row, offset(indices), offset(coefficients)) @@ -1933,9 +1967,7 @@ function MOI.modify( c::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}, chg::MOI.ScalarCoefficientChange{Float64} ) - glp_set_obj_coef( - model, _info(model, chg.variable).column, chg.new_coefficient - ) + glp_set_obj_coef(model, column(model, chg.variable), chg.new_coefficient) return end @@ -1979,8 +2011,8 @@ function MOI.get( c::MOI.ConstraintIndex{MOI.SingleVariable, S}, ) where {S <: _SCALAR_SETS} _throw_if_optimize_in_progress(model, attr) - column = _info(model, c).column - vbasis = glp_get_col_stat(model, column) + col = column(model, c) + vbasis = glp_get_col_stat(model, col) if vbasis == GLP_BS return MOI.BASIC elseif vbasis == GLP_NL diff --git a/src/MOI_wrapper/infeasibility_certificates.jl b/src/MOI_wrapper/infeasibility_certificates.jl index 9a5e47a..eec2a71 100644 --- a/src/MOI_wrapper/infeasibility_certificates.jl +++ b/src/MOI_wrapper/infeasibility_certificates.jl @@ -1,5 +1,6 @@ -# The functions getinfeasibilityray and getunboundedray are adapted from code -# taken from the LEMON C++ optimization library. This is the copyright notice: +# The functions _get_infeasibility_ray and _get_unbounded_ray are adapted from +# code taken from the LEMON C++ optimization library. This is the copyright +# notice: # ### Copyright (C) 2003-2010 ### Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport @@ -16,112 +17,84 @@ """ _get_infeasibility_ray(model::Optimizer, ray::Vector{Float64}) -Get the Farkas certificate of primal infeasiblity. +Compute a Farkas certificate of primal infeasibility (unbounded dual ray) and +store in `ray`. Returns `true` if successful, and `false` if a ray cannot be +found. -Can only be called when glp_simplex is used as the solver. +Assumes `ray` has been initialized to all `0.0`s. """ function _get_infeasibility_ray(model::Optimizer, ray::Vector{Float64}) - num_rows = glp_get_num_rows(model) - @assert length(ray) == num_rows - ur = glp_get_unbnd_ray(model) - if ur != 0 - if ur <= num_rows - k = ur - status = glp_get_row_stat(model, k) - bind = glp_get_row_bind(model, k) - primal = glp_get_row_prim(model, k) - upper_bound = glp_get_row_ub(model, k) - else - k = ur - num_rows - status = glp_get_col_stat(model, k) - bind = glp_get_col_bind(model, k) - primal = glp_get_col_prim(model, k) - upper_bound = glp_get_col_ub(model, k) - end - if status != GLP_BS - error("unbounded ray is primal (use getunboundedray)") - end - ray[bind] = (primal > upper_bound) ? -1 : 1 - glp_btran(model, ray) + m = glp_get_num_rows(model) + n = glp_get_num_cols(model) + @assert length(ray) == m + # Solve with dual simplex to find unbounded ray. + param = glp_smcp() + glp_init_smcp(param) + param.msg_lev = GLP_MSG_ERR + param.meth = GLP_DUAL + status = glp_simplex(model, param) + if status != 0 || glp_get_status(model) != GLP_NOFEAS + return false # Something went wrong finding an unbounded ray. + end + unbounded_index = glp_get_unbnd_ray(model) + if unbounded_index == 0 + return false # Something went wrong finding an unbounded ray. + end + primal = if unbounded_index <= m + glp_get_row_prim(model, unbounded_index) else - eps = 1e-7 - # We need to factorize here because sometimes GLPK will prove - # infeasibility before it has a factorized basis in memory. - glp_factorize(model) - for row in 1:num_rows - idx = glp_get_bhead(model, row) - if idx <= num_rows - k = idx - primal = glp_get_row_prim(model, k) - upper_bound = glp_get_row_ub(model, k) - lower_bound = glp_get_row_lb(model, k) - else - k = idx - num_rows - primal = glp_get_col_prim(model, k) - upper_bound = glp_get_col_ub(model, k) - lower_bound = glp_get_col_lb(model, k) - end - if primal > upper_bound + eps - ray[row] = -1 - elseif primal < lower_bound - eps - ray[row] = 1 - else - continue # ray[row] == 0 - end - if idx <= num_rows - ray[row] *= glp_get_rii(model, k) - else - ray[row] /= glp_get_sjj(model, k) - end - end - glp_btran(model, ray) - for row in 1:num_rows - ray[row] /= glp_get_rii(model, row) + glp_get_col_prim(model, unbounded_index - m) + end + scale = xor(glp_get_obj_dir(model) == GLP_MAX, primal > 0) ? -1 : 1 + if unbounded_index <= m + ray[unbounded_index] = scale + end + nnz = m + n + vind = Vector{Cint}(undef, nnz) + vval = Vector{Cdouble}(undef, nnz) + len = glp_eval_tab_row(model, unbounded_index, offset(vind), offset(vval)) + for i = 1:len + if vind[i] <= m + ray[vind[i]] = -scale * vval[i] end end + return true end """ - _get_unbounded_ray(model::Optimizer, ray::Vector{Float64}) + _get_unbounded_ray(model::Optimizer, ray::Vector{Float64})::Bool -Get the certificate of primal unboundedness. +Compute an unbounded primal ray and store in `ray`. Returns `true` if +successful, and `false` if a ray cannot be found. -Can only be called when glp_simplex is used as the solver. +Assumes the primal has been solved with primal simplex and is proven unbounded. +Assumes `ray` has been initialized to all `0.0`s. """ function _get_unbounded_ray(model::Optimizer, ray::Vector{Float64}) - num_rows = glp_get_num_rows(model) + m = glp_get_num_rows(model) n = glp_get_num_cols(model) @assert length(ray) == n - ur = glp_get_unbnd_ray(model) - if ur <= num_rows - k = ur - status = glp_get_row_stat(model, k) - dual = glp_get_row_dual(model, k) + unbounded_index = glp_get_unbnd_ray(model) + if unbounded_index == 0 + return false # Something went wrong finding an unbounded ray. + end + dual = if unbounded_index <= m + glp_get_row_dual(model, unbounded_index) else - k = ur - num_rows - status = glp_get_col_stat(model, k) - dual = glp_get_col_dual(model, k) - ray[k] = 1 + glp_get_col_dual(model, unbounded_index - m) end - if status == GLP_BS - error("unbounded ray is dual (use _get_infeasibility_ray)") + scale = xor(glp_get_obj_dir(model) == GLP_MAX, dual > 0) ? -1 : 1 + if unbounded_index > m + ray[unbounded_index - m] = scale end - nnz = n + num_rows - indices, coefficients = zeros(Cint, nnz), zeros(Cdouble, nnz) - len = glp_eval_tab_col( - model, - nnz, - pointer(indices) - sizeof(Cint), - pointer(coefficients) - sizeof(Cdouble), - ) - splice!(indices, (len+1):nnz) - splice!(coefficients, (len+1):nnz) - for (row, coef) in zip(indices, coefficients) - if row > num_rows - ray[row - num_rows] = coef + nnz = m + n + vind = Vector{Cint}(undef, nnz) + vval = Vector{Cdouble}(undef, nnz) + len = glp_eval_tab_col(model, unbounded_index, offset(vind), offset(vval)) + for i = 1:len + if vind[i] > m + ray[vind[i] - m] = scale * vval[i] end end - if xor(glp_get_obj_dir(model) == GLP_MAX, dual > 0) - ray *= -1 - end + return true end diff --git a/test/MOI_wrapper.jl b/test/MOI_wrapper.jl index 3870e1b..29d1ea3 100644 --- a/test/MOI_wrapper.jl +++ b/test/MOI_wrapper.jl @@ -44,15 +44,10 @@ end @testset "Linear tests" begin @testset "Default Solver" begin MOIT.contlineartest(OPTIMIZER, MOIT.TestConfig(basis = true), [ - # This requires an infeasiblity certificate for a variable bound. - "linear12", # VariablePrimalStart not supported. - "partial_start" + "partial_start", ]) end - @testset "No certificate" begin - MOIT.linear12test(OPTIMIZER, MOIT.TestConfig(infeas_certificates=false)) - end end @testset "Linear Conic tests" begin @@ -337,14 +332,23 @@ end [MOI.ScalarAffineTerm(1.0, x[1]), MOI.ScalarAffineTerm(1.0, x[2])], 0.0 ), - MOI.EqualTo(1.0) + MOI.LessThan(1.0) + ) + c2 = MOI.add_constraint(model, MOI.SingleVariable(x[1]), MOI.EqualTo(1.0)) + c3 = MOI.add_constraint(model, MOI.SingleVariable(x[2]), MOI.EqualTo(1.0)) + MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE) + MOI.set( + model, + MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), + MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, x), 0.0) ) - MOI.add_constraint(model, MOI.SingleVariable(x[1]), MOI.EqualTo(1.0)) - MOI.add_constraint(model, MOI.SingleVariable(x[2]), MOI.EqualTo(1.0)) MOI.optimize!(model) @test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE @test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE - @test MOI.get(model, MOI.ConstraintDual(), c1) == -1 + cd1 = MOI.get(model, MOI.ConstraintDual(), c1) + @test cd1 <= 1e-6 + @test MOI.get(model, MOI.ConstraintDual(), c2) ≈ -cd1 atol=1e-6 + @test MOI.get(model, MOI.ConstraintDual(), c3) ≈ -cd1 atol=1e-6 end @testset "Default parameters" begin @@ -485,3 +489,125 @@ end names = MOI.get.(dest, MOI.VariableName(), v) @test names == ["a", "b", "c", "d"] end + +@testset "dual_farkas" begin + @testset "test_farkas_dual_min" begin + model = GLPK.Optimizer() + x = MOI.add_variables(model, 2) + MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE) + MOI.set( + model, + MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), + MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [x[1]]), 0.0), + ) + clb = MOI.add_constraint.( + model, MOI.SingleVariable.(x), MOI.GreaterThan(0.0) + ) + c = MOI.add_constraint( + model, + MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 1.0], x), 0.0), + MOI.LessThan(-1.0), + ) + MOI.optimize!(model) + @test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE + @test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE + clb_dual = MOI.get.(model, MOI.ConstraintDual(), clb) + c_dual = MOI.get(model, MOI.ConstraintDual(), c) + @show clb_dual, c_dual + @test clb_dual[1] > -1e-6 + @test clb_dual[2] > -1e-6 + @test c_dual[1] < 1e-6 + @test clb_dual[1] ≈ -2 * c_dual atol = 1e-6 + @test clb_dual[2] ≈ -c_dual atol = 1e-6 + end + + @testset "test_farkas_dual_max" begin + model = GLPK.Optimizer() + x = MOI.add_variables(model, 2) + MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE) + MOI.set( + model, + MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), + MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [x[1]]), 0.0), + ) + clb = MOI.add_constraint.( + model, MOI.SingleVariable.(x), MOI.GreaterThan(0.0) + ) + c = MOI.add_constraint( + model, + MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 1.0], x), 0.0), + MOI.LessThan(-1.0), + ) + MOI.optimize!(model) + @test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE + @test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE + clb_dual = MOI.get.(model, MOI.ConstraintDual(), clb) + c_dual = MOI.get(model, MOI.ConstraintDual(), c) + @show clb_dual, c_dual + @test clb_dual[1] > -1e-6 + @test clb_dual[2] > -1e-6 + @test c_dual[1] < 1e-6 + @test clb_dual[1] ≈ -2 * c_dual atol = 1e-6 + @test clb_dual[2] ≈ -c_dual atol = 1e-6 + end + + @testset "test_farkas_dual_min_ii" begin + model = GLPK.Optimizer() + x = MOI.add_variables(model, 2) + MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE) + MOI.set( + model, + MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), + MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(-1.0, x[1])], 0.0), + ) + clb = MOI.add_constraint.( + model, MOI.SingleVariable.(x), MOI.LessThan(0.0) + ) + c = MOI.add_constraint( + model, + MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-2.0, -1.0], x), 0.0), + MOI.LessThan(-1.0), + ) + MOI.optimize!(model) + @test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE + @test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE + clb_dual = MOI.get.(model, MOI.ConstraintDual(), clb) + c_dual = MOI.get(model, MOI.ConstraintDual(), c) + @show clb_dual, c_dual + @test clb_dual[1] < 1e-6 + @test clb_dual[2] < 1e-6 + @test c_dual[1] < 1e-6 + @test clb_dual[1] ≈ 2 * c_dual atol = 1e-6 + @test clb_dual[2] ≈ c_dual atol = 1e-6 + end + + @testset "test_farkas_dual_max_ii" begin + model = GLPK.Optimizer() + x = MOI.add_variables(model, 2) + MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE) + MOI.set( + model, + MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), + MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(-1.0, x[1])], 0.0), + ) + clb = MOI.add_constraint.( + model, MOI.SingleVariable.(x), MOI.LessThan(0.0) + ) + c = MOI.add_constraint( + model, + MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-2.0, -1.0], x), 0.0), + MOI.LessThan(-1.0), + ) + MOI.optimize!(model) + @test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE + @test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE + clb_dual = MOI.get.(model, MOI.ConstraintDual(), clb) + c_dual = MOI.get(model, MOI.ConstraintDual(), c) + @show clb_dual, c_dual + @test clb_dual[1] < 1e-6 + @test clb_dual[2] < 1e-6 + @test c_dual[1] < 1e-6 + @test clb_dual[1] ≈ 2 * c_dual atol = 1e-6 + @test clb_dual[2] ≈ c_dual atol = 1e-6 + end +end From 5f73763f55786b9165e3e46849710cb62e2a3402 Mon Sep 17 00:00:00 2001 From: odow Date: Tue, 3 Nov 2020 16:47:42 +1300 Subject: [PATCH 2/2] Simplify column passing --- src/MOI_wrapper/MOI_wrapper.jl | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/src/MOI_wrapper/MOI_wrapper.jl b/src/MOI_wrapper/MOI_wrapper.jl index a26f1c7..3a0c8a5 100644 --- a/src/MOI_wrapper/MOI_wrapper.jl +++ b/src/MOI_wrapper/MOI_wrapper.jl @@ -1659,11 +1659,11 @@ function _dual_multiplier(model::Optimizer) return MOI.get(model, MOI.ObjectiveSense()) == MOI.MIN_SENSE ? 1.0 : -1.0 end -function _farkas_variable_dual(model::Optimizer, col::Cint) +function _farkas_variable_dual(model::Optimizer, col::Int) nnz = glp_get_mat_col(model, col, C_NULL, C_NULL) vind = Vector{Cint}(undef, nnz) vval = Vector{Cdouble}(undef, nnz) - nnz = glp_get_mat_col(model, col, offset(vind), offset(vval)) + nnz = glp_get_mat_col(model, Cint(col), offset(vind), offset(vval)) return sum( model.infeasibility_cert[row] * val for (row, val) in zip(vind, vval) ) @@ -1676,11 +1676,11 @@ function MOI.get( ) _throw_if_optimize_in_progress(model, attr) MOI.check_result_index_bounds(model, attr) + col = column(model, c) if model.infeasibility_cert !== nothing - dual = _farkas_variable_dual(model, Cint(column(model, c))) + dual = _farkas_variable_dual(model, col) return min(dual, 0.0) end - col = column(model, c) reduced_cost = if model.method == SIMPLEX || model.method == EXACT glp_get_col_dual(model, col) else @@ -1713,11 +1713,11 @@ function MOI.get( ) _throw_if_optimize_in_progress(model, attr) MOI.check_result_index_bounds(model, attr) + col = column(model, c) if model.infeasibility_cert !== nothing - dual = _farkas_variable_dual(model, Cint(column(model, c))) + dual = _farkas_variable_dual(model, col) return max(dual, 0.0) end - col = column(model, c) reduced_cost = if model.method == SIMPLEX || model.method == EXACT glp_get_col_dual(model, col) else @@ -1750,10 +1750,11 @@ function MOI.get( ) where {S <: Union{MOI.EqualTo, MOI.Interval}} _throw_if_optimize_in_progress(model, attr) MOI.check_result_index_bounds(model, attr) + col = column(model, c) if model.infeasibility_cert !== nothing - return _farkas_variable_dual(model, Cint(column(model, c))) + return _farkas_variable_dual(model, col) end - return _get_col_dual(model, column(model, c)) + return _get_col_dual(model, col) end function MOI.get(