R/248_reductions_solvers_conic_solvers_xpress_conif.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: reductions/solvers/conic_solvers/xpress_conif.py
## XPRESS conic solver interface for LP, SOCP, MI-LP, MI-SOCP problems
##
## Uses ConicSolver pipeline (ConeMatrixStuffing -> format_constraints).
## MIP_CAPABLE = TRUE, SUPPORTED_CONSTRAINTS = {Zero, NonNeg, SOC}.
##
## Key differences from Python xpress_conif.py:
##   - R xpress uses xprs_loadproblemdata() + addrows() + chgmcoef() + addqmatrix()
##   - xprs_newcol() returns 0-based column index (needed for addqmatrix)
##   - Controls set via getcontrolinfo() type dispatch
##   - Duals scoped to first nrows_linear rows via explicit getduals(first, last)
##   - INPUTROWS updates after addrows() — cannot rely on it for dual scoping


## -- XPRESS status map --------------------------------------------------------
## Empirically verified (2026-03-26, xpress 46.01 / Xpress 9.8):
##   xpress:::SOLSTATUS_NOTFOUND   = 0
##   xpress:::SOLSTATUS_OPTIMAL    = 1
##   xpress:::SOLSTATUS_FEASIBLE   = 2
##   xpress:::SOLSTATUS_INFEASIBLE = 3  (NOT unbounded!)
##   xpress:::SOLSTATUS_UNBOUNDED  = 4  (NOT infeasible!)

XPRESS_STATUS_MAP <- list(
  "0" = SOLVER_ERROR,
  "1" = OPTIMAL,
  "2" = OPTIMAL_INACCURATE,
  "3" = INFEASIBLE,
  "4" = UNBOUNDED
)

## -- Helper: set XPRESS controls from named list ------------------------------
## Maps string control names to integer IDs via getcontrolinfo(), then
## dispatches to setintcontrol() (type=1) or setdblcontrol() (type=3).
.xpress_set_controls <- function(prob, opts) {
  for (nm in names(opts)) {
    info <- tryCatch(
      xpress::getcontrolinfo(prob, nm),
      error = function(e) list(id = 0L, type = 0L)
    )
    if (info$id == 0L) next
    if (info$type == 1L) {
      xpress::setintcontrol(prob, info$id, as.integer(opts[[nm]]))
    } else if (info$type == 3L) {
      xpress::setdblcontrol(prob, info$id, as.numeric(opts[[nm]]))
    }
  }
}


# -- XPRESS_Conic_Solver class ------------------------------------------------
## CVXPY SOURCE: xpress_conif.py class XPRESS
## Supports Zero + NonNeg + SOC constraints.
## SOC encoded via auxiliary variables + addqmatrix (quadratic row constraints).

#' @keywords internal
XPRESS_Conic_Solver <- new_class("XPRESS_Conic_Solver", parent = ConicSolver,
  package = "CVXR",
  constructor = function() {
    new_object(S7_object(),
      .cache = new.env(parent = emptyenv()),
      MIP_CAPABLE = TRUE,
      BOUNDED_VARIABLES = FALSE,
      SUPPORTED_CONSTRAINTS = list(Zero, NonNeg, SOC),
      EXP_CONE_ORDER = NULL,
      REQUIRES_CONSTR = FALSE
    )
  }
)

method(solver_name, XPRESS_Conic_Solver) <- function(x) XPRESS_SOLVER


# -- reduction_accepts ---------------------------------------------------------
## Override: XPRESS conic handles Zero + NonNeg + SOC only.

method(reduction_accepts, XPRESS_Conic_Solver) <- function(x, problem, ...) {
  if (!is.list(problem) || is.null(problem[["constraints"]])) return(FALSE)
  constrs <- problem[["constraints"]]
  all(vapply(constrs, function(c) {
    S7_inherits(c, Zero) || S7_inherits(c, NonNeg) || S7_inherits(c, SOC)
  }, logical(1L)))
}

# -- reduction_apply -----------------------------------------------------------
## Override to add MIP index pass-through (same pattern as Gurobi).

