R/241_reductions_solvers_conic_solvers_highs_conif.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: reductions/solvers/conic_solvers/highs_conif.py
## HiGHS conic solver interface for LP/MILP problems
##
## HiGHS conic path: for LP and MILP problems (Zero + NonNeg only).
## Uses the standard ConicSolver pipeline (format_constraints -> negate A).
## MIP-capable: supports boolean and integer variables.
##
## QP problems are handled by HiGHS_QP_Solver (qp_solvers/highs_qp_solver.R).


# -- HiGHS status map (shared with HiGHS_QP_Solver) --------------------------
## CVXPY SOURCE: highs_conif.py lines 57-69
## R highs returns integer status codes (unlike Python highspy enum names).

HIGHS_STATUS_MAP <- list(
  "7"  = OPTIMAL,                    # Optimal
  "8"  = INFEASIBLE,                 # Infeasible
  "9"  = INFEASIBLE_OR_UNBOUNDED,    # Unbounded or Infeasible
  "10" = UNBOUNDED,                  # Unbounded
  "11" = USER_LIMIT,                 # Objective Bound
  "12" = USER_LIMIT,                 # Objective Target
  "13" = USER_LIMIT,                 # Time limit
  "14" = USER_LIMIT,                 # Iteration limit
  "15" = USER_LIMIT                  # Solution limit
)
## Codes 0-6, 16 -> SOLVER_ERROR (default fallback)

# -- HiGHS_Conic_Solver class -------------------------------------------------
## CVXPY SOURCE: highs_conif.py class HIGHS(ConicSolver)
## Only supports Zero and NonNeg constraints (LP/MILP only, no SOC).
## MIP_CAPABLE = TRUE (supports MILP).

HiGHS_Conic_Solver <- new_class("HiGHS_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),
      EXP_CONE_ORDER = NULL,
      REQUIRES_CONSTR = FALSE
    )
  }
)

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

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

method(reduction_accepts, HiGHS_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)
  }, logical(1L)))
}

# -- solve_via_data ------------------------------------------------------------
## CVXPY SOURCE: highs_conif.py solve_via_data()
##
## Receives conic data from ConicSolver.apply():
## A (negated), b, c, dims -- format: -A*x + s = b, s in K
##
## For HiGHS, undo the conic negation since HiGHS expects lhs <= Ax <= rhs:
##   Zero rows: A_z*x = b_z -> lhs = rhs = b_z (A is already un-negated for Zero)
##   NonNeg rows: solver has -A_raw*x + s = b, s >= 0 -> A_raw*x <= b
##     -> A_gurobi = solver_data$A (which is -formatted_A = A_raw for Zero,
##       = -A_raw for NonNeg). For HiGHS: negate NonNeg rows back.
##
## Actually simpler: use the raw ConeMatrixStuffing data from the conic base.
## ConicSolver.apply() gives us solver_data$A = -formatted_A and solver_data$B = formatted_b.
## For Zero: formatted_A = -A_raw, so solver_data$A = A_raw, b = -b_raw
## For NonNeg: formatted_A = A_raw, so solver_data$A = -A_raw, b = b_raw
## HiGHS needs lhs <= Ax <= rhs where A = original:
##   Zero: A_raw*x = -b_raw -> but -b_raw != what we want. Let's think differently.
##
## From ConeMatrixStuffing: A_raw*x + b_raw = 0 (Zero), A_raw*x + b_raw >= 0 (NonNeg)
## format_constraints: Zero -> -I block -> formatted_A[Zero] = -A_raw, formatted_b[Zero] = -b_raw
##                     NonNeg -> I block -> formatted_A[NonNeg] = A_raw, formatted_b[NonNeg] = b_raw
## ConicSolver.apply: solver_data$A = -formatted_A, solver_data$B = formatted_b
##   Zero: solver_A = A_raw, solver_b = -b_raw -> A_raw*x + s = -b_raw, s=0 -> A_raw*x = -b_raw
##   NonNeg: solver_A = -A_raw, solver_b = b_raw -> -A_raw*x + s = b_raw, s>=0 -> A_raw*x <= b_raw
## HiGHS: lhs <= Ax <= rhs
##   Zero: A=A_raw, lhs=rhs=-b_raw -> already from solver_A, solver_b -> A=solver_A, lhs=rhs=solver_b
##   NonNeg: A=A_raw=-solver_A, lhs=-Inf, rhs=b_raw=solver_b -> A=-solver_A, rhs=solver_b
##
## So: for the combined matrix, negate NonNeg rows of solver_A to get original A,
##     and use solver_b as the rhs.

