R/239_reductions_solvers_conic_solvers_mosek_conif.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: reductions/solvers/conic_solvers/mosek_conif.py
## MOSEK solver interface via Rmosek
##
## MOSEK is a commercial conic solver with native support for LP, QP,
## SOCP, ExpCone, PowCone, and SDP. Uses the standard ConicSolver pipeline
## (ConeMatrixStuffing -> format_constraints -> ConicSolver.reduction_apply).
##
## Key design:
## - Zero constraints -> prob$A + prob$bc (equality bounds, blc = buc = b)
## - NonNeg constraints -> prob$A + prob$bc (inequality bounds, blc = -Inf, buc = b)
## - Conic cones (SOC, PSD, ExpCone, PowCone3D) -> ACC (prob$F + prob$g + prob$cones)
## - PSD via SVEC_PSD_CONE in ACC (not bar variables)
## - SCS's lower-tri column-major with sqrt2 scaling matches MOSEK's SVEC_PSD_CONE
## - ExpCone order: MOSEK PEXP(x,y,z) where x >= y*exp(z/y)
##   vs CVXR (x,y,z) where z >= y*exp(x/y), so EXP_CONE_ORDER = c(2,1,0)
## - QP via prob$qobj lower-tri triplets
## - Rmosek SVEC_PSD_CONE bug workaround: dummy bardim + barc when PSD present


# -- MOSEK status map ----------------------------------------------
## Maps MOSEK solution status strings to CVXR status constants.

MOSEK_STATUS_MAP <- list(
  "OPTIMAL"                  = OPTIMAL,
  "PRIMAL_INFEASIBLE_CER"    = INFEASIBLE,
  "DUAL_INFEASIBLE_CER"      = UNBOUNDED,
  "NEAR_OPTIMAL"             = OPTIMAL_INACCURATE,
  "PRIMAL_ILLPOSED_CER"      = INFEASIBLE,
  "DUAL_ILLPOSED_CER"        = UNBOUNDED,
  "PRIMAL_AND_DUAL_FEASIBLE" = OPTIMAL_INACCURATE,
  "UNKNOWN"                  = SOLVER_ERROR
)

# -- mosek_tri_to_full: expand lower-tri svec to full matrix -------
## Same convention as SCS: lower-tri column-major with sqrt2 off-diag scaling.
## Scales off-diagonal by 1/sqrt(2) to recover the actual matrix entries.
## Returns a numeric vector of length n*n (column-major).

mosek_tri_to_full <- function(lower_tri, n) {
  scs_tri_to_full(lower_tri, n)
}

# -- dims_to_mosek_cones: ConeDims -> prob$cones matrix -------------
## Returns a 3-row matrix with one column per cone block.
## Row 1: type string, Row 2: dimension, Row 3: conepar.
## Also returns the total number of ACC rows (for building F and g).

dims_to_mosek_cones <- function(cone_dims) {
  cols <- list()

  ## NOTE: NonNeg is handled via prob$A/bc bounds, NOT ACC cones.

  ## SOC -> QUAD
  if (length(cone_dims@soc) > 0L) {
    for (q in cone_dims@soc) {
      cols <- c(cols, list(list("QUAD", q, NULL)))
    }
  }

  ## PSD -> SVEC_PSD_CONE
  if (length(cone_dims@psd) > 0L) {
    for (d in cone_dims@psd) {
      svec_dim <- (d * (d + 1L)) %/% 2L
      cols <- c(cols, list(list("SVEC_PSD_CONE", svec_dim, NULL)))
    }
  }

  ## ExpCone -> PEXP (3 elements each)
  if (cone_dims@exp > 0L) {
    for (i in seq_len(cone_dims@exp)) {
      cols <- c(cols, list(list("PEXP", 3L, NULL)))
    }
  }

  ## PowCone3D -> PPOW (3 elements each, with alpha parameters)
  ## MOSEK PPOW expects conepar = c(alpha, 1-alpha) for x1^alpha * x2^(1-alpha) >= |x3|
  ## CVXR stores single alpha in cone_dims@p3d
  if (length(cone_dims@p3d) > 0L) {
    for (alpha in cone_dims@p3d) {
      cols <- c(cols, list(list("PPOW", 3L, c(alpha, 1 - alpha))))
    }
  }

  if (length(cols) == 0L) return(NULL)

  ## Build 3-row matrix
  cones_mat <- matrix(list(), nrow = 3, ncol = length(cols))
  rownames(cones_mat) <- c("type", "dim", "conepar")
  for (j in seq_along(cols)) {
    cones_mat[1, j] <- cols[[j]][[1L]]
    cones_mat[2, j] <- cols[[j]][[2L]]
    cones_mat[3, j] <- list(cols[[j]][[3L]])
  }
  cones_mat
}

