R/288_reductions_solvers_qp_solvers_osqp_qpif.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: reductions/solvers/qp_solvers/osqp_qpif.py
## OSQP QP solver interface for QP/LP problems
##
## OSQP solves: minimize 0.5 x'Px + q'x  s.t. l <= Ax <= u
## Accepts ONLY Zero (equality) and NonNeg (inequality) constraints.
## Inherits from QpSolver -- uses QpSolver.apply() for sign-correct
## A_eq/b_eq/F_ineq/g_ineq data, then stacks into OSQP format.
##
## Requires osqp >= 1.0.0 (S7 API with @Solve/@Update/@WarmStart,
## polishing/warm_starting settings, positive status codes).


# -- OSQP status map (v1.0+) --------------------------------------------------
## CVXPY SOURCE: osqp_qpif.py lines 28-37
## osqp 1.0.0 uses positive integer status codes.

OSQP_STATUS_MAP <- list(
  "1"  = OPTIMAL,                    # OSQP_SOLVED
  "2"  = OPTIMAL_INACCURATE,         # OSQP_SOLVED_INACCURATE
  "3"  = INFEASIBLE,                 # OSQP_PRIMAL_INFEASIBLE
  "4"  = INFEASIBLE_INACCURATE,      # OSQP_PRIMAL_INFEASIBLE_INACCURATE
  "5"  = UNBOUNDED,                  # OSQP_DUAL_INFEASIBLE
  "6"  = UNBOUNDED_INACCURATE,       # OSQP_DUAL_INFEASIBLE_INACCURATE
  "7"  = USER_LIMIT,                 # OSQP_MAX_ITER_REACHED
  "8"  = USER_LIMIT,                 # OSQP_TIME_LIMIT_REACHED
  "9"  = SOLVER_ERROR,               # OSQP_NON_CVX (was -10 in pre-v1)
  "10" = SOLVER_ERROR,               # OSQP_SIGINT
  "11" = SOLVER_ERROR                # OSQP_UNSOLVED
)

# -- OSQP_QP_Solver class -----------------------------------------------------
## CVXPY SOURCE: osqp_qpif.py lines 13-169
## QP path: QpSolver.apply() produces A_eq, b_eq, F_ineq, g_ineq.
## OSQP needs: l <= [A_eq; F_ineq] * x <= u

OSQP_QP_Solver <- new_class("OSQP_QP_Solver", parent = QpSolver, package = "CVXR",
  constructor = function() {
    if (FALSE) new_object(S7_object())  ## S7 static-check guard
    .fast_new(OSQP_QP_Solver, S7_object(),
      .cache = new.env(parent = emptyenv()),
      MIP_CAPABLE = FALSE,
      BOUNDED_VARIABLES = FALSE,
      PSD_TRIANGLE_KIND = NA_character_,
      PSD_SQRT2_SCALING = NA,
      SUPPORTED_CONSTRAINTS = list(Zero, NonNeg),
      REQUIRES_CONSTR = FALSE
    )
  }
)

method(solver_name, OSQP_QP_Solver) <- function(x) OSQP_SOLVER

# -- .osqp_stack_data ----------------------------------------------------------
## Stack QP data into OSQP format: A = [A_eq; F_ineq], l/u bounds.
## Returns list(P, q, A, l, u, len_eq).

.osqp_stack_data <- function(data) {
  q_vec <- data[["q"]]
  nvars <- length(q_vec)

  ## P matrix (quadratic objective) -- upper-triangular for OSQP.
  ## CVXPY SOURCE: osqp_qpif.py line 151 -- upstream hands the full P to
  ## `solver.setup(P, q, A, lA, uA, ...)` and takes sp.triu(P) only for the
  ## warm-path `Px` update (:126). OSQP reads the upper triangle either way, so
  ## passing triu(P) on both paths here is the same problem.
  P <- if (!is.null(data[[SD_P]])) Matrix::triu(data[[SD_P]]) else NULL

  ## Stack A_eq and F_ineq into combined OSQP constraint matrix
  A_eq <- data[["A_eq"]]
  b_eq <- data[["b_eq"]]
  F_ineq <- data[["F_ineq"]]
  g_ineq <- data[["g_ineq"]]

  len_eq <- nrow(A_eq)
  len_ineq <- nrow(F_ineq)

  if (len_eq > 0L && len_ineq > 0L) {
    A <- rbind(A_eq, F_ineq)
    l <- c(b_eq, rep(-Inf, len_ineq))
    u <- c(b_eq, g_ineq)
  } else if (len_eq > 0L) {
    A <- A_eq
    l <- b_eq
    u <- b_eq
  } else if (len_ineq > 0L) {
    A <- F_ineq
    l <- rep(-Inf, len_ineq)
    u <- g_ineq
  } else {
    A <- Matrix::sparseMatrix(i = integer(0), j = integer(0),
                               x = numeric(0), dims = c(0L, nvars))
    l <- numeric(0)
    u <- numeric(0)
  }

  ## Ensure A is dgCMatrix for osqp
  if (!inherits(A, "dgCMatrix")) {
    A <- methods::as(A, "dgCMatrix")
  }
  ## P also needs to be dgCMatrix if non-NULL
  if (!is.null(P) && !inherits(P, "dgCMatrix")) {
    P <- methods::as(P, "generalMatrix")
  }

  list(P = P, q = q_vec, A = A, l = l, u = u, len_eq = len_eq)
}

