R/279_reductions_solvers_conic_solvers_cplex_conif.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: reductions/solvers/conic_solvers/cplex_conif.py
## CPLEX conic solver interface via Rcplex.
##
## CVXPY 1.9 parity notes:
##  - Mirrors cplex_conif.py for LP/SOCP/MI-LP/MI-SOCP: Zero and NonNeg are
##    linear rows; SOC blocks are represented by auxiliary variables plus one
##    convex quadratic constraint per cone.
##  - SOC duals come from the duals of each cone's linking rows, which
##    `Rcplex_solve_QCP()` returns with the other row duals in `extra$lambda`,
##    exactly as upstream reads them from `get_dual_values()`.
##  - Upstream forces `preprocessing.qcpduals` (cplex_conif.py:369-373);
##    Rcplex has no control for it, so CPLEX's default applies, and if CPLEX
##    computes no row duals `extra$lambda` is NA and no duals are returned.


# -- 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()

## CPX status constants, transcribed from
## CPLEX_Studio2211/cplex/include/ilcplex/cpxconst.h. Named, not inlined,
## because the flat table these replace carried COMMENTS that were wrong -- and
## a wrong comment on a magic number is invisible. Checked against the header:
## the old table labelled 10 "abort_user" (it is ABORT_IT_LIM), 11
## "abort_iteration_limit" (ABORT_TIME_LIM), 12 "abort_time_limit"
## (ABORT_OBJ_LIM), 13 "abort_dettime_limit" (ABORT_USER), 6 "feasible"
## (NUM_BEST; FEASIBLE is 23) and 115 "MIP_infeasible_or_unbounded"
## (MIP_OPTIMAL_INFEAS; MIP INForUNBD is 119). The labels were shifted by one
## through the 10-13 block, and the mappings followed the labels.
CPX_STAT_OPTIMAL              <- 1L
CPX_STAT_UNBOUNDED            <- 2L
CPX_STAT_INFEASIBLE           <- 3L
CPX_STAT_INForUNBD            <- 4L
CPX_STAT_OPTIMAL_INFEAS       <- 5L
CPX_STAT_NUM_BEST             <- 6L
CPX_STAT_ABORT_IT_LIM         <- 10L
CPX_STAT_ABORT_TIME_LIM       <- 11L
CPX_STAT_ABORT_OBJ_LIM        <- 12L
CPX_STAT_ABORT_USER           <- 13L
CPX_STAT_ABORT_PRIM_OBJ_LIM   <- 21L
CPX_STAT_ABORT_DUAL_OBJ_LIM   <- 22L
CPX_STAT_FEASIBLE             <- 23L
CPX_STAT_FIRSTORDER           <- 24L
CPX_STAT_ABORT_DETTIME_LIM    <- 25L
CPX_STAT_BENDERS_NUM_BEST     <- 41L
CPXMIP_OPTIMAL                <- 101L
CPXMIP_OPTIMAL_TOL            <- 102L
CPXMIP_INFEASIBLE             <- 103L
CPXMIP_SOL_LIM                <- 104L
CPXMIP_NODE_LIM_FEAS          <- 105L
CPXMIP_NODE_LIM_INFEAS        <- 106L
CPXMIP_TIME_LIM_FEAS          <- 107L
CPXMIP_TIME_LIM_INFEAS        <- 108L
CPXMIP_FAIL_FEAS              <- 109L
CPXMIP_FAIL_INFEAS            <- 110L
CPXMIP_MEM_LIM_FEAS           <- 111L
CPXMIP_MEM_LIM_INFEAS         <- 112L
CPXMIP_ABORT_FEAS             <- 113L
CPXMIP_ABORT_INFEAS           <- 114L
CPXMIP_OPTIMAL_INFEAS         <- 115L
CPXMIP_FAIL_FEAS_NO_TREE      <- 116L
CPXMIP_FAIL_INFEAS_NO_TREE    <- 117L
CPXMIP_UNBOUNDED              <- 118L
CPXMIP_INForUNBD              <- 119L
CPXMIP_FEASIBLE_RELAXED_SUM   <- 120L
CPXMIP_OPTIMAL_RELAXED_SUM    <- 121L
CPXMIP_FEASIBLE_RELAXED_INF   <- 122L
CPXMIP_OPTIMAL_RELAXED_INF    <- 123L
CPXMIP_FEASIBLE_RELAXED_QUAD  <- 124L
CPXMIP_OPTIMAL_RELAXED_QUAD   <- 125L
CPXMIP_ABORT_RELAXED          <- 126L
CPXMIP_FEASIBLE               <- 127L
CPXMIP_POPULATESOL_LIM        <- 128L
CPXMIP_OPTIMAL_POPULATED      <- 129L
CPXMIP_OPTIMAL_POPULATED_TOL  <- 130L
CPXMIP_DETTIME_LIM_FEAS       <- 131L
CPXMIP_DETTIME_LIM_INFEAS     <- 132L
CPX_STAT_FEASIBLE_RELAXED_SUM  <- 14L
CPX_STAT_OPTIMAL_RELAXED_SUM   <- 15L
CPX_STAT_FEASIBLE_RELAXED_INF  <- 16L
CPX_STAT_OPTIMAL_RELAXED_INF   <- 17L
CPX_STAT_FEASIBLE_RELAXED_QUAD <- 18L
CPX_STAT_OPTIMAL_RELAXED_QUAD  <- 19L


