R/238_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
## Clarabel solver interface
##
## 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
)

# -- triu_to_full (Clarabel): expand upper-tri svec to full matrix -
## CVXPY SOURCE: clarabel_conif.py lines 65-96
## Upper triangle version. Scales off-diagonal by 1/sqrt(2).

clarabel_triu_to_full <- function(upper_tri, n) {
  full <- matrix(0, n, n)
  ## Clarabel uses upper-tri in a way that maps to numpy tril_indices
  ## (confusing but correct per CVXPY comment)
  full[lower.tri(full, diag = TRUE)] <- upper_tri
  full <- full + t(full)
  diag(full) <- diag(full) / 2
  ## Scale off-diagonal by 1/sqrt(2)
  off_lower <- lower.tri(full)
  off_upper <- upper.tri(full)
  full[off_lower] <- full[off_lower] / sqrt(2)
  full[off_upper] <- full[off_upper] / sqrt(2)
  as.vector(full)  # column-major
}

# -- Clarabel psd_format_mat --------------------------------------
## CVXPY SOURCE: clarabel_conif.py lines 190-225
## Upper-triangular extraction with sqrt(2) off-diagonal scaling
## + symmetrization.

clarabel_psd_format_mat_fn <- function(constr) {
  n <- constr@args[[1L]]@shape[1L]
  entries <- (n * (n + 1L)) %/% 2L

  ## Row indices: 0 to entries-1
  row_arr <- seq(0L, length.out = entries)

  ## Column indices: upper triangle positions as column-major flat indices
  idx <- which(upper.tri(matrix(0, n, n), diag = TRUE), arr.ind = TRUE)
  ## Convert to 0-based column-major flat indices (already sorted)
  col_arr <- (idx[, 2L] - 1L) * n + (idx[, 1L] - 1L)

  ## Value array: sqrt(2) for off-diagonal, 1.0 for diagonal
  val_mat <- matrix(0, n, n)
  val_mat[upper.tri(val_mat, diag = TRUE)] <- sqrt(2)
  diag(val_mat) <- 1.0
  val_arr <- as.vector(val_mat)  # column-major
  val_arr <- val_arr[val_arr != 0]  # nonzero entries

  scaled_upper_tri <- Matrix::sparseMatrix(
    i = row_arr + 1L, j = col_arr + 1L, x = val_arr,
    dims = c(entries, as.integer(n * n))
  )

  ## Symmetrization matrix: (M + M^T) / 2
  nn <- as.integer(n * n)
  K <- matrix(seq_len(nn) - 1L, nrow = n, ncol = n, byrow = TRUE)  # C-order
  row_symm <- c(seq_len(nn) - 1L, as.vector(K)) + 1L      # 1-based
  col_symm <- c(seq_len(nn) - 1L, as.vector(t(K))) + 1L   # 1-based
  val_symm <- rep(0.5, 2L * nn)

  symm_matrix <- Matrix::sparseMatrix(
    i = row_symm, j = col_symm, x = val_symm,
    dims = c(nn, nn)
  )

  scaled_upper_tri %*% symm_matrix
}

# -- 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() {
    new_object(S7_object(),
      .cache = new.env(parent = emptyenv()),
      MIP_CAPABLE = FALSE,
      BOUNDED_VARIABLES = FALSE,
      SUPPORTED_CONSTRAINTS = list(Zero, NonNeg, SOC, ExpCone,
                                   PowCone3D, PSD, PowConeND),
      EXP_CONE_ORDER = c(0L, 1L, 2L),
      REQUIRES_CONSTR = FALSE
    )
  }
)

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

## Override PSD format for Clarabel (upper-triangular + sqrt(2))
method(solver_psd_format_mat, Clarabel_Solver) <- function(solver, constr) {
  clarabel_psd_format_mat_fn(constr)
}

# -- Clarabel extract_dual_value ----------------------------------
## CVXPY SOURCE: clarabel_conif.py lines 227-243

clarabel_extract_dual_value <- function(result_vec, offset, constraint) {
  if (S7_inherits(constraint, PSD)) {
    dim <- constraint@args[[1L]]@shape[1L]
    upper_tri_dim <- (dim * (dim + 1L)) %/% 2L
    new_offset <- offset + upper_tri_dim
    upper_tri <- result_vec[(offset + 1L):(offset + upper_tri_dim)]
    full <- clarabel_triu_to_full(upper_tri, dim)
    list(value = full, new_offset = new_offset)
  } else {
    extract_dual_value(result_vec, offset, constraint)
  }
}

# -- 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"
  }
  status <- CLARABEL_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

  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

    ## Dual variables: split at zero cone boundary
    z <- solution$z
    zero_dim <- inverse_data[[SD_DIMS]]@zero
    if (zero_dim > 0L) {
      eq_dual <- get_dual_values(
        z[seq_len(zero_dim)],
        clarabel_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)],
        clarabel_extract_dual_value,
        inverse_data[[SOLVER_NEQ_CONSTR]]
      )
    } else {
      ineq_dual <- list()
    }
    dual_vars <- c(eq_dual, ineq_dual)
    Solution(status, opt_val, primal_vars, dual_vars, attr_list)
  } else {
    failure_solution(status, attr_list)
  }
}

# -- 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(), ...) {
  if (!requireNamespace("clarabel", quietly = TRUE)) {
    cli_abort("Package {.pkg clarabel} is required but not installed.")
  }

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

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

  result
}

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