R/258_reductions_dcp2cone_cone_matrix_stuffing.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: reductions/dcp2cone/cone_matrix_stuffing.py
## ConeMatrixStuffing -- convert affine expressions to sparse matrices
## ConeDims -- summary of cone dimensions


# ==================================================================
# ConeDims -- summary of cone dimensions present in constraints
# ==================================================================
## CVXPY SOURCE: cone_matrix_stuffing.py lines 57-141

ConeDims <- new_class("ConeDims", package = "CVXR",
  properties = list(
    zero   = class_integer,
    nonneg = class_integer,
    exp    = class_integer,
    soc    = class_integer,     # vector of SOC cone dims
    psd    = class_integer,     # vector of PSD sizes (n for n x n)
    p3d    = class_numeric,     # vector of PowCone3D alphas
    pnd    = class_list         # list of numeric vectors for PowConeND
  ),
  constructor = function(constr_map) {
    if (FALSE) new_object(S7_object())  ## S7 static-check guard
    ## Zero cone dimension
    zero_constrs <- constr_map[["Zero"]]
    zero <- if (length(zero_constrs) == 0L) 0L
            else as.integer(sum(vapply(zero_constrs, constr_size, numeric(1L))))

    ## NonNeg cone dimension
    nn_constrs <- constr_map[["NonNeg"]]
    nonneg <- if (length(nn_constrs) == 0L) 0L
              else as.integer(sum(vapply(nn_constrs, constr_size, numeric(1L))))

    ## Exponential cone count
    exp_constrs <- constr_map[["ExpCone"]]
    exp_count <- if (length(exp_constrs) == 0L) 0L
                 else as.integer(sum(vapply(exp_constrs, num_cones, numeric(1L))))

    ## SOC dimensions (one entry per individual cone)
    soc_constrs <- constr_map[["SOC"]]
    soc <- if (length(soc_constrs) == 0L) integer(0)
           else as.integer(unlist(lapply(soc_constrs, cone_sizes)))

    ## PSD dimensions (the matrix side length n, whether the constraint is a
    ## full-matrix PSD or the packed SvecPSD the solver actually receives).
    ## CVXPY SOURCE: cone_matrix_stuffing.py:92-93
    psd_constrs <- c(constr_map[["PSD"]], constr_map[["SvecPSD"]])
    psd <- if (length(psd_constrs) == 0L) integer(0)
           else as.integer(unlist(lapply(psd_constrs, cone_sizes)))

    ## PowCone3D alphas
    p3d_constrs <- constr_map[["PowCone3D"]]
    p3d <- if (length(p3d_constrs) == 0L) numeric(0)
           else unlist(lapply(p3d_constrs, function(c) as.numeric(value(c@alpha))))

    ## PowConeND alphas
    pnd_constrs <- constr_map[["PowConeND"]]
    pnd <- if (length(pnd_constrs) == 0L) list()
           else {
             alpha_chunks <- vector("list", length(pnd_constrs))
             for (ci in seq_along(pnd_constrs)) {
               a <- value(pnd_constrs[[ci]]@alpha)
               if (is.matrix(a)) {
                 ## Each column is one cone's alpha vector
                 alpha_chunks[[ci]] <- lapply(seq_len(ncol(a)), function(j) a[, j])
               } else {
                 alpha_chunks[[ci]] <- list(as.numeric(a))
               }
             }
             unlist(alpha_chunks, recursive = FALSE)
           }

    .fast_new(ConeDims, S7_object(),
      zero = zero, nonneg = nonneg, exp = exp_count,
      soc = soc, psd = psd, p3d = p3d, pnd = pnd)
  }
)

method(print, ConeDims) <- function(x, ...) {
  cat(sprintf("%d equalities, %d inequalities, %d exponential cones\n",
    x@zero, x@nonneg, x@exp))
  cat(sprintf("SOC: %s, PSD: %s\n",
    paste(x@soc, collapse = ","),
    paste(x@psd, collapse = ",")))
  invisible(x)
}

# -- dims_to_solver_dict -------------------------------------------
## CVXPY SOURCE: conic_solver.py lines 87-97

