R/141_problems_problem.R

In CVXR: Disciplined Convex Optimization

#####
## DO NOT EDIT THIS FILE!! EDIT THE SOURCE INSTEAD: rsrc_tree/problems/problem.R
#####

## CVXPY SOURCE: problems/problem.py
## Problem -- optimization problem with objective and constraints
##
## CVXPY 1.9 parity notes:
##  - `is_dnlp` mirrors problem.py:295-299: objective and all constraints must
##    satisfy the DNLP predicate.
##  - CVXR exposes NLP solving through `psolve(..., nlp = TRUE)` and the R NLP
##    solver chain rather than CVXPY's Python `solve_nlp` helper surface.

# -- Problem class -------------------------------------------------
## CVXPY SOURCE: problem.py lines 156-193

#' Create an Optimization Problem
#'
#' Constructs a convex optimization problem from an objective and a list of
#' constraints. Use \code{\link{psolve}} to solve the problem.
#'
#' @param objective A \code{\link{Minimize}} or \code{\link{Maximize}} object.
#' @param constraints A list of \code{Constraint} objects (e.g., created by
#'   \code{==}, \code{<=}, \code{>=} operators on expressions). Defaults to
#'   an empty list (unconstrained).
#' @returns A \code{Problem} object.
#'
#' @section Known limitations:
#' \itemize{
#'   \item Problems must contain at least one \code{\link{Variable}}.
#'     Zero-variable problems (e.g., minimizing a constant) will cause an
#'     internal error in the reduction pipeline.
#' }
#'
#' @examples
#' x <- Variable(2)
#' prob <- Problem(Minimize(sum_entries(x)), list(x >= 1))
#'
#' @export
Problem <- new_class("Problem", package = "CVXR",
  properties = list(
    objective   = class_any,
    constraints = class_list,
    .cache      = class_environment
  ),
  constructor = function(objective, constraints = list()) {
    if (FALSE) new_object(S7_object())  ## S7 static-check guard
    ## Validate objective type
    ## CVXPY SOURCE: problem.py lines 163-164
    if (!.s7_is(objective, Objective)) {
      cli_abort("Problem objective must be {.cls Minimize} or {.cls Maximize}.")
    }
    ## Validate constraints list
    ## CVXPY SOURCE: problem.py lines 127-136 (_validate_constraint)
    for (i in seq_along(constraints)) {
      ci <- constraints[[i]]
      if (isTRUE(ci)) {
        ## TRUE -> trivially satisfied: 0 <= 1
        constraints[[i]] <- Inequality(Constant(0), Constant(1))
      } else if (identical(ci, FALSE)) {
        ## FALSE -> infeasible: 1 <= 0
        constraints[[i]] <- Inequality(Constant(1), Constant(0))
      } else if (!.s7_is(ci, Constraint)) {
        cli_abort("Element {i} of constraints is not a {.cls Constraint} object.")
      }
    }
    .fast_new(Problem, S7_object(),
      objective   = objective,
      constraints = constraints,
      .cache      = new.env(parent = emptyenv())
    )
  }
)

## Lazy-init solver cache on Problem's .cache environment.
## Uses an env (reference semantics) so solver writes propagate back.
## Matches CVXPY's problem._solver_cache dict.
.get_solver_cache <- function(problem) {
  if (is.null(problem@.cache$solver_cache)) {
    problem@.cache$solver_cache <- new.env(hash = TRUE, parent = emptyenv())
  }
  problem@.cache$solver_cache
}

# -- is_dcp --------------------------------------------------------
## CVXPY SOURCE: problem.py lines 273-294
## DCP if objective and all constraints are DCP.

method(is_dcp, Problem) <- function(x) {
  key <- .dpp_key("is_dcp")
  cached <- cache_get(x, key)
  if (!cache_miss(cached)) return(cached)
  result <- is_dcp(x@objective) &&
    all(vapply(x@constraints, is_dcp, logical(1)))
  cache_set(x, key, result)
  result
}

## is_dqcp: DQCP if objective and all constraints are DQCP
## CVXPY SOURCE: problem.py lines 334-338
method(is_dqcp, Problem) <- function(x) {
  key <- "is_dqcp"
  cached <- cache_get(x, key)
  if (!cache_miss(cached)) return(cached)
  result <- is_dqcp(x@objective) &&
    all(vapply(x@constraints, is_dqcp, logical(1)))
  cache_set(x, key, result)
  result
}

## is_dnlp: DNLP if objective and all constraints are DNLP
## CVXPY SOURCE: problem.py lines 294-299 (@perf.compute_once)
method(is_dnlp, Problem) <- function(x) {
  key <- "is_dnlp"
  cached <- cache_get(x, key)
  if (!cache_miss(cached)) return(cached)
  result <- is_dnlp(x@objective) &&
    all(vapply(x@constraints, is_dnlp, logical(1)))
  cache_set(x, key, result)
  result
}

## is_dgp: DGP if objective and all constraints are DGP
## CVXPY SOURCE: problem.py lines 310-331
method(is_dgp, Problem) <- function(x) {
  key <- "is_dgp"
  cached <- cache_get(x, key)
  if (!cache_miss(cached)) return(cached)
  result <- is_dgp(x@objective) &&
    all(vapply(x@constraints, is_dgp, logical(1)))
  cache_set(x, key, result)
  result
}

## is_dpp: DPP compliance (objective + all constraints DPP)
## CVXPY SOURCE: problem.py lines 341-369
## Extended: also checks atom-level is_dpp on all sub-expressions.
## This catches atoms like Kron/Conv whose C++ handlers can't process
## PARAM LinOp nodes even though the expression is structurally DCP.
method(is_dpp, Problem) <- function(x, context = "dcp") {
  .problem_is_dpp(x, context, quad_form_dpp = NULL)
}

## CVXPY v1.9.0 #3142: Problem.is_dpp(context, quad_form_dpp).
## quad_form_dpp = "qp" enters quad_form_dpp_scope for the OBJECTIVE ONLY, so a
## parametric-P quad_form in the objective is DPP (a quad-obj solver maps P
## linearly to data), while a parametric-P quad_form in a CONSTRAINT stays
## non-DPP -- the conic canonicalizer would bake in numeric Cholesky factors,
## which a re-solve with a changed P would silently reuse (stale). The solving
## chain passes "qp" iff the chosen solver supports_quad_obj().
## The generic is_dpp signature is unchanged; this parameterized helper is
## called directly by the chain.
.problem_is_dpp <- function(x, context = "dcp", quad_form_dpp = NULL) {
  dgp <- identical(tolower(context), "dgp")
  qp_relax <- !dgp && identical(quad_form_dpp, "qp")
  ## Run a thunk inside quad_form_dpp_scope only on the objective-checking path.
  in_qf_scope <- function(thunk) {
    if (qp_relax) with_quad_form_dpp_scope(thunk()) else thunk()
  }

  ## First: structural DCP/DGP compliance in DPP scope. The objective is
  ## additionally checked in quad_form_dpp_scope when qp_relax; constraints
  ## are NOT (parametric P in a constraint must stay non-DPP).
  base_ok <- with_dpp_scope({
    obj_ok <- in_qf_scope(function() {
      if (dgp) is_dgp(x@objective) else is_dcp(x@objective)
    })
    if (!obj_ok) {
      FALSE
    } else if (dgp) {
      all(vapply(x@constraints, is_dgp, logical(1)))
    } else {
      all(vapply(x@constraints, is_dcp, logical(1)))
    }
  })
  if (!base_ok) return(FALSE)
  ## Second: all atoms in the tree must be DPP-compatible (catches non-affine /
  ## non-log-log-affine variable bounds, parametric atom restrictions, etc.).
  ## The objective's atom-level check also runs in quad_form_dpp_scope.
  obj_atom_ok <- in_qf_scope(function() is_dpp(x@objective@args[[1L]], context))
  if (!obj_atom_ok) return(FALSE)
  all(vapply(x@constraints,
             function(c) .all_args(c, function(a) is_dpp(a, context)), logical(1)))
}