# -- .cplex_normalize_solstat ------------------------------------------------
## CVXPY SOURCE: cplex_conif.py:75-144 (`_handle_solve_status`) -- collapse the
## MIP status codes onto their LP equivalents so the ladder below only has to
## reason about one vocabulary.
.cplex_normalize_solstat <- function(solstat) {
  if (solstat == CPXMIP_OPTIMAL) return(CPX_STAT_OPTIMAL)
  if (solstat == CPXMIP_INFEASIBLE) return(CPX_STAT_INFEASIBLE)
  if (solstat %in% c(CPXMIP_TIME_LIM_FEAS, CPXMIP_TIME_LIM_INFEAS))
    return(CPX_STAT_ABORT_TIME_LIM)
  if (solstat %in% c(CPXMIP_DETTIME_LIM_FEAS, CPXMIP_DETTIME_LIM_INFEAS))
    return(CPX_STAT_ABORT_DETTIME_LIM)
  if (solstat %in% c(CPXMIP_ABORT_FEAS, CPXMIP_ABORT_INFEAS))
    return(CPX_STAT_ABORT_USER)
  if (solstat == CPXMIP_OPTIMAL_INFEAS) return(CPX_STAT_OPTIMAL_INFEAS)
  if (solstat == CPXMIP_INForUNBD) return(CPX_STAT_INForUNBD)
  if (solstat == CPXMIP_UNBOUNDED) return(CPX_STAT_UNBOUNDED)
  ## feasopt results. CVXR never calls feasopt, so reaching one means the model
  ## was built by something other than this interface. Upstream raises
  ## AssertionError (cplex_conif.py:110-123); same here.
  if (solstat %in% c(CPX_STAT_FEASIBLE_RELAXED_SUM, CPXMIP_FEASIBLE_RELAXED_SUM,
                     CPX_STAT_OPTIMAL_RELAXED_SUM, CPXMIP_OPTIMAL_RELAXED_SUM,
                     CPX_STAT_FEASIBLE_RELAXED_INF, CPXMIP_FEASIBLE_RELAXED_INF,
                     CPX_STAT_OPTIMAL_RELAXED_INF, CPXMIP_OPTIMAL_RELAXED_INF,
                     CPX_STAT_FEASIBLE_RELAXED_QUAD, CPXMIP_FEASIBLE_RELAXED_QUAD,
                     CPX_STAT_OPTIMAL_RELAXED_QUAD, CPXMIP_OPTIMAL_RELAXED_QUAD)) {
    cli_abort("feasopt status encountered: {.val {solstat}}.")
  }
  ## Conflict-refiner results, 30-39. Same reasoning (cplex_conif.py:124-135).
  if (solstat >= 30L && solstat <= 39L) {
    cli_abort("conflict refiner status encountered: {.val {solstat}}.")
  }
  if (solstat %in% c(CPX_STAT_FEASIBLE, CPXMIP_FEASIBLE)) return(CPX_STAT_FEASIBLE)
  if (solstat == CPX_STAT_BENDERS_NUM_BEST) return(CPX_STAT_NUM_BEST)
  solstat
}


