##
#   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#   File:      pinfeas.py
#
#   Purpose: Demonstrates how to fetch a primal infeasibility certificate
#            for a linear problem
#
##

from mosek.fusion import *
import sys

# Analyzes and prints infeasibility certificate for a single object,
# which can be a variable or constraint
def analyzeCertificate(name,          # name of the analyzed object
                       size,          # size of the object
                       duals,         # double[], actual dual values
                       eps):          # tolerance determining when a dual value is considered important
    for i in range(0, size):
        if abs(duals[i]) > eps:
            print(f"{name}[{i}], dual = {duals[i]}")

# Construct the sample model from the example in the manual
sMat = Matrix.sparse(3, 7, [0,0,1,1,2,2,2],
                           [0,1,2,3,4,5,6],
                           [1.0]*7)
sBound = [200, 1000, 1000]
dMat = Matrix.sparse(4, 7, [0,0,1,2,2,3,3],
                           [0,4,1,2,5,3,6],
                           [1.0]*7)
dBound = [1100, 200, 500, 500]
c = [1, 2, 5, 2, 1, 2, 1]

M = Model("pinfeas")
x = M.variable("x", 7, Domain.greaterThan(0))
s = M.constraint("s", Expr.mul(sMat, x), Domain.lessThan(sBound))
d = M.constraint("d", Expr.mul(dMat, x), Domain.equalsTo(dBound))
M.objective(ObjectiveSense.Minimize, Expr.dot(c,x))

# Set useful debugging tools and solve the model
M.setLogHandler(sys.stdout)
M.writeTask("pinfeas.ptf")
M.solve()

# Check problem status
if M.getProblemStatus() == ProblemStatus.PrimalInfeasible:
    # Set the tolerance at which we consider a dual value as essential
    eps = 1e-7

    # We want to retrieve infeasibility certificates
    M.acceptedSolutionStatus(AccSolutionStatus.Certificate)

    # Go through variable bounds
    print("Variable bounds important for infeasibility: ")
    analyzeCertificate(name = "x", size = x.getSize(), duals = x.dual(), eps = eps)

    # Go through constraint bounds
    print("Constraint bounds important for infeasibility: ")
    analyzeCertificate(name = "s", size = s.getSize(), duals = s.dual(), eps = eps)
    analyzeCertificate(name = "d", size = d.getSize(), duals = d.dual(), eps = eps)

else:
    print("The problem is not primal infeasible, no certificate to show")