## See hdrm R package, thank you to Patrick Breheny for providing useful simulation code

genS <- function (p, rho, corr) {
  if (corr == "exchangeable") {
    S <- matrix(rho, p, p) + (1 - rho) * diag(p)
  }
  else if (corr == "autoregressive") {
    RHO <- matrix(rho^(0:(p - 1)), p, p, byrow = TRUE)
    S <- Matrix::bandSparse(p, k = 0:(p - 1), diagonals = RHO, 
                            symmetric = TRUE)
  }
  as.matrix(S)
}

genX <- function (n, p, S) {
  R <- chol(S)
  as.matrix(matrix(rnorm(n * p), n, p) %*% R)
}

genY <- function (eta, family = c("gaussian", "binomial"), sigma = 1) {
  family = match.arg(family)
  n <- length(eta)
  if (family == "gaussian") 
    y <- rnorm(n, mean = eta, sd = sigma)
  else if (family == "binomial") {
    pi. <- exp(eta)/(1 + exp(eta))
    pi.[eta > log(0.9999/1e-04)] <- 1
    pi.[eta < log(1e-04/0.9999)] <- 0
    y <- rbinom(n, 1, pi.)
  }
  return(y)
}


## See hdrm R package
genData <- function (n, p, p1 = floor(p/2), beta, family = c("gaussian", "binomial"), 
                     SNR = 1, signal = c("homogeneous", "heterogeneous"), 
                     corr = c("exchangeable", "autoregressive"), rho = 0) {
  family <- match.arg(family)
  signal <- match.arg(signal)
  corr <- match.arg(corr)
  S <- genS(p, rho, corr)
  X <- genX(n, p, S)
  if (missing(beta) || length(beta) == 1) {
    j <- 1:p
    s <- c(-1, 1)[j%%2 + 1]
    b <- (j <= p1)
    if (missing(beta)) {
      if (signal == "heterogeneous") 
        b <- b * rev(j)
      b <- b * s
      beta <- b * sqrt(SNR)/sqrt(drop(crossprod(b, S) %*% b))
    }
    else {
      beta <- b * s * beta
    }
  }
  y <- genY(X %*% beta, family = family, sigma = 1)
  w <- 1 + floor(log10(p))
  vlab <- paste0("V", formatC(1:p, format = "d", width = w, flag = "0"))
  colnames(X) <- names(beta) <- vlab
  list(X = X, y = y, beta = beta, family = family)
}
