R/226_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()
.lower_problem <- function(problem) {
  lazy_fns <- problem@.cache$lazy_constraints
  if (is.null(lazy_fns)) lazy_fns <- list()
  lazy_constrs <- list()
  for (f in lazy_fns) {
    result <- f()
    if (identical(result, FALSE)) return(NULL)  # infeasible
    if (identical(result, TRUE)) next            # trivially satisfied
    lazy_constrs <- c(lazy_constrs, list(result))
  }
  ## Filter TRUE/FALSE from real constraints (can arise from infinity check
  ## in .dqcp_canonicalize_constraint)
  real_constrs <- list()
  for (con in problem@constraints) {
    if (identical(con, TRUE)) next
    if (identical(con, FALSE)) return(NULL)  # infeasible
    real_constrs <- c(real_constrs, list(con))
  }
  Problem(Minimize(0),
          c(real_constrs, 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
        high <- high * 2
        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
        low <- low * 2
        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

  for (i in seq_len(max_iters)) {
    if (!is.null(soln) && (high - low) <= eps) {
      return(list(soln = soln, low = low, high = high))
    }
    query_pt <- (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)
    } 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)
    } else {
      ## Solver failed -- at boundary, this is typically numerical
      ## degeneracy. If we already have a feasible solution, treat
      ## as infeasible and continue narrowing. Otherwise abort.
      if (is.null(soln)) {
        if (verbose) cli_inform(c("!" = "Aborting; the solver failed."))
        cli_abort("Solver failed with status {.val {status(lowered)}}")
      }
      if (verbose && ((i - 1L) %% verbose_freq == 0L)) {
        cli_inform(c("!" = "(iteration {i - 1L}) solver error; treating as infeasible."))
      }
      low <- tighten_lower(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
}