# -- .cplex_get_status -------------------------------------------------------
## CVXPY SOURCE: cplex_conif.py:147-204 (`get_status`).
##
## This replaces a FLAT TABLE, and the difference is not stylistic: upstream
## reads two secondary conditions -- `model.solution.is_primal_feasible()` and
## `is_dual_feasible()` -- that decide the answer for eleven status codes. A
## table cannot express "iteration limit, but a feasible point was found", so
## CVXR answered it wrongly. Measured, `psolve(prob, solver = "CPLEX",
## itlim = 1)` on a 60-variable LP: CPLEX returns status 10 WITH a feasible
## primal and duals, and CVXR reported
##     Solver "CPLEX" failed.
## discarding a usable solution. CVXPY returns `optimal_inaccurate` with it.
##
## `pfeas` / `dfeas`: Rcplex exposes no is_primal_feasible()/is_dual_feasible(),
## so they are read off the returned solution -- a primal iterate is present iff
## `xopt` is non-NULL and non-NA (verified: itlim = 1 yields a full xopt, a
## tilim so small the solve never starts yields NA), and duals likewise from
## `extra$lambda`. Upstream notes at cplex_conif.py:150 that dfeas is always
## false for a MIP; CVXR passes `dfeas = FALSE` for the MIP path rather than
## inferring it.
.cplex_get_status <- function(solstat, pfeas, dfeas) {
  if (is.null(solstat) || length(solstat) != 1L || is.na(solstat)) {
    return(SOLVER_ERROR)
  }
  solstat <- .cplex_normalize_solstat(as.integer(solstat))

  if (solstat %in% c(CPXMIP_NODE_LIM_INFEAS, CPXMIP_FAIL_INFEAS,
                     CPXMIP_MEM_LIM_INFEAS, CPXMIP_FAIL_INFEAS_NO_TREE,
                     CPX_STAT_NUM_BEST)) {
    return(SOLVER_ERROR)
  }
  if (solstat %in% c(CPX_STAT_ABORT_USER, CPX_STAT_ABORT_IT_LIM,
                     CPX_STAT_ABORT_TIME_LIM, CPX_STAT_ABORT_DETTIME_LIM,
                     CPX_STAT_ABORT_OBJ_LIM, CPX_STAT_ABORT_PRIM_OBJ_LIM,
                     CPX_STAT_ABORT_DUAL_OBJ_LIM, CPXMIP_ABORT_RELAXED,
                     CPX_STAT_FIRSTORDER)) {
    return(if (pfeas) OPTIMAL_INACCURATE else SOLVER_ERROR)
  }
  if (solstat %in% c(CPXMIP_NODE_LIM_FEAS, CPXMIP_SOL_LIM,
                     CPXMIP_POPULATESOL_LIM, CPXMIP_FAIL_FEAS,
                     CPXMIP_MEM_LIM_FEAS, CPXMIP_FAIL_FEAS_NO_TREE,
                     CPX_STAT_FEASIBLE)) {
    return(if (dfeas) OPTIMAL else OPTIMAL_INACCURATE)
  }
  if (solstat %in% c(CPX_STAT_OPTIMAL, CPXMIP_OPTIMAL_TOL,
                     CPX_STAT_OPTIMAL_INFEAS, CPXMIP_OPTIMAL_POPULATED,
                     CPXMIP_OPTIMAL_POPULATED_TOL)) {
    return(OPTIMAL)
  }
  if (solstat %in% c(CPX_STAT_INFEASIBLE, CPXMIP_OPTIMAL_RELAXED_SUM,
                     CPXMIP_OPTIMAL_RELAXED_INF, CPXMIP_OPTIMAL_RELAXED_QUAD)) {
    return(INFEASIBLE)
  }
  if (solstat %in% c(CPXMIP_FEASIBLE_RELAXED_QUAD, CPXMIP_FEASIBLE_RELAXED_INF,
                     CPXMIP_FEASIBLE_RELAXED_SUM)) {
    return(SOLVER_ERROR)
  }
  if (solstat == CPX_STAT_UNBOUNDED) return(UNBOUNDED)
  if (solstat == CPX_STAT_INForUNBD) return(INFEASIBLE_OR_UNBOUNDED)
  SOLVER_ERROR
}


# -- .cplex_solution_feasibility ---------------------------------------------
## The pfeas / dfeas pair for an Rcplex result. `is_mip` forces dfeas FALSE,
## mirroring the note at cplex_conif.py:150.
.cplex_solution_feasibility <- function(solution, is_mip = FALSE) {
  xopt <- solution$xopt
  lambda <- solution$extra$lambda
  list(
    pfeas = !is.null(xopt) && length(xopt) > 0L && !anyNA(xopt),
    dfeas = !isTRUE(is_mip) && !is.null(lambda) && length(lambda) > 0L &&
            !anyNA(lambda)
  )
}

