R/045_atoms_affine_binary_operators.R

In CVXR: Disciplined Convex Optimization

#####
## DO NOT EDIT THIS FILE!! EDIT THE SOURCE INSTEAD: rsrc_tree/atoms/affine/binary_operators.R
#####

## CVXPY SOURCE: atoms/affine/binary_operators.py
## BinaryOperator, MulExpression, Multiply, DivExpression


# ===================================================================
# BinaryOperator -- base class for binary operations (other than add)
# ===================================================================

BinaryOperator <- new_class("BinaryOperator", parent = AffAtom, package = "CVXR",
  constructor = function(lh_exp, rh_exp, shape) {
    ## NEW_OBJECT GUARD: satisfy S7's check_S7_constructor static check.
    if (FALSE) new_object(S7_object())
    lh_exp <- as_expr(lh_exp)
    rh_exp <- as_expr(rh_exp)
    shape <- validate_shape(shape)
    obj <- .fast_new(BinaryOperator, S7_object(),
      id    = next_expr_id(),
      .cache = new.env(parent = emptyenv()),
      args  = list(lh_exp, rh_exp),
      shape = shape
    )
    validate_arguments(obj)
    obj
  }
)

# -- sign: multiply sign rules ---------------------------------------
## CVXPY SOURCE: binary_operators.py lines 92-95

method(sign_from_args, BinaryOperator) <- function(x) {
  mul_sign(.args(x)[[1L]], .args(x)[[2L]])
}

# -- Complex propagation ---------------------------------------------
## CVXPY SOURCE: binary_operators.py lines 97-107

method(is_imag, BinaryOperator) <- function(x) {
  (is_imag(.args(x)[[1L]]) && is_real(.args(x)[[2L]])) ||
    (is_real(.args(x)[[1L]]) && is_imag(.args(x)[[2L]]))
}

method(is_complex, BinaryOperator) <- function(x) {
  (is_complex(.args(x)[[1L]]) || is_complex(.args(x)[[2L]])) &&
    !(is_imag(.args(x)[[1L]]) && is_imag(.args(x)[[2L]]))
}


# ===================================================================
# MulExpression -- matrix multiplication (lhs %*% rhs)
# ===================================================================

MulExpression <- new_class("MulExpression", parent = BinaryOperator, package = "CVXR",
  constructor = function(lh_exp, rh_exp) {
    if (FALSE) new_object(S7_object())  ## S7 static-check guard
    lh_exp <- as_expr(lh_exp)
    rh_exp <- as_expr(rh_exp)
    shape <- mul_shapes(.shape(lh_exp), .shape(rh_exp))
    obj <- .fast_new(MulExpression, S7_object(),
      id    = next_expr_id(),
      .cache = new.env(parent = emptyenv()),
      args  = list(lh_exp, rh_exp),
      shape = shape
    )
    validate_arguments(obj)
    obj
  }
)

# -- shape_from_args -------------------------------------------------
## CVXPY SOURCE: binary_operators.py lines 188-191

method(shape_from_args, MulExpression) <- function(x) {
  mul_shapes(.arg_shape(x), .arg_shape(x, 2L))
}

# -- Convexity: requires one constant arg (with DPP extension) --------
## CVXPY SOURCE: binary_operators.py lines 193-214

method(is_atom_convex, MulExpression) <- function(x) {
  lhs <- .args(x)[[1L]]
  rhs <- .args(x)[[2L]]
  if (is_constant(lhs) || is_constant(rhs)) return(TRUE)
  ## DPP rule: product is DPP-convex if one arg is param-affine
  ## and the other is parameter-free.
  if (dpp_scope_active()) {
    return((is_param_affine(lhs) && is_param_free(rhs)) ||
           (is_param_affine(rhs) && is_param_free(lhs)))
  }
  FALSE
}

method(is_atom_concave, MulExpression) <- function(x) {
  is_atom_convex(x)
}

