R/240_reductions_solvers_conic_solvers_gurobi_conif.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: reductions/solvers/conic_solvers/gurobi_conif.py
## Gurobi conic solver interface via R gurobi package
##
## Handles LP, SOCP, MI-LP, MI-SOCP via the conic path.
## QP/MIQP problems are handled by Gurobi_QP_Solver (qp_solvers/gurobi_qp_solver.R).
## Uses the standard ConicSolver pipeline (ConeMatrixStuffing ->
## format_constraints -> ConicSolver.reduction_apply).
##
## Key design:
## - Zero + NonNeg constraints -> model$A with sense "="/"<"
## - SOC constraints -> auxiliary variables + model$quadcon
## - QP via model$Q (lower-tri, NO 1/2 factor -- Gurobi uses x'Qx + c'x)
## - MIP via model$vtype ("B"=binary, "I"=integer, "C"=continuous)
## - R gurobi returns STRING status codes (e.g., "OPTIMAL")
## - Dual sign: negate all linear duals (matching CVXPY gurobi_conif.py)


# -- Gurobi status map ------------------------------------------
## CVXPY SOURCE: gurobi_conif.py STATUS_MAP
## R gurobi returns string status codes (not integer like gurobipy).

GUROBI_STATUS_MAP <- list(
  "OPTIMAL"         = OPTIMAL,
  "INFEASIBLE"      = INFEASIBLE,
  "INF_OR_UNBD"     = INFEASIBLE_OR_UNBOUNDED,
  "UNBOUNDED"       = UNBOUNDED,
  "NUMERIC"         = SOLVER_ERROR,
  "ITERATION_LIMIT" = USER_LIMIT,
  "NODE_LIMIT"      = USER_LIMIT,
  "TIME_LIMIT"      = USER_LIMIT,
  "SOLUTION_LIMIT"  = USER_LIMIT,
  "INTERRUPTED"     = USER_LIMIT,
  "SUBOPTIMAL"      = OPTIMAL_INACCURATE,
  "INPROGRESS"      = USER_LIMIT,
  "USER_OBJ_LIMIT"  = USER_LIMIT,
  "WORK_LIMIT"      = USER_LIMIT,
  "MEM_LIMIT"       = USER_LIMIT
)

# -- Gurobi_Conic_Solver class ----------------------------------------
## Supports Zero, NonNeg, SOC constraints.
## SOC is handled via quadratic constraints (model$quadcon).

#' @keywords internal
Gurobi_Conic_Solver <- new_class("Gurobi_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, Gurobi_Conic_Solver) <- function(x) GUROBI_SOLVER

# -- reduction_accepts --------------------------------------------
## Override: Gurobi conic only handles Zero + NonNeg + SOC.
## Rejects ExpCone, PSD, PowCone3D constraints.

method(reduction_accepts, Gurobi_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 and SOC tracking.
## Uses standard conic path via format_constraints for
## Zero/NonNeg/SOC restructuring.

method(reduction_apply, Gurobi_Conic_Solver) <- function(x, problem, ...) {
  data <- problem  ## data list from ConeMatrixStuffing
  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

  ## Split into eq (Zero) and ineq (NonNeg, SOC)
  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"]]
  )

  ## Format constraints: restructure A and b
  formatted <- format_constraints(
    constraints, data[[SD_A]], data[[SD_B]],
    exp_cone_order = x@EXP_CONE_ORDER
  )

  ## Build solver data dict
  ## Convention: solver sees A*x + s = b, s in K
  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[[SD_P]])) solver_data[[SD_P]] <- data[[SD_P]]

  ## Pass through MIP indices if present
  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 ----------------------------------------------