# -- .cplex_set_parameters -----------------------------------------------------
## CVXPY SOURCE: cplex_conif.py:42-62 (`set_parameters`), shared by both
## interfaces exactly as cplex_qpif.py:7-11 imports it from here. Upstream
## consumes `cplex_params` (a dict of CPLEX parameter paths, ValueError on a
## name the model does not have, :45-55), `cplex_filename` (`model.write`,
## :56-60) and raises `ValueError("invalid keyword-argument ...")` on anything
## else (:61-62). `reoptimize` is popped by the callers before the call
## (cplex_conif.py:376, cplex_qpif.py:176).
##
## Rcplex takes a FLAT control list, not CPLEX's dotted parameter tree
## (completeness_known.tsv: `set_parameters` n/a-for-r for the carrier), so a
## flat option IS the R-native spelling of a `cplex_params` entry and both are
## routed into `control`. The names Rcplex accepts are those of its default
## list (Rcplex 0.3.8 R/control.R `check.Rcplex.control`; the C dispatch in
## src/Rcplex_params.c warns "Unknown CPLEX parameter %s. Ignoring it." and
## DROPS any other name). Before 9156 every option was forwarded and that
## warning was the only signal, so `psolve(prob, solver = "CPLEX",
## reoptimize = TRUE)` silently solved without the re-solve. Now an unknown
## name aborts, as upstream and as the MOSEK twin (mosek_conif.R).
##
## `cplex_filename` cannot be honored: Rcplex exposes no model-writing entry
## point (no CPXwriteprob binding in src/). It aborts with that reason rather
## than being dropped.
.RCPLEX_CONTROL_NAMES <- c(
  "trace", "method", "preind", "aggind", "itlim", "epagap", "epgap", "tilim",
  "disjcuts", "mipemphasis", "cliques", "nodesel", "probe", "varsel",
  "flowcovers", "solnpoolagap", "solnpoolgap", "solnpoolintensity",
  "maxcalls", "round"
)

.cplex_set_parameters <- function(solver_opts, verbose = FALSE) {
  known <- .RCPLEX_CONTROL_NAMES   # cli literals may not start with a dot
  control <- list(trace = if (verbose) 1L else 0L)
  reoptimize <- isTRUE(solver_opts[["reoptimize"]])
  solver_opts[["reoptimize"]] <- NULL

  cplex_params <- solver_opts[["cplex_params"]]
  solver_opts[["cplex_params"]] <- NULL
  for (param in names(cplex_params)) {
    if (!(param %in% .RCPLEX_CONTROL_NAMES)) {
      ## CVXPY SOURCE: cplex_conif.py:51-54
      cli_abort(c(
        "Invalid CPLEX parameter, value pair ({.val {param}}, {.val {cplex_params[[param]]}}).",
        "i" = "Rcplex accepts {.val {known}}."
      ))
    }
    control[[param]] <- cplex_params[[param]]
  }

  cplex_filename <- solver_opts[["cplex_filename"]]
  solver_opts[["cplex_filename"]] <- NULL
  if (!is.null(cplex_filename) && !isFALSE(cplex_filename) &&
      !identical(cplex_filename, "")) {
    cli_abort(c(
      "{.arg cplex_filename} is not supported for solver {.val CPLEX}.",
      "i" = "Rcplex exposes no model-writing API."
    ))
  }

  for (opt_name in names(solver_opts)) {
    if (!(opt_name %in% .RCPLEX_CONTROL_NAMES)) {
      ## CVXPY SOURCE: cplex_conif.py:61-62
      cli_abort(c(
        "Invalid option {.arg {opt_name}} for solver {.val CPLEX}.",
        "i" = "Rcplex control parameters are {.val {known}}; {.arg cplex_params} may carry them as a named list, and {.arg reoptimize} re-solves with presolve off on an ambiguous status."
      ))
    }
    control[[opt_name]] <- solver_opts[[opt_name]]
  }

  list(control = control, reoptimize = reoptimize)
}