# -- Monotonicity ----------------------------------------------------
## CVXPY SOURCE: binary_operators.py lines 231-239
## idx is 1-based (R convention): is_incr(1) checks args[[2]], is_incr(2) checks args[[1]]

method(is_incr, MulExpression) <- function(x, idx, ...) {
  ## self.args[1-idx] in Python (0-based) -> args[[3L - idx]] in R (1-based)
  is_nonneg(.args(x)[[3L - idx]])
}

method(is_decr, MulExpression) <- function(x, idx, ...) {
  is_nonpos(.args(x)[[3L - idx]])
}

# -- numeric_value ---------------------------------------------------
## CVXPY SOURCE: binary_operators.py lines 180-186

method(numeric_value, MulExpression) <- function(x, values, ...) {
  lhs <- values[[1L]]
  rhs <- values[[2L]]
  ## Handle scalar multiplication
  if (length(lhs) == 1L) return(drop(lhs) * rhs)
  if (length(rhs) == 1L) return(lhs * drop(rhs))
  ## Matrix multiplication
  lhs %*% rhs
}

# -- graph_implementation --------------------------------------------
## CVXPY SOURCE: binary_operators.py lines 295-323

method(graph_implementation, MulExpression) <- function(x, arg_objs, shape, data = NULL, ...) {
  lhs <- arg_objs[[1L]]
  rhs <- arg_objs[[2L]]
  if (is_constant(.args(x)[[1L]])) {
    list(mul_expr_linop(lhs, rhs, shape), list())
  } else if (is_constant(.args(x)[[2L]])) {
    list(rmul_expr_linop(lhs, rhs, shape), list())
  } else {
    cli_abort("Product of two non-constant expressions is not DCP.", class = "DCPError")
  }
}

# -- .grad: matmul Jacobian via Kronecker products ----------------
## CVXPY SOURCE: binary_operators.py:241-294 (MulExpression._grad).
## For C = X %*% Y with X (m, k) and Y (k, n):
##   grad_X = kron(Y, I_m)         shape (m*k, m*n)
##   grad_Y = kron(I_n, X)^T       shape (k*n, m*n)
## (CVXPY convention: grad[i, j] = d(output[j]) / d(input[i]); F-order.)
##
## CVXPY's _grad branches on `is_constant` and falls back to AffAtom's
## canon-interface coefficient extraction in that case. Until our
## .grad(AffAtom) lands, we use the explicit Kronecker formulas
## unconditionally — same math, works for both constant and
## non-constant args; the chain rule walker handles either side.
method(.grad, MulExpression) <- function(x, values, ...) {
  X <- values[[1L]]; Y <- values[[2L]]
  x_shape <- .arg_shape(x)
  y_shape <- .arg_shape(x, 2L)

  ## Promote to 2D matrices for consistent Kronecker computation.
  ## CVXR is 2D-only so X, Y are already matrices in the expected sense.
  X_mat <- as.matrix(X)
  Y_mat <- as.matrix(Y)

  m <- as.integer(x_shape[1L])
  n <- as.integer(y_shape[2L])

  I_m <- Matrix::Diagonal(m)
  I_n <- Matrix::Diagonal(n)

  ## DX = kron(Y, I_m); DY = kron(I_n, X)^T.
  ## Use Matrix::kronecker (operator %x%) for sparse output.
  DX <- methods::kronecker(Y_mat, I_m, FUN = "*", make.dimnames = FALSE)
  DY <- methods::kronecker(I_n,   X_mat, FUN = "*", make.dimnames = FALSE)
  DY <- t(DY)

  ## Coerce to CSC.
  list(as(DX, "CsparseMatrix"), as(DY, "CsparseMatrix"))
}

# -- expr_name -------------------------------------------------------

## CVXPY SOURCE: binary_operators.py lines 221-229
method(is_atom_log_log_convex, MulExpression) <- function(x) TRUE
method(is_atom_log_log_concave, MulExpression) <- function(x) FALSE

