

benchmark_blasso_bayeslm <- function(
    vy,
    mX,
    lambda_init,
    sigma2_init,
    a, b, u, v,
    nburn, nsamples,
    trials,
    beta_inds=NA)
{

  if (!requireNamespace("bayeslm", quietly = TRUE)) {
    stop(
      "The 'bayeslm' package is required for this benchmark.\n",
      "Install it with install.packages('bayeslm').",
      call. = FALSE
    )
  }

  if (!is.numeric(nsamples) || !is.numeric(nburn) || length(nsamples) != 1L || length(nburn) != 1L) {
    stop("'nsamples' and 'nburn' must be numeric scalars.", call. = FALSE)
  }
  if (nsamples <= nburn) {
    stop("'nsamples' must be greater than 'nburn'.", call. = FALSE)
  }
  if (!is.numeric(trials) || length(trials) != 1L || trials < 1) {
    stop("'trials' must be a positive scalar.", call. = FALSE)
  }

  mX <- as.matrix(mX)
  p <- ncol(mX)

  block_vec = rep(1, p)

  # Run the Gibbs sampler trials times
  mStat = c()
  for (i in seq_len(trials)) {
    time_val = system.time({
      res_mcmc = bayeslm::bayeslm(Y=vy, X=mX, prior = 'laplace', icept = FALSE,
                         block_vec = block_vec, N = nsamples, burnin=nburn,
                         singular=TRUE)
    })[3]
    # print(time_val)

    inds_use = 1:(nsamples - nburn)

    # Calculate summary statistics of efficiencies and mixing rates
    stats = FastGibbsSamplers:::mcmc_stats(res_mcmc$beta, res_mcmc$sigma, res_mcmc$vglobal, time_val, inds_use)
    # print(stats)

    mStat = rbind(mStat,stats)
  }

  colname_vals = c("eff_beta", "mix_beta", "eff_sigma2", "mix_sigma2", "eff_lambda2", "mix_lambda2", "time")

  #rownames(mStat) <- NA
  colnames(mStat) <- colname_vals


  # Optional plotting
  if (isTRUE(plot)) {
    # Plot the acf for sigma2 and lambda2
    acf(res_mcmc_last$sigma[inds_use])
    acf(res_mcmc_last$vglobal[inds_use])

    # Trace plots for sigma2 and lambda2
    plot(res_mcmc_last$sigma[inds_use], type = "l")
    plot(res_mcmc_last$vglobal[inds_use], type = "l")
  }


  # Calculate the rhat values
  # Optional rhat (posterior)
  rhat_sigma2 <- NA_real_
  rhat_lambda2 <- NA_real_
  if (isTRUE(compute_rhat)) {
    if (!requireNamespace("posterior", quietly = TRUE)) {
      message("Skipping R-hat: package 'posterior' is not installed.")
    } else {
      rhat_sigma2  <- posterior::rhat(res_mcmc_last$sigma[inds_use])
      rhat_lambda2 <- posterior::rhat(res_mcmc_last$vglobal[inds_use])
    }
  }


  if (any(is.na(beta_inds))) {
    vbeta_hat = apply(res_mcmc$beta[inds_use,],2,mean)
    beta_inds = order(abs(vbeta_hat),decreasing = TRUE)[1:10]
  }

  ldens_beta <- list()
  for (J in 1:length(beta_inds)) {
    j = beta_inds[J]
    dens  = density(res_mcmc$beta[inds_use,j])
    ldens_beta[[J]] = dens
  }

  dens_sigma2  = density(res_mcmc$sigma[inds_use])
  dens_lambda2 = density(res_mcmc$vglobal[inds_use])


  return(list(
    ldens_beta = ldens_beta,
    dens_sigma2 = dens_sigma2,
    dens_lambda2 = dens_lambda2,
    mStat = mStat,
    beta_inds = beta_inds,
    rhat_sigma2 = rhat_sigma2,
    rhat_lambda2 = rhat_lambda2
  ))
}
