R/294_problems_problem_form.R

In CVXR: Disciplined Convex Optimization

#####
## DO NOT EDIT THIS FILE!! EDIT THE SOURCE INSTEAD: rsrc_tree/problems/problem_form.R
#####

## CVXPY SOURCE: problems/problem_form.py
## PARTIAL PORT: in: _objective_cone_atoms, _expr_cone_atoms, ProblemForm
##   (has_quadratic_objective, cones, _compute_cones, _compute_gp_cones,
##   is_mixed_integer, is_lp, has_constraints), pick_default_solver,
##   make_problem_form.
##   out: the SuppFuncAtom branch of _compute_cones (problem_form.py:230-243),
##   because suppfunc is not implemented (atoms/suppfunc.R); the
##   build_non_disciplined_error_msg detail in make_problem_form, because CVXR
##   keeps its own DCP/DGP messages; APPROX_CONE_CONVERSIONS keys
##   (problem_form.py:265-268), because CVXR has no such table
##   (reductions/solvers/solver.R, expand_cones).
##
## ProblemForm is the structural analysis of a problem BEFORE canonicalization:
## which cones it will need, whether the objective is quadratic, whether it is
## mixed-integer. Solver selection (`.resolve_solver`, solving_chain.R) asks it
## questions, so the solver is chosen first and the QP-or-conic interface
## follows from the solver, as in CVXPY 1.9 (ADR D_19.8, superseding the
## D_19.1 cluster-1 deferral).
##
## Atom types are compared by EXACT class, as upstream's `atom in SOC_ATOMS`
## compares types: `class(x)[1L]` of an S7 object is its fully qualified name
## ("CVXR::PnormApprox"). Upstream's `ptp` (NONPOS_ATOMS) and `std` (SOC_ATOMS)
## are functions in CVXPY 1.9.3, not classes, so those entries never match
## there and are omitted here.

# -- atom-type sets ------------------------------------------------
## CVXPY SOURCE: atoms/__init__.py:120-192. Built on first use, since the
## class objects must exist by then.

.pf_sets <- new.env(hash = TRUE, parent = emptyenv())

.pf_fqn <- function(cls) paste0("CVXR::", attr(cls, "name", exact = TRUE))

.pf_atom_sets <- function() {
  if (isTRUE(.pf_sets$built)) return(.pf_sets)
  set <- function(...) vapply(list(...), .pf_fqn, character(1L))
  .pf_sets$SOC_ATOMS <- set(GeoMeanApprox, PnormApprox, QuadForm, QuadOverLin,
                            PowerApprox, Huber)
  .pf_sets$POWCONE_ATOMS <- set(Pnorm, Power)
  .pf_sets$POWCONE_ND_ATOMS <- set(GeoMean)
  ## von_neumann_entr is NOT TARGETED (CLAUDE.md "Remaining Work").
  .pf_sets$EXP_ATOMS <- set(LogSumExp, LogDet, Entr, Exp, KlDiv, RelEntr, Log,
                            Log1p, Logistic, Xexp)
  ## quantum_rel_entr and von_neumann_entr are NOT TARGETED.
  .pf_sets$PSD_ATOMS <- set(LambdaMax, LambdaSumLargest, ConditionNumber,
                            LogDet, MatrixFrac, NormNuc, SigmaMax, TrInv)
  .pf_sets$NONPOS_ATOMS <- set(Norm1, Abs, Huber)
  .pf_sets$GP_EXP_ATOMS <- set(AddExpression, Exp, EyeMinusInv, Log,
                               MulExpression, Norm1, OneMInusPos, Pnorm,
                               PnormApprox, QuadForm, SumEntries, Trace, Xexp)
  .pf_sets$GP_NONPOS_ATOMS <- set(NormInf, MaxEntries, MinEntries, Maximum,
                                  Minimum, PfEigenvalue)
  .pf_sets$built <- TRUE
  .pf_sets
}

