#' @examples
#' nrows <- c(3, 2, 5)
#' ncols <- c(2, 4, 8)
#'
#' A <- Reduce(kronecker, Map(matrix, Map(rnorm, nrows * ncols), nrows))
#'
#' all.equal(A, Reduce(kronecker, approx.kron(A, nrows, ncols)))
#'
#' @export
approx.kron <- function(A, nrows, ncols) {

    # rearrange kronecker product `A` into outer product `R`
    dim(A) <- c(rev(nrows), rev(ncols))
    axis.perm <- as.vector(t(matrix(seq_along(dim(A)), ncol = 2L))[, rev(seq_along(nrows))])
    R <- aperm(A, axis.perm, resize = FALSE)
    dim(R) <- nrows * ncols

    # scaling factor for every product component
    s <- sum(A^2)^(1 / length(dim(A)))

    # higher order SVD on `R` with one singular vector
    Map(function(mode, nr, nc) {
        s * `dim<-`(La.svd(mcrossprod(R, mode = mode), 1L, 0L)$u, c(nr, nc))
    }, seq_along(nrows), nrows, ncols)

}


#' Approximates Kronecker Product decomposition.
#'
#' Approximates the matrices `A` and `B` such that
#'      C = A %x% B
#' with `%x%` the kronecker product of the matrixes `A` and `B`
#' of dimensions `dimA` and `dimB` respectively.
#'
#' @param C desired kronecker product result.
#' @param dimA length 2 vector of dimensions of \code{A}.
#' @param dimB length 2 vector of dimensions of \code{B}.
#'
#' @return list with attributes `A` and `B`.
#'
#' @examples
#' A <- matrix(seq(14), 7, 2)
#' B <- matrix(c(1, 0), 3, 4)
#' C <- kronecker(A, B) # the same as 'C <- A %x% B'
#' approx.kronecker(C, dim(A), dim(B))
#'
#' @seealso C.F. Van Loan / Journal of Computational and Applied Mathematics
#'          123 (2000) 85-100 (pp. 93-95)
#'
#' @export
approx.kronecker <- function(C, dimA, dimB = dim(C) / dimA) {

    dim(C) <- c(dimB[1L], dimA[1L], dimB[2L], dimA[2L])
    R <- aperm(C, c(2L, 4L, 1L, 3L))
    dim(R) <- c(prod(dimA), prod(dimB))

    if (requireNamespace("RSpectra", quietly = TRUE)) {
        svdR <- RSpectra::svds(R, 1L)
    } else {
        svdR <- svd(R, 1L, 1L)
    }

    list(
        A = array(sqrt(svdR$d[1]) * svdR$u, dimA),
        B = array(sqrt(svdR$d[1]) * svdR$v, dimB)
    )
}

#' Kronecker Product Decomposition.
#'
#' Computes the components summation components `A_i`, `B_i` of a sum of
#' Kronecker products
#'      C = sum_i A_i %x% B_i
#' with the minimal estimated number of summands.
#'
#' @param C desired kronecker product result.
#' @param dimA length 2 vector of dimensions of \code{A}.
#' @param dimB length 2 vector of dimensions of \code{B}.
#' @param tol tolerance of singular values of \code{C} to determin the number of
#'  needed summands.
#'
#' @return list of lenghth with estimated number of summation components, each
#'  entry consists of a list with named entries \code{"A"} and \code{"B"} of
#'  dimensions \code{dimA} and \code{dimB}.
#'
#' @examples
#' As <- replicate(3, matrix(rnorm(2 * 7), 2), simplify = FALSE)
#' Bs <- replicate(3, matrix(rnorm(5 * 3), 5), simplify = FALSE)
#' C <- Reduce(`+`, Map(kronecker, As, Bs))
#'
#' decomposition <- decompose.kronecker(C, c(2, 7))
#'
#' reconstruction <- Reduce(`+`, Map(function(summand) {
#'     kronecker(summand[[1]], summand[[2]])
#' }, decomposition), array(0, dim(C)))
#'
#' stopifnot(all.equal(C, reconstruction))
#'
#' @export
decompose.kronecker <- function(C, dimA, dimB = dim(C) / dimA, tol = 1e-7) {
    # convert the equivalent outer product
    dim(C) <- c(dimB[1L], dimA[1L], dimB[2L], dimA[2L])
    C <- aperm(C, c(2L, 4L, 1L, 3L), resize = FALSE)
    dim(C) <- c(prod(dimA), prod(dimB))

    # Singular Valued Decomposition
    svdC <- La.svd(C)

    # Sum of Kronecker Components
    lapply(seq_len(sum(svdC$d > tol)), function(i) list(
        A = matrix(svdC$d[i] * svdC$u[, i], dimA),
        B = matrix(svdC$vt[i, ], dimB)
    ))
}

### Given that C is a Kronecker product this is a fast method but a bit
### unreliable in full generality.
# decompose.kronecker <- function(C, dimA, dimB = dim(C) / dimA) {
#     dim(C) <- c(dimB[1L], dimA[1L], dimB[2L], dimA[2L])
#     R <- aperm(C, c(2L, 4L, 1L, 3L))
#     dim(R) <- c(prod(dimA), prod(dimB))
#     max.index <- which.max(abs(R))
#     row.index <- ((max.index - 1L) %% nrow(R)) + 1L
#     col.index <- ((max.index - 1L) %/% nrow(R)) + 1L
#     max.elem <- if (abs(R[max.index]) > .Machine$double.eps) R[max.index] else 1
#     list(
#         A = array(R[, col.index], dimA),
#         B = array(R[row.index, ] / max.elem, dimB)
#     )
# }


# kron <- function(A, B) {
#     perm <- as.vector(t(matrix(seq_len(2 * length(dim(A))), ncol = 2)[, 2:1]))
#     K <- aperm(outer(A, B), perm)
#     dim(K) <- dim(A) * dim(B)
#     K
# }

# kronperm <- function(A) {
#     # force A to have even number of dimensions
#     dim(A) <- c(dim(A), rep(1L, length(dim(A)) %% 2L))
#     # compute axis permutation
#     perm <- as.vector(t(matrix(seq_along(dim(A)), ncol = 2)[, 2:1]))
#     # permute elements of A
#     K <- aperm(A, perm, resize = FALSE)
#     # collapse/set dimensions
#     dim(K) <- head(dim(A), length(dim(A)) / 2) * tail(dim(A), length(dim(A)) / 2)
#     K
# }

# p <- c(2, 3, 5)
# q <- c(3, 4, 7)
# A <- array(rnorm(prod(p)), p)
# B <- array(rnorm(prod(q)), q)
# all.equal(kronperm(outer(A, B)), kronecker(A, B))
# all.equal(kron(A, B), kronecker(A, B))


# dA <- c(2, 3, 5)
# dB <- c(3, 4, 7)
# A <- array(rnorm(prod(dA)), dA)
# B <- array(rnorm(prod(dB)), dB)
# X <- outer(A, B)

# r <- length(dim(X)) / 2
# I <- t(do.call(expand.grid, Map(seq_len, head(dim(X), r) * tail(dim(X), r))))
# K <- apply(rbind(
#     (I - 1) %/% tail(dim(X), r) + 1,
#     (I - 1) %% tail(dim(X), r) + 1
# ), 2, function(i) X[sum(c(1, cumprod(head(dim(X), -1))) * (i - 1)) + 1])
# dim(K) <- head(dim(X), r) * tail(dim(X), r)

# all.equal(kronecker(A, B), K)