#' Stability Score
#'
#' Computing entry wise stability score for selection
#'
#' @param X The data matrix to operate on
#' @param K Number of k-folds (default 25)
#' @param subsample_fraction Hold out fraction (default 0.7)
#'
#' @return Score array
#' @export
stability_score <- function(X, K = 25, subsample_fraction = 0.7){
	n <- nrow(X)
	p <- ncol(X)
	rho <- vector('list', K)
	for(r in seq_len(K)){
		idx <- sample(n, size = floor(subsample_fraction * n))
		rho[[r]] <- cor(X[idx, ])
	}
	R <- simplify2array(rho)
	theta <- acos(pmax(pmin(R, 1), -1))
	delta <- theta - pi/2
	delta_sd <- apply(delta, c(1, 2), sd)
	delta_sd <- pmax(delta_sd, quantile(delta_sd[lower.tri(delta_sd)], 0.05))
	score <- 1 / delta_sd
	W <- score / median(score[lower.tri(score)])
	return(W[lower.tri(W)])
}