## The children an atom walk visits. CVXPY's PartialProblem has
## `args = [prob]`, so upstream's walkers reach the inner problem; CVXR's
## PartialProblem is an argument-less node with the inner problem in
## `@inner_problem` (transforms/partial_optimize.R), so descend explicitly.
## Missing this predicted no cones, and HIGHS returned "infeasible" on a
## feasible SOC problem (9202). Only zero-args nodes, i.e. leaves, pay the one
## string compare. STOPGAP: revisit if PartialProblem gets real args (option
## A, deferred; ledger divergence + handoff 2026-10-07 record it).
.pf_children <- function(expr) {
  kids <- expr@args
  if (length(kids) == 0L && class(expr)[1L] == "CVXR::PartialProblem") {
    inner <- expr@inner_problem
    kids <- c(list(inner@objective@args[[1L]]), inner@constraints)
  }
  kids
}

## CVXPY SOURCE: atom.py:567-573, constraint.py (`atoms`) -- the class of
## every node in the trees, leaves included; membership in the atom sets
## above keeps only atoms, so leaves and constraints never count.
## Level by level: each level is built in one unlist(), so the order is the
## breadth-first order of a single queue over all roots, without growing a
## queue (and a name vector) one element at a time, which reallocated
## repeatedly (+1.2 MB and 15 ms per solve at 1000 constraints, 9212).
.pf_atoms <- function(roots) {
  level <- roots
  out <- character(0)
  while (length(level) > 0L) {
    out <- c(out, unique(vapply(level, function(e) class(e)[1L], character(1L),
                                USE.NAMES = FALSE)))
    level <- unlist(lapply(level, .pf_children), recursive = FALSE,
                    use.names = FALSE)
  }
  unique(out)
}

# -- _objective_cone_atoms -----------------------------------------
## CVXPY SOURCE: problem_form.py:61-107
## Atom types that still need conic canonicalization under quad_obj = TRUE.
## Mirrors the Dcp2Cone walk: an atom with a quadratic canonicalization in the
## objective's affine head goes down the QP path and produces no cone.
## "type in CANON_METHODS" is `has_dcp_canon` and "type in
## QUAD_CANON_METHODS" is `has_quad_canon`, the predicates Dcp2Cone itself
## dispatches on (reductions/dcp2cone/dcp2cone.R).
.pf_objective_cone_atoms <- function(expr, affine_above = TRUE,
                                     eval_params = FALSE) {
  if (.s7_is(expr, Leaf)) return(character(0))
  ## Constraints (e.g. Equality) may appear as args of atoms like indicator.
  if (.s7_is(expr, Constraint)) {
    result <- character(0)
    for (arg in expr@args) {
      result <- c(result, .pf_objective_cone_atoms(arg, FALSE, eval_params))
    }
    return(result)
  }
  if (is_constant(expr) && (eval_params || length(parameters(expr)) == 0L)) {
    return(character(0))
  }
  canon <- has_dcp_canon(expr)
  child_affine_above <- !canon && affine_above
  result <- character(0)
  for (arg in .pf_children(expr)) {
    result <- c(result,
                .pf_objective_cone_atoms(arg, child_affine_above, eval_params))
  }
  key <- class(expr)[1L]
  if (affine_above && has_quad_canon(expr)) {
    if (.s7_is(expr, Power) && !.quadratic_power(expr)) {
      ## Power with a non-quadratic exponent falls back to cone canon.
      result <- c(result, key)
    } else if (key == .pf_fqn(QuadOverLin) && !is_constant(expr@args[[2L]])) {
      ## quad_over_lin needs a constant denominator for the QP path.
      result <- c(result, key)
    }
    ## Otherwise the QP path handles it: no cone needed.
  } else if (canon) {
    result <- c(result, key)
  }
  result
}

# -- _expr_cone_atoms ----------------------------------------------
## CVXPY SOURCE: problem_form.py:110-136
## Atom types needing conic canonicalization, without the QP-path filtering;
## parameter-only subtrees count as constant when `eval_params` is TRUE.
.pf_expr_cone_atoms <- function(expr, eval_params = FALSE) {
  if (.s7_is(expr, Leaf)) return(character(0))
  if (.s7_is(expr, Constraint)) {
    result <- character(0)
    for (arg in expr@args) result <- c(result, .pf_expr_cone_atoms(arg, eval_params))
    return(result)
  }
  if (is_constant(expr) && (eval_params || length(parameters(expr)) == 0L)) {
    return(character(0))
  }
  result <- character(0)
  for (arg in .pf_children(expr)) {
    result <- c(result, .pf_expr_cone_atoms(arg, eval_params))
  }
  if (has_dcp_canon(expr)) result <- c(result, class(expr)[1L])
  result
}