dims_to_solver_dict <- function(cone_dims) {
  list(
    f   = cone_dims@zero,
    l   = cone_dims@nonneg,
    q   = cone_dims@soc,
    ep  = cone_dims@exp,
    s   = cone_dims@psd,
    p   = cone_dims@p3d,
    pnd = cone_dims@pnd
  )
}


# -- extract_mip_idx -----------------------------------------------
## CVXPY SOURCE: matrix_stuffing.py lines 79-96
## Maps per-variable boolean/integer flags to global flattened variable indices.

.extract_lower_bounds <- function(vars, x_length) {
  has_bounds <- any(vapply(vars, function(v) {
    b <- .attributes(v)$bounds
    (!is.null(b) && is.list(b) && !all(is.infinite(as.numeric(b[[1L]])))) ||
      is_nonneg(v)
  }, logical(1L)))
  if (!has_bounds) return(NULL)

  lower <- rep(-Inf, x_length)
  offset <- 0L
  for (v in vars) {
    sz <- expr_size(v)
    idx <- (offset + 1L):(offset + sz)
    b <- .attributes(v)$bounds
    if (!is.null(b) && is.list(b)) {
      if (.s7_is(b[[1L]], Expression) || inherits(b[[1L]], "Matrix")) {
        cli_abort("Sparse or parametric bounds should not reach matrix stuffing.")
      }
      lower[idx] <- as.numeric(array(b[[1L]], dim = .shape(v)))
    }
    if (is_nonneg(v)) lower[idx] <- pmax(lower[idx], 0)
    offset <- offset + sz
  }
  lower
}

.extract_upper_bounds <- function(vars, x_length) {
  has_bounds <- any(vapply(vars, function(v) {
    b <- .attributes(v)$bounds
    (!is.null(b) && is.list(b) && !all(is.infinite(as.numeric(b[[2L]])))) ||
      is_nonpos(v)
  }, logical(1L)))
  if (!has_bounds) return(NULL)

  upper <- rep(Inf, x_length)
  offset <- 0L
  for (v in vars) {
    sz <- expr_size(v)
    idx <- (offset + 1L):(offset + sz)
    b <- .attributes(v)$bounds
    if (!is.null(b) && is.list(b)) {
      if (.s7_is(b[[2L]], Expression) || inherits(b[[2L]], "Matrix")) {
        cli_abort("Sparse or parametric bounds should not reach matrix stuffing.")
      }
      upper[idx] <- as.numeric(array(b[[2L]], dim = .shape(v)))
    }
    if (is_nonpos(v)) upper[idx] <- pmin(upper[idx], 0)
    offset <- offset + sz
  }
  upper
}

.has_parametric_bounds <- function(vars) {
  any(vapply(vars, function(v) {
    b <- .attributes(v)$bounds
    !is.null(b) && is.list(b) &&
      any(vapply(b, function(x) .s7_is(x, Expression), logical(1L)))
  }, logical(1L)))
}

.bound_expr_for_tensor <- function(v, which) {
  bound_idx <- if (identical(which, "lower")) 1L else 2L
  default_val <- if (identical(which, "lower")) -Inf else Inf
  sz <- expr_size(v)
  b <- .attributes(v)$bounds

  if (identical(which, "lower") && is_nonneg(v)) {
    return(Constant(matrix(0, sz, 1L)))
  }
  if (identical(which, "upper") && is_nonpos(v)) {
    return(Constant(matrix(0, sz, 1L)))
  }

  if (!is.null(b) && is.list(b)) {
    slot <- b[[bound_idx]]
    if (.s7_is(slot, Expression)) {
      if (expr_is_scalar(slot) && sz > 1L) {
        slot <- cvxr_promote(slot, .shape(v))
      }
      if (!identical(as.integer(.shape(slot)), as.integer(c(sz, 1L)))) {
        slot <- reshape_expr(slot, c(sz, 1L), order = "F")
      }
      return(slot)
    }
    if (inherits(slot, "Matrix")) {
      cli_abort("Sparse bounds should not reach matrix stuffing.")
    }
    return(Constant(matrix(as.numeric(array(slot, dim = .shape(v))), sz, 1L)))
  }

  Constant(matrix(default_val, sz, 1L))
}