method(expr_name, MulExpression) <- function(x) {
  .binop_name(x, "%*%")
}


# ===================================================================
# Multiply -- elementwise multiplication (lhs * rhs)
# ===================================================================

Multiply <- new_class("Multiply", parent = MulExpression, package = "CVXR",
  constructor = function(lh_exp, rh_exp) {
    if (FALSE) new_object(S7_object())  ## S7 static-check guard
    lh_exp <- as_expr(lh_exp)
    rh_exp <- as_expr(rh_exp)
    ## Broadcast scalars
    bcast <- broadcast_args(lh_exp, rh_exp)
    lh_exp <- bcast[[1L]]
    rh_exp <- bcast[[2L]]
    ## Shape from broadcasting
    shape <- sum_shapes(list(.shape(lh_exp), .shape(rh_exp)))
    obj <- .fast_new(Multiply, S7_object(),
      id    = next_expr_id(),
      .cache = new.env(parent = emptyenv()),
      args  = list(lh_exp, rh_exp),
      shape = shape
    )
    ## Multiply allows complex (inherits AffAtom validate_arguments)
    obj
  }
)

# -- shape_from_args -------------------------------------------------
## CVXPY SOURCE: binary_operators.py lines 371-373

method(shape_from_args, Multiply) <- function(x) {
  sum_shapes(list(.arg_shape(x), .arg_shape(x, 2L)))
}

# -- validate_arguments ----------------------------------------------
## CVXPY SOURCE: binary_operators.py lines 364-369
## Broadcast compatibility already checked by sum_shapes in constructor.

method(validate_arguments, Multiply) <- function(x) {
  invisible(NULL)
}

# -- numeric_value: elementwise --------------------------------------
## CVXPY SOURCE: binary_operators.py lines 355-362

method(numeric_value, Multiply) <- function(x, values, ...) {
  lhs <- values[[1L]]
  rhs <- values[[2L]]
  ## R's * handles sparse matrices correctly
  lhs * rhs
}

# -- symmetry --------------------------------------------------------
## CVXPY v1.9.0 #3142: elementwise multiply of symmetric matrices is
## symmetric. Needed so quad_form(x, 2*P) accepts 2*P (a Multiply) as
## param-affine + symmetric when P is a PSD/symmetric Parameter.
## CVXPY SOURCE: binary_operators.py multiply.is_symmetric (lines 403-406).
method(is_symmetric, Multiply) <- function(x) {
  .shape(x)[1L] == .shape(x)[2L] &&
    all(vapply(.args(x), is_symmetric, logical(1L)))
}

# -- PSD/NSD ---------------------------------------------------------
## CVXPY SOURCE: binary_operators.py lines 375-385

method(is_psd, Multiply) <- function(x) {
  (is_psd(.args(x)[[1L]]) && is_psd(.args(x)[[2L]])) ||
    (is_nsd(.args(x)[[1L]]) && is_nsd(.args(x)[[2L]]))
}

method(is_nsd, Multiply) <- function(x) {
  (is_psd(.args(x)[[1L]]) && is_nsd(.args(x)[[2L]])) ||
    (is_nsd(.args(x)[[1L]]) && is_psd(.args(x)[[2L]]))
}

## CVXPY v1.8.2: is_hermitian for elementwise multiply
## Enables trace(A@B) optimization to detect Hermitian products.
method(is_hermitian, Multiply) <- function(x) {
  is_hermitian(.args(x)[[1L]]) && is_hermitian(.args(x)[[2L]])
}

# -- Quasiconvexity --------------------------------------------------
## CVXPY SOURCE: binary_operators.py lines 343-353
method(is_atom_quasiconvex, Multiply) <- function(x) {
  (is_constant(.args(x)[[1L]]) || is_constant(.args(x)[[2L]])) ||
    (is_nonneg(.args(x)[[1L]]) && is_nonpos(.args(x)[[2L]])) ||
    (is_nonpos(.args(x)[[1L]]) && is_nonneg(.args(x)[[2L]]))
}

