R/128_constraints_power.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: constraints/power.py
## PowCone3D -- 3D Power Cone constraint
## PowConeND -- N-dimensional Power Cone constraint

# -- Tolerance for PowConeND alpha sum check ----------------------
## CVXPY SOURCE: power.py line 182
.POWCONE_TOL <- 1e-6

# ==================================================================
# PowCone3D
# ==================================================================

#' Create a 3D Power Cone Constraint
#'
#' Constrains \eqn{(x, y, z)} to lie in the 3D power cone:
#' \deqn{x^\alpha \cdot y^{1-\alpha} \ge |z|, \quad x \ge 0, \; y \ge 0}
#'
#' @param x_expr A CVXR expression.
#' @param y_expr A CVXR expression.
#' @param z_expr A CVXR expression.
#' @param alpha A CVXR expression or numeric value in \eqn{(0, 1)}.
#' @param constr_id Optional integer constraint ID.
#' @returns A \code{PowCone3D} constraint object.
#' @seealso \code{\link{PowConeND}}
#' @export
PowCone3D <- new_class("PowCone3D", parent = Cone, package = "CVXR",
  properties = list(
    .x    = new_property(class = Expression),
    .y    = new_property(class = Expression),
    .z    = new_property(class = Expression),
    alpha = new_property(class = Expression)
  ),
  constructor = function(x_expr, y_expr, z_expr, alpha, constr_id = NULL) {
    x_expr <- as_expr(x_expr)
    y_expr <- as_expr(y_expr)
    z_expr <- as_expr(z_expr)

    ## CVXPY SOURCE: power.py lines 48-50
    for (val in list(x_expr, y_expr, z_expr)) {
      if (!(is_affine(val) && is_real(val))) {
        cli_abort("All arguments must be affine and real.")
      }
    }

    ## CVXPY SOURCE: power.py lines 51-61 -- alpha promotion
    alpha <- as_expr(alpha)
    alpha_promoted_to_vec <- FALSE
    if (expr_is_scalar(alpha)) {
      if (!expr_is_scalar(x_expr)) {
        alpha <- cvxr_promote(alpha, x_expr@shape)
      } else {
        alpha_promoted_to_vec <- TRUE
        alpha <- cvxr_promote(alpha, c(1L, 1L))
      }
    }

    ## CVXPY SOURCE: power.py lines 63-65
    alpha_val <- value(alpha)
    if (any(alpha_val <= 0) || any(alpha_val >= 1)) {
      cli_abort("Argument {.arg alpha} must have entries in the open interval (0, 1).")
    }

    ## CVXPY SOURCE: power.py lines 66-73 -- shape compatibility
    if (alpha_promoted_to_vec) {
      arg_shapes <- list(x_expr@shape, y_expr@shape, z_expr@shape)
    } else {
      arg_shapes <- list(x_expr@shape, y_expr@shape, z_expr@shape, alpha@shape)
    }
    ref <- arg_shapes[[1L]]
    for (s in arg_shapes[-1L]) {
      if (!identical(ref, s)) {
        cli_abort(
          "All arguments must have the same shapes. Provided shapes: {.val {paste(lapply(arg_shapes, paste, collapse='x'), collapse=', ')}}."
        )
      }
    }

    if (is.null(constr_id)) constr_id <- next_expr_id()

    ## CVXPY SOURCE: power.py lines 127-131 -- shape override
    ## shape = (3,) + self.x.shape
    cone_shape <- c(3L, as.integer(prod(x_expr@shape)))

    args <- list(x_expr, y_expr, z_expr)
    dvars <- lapply(args, function(a) Variable(a@shape))

    new_object(S7_object(),
      id    = as.integer(constr_id),
      .cache = new.env(parent = emptyenv()),
      args  = args,
      dual_variables = dvars,
      shape = cone_shape,
      .label = "",
      .x    = x_expr,
      .y    = y_expr,
      .z    = z_expr,
      alpha = alpha
    )
  }
)

# -- expr_name ----------------------------------------------------
## CVXPY SOURCE: power.py lines 76-77

method(expr_name, PowCone3D) <- function(x) {
  sprintf("PowCone3D(%s, %s, %s; %s)",
    expr_name(x@.x), expr_name(x@.y), expr_name(x@.z), expr_name(x@alpha))
}

# -- is_dcp -------------------------------------------------------
## CVXPY SOURCE: power.py lines 113-119

method(is_dcp, PowCone3D) <- function(x) {
  .all_args(x, is_affine)
}