# -- solve_via_data ------------------------------------------------------------
## CVXPY SOURCE: osqp_qpif.py lines 92-159
## Receives QP data from QpSolver.apply(): P, q, A_eq, b_eq, F_ineq, g_ineq
## Implements warm-start: caches OSQP model, detects data changes, updates
## only what changed, and warm-starts from the previous solution.

method(solve_via_data, OSQP_QP_Solver) <- function(x, data, warm_start = FALSE, verbose = FALSE,
                                                     solver_opts = list(), ...) {
  .require_solver_package(OSQP_SOLVER)

  dots <- list(...)
  solver_cache <- dots[["solver_cache"]]

  ## Stack QP data into OSQP format
  stacked <- .osqp_stack_data(data)
  P <- stacked$P
  q_vec <- stacked$q
  A <- stacked$A
  l <- stacked$l
  u <- stacked$u
  len_eq <- stacked$len_eq

  ## CVXPY defaults (tighter than OSQP defaults)
  ## CVXPY SOURCE: osqp_qpif.py lines 111-114
  solver_opts[["eps_abs"]]  <- solver_opts[["eps_abs"]]  %||% 1e-5
  solver_opts[["eps_rel"]]  <- solver_opts[["eps_rel"]]  %||% 1e-5
  solver_opts[["max_iter"]] <- as.integer(solver_opts[["max_iter"]] %||% 10000L)

  cache_key <- OSQP_SOLVER
  used_warm <- FALSE

  ## -- Warm path ----------------------------------------------------
  ## CVXPY SOURCE: osqp_qpif.py lines 117-144
  if (warm_start && !is.null(solver_cache) && exists(cache_key, envir = solver_cache)) {
    cached <- get(cache_key, envir = solver_cache)
    old_model  <- cached$model
    old_data   <- cached$data
    old_result <- cached$result

    ## Dimension check: if n or m changed, fall through to cold path.
    ##
    ## Sparsity-pattern check (R-specific guard, as the PIQP interface does at
    ## piqp_qpif.R:143-167): an entry whose parameter is currently 0 leaves A
    ## (.dpp_contract_reshape), and `Update(Ax = ...)` writes values into the
    ## model's old pattern. With an equal nnz that was a silent wrong answer
    ## (measured: p = (1, 0) -> (0, 1) in p1 x1 + p2 x2 <= 1 gave 4.75 for
    ## 5.74, status optimal); with a changed nnz osqp_update_data_mat failed
    ## with error code 2 instead of the cold path taking over. CVXPY's data
    ## comparison is the same as the one below (osqp_qpif.py:124-132); it is
    ## safe there because its A keeps a fixed pattern (ReducedMat,
    ## reductions/utilities.py:113-140).
    same_pattern <- function(new, old) {
      identical(dim(new), dim(old)) && length(new@x) == length(old@x) &&
        identical(new@i, old@i) && identical(new@p, old@p)
    }
    P_pattern_ok <- is.null(P) || is.null(old_data$P) ||
      same_pattern(P, old_data$P)
    if (length(q_vec) == length(old_data$q) && nrow(A) == nrow(old_data$A) &&
        same_pattern(A, old_data$A) && P_pattern_ok) {
      new_args <- list()

      ## Compare q, l, u -- only send changed vectors
      if (!identical(q_vec, old_data$q)) new_args$q <- q_vec
      if (!identical(l, old_data$l))     new_args$l <- l
      if (!identical(u, old_data$u))     new_args$u <- u

      ## Compare P@x (non-zero values) -- if changed, send Px values
      factorizing <- FALSE
      if (!is.null(P) && !is.null(old_data$P)) {
        if (!identical(P@x, old_data$P@x)) {
          new_args$Px <- P@x
          factorizing <- TRUE
        }
      }

      ## Compare A@x (non-zero values) -- if changed, send Ax values
      if (!identical(A@x, old_data$A@x)) {
        new_args$Ax <- A@x
        factorizing <- TRUE
      }

      ## Send incremental updates to OSQP model
      if (length(new_args) > 0L) {
        do.call(old_model@Update, new_args)
      }

      ## Warm-start from previous solution if it was optimal
      old_status <- OSQP_STATUS_MAP[[as.character(old_result$info$status_val)]]
      if (!is.null(old_status) && old_status == OPTIMAL) {
        old_model@WarmStart(x = old_result$x, y = old_result$y)
      }

      ## Polish if factorizing (unless user specified)
      solver_opts[["polishing"]] <- solver_opts[["polishing"]] %||% factorizing

      ## Update solver settings
      pars <- c(list(verbose = verbose), solver_opts)
      old_model@UpdateSettings(newpars = pars)

      model <- old_model
      used_warm <- TRUE
    }
  }

  ## -- Cold path ----------------------------------------------------
  ## CVXPY SOURCE: osqp_qpif.py lines 145-153
  if (!used_warm) {
    solver_opts[["polishing"]] <- solver_opts[["polishing"]] %||% TRUE
    pars <- c(list(verbose = verbose, warm_starting = TRUE), solver_opts)
    ## CVXPY SOURCE: osqp_qpif.py:150-153 -- `try: solver.setup(...) except
    ## Exception as e: raise SolverError(e)`. A setup failure (e.g. the KKT
    ## factorization rejecting a non-convex P) is a SolverError, so a
    ## solver_path attempt falls through to the next solver
    ## (claims-audit decision 19, CVXR 1.9.2.9151).
    model <- tryCatch(
      osqp::osqp(P = P, q = q_vec, A = A, l = l, u = u, pars = pars),
      error = function(e) {
        ## The R osqp message already reads "OSQP setup failed with error
        ## code N: ..."; pass it through as CVXPY passes `e`.
        msg <- conditionMessage(e)
        cli_abort("{msg}", class = "SolverError", parent = e)
      }
    )
  }

  ## -- Solve --------------------------------------------------------
  result <- model@Solve()

  ## -- Cache for future warm-starts ---------------------------------
  ## CVXPY SOURCE: osqp_qpif.py lines 157-158
  if (!is.null(solver_cache)) {
    assign(cache_key, list(model = model, data = stacked, result = result),
           envir = solver_cache)
  }

  ## Store len_eq for dual splitting in reduction_invert
  result$.len_eq <- len_eq
  result
}