## Converts CVXR solver data (A, b, c, dims, P) to Gurobi model
## and calls gurobi::gurobi().
##
## SOC constraints are translated to Gurobi quadratic constraints:
## For each SOC block (t, x1..xk): add auxiliary vars, link via
## equality constraints, then add x1^2+...+xk^2 <= t^2.

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

  dots <- list(...)
  solver_cache <- dots[["solver_cache"]]

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

  ## -- Row splitting ----------------------------------------------
  ## After format_constraints, rows are ordered:
  ## [Zero (eq)] [NonNeg (ineq)] [SOC blocks]
  zero_dim   <- dims@zero
  nonneg_dim <- dims@nonneg
  linear_dim <- zero_dim + nonneg_dim
  soc_total  <- sum(dims@soc)
  total_rows <- linear_dim + soc_total

  ## -- Linear constraints -> model$A, sense, rhs ------------------
  ## After conic path negation: solver_data$A = -formatted_A
  ## For Zero:   -(-A)*x = b  ->  A*x = b  (sense "=")
  ## For NonNeg: -(A)*x = b   ->  -A*x <= b (sense "<")
  ##
  ## But format_constraints already handles:
  ##   Zero -> -A (negated), NonNeg -> A (identity)
  ## Then solver_data$A = -formatted_A, so:
  ##   Zero rows of solver_data$A = A (positive)
  ##   NonNeg rows of solver_data$A = -A (negated)
  ##
  ## Gurobi format: A*x sense rhs
  ## For Zero: A_z * x = b_z
  ## For NonNeg: solver sees -A*x + s = b, s >= 0 -> -A*x <= b -> A_gurobi*x <= b
  ##   solver_data$A for NonNeg rows = -A_raw, so A_gurobi = solver_data$A, sense = "<"

  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)] <- "="
    }
    if (nonneg_dim > 0L) {
      sense[(zero_dim + 1L):linear_dim] <- "<"
    }
  } 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 constraints -> quadcon ----------------------------------
  ## For each SOC block of size k = dims@soc[i]:
  ##   SOC rows in A: k rows, first is t, rest are x1..x_{k-1}
  ##   Create k auxiliary variables (aux_t, aux_x1, ..., aux_{k-1})
  ##   Equality constraint: aux_i = b_soc[i] - A_soc[i,:] * x  (note: A already has
  ##     the conic negation applied, so aux_i = b_soc[i] + A_orig_soc[i,:] * x)
  ##   Quadratic constraint: aux_x1^2 + ... + aux_{k-1}^2 <= aux_t^2
  ##
  ## We add the auxiliary variable equalities as extra rows in the linear
  ## constraint matrix, and store quadratic constraints in model$quadcon.

  n_soc_aux <- soc_total  # total auxiliary variables needed
  n_total_vars <- n_orig + n_soc_aux

  quadcon_list <- list()
  soc_eq_A_rows <- list()
  soc_eq_rhs <- numeric(0)
  soc_eq_sense <- character(0)
  aux_offset <- 0L  # offset into auxiliary variables

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

      ## Build equality constraints: aux_i = b_soc[i] - A_soc[i,:] * x
      ## (A_soc = solver_data$A = -formatted_A; aux = slack = b - A*x)
      ## Equivalently: -A_soc[i,:] * x - aux_i = -b_soc[i]
      ## In expanded form: [-A_soc | -I_aux] [x; aux]' = -b_soc
      for (i in seq_len(k)) {
        ## Row with NEGATED original vars from A_soc and -1 for the aux var
        aux_col_idx <- n_orig + aux_offset + i  # 1-based column index
        row_orig <- A_soc[i, , drop = FALSE]  # 1 x n_orig sparse

        ## Extend to n_total_vars columns by adding the -1 for aux variable
        aux_part <- Matrix::sparseMatrix(
          i = 1L, j = aux_col_idx,
          x = -1.0, dims = c(1L, n_total_vars)
        )
        ## Combine: copy NEGATED orig columns into full-width row
        full_row <- Matrix::sparseMatrix(
          i = integer(0), j = integer(0), x = numeric(0),
          dims = c(1L, n_total_vars)
        )
        full_row[1, seq_len(n_orig)] <- -row_orig[1, ]  # NEGATE
        full_row <- full_row + aux_part

        soc_eq_A_rows <- c(soc_eq_A_rows, list(full_row))
        soc_eq_rhs <- c(soc_eq_rhs, -b_soc[i])
        soc_eq_sense <- c(soc_eq_sense, "=")
      }

      ## Quadratic constraint: sum(aux_x_i^2) <= aux_t^2
      ## => aux_x_1^2 + ... + aux_x_{k-1}^2 - aux_t^2 <= 0
      ## Qc matrix: diagonal with +1 for body vars, -1 for head var
      aux_base <- n_orig + aux_offset  # 0-based offset
      qc_i <- integer(k)
      qc_j <- integer(k)
      qc_v <- numeric(k)
      ## Head variable (t): index aux_base + 1 (1-based), coeff = -1
      qc_i[1] <- aux_base + 1L
      qc_j[1] <- aux_base + 1L
      qc_v[1] <- -1.0
      ## Body variables (x): indices aux_base + 2 to aux_base + k, coeff = +1
      for (ii in 2:k) {
        qc_i[ii] <- aux_base + ii
        qc_j[ii] <- aux_base + ii
        qc_v[ii] <- 1.0
      }
      Qc <- Matrix::sparseMatrix(
        i = qc_i, j = qc_j, x = qc_v,
        dims = c(n_total_vars, n_total_vars)
      )
      quadcon_list <- c(quadcon_list, list(list(
        Qc = Qc,
        q = rep(0, n_total_vars),
        rhs = 0,
        sense = "<"
      )))

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

  ## -- Build Gurobi model -----------------------------------------
  model <- list()
  model$modelsense <- "min"

  ## Objective: only original variables contribute
  model$obj <- c(c_vec, rep(0, n_soc_aux))

  ## Combine linear constraints: original + SOC equality
  if (length(soc_eq_A_rows) > 0L) {
    soc_eq_A <- do.call(rbind, soc_eq_A_rows)
    ## Expand A_linear to n_total_vars columns
    if (linear_dim > 0L) {
      A_linear_expanded <- cbind(
        A_linear,
        Matrix::sparseMatrix(i = integer(0), j = integer(0), x = numeric(0),
                              dims = c(linear_dim, n_soc_aux))
      )
    } else {
      A_linear_expanded <- Matrix::sparseMatrix(
        i = integer(0), j = integer(0), x = numeric(0),
        dims = c(0L, n_total_vars)
      )
    }
    model$A <- rbind(A_linear_expanded, soc_eq_A)
    model$rhs <- c(b_linear, soc_eq_rhs)
    model$sense <- c(sense, soc_eq_sense)
  } else {
    if (linear_dim > 0L) {
      model$A <- A_linear
    } else {
      model$A <- Matrix::sparseMatrix(i = integer(0), j = integer(0),
                                       x = numeric(0), dims = c(0L, n_orig))
    }
    model$rhs <- b_linear
    model$sense <- sense
  }

  ## Quadratic constraints
  if (length(quadcon_list) > 0L) {
    model$quadcon <- quadcon_list
  }

  ## -- Quadratic objective ----------------------------------------
  ## ConeMatrixStuffing produces P for 0.5 x'Px + c'x
  ## Gurobi minimizes x'Qx + c'x (NO 1/2 factor)
  ## So Q = P / 2, but we pass lower triangle
  ## R gurobi reads lower triangle of Q for quadratic form
  if (!is.null(data[[SD_P]])) {
    P <- data[[SD_P]]
    Q_half <- P / 2  # 0.5 * P -> Q for Gurobi's x'Qx convention
    ## Use lower triangle (R gurobi reads lower-tri)
    Q_half <- Matrix::tril(Q_half)
    ## Expand to n_total_vars if we have SOC auxiliary vars
    if (n_soc_aux > 0L) {
      Q_expanded <- Matrix::bdiag(
        Q_half,
        Matrix::sparseMatrix(i = integer(0), j = integer(0), x = numeric(0),
                              dims = c(n_soc_aux, n_soc_aux))
      )
      model$Q <- Q_expanded
    } else {
      model$Q <- Q_half
    }
  }

  ## -- Variable bounds --------------------------------------------
  ## CRITICAL: R gurobi defaults lb to 0, NOT -Inf!
  lb <- rep(-Inf, n_total_vars)
  ub <- rep(Inf, n_total_vars)

  ## SOC head variables (t) must be >= 0
  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
    }
  }

  ## -- Variable types for MIP -------------------------------------
  vtype <- rep("C", n_total_vars)

  bool_idx <- data[["bool_idx"]]
  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)
    }
  }

  int_idx <- data[["int_idx"]]
  if (length(int_idx) > 0L) {
    for (idx in int_idx) {
      vtype[idx] <- "I"
    }
  }

  model$lb <- lb
  model$ub <- ub
  model$vtype <- vtype

  ## -- Solver parameters ------------------------------------------
  params <- list(OutputFlag = if (verbose) 1L else 0L)
  params$QCPDual <- 1L  # Enable SOC dual computation

  ## Apply user solver options (skip 'reoptimize' which is our flag)
  skip_opts <- c("reoptimize")
  for (opt_name in setdiff(names(solver_opts), skip_opts)) {
    params[[opt_name]] <- solver_opts[[opt_name]]
  }

  ## -- Warm-start: set initial point from previous solution -------
  ## CVXPY SOURCE: gurobi_conif.py lines 207-215
  cache_key <- GUROBI_SOLVER
  if (warm_start && !is.null(solver_cache) && exists(cache_key, envir = solver_cache)) {
    cached <- get(cache_key, envir = solver_cache)
    cached_status <- GUROBI_STATUS_MAP[[cached$status]]
    if (!is.null(cached_status) && cached_status %in% SOLUTION_PRESENT &&
        !is.null(cached$x) && length(cached$x) == n_total_vars) {
      model$start <- cached$x
    }
  }

  ## -- Call Gurobi ------------------------------------------------
  result <- gurobi::gurobi(model, params)

  ## -- Reoptimize on INF_OR_UNBD ----------------------------------
  if (identical(result$status, "INF_OR_UNBD") &&
      isTRUE(solver_opts[["reoptimize"]])) {
    params$DualReductions <- 0L
    result <- gurobi::gurobi(model, params)
  }

  ## -- Cache for future warm-starts -------------------------------
  ## CVXPY SOURCE: gurobi_conif.py lines 306-307
  if (!is.null(solver_cache)) {
    assign(cache_key, result, envir = solver_cache)
  }

  ## Store metadata for reduction_invert
  result$.n_orig <- n_orig
  result$.n_linear <- linear_dim
  result$.n_soc_eq <- length(soc_eq_rhs)

  result
}

