R/271_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) {
    if (FALSE) new_object(S7_object())  ## S7 static-check guard
    .fast_new(ConicSolver, S7_object(),
      .cache = new.env(parent = emptyenv()),
      MIP_CAPABLE = MIP_CAPABLE,
      BOUNDED_VARIABLES = FALSE,
      PSD_TRIANGLE_KIND = NA_character_,
      PSD_SQRT2_SCALING = NA,
      SUPPORTED_CONSTRAINTS = supported,
      EXP_CONE_ORDER = exp_cone_order,
      REQUIRES_CONSTR = requires_constr
    )
  }
)

## CVXPY v1.8.2: supports_quad_obj — default FALSE for conic solvers.
## Solvers that can handle a quadratic objective (P matrix) in the conic
## path override this to TRUE (Clarabel, SCS >= 3.0).
method(supports_quad_obj, ConicSolver) <- function(x) FALSE

# -- 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)
  )
}

# -- format_constraints --------------------------------------------
## CVXPY SOURCE: conic_solver.py lines 169-320
## Builds a block-diagonal restructuring matrix and applies to A and b.
##
## Returns list(A = formatted_A, b = formatted_b, param_prog = updated_pp).
##
## If `param_prog` is non-NULL and its `formatted` slot is FALSE, the
## same row permutation is also applied to `param_prog@A_tensor` (via
## the I_{x_length+1} ⊗ restruct_mat Kronecker pattern that CVXPY uses
## at conic_solver.py:290-307), and the returned `param_prog` is marked
## `formatted = TRUE`.  If `param_prog@formatted` is already TRUE, the
## function is a no-op on both the numeric A,b and the param_prog --
## this is what makes the DPP fast-path re-solve safe against double
## formatting.

format_constraints <- function(constraints, A, b, exp_cone_order,
                               param_prog    = NULL) {
  ## Fast-path / DPP re-solve: if the caller's param_prog has already
  ## been formatted on a previous solve, A and b coming in are already
  ## restructured -- pass everything through unchanged.
  if (!is.null(param_prog) && isTRUE(param_prog@formatted)) {
    return(list(A = A, b = b, param_prog = param_prog))
  }
  if (length(constraints) == 0L) {
    return(list(A = A, b = b, param_prog = param_prog))
  }

  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_is(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_is(constr, NonNeg)) {
      ## Identity: Ax + b >= 0
      n <- constr_size(constr)
      restruct_blocks[[ci]] <- Matrix::Diagonal(n)

    } else if (.s7_is(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_is(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_is(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_is(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_is(constr, PSD) || .s7_is(constr, SvecPSD)) {
      ## CVXPY SOURCE: conic_solver.py:267-270 -- identity for BOTH.  A PSD
      ## constraint reaching a solver unconverted means that solver takes full
      ## matrices (CVXOPT); an SvecPSD has already been packed by
      ## `PSDToSvecPSD` (reductions/cone2cone/exact.R).  Neither is repacked
      ## here -- that is what the svec packing moving into the reduction bought.
      restruct_blocks[[ci]] <- Matrix::Diagonal(constr_size(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 to the numeric A and b (existing behaviour).
  formatted_A <- restruct_mat %*% A
  formatted_b <- as.numeric(restruct_mat %*% b)

  ## CVXPY SOURCE: conic_solver.py:290-307 — apply the SAME row
  ## permutation to the parametric A_tensor so the cached tensor agrees
  ## with the solver's cone-interleaved row ordering.  Pre-fix, the
  ## numeric A,b were permuted but A_tensor was left in pre-format
  ## (block) layout, so `apply_param_jac` in the backward pass produced
  ## gradients in a coordinate system different from the one diffcp's
  ## DT returned -- contributions leaked between parameters when one
  ## parameter appeared in both log-transformed (alpha^2) and
  ## non-log-transformed (Constant^alpha) slots of the same DGP problem.
  ##
  ## The mapping in F-order flat coordinates is
  ##   I_{x_length + 1} ⊗ restruct_mat
  ## applied on the left of A_tensor (rows).
  updated_pp <- param_prog
  if (!is.null(param_prog) && !is.null(param_prog@A_tensor) &&
      !isTRUE(param_prog@formatted)) {
    n_cols <- param_prog@x_length + 1L
    P_kron <- Matrix::Diagonal(n_cols) %x% restruct_mat
    updated_pp <- ParamConeProg(
      c_tensor       = param_prog@c_tensor,
      A_tensor       = P_kron %*% param_prog@A_tensor,
      x_length       = param_prog@x_length,
      x_id           = param_prog@x_id,
      parameters     = param_prog@parameters,
      param_id_to_col= param_prog@param_id_to_col,
      param_to_size  = param_prog@param_to_size,
      variables      = param_prog@variables,
      var_id_to_col  = param_prog@var_id_to_col,
      constraints    = param_prog@constraints,
      cone_dims      = param_prog@cone_dims,
      P_tensor       = param_prog@P_tensor,
      lower_bounds   = param_prog@lower_bounds,
      upper_bounds   = param_prog@upper_bounds,
      lb_tensor      = param_prog@lb_tensor,
      ub_tensor      = param_prog@ub_tensor,
      formatted      = TRUE
    )
  }

  list(A = formatted_A, b = formatted_b, param_prog = updated_pp)
}

# -- 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"]]
  ## CVXPY SOURCE: conic_solver.py:358-362
  inv_data[[SOLVER_NEQ_CONSTR]] <- c(
    constr_map[["NonNeg"]], constr_map[["SOC"]],
    constr_map[["PSD"]], constr_map[["SvecPSD"]], constr_map[["ExpCone"]],
    constr_map[["PowCone3D"]], constr_map[["PowConeND"]]
  )

  ## Format constraints: restructure A and b.  Thread `param_prog` so
  ## its `A_tensor` rows get the same permutation, then propagate the
  ## formatted param_prog forward so the Problem-level cache picks up
  ## the post-format version.  See `format_constraints` for details.
  formatted <- format_constraints(
    constraints, data[[SD_A]], data[[SD_B]],
    exp_cone_order = x@EXP_CONE_ORDER,
    param_prog     = data[[SD_PARAM_PROB]]
  )

  ## 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
  solver_data[[LOWER_BOUNDS]] <- data[[LOWER_BOUNDS]]
  solver_data[[UPPER_BOUNDS]] <- data[[UPPER_BOUNDS]]
  ## Carry the formatted (or pass-through) param_prog so it can be
  ## cached for DPP fast-path re-solves and used by backward()/derivative().
  solver_data[[SD_PARAM_PROB]] <- formatted$param_prog
  ## 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)
}

# -- 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)
}