# -- ProblemForm ---------------------------------------------------
## CVXPY SOURCE: problem_form.py:139-330
## Lazily computed, cached structural properties of a problem.

ProblemForm <- new_class("ProblemForm", package = "CVXR",
  properties = list(
    problem     = class_any,
    gp          = class_logical,
    eval_params = class_logical,
    .cache      = class_environment
  ),
  constructor = function(problem, gp = FALSE, eval_params = FALSE) {
    if (FALSE) new_object(S7_object())  ## S7 static-check guard
    .fast_new(ProblemForm, S7_object(),
      problem     = problem,
      gp          = gp,
      eval_params = eval_params,
      .cache      = new.env(parent = emptyenv())
    )
  }
)

## CVXPY SOURCE: problem_form.py:157-163
.pf_has_quadratic_objective <- function(pf) {
  ans <- pf@.cache$has_quadratic_objective
  if (is.null(ans)) {
    ans <- has_quadratic_term(pf@problem@objective@args[[1L]])
    pf@.cache$has_quadratic_objective <- ans
  }
  ans
}

## CVXPY SOURCE: problem_form.py:165-179
## The cone constraint types required after canonicalization, as a list of
## S7 class objects. `quad_obj = TRUE` drops the objective atoms the QP path
## handles; FALSE gives the conservative full set.
.pf_cones <- function(pf, quad_obj = FALSE) {
  if (is.null(pf@.cache$cones_quad)) .pf_compute_cones(pf)
  if (isTRUE(quad_obj)) pf@.cache$cones_quad else pf@.cache$cones_full
}

## Is `cls` in a list of cone classes? By name, as `.cone_has` (solver.R).
.pf_has_cone <- function(cones, cls) .cone_has(cones, cls)

## CVXPY SOURCE: problem_form.py:181-282
.pf_compute_cones <- function(pf) {
  if (isTRUE(pf@gp)) return(.pf_compute_gp_cones(pf))
  problem <- pf@problem
  eval_params <- isTRUE(pf@eval_params)
  sets <- .pf_atom_sets()

  ## With eval_params, constraints with no variables become trivial once
  ## EvalParams runs, so they are left out of the analysis.
  relevant <- if (eval_params) {
    Filter(function(c) length(variables(c)) > 0L, problem@constraints)
  } else {
    problem@constraints
  }
  constr_types <- unique(vapply(relevant, function(c) class(c)[1L],
                                character(1L), USE.NAMES = FALSE))

  constr_atoms <- if (eval_params) {
    unlist(lapply(relevant, .pf_expr_cone_atoms, eval_params = TRUE),
           use.names = FALSE)
  } else {
    .pf_atoms(relevant)
  }

  obj_expr <- problem@objective@args[[1L]]
  base_obj_atoms <- if (.pf_has_quadratic_objective(pf)) {
    .pf_objective_cone_atoms(obj_expr, eval_params = eval_params)
  } else if (eval_params) {
    .pf_expr_cone_atoms(obj_expr, eval_params = TRUE)
  } else {
    .pf_atoms(list(obj_expr))
  }
  atoms <- unique(c(base_obj_atoms, constr_atoms))

  has_type <- function(cls) .pf_fqn(cls) %in% constr_types
  cones <- list()
  if (has_type(SOC) || any(atoms %in% sets$SOC_ATOMS))
    cones[[length(cones) + 1L]] <- SOC
  if (has_type(ExpCone) || any(atoms %in% sets$EXP_ATOMS))
    cones[[length(cones) + 1L]] <- ExpCone
  if (has_type(Inequality) || has_type(NonPos) || has_type(NonNeg) ||
      any(atoms %in% sets$NONPOS_ATOMS))
    cones[[length(cones) + 1L]] <- NonNeg
  if (has_type(NonPos))
    cones[[length(cones) + 1L]] <- NonPos
  if (has_type(Equality) || has_type(Zero))
    cones[[length(cones) + 1L]] <- Zero
  if (has_type(PSD) || any(atoms %in% sets$PSD_ATOMS) ||
      any(vapply(variables(problem), function(v) is_psd(v) || is_nsd(v),
                 logical(1L))))
    cones[[length(cones) + 1L]] <- PSD
  if (has_type(PowCone3D) || any(atoms %in% sets$POWCONE_ATOMS))
    cones[[length(cones) + 1L]] <- PowCone3D
  if (has_type(PowConeND) || any(atoms %in% sets$POWCONE_ND_ATOMS))
    cones[[length(cones) + 1L]] <- PowConeND

  ## The QP-filtered set is the base.
  pf@.cache$cones_quad <- cones
  pf@.cache$cones_full <- cones
  ## Full set: every QUAD_CANON_METHODS atom canonicalizes to SOC, so SOC is
  ## the only cone the unfiltered objective can add.
  if (!.pf_has_cone(cones, SOC) && .pf_has_quadratic_objective(pf)) {
    if (any(.pf_atoms(list(obj_expr)) %in% sets$SOC_ATOMS)) {
      pf@.cache$cones_full <- c(cones, list(SOC))
    }
  }
  invisible(NULL)
}

