#' Simulate Multivariate Normal Data
#'
#' Given a set of parameters controlling the properties
#' of the true covariance matrix, simulate data
#' with the intention of evaluating the performance
#' of covariance or correlation matrix estimators.
#'
#' @param n Number of observations to generate.
#' @param p Number of features to simulate.
#' @param rho The AR(1) decay parameter (default 0.8). Setting rho = 0 uses the identity matrix as the true covariance (0 <= rho < 1).
#' @param contamination The proportions of rows to contaminate as outlier observations (default 0, 0 <= contamination <= 1).
#' @param block_fraction The fraction of the covariance matrix to form as AR(1) structure (default 0.5, 0 < block_fraction <= 1).
#' @param t_df The degrees of freedom to use for scaling observations into heavy tails (default Inf, 0 <= t_df < 2).
#'
#' @return A list with X, S and a vector identifying the contaminated rows.
#'
#' @examples
#' X_obj <- generate_data(n = 10, p = 10)
#' X <- X_obj$X
#' S_true <- X_obj$S
#'
#' @importFrom MASS mvrnorm
#'
#' @export
simulate_data <- function(
	n,
	p,
	rho = 0.8,
	contamination = 0,
	block_fraction = 0.5,
	t_df = Inf,
	exact_contamination = FALSE
){
	if(any(
		rho < 0,
		rho >= 1,
		contamination < 0,
		contamination > 1,
		block_fraction <= 0,
		block_fraction > 1,
		t_df <= 2
	)){
		stop('Incompatible settings')
	}

	# True correlation matrix
	Sigma <- diag(p)
	if(rho > 0){
		b <- floor(block_fraction * p)
		if(b >= 2){
			idx <- seq_len(b)
			Sigma[idx, idx] <- rho ^ abs(outer(idx, idx, `-`))
		}
	}

	# Regular clean observations
	Z <- MASS::mvrnorm(n = n, mu = numeric(p), Sigma = Sigma)
	if(is.finite(t_df)){
		w <- rchisq(n, df = t_df)
		X <- Z / sqrt(w / t_df)
	}else{
		X <- Z
	}

	# Contamination - by row
	outlier_rows <- integer(0)
	if(contamination > 0){
		if(exact_contamination){
			outlier_rows <- sample.int(n, ceiling(contamination * n))
		}else{
			outlier_rows <- which(runif(n) < contamination)
		}

		if(length(outlier_rows) > 0){
			X[outlier_rows, ] <- matrix(
				rchisq(length(outlier_rows) * p, df = 1),
				nrow = length(outlier_rows),
				ncol = p
			)
		}
	}
	list(
		X = X,
		S = Sigma,
		contaminated = sort(outlier_rows)
	)
}