R/260_reductions_cone2cone_exact.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: reductions/cone2cone/exact.py
## PARTIAL PORT: in: PSDToSvecPSD + ExactCone2Cone (apply/invert, recovery
##   chains). out: NonPosConversion -- implemented at reductions/utilities.R
##   (`nonpos2nonneg`); PowNDConversion -- implemented at the dcp2cone
##   canonicalizers, and moving it here needs cone_tree.R plus tree-based dual
##   recovery; SOCConversion -- not needed (no CVXR solver is PSD-only).
##   Sequencing: notes/dual_value_shape_port_plan.md.
##
## Exact cone-to-cone conversions: rewrites that replace one cone with a set of
## simpler cones the target solver actually accepts, preserving equivalence.
## New at CVXPY 1.9.0 (PR #3079/#3223/#3268); ported here at 1.9.2 under ADR
## D_19.6, which reopened the D_19.1 cluster-2 defer.
##
## Each conversion is a list with
##   source        -- the cone class converted
##   targets       -- the cone classes it maps to
##   canonicalize(con, args, solver_context) -> list(canon_constr, aux_constrs)
##   recover_dual(cons, dual_var, inverse_data, dvars) [optional]
##       -- map a solver dual on the converted cone back to the original's.
## (`apply_hook` exists upstream for the PowConeND tree and packed-SOC cases;
## neither conversion is ported yet, so no hook is needed here. Add it with the
## conversion, not before.)
##
## PORTED SO FAR: PSD -> SvecPSD. The other three upstream conversions are NOT
## yet routed through this reduction:
##   - NonPos -> NonNeg lives in reductions/utilities.R (`nonpos2nonneg`);
##   - PowConeND -> PowCone3D lives in the dcp2cone canonicalizers and would
##     need cone_tree.R (LeafNode/SplitNode/SingleVarNode) plus its tree-based
##     dual recovery;
##   - SOC -> PSD is unneeded: every CVXR conic solver that lacks PSD also
##     lacks nothing SOC-shaped, and no solver here is PSD-only.
## Moving the first two here is follow-on work, tracked in
## notes/dual_value_shape_port_plan.md; until then they are deviations that are
## recorded rather than hidden.


# -- EXACT_CONE_CONVERSIONS ----------------------------------------
## CVXPY SOURCE: exact.py lines 61-66 -- {source_cone: {target_cones}}, which
## must form a DAG.  Consulted by `.expand_cones` (reductions/solvers/solver.R)
## when deciding whether a solver that does not accept a cone directly can
## still take it after conversion.
##
## Represented as a list of (source, targets) pairs because S7 class objects
## cannot be list names.

EXACT_CONE_CONVERSIONS <- list(
  list(source = PSD, targets = list(SvecPSD))
)


# -- PSDToSvecPSD --------------------------------------------------
## CVXPY SOURCE: exact.py lines 420-446
##
## PSD -> SvecPSD via scaled vectorization.  Which triangle, and whether the
## off-diagonal is sqrt(2)-scaled, come from the solver context -- the same two
## values the solver declares as PSD_TRIANGLE_KIND / PSD_SQRT2_SCALING.

PSDToSvecPSD <- list(
  source  = PSD,
  targets = list(SvecPSD),

  ## CVXPY SOURCE: exact.py lines 430-437
  canonicalize = function(con, args, solver_context = NULL) {
    X <- args[[1L]]
    n <- .shape(X)[2L]
    M <- psd_format_mat(con, solver_context@psd_triangle_kind,
                        solver_context@psd_sqrt2_scaling)
    ## `vec(X)` is CVXPY's `vec(X, order='F')`: R is column-major throughout.
    svec_expr <- Constant(M) %*% vec(X)
    list(SvecPSD(svec_expr, n = n), list())
  },

  ## CVXPY SOURCE: exact.py lines 439-446
  recover_dual = function(cons, dual_var, inverse_data, dvars) {
    ctx <- get0("solver_context", envir = inverse_data@.extra, ifnotfound = NULL)
    n <- .arg_shape(cons)[2L]
    ## Upstream reshapes `tri_to_full`'s output to the argument's shape; CVXR's
    ## `tri_to_full` already returns the (n, n) matrix, and a PSD constraint's
    ## argument is square by construction (psd.R constructor).
    tri_to_full(dual_var, n, ctx@psd_triangle_kind, ctx@psd_sqrt2_scaling)
  }
)

## All conversions this reduction knows how to perform.
## CVXPY SOURCE: exact.py line 451
.EXACT_CONVERSIONS <- list(PSDToSvecPSD)


