#' Sliced Inverse Regression
#'
#' @export
SIR <- function(X, y, d, nr.slices = 10L, slice.method = c("cut", "ecdf")) {

    if (!(is.factor(y) || is.integer(y))) {
        slice.method <- match.arg(slice.method)
        if (slice.method == "ecdf") {
            y <- cut(ecdf(y)(y), nr.slices)
        } else {
            y <- cut(y, nr.slices)
            # ensure there are no empty slices
            if (any(table(y) == 0)) {
                y <- as.factor(as.integer(y))
            }
        }
    }

    # Center `X`
    Z <- scale(X, scale = FALSE)

    # Split `Z` into slices determined by `y`
    slices <- Map(function(i) Z[i, , drop = FALSE], split(seq_along(y), y))

    # Sizes and Means for each slice
    slice.sizes <- mapply(nrow, slices)
    slice.means <- Map(colMeans, slices)

    # Inbetween slice covariances
    sCov <- Reduce(`+`, Map(function(mean_s, n_s) {
        n_s * tcrossprod(mean_s)
    }, slice.means, slice.sizes)) / nrow(X)

    # Compute EDR directions
    La.svd(sCov, d, 0L)$u
}