R/292_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)
delayedAssign(XPRESS_SOLVER, XPRESS_QP_Solver(),  assign.env = SOLVER_MAP_QP)

QP_SOLVER_PREFERENCE <- c(OSQP_SOLVER, GUROBI_SOLVER, XPRESS_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(CPLEX_SOLVER,    CPLEX_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)
delayedAssign(SCIP_SOLVER,     SCIP_Solver(),           assign.env = SOLVER_MAP_CONIC)
delayedAssign(XPRESS_SOLVER,   XPRESS_Conic_Solver(),   assign.env = SOLVER_MAP_CONIC)
delayedAssign(DIFFCP_SOLVER,   DIFFCP_Solver(),         assign.env = SOLVER_MAP_CONIC)

CONIC_SOLVER_PREFERENCE <- c(MOSEK_SOLVER, CLARABEL_SOLVER, SCS_SOLVER,
                              ECOS_SOLVER, GUROBI_SOLVER, CPLEX_SOLVER,
                              XPRESS_SOLVER,
                              SCIP_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",
  SCIP     = "scip",
  XPRESS   = "xpress"
)

## Check if a solver package is installed and not excluded
## Is this solver's R package physically usable?
## (Real install with version >= required minimum, not a stub.)
## PERFORMANCE (2026-08-13). Memoised for the session. Which packages are
## installed cannot change under a running R process in any way this predicate
## should react to, and the MOSEK branch is expensive: `packageVersion()` reads
## and parses Rmosek's DESCRIPTION off disk on every chain construction
## (measured 85us) for users who have Rmosek installed. CVXPY resolves the same
## question exactly once per process -- `INSTALLED_SOLVERS = installed_solvers()`
## is evaluated at import time (defines.py:129) -- so caching here is closer to
## CVXPY's behavior than probing was, not further from it. The exclusion set is
## deliberately NOT part of the key: it is applied by the caller
## `.solver_package_available()`, which stays live.
.solver_package_usable <- function(solver_name) {
  cached <- .cvxr_env$solver_usable[[solver_name]]
  if (!is.null(cached)) return(cached)
  ans <- .solver_package_usable_uncached(solver_name)
  if (is.null(.cvxr_env$solver_usable)) .cvxr_env$solver_usable <- list()
  .cvxr_env$solver_usable[[solver_name]] <- ans
  ans
}

.solver_package_usable_uncached <- function(solver_name) {
  ## Rmosek on CRAN is an ancient stub; only versions >= 10 are usable.
  if (solver_name == "MOSEK") {
    if (!requireNamespace("Rmosek", quietly = TRUE)) return(FALSE)
    return(utils::packageVersion("Rmosek") >= "10")
  }
  pkg <- .SOLVER_PACKAGES[[solver_name]]
  if (is.null(pkg)) return(FALSE)
  requireNamespace(pkg, quietly = TRUE)
}

## Is this solver available for use? (Usable AND not excluded.)
.solver_package_available <- function(solver_name) {
  if (solver_name %in% .cvxr_env$excluded_solvers) return(FALSE)
  .solver_package_usable(solver_name)
}

## Runtime guard: abort if the named solver's package isn't usable.
.require_solver_package <- function(solver_name) {
  if (.solver_package_usable(solver_name)) return(invisible(TRUE))
  pkg <- .SOLVER_PACKAGES[[solver_name]]
  cli_abort("Solver {.code {solver_name}} unavailable: package {.pkg {pkg}} is missing or unusable.")
}

# -- 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()) {
    if (FALSE) new_object(S7_object())  ## S7 static-check guard
    solver_inst <- if (length(reductions) > 0L) {
      reductions[[length(reductions)]]
    } else {
      NULL
    }
    .fast_new(SolvingChain, 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_is(con, Equality) || .s7_is(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)
}

# -- .preserve_variable_bounds() ------------------------------------
## CVXPY SOURCE: solving_chain.py:224-225 wires
## `CvxAttr2Constr(reduce_bounds = not solver_instance.BOUNDED_VARIABLES)`.
##
## CVXPY SOURCE: matrix_stuffing.py:48-52,84-88 rejects sparse bounds that
## reach matrix stuffing, and cone_matrix_stuffing.py:423-431 handles
## parametric bounds through dedicated tensors. CVXR supports that tensor path
## for its 2D scalar-or-exact-shape bound model; sparse bounds still lower to
## constraints.
.preserve_variable_bounds <- function(problem, solver_inst) {
  if (!isTRUE(solver_inst@BOUNDED_VARIABLES)) return(FALSE)
  for (v in variables(problem)) {
    b <- v@attributes$bounds
    if (is.null(b) || !is.list(b)) next
    if (any(vapply(b, function(x) inherits(x, "Matrix"), logical(1L)))) {
      return(FALSE)
    }
  }
  TRUE
}

# -- .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)
}

