R/059_atoms_elementwise_power.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: atoms/elementwise/power.py
## Power -- elementwise power x^p
##
## p is stored as a property (NOT as an arg). DCP curvature depends on p.
## Uses is_power2() from power_tools.R.


Power <- new_class("Power", parent = Elementwise, package = "CVXR",
  properties = list(
    p         = new_property(class = class_any),       # Constant or numeric
    p_used    = new_property(class = class_any),       # float value used for DCP
    max_denom = new_property(class = class_integer),
    p_orig    = new_property(class = class_any)        # original p value
  ),
  constructor = function(x, p, max_denom = 1024L, id = NULL) {
    if (is.null(id)) id <- next_expr_id()
    x <- as_expr(x)
    shape <- x@shape

    ## Store p as Constant if numeric
    p_orig <- p
    if (is.numeric(p) && !S7_inherits(p, Expression)) {
      p_const <- Constant(p)
    } else if (S7_inherits(p, Expression)) {
      p_const <- p
    } else {
      cli_abort("Exponent {.arg p} must be numeric or an Expression.")
    }

    ## Compute p_used (float value for DCP)
    if (S7_inherits(p_const, Constant)) {
      p_used <- as.numeric(value(p_const))
    } else {
      p_used <- NULL
    }

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

# -- sign ---------------------------------------------------------
## CVXPY: power.py lines 181-189
method(sign_from_args, Power) <- function(x) {
  pval <- x@p_used
  if (!is.null(pval) && pval == 1) {
    ## Same as input
    list(is_nonneg = is_nonneg(x@args[[1L]]),
         is_nonpos = is_nonpos(x@args[[1L]]))
  } else {
    ## Always nonneg
    list(is_nonneg = TRUE, is_nonpos = FALSE)
  }
}

# -- curvature ----------------------------------------------------
## CVXPY: power.py lines 191-206
method(is_atom_convex, Power) <- function(x) {
  pval <- x@p_used
  if (is.null(pval)) return(FALSE)  # param exponent: not DCP
  pval <= 0 || pval >= 1
}

method(is_atom_concave, Power) <- function(x) {
  pval <- x@p_used
  if (is.null(pval)) return(FALSE)
  pval >= 0 && pval <= 1
}

# -- monotonicity -------------------------------------------------
## CVXPY: power.py lines 258-290
method(is_incr, Power) <- function(x, idx, ...) {
  pval <- x@p_used
  if (is.null(pval)) return(FALSE)
  if (pval >= 0 && pval <= 1) return(TRUE)
  if (pval > 1) {
    if (is_power2(pval)) {
      return(is_nonneg(x@args[[idx]]))
    } else {
      return(TRUE)
    }
  }
  FALSE  # p < 0
}

method(is_decr, Power) <- function(x, idx, ...) {
  pval <- x@p_used
  if (is.null(pval)) return(FALSE)
  if (pval <= 0) return(TRUE)
  if (pval > 1 && is_power2(pval)) {
    return(is_nonpos(x@args[[idx]]))
  }
  FALSE
}

# -- is_constant: p == 0 makes it constant ------------------------
method(is_constant, Power) <- function(x) {
  pval <- x@p_used
  if (!is.null(pval) && pval == 0) return(TRUE)
  ## Fall through to default (all args constant)
  (0L %in% x@shape) || .all_args(x, is_constant)
}

# -- quadratic/PWL ------------------------------------------------
## CVXPY: power.py lines 292-337
method(is_quadratic, Power) <- function(x) {
  pval <- x@p_used
  if (is.null(pval)) return(FALSE)
  if (pval == 0) return(TRUE)
  if (pval == 1) return(is_quadratic(x@args[[1L]]))
  if (pval == 2) return(is_affine(x@args[[1L]]))
  is_constant(x@args[[1L]])
}

method(has_quadratic_term, Power) <- function(x) {
  pval <- x@p_used
  if (is.null(pval)) return(FALSE)
  if (pval == 1) return(has_quadratic_term(x@args[[1L]]))
  if (pval == 2) return(TRUE)
  FALSE
}

method(is_qpwa, Power) <- function(x) {
  pval <- x@p_used
  if (is.null(pval)) return(FALSE)
  if (pval == 0) return(TRUE)
  if (pval == 1) return(is_qpwa(x@args[[1L]]))
  if (pval == 2) return(is_pwl(x@args[[1L]]))
  is_constant(x@args[[1L]])
}

# -- domain -------------------------------------------------------
## CVXPY: power.py lines 376-390
method(atom_domain, Power) <- function(x) {
  pval <- x@p_used
  if (is.null(pval)) return(list())
  if ((pval < 1 && pval != 0) || (pval > 1 && !is_power2(pval))) {
    return(list(x@args[[1L]] >= 0))
  }
  list()
}

# -- get_data -----------------------------------------------------
method(get_data, Power) <- function(x) {
  list(x@p_orig, x@max_denom)
}

# -- name ---------------------------------------------------------
method(expr_name, Power) <- function(x) {
  pval <- if (!is.null(x@p_used)) x@p_used else "?"
  sprintf("Power(%s, %s)", expr_name(x@args[[1L]]), pval)
}

# -- numeric ------------------------------------------------------
method(numeric_value, Power) <- function(x, values, ...) {
  pval <- if (!is.null(x@p_used)) as.numeric(x@p_used) else as.numeric(value(x@p))
  values[[1L]]^pval
}

# -- graph_implementation: stub -----------------------------------
## CVXPY SOURCE: power.py lines 226-251
method(is_atom_log_log_convex, Power) <- function(x) {
  if (dpp_scope_active()) {
    ## DPP rule: power x^p is NOT log-log convex if BOTH x and p have parameters
    arg_x <- x@args[[1L]]
    p_expr <- x@p
    return(!(length(parameters(arg_x)) > 0L && length(parameters(p_expr)) > 0L))
  }
  TRUE
}

method(is_atom_log_log_concave, Power) <- function(x) {
  is_atom_log_log_convex(x)
}

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

# ===================================================================
# PowerApprox -- SOC-based rational approximation of Power
# ===================================================================
## CVXPY SOURCE: atoms/elementwise/power.py lines 421-449
## Subclass of Power. Overrides p_used with rational approximation
## and adds w (dyadic weights) for gm_constrs-based canonicalization.
## The factory function power() dispatches to PowerApprox when approx=TRUE.

PowerApprox <- new_class("PowerApprox", parent = Power, package = "CVXR",
  properties = list(
    w            = new_property(class = class_any),   # bigq weight vector or NULL
    approx_error = new_property(class = class_numeric)
  ),
  constructor = function(x, p, max_denom = 1024L, id = NULL) {
    if (is.null(id)) id <- next_expr_id()
    x <- as_expr(x)
    shape <- x@shape

    ## Store p as Constant if numeric (same as Power)
    p_orig <- p
    if (is.numeric(p) && !S7_inherits(p, Expression)) {
      p_const <- Constant(p)
    } else if (S7_inherits(p, Expression)) {
      p_const <- p
    } else {
      cli_abort("Exponent {.arg p} must be numeric or an Expression.")
    }

    ## Compute p_used (float value for DCP)
    if (S7_inherits(p_const, Constant)) {
      p_used <- as.numeric(value(p_const))
    } else {
      p_used <- NULL
    }

    ## Rational approximation -- override p_used and compute w
    w <- NULL
    approx_error <- 0.0
    if (!is.null(p_used)) {
      p_val <- p_used
      if (p_val > 1) {
        result <- pow_high(p_val, max_denom, approx = TRUE)
        p_used <- result[[1L]]
        w <- result[[2L]]
      } else if (p_val > 0 && p_val < 1) {
        result <- pow_mid(p_val, max_denom, approx = TRUE)
        p_used <- result[[1L]]
        w <- result[[2L]]
      } else if (p_val < 0) {
        result <- pow_neg(p_val, max_denom, approx = TRUE)
        p_used <- result[[1L]]
        w <- result[[2L]]
      }
      ## p == 0 or p == 1: no approximation needed, w stays NULL
      approx_error <- as.numeric(abs(as.numeric(p_used) - as.numeric(value(p_const))))
    }

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

# -- Factory function ---------------------------------------------
#' Create a Power atom
#'
#' @param x An Expression
#' @param p Numeric exponent
#' @param max_denom Maximum denominator for rational approximation
#' @param approx If TRUE (default), use SOC approximation. If FALSE, use exact power cone.
#' @returns A Power or PowerApprox atom
#' @note \code{sqrt(x)} on a CVXR expression dispatches to
#'   \code{Power(x, 0.5)} via the Math group generic.
#'   See \code{\link{math_atoms}} for all standard R function dispatch.
#' @export
power <- function(x, p, max_denom = 1024L, approx = TRUE) {
  if (approx) {
    PowerApprox(x, p, max_denom)
  } else {
    Power(x, p, max_denom)
  }
}