## is_dgp_dpp: DGP compliance in DPP scope
## CVXPY SOURCE: problem.py is_dpp(context='dgp') -> is_dgp(dpp=True)
## Checks log-log convexity/concavity with Parameters treated as positive/affine.
.is_dgp_dpp <- function(problem) {
  with_dpp_scope({
    is_dgp(problem@objective) &&
      all(vapply(problem@constraints, is_dgp, logical(1)))
  })
}

# -- is_qp ---------------------------------------------------------
## CVXPY SOURCE: problem.py lines 371-393
## QP if DCP, all inequality constraints are PWL, no conic constraints,
## and objective is QPWA (quadratic or piecewise affine).

method(is_qp, Problem) <- function(x) {
  if (!is_dcp(x)) return(FALSE)
  for (con in x@constraints) {
    if (.s7_is(con, Inequality) || .s7_is(con, NonPos) ||
        .s7_is(con, NonNeg)) {
      ## Get the constraint expression
      con_expr <- if (.s7_is(con, Inequality)) con@.expr else con@args[[1L]]
      if (!is_pwl(con_expr)) return(FALSE)
    } else if (!.s7_is(con, Equality) && !.s7_is(con, Zero)) {
      return(FALSE)
    }
  }
  ## CVXPY SOURCE: problem.py lines 390-392
  ## Reject PSD/NSD/hermitian variables (these require SDP, not QP)
  for (v in variables(x)) {
    a <- v@attributes
    if (isTRUE(a$PSD) || isTRUE(a$NSD) || isTRUE(a$hermitian)) return(FALSE)
  }
  is_qpwa(x@objective@args[[1L]])
}

# -- is_lp ---------------------------------------------------------
## CVXPY SOURCE: problem.py lines 395-422
## LP if QP and objective is also PWL (linear, not quadratic).

method(is_lp, Problem) <- function(x) {
  is_qp(x) && is_pwl(x@objective@args[[1L]])
}

# -- is_mixed_integer ----------------------------------------------
## CVXPY SOURCE: problem.py lines 473-488
## MIP if any variable has boolean or integer attribute.

#' Check if a Problem is Mixed-Integer
#'
#' Returns \code{TRUE} if any variable in the problem has a
#' \code{boolean} or \code{integer} attribute.
#'
#' @param problem A \code{\link{Problem}} object.
#' @returns Logical scalar.
#' @export
is_mixed_integer <- function(problem) {
  cached <- cache_get(problem, "is_mixed_integer")
  if (!cache_miss(cached)) return(cached)
  ## CVXPY SOURCE: problem.py:451-453 -- `any(v.attributes['boolean'] or
  ## v.attributes['integer'] ...)`, i.e. TRUTHINESS: a non-empty index list
  ## counts. `isTRUE()` is FALSE for a vector, so a PARTIAL boolean/integer
  ## attribute used to leave the problem looking continuous, and it was then
  ## solved as its relaxation and reported `optimal`.
  result <- any(vapply(variables(problem), function(v) {
    length(v@.boolean_idx) > 0L || length(v@.integer_idx) > 0L
  }, logical(1L)))
  cache_set(problem, "is_mixed_integer", result)
  result
}

# -- variables -----------------------------------------------------
## CVXPY SOURCE: problem.py lines 425-436

method(variables, Problem) <- function(x) {
  cached <- cache_get(x, "variables")
  if (!cache_miss(cached)) return(cached)
  nc <- length(x@constraints)
  parts <- vector("list", nc + 1L)
  parts[[1L]] <- variables(x@objective)
  for (i in seq_len(nc)) {
    parts[[i + 1L]] <- variables(x@constraints[[i]])
  }
  result <- unique_list(unlist(parts, recursive = FALSE))
  cache_set(x, "variables", result)
  result
}

# -- parameters ----------------------------------------------------
## CVXPY SOURCE: problem.py lines 438-450

method(parameters, Problem) <- function(x) {
  cached <- cache_get(x, "parameters")
  if (!cache_miss(cached)) return(cached)
  nc <- length(x@constraints)
  parts <- vector("list", nc + 1L)
  parts[[1L]] <- parameters(x@objective)
  for (i in seq_len(nc)) {
    parts[[i + 1L]] <- parameters(x@constraints[[i]])
  }
  result <- unique_list(unlist(parts, recursive = FALSE))
  cache_set(x, "parameters", result)
  result
}

# -- constants -----------------------------------------------------
## CVXPY SOURCE: problem.py lines 452-468

method(constants, Problem) <- function(x) {
  cached <- cache_get(x, "constants")
  if (!cache_miss(cached)) return(cached)
  consts <- constants(x@objective)
  for (con in x@constraints) {
    consts <- c(consts, constants(con))
  }
  result <- unique_list(consts)
  cache_set(x, "constants", result)
  result
}

# -- param_dict / var_dict / size_metrics --------------------------
## CVXPY SOURCE: problem.py lines 260-264 (param_dict),
##               problem.py lines 267-271 (var_dict),
##               problem.py lines 486-490 (size_metrics factory),
##               problem.py lines 1690-1752 (SizeMetrics class).

method(param_dict, Problem) <- function(x) {
  ps <- parameters(x)
  if (length(ps) == 0L) return(list())
  nms <- vapply(ps, expr_name, character(1L))
  setNames(ps, nms)
}

method(var_dict, Problem) <- function(x) {
  vs <- variables(x)
  if (length(vs) == 0L) return(list())
  nms <- vapply(vs, expr_name, character(1L))
  setNames(vs, nms)
}

# -- SizeMetrics class ---------------------------------------------
## CVXPY SOURCE: problem.py lines 1690-1752

#' Problem Size Metrics
#'
#' Reports scalar-counts and data-dimension metrics for a
#' \code{\link{Problem}}.  Constructed by \code{\link{size_metrics}};
#' end users normally call \code{size_metrics(prob)} rather than this
#' constructor directly.
#'
#' @param num_scalar_variables Total scalar entries across all
#'   variables in the problem.
#' @param num_scalar_data Total scalar entries across all constants
#'   and parameters.
#' @param num_scalar_eq_constr Total scalar entries in equality
#'   (\code{Equality}, \code{Zero}) constraints.
#' @param num_scalar_leq_constr Total scalar entries in inequality
#'   (\code{Inequality}, \code{NonNeg}, \code{NonPos}) constraints.
#' @param max_data_dimension Largest single dimension across any data
#'   block (constant or parameter).
#' @param max_big_small_squared Maximum of \code{big * small^2} over
#'   all data blocks, where \code{big}/\code{small} are the larger/
#'   smaller dimension of the block.
#' @returns A \code{SizeMetrics} object.
#' @export
SizeMetrics <- new_class("SizeMetrics", package = "CVXR",
  properties = list(
    num_scalar_variables   = new_property(class = class_integer, default = 0L),
    num_scalar_data        = new_property(class = class_integer, default = 0L),
    num_scalar_eq_constr   = new_property(class = class_integer, default = 0L),
    num_scalar_leq_constr  = new_property(class = class_integer, default = 0L),
    max_data_dimension     = new_property(class = class_integer, default = 0L),
    max_big_small_squared  = new_property(class = class_double,  default = 0)
  )
)

method(size_metrics, Problem) <- function(x) {
  cached <- cache_get(x, "size_metrics")
  if (!cache_miss(cached)) return(cached)

  ## num_scalar_variables
  vars <- variables(x)
  n_vars <- 0L
  for (v in vars) n_vars <- n_vars + as.integer(prod(v@shape))

  ## num_scalar_data, max_data_dimension, max_big_small_squared
  data_blocks <- c(constants(x), parameters(x))
  n_data <- 0L
  max_dim <- 0L
  max_bss <- 0
  for (d in data_blocks) {
    sh <- d@shape
    n_data <- n_data + as.integer(prod(sh))
    big   <- max(sh)
    small <- min(sh)
    if (big > max_dim) max_dim <- as.integer(big)
    bss <- as.numeric(big) * (as.numeric(small)^2)
    if (bss > max_bss) max_bss <- bss
  }

  ## num_scalar_eq_constr / num_scalar_leq_constr
  n_eq <- 0L; n_leq <- 0L
  for (con in x@constraints) {
    sz <- as.integer(prod(con@shape))
    if (.s7_is(con, Equality) || .s7_is(con, Zero)) {
      n_eq <- n_eq + sz
    } else if (.s7_is(con, Inequality) ||
               .s7_is(con, NonNeg) ||
               .s7_is(con, NonPos)) {
      n_leq <- n_leq + sz
    }
  }

  result <- SizeMetrics(
    num_scalar_variables  = n_vars,
    num_scalar_data       = n_data,
    num_scalar_eq_constr  = n_eq,
    num_scalar_leq_constr = n_leq,
    max_data_dimension    = max_dim,
    max_big_small_squared = max_bss
  )
  cache_set(x, "size_metrics", result)
  result
}

