2017-03-13 14:49:55 +01:00
|
|
|
using JuMP
|
|
|
|
|
2018-08-20 03:50:03 +02:00
|
|
|
function constraints(pm::Matrix{I}, total_length=maximum(pm)) where {I<:Integer}
|
|
|
|
cnstrs = [Vector{I}() for _ in 1:total_length]
|
|
|
|
for i in eachindex(pm)
|
|
|
|
push!(cnstrs[pm[i]], i)
|
|
|
|
end
|
|
|
|
return cnstrs
|
|
|
|
end
|
|
|
|
|
2018-08-20 03:50:16 +02:00
|
|
|
function constraint(pm::Matrix{I}, k) where {I<:Integer}
|
|
|
|
cnstr = Vector{I}()
|
|
|
|
for i in eachindex(pm)
|
|
|
|
if pm[i] == k
|
|
|
|
push!(cnstr, i)
|
2017-03-13 14:49:55 +01:00
|
|
|
end
|
|
|
|
end
|
2018-08-20 03:50:16 +02:00
|
|
|
return cnstr
|
2017-03-13 14:49:55 +01:00
|
|
|
end
|
|
|
|
|
|
|
|
|
2018-08-20 03:46:44 +02:00
|
|
|
function SOS_problem(X::GroupRingElem, orderunit::GroupRingElem; upper_bound=Inf)
|
|
|
|
N = size(parent(X).pm, 1)
|
|
|
|
matrix_constraints = PropertyT.constraints(parent(X).pm)
|
2017-03-13 14:49:55 +01:00
|
|
|
m = JuMP.Model();
|
2018-08-20 03:46:44 +02:00
|
|
|
|
2017-06-05 13:55:24 +02:00
|
|
|
JuMP.@variable(m, P[1:N, 1:N])
|
2017-04-01 15:21:57 +02:00
|
|
|
JuMP.@SDconstraint(m, P >= 0)
|
|
|
|
JuMP.@constraint(m, sum(P[i] for i in eachindex(P)) == 0)
|
2017-03-20 21:41:12 +01:00
|
|
|
|
2018-08-20 03:46:44 +02:00
|
|
|
JuMP.@variable(m, λ)
|
2017-03-13 14:49:55 +01:00
|
|
|
if upper_bound < Inf
|
2018-08-20 03:46:44 +02:00
|
|
|
JuMP.@constraint(m, λ <= upper_bound)
|
2017-03-13 14:49:55 +01:00
|
|
|
end
|
|
|
|
|
2018-08-20 03:46:44 +02:00
|
|
|
for (cnstr, x, u) in zip(matrix_constraints, X.coeffs, orderunit.coeffs)
|
|
|
|
JuMP.@constraint(m, sum(P[cnstr]) == x - λ*u)
|
2017-03-13 14:49:55 +01:00
|
|
|
end
|
2017-03-20 21:41:12 +01:00
|
|
|
|
2017-04-01 15:21:57 +02:00
|
|
|
JuMP.@objective(m, Max, λ)
|
2017-06-04 20:13:27 +02:00
|
|
|
return m, λ, P
|
2017-03-13 14:49:55 +01:00
|
|
|
end
|
|
|
|
|
2018-08-20 03:45:50 +02:00
|
|
|
|
|
|
|
end
|
|
|
|
|
|
|
|
|
|
|
|
end
|
|
|
|
|
2018-09-05 09:14:50 +02:00
|
|
|
###############################################################################
|
|
|
|
#
|
|
|
|
# Low-level solve
|
|
|
|
#
|
|
|
|
###############################################################################
|
|
|
|
using MathProgBase
|
|
|
|
|
|
|
|
function solve(m::JuMP.Model, varλ::JuMP.Variable, varP; warmstart=nothing)
|
2017-03-13 14:49:55 +01:00
|
|
|
|
2018-01-02 02:52:45 +01:00
|
|
|
traits = JuMP.ProblemTraits(m, relaxation=false)
|
2017-03-16 09:35:32 +01:00
|
|
|
|
2018-01-02 02:52:45 +01:00
|
|
|
JuMP.build(m, traits=traits)
|
|
|
|
if warmstart != nothing
|
|
|
|
p_sol, d_sol, s = warmstart
|
2018-09-05 09:14:50 +02:00
|
|
|
MathProgBase.SolverInterface.setwarmstart!(m.internalModel, p_sol;
|
|
|
|
dual_sol=d_sol, slack=s);
|
2018-01-02 02:52:45 +01:00
|
|
|
end
|
2017-03-16 09:35:32 +01:00
|
|
|
|
2018-01-02 02:52:45 +01:00
|
|
|
MathProgBase.optimize!(m.internalModel)
|
2017-03-13 14:49:55 +01:00
|
|
|
|
2018-01-02 02:52:45 +01:00
|
|
|
λ = MathProgBase.getobjval(m.internalModel)
|
|
|
|
|
|
|
|
warmstart = (m.internalModel.primal_sol, m.internalModel.dual_sol,
|
|
|
|
m.internalModel.slack)
|
|
|
|
|
|
|
|
fillfrominternal!(m, traits)
|
|
|
|
|
|
|
|
P = JuMP.getvalue(varP)
|
|
|
|
λ = JuMP.getvalue(varλ)
|
2017-03-13 14:49:55 +01:00
|
|
|
|
2018-01-02 02:52:45 +01:00
|
|
|
return λ, P, warmstart
|
2017-03-13 14:49:55 +01:00
|
|
|
end
|
2018-01-01 23:59:31 +01:00
|
|
|
|
2018-09-05 09:14:50 +02:00
|
|
|
function solve_logged(model::JuMP.Model, varλ::JuMP.Variable, varP, warmstart=nothing; solverlog::String=tempname()*".log")
|
|
|
|
|
|
|
|
function f()
|
|
|
|
Base.Libc.flush_cstdio()
|
|
|
|
λ, P, ws = PropertyT.solve(model, varλ, varP, warmstart=warmstart)
|
|
|
|
Base.Libc.flush_cstdio()
|
|
|
|
return λ, P, ws
|
|
|
|
end
|
|
|
|
|
|
|
|
log = open(solverlog, "a+")
|
|
|
|
λ, P, warmstart = redirect_stdout(f, log)
|
|
|
|
close(log)
|
|
|
|
|
|
|
|
return λ, P, warmstart
|
|
|
|
end
|
2018-01-01 23:59:31 +01:00
|
|
|
function fillfrominternal!(m::JuMP.Model, traits)
|
|
|
|
# Copied from JuMP/src/solvers.jl:178
|
|
|
|
|
|
|
|
stat::Symbol = MathProgBase.status(m.internalModel)
|
|
|
|
|
|
|
|
numRows, numCols = length(m.linconstr), m.numCols
|
|
|
|
m.objBound = NaN
|
|
|
|
m.objVal = NaN
|
|
|
|
m.colVal = fill(NaN, numCols)
|
|
|
|
m.linconstrDuals = Array{Float64}(0)
|
|
|
|
|
|
|
|
discrete = (traits.int || traits.sos)
|
|
|
|
|
|
|
|
if stat == :Optimal
|
|
|
|
# If we think dual information might be available, try to get it
|
|
|
|
# If not, return an array of the correct length
|
|
|
|
if discrete
|
|
|
|
m.redCosts = fill(NaN, numCols)
|
|
|
|
m.linconstrDuals = fill(NaN, numRows)
|
|
|
|
else
|
|
|
|
if !traits.conic
|
|
|
|
m.redCosts = try
|
|
|
|
MathProgBase.getreducedcosts(m.internalModel)[1:numCols]
|
|
|
|
catch
|
|
|
|
fill(NaN, numCols)
|
|
|
|
end
|
|
|
|
|
|
|
|
m.linconstrDuals = try
|
|
|
|
MathProgBase.getconstrduals(m.internalModel)[1:numRows]
|
|
|
|
catch
|
|
|
|
fill(NaN, numRows)
|
|
|
|
end
|
|
|
|
elseif !traits.qp && !traits.qc
|
|
|
|
JuMP.fillConicDuals(m)
|
|
|
|
end
|
|
|
|
end
|
|
|
|
else
|
|
|
|
# Problem was not solved to optimality, attempt to extract useful
|
|
|
|
# information anyway
|
|
|
|
|
|
|
|
if traits.lin
|
|
|
|
if stat == :Infeasible
|
|
|
|
m.linconstrDuals = try
|
|
|
|
infray = MathProgBase.getinfeasibilityray(m.internalModel)
|
|
|
|
@assert length(infray) == numRows
|
|
|
|
infray
|
|
|
|
catch
|
|
|
|
suppress_warnings || warn("Infeasibility ray (Farkas proof) not available")
|
|
|
|
fill(NaN, numRows)
|
|
|
|
end
|
|
|
|
elseif stat == :Unbounded
|
|
|
|
m.colVal = try
|
|
|
|
unbdray = MathProgBase.getunboundedray(m.internalModel)
|
|
|
|
@assert length(unbdray) == numCols
|
|
|
|
unbdray
|
|
|
|
catch
|
|
|
|
suppress_warnings || warn("Unbounded ray not available")
|
|
|
|
fill(NaN, numCols)
|
|
|
|
end
|
|
|
|
end
|
|
|
|
end
|
|
|
|
# conic duals (currently, SOC and SDP only)
|
|
|
|
if !discrete && traits.conic && !traits.qp && !traits.qc
|
|
|
|
if stat == :Infeasible
|
|
|
|
JuMP.fillConicDuals(m)
|
|
|
|
end
|
|
|
|
end
|
|
|
|
end
|
|
|
|
|
|
|
|
# If the problem was solved, or if it terminated prematurely, try
|
|
|
|
# to extract a solution anyway. This commonly occurs when a time
|
|
|
|
# limit or tolerance is set (:UserLimit)
|
|
|
|
if !(stat == :Infeasible || stat == :Unbounded)
|
|
|
|
try
|
|
|
|
# Do a separate try since getobjval could work while getobjbound does not and vice versa
|
|
|
|
objBound = MathProgBase.getobjbound(m.internalModel) + m.obj.aff.constant
|
|
|
|
m.objBound = objBound
|
|
|
|
end
|
|
|
|
try
|
|
|
|
objVal = MathProgBase.getobjval(m.internalModel) + m.obj.aff.constant
|
|
|
|
colVal = MathProgBase.getsolution(m.internalModel)[1:numCols]
|
|
|
|
# Rescale off-diagonal terms of SDP variables
|
|
|
|
if traits.sdp
|
|
|
|
offdiagvars = JuMP.offdiagsdpvars(m)
|
|
|
|
colVal[offdiagvars] /= sqrt(2)
|
|
|
|
end
|
|
|
|
# Don't corrupt the answers if one of the above two calls fails
|
|
|
|
m.objVal = objVal
|
|
|
|
m.colVal = colVal
|
|
|
|
end
|
|
|
|
end
|
|
|
|
|
|
|
|
return stat
|
|
|
|
end
|