R/278_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).


# -- LP-format column-name validation ----------------------------------------
## CVXPY SOURCE: highs_conif.py lines 31-48
## Validates a variable/column name against HiGHS LP-file format rules; used
## when writing a model to an .lp file.  A pure string check -- no dependency
## on the highs package.  Pattern is byte-identical to CVXPY's
## VALID_COLUMN_NAME_PATTERN (Python's `{,254}` written as PCRE `{0,254}`).

VALID_COLUMN_NAME_PATTERN <- paste0(
  "^(?!st$|bounds$|min$|max$|bin$|binary$|gen$|semi$|end$)",
  "[a-df-zA-DF-Z\"!#$%&/}{,;?@_\u2018\u2019'`|~]{1}",
  "[a-zA-Z0-9\"!#$%&/}{,;?@_\u2018\u2019'`|~.=()<>[\\]]{0,254}$"
)

INVALID_COLUMN_NAME_MESSAGE_TEMPLATE <- paste0(
  "Invalid column name: {name}",
  "\nA column name must:",
  "\n- not be equal to one of the keywords: st, bounds, min, max, bin, binary, gen, semi or end",
  "\n- not begin with a number, the letter e or E or any of the following characters: .=()<>[]",
  "\n- be alphanumeric (a-z, A-Z, 0-9) or one of these symbols: \"!#$%&/}{,;?@_\u2018\u2019'`|~.=()<>[]",
  "\n- be no longer than 255 characters."
)

validate_column_name <- function(name) {
  ## CVXPY SOURCE: highs_conif.py validate_column_name (lines 45-48)
  if (!grepl(VALID_COLUMN_NAME_PATTERN, name, perl = TRUE)) {
    msg <- gsub("{name}", name, INVALID_COLUMN_NAME_MESSAGE_TEMPLATE, fixed = TRUE)
    ## Escape literal braces so cli/glue renders them verbatim.
    msg <- gsub("}", "}}", gsub("{", "{{", msg, fixed = TRUE), fixed = TRUE)
    cli::cli_abort(msg)
  }
  invisible(NULL)
}


