R/142_problems_param_prob.R

In CVXR: Disciplined Convex Optimization

#####
## DO NOT EDIT THIS FILE!! EDIT THE SOURCE INSTEAD: rsrc_tree/problems/param_prob.R
#####

## CVXPY SOURCE: problems/param_prob.py + reductions/dcp2cone/cone_matrix_stuffing.py
## ParamConeProg -- parameterized cone program for DPP fast re-solve
##
## Stores sparse tensor representation of the problem data as a function
## of parameters. apply_parameters() multiplies tensor @ param_vec to
## recover numeric A, b, c (and optionally P for QP).


# -- get_parameter_vector -----------------------------------------------
## CVXPY SOURCE: canonInterface.py lines 28-62
## Builds a flattened parameter vector from current parameter values.
## Length = total_param_size (sum of all param sizes including CONSTANT_ID).
## CONSTANT_ID gets value 1.0 at its column (unless zero_offset = TRUE).

get_parameter_vector <- function(param_to_size, param_id_to_col,
                                  parameters, zero_offset = FALSE,
                                  id_to_param_value = NULL) {
  total <- sum(as.integer(unlist(param_to_size)))
  param_vec <- numeric(total)
  for (pid_str in names(param_id_to_col)) {
    col <- param_id_to_col[[pid_str]]  # 0-based offset
    pid <- as.integer(pid_str)
    if (pid == LINOP_CONSTANT_ID) {
      if (!zero_offset) param_vec[col + 1L] <- 1.0
    } else {
      sz <- param_to_size[[pid_str]]
      ## Override path: if a per-param-id dict is supplied, prefer
      ## those values over the live Parameter objects.  Used by
      ## Problem$derivative() to plug perturbations into the same
      ## tensor multiply that produces (c, A, b) in the forward
      ## direction.
      if (!is.null(id_to_param_value) &&
          !is.null(id_to_param_value[[pid_str]])) {
        param_vec[(col + 1L):(col + sz)] <-
          as.vector(id_to_param_value[[pid_str]])
        next
      }
      ## Find matching Parameter object
      for (p in parameters) {
        if (p@id == pid) {
          val <- value(p)
          if (is.null(val)) {
            cli_abort(c(
              "Problem contains an unspecified parameter.",
              "i" = "Set {.code value({name(p)}) <- <value>} before solving."
            ), class = "ParameterError")   ## CVXPY: raise ParameterError, eval_params.py:15
          }
          param_vec[(col + 1L):(col + sz)] <- as.vector(val)  # column-major
          break
        }
      }
    }
  }
  param_vec
}

# -- ParamConeProg class ------------------------------------------------
## CVXPY SOURCE: cone_matrix_stuffing.py lines 145-318
## Stores tensor representation and applies parameters to produce data.

## CVXPY SOURCE: cone_matrix_stuffing.py:158 — `formatted: bool = False`
## flag indicating whether `A_tensor` has been row-permuted into solver
## cone-interleaved layout by `ConicSolver.format_constraints`.
ParamConeProg <- new_class("ParamConeProg", package = "CVXR",
  properties = list(
    c_tensor     = class_any,  # sparse Matrix: objective tensor (x_length+1, param_size)
    A_tensor     = class_any,  # sparse Matrix: constraint tensor (constr*(x_length+1), param_size)
    P_tensor     = class_any,  # sparse Matrix or NULL: QP tensor (x_length*x_length, param_size)
    x_length     = class_integer,
    x_id         = class_integer,  # variable id for solver inverse data
    parameters   = class_list,
    param_id_to_col = class_list,
    param_to_size   = class_list,
    variables    = class_list,
    var_id_to_col = class_list,
    constraints  = class_list,
    cone_dims    = class_any,
    lower_bounds = class_any,
    upper_bounds = class_any,
    lb_tensor    = class_any,
    ub_tensor    = class_any,
    formatted    = class_logical   # FALSE = pre-format (block layout from
                                   # ConeMatrixStuffing); TRUE = post-format
                                   # (cone-interleaved by format_constraints).
  ),
  constructor = function(c_tensor, A_tensor, x_length, x_id, parameters,
                          param_id_to_col, param_to_size,
                          variables, var_id_to_col, constraints,
                          cone_dims, P_tensor = NULL,
                          lower_bounds = NULL, upper_bounds = NULL,
                          lb_tensor = NULL, ub_tensor = NULL,
                          formatted = FALSE) {
    if (FALSE) new_object(S7_object())  ## S7 static-check guard
    .fast_new(ParamConeProg, S7_object(),
      c_tensor = c_tensor,
      A_tensor = A_tensor,
      P_tensor = P_tensor,
      x_length = as.integer(x_length),
      x_id = as.integer(x_id),
      parameters = parameters,
      param_id_to_col = param_id_to_col,
      param_to_size = param_to_size,
      variables = variables,
      var_id_to_col = var_id_to_col,
      constraints = constraints,
      cone_dims = cone_dims,
      lower_bounds = lower_bounds,
      upper_bounds = upper_bounds,
      lb_tensor = lb_tensor,
      ub_tensor = ub_tensor,
      formatted = as.logical(formatted)
    )
  }
)

