R/231_reductions_cvx_attr2constr.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: reductions/cvx_attr2constr.py
## CvxAttr2Constr -- expand convex variable attributes into constraints
##
## Phase 5b scope: handle nonneg, nonpos, PSD, NSD, symmetric.
## Deferred: diag, sparsity, bounds.


# -- upper_tri_to_full: sparse matrix to fill symmetric from upper-tri --
## CVXPY SOURCE: atoms/affine/upper_tri.py lines 147-176
## Returns an (n^2 x n*(n+1)/2) sparse matrix A such that
## (A %*% v) reshaped as (n, n) is symmetric.

upper_tri_to_full <- function(n) {
  entries <- (n * (n + 1L)) %/% 2L

  ## Upper-triangular row/col indices (including diagonal), row-major order
  ## CVXPY uses np.triu_indices(n) which is row-major
  idx <- which(upper.tri(matrix(0, n, n), diag = TRUE), arr.ind = TRUE)
  idx <- idx[order(idx[, 1L], idx[, 2L]), , drop = FALSE]
  rows <- idx[, 1L] - 1L  # 0-based
  cols <- idx[, 2L] - 1L  # 0-based

  ## Mask for off-diagonal entries (i != j)
  mask <- rows != cols

  ## Row indices in the n*n flattened (column-major for R) representation
  ## CVXPY uses rows * n + cols (row-major), but since we'll reshape with

  ## matrix(, n, n) in R (column-major), use cols * n + rows
  row_idx <- c(cols * n + rows, rows[mask] * n + cols[mask]) + 1L  # 1-based
  col_idx <- c(seq_len(entries), which(mask)) # 1-based
  values <- rep(1.0, length(col_idx))

  Matrix::sparseMatrix(i = row_idx, j = col_idx, x = values,
                       dims = c(as.integer(n * n), entries))
}


# -- CvxAttr2Constr reduction class -------------------------------
## CVXPY SOURCE: cvx_attr2constr.py lines 105-215

CvxAttr2Constr <- new_class("CvxAttr2Constr", parent = Reduction,
  package = "CVXR",
  constructor = function() {
    new_object(S7_object(),
      .cache = new.env(parent = emptyenv())
    )
  }
)

method(reduction_accepts, CvxAttr2Constr) <- function(x, problem, ...) TRUE

