#' Mode wise Covariance Estimates
#'
#' Estimates Covariances \eqn{\Sigma_k}{Sigma_k} for each mode \eqn{k}.
#' This is equivalent to assuming a Kronecker structured Covariance
#'
#' \deqn{\Sigma = \Sigma_r\otimes ... \otimes\Sigma_1}{%
#'       Sigma = Sigma_r %x% ... %x% Sigma_1}
#'
#' where \eqn{\Sigma}{Sigma} is the Covariance \eqn{cov(vec(X))} of the
#' vectorized variables. This function estimates the Kronecerk components
#' \eqn{\Sigma_k}{Sigma_k}.
#'
#' @param X multi-dimensional array
#' @param sample.axis observation axis index
#'
#' @export
mcov <- function(X, sample.axis = 1L) {
    # observation modes (axis indices)
    modes <- seq_along(dim(X))[-sample.axis]
    # observation dimensions
    p <- dim(X)[modes]
    # observation tensor order
    r <- length(p)

    # ensure observations are on the last mode
    if (sample.axis != r + 1L) {
        X <- aperm(X, c(modes, sample.axis))
    }
    # centering: Z = X - E[X]
    Z <- X - c(rowMeans(X, dims = r))

    # estimes (unscaled) covariances for each mode
    Sigmas <- .mapply(mcrossprod, list(mode = seq_len(r)), MoreArgs = list(Z))
    # scale by per mode "sample" size
    Sigmas <- .mapply(`*`, list(Sigmas, p / prod(dim(X))), NULL)

    # estimate trace of Kronecker product of covariances
    tr.est <- prod(p) * mean(Z^2)

    # as well as the current trace of the unscaled covariances
    tr.Sigmas <- prod(unlist(.mapply(function(S) sum(diag(S)), list(Sigmas), NULL)))

    # Scale each mode Covariance to match the estimated Kronecker product scale
    .mapply(`*`, list(Sigmas), MoreArgs = list((tr.est / tr.Sigmas)^(1 / r)))
}