R/130_zzz_R_specific_operators.R

In CVXR: Disciplined Convex Optimization

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

## R-SPECIFIC: Operator dispatch for CVXR expressions
##
## R's operator dispatch differs from Python's __add__/__radd__ etc.
## We use two mechanisms:
##   1. S3 Ops group handler for +, -, *, / (handles both unary and binary)
##   2. S7 method() for %*% (not part of Ops group)
##
## Why S3 Ops instead of S7 method() for arithmetic:
##   S7 `method("-", list(Expression, class_any))` creates `-.CVXR::Expression`
##   which takes precedence over `Ops.CVXR::Expression` for BOTH unary and
##   binary `-`. The S7 method would fail for unary because e2 is missing.
##   Using S3 Ops handler for all arithmetic avoids this limitation.


# -- S3 Ops group handler --------------------------------------------
# Handles +, -, *, / for Expression objects (both unary and binary).
# Registered via registerS3method() in .onLoad().

.cvxr_Ops_handler <- function(e1, e2) {
  ## Unary operations: -x, +x, !x
  if (missing(e2)) {
    switch(.Generic,
      "-" = return(NegExpression(e1)),
      "+" = return(e1),
      "!" = return(Not(e1)),
      cli_abort("Unary {.val {(.Generic)}} not defined for CVXR expressions.")
    )
  }

  ## Binary: coerce non-Expression operands to Constant
  e1 <- as_expr(e1)
  e2 <- as_expr(e2)
  switch(.Generic,
    ## Arithmetic
    "+"  = { bcast <- broadcast_args(e1, e2); AddExpression(list(bcast[[1L]], bcast[[2L]])) },
    "-"  = { bcast <- broadcast_args(e1, NegExpression(e2)); AddExpression(list(bcast[[1L]], bcast[[2L]])) },
    "*"  = Multiply(e1, e2),
    "/"  = DivExpression(e1, e2),
    "^"  = .power_atom(e1, e2),
    ## Comparison -> Constraint objects
    "==" = Equality(e1, e2),
    "<=" = Inequality(e1, e2),         # lhs <= rhs
    ">=" = Inequality(e2, e1),         # lhs >= rhs  ==>  rhs <= lhs
    ## Boolean logic
    "&"  = And(e1, e2),
    "|"  = Or(e1, e2),
    cli_abort("Operator {.val {(.Generic)}} not yet implemented for CVXR expressions.")
  )
}

# -- chooseOpsMethod -------------------------------------------------
# When both operands have Ops methods (e.g., Matrix + Expression),
# this ensures CVXR's handler wins.  Requires R >= 4.3.

.cvxr_chooseOpsMethod <- function(x, y, mx, my, cl, reverse) TRUE

# -- %*% dispatch via S7 --------------------------------------------
# %*% is NOT in the Ops group; S7 method() works directly.

## CVXPY: Expression.__matmul__ / __rmatmul__ (expression.py:842-843,
## 894-895) raise "Scalar operands are not allowed, use '*' instead" when
## either operand has shape (). R has no 0-d value: a bare length-1 vector is
## also R's 1-vector, which base R's `%*%` promotes to 1 x 1 and accepts
## whenever the other operand's inner dimension is 1 (?matmult), so
## `x(n x 1) %*% 2` and `2 %*% y(1 x n)` are valid R and were valid in CVXR
## 1.9.2. Reject only what base R rejects too, with CVXPY's message; a
## conformable bare scalar goes on as a 1 x 1 constant through the
## inner-dimension check in mul_shapes_promote(), as a (1, 1) array does in
## CVXPY. (1.9.2.9500 rejected every bare scalar, which broke a CRAN reverse
## dependency; recorded deviation from CVXPY, 2026-10-17.)
.check_matmul_scalar <- function(v, other, side = c("right", "left")) {
  if (is.null(dim(v)) && length(v) == 1L) {
    side <- match.arg(side)
    inner <- if (side == "right") other@shape[2L] else other@shape[1L]
    if (inner != 1L) {
      cli_abort("Scalar operands are not allowed, use {.code *} instead")
    }
  }
  invisible(NULL)
}

method(`%*%`, list(Expression, class_any)) <- function(x, y) {
  .check_matmul_scalar(y, x, "right")
  MulExpression(x, as_expr(y))
}

method(`%*%`, list(class_any, Expression)) <- function(x, y) {
  .check_matmul_scalar(x, y, "left")
  MulExpression(as_expr(x), y)
}

## Resolve ambiguity when both operands are Expressions
method(`%*%`, list(Expression, Expression)) <- function(x, y) {
  MulExpression(x, y)
}

# -- [ indexing via S7 ------------------------------------------------
# S7 method for `[` on Expression objects.
# Handles:
#   x[i, j]  -- both row and column specified
#   x[i, ]   -- row only (j is missing -> all columns)
#   x[, j]   -- column only (i is missing -> all rows)
#   x[i]     -- single index for vectors; for matrices, selects rows