method(is_atom_quasiconcave, Multiply) <- function(x) {
  (is_constant(.args(x)[[1L]]) || is_constant(.args(x)[[2L]])) ||
    (is_nonneg(.args(x)[[1L]]) && is_nonneg(.args(x)[[2L]])) ||
    (is_nonpos(.args(x)[[1L]]) && is_nonpos(.args(x)[[2L]]))
}

# -- graph_implementation --------------------------------------------
## CVXPY SOURCE: binary_operators.py lines 416-444

method(graph_implementation, Multiply) <- function(x, arg_objs, shape, data = NULL, ...) {
  lhs <- arg_objs[[1L]]
  rhs <- arg_objs[[2L]]
  if (is_constant(.args(x)[[1L]])) {
    list(multiply_linop(lhs, rhs), list())
  } else if (is_constant(.args(x)[[2L]])) {
    list(multiply_linop(rhs, lhs), list())
  } else {
    cli_abort("Product of two non-constant expressions is not DCP.", class = "DCPError")
  }
}

# -- .grad: elementwise product diagonal Jacobians ----------------
## CVXPY SOURCE: binary_operators.py:392-419 (multiply._grad).
## For Z = multiply(X, Y), Z[i] = X[i] * Y[i] elementwise:
##   grad_X = diag(Y_flat),  grad_Y = diag(X_flat)
## Both flat in column-major (F-order). Same math whether or not one
## arg is constant; we use the explicit form unconditionally.
method(.grad, Multiply) <- function(x, values, ...) {
  x_flat <- as.numeric(values[[1L]])
  y_flat <- as.numeric(values[[2L]])
  ## Diagonals must have the *same* length on both sides — they're
  ## defined over the (broadcast) output. CVXPY relies on numpy's
  ## broadcasting in self.numeric having already produced matching
  ## flat lengths; here CVXR's matrix shapes are assumed pre-aligned
  ## by the constructor.
  list(
    Matrix::sparseMatrix(
      i = seq_along(y_flat), j = seq_along(y_flat),
      x = y_flat, dims = c(length(y_flat), length(y_flat)),
      repr = "C"
    ),
    Matrix::sparseMatrix(
      i = seq_along(x_flat), j = seq_along(x_flat),
      x = x_flat, dims = c(length(x_flat), length(x_flat)),
      repr = "C"
    )
  )
}

# -- expr_name -------------------------------------------------------

## CVXPY SOURCE: binary_operators.py lines 335-341
method(is_atom_log_log_convex, Multiply) <- function(x) TRUE
method(is_atom_log_log_concave, Multiply) <- function(x) TRUE

method(expr_name, Multiply) <- function(x) {
  .binop_name(x, "*")
}


# ===================================================================
# DivExpression -- division (lhs / rhs)
# ===================================================================

DivExpression <- new_class("DivExpression", parent = BinaryOperator, package = "CVXR",
  constructor = function(lh_exp, rh_exp) {
    if (FALSE) new_object(S7_object())  ## S7 static-check guard
    lh_exp <- as_expr(lh_exp)
    rh_exp <- as_expr(rh_exp)
    ## Broadcast scalars
    bcast <- broadcast_args(lh_exp, rh_exp)
    lh_exp <- bcast[[1L]]
    rh_exp <- bcast[[2L]]
    ## Shape is numerator shape
    shape <- .shape(lh_exp)
    obj <- .fast_new(DivExpression, S7_object(),
      id    = next_expr_id(),
      .cache = new.env(parent = emptyenv()),
      args  = list(lh_exp, rh_exp),
      shape = shape
    )
    validate_arguments(obj)
    obj
  }
)

# -- shape_from_args -------------------------------------------------
## CVXPY SOURCE: binary_operators.py lines 478-481

method(shape_from_args, DivExpression) <- function(x) .arg_shape(x)

# -- Convexity: requires constant denominator ------------------------
## CVXPY SOURCE: binary_operators.py lines 483-490

