R/065_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))
    } 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:208-216 — 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 defaults to 1.0 when unset
## (matches CVXPY's `float(self.p.value) if ... else 1.0`).
method(bounds_from_args, Power) <- function(x) {
  b <- get_bounds(.args(x)[[1L]])
  ## as.numeric() mirrors CVXPY's `float(self.p.value)` (power.py:209), whose
  ## power_bounds() is typed `p: float`.  It is load-bearing, not cosmetic:
  ## PowerApprox stores `p_used` as an exact <bigq> rational, and `^.bigq`
  ## REFUSES a non-integer exponent ("<bigq> ^ <non-int> is not rational"), so
  ## without the cast get_bounds() errors on every fractional power.  Reachable
  ## in released 1.9.1 via the documented `power(x, 0.5)` default; it stayed
  ## hidden for `sqrt()` only while sqrt wrongly bypassed PowerApprox.
  p <- if (is.null(x@p_used)) 1.0 else as.numeric(x@p_used)
  power_bounds(b[[1L]], b[[2L]], p)
}

# -- 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(.args(x)[[1L]]),
         is_nonpos = is_nonpos(.args(x)[[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
}

## 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: 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(.args(x)[[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(.args(x)[[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% .shape(x)) || .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(.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: 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(.args(x)[[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(.args(x)[[1L]]), pval)
}

## CVXPY SOURCE: atoms/elementwise/power.py Power.format_labeled.
method(format_labeled, Power) <- function(x) {
  lbl <- label(x); if (!is.null(lbl)) return(lbl)
  pval <- if (!is.null(x@p_used)) x@p_used else "?"
  sprintf("Power(%s, %s)", format_labeled(.args(x)[[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 <- .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:344-374 (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:373). 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 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 (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))
    } 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
#'
#' @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
power <- 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)
  }
}