# -- .cplex_run ----------------------------------------------------------------
## CVXPY SOURCE: cplex_conif.py:379-396 and cplex_qpif.py:179-199 -- the solve
## itself, timed with `model.get_time()` (wall clock) into SOLVE_TIME, and the
## `reoptimize` re-solve: when the first status is INFEASIBLE_OR_UNBOUNDED and
## the flag is set, presolve is switched off (`preprocessing.presolve.set(0)`,
## Rcplex control `preind = 0L`) and the model is solved again, the second
## run's time added to the first. Shared by both interfaces; `call_fn` is the
## Rcplex call closed over the problem data, taking the control list.
##
## Rcplex returns no timing (Rcplex.Rd Value: xopt, obj, status, extra), so
## the wall clock around the call stands in for `model.get_time()`, which is
## also wall clock by CPLEX's default clocktype. Rcplex prints "CPLEX
## environment opened" to stderr on first call (C-level init message); it is
## captured when not verbose.
.cplex_run <- function(call_fn, control, verbose, reoptimize, is_mip) {
  run_once <- function(control) {
    t0 <- proc.time()[["elapsed"]]
    res <- tryCatch({
      if (!verbose) {
        r <- NULL
        utils::capture.output(r <- call_fn(control), type = "message")
        r
      } else {
        call_fn(control)
      }
    },
      error = function(e) {
        list(xopt = NA, obj = NA, status = -1L,
             extra = list(lambda = NA, slack = NA),
             .error = conditionMessage(e))
      }
    )
    ## upstream records cputime only when solve() returned (:182-185)
    if (is.null(res$.error)) res$.solve_time <- proc.time()[["elapsed"]] - t0
    res
  }

  result <- run_once(control)
  if (reoptimize) {
    feas <- .cplex_solution_feasibility(result, is_mip = is_mip)
    ambiguous <- identical(
      .cplex_get_status(result$status, feas$pfeas, feas$dfeas),
      INFEASIBLE_OR_UNBOUNDED
    )
    if (ambiguous) {
      first_time <- result$.solve_time
      control$preind <- 0L
      result <- run_once(control)
      if (!is.null(result$.solve_time)) {
        result$.solve_time <- result$.solve_time + first_time
      }
    }
  }
  result
}

# -- CPLEX_Conic_Solver class ----------------------------------------
## CVXPY SOURCE: cplex_conif.py lines 207-215

#' @keywords internal
CPLEX_Conic_Solver <- new_class("CPLEX_Conic_Solver", parent = ConicSolver,
  package = "CVXR",
  constructor = function() {
    if (FALSE) new_object(S7_object())  ## S7 static-check guard
    .fast_new(CPLEX_Conic_Solver, S7_object(),
      .cache = new.env(parent = emptyenv()),
      MIP_CAPABLE = TRUE,
      BOUNDED_VARIABLES = TRUE,
      PSD_TRIANGLE_KIND = NA_character_,
      PSD_SQRT2_SCALING = NA,
      SUPPORTED_CONSTRAINTS = list(Zero, NonNeg, SOC),
      EXP_CONE_ORDER = NULL,
      REQUIRES_CONSTR = FALSE
    )
  }
)

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

# -- reduction_accepts ------------------------------------------------

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

# -- reduction_apply --------------------------------------------------
## CVXPY SOURCE: cplex_conif.py apply() inherits ConicSolver structure.

