Add the block-based SDA workflow. #3197
Conversation
…tead of the foreach function for the parallelization.
2) Allow pre-specified prior of Q imported from outside. 3) Allow NA values for any site that doesn't have any observation. 4) Solve the one observation or one site problem. 5) Solved the aqq & bqq length issue.
There was a problem hiding this comment.
Shall we add "aqq.Init" and "bqq.Init" under the tag "<state.data.assimilation>"?
There was a problem hiding this comment.
Fixed. You can now add those two arguments inside this script.
| Profiling = FALSE, | ||
| OutlierDetection=FALSE, | ||
| parallel_qsub = TRUE, | ||
| free_run = FALSE, |
There was a problem hiding this comment.
So we probably need to remove "free_run" here? And "SDA_runner.R" should be updated too.
|
|
||
| * pre_enkf_params - (optional) Used for carrying out SDA with pre-existed enkf.params, in which the Pf, aqq, and bqq can be used for the analysis step. Defualt is NULL. | ||
|
|
||
| * ensemble.samples - (optional) Pass ensemble.samples from outside to avoid GitHub check issues. Defualt is NULL. |
There was a problem hiding this comment.
- check spelling on "default" throughout.
- explain what ensemble.samples can be used for productively! This argument doesn't just exist to avoid a github check. Indeed, there's a legitimate question about whether pre_enkf_params and ensemble.samples should just be part of the
restartlist rather than their own arguments (e.g. are there scenarios where you'd use them outside of a restart run?)
| #### **sda.enkf.multisite.R Arguments** | ||
| * settings - (required) [State Data Assimilation Tags Example] settings object | ||
|
|
||
| * obs.mean - (required) List of sites named by site ids, which contains dataframe for observation means, named with observation datetime. |
There was a problem hiding this comment.
I don't know if it's documented somewhere else, but this explanation of how obs.mean and obs.cov have to be structured is clearly insufficient. Also, lets make sure the documentation here matches the ROxygen documentation of the function itself.
There was a problem hiding this comment.
Added documentation about why we need observations, how to generate them, and what are the basic formats for them.
|
|
||
| * control - (optional) List of flags controlling the behaviour of the SDA. Default is as follows: |
There was a problem hiding this comment.
If the list below isn't documented here, point the reader to where they can find an explanation of what all the flags below do (e.g. make sure they're documented in the function's Roxygen
| @@ -324,34 +363,124 @@ $`2010/12/31`$`1000000651` | |||
| [1,] 15.2821691 0.513584319 | |||
| [2,] 0.1213583 0.001162113 | |||
| ``` | |||
| An example of multi-settings pecan xml file also may look like below: | |||
| #### Anlysis SDA workflow | |||
| Before running the SDA analysis functions, the ensemble forecast results have to be generated, and arguments such as H matrix, MCMC arguments, and multi-site Y and R (by `Construct.R` function) have to be generated as well. If you want to proceed the block-based SDA workflow for a free run (estimate the process error to the free run), you need to first specify the `state.data.assimilation` as follows. Beyond that, if any of your sites has zero observation, you should flag the `by.site` as TRUE. | |||
There was a problem hiding this comment.
The sentence "If you want to proceed the block-based SDA workflow for a free run (estimate the process error to the free run)" is confusing since you haven't said what the block-block based SDA is, what a free run is, or what each of the tags below actually does. e.g. is by.site what turns on the block-based version? And if so, why isn't it called by.block and why would you want it turned on if sites have zero observations? What does turning on the free.run tag cause to happen?
| ##' @param pft.path physical path to the pft.csv file. | ||
| ##' @param by criteria, it supports by variable, site, pft, all, and single Q. | ||
| ##' | ||
| ##' @description This function is an upgrade to the Construct.H.multisite function which provides the index by different criteria. |
There was a problem hiding this comment.
If it's an upgrade, why not just fix the previous function? The whole point of a function is that if the inputs and outputs are defined then anything internal should be malleable with affecting the rest of the code.
There was a problem hiding this comment.
I haven't cleaned up this function. Only partial features work now (but this is enough for the current SDA workflow). I will make another PR to build the H matrix under different situations.
| #I think the blocked nimble has to be implemented and used instead of a long vector sampling. | ||
| #1) due to the convergence of X.mod. | ||
| #2) temporal efficiency. MVN sampling over a large cov matrix can be time consuming. | ||
| #4) this data structure design allows us to implement the MCMC sampling parallely. |
There was a problem hiding this comment.
Is this model actually different or is it just called by block instead of overall? I'm not sure I 100% understand why we need a new nimble model.
| @@ -19,6 +19,21 @@ SDA_OBS_Assembler <- function(settings){ | |||
| #extract Obs_Prep object from settings. | |||
| Obs_Prep <- settings$state.data.assimilation$Obs_Prep | |||
|
|
|||
| #check if we want to proceed the free run without any observations. | |||
| if (as.logical(settings$state.data.assimilation$free.run)) { | |||
| #calcualte time points. | |||
| @@ -362,7 +355,7 @@ sda.enkf.multisite <- function(settings, | |||
| ) | |||
| }) | |||
| ###------------------------------------------------------------------------------------------------### | |||
| ### loop over time ### | |||
| ### w over time ### | |||
There was a problem hiding this comment.
not sure I'm following the change in comment
| @@ -20,6 +20,8 @@ pool_ic_netcdf2list <- function(nc.path){ | |||
| for(varname in names(vals)){ | |||
| vals[[varname]] <- ncdf4::ncvar_get(IC.nc,varname) | |||
| } | |||
| ncdf4::nc_close(IC.nc) | |||
| on.exit() | |||
There was a problem hiding this comment.
These two lines should be combined into one line and moved up to right after the nc_open. That ensures that the nc file is closed when the function exits, even if the function exits with an error. See for example
pecan/models/ed/R/model2netcdf.ED2.R
Line 1005 in 91531ae
824bcfa
Description
This PR added the highly modulized block-based SDA workflow, which has complete
try-catchdetections that will facilitate the user-end debugging process. This workflow is also highly parallelized, which will increase the speed a lot.Motivation and Context
The multisite SDA workflow has been problematic due to the high-dimension Bayesian sampling and speed. By using the block-based SDA workflow, people can maintain the multi-site structure and, meanwhile, the covariance structure. For the Wishart case of the MCMC sampling, it's been problematic due to the changes in observations, which is now been solved by employing the global
aqqthat is defined by the number of state variables.Beyond that, this PR has the following updates:
qsub_parallelfunction is updated and more robust.nc_closeandon.exit()paired with eachnc_openfunction for theprepare_poolsfunction.Review Time Estimate
Types of changes
Checklist: