R/252_reductions_solvers_solving_chain.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: reductions/solvers/solving_chain.py
## SolvingChain -- a Chain with a terminal solver reduction
##
## Dual-interface architecture matching CVXPY:
## - SOLVER_MAP_CONIC: conic path solvers (ConicSolver subclasses)
## - SOLVER_MAP_QP: QP path solvers (QpSolver subclasses)
## - _solve_as_qp(): routes QP problems to QP solvers
## - construct_solving_chain(): builds the full reduction chain


# -- Solver registries --------------------------------------------
## Two registries matching CVXPY defines.py:
## SOLVER_MAP_CONIC for conic path, SOLVER_MAP_QP for QP path.

## QP solver preference order (matches CVXPY defines.py QP_SOLVERS)
## Lazy-init via delayedAssign: solver objects constructed on first access.
SOLVER_MAP_QP <- new.env(hash = TRUE, parent = emptyenv())
delayedAssign(OSQP_SOLVER,   OSQP_QP_Solver(),   assign.env = SOLVER_MAP_QP)
delayedAssign(GUROBI_SOLVER, Gurobi_QP_Solver(),  assign.env = SOLVER_MAP_QP)
delayedAssign(CPLEX_SOLVER,  CPLEX_QP_Solver(),   assign.env = SOLVER_MAP_QP)
delayedAssign(HIGHS_SOLVER,  HiGHS_QP_Solver(),   assign.env = SOLVER_MAP_QP)
delayedAssign(PIQP_SOLVER,   PIQP_QP_Solver(),    assign.env = SOLVER_MAP_QP)

QP_SOLVER_PREFERENCE <- c(OSQP_SOLVER, GUROBI_SOLVER, CPLEX_SOLVER, HIGHS_SOLVER, PIQP_SOLVER)

## Conic solver preference order (matches CVXPY defines.py CONIC_SOLVERS)
## Lazy-init via delayedAssign: solver objects constructed on first access.
SOLVER_MAP_CONIC <- new.env(hash = TRUE, parent = emptyenv())
delayedAssign(CLARABEL_SOLVER, Clarabel_Solver(),     assign.env = SOLVER_MAP_CONIC)
delayedAssign(SCS_SOLVER,      SCS_Solver(),           assign.env = SOLVER_MAP_CONIC)
delayedAssign(MOSEK_SOLVER,    Mosek_Solver(),         assign.env = SOLVER_MAP_CONIC)
delayedAssign(GUROBI_SOLVER,   Gurobi_Conic_Solver(),  assign.env = SOLVER_MAP_CONIC)
delayedAssign(HIGHS_SOLVER,    HiGHS_Conic_Solver(),   assign.env = SOLVER_MAP_CONIC)
delayedAssign(GLPK_SOLVER,     GLPK_Solver(),          assign.env = SOLVER_MAP_CONIC)
delayedAssign(GLPK_MI_SOLVER,  GLPK_MI_Solver(),       assign.env = SOLVER_MAP_CONIC)
delayedAssign(ECOS_SOLVER,     ECOS_Solver(),           assign.env = SOLVER_MAP_CONIC)
delayedAssign(ECOS_BB_SOLVER,  ECOS_BB_Solver(),        assign.env = SOLVER_MAP_CONIC)
delayedAssign(CVXOPT_SOLVER,   CVXOPT_Solver(),         assign.env = SOLVER_MAP_CONIC)

CONIC_SOLVER_PREFERENCE <- c(MOSEK_SOLVER, CLARABEL_SOLVER, SCS_SOLVER,
                              ECOS_SOLVER, GUROBI_SOLVER, GLPK_SOLVER,
                              GLPK_MI_SOLVER, CVXOPT_SOLVER, HIGHS_SOLVER,
                              ECOS_BB_SOLVER)

