R/265_reductions_solvers_nlp_solvers_ipopt_nlpif.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: reductions/solvers/nlp_solvers/ipopt_nlpif.py
## IPOPT(NLPsolver) -- wraps the lightweight R `ipopt` package, fed by the
## shared diff_engine Oracles built in nlp_solver.R. Mirrors cvxpy
## ipopt_nlpif.py:25-211.


# ==================================================================
# IPOPT status map (cvxpy ipopt_nlpif.py:34-62)
# ==================================================================
IPOPT_STATUS_MAP <- list(
  `0`    = OPTIMAL,
  `1`    = OPTIMAL_INACCURATE,
  `6`    = OPTIMAL,
  `2`    = INFEASIBLE,
  `4`    = UNBOUNDED,
  `3`    = SOLVER_ERROR,
  `-2`   = SOLVER_ERROR,
  `-3`   = SOLVER_ERROR,
  `-13`  = SOLVER_ERROR,
  `-100` = SOLVER_ERROR,
  `-101` = SOLVER_ERROR,
  `-199` = SOLVER_ERROR,
  `5`    = USER_LIMIT,
  `-1`   = USER_LIMIT,
  `-4`   = USER_LIMIT,
  `-5`   = USER_LIMIT,
  `-102` = USER_LIMIT,
  `-10`  = SOLVER_ERROR,
  `-11`  = SOLVER_ERROR,
  `-12`  = SOLVER_ERROR
)


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

method(solver_name, IPOPT_NLP_Solver) <- function(x) IPOPT_SOLVER

method(reduction_invert, IPOPT_NLP_Solver) <- function(x, solution, inverse_data, ...) {
  attr_list <- list()
  attr_list[[RK_NUM_ITERS]] <- solution$num_iters %||% NA_integer_
  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)
  }

  status <- IPOPT_STATUS_MAP[[as.character(solution$status_code)]] %||% SOLVER_ERROR
  if (status %in% SOLUTION_PRESENT) {
    primal_val <- solution$objective
    opt_val <- primal_val + inverse_data@.extra$offset
    primal_vars <- list()
    x_opt <- solution$solution
    for (id in names(inverse_data@var_offsets)) {
      offset <- inverse_data@var_offsets[[id]]
      shape  <- inverse_data@var_shapes[[id]]
      size   <- prod(shape)
      primal_vars[[id]] <- matrix(x_opt[(offset + 1L):(offset + size)],
                                  nrow = shape[1L], ncol = shape[2L])
    }
    ## CVXPY's IPOPT interface returns no dual variables.
    Solution(status = status, opt_val = opt_val,
             primal_vars = primal_vars, dual_vars = list(), attr = attr_list)
  } else {
    failure_solution(status, attr_list)
  }
}

method(solve_via_data, IPOPT_NLP_Solver) <- function(x, data, warm_start = FALSE,
                                                     verbose = FALSE,
                                                     solver_opts = list(), ...) {
  solver_cache <- list(...)[["solver_cache"]]
  if (!requireNamespace("ipopt", quietly = TRUE)) {
    cli_abort(c(
      "NLP solver {.val IPOPT} unavailable: package {.pkg ipopt} is not installed.",
      "i" = "Install it from {.url https://bnaras.github.io/ipopt/}."
    ))
  }

  bounds <- data[["_bounds"]]

  opts <- if (length(solver_opts) > 0L) as.list(solver_opts) else list()
  hessian_approx <- opts[["hessian_approximation"]] %||% "exact"
  use_hessian <- identical(hessian_approx, "exact")

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

  x0 <- data[["x0"]]
  cl <- data[["cl"]]
  js <- oracles$jacobianstructure()
  hs <- oracles$hessianstructure()

  structure_matrix <- function(sp) {
    if (length(sp$rows) == 0L) matrix(integer(), ncol = 2)
    else cbind(as.integer(sp$rows) + 1L, as.integer(sp$cols) + 1L)
  }

  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, obj_factor, lambda) oracles$hessian(u, lambda, obj_factor)

  default_options <- list(
    mu_strategy = "adaptive",
    tol = 1e-7,
    bound_relax_factor = 0.0,
    hessian_approximation = "exact",
    derivative_test = "none",
    least_square_init_duals = "no"
  )
  default_options[names(opts)] <- opts
  if (!verbose && is.null(default_options[["print_level"]])) {
    default_options[["print_level"]] <- 3L
  }

  ## Ipopt may query derivative callbacks before a normal objective/constraint
  ## forward pass, so prime the diff engine exactly as CVXPY does.
  oracles$objective(x0)
  if (length(cl) > 0L) oracles$constraints(x0)

  ans <- ipopt::ipopt_solve(
    x0 = x0,
    lower = data[["lb"]],
    upper = data[["ub"]],
    constraint_lower = cl,
    constraint_upper = data[["cu"]],
    eval_f = obj_cb,
    eval_grad_f = grad_cb,
    eval_g = cons_cb,
    eval_jac_g = jac_cb,
    jacobian_structure = structure_matrix(js),
    eval_h = hess_cb,
    hessian_structure = structure_matrix(hs),
    options = default_options
  )
  ans$num_iters <- NA_integer_
  ans
}