# -- cone-scan lookup tables ------------------------------------
## PERFORMANCE (2026-08-13). The scan below walks every node of every constraint
## tree plus the objective.  Written as a chain of ~30 sequential `.s7_is()`
## tests per node it cost 21-25% of total solve time and was measured to account
## for 83-100% of the whole CVXR 1.9.2 regression against CRAN 1.9.1
## (notes/session_handoff_2026-08-12_perf_regression.md).  Two facts made it
## expensive: descending into Constraints as well as Expressions multiplied the
## node count 6x (large_lp) to 12x (socp) -- a REQUIRED correctness fix, see the
## comment in the walk below -- and each node paid up to 30 `inherits()` calls.
##
## These tables replace the if-chain with ONE lookup keyed by `class(x)[1L]`,
## which for an S7 object already IS the fully-qualified name ("CVXR::Pnorm").
## They are derived FROM the priority-ordered rule list by walking each class's
## real ancestry, so inheritance semantics are preserved exactly rather than
## re-stated: a subclass lands on the first rule it inherits from, in the same
## order the if-chain tested.  Built lazily on first use because the class
## objects must all exist by then.

.cone_lookup  <- new.env(hash = TRUE, parent = emptyenv())  # fqn -> cone class
.cone_descend <- new.env(hash = TRUE, parent = emptyenv())  # fqn -> TRUE
.cone_maps_built <- new.env(hash = TRUE, parent = emptyenv())

## Does S7 class `cls` have `target` among its ancestors (or equal it)?
.s7_class_inherits <- function(cls, target) {
  while (inherits(cls, "S7_class")) {
    if (identical(cls, target)) return(TRUE)
    cls <- S7::prop(cls, "parent")
  }
  FALSE
}

.build_cone_maps <- function() {
  if (isTRUE(.cone_maps_built$done)) return(invisible(NULL))
  ## PRIORITY ORDER IS LOAD-BEARING and must match the historical if-chain:
  ## Approx subclasses before their exact parents (PowerApprox -> SOC must win
  ## over Power -> PowCone3D).  Classifications follow CVXPY's SOC_ATOMS,
  ## PSD_ATOMS, EXP_ATOMS, POWCONE_ATOMS, POWCONE_ND_ATOMS
  ## (cvxpy/atoms/__init__.py).
  rules <- list(
    list(SOC, SOC), list(PSD, PSD), list(ExpCone, ExpCone),
    list(PowCone3D, PowCone3D), list(PowConeND, PowConeND),
    list(PnormApprox, SOC), list(PowerApprox, SOC), list(GeoMeanApprox, SOC),
    list(QuadOverLin, SOC), list(QuadForm, SOC), list(SymbolicQuadForm, SOC),
    list(Huber, SOC),
    list(Pnorm, PowCone3D), list(Power, PowCone3D),
    list(GeoMean, PowConeND),
    list(SigmaMax, PSD), list(NormNuc, PSD),
    list(MatrixFrac, PSD), list(TrInv, PSD),
    list(LambdaMax, PSD), list(LambdaSumLargest, PSD),
    list(LogDet, PSD), list(ConditionNumber, PSD),
    list(LogSumExp, ExpCone), list(Exp, ExpCone), list(Log, ExpCone),
    list(Entr, ExpCone), list(KlDiv, ExpCone), list(RelEntr, ExpCone),
    list(Logistic, ExpCone), list(Xexp, ExpCone), list(Log1p, ExpCone)
  )
  ns <- asNamespace("CVXR")
  for (nm in ls(ns, all.names = TRUE)) {
    obj <- get0(nm, envir = ns, inherits = FALSE)
    if (!inherits(obj, "S7_class")) next
    fqn <- paste0(attr(obj, "package"), "::", attr(obj, "name"))
    for (r in rules) {
      if (.s7_class_inherits(obj, r[[1L]])) {
        assign(fqn, r[[2L]], envir = .cone_lookup)
        break
      }
    }
    ## Only Expressions and Constraints have `@args` worth descending into.
    ## Precomputing this replaces one-or-two `.s7_is()` calls PER NODE -- for
    ## large_lp that alone was 3k-6k `inherits()` calls per scan.
    if (.s7_class_inherits(obj, Expression) || .s7_class_inherits(obj, Constraint))
      assign(fqn, TRUE, envir = .cone_descend)
  }
  .cone_maps_built$done <- TRUE
  invisible(NULL)
}

