R/250_reductions_dqcp2dcp_dqcp2dcp.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: reductions/dqcp2dcp/dqcp2dcp.py
## Dqcp2Dcp -- reduce DQCP problems to parameterized DCP for bisection

## -- Helper: separate constraints into lazy, real, filtering TRUE/FALSE --
## CVXPY SOURCE: dqcp2dcp.py _get_lazy_and_real_constraints()
## TRUE (redundant) constraints are dropped, FALSE (infeasible) constraints
## are replaced with an explicitly infeasible DCP constraint.
.dqcp_separate_constraints <- function(constraints) {
  n <- length(constraints)
  lazy_chunks <- vector("list", n)
  real_chunks <- vector("list", n)
  nl <- 0L; nr <- 0L
  for (i in seq_len(n)) {
    con <- constraints[[i]]
    if (identical(con, TRUE)) next  # redundant -- always satisfied
    if (identical(con, FALSE)) {
      ## Infeasible: add a contradiction constraint (1 <= 0)
      nr <- nr + 1L
      real_chunks[[nr]] <- list(Inequality(Constant(1), Constant(0)))
      next
    }
    if (is.function(con)) {
      nl <- nl + 1L
      lazy_chunks[[nl]] <- list(con)
    } else {
      nr <- nr + 1L
      real_chunks[[nr]] <- list(con)
    }
  }
  lazy <- unlist(lazy_chunks[seq_len(nl)], recursive = FALSE)
  real <- unlist(real_chunks[seq_len(nr)], recursive = FALSE)
  if (is.null(lazy)) lazy <- list()
  if (is.null(real)) real <- list()
  list(lazy = lazy, real = real)
}

# -- Dqcp2Dcp class -----------------------------------------------
## Inherits from Canonicalization but uses its own _canonicalize_constraint
## CVXPY SOURCE: dqcp2dcp.py class Dqcp2Dcp(Canonicalization)

Dqcp2Dcp <- new_class("Dqcp2Dcp", parent = Canonicalization, package = "CVXR",
  constructor = function() {
    if (FALSE) new_object(S7_object())  ## S7 static-check guard
    .fast_new(Dqcp2Dcp, S7_object(),
      .cache = new.env(parent = emptyenv())
    )
  }
)

## accepts: minimization DQCP problem
## CVXPY SOURCE: dqcp2dcp.py accepts()
method(reduction_accepts, Dqcp2Dcp) <- function(x, problem, ...) {
  .s7_is(problem@objective, Minimize) && is_dqcp(problem)
}

## invert: map solution back to original variables
## CVXPY SOURCE: dqcp2dcp.py invert()
method(reduction_invert, Dqcp2Dcp) <- function(x, solution, inverse_data, ...) {
  pvars <- list()
  for (vid in names(inverse_data@id_map)) {
    if (vid %in% names(solution@primal_vars)) {
      pvars[[vid]] <- solution@primal_vars[[vid]]
    } else {
      pvars[[vid]] <- 0.0
    }
  }
  Solution(solution@status, solution@opt_val, pvars, list(), solution@attr)
}

## apply: the main DQCP reduction
## CVXPY SOURCE: dqcp2dcp.py apply()
method(reduction_apply, Dqcp2Dcp) <- function(x, problem, ...) {
  ## Canonicalize all constraints -- collect chunks, flatten once
  n_cons <- length(problem@constraints)
  constr_chunks <- vector("list", n_cons)
  for (i in seq_len(n_cons)) {
    constr_chunks[[i]] <- .dqcp_canonicalize_constraint(problem@constraints[[i]])
  }
  constraints <- unlist(constr_chunks, recursive = FALSE)
  if (is.null(constraints)) constraints <- list()

  ## Separate lazy (callable), real, and boolean (TRUE/FALSE) constraints
  ## CVXPY SOURCE: dqcp2dcp.py _get_lazy_and_real_constraints()
  ## TRUE = redundant (always satisfied), FALSE = infeasible
  sep1 <- .dqcp_separate_constraints(constraints)
  lazy <- sep1$lazy
  real <- sep1$real

  ## Build feasibility problem
  feas_problem <- Problem(Minimize(0), real)
  feas_problem@.cache$lazy_constraints <- lazy

  ## Create parameter t bounded by objective sign
  objective <- .args(problem@objective)[[1L]]
  if (is_nonneg(objective)) {
    t_param <- Parameter(nonneg = TRUE)
  } else if (is_nonpos(objective)) {
    t_param <- Parameter(nonpos = TRUE)
  } else {
    t_param <- Parameter()
  }

  ## Add objective <= t constraint
  obj_constr <- .dqcp_canonicalize_constraint(Inequality(objective, t_param))
  constraints <- c(constraints, obj_constr)

  ## Re-separate
  sep2 <- .dqcp_separate_constraints(constraints)
  lazy2 <- sep2$lazy
  real2 <- sep2$real

  ## Build parameterized problem
  param_problem <- Problem(Minimize(0), real2)
  param_problem@.cache$lazy_constraints <- lazy2

  ## Store bisection data
  tighten <- .tighten_fns(objective)
  param_problem@.cache$bisection_data <- list(
    feas_problem = feas_problem,
    param = t_param,
    tighten_lower = tighten$lower,
    tighten_upper = tighten$upper
  )

  list(param_problem, InverseData(problem))
}

