R/255_reductions_solvers_bisection.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: reductions/solvers/bisection.py
## Bisection solver for DQCP problems

## -- Helper: lower problem (evaluate lazy constraints) ----------
## CVXPY SOURCE: bisection.py _lower_problem():
##   return Problem(Minimize(0),
##                  problem.constraints + [c() for c in problem._lazy_constraints])
##
## CVXPY just concatenates -- TRUE / FALSE constraints flow into the
## Problem() constructor, which (mirrored on R side in problem.R:53-58)
## coerces TRUE -> `Inequality(Constant(0), Constant(1))` and FALSE ->
## `Inequality(Constant(1), Constant(0))`.  The inner solver then
## reports INFEASIBLE if any FALSE was present, and the bisection loop
## narrows via .bisect_infeasible.
##
## An earlier R-side optimization shortcut FALSE -> return(NULL) and
## skipped TRUE outright, which (a) diverged from CVXPY and (b) left a
## latent crash in .bisect_solve(NULL, ...).  Removed in v9512 so the
## R path is byte-equivalent to CVXPY's.
.lower_problem <- function(problem) {
  lazy_fns <- problem@.cache$lazy_constraints
  if (is.null(lazy_fns)) lazy_fns <- list()
  ## R DEVIATION (ADR D_PERF.9, user-approved 2026-08-22): with no lazy
  ## constraints there is nothing to re-evaluate -- Dqcp2Dcp built this
  ## problem as Problem(Minimize(0), real), so it IS its own lowering.
  ## Return it unchanged: repeated solves of the SAME object hit the
  ## parameterized-problem cache instead of recanonicalizing per iteration
  ## (the emitted problem is DPP -- the bisection parameter enters as
  ## param * variable).  CVXPY reconstructs unconditionally
  ## (bisection.py:27-29) and carries a TODO (bisection.py:51-55) wishing
  ## it did not.  When lazy constraints ARE present the reconstruction is
  ## load-bearing (sign-branching callables) and we keep it.
  if (length(lazy_fns) == 0L) return(problem)
  lazy_constrs <- lapply(lazy_fns, function(f) f())
  Problem(Minimize(0), c(problem@constraints, lazy_constrs))
}

## -- Helper: solve sub-problem ----------------------------------
## CVXPY SOURCE: bisection.py _solve()
## CVXPY SOURCE: bisection.py:32-43 -- only DPP warnings are silenced (:35),
## every other sub-solve warning ("Solution may be inaccurate", ...) reaches the
## user; a SolverError is caught, warned about and recorded as SOLVER_ERROR
## (:38-43); any other error propagates. Upstream filters by message
## (`.*DPP.*`); CVXR's one DPP warning carries class DPPWarning
## (solving_chain.R), so it is filtered by class. Through 1.9.2.9170 this
## suppressed every warning and swallowed every error silently (claims audit
## C0723 side findings 1 and 2, decision 20 side finding 1).
.bisect_solve <- function(problem, solver) {
  tryCatch(
    withCallingHandlers(
      psolve(problem, solver = solver),
      DPPWarning = function(w) invokeRestart("muffleWarning")),
    SolverError = function(e) {
      cli_warn("Solver failed, attempting to continue bisection.")
      ## Directly set cache (problem_unpack rejects non-solution status)
      problem@.cache$status <- SOLVER_ERROR
      problem@.cache$solution <- failure_solution(SOLVER_ERROR)
      invisible(NULL)
    }
  )
}

## -- Helper: check infeasible -----------------------------------
## CVXPY SOURCE: bisection.py _infeasible()
.bisect_infeasible <- function(problem) {
  is.null(problem) || status(problem) %in% c(INFEASIBLE, INFEASIBLE_INACCURATE)
}

