R/242_reductions_solvers_conic_solvers_glpk_conif.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: reductions/solvers/conic_solvers/glpk_conif.py
## GLPK conic solver interface for LP problems (continuous only).
##
## GLPK supports Zero + NonNeg cones (LP only, no SOC/ExpCone/PSD).
## Uses the standard ConicSolver pipeline (format_constraints -> negate A).
## NOT MIP-capable; see glpk_mi_conif.R for MILP support.
##
## R interface: Rglpk::Rglpk_solve_LP()
## CRITICAL: Rglpk defaults to bounds [0, Inf) -- MUST override to (-Inf, Inf).


# -- GLPK status map (native codes from glpk.h) --------------------------------
## With canonicalize_status = FALSE, Rglpk returns GLPK's native status codes:
##   1 = GLP_UNDEF  (solution undefined)
##   2 = GLP_FEAS   (feasible but not optimal)
##   3 = GLP_INFEAS (integer infeasible -- MIP only)
##   4 = GLP_NOFEAS (no feasible solution exists)
##   5 = GLP_OPT    (optimal solution found)
##   6 = GLP_UNBND  (problem is unbounded)

GLPK_STATUS_MAP <- list(
  "5" = OPTIMAL,                    # GLP_OPT
  "2" = OPTIMAL_INACCURATE,         # GLP_FEAS (feasible, not proven optimal)
  "4" = INFEASIBLE,                 # GLP_NOFEAS
  "3" = INFEASIBLE,                 # GLP_INFEAS (MIP infeasible)
  "6" = UNBOUNDED,                  # GLP_UNBND
  "1" = SOLVER_ERROR                # GLP_UNDEF
)

# -- Helper: dgCMatrix -> slam simple_triplet_matrix ----------------------------
## Rglpk requires slam::simple_triplet_matrix format for sparse matrices.

.dgC_to_stm <- function(x) {
  if (!inherits(x, "dgCMatrix")) {
    x <- methods::as(x, "dgCMatrix")
  }
  ## dgCMatrix is CSC format: p (col pointers), i (row indices 0-based), x (values)
  dims <- dim(x)
  ## Convert to triplet form
  cols <- rep(seq_len(dims[2L]), diff(x@p))
  rows <- x@i + 1L  # 0-based to 1-based
  slam::simple_triplet_matrix(
    i = rows, j = cols, v = x@x,
    nrow = dims[1L], ncol = dims[2L]
  )
}

# -- GLPK_Solver class ----------------------------------------------------------
## CVXPY SOURCE: glpk_conif.py class GLPK(CVXOPT)
## LP only (MIP_CAPABLE = FALSE). Supports Zero + NonNeg constraints.
## In CVXPY, GLPK inherits from CVXOPT; here we inherit directly from
## ConicSolver since we use Rglpk (not CVXOPT).

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

method(solver_name, GLPK_Solver) <- function(x) GLPK_SOLVER

# -- reduction_accepts -----------------------------------------------------------
## Override: GLPK only handles Zero + NonNeg (LP problems).

method(reduction_accepts, GLPK_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: glpk_conif.py solve_via_data() (adapted for Rglpk)
##
## Conic data from ConicSolver.apply():
##   A_solver * x + s = b_solver, s in K
##   Zero rows: s = 0 -> A_solver * x = b_solver (equality)
##   NonNeg rows: s >= 0 -> A_solver * x <= b_solver (inequality)
##
## Maps to Rglpk_solve_LP(obj, mat, dir, rhs, bounds, types, control).
## CRITICAL: bounds must be (-Inf, Inf) since conic pipeline handles all bounds.

method(solve_via_data, GLPK_Solver) <- function(x, data, warm_start = FALSE, verbose = FALSE,
                                                  solver_opts = list(), ...) {
  if (!requireNamespace("Rglpk", quietly = TRUE)) {
    cli_abort("Package {.pkg Rglpk} 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 constraint matrix and direction vector
  if (total_rows > 0L) {
    ## Convert sparse matrix to slam format
    mat <- .dgC_to_stm(A_solver)

    ## Direction: Zero rows -> "==", NonNeg rows -> "<="
    dir <- character(total_rows)
    if (zero_dim > 0L) {
      dir[seq_len(zero_dim)] <- "=="
    }
    if (nonneg_dim > 0L) {
      dir[(zero_dim + 1L):(zero_dim + nonneg_dim)] <- "<="
    }

    rhs <- b_solver
  } else {
    ## No constraints -- create a trivial constraint to avoid Rglpk error
    mat <- slam::simple_triplet_matrix(
      i = integer(0), j = integer(0), v = numeric(0),
      nrow = 0L, ncol = nvars
    )
    dir <- character(0)
    rhs <- numeric(0)
  }

  ## Variable bounds: MUST be (-Inf, Inf) -- conic pipeline handles bounds
  bounds <- list(
    lower = list(ind = seq_len(nvars), val = rep(-Inf, nvars)),
    upper = list(ind = seq_len(nvars), val = rep(Inf, nvars))
  )

  ## Variable types (all continuous for LP)
  types <- NULL  # NULL = all continuous in Rglpk

  ## Build control list
  ctrl <- list(canonicalize_status = FALSE, verbose = verbose)
  for (opt_name in names(solver_opts)) {
    ctrl[[opt_name]] <- solver_opts[[opt_name]]
  }

  ## Call Rglpk
  result <- Rglpk::Rglpk_solve_LP(
    obj     = c_vec,
    mat     = mat,
    dir     = dir,
    rhs     = rhs,
    bounds  = bounds,
    types   = types,
    max     = FALSE,
    control = ctrl
  )

  result
}

# -- reduction_invert ------------------------------------------------------------
## Parse Rglpk result, extract primal + dual values.
## Dual sign: Rglpk returns standard LP duals (pi), but conic convention
## uses y = -pi. Negate ALL duals to match Clarabel/SCS convention.

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

  status_code <- solution$status
  status <- GLPK_STATUS_MAP[[as.character(status_code)]]
  if (is.null(status)) status <- SOLVER_ERROR

  if (status %in% SOLUTION_PRESENT) {
    opt_val <- solution$optimum + inverse_data[[SD_OFFSET]]

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

    ## Extract dual values from auxiliary$dual
    ## Negate: Rglpk returns LP duals (pi); conic convention needs y = -pi
    raw_dual <- solution$auxiliary$dual
    if (!is.null(raw_dual) && !anyNA(raw_dual)) {
      raw_dual <- -raw_dual
      eq_dual <- if (inverse_data[[SD_DIMS]]@zero > 0L) {
        zero_dim <- inverse_data[[SD_DIMS]]@zero
        get_dual_values(
          raw_dual[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(
          raw_dual[(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, GLPK_Solver) <- function(x, ...) {
  cat("GLPK_Solver()\n")
  invisible(x)
}