R/066_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
## At 1.9.0: power() factory const-base / variable-exponent dispatch
## (#3180 _pow_const_base: power(b, x) and b^x for positive constant b =>
## exp(x*log(b))), bounds_from_args (power_bounds), and is_atom_smooth all ported.
## 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 (FALSE) new_object(S7_object())  ## S7 static-check guard
    if (is.null(id)) id <- next_expr_id()
    x <- as_expr(x)
    shape <- .shape(x)

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

    ## CVXPY SOURCE: power.py Power.__init__ -- the exponent must be a Constant
    ## or a Parameter (e.g. y^x with a Variable exponent is rejected; the
    ## positive-constant-base / variable-exponent case is handled upstream by
    ## the power() factory via .pow_const_base).
    if (!(.s7_is(p_const, Constant) || .s7_is(p_const, Parameter))) {
      cli_abort("The exponent {.arg p} must be either a Constant or a Parameter.")
    }

    ## Compute p_used (float value for DCP)
    if (.s7_is(p_const, Constant)) {
      p_used <- as.numeric(value(p_const))
      ## power.py:183 float(self.p.value) rejects a non-scalar exponent.
      if (length(p_used) != 1L)
        cli_abort("The exponent must be a scalar, not of length {length(p_used)}.")
    } else {
      p_used <- NULL
    }

    obj <- .fast_new(Power, 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
  }
)

# -- parameters ---------------------------------------------------
## CVXPY: power.py:233-249 — override needed because `p` is stored as
## a property, not in `args`, so the default Canonical traversal misses
## any Parameter held in `p` (e.g., a parametric exponent).
method(parameters, Power) <- function(x) {
  unique_list(c(
    unlist(lapply(.args(x), parameters), recursive = FALSE),
    parameters(x@p)
  ))
}

# -- bounds: x^p (#3080) ------------------------------------------
## CVXPY SOURCE: elementwise/power.py:206-210.
## `p = float(self.p.value) if self.p.value is not None else 1.0`: the exponent
## is the ORIGINAL `p` (`x@p`, a Constant or a Parameter), read at call time,
## so a Parameter exponent uses its current value and only an UNSET Parameter
## falls back to 1.0.  Not `p_used`: for PowerApprox that is the rational
## approximation, which CVXPY uses only for monotonicity, quadratic/PWL and
## the gradient (claims audit C0187, 1.9.2.9158).  as.numeric() mirrors the
## `float()` cast; power_bounds() is typed `p: float`.
method(bounds_from_args, Power) <- function(x) {
  b <- get_bounds(.args(x)[[1L]])
  pv <- value(x@p)
  p <- if (is.null(pv)) 1.0 else as.numeric(pv)
  power_bounds(b[[1L]], b[[2L]], p)
}

# -- sign ---------------------------------------------------------
## CVXPY SOURCE: power.py:196-204 -- tests the exponent as written
## (`self.p.value == 1`; a Parameter's current value), not p_used: a
## PowerApprox whose p rounds to 1 is still nonneg.
method(sign_from_args, Power) <- function(x) {
  pval <- value(x@p)
  if (!is.null(pval) && as.numeric(pval) == 1) {
    ## Same as input
    list(is_nonneg = is_nonneg(.args(x)[[1L]]),
         is_nonpos = is_nonpos(.args(x)[[1L]]))
  } else {
    ## Always nonneg
    list(is_nonneg = TRUE, is_nonpos = FALSE)
  }
}

# -- curvature ----------------------------------------------------
## CVXPY SOURCE: power.py:212-227 -- `_is_const(self.p)` and the exponent as
## written (`self.p.value`), not p_used.
method(is_atom_convex, Power) <- function(x) {
  if (!.s7_is(x@p, Constant)) return(FALSE)  # param exponent: not DCP
  pval <- as.numeric(value(x@p))
  pval <= 0 || pval >= 1
}

method(is_atom_concave, Power) <- function(x) {
  if (!.s7_is(x@p, Constant)) return(FALSE)
  pval <- as.numeric(value(x@p))
  pval >= 0 && pval <= 1
}

## CVXPY power.py:229-231: smooth iff the exponent is a constant (p_used set);
## a parametric exponent (p_used NULL) is not smooth. Mirrors _is_const(self.p).
method(is_atom_smooth, Power) <- function(x) !is.null(x@p_used)

