R/014_utilities_bounds.R

In CVXR: Disciplined Convex Optimization

#####
## DO NOT EDIT THIS FILE!! EDIT THE SOURCE INSTEAD: rsrc_tree/utilities/bounds.R
#####

## CVXPY SOURCE: utilities/bounds.py
## Interval-arithmetic bounds propagation (CVXPY 1.9.0, #3080). Each helper
## takes interval bounds (lower, upper) of operand(s) and returns the interval
## bounds of an operation's result. Atoms call these from `bounds_from_args`;
## the recursion bottoms out at Leaf$get_bounds.
##
## Representation: a "Bounds" is `list(lb, ub)` where lb/ub are real R matrices
## matching the expression shape (column-major, R's native order = numpy 'F').
## Scalars are length-1. This mirrors bounds.py's `Bounds = (np.ndarray, np.ndarray)`.
##
## Axis convention: the reduction helpers (sum/max/min/norm) take a NUMPY-style
## `axis` (NULL = all, 0 = down rows -> per-column, 1 = across cols -> per-row),
## exactly as bounds.py does. The calling atom converts CVXR's 1-based axis to
## this convention before calling (mirrors "axis compiled away" -- ADR D_CONV.1).
##
## R note: sparse Matrix inputs are accepted where the CVXR code can preserve them,
## but helpers that need elementwise inspection may densify via `.bnd_dense`.


# -- internal helpers ---------------------------------------------

## Elementwise select preserving the shape of `cond` (mirrors np.where).
## `yes`/`no` are scalars or same-shape as cond; both are evaluated (like numpy),
## so wrap domain-violating computations (log of <=0, etc.) in suppressWarnings.
.bnd_where <- function(cond, yes, no) {
  out <- ifelse(cond, yes, no)         # ifelse keeps cond's dim attribute
  dim(out) <- dim(cond)
  out
}

## Densify a possibly-sparse Matrix bound to a base-R matrix (dense path only).
.bnd_dense <- function(arr) {
  if (inherits(arr, "Matrix")) as.matrix(arr) else arr
}

.bounds_all_isinf <- function(arr) {
  arr <- .bnd_dense(arr)
  all(is.infinite(arr))
}

.bounds_any_isnan <- function(arr) {
  arr <- .bnd_dense(arr)
  any(is.nan(arr))
}

.bounds_all_zero_or_inf <- function(arr) {
  arr <- .bnd_dense(arr)
  all(arr == 0 | is.infinite(arr))
}

.bounds_ensure_dense <- function(arr, shape) {
  arr <- .bnd_dense(arr)
  if (length(arr) == 1L) {
    array(as.numeric(arr), dim = shape)
  } else {
    dim(arr) <- shape
    arr
  }
}

## CVXPY SOURCE: bounds.py:809-856
## Return auxiliary-variable bounds only when the chosen solver consumes native
## variable bounds and the expression's interval bounds carry useful information.
get_expr_bounds_if_supported <- function(expr, solver_context) {
  if (is.null(solver_context) || !isTRUE(solver_context@solver_supports_bounds)) {
    return(NULL)
  }
  b <- get_bounds(expr)
  lb <- b[[1L]]
  ub <- b[[2L]]

  if (.bounds_all_isinf(lb) && .bounds_all_isinf(ub)) return(NULL)
  if (.bounds_any_isnan(lb) || .bounds_any_isnan(ub)) return(NULL)

  if (.bounds_all_zero_or_inf(lb) && .bounds_all_isinf(ub) && is_nonneg(expr)) {
    return(NULL)
  }
  if (.bounds_all_isinf(lb) && .bounds_all_zero_or_inf(ub) && is_nonpos(expr)) {
    return(NULL)
  }

  list(.bounds_ensure_dense(lb, .shape(expr)), .bounds_ensure_dense(ub, .shape(expr)))
}

# -- Bounds constructors ------------------------------------------

## CVXPY SOURCE: bounds.py:110-125
unbounded <- function(shape) {
  list(array(-Inf, dim = shape), array(Inf, dim = shape))
}

## CVXPY SOURCE: bounds.py:128-151 (broadcast views -> plain filled arrays in R)
uniform_bounds <- function(shape, lb, ub) {
  list(array(lb, dim = shape), array(ub, dim = shape))
}

