R/251_reductions_solvers_qp_solvers_piqp_qpif.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: reductions/solvers/qp_solvers/piqp_qpif.py
## PIQP QP solver interface for QP/LP problems
##
## PIQP solves: minimize 0.5 x'Px + c'x
##              s.t.  A x = b           (equality)
##                    h_l <= G x <= h_u (inequality, double-sided)
##                    x_l <= x <= x_u   (box constraints)
##
## 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 converts to PIQP format.
##
## R piqp package uses integer status codes (not string enums):
##   1 = solved, -1 = max_iter, -2 = primal_infeasible,
##   -3 = dual_infeasible, -8 = numerical_error


# -- PIQP status map -----------------------------------------------------------
## CVXPY SOURCE: piqp_qpif.py lines 33-37
## R piqp uses integer status codes; CVXPY Python piqp uses string enums.

PIQP_STATUS_MAP <- list(
  "1"   = OPTIMAL,            # PIQP_SOLVED
  "-1"  = USER_LIMIT,         # PIQP_MAX_ITER_REACHED
  "-2"  = INFEASIBLE,         # PIQP_PRIMAL_INFEASIBLE
  "-3"  = UNBOUNDED,          # PIQP_DUAL_INFEASIBLE
  "-8"  = SOLVER_ERROR,       # PIQP_NUMERICAL_ERROR
  "-9"  = SOLVER_ERROR,       # PIQP_UNSOLVED
  "-10" = SOLVER_ERROR        # PIQP_INVALID_SETTINGS
)

# -- PIQP_QP_Solver class ------------------------------------------------------
## CVXPY SOURCE: piqp_qpif.py line 28