# -- ExactCone2Cone ------------------------------------------------
## CVXPY SOURCE: exact.py lines 449-476
##
## `target_cones` restricts the reduction to the conversions the chain actually
## asked for (computed by `.expand_cones`); NULL means "all of them".
##
## Upstream subclasses `Canonicalization` to reuse `canonicalize_tree`.  CVXR's
## `Canonicalization` walks EXPRESSION trees through S7 generics and has no
## canon_methods dict, and by this point in the chain every constraint argument
## is already affine -- so the expression walk would be a no-op.  This is a
## `Reduction` that converts the constraints directly, which is what
## `canonicalize_tree` reduces to for affine arguments.

ExactCone2Cone <- new_class("ExactCone2Cone", parent = Reduction,
  package = "CVXR",
  properties = list(
    target_cones   = class_list,
    solver_context = class_any
  ),
  constructor = function(target_cones = list(), solver_context = NULL) {
    if (FALSE) new_object(S7_object())  ## S7 static-check guard
    .fast_new(ExactCone2Cone, S7_object(),
      .cache         = new.env(parent = emptyenv()),
      target_cones   = target_cones,
      solver_context = solver_context
    )
  }
)

## The conversion applicable to `con`, or NULL.
## CVXPY SOURCE: exact.py lines 457-460 (the `target_cones` filter) and the
## `type(current) in self.canon_methods` test at line 486.
.exact_conversion_for <- function(x, con) {
  for (conv in .EXACT_CONVERSIONS) {
    if (!.s7_is(con, conv$source)) next
    if (length(x@target_cones) > 0L &&
        !any(vapply(x@target_cones, identical, logical(1L), conv$source))) next
    return(conv)
  }
  NULL
}

## CVXPY SOURCE: exact.py lines 478-503
## Follow the conversion chain transitively, recording the constraints passed
## through so `invert` can unwind them in reverse.
.convert_constraint <- function(x, constraint, inverse_data, canon_constraints) {
  chain <- list()
  current <- constraint

  repeat {
    conv <- .exact_conversion_for(x, current)
    if (is.null(conv)) break
    chain[[length(chain) + 1L]] <- current
    res <- conv$canonicalize(current, .args(current), x@solver_context)
    canon_constr <- res[[1L]]
    aux_constr <- res[[2L]]

    ## Auxiliary constraints convert recursively too.
    for (aux in aux_constr) {
      canon_constraints <- .convert_constraint(x, aux, inverse_data,
                                               canon_constraints)
    }
    current <- canon_constr
  }

  canon_constraints[[length(canon_constraints) + 1L]] <- current
  assign(as.character(.id(constraint)), .id(current),
         envir = inverse_data@cons_id_map)
  if (length(chain) > 0L) {
    chains <- get0("recovery_chains", envir = inverse_data@.extra,
                   ifnotfound = list())
    chains[[as.character(constraint@id)]] <- chain
    assign("recovery_chains", chains, envir = inverse_data@.extra)
  }
  canon_constraints
}

## CVXPY SOURCE: exact.py lines 505-520
method(reduction_apply, ExactCone2Cone) <- function(x, problem, ...) {
  inverse_data <- InverseData(problem)
  assign("solver_context", x@solver_context, envir = inverse_data@.extra)

  canon_constraints <- list()
  for (constraint in problem@constraints) {
    canon_constraints <- .convert_constraint(x, constraint, inverse_data,
                                             canon_constraints)
  }

  list(Problem(problem@objective, canon_constraints), inverse_data)
}

## CVXPY SOURCE: exact.py lines 522-551
method(reduction_invert, ExactCone2Cone) <- function(x, solution, inverse_data, ...) {
  ## Re-key the duals from converted-constraint ids onto original ids.
  ## O(n) via one vectorized match; see `.remap_by_id_map`.
  dvars <- .remap_by_id_map(solution@dual_vars, inverse_data@cons_id_map)

  if (length(dvars) == 0L) {
    return(Solution(status      = solution@status,
                    opt_val     = solution@opt_val,
                    primal_vars = solution@primal_vars,
                    dual_vars   = dvars,
                    attr        = solution@attr))
  }

  ## Unwind each recovery chain in reverse, so transitive conversions come
  ## apart in the order they were applied.
  chains <- get0("recovery_chains", envir = inverse_data@.extra,
                 ifnotfound = list())
  for (orig_id in names(chains)) {
    if (is.null(dvars[[orig_id]])) next
    dual <- dvars[[orig_id]]
    chain <- chains[[orig_id]]
    for (i in rev(seq_along(chain))) {
      cons <- chain[[i]]
      conv <- .exact_conversion_for(x, cons)
      if (!is.null(conv) && !is.null(conv$recover_dual)) {
        dual <- conv$recover_dual(cons, dual, inverse_data, dvars)
      }
    }
    dvars[[orig_id]] <- dual
  }

  Solution(status      = solution@status,
           opt_val     = solution@opt_val,
           primal_vars = solution@primal_vars,
           dual_vars   = dvars,
           attr        = solution@attr)
}