# -- is_dgp -------------------------------------------------------
## CVXPY SOURCE: power.py lines 121-122

method(is_dgp, PowCone3D) <- function(x) FALSE

# -- get_data -----------------------------------------------------
## CVXPY SOURCE: power.py lines 94-95

method(get_data, PowCone3D) <- function(x) {
  list(x@alpha, x@id)
}

# -- num_cones ----------------------------------------------------
## CVXPY SOURCE: power.py lines 107-108

method(num_cones, PowCone3D) <- function(x) {
  expr_size(x@.x)
}

# -- cone_sizes ---------------------------------------------------
## CVXPY SOURCE: power.py lines 110-111

method(cone_sizes, PowCone3D) <- function(x) {
  rep(3L, num_cones(x))
}

# -- constr_size --------------------------------------------------
## CVXPY SOURCE: power.py lines 103-105

method(constr_size, PowCone3D) <- function(x) {
  3L * num_cones(x)
}

# -- residual -----------------------------------------------------
## CVXPY SOURCE: power.py lines 79-92
## TODO: projection should be implemented directly

method(residual, PowCone3D) <- function(x) {
  xv <- value(x@.x)
  yv <- value(x@.y)
  zv <- value(x@.z)
  if (is.null(xv) || is.null(yv) || is.null(zv)) return(NULL)

  alpha_val <- as.numeric(value(x@alpha))
  xv <- as.numeric(xv)
  yv <- as.numeric(yv)
  zv <- as.numeric(zv)

  ## Feasibility: x^alpha * y^(1-alpha) >= |z|, x >= 0, y >= 0
  n <- length(xv)
  resid <- numeric(n)
  for (i in seq_len(n)) {
    if (xv[i] >= 0 && yv[i] >= 0) {
      lhs <- xv[i]^alpha_val[i] * yv[i]^(1 - alpha_val[i])
      resid[i] <- max(0, abs(zv[i]) - lhs)
    } else {
      resid[i] <- abs(min(0, xv[i])) + abs(min(0, yv[i])) + abs(zv[i])
    }
  }
  if (n == 1L) resid[1L] else resid
}

# -- save_dual_value ----------------------------------------------
## CVXPY SOURCE: power.py lines 133-142
## NOTE: Different reshape order than ExpCone -- (3, -1) not (-1, 3)
## Uses C-order reshape (CR-1)

method(save_dual_value, PowCone3D) <- function(x, val) {
  nc <- num_cones(x)
  ## CVXPY: np.reshape(value, (3, -1)) in C-order
  val_mat <- .reshape_c_order(val, 3L, nc)
  dv0 <- matrix(val_mat[1L, ], nrow = x@.x@shape[1L], ncol = x@.x@shape[2L])
  dv1 <- matrix(val_mat[2L, ], nrow = x@.y@shape[1L], ncol = x@.y@shape[2L])
  dv2 <- matrix(val_mat[3L, ], nrow = x@.z@shape[1L], ncol = x@.z@shape[2L])
  value(x@dual_variables[[1L]]) <- dv0
  value(x@dual_variables[[2L]]) <- dv1
  value(x@dual_variables[[3L]]) <- dv2
  invisible(x)
}

# -- dual_cone ----------------------------------------------------
## CVXPY SOURCE: power.py lines 144-157
## dual = PowCone3D(x/alpha, y/(1-alpha), z, alpha)
## NOTE: CVXPY has a missing `return` on the `args is None` branch

method(dual_cone, PowCone3D) <- function(x, ...) {
  args <- list(...)
  if (length(args) == 0L) {
    PowCone3D(x@dual_variables[[1L]] / x@alpha,
              x@dual_variables[[2L]] / (1 - x@alpha),
              x@dual_variables[[3L]],
              x@alpha)
  } else {
    PowCone3D(args[[1L]] / x@alpha,
              args[[2L]] / (1 - x@alpha),
              args[[3L]],
              x@alpha)
  }
}

# -- is_complex / is_imag ----------------------------------------
## CVXPY SOURCE: power.py lines 97-101

method(is_complex, PowCone3D) <- function(x) FALSE
method(is_imag, PowCone3D) <- function(x) FALSE


# ==================================================================
# PowConeND
# ==================================================================