PIQP_QP_Solver <- new_class("PIQP_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, PIQP_QP_Solver) <- function(x) PIQP_SOLVER

# -- solve_via_data -------------------------------------------------------------
## CVXPY SOURCE: piqp_qpif.py lines 80-178
## Receives QP data from QpSolver.apply(): P, q, A_eq, b_eq, F_ineq, g_ineq
## Converts to PIQP format and calls piqp model-based API.
## Implements warm-start: caches piqp model, detects data changes, updates
## only what changed.

method(solve_via_data, PIQP_QP_Solver) <- function(x, data, warm_start = FALSE, verbose = FALSE,
                                                     solver_opts = list(), ...) {
  if (!requireNamespace("piqp", quietly = TRUE)) {
    cli_abort("Package {.pkg piqp} is required but not installed.")
  }

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

  ## Extract QP data
  q_vec  <- data[["q"]]
  nvars  <- length(q_vec)
  P      <- data[[SD_P]]
  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)

  ## Ensure sparse matrices are dgCMatrix for piqp
  if (!is.null(P) && !inherits(P, "dgCMatrix")) {
    P <- methods::as(P, "dgCMatrix")
  }
  if (!inherits(A_eq, "dgCMatrix")) {
    A_eq <- methods::as(A_eq, "dgCMatrix")
  }
  if (!inherits(F_ineq, "dgCMatrix")) {
    F_ineq <- methods::as(F_ineq, "dgCMatrix")
  }

  ## Handle NULL P (LP case) -- zero matrix
  if (is.null(P)) {
    P <- Matrix::sparseMatrix(i = integer(0), j = integer(0),
                               dims = c(nvars, nvars))
    P <- methods::as(P, "dgCMatrix")
  }

  ## Backend selection (default: sparse)
  ## CVXPY SOURCE: piqp_qpif.py line 88
  backend <- solver_opts[["backend"]] %||% "sparse"
  solver_opts[["backend"]] <- NULL

  ## Apply standard parameter mappings
  solver_opts[["eps_abs"]]  <- solver_opts[["eps_abs"]]  %||% 1e-8
  solver_opts[["eps_rel"]]  <- solver_opts[["eps_rel"]]  %||% 1e-9
  solver_opts[["max_iter"]] <- as.integer(solver_opts[["max_iter"]] %||% 250L)

  ## Verbose setting
  solver_opts[["verbose"]] <- verbose

  cache_key   <- PIQP_SOLVER
  used_warm   <- FALSE
  structure_changed <- TRUE

  ## -- Warm path ----------------------------------------------------
  ## CVXPY SOURCE: piqp_qpif.py lines 105-142
  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

    structure_changed <- FALSE
    new_args <- list()

    ## Check vector changes: q (c), b_eq (b), g_ineq (h_u)
    if (!identical(q_vec, old_data$q))       new_args[["c"]]   <- q_vec
    if (!identical(b_eq, old_data$b_eq))     new_args[["b"]]   <- b_eq
    if (!identical(g_ineq, old_data$g_ineq)) new_args[["h_u"]] <- g_ineq

    ## Check P matrix sparsity pattern change (structure change)
    if (!identical(dim(P), dim(old_data$P)) ||
        length(P@x) != length(old_data$P@x) ||
        !identical(P@i, old_data$P@i) || !identical(P@p, old_data$P@p)) {
      structure_changed <- TRUE
    } else if (!identical(P@x, old_data$P@x)) {
      new_args[["P"]] <- P
    }

    ## Check A_eq matrix sparsity pattern change
    if (!identical(dim(A_eq), dim(old_data$A_eq)) ||
        length(A_eq@x) != length(old_data$A_eq@x) ||
        !identical(A_eq@i, old_data$A_eq@i) || !identical(A_eq@p, old_data$A_eq@p)) {
      structure_changed <- TRUE
    } else if (!identical(A_eq@x, old_data$A_eq@x)) {
      new_args[["A"]] <- A_eq
    }

    ## Check F_ineq matrix sparsity pattern change
    if (!identical(dim(F_ineq), dim(old_data$F_ineq)) ||
        length(F_ineq@x) != length(old_data$F_ineq@x) ||
        !identical(F_ineq@i, old_data$F_ineq@i) || !identical(F_ineq@p, old_data$F_ineq@p)) {
      structure_changed <- TRUE
    } else if (!identical(F_ineq@x, old_data$F_ineq@x)) {
      new_args[["G"]] <- F_ineq
    }

    if (!structure_changed && length(new_args) > 0L) {
      piqp::update_settings(old_model, solver_opts)
      do.call(update, c(list(old_model), new_args))
      model     <- old_model
      used_warm <- TRUE
    } else if (!structure_changed) {
      ## No changes at all -- just re-solve
      piqp::update_settings(old_model, solver_opts)
      model     <- old_model
      used_warm <- TRUE
    }
  }

  ## -- Cold path ----------------------------------------------------
  ## CVXPY SOURCE: piqp_qpif.py lines 144-169
  if (!used_warm) {
    model <- piqp::piqp(
      P        = P,
      c        = q_vec,
      A        = A_eq,
      b        = b_eq,
      G        = F_ineq,
      h_u      = g_ineq,
      settings = solver_opts,
      backend  = backend
    )
  }

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

  ## -- Cache for future warm-starts ---------------------------------
  if (!is.null(solver_cache)) {
    assign(cache_key, list(
      model = model,
      data  = list(P = P, q = q_vec, A_eq = A_eq, b_eq = b_eq,
                   F_ineq = F_ineq, g_ineq = g_ineq)
    ), envir = solver_cache)
  }

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

# -- reduction_invert -----------------------------------------------------------
## CVXPY SOURCE: piqp_qpif.py lines 46-78
## Dual sign: use raw y and z_u (NO negation) -- QpSolver.apply() sign flip
## handles it, same as OSQP.

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

  ## Map status via integer status code
  ## CVXPY SOURCE: piqp_qpif.py line 51
  status_val <- solution$status
  status <- PIQP_STATUS_MAP[[as.character(status_val)]]
  if (is.null(status)) status <- SOLVER_ERROR

  ## Timing and iteration info
  ## CVXPY SOURCE: piqp_qpif.py line 47
  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$primal_obj + inverse_data[[SD_OFFSET]]

    ## Primal variables
    primal_vars <- list()
    primal_vars[[as.character(inverse_data[[SOLVER_VAR_ID]])]] <- solution$x

    ## Dual variables: PIQP returns y (eq duals) and z_u (ineq duals)
    ## CVXPY SOURCE: piqp_qpif.py lines 62-73
    ## Use raw duals, NO negation -- same convention as OSQP.
    len_eq   <- solution$.len_eq
    len_ineq <- solution$.len_ineq

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

    ineq_dual <- if (len_ineq > 0L) {
      get_dual_values(
        solution$z_u,
        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, PIQP_QP_Solver) <- function(x, ...) {
  cat("PIQP_QP_Solver()\n")
  invisible(x)
}