# -- value ---------------------------------------------------------
## CVXPY SOURCE: problem.py lines 216-230
## Returns the objective value from last solve (or NULL).

method(value, Problem) <- function(x) {
  v <- x@.cache$value
  if (is.null(v)) return(NULL)
  scalar_value(v)
}

# -- status accessor ----------------------------------------------

#' Get the Solution Status of a Problem
#'
#' Returns the status string from the most recent solve, such as
#' \code{"optimal"}, \code{"infeasible"}, or \code{"unbounded"}.
#'
#' @param x A \code{\link{Problem}} object.
#' @returns Character string, or \code{NULL} if the problem has not been
#'   solved.
#' @seealso \code{\link{OPTIMAL}}, \code{\link{INFEASIBLE}},
#'   \code{\link{UNBOUNDED}}
#' @export
status <- function(x) {
  if (!.s7_is(x, Problem)) {
    cli_abort("{.fn status} requires a {.cls Problem} object.")
  }
  x@.cache$status
}

#' Get the Solution Status of a Problem (deprecated)
#'
#' `r lifecycle::badge("deprecated")`
#'
#' Use \code{\link{status}} instead.
#'
#' @param x A \code{\link{Problem}} object.
#' @returns Character string, or \code{NULL} if the problem has not been
#'   solved.
#' @seealso \code{\link{status}}
#' @export
problem_status <- function(x) {
  cli_warn("{.fn problem_status} is deprecated. Use {.fn status} instead.",
           .frequency = "once", .frequency_id = "cvxr_problem_status_deprecated")
  status(x)
}

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

method(print, Problem) <- function(x, ...) {
  cat(sprintf("Problem(%s, %d constraints)\n",
    x@objective@.name, length(x@constraints)))
  invisible(x)
}

# -- Compilation caching ------------------------------------------
## CVXPY SOURCE: problem.py Cache class (lines 107-125)
## CVXPY v1.8.2: cache key includes use_quad_obj (from solver_opts).

.compile <- function(problem, solver = NULL, gp = FALSE, opts = solver_opts(),
                     enforce_dpp = FALSE, ignore_dpp = FALSE) {
  ## CVXPY SOURCE: problem.py lines 816-818 --
  ##     # Invalid DPP setting.
  ##     # Must be checked here to avoid cache issues.
  ##     if enforce_dpp and ignore_dpp:
  ##         raise DPPError("Cannot set enforce_dpp = True and ignore_dpp = True.")
  ##
  ## The two flags are contradictory: enforce_dpp says "refuse to solve unless
  ## this is DPP", ignore_dpp says "pretend it is not DPP". Upstream rejects the
  ## pair UNCONDITIONALLY -- verified against CVXPY 1.9.2 for a DPP, a non-DPP
  ## and a parameter-free problem.
  ##
  ## Without this check CVXR gave two wrong answers. `ignore_dpp` forces
  ## `dpp_ok <- FALSE` in construct_solving_chain(), so `enforce_dpp` then
  ## aborted with "Problem does not follow DPP rules" -- a FALSE STATEMENT when
  ## the problem is in fact DPP -- and for a PARAMETER-FREE problem the abort is
  ## gated on has_params, so the contradictory pair was silently ACCEPTED and the
  ## problem solved. Wrong message in two cases, wrong outcome in the third.
  ##
  ## Position matters and is upstream's own reason: BEFORE the cache key, so a
  ## rejected call cannot interact with the cache at all (cf. the ordering bug
  ## fixed just below).
  if (isTRUE(enforce_dpp) && isTRUE(ignore_dpp)) {
    cli_abort(c(
      "Cannot set {.code enforce_dpp = TRUE} and {.code ignore_dpp = TRUE}.",
      "i" = "{.code enforce_dpp} requires the problem to be DPP; {.code ignore_dpp} treats it as though it were not."
    ), class = "DPPError")   ## CVXPY: raise DPPError, problem.py:818
  }
  ## CVXPY SOURCE: problem.py lines 786-806, and the cache block at 831-841:
  ##     if key != self._cache.key:
  ##         self._cache.invalidate()          # clear FIRST
  ##         solving_chain = self._construct_chain(...)   # may raise
  ##         self._cache.key = key             # record only on success
  ##         self._cache.solving_chain = solving_chain
  ##
  ## The two orderings matter and CVXR had both wrong. It recorded the key
  ## BEFORE building, and cleared nothing, so an abort from
  ## construct_solving_chain() -- non-DCP, non-DGP under gp = TRUE, an
  ## unavailable solver, enforce_dpp -- left the new key paired with the
  ## PREVIOUS chain. The next call with the same arguments then matched the
  ## cache and returned that stale chain:
  ##
  ##     prob <- Problem(Minimize(sum_entries(x)), list(x >= 1))  # DCP, not DGP
  ##     psolve(prob, solver = "CLARABEL")   # 2
  ##     psolve(prob, gp = TRUE)             # aborts, "not DGP compliant"
  ##     psolve(prob, gp = TRUE)             # 2  <- SILENTLY solved as DCP
  ##
  ## With no prior successful compile the same path yields a NULL chain and
  ## "no applicable method for `@` applied to an object of class NULL".
  ## Found while writing tests for the completeness audit; see
  ## notes/session_handoff_2026-08-17_completeness_audit_and_fixes.md section 5.1.
  cache_key <- list(solver, gp, opts$use_quad_obj, enforce_dpp, ignore_dpp)
  if (!identical(cache_key, problem@.cache$compile_key)) {
    ## Invalidate first: if the build aborts, the cache must be EMPTY, not
    ## stale. Anything that survives here would be paired with the wrong key.
    problem@.cache$compile_key <- NULL
    problem@.cache$compile_chain <- NULL
    problem@.cache$param_prog <- NULL
    problem@.cache$compile_inverse_data <- NULL
    ## Clear solver-specific cache when chain changes
    problem@.cache$solver_cache <- NULL
    ## Build into a local; record the key only once it has succeeded.
    chain <- construct_solving_chain(problem, solver,
                                     gp = gp, opts = opts,
                                     enforce_dpp = enforce_dpp,
                                     ignore_dpp = ignore_dpp)
    problem@.cache$compile_key <- cache_key
    problem@.cache$compile_chain <- chain
  }
  problem@.cache$compile_chain
}

# -- problem_data --------------------------------------------------
## CVXPY SOURCE: problem.py get_problem_data() (simplified)
## Applies the solving chain and returns solver-ready data.
##
## Returns list(data, chain, inverse_data) where:
##   data: solver-ready data (A, b, c, cone_dims)
##   chain: SolvingChain used
##   inverse_data: list of inverse data from each reduction

method(problem_data, Problem) <- function(x, solver = NULL, gp = FALSE,
                                          enforce_dpp = FALSE, ignore_dpp = FALSE, ...) {
  opts <- solver_opts(...)
  chain <- .compile(x, solver, gp = gp, opts = opts,
                    enforce_dpp = enforce_dpp, ignore_dpp = ignore_dpp)
  result <- reduction_apply(chain, x)
  list(data = result[[1L]], chain = chain, inverse_data = result[[2L]])
}

## Deprecated wrapper -- delegates to problem_data()
method(get_problem_data, Problem) <- function(x, solver = NULL, gp = FALSE,
                                              enforce_dpp = FALSE, ignore_dpp = FALSE, ...) {
  cli_warn("{.fn get_problem_data} is deprecated. Use {.fn problem_data} instead.",
           .frequency = "once", .frequency_id = "cvxr_get_problem_data_deprecated")
  problem_data(x, solver = solver, gp = gp,
               enforce_dpp = enforce_dpp, ignore_dpp = ignore_dpp, ...)
}

