###
##  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
##
##  File :      solutionquality.R
##
##  Purpose :   To demonstrate how to examine the quality of a solution. 
###
library("Rmosek")

solutionquality <- function(filename)
{
    # Read input file 
    r <- mosek_read(filename, list(verbose=1, usesol=FALSE, useparam=TRUE))
    stopifnot(identical(r$response$code, 0))
    prob <- r$prob

    # Ensure that the interior-point solution is available
    prob$iparam <- list(REQUEST_INTPNT=1)

    # Perform the optimization.
    r <- mosek(prob, list(verbose=10, soldetail=2))
    stopifnot(identical(r$response$code, 0))

    # Examine interior-point solution quality.
    isol <- r$sol$itr
    if( is.element(isol$solsta, list("OPTIMAL")) )
    {
      pobj <- isol$pobjval
      dobj <- isol$dobjval

      abs_obj_gap <- abs(dobj-pobj)
      rel_obj_gap <- abs_obj_gap/(1.0 + min(abs(pobj),abs(dobj)))

      max_primal_viol <- max( isol$maxinfeas$pcon, 
                              isol$maxinfeas$pbound,
                              isol$maxinfeas$pbarvar,
                              isol$maxinfeas$pcone )

      max_dual_viol <- max( isol$maxinfeas$dcon,
                            isol$maxinfeas$dbound,
                            isol$maxinfeas$dbarvar,
                            isol$maxinfeas$dcone )

      ## Assume the application needs the solution to be within
      ## 1e-6 of optimality in an absolute sense. Another approach
      ## would be looking at the relative objective gap 

      cat("Customized solution information.\n")
      cat("  Absolute objective gap: ",abs_obj_gap,"\n")
      cat("  Relative objective gap: ",rel_obj_gap,"\n")
      cat("  Max primal violation  : ",max_primal_viol,"\n")
      cat("  Max dual violation    : ",max_dual_viol,"\n")

      accepted <- TRUE

      if( rel_obj_gap > 1e-6 ) {
        print("Warning: The relative objective gap is LARGE.")
        accepted <- FALSE
      }

      ## We will accept a primal and dual infeasibility of 1e-6. 
      ## These numbers should chosen problem dependent.

      if( max_primal_viol > 1e-6 ) {
        print("Warning: Primal violation is too LARGE")
        accepted <- FALSE
      }

      if( max_dual_viol > 1e-6 ) {
        print("Warning: Dual violation is too LARGE.")
        accepted <- FALSE
      }

      if( accepted ) {
        print("Optimal primal solution")
      }

    } else if( is.element(isol$solsta, list("DUAL_INFEASIBLE_CER","PRIMAL_INFEASIBLE_CER")) ) {
        cat("Primal or dual infeasibility certificate found.\n")

    } else {
        cat("Other solution status.\n")
    }
}

# Run only if called directly from the command-line
if( sys.nframe() == 0L )
{
    args <- commandArgs(trailingOnly=TRUE)
    if( !file.exists(args[1]) ) {
        stop("Expected problem file as input argument.")
    }

    solutionquality(args[1])
}