#' A simple higher order (multi-way) canonical correlation analysis.
#'
#' @param X multi-dimensional array
#' @param Y multi-dimensional array with the same nr. of dimensions and equal
#'  sample axis to `X`.
#' @param sample.axis integer indicationg which axis enumerates observations
#'
#' @export
HOCCA <- function(X, Y, sample.axis = length(dim(X)), centerX = TRUE, centerY = TRUE) {

    # ensure sample axis is the last axis
    if (!missing(sample.axis)) {
        modes <- seq_along(dim(X))[-sample.axis]
        X <- aperm(X, c(modes, sample.axis))
        Y <- aperm(Y, c(modes, sample.axis))
    }
    modes <- seq_len(length(dim(X)) - 1L)
    dimX <- head(dim(X), -1L)
    dimF <- head(dim(F), -1L)
    sample.size <- tail(dim(X), 1L)

    # center `X` and `Y`
    if (centerX) {
        X <- X - as.vector(rowMeans(X, dims = length(dim(X)) - 1L))
    }
    if (centerY) {
        Y <- Y - as.vector(rowMeans(Y, dims = length(dim(Y)) - 1L))
    }

    # estimate marginal covariance matrices
    CovXX <- Map(function(mode) mcrossprod(X, mode = mode) / prod(dim(X)[-mode]), modes)
    CovYY <- Map(function(mode) mcrossprod(Y, mode = mode) / prod(dim(Y)[-mode]), modes)
    # and the "covariance tensor"
    CovXY <- array(tcrossprod(mat(X, modes), mat(Y, modes)) / sample.size, c(dimX, dimF))

    # Compute standardized X and Y correlation tensor
    SCovXY <- mlm(CovXY, Map(matpow, c(CovXX, CovYY), -1 / 2))

    # mode-wise canonical correlation directions
    hosvd <- HOSVD(SCovXY, nu = rep(pmin(dimX, dimF), 2L))
    dirsX <- hosvd$Us[modes]
    dirsY <- hosvd$Us[modes + length(modes)]

    list(dirsX = dirsX, dirsY = dirsY)
}


# c(dirsX, dirsY) %<-% HOCCA(X, F)
# B.hocca <- Reduce(kronecker, rev(Map(tcrossprod, dirsX, dirsY)))
# dist.subspace(B.true, B.hocca, normalize = TRUE)
# dist.subspace(B.true, Reduce(kronecker, rev(dirsX)), normalize = TRUE)

# cca <- cancor(mat(X, 4), mat(F, 4))
# B.cca <- tcrossprod(cca$xcoef[, prod(dim(X)[-4])], cca$xcoef[, prod(dim(F)[-4])])
# dist.subspace(B.true, cca$xcoef[, prod(dim(X)[-4])], normalize = TRUE)