R/273_reductions_solvers_conic_solvers_clarabel_conif.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: reductions/solvers/conic_solvers/clarabel_conif.py
## PARTIAL PORT: in: cone conversion, status map, invert, solve_via_data,
##   warm start, supports_quad_obj, EXP_CONE_ORDER, and the `accept_unknown`
##   option (clarabel_conif.py:87, 126-134, 176-181) -- all verified
##   entry-for-entry against 1.9.2.
##   out: upstream's `parse_solver_opts()` exists as a separate static method;
##   CVXR inlines the equivalent in solve_via_data, so there is no R function
##   of that name to call. Behavior is ported, structure is not.
##
##   DELIBERATE MECHANISM DEVIATION for `accept_unknown`: upstream reads the
##   flag from `inverse_data.solver_options`, which `Solver.solve()`
##   (solver.py:171-177) can attach because apply/solve/invert happen inside
##   one call. CVXR's solve API is decomposed and PUBLIC -- problem_data() /
##   solve_via_data() / problem_unpack_results() -- so invert only ever sees
##   compile-time inverse_data. The flag therefore rides on the solution,
##   the sole channel between the halves. Same information, same guards, no
##   public signature change.
##
## THE `accept_unknown` GAP DATED TO THE ORIGINAL 1.8.2 PORT and is recorded
## here because the reason it survived is reusable: ACCEPT_UNKNOWN is present
## in 1.8.2, 1.9.0, 1.9.1 and 1.9.2 alike, so it never appeared in any
## upstream DELTA -- and all three parity passes were delta-driven. A gap
## present at the first port is structurally invisible to a diff. Closing the
## rest needs a COMPLETENESS audit of the solver interfaces, not another
## delta. Without this sentinel a half-ported file is indistinguishable from
## a complete one, which is what constraint 15(h) exists to prevent.
##
## Clarabel uses upper-triangular svec format for PSD constraints.
## Convention: A*x + s = b, s in K


# -- Clarabel status map -------------------------------------------
## CVXPY SOURCE: clarabel_conif.py lines 158-169

CLARABEL_STATUS_MAP <- list(
  "Solved"               = OPTIMAL,
  "PrimalInfeasible"     = INFEASIBLE,
  "DualInfeasible"       = UNBOUNDED,
  "AlmostSolved"         = OPTIMAL_INACCURATE,
  "AlmostPrimalInfeasible" = INFEASIBLE_INACCURATE,
  "AlmostDualInfeasible" = UNBOUNDED_INACCURATE,
  "MaxIterations"        = USER_LIMIT,
  "MaxTime"              = USER_LIMIT,
  "NumericalError"       = SOLVER_ERROR,
  "InsufficientProgress" = SOLVER_ERROR,
  ## CVXPY v1.9.0 fix: #3180 -- map the "Unsolved" status to SOLVER_ERROR
  ## instead of falling through to a KeyError.
  "Unsolved"             = SOLVER_ERROR
)

## CVXPY SOURCE: clarabel_conif.py:87 -- `ACCEPT_UNKNOWN = "accept_unknown"`.
## Opt-in: when set, and the solver returned BOTH a primal and a dual iterate,
## `InsufficientProgress` is read as OPTIMAL_INACCURATE rather than
## SOLVER_ERROR.  Clarabel reports InsufficientProgress when it stalls, which
## on a near-degenerate problem can still leave a usable iterate; the default
## remains to reject it, and the caller must ask for the weaker guarantee.
##
## SCOPE -- this reaches a DIRECT solve only. DQCP bisection subproblems go
## through bisection.py:32-42 `_solve(problem, solver)`, which passes the
## solver and NOTHING ELSE; upstream forwards no solver options there, and
## neither does CVXR's .bisect_solve(). So `accept_unknown` cannot rescue a
## stall inside a bisection in either library. That is parity, not a gap --
## do not "fix" .bisect_solve to forward options unless upstream does.
CLARABEL_ACCEPT_UNKNOWN <- "accept_unknown"

# -- dims_to_solver_dict_clarabel --------------------------------
## Converts ConeDims to R clarabel's named-list cone format.
## R clarabel uses: z (zero), l (nonneg), q (soc), s (psd),
## ep (exp), p (power), gp (generalized power).

