R/248_reductions_solvers_qp_solvers_highs_qpif.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: reductions/solvers/qp_solvers/highs_qpif.py
## HiGHS QP solver interface for LP/QP problems
##
## HiGHS solves: minimize 0.5 x'Qx + L'x  s.t. lhs <= Ax <= rhs,
##               lower <= x <= upper
## Accepts ONLY Zero (equality) and NonNeg (inequality) constraints.
## Inherits from QpSolver -- uses QpSolver.apply() for sign-correct data.
##
## Key differences from OSQP:
##   - Q matrix: full symmetric dgCMatrix (NOT upper-triangle)
##   - Dual variables: NEGATE ALL row_duals (matching CVXPY highs_qpif.py line 97)
##   - Status codes: integer codes from R highs (not string enum names)
##   - NOT MIP capable in QP path (MIQP not supported -- CVXPY highs_qpif.py line 37)


# -- HiGHS status map ---------------------------------------------------------
## Reuse HIGHS_STATUS_MAP from highs_conic_solver.R (defined there to avoid
## duplicate, since both files are loaded)

# -- HiGHS_QP_Solver class ----------------------------------------------------
## CVXPY SOURCE: highs_qpif.py lines 30-130

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

method(solver_name, HiGHS_QP_Solver) <- function(x) HIGHS_SOLVER

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

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

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

  ## Q matrix (quadratic objective) -- full symmetric dgCMatrix for HiGHS
  ## DIFFERS from OSQP (upper-triangle) and Clarabel/SCS (dsCMatrix)
  if (!is.null(data[[SD_P]])) {
    Q <- methods::as(methods::as(data[[SD_P]], "generalMatrix"), "CsparseMatrix")
  } else {
    Q <- NULL
  }

  ## Stack A_eq and F_ineq into combined HiGHS constraint matrix
  ## Bounds: lhs = [b_eq, -inf*ones], rhs = [b_eq, g_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)

  if (len_eq > 0L && len_ineq > 0L) {
    A <- rbind(A_eq, F_ineq)
    lhs <- c(b_eq, rep(-Inf, len_ineq))
    rhs <- c(b_eq, g_ineq)
  } else if (len_eq > 0L) {
    A <- A_eq
    lhs <- b_eq
    rhs <- b_eq
  } else if (len_ineq > 0L) {
    A <- F_ineq
    lhs <- rep(-Inf, len_ineq)
    rhs <- g_ineq
  } else {
    A <- Matrix::sparseMatrix(i = integer(0), j = integer(0),
                               dims = c(0L, nvars))
    lhs <- numeric(0)
    rhs <- numeric(0)
  }

  ## Ensure A is dgCMatrix for highs
  if (!is.null(A) && !inherits(A, "dgCMatrix")) {
    A <- methods::as(A, "dgCMatrix")
  }

  ## Variable bounds: default to (-Inf, Inf) -- no MIP handling in QP path
  lower <- rep(-Inf, nvars)
  upper <- rep(Inf, nvars)

  ## All continuous (QP path is NOT MIP-capable)
  types <- rep(1L, nvars)

  ## Build HiGHS control
  ctrl <- highs::highs_control()
  ctrl$log_to_console <- verbose

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

  ## Call HiGHS
  result <- highs::highs_solve(
    Q       = Q,
    L       = L_vec,
    lower   = lower,
    upper   = upper,
    A       = A,
    lhs     = lhs,
    rhs     = rhs,
    types   = types,
    maximum = FALSE,
    offset  = 0,
    control = ctrl
  )

  ## Store len_eq for dual splitting
  result$.len_eq <- len_eq
  result
}

# -- reduction_invert ----------------------------------------------------------
## CVXPY SOURCE: highs_qpif.py lines 76-120
## Dual sign: negate ALL row_duals -- matching CVXPY highs_qpif.py line 97.

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

  ## Map status via integer status code
  status_code <- solution$status
  status <- HIGHS_STATUS_MAP[[as.character(status_code)]]
  if (is.null(status)) status <- SOLVER_ERROR

  ## Timing and iteration info
  info <- solution$info
  if (!is.null(info)) {
    num_iters <- (info$simplex_iteration_count %||% 0L) +
                 (info$ipm_iteration_count %||% 0L) +
                 (info$qp_iteration_count %||% 0L) +
                 (info$crossover_iteration_count %||% 0L)
    if (num_iters > 0L) attr_list[[RK_NUM_ITERS]] <- num_iters
  }

  if (status %in% SOLUTION_PRESENT) {
    ## Objective value
    opt_val <- solution$objective_value + inverse_data[[SD_OFFSET]]

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

    ## Dual variables: negate ALL row_duals
    ## CVXPY SOURCE: highs_qpif.py line 97: y = -np.array(results["solution"].row_dual)
    if (!is.null(solution$solver_msg) &&
        isTRUE(solution$solver_msg$dual_valid)) {
      raw_dual <- solution$solver_msg$row_dual
      y <- -raw_dual  # negate ALL
      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, HiGHS_QP_Solver) <- function(x, ...) {
  cat("HiGHS_QP_Solver()\n")
  invisible(x)
}