R/263_reductions_solvers_nlp_solvers_nlp_solver.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: reductions/solvers/nlp_solvers/nlp_solver.py
## NLPsolver (base NLP solver) + Bounds (problem -> NLP standard form bounds) +
## Oracles (value/derivative oracle wrapping the diff_engine C_problem).
##
## The runnable backend (Uno) is uno_nlpif.R; the solve dispatch is
## nlp_solving_chain.R. This file provides the shared building blocks.


# ==================================================================
# NLPsolver
# ==================================================================
## CVXPY SOURCE: nlp_solver.py:38-72 (NLPsolver(Solver)).
## Class attrs: REQUIRES_CONSTR=False (default), MIP_CAPABLE=False,
## BOUNDED_VARIABLES=True (the NLP solver consumes variable bounds directly, so
## CvxAttr2Constr(reduce_bounds = not BOUNDED_VARIABLES) leaves them intact).

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

## accepts: only disciplined nonlinear programs (CVXPY nlp_solver.py:46-50)
method(reduction_accepts, NLPsolver) <- function(x, problem, ...) {
  is_dnlp(problem)
}

## apply: build NLP problem data (CVXPY nlp_solver.py:52-72).
##   minimize    f(x)
##   subject to  g^l <= g(x) <= g^u,   x^l <= x <= x^u
method(reduction_apply, NLPsolver) <- function(x, problem, ...) {
  bounds <- Bounds(problem)
  inverse_data <- InverseData(bounds@new_problem)
  inverse_data@.extra$offset <- 0.0
  ## Per-constraint dual metadata, in the lowered-constraint order the solver
  ## sees (= the order of the constraint-dual vector). Lets the solver's invert
  ## map NLP constraint duals back to the original constraints by id. The id is
  ## preserved by lower_equality / lower_ineq_to_nonneg, and remapped upstream
  ## through Dnlp2Smooth / CvxAttr2Constr via their cons_id_maps.
  ## is_eq marks Zero (from equality, dual sign kept) vs NonNeg (from
  ## inequality, dual sign flipped) -- see reduction_invert(UNO_NLP_Solver).
  inverse_data@.extra$nlp_dual_info <- lapply(
    bounds@new_problem@constraints,
    function(con) list(
      id    = as.character(con@id),
      size  = expr_size(con),
      shape = con@args[[1L]]@shape,
      is_eq = .s7_is(con, Zero)
    )
  )
  data <- list(
    problem  = bounds@new_problem,
    cl       = bounds@cl,
    cu       = bounds@cu,
    lb       = bounds@lb,
    ub       = bounds@ub,
    x0       = bounds@x0,
    `_bounds` = bounds   # kept for deferred Oracles creation in solve_via_data
  )
  list(data, inverse_data)
}


# ==================================================================
# Bounds
# ==================================================================
## CVXPY SOURCE: nlp_solver.py:74-148 (Bounds).
## Lowers the problem to NLP standard form: equalities -> Zero (cl=cu=0),
## inequalities -> NonNeg (cl=0, cu=Inf); variable bounds from get_bounds();
## the initial point from variable values.

Bounds <- new_class("Bounds", package = "CVXR",
  properties = list(
    problem     = class_any,
    new_problem = class_any,
    cl          = class_numeric,
    cu          = class_numeric,
    lb          = class_numeric,
    ub          = class_numeric,
    x0          = class_numeric
  ),
  constructor = function(problem) {
    if (FALSE) new_object(S7_object())  ## S7 static-check guard
    main_var <- variables(problem)
    cb <- .nlp_constraint_bounds(problem)
    vb <- .nlp_variable_bounds(main_var)
    x0 <- .nlp_initial_point(main_var)
    .fast_new(Bounds, S7_object(),
      problem     = problem,
      new_problem = cb$new_problem,
      cl          = cb$cl,
      cu          = cb$cu,
      lb          = vb$lb,
      ub          = vb$ub,
      x0          = x0
    )
  }
)