dims_to_solver_dict_clarabel <- function(cone_dims) {
  cones <- list()
  if (cone_dims@zero > 0L)
    cones[["z"]] <- as.integer(cone_dims@zero)
  if (cone_dims@nonneg > 0L)
    cones[["l"]] <- as.integer(cone_dims@nonneg)
  if (length(cone_dims@soc) > 0L)
    cones[["q"]] <- as.integer(cone_dims@soc)
  if (length(cone_dims@psd) > 0L)
    cones[["s"]] <- as.integer(cone_dims@psd)
  if (cone_dims@exp > 0L)
    cones[["ep"]] <- as.integer(cone_dims@exp)
  if (length(cone_dims@p3d) > 0L)
    cones[["p"]] <- cone_dims@p3d
  if (length(cone_dims@pnd) > 0L) {
    ## Generalized power cones: each needs a unique name (gp1, gp2, ...)
    ## and strict_cone_order = FALSE when calling clarabel_solver.
    for (i in seq_along(cone_dims@pnd)) {
      cones[[paste0("gp", i)]] <- list(a = cone_dims@pnd[[i]], n = 1L)
    }
  }
  cones
}

# -- Clarabel_Solver class ----------------------------------------
## CVXPY SOURCE: clarabel_conif.py lines 136-173

Clarabel_Solver <- new_class("Clarabel_Solver", parent = ConicSolver,
  package = "CVXR",
  constructor = function() {
    if (FALSE) new_object(S7_object())  ## S7 static-check guard
    .fast_new(Clarabel_Solver, S7_object(),
      .cache = new.env(parent = emptyenv()),
      MIP_CAPABLE = FALSE,
      BOUNDED_VARIABLES = FALSE,
      ## CVXPY SOURCE: clarabel_conif.py lines 73-74
      PSD_TRIANGLE_KIND = TriangleKind$UPPER,
      PSD_SQRT2_SCALING = TRUE,
      ## CVXPY SOURCE: clarabel_conif.py:71-72 -- SvecPSD, not PSD: the
      ## packing is done by PSDToSvecPSD before the solver sees it.
      SUPPORTED_CONSTRAINTS = list(Zero, NonNeg, SOC, ExpCone,
                                   PowCone3D, SvecPSD, PowConeND),
      EXP_CONE_ORDER = c(0L, 1L, 2L),
      REQUIRES_CONSTR = FALSE
    )
  }
)

method(solver_name, Clarabel_Solver) <- function(x) CLARABEL_SOLVER

## CVXPY v1.8.2: Clarabel supports quadratic objective with any conic constraints
method(supports_quad_obj, Clarabel_Solver) <- function(x) TRUE

# -- Clarabel invert ----------------------------------------------
## CVXPY SOURCE: clarabel_conif.py lines 245-284

