R/236_reductions_solvers_conic_solvers_conic_solver.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: reductions/solvers/conic_solvers/conic_solver.py
## ConicSolver -- conic solver base with format_constraints
##
## Non-parametric simplification: works directly with A, b matrices
## rather than parameter tensors.


# -- ConicSolver class --------------------------------------------
## CVXPY SOURCE: conic_solver.py lines 100-401

ConicSolver <- new_class("ConicSolver", parent = Solver, package = "CVXR",
  properties = list(
    SUPPORTED_CONSTRAINTS = class_list,
    EXP_CONE_ORDER = class_any,
    REQUIRES_CONSTR = class_logical
  ),
  constructor = function(supported = list(Zero, NonNeg),
                         exp_cone_order = NULL,
                         requires_constr = FALSE,
                         MIP_CAPABLE = FALSE) {
    new_object(S7_object(),
      .cache = new.env(parent = emptyenv()),
      MIP_CAPABLE = MIP_CAPABLE,
      BOUNDED_VARIABLES = FALSE,
      SUPPORTED_CONSTRAINTS = supported,
      EXP_CONE_ORDER = exp_cone_order,
      REQUIRES_CONSTR = requires_constr
    )
  }
)

# -- get_spacing_matrix --------------------------------------------
## CVXPY SOURCE: conic_solver.py lines 132-160
## Static helper that spaces out rows with interleaving.
##
## Returns a sparse matrix of shape (shape[1], shape[2]) with ones
## placed at regular intervals: num_blocks blocks of `streak` consecutive
## ones, separated by `spacing` zero rows, starting at row `offset`.
##
## All indices are 0-based internally, converted to 1-based for sparseMatrix.

get_spacing_matrix <- function(shape, spacing, streak, num_blocks, offset) {
  num_values <- num_blocks * streak
  streak_plus_spacing <- streak + spacing

  ## Row indices: for each block, take the first `streak` rows of a

  ## (streak + spacing)-sized chunk, then shift by offset
  block_starts <- seq(0L, by = as.integer(streak_plus_spacing),
                      length.out = num_blocks)
  row_arr <- integer(num_values)
  k <- 1L
  for (b in seq_len(num_blocks)) {
    start <- block_starts[b] + offset
    for (s in seq_len(streak)) {
      row_arr[k] <- start + s - 1L  # 0-based
      k <- k + 1L
    }
  }

  col_arr <- seq(0L, length.out = num_values)  # 0-based

  Matrix::sparseMatrix(
    i = row_arr + 1L, j = col_arr + 1L,
    x = rep(1.0, num_values),
    dims = as.integer(shape)
  )
}

# -- psd_format_mat (default: identity) ----------------------------
## CVXPY SOURCE: conic_solver.py lines 162-167
## Subclasses (SCS, Clarabel) override with triangle extraction + scaling.

psd_format_mat <- function(constr) {
  Matrix::Diagonal(constr_size(constr))
}

# -- format_constraints --------------------------------------------
## CVXPY SOURCE: conic_solver.py lines 169-320
## Builds a block-diagonal restructuring matrix and applies to A and b.
##
## Non-parametric simplification: directly multiply sparse matrices
## instead of using LinearOperator abstractions.
##
## Returns list(A = formatted_A, b = formatted_b)