# -- Mosek_Solver class --------------------------------------------

#' @keywords internal
Mosek_Solver <- new_class("Mosek_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),
      EXP_CONE_ORDER = c(2L, 1L, 0L),
      REQUIRES_CONSTR = FALSE
    )
  }
)

method(solver_name, Mosek_Solver) <- function(x) MOSEK_SOLVER

## PSD format: reuse SCS's lower-tri column-major with sqrt2 scaling
## (matches MOSEK's SVEC_PSD_CONE convention exactly)
method(solver_psd_format_mat, Mosek_Solver) <- function(solver, constr) {
  scs_psd_format_mat_fn(constr)
}

# -- mosek_extract_dual_value --------------------------------------
## PSD constraints use lower-tri svec format; expand to full.
## Same logic as SCS.

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

# -- MOSEK solve_via_data ------------------------------------------
## Converts CVXR solver data (A, b, c, dims, P) to MOSEK prob list
## and calls Rmosek::mosek().

method(solve_via_data, Mosek_Solver) <- function(x, data, warm_start = FALSE, verbose = FALSE,
                                                  solver_opts = list(), ...) {
  if (!requireNamespace("Rmosek", quietly = TRUE)) {
    cli_abort("Package {.pkg Rmosek} is required but not installed.")
  }

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

  A_full <- data[[SD_A]]
  b_full <- data[[SD_B]]
  c_vec  <- data[[SD_C]]
  dims   <- data[[SD_DIMS]]
  n      <- length(c_vec)  # number of variables

  ## -- Build MOSEK prob list --------------------------------------
  prob <- list()
  prob$sense <- "min"
  prob$c     <- c_vec
  prob$bx    <- rbind(blx = rep(-Inf, n), bux = rep(Inf, n))

  ## -- Row splitting ----------------------------------------------
  ## Solver convention: data$A * x + s = data$b, s in K
  ##
  ## MOSEK row assignment:
  ##   Zero rows     -> prob$A/bc with blc = buc = b  (equality)
  ##   NonNeg rows   -> prob$A/bc with blc = b, buc = Inf (inequality)
  ##   Conic rows    -> ACC: F*x + g in K  where F = -A, g = b
  ##
  ## NonNeg goes to prob$A/bc (not ACC) because MOSEK does not allow

  ## mixing quadratic objectives (qobj) with conic ACC constraints.

  zero_dim   <- dims@zero
  nonneg_dim <- dims@nonneg
  linear_dim <- zero_dim + nonneg_dim

  ## Linear portion (Zero + NonNeg) -> prob$A / prob$bc
  if (linear_dim > 0L) {
    A_linear <- A_full[seq_len(linear_dim), , drop = FALSE]
    b_linear <- b_full[seq_len(linear_dim)]
    prob$A <- A_linear

    ## Bounds: Zero rows get equality, NonNeg rows get one-sided
    ## Zero: A*x + s = b, s = 0  ->  A*x = b  ->  blc = buc = b
    ## NonNeg: A*x + s = b, s >= 0  ->  A*x <= b  ->  blc = -Inf, buc = b
    blc <- rep(-Inf, linear_dim)
    buc <- b_linear
    if (zero_dim > 0L) {
      blc[seq_len(zero_dim)] <- b_linear[seq_len(zero_dim)]
    }
    prob$bc <- rbind(blc = blc, buc = buc)
  } else {
    prob$A  <- Matrix::sparseMatrix(i = integer(0), j = integer(0),
                                    x = numeric(0), dims = c(0L, n))
    prob$bc <- rbind(blc = numeric(0), buc = numeric(0))
  }

  ## -- ACC: conic rows (SOC, PSD, ExpCone, PowCone3D) ------------
  ## ACC: F*x + g in K  where s = b - A*x in K  ->  F = -A, g = b
  total_rows <- nrow(A_full)
  conic_rows <- if (linear_dim < total_rows) (linear_dim + 1L):total_rows else integer(0)

  if (length(conic_rows) > 0L) {
    prob$F <- -A_full[conic_rows, , drop = FALSE]
    prob$g <- b_full[conic_rows]

    ## Build cone specifications (excluding NonNeg which is in prob$A/bc)
    cones_mat <- dims_to_mosek_cones(dims)
    if (!is.null(cones_mat)) {
      prob$cones <- cones_mat
    }
  }

  ## -- SVEC_PSD_CONE bug workaround ------------------------------
  ## Rmosek 11.1.1 BETA crashes if SVEC_PSD_CONE is used but prob$barc
  ## is not defined. Add dummy 1x1 bar variable with zero objective.
  if (length(dims@psd) > 0L) {
    prob$bardim <- c(1L)
    prob$barc <- list(j = 1L, k = 1L, l = 1L, v = 0)
  }

  ## -- QP: quadratic objective via prob$qobj ----------------------
  if (!is.null(data[[SD_P]])) {
    P <- data[[SD_P]]
    ## Extract lower-triangle triplets (MOSEK wants i >= j, 1-based)
    P_tril <- Matrix::tril(P)
    P_tril <- methods::as(P_tril, "TsparseMatrix")
    if (length(P_tril@x) > 0L) {
      prob$qobj <- list(
        i = P_tril@i + 1L,   # 1-based
        j = P_tril@j + 1L,   # 1-based
        v = P_tril@x
      )
    }
  }

  ## -- Solver options ---------------------------------------------
  opts <- list(verbose = if (verbose) 1L else 0L, soldetail = 1L)
  for (opt_name in names(solver_opts)) {
    opts[[opt_name]] <- solver_opts[[opt_name]]
  }

  ## -- Warm-start: pass previous solution as initial point -------
  ## Rmosek supports warm-start via prob$sol: pass the previous
  ## solution structure and MOSEK uses it as initial point.
  ## opts$usesol = TRUE is the Rmosek default.
  ## NOTE: CVXPY does NOT implement MOSEK warm-start -- this is an
  ## R-specific enhancement using native Rmosek capabilities.
  cache_key <- MOSEK_SOLVER
  if (warm_start && !is.null(solver_cache) && exists(cache_key, envir = solver_cache)) {
    cached <- get(cache_key, envir = solver_cache)
    cached_sol <- cached$sol$itr
    if (!is.null(cached_sol) &&
        !is.null(MOSEK_STATUS_MAP[[cached_sol$solsta]]) &&
        MOSEK_STATUS_MAP[[cached_sol$solsta]] %in% SOLUTION_PRESENT &&
        !is.null(cached_sol$xx) && length(cached_sol$xx) == n) {
      prob$sol <- cached$sol
    }
  }

  ## -- Call MOSEK -------------------------------------------------
  result <- Rmosek::mosek(prob, opts)

  ## -- Cache for future warm-starts ------------------------------
  if (!is.null(solver_cache) && result$response$code == 0L &&
      !is.null(result$sol$itr)) {
    assign(cache_key, result, envir = solver_cache)
  }

  result
}