# ==================================================================
# SolverStats -- miscellaneous solver output
# ==================================================================
## CVXPY SOURCE: problem.py lines 1637-1687

SolverStats <- new_class("SolverStats", package = "CVXR",
  properties = list(
    solver_name = class_character,
    solve_time  = class_any,     # numeric or NULL
    setup_time  = class_any,     # numeric or NULL
    num_iters   = class_any,     # integer or NULL
    extra_stats = class_any      # list or NULL
  ),
  constructor = function(solver_name, solve_time = NULL, setup_time = NULL,
                         num_iters = NULL, extra_stats = NULL) {
    if (FALSE) new_object(S7_object())  ## S7 static-check guard
    .fast_new(SolverStats, S7_object(),
      solver_name = solver_name,
      solve_time  = solve_time,
      setup_time  = setup_time,
      num_iters   = num_iters,
      extra_stats = extra_stats
    )
  }
)

method(print, SolverStats) <- function(x, ...) {
  cat(sprintf("SolverStats(solver=%s, time=%.4f, iters=%s)\n",
    x@solver_name,
    if (is.null(x@solve_time)) NA_real_ else x@solve_time,
    if (is.null(x@num_iters)) "NA" else as.character(x@num_iters)))
  invisible(x)
}

## Factory: construct from Solution attr dict
## CVXPY SOURCE: problem.py SolverStats.from_dict()
solver_stats_from_dict <- function(attr, solver_name) {
  SolverStats(
    solver_name = solver_name,
    solve_time  = attr[[RK_SOLVE_TIME]],
    setup_time  = attr[[RK_SETUP_TIME]],
    num_iters   = attr[[RK_NUM_ITERS]],
    extra_stats = attr[[RK_EXTRA_STATS]]
  )
}

# ==================================================================
# problem_unpack -- apply solution to variables/constraints
# ==================================================================
## CVXPY SOURCE: problem.py lines 1482-1518

problem_unpack <- function(problem, solution) {
  ## `solution@dual_vars[[as.character(con@id)]]` is a LINEAR SCAN of the names
  ## of an n-entry list, so looking one up per constraint was O(n^2). Resolve
  ## all n at once with a single vectorized `match()` (one C call over a hash
  ## table) and index positionally. Same lookups, same order, same values.
  ## See notes/string_key_hashing_sweep_2026-08-13.md and ADR D_PERF.7.
  dual_idx <- match(vapply(problem@constraints,
                           function(con) as.character(con@id), character(1)),
                    names(solution@dual_vars))
  dual_at <- function(i) {
    if (is.na(dual_idx[i])) NULL else solution@dual_vars[[dual_idx[i]]]
  }
  if (solution@status %in% SOLUTION_PRESENT) {
    for (v in variables(problem)) {
      save_leaf_value(v, solution@primal_vars[[as.character(v@id)]])
    }
    for (i in seq_along(problem@constraints)) {
      con <- problem@constraints[[i]]
      dual_val <- dual_at(i)
      if (!is.null(dual_val)) {
        save_dual_value(con, dual_val)
      }
    }
    ## Store objective value
    problem@.cache$value <- value(problem@objective)
  } else if (solution@status %in% INF_OR_UNB) {
    for (v in variables(problem)) {
      save_leaf_value(v, NULL)
    }
    ## CVXPY v1.9.0 fix: #3197 -- if the solver returned an infeasibility
    ## certificate for a constraint, expose it as that constraint's dual value;
    ## otherwise clear the dual variables. Mirrors problem.py unpack INF_OR_UNB.
    for (i in seq_along(problem@constraints)) {
      con <- problem@constraints[[i]]
      dual_val <- dual_at(i)
      if (!is.null(dual_val)) {
        save_dual_value(con, dual_val)
      } else {
        for (dv in con@dual_variables) {
          save_leaf_value(dv, NULL)
        }
      }
    }
    problem@.cache$value <- solution@opt_val
  } else {
    cli_abort("Cannot unpack invalid solution: {solution@status}")
  }

  problem@.cache$status <- solution@status
  problem@.cache$solution <- solution
  invisible(problem)
}

# ==================================================================
# problem_unpack_results -- invert through chain, then unpack
# ==================================================================
## CVXPY SOURCE: problem.py lines 1520-1558

#' Unpack Solver Results into a Problem
#'
#' Inverts the reduction chain and unpacks the raw solver solution into
#' the original problem's variables and constraints. This is step 3 of
#' the decomposed solve pipeline:
#' \enumerate{
#'   \item \code{\link{problem_data}()} -- compile the problem
#'   \item \code{\link{solve_via_data}(chain, data)} -- call the solver
#'   \item \code{problem_unpack_results()} -- invert and unpack
#' }
#'
#' After calling this function, variable values are available via
#' \code{\link{value}()} and constraint duals via \code{\link{dual_value}()}.
#'
#' @param problem A \code{\link{Problem}} object.
#' @param solution The raw solver result from \code{\link{solve_via_data}()}.
#' @param chain The \code{SolvingChain} from \code{\link{problem_data}()}.
#' @param inverse_data The inverse data list from \code{\link{problem_data}()}.
#' @returns The problem object (invisibly), with solution unpacked.
#'
#' @seealso \code{\link{problem_data}}, \code{\link{solve_via_data}}
#' @export
problem_unpack_results <- function(problem, solution, chain, inverse_data) {
  solution <- reduction_invert(chain, solution, inverse_data)

  if (solution@status %in% INACCURATE_STATUS) {
    cli_warn(
      "Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with {.code verbose = TRUE} for more information."
    )
  }
  if (solution@status == INFEASIBLE_OR_UNBOUNDED) {
    cli_warn(
      "The problem is either infeasible or unbounded, but the solver cannot tell which. Disable any solver-specific presolve methods and re-solve."
    )
  }
  if (solution@status %in% ERROR_STATUS) {
    cli_abort(
      "Solver {.val {solver_name(chain@solver)}} failed. Try another solver, or solve with {.code verbose = TRUE} for more information."
    )
  }

  problem_unpack(problem, solution)
  problem@.cache$solver_stats <- solver_stats_from_dict(
    problem@.cache$solution@attr,
    solver_name(chain@solver)
  )
  invisible(problem)
}

# ==================================================================
# Problem data validation -- catch NaN/Inf before solver
# ==================================================================

.check_finite <- function(val, label, inf_allowed = FALSE) {
  if (is.null(val) || length(val) == 0L) return(invisible(NULL))
  ## Sparse matrices: check @x slot (non-zero entries only -- O(nnz) not O(n*m))
  nums <- if (inherits(val, "sparseMatrix")) val@x else as.numeric(val)
  ## `anyNA` is TRUE for NaN, which is what upstream's np.isnan catches.
  if (anyNA(nums))
    cli_abort("Problem data {.val {label}} contains NaN values.")
  if (!inf_allowed && any(is.infinite(nums)))
    cli_abort("Problem data {.val {label}} contains Inf values.")
}

## CVXPY SOURCE: solving_chain.py:412-449 (_validate_problem_data).
##
## Two things were missing. First, upstream checks NINE keys --
## [P, Q, C, A, B, F, G, LOWER_BOUNDS, UPPER_BOUNDS] (:428-429) -- where CVXR
## checked four, so a NaN in the QP inequality data or in a bounds vector went
## through unnoticed. Second, and user-visible: upstream defines
##     inf_allowed_keys = {s.B, s.G, s.LOWER_BOUNDS, s.UPPER_BOUNDS}   (:427)
## and NaN-checks those keys only, with the reason stated at :419-420 --
## "users sometimes use inf for unbounded constraints/variables". CVXR ran every
## key through an indiscriminate NaN-or-Inf check, so two problems CVXPY solves
## errored out here:
##     Minimize(sum(x)), [x >= 1, x <= Inf]   CVXPY 2.0    CVXR "b contains Inf"
##     Minimize(sum(y)), [y >= c(1, -Inf)]    CVXPY -Inf   CVXR "b contains Inf"
##
## Key-name mapping: CVXR's QP path spells the four inequality/equality slots
## `A_eq` / `b_eq` / `F_ineq` / `g_ineq` where upstream uses `A` / `b` / `F` /
## `G` (qp_solver.R:114-121 vs qp_solver.py:139-148), and the conic path uses
## `A` / `b`. Both spellings are listed so the same rule covers both paths.
## Upstream's `s.G` is the inequality RHS *vector* (`F x <= G`), i.e. CVXR's
## `g_ineq` -- not a matrix -- which is why it is Inf-allowed.
.validate_problem_data <- function(data) {
  ## Skip validation for non-dict data (e.g., ConstantSolver returns Problem)
  if (!is.list(data)) return(invisible(NULL))

  ## Upstream s.B, s.G, s.LOWER_BOUNDS, s.UPPER_BOUNDS, in CVXR's spellings.
  inf_allowed <- c(SD_B, "b_eq", "g_ineq", "G", LOWER_BOUNDS, UPPER_BOUNDS)
  keys_to_check <- c(SD_P, SD_C, "q", SD_A, SD_B, "A_eq", "b_eq",
                     "F", "F_ineq", "G", "g_ineq",
                     LOWER_BOUNDS, UPPER_BOUNDS)

  for (key in keys_to_check) {
    .check_finite(data[[key]], key, inf_allowed = key %in% inf_allowed)
  }
  invisible(NULL)
}