# -- .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.
##
## MEMOISED on the problem's own cache: chain construction calls this three
## times for one `psolve()` on an immutable problem.  Cone types depend on the
## expression TREE only, never on parameter values, so a changed parameter does
## not invalidate the answer.

.required_cone_types <- function(problem) {
  cached <- problem@.cache$required_cone_types
  if (!is.null(cached)) return(cached)
  result <- .required_cone_types_uncached(problem)
  problem@.cache$required_cone_types <- result
  result
}

.required_cone_types_uncached <- 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.
  ## The former ~30-test `.s7_is()` chain now lives in `.build_cone_maps()` as a
  ## priority-ordered rule list compiled into `.cone_lookup`; see the block above
  ## for why. Semantics are unchanged -- the map is derived from those same rules
  ## applied over each class's real ancestry.
  .build_cone_maps()
  lookup  <- .cone_lookup
  descend <- .cone_descend

  ## 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
      ## One env lookup on the fully-qualified class name, which for an S7
      ## object IS class(e)[1L] -- replaces up to 30 inherits() calls.
      key <- class(e)[1L]
      cone <- lookup[[key]]
      if (!is.null(cone)) .add_cone(cone)
      ## Descend through Constraints as well as Expressions.  CVXPY scans
      ## `problem.atoms()` (problem_form.py:243-262), which covers atoms inside
      ## CONSTRAINTS, not just the objective.  Restricting the walk to
      ## Expressions meant a constraint's args were never enqueued, so
      ## `lambda_max(X) <= t` reported NO required cones at all.  That was
      ## invisible while it only affected solver preference; once D_19.6 made
      ## the cone set decide whether PSD -> SvecPSD is scheduled, it became a
      ## hard solver failure ("Constraint dimensions inconsistent with cones").
      ##
      ## `descend` is the precomputed Expression-or-Constraint test. Leaves
      ## (Variable, Constant, Parameter) miss it and are never expanded, which
      ## is also what keeps them off the queue in the first place.
      if (!is.null(descend[[key]])) {
        kids <- e@args
        nk <- length(kids)
        if (nk > 0L) {
          nq <- length(queue)
          for (ki in seq_len(nk)) queue[[nq + ki]] <- kids[[ki]]
        }
      }
    }
  }

  ## CVXPY SOURCE: problem_form.py:256-258 -- a PSD/NSD VARIABLE ATTRIBUTE needs
  ## the PSD cone just as a PSD constraint does; `CvxAttr2Constr` turns it into
  ## one later in the chain, long after this scan. Missing it meant both a
  ## non-PSD solver could be selected AND (post-D_19.6) the PSD -> SvecPSD
  ## conversion would not be scheduled, so the packed cone dims disagreed with
  ## the stuffed rows ("cone dimensions 4 not equal to num rows in A = m = 5").
  if (any(vapply(variables(problem), function(v) is_psd(v) || is_nsd(v),
                 logical(1L)))) {
    .add_cone(PSD)
  }

  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
  ## CVXPY SOURCE: solving_chain.py:214-216 -- a cone the solver does not take
  ## directly still counts as supported when an EXACT conversion reaches one it
  ## does (PSD -> SvecPSD for every packed-triangle solver). `expand_cones`
  ## rewrites the requirement set; what survives it must be natively supported.
  expanded <- expand_cones(required_cones, supported)
  missing <- !vapply(expanded$cones, function(cone) {
    any(vapply(supported, identical, logical(1L), cone))
  }, logical(1L), USE.NAMES = FALSE)
  list(ok = !any(missing), unsupported = expanded$cones[missing],
       exact_targets = expanded$exact_targets)
}