#' Create an N-Dimensional Power Cone Constraint
#'
#' Constrains \eqn{(W, z)} to lie in the N-dimensional power cone:
#' \deqn{\prod W_i^{\alpha_i} \ge |z|, \quad W \ge 0}
#' where \eqn{\alpha_i > 0} and \eqn{\sum \alpha_i = 1}.
#'
#' @param W A CVXR expression (vector or matrix).
#' @param z A CVXR expression (scalar or vector).
#' @param alpha A CVXR expression with positive entries summing to 1
#'   along the specified axis.
#' @param axis Integer, 2 (default, column-wise) or 1 (row-wise).
#' @param constr_id Optional integer constraint ID.
#' @returns A \code{PowConeND} constraint object.
#'
#' @section Known limitations:
#' The R \pkg{clarabel} solver does not currently support the PowConeND cone
#' specification. Problems involving PowConeND (e.g., exact geometric mean
#' with more than 2 arguments) should use SCS or MOSEK as the solver, or
#' use approximation-based atoms (e.g., \code{geo_mean(x, approx = TRUE)}).
#'
#' @seealso \code{\link{PowCone3D}}
#' @export
PowConeND <- new_class("PowConeND", parent = Cone, package = "CVXR",
  properties = list(
    .W    = new_property(class = Expression),
    .z    = new_property(class = Expression),
    alpha = new_property(class = Expression),
    axis  = new_property(class = class_integer, default = 2L)
  ),
  constructor = function(W, z, alpha, axis = 2L, constr_id = NULL) {
    W <- as_expr(W)
    z <- as_expr(z)
    axis <- as.integer(axis)

    ## CVXPY SOURCE: power.py lines 186-189
    if (!(is_real(W) && is_affine(W))) {
      cli_abort("Invalid first argument; {.arg W} must be affine and real.")
    }

    ## CVXPY SOURCE: power.py lines 190-194
    ## z must be <= 1D, affine, real
    ## In our R shapes (always 2D), "1D" means min(shape) == 1
    if (!expr_is_vector(z) || !(is_real(z) && is_affine(z))) {
      cli_abort("Invalid second argument. {.arg z} must be affine, real, and at most 1D.")
    }

    ## CVXPY SOURCE: power.py lines 196-201 -- dimension compatibility
    W_shape <- W@shape
    z_size <- expr_size(z)
    if ((expr_is_vector(W) && z_size > 1L) ||
        (expr_is_matrix(W) && z_size != W_shape[axis]) ||
        (expr_is_vector(W) && axis == 1L)) {
      cli_abort(
        "Argument dimensions {.val {W_shape}} and {.val {z@shape}}, with axis={axis}, are incompatible."
      )
    }

    ## CVXPY SOURCE: power.py lines 202-204
    if (expr_is_matrix(W) && W_shape[3L - axis] <= 1L) {
      cli_abort("PowConeND requires left-hand-side to have at least two terms.")
    }

    ## CVXPY SOURCE: power.py lines 205-208
    alpha <- as_expr(alpha)
    if (!identical(alpha@shape, W_shape)) {
      cli_abort(
        "Argument dimensions {.val {W_shape}} and {.val {alpha@shape}} are not equal."
      )
    }

    ## CVXPY SOURCE: power.py lines 209-212
    alpha_val <- value(alpha)
    if (any(as.numeric(alpha_val) <= 0)) {
      cli_abort("Argument {.arg alpha} must be entry-wise positive.")
    }
    alpha_sums <- apply(as.matrix(alpha_val), axis, sum)
    if (any(abs(1 - alpha_sums) > .POWCONE_TOL)) {
      cli_abort("Argument {.arg alpha} must sum to 1 along axis {axis}.")
    }

    if (is.null(constr_id)) constr_id <- next_expr_id()

    ## CVXPY SOURCE: power.py lines 233-247 -- shape override
    if (expr_is_vector(W)) {
      m <- W_shape[1L]
      ncol_val <- 1L
    } else {
      m <- W_shape[3L - axis]    ## cone dimension
      ncol_val <- W_shape[axis]  ## number of cones
    }
    cone_shape <- c(m + 1L, ncol_val)

    ## CVXPY SOURCE: power.py lines 217-219 -- flatten 0-d z
    if (expr_is_scalar(z)) {
      z <- reshape_expr(z, c(1L, 1L))
    }

    args <- list(W, z)
    dvars <- lapply(args, function(a) Variable(a@shape))

    new_object(S7_object(),
      id    = as.integer(constr_id),
      .cache = new.env(parent = emptyenv()),
      args  = args,
      dual_variables = dvars,
      shape = cone_shape,
      .label = "",
      .W    = W,
      .z    = z,
      alpha = alpha,
      axis  = axis
    )
  }
)