format_constraints <- function(constraints, A, b, exp_cone_order,
                               psd_format_fn = psd_format_mat) {
  if (length(constraints) == 0L) {
    return(list(A = A, b = b))
  }

  restruct_blocks <- vector("list", length(constraints))

  for (ci in seq_along(constraints)) {
    constr <- constraints[[ci]]
    total_height <- sum(vapply(constr@args, expr_size, integer(1L)))

    if (S7_inherits(constr, Zero)) {
      ## Negate: Ax + b = 0 -> -Ax - b = 0 (solver: Ax + s = b, s = 0)
      n <- constr_size(constr)
      restruct_blocks[[ci]] <- -Matrix::Diagonal(n)

    } else if (S7_inherits(constr, NonNeg)) {
      ## Identity: Ax + b >= 0
      n <- constr_size(constr)
      restruct_blocks[[ci]] <- Matrix::Diagonal(n)

    } else if (S7_inherits(constr, SOC)) {
      ## Interleave t rows and X rows for each cone
      ## SOC axis must be 0 (lowered by ConeMatrixStuffing)
      x_dim <- constr@args[[2L]]@shape[1L]

      t_spacer <- get_spacing_matrix(
        shape = c(total_height, expr_size(constr@args[[1L]])),
        spacing = x_dim,
        streak = 1L,
        num_blocks = expr_size(constr@args[[1L]]),
        offset = 0L
      )
      X_spacer <- get_spacing_matrix(
        shape = c(total_height, expr_size(constr@args[[2L]])),
        spacing = 1L,
        streak = x_dim,
        num_blocks = expr_size(constr@args[[1L]]),
        offset = 1L
      )
      restruct_blocks[[ci]] <- cbind(t_spacer, X_spacer)

    } else if (S7_inherits(constr, ExpCone)) {
      ## 3-way interleave of x, y, z args
      n_eargs <- length(constr@args)
      arg_mats <- vector("list", n_eargs)
      for (i in seq_len(n_eargs)) {
        arg <- constr@args[[i]]
        arg_mats[[i]] <- get_spacing_matrix(
          shape = c(total_height, expr_size(arg)),
          spacing = length(exp_cone_order) - 1L,
          streak = 1L,
          num_blocks = expr_size(arg),
          offset = exp_cone_order[i]
        )
      }
      restruct_blocks[[ci]] <- do.call(cbind, arg_mats)

    } else if (S7_inherits(constr, PowCone3D)) {
      ## 3-way interleave of x, y, z args (like ExpCone but offset = i-1)
      n_pargs <- length(constr@args)
      arg_mats <- vector("list", n_pargs)
      for (i in seq_len(n_pargs)) {
        arg <- constr@args[[i]]
        arg_mats[[i]] <- get_spacing_matrix(
          shape = c(total_height, expr_size(arg)),
          spacing = 2L,
          streak = 1L,
          num_blocks = expr_size(arg),
          offset = i - 1L
        )
      }
      restruct_blocks[[ci]] <- do.call(cbind, arg_mats)

    } else if (S7_inherits(constr, PowConeND)) {
      ## CVXPY SOURCE: conic_solver.py lines 255-277
      w_arg <- constr@args[[1L]]
      if (length(w_arg@shape) <= 1L) {
        m <- w_arg@shape[1L]
        n <- 1L
      } else {
        m <- w_arg@shape[1L]
        n <- w_arg@shape[2L]
      }
      arg_mats <- vector("list", n + 1L)
      for (j in seq_len(n)) {
        arg_mats[[j]] <- get_spacing_matrix(
          shape = c(total_height, m),
          spacing = 0L,
          streak = 1L,
          num_blocks = m,
          offset = (m + 1L) * (j - 1L)
        )
      }
      ## Hypo columns
      arg_mats[[n + 1L]] <- get_spacing_matrix(
        shape = c(total_height, n),
        spacing = m,
        streak = 1L,
        num_blocks = n,
        offset = m
      )
      restruct_blocks[[ci]] <- do.call(cbind, arg_mats)

    } else if (S7_inherits(constr, PSD)) {
      ## Solver-specific PSD format matrix
      restruct_blocks[[ci]] <- psd_format_fn(constr)

    } else {
      cli_abort("Unsupported constraint type: {.cls {short_class_name(constr)}}.")
    }
  }

  ## Build block-diagonal restructuring matrix
  restruct_mat <- Matrix::bdiag(restruct_blocks)

  ## Apply restructuring
  formatted_A <- restruct_mat %*% A
  formatted_b <- as.numeric(restruct_mat %*% b)

  list(A = formatted_A, b = formatted_b)
}

# -- reduction_accepts ---------------------------------------------
## CVXPY SOURCE: conic_solver.py lines 124-130
## Checks if solver supports all constraint types in the data.

