R/247_reductions_solvers_conic_solvers_scip_conif.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: reductions/solvers/conic_solvers/scip_conif.py
## SCIP conic solver interface via R scip package
##
## Handles LP, SOCP, MI-LP, MI-SOCP via the conic path.
## Uses the standard ConicSolver pipeline (ConeMatrixStuffing ->
## format_constraints -> ConicSolver.reduction_apply).
##
## Key design (mirrors CVXPY scip_conif.py):
## - Zero + NonNeg constraints -> linear constraints via scip_add_linear_cons
## - SOC constraints -> auxiliary variables + quadratic constraint (x'x <= t^2)
## - MIP via vtype "B"=binary, "I"=integer, "C"=continuous
## - Duals only available for pure LP (no MIP, no SOC)
## - Dual sign: negate all linear duals (matching CVXPY)


# -- SCIP_Solver class ------------------------------------------------
## CVXPY SOURCE: scip_conif.py class SCIP(ConicSolver)

#' @keywords internal
SCIP_Solver <- new_class("SCIP_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, SOC),
      EXP_CONE_ORDER = NULL,
      REQUIRES_CONSTR = FALSE
    )
  }
)

method(solver_name, SCIP_Solver) <- function(x) SCIP_SOLVER

# -- reduction_apply ---------------------------------------------------
## CVXPY SOURCE: scip_conif.py apply() lines 97-135
## Override to add MIP index pass-through and SOC tracking.

method(reduction_apply, SCIP_Solver) <- function(x, problem, ...) {
  data <- problem  ## data list from ConeMatrixStuffing
  inv_data <- list()

  inv_data[[SOLVER_VAR_ID]] <- data[["x_id"]]

  constraints <- data[["constraints"]]
  cone_dims <- data[[SD_DIMS]]
  inv_data[[SD_DIMS]] <- cone_dims

  ## Split into eq (Zero) and ineq (NonNeg, SOC)
  constr_map <- group_constraints(constraints)
  inv_data[[SOLVER_EQ_CONSTR]] <- constr_map[["Zero"]]
  inv_data[[SOLVER_NEQ_CONSTR]] <- c(
    constr_map[["NonNeg"]], constr_map[["SOC"]]
  )

  ## Format constraints: restructure A and b
  formatted <- format_constraints(
    constraints, data[[SD_A]], data[[SD_B]],
    exp_cone_order = x@EXP_CONE_ORDER
  )

  ## Build solver data dict
  ## Convention: solver sees A*x + s = b, s in K
  solver_data <- list()
  solver_data[[SD_C]] <- data[[SD_C]]
  solver_data[[SD_A]] <- -formatted$A
  solver_data[[SD_B]] <- formatted$b
  solver_data[[SD_DIMS]] <- cone_dims

  ## Pass through MIP indices if present
  if (!is.null(data[["bool_idx"]])) solver_data[["bool_idx"]] <- data[["bool_idx"]]
  if (!is.null(data[["int_idx"]]))  solver_data[["int_idx"]]  <- data[["int_idx"]]

  inv_data[[SD_OFFSET]] <- data[[SD_OFFSET]]
  inv_data[["is_mip"]] <- (length(data[["bool_idx"]] %||% integer(0)) > 0L ||
                            length(data[["int_idx"]]  %||% integer(0)) > 0L)

  list(solver_data, inv_data)
}

