Improved the current SDA workflow to reach the North American runs with 6400 sites.

docker/depends/pecan_package_dependencies.csv

@@ -43,6 +43,7 @@
 "doParallel","*","modules/data.atmosphere","Suggests",FALSE
 "doParallel","*","modules/data.remote","Imports",FALSE
 "doSNOW","*","base/remote","Suggests",FALSE
+"doSNOW","*","modules/assim.sequential","Suggests",FALSE
 "dplR","*","modules/data.land","Imports",FALSE
 "dplyr","*","base/qaqc","Imports",FALSE
 "dplyr","*","base/remote","Imports",FALSE
@@ -62,6 +63,7 @@
 "ellipse","*","modules/assim.batch","Imports",FALSE
 "emdbook","*","modules/assim.sequential","Suggests",FALSE
 "foreach","*","base/remote","Imports",FALSE
+"foreach","*","modules/assim.sequential","Suggests",FALSE
 "foreach","*","modules/data.atmosphere","Suggests",FALSE
 "foreach","*","modules/data.remote","Imports",FALSE
 "fs","*","base/db","Imports",FALSE

models/template/inst/analysis_qsub.job

@@ -0,0 +1,5 @@
+#!/bin/bash -l
+module load R/4.1.2
+echo "require (PEcAnAssimSequential)
+    qsub_analysis('@FOLDER_PATH@')
+    " | R --no-save

modules/assim.sequential/DESCRIPTION

@@ -32,6 +32,8 @@ Imports:
     stringr
 Suggests:
     corrplot,
+    doSNOW,
+    foreach,
     ggrepel,
     emdbook,
     glue,

modules/assim.sequential/NAMESPACE

@@ -44,6 +44,8 @@ export(post.analysis.multisite.ggplot)
 export(postana.bias.plotting.sda)
 export(postana.bias.plotting.sda.corr)
 export(postana.timeser.plotting.sda)
+export(qsub_analysis)
+export(qsub_analysis_submission)
 export(rescaling_stateVars)
 export(rwtmnorm)
 export(sample_met)
@@ -61,6 +63,7 @@ import(furrr)
 import(lubridate)
 import(nimble)
 importFrom(dplyr,"%>%")
+importFrom(foreach,"%dopar%")
 importFrom(lubridate,"%m+%")
 importFrom(magrittr,"%>%")
 importFrom(rlang,.data)

modules/assim.sequential/R/Analysis_sda_block.R

