#' @rdname matProj
#' @export
projSym <- function(A) 0.5 * (A + t(A))
#' @rdname matProj
#' @export
projDiag <- function(A) diag(diag(A))
#' @rdname matProj
#' @export
.projBand <- function(dims, low, high) {
    diag.index <- .row(dims) - .col(dims)
    mask <- (diag.index <= low) & (-high <= diag.index)
    function(A) A * mask
}
#' @rdname matProj
#' @export
.projSymBand <- function(dims, low, high) {
    diag.index <- .row(dims) - .col(dims)
    mask <- (diag.index <= low) & (-high <= diag.index)
    function(A) projSym(A) * mask
}
#' @rdname matProj
#' @export
.projPSD <- function(sym = FALSE) {
    if (sym) {
        function(A) {
            eig <- eigen(A, symmetric = TRUE)
            eig$vectors %*% (pmax(0, eig$values) * t(eig$vectors))
        }
    } else {
        function(A) {
            eig <- eigen(0.5 * (A + t(A)), symmetric = TRUE)
            eig$vectors %*% (pmax(0, eig$values) * t(eig$vectors))
        }
    }
}
#' @rdname matProj
#' @export
.projRank <- function(rank) {
    force(rank)
    function(A) {
        rank <- min(dim(A), rank)
        svdA <- La.svd(A, rank, rank)
        svdA$u %*% (svdA$d[seq_len(rank)] * svdA$vt)
    }
}
#' @rdname matProj
#' @export
.projSymRank <- function(rank) {
    force(rank)
    function(A) {
        rank <- min(dim(A), rank)
        svdA <- La.svd(0.5 * (A + t(A)), rank, rank)
        svdA$u %*% (svdA$d[seq_len(rank)] * svdA$vt)
    }
}
#' @rdname matProj
#' @export
projStiefel <- function(A) {
    # Using a polar decomposition of `A = Q P` via SVD `A = U D V^T`. Compaired
    # to a QR decomposition the polar decomposition is unique, making it "stabel".
    svdA <- La.svd(A)
    svdA$u %*% svdA$vt  # = Q
}
# .projKron <- function(dims) {
#     ... # TODO: Implement this!
# }

#' @rdname matProj
#' @export
.projMaskedMean <- function(mask) {
    force(mask)
    function(A) {
        `[<-`(matrix(0, nrow(A), ncol(A)), mask, mean(A[mask]))
    }
}