## Helper: canonicalize a single DQCP constraint
## CVXPY SOURCE: dqcp2dcp.py _canonicalize_constraint()
.dqcp_canonicalize_constraint <- function(constr) {
  ## DCP constraints: standard canonicalization
  ## CVXPY SOURCE: dqcp2dcp.py:167
  ##   canon_constr, aux_constr = self.canonicalize_tree(
  ##       constr, canonicalize_params=False)
  ## canonicalize_params = FALSE keeps the bisection parameter's inverted
  ## sublevel bound as a constant expression instead of expanding it into
  ## epigraph variables; see the note in .canonicalize_expr().
  if (is_dcp(constr)) {
    result <- .canonicalize_tree(constr, canonicalize_params = FALSE)
    canon_constr <- result[[1L]]
    aux_constr <- result[[2L]]
    return(c(list(canon_constr), aux_constr))
  }

  ## Must be Inequality for DQCP non-DCP constraints
  if (!.s7_is(constr, Inequality)) {
    cli_abort("Non-DCP constraint must be {.cls Inequality} in DQCP.")
  }

  lhs <- .args(constr)[[1L]]
  rhs <- .args(constr)[[2L]]

  ## CVXPY SOURCE: dqcp2dcp.py lines 174-182
  ## Taking inverses can yield +/- infinity; handle here.
  ## lhs <= rhs: if lhs == -Inf or rhs == Inf, constraint is redundant (TRUE).
  ## If lhs == Inf or rhs == -Inf, constraint is infeasible (FALSE).
  lhs_val <- tryCatch(value(lhs), error = function(e) NULL)
  rhs_val <- tryCatch(value(rhs), error = function(e) NULL)
  if (!is.null(lhs_val) && is.numeric(lhs_val) &&
      all(is.infinite(lhs_val) & lhs_val < 0)) {
    return(list(TRUE))
  }
  if (!is.null(rhs_val) && is.numeric(rhs_val) &&
      all(is.infinite(rhs_val) & rhs_val > 0)) {
    return(list(TRUE))
  }
  if (!is.null(lhs_val) && is.numeric(lhs_val) &&
      any(is.infinite(lhs_val) & lhs_val > 0)) {
    return(list(FALSE))
  }
  if (!is.null(rhs_val) && is.numeric(rhs_val) &&
      any(is.infinite(rhs_val) & rhs_val < 0)) {
    return(list(FALSE))
  }

  ## Short-circuit zero-valued expressions
  if (is_zero(lhs)) {
    return(.dqcp_canonicalize_constraint(Inequality(Constant(0), rhs)))
  }
  if (is_zero(rhs)) {
    return(.dqcp_canonicalize_constraint(Inequality(lhs, Constant(0))))
  }

  ## Case 1: quasiconvex LHS (not convex), constant RHS
  if (is_quasiconvex(lhs) && !is_convex(lhs)) {
    if (.dqcp_invertible(lhs)) {
      rhs_new <- .dqcp_inverse(lhs)(rhs)
      idx <- .non_const_idx(lhs)[1L]
      expr <- .args(lhs)[[idx]]
      if (is_incr(lhs, idx)) {
        return(.dqcp_canonicalize_constraint(Inequality(expr, rhs_new)))
      }
      ## decreasing
      return(.dqcp_canonicalize_constraint(Inequality(rhs_new, expr)))
    } else if (.s7_is(lhs, Maximum) || .s7_is(lhs, MaxEntries)) {
      ## Lower: each arg <= rhs
      result <- list()
      for (arg in .args(lhs)) {
        result <- c(result, .dqcp_canonicalize_constraint(Inequality(arg, rhs)))
      }
      return(result)
    } else {
      ## Sublevel set
      canon_args <- .dqcp_canon_args(lhs)
      lhs_copy <- expr_copy(lhs, args = canon_args[[1L]])
      sublevel <- .dqcp_sublevel(lhs_copy, t = rhs)
      return(c(sublevel, canon_args[[2L]]))
    }
  }

  ## Case 2: constant LHS, quasiconcave RHS
  ## CVXPY SOURCE: dqcp2dcp.py:216-218 -- this is the plain ELSE of case 1, and
  ## the only conditions on it are the two asserts
  ##     assert rhs.is_quasiconcave()
  ##     assert lhs.is_constant()
  ## CVXR made it a second `if` with an extra `&& !is_concave(rhs)` conjunct,
  ## so an RHS that is BOTH quasiconcave and concave fell past it to the
  ## cli_abort below. That is reachable through the recursive rewrites: case 1
  ## applies an inverse and re-enters with a CONSTANT lhs, for which
  ## `is_quasiconvex && !is_convex` is FALSE, and if the surviving rhs is
  ## concave the extra conjunct then rejected a constraint upstream lowers.
  ## Upstream's two asserts. Everything below is unconditional, as there, so
  ## the old trailing cli_abort is gone -- every branch returns.
  if (!is_quasiconcave(rhs)) {
    cli_abort("Cannot canonicalize constraint in DQCP context.")
  }
  if (.dqcp_invertible(rhs)) {
    lhs_new <- .dqcp_inverse(rhs)(lhs)
    idx <- .non_const_idx(rhs)[1L]
    expr <- .args(rhs)[[idx]]
    if (is_incr(rhs, idx)) {
      return(.dqcp_canonicalize_constraint(Inequality(lhs_new, expr)))
    }
    return(.dqcp_canonicalize_constraint(Inequality(expr, lhs_new)))
  } else if (.s7_is(rhs, Minimum) || .s7_is(rhs, MinEntries)) {
    result <- list()
    for (arg in .args(rhs)) {
      result <- c(result, .dqcp_canonicalize_constraint(Inequality(lhs, arg)))
    }
    return(result)
  } else {
    canon_args <- .dqcp_canon_args(rhs)
    rhs_copy <- expr_copy(rhs, args = canon_args[[1L]])
    superlevel <- .dqcp_superlevel(rhs_copy, t = lhs)
    return(c(superlevel, canon_args[[2L]]))
  }
}