# ==================================================================
# .psolve_via_solver_path -- internal fallback-chain helper
# ==================================================================
## CVXPY SOURCE: problem.py lines 506-552 (_solve_solver_path).
##
## solver_path forms accepted:
##   - character vector: c("OSQP", "CLARABEL")
##   - list of length-1 character names: list("OSQP", "CLARABEL")
##   - list of mixed character / length-2 list entries with per-entry opts:
##       list(list("OSQP", list(max_iter = 1)), "CLARABEL")
##
## Per-CVXPY semantics, user kwargs (passed through ...) take precedence
## over entry-specific opts on collision.  An error during a given
## solver's invocation is caught and the next solver is tried; if every
## solver fails, a SolverError-classed condition is raised with the
## per-solver error messages.

## Sentinel value the tryCatch handler uses to signal "this solver
## failed; try the next one".  Not exported.
.PSOLVE_FAILED <- structure(list(), class = ".psolve_failed")

.psolve_via_solver_path <- function(problem, solver_path,
                                     gp, qcp, verbose, warm_start,
                                     requires_grad, ...) {
  user_kwargs <- list(...)

  ## Accept bare character vector form for ergonomics.
  if (is.character(solver_path)) {
    solver_path <- as.list(solver_path)
  }
  ENTRY_MSG <- paste0(
    "Each {.arg solver_path} entry must be a length-1 character solver name ",
    "or a length-2 list of {.code list(name, opts)}."
  )
  if (!is.list(solver_path)) {
    cli_abort(c("{.arg solver_path} must be a list or character vector.",
                "i" = ENTRY_MSG))
  }
  if (length(solver_path) == 0L) {
    cli_abort("{.arg solver_path} must contain at least one solver.")
  }

  ## Parse + validate every entry up-front, so a malformed `solver_path`
  ## errors out before any compile/solve is attempted (CVXPY's
  ## solvers_invalid_inner_input cases all hit this branch).
  parsed <- lapply(solver_path, function(entry) {
    if (is.character(entry) && length(entry) == 1L && !is.na(entry)) {
      list(name = entry, opts = list())
    } else if (is.list(entry) && length(entry) == 2L &&
               is.character(entry[[1L]]) && length(entry[[1L]]) == 1L &&
               !is.na(entry[[1L]]) && is.list(entry[[2L]])) {
      list(name = entry[[1L]], opts = entry[[2L]])
    } else {
      cli_abort(c(ENTRY_MSG,
                  "x" = "Got entry of class {.cls {paste(class(entry), collapse = '/')}}."))
    }
  })

  errors <- list()
  for (e in parsed) {
    name <- toupper(e$name)
    merged_kwargs <- utils::modifyList(e$opts, user_kwargs)
    result <- tryCatch(
      do.call(psolve, c(
        list(problem = problem, solver = name,
             gp = gp, qcp = qcp, verbose = verbose,
             warm_start = warm_start, requires_grad = requires_grad),
        merged_kwargs
      )),
      error = function(err) {
        errors[[name]] <<- conditionMessage(err)
        .PSOLVE_FAILED
      }
    )
    if (!identical(result, .PSOLVE_FAILED)) {
      ## CVXPY v1.9.0 fix: #3324 -- a solver_path entry "succeeds" only on a
      ## clean OPTIMAL status; any other status (infeasible / unbounded /
      ## inaccurate / solver_error) falls through to the next solver, matching
      ## Problem._solve_solver_path (problem.py:578-586).
      if (identical(status(problem), OPTIMAL)) {
        return(result)
      }
      errors[[name]] <- sprintf("non-optimal status: %s", status(problem))
    }
  }
  cli_abort(
    c("All solvers in {.arg solver_path} failed.",
      stats::setNames(
        sprintf("%s: %s", names(errors), unlist(errors)),
        rep("x", length(errors)))),
    class = "SolverError"
  )
}

# ==================================================================
# psolve -- main solve entry point
# ==================================================================
## CVXPY SOURCE: problem.py _solve() (simplified, non-parametric)

