R/005_zzz_R_specific_utility.R

In CVXR: Disciplined Convex Optimization

#####
## DO NOT EDIT THIS FILE!! EDIT THE SOURCE INSTEAD: rsrc_tree/zzz_R_specific/utility.R
#####

## CVXPY SOURCE: (R-specific utility functions)
## Caching helpers and common utilities

# -- Caching -----------------------------------------------------------

## Sentinel object for cache misses.  Using a dedicated environment avoids
## the NULL-sentinel pitfall: if a computation legitimately returns NULL or
## FALSE, a NULL default would cause re-computation on every access.
.NOT_CACHED <- new.env(parent = emptyenv())

#' Get a cached value from an expression's cache environment
#' @param x An expression with a `.cache` property
#' @param key Character key
#' @returns The cached value, or `.NOT_CACHED` sentinel on miss
#' @noRd
cache_get <- function(x, key) {
  cache <- x@.cache
  if (exists(key, envir = cache, inherits = FALSE)) {
    get(key, envir = cache, inherits = FALSE)
  } else {
    .NOT_CACHED
  }
}

#' Test whether a value is the cache-miss sentinel
#' @param val Value returned by `cache_get`
#' @returns Logical
#' @noRd
cache_miss <- function(val) identical(val, .NOT_CACHED)

#' Set a cached value in an expression's cache environment
#' @param x An expression with a `.cache` property
#' @param key Character key
#' @param val Value to cache
#' @returns Invisible NULL (side effect: modifies cache)
#' @noRd
cache_set <- function(x, key, val) {
  assign(key, val, envir = x@.cache)
  invisible(NULL)
}

#' Check if a key exists in the cache
#' @param x An expression with a `.cache` property
#' @param key Character key
#' @returns Logical
#' @noRd
cache_has <- function(x, key) {
  exists(key, envir = x@.cache, inherits = FALSE)
}

#' Clear all cached values for an expression
#' @param x An expression with a `.cache` property
#' @returns Invisible NULL
#' @noRd
cache_clear <- function(x) {
  rm(list = ls(x@.cache, all.names = TRUE), envir = x@.cache)
  invisible(NULL)
}

# -- Argument predicate helpers ----------------------------------------
## Shorthand for the recurring vapply(x@args, pred, logical(1)) pattern.
## @noRd
.all_args <- function(x, pred) all(vapply(x@args, pred, logical(1)))
.any_args <- function(x, pred) any(vapply(x@args, pred, logical(1)))

# -- C-order reshape helper --------------------------------------------
## CVXPY's np.reshape() defaults to C-order (row-major). R's matrix()
## defaults to column-major (Fortran order). This helper bridges the gap.
## Used by save_dual_value methods on cone constraints.

#' Reshape vector into matrix using C-order (row-major)
#' @param x Numeric vector
#' @param nrow Number of rows
#' @param ncol Number of columns
#' @returns Matrix with C-order (row-major) fill
#' @noRd
.reshape_c_order <- function(x, nrow, ncol) {
  matrix(x, nrow = nrow, ncol = ncol, byrow = TRUE)
}

# -- Type checking helpers ---------------------------------------------

#' Convert a value to an expression (promoting scalars/matrices to Constant)
#' @param x A value: numeric, matrix, or Expression
#' @returns An Expression object
#' @noRd
#' @note Expression and Constant classes defined in Phase 1
as_expr <- function(x) {
  # Expression and Constant will be available after Phase 1
  if (inherits(x, "CVXR::Expression")) {
    x
  } else if (is.numeric(x) || is.complex(x) || is.logical(x) || inherits(x, "Matrix") || inherits(x, "sparseVector")) {
    Constant(x)
  } else {
    cli_abort("Cannot convert object of class {.cls {class(x)}} to a CVXR Expression.")
  }
}