method(is_atom_convex, DivExpression) <- function(x) {
  is_constant(.args(x)[[2L]])
}

method(is_atom_concave, DivExpression) <- function(x) {
  is_atom_convex(x)
}

# -- Quasiconvexity --------------------------------------------------
## CVXPY SOURCE: binary_operators.py lines 502-506
method(is_atom_quasiconvex, DivExpression) <- function(x) {
  is_nonneg(.args(x)[[2L]]) || is_nonpos(.args(x)[[2L]])
}

method(is_atom_quasiconcave, DivExpression) <- function(x) {
  is_atom_quasiconvex(x)
}

# -- Monotonicity ----------------------------------------------------
## CVXPY SOURCE: binary_operators.py lines 508-522
## idx is 1-based (R convention)

method(is_incr, DivExpression) <- function(x, idx, ...) {
  if (idx == 1L) {
    ## d(lhs/rhs)/d(lhs) > 0 when rhs > 0
    is_nonneg(.args(x)[[2L]])
  } else {
    ## d(lhs/rhs)/d(rhs) > 0 when lhs < 0 (nonpositive)
    is_nonpos(.args(x)[[1L]])
  }
}

method(is_decr, DivExpression) <- function(x, idx, ...) {
  if (idx == 1L) {
    ## d(lhs/rhs)/d(lhs) < 0 when rhs < 0
    is_nonpos(.args(x)[[2L]])
  } else {
    ## d(lhs/rhs)/d(rhs) < 0 when lhs > 0 (nonneg)
    is_nonneg(.args(x)[[1L]])
  }
}

# -- numeric_value ---------------------------------------------------
## CVXPY SOURCE: binary_operators.py lines 460-466

method(numeric_value, DivExpression) <- function(x, values, ...) {
  lhs <- if (inherits(values[[1L]], "sparseMatrix")) as.matrix(values[[1L]]) else values[[1L]]
  rhs <- if (inherits(values[[2L]], "sparseMatrix")) as.matrix(values[[2L]]) else values[[2L]]
  lhs / rhs
}

# -- graph_implementation --------------------------------------------
## CVXPY SOURCE: binary_operators.py lines 524-543

method(graph_implementation, DivExpression) <- function(x, arg_objs, shape, data = NULL, ...) {
  list(div_expr_linop(arg_objs[[1L]], arg_objs[[2L]]), list())
}

# -- expr_name -------------------------------------------------------

## CVXPY SOURCE: binary_operators.py lines 492-500
method(is_atom_log_log_convex, DivExpression) <- function(x) TRUE
method(is_atom_log_log_concave, DivExpression) <- function(x) TRUE

method(expr_name, DivExpression) <- function(x) {
  .binop_name(x, "/")
}


# ===================================================================
# Helper for binary operator names
# ===================================================================

.binop_name <- function(x, op_name) {
  lhs_name <- expr_name(.args(x)[[1L]])
  rhs_name <- expr_name(.args(x)[[2L]])
  ## Parenthesize AddExpression and DivExpression args
  if (.s7_is(.args(x)[[1L]], AddExpression) ||
      .s7_is(.args(x)[[1L]], DivExpression)) {
    lhs_name <- sprintf("(%s)", lhs_name)
  }
  if (.s7_is(.args(x)[[2L]], AddExpression) ||
      .s7_is(.args(x)[[2L]], DivExpression)) {
    rhs_name <- sprintf("(%s)", rhs_name)
  }
  sprintf("%s %s %s", lhs_name, op_name, rhs_name)
}

## Label-aware twin of .binop_name; mirrors precedence rules but
## recurses with format_labeled on each operand.
.binop_format_labeled <- function(x, op_name) {
  lhs_name <- format_labeled(.args(x)[[1L]])
  rhs_name <- format_labeled(.args(x)[[2L]])
  if (.s7_is(.args(x)[[1L]], AddExpression) ||
      .s7_is(.args(x)[[1L]], DivExpression)) {
    lhs_name <- sprintf("(%s)", lhs_name)
  }
  if (.s7_is(.args(x)[[2L]], AddExpression) ||
      .s7_is(.args(x)[[2L]], DivExpression)) {
    rhs_name <- sprintf("(%s)", rhs_name)
  }
  sprintf("%s %s %s", lhs_name, op_name, rhs_name)
}