## CVXPY SOURCE: bounds.py:154-169
scalar_bounds <- function(lb, ub) {
  list(lb, ub)
}

# -- affine / elementwise binary ----------------------------------

## CVXPY SOURCE: bounds.py:172-188 (x + y)
add_bounds <- function(lb1, ub1, lb2, ub2) {
  list(lb1 + lb2, ub1 + ub2)
}

## CVXPY SOURCE: bounds.py:233-246 (-x)
neg_bounds <- function(lb, ub) {
  list(-ub, -lb)
}

## CVXPY SOURCE: bounds.py:249-298 (x * y, elementwise, interval arithmetic).
## NaN from 0*Inf -> 0 (interval arithmetic: 0 * anything = 0).
mul_bounds <- function(lb1, ub1, lb2, ub2) {
  .mul <- function(a, b) {
    p <- suppressWarnings(a * b)
    p[is.nan(p)] <- 0
    p
  }
  p1 <- .mul(lb1, lb2)
  p2 <- .mul(lb1, ub2)
  p3 <- .mul(ub1, lb2)
  p4 <- .mul(ub1, ub2)
  list(pmin(p1, p2, p3, p4), pmax(p1, p2, p3, p4))
}

## CVXPY SOURCE: bounds.py:301-332 (x / y). Divisor interval containing 0 -> unbounded.
div_bounds <- function(lb1, ub1, lb2, ub2) {
  lb2 <- .bnd_dense(lb2)
  ub2 <- .bnd_dense(ub2)
  contains_zero <- (lb2 <= 0) & (ub2 >= 0)
  inv_lb <- .bnd_where(contains_zero, -Inf, 1.0 / ub2)
  inv_ub <- .bnd_where(contains_zero,  Inf, 1.0 / lb2)
  mul_bounds(lb1, ub1, inv_lb, inv_ub)
}

## CVXPY SOURCE: bounds.py:335-353 (c * x, scalar c)
scale_bounds <- function(lb, ub, c) {
  if (c >= 0) list(c * lb, c * ub) else list(c * ub, c * lb)
}

## CVXPY SOURCE: bounds.py:356-381 (|x|)
abs_bounds <- function(lb, ub) {
  spans_zero <- (lb <= 0) & (ub >= 0)
  entirely_positive <- lb >= 0
  entirely_negative <- ub <= 0
  new_lb <- .bnd_where(spans_zero, 0.0,
              .bnd_where(entirely_positive, lb, -ub))
  new_ub <- .bnd_where(entirely_positive, ub,
              .bnd_where(entirely_negative, -lb, pmax(-lb, ub)))
  list(new_lb, new_ub)
}

# -- elementwise max/min over several operands --------------------

## CVXPY SOURCE: bounds.py:384-403 (max(x1, x2, ...))
maximum_bounds <- function(bounds_list) {
  lb <- bounds_list[[1L]][[1L]]
  ub <- bounds_list[[1L]][[2L]]
  for (b in bounds_list[-1L]) {
    lb <- pmax(lb, b[[1L]])
    ub <- pmax(ub, b[[2L]])
  }
  list(lb, ub)
}

## CVXPY SOURCE: bounds.py:406-425 (min(x1, x2, ...))
minimum_bounds <- function(bounds_list) {
  lb <- bounds_list[[1L]][[1L]]
  ub <- bounds_list[[1L]][[2L]]
  for (b in bounds_list[-1L]) {
    lb <- pmin(lb, b[[1L]])
    ub <- pmin(ub, b[[2L]])
  }
  list(lb, ub)
}

# -- axis reductions (numpy-style axis: NULL/0/1, optional keepdims) ----

## Reduce a matrix `m` with reducer `fun` (e.g. sum/max/min) along a numpy axis.
.bnd_reduce <- function(m, fun, axis, keepdims) {
  if (!is.matrix(m)) m <- as.matrix(m)
  if (is.null(axis)) {
    val <- fun(m)
    return(if (keepdims) matrix(val, 1L, 1L) else val)
  }
  if (axis == 0L) {                      # down rows -> one value per column
    val <- apply(m, 2L, fun)
    if (keepdims) matrix(val, nrow = 1L) else val
  } else {                               # axis == 1: across cols -> per row
    val <- apply(m, 1L, fun)
    if (keepdims) matrix(val, ncol = 1L) else val
  }
}