#' Convert a value to a CVXR Expression
#'
#' Wraps numeric vectors, matrices, and Matrix package objects as CVXR
#' [Constant] objects. Values that are already CVXR expressions are returned
#' unchanged.
#'
#' @section Matrix package interoperability:
#' Objects from the \pkg{Matrix} package (`dgCMatrix`, `dgeMatrix`,
#' `ddiMatrix`, `sparseVector`, etc.) are S4 classes.
#' Because S4 dispatch preempts S7/S3 dispatch, **raw Matrix objects cannot be
#' used directly with CVXR operators** (`+`, `-`, `*`, `/`, `%*%`, `>=`, `==`,
#' etc.).
#'
#' Use `as_cvxr_expr()` to wrap a Matrix object as a CVXR [Constant] before
#' combining it with CVXR variables or expressions.  This preserves sparsity
#' (unlike [as.matrix()], which densifies).
#'
#' Base R `matrix` and `numeric` objects work natively with CVXR operators ---
#' no wrapping is needed.
#'
#' @param x A numeric vector, matrix, [Matrix::Matrix-class] object,
#'   [Matrix::sparseVector-class] object, or CVXR expression.
#' @return A CVXR expression (either the input unchanged or wrapped in
#'   [Constant]).
#' @examples
#' x <- Variable(3)
#'
#' ## Sparse Matrix needs as_cvxr_expr() for CVXR operator dispatch:
#' A <- Matrix::sparseMatrix(i = 1:3, j = 1:3, x = 1.0)
#' expr <- as_cvxr_expr(A) %*% x
#'
#' ## All operators work with wrapped Matrix objects:
#' y <- Variable(c(3, 3))
#' expr2 <- as_cvxr_expr(A) + y
#' constr <- as_cvxr_expr(A) >= y
#'
#' ## Base R matrix works natively (no wrapping needed):
#' D <- matrix(1:9, 3, 3)
#' expr3 <- D %*% x
#' @export
as_cvxr_expr <- function(x) as_expr(x)

# -- Numeric broadcasting (R-specific) ---------------------------------
## R's matrix arithmetic does NOT broadcast like numpy.
## e.g., matrix(1:2, 2, 1) + matrix(1, 1, 1) -> "non-conformable arrays"
## This helper broadcasts a value to a target shape for numeric evaluation.

#' Broadcast a numeric value to a target shape
#' @param val A numeric matrix or vector
#' @param target_shape Integer(2) target shape
#' @returns A matrix with dimensions matching target_shape
#' @noRd
.broadcast_numeric <- function(val, target_shape) {
  if (!is.matrix(val)) val <- matrix(val, nrow = length(val), ncol = 1L)
  vdim <- dim(val)
  if (identical(vdim, as.integer(target_shape))) return(val)
  ## Scalar -> full matrix
  if (vdim[1L] == 1L && vdim[2L] == 1L) {
    return(matrix(val[1L, 1L], target_shape[1L], target_shape[2L]))
  }
  ## Column (n,1) -> (n,m)
  if (vdim[1L] == target_shape[1L] && vdim[2L] == 1L && target_shape[2L] > 1L) {
    return(matrix(val[, 1L], target_shape[1L], target_shape[2L]))
  }
  ## Row (1,m) -> (n,m)
  if (vdim[1L] == 1L && vdim[2L] == target_shape[2L] && target_shape[1L] > 1L) {
    return(matrix(val[1L, ], target_shape[1L], target_shape[2L], byrow = TRUE))
  }
  val
}

# -- Shape utilities ---------------------------------------------------

#' Validate and normalize a shape to integer vector of length 2
#' @param shape Shape specification (integer vector, single integer, or NULL)
#' @returns Integer vector of length 2
#' @noRd
validate_shape <- function(shape) {
  if (is.null(shape)) {
    return(c(1L, 1L))
  }
  shape <- as.integer(shape)
  if (length(shape) == 1L) {
    shape <- c(shape, 1L)
  }
  if (length(shape) != 2L) {
    cli_abort("Shape must be a vector of length 1 or 2.")
  }
  if (any(shape <= 0L)) {
    cli_abort("Shape dimensions must be positive.")
  }
  shape
}

#' Check if a shape represents a scalar
#' @param shape Integer vector of length 2
#' @returns Logical
#' @noRd
is_scalar_shape <- function(shape) {
  all(shape == c(1L, 1L))
}

# -- Dedup utility -----------------------------------------------------
## CVXPY SOURCE: utilities/deterministic.py::unique_list

#' Deduplicate a list of expression objects by their \code{@id}
#' @param lst A list of objects with an \code{@id} property
#' @returns A deduplicated list preserving first-occurrence order
#' @noRd
unique_list <- function(lst) {
  if (length(lst) == 0L) return(list())
  ## Use environment as hash set for O(1) lookup instead of O(n) %in%
  seen_env <- new.env(hash = TRUE, parent = emptyenv())
  result <- vector("list", length(lst))
  n <- 0L
  for (item in lst) {
    key <- as.character(item@id)
    if (!exists(key, envir = seen_env, inherits = FALSE)) {
      assign(key, TRUE, envir = seen_env)
      n <- n + 1L
      result[[n]] <- item
    }
  }
  result[seq_len(n)]
}

# -- Shape query helpers (work on Expression objects) ------------------
## CVXPY SOURCE: expressions/expression.py