## Package names for solver availability checks
.SOLVER_PACKAGES <- list(
  CLARABEL = "clarabel",
  SCS      = "scs",
  OSQP     = "osqp",
  HIGHS    = "highs",
  MOSEK    = "Rmosek",
  GUROBI   = "gurobi",
  GLPK     = "Rglpk",
  GLPK_MI  = "Rglpk",
  ECOS     = "ECOSolveR",
  ECOS_BB  = "ECOSolveR",
  CPLEX    = "Rcplex",
  CVXOPT   = "cccp",
  PIQP     = "piqp"
)

## Check if a solver package is installed and not excluded
.solver_package_available <- function(solver_name) {
  if (solver_name %in% .cvxr_env$excluded_solvers) return(FALSE)
  pkg <- .SOLVER_PACKAGES[[solver_name]]
  if (is.null(pkg)) return(FALSE)
  requireNamespace(pkg, quietly = TRUE)
}

# -- SolvingChain class -------------------------------------------
## CVXPY SOURCE: solving_chain.py lines 20-80

SolvingChain <- new_class("SolvingChain", parent = Chain, package = "CVXR",
  properties = list(
    solver = class_any  # the terminal solver instance
  ),
  constructor = function(reductions = list()) {
    solver_inst <- if (length(reductions) > 0L) {
      reductions[[length(reductions)]]
    } else {
      NULL
    }
    new_object(S7_object(),
      .cache     = new.env(parent = emptyenv()),
      reductions = reductions,
      solver     = solver_inst
    )
  }
)

method(print, SolvingChain) <- function(x, ...) {
  names <- vapply(x@reductions, function(r) short_class_name(r),
                  character(1L))
  cat(sprintf("SolvingChain(%s)\n", paste(names, collapse = " -> ")))
  invisible(x)
}

# -- _is_lp() -----------------------------------------------------
## CVXPY SOURCE: solving_chain.py lines 89-98
## Checks if problem is a linear program.

.is_lp <- function(problem) {
  ## Check constraints: all must be Equality/Zero or have PWL args
  for (con in problem@constraints) {
    if (!(S7_inherits(con, Equality) || S7_inherits(con, Zero))) {
      if (!is_pwl(con@args[[1L]])) return(FALSE)
    }
  }
  ## Check variables: no PSD/NSD
  for (v in variables(problem)) {
    attrs <- v@attributes
    if (isTRUE(attrs$PSD) || isTRUE(attrs$NSD)) return(FALSE)
  }
  ## Objective must be PWL and problem DCP
  is_dcp(problem) && is_pwl(problem@objective@args[[1L]])
}

# -- _solve_as_qp() ----------------------------------------------
## CVXPY SOURCE: solving_chain.py lines 101-113
## Decides if we should use the QP path.

.solve_as_qp <- function(problem, candidates) {
  ## LP: prefer conic path (OSQP is slow at LPs)
  if (.is_lp(problem) &&
      length(setdiff(candidates$conic_solvers, candidates$qp_solvers)) > 0L) {
    return(FALSE)
  }
  ## QP: use QP path if QP solvers available and problem is QP
  length(candidates$qp_solvers) > 0L && is_qp(problem)
}

# -- .build_candidates() -----------------------------------------
## CVXPY SOURCE: solving_chain.py construct() (candidate selection logic)
## Builds candidate lists of available solvers for both QP and conic paths.