# -- expr_name ----------------------------------------------------
## CVXPY SOURCE: power.py lines 221-222

method(expr_name, PowConeND) <- function(x) {
  sprintf("PowConeND(%s, %s; %s)",
    expr_name(x@.W), expr_name(x@.z), expr_name(x@alpha))
}

# -- is_dcp -------------------------------------------------------
## CVXPY SOURCE: power.py lines 276-284
## Returns TRUE always (non-DPP mode)

method(is_dcp, PowConeND) <- function(x) TRUE

# -- is_dgp -------------------------------------------------------
## CVXPY SOURCE: power.py lines 286-287

method(is_dgp, PowConeND) <- function(x) FALSE

# -- get_data -----------------------------------------------------
## CVXPY SOURCE: power.py lines 230-231

method(get_data, PowConeND) <- function(x) {
  list(x@alpha, x@axis, x@id)
}

# -- num_cones ----------------------------------------------------
## CVXPY SOURCE: power.py lines 264-265

method(num_cones, PowConeND) <- function(x) {
  expr_size(x@.z)
}

# -- cone_sizes ---------------------------------------------------
## CVXPY SOURCE: power.py lines 272-274

method(cone_sizes, PowConeND) <- function(x) {
  cone_size <- 1L + x@args[[1L]]@shape[3L - x@axis]
  rep(cone_size, num_cones(x))
}

# -- constr_size --------------------------------------------------
## CVXPY SOURCE: power.py lines 267-270

method(constr_size, PowConeND) <- function(x) {
  cone_size <- 1L + x@args[[1L]]@shape[3L - x@axis]
  cone_size * num_cones(x)
}

# -- residual -----------------------------------------------------
## CVXPY SOURCE: power.py lines 250-262
## TODO: projection should be implemented directly

method(residual, PowConeND) <- function(x) {
  W_val <- value(x@.W)
  z_val <- value(x@.z)
  if (is.null(W_val) || is.null(z_val)) return(NULL)

  alpha_val <- as.matrix(value(x@alpha))
  W_val <- as.matrix(W_val)
  z_val <- as.numeric(z_val)

  ## Transpose if axis == 1 to normalize to column-wise
  if (x@axis == 1L) {
    W_val <- t(W_val)
    alpha_val <- t(alpha_val)
  }

  n_cones <- length(z_val)
  resid <- numeric(n_cones)
  for (j in seq_len(n_cones)) {
    w_col <- W_val[, j]
    a_col <- alpha_val[, j]
    if (all(w_col >= 0)) {
      lhs <- prod(w_col^a_col)
      resid[j] <- max(0, abs(z_val[j]) - lhs)
    } else {
      resid[j] <- sum(abs(pmin(0, w_col))) + abs(z_val[j])
    }
  }
  if (n_cones == 1L) resid[1L] else resid
}

# -- save_dual_value ----------------------------------------------
## CVXPY SOURCE: power.py lines 292-305
## Value has shape (m+1, k); first m rows are W duals, last row is z duals.

method(save_dual_value, PowConeND) <- function(x, val) {
  if (!is.matrix(val)) {
    val <- matrix(val, nrow = x@shape[1L], ncol = x@shape[2L])
  }
  dW <- val[seq_len(nrow(val) - 1L), , drop = FALSE]
  dz <- val[nrow(val), ]

  if (x@axis == 1L) dW <- t(dW)

  ## Drop trailing dimension if single column
  if (ncol(dW) == 1L) dW <- as.numeric(dW)

  value(x@dual_variables[[1L]]) <- dW
  value(x@dual_variables[[2L]]) <- dz
  invisible(x)
}

# -- dual_cone ----------------------------------------------------
## CVXPY SOURCE: power.py lines 307-320
## dual = PowConeND(W/alpha, z, alpha, axis)

method(dual_cone, PowConeND) <- function(x, ...) {
  args <- list(...)
  if (length(args) == 0L) {
    PowConeND(x@dual_variables[[1L]] / x@alpha,
              x@dual_variables[[2L]],
              x@alpha, axis = x@axis)
  } else {
    PowConeND(args[[1L]] / x@alpha,
              args[[2L]],
              x@alpha, axis = x@axis)
  }
}

# -- is_complex / is_imag ----------------------------------------
## CVXPY SOURCE: power.py lines 224-228

method(is_complex, PowConeND) <- function(x) FALSE
method(is_imag, PowConeND) <- function(x) FALSE