# -- reduction_invert ----------------------------------------------------------
## CVXPY SOURCE: osqp_qpif.py lines 61-90
## Dual sign: use raw y (NO negation) -- QpSolver.apply() sign flip handles it.

method(reduction_invert, OSQP_QP_Solver) <- function(x, solution, inverse_data, ...) {
  attr_list <- list()

  ## Map status via integer status_val (v1.0 positive codes)
  ## CVXPY SOURCE: osqp_qpif.py lines 69-70
  status_val <- solution$info$status_val
  status <- OSQP_STATUS_MAP[[as.character(status_val)]]
  if (is.null(status)) status <- SOLVER_ERROR

  ## Timing and the raw result
  ## CVXPY SOURCE: osqp_qpif.py:65-66 -- SOLVE_TIME and EXTRA_STATS = solution.
  ## `.len_eq` is CVXR's carrier for the dual split, not OSQP output.
  if (!is.null(solution$info$run_time))
    attr_list[[RK_SOLVE_TIME]] <- solution$info$run_time
  attr_list[[RK_EXTRA_STATS]] <- .drop_private_fields(solution)

  ## Split a [eq; ineq] dual vector keyed by constraint IDs.
  ## CVXPY SOURCE: osqp_qpif.py extract_duals_from_solution_attribute()
  ## (added in CVXPY v1.9.0 #3228). solution_attribute is solution$y on
  ## success and solution$prim_inf_cert on infeasibility.
  ## Dual sign: use raw values, NO negation -- QpSolver.apply() negated F.
  len_eq <- solution$.len_eq
  .extract_duals <- function(solution_attribute) {
    eq_dual <- if (len_eq > 0L) {
      get_dual_values(
        solution_attribute[seq_len(len_eq)],
        extract_dual_value,
        inverse_data[[SOLVER_EQ_CONSTR]]
      )
    } else {
      list()
    }
    ineq_dual <- if (len_eq < length(solution_attribute)) {
      get_dual_values(
        solution_attribute[(len_eq + 1L):length(solution_attribute)],
        extract_dual_value,
        inverse_data[[SOLVER_NEQ_CONSTR]]
      )
    } else {
      list()
    }
    c(eq_dual, ineq_dual)
  }

  if (status %in% SOLUTION_PRESENT) {
    ## Objective value
    opt_val <- solution$info$obj_val + inverse_data[[SD_OFFSET]]

    ## Primal variables
    primal_vars <- list()
    primal_vars[[as.character(inverse_data[[SOLVER_VAR_ID]])]] <- solution$x
    ## CVXPY SOURCE: osqp_qpif.py:78 -- NUM_ITERS only with a solution.
    if (!is.null(solution$info$iter))
      attr_list[[RK_NUM_ITERS]] <- solution$info$iter

    ## OSQP returns y for [eq; ineq] combined
    dual_vars <- .extract_duals(solution$y)
    Solution(status, opt_val, primal_vars, dual_vars, attr_list)
  } else {
    ## CVXPY v1.9.0 fix: #3228 propagate the infeasibility certificate.
    ## y does NOT hold a certificate when infeasible; OSQP exposes it
    ## separately as prim_inf_cert (also ordered [eq; ineq]).
    dual_vars <- if (!is.null(solution$prim_inf_cert)) {
      .extract_duals(solution$prim_inf_cert)
    } else {
      list()
    }
    failure_solution(status, attr_list, dual_vars)
  }
}

# -- print ---------------------------------------------------------------------

method(print, OSQP_QP_Solver) <- function(x, ...) {
  cat("OSQP_QP_Solver()\n")
  invisible(x)
}