R/274_reductions_solvers_conic_solvers_diffcp_conif.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: cvxpy/reductions/solvers/conic_solvers/diffcp_conif.py
##
## DIFFCP solver wrapper.  Inherits canonicalisation from SCS (so the
## R-level apply / reduction_apply machinery is reused) and overrides
## `solve_via_data` to call the R `diffcp::solve_and_derivative`,
## which returns the optimum (x, y, s) together with the forward and
## adjoint derivative closures D / DT.

## NOTE: diffcp is an optional dependency (Enhances, not Imports), so it is
## deliberately NOT imported into the NAMESPACE.  solve_and_derivative is
## called fully qualified as diffcp::solve_and_derivative under a
## requireNamespace() guard (see solve_via_data below).
##
## CVXPY's DIFFCP class supports both SCS and Clarabel as the inner
## forward solver.  Mirroring that, the R version accepts a
## `solve_method` solver-opts entry; defaults to Clarabel (matches
## diffcp R's default).

# -- DIFFCP status map --------------------------------------------
## CVXPY SOURCE: diffcp_conif.py:32-41
DIFFCP_STATUS_MAP <- list(
  "Solved"               = OPTIMAL,
  "Solved/Inaccurate"    = OPTIMAL_INACCURATE,
  "Optimal Inaccurate"   = OPTIMAL_INACCURATE,    # Clarabel "AlmostSolved"
  "Unbounded"            = UNBOUNDED,
  "Unbounded/Inaccurate" = UNBOUNDED_INACCURATE,
  "Unbounded Inaccurate" = UNBOUNDED_INACCURATE,
  "Infeasible"           = INFEASIBLE,
  "Infeasible/Inaccurate"= INFEASIBLE_INACCURATE,
  "Infeasible Inaccurate"= INFEASIBLE_INACCURATE,
  "Failure"              = SOLVER_ERROR,
  "Indeterminate"        = SOLVER_ERROR,
  "Interrupted"          = SOLVER_ERROR
)

# -- DIFFCP_Solver class ------------------------------------------
## CVXPY SOURCE: diffcp_conif.py:28-46

DIFFCP_Solver <- new_class("DIFFCP_Solver", parent = SCS_Solver,
                           package = "CVXR",
  constructor = function() {
    if (FALSE) new_object(S7_object())  ## S7 static-check guard
    .fast_new(DIFFCP_Solver, S7_object(),
      .cache = new.env(parent = emptyenv()),
      MIP_CAPABLE = FALSE,
      BOUNDED_VARIABLES = FALSE,
      ## Same svec format as SCS, whose interface this class extends and whose
      ## the shared `extract_dual_value` it reuses below.  MUST be stated explicitly:
      ## the PSD format used to be an S7 method override that DIFFCP inherited
      ## from SCS_Solver, but it is now a property, and `.fast_new` gives a
      ## subclass NOTHING it does not name (constraint 17).
      PSD_TRIANGLE_KIND = TriangleKind$LOWER,
      PSD_SQRT2_SCALING = TRUE,
      ## Same cone set as SCS, whose data format DIFFCP reuses.
      SUPPORTED_CONSTRAINTS = list(Zero, NonNeg, SOC, ExpCone,
                                   SvecPSD, PowCone3D),
      EXP_CONE_ORDER = c(0L, 1L, 2L),
      REQUIRES_CONSTR = TRUE
    )
  }
)

method(solver_name, DIFFCP_Solver) <- function(x) DIFFCP_SOLVER

## CVXPY SOURCE: diffcp_conif.py:60-63 — DIFFCP itself does not
## support quadratic objectives via dense / lsqr derivative modes;
## the R diffcp port mirrors that constraint.
method(supports_quad_obj, DIFFCP_Solver) <- function(x) FALSE