## CVXPY SOURCE: atoms/affine/binary_operators.py BinaryOperator.format_labeled.
method(format_labeled, MulExpression) <- function(x) {
  lbl <- label(x); if (!is.null(lbl)) return(lbl)
  .binop_format_labeled(x, "%*%")
}
method(format_labeled, Multiply) <- function(x) {
  lbl <- label(x); if (!is.null(lbl)) return(lbl)
  .binop_format_labeled(x, "*")
}
method(format_labeled, DivExpression) <- function(x) {
  lbl <- label(x); if (!is.null(lbl)) return(lbl)
  .binop_format_labeled(x, "/")
}

# -- bounds_from_args (#3080) -----------------------------------------
## Defined after the classes above so each is in scope at registration.

## CVXPY SOURCE: binary_operators.py:237-247 (MulExpression). Scalar operand ->
## elementwise mul_bounds; otherwise matrix-product matmul_bounds.
method(bounds_from_args, MulExpression) <- function(x) {
  b1 <- get_bounds(.args(x)[[1L]]); b2 <- get_bounds(.args(x)[[2L]])
  if (length(b1[[1L]]) == 1L || length(b2[[1L]]) == 1L) {
    mul_bounds(b1[[1L]], b1[[2L]], b2[[1L]], b2[[2L]])
  } else {
    matmul_bounds(b1[[1L]], b1[[2L]], b2[[1L]], b2[[2L]])
  }
}

## CVXPY SOURCE: binary_operators.py:361-369 (multiply -- always elementwise).
method(bounds_from_args, Multiply) <- function(x) {
  b1 <- get_bounds(.args(x)[[1L]]); b2 <- get_bounds(.args(x)[[2L]])
  mul_bounds(b1[[1L]], b1[[2L]], b2[[1L]], b2[[2L]])
}

## CVXPY SOURCE: binary_operators.py:546-550 (DivExpression).
method(bounds_from_args, DivExpression) <- function(x) {
  b1 <- get_bounds(.args(x)[[1L]]); b2 <- get_bounds(.args(x)[[2L]])
  div_bounds(b1[[1L]], b1[[2L]], b2[[1L]], b2[[2L]])
}

# =====================================================================
# Module-level products: vdot / scalar_product / outer
# =====================================================================
## These lived in zzz_R_specific/convenience.R until 1.9.1.9040, under a header
## claiming they had "no direct CVXPY counterpart". That was false -- all three
## are defined right here in binary_operators.py -- and it is very likely why
## vdot's missing conjugation went unreviewed for the life of the project.
## Moved to their real home (constraint 15f: zzz_R_specific is for files with NO
## counterpart). Found by the completeness audit; see
## notes/audit/completeness_ledger.md finding 2.

