#' Random sampling from a tensor (multi-array) normal distribution
#'
#' @examples
#' n <- 1000
#' Sigma.1 <- 0.5^abs(outer(1:5, 1:5, "-"))
#' Sigma.2 <- diag(1:4)
#' X <- rtensornorm(n, 0, Sigma.1, Sigma.2)
#'
#' @export
rtensornorm <- function(n, mean, ..., sample.axis) {
    # get covariance matrices
    cov <- list(...)

    # get sample shape (dimensions of an observation)
    shape <- unlist(Map(nrow, cov))

    # result tensor dimensions
    dims <- c(shape, n)

    # validate mean dimensions
    if (!missing(mean)) {
        stopifnot(is.numeric(mean))
        stopifnot(is.array(mean) || length(mean) == 1)
        stopifnot(length(shape) == length(cov))
        stopifnot(all(shape == dim(mean)))
    }

    # sample i.i.d. normal entries
    X <- array(rnorm(prod(dims)), dim = dims)

    # transform from standard normal to tensor normal with given covariances
    for (i in seq_along(cov)) {
        X <- ttm(X, matpow(cov[[i]], 1 / 2), i)
    }

    # add mean (using recycling, observations on last mode)
    X <- X + mean

    # permute axis for indeing observations on sample mode (permute first axis
    # with sample mode axis)
    if (!missing(sample.axis)) {
        axis <- seq_len(length(dims) - 1)
        start <- seq_len(sample.axis - 1)
        end <- seq_len(length(dims) - sample.axis) + sample.axis - 1
        X <- aperm(X, c(axis[start], length(dims), axis[end]))
    }

    X
}