.extract_bounds_tensor <- function(vars, x_length, param_to_size, param_id_map,
                                   which = c("lower", "upper")) {
  which <- match.arg(which)
  op_list <- lapply(vars, function(v) {
    canonical_form(.bound_expr_for_tensor(v, which))[[1L]]
  })
  get_problem_matrix_tensor(
    op_list,
    id_to_col = integer(0),
    var_length = 0L,
    param_to_size = param_to_size,
    param_id_map = param_id_map
  )
}

.extract_mip_idx <- function(problem, inverse_data) {
  vars <- variables(problem)
  bool_chunks <- vector("list", length(vars))
  int_chunks  <- vector("list", length(vars))

  for (i in seq_along(vars)) {
    var <- vars[[i]]
    vid <- as.character(.id(var))
    offset <- inverse_data@var_offsets[[vid]]
    sz <- expr_size(var)
    ## CVXPY SOURCE: matrix_stuffing.py:107-119 -- shift each variable's OWN
    ## index set by its offset in the stuffed vector. Upstream ravels a numpy
    ## multi-index; CVXR's canonical form is already a flat 1-based vector
    ## (`.mip_idx`, expressions/leaf.R), so the shift is all that is left.
    ##
    ## Two INDEPENDENT ifs, matching upstream: a variable may carry a boolean
    ## index list AND an integer one (upstream's own test_bool_int_variable
    ## does exactly that). The `else if` this replaces silently dropped the
    ## integer set in that case.
    if (length(var@.boolean_idx) > 0L) {
      bool_chunks[[i]] <- offset + var@.boolean_idx
    }
    if (length(var@.integer_idx) > 0L) {
      int_chunks[[i]] <- offset + var@.integer_idx
    }
  }

  bool_idx <- unlist(bool_chunks)
  int_idx  <- unlist(int_chunks)

  ## INVARIANT: the two sets are disjoint, boolean winning.
  ##
  ## Boolean is a SUBSET of integer -- {0,1} is the integers intersected with
  ## [0,1] -- so an entry named by both is not a conflict, and applying both
  ## constraints yields the boolean one. Upstream never resolves this: each
  ## solver marks integrality from the union and then applies the [0,1]
  ## restriction for the boolean list ON TOP (scipy_conif.py:132-141,
  ## cbc_conif.py:154-160), so boolean wins there by write order.
  ##
  ## CVXR's eight MIP interfaces happen to apply boolean FIRST and integer
  ## second, so inheriting upstream's overlap would silently give the opposite
  ## answer -- and would leave the semantics resting on the statement order in
  ## six separate files, which is the same duplicated-knowledge shape that
  ## produced the transposed PSD duals. Enforcing disjointness once here is
  ## order-independent and behaviourally identical to upstream.
  if (!is.null(bool_idx) && !is.null(int_idx)) {
    int_idx <- setdiff(int_idx, bool_idx)
  }
  list(bool_idx = if (is.null(bool_idx)) integer(0) else bool_idx,
       int_idx  = if (is.null(int_idx))  integer(0) else int_idx)
}

# ==================================================================
# ConeMatrixStuffing -- reduction from affine expressions to matrices
# ==================================================================
## CVXPY SOURCE: cone_matrix_stuffing.py lines 321-470

ConeMatrixStuffing <- new_class("ConeMatrixStuffing", parent = Reduction,
  package = "CVXR",
  properties = list(
    quad_obj = class_logical
  ),
  constructor = function(quad_obj = FALSE) {
    if (FALSE) new_object(S7_object())  ## S7 static-check guard
    .fast_new(ConeMatrixStuffing, S7_object(),
      .cache = new.env(parent = emptyenv()),
      quad_obj = quad_obj
    )
  }
)