method(reduction_apply, CPLEX_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
  solver_data[[LOWER_BOUNDS]] <- data[[LOWER_BOUNDS]]
  solver_data[[UPPER_BOUNDS]] <- data[[UPPER_BOUNDS]]
  if (!is.null(data[[SD_P]])) solver_data[[SD_P]] <- data[[SD_P]]

  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: cplex_conif.py lines 289-396 and add_model_soc_constr()
## lines 440-522. Rcplex takes one sparse `Amat` plus a QCP `QC` list:
##   min 0.5 x'Qx + c'x
##   s.t. Amat x {sense} bvec
##        q_i'x + x'Q_i x {L/G} r_i

method(solve_via_data, CPLEX_Conic_Solver) <- function(x, data, warm_start = FALSE, verbose = FALSE,
                                                        solver_opts = list(), ...) {
  .require_solver_package(CPLEX_SOLVER)

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

  zero_dim <- dims@zero
  nonneg_dim <- dims@nonneg
  linear_dim <- zero_dim + nonneg_dim
  soc_total <- sum(dims@soc)
  n_soc_aux <- soc_total
  n_total <- n_orig + n_soc_aux

  ## Linear constraints are ordered Zero, NonNeg, then SOC rows.
  if (linear_dim > 0L) {
    A_linear <- A[seq_len(linear_dim), , drop = FALSE]
    b_linear <- b[seq_len(linear_dim)]
    sense <- character(linear_dim)
    if (zero_dim > 0L) sense[seq_len(zero_dim)] <- "E"
    if (nonneg_dim > 0L) sense[(zero_dim + 1L):linear_dim] <- "L"
  } else {
    A_linear <- Matrix::sparseMatrix(i = integer(0), j = integer(0),
                                      x = numeric(0), dims = c(0L, n_orig))
    b_linear <- numeric(0)
    sense <- character(0)
  }

  soc_eq_rows <- list()
  soc_eq_rhs <- numeric(0)
  qc_Q <- list()
  qc_L <- list()
  qc_dir <- character(0)
  qc_b <- numeric(0)
  aux_offset <- 0L

  if (soc_total > 0L) {
    soc_start <- linear_dim
    for (k in dims@soc) {
      soc_rows <- (soc_start + 1L):(soc_start + k)
      A_soc <- A[soc_rows, , drop = FALSE]
      b_soc <- b[soc_rows]

      ## CVXPY SOURCE: cplex_conif.py:497-510 -- row A_i x + aux_i = b_i,
      ## sense "E". Its sign fixes the sign of the row's dual, which invert
      ## reads back as the SOC dual.
      for (i in seq_len(k)) {
        aux_col <- n_orig + aux_offset + i
        full_row <- Matrix::sparseMatrix(i = integer(0), j = integer(0),
                                          x = numeric(0), dims = c(1L, n_total))
        full_row[1L, seq_len(n_orig)] <- A_soc[i, , drop = TRUE]
        full_row[1L, aux_col] <- 1.0
        soc_eq_rows[[length(soc_eq_rows) + 1L]] <- full_row
        soc_eq_rhs <- c(soc_eq_rhs, b_soc[i])
      }

      ## Quadratic SOC row: sum(aux_tail^2) - aux_head^2 <= 0.
      aux_base <- n_orig + aux_offset
      q_idx <- aux_base + seq_len(k)
      q_val <- c(-1.0, rep(1.0, k - 1L))
      Qc <- Matrix::sparseMatrix(i = q_idx, j = q_idx, x = q_val,
                                  dims = c(n_total, n_total))
      qc_Q[[length(qc_Q) + 1L]] <- Qc
      qc_L[[length(qc_L) + 1L]] <- rep(0, n_total)
      qc_dir <- c(qc_dir, "L")
      qc_b <- c(qc_b, 0)

      aux_offset <- aux_offset + k
      soc_start <- soc_start + k
    }
  }

  ## Expand original linear rows to include SOC auxiliaries.
  if (linear_dim > 0L && n_soc_aux > 0L) {
    A_linear <- cbind(
      A_linear,
      Matrix::sparseMatrix(i = integer(0), j = integer(0), x = numeric(0),
                            dims = c(linear_dim, n_soc_aux))
    )
  }

  if (length(soc_eq_rows) > 0L) {
    A_soc_eq <- do.call(rbind, soc_eq_rows)
    Amat <- if (linear_dim > 0L) rbind(A_linear, A_soc_eq) else A_soc_eq
    bvec <- c(b_linear, soc_eq_rhs)
    sense <- c(sense, rep("E", length(soc_eq_rhs)))
  } else {
    Amat <- A_linear
    bvec <- b_linear
  }

  ## Rows whose duals invert maps: Zero, NonNeg, then the SOC rows.
  n_user_linear <- linear_dim + soc_total
  if (nrow(Amat) == 0L) {
    ## Rcplex LP/QP paths expect at least one row; keep it out of dual mapping.
    Amat <- Matrix::sparseMatrix(i = integer(0), j = integer(0),
                                  x = numeric(0), dims = c(1L, n_total))
    bvec <- 0
    sense <- "L"
  }

  ## Objective: Rcplex uses 0.5 x'Qx + c'x, matching CVXR's P convention.
  c_full <- c(c_vec, rep(0, n_soc_aux))
  Qmat <- data[[SD_P]]
  if (!is.null(Qmat)) {
    if (n_soc_aux > 0L) {
      Qmat <- Matrix::bdiag(
        Qmat,
        Matrix::sparseMatrix(i = integer(0), j = integer(0), x = numeric(0),
                              dims = c(n_soc_aux, n_soc_aux))
      )
    }
    if (!inherits(Qmat, "dgCMatrix")) {
      Qmat <- methods::as(methods::as(Qmat, "generalMatrix"), "CsparseMatrix")
    }
  }

  lb_orig <- data[[LOWER_BOUNDS]] %||% rep(-Inf, n_orig)
  ub_orig <- data[[UPPER_BOUNDS]] %||% rep(Inf, n_orig)
  lb <- c(lb_orig, rep(-Inf, n_soc_aux))
  ub <- c(ub_orig, rep(Inf, n_soc_aux))

  if (soc_total > 0L) {
    soc_aux_offset <- 0L
    for (k in dims@soc) {
      head_idx <- n_orig + soc_aux_offset + 1L
      lb[head_idx] <- 0
      soc_aux_offset <- soc_aux_offset + k
    }
  }

  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", n_total)
    if (length(bool_idx) > 0L) {
      for (idx in bool_idx) {
        vtype[idx] <- "B"
        lb[idx] <- max(lb[idx], 0)
        ub[idx] <- min(ub[idx], 1)
      }
    }
    if (length(int_idx) > 0L) {
      for (idx in int_idx) vtype[idx] <- "I"
    }
  }

  ## CVXPY SOURCE: cplex_conif.py:375-377 -- reoptimize popped, then
  ## set_parameters(); the solve and the re-solve are :379-396 (`.cplex_run`).
  params <- .cplex_set_parameters(solver_opts, verbose)

  QC <- if (length(qc_Q) > 0L) {
    list(QC = list(Q = qc_Q, L = qc_L), dir = qc_dir, b = qc_b)
  } else {
    NULL
  }

  .cplex_call <- function(control) {
    Rcplex::Rcplex_solve_QCP(
      cvec = as.double(c_full),
      Amat = Amat,
      bvec = as.double(bvec),
      Qmat = Qmat,
      QC = QC,
      lb = as.double(lb),
      ub = as.double(ub),
      control = control,
      objsense = "min",
      sense = sense,
      vtype = vtype
    )
  }

  result <- .cplex_run(.cplex_call, params$control, verbose,
                       params$reoptimize, is_mip)

  result$.n_orig <- n_orig
  result$.n_user_linear <- n_user_linear
  result$.is_mip <- is_mip
  result
}

