R/058_atoms_affine_conv.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: atoms/affine/conv.py
## Convolve -- 1D discrete convolution (one arg must be constant)


Convolve <- new_class("Convolve", parent = AffAtom, package = "CVXR",
  constructor = function(a, b, id = NULL) {
    if (FALSE) new_object(S7_object())  ## S7 static-check guard
    if (is.null(id)) id <- next_expr_id()
    a <- as_expr(a)
    b <- as_expr(b)
    ## At least one must be constant
    if (!is_constant(a) && !is_constant(b)) {
      cli_abort("At least one argument to {.fn conv} must be constant.")
    }
    ## Both must be vectors (n, 1)
    if (.shape(a)[2L] != 1L || .shape(b)[2L] != 1L) {
      cli_abort("{.fn conv} requires vector inputs (column vectors).")
    }
    ## Output length = len(a) + len(b) - 1
    out_len <- .shape(a)[1L] + .shape(b)[1L] - 1L
    shape <- c(out_len, 1L)

    obj <- .fast_new(Convolve, S7_object(),
      id    = as.integer(id),
      .cache = new.env(parent = emptyenv()),
      args  = list(a, b),
      shape = shape
    )
    obj
  }
)

method(shape_from_args, Convolve) <- function(x) {
  c(.arg_shape(x)[1L] + .arg_shape(x, 2L)[1L] - 1L, 1L)
}

# -- sign: same as multiplication (CVXPY conv.py lines 100-103) ---
method(sign_from_args, Convolve) <- function(x) {
  mul_sign(.args(x)[[1L]], .args(x)[[2L]])
}

# -- monotonicity (CVXPY conv.py lines 105-113) ------------------
## Monotonicity follows the sign of the CONSTANT KERNEL, whichever argument
## holds it.  Upstream reads args[0] unconditionally, which is safe there only
## because validate_arguments (conv.py:64-70) rejects a non-constant first
## argument.  CVXR deliberately accepts the constant in either position (see the
## constructor above, and graph_implementation / is_dpp below, which both locate
## it), so these two must locate it as well -- exactly as Kron does at
## kron.R:38-45 for the atom where upstream itself allows either order.
##
## Reading args[[1]] unconditionally made `conv(convex_expr, negative_kernel)`
## report is_incr = TRUE (the sign of the VARIABLE), so DCP called a concave
## composition convex and accepted it for minimization.  Reproducer:
## notes/audit/repro/01_conv_monotonicity.{py,R}.
.conv_cst_loc <- function(x) if (is_constant(.args(x)[[1L]])) 1L else 2L

method(is_incr, Convolve) <- function(x, idx, ...) {
  is_nonneg(.args(x)[[.conv_cst_loc(x)]])
}
method(is_decr, Convolve) <- function(x, idx, ...) {
  is_nonpos(.args(x)[[.conv_cst_loc(x)]])
}

# -- log-log: affine (CVXPY conv.py) -----------------------------
method(is_atom_log_log_convex, Convolve) <- function(x) TRUE
method(is_atom_log_log_concave, Convolve) <- function(x) TRUE

method(numeric_value, Convolve) <- function(x, values, ...) {
  a <- as.vector(values[[1L]])
  b <- as.vector(values[[2L]])
  ## R's convolve() uses Conj(FFT(y)), which gives cross-correlation
  ## instead of polynomial multiplication for complex inputs. Use FFT
  ## without conjugation for complex, R's convolve for real.
  if (is.complex(a) || is.complex(b)) {
    n <- length(a) + length(b) - 1L
    fa <- fft(c(a, rep(0+0i, n - length(a))))
    fb <- fft(c(b, rep(0+0i, n - length(b))))
    result <- fft(fa * fb, inverse = TRUE) / n
    ## Clean up near-zero imaginary parts if both inputs were real-ish
    matrix(result, ncol = 1L)
  } else {
    matrix(stats::convolve(a, rev(b), type = "open"), ncol = 1L)
  }
}

# -- DPP: conv is NOT DPP when the kernel arg is parametric ------
## Same as Kron -- C++ get_conv_mat uses get_constant_data and cannot
## handle PARAM LinOp nodes. Return FALSE when the kernel has parameters.
method(is_dpp, Convolve) <- function(x, context = "dcp") {
  cst_idx <- if (is_constant(.args(x)[[1L]])) 1L else 2L
  if (length(parameters(.args(x)[[cst_idx]])) > 0L) return(FALSE)
  if (identical(tolower(context), "dgp")) with_dpp_scope(is_dgp(x))
  else with_dpp_scope(is_dcp(x))
}

method(graph_implementation, Convolve) <- function(x, arg_objs, shape, data = NULL, ...) {
  if (is_constant(.args(x)[[1L]])) {
    list(conv_linop(arg_objs[[1L]], arg_objs[[2L]], shape), list())
  } else {
    list(conv_linop(arg_objs[[2L]], arg_objs[[1L]], shape), list())
  }
}

#' 1D discrete convolution
#'
#' @param a An Expression (vector, one must be constant)
#' @param b An Expression (vector)
#' @returns A Convolve atom
#' @export
conv <- function(a, b) {
  Convolve(a, b)
}

#' 1D discrete convolution (numpy-style)
#'
#' \code{convolve()} is CVXPY 1.9's preferred name for \code{\link{conv}}
#' (CVXPY's \code{conv} class is deprecated in favor of \code{convolve}). For a
#' CVXR \code{Expression} argument it builds a \code{Convolve} atom; for plain
#' numeric input it computes the same numpy-style convolution the atom does
#' (\code{numpy.convolve(a, b) == stats::convolve(a, rev(b), type = "open")}),
#' so the numeric and Expression paths agree.
#'
#' CVXR masks \code{stats::convolve} when attached (R reports this on load). If
#' you specifically want \pkg{stats}' circular cross-correlation or its other
#' options, call \code{\link[stats]{convolve}} directly.
#'
#' @param a,b Expressions or numeric vectors; at least one must be constant
#'   when building an atom.
#' @returns A \code{Convolve} atom (for expressions) or a numeric vector.
#' @export
convolve <- function(a, b) {
  if (inherits(a, "CVXR::Expression") || inherits(b, "CVXR::Expression")) {
    Convolve(a, b)
  } else {
    ## numpy-style ("open") convolution, identical to what the Convolve atom
    ## computes: numpy.convolve(a, b) == stats::convolve(a, rev(b), type = "open")
    stats::convolve(a, rev(b), type = "open")
  }
}