.build_candidates <- function(problem, solver = NULL) {
  is_mip <- is_mixed_integer(problem)

  if (!is.null(solver)) {
    ## User specified a solver -- find it in either map
    qp_solvers <- character(0)
    conic_solvers <- character(0)
    if (!is.null(SOLVER_MAP_QP[[solver]])) {
      inst <- SOLVER_MAP_QP[[solver]]
      if (!is_mip || inst@MIP_CAPABLE) {
        qp_solvers <- solver
      }
    }
    if (!is.null(SOLVER_MAP_CONIC[[solver]])) {
      inst <- SOLVER_MAP_CONIC[[solver]]
      if (!is_mip || inst@MIP_CAPABLE) {
        conic_solvers <- solver
      }
    }
    if (length(qp_solvers) == 0L && length(conic_solvers) == 0L) {
      if (is_mip) {
        ## Suggest MIP-capable solvers (check both maps for MIP_CAPABLE)
        all_s <- unique(c(QP_SOLVER_PREFERENCE, CONIC_SOLVER_PREFERENCE))
        mip_solvers <- character(length(all_s))
        mi <- 0L
        for (s in all_s) {
          if (!.solver_package_available(s)) next
          qp_inst <- SOLVER_MAP_QP[[s]]
          conic_inst <- SOLVER_MAP_CONIC[[s]]
          if ((!is.null(qp_inst) && qp_inst@MIP_CAPABLE) ||
              (!is.null(conic_inst) && conic_inst@MIP_CAPABLE)) {
            mi <- mi + 1L
            mip_solvers[[mi]] <- s
          }
        }
        mip_solvers <- mip_solvers[seq_len(mi)]
        suggest <- if (length(mip_solvers) > 0L) {
          paste0("Use ", paste0("{.val ", mip_solvers, "}", collapse = " or "), " instead.")
        } else {
          "Install a MIP-capable solver such as {.val HIGHS} or {.val GUROBI}."
        }
        cli_abort(c(
          "Solver {.val {solver}} does not support mixed-integer problems.",
          "i" = suggest
        ))
      } else {
        cli_abort("Solver {.val {solver}} cannot handle this problem type.")
      }
    }
  } else {
    ## Auto-select: filter by installed + capable
    qp_solvers <- character(length(QP_SOLVER_PREFERENCE))
    qi <- 0L
    for (s in QP_SOLVER_PREFERENCE) {
      inst <- SOLVER_MAP_QP[[s]]
      if (.solver_package_available(s) && (!is_mip || inst@MIP_CAPABLE)) {
        qi <- qi + 1L
        qp_solvers[[qi]] <- s
      }
    }
    qp_solvers <- qp_solvers[seq_len(qi)]
    conic_solvers <- character(length(CONIC_SOLVER_PREFERENCE))
    ci <- 0L
    for (s in CONIC_SOLVER_PREFERENCE) {
      inst <- SOLVER_MAP_CONIC[[s]]
      if (.solver_package_available(s) && (!is_mip || inst@MIP_CAPABLE)) {
        ci <- ci + 1L
        conic_solvers[[ci]] <- s
      }
    }
    conic_solvers <- conic_solvers[seq_len(ci)]
    if (length(qp_solvers) == 0L && length(conic_solvers) == 0L) {
      cli_abort("No installed solver can handle this problem.")
    }
  }

  list(qp_solvers = qp_solvers, conic_solvers = conic_solvers)
}

# -- .required_cone_types() -------------------------------------
## Predict which cone types a problem will require after Dcp2Cone.
## This is a heuristic used for solver selection -- exact cones are
## determined after canonicalization, but we can predict from atoms.

