#' Per mode (axis) MLE
#'
#' @export
kpir.mle <- function(X, Fy, max.iter = 500L, sample.axis = 1L,
    logger = NULL #, eps = .Machine$double.eps
) {
    ### Step 0: Setup/Initialization
    if (!is.array(Fy)) {
        # scalar response case (add new axis of size 1)
        dim(Fy) <- ifelse(seq_along(dim(X)) == sample.axis, dim(X)[sample.axis], 1L)
    }
    # Check dimensions and matching of axis (tensor order)
    stopifnot(exprs = {
        length(dim(X)) == length(dim(Fy))
        dim(X)[sample.axis] == dim(Fy)[sample.axis]
    })
    # warn about model constraints
    if (any(dim(Fy)[-sample.axis] >= dim(X)[-sample.axis])) {
        warning("Degenerate case 'any(dim(Fy)[-sample.axis] >= dim(X)[-sample.axis])'")
    }

    # extract dimensions (for easier handling as local variables)
    modes <- seq_along(dim(X))[-sample.axis]    # predictor axis indices
    n <- dim(X)[sample.axis]                    # observation count (scalar)
    p <- dim(X)[-sample.axis]                   # predictor dimensions (vector)
    # q <- dim(Fy)[-sample.axis]                  # response dimensions (vector)
    # r <- length(dim(X)) - 1L                    # tensor order (scalar)

    # Means for X and Fy (a.k.a. sum elements over the sample axis)
    meanX <- apply(X, modes, mean, simplify = TRUE)
    meanFy <- apply(Fy, modes, mean, simplify = TRUE)
    # Center both X and Fy
    X <- sweep(X, modes, meanX)
    Fy <- sweep(Fy, modes, meanFy)


    ### Step 1: Initial value estimation
    alphas <- Map(function(mode, ncol) {
        La.svd(mcrossprod(X, mode = mode), ncol)$u
    }, modes, dim(Fy)[modes])
    # Residuals
    R <- X - mlm(Fy, alphas, modes = modes)
    # Covariance Moment estimates
    Deltas <- Map(mcrossprod, list(R), mode = modes)
    Deltas <- Map(function(Delta, j) (n * prod(p[-j]))^(-1) * Delta,
        Deltas, seq_along(Deltas))

    # Call history callback (logger) before the first iteration
    if (is.function(logger)) { do.call(logger, c(0L, NA, alphas, Deltas)) }


    ### Step 2: Alternating estimate updates
    for (iter in seq_len(max.iter)) {
        # Compute covariance inverses
        Deltas.inv <- Map(solve, Deltas)

        # "standardize" X
        Z <- mlm(X, Deltas.inv, modes = modes)

        # Compute new alpha estimates
        alphas <- Map(function(j) {
            # MLE estimate for alpha_j | alpha_k, Delta_l for all k != j and l
            FF <- mlm(Fy, alphas[-j], modes = modes[-j])
            Deltas[[j]] %*% t(solve(
                t(mcrossprod(mlm(FF, Deltas.inv[-j], modes = modes[-j]), FF, mode = modes[j])),
                t(mcrossprod(Z, FF, mode = modes[j]))))
        }, seq_along(modes))

        # update residuals
        R <- X - mlm(Fy, alphas, modes = modes)

        # next Delta estimates
        Deltas <- Map(function(j) {
            # MLE estimate for Delta_j | Delta_k, alpha_l for all k != j and l
            (n * prod(p[-j]))^(-1) * mcrossprod(
                mlm(R, Deltas[-j], modes = modes[-j]), R, mode = modes[j])
        }, seq_along(modes))


        # TODO: break condition!!!


        # Call history callback
        if (is.function(logger)) { do.call(logger, c(iter, NA, alphas, Deltas)) }
    }

    list(alphas = structure(alphas, names = as.character(modes)),
         Deltas = structure(Deltas, names = as.character(modes)),
         meanX = meanX, meanFy = meanFy)
}