R/046_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 247-250

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 252-273

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 302-310
## 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 239-245

method(numeric_value, MulExpression) <- function(x, values, ...) {
  lhs <- values[[1L]]
  rhs <- values[[2L]]
  ## Scalar multiplication (CVXPY: a `()` shape; an R scalar value is 1 x 1,
  ## ADR D_AP.3). as.vector(), not drop(): base::drop leaves a 1 x 1 Matrix
  ## object as it is, since Matrix's drop is an S4 generic this namespace
  ## does not see. A dense Matrix result becomes a base matrix, as scipy
  ## returns an ndarray; sparse stays sparse.
  if (length(lhs) == 1L) return(.dense_Matrix_to_base(as.vector(lhs) * rhs))
  if (length(rhs) == 1L) return(.dense_Matrix_to_base(lhs * as.vector(rhs)))
  ## Matrix multiplication
  .dense_Matrix_to_base(lhs %*% rhs)
}

# -- graph_implementation --------------------------------------------
## CVXPY SOURCE: binary_operators.py lines 366-394

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:312-364 (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` (binary_operators.py:329-330) and
## falls back to AffAtom._grad's canon-interface coefficient extraction when
## either argument is constant, using the Kronecker formulas only for the
## bilinear case. CVXR uses the Kronecker formulas unconditionally: for a
## constant argument they are the exact Jacobian of the remaining linear map,
## and the constant argument's own block is dropped by the chain rule in
## grad(Atom) (a Constant has no variables), so the values agree with upstream
## on every input. .grad(AffAtom) exists (affine_atom.R) and is not needed
## here. Recorded as claims audit C0025 (2026-10-11, comment-only).
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 280-288
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 452-454

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 445-450
## 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 436-443

method(numeric_value, Multiply) <- function(x, values, ...) {
  lhs <- values[[1L]]
  rhs <- values[[2L]]
  ## Matrix's `*` handles sparse operands; sparse * dense is a dgeMatrix,
  ## which becomes a base matrix (the value path keeps base or dgCMatrix),
  ## sparse * sparse stays sparse.
  .dense_Matrix_to_base(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 456-459).
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 461-471

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 424-434
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 507-535

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:478-505 (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 406-412
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)
    orig_shapes <- list(.shape(lh_exp), .shape(rh_exp))
    ## Broadcast scalars
    bcast <- broadcast_args(lh_exp, rh_exp)
    lh_exp <- bcast[[1L]]
    rh_exp <- bcast[[2L]]
    ## CVXPY Expression.__div__ (expression.py:862-870): after broadcasting,
    ## the operands must be a scalar or of equal shape. CVXPY checks this in
    ## the operator; here the constructor is the single entry point for `/`,
    ## log2/log10 and the DQCP inverse, so the check lives here.
    ls <- .shape(lh_exp)
    rs <- .shape(rh_exp)
    if (!(expr_is_scalar(lh_exp) || expr_is_scalar(rh_exp) || identical(ls, rs))) {
      hint <- if (!identical(orig_shapes, list(ls, rs))) {
        c(i = paste0("The operands were (", orig_shapes[[1L]][1L], ", ",
                     orig_shapes[[1L]][2L], ") and (", orig_shapes[[2L]][1L],
                     ", ", orig_shapes[[2L]][2L], ") before broadcasting."))
      }
      cli_abort(c(
        "Incompatible shapes for division (({ls[1L]}, {ls[2L]}) / ({rs[1L]}, {rs[2L]}))",
        hint
      ))
    }
    ## Shape is numerator shape
    shape <- ls
    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 569-572

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

# -- Quadratic predicates: require a constant divisor -----------------
## CVXPY SOURCE: binary_operators.py lines 559-567

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

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

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

# -- Convexity: requires constant denominator ------------------------
## CVXPY SOURCE: binary_operators.py lines 574-581

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 593-597
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 605-619
## 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 551-557

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 621-640

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 583-591
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:290-300 (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:414-422 (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:599-603 (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:669-672
  ##     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:138 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:675-679
  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:682-709 (`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
}