## CVXPY SOURCE: bounds.py:203-230 (sum reduction)
sum_bounds <- function(lb, ub, axis = NULL, keepdims = FALSE) {
  list(.bnd_reduce(lb, sum, axis, keepdims),
       .bnd_reduce(ub, sum, axis, keepdims))
}

## CVXPY SOURCE: bounds.py:428-449 (max reduction)
max_reduction_bounds <- function(lb, ub, axis = NULL, keepdims = FALSE) {
  list(.bnd_reduce(lb, max, axis, keepdims),
       .bnd_reduce(ub, max, axis, keepdims))
}

## CVXPY SOURCE: bounds.py:452-473 (min reduction)
min_reduction_bounds <- function(lb, ub, axis = NULL, keepdims = FALSE) {
  list(.bnd_reduce(lb, min, axis, keepdims),
       .bnd_reduce(ub, min, axis, keepdims))
}

# -- elementwise nonlinear ----------------------------------------

## CVXPY SOURCE: bounds.py:476-549 (x^p)
power_bounds <- function(lb, ub, p) {
  if (p == 0) {
    return(list(array(1, dim = dim(lb) %||% length(lb)),
                array(1, dim = dim(ub) %||% length(ub))))
  }
  if (p > 0) {
    if (p == as.integer(p) && as.integer(p) %% 2L == 0L) {
      ## even integer power
      spans_zero <- (lb <= 0) & (ub >= 0)
      entirely_positive <- lb >= 0
      lb_power <- abs(lb)^p
      ub_power <- abs(ub)^p
      new_lb <- .bnd_where(spans_zero, 0.0,
                  .bnd_where(entirely_positive, lb_power, ub_power))
      new_ub <- .bnd_where(spans_zero, pmax(lb_power, ub_power),
                  .bnd_where(entirely_positive, ub_power, lb_power))
      list(new_lb, new_ub)
    } else if (p == as.integer(p)) {
      ## odd integer power: monotonic
      list(lb^p, ub^p)
    } else {
      ## non-integer positive power: requires x >= 0 for a real result
      valid <- lb >= 0
      new_lb <- .bnd_where(valid, suppressWarnings(lb^p), -Inf)
      new_ub <- .bnd_where(valid, suppressWarnings(ub^p),  Inf)
      list(new_lb, new_ub)
    }
  } else {
    ## negative power
    spans_zero <- (lb <= 0) & (ub >= 0)
    entirely_positive <- lb > 0
    if (p == as.integer(p) && as.integer(-p) %% 2L == 1L) {
      ## negative odd power: monotonically decreasing
      new_lb <- .bnd_where(spans_zero, -Inf, pmin(lb^p, ub^p))
      new_ub <- .bnd_where(spans_zero,  Inf, pmax(lb^p, ub^p))
    } else {
      ## negative even / non-integer: only for positive x
      new_lb <- .bnd_where(entirely_positive, suppressWarnings(ub^p), -Inf)
      new_ub <- .bnd_where(entirely_positive, suppressWarnings(lb^p),  Inf)
    }
    list(new_lb, new_ub)
  }
}

## CVXPY SOURCE: bounds.py:552-567 (exp x, monotone increasing)
exp_bounds <- function(lb, ub) {
  list(exp(lb), exp(ub))
}

## CVXPY SOURCE: bounds.py:570-589 (log x, defined on x > 0)
log_bounds <- function(lb, ub) {
  new_lb <- .bnd_where(lb > 0, suppressWarnings(log(lb)), -Inf)
  new_ub <- .bnd_where(ub > 0, suppressWarnings(log(ub)),  Inf)
  list(new_lb, new_ub)
}

## CVXPY SOURCE: bounds.py:592-610 (sqrt x, defined on x >= 0)
sqrt_bounds <- function(lb, ub) {
  new_lb <- .bnd_where(lb >= 0, suppressWarnings(sqrt(lb)), -Inf)
  new_ub <- .bnd_where(ub >= 0, suppressWarnings(sqrt(ub)),  Inf)
  list(new_lb, new_ub)
}

# -- norms (compose abs + reduction) ------------------------------

