#' Half vectorization of a matrix (lower part)
#' @export
vech <- function(A) A[lower.tri(A, diag = TRUE)]

#' @export
vech.index <- function(p) which(.row(c(p, p)) >= .col(c(p, p)))

#' @export
vech.pinv.index <- function(p) {
    index <- matrix(NA_integer_, p, p)
    index[lower.tri(index, diag = TRUE)] <- seq_len(p * (p + 1L) / 2L)
    index[upper.tri(index)] <- t(index)[upper.tri(index)]

    index
}

#' pseudo inserse of the half vectorization
#'
#' @examples
#' # only valid for symmetric matrices
#' A <- matrix(rnorm(4^2), 4)
#' A <- A + t(A)
#' all.equal(A, vech.pinv(vech(A)))
#'
#' @export
vech.pinv <- function(a) {
    # determin original dimensions
    p <- as.integer((sqrt(8 * length(a) + 1) - 1) / 2)
    stopifnot(p * (p + 1L) == 2L * length(a))

    # de-vectorized matrix
    matrix(a[vech.pinv.index(p)], p, p)
}


# vech <- function(A) {
#     stopifnot(all(nrow(A) == dim(A)))

#     axis <- seq_len(nrow(A))
#     grid <- expand.grid(c(list(rev(axis)), rep(list(axis), length(dim(A)) - 1)))

#     A[rowSums(grid - 1) < nrow(A)]
# }

# p <- 4
# X <- matrix(rnorm(p^2), p)

# vech(outer(X, X))

# vech(outer(c(X), c(X)))

# sort(unique(c(sym(outer(X, X)))))
# sort(vech(sym(outer(X, X))))

# # (I <- matrix(c(
# #     1, 1, 1,
# #     2, 1, 1,
# #     3, 1, 1,
# #     2, 2, 1,
# #     3, 2, 1,
# #     3, 3, 1,
# #     2, 1, 2,
# #     3, 1, 2,
# #     3, 2, 2,
# #     3, 1, 3
# # ), ncol = 3, byrow = TRUE))

# # ((I - 1) %*% nrow(A)^(seq_along(dim(A)) - 1) + 1)

# # p <- 3
# # ord <- 4
# # A <- array(seq_len(p^ord), rep(p, ord))

# # axis <- seq_len(nrow(A))
# # grid <- expand.grid(c(list(rev(axis)), rep(list(axis), length(dim(A)) - 1)))
# # array(rowSums(grid - 1), dim(A))

# # A[rowSums(grid - 1) < nrow(A)]



# # apply(indices, 1, function(i) do.call(`[`, c(list(A), as.list(i))))