# -- reduction_invert ----------------------------------------------
## Parses Gurobi result format into a CVXR Solution.

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

  ## Map Gurobi string status to CVXR status
  status <- GUROBI_STATUS_MAP[[solution$status]]
  if (is.null(status)) status <- SOLVER_ERROR

  ## Post-adjustment: SOLVER_ERROR with solution -> OPTIMAL_INACCURATE
  if (status == SOLVER_ERROR && !is.null(solution$x)) {
    status <- OPTIMAL_INACCURATE
  }
  ## USER_LIMIT with no solution -> INFEASIBLE_INACCURATE
  if (status == USER_LIMIT && is.null(solution$x)) {
    status <- INFEASIBLE_INACCURATE
  }

  ## Timing and iteration info
  if (!is.null(solution$runtime))
    attr_list[[RK_SOLVE_TIME]] <- solution$runtime
  bar_iter <- solution$baritercount %||% 0L
  simplex_iter <- solution$itercount %||% 0L
  total_iter <- bar_iter + simplex_iter
  if (total_iter > 0L)
    attr_list[[RK_NUM_ITERS]] <- total_iter

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

    ## Primal variables: only original vars (exclude SOC auxiliaries)
    primal_vars <- list()
    primal_vars[[as.character(inverse_data[[SOLVER_VAR_ID]])]] <-
      solution$x[seq_len(n_orig)]

    ## -- Dual variables -------------------------------------------
    ## Skip duals for MIP problems
    is_mip <- isTRUE(inverse_data[["is_mip"]])

    if (!is_mip && !is.null(solution$pi)) {
      ## Gurobi returns pi for linear constraints (eq + ineq + soc_eq)
      ## CVXPY negates all duals: solution["y"] = -np.array(vals)
      ## We do the same: negate all linear constraint duals
      all_pi <- -solution$pi

      ## First n_linear entries are for Zero + NonNeg constraints
      ## Remaining entries are for SOC auxiliary equality constraints
      linear_pi <- if (n_linear > 0L) {
        all_pi[seq_len(n_linear)]
      } else {
        numeric(0)
      }

      dims <- inverse_data[[SD_DIMS]]
      zero_dim <- dims@zero

      ## Equality (zero cone) duals
      eq_dual <- if (zero_dim > 0L) {
        get_dual_values(
          linear_pi[seq_len(zero_dim)],
          extract_dual_value,
          inverse_data[[SOLVER_EQ_CONSTR]]
        )
      } else {
        list()
      }

      ## NonNeg duals
      nonneg_dim <- dims@nonneg
      nonneg_pi <- if (nonneg_dim > 0L) {
        linear_pi[(zero_dim + 1L):(zero_dim + nonneg_dim)]
      } else {
        numeric(0)
      }

      ## SOC duals from quadratic constraint duals (qcpi) +
      ## auxiliary equality constraint duals
      ## For SOC, we need the auxiliary eq duals and the QC duals
      ## The ineq_dual vector should cover NonNeg + SOC constraints
      n_soc_eq <- solution$.n_soc_eq
      soc_eq_pi <- if (n_soc_eq > 0L) {
        all_pi[(n_linear + 1L):(n_linear + n_soc_eq)]
      } else {
        numeric(0)
      }

      ## QC duals (one per SOC constraint)
      qcpi <- if (!is.null(solution$qcpi)) {
        -solution$qcpi  # negate to match convention
      } else {
        numeric(0)
      }

      ## Build ineq_dual vector: [nonneg_pi, soc_duals]
      ## For SOC constraints, the dual is the concatenation of
      ## auxiliary eq duals for that cone + the QC dual
      ## This matches how CVXPY builds ineq_dual from Pi + QCPi:
      ## vals = Pi[ineq_constrs] + Pi[new_leq_constrs] + QCPi[soc_constrs]
      soc_duals <- c(soc_eq_pi, qcpi)
      full_ineq_vec <- c(nonneg_pi, soc_duals)

      ineq_dual <- if (length(full_ineq_vec) > 0L &&
                       length(inverse_data[[SOLVER_NEQ_CONSTR]]) > 0L) {
        get_dual_values(
          full_ineq_vec,
          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)
  }
}

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