#' Generale a matrix of all permutations of `n` elements
permutations <- function(n) {
    if (n <= 0) {
        matrix(nrow = 0, ncol = 0)
    } else if (n == 1) {
        matrix(1)
    } else {
        sub.perm <- permutations(n - 1)
        p <- nrow(sub.perm)
        A <- matrix(NA, n * p, n)
        for (i in 1:n) {
            A[(i - 1) * p + 1:p, ] <- cbind(i, sub.perm + (sub.perm >= i))
        }
        A
    }
}


#' General symmetrization opperation for tensors (arrays) of equal dimensions
#'
#' @param A array of dimensions c(p, ..., p)
#'
#' @returns array of same dimensions as `A`
#'
#' @export
tsym <- function(A) {
    stopifnot(all(dim(A) == nrow(A)))

    if (is.matrix(A)) {
        return(0.5 * (A + t(A)))
    }

    axis.perm <- permutations(length(dim(A)))

    S <- array(0, dim(A))
    for (i in seq_len(nrow(axis.perm))) {
        S <- S + aperm(A, axis.perm[i, ])
    }

    S / nrow(axis.perm)
}



#' Genralized (pseudo) symmetrication for generel multi-dimensional arrays
sym <- function(A, FUN = `+`, scale = factorial(length(dim(A)))) {
    FUN <- match.fun(FUN)

    if (is.matrix(A) && (nrow(A) == ncol(A))) {
        A <- FUN(A, t(A))
        return(if (is.numeric(scale)) A / scale else A)
    }

    A.copy <- A
    perm <- seq_along(dim(A))
    while (length(pivot <- which(diff(perm) > 0))) {
        pivot <- max(pivot)
        successor <- max(which(perm[seq_along(perm) > pivot] > perm[pivot])) + pivot
        perm[c(pivot, successor)] <- perm[c(successor, pivot)]
        suffix <- seq(pivot + 1, length(perm))
        perm <- c(perm[-suffix], perm[rev(suffix)])

        modes <- which(perm != seq_along(perm))
        sub.dimA <- dim(A)
        sub.dimA[modes] <- min(dim(A)[modes])
        sub.indices <- Map(seq_len, sub.dimA)
        sub.selection <- do.call(`[`, c(list(A.copy), sub.indices, drop = FALSE))


        sub.assign <- do.call(call, c(list("[<-", quote(quote(A))), sub.indices,
            quote(quote(FUN(
                do.call(`[`, c(list(A), sub.indices)),
                aperm(do.call(`[`, c(list(A.copy), sub.indices, drop = FALSE)), perm)
            )))
        ))

        A <- eval(sub.assign)
    }

    if (is.numeric(scale)) A / scale else A
}