## CVXPY SOURCE: problem_form.py:284-314
## DGP atoms map to other DCP atoms than their DCP canonicalization does
## (norm1 needs NonNeg in DCP, ExpCone in GP). No QP path for GP.
.pf_compute_gp_cones <- function(pf) {
  problem <- pf@problem
  sets <- .pf_atom_sets()
  constr_types <- unique(vapply(problem@constraints, function(c) class(c)[1L],
                                character(1L), USE.NAMES = FALSE))
  has_type <- function(cls) .pf_fqn(cls) %in% constr_types
  atoms <- .pf_atoms(c(list(problem@objective@args[[1L]]), problem@constraints))
  cones <- list()
  if (any(atoms %in% sets$GP_EXP_ATOMS))
    cones[[length(cones) + 1L]] <- ExpCone
  if (has_type(Inequality) || has_type(NonPos) || has_type(NonNeg) ||
      any(atoms %in% sets$GP_NONPOS_ATOMS))
    cones[[length(cones) + 1L]] <- NonNeg
  if (has_type(Equality) || has_type(Zero))
    cones[[length(cones) + 1L]] <- Zero
  if (has_type(PSD) ||
      any(vapply(variables(problem), function(v) is_psd(v) || is_nsd(v),
                 logical(1L))))
    cones[[length(cones) + 1L]] <- PSD
  pf@.cache$cones_quad <- cones
  pf@.cache$cones_full <- cones
  invisible(NULL)
}

## CVXPY SOURCE: problem_form.py:316-320
.pf_is_mixed_integer <- function(pf) {
  ans <- pf@.cache$is_mixed_integer
  if (is.null(ans)) {
    ans <- is_mixed_integer(pf@problem)
    pf@.cache$is_mixed_integer <- ans
  }
  ans
}

## CVXPY SOURCE: problem_form.py:333 (`_QP_CONES`)
.pf_within_qp_cones <- function(cones) {
  all(vapply(cones, function(cone) .cone_name(cone) %in% c("Zero", "NonNeg"),
             logical(1L)))
}

## CVXPY SOURCE: problem_form.py:322-324
.pf_is_lp <- function(pf) {
  .pf_within_qp_cones(.pf_cones(pf)) && !.pf_has_quadratic_objective(pf)
}

## CVXPY SOURCE: problem_form.py:326-330
.pf_has_constraints <- function(pf) {
  problem <- pf@problem
  length(.pf_cones(pf)) > 0L ||
    length(problem@constraints) > 0L ||
    any(vapply(variables(problem), function(v) length(domain(v)) > 0L,
               logical(1L)))
}