method(reduction_invert, Clarabel_Solver) <- function(x, solution, inverse_data, ...) {
  attr_list <- list()

  ## R clarabel returns integer status; map to string via solver_status_descriptions()
  status_int <- solution$status
  status_names <- names(clarabel::solver_status_descriptions())
  status_str <- if (status_int >= 1L && status_int <= length(status_names)) {
    status_names[status_int]
  } else {
    "Unknown"
  }
  ## CVXPY SOURCE: clarabel_conif.py:126-134 -- take a COPY of the map and
  ## re-point one entry when `accept_unknown` was requested AND the solver
  ## returned both x and z.  Both guards matter: without a primal and a dual
  ## iterate there is nothing to hand back, so the status would be promoted on
  ## an empty solution.
  ##
  ## Upstream reads the flag off `inverse_data.solver_options`, which its
  ## `Solver.solve()` (solver.py:171-177) can attach because it runs apply,
  ## solve and invert in one call.  CVXR's solve API is DECOMPOSED and public
  ## -- problem_data() / solve_via_data() / problem_unpack_results() -- so
  ## invert receives compile-time inverse_data and never sees solver options.
  ## The flag therefore travels on the solution, the only channel between the
  ## two halves; solve_via_data() attaches it below.  Same information, same
  ## precondition, different carrier, and it costs no public signature change.
  status_map <- CLARABEL_STATUS_MAP
  if (isTRUE(solution[[CLARABEL_ACCEPT_UNKNOWN]]) &&
      !is.null(solution$x) && !is.null(solution$z)) {
    status_map[["InsufficientProgress"]] <- OPTIMAL_INACCURATE
  }
  status <- status_map[[status_str]]
  if (is.null(status)) status <- SOLVER_ERROR

  if (!is.null(solution$solve_time))
    attr_list[[RK_SOLVE_TIME]] <- solution$solve_time
  if (!is.null(solution$iterations))
    attr_list[[RK_NUM_ITERS]] <- solution$iterations

  ## Dual variables: recover whenever the solver returns z -- even for an
  ## infeasibility certificate (the dual ray), not just optimal solutions.
  ## CVXPY v1.9.0 fix: #3228 -- propagate infeasibility certificate for CLARABEL.
  ## Split z at the zero-cone boundary into equality / inequality duals.
  dual_vars <- list()
  if (!is.null(solution$z)) {
    z <- solution$z
    zero_dim <- inverse_data[[SD_DIMS]]@zero
    if (zero_dim > 0L) {
      eq_dual <- get_dual_values(
        z[seq_len(zero_dim)],
        extract_dual_value,
        inverse_data[[SOLVER_EQ_CONSTR]]
      )
    } else {
      eq_dual <- list()
    }
    if (zero_dim < length(z)) {
      ineq_dual <- get_dual_values(
        z[(zero_dim + 1L):length(z)],
        extract_dual_value,
        inverse_data[[SOLVER_NEQ_CONSTR]]
      )
    } else {
      ineq_dual <- list()
    }
    dual_vars <- c(eq_dual, ineq_dual)
  }

  if (status %in% SOLUTION_PRESENT) {
    primal_val <- solution$obj_val
    opt_val <- primal_val + inverse_data[[SD_OFFSET]]
    primal_vars <- list()
    primal_vars[[as.character(inverse_data[[SOLVER_VAR_ID]])]] <- solution$x
    Solution(status, opt_val, primal_vars, dual_vars, attr_list)
  } else {
    failure_solution(status, attr_list, dual_vars)
  }
}

# -- Clarabel solve_via_data --------------------------------------
## CVXPY SOURCE: clarabel_conif.py lines 313-388
## Warm-start: persistent solver via clarabel_solver() + solver_update().
## Follows the same cache pattern as OSQP (osqp_qpif.R).