# -- monotonicity -------------------------------------------------
## CVXPY SOURCE: power.py:283-315 -- a Parameter exponent is monotone by its
## sign (286-287, 303-304); a constant one by p_used.
method(is_incr, Power) <- function(x, idx, ...) {
  if (!.s7_is(x@p, Constant)) {
    return(is_nonneg(x@p) && is_nonneg(.args(x)[[idx]]))
  }
  pval <- x@p_used
  if (pval >= 0 && pval <= 1) return(TRUE)
  if (pval > 1) {
    if (is_power2(pval)) {
      return(is_nonneg(.args(x)[[idx]]))
    } else {
      return(TRUE)
    }
  }
  FALSE  # p < 0
}

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

# -- is_constant: p == 0 makes it constant ------------------------
## CVXPY SOURCE: power.py:278-281 -- the exponent as written, not p_used.
method(is_constant, Power) <- function(x) {
  if (.s7_is(x@p, Constant) && as.numeric(value(x@p)) == 0) return(TRUE)
  ## Fall through to default (all args constant)
  (0L %in% .shape(x)) || .all_args(x, is_constant)
}

# -- quadratic/PWL ------------------------------------------------
## CVXPY: power.py lines 317-362
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(.args(x)[[1L]]))
  if (pval == 2) return(is_affine(.args(x)[[1L]]))
  is_constant(.args(x)[[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(.args(x)[[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(.args(x)[[1L]]))
  if (pval == 2) return(is_pwl(.args(x)[[1L]]))
  is_constant(.args(x)[[1L]])
}

# -- domain -------------------------------------------------------
## CVXPY SOURCE: power.py:401-415 -- the exponent as written (`_p_orig`, or
## `self.p.value` when it was given as an Expression), not p_used.
method(atom_domain, Power) <- function(x) {
  pval <- if (!.s7_is(x@p_orig, Expression)) x@p_orig else value(x@p)
  if (is.null(pval)) {
    cli_abort("Cannot compute domain of parametrized power when parameter value is unspecified.")
  }
  pval <- as.numeric(pval)
  if ((pval < 1 && pval != 0) || (pval > 1 && !is_power2(pval))) {
    return(list(.args(x)[[1L]] >= 0))
  }
  list()
}

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

# -- name ---------------------------------------------------------
## CVXPY SOURCE: power.py:437-443 -- `type(self).__name__` and the exponent
## as written (`self.p.value`), not p_used. An unset Parameter prints "?"
## where Python prints None.
.power_p_label <- function(x) {
  pv <- value(x@p)
  if (is.null(pv)) "?" else format(as.numeric(pv), digits = 15L)
}

method(expr_name, Power) <- function(x) {
  sprintf("%s(%s, %s)", short_class_name(x), expr_name(.args(x)[[1L]]),
          .power_p_label(x))
}

method(format_labeled, Power) <- function(x) {
  lbl <- label(x); if (!is.null(lbl)) return(lbl)
  sprintf("%s(%s, %s)", short_class_name(x), format_labeled(.args(x)[[1L]]),
          .power_p_label(x))
}

# -- numeric ------------------------------------------------------
## CVXPY SOURCE: power.py:192-194 -- `np.power(values[0], float(self.p.value))`,
## the exponent as written, not p_used.
method(numeric_value, Power) <- function(x, values, ...) {
  values[[1L]]^as.numeric(value(x@p))
}

# -- graph_implementation: stub -----------------------------------
## CVXPY SOURCE: power.py lines 251-276
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 <- .args(x)[[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.")
}

# -- .grad: per-atom subgradient ----------------------------------
## CVXPY SOURCE: atoms/elementwise/power.py:369-399 (Power._grad).
## d/dx (x^p) = p * x^(p-1)  (elementwise diagonal).
## Edge cases:
##   p == 0          -> all zeros (constant atom).
##   p not power-of-2, x <= 0  -> outside domain, return list(NULL).
method(.grad, Power) <- function(x, values, ...) {
  rows <- as.integer(prod(.arg_shape(x)))
  cols <- as.integer(prod(.shape(x)))

  if (!is.null(x@p_used)) {
    p <- x@p_used
  } else if (!is.null(value(x@p))) {
    p <- as.numeric(value(x@p))
  } else {
    cli_abort(c(
      "Cannot compute {.fn .grad} for {.cls Power} when {.arg p} is unset.",
      "i" = "Set the parameter value before evaluating the gradient."
    ))
  }

  if (p == 0) {
    return(list(Matrix::sparseMatrix(
      i = integer(0), j = integer(0),
      x = numeric(0),
      dims = c(rows, cols),
      repr = "C"
    )))
  }

  v <- values[[1L]]
  if (!is_power2(p) && min(v) <= 0) {
    return(list(NULL))
  }

  ## Cast p to plain numeric only here, mirroring CVXPY's
  ## `float(p)*np.power(values[0], float(p)-1)` (power.py:398). p may be a
  ## bigq from PowerApprox (analogue of Python's Fraction); this cast
  ## avoids gmp::`^.bigq` and matches CVXPY's `float(p)` cast point.
  p_num <- as.numeric(p)
  grad_vals <- p_num * v^(p_num - 1)
  list(.elemwise_grad_to_diag(grad_vals, rows, cols))
}

# ===================================================================
# PowerApprox -- SOC-based rational approximation of Power
# ===================================================================
## CVXPY SOURCE: atoms/elementwise/power.py lines 446-473
## 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 (FALSE) new_object(S7_object())  ## S7 static-check guard
    if (is.null(id)) id <- next_expr_id()
    x <- as_expr(x)
    shape <- .shape(x)

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

    ## CVXPY SOURCE: power.py Power.__init__ -- the exponent must be a Constant
    ## or a Parameter (mirrors Power; const-base / variable-exponent handled by
    ## the power() factory via .pow_const_base).
    if (!(.s7_is(p_const, Constant) || .s7_is(p_const, Parameter))) {
      cli_abort("The exponent {.arg p} must be either a Constant or a Parameter.")
    }

    ## Compute p_used (float value for DCP)
    if (.s7_is(p_const, Constant)) {
      p_used <- as.numeric(value(p_const))
      ## power.py:183 float(self.p.value) rejects a non-scalar exponent.
      if (length(p_used) != 1L)
        cli_abort("The exponent must be a scalar, not of length {length(p_used)}.")
    } 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 <- .fast_new(PowerApprox, 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
#'
#' CVXPY's \code{power(x, p)}: \code{x} raised to the power \code{p}, for
#' any argument (a plain number becomes a constant). \code{cvxr_power()} is
#' the preferred name. \code{power()} is the same function, kept because
#' existing code uses it; it masks \code{\link[stats]{power}}, the unrelated
#' \code{glm()} link constructor, so with CVXR attached write
#' \code{stats::power()} in a \code{glm()} family. \code{x^p} is the same
#' atom for an expression \code{x}.
#'
#' @param x An Expression (the base), OR a positive constant if \code{p} is a
#'   variable (the identity \code{b^x = exp(x * log(b))} is used).
#' @param p Numeric exponent, Parameter, or Expression. If \code{p} is a
#'   non-constant Expression and \code{x} is a positive constant, dispatches to
#'   \code{exp(p * log(x))}.
#' @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, or an exp expression for the const-base case.
#' @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
cvxr_power <- function(x, p, max_denom = 1024L, approx = TRUE) {
  .power_atom(x, p, max_denom, approx)
}

#' @rdname cvxr_power
#' @export
power <- cvxr_power

## CVXPY SOURCE: power.py:31 -- def power(x, p, max_denom=1024, approx=True).
## Internal callers use this directly.
.power_atom <- function(x, p, max_denom = 1024L, approx = TRUE) {
  ## CVXPY v1.9.0 fix: #3180 (power.py:58-62) -- b^x where b is constant and x
  ## is a (non-constant) variable: dispatch to exp(x * log(b)).
  x_expr <- as_expr(x)
  p_expr <- as_expr(p)
  if (is_constant(x_expr) && !is_constant(p_expr)) {
    return(.pow_const_base(x_expr, p_expr))
  }
  if (approx) {
    PowerApprox(x, p, max_denom)
  } else {
    Power(x, p, max_denom)
  }
}