R/247_reductions_dqcp2dcp_inverse.R

In CVXR: Disciplined Convex Optimization

#####
## DO NOT EDIT THIS FILE!! EDIT THE SOURCE INSTEAD: rsrc_tree/reductions/dqcp2dcp/inverse.R
#####

## CVXPY SOURCE: reductions/dqcp2dcp/inverse.py
## Invertibility checks and inverse functions for DQCP atoms

## Compute the inverse of an invertible expression
## Returns a closure function(t) -> expression
## CVXPY SOURCE: inverse.py inverse()
.dqcp_inverse <- function(expr) {
  if (.s7_is(expr, Ceil)) {
    return(function(t) Floor(t))
  } else if (.s7_is(expr, Floor)) {
    return(function(t) Ceil(t))
  } else if (.s7_is(expr, NegExpression)) {
    return(function(t) -t)
  } else if (.s7_is(expr, Exp)) {
    return(function(t) {
      if (is_nonneg(t)) Log(t) else -Inf
    })
  } else if (.s7_is(expr, Log)) {
    return(function(t) Exp(t))
  } else if (.s7_is(expr, Log1p)) {
    return(function(t) Exp(t) - 1)
  } else if (.s7_is(expr, Logistic)) {
    return(function(t) {
      if (is_nonneg(t)) Log(Exp(t) - 1) else -Inf
    })
  } else if (.s7_is(expr, Power)) {
    p_val <- expr@p_orig
    return(function(t) {
      if (p_val == 1) return(t)
      ## CVXPY v1.9.0 fix: #3180 power_inv sentinel -- the infeasible
      ## sentinel depends on p's sign: p > 0 (increasing) needs -Inf,
      ## p < 0 (decreasing) needs +Inf, since the DQCP code flips the
      ## constraint direction for decreasing functions. Returning +Inf
      ## unconditionally made power(ceil(x), 2) <= -5 report OPTIMAL
      ## instead of INFEASIBLE.
      if (is_nonneg(t)) return(power(t, 1 / p_val))
      if (p_val > 0) -Inf else Inf
    })
  } else if (.s7_is(expr, Multiply)) {
    ## Multiply (ELEMENTWISE), not MulExpression. Upstream keys this branch on
    ## `type(expr) == atoms.multiply` (inverse.py:56) -- an exact-type test, so
    ## a matrix-product MulExpression reaches none of the branches and falls
    ## through to `raise ValueError`. CVXR's Multiply extends MulExpression
    ## (binary_operators.R:190), so testing the parent swept matrix products in
    ## and handed back t / A -- ELEMENTWISE division by the matrix, where the
    ## true inverse needs a solve. See notes/audit/completeness_ledger.md
    ## finding 6 and notes/audit/repro/03_dqcp_matmul.{py,R}.
    if (is_constant(.args(expr)[[1L]])) {
      const <- .args(expr)[[1L]]
    } else {
      const <- .args(expr)[[2L]]
    }
    return(function(t) DivExpression(t, const))
  } else if (.s7_is(expr, DivExpression)) {
    if (is_constant(.args(expr)[[1L]])) {
      const <- .args(expr)[[1L]]
      return(function(t) DivExpression(const, t))
    } else {
      const <- .args(expr)[[2L]]
      return(function(t) MulExpression(const, t))
    }
  } else if (.s7_is(expr, AddExpression)) {
    if (is_constant(.args(expr)[[1L]])) {
      const <- .args(expr)[[1L]]
    } else {
      const <- .args(expr)[[2L]]
    }
    return(function(t) t - const)
  } else if (.s7_is(expr, Abs)) {
    arg <- .args(expr)[[1L]]
    if (is_nonneg(arg)) {
      return(function(t) t)
    } else if (is_nonpos(arg)) {
      return(function(t) -t)
    } else {
      cli_abort("Sign of argument to {.cls Abs} must be known for DQCP inversion.")
    }
  } else if (.s7_is(expr, SumEntries) || .s7_is(expr, Cumsum)) {
    return(function(t) t)
  } else {
    cls <- class(expr)[[1L]]
    cli_abort("Cannot compute inverse of {.cls {cls}}.")
  }
}

## Check if an expression is invertible
## CVXPY SOURCE: inverse.py invertible()
.dqcp_invertible <- function(expr) {
  ## Multiply (ELEMENTWISE), not MulExpression -- upstream tests
  ## `isinstance(expr, atoms.multiply)` (inverse.py:93), and multiply is the
  ## NARROW class. Multiply is the only subclass of MulExpression in CVXR, so
  ## this is exactly upstream's semantics. See the note in .dqcp_inverse above.
  if (.s7_is(expr, Multiply) ||
      .s7_is(expr, DivExpression) ||
      .s7_is(expr, AddExpression)) {
    return(length(.non_const_idx(expr)) == 1L)
  } else if (.s7_is(expr, SumEntries) || .s7_is(expr, Cumsum)) {
    return(.is_real_fn(expr))
  } else {
    ## Always-invertible atom types
    ## CVXPY SOURCE: inverse.py INVERTIBLE set
    return(
      .s7_is(expr, Ceil) ||
      .s7_is(expr, Floor) ||
      .s7_is(expr, NegExpression) ||
      .s7_is(expr, Exp) ||
      .s7_is(expr, Log) ||
      .s7_is(expr, Log1p) ||
      .s7_is(expr, Logistic) ||
      .s7_is(expr, Power) ||
      .s7_is(expr, Abs)
    )
  }
}