R/132_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

# -- 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()) {
    ## Validate objective type
    ## CVXPY SOURCE: problem.py lines 163-164
    if (!S7_inherits(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_inherits(ci, Constraint)) {
        cli_abort("Element {i} of constraints is not a {.cls Constraint} object.")
      }
    }
    new_object(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_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) {
  ## First: structural DCP compliance in DPP scope
  dcp_ok <- with_dpp_scope({
    is_dcp(x@objective) &&
      all(vapply(x@constraints, is_dcp, logical(1)))
  })
  if (!dcp_ok) return(FALSE)
  ## Second: all atoms in the tree must be DPP-compatible
  if (!is_dpp(x@objective@args[[1L]])) return(FALSE)
  all(vapply(x@constraints, function(c) .all_args(c, is_dpp), 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_inherits(con, Inequality) || S7_inherits(con, NonPos) ||
        S7_inherits(con, NonNeg)) {
      ## Get the constraint expression
      con_expr <- if (S7_inherits(con, Inequality)) con@.expr else con@args[[1L]]
      if (!is_pwl(con_expr)) return(FALSE)
    } else if (!S7_inherits(con, Equality) && !S7_inherits(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)
  result <- any(vapply(variables(problem), function(v) {
    isTRUE(v@attributes$boolean) || isTRUE(v@attributes$integer)
  }, 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
}

# -- 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_inherits(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)
## Cache key: (solver, gp) -- matching CVXPY's (solver, gp, ignore_dpp, use_quad_obj).

.compile <- function(problem, solver = NULL, gp = FALSE) {
  ## CVXPY SOURCE: problem.py lines 786-806
  cache_key <- list(solver, gp)  # (solver, gp)
  if (!identical(cache_key, problem@.cache$compile_key)) {
    problem@.cache$compile_key <- cache_key
    problem@.cache$compile_chain <- construct_solving_chain(problem, solver,
                                                            gp = gp)
    problem@.cache$param_prog <- NULL
    problem@.cache$compile_inverse_data <- NULL
    ## Clear solver-specific cache when chain changes
    problem@.cache$solver_cache <- NULL
  }
  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, ...) {
  chain <- .compile(x, solver)
  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, ...) {
  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, ...)
}

# ==================================================================
# 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) {
    new_object(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) {
  if (solution@status %in% SOLUTION_PRESENT) {
    for (v in variables(problem)) {
      save_leaf_value(v, solution@primal_vars[[as.character(v@id)]])
    }
    for (con in problem@constraints) {
      dual_val <- solution@dual_vars[[as.character(con@id)]]
      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)
    }
    for (con in problem@constraints) {
      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) {
  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)
  if (anyNA(nums))
    cli_abort("Problem data {.val {label}} contains NaN values.")
  if (any(is.infinite(nums)))
    cli_abort("Problem data {.val {label}} contains Inf values.")
}

.validate_problem_data <- function(data) {
  ## CVXPY SOURCE: solving_chain.py lines 414-416
  ## Skip validation for non-dict data (e.g., ConstantSolver returns Problem)
  if (!is.list(data)) return(invisible(NULL))
  .check_finite(data[[SD_A]], SD_A)
  .check_finite(data[[SD_B]], SD_B)
  .check_finite(data[[SD_C]], SD_C)
  .check_finite(data[[SD_P]], SD_P)
}

# ==================================================================
# 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 feastol Numeric or \code{NULL}; feasibility tolerance. Mapped to
#'   the solver-specific parameter name (e.g., \code{tol_feas} for Clarabel,
#'   \code{eps_prim_inf} for OSQP). If \code{NULL} (default), the solver's
#'   own default is used.
#' @param reltol Numeric or \code{NULL}; relative tolerance. Mapped to the
#'   solver-specific parameter name (e.g., \code{tol_gap_rel} for Clarabel,
#'   \code{eps_rel} for OSQP/SCS). If \code{NULL} (default), the solver's
#'   own default is used.
#' @param abstol Numeric or \code{NULL}; absolute tolerance. Mapped to the
#'   solver-specific parameter name (e.g., \code{tol_gap_abs} for Clarabel,
#'   \code{eps_abs} for OSQP/SCS). If \code{NULL} (default), the solver's
#'   own default is used.
#' @param num_iter Integer or \code{NULL}; maximum number of solver
#'   iterations. Mapped to the solver-specific parameter name (e.g.,
#'   \code{max_iter} for Clarabel/OSQP, \code{max_iters} for SCS). If
#'   \code{NULL} (default), the solver's own default is used.
#' @param ... Additional solver-specific options passed directly to the
#'   solver. If a solver-native parameter name conflicts with a standard
#'   parameter (e.g., both \code{feastol = 1e-3} and \code{tol_feas = 1e-6}
#'   are given), the solver-native name in \code{...} takes priority.
#'   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,
                   feastol = NULL, reltol = NULL, abstol = NULL,
                   num_iter = NULL, ...) {
  if (!S7_inherits(problem, Problem)) {
    cli_abort("{.fn psolve} requires a {.cls Problem} object.")
  }

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

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

  pkg_ver <- utils::packageVersion("CVXR")

  ## -- Verbose header ----------------------------------------------
  if (verbose) {
    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})")
  }

  ## -- Compilation (chain construction + reductions) ---------------
  t0 <- proc.time()
  chain <- .compile(problem, solver, gp = gp)

  ## 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
    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
    ## 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 (before solver transforms data)
    ## 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_inherits(r, EvalParams), logical(1L)))
    if (safe_to_cache) {
      problem@.cache$param_prog <- 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]]
    }
    ## 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]]))
  }

  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_opts <- .apply_std_params(solver_nm, list(...),
                                   feastol, reltol, abstol, num_iter)
  raw_result <- solve_via_data(chain, data, warm_start, verbose, solver_opts,
                               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_inherits(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_inherits(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)
}

## DEFERRED: Derivative/sensitivity API (problem.py lines 1258-1470)
## These methods enable differentiating through the solution map:
##
## backward() -- compute gradients of parameters w.r.t. objective via diffcp
## derivative() -- apply forward derivatives (Jacobian-vector products)
## requires_grad parameter in psolve() -- triggers diffcp solver path
##
## Deferred because the core dependency (diffcp) has no R equivalent.
## See notes/derivative_api_deferred.md for rationale and future path.

# -- 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_inherits(e1, Problem) || !S7_inherits(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_inherits(e1, Problem) || !S7_inherits(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_inherits(e1, Problem) && is.numeric(e2)) {
        Problem(.mul_objective(e1@objective, e2), e1@constraints)
      } else if (is.numeric(e1) && S7_inherits(e2, Problem)) {
        Problem(.mul_objective(e2@objective, e1), e2@constraints)
      } else {
        cli_abort("Problem can only be multiplied by a numeric scalar.")
      }
    },
    "/" = {
      if (!S7_inherits(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.")
  )
}