# -- DIFFCP reduction_apply ---------------------------------------
## CVXPY SOURCE: diffcp_conif.py:65-78 -- DIFFCP overrides `apply` for exactly
## one reason, and states it in a comment:
##
##     # Keep zeros in A that are affected by parameters
##     c, d, A, b = problem.apply_parameters(keep_zeros=True)
##
## Without it, an entry of `A` driven by a parameter whose CURRENT VALUE IS
## ZERO multiplies out to zero and drops out of the sparsity pattern, so diffcp
## has no entry to differentiate with respect to and the gradient for that
## parameter comes back as exactly 0 -- no error, no warning. Measured on
## `max x1 + 2*x2  s.t.  p*x1 + x2 <= 1, 0 <= x1 <= 3, 0 <= x2 <= 5`
## (closed form x1 = 3, x2 = 1 - 3p, so d/dp of sum(x) is -3 everywhere):
##       p = 0     CVXR 0   CVXPY -3.000000
##       p = 0.1   CVXR -3  CVXPY -3.000004
## CVXR inherited SCS's reduction_apply and so could not ask for it.
##
## Recomputing `A` here rather than at the stuffing step keeps the request
## where upstream puts it -- on the solver that needs it -- and works on both of
## CVXR's paths, because `data[[SD_PARAM_PROB]]` and `data[[SD_A]]` are always
## in the same layout: pre-format on the first compile (ConicSolver's
## `format_constraints` then permutes both, and a row permutation preserves
## explicit zeros), post-format on the DPP fast path.
method(reduction_apply, DIFFCP_Solver) <- function(x, problem, ...) {
  data <- problem   ## the data list from ConeMatrixStuffing, as for ConicSolver
  pp <- data[[SD_PARAM_PROB]]
  if (!is.null(pp) && !is.null(pp@A_tensor)) {
    ap <- apply_parameters(pp, quad_obj = FALSE, keep_zeros = TRUE)
    data[[SD_A]] <- ap$A
    data[[SD_B]] <- ap$b
  }
  ## Delegate to the conic machinery (S7 has no callNextMethod; call the
  ## parent's method object directly).
  method(reduction_apply, ConicSolver)(x, data, ...)
}

# -- DIFFCP solve_via_data ----------------------------------------
## CVXPY SOURCE: diffcp_conif.py:129-184
##
## Calls `diffcp::solve_and_derivative` and returns the raw result
## dictionary (x, y, s, D, DT, info, solve_method) for `reduction_invert`
## to pick apart.  The D and DT closures live on the result so the
## caller (Problem$backward / Problem$derivative) can pluck them.
method(solve_via_data, DIFFCP_Solver) <- function(x, data, warm_start = FALSE,
                                                  verbose = FALSE,
                                                  solver_opts = list(), ...) {
  if (!requireNamespace("diffcp", quietly = TRUE)) {
    cli_abort("Package {.pkg diffcp} is required for {.code requires_grad = TRUE}.")
  }

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

  ## CVXPY scs_conif's apply yields A, b, c, dims; we reuse the same
  ## (data layout matches SCS).
  A <- data[[SD_A]]
  b <- data[[SD_B]]
  c <- data[[SD_C]]
  cone_dict <- dims_to_solver_dict_scs(data[[SD_DIMS]])

  ## solve_method: SCS or Clarabel.  Default to Clarabel (the user
  ## "SCS not preferred" preference; also matches diffcp R's default).
  solve_method <- solver_opts[["solve_method"]] %||% CLARABEL_SOLVER

  ## Tolerance defaults match Python diffcp_conif.
  if (toupper(solve_method) == SCS_SOLVER) {
    if (is.null(solver_opts[["eps_abs"]])) solver_opts[["eps_abs"]] <- 1e-5
    if (is.null(solver_opts[["eps_rel"]])) solver_opts[["eps_rel"]] <- 1e-5
  }

  warm_start_arg <- NULL
  if (toupper(solve_method) == SCS_SOLVER && warm_start &&
      !is.null(solver_cache) && exists(DIFFCP_SOLVER, envir = solver_cache)) {
    cached <- get(DIFFCP_SOLVER, envir = solver_cache)
    warm_start_arg <- list(cached$x, cached$y, cached$s)
  }

  ## Strip our own solve_method/mode from solver_opts before forwarding
  ## the rest as solver-control kwargs.
  forward_opts <- solver_opts
  forward_opts[["solve_method"]] <- NULL
  mode <- forward_opts[["mode"]] %||% "lsqr"
  forward_opts[["mode"]] <- NULL

  start <- Sys.time()
  res <- tryCatch(
    do.call(diffcp::solve_and_derivative,
            c(list(A, b, c, cone_dict,
                   solve_method = solve_method,
                   mode = mode,
                   warm_start = warm_start_arg),
              forward_opts)),
    error = function(e) list(error = conditionMessage(e))
  )
  end <- Sys.time()

  if (!is.null(res$error)) {
    return(list(
      info = list(status = "Failure"),
      solve_method = solve_method,
      error = res$error,
      TOT_TIME = as.numeric(difftime(end, start, units = "secs"))
    ))
  }

  res$solve_method <- solve_method
  res$TOT_TIME <- as.numeric(difftime(end, start, units = "secs"))

  if (!is.null(solver_cache)) {
    assign(DIFFCP_SOLVER, res, envir = solver_cache)
  }
  res
}

