R/293_reductions_solvers_nlp_solving_chain.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: reductions/solvers/nlp_solving_chain.py
## NLP solving entry point: builds the DNLP reduction chain and drives the solve.
## Mirrors cvxpy nlp_solving_chain.py:27-238.
##
## Chain order (cvxpy nlp_solving_chain.py:56-64):
##   [FlipObjective?] -> CvxAttr2Constr(reduce_bounds = !BOUNDED_VARIABLES)
##                    -> Dnlp2Smooth -> <NLP solver>
##
## Both the single-shot path and best-of-N random restarts (`best_of`, with the
## `_set_random_nlp_initial_point` machinery + `sample_bounds`) are ported --
## cvxpy nlp_solving_chain.py:69-238.


# -- .build_nlp_chain ---------------------------------------------
## CVXPY SOURCE: nlp_solving_chain.py:27-66.
## Returns list(chain = SolvingChain, kwargs = <possibly variant-updated>).
## Solver selection may merge variant kwargs (e.g. uno_ipm -> preset/linear_solver).

.build_nlp_chain <- function(problem, solver, kwargs) {
  if (is.null(solver)) {
    ## Pick the first available NLP solver in preference order.
    solver_instance <- NULL
    for (name in NLP_SOLVER_PREFERENCE) {
      if (.nlp_solver_available(name)) {
        solver_instance <- SOLVER_MAP_NLP[[name]]
        break
      }
    }
    if (is.null(solver_instance)) {
      cli_abort(c(
        "No NLP solver is installed.",
        "i" = 'Install {.pkg Uno} with {.run install.packages("Uno")},',
        "i" = "or {.pkg ipopt} from {.url https://bnaras.github.io/ipopt/}."
      ))
    }
  } else if (exists(solver, envir = SOLVER_MAP_NLP, inherits = FALSE)) {
    ## Base solver name (psolve already upper-cased it; map keys are upper-case).
    solver_instance <- SOLVER_MAP_NLP[[solver]]
  } else if (!is.null(NLP_SOLVER_VARIANTS[[tolower(solver)]])) {
    ## Variant name (e.g. uno_ipm / knitro_sqp): merge its kwargs, use the base.
    variant <- NLP_SOLVER_VARIANTS[[tolower(solver)]]
    kwargs <- utils::modifyList(kwargs, variant$kwargs)
    solver_instance <- SOLVER_MAP_NLP[[variant$base]]
  } else {
    cli_abort("Solver {.val {solver}} is not supported for NLP problems.")
  }

  ## DNLP derivatives are supplied by CVXR's sparsediff-backed diff engine
  ## (sparsediff is independent of the NLP solver -- a user may have Uno/ipopt
  ## but not sparsediff). Guard here so the failure is informative, not a bare
  ## "no package called 'sparsediff'" deep inside the diff engine.
  if (!requireNamespace("sparsediff", quietly = TRUE)) {
    cli_abort(c(
      "DNLP requires the {.pkg sparsediff} package for derivatives.",
      "i" = 'Install it with {.run install.packages("sparsediff")}.'
    ))
  }

  reductions <- if (.s7_is(problem@objective, Maximize)) {
    list(FlipObjective())
  } else {
    list()
  }
  reductions <- c(reductions, list(
    CvxAttr2Constr(reduce_bounds = !solver_instance@BOUNDED_VARIABLES),
    Dnlp2Smooth(),
    solver_instance
  ))

  list(chain = SolvingChain(reductions = reductions), kwargs = kwargs)
}


# -- .set_nlp_initial_point ---------------------------------------
## CVXPY SOURCE: nlp_solving_chain.py:69-93.
## For each variable without a user-set value, construct an initial point from
## get_bounds() (which folds in sign attributes): midpoint if both bounds finite,
## one unit inside a single finite bound, else zero. Stored via save_leaf_value
## (no validation/projection -- matches CVXPY's var.save_value).

.set_nlp_initial_point <- function(problem) {
  for (var in variables(problem)) {
    if (!is.null(value(var))) next

    gb <- get_bounds(var)
    lb <- gb[[1L]]
    ub <- gb[[2L]]

    lb_fin <- is.finite(lb)
    ub_fin <- is.finite(ub)

    init <- numeric(length(lb))
    both    <- lb_fin & ub_fin
    lb_only <- lb_fin & !ub_fin
    ub_only <- !lb_fin & ub_fin
    init[both]    <- 0.5 * (lb[both] + ub[both])
    init[lb_only] <- lb[lb_only] + 1.0
    init[ub_only] <- ub[ub_only] - 1.0

    save_leaf_value(var, matrix(init, nrow = var@shape[1L], ncol = var@shape[2L]))
  }
  invisible(problem)
}


# -- .set_random_nlp_initial_point --------------------------------
## CVXPY SOURCE: nlp_solving_chain.py:96-153.
## Random initial point for best_of restarts. A variable is randomized if its
## sample_bounds are set, OR its value is unset and it has finite variable
## bounds. Variables with a user-set value (and no sample_bounds) are restored
## to that value each run. `user_initials` is an environment keyed by var id;
## a stored NULL is the sentinel "always randomize this variable".