## -- Find bisection interval -----------------------------------
## CVXPY SOURCE: bisection.py _find_bisection_interval()
.find_bisection_interval <- function(problem, t, solver = NULL,
                                     low = NULL, high = NULL,
                                     max_iters = 100L) {
  if (is.null(low))  low  <- if (is_nonneg(t)) 0 else -1
  if (is.null(high)) high <- if (is_nonpos(t)) 0 else 1

  infeasible_low <- is_nonneg(t)
  feasible_high  <- is_nonpos(t)

  for (iter in seq_len(max_iters)) {
    if (!feasible_high) {
      value(t) <- high
      lowered <- .lower_problem(problem)
      .bisect_solve(lowered, solver)
      if (.bisect_infeasible(lowered)) {
        low <- high
        ## CVXPY SOURCE: bisection.py:73 -- `high * 2 if high != 0 else 1`.
        ## Doubling cannot move an endpoint that sits exactly at 0, so the
        ## bracket never grows and the loop burns all max_iters before failing
        ## with "Unable to find suitable interval". Reachable whenever the
        ## caller supplies an explicit `high = 0` (the default only lands on 0
        ## when `is_nonpos(t)`, which sets feasible_high and skips this branch).
        high <- if (high != 0) high * 2 else 1
        next
      } else if (status(lowered) %in% SOLUTION_PRESENT) {
        feasible_high <- TRUE
      } else {
        ## CVXPY SOURCE: bisection.py:78-79 -- raise error.SolverError.
        cli_abort("Solver failed with status {status(lowered)}", class = "SolverError")
      }
    }

    if (!infeasible_low) {
      value(t) <- low
      lowered <- .lower_problem(problem)
      .bisect_solve(lowered, solver)
      if (.bisect_infeasible(lowered)) {
        infeasible_low <- TRUE
      } else if (status(lowered) %in% SOLUTION_PRESENT) {
        high <- low
        ## CVXPY SOURCE: bisection.py:89 -- `low * 2 if low != 0 else -1`.
        ## Mirror of the `high` guard above.
        low <- if (low != 0) low * 2 else -1
        next
      } else {
        ## CVXPY SOURCE: bisection.py:92-93 -- raise error.SolverError.
        cli_abort("Solver failed with status {status(lowered)}", class = "SolverError")
      }
    }

    if (infeasible_low && feasible_high) {
      return(c(low, high))
    }
  }

  ## CVXPY SOURCE: bisection.py:98-100 -- raise error.SolverError, in these
  ## words (the doubled period is upstream's). Through 1.9.2.9170 the three
  ## interval-search aborts were unclassed (claims audit C0723 side finding 3).
  cli_abort(paste0("Unable to find suitable interval for bisection; ",
                   "your problem may be unbounded.."),
            class = "SolverError")
}

## -- Core bisection loop ---------------------------------------
## CVXPY SOURCE: bisection.py _bisect()
.bisect_loop <- function(problem, solver, t, low, high,
                         tighten_lower, tighten_upper,
                         eps = 1e-6, verbose = FALSE,
                         max_iters = 100L) {
  verbose_freq <- 5L
  soln <- NULL
  ## Lower end of the subinterval the query point is drawn from.  Solver
  ## FAILURES move it toward `high`, the feasible end of the bracket, so a
  ## failing query point is never retried verbatim -- while `low` itself, the
  ## reported bound, stays put.  CVXPY SOURCE: bisection.py:109-112 (1.9.2).
  query_low <- low

  for (i in seq_len(max_iters)) {
    if (!is.null(soln) && (high - low) <= eps) {
      return(list(soln = soln, low = low, high = high))
    }
    query_pt <- (query_low + high) / 2.0
    if (verbose && ((i - 1L) %% verbose_freq == 0L)) {
      cli_inform(c(
        "i" = "(iteration {i - 1L}) lower bound: {format(low, digits = 6)}",
        "i" = "(iteration {i - 1L}) upper bound: {format(high, digits = 6)}",
        "i" = "(iteration {i - 1L}) query point: {format(query_pt, digits = 6)}"
      ))
    }
    value(t) <- query_pt
    lowered <- .lower_problem(problem)
    .bisect_solve(lowered, solver)

    if (.bisect_infeasible(lowered)) {
      if (verbose && ((i - 1L) %% verbose_freq == 0L)) {
        cli_inform(c("i" = "(iteration {i - 1L}) query was infeasible."))
      }
      low <- tighten_lower(query_pt)
      query_low <- low
    } else if (status(lowered) %in% SOLUTION_PRESENT) {
      if (verbose && ((i - 1L) %% verbose_freq == 0L)) {
        cli_inform(c("i" = "(iteration {i - 1L}) query was feasible."))
      }
      soln <- solution(lowered)
      high <- tighten_upper(query_pt)
      query_low <- low
    } else {
      ## Solver FAILED -- neither feasible nor proven infeasible.  This is
      ## typically numerical degeneracy near the optimum of a badly scaled
      ## problem.  CVXPY SOURCE: bisection.py:142-156 (CVXPY 1.9.3).
      ##
      ## CVXR previously did `low <- tighten_lower(query_pt)` here, i.e. it
      ## treated the failure as infeasibility and MOVED THE LOWER BOUND on a
      ## point that was never proven infeasible -- so the reported `low`
      ## could be wrong.  Upstream moves only `query_low`, which perturbs the
      ## next query toward the feasible end without disturbing the bracket.
      ## (Through 1.9.1 upstream instead re-solved the identical failing
      ## point until it hit "Max iters hit during bisection"; CVXR never had
      ## that infinite-retry bug, but avoided it the unsound way.)
      ##
      ## Whether a solution has been seen yet plays no part here: upstream
      ## consults `soln` only in the termination test above.  Through CVXR
      ## 1.9.2.9153 this branch aborted when `soln` was NULL, so a failure
      ## before any feasible query (as on CVXPY's #3407 regression problem)
      ## errored where CVXPY perturbs and goes on to solve (claims audit C0723,
      ## decision 20, CVXR 1.9.2.9154).
      ## CVXPY SOURCE: bisection.py:143-146 -- raise error.SolverError.
      if (query_pt == high) {
        cli_abort(paste0("Solver failed at every query point between the ",
                         "current bisection bounds; unable to make progress."),
                  class = "SolverError")
      }
      ## CVXPY SOURCE: bisection.py:147-151 -- one RuntimeWarning per failed
      ## query, in these words. A CVXPY user sees every one of them: the
      ## iteration number makes each message distinct, so Python's default
      ## once-per-location filter does not collapse them (measured on 1.9.3:
      ## three failures, three warnings, again on a second solve in the same
      ## process). Through 1.9.2.9166 CVXR warned once per bisection on the
      ## contrary premise (claims audit C0723).
      cli_warn(paste0("Solver failed at iteration {i - 1L}, querying a point ",
                      "closer to the feasible end of the interval."))
      if (verbose) {
        cli_inform(c("!" = "(iteration {i - 1L}) solver failed, perturbing query point."))
      }
      query_low <- query_pt
    }
  }

  ## CVXPY SOURCE: bisection.py:157 -- raise error.SolverError.
  cli_abort("Max iters hit during bisection.", class = "SolverError")
}

