#' Subspace distance
#'
#' @param As,Bs Basis matrices or list of Kronecker components of basis matrices
#' as representations of elements of the Grassmann manifold.
#' @param is.ortho Boolean to specify if \eqn{A} and \eqn{B} are semi-orthogonal.
#'   If false, the projection matrices are computed as
#' \deqn{P_A = A (A' A)^{-1} A'}
#'  otherwise just \eqn{P_A = A A'} since \eqn{A' A} is the identity.
#' @param normalize Boolean to specify if the distance shall be normalized.
#'  Meaning, the maximal distance scaled to be \eqn{1} independent of dimensions.
#'
#' @seealso
#'  K. Ye and L.-H. Lim (2016) "Schubert varieties and distances between
#'  subspaces of different dimensions" <arXiv:1407.0900>
#'
#' @export
dist.subspace <- function (As, Bs, is.ortho = FALSE, normalize = FALSE,
    tol = sqrt(.Machine$double.eps)
) {
    As <- if (is.list(As)) Map(as.matrix, As) else list(as.matrix(As))
    Bs <- if (is.list(Bs)) Map(as.matrix, Bs) else list(as.matrix(Bs))

    if (!is.ortho) {
        As <- Map(function(A) {
            qrA <- qr(A, tol)
            if (qrA$rank < ncol(A)) {
                qr.Q(qrA)[, abs(diag(qr.R(qrA))) > tol, drop = FALSE]
            } else {
                qr.Q(qrA)
            }
        }, As)
        Bs <- Map(function(B) {
            qrB <- qr(B, tol)
            if (qrB$rank < ncol(B)) {
                qr.Q(qrB)[, abs(diag(qr.R(qrB))) > tol, drop = FALSE]
            } else {
                qr.Q(qrB)
            }
        }, Bs)
    }

    rankA <- prod(sapply(As, ncol))
    rankB <- prod(sapply(Bs, ncol))

    c <- if (normalize) {
        rankSum <- rankA + rankB
        sqrt(1 / max(1, min(rankSum, 2 * prod(sapply(As, nrow)) - rankSum)))
    } else {
        1
    }

    s <- prod(mapply(function(A, B) sum(crossprod(A, B)^2), As, Bs))
    c * sqrt(max(0, rankA + rankB - 2 * s))
}