R/059_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 99-102) ---
method(sign_from_args, Convolve) <- function(x) {
  mul_sign(.args(x)[[1L]], .args(x)[[2L]])
}

# -- monotonicity (CVXPY conv.py lines 104-112) ------------------
## 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_atom_convex below, which
## both locate it), so these two must locate it as well -- exactly as Kron does
## with .kron_cst_loc 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}.
##
## Located by "has no variables" rather than is_constant() so the answer is the
## same inside the DPP scope, where a Parameter is not constant; is_atom_convex
## relies on that (decision 29).
.conv_cst_loc <- function(x) {
  if (length(variables(.args(x)[[1L]])) == 0L) 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)]])
}

## No log-log hooks: conv.py defines none, so Atom's FALSE applies
## (atom.py:195-202). A TRUE/TRUE pair here made conv(c, x) log-log affine
## and accepted GPs CVXPY rejects (sweep item 5, finding 5.1).

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)
  }
}

# -- curvature: affine, and under the DPP scope the DPP rule ----
## CVXPY SOURCE: atoms/affine/conv.py lines 72-86 (convolve: lines 173-187)
##     def is_atom_convex(self) -> bool:
##         # TODO support a non-constant first argument.
##         if u.scopes.dpp_scope_active():
##             # conv is not DPP if any parameters are present.
##             c = self.args[0]
##             return c.is_constant() and is_param_free(c)
##         else:
##             return self.args[0].is_constant()
## The C++ get_conv_mat handler reads the kernel through get_constant_data(),
## which cannot take a PARAM LinOp, so under the DPP scope the kernel must be
## parameter-free. Upstream tests args[0] because it only accepts the kernel
## there; CVXR locates it (.conv_cst_loc), the same rule, parity measured
## against CVXPY 1.9.3 (completeness TSV rows conv.is_atom_convex). Outside the
## scope the constructor (conv.R:11-14) has already enforced the `else` branch.
## Until 1.9.2.9163 this rule sat in a method(is_dpp, Convolve) reached by an
## argument recursion (decision 29).
method(is_atom_convex, Convolve) <- function(x) {
  if (dpp_scope_active()) {
    return(is_param_free(.args(x)[[.conv_cst_loc(x)]]))
  }
  TRUE
}
method(is_atom_concave, Convolve) <- function(x) is_atom_convex(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, the
#' numpy-style convolution
#' (\code{numpy.convolve(x, y) == stats::convolve(x, rev(y), type = "open")}).
#' If neither argument is an Expression the call is forwarded as written to
#' \code{\link[stats]{convolve}}, so the result is exactly \pkg{stats}'
#' (circular by default).
#'
#' @param x,y Expressions or numeric vectors; at least one must be constant
#'   when building an atom.
#' @param conj,type As in \code{\link[stats]{convolve}}; only for plain
#'   numeric input, an error with an Expression.
#' @returns A \code{Convolve} atom (for expressions) or a numeric vector.
#' @export
convolve <- function(x, y, conj = TRUE, type = c("circular", "open", "filter")) {
  if (inherits(x, "CVXR::Expression") || inherits(y, "CVXR::Expression")) {
    if (!missing(conj) || !missing(type))
      cli_abort("{.arg conj} and {.arg type} are {.fn stats::convolve} options, only for plain numeric input; an Expression gets the numpy-style convolution.")
    Convolve(x, y)
  } else {
    ## Forward the call as written: stats' defaults and argument matching apply.
    cl <- match.call()
    cl[[1L]] <- quote(stats::convolve)
    eval(cl, parent.frame())
  }
}