method(`[`, Expression) <- function(x, i, j, ..., drop = FALSE) {
  ## Build key: list(row_indices, col_indices)
  ## NULL means "all" for that dimension
  if (missing(i) && missing(j)) {
    ## x[] -- return as-is
    return(x)
  }
  ## A user NULL selects nothing in base R, but NULL in a key means "all".
  if ((!missing(i) && is.null(i)) || (!missing(j) && is.null(j))) {
    cli_abort("A {.code NULL} subscript is not supported on an expression.")
  }
  ## Determine nargs: were we called as x[i] or x[i, j]?
  ## R passes nargs() but we detect via sys.call
  cl <- sys.call()
  has_comma <- length(cl) > 3L || (!missing(j))
  ## Alternative detection: if called as x[i], j won't be in the call at all
  ## We check nargs() to distinguish x[i] from x[i,]
  if (!has_comma && !missing(i) && missing(j)) {
    ## x[i] -- single subscript
    ## Step 1: Detect special indexing FIRST (before vector-shape branches)
    ## 2-col matrix and logical matrix work on vectors AND matrices
    if (is.matrix(i) && ncol(i) == 2L && (is.integer(i) || is.double(i))) {
      ## 2-column matrix indexing: each row is a (row, col) pair
      return(SpecialIndex(x, i))
    }
    if (is.matrix(i) && is.logical(i)) {
      ## Logical matrix indexing: select where TRUE
      return(SpecialIndex(x, i))
    }
    ## Step 2: For non-matrix index types, vectors use standard Index
    if (x@shape[2L] == 1L) {
      key <- list(i, NULL)
    } else if (x@shape[1L] == 1L) {
      key <- list(NULL, i)
    } else {
      ## Matrix with scalar/vector single subscript
      if (is.logical(i) && !is.matrix(i)) {
        ## Logical vector on matrix: convert to logical matrix
        if (length(i) != prod(x@shape))
          cli_abort("Logical index length ({length(i)}) must match expression size ({prod(x@shape)}).")
        i_mat <- matrix(i, nrow = x@shape[1L], ncol = x@shape[2L])
        return(SpecialIndex(x, i_mat))
      }
      if ((is.integer(i) || is.double(i)) && !is.matrix(i)) {
        ## Linear integer indexing on matrix (column-major)
        if (is.double(i)) i <- as.integer(i)
        return(SpecialIndex(x, i))
      }
      cli_abort("Unsupported single-index type for matrix expressions.")
    }
  } else {
    ## x[i, j] or x[i, ] or x[, j]
    key <- list(
      if (missing(i)) NULL else i,
      if (missing(j)) NULL else j
    )
  }
  Index(x, key)
}

# -- S3 Math group handler -----------------------------------------
# Dispatches base R Math generics (abs, exp, log, sqrt, etc.) to CVXR atoms.
# Registered via registerS3method() in .onLoad().
# NOTE: H1 -- cumsum, cummax, cumprod, cummin are Math group, NOT Summary
# NOTE: H2 -- sign() routes to the Sign atom (DQCP).  expr_sign() is
#             the distinct *curvature-query* generic.
# NOTE: H3 -- log2, log10 routed through Log atom; log1p is the Log1p atom

.cvxr_Math_handler <- function(x, ...) {
  switch(.Generic,
    "abs"     = Abs(x),
    "exp"     = Exp(x),
    "log"     = if (...length() == 0L) Log(x) else .log_base(x, ...),
    ## NOTE: H4 -- `power()`, NOT `Power()`.  CVXPY SOURCE:
    ##        atoms/elementwise/sqrt.py:24 `return power(x, Fraction(1, 2))`,
    ##        i.e. the wrapper, whose `approx = TRUE` default canonicalizes to
    ##        SOC via PowerApprox.  Constructing `Power` directly takes the
    ##        exact path and emits PowCone3D, which no SOC-only solver (ECOS,
    ##        and every MI-SOCP interface) can accept -- so `sqrt()` alone made
    ##        those solvers reject the problem outright, and SCS return a
    ##        silently wrong answer on the DQCP form.  p = 1/2 is exactly
    ##        representable, so the SOC form is exact, not approximate.
    "sqrt"    = .power_atom(x, 0.5),
    "cumsum"  = Cumsum(x),
    "sign"    = Sign(x),
    "log2"    = DivExpression(Log(x), Constant(log(2))),
    "log10"   = DivExpression(Log(x), Constant(log(10))),
    "log1p"   = Log1p(x),
    "cummax"  = Cummax(x),
    "cumprod" = Cumprod(x),
    "ceiling"  = Ceil(x),
    "floor"    = Floor(x),
    ## Trigonometric / hyperbolic (CVXPY 1.9.0 DNLP atoms; trig.R, hyperbolic.R).
    ## NOTE: no cosh atom exists upstream -> cosh(expr) falls through to abort.
    "sin"     = Sin(x),
    "cos"     = Cos(x),
    "tan"     = Tan(x),
    "sinh"    = Sinh(x),
    "tanh"    = Tanh(x),
    "asinh"   = Asinh(x),
    "atanh"   = Atanh(x),
    "cummin"  = cli_abort("{.fn cummin} is not yet supported for CVXR expressions."),
    cli_abort("{.fn {(.Generic)}} is not supported for CVXR expressions.")
  )
}

