R/264_reductions_solvers_nlp_solvers_uno_nlpif.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: reductions/solvers/nlp_solvers/uno_nlpif.py
## UNO(NLPsolver) -- the runnable NLP backend, wrapping the Uno solver via the
## `Uno` R package (`uno_solve()`), fed by the diff_engine `Oracles` built in
## nlp_solver.R. Mirrors cvxpy uno_nlpif.py:25-301.
##
## Status convention (R deviation, lossless): the `Uno` package returns
## `optimization_status` / `solution_status` as NAMED integers, e.g.
## `c(SUCCESS = 0L)` -- the integer is Uno's enum code and the name is its
## canonical label. CVXPY reads `str(result.optimization_status)` (the label);
## here the label is `names(status)`, so STATUS_MAP is keyed by that label,
## byte-for-byte the same keys as cvxpy uno_nlpif.py:36-56.


# ==================================================================
# UNO status map (cvxpy uno_nlpif.py:36-56)
# ==================================================================
## Keyed by Uno's canonical status label (names() of the returned named int).
UNO_STATUS_MAP <- list(
  ## Success cases (optimization_status)
  SUCCESS                     = OPTIMAL,
  ## Limit cases
  ITERATION_LIMIT             = USER_LIMIT,
  TIME_LIMIT                  = USER_LIMIT,
  ## Error cases
  EVALUATION_ERROR            = SOLVER_ERROR,
  ALGORITHMIC_ERROR           = SOLVER_ERROR,
  ## Solution status cases
  FEASIBLE_KKT_POINT          = OPTIMAL,
  FEASIBLE_FJ_POINT           = OPTIMAL_INACCURATE,
  FEASIBLE_SMALL_STEP         = OPTIMAL_INACCURATE,
  INFEASIBLE_STATIONARY_POINT = INFEASIBLE,
  INFEASIBLE_SMALL_STEP       = INFEASIBLE,
  UNBOUNDED                   = UNBOUNDED,
  NOT_OPTIMAL                 = SOLVER_ERROR
)


# ==================================================================
# UNO class
# ==================================================================
## CVXPY SOURCE: uno_nlpif.py:25-34 (class UNO(NLPsolver)).

UNO_NLP_Solver <- new_class("UNO_NLP_Solver", parent = NLPsolver, package = "CVXR",
  constructor = function() {
    if (FALSE) new_object(S7_object())  ## S7 static-check guard
    .fast_new(UNO_NLP_Solver, S7_object(),
      .cache = new.env(parent = emptyenv()),
      MIP_CAPABLE = FALSE,
      BOUNDED_VARIABLES = TRUE,
      PSD_TRIANGLE_KIND = NA_character_,
      PSD_SQRT2_SCALING = NA
    )
  }
)

## name (cvxpy uno_nlpif.py:58-62)
method(solver_name, UNO_NLP_Solver) <- function(x) UNO_SOLVER


# ==================================================================
# invert (cvxpy uno_nlpif.py:70-96)
# ==================================================================
## Maps the raw Uno result back to a CVXR Solution. Status handling mirrors
## CVXPY (prefer optimization_status, fall back to solution_status); on a
## present solution, x is reshaped column-major (order F) per
## var_offsets+var_shapes.
##
## CVXR EXTENSION (beyond CVXPY, which returns an empty dual map at
## uno_nlpif.py:94): recover constraint duals from Uno's `constraint_dual`.
## The vector is in lowered-constraint order (nlp_dual_info, from NLPsolver.apply).
## Sign convention, verified against the conic path and KKT stationarity:
##   Zero   (from an equality):   CVXR dual = +constraint_dual
##   NonNeg (from an inequality): CVXR dual = -constraint_dual
## Keyed by the (id-preserving) lowered-constraint id; the chain inverts
## (Dnlp2Smooth / CvxAttr2Constr) remap those ids back to the user constraints.

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

  ## Status: prefer optimization_status, else solution_status (uno_nlpif.py:78).
  ## The value is a named int; its name() is the canonical label STATUS_MAP keys on.
  status_int <- solution$optimization_status
  if (is.null(status_int)) status_int <- solution$solution_status
  status_key <- names(status_int)
  status <- if (!is.null(status_key)) UNO_STATUS_MAP[[status_key]] else NULL
  if (is.null(status)) status <- SOLVER_ERROR

  attr_list[[RK_NUM_ITERS]] <- solution$iterations %||% 0L
  if (!is.null(solution$cpu_time)) attr_list[[RK_SOLVE_TIME]] <- solution$cpu_time

  ## best_of surfaces the per-run objectives via solver_stats extra_stats.
  ## CVXPY wires this in ipopt_nlpif.py:89-91; UNO mirrors that pattern.
  if (!is.null(solution$all_objs_from_best_of)) {
    attr_list[[RK_EXTRA_STATS]] <-
      list(all_objs_from_best_of = solution$all_objs_from_best_of)
  }

  if (status %in% SOLUTION_PRESENT) {
    primal_val <- solution$objective
    opt_val <- primal_val + inverse_data@.extra$offset
    primal_vars <- list()
    x_opt <- solution$primal
    for (id in names(inverse_data@var_offsets)) {
      offset <- inverse_data@var_offsets[[id]]   # 0-based
      shape  <- inverse_data@var_shapes[[id]]
      size   <- prod(shape)
      ## column-major (numpy order='F') reshape
      primal_vars[[id]] <- matrix(x_opt[(offset + 1L):(offset + size)],
                                  nrow = shape[1L], ncol = shape[2L])
    }

    ## Constraint duals: slice constraint_dual into per-constraint blocks (in
    ## lowered order), apply the Zero/NonNeg sign convention, key by id.
    dual_vars <- list()
    cd <- solution$constraint_dual
    dinfo <- inverse_data@.extra$nlp_dual_info
    if (!is.null(cd) && !is.null(dinfo)) {
      pos <- 0L
      for (info in dinfo) {
        sz <- info$size
        block <- cd[(pos + 1L):(pos + sz)]
        pos <- pos + sz
        dval <- if (isTRUE(info$is_eq)) block else -block
        dual_vars[[info$id]] <- matrix(dval, nrow = info$shape[1L],
                                       ncol = info$shape[2L])
      }
    }

    Solution(status = status, opt_val = opt_val,
             primal_vars = primal_vars, dual_vars = dual_vars, attr = attr_list)
  } else {
    failure_solution(status, attr_list)
  }
}