# -- DPP tensor contraction helper --------------------------------------
## Contract a DPP parameter tensor with the parameter vector and reshape the
## resulting column-major vec into an (nrow x ncol) sparse matrix, staying
## SPARSE throughout -- the dense nrow*ncol form is never materialized.
##
## `tensor %*% pv_sparse` (pv as a sparse 1-column Matrix) keeps the result a
## 1-column dgCMatrix; its @i slot is exactly the 0-based column-major flat
## index of each stored value, which scatters directly into (row, col) via
## sparseMatrix() -> a general dgCMatrix (no ddiMatrix/dsCMatrix surprises that
## some solver interfaces, e.g. XPRESS, cannot coerce).
##
## `zero_rows` (1-based flat column-major indices, from A_mapping_nonzero_rows)
## are appended with value 0 so they survive as EXPLICIT ZEROS in the result.
## CVXPY SOURCE: canonInterface.py:240-247, which does the same by appending
## `np.zeros(nonzero_rows.size)` to the COO triplets before rebuilding the CSC.
## The product above DROPS them: an entry whose parameter is currently 0
## multiplies out to 0 and leaves the pattern, which is why a DIFFCP gradient
## was exactly 0 at a parameter value of 0. `Matrix::sparseMatrix()` sums
## duplicate (i, j) pairs, so appending a 0 at a position that is already
## present is a no-op, and it RETAINS a stored zero where it is not -- verified,
## along with the fact that the column subsetting, negation and row permutation
## downstream all preserve explicit zeros.
.dpp_contract_reshape <- function(tensor, pv_sparse, nrow, ncol,
                                  zero_rows = NULL) {
  flat <- methods::as(tensor %*% pv_sparse, "CsparseMatrix")  # (nrow*ncol, 1)
  k <- flat@i                                                  # 0-based flat idx
  x <- flat@x
  if (length(zero_rows) > 0L) {
    k <- c(k, as.integer(zero_rows) - 1L)
    x <- c(x, numeric(length(zero_rows)))
  }
  Matrix::sparseMatrix(
    i = (k %% nrow) + 1L,
    j = (k %/% nrow) + 1L,
    x = x,
    dims = c(nrow, ncol))
}

# -- apply_parameters ---------------------------------------------------
## CVXPY SOURCE: cone_matrix_stuffing.py lines 206-240
## Applies current parameter values to tensors, producing A, b, c, (P).
##
## Returns list(c, d, A, b) or list(P, c, d, A, b) if quad_obj.
## c: numeric(x_length), d: numeric(1), A: sparse (m, x_length), b: numeric(m)