## Constraint bounds + lowered problem (CVXPY get_constraint_bounds).
.nlp_constraint_bounds <- function(problem) {
  lower <- numeric(0)
  upper <- numeric(0)
  new_constr <- list()
  for (con in problem@constraints) {
    sz <- expr_size(con)
    if (.s7_is(con, Equality)) {
      lower <- c(lower, rep(0, sz))
      upper <- c(upper, rep(0, sz))
      new_constr <- c(new_constr, list(lower_equality(con)))
    } else if (.s7_is(con, Inequality)) {
      lower <- c(lower, rep(0, sz))
      upper <- c(upper, rep(Inf, sz))
      new_constr <- c(new_constr, list(lower_ineq_to_nonneg(con)))
    } else {
      cli_abort("NLP Bounds: unsupported constraint type {.cls {class(con)[[1L]]}}.")
    }
  }
  list(new_problem = Problem(problem@objective, new_constr), cl = lower, cu = upper)
}

## Variable bounds via get_bounds() (CVXPY get_variable_bounds). get_bounds()
## returns dense (lb, ub) arrays of length prod(shape), column-major.
.nlp_variable_bounds <- function(main_var) {
  lb <- numeric(0)
  ub <- numeric(0)
  for (v in main_var) {
    gb <- get_bounds(v)
    lb <- c(lb, as.numeric(gb[[1L]]))
    ub <- c(ub, as.numeric(gb[[2L]]))
  }
  list(lb = lb, ub = ub)
}

## Initial point from variable values (CVXPY construct_initial_point).
.nlp_initial_point <- function(main_var) {
  x0 <- numeric(0)
  for (v in main_var) {
    val <- value(v)
    if (is.null(val)) {
      cli_abort(c(
        "Variable {.val {expr_name(v)}} has no value.",
        "i" = "NLP solvers need an initial point; set values on all variables."
      ))
    }
    x0 <- c(x0, as.numeric(val))   # column-major flatten
  }
  x0
}


# ==================================================================
# Oracles
# ==================================================================
## CVXPY SOURCE: nlp_solver.py:150-252 (Oracles).
## Wraps the diff_engine C_problem and exposes value/derivative oracles.
## Returned as an environment (reference semantics, mutable sparsity cache) --
## the same pattern as C_problem itself.

.nlp_oracles <- function(problem, verbose = FALSE, use_hessian = TRUE) {
  cp <- .de_C_problem(problem, verbose = verbose)
  cp$init_jacobian_coo()
  if (use_hessian) {
    cp$init_hessian_coo_lower_tri()
  }
  jac_structure <- NULL
  hess_structure <- NULL

  orc <- new.env(parent = emptyenv())
  orc$c_problem   <- cp
  orc$use_hessian <- use_hessian

  orc$objective   <- function(x) cp$objective_forward(x)
  orc$gradient    <- function(x) cp$gradient()
  orc$constraints <- function(x) cp$constraint_forward(x)
  orc$jacobian    <- function(x) cp$jacobian_values()
  orc$jacobianstructure <- function() {
    if (is.null(jac_structure)) jac_structure <<- cp$jacobian_sparsity()
    jac_structure
  }
  orc$hessian <- function(x, duals, obj_factor) {
    if (!use_hessian) {
      cli_abort("Hessian oracle called but use_hessian is FALSE (bug).")
    }
    cp$hessian_values(obj_factor, duals)
  }
  orc$hessianstructure <- function() {
    if (!use_hessian) return(list(rows = integer(0), cols = integer(0)))
    if (is.null(hess_structure)) hess_structure <<- cp$hessian_sparsity()
    hess_structure
  }
  orc$update_params <- function(problem) {
    params <- parameters(problem)
    if (length(params) == 0L) {
      cli_abort("update_params called but problem has no parameters (bug).")
    }
    theta <- unlist(lapply(params, function(p) as.numeric(value(p))))
    cp$update_params(theta)
  }
  orc
}