#' Total number of elements in an expression shape
#' @param x An expression (with \code{@shape})
#' @returns Integer
#' @noRd
expr_size <- function(x) as.integer(prod(x@shape))

#' Number of dimensions in an expression shape
#' @param x An expression (with \code{@shape})
#' @returns Integer
#' @noRd
expr_ndim <- function(x) length(x@shape)

#' Is the expression a scalar (all shape dims are 1)?
#' @param x An expression
#' @returns Logical
#' @noRd
expr_is_scalar <- function(x) all(x@shape == 1L)

#' Is the expression a column or row vector?
#' @param x An expression
#' @returns Logical
#' @noRd
expr_is_vector <- function(x) {
  nd <- length(x@shape)
  nd <= 1L || (nd == 2L && min(x@shape) == 1L)
}

#' Is the expression a matrix (both dims > 1)?
#' @param x An expression
#' @returns Logical
#' @noRd
expr_is_matrix <- function(x) {
  nd <- length(x@shape)
  nd == 2L && x@shape[1L] > 1L && x@shape[2L] > 1L
}

# -- Curvature string helper ------------------------------------------
## CVXPY SOURCE: expressions/expression.py::curvature property

#' Get the curvature string for an expression
#' @param x An expression object
#' @returns Character: "CONSTANT", "AFFINE", "CONVEX", "CONCAVE", or "UNKNOWN"
#' @noRd
expr_curvature <- function(x) {
  if (is_constant(x)) CONSTANT_CURV
  else if (is_affine(x)) AFFINE
  else if (is_convex(x)) CONVEX
  else if (is_concave(x)) CONCAVE
  else UNKNOWN_CURVATURE
}

# -- Sign string helper -----------------------------------------------
## CVXPY SOURCE: expressions/expression.py::sign property

#' Get the sign string for an expression
#' @param x An expression object
#' @returns Character: "ZERO", "POSITIVE", "NEGATIVE", or "UNKNOWN"
#' @noRd
expr_sign_str <- function(x) {
  if (is_zero(x)) ZERO_SIGN
  else if (is_nonneg(x)) NONNEG_SIGN
  else if (is_nonpos(x)) NONPOS_SIGN
  else UNKNOWN_SIGN
}

# -- Safe eigenvalue decomposition for symmetric/Hermitian matrices ----
## Apple Accelerate's zheev (complex Hermitian eigenvalue) segfaults on
## macOS ARM64 when R is linked against vecLib.  R's bundled reference
## LAPACK is not affected.  This helper routes around the bug:
##   - Real matrix:        eigen(A, symmetric=TRUE) -- uses dsyev, safe.
##   - Complex, Im all 0:  Re(A) then dsyev -- avoids zheev entirely.
##   - Truly complex:      eigen(A, symmetric=FALSE) + sort -- uses zgeev.
## Mirrors numpy.linalg.eigvalsh() semantics.

#' Safe eigenvalues/vectors of a symmetric/Hermitian matrix
#'
#' Drop-in replacement for \code{eigen(A, symmetric = TRUE)} that avoids
#' the \code{zheev} segfault in Apple Accelerate on macOS ARM64.
#'
#' @param A A symmetric (real) or Hermitian (complex) matrix.
#' @param only_values If \code{TRUE}, return only eigenvalues (faster).
#' @returns When \code{only_values = TRUE}, a list with \code{$values}
#'   (numeric, decreasing order).
#'   When \code{only_values = FALSE}, a list with \code{$values} and
#'   \code{$vectors}, like \code{eigen()}.
#' @noRd
.eigvalsh <- function(A, only_values = TRUE) {
  if (!is.complex(A)) {
    ## Real matrix -- dsyev is safe
    return(eigen(A, symmetric = TRUE, only.values = only_values))
  }
  ## Complex matrix -- check if imaginary part is all zero
  if (all(Im(A) == 0)) {
    ## Strip +0i to route through dsyev instead of zheev
    return(eigen(Re(A), symmetric = TRUE, only.values = only_values))
  }
  ## Truly complex Hermitian -- use zgeev (symmetric=FALSE) to avoid zheev
  eig <- eigen(A, symmetric = FALSE, only.values = only_values)
  ## Eigenvalues of a Hermitian matrix are guaranteed real
  vals <- Re(eig$values)
  ## symmetric=FALSE returns eigenvalues sorted by decreasing modulus;
  ## re-sort by decreasing real value to match symmetric=TRUE convention
  ord <- order(vals, decreasing = TRUE)
  eig$values <- vals[ord]
  if (!only_values && !is.null(eig$vectors)) {
    eig$vectors <- eig$vectors[, ord, drop = FALSE]
  }
  eig
}