## base R's log(x, base); numpy's log has no base. log(base) folds to one
## constant, and a base below 1 divides by a negative, flipping curvature.
.log_base <- function(x, base) {
  if (!is.numeric(base) || length(base) != 1L || is.na(base) ||
      base <= 0 || base == 1) {
    cli_abort("{.arg base} must be a positive number other than 1.")
  }
  DivExpression(Log(x), Constant(log(base)))
}

# -- S3 Summary group handler -------------------------------------
# Dispatches base R Summary generics (sum, max, min, etc.) to CVXR atoms.
# Registered via registerS3method() in .onLoad().

.cvxr_Summary_handler <- function(..., na.rm = FALSE) {
  args <- list(...)
  nms <- names(args)
  ## CVXPY: sum/max/min/prod(x, axis, keepdims) reduce along axis
  ## (atoms/affine/sum.py:166-175, max.py:39, min.py:39, prod.py:141); any
  ## other keyword is a TypeError there. R's variadic convention differs: a
  ## name on a `...` entry labels an operand (sum(a = 1:3) is 6,
  ## do.call(sum, split(...))), and `na.rm` is the formal. So `axis` and
  ## `keepdims` are the options and every other entry, named or not, is an
  ## operand whose name is dropped. Recorded deviation, better-in-r.
  is_opt <- nms %in% c("axis", "keepdims") &
    .Generic %in% c("sum", "max", "min", "prod")
  if (any(is_opt)) {
    axis <- args[["axis"]]
    keepdims <- if (is.null(args[["keepdims"]])) FALSE else args[["keepdims"]]
    operands <- args[!is_opt]
    if (length(operands) != 1L) {
      cli_abort("{.arg axis} and {.arg keepdims} in {.fn {(.Generic)}} take a single expression.")
    }
    x <- operands[[1L]]
    return(switch(.Generic,
      "sum"  = SumEntries(x, axis = axis, keepdims = keepdims),
      "max"  = MaxEntries(x, axis = axis, keepdims = keepdims),
      "min"  = MinEntries(x, axis = axis, keepdims = keepdims),
      "prod" = Prod(x, axis = axis, keepdims = keepdims)
    ))
  }
  args <- unname(args)
  ## base R: sum(x, y, ...) is the sum of all the values present in all the
  ## arguments (?sum), max/min the largest/smallest of them (?max), prod their
  ## product; operands need not conform and the result has length 1. The
  ## elementwise forms are pmax/pmin, here max_elemwise/min_elemwise. So each
  ## operand is reduced on its own and the scalars are combined. (Until 1.9.3,
  ## sum(x, y) was sum(x + y) with broadcasting, prod(x, y) was prod(x * y),
  ## and max/min(x, y) were the elementwise Maximum/Minimum.)
  reduce <- switch(.Generic,
    "sum"   = SumEntries,
    "max"   = MaxEntries,
    "min"   = MinEntries,
    "prod"  = Prod,
    "range" = cli_abort("{.fn range} is not supported for CVXR expressions."),
    "any"   = cli_abort("{.fn any} is not applicable to CVXR expressions."),
    "all"   = cli_abort("{.fn all} is not applicable to CVXR expressions."),
    cli_abort("{.fn {(.Generic)}} is not supported for CVXR expressions.")
  )
  if (length(args) == 1L) {
    return(reduce(args[[1L]]))
  }
  parts <- lapply(args, function(a) reduce(as_expr(a)))
  switch(.Generic,
    "sum"   = Reduce(`+`, parts),
    "prod"  = Reduce(`*`, parts),
    "max"   = do.call(Maximum, parts),
    "min"   = do.call(Minimum, parts)
  )
}

# -- S3 Complex group handler -----------------------------------------
## Dispatches Re(), Im(), Conj() to CVXR atom constructors.
## R's Complex group generic includes: Arg, Conj, Im, Mod, Re.
## We support the three that map to CVXPY atoms.

.cvxr_Complex_handler <- function(z) {
  switch(.Generic,
    "Re"   = Real_(z),
    "Im"   = Imag_(z),
    "Conj" = Conj_(z),
    "Mod"  = Abs(z),
    "Arg"  = cli_abort("{.fn Arg} is not supported for CVXR expressions."),
    cli_abort("{.fn {(.Generic)}} is not supported for CVXR expressions.")
  )
}

# -- t() and mean() ----------------------------------------------------
## These are S3 generics (UseMethod) in base R, promoted to S4 by Matrix.
## S7 method() demotes them from S4 standardGeneric back to S3, breaking
## Matrix's S4 dispatch (t(dgeMatrix), mean(dgCMatrix) fail).
## Fix: use registerS3method() in .onLoad() instead — see globals.R.
## See conflicts/03_s7_method_demotes_s4.R and conflicts/04_registerS3method_safe.R.