.required_cone_types <- function(problem) {
  ## Returns a list of S7 class objects for the "advanced" cone types required.
  ## Base types (Zero, NonNeg) are always supported -- not returned.
  cones <- list()
  cone_seen <- new.env(hash = TRUE, parent = emptyenv())

  .add_cone <- function(cls) {
    nm <- cls@name
    if (!exists(nm, envir = cone_seen, inherits = FALSE)) {
      assign(nm, TRUE, envir = cone_seen)
      cones[[length(cones) + 1L]] <<- cls
    }
  }

  ## Walk constraint tree to detect atom types that require specific cones.
  ## Classifications match CVXPY's SOC_ATOMS, PSD_ATOMS, EXP_ATOMS,
  ## POWCONE_ATOMS, POWCONE_ND_ATOMS (cvxpy/atoms/__init__.py).
  ## Approx subclasses checked BEFORE exact parents (S7_inherits is TRUE for both).
  .check_expr <- function(expr) {
    ## Direct cone constraint types
    if (S7_inherits(expr, SOC))       return(SOC)
    if (S7_inherits(expr, PSD))       return(PSD)
    if (S7_inherits(expr, ExpCone))   return(ExpCone)
    if (S7_inherits(expr, PowCone3D)) return(PowCone3D)
    if (S7_inherits(expr, PowConeND)) return(PowConeND)
    ## SOC_ATOMS: Approx variants (SOC via rational approx), QuadForm, QuadOverLin, Huber
    if (S7_inherits(expr, PnormApprox))  return(SOC)
    if (S7_inherits(expr, PowerApprox))  return(SOC)
    if (S7_inherits(expr, GeoMeanApprox)) return(SOC)
    if (S7_inherits(expr, QuadOverLin) || S7_inherits(expr, QuadForm) ||
        S7_inherits(expr, SymbolicQuadForm)) return(SOC)
    if (S7_inherits(expr, Huber)) return(SOC)
    ## POWCONE_ATOMS: exact Pnorm, Power (PowCone3D)
    if (S7_inherits(expr, Pnorm))  return(PowCone3D)
    if (S7_inherits(expr, Power))  return(PowCone3D)
    ## POWCONE_ND_ATOMS: exact GeoMean (PowConeND)
    if (S7_inherits(expr, GeoMean)) return(PowConeND)
    ## PSD_ATOMS
    if (S7_inherits(expr, SigmaMax) || S7_inherits(expr, NormNuc)) return(PSD)
    if (S7_inherits(expr, MatrixFrac) || S7_inherits(expr, TrInv)) return(PSD)
    if (S7_inherits(expr, LambdaMax) || S7_inherits(expr, LambdaSumLargest)) return(PSD)
    if (S7_inherits(expr, LogDet)) return(PSD)
    if (S7_inherits(expr, ConditionNumber)) return(PSD)
    ## EXP_ATOMS
    if (S7_inherits(expr, LogSumExp) || S7_inherits(expr, Exp) ||
        S7_inherits(expr, Log) || S7_inherits(expr, Entr) ||
        S7_inherits(expr, KlDiv) || S7_inherits(expr, RelEntr) ||
        S7_inherits(expr, Logistic) || S7_inherits(expr, Xexp) ||
        S7_inherits(expr, Log1p)) return(ExpCone)
    NULL
  }

  ## Walk all constraints and the objective
  all_exprs <- c(problem@constraints, list(problem@objective@args[[1L]]))
  for (expr_root in all_exprs) {
    ## BFS over expression tree -- index-based to avoid O(n^2) c()/queue[-1L]
    queue <- list(expr_root)
    qi <- 1L
    while (qi <= length(queue)) {
      e <- queue[[qi]]
      qi <- qi + 1L
      cone <- .check_expr(e)
      if (!is.null(cone)) .add_cone(cone)
      if (S7_inherits(e, Expression) && length(e@args) > 0L) {
        for (a in e@args) queue[[length(queue) + 1L]] <- a
      }
    }
  }

  cones  # list of S7 class objects (may be empty)
}

# -- .solver_supports_cones() ----------------------------------
## Check if a conic solver supports the required cone types.

.solver_supports_cones <- function(solver_inst, required_cones) {
  supported <- solver_inst@SUPPORTED_CONSTRAINTS
  unsupported <- Filter(function(cone) {
    !any(vapply(supported, identical, logical(1L), cone))
  }, required_cones)
  list(ok = length(unsupported) == 0L, unsupported = unsupported)
}

# -- construct_solving_chain --------------------------------------
## CVXPY SOURCE: solving_chain.py construct() (simplified)
## Builds the reduction chain for a given problem and solver.