#' Vector dot product (inner product)
#'
#' @description
#' The standard inner product of `x` and `y`: both are flattened, multiplied
#' elementwise, and summed. **Conjugate-linear in `x`** — `x` is conjugated
#' before multiplying, matching `numpy.vdot()` and CVXPY's `vdot()`. For real
#' arguments conjugation is the identity, so this is the ordinary dot product.
#'
#' Either argument may also be a (possibly nested) list, which is flattened
#' before the product; the two are flattened independently, so their nesting
#' need not match.
#'
#' @param x An Expression, numeric value, or nested list thereof. The
#'   conjugate-linear argument.
#' @param y An Expression, numeric value, or nested list thereof. The linear
#'   argument.
#' @returns A scalar Expression representing `sum(Conj(x) * y)`.
#' @seealso [scalar_product()], [cvxr_outer()], [deep_flatten()]
#' @examples
#' x <- Variable(3)
#' vdot(x, c(1, 2, 3))
#'
#' a <- Variable(); b <- Variable()
#' vdot(list(a, b), c(1, 2))
#' @export
vdot <- function(x, y) {
  ## CVXPY SOURCE: binary_operators.py:616-619
  ##     x = deep_flatten(x); y = deep_flatten(y)
  ##     prod = multiply(conj(x), y)
  ##     return cvxpy_sum(prod)
  ##
  ## The Conj_ wrap is UNCONDITIONAL, as upstream's is. It is tempting to skip
  ## it when is_real(x) -- expression.R:145 uses exactly that short-circuit for
  ## the conjugate transpose -- but measured against CVXPY 1.9.2, conj() is
  ## non-monotone (conj.py:41-49) and so drops a real CONVEX argument to UNKNOWN
  ## curvature: cp.vdot(cp.square(z), c) is not DCP upstream. Short-circuiting
  ## would leave CVXR strictly more permissive than CVXPY, which is the exact
  ## shape of the conv() defect this audit just fixed. Divergence here would
  ## need to be a deliberate, recorded decision, not a side effect.
  ##
  ## Affine arguments -- the overwhelmingly common case -- are unaffected: an
  ## AffAtom over an affine argument stays affine regardless of monotonicity.
  ##
  ## Both arguments go through deep_flatten (reshape.R), so a nested list works
  ## as it does upstream: vdot(list(a, b), c(1, 2)) is a * 1 + b * 2. The two
  ## arguments are flattened INDEPENDENTLY -- their nesting need not match.
  sum_entries(Conj_(deep_flatten(x)) * deep_flatten(y))
}

#' Scalar product (alias for vdot)
#'
#' @inheritParams vdot
#' @returns A scalar Expression representing `sum(Conj(x) * y)`.
#' @seealso [vdot()]
#' @export
scalar_product <- function(x, y) {
  ## CVXPY SOURCE: binary_operators.py:622-626
  vdot(x, y)
}

#' Outer product of two vectors
#'
#' @description
#' Computes the outer product `x %*% t(y)`. Both inputs must be vectors.
#' Named `cvxr_outer()` because base R already exports `outer()`.
#'
#' @param x An Expression or numeric value (vector).
#' @param y An Expression or numeric value (vector).
#' @returns An Expression of shape (length(x), length(y)).
#' @seealso [vdot()]
#' @export
cvxr_outer <- function(x, y) {
  ## CVXPY SOURCE: binary_operators.py:629-660 (`outer`; renamed here because
  ## base::outer exists -- constraint 15e).
  ##
  ## NOTE: upstream's docstring says the inputs may be a "nested list thereof"
  ## and are "flattened if not already a vector". Its BODY does neither -- it
  ## calls cast_to_const and then rejects anything with ndim > 1. Measured
  ## against cvxpy 1.9.2: outer([[1],[2]], [3,4]) and outer(<2x2 matrix>, ...)
  ## both raise "x must be a 1-d array." So the vector check below is the
  ## faithful port, and deep_flatten deliberately is NOT used here. Following
  ## the docstring instead of the behavior would have put a divergence in.
  x <- as_expr(x)
  y <- as_expr(y)
  ## Validate: both must be vectors (one dimension == 1)
  if (min(x@shape) != 1L) {
    cli_abort("{.fn cvxr_outer}: {.arg x} must be a vector, got shape {x@shape[1]}x{x@shape[2]}.")
  }
  if (min(y@shape) != 1L) {
    cli_abort("{.fn cvxr_outer}: {.arg y} must be a vector, got shape {y@shape[1]}x{y@shape[2]}.")
  }
  ## Reshape to column and row vectors, then matrix multiply
  n <- expr_size(x)
  m <- expr_size(y)
  x_col <- reshape_expr(x, c(n, 1L))
  y_row <- reshape_expr(y, c(1L, m))
  x_col %*% y_row
}