# -- reduction_invert -------------------------------------------------

method(reduction_invert, CPLEX_Conic_Solver) <- function(x, solution, inverse_data, ...) {
  attr_list <- list()
  attr_list[[RK_EXTRA_STATS]] <- solution
  ## CVXPY SOURCE: cplex_conif.py:260-261 -- SOLVE_TIME when the solve ran.
  if (!is.null(solution$.solve_time)) {
    attr_list[[RK_SOLVE_TIME]] <- solution$.solve_time
  }

  ## CVXPY SOURCE: cplex_conif.py:147-204 -- the decision procedure, not a
  ## table: eleven status codes are resolved by whether a primal (or dual)
  ## iterate is present. See `.cplex_get_status`.
  feas <- .cplex_solution_feasibility(solution,
                                      is_mip = isTRUE(inverse_data[["is_mip"]]))
  status <- .cplex_get_status(solution$status, feas$pfeas, feas$dfeas)

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

    primal_vars <- list()
    primal_vars[[as.character(inverse_data[[SOLVER_VAR_ID]])]] <-
      solution$xopt[seq_len(solution$.n_orig)]

    dual_vars <- list()
    is_mip <- isTRUE(solution$.is_mip)

    ## CVXPY SOURCE: cplex_conif.py:274-283 -- y = -(row duals), in the order
    ## the rows were added (Zero, NonNeg, then each SOC's rows), mapped onto
    ## EQ_CONSTR followed by NEQ_CONSTR (NonNeg, then SOC). Rcplex returns the
    ## row duals as `extra$lambda` (CPXgetpi), NA when CPLEX computed none.
    if (!is_mip && !is.null(solution$extra$lambda) &&
        !any(is.na(solution$extra$lambda))) {
      y <- -solution$extra$lambda
      if (solution$.n_user_linear < length(y)) {
        y <- y[seq_len(solution$.n_user_linear)]
      }
      dual_vars <- get_dual_values(
        y,
        extract_dual_value,
        c(inverse_data[[SOLVER_EQ_CONSTR]], inverse_data[[SOLVER_NEQ_CONSTR]])
      )
    }

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

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