## accepts: affine (or quadratic when quad_obj) objective (Minimize),
## no convex attributes, all constraint args affine
## CVXPY SOURCE: cone_matrix_stuffing.py lines 335-342
method(reduction_accepts, ConeMatrixStuffing) <- function(x, problem, ...) {
  is_min <- .s7_is(problem@objective, Minimize)
  obj_expr <- .args(problem@objective)[[1L]]
  valid_obj <- if (x@quad_obj) {
    is_affine(obj_expr) || is_quadratic(obj_expr)
  } else {
    is_affine(obj_expr)
  }
  no_cvx_attr <- length(convex_attributes(variables(problem))) == 0L
  aff_con <- are_args_affine(problem@constraints)
  is_min && valid_obj && no_cvx_attr && aff_con
}

## apply: lower constraints and extract A, b, c matrices
## CVXPY SOURCE: cone_matrix_stuffing.py lines 360-430
method(reduction_apply, ConeMatrixStuffing) <- function(x, problem, ...) {
  inverse_data <- InverseData(problem)

  ## Step 1: Lower constraints -- pre-allocate
  n_cons_raw <- length(problem@constraints)
  cons <- vector("list", n_cons_raw)
  for (i in seq_len(n_cons_raw)) {
    con <- problem@constraints[[i]]
    if (.s7_is(con, Equality)) {
      con <- lower_equality(con)
    } else if (.s7_is(con, Inequality)) {
      con <- lower_ineq_to_nonneg(con)
    } else if (.s7_is(con, NonPos)) {
      con <- nonpos2nonneg(con)
    } else if (.s7_is(con, SOC) && con@axis == 1L) {
      ## Transpose X for axis=1 -> axis=2
      con <- SOC(.args(con)[[1L]], t(.args(con)[[2L]]), axis = 2L,
                 constr_id = .id(con))
    } else if (.s7_is(con, PowConeND) && con@axis == 1L) {
      ## CVXPY SOURCE: cone_matrix_stuffing.py lines 382-388
      ## Transpose W and alpha for axis=1 -> axis=2
      con <- PowConeND(t(con@.W), con@.z,
                        t(con@alpha), axis = 2L,
                        constr_id = .id(con))
    }
    ## ExpCone/PowCone3D flattening: only needed if multidimensional
    ## For Phase 5b, our canonicalizers produce 1D args; skip flattening for now
    cons[[i]] <- con
  }

  ## Step 2: Group and order constraints
  ## CVXPY order: Zero, NonNeg, SOC, PSD, ExpCone, PowCone3D, PowConeND
  constr_map <- group_constraints(cons)
  ## CVXPY SOURCE: cone_matrix_stuffing.py:408-412 -- SvecPSD occupies the same
  ## slot in the cone ordering as PSD, immediately after it.
  ordered_cons <- c(constr_map[["Zero"]], constr_map[["NonNeg"]],
                    constr_map[["SOC"]], constr_map[["PSD"]],
                    constr_map[["SvecPSD"]],
                    constr_map[["ExpCone"]], constr_map[["PowCone3D"]],
                    constr_map[["PowConeND"]])

  ## Step 3: Store constraint ID mapping (identity -- already lowered)
  ## CVXPY SOURCE: cone_matrix_stuffing.py line 406
  id2cons <- new.env(hash = TRUE, parent = emptyenv())
  for (con in ordered_cons) {
    assign(as.character(.id(con)), .id(con), envir = inverse_data@cons_id_map)
    assign(as.character(.id(con)), con, envir = id2cons)
  }

  ## Step 4: Create CoeffExtractor and extract objective
  extractor <- CoeffExtractor(inverse_data)

  ## Create the flattened variable
  x_var <- Variable(c(extractor@x_length, 1L))

  ## Detect DPP path: parameters exist and no EvalParams was applied
  ## (i.e., parameters are still present in the expressions).
  has_params <- length(parameters(problem)) > 0L

  ## Extract objective
  P_mat <- NULL
  c_tensor <- NULL
  P_tensor <- NULL
  pv <- NULL  ## parameter vector, lazily computed when has_params
  if (x@quad_obj) {
    if (has_params) {
      ## CVXPY v1.9.0 #3142: DPP tensor path for a PARAMETRIC quadratic
      ## objective. Extract P and the linear term as parameter tensors so a
      ## re-solve rebuilds them from new parameter values (no baking, no
      ## staleness). The current numeric P/q/d for THIS solve come from
      ## contracting the tensors with the current parameter vector.
      qt <- coeff_quad_form_tensor(extractor, .args(problem@objective)[[1L]])
      P_tensor <- qt$P_tensor
      c_tensor <- qt$c_tensor
      pv <- get_parameter_vector(extractor@param_to_size,
                                  extractor@param_id_map,
                                  parameters(problem))
      pv_sparse <- Matrix::Matrix(pv, ncol = 1L, sparse = TRUE)
      ## Reshape the contracted vec(P) to (x_length, x_length) without
      ## densifying (-> general dgCMatrix; see .dpp_contract_reshape).
      P_mat <- .dpp_contract_reshape(P_tensor, pv_sparse,
                                     extractor@x_length, extractor@x_length)
      c_flat <- as.numeric(c_tensor %*% pv)
      c_vec <- c_flat[seq_len(extractor@x_length)]
      d_offset <- c_flat[extractor@x_length + 1L]
    } else {
      ## QP path: extract 0.5*x'*P*x + q'*x + d (numeric, P baked)
      quad_result <- coeff_quad_form(extractor, .args(problem@objective)[[1L]])
      P_mat <- quad_result$P          # already 2x scaled for solver
      c_vec <- quad_result$q          # linear term
      d_offset <- quad_result$offset  # constant
    }
  } else {
    if (has_params) {
      ## DPP tensor path: extract c_tensor (x_length+1, param_size)
      c_tensor <- coeff_affine_tensor(extractor, list(.args(problem@objective)[[1L]]))
      ## Also get current numeric values via apply
      pv <- get_parameter_vector(extractor@param_to_size,
                                  extractor@param_id_map,
                                  parameters(problem))
      c_flat <- as.numeric(c_tensor %*% pv)
      c_vec <- c_flat[seq_len(extractor@x_length)]
      d_offset <- c_flat[extractor@x_length + 1L]
    } else {
      ## LP/conic path: extract c'*x + d
      obj_result <- coeff_affine(extractor, list(.args(problem@objective)[[1L]]))
      c_vec <- as.numeric(obj_result$A)  # x_length vector
      d_offset <- as.numeric(obj_result$b)  # scalar offset
    }
  }

  ## Step 5: Batch constraint args and extract A, b
  expr_chunks <- vector("list", length(ordered_cons))
  for (i in seq_along(ordered_cons)) {
    expr_chunks[[i]] <- .args(ordered_cons[[i]])
  }
  expr_list <- unlist(expr_chunks, recursive = FALSE)
  if (is.null(expr_list)) expr_list <- list()

  A_tensor <- NULL
  if (length(expr_list) > 0L) {
    if (has_params) {
      ## DPP tensor path: extract A_tensor
      A_tensor <- coeff_affine_tensor(extractor, expr_list)
      ## Also get current numeric values via apply
      if (is.null(pv)) {
        pv <- get_parameter_vector(extractor@param_to_size,
                                    extractor@param_id_map,
                                    parameters(problem))
      }
      ## Reshape the contracted column-major vec to (n_rows, x_length + 1)
      ## without densifying; first x_length columns = A, last column = b.
      pv_sparse <- Matrix::Matrix(pv, ncol = 1L, sparse = TRUE)
      n_cols <- extractor@x_length + 1L
      n_rows <- nrow(A_tensor) %/% n_cols
      Ab <- .dpp_contract_reshape(A_tensor, pv_sparse, n_rows, n_cols)
      A <- Ab[, seq_len(extractor@x_length), drop = FALSE]
      b <- as.numeric(Ab[, n_cols])
    } else {
      con_result <- coeff_affine(extractor, expr_list)
      A <- con_result$A
      b <- con_result$b
    }
  } else {
    A <- Matrix::sparseMatrix(i = integer(0), j = integer(0), x = numeric(0),
                              dims = c(0L, extractor@x_length))
    b <- numeric(0)
  }

  ## Step 6: Build ConeDims
  cone_dims <- ConeDims(constr_map)

  ## Step 7: Store extra data for inversion
  inverse_data@.extra <- new.env(hash = TRUE, parent = emptyenv())
  inverse_data@.extra$minimize <- TRUE
  inverse_data@.extra$constraints <- ordered_cons
  inverse_data@.extra$id2cons <- id2cons

  ## Native variable bounds. CVXPY builds bound tensors when bounds depend on
  ## Parameters; CVXR mirrors that for 2D scalar-or-exact-shape bounds.
  vars_ <- variables(problem)
  if (.has_parametric_bounds(vars_)) {
    if (is.null(pv)) {
      pv <- get_parameter_vector(extractor@param_to_size,
                                  extractor@param_id_map,
                                  parameters(problem))
    }
    lb_tensor <- .extract_bounds_tensor(
      vars_, extractor@x_length, extractor@param_to_size,
      extractor@param_id_map, which = "lower"
    )
    ub_tensor <- .extract_bounds_tensor(
      vars_, extractor@x_length, extractor@param_to_size,
      extractor@param_id_map, which = "upper"
    )
    lower_bounds <- as.numeric(lb_tensor %*% pv)
    upper_bounds <- as.numeric(ub_tensor %*% pv)
  } else {
    lb_tensor <- NULL
    ub_tensor <- NULL
    lower_bounds <- .extract_lower_bounds(vars_, extractor@x_length)
    upper_bounds <- .extract_upper_bounds(vars_, extractor@x_length)
  }

  ## Return data as named list + inverse data
  data <- list()
  data[[SD_C]]      <- c_vec
  data[[SD_OFFSET]]  <- d_offset
  data[[SD_A]]      <- A
  data[[SD_B]]      <- b
  data[[SD_DIMS]]   <- cone_dims
  data[[LOWER_BOUNDS]] <- lower_bounds
  data[[UPPER_BOUNDS]] <- upper_bounds
  data[["x_id"]]    <- .id(x_var)
  data[["constraints"]] <- ordered_cons
  if (!is.null(P_mat)) data[[SD_P]] <- P_mat

  ## Step 8: Build ParamConeProg for DPP caching
  if (has_params && !is.null(c_tensor)) {
    if (is.null(A_tensor)) {
      A_tensor <- Matrix::sparseMatrix(
        i = integer(0), j = integer(0), x = numeric(0),
        dims = c(0L, sum(as.integer(unlist(extractor@param_to_size))))
      )
    }
    data[[SD_PARAM_PROB]] <- ParamConeProg(
      c_tensor = c_tensor,
      A_tensor = A_tensor,
      x_length = extractor@x_length,
      x_id = .id(x_var),
      parameters = parameters(problem),
      param_id_to_col = extractor@param_id_map,
      param_to_size = extractor@param_to_size,
      variables = variables(problem),
      var_id_to_col = inverse_data@var_offsets,
      constraints = ordered_cons,
      cone_dims = cone_dims,
      P_tensor = P_tensor,
      lower_bounds = data[[LOWER_BOUNDS]],
      upper_bounds = data[[UPPER_BOUNDS]],
      lb_tensor = lb_tensor,
      ub_tensor = ub_tensor
    )
  }

  ## Extract MIP indices for boolean/integer variables
  mip_idx <- .extract_mip_idx(problem, inverse_data)
  data[[SD_BOOL_IDX]] <- mip_idx$bool_idx
  data[[SD_INT_IDX]]  <- mip_idx$int_idx

  list(data, inverse_data)
}

