From e94662e615e305a9ac2ef3e9d930b51bd08c0916 Mon Sep 17 00:00:00 2001
From: kalmarek <kalmar@amu.edu.pl>
Date: Fri, 1 Dec 2017 17:06:20 +0100
Subject: [PATCH] fill the JuMP model variables from its internal model

---
 src/PropertyT.jl | 98 ++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 98 insertions(+)

diff --git a/src/PropertyT.jl b/src/PropertyT.jl
index a8ba76d..813283a 100644
--- a/src/PropertyT.jl
+++ b/src/PropertyT.jl
@@ -161,6 +161,101 @@ function λandP(name::String, SDP_problem::JuMP.Model, varλ, varP)
 
 end
 
+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
+
 function compute_λandP(m, varλ, varP; warmstart=nothing)
     λ = 0.0
     P = nothing
@@ -183,6 +278,9 @@ function compute_λandP(m, varλ, varP; warmstart=nothing)
 
     warmstart = (m.internalModel.primal_sol, m.internalModel.dual_sol,
           m.internalModel.slack)
+
+    fillfrominternal!(m, traits)
+
     P = JuMP.getvalue(varP)
     λ = JuMP.getvalue(varλ)