# -- DIFFCP invert ------------------------------------------------
## CVXPY SOURCE: diffcp_conif.py:79-127
##
## We dispatch back to SCS_Solver's invert for the common bookkeeping
## (status, primal/dual var packaging) but tag the Solution with the
## raw diffcp output (x, y, s, D, DT) so Problem$backward /
## Problem$derivative can find them.
method(reduction_invert, DIFFCP_Solver) <- function(x, solution, inverse_data, ...) {
  ## CVXPY SOURCE: diffcp_conif.py:79-99 -- status mapping is
  ## *solve_method-aware*. This MUST mirror CVXPY exactly; do not collapse
  ## it back to a single string lookup (see the long note below + ADR
  ## D_DIFFCP.1 in notes/decisions.md / notes/diffcp_status_mapping.md).
  ##
  ## WHY two branches (the trap this guards against):
  ##   The R `diffcp` package emits DIFFERENT status shapes per solver,
  ##   faithfully mirroring Python diffcp:
  ##     * SCS path      -> info$status     = R scs string, LOWERCASE
  ##                        ("solved", "infeasible", ...), PLUS the
  ##                        integer info$status_val (1, -2, ...).
  ##     * Clarabel path -> info$status     = canonical SCS-convention
  ##                        STRING ("Solved", "Infeasible", ...), already
  ##                        normalized by R diffcp's CLARABEL2SCS map;
  ##                        NO status_val field.
  ##   CVXPY resolves this by branching on solve_method: SCS is mapped by
  ##   the INTEGER status_val (diffcp_conif.py:86/90), Clarabel/ECOS by the
  ##   STRING (py:94/98). A previous version of this method took a shortcut
  ##   -- one string lookup in DIFFCP_STATUS_MAP for both paths -- whose
  ##   keys are capitalized ("Solved"). That silently mapped the SCS path's
  ##   lowercase "solved" to NULL -> SOLVER_ERROR, so every SCS-backed
  ##   differentiable solve "failed" even though diffcp had solved it. The
  ##   bug hid because the default solve_method is Clarabel, so the SCS
  ##   branch was dead code until a test forced solve_method="SCS".
  ##   Guard: scripts/status_map_audit.R section (C) exercises BOTH paths.
  info <- solution[["info"]]
  solve_method <- solution[["solve_method"]] %||% CLARABEL_SOLVER
  if (toupper(solve_method) == SCS_SOLVER) {
    ## CVXPY SOURCE: diffcp_conif.py:90-92 (scs >= 3.0 branch).
    status <- SCS_STATUS_MAP[[as.character(info[["status_val"]])]]
    solve_time <- info[["solve_time"]]   # R scs field name
    setup_time <- info[["setup_time"]]
  } else {
    ## CVXPY SOURCE: diffcp_conif.py:97-99 (Clarabel branch). The string
    ## here is already in canonical SCS convention (normalized by R diffcp).
    status <- DIFFCP_STATUS_MAP[[as.character(info[["status"]] %||% "Failure")]]
    solve_time <- info[["solveTime"]]    # R diffcp Clarabel field name
    setup_time <- info[["setupTime"]]
  }
  if (is.null(status)) status <- SOLVER_ERROR

  attr_list <- list()
  attr_list[[RK_NUM_ITERS]] <- info[["iter"]] %||% NA_integer_
  if (!is.null(solve_time)) {
    attr_list[[RK_SOLVE_TIME]] <- solve_time
  }
  if (!is.null(setup_time)) {
    attr_list[[RK_SETUP_TIME]] <- setup_time
  }
  ## Stash the raw diffcp output so Problem$backward / Problem$derivative
  ## can recover D, DT, x, y, s without going through invert.
  attr_list[["diffcp_raw"]] <- solution

  if (status %in% SOLUTION_PRESENT) {
    primal_val <- solution[["info"]][["pobj"]]
    if (is.null(primal_val)) primal_val <- sum(inverse_data[[SD_C]] * solution$x)
    opt_val <- primal_val + inverse_data[[SD_OFFSET]]

    primal_vars <- list()
    primal_vars[[as.character(inverse_data[[SOLVER_VAR_ID]])]] <- solution[["x"]]

    y <- solution[["y"]]
    zero_dim <- inverse_data[[SD_DIMS]]@zero
    if (zero_dim > 0L) {
      eq_dual <- get_dual_values(
        y[seq_len(zero_dim)],
        extract_dual_value,
        inverse_data[[SOLVER_EQ_CONSTR]]
      )
    } else {
      eq_dual <- list()
    }
    if (zero_dim < length(y)) {
      ineq_dual <- get_dual_values(
        y[(zero_dim + 1L):length(y)],
        extract_dual_value,
        inverse_data[[SOLVER_NEQ_CONSTR]]
      )
    } else {
      ineq_dual <- list()
    }
    Solution(status, opt_val, primal_vars, c(eq_dual, ineq_dual), attr_list)
  } else {
    failure_solution(status, attr_list)
  }
}

method(print, DIFFCP_Solver) <- function(x, ...) {
  cat("DIFFCP_Solver()\n")
  invisible(x)
}