## Helper: canonicalize args preserving sign attributes
## CVXPY SOURCE: dqcp2dcp.py _canon_args()
.dqcp_canon_args <- function(expr) {
  n_args <- length(.args(expr))
  canon_args <- vector("list", n_args)
  constr_chunks <- vector("list", n_args)
  for (i in seq_len(n_args)) {
    result <- .dqcp_canonicalize_tree(.args(expr)[[i]])
    canon_arg <- result[[1L]]
    ## Preserve nonneg/nonpos attributes
    if (.s7_is(canon_arg, Variable)) {
      if (is_nonneg(.args(expr)[[i]])) canon_arg@attributes[["nonneg"]] <- TRUE
      else if (is_nonpos(.args(expr)[[i]])) canon_arg@attributes[["nonpos"]] <- TRUE
    }
    canon_args[[i]] <- canon_arg
    constr_chunks[[i]] <- result[[2L]]
  }
  constrs <- unlist(constr_chunks, recursive = FALSE)
  if (is.null(constrs)) constrs <- list()
  list(canon_args, constrs)
}

## Helper: tree walk for DCP sub-expressions
## CVXPY SOURCE: dqcp2dcp.py:122-127
##   def _canonicalize_tree(self, expr):
##       canon_args, constrs = self._canon_args(expr)
##       canon_expr, c = self.canonicalize_expr(
##           expr, canon_args, canonicalize_params=False)
##       constrs += c
##       return canon_expr, constrs
##
## The argument walk MUST go through .dqcp_canon_args(), not an inlined loop:
## that helper re-applies each argument's nonneg/nonpos sign to the Variable it
## canonicalized to (upstream `_canon_args`, dqcp2dcp.py:129-144). CVXR inlined
## the loop here and dropped the sign, so a `Variable(pos = TRUE)` lost its
## nonnegativity inside a DQCP lowering: the stuffed subproblem came out with
## one fewer nonneg row than CVXPY's (measured on
## `minimize sqrt(inv_pos(x)), 0.1 <= x <= 4`: nonneg 3 vs 4, A 6x2 vs 7x2).
## The missing sign leaves the bisection's near-boundary subproblems degenerate,
## so Clarabel returned InsufficientProgress -> SOLVER_ERROR where CVXPY got
## AlmostPrimalInfeasible -> INFEASIBLE_INACCURATE; the loop then perturbed the
## query point instead of tightening the bracket and burned all 100 iterations.
.dqcp_canonicalize_tree <- function(expr) {
  arg_result <- .dqcp_canon_args(expr)
  canon_args <- arg_result[[1L]]
  node_result <- .canonicalize_expr(expr, canon_args, canonicalize_params = FALSE)
  constrs <- c(arg_result[[2L]], node_result[[2L]])
  if (is.null(constrs)) constrs <- list()
  list(node_result[[1L]], constrs)
}