## ----setup, echo = FALSE, eval=TRUE, message=FALSE----------------------------
library(BiocStyle)
knitr::opts_knit$set(
  upload.fun = NULL,
  base.url = NULL) # use local files for images
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#"
)
# Override BiocStyle plot hook to avoid css_align issues
knitr::knit_hooks$set(plot = function(x, options) {
  paste0('![', basename(x), '](', x, ')')
})
runchunks = if (Sys.getenv("FORCE_VIGNETTE_REBUILD", "FALSE") == "TRUE") TRUE else FALSE

cache_file <- '../vignette_cache/pulsar-parallel.RData'
if (!runchunks && file.exists(cache_file)) {
  load(cache_file)
  # If we loaded trimmed objects, reassign them to original names
  if (exists("se1.trimmed")) {
    se1 <- se1.trimmed
    rm(se1.trimmed)
  }
  if (exists("se2.trimmed")) {
    se2 <- se2.trimmed
    rm(se2.trimmed)
  }
  if (exists("se3.trimmed")) {
    se3 <- se3.trimmed
    rm(se3.trimmed)
  }
  if (exists("se4.trimmed")) {
    se4 <- se4.trimmed
    rm(se4.trimmed)
  }
  if (exists("se.serial.trimmed")) {
    se.serial <- se.serial.trimmed
    rm(se.serial.trimmed)
  }

} else {
  if (!runchunks) {
    message("Cache file pulsar-parallel.RData not found - building from scratch")
  }
  runchunks <- TRUE
}
saveout   = runchunks


## ----eval=TRUE----------------------------------------------------------------
library(SpiecEasi)
data(amgut1.filt)

## ----eval=runchunks, message=FALSE, warning=FALSE-----------------------------
# 
# ## Default settings ##
# pargs1 <- list(rep.num=50, seed=10010)
# t1 <- system.time(
# se1 <- spiec.easi(amgut1.filt, method='mb', lambda.min.ratio=1e-3, nlambda=30,
#               sel.criterion='stars', pulsar.select=TRUE, pulsar.params=pargs1)
# )
# 
# ## Platform-specific or envrionment-specific parallel processing ##
# if (.Platform$OS.type == "windows" || Sys.getenv("CI") == "true" || nzchar(Sys.getenv("GITHUB_ACTIONS"))) {
#   # On Windows, use snow cluster or run serial
#   pargs2 <- list(rep.num=50, seed=10010, ncores=1)  # Serial for Windows
#   cat("Running on Windows - using serial processing\n")
# } else {
#   # On Unix-like systems, use multicore
#   n_cores <- min(2, parallel::detectCores())
#   pargs2 <- list(rep.num=50, seed=10010, ncores=n_cores)
#   cat("Running on Unix-like system - using parallel processing\n")
# }
# 
# t2 <- system.time(
# se2 <- spiec.easi(amgut1.filt, method='mb', lambda.min.ratio=1e-3, nlambda=30,
#               sel.criterion='stars', pulsar.select=TRUE, pulsar.params=pargs2)
# )

## ----message=FALSE, warning=FALSE, eval=runchunks-----------------------------
# t3 <- system.time(
# se3 <- spiec.easi(amgut1.filt, method='mb', lambda.min.ratio=1e-3, nlambda=30,
#                sel.criterion='bstars', pulsar.select=TRUE, pulsar.params=pargs1)
# )
# t4 <- system.time(
# se4 <- spiec.easi(amgut1.filt, method='mb', lambda.min.ratio=1e-3, nlambda=30,
#                sel.criterion='bstars', pulsar.select=TRUE, pulsar.params=pargs2)
# )

## ----eval=TRUE----------------------------------------------------------------
## serial vs parallel
identical(getRefit(se1), getRefit(se2))
t1[3] > t2[3]
## stars vs bstars
identical(getRefit(se1), getRefit(se3))
t1[3] > t3[3]
identical(getRefit(se2), getRefit(se4))
t2[3] > t4[3]

## ----eval=runchunks-----------------------------------------------------------
# # Simple serial processing for Windows
# pargs.serial <- list(rep.num=50, seed=10010, ncores=1)
# se.serial <- spiec.easi(amgut1.filt, method='mb', lambda.min.ratio=1e-3, nlambda=30,
#                         sel.criterion='stars', pulsar.select=TRUE, pulsar.params=pargs.serial)

## ----eval=FALSE---------------------------------------------------------------
# # Batch mode works on all platforms
# bargs <- list(rep.num=50, seed=10010, conffile="parallel")
# se.batch <- spiec.easi(amgut1.filt, method='mb', lambda.min.ratio=1e-3, nlambda=30,
#                        sel.criterion='stars', pulsar.select='batch', pulsar.params=bargs)

## ----eval=FALSE---------------------------------------------------------------
# ## bargs <- list(rep.num=50, seed=10010, conffile="path/to/conf.txt")
# bargs <- list(rep.num=50, seed=10010, conffile="parallel")
# ## See the config file stores:
# pulsar::findConfFile('parallel')
# 
# ## uncomment line below to turn off batchtools reporting
# # options(batchtools.verbose=FALSE)
# se5 <- spiec.easi(amgut1.filt, method='mb', lambda.min.ratio=1e-3, nlambda=30,
#             sel.criterion='stars', pulsar.select='batch', pulsar.params=bargs)

## ----eval=TRUE----------------------------------------------------------------
# Compare timing
cat("Serial StARS:", t1[3], "seconds\n")
cat("Platform-specific StARS:", t2[3], "seconds\n")
cat("Serial B-StARS:", t3[3], "seconds\n")
cat("Platform-specific B-StARS:", t4[3], "seconds\n")

# Speedup factors (only meaningful on Unix-like systems)
if (.Platform$OS.type != "windows") {
  cat("Parallel speedup (StARS):", t1[3]/t2[3], "\n")
  cat("B-StARS speedup (serial):", t1[3]/t3[3], "\n")
  cat("B-StARS speedup (parallel):", t2[3]/t4[3], "\n")
}

## -----------------------------------------------------------------------------
sessionInfo()

## ----echo = FALSE, eval=TRUE, message=FALSE-----------------------------------
# Save objects if they exist
if (exists("se1") && exists("se2")) {
  cache_file <- '../vignette_cache/pulsar-parallel.RData'
  tryCatch({
    # Load trimming function and trim objects to reduce size
    source('../vignette_cache/trim_spiec_easi.R')
    se1.trimmed <- trim_spiec_easi(se1)
    se2.trimmed <- trim_spiec_easi(se2)
    se3.trimmed <- trim_spiec_easi(se3)
    se4.trimmed <- trim_spiec_easi(se4)
    se.serial.trimmed <- trim_spiec_easi(se.serial)
    
    # Save trimmed objects and timing data
    save(se1.trimmed, se2.trimmed, se3.trimmed, se4.trimmed, se.serial.trimmed, t1, t2, t3, t4, pargs1, pargs2, file=cache_file)
    message("Saved trimmed cache file: ", cache_file, " in directory: ", getwd())
  }, error = function(e) {
    message("Failed to save cache file: ", e$message)
  })
}