## `keep_zeros` mirrors the upstream argument of the same name
## (cone_matrix_stuffing.py:213-214 -> ReducedMat.cache, utilities.py:138-140):
## when TRUE, every entry of `A` that a parameter CAN drive is kept in the
## sparsity pattern as an explicit zero, even where the current parameter value
## makes it numerically zero. Exactly one caller sets it -- the DIFFCP
## interface (diffcp_conif.py:72) -- because the derivative of the solution map
## is taken with respect to the entries of A, and an entry that is not in the
## pattern has no derivative at all.
##
## Measured, on max x1 + 2*x2 s.t. p*x1 + x2 <= 1, 0 <= x1 <= 3, 0 <= x2 <= 5:
##                    CVXR before      CVXR after      CVXPY 1.9.2
##      p = 0         0                -3              -3.000000
##      p = 0.1       -3               -3              -3.000004
## (seeding gradient(x) = sum(x); the objective seed gives -6 throughout.)
## A silently zero gradient, with no error and no warning -- gradient descent
## started at p = 0 never moves.
apply_parameters <- function(param_prog, quad_obj = FALSE,
                             id_to_param_value = NULL,
                             zero_offset = FALSE,
                             keep_zeros = FALSE) {
  pv <- get_parameter_vector(param_prog@param_to_size,
                              param_prog@param_id_to_col,
                              param_prog@parameters,
                              zero_offset = zero_offset,
                              id_to_param_value = id_to_param_value)
  xl <- param_prog@x_length

  ## Objective: c_tensor is (x_length + 1, param_size) sparse
  ## c_tensor %*% pv -> (x_length + 1) vector
  ## First x_length entries = c, last entry = d (constant offset)
  if (is.null(param_prog@c_tensor))
    cli_abort("DPP fast path requires c_tensor. This is a bug -- please report.")
  ## c and the bounds below are genuine dense vectors (one entry per variable),
  ## so the dense product is appropriate. A and P are matrices and stay sparse.
  pv_sparse <- Matrix::Matrix(pv, ncol = 1L, sparse = TRUE)
  c_flat <- as.numeric(param_prog@c_tensor %*% pv)
  c_vec <- c_flat[seq_len(xl)]
  d_offset <- c_flat[xl + 1L]

  ## Constraints: A_tensor is (constr * (x_length + 1), param_size) sparse.
  ## Reshape the contracted column-major vec to (constr, x_length + 1) WITHOUT
  ## densifying; first x_length columns = A, last column = b.
  n_cols <- xl + 1L
  n_rows <- nrow(param_prog@A_tensor) %/% n_cols
  if (n_rows > 0L) {
    zero_rows <- if (isTRUE(keep_zeros)) {
      A_mapping_nonzero_rows(param_prog@A_tensor, xl,
                             const_col = param_prog@param_id_to_col[["-1"]])
    } else NULL
    Ab <- .dpp_contract_reshape(param_prog@A_tensor, pv_sparse, n_rows, n_cols,
                                zero_rows = zero_rows)
    A <- Ab[, seq_len(xl), drop = FALSE]
    b <- as.numeric(Ab[, n_cols])
  } else {
    A <- Matrix::sparseMatrix(i = integer(0), j = integer(0), x = numeric(0),
                               dims = c(0L, xl))
    b <- numeric(0)
  }

  lower_bounds <- param_prog@lower_bounds
  upper_bounds <- param_prog@upper_bounds
  if (!is.null(param_prog@lb_tensor)) {
    lower_bounds <- as.numeric(param_prog@lb_tensor %*% pv)
  }
  if (!is.null(param_prog@ub_tensor)) {
    upper_bounds <- as.numeric(param_prog@ub_tensor %*% pv)
  }

  if (quad_obj && !is.null(param_prog@P_tensor)) {
    ## P_tensor is (x_length * x_length, param_size) sparse; reshape the
    ## contracted column-major vec to (x_length, x_length) without densifying.
    P_mat <- .dpp_contract_reshape(param_prog@P_tensor, pv_sparse, xl, xl)
    list(P = P_mat, c = c_vec, d = d_offset, A = A, b = b,
         lower_bounds = lower_bounds, upper_bounds = upper_bounds)
  } else {
    list(c = c_vec, d = d_offset, A = A, b = b,
         lower_bounds = lower_bounds, upper_bounds = upper_bounds)
  }
}

# -- Derivative-API helpers (Phase 4.3) ---------------------------
## CVXPY SOURCE: cone_matrix_stuffing.py:242-318
##
## Three pure-R linear-algebra functions used by Problem$backward()
## and Problem$derivative() in Phase 4.5.
##
## They depend on the same (c_tensor, A_tensor, x_length, var/param
## column maps) that apply_parameters() above already consumes.