# -- 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() {
    if (FALSE) new_object(S7_object())  ## S7 static-check guard
    .fast_new(HiGHS_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),
      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_is(c, Zero) || .s7_is(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(), ...) {
  .require_solver_package(HIGHS_SOLVER)

  ## CVXPY SOURCE: highs_conif.py:201 (warm_start, solver_cache plumbed via ...).
  dots <- list(...)
  solver_cache <- dots[["solver_cache"]]

  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 <- data[[LOWER_BOUNDS]] %||% rep(-Inf, nvars)
  upper <- data[[UPPER_BOUNDS]] %||% 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
  ## CVXPY SOURCE: highs_conif.py:296 sets only log_to_console. R-SPECIFIC: the
  ## R `highs` 1.14 persistent-solver API installs a log callback whose console
  ## output is gated by `output_flag` (master switch, default TRUE), so
  ## log_to_console alone no longer silences HiGHS. Mirror CVXPY's intent
  ## (quiet unless verbose) by also setting output_flag.
  ctrl$output_flag <- verbose

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

  ## CVXPY SOURCE: highs_conif.py:250-338.
  ## Persistent-solver flow (highs_model -> hi_new_solver -> optional
  ## hi_solver_set_solution -> hi_solver_run -> getter calls) replaces
  ## the one-shot highs::highs_solve() call.  Necessary so that warm-start
  ## can feed the prior solution into the new solver instance via
  ## hi_solver_set_solution(), mirroring CVXPY's `solver.setSolution()` at
  ## highs_conif.py:320.  Requires highs >= 1.14 (persistent-solver API).
  model <- highs::highs_model(
    Q       = Q,
    L       = c_vec,
    lower   = lower,
    upper   = upper,
    A       = A_highs,
    lhs     = lhs,
    rhs     = rhs,
    types   = types,
    maximum = FALSE,
    offset  = 0
  )
  solver <- highs::hi_new_solver(model)
  highs::hi_solver_set_options(solver, ctrl)

  ## CVXPY SOURCE: highs_conif.py:316-320 (warm-start primal/dual feed-in).
  ## If we have a cached solution from a prior solve and its status was
  ## SOLUTION_PRESENT and the dimensions match, hand it to HiGHS as the
  ## starting point.  Any failure falls through to a cold solve.
  cache_key <- HIGHS_SOLVER
  if (warm_start && !is.null(solver_cache) &&
      exists(cache_key, envir = solver_cache)) {
    cached <- get(cache_key, envir = solver_cache)
    old_status <- HIGHS_STATUS_MAP[[as.character(cached$result$status)]]
    if (!is.null(old_status) && old_status %in% SOLUTION_PRESENT &&
        length(cached$result$solver_msg$col_value) == nvars &&
        length(cached$result$solver_msg$row_value) == nrow(A_highs)) {
      tryCatch({
        prior <- cached$result$solver_msg
        highs::hi_solver_set_solution(
          solver,
          col_value   = prior$col_value,
          row_value   = prior$row_value,
          col_dual    = prior$col_dual,
          row_dual    = prior$row_dual,
          value_valid = isTRUE(prior$value_valid),
          dual_valid  = isTRUE(prior$dual_valid)
        )
      }, error = function(e) NULL)
    }
  }

  ## CVXPY SOURCE: highs_conif.py:323-331 (run + collect result fields).
  ## Result shape matches the prior highs::highs_solve() return so that
  ## reduction_invert() below is unchanged.
  highs::hi_solver_run(solver)
  solution <- highs::hi_solver_get_solution(solver)
  info <- highs::hi_solver_info(solver)
  result <- list(
    primal_solution = solution[["col_value"]],
    objective_value = info[["objective_function_value"]],
    status          = highs::hi_solver_status(solver),
    status_message  = highs::hi_solver_status_message(solver),
    solver_msg      = solution,
    info            = info
  )

  ## CVXPY SOURCE: highs_conif.py:348-349
  ##   if results["model_status"] == "kInfeasible":
  ##       results["dual_ray"] = solver.getDualRay()
  ## The dual (Farkas) ray certifies primal infeasibility; reduction_invert()
  ## maps it onto the constraints' dual_value.  Requires highs >= 1.14 for
  ## hi_solver_get_dual_ray(), which DESCRIPTION already mandates.
  if (identical(HIGHS_STATUS_MAP[[as.character(result$status)]], INFEASIBLE)) {
    result$dual_ray <- tryCatch(highs::hi_solver_get_dual_ray(solver),
                                error = function(e) NULL)
  }

  ## CVXPY SOURCE: highs_conif.py:335-336 (cache for next warm-start).
  if (!is.null(solver_cache)) {
    assign(cache_key,
           list(solver = solver, result = result),
           envir = solver_cache)
  }

  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
  ## DELIBERATE DIVERGENCE, kept on purpose: upstream falls back to `s.UNKNOWN`
  ## (highs_conif.py:170, highs_qpif.py:82), which is NOT a solver status -- it
  ## is the curvature/sign sentinel at settings.py:202 and is in none of
  ## SOLUTION_PRESENT / INF_OR_UNB / INACCURATE / ERROR. A solve that returns it
  ## therefore misses the `s.ERROR` branch of unpack_results (problem.py:1452)
  ## and dies in unpack() with "Cannot unpack invalid solution" instead of a
  ## SolverError naming the solver. Eight reachable HighsModelStatus values hit
  ## that path, including kMemoryLimit and both interrupt statuses.
  ##
  ## Upstream's own warm-start paths (highs_conif.py:334, highs_qpif.py:233)
  ## use SOLVER_ERROR on the same map, as does every other solver interface, so
  ## this reads as a slip rather than a design choice.
  ##
  ## Reported to CVXPY -- draft and full verification log in
  ## notes/upstream_report_cvxpy_highs_unknown.md. Revisit if upstream changes it.
  status <- .status_from_map(HIGHS_STATUS_MAP, status_code)

  ## 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 {
    ## CVXPY SOURCE: highs_conif.py:203-209
    ##   if status == s.INFEASIBLE:
    ##       dual_ray = -np.array(results["dual_ray"][2])
    ##       dual_vars = utilities.get_dual_values(
    ##           dual_ray, utilities.extract_dual_value,
    ##           inverse_data[HIGHS.EQ_CONSTR] + inverse_data[HIGHS.NEQ_CONSTR])
    ##   sol = failure_solution(status, attr, dual_vars)
    ##
    ## R-SPECIFIC guard (the two accessors disagree on the no-ray case):
    ## hi_solver_get_dual_ray() sets has_dual_ray = FALSE and dual_ray = NULL
    ## when HiGHS reports infeasible without producing a ray -- measured for
    ## solver = "ipm" and "pdlp", and for column-bound infeasibility, on BOTH
    ## highs 1.14.0-2 and 1.15.1.  highspy returns a ZERO VECTOR instead, which
    ## is why CVXPY's unguarded [2] does not error there: it silently reports an
    ## all-zero "certificate" satisfying y >= 0 and A'y = 0 but with b'y = 0,
    ## which certifies nothing.  Unguarded in R this would instead be a
    ## zero-length y and a length mismatch in get_dual_values().  Degrade to no
    ## duals -- the pre-existing behavior -- rather than error or fabricate.
    dual_vars <- list()
    ray <- solution$dual_ray
    if (status == INFEASIBLE && !is.null(ray) && isTRUE(ray$has_dual_ray)) {
      y <- -as.numeric(ray$dual_ray)
      zero_dim <- inverse_data[[SD_DIMS]]@zero
      nonneg_dim <- inverse_data[[SD_DIMS]]@nonneg
      if (length(y) == zero_dim + nonneg_dim) {
        eq_dual <- if (zero_dim > 0L) {
          get_dual_values(y[seq_len(zero_dim)], extract_dual_value,
                          inverse_data[[SOLVER_EQ_CONSTR]])
        } else {
          list()
        }
        ineq_dual <- if (nonneg_dim > 0L) {
          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)
      }
    }
    failure_solution(status, attr_list, dual_vars)
  }
}

# -- print ---------------------------------------------------------------------

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