method(solve_via_data, Clarabel_Solver) <- function(x, data, warm_start = FALSE, verbose = FALSE,
                                                    solver_opts = list(), ...) {
  .require_solver_package(CLARABEL_SOLVER)

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

  A <- data[[SD_A]]
  b <- data[[SD_B]]
  q <- data[[SD_C]]
  cones <- dims_to_solver_dict_clarabel(data[[SD_DIMS]])

  ## Build P (quadratic objective)
  ## CVXPY: clarabel_conif.py line 343 -- P = sp.triu(P).tocsc()
  ## R clarabel expects a dsCMatrix (symmetric sparse).
  ## forceSymmetric(triu(P)) stores upper triangle as symmetric.
  nvars <- length(q)
  if (!is.null(data[[SD_P]])) {
    P <- Matrix::forceSymmetric(Matrix::triu(data[[SD_P]]), uplo = "U")
  } else {
    P <- Matrix::sparseMatrix(i = integer(0), j = integer(0), x = numeric(0),
                              dims = c(nvars, nvars))
  }

  ## Parse settings
  settings <- clarabel::clarabel_control()
  settings$verbose <- verbose
  ## BACK-COMPAT SHIM for R clarabel <= 0.11.2, where `reduced_tol_infeas_abs`
  ## defaulted to 5e-5.  Clarabel.rs has used 5e-12 since v0.10.0; the R
  ## package had fallen behind and corrected it in 0.11.3 ("reduced_tol_infeas
  ## _abs is now 5e-12, not 5e-5, matching upstream", clarabel 0.11.3 NEWS).
  ## On 0.11.3+ this line is a NO-OP.  DESCRIPTION allows clarabel (>= 0.11)
  ## and CRAN still ships 0.11.2, so the shim stays until the floor can move.
  ##
  ## Why it matters: the value is a threshold an infeasibility certificate must
  ## EXCEED, so the larger 5e-5 is the STRICTER one and suppresses
  ## `AlmostPrimalInfeasible` entirely -- a near-degenerate LP that Clarabel
  ## could certify as almost-infeasible returns the inconclusive NumericalError
  ## instead.  DQCP bisection needs a conclusive answer per query to narrow its
  ## bracket, so it stalls.  Measured on the subproblem family
  ## a>=1, b<=3, a - t*b <= t-1, over 200 points of t in [0.5-3e-5, 0.5-1e-8]:
  ##      5e-5   ->   0 AlmostPrimalInfeasible, 198/200 inconclusive
  ##      5e-12  ->  69 AlmostPrimalInfeasible
  ## (verified side by side: clarabel 0.11.2 in a temp library vs 0.11.3.)
  settings$reduced_tol_infeas_abs <- 5e-12

  ## CVXPY SOURCE: clarabel_conif.py:176-181 -- `accept_unknown` is a CVXPY
  ## option, not a Clarabel setting, so it is REMOVED from the keys before the
  ## rest are applied.  Leaving it in would push an unknown field into
  ## clarabel_control(), which R accepts silently -- the option would then be
  ## quietly ignored rather than honored, which is the worst of both.
  ## (`use_quad_obj`, which upstream strips here too, never reaches this
  ## function in CVXR: it lives on solver_opts() and is consumed during chain
  ## construction.)
  accept_unknown <- isTRUE(solver_opts[[CLARABEL_ACCEPT_UNKNOWN]])
  solver_opts[[CLARABEL_ACCEPT_UNKNOWN]] <- NULL

  for (opt_name in names(solver_opts)) {
    settings[[opt_name]] <- solver_opts[[opt_name]]
  }

  cache_key <- CLARABEL_SOLVER
  used_warm <- FALSE

  ## -- Warm path --------------------------------------------------
  if (warm_start && !is.null(solver_cache) &&
      exists(cache_key, envir = solver_cache)) {
    cached <- get(cache_key, envir = solver_cache)
    old_solver <- cached$solver
    old_data   <- cached$data

    ## Dimension check: structure must match
    if (length(q) == length(old_data$q) &&
        nrow(A) == nrow(old_data$A) &&
        ncol(A) == ncol(old_data$A)) {

      ## Determine what changed and build update args
      new_P <- NULL
      new_q <- NULL
      new_A <- NULL
      new_b <- NULL

      if (!is.null(data[[SD_P]]) && !is.null(old_data$P)) {
        if (!identical(P@x, old_data$P@x)) new_P <- P
      }
      if (!identical(q, old_data$q)) new_q <- q
      if (!identical(A@x, old_data$A@x)) new_A <- A
      if (!identical(b, old_data$b)) new_b <- b

      ## Send incremental updates
      ## CVXPY v1.8.2 fix: wrap in tryCatch — if sparsity pattern changed,
      ## solver_update() fails; fall back to cold path re-initialization.
      tryCatch({
        if (!is.null(new_P) || !is.null(new_q) ||
            !is.null(new_A) || !is.null(new_b)) {
          clarabel::solver_update(old_solver,
                                  P = new_P, q = new_q,
                                  A = new_A, b = new_b)
        }
        result <- clarabel::solver_solve(old_solver)
        used_warm <- TRUE
      }, error = function(e) {
        ## Sparsity pattern or dimensions changed; cold path will handle it
      })
    }
  }

  ## -- Cold path --------------------------------------------------
  if (!used_warm) {
    ## For warm-start-capable solver, disable features that block updates
    if (warm_start) {
      settings$presolve_enable <- FALSE
      settings$chordal_decomposition_enable <- FALSE
      settings$input_sparse_dropzeros <- FALSE
    }

    ## Use strict_cone_order = FALSE when generalized power cones are present,
    ## since each gp cone needs a unique name (gp1, gp2, ...).
    has_gp <- any(grepl("^gp", names(cones)))
    solver_obj <- clarabel::clarabel_solver(
      A = A, b = b, q = q, P = P, cones = cones, control = settings,
      strict_cone_order = !has_gp
    )
    result <- clarabel::solver_solve(solver_obj)
  }

  ## -- Cache for future warm starts -------------------------------
  if (!is.null(solver_cache)) {
    solver_to_cache <- if (used_warm) old_solver else solver_obj
    assign(cache_key,
           list(solver = solver_to_cache,
                data = list(P = P, q = q, A = A, b = b)),
           envir = solver_cache)
  }

  ## Carry the flag to reduction_invert(), which is where the status is read.
  ## Attached AFTER both the warm and cold paths so neither can miss it, and
  ## after caching so the cached entry stays a pure solver result.
  result[[CLARABEL_ACCEPT_UNKNOWN]] <- accept_unknown

  result
}

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