R/006_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
#####

## R-SPECIFIC: 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())

# -- Fast constructor for trusted internal node construction -----------------
## Bypasses S7's `new_object()` for the hot-path Expression-tree construction.
##
## ===========================================================================
## AUTHORITY: this is the canonical reference for `.fast_new`.  Every CVXR
## constructor MUST use it (CLAUDE.md constraint 17).  Background: ADR
## D_PERF.1 in notes/decisions.md, migration status in
## notes/lever2_migration_status.md.
##
## USAGE PATTERN (every constructor, no exceptions):
##
##   constructor = function(...) {
##     if (FALSE) new_object(S7_object())  ## S7 static-check guard
##     ...                                  ## arg coercion / validation
##     .fast_new(ThisClassName, S7_object(),
##       id     = next_expr_id(),
##       .cache = new.env(parent = emptyenv()),
##       args   = <args>,
##       shape  = <shape>,
##       <other named properties>
##     )
##   }
##
## The `if (FALSE) new_object(...)` line satisfies S7's
## `check_S7_constructor` static AST check (which scans for a literal
## `new_object` call in the body) without executing it at runtime.
## ===========================================================================
##
## new_object() pays for: sys.function(-1) call-stack walk to find class,
## abstract-class check, named-args validation, has_setter dispatch loop,
## attribute deduplication, full recursive validate().  Hot-path microbench:
## 67 us per Variable(10) vs 3.4 us via .fast_new -- 20x speedup.
##
## CRAN safety: this helper uses ONLY exported S7 API (`S7_object` sentinel,
## the `@parent` / `@name` / `@package` accessors -- all part of the
## documented `S7_class` slot interface).  No `:::` or `getNamespace()` reach
## into S7 internals.
##
## Trade-off vs new_object(): NO property validation -- relies on the
## constructor body to type-correct its inputs.  All CVXR constructors
## already do this via explicit `validate_*()` / `cli_abort()` calls; the
## audit confirms zero classes in rsrc_tree use S7 property setters or
## validators (see notes/lever2_migration_status.md).
##
## Class dispatch vector ("CVXR::Foo", "CVXR::Bar", ..., "S7_object") is
## computed once per class via a parent-chain walk and cached in a closure-
## captured environment.  Subsequent constructions reuse the cached vector.
.fast_new <- local({
  cache <- new.env(hash = TRUE, parent = emptyenv())

  ## Reproduces S7's `class_dispatch()` output using only public accessors.
  ## Walks the class's parent chain, terminating at the `S7_object` sentinel.
  build_dispatch <- function(cls) {
    out <- character(0L)
    cur <- cls
    while (inherits(cur, "S7_class") && !identical(cur, S7::S7_object)) {
      out <- c(out, paste0(cur@package, "::", cur@name))
      cur <- cur@parent
    }
    c(out, "S7_object")
  }

  function(.class, .parent, ...) {
    cn <- .class@name
    cd <- get0(cn, envir = cache, ifnotfound = NULL)
    if (is.null(cd)) {
      cd <- build_dispatch(.class)
      assign(cn, cd, envir = cache)
    }
    obj <- .parent
    attributes(obj) <- c(list(class = cd, S7_class = .class), list(...))
    obj
  }
})