method(reduction_apply, XPRESS_Conic_Solver) <- function(x, problem, ...) {
  data <- problem
  inv_data <- list()

  inv_data[[SOLVER_VAR_ID]] <- data[["x_id"]]

  constraints <- data[["constraints"]]
  cone_dims <- data[[SD_DIMS]]
  inv_data[[SD_DIMS]] <- cone_dims

  constr_map <- group_constraints(constraints)
  inv_data[[SOLVER_EQ_CONSTR]] <- constr_map[["Zero"]]
  inv_data[[SOLVER_NEQ_CONSTR]] <- c(
    constr_map[["NonNeg"]], constr_map[["SOC"]]
  )

  formatted <- format_constraints(
    constraints, data[[SD_A]], data[[SD_B]],
    exp_cone_order = x@EXP_CONE_ORDER
  )

  solver_data <- list()
  solver_data[[SD_C]] <- data[[SD_C]]
  solver_data[[SD_A]] <- -formatted$A
  solver_data[[SD_B]] <- formatted$b
  solver_data[[SD_DIMS]] <- cone_dims

  if (!is.null(data[["bool_idx"]])) solver_data[["bool_idx"]] <- data[["bool_idx"]]
  if (!is.null(data[["int_idx"]]))  solver_data[["int_idx"]]  <- data[["int_idx"]]

  inv_data[[SD_OFFSET]] <- data[[SD_OFFSET]]
  inv_data[["is_mip"]] <- (length(data[["bool_idx"]] %||% integer(0)) > 0L ||
                            length(data[["int_idx"]]  %||% integer(0)) > 0L)

  list(solver_data, inv_data)
}