## CVXPY SOURCE: bounds.py:613-633 (1-norm = sum(|x|))
norm1_bounds <- function(lb, ub, axis = NULL, keepdims = FALSE) {
  ab <- abs_bounds(lb, ub)
  sum_bounds(ab[[1L]], ab[[2L]], axis = axis, keepdims = keepdims)
}

## CVXPY SOURCE: bounds.py:636-656 (inf-norm = max(|x|))
norm_inf_bounds <- function(lb, ub, axis = NULL, keepdims = FALSE) {
  ab <- abs_bounds(lb, ub)
  max_reduction_bounds(ab[[1L]], ab[[2L]], axis = axis, keepdims = keepdims)
}

# -- structural ---------------------------------------------------

## CVXPY SOURCE: bounds.py:659-685 (broadcast to a target shape)
broadcast_bounds <- function(lb, ub, target_shape) {
  list(array(.bnd_dense(lb), dim = target_shape),
       array(.bnd_dense(ub), dim = target_shape))
}

## CVXPY SOURCE: bounds.py:688-707 (reshape; order 'F' = column-major = R native)
reshape_bounds <- function(lb, ub, new_shape, order = "F") {
  .rs <- function(a) {
    a <- .bnd_dense(a)
    if (order == "F") {
      ## column-major (R native): preserve the column-major element order
      array(as.numeric(a), dim = new_shape)
    } else {
      ## 'C' (row-major, 2D): read `a` row-major (as.numeric(t(a))), then fill
      ## the target row-major (byrow). Verified against numpy reshape(order='C').
      if (!is.matrix(a)) a <- as.matrix(a)
      matrix(as.numeric(t(a)), nrow = new_shape[1L], ncol = new_shape[2L],
             byrow = TRUE)
    }
  }
  list(.rs(lb), .rs(ub))
}

## CVXPY SOURCE: bounds.py:710-726 (transpose)
transpose_bounds <- function(lb, ub) {
  list(t(.bnd_dense(lb)), t(.bnd_dense(ub)))
}

## CVXPY SOURCE: bounds.py:729-744 (index/slice; `key` is an R index expression
## applied via `[`. Caller passes drop = FALSE semantics as needed.)
index_bounds <- function(lb, ub, key_fn) {
  list(key_fn(.bnd_dense(lb)), key_fn(.bnd_dense(ub)))
}

## CVXPY SOURCE: bounds.py:747-806 (matrix product x %*% y).
## Exact split formula when one operand is a point (lb == ub); otherwise unbounded.
matmul_bounds <- function(lb1, ub1, lb2, ub2) {
  lb1 <- .bnd_dense(lb1); ub1 <- .bnd_dense(ub1)
  lb2 <- .bnd_dense(lb2); ub2 <- .bnd_dense(ub2)
  lhs_point <- isTRUE(all.equal(lb1, ub1)) && identical(dim(lb1), dim(ub1))
  rhs_point <- isTRUE(all.equal(lb2, ub2)) && identical(dim(lb2), dim(ub2))

  if (lhs_point) {
    a_pos <- pmax(lb1, 0); a_neg <- pmin(lb1, 0)
    return(list(a_pos %*% lb2 + a_neg %*% ub2,
                a_pos %*% ub2 + a_neg %*% lb2))
  }
  if (rhs_point) {
    b_pos <- pmax(lb2, 0); b_neg <- pmin(lb2, 0)
    return(list(lb1 %*% b_pos + ub1 %*% b_neg,
                ub1 %*% b_pos + lb1 %*% b_neg))
  }
  ## both intervals: no efficient exact formula
  unbounded(c(nrow(as.matrix(lb1)), ncol(as.matrix(lb2))))
}

# -- sign refinement ----------------------------------------------

## CVXPY SOURCE: bounds.py:879-901 (tighten with known sign)
refine_bounds_from_sign <- function(lb, ub, is_nonneg, is_nonpos) {
  if (isTRUE(is_nonneg)) lb <- pmax(lb, 0)
  if (isTRUE(is_nonpos)) ub <- pmin(ub, 0)
  list(lb, ub)
}

## CVXPY SOURCE: bounds.py:859-876 (sparsity-pattern equality; dense-coords here)
coords_equal <- function(coords1, coords2) {
  if (length(coords1) != length(coords2)) return(FALSE)
  all(vapply(seq_along(coords1),
             function(i) isTRUE(all.equal(coords1[[i]], coords2[[i]])),
             logical(1L)))
}