R/015_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, where the
## upstream helper densifies too (`_ensure_dense`, bounds.py:24-50; the dense
## arithmetic of div/broadcast). The predicates below and `matmul_bounds` stay
## sparse-aware, as bounds.py:79-107 and :747-806 do (sweep item 4, pass B).
.bnd_dense <- function(arr) {
  if (inherits(arr, "Matrix")) as.matrix(arr) else arr
}

## is.infinite / is.nan / == dispatch on a Matrix object: a structural zero is
## neither Inf nor NaN and is 0, which is what bounds.py:79-107 reads off `.data`.
.bounds_all_isinf <- function(arr) all(is.infinite(arr))

.bounds_any_isnan <- function(arr) any(is.nan(arr))

.bounds_all_zero_or_inf <- function(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
  }
}

## 2-D shape of a bound. A dim-less vector is a column (CVXR's shape for a
## 1-D leaf is (n, 1)) unless `ref`, the partner's dim, says it can only be a row.
.bnd_shape <- function(arr, ref = NULL) {
  d <- dim(arr)
  if (!is.null(d)) return(as.integer(d))
  n <- length(arr)
  if (n == 1L) return(c(1L, 1L))
  if (!is.null(ref) && ref[2L] == n && ref[1L] != n) return(c(1L, n))
  c(n, 1L)
}

## numpy-style broadcast of ONE bound to `shape` (np.broadcast_to for the 2-D
## case, bounds.py:677-683): a (1, n) row repeats down the rows, an (m, 1) column
## repeats across the columns, a scalar fills. R's own recycling (array(), pmin,
## pmax) is column-major and PERMUTES a row operand, so it must not be used here.
## A sparse bound that already has the target shape is kept as is (bounds.py:679-680).
.bnd_broadcast_to <- function(arr, shape) {
  shape <- as.integer(shape)
  if (inherits(arr, "Matrix")) {
    if (identical(as.integer(dim(arr)), shape)) return(arr)
    arr <- as.matrix(arr)
  }
  src <- .bnd_shape(arr, shape)
  if (identical(src, shape)) {
    if (is.null(dim(arr))) dim(arr) <- shape
    return(arr)
  }
  if (!all(src == shape | src == 1L)) {
    cli_abort("Cannot broadcast bounds of shape ({src[1L]}, {src[2L]}) to ({shape[1L]}, {shape[2L]}).")
  }
  out <- matrix(as.numeric(arr), nrow = src[1L], ncol = src[2L])
  if (src[1L] == 1L && shape[1L] != 1L) out <- out[rep(1L, shape[1L]), , drop = FALSE]
  if (src[2L] == 1L && shape[2L] != 1L) out <- out[, rep(1L, shape[2L]), drop = FALSE]
  out
}

## Bring two bounds to their common numpy broadcast shape before an elementwise
## fold (what np.maximum / np.minimum do implicitly, bounds.py:398-402, 420-424).
## Two dim-less vectors of equal length are left alone.
.bnd_broadcast_pair <- function(a, b) {
  if (is.null(dim(a)) && is.null(dim(b)) && length(a) == length(b)) return(list(a, b))
  shape <- .broadcast_2d(.bnd_shape(a, dim(b)), .bnd_shape(b, dim(a)))
  list(.bnd_broadcast_to(a, shape), .bnd_broadcast_to(b, shape))
}

## 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]) {
    ## np.maximum broadcasts; pmax recycles column-major (permutes a row operand)
    p <- .bnd_broadcast_pair(lb, b[[1L]])
    lb <- pmax(p[[1L]], p[[2L]])
    p <- .bnd_broadcast_pair(ub, b[[2L]])
    ub <- pmax(p[[1L]], p[[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]) {
    ## np.minimum broadcasts; pmin recycles column-major (permutes a row operand)
    p <- .bnd_broadcast_pair(lb, b[[1L]])
    lb <- pmin(p[[1L]], p[[2L]])
    p <- .bnd_broadcast_pair(ub, b[[2L]])
    ub <- pmin(p[[1L]], p[[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) {
  ## np.broadcast_to, not array(): array() recycles column-major and permutes a row
  list(.bnd_broadcast_to(lb, target_shape), .bnd_broadcast_to(ub, 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:774-778 -- `_is_point`: EXACT equality (np.array_equal;
## two sparse operands `(lb - ub).nnz == 0`). The split formula in
## matmul_bounds is only valid for lb == ub, so this must not be a tolerance
## test: `isTRUE(all.equal(lb, ub))` (relative 1.5e-8 on the differing entries)
## called a [1, 1 + 1e-9] operand a point and returned finite bounds for two
## narrow intervals where CVXPY returns unbounded. NaN entries are never equal,
## as in numpy (a stored NaN difference counts as nonzero). Claims-audit C1002,
## decision 14, CVXR 1.9.2.9152.
.bnd_is_point <- function(lb, ub) {
  if (!identical(dim(lb), dim(ub)) || length(lb) != length(ub)) return(FALSE)
  if (inherits(lb, "Matrix") && inherits(ub, "Matrix")) {
    return(Matrix::nnzero(lb - ub, na.counted = TRUE) == 0L)
  }
  isTRUE(all(lb == ub))
}

## CVXPY SOURCE: bounds.py:747-806 (matrix product x %*% y).
## Exact split formula when one operand is a point (lb == ub); otherwise
## unbounded. The point side may be a sparse constant's value and stays
## sparse through pmax/pmin (`_safe_maximum`/`_safe_minimum`, bounds.py:53-76);
## the interval side is dense, and so is the product (a dgeMatrix when one
## factor is sparse, returned as a base matrix). Sweep item 4, pass B.
matmul_bounds <- function(lb1, ub1, lb2, ub2) {
  lhs_point <- .bnd_is_point(lb1, ub1)
  rhs_point <- .bnd_is_point(lb2, ub2)

  if (lhs_point) {
    lb2 <- .bnd_dense(lb2); ub2 <- .bnd_dense(ub2)
    a_pos <- pmax(lb1, 0); a_neg <- pmin(lb1, 0)
    return(list(.dense_Matrix_to_base(a_pos %*% lb2 + a_neg %*% ub2),
                .dense_Matrix_to_base(a_pos %*% ub2 + a_neg %*% lb2)))
  }
  if (rhs_point) {
    lb1 <- .bnd_dense(lb1); ub1 <- .bnd_dense(ub1)
    b_pos <- pmax(lb2, 0); b_neg <- pmin(lb2, 0)
    return(list(.dense_Matrix_to_base(lb1 %*% b_pos + ub1 %*% b_neg),
                .dense_Matrix_to_base(ub1 %*% b_pos + lb1 %*% b_neg)))
  }
  ## both intervals: no efficient exact formula
  unbounded(c(NROW(lb1), NCOL(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).
## Exact per-coordinate equality (np.array_equal): the previous all.equal() test
## called the index vectors 1e8 and 1e8 + 1 equal. Claims-audit C1004, decision 14.
coords_equal <- function(coords1, coords2) {
  if (length(coords1) != length(coords2)) return(FALSE)
  all(vapply(seq_along(coords1),
             function(i) {
               c1 <- coords1[[i]]
               c2 <- coords2[[i]]
               length(c1) == length(c2) && isTRUE(all(c1 == c2))
             },
             logical(1L)))
}