Skip to content

Commit

Permalink
Add Farkas duals for variable bounds (#105)
Browse files Browse the repository at this point in the history
  • Loading branch information
@odow
odow authored Nov 8, 2020
1 parent cbb8000 commit 312e9c0
Show file tree
Hide file tree
Showing 2 changed files with 241 additions and 55 deletions.
100 changes: 50 additions & 50 deletions src/MOI_wrapper/MOI_wrapper.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -492,15 +492,15 @@ function _unsafe_wrap_clp_array(
model::Optimizer,
f::Function,
n::Integer,
indices;
indices = nothing;
own::Bool = false
)
p = f(model)
if p == C_NULL
return map(x -> NaN, indices)
end
x = unsafe_wrap(Array, p, (n,); own = own)
return x[indices]
return indices === nothing ? x : x[indices]
end

function MOI.get(
Expand Down Expand Up @@ -588,80 +588,81 @@ end
# If sense is maximize, we negate all the duals to follow MOI conventions
# Feasibility problems are treated as a minimization

"""
_farkas_variable_dual(model::Optimizer, col::Cint)
Return a Farkas dual associated with the variable bounds of `col`.
Compute the Farkas dual as:
ā * x = λ' * A * x <= λ' * b = -β + sum(āᵢ * Uᵢ | āᵢ < 0) + sum(āᵢ * Lᵢ | āᵢ > 0)
The Farkas dual of the variable is ā, and it applies to the upper bound if ā < 0,
and it applies to the lower bound if ā > 0.
"""
function _farkas_variable_dual(model::Optimizer, col::Integer)
m, n = Clp_getNumRows(model), Clp_getNumCols(model)
nnz = Clp_getNumElements(model)
vbeg = _unsafe_wrap_clp_array(model, Clp_getVectorStarts, n)
vlen = _unsafe_wrap_clp_array(model, Clp_getVectorLengths, n)
indices = vbeg[col] .+ (1:vlen[col])
vind = _unsafe_wrap_clp_array(model, Clp_getIndices, nnz, indices)
vval = _unsafe_wrap_clp_array(model, Clp_getElements, nnz, indices)
# We need to claim ownership of the pointer returned by Clp_infeasibilityRay.
λ = _unsafe_wrap_clp_array(model, Clp_infeasibilityRay, m; own = true)
return sum(v * λ[i + 1] for (i, v) in zip(vind, vval))
end

function MOI.get(
model::Optimizer,
attr::MOI.ConstraintDual,
c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, <:SCALAR_SETS}
c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, <:SCALAR_SETS},
)
MOI.check_result_index_bounds(model, attr)
sense = (Clp_getObjSense(model) == -1) ? -1 : 1
n = Clp_getNumRows(model)
dual_status = MOI.get(model, MOI.DualStatus())
sense = Clp_getObjSense(model)
if dual_status == MOI.FEASIBLE_POINT
dsol = _unsafe_wrap_clp_array(
model,
Clp_getRowPrice,
Clp_getNumRows(model),
c.value,
)
dsol = _unsafe_wrap_clp_array(model, Clp_getRowPrice, n, c.value)
return sense * dsol
elseif dual_status == MOI.INFEASIBILITY_CERTIFICATE
# We claim ownership of the pointer returned by Clp_infeasibilityRay.
dsol = _unsafe_wrap_clp_array(
model,
Clp_infeasibilityRay,
Clp_getNumRows(model),
c.value;
own = true,
return -_unsafe_wrap_clp_array(
model, Clp_infeasibilityRay, n, c.value; own = true,
)
return -sense * dsol
else
error("Dual solution not available")
end
end

# TODO: what happens if problem is unbounded / infeasible?
function MOI.get(
model::Optimizer,
attr::MOI.ConstraintDual,
c::MOI.ConstraintIndex{MOI.SingleVariable, MOI.LessThan{Float64}}
c::MOI.ConstraintIndex{MOI.SingleVariable, MOI.LessThan{Float64}},
)
MOI.check_result_index_bounds(model, attr)
if MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
return min(0.0, _farkas_variable_dual(model, c.value))
end
rc = _unsafe_wrap_clp_array(
model,
Clp_getReducedCost,
Clp_getNumCols(model),
c.value,
model, Clp_getReducedCost, Clp_getNumCols(model), c.value,
)
sense = (Clp_getObjSense(model) == -1) ? -1 : 1
if sense == 1 && rc <= 0.0
return rc
elseif sense == -1 && rc >= 0.0
return -rc
else
return 0.0
end
return min(0.0, Clp_getObjSense(model) * rc)
end

function MOI.get(
model::Optimizer,
attr::MOI.ConstraintDual,
c::MOI.ConstraintIndex{MOI.SingleVariable, MOI.GreaterThan{Float64}}
c::MOI.ConstraintIndex{MOI.SingleVariable, MOI.GreaterThan{Float64}},
)
MOI.check_result_index_bounds(model, attr)
if MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
return max(0.0, _farkas_variable_dual(model, c.value))
end
rc = _unsafe_wrap_clp_array(
model,
Clp_getReducedCost,
Clp_getNumCols(model),
c.value,
model, Clp_getReducedCost, Clp_getNumCols(model), c.value,
)
sense = (Clp_getObjSense(model) == -1) ? -1 : 1
if sense == 1 && rc >= 0.0
return rc
elseif sense == -1 && rc <= 0.0
return -rc
else
return 0.0
end
return max(0.0, Clp_getObjSense(model) * rc)
end

function MOI.get(
Expand All @@ -673,12 +674,11 @@ function MOI.get(
}
)
MOI.check_result_index_bounds(model, attr)
sense = (Clp_getObjSense(model) == -1) ? -1 : 1
if MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
return _farkas_variable_dual(model, c.value)
end
rc = _unsafe_wrap_clp_array(
model,
Clp_getReducedCost,
Clp_getNumCols(model),
c.value,
model, Clp_getReducedCost, Clp_getNumCols(model), c.value,
)
return sense * rc
return Clp_getObjSense(model) * rc
end
196 changes: 191 additions & 5 deletions test/MOI_wrapper.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,12 +57,8 @@ end

function test_contlinear()
MOI.Test.contlineartest(BRIDGED, CONFIG, [
# The linear12 test requires the InfeasibilityCertificate for variable
# bounds. These are available through C++, but not via the C interface.
"linear12",

# MOI.VariablePrimalStart not supported.
"partial_start"
"partial_start",
])
end

Expand Down Expand Up @@ -155,6 +151,196 @@ function test_options_after_empty!()
@test MOI.get(model, MOI.Silent()) == true
end

function test_farkas_dual_min()
MOI.empty!(BRIDGED)
model = BRIDGED
MOI.set(model, MOI.Silent(), true)
x = MOI.add_variables(model, 2)
MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
MOI.set(
model,
MOI.ObjectiveFunction{MOI.SingleVariable}(),
MOI.SingleVariable(x[1]),
)
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

function test_farkas_dual_min_interval()
MOI.empty!(BRIDGED)
model = BRIDGED
MOI.set(model, MOI.Silent(), true)
x = MOI.add_variables(model, 2)
MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
MOI.set(
model,
MOI.ObjectiveFunction{MOI.SingleVariable}(),
MOI.SingleVariable(x[1]),
)
clb = MOI.add_constraint.(
model, MOI.SingleVariable.(x), MOI.Interval(0.0, 10.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

function test_farkas_dual_min_equalto()
MOI.empty!(BRIDGED)
model = BRIDGED
MOI.set(model, MOI.Silent(), true)
x = MOI.add_variables(model, 2)
MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
MOI.set(
model,
MOI.ObjectiveFunction{MOI.SingleVariable}(),
MOI.SingleVariable(x[1]),
)
clb = MOI.add_constraint.(model, MOI.SingleVariable.(x), MOI.EqualTo(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

function test_farkas_dual_min_ii()
MOI.empty!(BRIDGED)
model = BRIDGED
MOI.set(model, MOI.Silent(), true)
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

function test_farkas_dual_max()
MOI.empty!(BRIDGED)
model = BRIDGED
MOI.set(model, MOI.Silent(), true)
x = MOI.add_variables(model, 2)
MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
MOI.set(
model,
MOI.ObjectiveFunction{MOI.SingleVariable}(),
MOI.SingleVariable(x[1]),
)
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

function test_farkas_dual_max_ii()
MOI.empty!(BRIDGED)
model = BRIDGED
MOI.set(model, MOI.Silent(), true)
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 # module TestMOIWrapper

runtests(TestMOIWrapper)

0 comments on commit 312e9c0

Please sign in to comment.