#' Adjoint of the parameter -> (c, d, A, b) tensor map
#'
#' Given derivatives `(delc, delA, delb)` of a downstream objective
#' with respect to the conic problem data, returns derivatives with
#' respect to each `Parameter` (keyed by parameter id).  Mirrors
#' `ParamConeProg.apply_param_jac` in `cone_matrix_stuffing.py:242-280`.
#'
#' Reusing the tensor identity that `apply_parameters` exploits in
#' the forward direction:
#'   `c = c_tensor[1:x_length, ] %*% p`  (omitting the offset row),
#'   `vec(A) // b = A_tensor %*% p`     (column-major; b is last block),
#' the adjoint is just the transpose of the same tensors applied to
#' the stacked (delc, vec(delA), delb).
#'
#' @param param_prog A `ParamConeProg`.
#' @param delc Numeric vector of length `x_length`.
#' @param delA Sparse matrix of shape `(m, x_length)` (same shape as
#'   the conic constraint matrix `A`).
#' @param delb Numeric vector of length `m`.
#' @param active_params Optional character vector of parameter ids to
#'   restrict the output to.  Default: all parameters.
#' @returns A named list mapping `as.character(param_id)` to a numeric
#'   array of the parameter's shape.
#' @keywords internal
apply_param_jac <- function(param_prog, delc, delA, delb,
                            active_params = NULL) {
  if (!is.null(param_prog@P_tensor)) {
    cli_abort("Cannot apply Jacobian when a quadratic objective is present.")
  }
  if (is.null(active_params)) {
    active_params <- vapply(param_prog@parameters,
                             function(p) as.character(id(p)),
                             character(1))
  }

  ## c contribution: delc %*% c_tensor[1:x_length, ]
  ## (drop the last row, which produced the constant offset d).
  xl <- param_prog@x_length
  c_tensor_no_off <- param_prog@c_tensor[seq_len(xl), , drop = FALSE]
  del_param_vec <- as.numeric(matrix(delc, nrow = 1L) %*% c_tensor_no_off)

  ## (A, b) contribution: A_tensor maps p -> column-major flatten of
  ## (cbind(A, b)) of shape (m, x_length + 1).  The adjoint is
  ## A_tensor^T %*% c(vec(delA, F-order), delb).
  delA_csc <- methods::as(delA, "CsparseMatrix")
  flat_delA <- as.numeric(delA_csc)   # column-major flatten
  delAb     <- c(flat_delA, as.numeric(delb))
  del_param_vec <- del_param_vec +
    as.numeric(matrix(delAb, nrow = 1L) %*% param_prog@A_tensor)

  ## Split into per-parameter chunks.
  out <- list()
  for (pid_str in names(param_prog@param_id_to_col)) {
    if (!(pid_str %in% active_params)) next
    col0 <- param_prog@param_id_to_col[[pid_str]]   # 0-based offset
    psize <- param_prog@param_to_size[[pid_str]]
    delta <- del_param_vec[(col0 + 1L):(col0 + psize)]
    ## Reshape to the parameter's shape.
    param <- NULL
    for (p in param_prog@parameters) {
      if (as.character(id(p)) == pid_str) { param <- p; break }
    }
    if (is.null(param)) next
    out[[pid_str]] <- array(as.numeric(delta), dim = param@shape)
  }
  out
}

#' Split a primal solution into per-variable arrays
#'
#' Mirrors `ParamConeProg.split_solution` in
#' `cone_matrix_stuffing.py:282-302`.
#'
#' @param param_prog A `ParamConeProg`.
#' @param sltn Numeric vector of length `x_length` -- the primal `x`
#'   from the conic forward solve.
#' @param active_vars Optional character vector of variable ids to
#'   restrict the output to.  Default: all variables.
#' @returns A named list mapping `as.character(var_id)` to a numeric
#'   array of the variable's shape.
#' @keywords internal
split_solution <- function(param_prog, sltn, active_vars = NULL) {
  if (is.null(active_vars)) {
    active_vars <- vapply(param_prog@variables,
                           function(v) as.character(id(v)),
                           character(1))
  }
  out <- list()
  for (vid_str in names(param_prog@var_id_to_col)) {
    if (!(vid_str %in% active_vars)) next
    col0 <- param_prog@var_id_to_col[[vid_str]]   # 0-based offset
    var <- NULL
    for (v in param_prog@variables) {
      if (as.character(id(v)) == vid_str) { var <- v; break }
    }
    if (is.null(var)) next
    vsize <- prod(var@shape)
    value <- sltn[(col0 + 1L):(col0 + vsize)]
    out[[vid_str]] <- array(as.numeric(value), dim = var@shape)
  }
  out
}

#' Adjoint of `split_solution`
#'
#' Given a named list mapping variable ids to their (shape-aligned)
#' deltas, packs them into a single length-`x_length` numeric vector
#' in `var_id_to_col` order.  Mirrors `ParamConeProg.split_adjoint`
#' in `cone_matrix_stuffing.py:304-318`.
#'
#' @param param_prog A `ParamConeProg`.
#' @param del_vars A named list of `as.character(var_id)` -> array.
#' @returns Numeric vector of length `x_length`.
#' @keywords internal
split_adjoint <- function(param_prog, del_vars) {
  var_vec <- numeric(param_prog@x_length)
  for (vid_str in names(del_vars)) {
    col0 <- param_prog@var_id_to_col[[vid_str]]
    if (is.null(col0)) next
    delta <- del_vars[[vid_str]]
    flat <- as.numeric(delta)   # column-major if matrix
    var_vec[(col0 + 1L):(col0 + length(flat))] <- flat
  }
  var_vec
}