## invert: split solution vector back into per-variable values
## CVXPY SOURCE: cone_matrix_stuffing.py lines 432-470
method(reduction_invert, ConeMatrixStuffing) <- function(x, solution, inverse_data, ...) {
  var_map <- inverse_data@var_offsets

  ## Flip sign of opt val if maximize
  opt_val <- solution@opt_val
  if (!(solution@status %in% ERROR_STATUS) &&
      !is.null(inverse_data@.extra$minimize) &&
      !inverse_data@.extra$minimize) {
    opt_val <- -opt_val
  }

  primal_vars <- list()
  dual_vars <- list()

  ## Remap dual variables. CVXPY v1.9.0 fix: #3197 -- hoisted above the
  ## failure short-circuit so an infeasibility-certificate dual (the dual ray)
  ## is remapped and returned for INF_OR_UNB solutions, not just optimal ones.
  if (length(solution@dual_vars) > 0L) {
    con_map_env <- inverse_data@cons_id_map
    id2cons <- inverse_data@.extra$id2cons
    old_ids <- ls(con_map_env, all.names = TRUE)
    ## One vectorized lookup instead of n linear scans of `names(dual_vars)`,
    ## and a preallocated result instead of one grown by name: O(n^2) -> O(n),
    ## measured 9.1 ms -> 3.4 ms at n = 1000. This is `.remap_by_id_map`
    ## (zzz_R_specific/utility.R) with the per-constraint reshape below kept
    ## inline, since it is what makes the dual's SHAPE contract true.
    ## See notes/string_key_hashing_sweep_2026-08-13.md and ADR D_PERF.7.
    new_ids <- as.character(unlist(mget(old_ids, envir = con_map_env),
                                   use.names = FALSE))
    dv_idx <- match(new_ids, names(solution@dual_vars))
    n_old <- length(old_ids)
    out_dv <- vector("list", n_old)
    for (i in seq_len(n_old)) {
      if (is.na(dv_idx[i])) next
      old_id <- old_ids[[i]]
      dv <- solution@dual_vars[[dv_idx[i]]]
      if (is.null(dv)) next

      ## CVXPY SOURCE: cone_matrix_stuffing.py:467-478 -- give the dual the
      ## SHAPE OF ITS CONSTRAINT, in F (column-major) order, which R's `matrix`
      ## is natively.  The exemption list transliterates verbatim; the
      ## `PSD && num_cones() > 1` clause is written out even though
      ## `num_cones(PSD)` is 1L (psd.R:67) so it can never fire, because the
      ## point is to read as CVXPY's does.
      ##
      ## This is what makes the roxygen at constraints/constraint.R:89-91 ("a
      ## numeric matrix ... or a list of numeric matrices") true.  Its absence
      ## was the root cause of both dual bugs fixed earlier on this branch:
      ## with no reshape here, every per-cone `save_dual_value` had to know the
      ## raw stuffed layout itself, and two of them got it wrong.
      con_obj <- get0(old_id, envir = id2cons, ifnotfound = NULL)
      shape <- if (is.null(con_obj)) NULL else .shape(con_obj)
      exempt <- is.null(shape) || prod(shape) == 1L ||
        .s7_is(con_obj, ExpCone) || .s7_is(con_obj, SOC) ||
        (.s7_is(con_obj, PSD) && num_cones(con_obj) > 1L)

      if (exempt) {
        out_dv[i] <- list(dv)
      } else if (length(dv) != prod(shape)) {
        ## numpy raises here; R would RECYCLE silently whenever one length
        ## divides the other, which is exactly the class of corruption this
        ## whole reshape exists to prevent.
        cli_abort(c(
          "Dual for constraint {old_id} has length {length(dv)}, but its shape is {.val {shape}}.",
          "i" = "This is an internal inconsistency between the solver's dual layout and the constraint."
        ))
      } else {
        out_dv[i] <- list(matrix(dv, nrow = shape[1L], ncol = shape[2L]))
      }
    }
    names(out_dv) <- old_ids
    filled <- !vapply(out_dv, is.null, logical(1))
    dual_vars <- out_dv[filled]
  }

  if (!(solution@status %in% SOLUTION_PRESENT)) {
    return(Solution(solution@status, opt_val, primal_vars, dual_vars,
                    solution@attr))
  }

  ## Split vectorized variable into components
  x_opt <- solution@primal_vars[[1L]]
  if (is.null(x_opt)) {
    ## Try to get the single primal var
    if (length(solution@primal_vars) > 0L) {
      x_opt <- solution@primal_vars[[1L]]
    }
  }

  if (!is.null(x_opt)) {
    for (var_id in names(var_map)) {
      offset <- var_map[[var_id]]
      shape <- inverse_data@var_shapes[[var_id]]
      sz <- prod(shape)
      val <- x_opt[(offset + 1L):(offset + sz)]
      primal_vars[[var_id]] <- matrix(val, nrow = shape[1L], ncol = shape[2L])
    }
  }

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