method(reduction_accepts, ConicSolver) <- function(x, problem, ...) {
  ## In our non-parametric path, "problem" is data from ConeMatrixStuffing
  if (!is.list(problem) || is.null(problem[["constraints"]])) return(FALSE)
  constrs <- problem[["constraints"]]
  if (length(constrs) == 0L && x@REQUIRES_CONSTR) return(FALSE)
  supported <- x@SUPPORTED_CONSTRAINTS
  all(vapply(constrs, function(c) {
    .inherits_any(c, supported)
  }, logical(1L)))
}

# -- reduction_apply (non-parametric) ------------------------------
## CVXPY SOURCE: conic_solver.py lines 344-401
## Non-parametric: receives data dict from ConeMatrixStuffing,
## restructures A and b, returns solver-ready data dict.

method(reduction_apply, ConicSolver) <- function(x, problem, ...) {
  data <- problem  ## "problem" is actually the data list from ConeMatrixStuffing
  inv_data <- list()

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

  ## Get constraints and cone dims
  constraints <- data[["constraints"]]
  cone_dims <- data[[SD_DIMS]]
  inv_data[[SD_DIMS]] <- cone_dims

  ## Split into eq (Zero) and ineq (all others)
  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"]],
    constr_map[["PSD"]], constr_map[["ExpCone"]],
    constr_map[["PowCone3D"]], constr_map[["PowConeND"]]
  )

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

  ## Build solver data dict
  ## Convention: solver sees A*x + s = b, s in K
  ## formatted_A encodes: for Zero: -A, for NonNeg: A, etc.
  ## We negate to get solver convention: data[A] = -formatted_A
  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 P matrix through for QP path
  if (!is.null(data[[SD_P]])) solver_data[[SD_P]] <- data[[SD_P]]

  ## Pass through MIP indices for MIP-capable conic solvers
  solver_data[["bool_idx"]] <- data[["bool_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)
}

# -- solver_psd_format_mat -----------------------------------------
## S7 generic for subclass-specific PSD format matrix.
## Default: identity (SCS overrides with lower-tri, Clarabel with upper-tri).

solver_psd_format_mat <- new_generic("solver_psd_format_mat", "solver",
  function(solver, constr) S7_dispatch())

method(solver_psd_format_mat, ConicSolver) <- function(solver, constr) {
  psd_format_mat(constr)
}

# -- reduction_invert ----------------------------------------------
## CVXPY SOURCE: conic_solver.py lines 322-342
## Builds Solution from solver result dictionary.
##
## This default implementation handles the standard result format
## where the solver returns a list with status, value, primal, eq_dual,
## ineq_dual keys. Subclasses override for solver-specific formats.

method(reduction_invert, ConicSolver) <- function(x, solution, inverse_data, ...) {
  status <- solution[[RK_STATUS]]

  if (status %in% SOLUTION_PRESENT) {
    opt_val <- solution[[RK_VALUE]]
    primal_vars <- list()
    primal_vars[[as.character(inverse_data[[SOLVER_VAR_ID]])]] <-
      solution[[RK_PRIMAL]]

    eq_dual <- get_dual_values(
      solution[[RK_EQ_DUAL]],
      extract_dual_value,
      inverse_data[[SOLVER_EQ_CONSTR]]
    )
    ineq_dual <- get_dual_values(
      solution[[RK_INEQ_DUAL]],
      extract_dual_value,
      inverse_data[[SOLVER_NEQ_CONSTR]]
    )
    dual_vars <- c(eq_dual, ineq_dual)
    Solution(status, opt_val, primal_vars, dual_vars, list())
  } else {
    failure_solution(status)
  }
}

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

method(print, ConicSolver) <- function(x, ...) {
  names <- vapply(x@SUPPORTED_CONSTRAINTS, function(cls) cls@name, character(1L))
  cat(sprintf("ConicSolver(supported=%s)\n",
    paste(names, collapse = ", ")))
  invisible(x)
}