# -- .pick_default_conic_solver() ----------------------------------
## CVXPY SOURCE: cvxpy/problems/problem_form.py::pick_default_solver (CVXPY 1.9
## ProblemForm refactor). CVXR keeps its solver resolution inline in
## solving_chain.R (ADR D_19.1), so this implements CVXPY 1.9's *structured*
## default-solver policy for the conic path. It is consulted only when the user
## did not name a solver; it returns a conic solver name, or NULL to fall back
## to the preference-order scan. (CVXR robustness extension: where CVXPY returns
## None and the caller errors, CVXR falls through to any cone-capable solver.)
##   1. commercial first; 2. MI-LP -> HIGHS, MI-other -> SCIP;
##   3. SDP (PSD) -> SCS; 4. LP / SOCP / Exp / Pow -> CLARABEL.
## QP -> OSQP is handled by the QP path (qp_solvers[1]); LP routes to the conic
## path and lands on CLARABEL via rule 4, matching CVXPY.
.COMMERCIAL_CONIC_SOLVERS <- c(MOSEK_SOLVER, GUROBI_SOLVER, CPLEX_SOLVER,
                               XPRESS_SOLVER)

.pick_default_conic_solver <- function(problem, conic_candidates, required_cones) {
  .ok <- function(name) {
    name %in% conic_candidates &&
      .solver_supports_cones(SOLVER_MAP_CONIC[[name]], required_cones)$ok
  }
  ## 1. Commercial solvers first -- unless that is switched off.
  ## CVXPY SOURCE: problem_form.py:356 (`if s.DEFAULT_TO_COMMERCIAL_SOLVERS:`),
  ## settings.py:114-123 (CVXPY 1.9.2, PR #3352). See settings.R for why CVXR
  ## reads the flag per call and also honors an R option.
  if (.default_to_commercial_solvers()) {
    for (name in .COMMERCIAL_CONIC_SOLVERS) if (.ok(name)) return(name)
  }
  ## 2. Mixed-integer: LP -> HIGHS, otherwise -> SCIP.
  if (is_mixed_integer(problem)) {
    if (.is_lp(problem) && .ok(HIGHS_SOLVER)) return(HIGHS_SOLVER)
    if (.ok(SCIP_SOLVER)) return(SCIP_SOLVER)
    return(NULL)
  }
  ## 3. SDP (PSD cone) -> SCS.
  if (any(vapply(required_cones, identical, logical(1L), PSD)) &&
      .ok(SCS_SOLVER)) {
    return(SCS_SOLVER)
  }
  ## 4. LP / SOCP / Exp / Pow -> CLARABEL.
  if (.ok(CLARABEL_SOLVER)) return(CLARABEL_SOLVER)
  NULL
}

# -- .select_conic_solver_name() -----------------------------------
## Single source of truth for which conic solver is chosen: the CVXPY 1.9
## structured default policy when no solver was named, otherwise (and as a
## fallback) the first cone-capable solver in preference order. Returns a name
## or NULL when nothing fits.
.select_conic_solver_name <- function(problem, conic_candidates, required_cones,
                                      solver) {
  if (is.null(solver)) {
    pick <- .pick_default_conic_solver(problem, conic_candidates, required_cones)
    if (!is.null(pick)) return(pick)
  }
  for (s in conic_candidates) {
    if (.solver_supports_cones(SOLVER_MAP_CONIC[[s]], required_cones)$ok) {
      return(s)
    }
  }
  NULL
}

