R/090_atoms_pnorm.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: atoms/pnorm.py
## Pnorm -- general p-norm, axis-aware


Pnorm <- new_class("Pnorm", parent = AxisAtom, package = "CVXR",
  properties = list(
    p         = new_property(class = class_numeric),
    max_denom = new_property(class = class_integer)
  ),
  constructor = function(x, p = 2, axis = NULL, keepdims = FALSE,
                         max_denom = 1024L, id = NULL) {
    if (is.null(id)) id <- next_expr_id()
    x <- as_expr(x)
    p <- as.numeric(p)
    if (!is.null(axis)) .validate_axis(axis, length(x@shape))
    shape <- .axis_shape(x@shape, axis, keepdims)

    obj <- new_object(S7_object(),
      id        = as.integer(id),
      .cache    = new.env(parent = emptyenv()),
      args      = list(x),
      shape     = shape,
      axis      = axis,
      keepdims  = keepdims,
      p         = p,
      max_denom = as.integer(max_denom)
    )
    validate_arguments(obj)
    obj
  }
)

# -- validate -----------------------------------------------------
method(validate_arguments, Pnorm) <- function(x) {
  if (x@p == 0) {
    cli_abort("{.fn p_norm} does not accept p=0.")
  }
  ## axis only supported for p=2 (CVXPY pnorm.py lines 157-158)
  if (!is.null(x@axis) && x@p != 2) {
    cli_abort("The {.arg axis} parameter is only supported for p=2.")
  }
  ## CVXPY: pnorm.py lines 160-161 -- complex not allowed for p < 1
  if (x@p < 1 && .any_args(x, is_complex)) {
    cli_abort("{.fn p_norm} cannot have complex {.arg x} for p < 1.")
  }
  invisible(NULL)
}

# -- sign: always nonneg ------------------------------------------
method(sign_from_args, Pnorm) <- function(x) {
  list(is_nonneg = TRUE, is_nonpos = FALSE)
}

# -- curvature ----------------------------------------------------
method(is_atom_convex, Pnorm) <- function(x) x@p >= 1
method(is_atom_concave, Pnorm) <- function(x) x@p < 1

# -- monotonicity -------------------------------------------------
method(is_incr, Pnorm) <- function(x, idx, ...) {
  x@p < 1 || (x@p >= 1 && is_nonneg(x@args[[1L]]))
}
method(is_decr, Pnorm) <- function(x, idx, ...) {
  x@p >= 1 && is_nonpos(x@args[[1L]])
}

# -- domain -------------------------------------------------------
method(atom_domain, Pnorm) <- function(x) {
  if (x@p < 1 && x@p != 0) {
    list(x@args[[1L]] >= 0)
  } else {
    list()
  }
}

# -- get_data -----------------------------------------------------
method(get_data, Pnorm) <- function(x) {
  list(x@p, x@axis, x@keepdims, x@max_denom)
}

# -- numeric ------------------------------------------------------
method(numeric_value, Pnorm) <- function(x, values, ...) {
  v <- values[[1L]]
  p <- x@p
  if (is.null(x@axis)) {
    ## Full p-norm -> scalar
    if (p == Inf) return(matrix(max(abs(v)), 1L, 1L))
    if (p == -Inf) return(matrix(min(abs(v)), 1L, 1L))
    matrix(sum(abs(v)^p)^(1/p), 1L, 1L)
  } else {
    ## Axis-wise p-norm
    margin <- if (x@axis == 2L) 2L else 1L
    res <- apply(v, margin, function(col) sum(abs(col)^p)^(1/p))
    if (x@keepdims) {
      if (x@axis == 2L) matrix(res, nrow = 1L) else matrix(res, ncol = 1L)
    } else {
      if (x@axis == 2L) matrix(res, nrow = 1L) else matrix(res, ncol = 1L)
    }
  }
}

# -- graph_implementation: stub -----------------------------------
## CVXPY SOURCE: pnorm.py lines 179-186
method(is_atom_log_log_convex, Pnorm) <- function(x) TRUE
method(is_atom_log_log_concave, Pnorm) <- function(x) FALSE

method(graph_implementation, Pnorm) <- function(x, arg_objs, shape, data = NULL, ...) {
  cli_abort("graph_implementation for {.cls Pnorm} not yet implemented.")
}

# ===================================================================
# PnormApprox -- SOC-based rational approximation of Pnorm
# ===================================================================
## CVXPY SOURCE: atoms/pnorm.py lines 272-289
## Subclass of Pnorm. Overrides p with rational approximation.
## The factory function p_norm() dispatches to PnormApprox when approx=TRUE.

PnormApprox <- new_class("PnormApprox", parent = Pnorm, package = "CVXR",
  properties = list(
    approx_error = new_property(class = class_numeric)
  ),
  constructor = function(x, p = 2, axis = NULL, keepdims = FALSE,
                         max_denom = 1024L, id = NULL) {
    if (is.null(id)) id <- next_expr_id()
    x <- as_expr(x)
    p_orig <- as.numeric(p)
    if (!is.null(axis)) .validate_axis(axis, length(x@shape))
    shape <- .axis_shape(x@shape, axis, keepdims)

    ## Rational approximation -- override p (convert to numeric for property)
    p_used <- p_orig
    if (p_orig < 0) {
      p_used <- as.numeric(pow_neg(p_orig, max_denom)[[1L]])
    } else if (p_orig > 0 && p_orig < 1) {
      p_used <- as.numeric(pow_mid(p_orig, max_denom)[[1L]])
    } else if (p_orig > 1) {
      p_used <- as.numeric(pow_high(p_orig, max_denom)[[1L]])
    }
    approx_error <- abs(p_used - p_orig)

    obj <- new_object(S7_object(),
      id           = as.integer(id),
      .cache       = new.env(parent = emptyenv()),
      args         = list(x),
      shape        = shape,
      axis         = axis,
      keepdims     = keepdims,
      p            = p_used,
      max_denom    = as.integer(max_denom),
      approx_error = approx_error
    )
    validate_arguments(obj)
    obj
  }
)

#' General p-norm of an expression
#'
#' @param x An Expression
#' @param p Numeric exponent (default 2)
#' @param axis NULL (all), 1 (row-wise), or 2 (column-wise)
#' @param keepdims Logical
#' @param max_denom Integer max denominator for rational approx
#' @param approx If TRUE (default), use SOC approximation. If FALSE, use exact power cone.
#' @returns A Pnorm, PnormApprox, Norm1, or NormInf atom
#' @export
p_norm <- function(x, p = 2, axis = NULL, keepdims = FALSE,
                   max_denom = 1024L, approx = TRUE) {
  ## Dispatch to specialized norms
  if (p == 1) return(Norm1(x, axis, keepdims))
  if (p == Inf) return(NormInf(x, axis, keepdims))
  if (approx) {
    PnormApprox(x, p, axis, keepdims, max_denom)
  } else {
    Pnorm(x, p, axis, keepdims, max_denom)
  }
}