# ==================================================================
# solve_via_data (cvxpy uno_nlpif.py:98-291)
# ==================================================================
## Builds the diff_engine Oracles (deferred from apply so we have `verbose`),
## then drives Uno via `uno_solve()`. The oracle convention is
## sigma*grad^2 f + sum_i lambda_i grad^2 g_i (positive Lagrangian sign), so we
## MUST pass lagrangian_sign = "positive" -- the wrong sign silently corrupts
## the Hessian's constraint terms for NONLINEAR constraints (uno_nlpif.py:237-240).
##
## `data` holds: "_bounds" (the Bounds object), "x0", "lb", "ub", "cl", "cu".

method(solve_via_data, UNO_NLP_Solver) <- function(x, data, warm_start = FALSE,
                                                   verbose = FALSE,
                                                   solver_opts = list(), ...) {
  solver_cache <- list(...)[["solver_cache"]]
  if (!requireNamespace("Uno", quietly = TRUE)) {
    cli_abort(c(
      "NLP solver {.val UNO} unavailable: package {.pkg Uno} is not installed.",
      "i" = 'Install it with {.run install.packages("Uno")}.'
    ))
  }

  bounds <- data[["_bounds"]]
  use_hessian <- TRUE  # Uno uses the exact Hessian (L-BFGS only when none given)

  ## Oracle creation, with optional reuse across solves (cvxpy uno_nlpif.py:138-146).
  ## solver_cache (when supplied) is an environment.
  if (is.null(solver_cache)) {
    oracles <- .nlp_oracles(bounds@new_problem, verbose = verbose,
                            use_hessian = use_hessian)
  } else if (exists("oracles", envir = solver_cache, inherits = FALSE)) {
    oracles <- get("oracles", envir = solver_cache, inherits = FALSE)
    if (length(parameters(bounds@new_problem)) > 0L) {
      oracles$update_params(bounds@new_problem)
    }
  } else {
    oracles <- .nlp_oracles(bounds@new_problem, verbose = verbose,
                            use_hessian = use_hessian)
    assign("oracles", oracles, envir = solver_cache)
  }

  ## Standard-form data (cvxpy uno_nlpif.py:149-156).
  x0 <- data[["x0"]]
  lb <- data[["lb"]]
  ub <- data[["ub"]]
  cl <- data[["cl"]]
  cu <- data[["cu"]]
  n <- length(x0)
  m <- length(cl)

  ## COO sparsity (0-based, the engine convention -- uno_solve base_indexing = 0).
  js <- oracles$jacobianstructure()
  hs <- oracles$hessianstructure()

  ## Callbacks. Note the Hessian arg order: uno_solve's hess(x, sigma, lambda)
  ## maps to oracles$hessian(x, duals = lambda, obj_factor = sigma).
  obj_cb  <- function(u) oracles$objective(u)
  grad_cb <- function(u) oracles$gradient(u)
  cons_cb <- function(u) oracles$constraints(u)
  jac_cb  <- function(u) oracles$jacobian(u)
  hess_cb <- function(u, sigma, lambda) oracles$hessian(u, lambda, sigma)

  ## solver_opts: pop "preset", the rest are forwarded as named Uno options.
  ## CVXR DEVIATES from CVXPY's "filtersqp" default (uno_nlpif.py:252): CVXPY's
  ## filtersqp uses BQPD for its QP subproblem (handles indefinite Hessians), but
  ## this R build is HiGHS-only -- BQPD's Fortran source is not bundled and its
  ## license is not CRAN-compatible. HiGHS solves only CONVEX QPs, so filtersqp
  ## aborts ("Algorithmic error") on the indefinite QP subproblems that nonconvex
  ## DNLPs generate. We therefore default to the interior-point "ipopt" preset,
  ## which factors the regularized KKT system with MUMPS (symmetric-indefinite)
  ## and is robust on convex AND nonconvex problems. Force SQP explicitly with
  ## preset = "filtersqp" (or solver = "uno_sqp").
  opts <- if (length(solver_opts) > 0L) as.list(solver_opts) else list()
  preset <- opts[["preset"]] %||% "ipopt"
  opts[["preset"]] <- NULL
  ## The "ipopt" preset needs a symmetric-indefinite linear solver; default it to
  ## MUMPS (built from rmumps) unless the user specified one.
  if (identical(preset, "ipopt") && is.null(opts[["linear_solver"]])) {
    opts[["linear_solver"]] <- "MUMPS"
  }

  Uno::uno_solve(
    n = n, lb = lb, ub = ub, sense = "minimize",
    obj = obj_cb, grad = grad_cb,
    m = m, cl = cl, cu = cu, cons = cons_cb,
    jac_rows = js$rows, jac_cols = js$cols, jac = jac_cb,
    hess_rows = hs$rows, hess_cols = hs$cols, hess = hess_cb,
    x0 = x0, preset = preset, base_indexing = 0L, verbose = verbose,
    options = opts, lagrangian_sign = "positive"
  )
}