R/254_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:44-54)
## 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:29-31) and carries a TODO (bisection.py:47-52) 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's problem.solve() does not throw on solver failure -- it sets
## problem.status.  Our psolve() throws via problem_unpack_results.
## Catch errors and store the status so status() works.
.bisect_solve <- function(problem, solver) {
  tryCatch(
    suppressWarnings(psolve(problem, solver = solver)),
    error = function(e) {
      ## 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:74 -- `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 {
        cli_abort("Solver failed with status {.val {status(lowered)}}")
      }
    }

    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:88 -- `low * 2 if low != 0 else -1`.
        ## Mirror of the `high` guard above.
        low <- if (low != 0) low * 2 else -1
        next
      } else {
        cli_abort("Solver failed with status {.val {status(lowered)}}")
      }
    }

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

  cli_abort(c(
    "Unable to find suitable interval for bisection.",
    "i" = "Your problem may be unbounded."
  ))
}

## -- 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, initial_soln = NULL) {
  verbose_freq <- 5L
  soln <- initial_soln
  ## 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
  warned_failure <- FALSE

  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:141-155 (CVXPY 1.9.2).
      ##
      ## 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.)
      if (is.null(soln)) {
        if (verbose) cli_inform(c("!" = "Aborting; the solver failed."))
        cli_abort("Solver failed with status {.val {status(lowered)}}")
      }
      if (query_pt == high) {
        cli_abort(paste0("Solver failed at every query point between the ",
                         "current bisection bounds; unable to make progress."))
      }
      ## Warn ONCE per bisection, not once per failed query.  Upstream calls
      ## warnings.warn() on every failure, but Python's default filter shows a
      ## given warning once per location, so a user sees one message; R's
      ## warning() fires every time.  Emitting per-iteration is the literal
      ## port and the wrong one -- measured, it turned 2 test-suite warnings
      ## into 167.  One warning per bisect() matches what upstream users see.
      if (!warned_failure) {
        cli_warn(paste0("Solver failed at iteration {i - 1L} of bisection; ",
                        "querying points closer to the feasible end of the ",
                        "interval. Results may be less accurate."))
        warned_failure <- TRUE
      }
      if (verbose) {
        cli_inform(c("!" = "(iteration {i - 1L}) solver failed, perturbing query point."))
      }
      query_low <- query_pt
    }
  }

  cli_abort("Max iterations hit during bisection.")
}

## -- 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)}"
    ))
  }

  ## When interval was auto-detected, solve at high to seed initial soln.
  ## The interval search guarantees high is feasible; we need soln
  ## so the bisection convergence check (soln not NULL && gap < eps)
  ## can succeed even when all midpoint queries are infeasible
  ## (this occurs when the optimum lies exactly at the boundary).
  ## Skip when user provides explicit bounds (CVXPY doesn't do this step).
  initial_soln <- NULL
  if (auto_interval) {
    value(t) <- high
    lowered_init <- .lower_problem(problem)
    .bisect_solve(lowered_init, solver)
    if (!.bisect_infeasible(lowered_init) &&
        status(lowered_init) %in% SOLUTION_PRESENT) {
      initial_soln <- solution(lowered_init)
    }
  }

  ## Run bisection
  result <- .bisect_loop(problem, solver, t, low, high,
                         tighten_lower, tighten_upper,
                         eps, verbose, max_iters, initial_soln)
  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
}