Add the block-based SDA workflow.

CHANGELOG.md

@@ -44,6 +44,7 @@ see if you need to change any of these:
 - Added new features of the SDA function including: 1) allow user-defined free-run mode; 
 2) allow user-defined parallel mode for the qsub submission; 3) allow user-defined email option to report the progress.
 - The analysis function now supports the parallelization of multi-chain MCMC sampling with the fully randomized inits function.
+- Added the new feature of the block-based SDA workflow, which supports the parallel computation.
 
 ### Fixed
 

base/remote/R/qsub_parallel.R

@@ -3,6 +3,7 @@
 #' @param settings pecan settings object
 #' @param files allow submit jobs based on job.sh file paths.
 #' @param prefix used for detecting if jobs are completed or not.
+#' @param sleep time (in second) that we wait each time for the jobs to be completed.
 #' @export
 #' @examples
 #' \dontrun{
@@ -11,7 +12,7 @@
 #' @author Dongchen Zhang
 #' 
 #' @importFrom foreach %dopar%
-qsub_parallel <- function(settings, files = NULL, prefix = "sipnet.out") {
+qsub_parallel <- function(settings, files = NULL, prefix = "sipnet.out", sleep = 10) {
   if("try-error" %in% class(try(find.package("doSNOW"), silent = T))){
     PEcAn.logger::logger.info("Package doSNOW is not installed! Please install it and rerun the function!")
     return(0)
@@ -91,19 +92,17 @@ qsub_parallel <- function(settings, files = NULL, prefix = "sipnet.out") {
   pb <- utils::txtProgressBar(min = 0, max = length(unlist(run_list)), style = 3)
   pbi <- 0
   folders <- file.path(settings$host$outdir, run_list)
-  while (length(folders) > 0) {
-    Sys.sleep(10)
-    completed_folders <- foreach::foreach(folder = folders, settings = rep(settings, length(run_list))) %dopar% {
+  completed_folders <- c()
+  while (length(completed_folders) < length(folders)) {
+    Sys.sleep(sleep)
+    completed_folders <- foreach::foreach(folder = folders) %dopar% {
       if(file.exists(file.path(folder, prefix))){
         return(folder)
       }
     }
-    if(length(unlist(completed_folders)) > 0){
-      Ind <- which(unlist(completed_folders) %in% folders)
-      folders <- folders[-Ind]
-      pbi <- pbi + length(completed_folders)
-      utils::setTxtProgressBar(pb, pbi)
-    }
+    completed_folders <- unlist(completed_folders)
+    pbi <- length(completed_folders)
+    utils::setTxtProgressBar(pb, pbi)
   }
   close(pb)
   parallel::stopCluster(cl)

book_source/03_topical_pages/03_pecan_xml.Rmd

@@ -750,13 +750,16 @@ The following tags can be used for state data assimilation. More detailed inform
 ```xml
 <state.data.assimilation>
 	<process.variance>TRUE</process.variance>
+	<aqq.Init>1</aqq.Init>
+	<bqq.Init>1</bqq.Init>
   <sample.parameters>FALSE</sample.parameters>
   <adjustment>TRUE</adjustment>
   <censored.data>FALSE</censored.data>
   <FullYearNC>TRUE</FullYearNC>
   <NC.Overwrite>FALSE</NC.Overwrite>
   <NC.Prefix>sipnet.out</NC.Prefix>
   <q.type>SINGLE</q.type>
+  <free.run>FALSE</free.run>
   <Localization.FUN>Local.support</Localization.FUN>
   <scalef>1</scalef>
   <chains>5</chains>
@@ -835,13 +838,16 @@ The following tags can be used for state data assimilation. More detailed inform
 ```
 
 * **process.variance** : [optional] TRUE/FLASE flag for if process variance should be estimated (TRUE) or not (FALSE). If TRUE, a generalized ensemble filter will be used. If FALSE, an ensemble Kalman filter will be used. Default is FALSE.
+* **aqq.Init** : [optional] The initial value of `aqq` used for estimate the Q distribution, the default value is 1 (note that, the `aqq.init` and `bqq.init` right now only work on the `VECTOR` q type, and we didn't account for the variabilities of them across sites or variables, meaning we initialize the `aqq` and `bqq` given single value).
+* **bqq.Init** : [optional] The initial value of `bqq` used for estimate the Q distribution, the default value is 1.
 * **sample.parameters** : [optional] TRUE/FLASE flag for if parameters should be sampled for each ensemble member or not. This allows for more spread in the intial conditions of the forecast.
 * **adjustment** : [optional] Bool variable decide if you want to adjust analysis results by the likelihood.
 * **censored.data** : [optional] Bool variable decide if you want to do MCMC sampling for the forecast ensemble space, the default is FALSE.
 * **FullYearNC** : [optional] Bool variable decide if you want to generate the full-year netcdf file when there is a overlap in time, the default is TRUE.
 * **NC.Overwrite** : [optional] Bool variable decide if you want to overwrite the previous netcdf file when there is a overlap in time, the default is FALSE.
 * **NC.Prefix** : [optional] The prefix for the generation of the full-year netcdf file, the default is sipnet.out.
 * **q.type** : [optional] The type of process variance that will be estimated, the default is SINGLE.
+* **free.run** : [optional] If it's a free run without any observations, the default is FALSE.
 * **Localization.FUN** : [optional] The localization function name for the localization operation, the default is Local.support.
 * **scalef** : [optional] The scale parameter used for the localization operation, the smaller the value is, the sites are more isolated.
 * **chains** : [optional] The number of chains needed to be estimated during the MCMC sampling process.

models/sipnet/R/write.configs.SIPNET.R

@@ -441,7 +441,12 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
     plant_wood_vars <- c("AbvGrndWood", "abvGrndWoodFrac", "coarseRootFrac", "fineRootFrac")
     if (all(plant_wood_vars %in% ic.names)) {
       # reconstruct total wood C
-      wood_total_C <- IC$AbvGrndWood / IC$abvGrndWoodFrac
+      if(IC$abvGrndWoodFrac < 0.05){
+        wood_total_C <- IC$AbvGrndWood
+      }else{
+        wood_total_C <- IC$AbvGrndWood / IC$abvGrndWoodFrac
+      }
+
       #Sanity check
       if (is.infinite(wood_total_C) | is.nan(wood_total_C) | wood_total_C < 0) {
         wood_total_C <- 0

modules/assim.sequential/NAMESPACE

@@ -11,6 +11,7 @@ export(EnKF.MultiSite)
 export(GEF)
 export(GEF.MultiSite)
 export(GEF.MultiSite.Nimble)
+export(GrabFillMatrix)
 export(Local.support)
 export(Obs.data.prepare.MultiSite)
 export(Prep_OBS_SDA)
@@ -25,13 +26,15 @@ export(alr)
 export(assessParams)
 export(block_matrix)
 export(conj_wt_wishart_sampler)
+export(construct_nimble_H)
 export(dwtmnorm)
 export(generate_colors_sda)
 export(get_ensemble_weights)
 export(hop_test)
 export(interactive.plotting.sda)
 export(inv.alr)
 export(load_data_paleon_sda)
+export(matrix_network)
 export(metSplit)
 export(obs_timestep2timepoint)
 export(outlier.detector.boxplot)
@@ -57,6 +60,7 @@ export(y_star_create)
 import(furrr)
 import(lubridate)
 import(nimble)
+importFrom(dplyr,"%>%")
 importFrom(lubridate,"%m+%")
 importFrom(magrittr,"%>%")
 importFrom(rlang,.data)