# -- solve_via_data ----------------------------------------------------
## CVXPY SOURCE: scip_conif.py solve_via_data() lines 171-385
## Converts CVXR solver data (A, b, c, dims) to SCIP model
## using the R scip package's model-building API, then solves.

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

  c_vec <- data[[SD_C]]
  b_vec <- data[[SD_B]]
  A_mat <- data[[SD_A]]
  dims  <- data[[SD_DIMS]]
  n_orig <- length(c_vec)

  ## Convert A to dgTMatrix for efficient (i,j,x) triplet access
  A_T <- as(A_mat, "TsparseMatrix")

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

  ## -- Create SCIP model -------------------------------------------
  model <- scip::scip_model()

  ## -- Create variables --------------------------------------------
  ## CVXPY SOURCE: scip_conif.py _create_variables() lines 202-216
  bool_idx <- data[["bool_idx"]]
  int_idx  <- data[["int_idx"]]
  bool_set <- if (length(bool_idx) > 0L) bool_idx else integer(0)
  int_set  <- if (length(int_idx) > 0L) int_idx else integer(0)

  vtype <- rep("C", n_orig)
  lb <- rep(-Inf, n_orig)
  ub <- rep(Inf, n_orig)

  if (length(bool_set) > 0L) {
    vtype[bool_set] <- "B"
    lb[bool_set] <- 0
    ub[bool_set] <- 1
  }
  if (length(int_set) > 0L) {
    vtype[int_set] <- "I"
  }

  scip::scip_add_vars(model,
    obj   = as.double(c_vec),
    lb    = lb,
    ub    = ub,
    vtype = vtype,
    names = paste0("x", seq_len(n_orig))
  )

  ## -- Add linear constraints (equality + inequality) ---------------
  ## CVXPY SOURCE: scip_conif.py _add_constraints() / add_model_lin_constr()
  ## Build a row-indexed list of (col, coef) pairs from triplet form.
  ## Only need rows up to linear_dim.
  if (linear_dim > 0L) {
    ## Filter triplets to linear rows (0-based i < linear_dim)
    lin_mask <- A_T@i < linear_dim
    tri_i <- A_T@i[lin_mask] + 1L  # 1-based row
    tri_j <- A_T@j[lin_mask] + 1L  # 1-based col (= variable index in scip)
    tri_x <- A_T@x[lin_mask]

    ## Group by row
    row_groups <- split(seq_along(tri_i), tri_i)

    for (row_str in names(row_groups)) {
      row <- as.integer(row_str)
      idx <- row_groups[[row_str]]
      vars  <- tri_j[idx]
      coefs <- tri_x[idx]

      if (row <= zero_dim) {
        ## Equality constraint: sum(coefs * x[vars]) == b[row]
        scip::scip_add_linear_cons(model, vars = vars, coefs = coefs,
                                   lhs = b_vec[row], rhs = b_vec[row])
      } else {
        ## Inequality constraint: sum(coefs * x[vars]) <= b[row]
        scip::scip_add_linear_cons(model, vars = vars, coefs = coefs,
                                   lhs = -Inf, rhs = b_vec[row])
      }
    }

    ## Handle empty rows (all-zero coefficients) — still need the constraint
    all_rows <- seq_len(linear_dim)
    nonempty <- as.integer(names(row_groups))
    empty_rows <- setdiff(all_rows, nonempty)
    for (row in empty_rows) {
      ## 0 == b or 0 <= b — add trivially satisfied/infeasible constraint
      ## Use a dummy variable reference (first var, coef 0)
      if (row <= zero_dim) {
        scip::scip_add_linear_cons(model, vars = 1L, coefs = 0,
                                   lhs = b_vec[row], rhs = b_vec[row])
      } else {
        scip::scip_add_linear_cons(model, vars = 1L, coefs = 0,
                                   lhs = -Inf, rhs = b_vec[row])
      }
    }
  }

  ## -- Add SOC constraints ------------------------------------------
  ## CVXPY SOURCE: scip_conif.py _add_constraints() / add_model_soc_constr()
  ## For each SOC block of size k:
  ##   Create k auxiliary variables (soc_t, soc_x1, ..., soc_{k-1})
  ##   Equality constraint: aux_i = b[i] - A[i,:] * x for each row i
  ##   Quadratic constraint: sum(aux_x_j^2) <= aux_t^2
  soc_dims <- dims@soc
  has_soc <- length(soc_dims) > 0L && sum(soc_dims) > 0L

  if (has_soc) {
    soc_start <- linear_dim  # 0-based row offset

    for (k in soc_dims) {
      soc_rows <- (soc_start + 1L):(soc_start + k)  # 1-based

      ## Add auxiliary variables for this SOC block
      ## First aux var (t) has lb=0; rest are unbounded
      soc_lb <- c(0, rep(-Inf, k - 1L))
      first_aux <- scip::scip_add_vars(model,
        obj   = rep(0, k),
        lb    = soc_lb,
        ub    = rep(Inf, k),
        vtype = rep("C", k),
        names = paste0("soc_", soc_start, "_", seq_len(k))
      )
      aux_indices <- seq(first_aux, first_aux + k - 1L)

      ## For each row in the SOC block: aux_i = b[row] - A[row,:] * x
      ## Equivalently: A[row,:] * x + aux_i = b[row]
      soc_mask <- A_T@i >= soc_start & A_T@i < (soc_start + k)
      soc_tri_i <- A_T@i[soc_mask] - soc_start  # 0-based within block
      soc_tri_j <- A_T@j[soc_mask] + 1L  # 1-based col
      soc_tri_x <- A_T@x[soc_mask]

      soc_row_groups <- split(seq_along(soc_tri_i), soc_tri_i)

      for (local_row in 0:(k - 1L)) {
        row_str <- as.character(local_row)
        aux_var <- aux_indices[local_row + 1L]
        global_row <- soc_start + local_row + 1L  # 1-based into b_vec

        if (row_str %in% names(soc_row_groups)) {
          idx <- soc_row_groups[[row_str]]
          vars  <- c(soc_tri_j[idx], aux_var)
          coefs <- c(soc_tri_x[idx], 1)
        } else {
          vars  <- aux_var
          coefs <- 1
        }
        ## A[row,:]*x + aux = b[row]
        scip::scip_add_linear_cons(model, vars = vars, coefs = coefs,
                                   lhs = b_vec[global_row],
                                   rhs = b_vec[global_row])
      }

      ## Quadratic constraint: sum(aux_x_j^2) <= aux_t^2
      ## Rewrite: sum_{j=2}^{k} aux_j^2 - aux_1^2 <= 0
      ## quadvars1, quadvars2 = same (diagonal), quadcoefs = +1 for body, -1 for head
      qv1 <- aux_indices
      qv2 <- aux_indices
      qc  <- c(-1, rep(1, k - 1L))

      scip::scip_add_quadratic_cons(model,
        quadvars1 = qv1, quadvars2 = qv2, quadcoefs = qc,
        lhs = -Inf, rhs = 0
      )

      soc_start <- soc_start + k
    }
  }

  ## -- Set parameters -----------------------------------------------
  ## CVXPY SOURCE: scip_conif.py _set_params() lines 270-316
  is_mip <- length(bool_set) > 0L || length(int_set) > 0L

  ## For pure LP (no MIP, no SOC), disable presolving/heuristics
  ## to enable dual extraction (matching CVXPY)
  if (!is_mip && !has_soc) {
    scip::scip_set_param(model, "presolving/maxrounds", 0L)
    scip::scip_set_param(model, "misc/usesymmetry", 0L)
  }

  ## Verbosity
  if (!verbose) {
    scip::scip_set_param(model, "display/verblevel", 0L)
  }

  ## User solver options — apply as SCIP native parameters
  ## Support scip_params sub-list (matching CVXPY convention)
  scip_params <- solver_opts[["scip_params"]]
  solver_opts[["scip_params"]] <- NULL

  for (nm in names(solver_opts)) {
    tryCatch(
      scip::scip_set_param(model, nm, solver_opts[[nm]]),
      error = function(e) {
        cli_abort(c(
          "Invalid SCIP solver parameter: {.val {nm}}.",
          "i" = conditionMessage(e)
        ))
      }
    )
  }
  if (!is.null(scip_params)) {
    for (nm in names(scip_params)) {
      tryCatch(
        scip::scip_set_param(model, nm, scip_params[[nm]]),
        error = function(e) {
          cli_abort(c(
            "Invalid SCIP parameter: {.val {nm}}.",
            "i" = conditionMessage(e)
          ))
        }
      )
    }
  }

  ## -- Solve --------------------------------------------------------
  ## CVXPY SOURCE: scip_conif.py _solve() lines 318-385
  tryCatch(
    scip::scip_optimize(model),
    error = function(e) {
      ## Non-fatal: log but continue to check status
    }
  )

  scip_status <- scip::scip_get_status(model)
  info <- scip::scip_get_info(model)
  nsols <- scip::scip_get_nsols(model)

  solution <- list(
    scip_status = scip_status,
    model = model
  )
  solution[[RK_SOLVE_TIME]] <- info$solve_time
  solution[[RK_NUM_ITERS]] <- info$iterations

  ## Extract primal solution if available
  if (nsols > 0L) {
    sol <- scip::scip_get_solution(model)
    solution[["primal"]] <- sol$x[seq_len(n_orig)]  # only original vars

    ## Objective value
    if (scip_status == "timelimit") {
      solution[["value"]] <- NaN
    } else {
      solution[["value"]] <- sol$objval
    }

    ## Dual values (only for pure LP — no MIP, no SOC)
    ## CVXPY SOURCE: scip_conif.py lines 352-369
    ## Note: R scip package currently does NOT expose getDualsolLinear.
    ## Duals are skipped for now (matching CVXPY's known inaccuracy caveat).
  }

  ## Map status
  status <- SCIP_STATUS_MAP[[scip_status]]
  if (is.null(status)) status <- SOLVER_ERROR

  ## Post-adjustment: SOLVER_ERROR with solutions -> OPTIMAL_INACCURATE
  if (status == SOLVER_ERROR && nsols > 0L) {
    status <- OPTIMAL_INACCURATE
  }

  ## Timelimit with no solutions -> SOLVER_ERROR
  if (scip_status == "timelimit" && nsols == 0L) {
    status <- SOLVER_ERROR
  }

  ## Timelimit with solutions -> OPTIMAL_INACCURATE
  if (scip_status == "timelimit" && nsols > 0L) {
    status <- OPTIMAL_INACCURATE
  }

  solution[["status"]] <- status
  solution[[".n_orig"]] <- n_orig

  solution
}

# -- reduction_invert --------------------------------------------------
## CVXPY SOURCE: scip_conif.py invert() lines 137-169

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

  status <- solution[["status"]]
  if (!is.null(solution[[RK_SOLVE_TIME]]))
    attr_list[[RK_SOLVE_TIME]] <- solution[[RK_SOLVE_TIME]]
  if (!is.null(solution[[RK_NUM_ITERS]]))
    attr_list[[RK_NUM_ITERS]] <- solution[[RK_NUM_ITERS]]
  attr_list[[RK_EXTRA_STATS]] <- solution

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

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

    ## Duals: skip for MIP and SOC (matching CVXPY)
    ## R scip package does not currently expose dual extraction API
    dual_vars <- list()

    Solution(status, opt_val, primal_vars, dual_vars, attr_list)
  } else {
    failure_solution(status, attr_list)
  }
}

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

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