# -- pick_default_solver -------------------------------------------
## CVXPY SOURCE: problem_form.py:336-387
## The default solver for the problem's structure: commercial first when that
## is switched on, then the open-source defaults. Returns a solver instance or
## NULL when no suitable installed solver is found (the caller falls back).
## "Installed" is `.solver_package_available`, so the solver-exclusion API
## (exclude_solvers()) makes a solver look uninstalled, as CVXPY's
## `is_installed()` would.
## `use_quad_obj` is R-only (recorded at `can_solve`, solver.R): with FALSE the
## QP -> OSQP rule is skipped, since the chain will hand OSQP the objective's
## SOC, and the problem falls through to the conic rules.
pick_default_solver <- function(problem_form, use_quad_obj = TRUE) {
  .get <- function(solver_map, name) {
    inst <- solver_map[[name]]
    if (!is.null(inst) && .solver_package_available(name)) inst else NULL
  }

  ## CVXPY SOURCE: problem_form.py:356 (`if s.DEFAULT_TO_COMMERCIAL_SOLVERS:`);
  ## see settings.R for why CVXR reads the flag per call.
  if (.default_to_commercial_solvers()) {
    for (solver_name in .COMMERCIAL_CONIC_SOLVERS) {
      solver <- .get(SOLVER_MAP_CONIC, solver_name)
      if (!is.null(solver) && can_solve(solver, problem_form, use_quad_obj)) return(solver)
    }
  }

  ## Mixed-integer: LP -> HIGHS, non-LP -> SCIP.
  if (.pf_is_mixed_integer(problem_form)) {
    if (.pf_is_lp(problem_form)) {
      solver <- .get(SOLVER_MAP_CONIC, HIGHS_SOLVER)
      if (!is.null(solver)) return(solver)
    }
    solver <- .get(SOLVER_MAP_CONIC, SCIP_SOLVER)
    if (!is.null(solver) && can_solve(solver, problem_form, use_quad_obj)) return(solver)
    return(NULL)
  }

  ## LP -> Clarabel.
  if (.pf_within_qp_cones(.pf_cones(problem_form))) {
    return(.get(SOLVER_MAP_CONIC, CLARABEL_SOLVER))
  }

  ## QP -> OSQP.
  if (isTRUE(use_quad_obj) && .pf_has_quadratic_objective(problem_form) &&
      .pf_within_qp_cones(.pf_cones(problem_form, quad_obj = TRUE))) {
    return(.get(SOLVER_MAP_QP, OSQP_SOLVER))
  }

  ## SDP -> SCS.
  if (.pf_has_cone(.pf_cones(problem_form), PSD)) {
    return(.get(SOLVER_MAP_CONIC, SCS_SOLVER))
  }

  ## Everything else (SOCP, ExpCone, PowCone, ...) -> Clarabel.
  .get(SOLVER_MAP_CONIC, CLARABEL_SOLVER)
}

# -- make_problem_form ---------------------------------------------
## CVXPY SOURCE: problem_form.py:390-435
## Check DCP (or DGP under gp = TRUE), then build the ProblemForm. The checks
## run BEFORE any solver is resolved, as upstream's do. The messages are
## CVXR's own (recorded divergence, ADR D_PORT.14).
make_problem_form <- function(problem, gp, ignore_dpp) {
  if (!gp && !is_dcp(problem)) {
    if (is_dgp(problem)) {
      cli_abort(c(
        "Problem is not DCP compliant.",
        "i" = "However, the problem is DGP. Try {.code psolve(problem, gp = TRUE)}."
      ))
    }
    if (is_dqcp(problem)) {
      cli_abort(c(
        "Problem is not DCP compliant.",
        "i" = "However, the problem is DQCP. Try {.code psolve(problem, qcp = TRUE)}."
      ))
    }
    cli_abort("Problem is not DCP compliant.")
  } else if (gp && !is_dgp(problem)) {
    cli_abort(c(
      "Problem is not DGP compliant.",
      "i" = "Remove {.code gp = TRUE} or reformulate as a geometric program."
    ), class = "DGPError")   ## CVXPY: raise DGPError, problem_form.py:432
  }
  eval_params <- isTRUE(ignore_dpp) && length(parameters(problem)) > 0L
  ProblemForm(problem, gp = gp, eval_params = eval_params)
}