construct_solving_chain <- function(problem, solver = NULL, gp = FALSE) {
  ## CVXPY SOURCE: solving_chain.py line 245-246
  ## Zero-variable problems: use ConstantSolver (evaluates constraints directly)
  if (length(variables(problem)) == 0L) {
    return(SolvingChain(reductions = list(ConstantSolver())))
  }

  reductions <- list()

  ## 0. DPP-aware parameter handling
  ## CVXPY SOURCE: solving_chain.py lines 250-267
  has_params <- length(parameters(problem)) > 0L
  obj_expr <- problem@objective@args[[1L]]
  ## DPP context depends on gp flag (CVXPY: dpp_context = 'dgp' if gp else 'dcp')
  dpp_ok <- if (gp) .is_dgp_dpp(problem) else is_dpp(problem)
  qp_with_params <- has_params && !gp &&
    is_quadratic(obj_expr) && !is_affine(obj_expr)
  ## Complex parameters need EvalParams before Complex2Real because:
  ## (1) Imag_/Real_ atoms lack graph_implementation (can't enter LinOp tree)
  ## (2) R's Matrix package doesn't support complex sparse matrices (zgeMatrix)
  ## (3) Complex DPP is not yet supported (all complex-dpp-parity tests skipped)
  has_complex_params <- has_params &&
    any(vapply(parameters(problem), function(p) is_complex(p) || is_imag(p),
               logical(1L)))
  if (has_params && (!dpp_ok || qp_with_params || has_complex_params)) {
    reductions <- c(reductions, list(EvalParams()))
  }

  ## 1. Complex2Real: if any leaf is complex, reduce to real
  if (complex2real_accepts(problem)) {
    reductions <- c(reductions, list(Complex2Real()))
  }

  ## 2. DGP path: insert Dgp2Dcp reduction (G7)
  ## Chain order: [EvalParams] -> [Complex2Real] -> [Dgp2Dcp] -> [FlipObjective]
  ##              -> [Dcp2Cone] -> [CvxAttr2Constr] -> [ConeMatrixStuffing] -> [Solver]
  if (gp) {
    if (!is_dgp(problem)) {
      cli_abort(c(
        "Problem is not DGP compliant.",
        "i" = "Remove {.code gp = TRUE} or reformulate as a geometric program."
      ))
    }
    reductions <- c(reductions, list(Dgp2Dcp()))
  } else if (!is_dcp(problem)) {
    ## Check if it might be DGP and suggest gp=TRUE
    if (is_dgp(problem)) {
      cli_abort(c(
        "Problem is not DCP compliant.",
        "i" = "However, the problem is DGP. Try {.code psolve(problem, gp = TRUE)}."
      ))
    }
    ## Check if it might be DQCP and suggest qcp=TRUE
    if (is_dqcp(problem)) {
      cli_abort(c(
        "Problem is not DCP compliant.",
        "i" = "However, the problem is DQCP. Try {.code psolve(problem, qcp = TRUE)}."
      ))
    }
    cli_abort("Problem is not DCP compliant.")
  }

  ## 3. Flip objective if Maximize -> Minimize
  if (S7_inherits(problem@objective, Maximize)) {
    reductions <- c(reductions, list(FlipObjective()))
  }

  ## 3a. FiniteSet -> MIP reduction (used by both QP and conic pathways)
  ## CVXPY SOURCE: solving_chain.py lines 178-182
  if (any(vapply(problem@constraints, function(c) S7_inherits(c, FiniteSet), logical(1)))) {
    reductions <- c(reductions, list(Valinvec2mixedint()))
  }

  ## 3b. Build candidate lists
  candidates <- .build_candidates(problem, solver)

  ## 3b. Early MIQP rejection: MIP + QP (not LP) but no MIP-capable QP solver
  ## Note: is_qp() returns TRUE for LP (LP is a subset of QP), so we must
  ## exclude LP using has_quadratic_term to only detect actual QP problems.
  if (is_mixed_integer(problem) && is_qp(problem) &&
      has_quadratic_term(problem@objective@args[[1L]])) {
    ## This is a true MIQP (mixed-integer + quadratic objective)
    has_miqp <- FALSE
    for (s in candidates$qp_solvers) {
      inst <- SOLVER_MAP_QP[[s]]
      if (inst@MIP_CAPABLE) { has_miqp <- TRUE; break }
    }
    if (!has_miqp) {
      miqp_solvers <- c()
      for (s in QP_SOLVER_PREFERENCE) {
        inst <- SOLVER_MAP_QP[[s]]
        if (inst@MIP_CAPABLE) miqp_solvers <- c(miqp_solvers, s)
      }
      suggest <- if (length(miqp_solvers) > 0L) {
        paste0("Try solver: ", paste(miqp_solvers, collapse = ", "), ".")
      } else {
        "No MIQP-capable solver is available. Use LP constraints with integer variables, or remove integer/boolean attributes for QP problems."
      }
      cli_abort(c(
        "Mixed-integer QP (MIQP) is not supported by the selected solver.",
        "i" = suggest
      ))
    }
  }

  ## 4. Route: QP path or conic path
  if (.solve_as_qp(problem, candidates)) {
    ## QP path
    solver_name_sel <- candidates$qp_solvers[1L]
    solver_inst <- SOLVER_MAP_QP[[solver_name_sel]]
    quad_obj <- has_quadratic_term(problem@objective@args[[1L]])
    reductions <- c(reductions, list(
      Dcp2Cone(quad_obj = quad_obj),
      CvxAttr2Constr(),
      ConeMatrixStuffing(quad_obj = quad_obj),
      solver_inst
    ))
  } else {
    ## Conic path -- find a solver that accepts the problem's cones
    solver_name_sel <- NULL
    solver_inst <- NULL
    required_cones <- .required_cone_types(problem)

    ## Try each conic solver in preference order
    for (s in candidates$conic_solvers) {
      inst <- SOLVER_MAP_CONIC[[s]]
      check <- .solver_supports_cones(inst, required_cones)
      if (check$ok) {
        solver_name_sel <- s
        solver_inst <- inst
        break
      }
    }

    if (is.null(solver_inst)) {
      ## Build informative error message
      if (length(required_cones) > 0L) {
        cone_str <- paste(vapply(required_cones, function(c) c@name, character(1L)),
                          collapse = ", ")
        if (!is.null(solver)) {
          cli_abort(c(
            "Solver {.val {solver}} does not support the required cone types: {cone_str}.",
            "i" = "This problem requires {cone_str} cones. Choose a solver that supports them."
          ))
        } else {
          cli_abort(c(
            "No installed solver supports the required cone types: {cone_str}.",
            "i" = "Install a solver that supports {cone_str} (e.g., Clarabel, SCS, or MOSEK)."
          ))
        }
      } else {
        cli_abort("No installed conic solver can handle this problem.")
      }
    }

    quad_obj <- has_quadratic_term(problem@objective@args[[1L]])
    reductions <- c(reductions, list(
      Dcp2Cone(quad_obj = quad_obj),
      CvxAttr2Constr(),
      ConeMatrixStuffing(quad_obj = quad_obj),
      solver_inst
    ))
  }

  SolvingChain(reductions = reductions)
}

# -- solve_via_data (SolvingChain) ---------------------------------
## Delegates to the terminal solver, managing the solver_cache.
## When `problem` is supplied, uses the problem's solver_cache
## (for warm-start state persistence across repeated solves).

#' @rdname solve_via_data
#' @name solve_via_data
#' @export
method(solve_via_data, SolvingChain) <- function(x, data, warm_start = FALSE,
                                                  verbose = FALSE,
                                                  solver_opts = list(), ...,
                                                  problem = NULL) {
  solver_cache <- if (!is.null(problem)) {
    .get_solver_cache(problem)
  } else {
    new.env(parent = emptyenv())
  }
  solve_via_data(x@solver, data, warm_start, verbose, solver_opts,
                 solver_cache = solver_cache)
}