.set_random_nlp_initial_point <- function(problem, run, user_initials) {
  vars <- variables(problem)

  ## Run 0: capture user-specified initial values (variable.py / nlp chain).
  if (run == 0L) {
    rm(list = ls(envir = user_initials, all.names = TRUE), envir = user_initials)
    for (v in vars) {
      id <- as.character(v@id)
      if (!is.null(sample_bounds(v))) {
        assign(id, NULL, envir = user_initials)        # sentinel: always randomize
      } else {
        assign(id, value(v), envir = user_initials)    # capture (possibly NULL)
      }
    }
  }

  for (v in vars) {
    id <- as.character(v@id)
    ui <- if (exists(id, envir = user_initials, inherits = FALSE)) {
      get(id, envir = user_initials, inherits = FALSE)
    } else {
      NULL
    }
    if (!is.null(ui)) {
      save_leaf_value(v, ui)            # restore user value; do not randomize
      next
    }

    ## Randomize within sample_bounds (preferred) or variable bounds (fallback).
    sb <- sample_bounds(v)
    if (is.null(sb)) sb <- get_bounds(v)
    low <- sb[[1L]]
    high <- sb[[2L]]
    if (!all(is.finite(low)) || !all(is.finite(high))) {
      cli_abort(c(
        "Variable {.val {expr_name(v)}} has non-finite sampling bounds.",
        "i" = "Set {.code sample_bounds(var) <- c(low, high)} or finite variable bounds for {.arg best_of}."
      ))
    }
    n <- prod(v@shape)
    ## low/high are scalars (broadcast) or length-n vectors; both recycle.
    init <- low + stats::runif(n) * (high - low)
    save_leaf_value(v, matrix(init, nrow = v@shape[1L], ncol = v@shape[2L]))
  }
  invisible(problem)
}


# -- solve_nlp ----------------------------------------------------
## CVXPY SOURCE: nlp_solving_chain.py:156-238 (single-shot path).
## Applies the DNLP chain, calls the terminal NLP solver, and unpacks back into
## the problem. Returns the optimal value.

solve_nlp <- function(problem, solver = NULL, warm_start = FALSE,
                      verbose = FALSE, ...) {
  kwargs <- list(...)
  built <- .build_nlp_chain(problem, solver, kwargs)
  nlp_chain <- built$chain
  kwargs <- built$kwargs

  has_best_of <- !is.null(kwargs[["best_of"]])

  ## Reuse cached Oracles across solves when the problem has parameters, or for
  ## best_of restarts (the C problem structure is identical across runs -- only
  ## the initial point changes). cvxpy nlp_solving_chain.py:182-185.
  solver_cache <- problem@.cache$nlp_solver_cache
  if (is.null(solver_cache) &&
      (length(parameters(problem)) > 0L || has_best_of)) {
    solver_cache <- new.env(parent = emptyenv())
    problem@.cache$nlp_solver_cache <- solver_cache
  }

  ## -- Single-shot path (cvxpy nlp_solving_chain.py:187-194) ------------
  if (!has_best_of) {
    .set_nlp_initial_point(problem)
    applied <- reduction_apply(nlp_chain, problem)
    data <- applied[[1L]]
    inverse_data <- applied[[2L]]
    raw <- solve_via_data(nlp_chain@solver, data, warm_start, verbose,
                          solver_opts = kwargs, solver_cache = solver_cache)
    problem_unpack_results(problem, raw, nlp_chain, inverse_data)
    return(value(problem))
  }

  ## -- best-of-N restarts (cvxpy nlp_solving_chain.py:196-238) ----------
  best_of <- kwargs[["best_of"]]
  kwargs[["best_of"]] <- NULL
  if (!is.numeric(best_of) || length(best_of) != 1L ||
      best_of < 1 || best_of != as.integer(best_of)) {
    cli_abort("{.arg best_of} must be a positive integer.")
  }
  best_of <- as.integer(best_of)

  best_obj <- Inf
  best_solution <- NULL
  all_objs <- numeric(best_of)
  user_initials <- new.env(parent = emptyenv())
  inverse_data <- NULL

  for (run in seq_len(best_of)) {
    ## CVXPY runs are 0-based; .set_random_nlp_initial_point captures user
    ## values on run 0.
    .set_random_nlp_initial_point(problem, run - 1L, user_initials)
    applied <- reduction_apply(nlp_chain, problem)
    data <- applied[[1L]]
    inverse_data <- applied[[2L]]
    solution <- solve_via_data(nlp_chain@solver, data, warm_start, verbose,
                               solver_opts = kwargs, solver_cache = solver_cache)
    ## Unpack to read the objective in original-problem space
    ## (+Inf for infeasible runs, -Inf for unbounded runs).
    problem_unpack_results(problem, solution, nlp_chain, inverse_data)
    obj_value <- value(problem)
    all_objs[run] <- obj_value

    ## Always keep the first run's solution, so an all-infeasible best_of still
    ## has a solution to unpack at the end (its status then propagates).
    if (is.null(best_solution) || obj_value < best_obj) {
      best_obj <- obj_value
      best_solution <- solution
    }
    if (verbose) {
      cli_inform("Run {run}/{best_of}: obj = {format(obj_value, digits = 6)} | best = {format(best_obj, digits = 6)}")
    }
  }

  ## Report all run objectives to the user in original-problem space.
  if (.s7_is(problem@objective, Maximize)) all_objs <- -all_objs
  best_solution[["all_objs_from_best_of"]] <- all_objs
  problem_unpack_results(problem, best_solution, nlp_chain, inverse_data)
  value(problem)
}