method(reduction_apply, CvxAttr2Constr) <- function(x, problem, ...) {
  vars_ <- variables(problem)
  if (length(convex_attributes(vars_)) == 0L) {
    return(list(problem, list()))
  }

  ## For each unique variable, create replacement and constraints
  id2new_var <- new.env(hash = TRUE, parent = emptyenv())
  id2new_obj <- new.env(hash = TRUE, parent = emptyenv())
  id2old_var <- new.env(hash = TRUE, parent = emptyenv())
  ## Collect attribute constraint chunks per variable (up to ~4 each: PSD/NSD + nonneg + nonpos + bounds)
  attr_constr_chunks <- vector("list", length(vars_))
  n_attr_chunks <- 0L

  for (var in vars_) {
    vid <- as.character(var@id)
    if (exists(vid, envir = id2new_var, inherits = FALSE)) next
    assign(vid, var, envir = id2old_var)

    attrs <- var@attributes
    has_sym <- !is.null(attrs$symmetric) && attrs$symmetric
    has_psd <- !is.null(attrs$PSD) && attrs$PSD
    has_nsd <- !is.null(attrs$NSD) && attrs$NSD
    has_nonneg <- !is.null(attrs$nonneg) && attrs$nonneg
    has_nonpos <- !is.null(attrs$nonpos) && attrs$nonpos
    has_bounds <- !is.null(attrs$bounds) && is.list(attrs$bounds)

    ## CVXPY approach: copy all attributes, then clear only the
    ## "reduction_attributes" (symmetric, PSD, NSD). Preserve boolean,
    ## integer, nonneg, nonpos, bounds on the new variable.
    ## reduction_attributes = c("symmetric", "PSD", "NSD")
    ## (diag/sparsity not implemented in R)
    new_attrs <- attrs
    needs_new_var <- FALSE
    for (key in c("symmetric", "PSD", "NSD")) {
      if (!is.null(new_attrs[[key]]) && new_attrs[[key]]) {
        needs_new_var <- TRUE
        new_attrs[[key]] <- FALSE
      }
    }
    ## Also mark as needing new var if nonneg/nonpos/bounds present
    if (has_nonneg || has_nonpos || has_bounds) needs_new_var <- TRUE

    if (has_sym || has_psd || has_nsd) {
      ## Create upper-triangular variable and fill to full matrix
      n <- var@shape[1L]
      tri_size <- (n * (n + 1L)) %/% 2L
      ## New variable with same id but stripped symmetric/PSD/NSD attrs
      ## Preserve boolean/integer on the upper-tri variable
      upper_tri_args <- list(c(tri_size, 1L), var_id = var@id)
      if (!is.null(new_attrs$boolean) && new_attrs$boolean)
        upper_tri_args$boolean <- TRUE
      if (!is.null(new_attrs$integer) && new_attrs$integer)
        upper_tri_args$integer <- TRUE
      upper_tri_var <- do.call(Variable, upper_tri_args)
      assign(vid, upper_tri_var, envir = id2new_var)

      ## Fill coefficient: n^2 x tri_size sparse matrix
      fill_coeff <- Constant(upper_tri_to_full(n))
      full_mat <- fill_coeff %*% upper_tri_var
      obj <- reshape_expr(full_mat, c(n, n))

      ## Map original var's python id to the new expression
      assign(as.character(var@id), obj, envir = id2new_obj)

      ## Add PSD/NSD constraint on the full matrix
      var_constrs <- list()
      if (has_psd) {
        var_constrs <- list(PSD(obj))
      } else if (has_nsd) {
        var_constrs <- list(PSD(-obj))
      }
    } else if (needs_new_var) {
      ## Create replacement variable preserving boolean/integer
      new_var_args <- list(var@shape, var_id = var@id)
      if (!is.null(new_attrs$boolean) && new_attrs$boolean)
        new_var_args$boolean <- TRUE
      if (!is.null(new_attrs$integer) && new_attrs$integer)
        new_var_args$integer <- TRUE
      obj <- do.call(Variable, new_var_args)
      assign(vid, obj, envir = id2new_var)
      assign(as.character(var@id), obj, envir = id2new_obj)
      var_constrs <- list()
    } else {
      ## No attribute to reduce -- keep the variable as-is
      obj <- var
      assign(vid, var, envir = id2new_var)
      assign(as.character(var@id), var, envir = id2new_obj)
      var_constrs <- list()
    }

    ## Add nonneg/nonpos/bounds constraints independently
    ## CVXPY SOURCE: leaf.py _bound_domain() -- each is a separate `if`
    obj <- get(as.character(var@id), envir = id2new_obj)
    if (has_nonneg) {
      var_constrs <- c(var_constrs, list(NonNeg(obj)))
    }
    if (has_nonpos) {
      var_constrs <- c(var_constrs, list(NonPos(obj)))
    }
    if (has_bounds) {
      lb <- attrs$bounds[[1L]]
      ub <- attrs$bounds[[2L]]
      if (any(is.finite(lb))) {
        var_constrs <- c(var_constrs, list(obj >= lb))
      }
      if (any(is.finite(ub))) {
        var_constrs <- c(var_constrs, list(obj <= ub))
      }
    }
    if (length(var_constrs) > 0L) {
      n_attr_chunks <- n_attr_chunks + 1L
      attr_constr_chunks[[n_attr_chunks]] <- var_constrs
    }
  }

  ## Substitute variables in the objective and constraints via tree_copy
  new_obj <- tree_copy(problem@objective, id_objects = id2new_obj)
  n_pcons <- length(problem@constraints)
  tree_copy_chunks <- vector("list", n_pcons)
  cons_id_map <- list()
  for (i in seq_len(n_pcons)) {
    con <- problem@constraints[[i]]
    new_con <- tree_copy(con, id_objects = id2new_obj)
    tree_copy_chunks[[i]] <- list(new_con)
    cons_id_map[[as.character(con@id)]] <- new_con@id
  }

  ## Combine all constraint chunks: attribute constraints + tree-copied constraints
  all_constr_chunks <- c(attr_constr_chunks[seq_len(n_attr_chunks)], tree_copy_chunks)
  new_constrs <- unlist(all_constr_chunks, recursive = FALSE)
  if (is.null(new_constrs)) new_constrs <- list()

  new_problem <- Problem(new_obj, new_constrs)
  inverse_data <- list(id2new_var = id2new_var,
                       id2old_var = id2old_var,
                       cons_id_map = cons_id_map)
  list(new_problem, inverse_data)
}

method(reduction_invert, CvxAttr2Constr) <- function(x, solution, inverse_data, ...) {
  if (length(inverse_data) == 0L) return(solution)

  id2new_var <- inverse_data$id2new_var
  id2old_var <- inverse_data$id2old_var
  cons_id_map <- inverse_data$cons_id_map

  pvars <- list()
  old_ids <- ls(id2old_var, all.names = TRUE)
  for (vid in old_ids) {
    old_var <- get(vid, envir = id2old_var)
    new_var <- get(vid, envir = id2new_var)
    new_vid <- as.character(new_var@id)
    if (!is.null(solution@primal_vars[[new_vid]])) {
      raw_val <- solution@primal_vars[[new_vid]]
      ## Recover full matrix for symmetric/PSD/NSD
      attrs <- old_var@attributes
      has_sym <- !is.null(attrs$symmetric) && attrs$symmetric
      has_psd <- !is.null(attrs$PSD) && attrs$PSD
      has_nsd <- !is.null(attrs$NSD) && attrs$NSD
      if (has_sym || has_psd || has_nsd) {
        n <- old_var@shape[1L]
        val <- numeric(n * n)
        ## Indices for upper triangle (row-major, 0-based)
        idx <- which(upper.tri(matrix(0, n, n), diag = TRUE), arr.ind = TRUE)
        idx <- idx[order(idx[, 1L], idx[, 2L]), , drop = FALSE]
        flat_val <- as.numeric(raw_val)
        ## Fill upper triangle and mirror
        full <- matrix(0, n, n)
        for (k in seq_len(nrow(idx))) {
          i <- idx[k, 1L]
          j <- idx[k, 2L]
          full[i, j] <- flat_val[k]
          full[j, i] <- flat_val[k]
        }
        pvars[[vid]] <- full
      } else {
        pvars[[vid]] <- raw_val
      }
    }
  }

  ## Remap dual variables
  dvars <- list()
  for (orig_id in names(cons_id_map)) {
    new_id <- as.character(cons_id_map[[orig_id]])
    if (!is.null(solution@dual_vars[[new_id]])) {
      dvars[[orig_id]] <- solution@dual_vars[[new_id]]
    }
  }

  Solution(solution@status, solution@opt_val, pvars, dvars, solution@attr)
}