#' Projection Distance of two matrices
#'
#' Defined as sine of the maximum principal angle between the column spaces
#' of the matrices
#'      max{ sin theta_i, i = 1, ..., min(d1, d2) }
#' In case of rank(A) = rank(B) this measure is equivalent to
#'      || A A' - B B' ||_2
#' where ||.||_2 is the spectral norm (max singular value).
#'
#' @param A,B matrices of size \eqn{p\times d_1} and \eqn{p\times d_2}.
#'
#' @export
dist.projection <- function(A, B, is.ortho = FALSE,
    tol = sqrt(.Machine$double.eps)
) {
    if (!is.matrix(A)) A <- as.matrix(A)
    if (!is.matrix(B)) B <- as.matrix(B)

    if (!is.ortho) {
        qrA <- qr(A, tol)
        rankA <- qrA$rank
        A <- qr.Q(qrA)[, seq_len(qrA$rank), drop = FALSE]
        qrB <- qr(B, tol)
        rankB <- qrB$rank
        B <- qr.Q(qrB)[, seq_len(qrB$rank), drop = FALSE]
    } else {
        rankA <- ncol(A)
        rankB <- ncol(B)
    }

    if (rankA == 0 || rankB == 0) {
        return(as.double(rankA != rankB))
    } else if (rankA == 1 && rankB == 1) {
        sigma.min <- min(abs(sum(A * B)), 1)
    } else {
        sigma.min <- min(La.svd(crossprod(A, B), 0, 0)$d, 1)
    }
    if (sigma.min < 0.5) {
        sin(acos(sigma.min))
    } else {
        cos(asin(sigma.min))
    }
}