# -- 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,
                                     opts = solver_opts(),
                                     enforce_dpp = FALSE, ignore_dpp = 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]]

  ## CVXPY v1.9.0 #3142: select the solver BEFORE the DPP check so we know
  ## whether the objective's quadratic part will be handled directly by a
  ## quad-objective solver (P maps linearly to data -> param-affine P is
  ## DPP). candidates are reused by the routing section below.
  candidates <- .build_candidates(problem, solver)
  will_solve_as_qp <- .solve_as_qp(problem, candidates) && opts$use_quad_obj
  quad_obj_capable <- if (will_solve_as_qp) {
    TRUE  # QP-path solvers handle quadratic objectives directly
  } else {
    ## Conic path: would the conic solver that gets selected support quad obj?
    ## Use the SAME selection the routing section uses so this probe matches the
    ## solver actually chosen (CVXPY 1.9 structured default policy).
    req_cones <- .required_cone_types(problem)
    sel <- .select_conic_solver_name(problem, candidates$conic_solvers,
                                     req_cones, solver)
    if (is.null(sel)) FALSE else supports_quad_obj(SOLVER_MAP_CONIC[[sel]])
  }
  ## "qp" relaxation applies only to a quadratic objective that a quad-obj
  ## solver will take directly (so parametric P in the OBJECTIVE is DPP;
  ## constraint quad_forms stay non-DPP -- see .problem_is_dpp).
  ## PERF (v1.9): gate on has_params. The relaxation only makes a PARAMETRIC P
  ## objective DPP, so it is meaningless without parameters. Without this gate a
  ## quadratic-objective, parameter-FREE problem (e.g. the Kalman cells) ran a
  ## redundant SECOND .problem_is_dpp full-tree curvature pass (dpp_ok_std +
  ## dpp_ok), a per-node S7-dispatch cost that scales with expression size
  ## (+10-19% on node-dense problems vs 1.8.2-1). See
  ## notes/v19_perf_regression_dpp_rootcause.md. #3142 is preserved: parametric
  ## quad_form(x,P) always has has_params = TRUE.
  quad_form_dpp <- if (!gp && opts$use_quad_obj && quad_obj_capable &&
                       has_quadratic_term(obj_expr) && has_params) {
    "qp"
  } else {
    NULL
  }

  ## DPP context depends on gp flag (CVXPY: dpp_context = 'dgp' if gp else 'dcp')
  ## relaxation_enabled = the quad_form ("qp") relaxation is what makes the
  ## problem DPP (it is NOT DPP under standard rules). Only THEN is it safe to
  ## drop the EvalParams baking and use the P_tensor path -- the parametric
  ## dependence is then genuinely a quad_form(x, P) the tensor path represents.
  ## A problem already DPP by standard rules but with a parametric non-quadratic
  ## atom (e.g. log_det(P)) must keep EvalParams: the QP tensor path cannot
  ## represent it, and a quad-obj solver would otherwise be handed conic data.
  relaxation_enabled <- FALSE
  if (gp) {
    dpp_ok <- .is_dgp_dpp(problem)
  } else if (!is.null(quad_form_dpp)) {
    dpp_ok_std <- .problem_is_dpp(problem, "dcp", NULL)
    dpp_ok     <- .problem_is_dpp(problem, "dcp", quad_form_dpp)  # >= dpp_ok_std
    relaxation_enabled <- dpp_ok && !dpp_ok_std
  } else {
    dpp_ok <- .problem_is_dpp(problem, "dcp", NULL)
  }
  ## CVXPY SOURCE: solving_chain.py:199-204 -- ignore_dpp forces the non-DPP
  ## path; enforce_dpp raises; otherwise a non-DPP parameterized problem WARNS
  ## instead of silently falling back:
  ##
  ##     if ignore_dpp or not is_dpp:
  ##         if not ignore_dpp and enforce_dpp: raise DPPError(DPP_ERROR_MSG)
  ##         if not ignore_dpp:                 warn(DPP_ERROR_MSG)
  ##         reductions = [EvalParams()] + reductions
  ##
  ## Both branches read the TRUE DPP status, and the ignore_dpp line below
  ## OVERWRITES `dpp_ok`, so capture it first.
  dpp_violated <- has_params && !dpp_ok
  if (isTRUE(ignore_dpp)) { dpp_ok <- FALSE; relaxation_enabled <- FALSE }
  ## The abort tests `dpp_violated`, NOT `dpp_ok`: the line above OVERWRITES
  ## `dpp_ok` to FALSE for `ignore_dpp`, so testing it here asks the wrong
  ## question. Upstream's guard is `if not ignore_dpp and enforce_dpp`
  ## (solving_chain.py:200) -- ignore_dpp short-circuits the enforce check --
  ## and `dpp_violated` is `has_params && !dpp_ok` captured before the
  ## overwrite, i.e. upstream's `not is_dpp`.
  ##
  ## No psolve() path reaches here with both flags set: `.compile()`
  ## (problems/problem.R:489) rejects the contradictory pair first, as upstream
  ## does at problem.py:817-818, and for the same stated reason (before the
  ## cache key). This matters only for a direct call to the chain constructor,
  ## where CVXR now behaves as CVXPY's does.
  if (isTRUE(enforce_dpp) && !isTRUE(ignore_dpp) && dpp_violated) {
    cli_abort(c(
      "Problem does not follow DPP rules but {.code enforce_dpp = TRUE}.",
      "i" = "Set {.code enforce_dpp = FALSE} to canonicalize it as a non-DPP problem."
    ), class = "DPPError")   ## CVXPY: raise DPPError, solving_chain.py:201
  }
  ## CVXPY v1.9.2 parity (audit finding #5): warn on the ordinary path. Without
  ## this a user whose re-solves are NOT being accelerated is never told -- the
  ## problem solves correctly, just without the DPP fast path, which is exactly
  ## the silent-underperformance case upstream's message exists to prevent.
  ## `ignore_dpp = TRUE` means the user already knows, so it stays silent.
  ## The abort above returns first, so at most one of the two ever fires.
  if (dpp_violated && !isTRUE(ignore_dpp)) {
    cli_warn(c(
      "Solving a parameterized problem that is not DPP.",
      "i" = "Subsequent solves will not be faster than the first one.",
      "i" = "Check with {.fn is_dpp}; see {.url https://www.cvxpy.org/tutorial/dpp/index.html}.",
      "i" = "Set {.code ignore_dpp = TRUE} to silence this."
    ), class = "DPPWarning")   ## CVXPY: warn(DPP_ERROR_MSG), solving_chain.py:203
  }
  qp_with_params <- has_params && !gp &&
    is_quadratic(obj_expr) && !is_affine(obj_expr)
  ## CVXPY v1.9.0 #3142: when the qp relaxation is what makes the parametric
  ## quadratic objective DPP (genuine quad_form(x, P)), do NOT bake P via
  ## EvalParams -- let it flow through cone_matrix_stuffing as a P_tensor so
  ## re-solves rebuild P from new parameter values (no staleness).
  if (relaxation_enabled) {
    qp_with_params <- FALSE
  }
  ## Deliberate feature gap: complex parameters use EvalParams before
  ## Complex2Real. Ordinary complex solves are supported, but the DPP fast path
  ## would need complex sparse parameter tensors; R's Matrix package does not
  ## provide the required complex sparse representation here. The
  ## complex-dpp-parity tests document the deferred CVXPY fast-path behavior.
  has_complex_params <- has_params &&
    any(vapply(parameters(problem), function(p) is_complex(p) || is_imag(p),
               logical(1L)))
  ## When the user requested requires_grad (forces solver = DIFFCP), the
  ## parameters MUST survive into the conic representation so backward
  ## and derivative can chain-rule through them; skip the QP-with-params
  ## EvalParams optimisation in that case.  The solver being DIFFCP
  ## already implies DPP-compliance was checked in psolve().
  is_diffcp_path <- !is.null(solver) && identical(toupper(solver), DIFFCP_SOLVER)
  if (has_params &&
      (!dpp_ok || (qp_with_params && !is_diffcp_path) || 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."
      ), class = "DGPError")   ## CVXPY: raise DGPError, intermediate_chain.py:81
    }
    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_is(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_is(c, FiniteSet), logical(1)))) {
    reductions <- c(reductions, list(Valinvec2mixedint()))
  }

  ## 3b. (candidates were built above, before the DPP check -- #3142)

  ## 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
  ## CVXPY v1.8.2: respect opts$use_quad_obj — when FALSE, skip QP path
  ## entirely so quadratic objectives go through conic decomposition
  ## (where quad_form_canon can catch indefinite P).
  ## PERF (v1.9): reuse will_solve_as_qp computed above (line ~378). #3142 moved
  ## the .solve_as_qp/.build_candidates evaluation before the DPP check; this
  ## routing test is the SAME expression, so recomputing it called
  ## .solve_as_qp -> .is_lp a second time per solve (an extra uncached
  ## variables(problem) walk + cached is_dcp/is_pwl traversal). Reuse the value.
  if (will_solve_as_qp) {
    ## 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]])
    ## CVXPY SOURCE: solving_chain.py:345 -- carries the chosen solver's
    ## supported-constraint set into Dcp2Cone so SOC-approx canonicalizers
    ## can emit a warning when the solver actually supports power cones.
    solver_context <- SolverInfo(
      solver_name = solver_name_sel,
      solver_supported_constraints = solver_inst@SUPPORTED_CONSTRAINTS,
      solver_supports_bounds = isTRUE(solver_inst@BOUNDED_VARIABLES),
      ## CVXPY SOURCE: solving_chain.py:173-174
      psd_triangle_kind = solver_inst@PSD_TRIANGLE_KIND,
      psd_sqrt2_scaling = solver_inst@PSD_SQRT2_SCALING
    )
    reductions <- c(reductions, list(
      Dcp2Cone(quad_obj = quad_obj, solver_context = solver_context),
      CvxAttr2Constr(reduce_bounds = !.preserve_variable_bounds(problem, solver_inst))
    ))
    ## CVXPY SOURCE: solving_chain.py:227-229 -- the exact cone conversions the
    ## chosen solver needs, computed by `expand_cones` from the problem's cones
    ## and the solver's SUPPORTED_CONSTRAINTS.
    exact_targets <- expand_cones(.required_cone_types(problem),
                                  solver_inst@SUPPORTED_CONSTRAINTS)$exact_targets
    if (length(exact_targets) > 0L) {
      reductions <- c(reductions, list(
        ExactCone2Cone(target_cones = exact_targets,
                       solver_context = solver_context)))
    }
    reductions <- c(reductions, list(
      ConeMatrixStuffing(quad_obj = quad_obj),
      solver_inst
    ))
  } else {
    ## Conic path -- choose a solver that accepts the problem's cones.
    ## CVXPY 1.9 structured default policy when no solver was named
    ## (.pick_default_conic_solver), else first cone-capable in preference order.
    required_cones <- .required_cone_types(problem)
    solver_name_sel <- .select_conic_solver_name(problem, candidates$conic_solvers,
                                                 required_cones, solver)
    solver_inst <- if (is.null(solver_name_sel)) NULL else
      SOLVER_MAP_CONIC[[solver_name_sel]]

    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.")
      }
    }

    ## CVXPY v1.8.2: quad_obj in conic path depends on use_quad_obj,
    ## solver capability, and whether objective has a quadratic term.
    quad_obj <- opts$use_quad_obj &&
                supports_quad_obj(solver_inst) &&
                has_quadratic_term(problem@objective@args[[1L]])
    ## CVXPY SOURCE: solving_chain.py:345
    solver_context <- SolverInfo(
      solver_name = solver_name_sel,
      solver_supported_constraints = solver_inst@SUPPORTED_CONSTRAINTS,
      solver_supports_bounds = isTRUE(solver_inst@BOUNDED_VARIABLES),
      ## CVXPY SOURCE: solving_chain.py:173-174
      psd_triangle_kind = solver_inst@PSD_TRIANGLE_KIND,
      psd_sqrt2_scaling = solver_inst@PSD_SQRT2_SCALING
    )
    reductions <- c(reductions, list(
      Dcp2Cone(quad_obj = quad_obj, solver_context = solver_context),
      CvxAttr2Constr(reduce_bounds = !.preserve_variable_bounds(problem, solver_inst))
    ))
    ## CVXPY SOURCE: solving_chain.py:227-229 -- the exact cone conversions the
    ## chosen solver needs, computed by `expand_cones` from the problem's cones
    ## and the solver's SUPPORTED_CONSTRAINTS.
    exact_targets <- expand_cones(.required_cone_types(problem),
                                  solver_inst@SUPPORTED_CONSTRAINTS)$exact_targets
    if (length(exact_targets) > 0L) {
      reductions <- c(reductions, list(
        ExactCone2Cone(target_cones = exact_targets,
                       solver_context = solver_context)))
    }
    reductions <- c(reductions, list(
      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).

## S7 method (registered via S7::methods_register() in .onLoad; not exported and
## not separately documented -- it carries an extra `problem` arg beyond the
## generic, which roxygen 8 would emit as a bare `solve_via_data(..., problem=)`
## usage that R's codoc cannot distinguish from the generic. The generic's
## `@param ...` documents `problem`. See the generic in generics.R.
method(solve_via_data, SolvingChain) <- function(x, data, warm_start = FALSE,
                                                  verbose = FALSE,
                                                  solver_opts = list(), ...,
                                                  problem = NULL) {
  ## CVXPY SOURCE: solving_chain.py:556-560 -- upstream calls
  ## `_validate_problem_data(data)` HERE as well as in `SolvingChain.solve`
  ## (:508-512), with the comment naming why: "These are the two possible entry
  ## points for executing the solving chain." CVXR had only the psolve() site
  ## (problem.R), so the documented decomposed API --
  ## problem_data() -> solve_via_data() -> problem_unpack_results() -- accepted
  ## poisoned data: a NaN planted in `c` reached Clarabel, which returned status
  ## 10 and garbage that propagated to the caller. CVXPY raises ValueError.
  .validate_problem_data(data)

  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)
}