# -- solve_via_data ------------------------------------------------------------
## CVXPY SOURCE: xpress_conif.py solve_via_data
## Builds XPRESS problem: linear constraints via xprs_loadproblemdata,
## SOC via auxiliary variables + addqmatrix.

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

  c_vec <- data[[SD_C]]
  A     <- data[[SD_A]]
  b     <- data[[SD_B]]
  dims  <- data[[SD_DIMS]]
  nvars <- length(c_vec)

  ## -- Dimension splitting ------------------------------------------------
  nrowsEQ  <- dims@zero
  nrowsLEQ <- dims@nonneg
  nrows_linear <- nrowsEQ + nrowsLEQ
  soc_dims <- dims@soc

  ## -- MIP detection ------------------------------------------------------
  bool_idx <- data[["bool_idx"]] %||% integer(0)
  int_idx  <- data[["int_idx"]]  %||% integer(0)
  is_mip   <- length(bool_idx) > 0L || length(int_idx) > 0L

  ## -- Variable types and bounds ------------------------------------------
  columntypes <- rep("C", nvars)
  lb <- rep(-Inf, nvars)
  ub <- rep(Inf, nvars)

  if (length(bool_idx) > 0L) {
    columntypes[bool_idx] <- "B"
    lb[bool_idx] <- pmax(lb[bool_idx], 0)
    ub[bool_idx] <- pmin(ub[bool_idx], 1)
  }
  if (length(int_idx) > 0L) {
    columntypes[int_idx] <- "I"
  }

  ## -- Build problem data for linear part ---------------------------------
  problemdata <- list(probname = "CVX_xpress_conic", objcoef = c_vec)

  if (nrows_linear > 0L) {
    A_lin <- A[seq_len(nrows_linear), , drop = FALSE]
    b_lin <- b[seq_len(nrows_linear)]
    problemdata$A <- A_lin
    problemdata$rhs <- b_lin
    problemdata$rowtype <- c(rep("E", nrowsEQ), rep("L", nrowsLEQ))
  } else {
    problemdata$A <- Matrix::sparseMatrix(i = integer(0), j = integer(0),
                                           x = numeric(0), dims = c(0L, nvars))
    problemdata$rhs <- numeric(0)
    problemdata$rowtype <- character(0)
  }

  problemdata$lb <- lb
  problemdata$ub <- ub
  if (is_mip) problemdata$columntypes <- columntypes

  ## -- Create problem and load --------------------------------------------
  prob <- xpress::createprob()
  on.exit(xpress::destroyprob(prob), add = TRUE)

  xpress::xprs_loadproblemdata(prob, problemdata = problemdata)

  ## -- SOC constraints ----------------------------------------------------
  ## For each SOC of dimension k:
  ##   (a) Add k auxiliary cone variables via xprs_newcol (0-based indices)
  ##   (b) Add k transformation rows (equality): A_soc * x_orig - cone_var = -b_soc
  ##       Equivalently: add rows with A_soc coefficients, then chgmcoef to add -1
  ##       for cone vars. But CVXPY uses +1 for cone vars (b - A*x = cone_var).
  ##       Since solver_data$A = -formatted_A, the transformation is:
  ##       A_soc_row * x_orig + 1*cone_var = b_soc  (cone_var = b - A*x)
  ##   (c) For k > 1: add quadratic constraint v0^2 - v1^2 - ... - v_{k-1}^2 >= 0

  if (length(soc_dims) > 0L) {
    currow <- nrows_linear  ## 0-based row offset into full A/b

    for (i_cone in seq_along(soc_dims)) {
      k <- soc_dims[i_cone]
      soc_rows <- (currow + 1L):(currow + k)
      A_soc <- A[soc_rows, , drop = FALSE]  ## k x nvars
      b_soc <- b[soc_rows]

      ## (a) Add k cone auxiliary variables
      ## First variable (head, t): lb=0 (SOC requires t >= 0)
      ## Rest: lb=-Inf
      cone_cols <- integer(k)
      cone_cols[1L] <- xpress::xprs_newcol(prob, 0, Inf, "C",
                                            paste0("cX", i_cone, "_0"), 0)
      if (k > 1L) {
        for (j in 2:k) {
          cone_cols[j] <- xpress::xprs_newcol(prob, -Inf, Inf, "C",
                                               paste0("cX", i_cone, "_", j - 1L), 0)
        }
      }

      ## (b) Add k transformation rows (equality)
      ## Row i: A_soc[i, :] * x_orig + cone_cols[i] = b_soc[i]
      ## First add the rows with original variable coefficients
      A_soc_csr <- methods::as(A_soc, "RsparseMatrix")
      xpress::addrows(prob,
        rowtype = rep("E", k),
        rhs     = b_soc,
        start   = A_soc_csr@p,
        colind  = A_soc_csr@j,
        rowcoef = A_soc_csr@x
      )

      ## Then add the identity block for cone variables via chgmcoef
      ## The transformation rows just added are the last k rows
      trans_row_start <- xpress::getintattrib(prob, xpress:::ROWS) - k
      xpress::chgmcoef(prob,
        rowind  = seq(trans_row_start, trans_row_start + k - 1L),
        colind  = cone_cols,
        rowcoef = rep(1, k)
      )

      ## (c) Quadratic SOC constraint (only if k > 1)
      ## v0^2 - v1^2 - ... - v_{k-1}^2 >= 0
      if (k > 1L) {
        xpress::addrows(prob, rowtype = "G", rhs = 0,
                         start = 0L, colind = integer(0), rowcoef = numeric(0))
        qc_row <- xpress::getintattrib(prob, xpress:::ROWS) - 1L
        xpress::addqmatrix(prob, qc_row,
          rowqcol1 = cone_cols,
          rowqcol2 = cone_cols,
          rowqcoef = c(1, rep(-1, k - 1L))
        )
      }

      currow <- currow + k
    }
  }

  ## -- Set solver controls ------------------------------------------------
  defaults <- list()
  if (!"BARGAPTARGET" %in% toupper(names(solver_opts)))
    defaults[["BARGAPTARGET"]] <- 1e-30
  if (!"FEASTOL" %in% toupper(names(solver_opts)))
    defaults[["FEASTOL"]] <- 1e-9
  .xpress_set_controls(prob, defaults)

  if (verbose) {
    .xpress_set_controls(prob, list(OUTPUTLOG = 1L, MIPLOG = 2L, LPLOG = 1L))
  } else {
    .xpress_set_controls(prob, list(OUTPUTLOG = 0L, MIPLOG = 0L, LPLOG = 0L))
  }

  user_opts <- solver_opts[!names(solver_opts) %in% c("write_mps")]
  if (length(user_opts) > 0L) .xpress_set_controls(prob, user_opts)

  if (!is.null(solver_opts[["write_mps"]])) {
    tryCatch(xpress::writeprob(prob, solver_opts[["write_mps"]]),
             error = function(e) NULL)
  }

  ## -- MIP warm-start -----------------------------------------------------
  if (warm_start && is_mip) {
    mip_indices <- c(bool_idx, int_idx) - 1L
    dots <- list(...)
    solver_cache <- dots[["solver_cache"]]
    if (!is.null(solver_cache) && exists(XPRESS_SOLVER, envir = solver_cache)) {
      cached <- get(XPRESS_SOLVER, envir = solver_cache)
      if (!is.null(cached$x) && length(cached$x) >= max(c(bool_idx, int_idx))) {
        warm_vals <- cached$x[c(bool_idx, int_idx)]
        tryCatch(
          xpress::addmipsol(prob, warm_vals, mip_indices, "warmstart"),
          error = function(e) NULL
        )
      }
    }
  }

  ## -- Optimize -----------------------------------------------------------
  tryCatch(
    xpress::xprs_optimize(prob),
    error = function(e) NULL
  )

  ## -- Extract results ----------------------------------------------------
  result <- list()

  sol_status <- xpress::getintattrib(prob, xpress:::SOLSTATUS)
  result$status <- as.character(sol_status)
  result$.nrows_linear <- nrows_linear
  result$.nrowsEQ <- nrowsEQ
  result$.is_mip <- is_mip
  result$.nvars <- nvars

  status <- XPRESS_STATUS_MAP[[result$status]]
  if (is.null(status)) status <- SOLVER_ERROR

  if (status %in% SOLUTION_PRESENT) {
    ## Primal: only first nvars entries (exclude cone auxiliary variables)
    sol <- xpress::getsolution(prob)$x
    result$x <- sol[seq_len(nvars)]

    result$objval <- if (is_mip) {
      xpress::getdblattrib(prob, xpress:::MIPOBJVAL)
    } else {
      xpress::getdblattrib(prob, xpress:::LPOBJVAL)
    }

    ## Duals: only linear constraint duals (not MIP, not SOC rows)
    ## Scope explicitly to first nrows_linear rows
    if (!is_mip && nrows_linear > 0L) {
      result$duals <- xpress::getduals(prob, first = 0L,
                                        last = nrows_linear - 1L)$duals
    }
  }

  tryCatch(
    result$solve_time <- xpress::getdblattrib(prob, xpress:::TIME),
    error = function(e) NULL
  )

  ## Cache for future warm-starts
  dots <- list(...)
  solver_cache <- dots[["solver_cache"]]
  if (!is.null(solver_cache)) {
    assign(XPRESS_SOLVER, result, envir = solver_cache)
  }

  result
}

# -- reduction_invert ----------------------------------------------------------
## CVXPY SOURCE: xpress_conif.py invert
## Dual sign: negate ALL linear duals (matching CVXPY)

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

  status <- XPRESS_STATUS_MAP[[solution$status]]
  if (is.null(status)) status <- SOLVER_ERROR

  if (!is.null(solution$solve_time))
    attr_list[[RK_SOLVE_TIME]] <- solution$solve_time

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

    ## Primal variables (already trimmed to nvars in solve_via_data)
    primal_vars <- list()
    primal_vars[[as.character(inverse_data[[SOLVER_VAR_ID]])]] <- solution$x

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

    if (!is_mip && !is.null(solution$duals)) {
      y <- -solution$duals
      nrowsEQ <- solution$.nrowsEQ

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

      ineq_dual <- if (nrowsEQ < length(y)) {
        get_dual_values(
          y[(nrowsEQ + 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, XPRESS_Conic_Solver) <- function(x, ...) {
  cat("XPRESS_Conic_Solver()\n")
  invisible(x)
}