# -- MOSEK reduction_invert ----------------------------------------
## Parses MOSEK-specific result format into a CVXR Solution.

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

  ## Check for MOSEK errors
  if (solution$response$code != 0L) {
    return(failure_solution(SOLVER_ERROR, attr_list))
  }

  ## Extract interior-point solution
  sol <- solution$sol$itr
  if (is.null(sol)) {
    return(failure_solution(SOLVER_ERROR, attr_list))
  }

  ## Map MOSEK status to CVXR status
  solsta <- sol$solsta
  status <- MOSEK_STATUS_MAP[[solsta]]
  if (is.null(status)) status <- SOLVER_ERROR

  if (status %in% SOLUTION_PRESENT) {
    dims <- inverse_data[[SD_DIMS]]
    opt_val <- sol$pobjval + inverse_data[[SD_OFFSET]]

    ## Primal variables from sol$xx
    primal_vars <- list()
    primal_vars[[as.character(inverse_data[[SOLVER_VAR_ID]])]] <- sol$xx

    ## -- Dual variables -------------------------------------------
    ## Rmosek 11.x does NOT return sol$y directly. Instead, linear
    ## constraint duals come from sol$slc (lower bound slack) and
    ## sol$suc (upper bound slack).
    ##
    ## Sign convention: MOSEK's y = slc - suc satisfies A^T y = c,
    ## but Clarabel/SCS use z satisfying A^T z = -c. So z = -(slc - suc)
    ## = suc - slc.  This maps both Zero and NonNeg duals correctly.
    ##
    ## ACC conic duals come from sol$doty (already in CVXR convention).
    ##
    ## CVXR inverse_data splits:
    ##   SOLVER_EQ_CONSTR  -> Zero constraints
    ##   SOLVER_NEQ_CONSTR -> [NonNeg, SOC, PSD, ExpCone, PowCone3D]
    ##
    ## We concatenate nonneg portion of linear_y with sol$doty to form
    ## the full ineq_dual vector.

    zero_dim   <- dims@zero
    nonneg_dim <- dims@nonneg
    linear_dim <- zero_dim + nonneg_dim

    ## Reconstruct dual vector from slc/suc (suc - slc = Clarabel convention)
    linear_y <- if (linear_dim > 0L && !is.null(sol$suc) &&
                    length(sol$suc) >= linear_dim) {
      sol$suc[seq_len(linear_dim)] - sol$slc[seq_len(linear_dim)]
    } else {
      numeric(0)
    }

    ## Equality (zero cone) duals
    if (zero_dim > 0L && length(linear_y) >= zero_dim) {
      eq_dual <- get_dual_values(
        linear_y[seq_len(zero_dim)],
        mosek_extract_dual_value,
        inverse_data[[SOLVER_EQ_CONSTR]]
      )
    } else {
      eq_dual <- list()
    }

    ## NonNeg duals from linear_y
    nonneg_duals <- if (nonneg_dim > 0L && length(linear_y) >= linear_dim) {
      linear_y[(zero_dim + 1L):linear_dim]
    } else {
      numeric(0)
    }

    conic_duals <- if (!is.null(sol$doty) && length(sol$doty) > 0L) {
      doty <- sol$doty

      ## Un-permute ExpCone dual blocks from MOSEK order [2,1,0] back to CVXR [0,1,2]
      if (dims@exp > 0L) {
        ## Calculate offset to ExpCone blocks within doty
        ## ACC order (no NonNeg): SOC, PSD, ExpCone, PowCone3D
        exp_offset <- sum(dims@soc)
        if (length(dims@psd) > 0L) {
          exp_offset <- exp_offset + sum(vapply(dims@psd, function(d) (d * (d + 1L)) %/% 2L, integer(1L)))
        }
        for (k in seq_len(dims@exp)) {
          base_idx <- exp_offset + (k - 1L) * 3L
          ## Swap indices 1 and 3 (MOSEK [z,y,x] -> CVXR [x,y,z])
          tmp <- doty[base_idx + 1L]
          doty[base_idx + 1L] <- doty[base_idx + 3L]
          doty[base_idx + 3L] <- tmp
        }
      }
      doty
    } else {
      numeric(0)
    }

    ## Concatenate for get_dual_values (order matches SOLVER_NEQ_CONSTR)
    full_ineq_vec <- c(nonneg_duals, conic_duals)
    if (length(full_ineq_vec) > 0L && length(inverse_data[[SOLVER_NEQ_CONSTR]]) > 0L) {
      ineq_dual <- get_dual_values(
        full_ineq_vec,
        mosek_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)
  }
}

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