R/249_reductions_solvers_qp_solvers_gurobi_qpif.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: reductions/solvers/qp_solvers/gurobi_qpif.py
## Gurobi QP solver interface for QP/MIQP problems
##
## Uses QpSolver.apply() for sign-correct A_eq/b_eq/F_ineq/g_ineq data.
## The ONLY QP solver that is MIP_CAPABLE (supports MIQP).
##
## Key differences from Gurobi_Conic_Solver:
##   - No SOC handling (quadratic constraints)
##   - Uses QP path (simpler, no auxiliary variables)
##   - Dual sign: negate ALL constraint duals (matching CVXPY gurobi_qpif.py line 121)


# -- Gurobi_QP_Solver class ---------------------------------------------------
## CVXPY SOURCE: gurobi_qpif.py lines 30-175

Gurobi_QP_Solver <- new_class("Gurobi_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, Gurobi_QP_Solver) <- function(x) GUROBI_SOLVER

# -- solve_via_data ------------------------------------------------------------
## CVXPY SOURCE: gurobi_qpif.py lines 140-175
## Receives QP data from QpSolver.apply(): P, q, A_eq, b_eq, F_ineq, g_ineq

method(solve_via_data, Gurobi_QP_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"]]

  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) {
      A_model <- rbind(A_eq, F_ineq)
      rhs_model <- c(b_eq, g_ineq)
      sense_model <- c(rep("=", len_eq), rep("<", len_ineq))
    } else if (len_eq > 0L) {
      A_model <- A_eq
      rhs_model <- b_eq
      sense_model <- rep("=", len_eq)
    } else {
      A_model <- F_ineq
      rhs_model <- g_ineq
      sense_model <- rep("<", len_ineq)
    }
  } else {
    A_model <- Matrix::sparseMatrix(i = integer(0), j = integer(0),
                                     x = numeric(0), dims = c(0L, nvars))
    rhs_model <- numeric(0)
    sense_model <- character(0)
  }

  ## Build Gurobi model
  model <- list()
  model$modelsense <- "min"
  model$obj <- q_vec
  model$A <- A_model
  model$rhs <- rhs_model
  model$sense <- sense_model

  ## Quadratic objective: P for 0.5 x'Px + q'x -> Q = P/2 for Gurobi's x'Qx + c'x
  ## R gurobi reads lower triangle
  if (!is.null(data[[SD_P]])) {
    P <- data[[SD_P]]
    Q_half <- P / 2
    model$Q <- Matrix::tril(Q_half)
  }

  ## Variable bounds
  lb <- rep(-Inf, nvars)
  ub <- rep(Inf, nvars)

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

  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)

  ## Apply user solver options (skip 'reoptimize')
  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_qpif.py lines 186-198
  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) == nvars) {
      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_qpif.py lines 239-240
  if (!is.null(solver_cache)) {
    assign(cache_key, result, envir = solver_cache)
  }

  ## Store metadata
  result$.len_eq <- len_eq
  result$.n_constrs <- n_constrs
  result$.is_mip <- length(bool_idx %||% integer(0)) > 0L ||
    length(int_idx %||% integer(0)) > 0L

  result
}

# -- reduction_invert ----------------------------------------------------------
## CVXPY SOURCE: gurobi_qpif.py lines 76-135
## Dual sign: negate ALL pi duals -- matching CVXPY gurobi_qpif.py line 121.

method(reduction_invert, Gurobi_QP_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
  if (status == SOLVER_ERROR && !is.null(solution$x)) {
    status <- OPTIMAL_INACCURATE
  }
  if (status == USER_LIMIT && is.null(solution$x)) {
    status <- INFEASIBLE_INACCURATE
  }

  ## Timing
  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) {
    opt_val <- solution$objval + inverse_data[[SD_OFFSET]]

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

    ## Dual variables: negate ALL pi duals
    ## CVXPY: y = -np.array([constraints_grb[i].Pi for i in range(m)])
    is_mip <- isTRUE(solution$.is_mip)

    if (!is_mip && !is.null(solution$pi)) {
      y <- -solution$pi
      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, Gurobi_QP_Solver) <- function(x, ...) {
  cat("Gurobi_QP_Solver()\n")
  invisible(x)
}