#' 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
## MEASUREMENT VARIANT B (2026-08-14, uncommitted -- see ADR D_PERF.8, status
## PROPOSED). Two changes, semantics IDENTICAL to what they replace: the same
## environment, the same sentinel, the same by-reference `.cache`.
##
##   `x@.cache`  -> `attr(x, ".cache", exact = TRUE)`   0.77us -> 0.12us
##   `exists()` + `get()` -> `get0()`                   one lookup, not two
##
## The `@` bypass is legitimate here ONLY because `.fast_new` writes every
## property as a plain attribute and CVXR has zero getters/validators; if that
## ever stops holding, these four functions silently read past the getter. The
## reasoning, the measurements and the hazard are in D_PERF.8 -- read it before
## copying this pattern anywhere else.
cache_get <- function(x, key) {
  ## was: cache <- x@.cache
  ##      if (exists(key, envir = cache, inherits = FALSE)) {
  ##        get(key, envir = cache, inherits = FALSE)
  ##      } else .NOT_CACHED
  get0(key, envir = attr(x, ".cache", exact = TRUE),   ## @-bypass: D_PERF.8
       inherits = FALSE, ifnotfound = .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) {
  ## was: assign(key, val, envir = x@.cache)
  assign(key, val, envir = attr(x, ".cache", exact = TRUE))  ## @-bypass: D_PERF.8
  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) {
  ## was: exists(key, envir = x@.cache, inherits = FALSE)
  ## @-bypass: D_PERF.8
  exists(key, envir = attr(x, ".cache", exact = TRUE), 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) {
  ## was: rm(list = ls(x@.cache, all.names = TRUE), envir = x@.cache)
  e <- attr(x, ".cache", exact = TRUE)   ## @-bypass: D_PERF.8
  rm(list = ls(e, all.names = TRUE), envir = e)
  invisible(NULL)
}

# -- Hot property accessors (ADR D_PERF.8) -----------------------------
## `x@prop` costs ~533ns; `attr(x, "prop", exact = TRUE)` returns the same value
## in ~82ns, because the excess is S3 dispatch (`@` is a base primitive generic
## whose S7 method sits in baseenv()'s S3 methods table, so every `CVXR::` entry
## in the class vector misses first) rather than S7's property machinery.
##
## These four properties are 81.5% of all `@` reads on a solve -- per `large_lp`
## solve: args 65,682, shape 25,185, id 12,541, attributes 6,048 of 134,206.
## A helper costs ~82ns more than a bare `attr()` and still keeps ~85% of the
## saving, while buying what a bare `attr()` at a thousand call sites cannot:
## the reasoning lives in ONE place, call sites stay readable, and the entire
## sweep reverts by redefining four functions.
##
##   min / median, 200k iterations:  x@args 533/697ns
##                                   attr(x,"args",exact=TRUE) ~0/82ns
##                                   .args(x) 82/164ns
##
## READS ONLY. An assignment target keeps `@`: `x@attributes[["k"]] <- v`
## expands through `@<-`, and the read form would need a `.attributes<-`
## replacement function that does not exist.
##
## SAFE ONLY FOR CVXR OBJECTS. The equivalence rests on `.fast_new` writing every
## property as a plain attribute (constraint 17) and CVXR having zero getters,
## setters and validators across its 162 classes. It does NOT extend to S7 class
## objects, to S4/Matrix objects, or to any class that later gains a getter.
## `test-attr-property-equivalence.R` enumerates every class and fails the day
## that stops holding, which is when these helpers must be reverted to `@`.
##
## `exact = TRUE` is mandatory: without it `attr()` partial-matches, so a lookup
## of "arg" would silently return `args`.
##
## Call sites were converted mechanically by a parser-driven rewriter
## (`scripts/rewrite_at_to_attr.R`), which locates `@` nodes in parse data, edits by
## exact source span, and verifies each file by comparing its rewritten AST
## against the original AST transformed independently.
## @noRd
.args       <- function(x) attr(x, "args", exact = TRUE)
.shape      <- function(x) attr(x, "shape", exact = TRUE)
.id         <- function(x) attr(x, "id", exact = TRUE)
.attributes <- function(x) attr(x, "attributes", exact = TRUE)

## "the shape of x's i-th argument" -- 89 sites, by far the most common
## composite read in the package.
##
## READABILITY, not just speed. `x@args[[1L]]@shape` reads left to right; the
## mechanical rewrite `.shape(.args(x)[[1L]])` has to be read inside out, which
## is worse than what it replaced. Naming the idiom is better than either form:
## `.arg_shape(x)` says what it means, and it is also one call rather than two.
## Where a mechanical sweep makes code less readable, the fix is usually that
## the sweep found the wrong abstraction, not that the sweep was wrong.
## @noRd
.arg_shape <- function(x, i = 1L)
  attr(attr(x, "args", exact = TRUE)[[i]], "shape", exact = TRUE)

# -- Dual re-keying ----------------------------------------------------
## Every `reduction_invert` that renames constraints has to move a solution's
## dual variables from the NEW constraint ids back onto the ORIGINAL ones.
## CVXPY writes that as a dict comprehension over `cons_id_map.items()`
## (canonicalization.py:77-83), which is O(n) because Python dict insert and
## lookup are O(1).
##
## Transliterated literally into R it is O(n^2) TWICE OVER, because an R named
## list is not a hash map: `l[[key]]` is a linear scan of `names(l)`, and
## `l[[key]] <- v` reallocates. Measured on a 1000-constraint problem, the loop
## in `reduction_invert(Canonicalization)` cost 10.8 ms where this costs 1.0 ms
## (11.1x), and the growth is quadratic -- at n = 4000 it is 216 ms vs 4.7 ms.
##
## The fix is not a trick: do the n lookups as ONE vectorized `match()`, which
## is a single C call over a hash table, and build the result by subsetting
## once instead of growing. Same complexity CVXPY has, in R's idiom.
## See notes/string_key_hashing_sweep_2026-08-13.md and ADR D_PERF.7.
##
## `id_map` is an environment (Canonicalization, ExactCone2Cone) or a named
## list (CvxAttr2Constr) mapping `as.character(old_id)` -> new id. Key ORDER is
## preserved exactly as before (`ls()` sorts; `names()` does not), since callers
## and tests may depend on it.
## @noRd
.remap_by_id_map <- function(dual_vars, id_map) {
  if (length(dual_vars) == 0L) return(list())
  is_env   <- is.environment(id_map)
  old_ids  <- if (is_env) ls(id_map, all.names = TRUE) else names(id_map)
  if (length(old_ids) == 0L) return(list())
  new_vals <- if (is_env) mget(old_ids, envir = id_map) else id_map[old_ids]
  new_ids  <- as.character(unlist(new_vals, use.names = FALSE))
  idx      <- match(new_ids, names(dual_vars))
  keep     <- !is.na(idx)
  if (!any(keep)) return(list())
  out <- dual_vars[idx[keep]]
  names(out) <- old_ids[keep]
  ## `l[[k]] <- NULL` deletes rather than stores, so the loop this replaces
  ## could never produce a NULL entry; drop any so the result is identical.
  null_out <- vapply(out, is.null, logical(1))
  if (any(null_out)) out <- out[!null_out]
  out
}

# -- Argument predicate helpers ----------------------------------------
## Shorthand for the recurring vapply(x@args, pred, logical(1)) pattern.
## @noRd
.all_args <- function(x, pred) all(vapply(.args(x), pred, logical(1)))
.any_args <- function(x, pred) any(vapply(.args(x), 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())
  ## PERFORMANCE (2026-08-13). This was an R-level loop over an environment used
  ## as a hash set: one `new.env()` per call plus `exists()`/`assign()` and an
  ## `as.character()` key per item. `variables()`/`parameters()`/`constants()`/
  ## `atoms()` re-dedup each subtree at EVERY node of the expression tree
  ## (utilities/canonical.R:24-44, matching CVXPY canonical.py:56-73, which
  ## deliberately does not cache them), so the call volume is large: 8540 calls
  ## per `large_lp` solve, and the env allocations alone were ~47% of that
  ## cell's 18k `new.env()` calls. Measured on 1000 items / 500 distinct:
  ## 2763us -> 950us (x2.9). Semantics are unchanged -- dedup is still by `@id`
  ## and still keeps first-occurrence order; `duplicated()` does both in C.
  ## `id` is declared `class_integer` (utilities/canonical.R:7), so `integer(1)`
  ## is the exact vapply type.
  ## A single-element list is already unique, and that is the COMMON case: most
  ## of the 8540 calls in a `large_lp` solve come from leaf and near-leaf nodes.
  ## Measured at n=1: 4.77us (old) / 3.08us (vapply+duplicated) / 0.23us (here).
  if (length(lst) == 1L) return(lst)
  ids <- vapply(lst, function(item) .id(item), integer(1L), USE.NAMES = FALSE)
  lst[!duplicated(ids)]
}

# -- 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(.shape(x)))

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

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

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

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