R/250_reductions_solvers_qp_solvers_cplex_qpif.R

In CVXR: Disciplined Convex Optimization

#####
## DO NOT EDIT THIS FILE!! EDIT THE SOURCE INSTEAD: rsrc_tree/reductions/solvers/qp_solvers/cplex_qpif.R
#####

## CVXPY SOURCE: reductions/solvers/qp_solvers/cplex_qpif.py
## CPLEX QP solver interface for LP/QP/MIP problems via Rcplex
##
## Uses QpSolver.apply() for sign-correct A_eq/b_eq/F_ineq/g_ineq data.
## MIP_CAPABLE: supports boolean and integer variables.
##
## Key conventions:
##   - Rcplex uses 0.5 x'Qx + c'x -- matches CVXR's QP path directly (no halving)
##   - Rcplex status: 1=optimal, 2=unbounded, 3=infeasible, 4=inf_or_unbd, 101=MIP_optimal, etc.
##   - Duals in extra$lambda (absent for MIP). CVXPY negates all: y = -lambda
##   - vtype = NULL for continuous (avoids CPLEX MIP mode); set to "B"/"I" vector only for MIP


# -- CPLEX status map ----------------------------------------------------------
## Maps Rcplex integer status codes to CVXR status strings.
## LP codes: 1-6, MIP codes: 101-108+
## CVXPY SOURCE: cplex_conif.py get_status() / cplex_qpif.py invert()

CPLEX_STATUS_MAP <- list(
  "1"   = OPTIMAL,                    # optimal
  "2"   = UNBOUNDED,                  # unbounded
  "3"   = INFEASIBLE,                 # infeasible
  "4"   = INFEASIBLE_OR_UNBOUNDED,    # infeasible_or_unbounded
  "5"   = OPTIMAL_INACCURATE,         # optimal_with_unscaled_infeasibilities
  "6"   = OPTIMAL_INACCURATE,         # feasible (non-optimal)
  "10"  = SOLVER_ERROR,               # abort_user
  "11"  = SOLVER_ERROR,               # abort_iteration_limit
  "12"  = USER_LIMIT,                 # abort_time_limit
  "13"  = USER_LIMIT,                 # abort_dettime_limit
  "101" = OPTIMAL,                    # MIP_optimal
  "102" = OPTIMAL_INACCURATE,         # MIP_optimal_tolerance
  "103" = INFEASIBLE,                 # MIP_infeasible
  "104" = USER_LIMIT,                 # MIP_sol_limit
  "107" = USER_LIMIT,                 # MIP_time_limit_feasible
  "108" = INFEASIBLE_INACCURATE,      # MIP_time_limit_infeasible
  "115" = INFEASIBLE_OR_UNBOUNDED,    # MIP_infeasible_or_unbounded
  "118" = UNBOUNDED                   # MIP_unbounded
)

# -- CPLEX_QP_Solver class -----------------------------------------------------
## CVXPY SOURCE: cplex_qpif.py lines 31-34

CPLEX_QP_Solver <- new_class("CPLEX_QP_Solver", parent = QpSolver,
  package = "CVXR",
  constructor = function() {
    new_object(S7_object(),
      .cache = new.env(parent = emptyenv()),
      MIP_CAPABLE = TRUE,
      BOUNDED_VARIABLES = FALSE,
      SUPPORTED_CONSTRAINTS = list(Zero, NonNeg),
      REQUIRES_CONSTR = FALSE
    )
  }
)

method(solver_name, CPLEX_QP_Solver) <- function(x) CPLEX_SOLVER

# -- solve_via_data ------------------------------------------------------------
## CVXPY SOURCE: cplex_qpif.py lines 91-189
## Receives QP data from QpSolver.apply(): P, q, A_eq, b_eq, F_ineq, g_ineq
## Rcplex convention: min 0.5 x'Qx + c'x s.t. Ax {sense} b, lb <= x <= ub