#' Solve a Convex Optimization Problem
#'
#' Solves the problem and returns the optimal objective value. After solving,
#' variable values can be retrieved with \code{\link{value}}, constraint
#' dual values with \code{\link{dual_value}}, and solver information with
#' \code{\link{solver_stats}}.
#'
#' @param problem A \code{\link{Problem}} object.
#' @param solver Character string naming the solver to use (e.g.,
#'   \code{"CLARABEL"}, \code{"SCS"}, \code{"OSQP"}, \code{"HIGHS"}),
#'   or \code{NULL} for automatic selection.
#' @param gp Logical; if \code{TRUE}, solve as a geometric program (DGP).
#' @param qcp Logical; if \code{TRUE}, solve as a quasiconvex program (DQCP)
#'   via bisection. Only needed for non-DCP DQCP problems.
#' @param verbose Logical; if \code{TRUE}, print solver output.
#' @param warm_start Logical; if \code{TRUE}, use the current variable
#'   values as a warm-start point for the solver.
#' @param requires_grad Logical; if \code{TRUE}, route the solve through
#'   the DIFFCP wrapper so \code{\link{backward}()} /
#'   \code{\link{derivative}()} can recover gradients.
#' @param nlp Logical; if \code{TRUE}, solve the problem as a disciplined
#'   nonlinear program (DNLP) using the NLP reduction chain and an NLP solver
#'   (e.g. \code{"UNO"}).  The problem must satisfy \code{\link{is_dnlp}()}.
#' @param enforce_dpp Logical; if \code{TRUE}, raise an error when a
#'   parametrized problem is not DPP instead of compiling it as non-DPP.
#' @param ignore_dpp Logical; if \code{TRUE}, treat a DPP problem as non-DPP
#'   (skip the DPP fast path).
#' @param solver_path Optional fallback chain.  A character vector of
#'   solver names or a list whose entries are either character names or
#'   length-2 \code{list(name, opts)} pairs.  Each solver is tried in
#'   sequence; the first that succeeds returns its result.  If every
#'   solver fails, a \code{SolverError}-classed condition is raised
#'   with the per-solver error messages.  Cannot be combined with
#'   \code{solver}.  Mirrors CVXPY's \code{solver_path} argument.
#' @param ... Solver options passed to \code{\link{solver_opts}()}.
#'   Includes chain-construction options (\code{use_quad_obj}), standard
#'   tolerances (\code{feastol}, \code{reltol}, \code{abstol}, \code{num_iter}),
#'   and solver-specific parameters (e.g., \code{eps_abs}, \code{scip_params}).
#'   See \code{\link{solver_opts}} for details.
#'   For DQCP problems (\code{qcp = TRUE}), additional arguments include
#'   \code{low}, \code{high}, \code{eps}, \code{max_iters}, and
#'   \code{max_iters_interval_search}.
#' @returns The optimal objective value (numeric scalar), or \code{Inf} /
#'   \code{-Inf} for infeasible / unbounded problems.
#'
#' @examples
#' x <- Variable()
#' prob <- Problem(Minimize(x), list(x >= 5))
#' result <- psolve(prob, solver = "CLARABEL")
#'
#' @seealso \code{\link{Problem}}, \code{\link{status}},
#'   \code{\link{solver_stats}}, \code{\link{solver_default_param}}
#' @export
psolve <- function(problem, solver = NULL, gp = FALSE, qcp = FALSE,
                   verbose = FALSE, warm_start = FALSE,
                   requires_grad = FALSE, nlp = FALSE,
                   enforce_dpp = FALSE, ignore_dpp = FALSE,
                   solver_path = NULL, ...) {
  if (!.s7_is(problem, Problem)) {
    cli_abort("{.fn psolve} requires a {.cls Problem} object.")
  }

  ## CVXPY SOURCE: problem.py lines 643-649 -- solver_path dispatch
  ## intercepts before any other routing.  Mutually exclusive with
  ## `solver`; tries each entry in turn until one succeeds.
  if (!is.null(solver_path)) {
    if (!is.null(solver)) {
      cli_abort(c(
        "Cannot specify both {.arg solver} and {.arg solver_path}.",
        "i" = "Use {.arg solver} for a single solver or {.arg solver_path} for a fallback list."
      ))
    }
    return(.psolve_via_solver_path(
      problem, solver_path,
      gp = gp, qcp = qcp, verbose = verbose,
      warm_start = warm_start, requires_grad = requires_grad,
      ...
    ))
  }

  ## Normalize solver name to uppercase (CVXPY accepts case-insensitive names)
  if (!is.null(solver)) {
    solver <- toupper(solver)
  }

  ## requires_grad routes the solve through the DIFFCP solver wrapper,
  ## which captures the raw (x, y, s, D, DT) on the solver cache so
  ## Problem$backward() / Problem$derivative() can recover them.
  ## CVXPY SOURCE: problem.py:1167-1183 (the requires_grad guard).
  if (requires_grad) {
    if (!is.null(solver) && solver != DIFFCP_SOLVER) {
      cli_abort(c(
        "When {.code requires_grad = TRUE} the solver must be {.val DIFFCP} (or unspecified).",
        "i" = "Got {.arg solver} = {.val {solver}}."
      ))
    }
    solver <- DIFFCP_SOLVER
    if (!is_dcp(problem) && !gp) {
      cli_abort(c(
        "{.code requires_grad = TRUE} requires a DCP or DGP problem.",
        "i" = "If the problem is DGP, pass {.code gp = TRUE}."
      ))
    }
  }

  ## -- Validate gp/qcp mutual exclusivity -------------------------
  ## CVXPY SOURCE: problem.py lines 1186-1187
  if (gp && qcp) {
    cli_abort("At most one of {.arg gp} and {.arg qcp} can be {.val TRUE}.")
  }

  ## -- DQCP path: bisection solver --------------------------------
  ## CVXPY SOURCE: problem.py lines 1188-1213
  if (qcp && !is_dcp(problem)) {
    if (!is_dqcp(problem)) {
      cli_abort(c(
        "The problem is not DQCP.",
        "i" = "Check that the objective is quasiconvex (Minimize) or quasiconcave (Maximize), and all constraints are DQCP."
      ), class = "DQCPError")   ## CVXPY: raise error.DQCPError, problem.py:1084
    }
    if (verbose) {
      pkg_ver <- utils::packageVersion("CVXR")
      cli_rule(center = "CVXR v{pkg_ver}")
      cli_inform(c("i" = "Reducing DQCP problem to a one-parameter family of DCP problems, for bisection."))
    }

    ## Build reduction chain
    reductions <- list(Dqcp2Dcp())
    extra_args <- list(...)

    if (.s7_is(problem@objective, Maximize)) {
      reductions <- c(list(FlipObjective()), reductions)
      ## FlipObjective negates the objective, so flip the bisection bounds.
      ## Must clear originals first: if user provides only high, the stale
      ## high must not remain (it may violate the flipped parameter's sign).
      low  <- extra_args[["low"]]
      high <- extra_args[["high"]]
      extra_args[["low"]]  <- NULL
      extra_args[["high"]] <- NULL
      if (!is.null(high)) extra_args[["low"]]  <- -high
      if (!is.null(low))  extra_args[["high"]] <- -low
    }

    dqcp_chain <- Chain(reductions = reductions)
    chain_result <- reduction_apply(dqcp_chain, problem)
    reduced_problem <- chain_result[[1L]]
    inverse_data    <- chain_result[[2L]]

    ## Call bisect with forwarded arguments
    bisect_args <- c(
      list(problem = reduced_problem, solver = solver, verbose = verbose),
      extra_args[intersect(names(extra_args),
                           c("low", "high", "eps", "max_iters",
                             "max_iters_interval_search"))]
    )
    soln <- do.call(bisect, bisect_args)

    ## Invert through chain and unpack
    soln <- reduction_invert(dqcp_chain, soln, inverse_data)
    problem_unpack(problem, soln)
    problem@.cache$status <- soln@status

    return(value(problem))
  }

  ## -- NLP path: disciplined nonlinear programming ------------------
  ## CVXPY SOURCE: problem.py:1109-1115. Routes DNLP problems through the NLP
  ## reduction chain + an NLP solver (solve_nlp); errors if nlp = TRUE but the
  ## problem is not DNLP.
  if (nlp && is_dnlp(problem)) {
    return(solve_nlp(problem, solver, warm_start, verbose, ...))
  } else if (nlp && !is_dnlp(problem)) {
    cli_abort(c(
      "The problem you specified is not DNLP.",
      "i" = "{.code nlp = TRUE} requires a disciplined nonlinear program (see {.fn is_dnlp})."
    ), class = "DNLPError")   ## CVXPY: raise error.DNLPError, problem.py:1115
  }

  ## -- Verbose header ----------------------------------------------
  ## PERFORMANCE (2026-08-13). `packageVersion()` was computed here on EVERY
  ## solve although its only use is the header below. It reads and parses
  ## DESCRIPTION off disk (packageDescription -> file.exists + read.dcf):
  ## measured 0.43ms per solve, 4.7% of the `qp` bench cell (A/B with the call
  ## stubbed, N=400, drift -1.3%). The DQCP branch at line ~900 already scoped
  ## it correctly; this one did not.
  if (verbose) {
    pkg_ver <- utils::packageVersion("CVXR")
    cli_rule(center = "CVXR v{pkg_ver}")
    nvars <- length(variables(problem))
    ncons <- length(problem@constraints)
    prob_type <- if (gp) "DGP"
      else if (is_lp(problem)) "LP"
      else if (is_qp(problem)) "QP"
      else if (is_dcp(problem)) "DCP"
      else "non-DCP"
    cli_alert_info("Problem: {nvars} variable{?s}, {ncons} constraint{?s} ({prob_type})")
  }

  ## -- Build solver opts and compile --------------------------------
  opts <- solver_opts(...)
  t0 <- proc.time()
  chain <- .compile(problem, solver, gp = gp, opts = opts,
                    enforce_dpp = enforce_dpp, ignore_dpp = ignore_dpp)

  ## DPP fast path: if we have a cached param_prog, skip full chain apply.
  ## On first solve, we split the chain into pre-solver + solver so we can
  ## intercept the ConeMatrixStuffing data dict (which contains SD_PARAM_PROB)
  ## before the solver's reduction_apply strips it.
  ## CVXPY SOURCE: problem.py lines 811-860
  n_red <- length(chain@reductions)

  if (!is.null(problem@.cache$param_prog)) {
    ## Fast path: re-apply parameters to cached tensor
    if (verbose) cli_alert_info("Using cached DPP tensor (fast path)")
    ## CVXPY SOURCE: problem.py lines 820-821
    ## Update parameter values for reductions that transform them
    ## (e.g., Dgp2Dcp applies log() to parameter values)
    for (red in chain@reductions) {
      update_parameters(red, problem)
    }
    pp <- problem@.cache$param_prog
    pp_result <- apply_parameters(pp, quad_obj = !is.null(pp@P_tensor))
    ## Rebuild ConeMatrixStuffing-format data dict
    cms_data <- list()
    cms_data[[SD_C]] <- pp_result$c
    cms_data[[SD_OFFSET]] <- pp_result$d
    cms_data[[SD_A]] <- pp_result$A
    cms_data[[SD_B]] <- pp_result$b
    cms_data[[LOWER_BOUNDS]] <- pp_result$lower_bounds
    cms_data[[UPPER_BOUNDS]] <- pp_result$upper_bounds
    if (!is.null(pp@P_tensor)) cms_data[[SD_P]] <- pp_result$P
    cms_data[[SD_DIMS]] <- pp@cone_dims
    cms_data[["constraints"]] <- pp@constraints
    cms_data[["x_id"]] <- pp@x_id
    cms_data[[SD_BOOL_IDX]] <- problem@.cache$compile_bool_idx
    cms_data[[SD_INT_IDX]] <- problem@.cache$compile_int_idx
    ## Carry the cached param_prog so the solver's reduction_apply can
    ## see `param_prog@formatted = TRUE` and short-circuit
    ## `format_constraints` (A and b from `apply_parameters` on the
    ## post-format A_tensor are already in cone-interleaved layout).
    cms_data[[SD_PARAM_PROB]] <- pp
    ## Apply solver's reduction_apply to format data for solver
    solver_result <- reduction_apply(chain@reductions[[n_red]], cms_data)
    data <- solver_result[[1L]]
    inverse_data <- c(problem@.cache$compile_inverse_data,
                      list(solver_result[[2L]]))
  } else {
    ## Full path: apply pre-solver reductions, then solver separately
    ## so we can intercept SD_PARAM_PROB before solver strips it.
    pre_data <- problem
    pre_inv_data <- list()
    for (i in seq_len(n_red - 1L)) {
      result <- reduction_apply(chain@reductions[[i]], pre_data)
      pre_data <- result[[1L]]
      pre_inv_data <- c(pre_inv_data, list(result[[2L]]))
    }
    ## Check if we can cache for DPP fast path
    ## CVXPY SOURCE: problem.py lines 840-860
    safe_to_cache <- is.list(pre_data) &&
      !is.null(pre_data[[SD_PARAM_PROB]]) &&
      !any(vapply(chain@reductions, function(r) .s7_is(r, EvalParams), logical(1L)))
    ## Apply solver's reduction_apply
    solver_result <- reduction_apply(chain@reductions[[n_red]], pre_data)
    data <- solver_result[[1L]]
    inverse_data <- c(pre_inv_data, list(solver_result[[2L]]))
    ## Cache the POST-format param_prog (its A_tensor has been row-permuted
    ## by ConicSolver.format_constraints to match the solver's cone-interleaved
    ## row ordering, with `formatted = TRUE` set).  Fall back to the pre-solver
    ## param_prog if the chosen solver did not produce one (e.g. QP path).
    if (safe_to_cache) {
      problem@.cache$param_prog <- data[[SD_PARAM_PROB]] %||%
                                   pre_data[[SD_PARAM_PROB]]
      problem@.cache$compile_inverse_data <- pre_inv_data
      problem@.cache$compile_bool_idx <- pre_data[[SD_BOOL_IDX]]
      problem@.cache$compile_int_idx <- pre_data[[SD_INT_IDX]]
    }
  }

  compile_elapsed <- (proc.time() - t0)[["elapsed"]]
  problem@.cache$compile_time <- compile_elapsed

  if (verbose) {
    solver_nm <- solver_name(chain@solver)
    chain_names <- vapply(chain@reductions, function(r) class(r)[1L], character(1L))
    cli_alert_info("Compilation: {.val {solver_nm}} via {paste(chain_names, collapse = ' -> ')}")
    cli_alert_info("Compile time: {round(compile_elapsed, 4)}s")
  }

  ## Validate solver data before passing to solver
  .validate_problem_data(data)

  ## -- Solver invocation -------------------------------------------
  if (verbose) cli_rule(center = "Numerical solver")
  t0 <- proc.time()
  solver_nm <- solver_name(chain@solver)
  solver_params <- .build_solver_params(solver_nm, opts)
  raw_result <- solve_via_data(chain, data, warm_start, verbose, solver_params,
                               problem = problem)
  solve_elapsed <- (proc.time() - t0)[["elapsed"]]
  problem@.cache$solve_time <- solve_elapsed

  ## Invert through chain and unpack
  problem_unpack_results(problem, raw_result, chain, inverse_data)

  ## -- Verbose summary ---------------------------------------------
  if (verbose) {
    cli_rule(center = "Summary")
    prob_status <- status(problem)
    opt_val <- value(problem)
    cli_alert_success("Status: {prob_status}")
    cli_alert_success("Optimal value: {format(opt_val, digits = 6)}")
    cli_alert_info("Compile time: {round(compile_elapsed, 4)}s")
    cli_alert_info("Solver time: {round(solve_elapsed, 4)}s")
  }

  ## Return optimal value
  value(problem)
}