@@ -58,9 +58,18 @@ analysis_sda_block <- function (settings, block.list.all, X, obs.mean, obs.cov,
 
   #parallel for loop over each block.
   PEcAn.logger::logger.info(paste0("Running MCMC ", "for ", length(block.list.all[[t]]), " blocks"))
-  if ("try-error" %in% class(try(block.list.all[[t]] <- furrr::future_map(block.list.all[[t]], MCMC_block_function, .progress = T)))) {
-    PEcAn.logger::logger.severe("Something wrong within the MCMC_block_function function.")
-    return(0)
+  if (!is.null(settings$state.data.assimilation$qsub_analysis)) {# qsub_analysis <- list(cores = 28)
+    if ("try-error" %in% class(try(block.list.all[[t]] <- qsub_analysis_submission(block.list = block.list.all[[t]], 
+                                                                                   outdir = settings$outdir, 
+                                                                                   cores = as.numeric(settings$state.data.assimilation$qsub_analysis$cores))))) {
+      PEcAn.logger::logger.severe("Something wrong within the MCMC_block_function function.")
+      return(0)
+    }
+  } else {
+    if ("try-error" %in% class(try(block.list.all[[t]] <- furrr::future_map(block.list.all[[t]], MCMC_block_function, .progress = T)))) {
+      PEcAn.logger::logger.severe("Something wrong within the MCMC_block_function function.")
+      return(0)
+    }
   }
   PEcAn.logger::logger.info("Completed!")
 
@@ -632,4 +641,117 @@ block.2.vector <- function (block.list, X, H) {
               Pf = Pf,
               mu.a = mu.a,
               Pa = Pa))
+}
+
+##' This function provides means to split large SDA analysis (MCMC) runs into separate `qsub` jobs.
+##' Including job creation, submission, and assemble.
+##' @title qsub_analysis_submission
+##' @param block.list list: MCMC configuration lists for the block SDA analysis.
+##' @param outdir character: SDA output path.
+##' @param job.per.folder numeric: number of jobs per folder.
+##' @param cores numeric: number of cpus used for parallel computaion. Default is NULL.
+##' @export
+##' 
+qsub_analysis_submission <- function(block.list, outdir, job.per.folder = 200, cores = NULL) {
+  L <- length(block.list)
+  # calculate proper folder number based on settings.
+  folder.num <- ceiling(L/job.per.folder)
+  # create folder.
+  # if we have previous outputs, remove them.
+  if (file.exists(file.path(outdir, "qsub_analysis"))) {
+    unlink(file.path(outdir, "qsub_analysis"), recursive = T)
+  }
+  # create new folder.
+  dir.create(file.path(outdir, "qsub_analysis"))
+  # loop over sub-folders.
+  folder.paths <- job.ids <- c()
+  PEcAn.logger::logger.info(paste("Submitting", folder.num, "jobs."))
+  for (i in 1:folder.num) {
+    # create folder for each set of job runs.
+    # calculate start and end index for the current folder.
+    head.num <- (i-1)*job.per.folder + 1
+    if (i*job.per.folder > L) {
+      tail.num <- L
+    } else {
+      tail.num <- i*job.per.folder
+    }
+    # naming and creating folder.
+    folder.name <- paste0("From_", head.num, "_to_", tail.num)
+    folder.path <- file.path(outdir, "qsub_analysis", folder.name)
+    folder.paths <- c(folder.paths, folder.path)
+    dir.create(folder.path)
+    # save corresponding block list to the folder.
+    blocks <- block.list[head.num:tail.num]
+    save(blocks, file = file.path(folder.path, "block.Rdata"))
+    # create job file.
+    jobsh <- readLines(con = system.file("analysis_qsub.job", package = "PEcAn.ModelName"), n = -1, warn=FALSE)
+    jobsh <- gsub("@FOLDER_PATH@", folder.path, jobsh)
+    writeLines(jobsh, con = file.path(folder.path, "job.sh"))
+    # qsub command.
+    qsub <- "qsub -l h_rt=48:00:00 -l buyin -pe omp 28 -V -N @NAME@ -o @STDOUT@ -e @STDERR@ -S /bin/bash"
+    qsub <- gsub("@NAME@", paste0("Job-", i), qsub)
+    qsub <- gsub("@STDOUT@", file.path(folder.path, "stdout.log"), qsub)
+    qsub <- gsub("@STDERR@", file.path(folder.path, "stderr.log"), qsub)
+    qsub <- strsplit(qsub, " (?=([^\"']*\"[^\"']*\")*[^\"']*$)", perl = TRUE)
+    cmd <- qsub[[1]]
+    out <- system2(cmd, file.path(folder.path, "job.sh"), stdout = TRUE, stderr = TRUE)
+    # grab job ids for future job completion detection.
+    job.ids <- c(job.ids, PEcAn.remote::qsub_get_jobid(
+      out = out[length(out)],
+      qsub.jobid = settings$host$qsub.jobid,
+      stop.on.error = TRUE))
+  }
+  # checking results.
+  PEcAn.logger::logger.info("Checking results.")
+  # if remaining number of jobs larger than 0.
+  while (length(job.ids) > 0) {
+    Sys.sleep(60)
+    completed_jobs <- job.ids %>% purrr::map(function(id) {
+      if (PEcAn.remote::qsub_run_finished(
+        run = id,
+        host = host,
+        qstat = qstat)) {
+        return(id)
+      }
+    }) %>% unlist()
+    job.ids <- job.ids[which(!job.ids %in% completed_jobs)]
+  }
+  # assemble results.
+  PEcAn.logger::logger.info("Assembling results.")
+  analysis <- c()
+  for (path in seq_along(folder.paths)) {
+    res_env <- new.env()
+    load(file.path(path, "results.Rdata"), envir = res_env)
+    analysis <- c(analysis, res_env$results)
+  }
+  return(analysis)
+}
+
+##' This function can help to execute `foreach` parallel MCMC sampling given generated MCMC configuration lists.
+##' @title qsub_analysis
+##' @param folder.path character: path where the `block.Rdata` file is stored.
+##' @param cores numeric: number of cpus used for parallel computaion. Default is NULL.
+##' @importFrom foreach %dopar%
+##' @export
+qsub_analysis <- function(folder.path, cores = NULL) {
+  # load file.
+  load(file.path(folder.path, "block.Rdata"))
+  # initialize parallel.
+  if (is.null(cores)) {
+    cores <- parallel::detectCores()
+  }
+  cl <- parallel::makeCluster(cores)
+  doSNOW::registerDoSNOW(cl)
+  # progress bar
+  pb <- utils::txtProgressBar(min=1, max=length(blocks), style=3)
+  progress <- function(n) utils::setTxtProgressBar(pb, n)
+  opts <- list(progress=progress)
+  # parallel computation.
+  l <- NULL # fix GitHub check issue.
+  results <- foreach::foreach(l = blocks, .packages=c("Kendall", "purrr"), .options.snow=opts) %dopar% {
+    MCMC_block_function(l)
+  }
+  # wrap results.
+  parallel::stopCluster(cl)
+  save(results, file = file.path(folder.path, "results.Rdata"))
 }