method(solve_via_data, CPLEX_QP_Solver) <- function(x, data, warm_start = FALSE, verbose = FALSE,
                                                      solver_opts = list(), ...) {
  if (!requireNamespace("Rcplex", quietly = TRUE)) {
    cli_abort("Package {.pkg Rcplex} is required but not installed.")
  }

  q_vec <- data[["q"]]
  nvars <- length(q_vec)

  ## Build constraint matrix: stack A_eq and F_ineq
  A_eq <- data[["A_eq"]]
  b_eq <- data[["b_eq"]]
  F_ineq <- data[["F_ineq"]]
  g_ineq <- data[["g_ineq"]]

  len_eq <- nrow(A_eq)
  len_ineq <- nrow(F_ineq)
  n_constrs <- len_eq + len_ineq

  if (n_constrs > 0L) {
    if (len_eq > 0L && len_ineq > 0L) {
      Amat <- rbind(A_eq, F_ineq)
      bvec <- c(b_eq, g_ineq)
      sense <- c(rep("E", len_eq), rep("L", len_ineq))
    } else if (len_eq > 0L) {
      Amat <- A_eq
      bvec <- b_eq
      sense <- rep("E", len_eq)
    } else {
      Amat <- F_ineq
      bvec <- g_ineq
      sense <- rep("L", len_ineq)
    }
  } else {
    ## No constraints: add a dummy row (Rcplex requires at least one constraint)
    Amat <- Matrix::sparseMatrix(i = integer(0), j = integer(0),
                                  x = numeric(0), dims = c(1L, nvars))
    bvec <- 0
    sense <- "L"
  }

  ## Ensure dgCMatrix for Rcplex
  if (!inherits(Amat, "dgCMatrix")) {
    Amat <- methods::as(Amat, "dgCMatrix")
  }

  ## Quadratic objective matrix
  ## Rcplex uses 0.5 x'Qx + c'x -- matches CVXR's P convention directly
  Qmat <- data[[SD_P]]
  if (!is.null(Qmat)) {
    if (!inherits(Qmat, "dgCMatrix")) {
      Qmat <- methods::as(methods::as(Qmat, "generalMatrix"), "CsparseMatrix")
    }
  }

  ## Variable types for MIP
  ## CRITICAL: vtype = NULL for continuous problems (avoids CPLEX MIP mode)
  ## Only set vtype when bool_idx or int_idx present
  ## CVXPY SOURCE: cplex_qpif.py lines 118-124
  bool_idx <- data[["bool_idx"]]
  int_idx <- data[["int_idx"]]
  is_mip <- length(bool_idx %||% integer(0)) > 0L ||
    length(int_idx %||% integer(0)) > 0L

  vtype <- NULL
  if (is_mip) {
    vtype <- rep("C", nvars)
    if (length(bool_idx) > 0L) {
      for (idx in bool_idx) {
        vtype[idx] <- "B"
      }
    }
    if (length(int_idx) > 0L) {
      for (idx in int_idx) {
        vtype[idx] <- "I"
      }
    }
  }

  ## Control parameters
  control <- list(trace = if (verbose) 1L else 0L)

  ## Apply user solver options
  for (opt_name in names(solver_opts)) {
    control[[opt_name]] <- solver_opts[[opt_name]]
  }

  ## Call Rcplex
  ## Rcplex prints "CPLEX environment opened" to stderr on first call
  ## (C-level init message). Suppress when not verbose via capture.output.
  .cplex_call <- function() {
    Rcplex::Rcplex(
      cvec = q_vec,
      Amat = Amat,
      bvec = bvec,
      Qmat = Qmat,
      lb = rep(-Inf, nvars),
      ub = rep(Inf, nvars),
      control = control,
      objsense = "min",
      sense = sense,
      vtype = vtype
    )
  }
  result <- tryCatch({
    if (!verbose) {
      res <- NULL
      capture.output(res <- .cplex_call(), type = "message")
      res
    } else {
      .cplex_call()
    }
  },
    error = function(e) {
      list(xopt = NA, obj = NA, status = -1L,
           extra = list(lambda = NA, slack = NA))
    }
  )

  ## Store metadata for reduction_invert
  result$.len_eq <- len_eq
  result$.n_constrs <- n_constrs
  result$.is_mip <- is_mip

  result
}

# -- reduction_invert ----------------------------------------------------------
## CVXPY SOURCE: cplex_qpif.py lines 43-89
## Dual sign: negate ALL lambda duals -- matching CVXPY cplex_qpif.py line 72.

method(reduction_invert, CPLEX_QP_Solver) <- function(x, solution, inverse_data, ...) {
  attr_list <- list()

  ## Map CPLEX integer status to CVXR status
  status <- CPLEX_STATUS_MAP[[as.character(solution$status)]]
  if (is.null(status)) status <- SOLVER_ERROR

  if (status %in% SOLUTION_PRESENT) {
    opt_val <- solution$obj + inverse_data[[SD_OFFSET]]

    ## Primal variables
    primal_vars <- list()
    primal_vars[[as.character(inverse_data[[SOLVER_VAR_ID]])]] <- solution$xopt

    ## Dual variables: negate ALL lambda duals
    ## CVXPY: y = -np.array(model.solution.get_dual_values())
    is_mip <- isTRUE(solution$.is_mip)

    if (!is_mip && !is.null(solution$extra$lambda) &&
        !any(is.na(solution$extra$lambda))) {
      y <- -solution$extra$lambda
      len_eq <- solution$.len_eq

      eq_dual <- if (len_eq > 0L) {
        get_dual_values(
          y[seq_len(len_eq)],
          extract_dual_value,
          inverse_data[[SOLVER_EQ_CONSTR]]
        )
      } else {
        list()
      }

      ineq_dual <- if (len_eq < length(y)) {
        get_dual_values(
          y[(len_eq + 1L):length(y)],
          extract_dual_value,
          inverse_data[[SOLVER_NEQ_CONSTR]]
        )
      } else {
        list()
      }

      dual_vars <- c(eq_dual, ineq_dual)
    } else {
      dual_vars <- list()
    }

    Solution(status, opt_val, primal_vars, dual_vars, attr_list)
  } else {
    failure_solution(status, attr_list)
  }
}

# -- print ---------------------------------------------------------------------

method(print, CPLEX_QP_Solver) <- function(x, ...) {
  cat("CPLEX_QP_Solver()\n")
  invisible(x)
}