# ==================================================================
# .make_cvxr_result -- backward-compatible result object
# ==================================================================
## Returns an S3 "cvxr_result" list mimicking old CVXR's solve() return.
## $value and $status work silently; $getValue() and $getDualValue()
## emit one-time deprecation warnings pointing to the new API.

.make_cvxr_result <- function(problem, solver_name) {
  result <- list(
    value  = value(problem),
    status = status(problem),
    solver = solver_name
  )

  ## Closure captures the problem environment
  result$getValue <- function(object) {
    cli_warn(
      c("{.fn getValue} is deprecated.",
        "i" = "Use {.code value(x)} after solving instead."),
      .frequency = "once",
      .frequency_id = "cvxr_getValue_deprecated"
    )
    value(object)
  }

  result$getDualValue <- function(object) {
    cli_warn(
      c("{.fn getDualValue} is deprecated.",
        "i" = "Use {.code dual_value(constraint)} after solving instead."),
      .frequency = "once",
      .frequency_id = "cvxr_getDualValue_deprecated"
    )
    dual_value(object)
  }

  class(result) <- "cvxr_result"
  result
}

#' @export
print.cvxr_result <- function(x, ...) {
  cat(sprintf("Solver: %s\nStatus: %s\nOptimal value: %s\n",
              x$solver, x$status, format(x$value, digits = 6)))
  invisible(x)
}

# -- Accessors ----------------------------------------------------

#' Get Solver Statistics
#'
#' Returns solver statistics from the most recent solve, including
#' solve time, setup time, and iteration count.
#'
#' @param x A \code{\link{Problem}} object.
#' @returns A \code{SolverStats} object, or \code{NULL} if the problem
#'   has not been solved.
#' @export
solver_stats <- function(x) {
  if (!.s7_is(x, Problem)) {
    cli_abort("{.fn solver_stats} requires a {.cls Problem} object.")
  }
  x@.cache$solver_stats
}

#' Get the Raw Solution Object
#'
#' Returns the raw \code{Solution} object from the most recent solve,
#' containing primal and dual variable values, status, and solver
#' attributes.
#'
#' @param x A \code{\link{Problem}} object.
#' @returns A \code{Solution} object, or \code{NULL} if the problem
#'   has not been solved.
#' @export
solution <- function(x) {
  if (!.s7_is(x, Problem)) {
    cli_abort("{.fn solution} requires a {.cls Problem} object.")
  }
  x@.cache$solution
}

#' Get the Raw Solution Object (deprecated)
#'
#' `r lifecycle::badge("deprecated")`
#'
#' Use \code{\link{solution}} instead.
#'
#' @param x A \code{\link{Problem}} object.
#' @returns A \code{Solution} object, or \code{NULL} if the problem
#'   has not been solved.
#' @seealso \code{\link{solution}}
#' @export
problem_solution <- function(x) {
  cli_warn("{.fn problem_solution} is deprecated. Use {.fn solution} instead.",
           .frequency = "once", .frequency_id = "cvxr_problem_solution_deprecated")
  solution(x)
}

# -- backward / derivative (Phase 4.5) ----------------------------
## CVXPY SOURCE: problem.py:1258-1471.