method(solve_via_data, HiGHS_Conic_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.")
  }

  c_vec <- data[[SD_C]]
  nvars <- length(c_vec)
  dims <- data[[SD_DIMS]]

  zero_dim <- dims@zero
  nonneg_dim <- dims@nonneg
  total_rows <- zero_dim + nonneg_dim

  A_solver <- data[[SD_A]]  # = -formatted_A
  b_solver <- data[[SD_B]]  # = formatted_b

  ## Build HiGHS constraint matrix and bounds
  ## From conic convention: A_solver * x + s = b_solver, s in K
  ## Zero rows (first zero_dim): A_solver * x + 0 = b_solver -> A*x = b (equality)
  ##   -> A_highs = A_solver, lhs = rhs = b_solver
  ## NonNeg rows (next nonneg_dim): A_solver * x + s = b_solver, s >= 0 -> A*x <= b
  ##   -> A_highs = A_solver, lhs = -Inf, rhs = b_solver

  if (total_rows > 0L) {
    A_highs <- A_solver
    lhs <- b_solver
    rhs <- b_solver

    if (nonneg_dim > 0L) {
      ineq_idx <- (zero_dim + 1L):(zero_dim + nonneg_dim)
      ## A_solver rows are already correct (no negation needed)
      lhs[ineq_idx] <- -Inf
    }
  } else {
    A_highs <- Matrix::sparseMatrix(i = integer(0), j = integer(0),
                                     dims = c(0L, nvars))
    lhs <- numeric(0)
    rhs <- numeric(0)
  }

  ## Ensure A is dgCMatrix
  if (!inherits(A_highs, "dgCMatrix")) {
    A_highs <- methods::as(A_highs, "dgCMatrix")
  }

  ## Q matrix (quadratic objective) -- should be NULL for LP/MILP
  ## MIQP rejection -- HiGHS does not support mixed-integer QP
  if (!is.null(data[[SD_P]])) {
    is_mip <- length(data[["bool_idx"]] %||% integer(0)) > 0L ||
              length(data[["int_idx"]] %||% integer(0)) > 0L
    if (is_mip) {
      cli_abort(c(
        "HiGHS does not support mixed-integer QP (MIQP).",
        "i" = "Use LP constraints with integer variables, or remove integer/boolean attributes for QP problems."
      ))
    }
  }
  Q <- if (!is.null(data[[SD_P]])) {
    methods::as(methods::as(data[[SD_P]], "generalMatrix"), "CsparseMatrix")
  } else {
    NULL
  }

  ## Variable bounds and types
  lower <- rep(-Inf, nvars)
  upper <- rep(Inf, nvars)
  types <- rep(1L, nvars)  # 1=continuous

  ## MIP variable types
  bool_idx <- data[["bool_idx"]]
  if (length(bool_idx) > 0L) {
    for (idx in bool_idx) {
      types[idx] <- 2L
      lower[idx] <- max(lower[idx], 0)
      upper[idx] <- min(upper[idx], 1)
    }
  }

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

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

  for (opt_name in names(solver_opts)) {
    ctrl[[opt_name]] <- solver_opts[[opt_name]]
  }

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

  result
}

# -- reduction_invert ----------------------------------------------------------
## CVXPY SOURCE: highs_conif.py lines 172-190
## Dual sign: negate ALL row_duals (matching CVXPY highs_conif.py line 182)

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

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

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

    ## Check if MIP (duals not valid for MIP)
    is_mip <- isTRUE(inverse_data[["is_mip"]])

    if (!is_mip && !is.null(solution$solver_msg) &&
        isTRUE(solution$solver_msg$dual_valid)) {
      ## Negate ALL row_duals -- matching CVXPY highs_conif.py line 182
      raw_dual <- solution$solver_msg$row_dual
      y <- -raw_dual

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

      ineq_dual <- if (inverse_data[[SD_DIMS]]@nonneg > 0L) {
        zero_dim <- inverse_data[[SD_DIMS]]@zero
        nonneg_dim <- inverse_data[[SD_DIMS]]@nonneg
        get_dual_values(
          y[(zero_dim + 1L):(zero_dim + nonneg_dim)],
          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_Conic_Solver) <- function(x, ...) {
  cat("HiGHS_Conic_Solver()\n")
  invisible(x)
}