R/246_reductions_solvers_conic_solvers_cvxopt_conif.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: reductions/solvers/conic_solvers/cvxopt_conif.py
## CVXOPT solver interface (backed by R cccp package)
##
## cccp implements the same interior-point algorithm as Python's cvxopt.
## Supports Zero, NonNeg, SOC, PSD cones. No ExpCone, no MIP.
## Convention: A*x + s = b, s in K


# -- CVXOPT status map ---------------------------------------------
## CVXPY SOURCE: cvxopt_conif.py lines 54-64

CVXOPT_STATUS_MAP <- list(
  "optimal"                          = OPTIMAL,
  "feasible"                         = OPTIMAL_INACCURATE,
  "infeasible problem"               = INFEASIBLE,
  "primal infeasible"                = INFEASIBLE,
  "LP relaxation is primal infeasible" = INFEASIBLE,
  "LP relaxation is dual infeasible" = UNBOUNDED,
  "unbounded"                        = UNBOUNDED,
  "dual infeasible"                  = UNBOUNDED,
  "unknown"                          = SOLVER_ERROR,
  "undefined"                        = SOLVER_ERROR,
  "solver_error"                     = SOLVER_ERROR
)

# -- CVXOPT_Solver class -------------------------------------------
## CVXPY SOURCE: cvxopt_conif.py lines 44-52

CVXOPT_Solver <- new_class("CVXOPT_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, SOC, PSD),
      EXP_CONE_ORDER = NULL,
      REQUIRES_CONSTR = FALSE
    )
  }
)

method(solver_name, CVXOPT_Solver) <- function(x) CVXOPT_SOLVER

# -- CVXOPT solve_via_data ----------------------------------------
## CVXPY SOURCE: cvxopt_conif.py lines 164-246
##
## cccp does NOT take a single (G, h, dims) matrix. It takes individual
## cone constraint objects: nnoc, socc, psdc. We split the formatted
## A, b by cone dimensions and construct cccp objects.

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

  dims <- data[[SD_DIMS]]
  A_s  <- data[[SD_A]]
  b_s  <- data[[SD_B]]
  c_obj <- as.numeric(data[[SD_C]])
  nvars <- length(c_obj)

  ## -- Split by cone dimensions ----------------------------------
  zero_dim   <- dims@zero
  nonneg_dim <- dims@nonneg
  soc_dims   <- dims@soc   # integer vector of SOC sizes
  psd_dims   <- dims@psd   # integer vector of PSD sizes (n for n x n)

  ## Equality constraints (Zero cone)
  A_eq <- NULL
  b_eq <- NULL
  if (zero_dim > 0L) {
    A_eq <- as.matrix(A_s[seq_len(zero_dim), , drop = FALSE])
    b_eq <- b_s[seq_len(zero_dim)]
  }

  ## Inequality cone constraints
  offset <- zero_dim
  cList <- list()

  ## NonNeg cone
  if (nonneg_dim > 0L) {
    rows <- offset + seq_len(nonneg_dim)
    G_nn <- as.matrix(A_s[rows, , drop = FALSE])
    h_nn <- b_s[rows]
    cList <- c(cList, list(cccp::nnoc(G = G_nn, h = h_nn)))
    offset <- offset + nonneg_dim
  }

  ## SOC cones
  for (q_i in soc_dims) {
    ## SOC format: first row is t (scalar bound), remaining is body
    ## Convention: A_s * x + s = b_s, s in SOC (s[0] >= ||s[1:]||)
    ## socc(F, g, d, f): ||F*x + g|| <= d'*x + f
    ##   F = -A_body, g = b_body, d = -a_t, f = b_t
    t_row <- offset + 1L
    a_t <- as.numeric(A_s[t_row, ])
    b_t <- b_s[t_row]

    if (q_i > 1L) {
      body_rows <- offset + 1L + seq_len(q_i - 1L)
      A_body <- as.matrix(A_s[body_rows, , drop = FALSE])
      b_body <- b_s[body_rows]
    } else {
      ## Degenerate SOC with no body (shouldn't happen but be safe)
      A_body <- matrix(0, nrow = 0L, ncol = nvars)
      b_body <- numeric(0)
    }

    cList <- c(cList, list(cccp::socc(
      F = -A_body,
      g = b_body,
      d = -a_t,
      f = b_t
    )))
    offset <- offset + q_i
  }

  ## PSD cones
  for (n_psd in psd_dims) {
    ## Identity psd_format_mat: n^2 rows, full vectorized
    ## Convention: mat(b_psd - A_psd * x) is PSD
    ## psdc(Flist, F0): F0 - sum_j x_j F_j is PSD
    n2 <- as.integer(n_psd) * as.integer(n_psd)
    rows <- offset + seq_len(n2)
    A_psd <- A_s[rows, , drop = FALSE]
    b_psd <- b_s[rows]

    F0 <- matrix(b_psd, n_psd, n_psd)
    Flist <- lapply(seq_len(nvars), function(j) {
      matrix(as.numeric(A_psd[, j]), n_psd, n_psd)
    })
    cList <- c(cList, list(cccp::psdc(Flist = Flist, F0 = F0)))
    offset <- offset + n2
  }

  ## -- Build control ---------------------------------------------
  ctrl_args <- list(trace = verbose)
  ## Map standard solver_opts to ctrl() args
  ctrl_params <- c("maxiters", "abstol", "reltol", "feastol", "stepadj", "beta")
  for (nm in names(solver_opts)) {
    if (nm %in% ctrl_params) {
      val <- solver_opts[[nm]]
      if (nm == "maxiters") val <- as.integer(val)
      ctrl_args[[nm]] <- val
    }
  }
  optctrl <- do.call(cccp::ctrl, ctrl_args)

  ## -- Call cccp -------------------------------------------------
  result <- tryCatch(
    cccp::cccp(q = c_obj, A = A_eq, b = b_eq, cList = cList, optctrl = optctrl),
    error = function(e) {
      list(.cvxr_error = TRUE, status = "unknown")
    }
  )

  result
}