#' Compute the gradient of a solution with respect to Parameters
#'
#' Differentiates through the solution map of `problem`: populates
#' the `gradient` slot of each `Parameter` with the sensitivity of
#' a scalar-valued function of the variables (defaulting to the
#' sum-of-x loss; override per variable by setting
#' `gradient(variable) <-` before calling) with respect to that
#' parameter.  Mirrors `cvxpy.Problem.backward()`.
#'
#' Must be called after [psolve()] with `requires_grad = TRUE`.
#'
#' @param problem A solved `Problem`.
#' @returns The `problem` (for piping); side-effect sets
#'   `gradient(param)` on each parameter.
#' @seealso [derivative()], [psolve()], [gradient()]
#' @export
backward <- function(problem) {
  if (!.s7_is(problem, Problem))
    cli_abort("{.fn backward} requires a {.cls Problem} object.")
  cache <- .get_solver_cache(problem)
  if (!exists(DIFFCP_SOLVER, envir = cache))
    cli_abort(c("{.fn backward} can only be called after",
                "i" = "{.code psolve(problem, requires_grad = TRUE)}."))
  diffcp_raw <- get(DIFFCP_SOLVER, envir = cache)
  status <- problem@.cache$status
  if (!is.null(status) && !(status %in% SOLUTION_PRESENT)) {
    cli_abort("Cannot backpropagate through an infeasible / unbounded problem.")
  }

  DT <- diffcp_raw$DT
  zeros_y <- numeric(length(diffcp_raw$y))
  zeros_s <- numeric(length(diffcp_raw$s))

  param_prog <- problem@.cache$param_prog
  reductions <- problem@.cache$compile_chain@reductions

  ## CVXPY SOURCE: problem.py backward() (#3147 part A -- dict-in/dict-out).
  ## Seed variable gradients in the OUTER (original) representation, keyed by
  ## var id, as SHAPED arrays. NULL on a Variable => all-ones ("sum-of-x loss").
  ## Keep-dims (ADR D_19.5 addendum 2): never flatten until the split_adjoint
  ## boundary -- a bare vector re-enables silent recycling.
  del_vars <- list()
  for (v in variables(problem)) {
    g <- gradient(v)
    if (is.null(g)) g <- array(1, dim = v@shape)
    del_vars[[as.character(v@id)]] <- g
  }
  ## Forward through the chain (outer -> inner): each reduction's var_backward
  ## transforms the whole dict at once.
  for (red in reductions) del_vars <- var_backward(red, del_vars)

  ## Flatten boundary -> conic adjoint (DT) -> parameter Jacobian.
  dx <- split_adjoint(param_prog, del_vars)
  dA_db_dc <- DT(dx, zeros_y, zeros_s)
  ## Note the sign on dA mirrors CVXPY (problem.py:1369: -dA).
  dparams <- apply_param_jac(param_prog,
                              dA_db_dc$dc, -dA_db_dc$dA, dA_db_dc$db)

  ## Reverse through the chain (inner -> outer): param_backward.
  for (i in rev(seq_along(reductions))) {
    dparams <- param_backward(reductions[[i]], dparams)
  }

  ## Write back, one gradient per original parameter (default zero).
  for (p in parameters(problem)) {
    grad <- dparams[[as.character(p@id)]]
    if (is.null(grad)) grad <- array(0, dim = p@shape)
    gradient(p) <- grad
  }

  invisible(problem)
}

#' Apply the derivative of the solution map to perturbations
#'
#' Forward-mode counterpart of [backward()]: reads `delta(param)` for
#' each parameter, applies the cone-program derivative, and writes
#' the predicted change in each variable's optimum to `delta(var)`.
#' Mirrors `cvxpy.Problem.derivative()`.
#'
#' Must be called after [psolve()] with `requires_grad = TRUE`.
#'
#' @inheritParams backward
#' @returns The `problem` (for piping); side-effect sets
#'   `delta(variable)` on each variable.
#' @seealso [backward()], [psolve()], [delta()]
#' @export
derivative <- function(problem) {
  if (!.s7_is(problem, Problem))
    cli_abort("{.fn derivative} requires a {.cls Problem} object.")
  cache <- .get_solver_cache(problem)
  if (!exists(DIFFCP_SOLVER, envir = cache))
    cli_abort(c("{.fn derivative} can only be called after",
                "i" = "{.code psolve(problem, requires_grad = TRUE)}."))
  diffcp_raw <- get(DIFFCP_SOLVER, envir = cache)
  status <- problem@.cache$status
  if (!is.null(status) && !(status %in% SOLUTION_PRESENT)) {
    cli_abort("Cannot apply derivative on an infeasible / unbounded problem.")
  }

  D <- diffcp_raw$D
  param_prog <- problem@.cache$param_prog
  reductions <- problem@.cache$compile_chain@reductions

  if (length(parameters(problem)) == 0L) {
    for (v in variables(problem)) delta(v) <- array(0, dim = v@shape)
    return(invisible(problem))
  }

  ## CVXPY SOURCE: problem.py derivative() (#3147 part A -- dict-in/dict-out).
  ## Seed parameter deltas in the OUTER representation, keyed by param id, as
  ## SHAPED arrays (0 where unperturbed -- every stuffed param must be covered,
  ## else get_parameter_vector falls back to the live value). Keep-dims until
  ## the get_parameter_vector boundary.
  param_deltas <- list()
  for (p in parameters(problem)) {
    d <- delta(p)
    if (is.null(d)) d <- array(0, dim = p@shape)
    param_deltas[[as.character(p@id)]] <- d
  }
  ## Forward through the chain (outer -> inner): param_forward.
  for (red in reductions) param_deltas <- param_forward(red, param_deltas)

  ## Re-run the parameter -> data tensor multiply with deltas as values and the
  ## constant offset zeroed; conic forward derivative (D).
  ap <- apply_parameters(param_prog,
                         id_to_param_value = param_deltas,
                         zero_offset = TRUE)
  dx_dy_ds <- D(-ap$A, ap$b, ap$c)

  ## Split, then reverse through the chain (inner -> outer): var_forward.
  dvars <- split_solution(param_prog, dx_dy_ds$dx)
  for (i in rev(seq_along(reductions))) {
    dvars <- var_forward(reductions[[i]], dvars)
  }

  ## Write back, one delta per original variable (default zero).
  for (v in variables(problem)) {
    dv <- dvars[[as.character(v@id)]]
    if (is.null(dv)) dv <- array(0, dim = v@shape)
    delta(v) <- dv
  }

  invisible(problem)
}

# -- Problem arithmetic ----------------------------------------------
## CVXPY SOURCE: problem.py lines 1593-1634

## Register S3 Ops for Problem so +, -, *, / work
## We use .onLoad registration (same pattern as Expression)

#' @keywords internal
.problem_Ops_handler <- function(e1, e2) {
  op <- .Generic
  unary <- (nargs() == 1L)

  if (unary && op == "-") {
    return(Problem(.negate_objective(e1@objective), e1@constraints))
  }
  if (unary) cli_abort("Unary {.val {op}} not supported on Problem objects.")

  switch(op,
    "+" = {
      if (is.numeric(e1) && e1 == 0) return(e2)
      if (is.numeric(e2) && e2 == 0) return(e1)
      if (!.s7_is(e1, Problem) || !.s7_is(e2, Problem))
        cli_abort("Can only add two {.cls Problem} objects.")
      Problem(.add_objectives(e1@objective, e2@objective),
              unique_list(c(e1@constraints, e2@constraints)))
    },
    "-" = {
      if (is.numeric(e1) && e1 == 0) return(Problem(.negate_objective(e2@objective), e2@constraints))
      if (!.s7_is(e1, Problem) || !.s7_is(e2, Problem))
        cli_abort("Can only subtract two {.cls Problem} objects.")
      Problem(.sub_objectives(e1@objective, e2@objective),
              unique_list(c(e1@constraints, e2@constraints)))
    },
    "*" = {
      if (.s7_is(e1, Problem) && is.numeric(e2)) {
        Problem(.mul_objective(e1@objective, e2), e1@constraints)
      } else if (is.numeric(e1) && .s7_is(e2, Problem)) {
        Problem(.mul_objective(e2@objective, e1), e2@constraints)
      } else {
        cli_abort("Problem can only be multiplied by a numeric scalar.")
      }
    },
    "/" = {
      if (!.s7_is(e1, Problem) || !is.numeric(e2))
        cli_abort("Problem can only be divided by a numeric scalar.")
      Problem(.div_objective(e1@objective, e2), e1@constraints)
    },
    cli_abort("Operator {.val {op}} not supported on Problem objects.")
  )
}