## -- Main bisection entry point --------------------------------
## CVXPY SOURCE: bisection.py bisect()
bisect <- function(problem, solver = NULL, low = NULL, high = NULL,
                   eps = 1e-6, verbose = FALSE,
                   max_iters = 100L,
                   max_iters_interval_search = 100L) {
  bdata <- problem@.cache$bisection_data
  if (is.null(bdata)) {
    cli_abort("{.fn bisect} only accepts problems emitted by {.cls Dqcp2Dcp}.")
  }

  feas_problem   <- bdata$feas_problem
  t              <- bdata$param
  tighten_lower  <- bdata$tighten_lower
  tighten_upper  <- bdata$tighten_upper

  if (verbose) {
    cli_rule(center = "DQCP bisection")
    cli_inform(c("i" = "Preparing to bisect problem"))
  }

  ## Check feasibility first
  lowered_feas <- .lower_problem(feas_problem)
  .bisect_solve(lowered_feas, solver)
  if (.bisect_infeasible(lowered_feas)) {
    if (verbose) cli_inform(c("!" = "Problem is infeasible."))
    return(failure_solution(INFEASIBLE))
  }

  ## Find bisection interval if not provided
  auto_interval <- is.null(low) || is.null(high)
  if (auto_interval) {
    if (verbose) cli_inform(c("i" = "Finding interval for bisection ..."))
    interval <- .find_bisection_interval(problem, t, solver, low, high,
                                         max_iters_interval_search)
    low  <- interval[1L]
    high <- interval[2L]
  }
  if (verbose) {
    cli_inform(c(
      "i" = "Initial lower bound: {format(low, digits = 6)}",
      "i" = "Initial upper bound: {format(high, digits = 6)}"
    ))
  }

  ## CVXPY SOURCE: bisection.py:212-214 -- straight into the loop, which starts
  ## from no solution (:108).  Through 1.9.2.9171 CVXR first solved once at
  ## `high` after an auto-detected interval and handed the loop that solution.
  ## When the optimum sits exactly at `high` that halved the sub-solves
  ## (25 vs 57) and let a small `max_iters` succeed where CVXPY raises "Max
  ## iters hit"; on every other problem it was one extra solve.  Removed for
  ## parity (claims audit decision 20 side finding 3, repro 32).
  result <- .bisect_loop(problem, solver, t, low, high,
                         tighten_lower, tighten_upper,
                         eps, verbose, max_iters)
  soln <- result$soln
  low  <- result$low
  high <- result$high

  ## Set optimal value to midpoint
  soln <- Solution(
    status      = soln@status,
    opt_val     = (low + high) / 2.0,
    primal_vars = soln@primal_vars,
    dual_vars   = soln@dual_vars,
    attr        = soln@attr
  )

  if (verbose) {
    cli_inform(c(
      "v" = "Bisection completed.",
      "i" = "Lower bound: {format(low, digits = 6)}",
      "i" = "Upper bound: {format(high, digits = 6)}"
    ))
    cli_rule()
  }

  soln
}