Skip to content

Add infeasibility certificate for variable bounds #155

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 2 commits into from
Nov 8, 2020
Merged

Add infeasibility certificate for variable bounds #155

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
ec41a47
Add farkas duals for variable bounds and refactor how we compute
odow Nov 3, 2020
5f73763
Simplify column passing
odow Nov 3, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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"
Expand Down
123 changes: 78 additions & 45 deletions src/MOI_wrapper/MOI_wrapper.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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}}
Expand Down Expand Up @@ -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

Expand All @@ -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

Expand All @@ -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

Expand All @@ -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

Expand All @@ -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

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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

Expand All @@ -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::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, Cint(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
col = column(model, c)
if model.infeasibility_cert !== nothing
dual = _farkas_variable_dual(model, col)
return min(dual, 0.0)
end
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
Expand All @@ -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
col = column(model, c)
if model.infeasibility_cert !== nothing
dual = _farkas_variable_dual(model, col)
return max(dual, 0.0)
end
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
Expand All @@ -1715,23 +1744,29 @@ 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)
col = column(model, c)
if model.infeasibility_cert !== nothing
return _farkas_variable_dual(model, col)
end
return _get_col_dual(model, col)
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
Expand Down Expand Up @@ -1911,7 +1946,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))
Expand All @@ -1933,9 +1968,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

Expand Down Expand Up @@ -1979,8 +2012,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
Expand Down
Loading