R/247_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 40-151
## 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() {
    new_object(S7_object(),
      .cache = new.env(parent = emptyenv()),
      MIP_CAPABLE = FALSE,
      BOUNDED_VARIABLES = FALSE,
      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 93
  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),
                               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 84-151
## 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(), ...) {
  if (!requireNamespace("osqp", quietly = TRUE)) {
    cli_abort("Package {.pkg osqp} (>= 1.0.0) is required but not installed.")
  }

  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 103-106
  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 108-136
  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
    if (length(q_vec) == length(old_data$q) && nrow(A) == nrow(old_data$A)) {
      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 137-145
  if (!used_warm) {
    solver_opts[["polishing"]] <- solver_opts[["polishing"]] %||% TRUE
    pars <- c(list(verbose = verbose, warm_starting = TRUE), solver_opts)
    model <- osqp::osqp(P = P, q = q_vec, A = A, l = l, u = u, pars = pars)
  }

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

  ## -- Cache for future warm-starts ---------------------------------
  ## CVXPY SOURCE: osqp_qpif.py lines 149-150
  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 46-82
## 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 48-55
  status_val <- solution$info$status_val
  status <- OSQP_STATUS_MAP[[as.character(status_val)]]
  if (is.null(status)) status <- SOLVER_ERROR

  ## Timing and iteration info
  if (!is.null(solution$info$run_time))
    attr_list[[RK_SOLVE_TIME]] <- solution$info$run_time
  if (!is.null(solution$info$iter))
    attr_list[[RK_NUM_ITERS]] <- solution$info$iter

  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

    ## Dual variables: OSQP returns y for [eq; ineq] combined
    ## CVXPY SOURCE: osqp_qpif.py -- use raw y, NO negation.
    ## QpSolver.apply() negated F, so OSQP's duals are already correct.
    y <- solution$y
    len_eq <- solution$.len_eq

    eq_dual <- if (len_eq > 0L) {
      get_dual_values(
        y[seq_len(len_eq)],
        extract_dual_value,
        inverse_data[[SOLVER_EQ_CONSTR]]
      )
    } else {
      list()
    }

    ineq_dual <- if (len_eq < length(y)) {
      get_dual_values(
        y[(len_eq + 1L):length(y)],
        extract_dual_value,
        inverse_data[[SOLVER_NEQ_CONSTR]]
      )
    } else {
      list()
    }

    dual_vars <- c(eq_dual, ineq_dual)
    Solution(status, opt_val, primal_vars, dual_vars, attr_list)
  } else {
    failure_solution(status, attr_list)
  }
}

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

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