# -- CVXOPT reduction_invert ---------------------------------------
## CVXPY SOURCE: cvxopt_conif.py lines 134-137, 211-246
## Map cccp result to CVXR Solution object.

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

  ## Handle solver errors caught by tryCatch
  if (is.list(solution) && isTRUE(solution$.cvxr_error)) {
    return(failure_solution(SOLVER_ERROR))
  }

  status_str <- cccp::getstatus(solution)
  status <- CVXOPT_STATUS_MAP[[status_str]]
  if (is.null(status)) status <- SOLVER_ERROR

  ## Iteration count
  niter <- solution$niter
  if (!is.null(niter)) attr_list[[RK_NUM_ITERS]] <- niter

  if (status %in% SOLUTION_PRESENT) {
    ## Primal values
    primal_vec <- as.numeric(solution$pdv$x)
    state <- cccp::getstate(solution)
    opt_val <- as.numeric(state[["pobj"]]) + inverse_data[[SD_OFFSET]]

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

    ## Equality duals
    y_raw <- as.numeric(solution$pdv$y)
    if (length(y_raw) > 0L && length(inverse_data[[SOLVER_EQ_CONSTR]]) > 0L) {
      eq_dual <- get_dual_values(
        y_raw, extract_dual_value,
        inverse_data[[SOLVER_EQ_CONSTR]]
      )
    } else {
      eq_dual <- list()
    }

    ## Inequality duals (stacked: nonneg, then SOC, then PSD)
    z_raw <- as.numeric(solution$pdv$z)
    if (length(z_raw) > 0L && length(inverse_data[[SOLVER_NEQ_CONSTR]]) > 0L) {
      ineq_dual <- get_dual_values(
        z_raw, extract_dual_value,
        inverse_data[[SOLVER_NEQ_CONSTR]]
      )
    } else {
      ineq_dual <- list()
    }

    dual_vars <- c(eq_dual, ineq_dual)
    Solution(status, opt_val, primal_vars, dual_vars, attr_list)
  } else {
    failure_solution(status, attr_list)
  }
}

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