Newer
Older
# Practical guide for processing large data sets at BSC using startR on HPCs
This guide includes explanations and practical examples for you to learn how to use startR to efficiently process large data sets in parallel on the BSC's HPCs (CTE-Power 9, Marenostrum 4, ...). See the main page of the [**startR**](README.md) project for a general overview of the features of startR, without actual guidance on how to use it.
If you would like to start using startR rightaway on the BSC infrastructure, you can directly go through the "Configuring startR" section, copy/paste the basic startR script example shown at the end of the "Motivation" section onto the text editor of your preference, adjust the paths and user names specified in the `Compute()` function, and run the code in an R session after loading the relevant modules.
What would you do if you had to apply a custom statistical analysis procedure to a 10TB climate data set? Probably, you would need to use a scripting language to write a procedure which is able to retrieve a subset of data from the file system (it would rarely be possible to handle all of it at once on a single node), encode the procedure in that language, and apply it carefully and efficiently to the data. Afterwards, you would need to think of and develop a mechanism to dispatch the job mutiple times in parallel to an HPC of your choice, each of the jobs processing a different subset of the data set. You could do this by hand, but ideally you would use EC-Flow or a similar general purpose workflow manager, which would orchestrate the work for you, and would allow you to monitor and control the progress, as well as keep an easy-to-understand record of what you did, to be reused in the future if needed. The mentioned solution, although it is the recommended way to go, is a demanding one and you could easily spend a few days until you get it running smoothly. Additionally, when developing the job script, you would be exposed to the difficulties of efficiently managing the data and applying the encoded procedure to it.
With the constant increase of resolution (in all possible dimensions) of weather and climate model output, and with the growing need for using computationally demanding analytical methodologies (e.g. bootstraping with thousands of repetitions), this kind of divide-and-conquer approach is becoming indispensable. While tools exist to simplify and automate this complex procedure, they usually require adapting your data to specific formats, migrating to specific database systems, or an advanced knowledge of computer sciences or of specific programming languages or frameworks.
startR is yet another tool which allows the R user to apply user-defined functions or procedures to collections of NetCDF files (see Note 1) as large as desired, transparently using computational resources in HPCs (see Note 2) to minimize the time to solution. Although it has been designed to require as few mandatory technical parameters as possible from the user, an experienced user can configure a number of additional parameters to adjust the execution. startR operates on, and provides as outputs, multidimensional arrays with named dimensions, a basic and widely used data structure in R, and this makes the framework more familiar to the general R user.
Although startR can be difficult to use if learnt from the documentation of its functions, it can also be used with little effort if re-using and adapting already existing startR scripts, such as the ones provided later in this guide. startR scripts are written in R, and are short (usually under 50 lines of code), concise, and easy to read.
Other things you can expect from startR:
- Combining data from multiple model executions or observational data sources.
- Extracting values for a certain time period, geographical location, etc., from a collection of NetCDF files.
- Obtaining data arrays with results of analytical procedures that are to be plotted or stored as RData or NetCDF for later use in the analysis workflow.
Things that are not supposed to be done with startR:
- Curating/homogenizing model output files or generating files to be stored under /esarchive following the department/community conventions. Although metadata is understood and used by startR, its handling is not 100% consistent yet.
If startR is suitable for your use case and decide to use it, you will then need to follow the configuration steps listed in the first section of this guide to make sure startR works on your workstation with the HPC of your choice.
Afterwards, you will need to understand and use five functions, all of them included in the startR package:
- **Start()**, for declaing the data sets to be processed
- **Step()** and **AddStep()**, for specifying the operation to be applied to the data
- **Compute()**, for specifying the HPC to be employed, the number of chunks and cores, and to trigger the computation
- **Collect()** and the **EC-Flow graphical user interface**, for monitoring of the progress and collection of results
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
Next, you can see an example of startR script performing an ensemble mean of a small data set on CTE-Power9, for you to get a broad picture of how the startR functions interact and the information that is represented in a startR script. Note that the `temp_dir` and `ecflow_suite_dir` parameters in the `Compute()` call are user-specific.
```r
library(startR)
repos <- '/esarchive/exp/ecmwf/system5_m1/6hourly/$var$/$var$_$sdate$.nc'
data <- Start(dat = repos,
var = 'tas',
sdate = '20180101',
ensemble = 'all',
time = 'all',
latitude = indices(1:40),
longitude = indices(1:40))
fun <- function(x) {
# Expected inputs:
# x: array with dimensions ('ensemble', 'time')
apply(x + 1, 2, mean)
}
step <- Step(fun,
target_dims = c('ensemble', 'time'),
output_dims = c('time'))
wf <- AddStep(data, step)
res <- Compute(wf,
chunks = list(latitude = 2,
longitude = 2),
threads_load = 2,
threads_compute = 4,
cluster = list(queue_host = 'p9login1.bsc.es',
queue_type = 'slurm',
temp_dir = '/gpfs/scratch/bsc32/bsc32473/startR_hpc/',
lib_dir = '/gpfs/projects/bsc32/share/R_libs/3.5/',
r_module = 'R/3.5.0-foss-2018b',
job_wallclock = '00:10:00',
cores_per_job = 4,
max_jobs = 4,
bidirectional = FALSE,
polling_period = 10
),
ecflow_suite_dir = '/home/Earth/nmanuben/startR_local/')
```
***Note 1****: The data files do not need to be migrated to a database system, nor have to comply any specific convention for their file names, name, number or order of the dimensions and variables, or distribution of files in the file system. Although only files in the NetCDF format are supported for now, plug-in functions (of approximately 200 lines of code) can be programmed to add support for additional file formats.*
***Note 2****: The HPCs startR is designed to run on, are understood as multi-core multi-node clusters. startR relies on a shared file system across all HPC nodes, and does not implement any kind of distributed storage system for now.*
At BSC, the only configuration step you need to follow is to set up passwordless connection with the HPC. You do not need to follow the complete list of deployment steps since all dependencies are already installed for you to use, but you can find the steps listed under the [**Deployment**](inst/doc/deployment.md) section.
Specifically, you need to set up passwordless, userless access from your machine to the HPC login node, and from the HPC login node to your machine if at all possible. In order to establish the connection in one of the directions, you need do the following:
1- generate an ssh pair of keys on the origin host if you do not have one, using `ssh-keygen -t rsa`
2- ssh to the destionation node and create a directory where to store the public key, using `ssh username@hostname_or_ip mkdir -p .ssh`. 'hostname_or_ip' refers to the host name or IP address of the login node of the selected HPC, and 'username' to your account name on the HPC, which may not coincide with the one in your workstation.
3- dump your public key on a new file under that folder, using `cat .ssh/id_rsa.pub | ssh username@hostname_or_ip 'cat >> .ssh/authorized_keys'`
4- adjust the permissions of the repository of keys, using `ssh username@hostname_or_ip "chmod 700 .ssh; chmod 640 .ssh/authorized_keys"`
5- if your username is different on your workstation and on the login node of the HPC, add an entry in the file .ssh/config in your workstation as follows:
```
Host short_name_of_the_host
HostName hostname_or_ip
User username
IdentityFile ~/.ssh/id_rsa
```
After following these steps for the connections in both directions (although from the HPC to the workstation might not be possible), you are good to go.
Do not forget adding the following lines in your .bashrc on CTE-Power if you are planning to run on CTE-Power:
```
if [[ $BSC_MACHINE == "power" ]] ; then
module unuse /apps/modules/modulefiles/applications
module use /gpfs/projects/bsc32/software/rhel/7.4/ppc64le/POWER9/modules/all/
fi
```
You can add the following lines in your .bashrc file on your workstation for convenience:
alias ctp='ssh -XY username@hostname_or_ip'
alias start='module load R CDO ecFlow'
```
## Using startR
If you have successfully gone through the configuration steps, you will just need to run the following commands in a terminal session and a fresh R session will pop up with the startR environment ready to use.
library(startR)
```
### Start()
In order for startR to recognize the data sets you want to process, you first need to declare them. The first step in the declaration of a data set is to build a special path string that encodes where all the involved NetCDF files to be processed are stored. It contains some wildcards in those parts of the path that vary across files. This special path string is also called "path pattern".
Before defining an example path pattern, let's introduce some target NetCDF files. In the esarchive, we can find the following files:
```
/esarchive/exp/ecmwf/system5_m1/6hourly/
|--tas/
| |--tas_19930101.nc
| |--tas_19930201.nc
| | ...
| |--tas_20171201.nc
|--tos/
|--tos_19930101.nc
|--tos_19930201.nc
| ...
|--tos_20171201.nc
```
A path pattern that can be used to encode the location of these files in a compact way is the following:
```r
repos <- '/esarchive/exp/ecmwf/system5_m1/6hourly/$var$/$var$_$sdate$.nc'
```
The wildcards used (the pieces wrapped between '$' symbols) can be given any names you like. They do not necessarily need to be 'var' or 'sdate' or match any specific key word (although in this case, as explained later, the 'var' name will trigger a special feature of `Start()`).
Once the path pattern is specified, a `Start()` call can be built, in which you need to provide, as parameters, the specific values of interest for each of the wildcards (also called outer dimensions, or file dimensions), as well as for each of the dimensions inside the NetCDF files (inner dimensions).
You can check in advance which dimensions are inside the NetCDF files by using e.g. easyNCDF on one of the files:
```r
easyNCDF::NcReadDims('/esarchive/exp/ecmwf/system5_m1/6hourly/tas/tas_19930101.nc',
var_names = 'tas')
This will show the names and lengths of the dimensions of the selected variable:
```r
var longitude latitude ensemble time
1 1296 640 25 860
***Note****: If you check the dimensions of that file with `ncdump -h`, you will realize the 'var' dimension is actually not defined inside. `NcReadDims()` and `Start()`, though, perceive the different variables inside a file as if stored along a virtual dimension called 'var'. You can ignore this for now and assume 'var' is simply a file dimension (since it appears as a wildcard in the path pattern). Read more on this in Note 1 at the end of this section.*
Once we know the dimension names, we have all the information we need to put the `Start()` call together:
```r
data <- Start(dat = repos,
# outer dimensions
var = 'tas',
sdate = '19930101',
# inner dimensions
ensemble = 'all',
time = 'all',
latitude = 'all',
longitude = 'all')
```
For each of the dimensions, the values or indices of interest (a.k.a. selectors) can be specified in three possible ways:
- Using one or more numeric indices, for example `time = indices(c(1, 3, 5))`, or `sdate = indices(3:5)`. In the latter case, the third, fourht and fifth start dates appearing in the file system in alphabetical order would be selected ('19930301', '19930401' and '19930501').
- Using one or more actual values, for example `sdate = values('19930101')`, or `ensemble = values(c('r1i1p1', 'r2i1p1'))`, or `latitude = values(c(10, 10.5, 11))`. The `values()` helper function can be omitted (as shown in the example). See Note 2 for details on how values are handled for inner dimensions.
- Using a list of two numeric values, for example `sdate = indices(list(5, 10))`. This will take all indices from the 5th to the 10th.
- Using a list of two actual values, for example `sdate = values(list('r1i1p1', 'r5i1p1'))` or `latitude = values(list(-45, 75))`. This will take all values, in order, placed between the two values specified (both ends included).
- Using the special keywords 'all', 'first' or 'last'.
The dimensions specified in the `Start()` call do not need to follow any specific order, not even the actual order in the path pattern or inside the file. The order, though, can have an impact on the performance of `Start()` as explained in Note 3.
Running the `Start()` call shown above will display some progress and information messages:
```r
* Exploring files... This will take a variable amount of time depending
* on the issued request and the performance of the file server...
* Detected dimension sizes:
* dat: 1
* var: 1
* sdate: 1
* ensemble: 25
* time: 860
* latitude: 640
* longitude: 1296
* Total size of involved data:
* 1 x 1 x 1 x 25 x 860 x 640 x 1296 x 8 bytes = 132.9 Gb
* Successfully discovered data dimensions.
Warning messages:
1: ! Warning: Parameter 'pattern_dims' not specified. Taking the first dimension,
! 'dat' as 'pattern_dims'.
2: ! Warning: Could not find any pattern dim with explicit data set descriptions (in
! the form of list of lists). Taking the first pattern dim, 'dat', as
! dimension with pattern specifications.
```
The warnings shown are normal, and could be avoided with a more wordy specification of the parameters to the `Start()` function.
The dimensions of the selected data set and the total size are shown. As you have probably noticed, this `Start()` call is very fast, even though several GB of data are involved (see Note 4 on the size of the data in R). This is because `Start()` is simply discovering the location and dimension of the involved data.
You can give a quick look to the collected metadata as follows:
```r
Class 'startR_header' length 9 Start(dat = "/esarchive/exp/ecmwf/system5_m1/6hourly/$var$/$var$_$sdate$.nc", var = "tas", sdate = "19930101", ensemble = "all", time = "all", latitude = "all", ...
..- attr(*, "Dimensions")= Named num [1:7] 1 1 1 25 860 ...
.. ..- attr(*, "names")= chr [1:7] "dat" "var" "sdate" "ensemble" ...
..- attr(*, "Variables")=List of 2
.. ..$ common: NULL
.. ..$ dat1 : NULL
..- attr(*, "ExpectedFiles")= chr [1, 1, 1] "/esarchive/exp/ecmwf/system5_m1/6hourly/tas/tas_19930101.nc"
..- attr(*, "FileSelectors")=List of 1
.. ..$ dat1:List of 3
.. .. ..$ dat :List of 1
.. .. .. ..$ : chr "dat1"
.. .. ..$ var :List of 1
.. .. .. ..$ : chr "tas"
.. .. ..$ sdate:List of 1
.. .. .. ..$ : chr "19930101"
..- attr(*, "PatternDim")= chr "dat"
```
The retrieved information can be accessed with the `attr()` function. For example:
```r
attr(data, 'FileSelectors')$dat1
```
```r
$dat
$dat[[1]]
[1] "dat1"
$var
$var[[1]]
[1] "tas"
$sdate
$sdate[[1]]
[1] "19930101"
```
If you are interested in actually loading the entire data set in your machine you can do so in two ways (***be careful****, loading the data involved in the `Start()` call in this example will most likely stall your machine. Try it with a smaller region or a subset of forecast time steps*):
- adding the parameter `retrieve = TRUE` in your `Start()` call.
- evaluating the object returned by `Start()`: `data_load <- eval(data)`
You may realize that this functionality is similar to the `Load()` function in the s2dverification package. In fact, `Start()` is more advanced and flexible, although `Load()` is more mature and consistent for loading typical seasonal to decadal forecasting data. `Load()` will be adapted in the future to use `Start()` internally.
There are no constrains for the number or names of the outer or inner dimensions used in a `Start()` call. In other words, `Start()` will handle NetCDF files with any number of dimensions with any name, as well as files distributed across folders in complex ways, since you can use customized wildcards in the path pattern.
There are a number of advanced parameters and features in `Start()` to handle heterogeneities across files involved in a `Start()` call, such as the `synonims` parameter, or to handle dimensions extending across multiple NetCDF files, such as the `*_across` parameter. See the documentation on the `Start()` function (https://earth.bsc.es/gitlab/es/startR/blob/master/vignettes/start.md) or in `?Start` for more information.
***Note 1 ****on the 'var' dimension*: as mentioned above in this section, `NcVarReader()` is showing as if a virtual dimension 'var' appears inside the file. The existence of this dimension is justified by the fact that, many times, NetCDF files contain more than one variable. The 'var' dimension should hence be considered a "inner" dimension. But, in our example, the dimension 'var' is also defined as a file dimension in the path pattern. So, following the logic of `Start()`, there would be two 'var' dimensions, one of them outer and the other inner, and we should consequently specify indices for each of them. However, as exception, they are automatically understood to be the same dimension, and the target variable name specified as index for the outer 'var' dimension is also re-used to select the variable inside the file. This is a feature triggered only by the 'var' dimension name and, if other dimension names appeared more than once as inner or outer dimensions, `Start()` would crash throw an error. The feature described here is useful for the very common case where file paths contain the variable name and that variable is the only climate variable inside the file. If this feature was not available, one could still define the data set as shown in the code snippet below, where there would be some redundancy in the `Start()` call and in the dimensions of the resulting array.
***Note 2 ****on providing values as selectors for inner dimensions*: when values are requested for a inner dimension, the corresponding numeric indices are automatically calculated by comparing the provided values with a variable inside the file with the same name as the dimension for which the values have been requested. In the last example where specific values are requested for the latitudes, the variable 'latitude' is automatically retrieved from the files. If the name of the variable does not coincide with the name of the dimension, the parameter `*_var` can be specified in the `Start()` call, as detailed in `?Start`.
***Note 3 ****on the order of dimensions*: neither the file dimensions nor the inner dimensions need to be specified in the same order as they appear in the path pattern or inside the NetCDF files, respectively. The resulting arrays returned by `Start()` will have the dimensions in the same order as requested in `Start()`, so changing the order in the call can potentially trigger automatic reordering of data, which is time consuming. But, depending on the use case, it may be a good idea to ask for a specific dimension order so that the data is properly arranged for making posterior calculations more efficient. Remember that the order of the dimensions in R is "big endian"; the values are consecutive along the first (left-most) dimension. In contrast, the order of dimensions in NetCDF files is "little endian". This means that if you want to respect the order of the data values in memory as stored in the NetCDF files, you should request the dimensions in your `Start()` call in the opposite order.
***Note 4 ****on the size of the data in R*: if you check the size of the involved file in the example `Start()` call used above ('/esarchive/exp/ecmwf/system5_m1/6hourly/tas/tas_19930101.nc'), you will realize it only weights 34GB. Why is the data reported to occupy 134GB then? This is due to two facts: by one side, NetCDF files are usually compressed, and their uncompressed size can be substantially greater. In this case, the uncompressed data would occupy about 72GB. Besides, the variable we are targetting in the example ('tas') is defined as a float variable inside the NetCDF file. This means each value is a 4-byte real number. However, R cannot represent 4-byte real numbers; it always takes 8 bytes to represent a real number. This is why, when float numbers are represented in R, the effective size of the data is doubled. In this case, the 72GB of uncompressed float numbers need to be represented using 132GB of RAM in R.
Once the data sources are declared, you can define the operation to be applied to them. The operation needs to be encapsulated in the form of an R function receiving one or more multidimensional arrays with named dimensions (plus additional helper parameters) and returning one or more multidimensional arrays, which should also have named dimensions. For example:
```r
fun <- function(x) {
r <- sqrt(sum(x ^ 2) / length(x))
for (i in 1:100) {
r <- sqrt(sum(x ^ 2) / length(x))
}
dim(r) <- c(time = 1)
r
}
```
As you see, the function only needs to operate on the essential dimensions, not on the whole set of dimensions of the data set. This example function receives only one multidimensional array (with dimensions 'ensemble' and 'time', although not expressed in the code), and returns one multidimensional array (with a single dimension 'time' of length 1). startR will automatically feed the function with subsets of data with only the essential dinmensions, but first, a startR "step" for the function has to be built with with the `Step()` function.
The `Step()` function requires, as parameters, the names of the dimensions of the input arrays expected by the function (in this example, a single array with the dimensions 'ensemble' and 'time'), as well as the names of the dimensions of the arrays the function returns (in this example, a single array with the dimension 'time'):
```r
step <- Step(fun = fun,
target_dims = c('ensemble', 'time'),
output_dims = c('time'))
```
The step function should ideally expect arrays with the dimensions in the same order as requested in the `Start()` call, and consequently the dimension names specified in the `Step()` function should appear in the same order. If a different order was specified, startR would reorder the subsets for the step function to receive them in the expected dimension order.
Functions that receive or return multiple multidimensional arrays are also supported. In such cases, lists of vectors of dimension names should be provided as `target_dims` or `output_dims`.
Since functions wrapped with the `Step()` function will potentially be called thousands of times, it is recommended to keep them as light as possible by, for example, avoiding calls to the `library()` function to load other packages or interacting with files on disk. See the documentation on the parameter `use_libraries` of the `Step()` function, or consider adding additional parameters to the step function with extra information.
Once the step is built, a workflow of steps can be assembled as follows:
```r
wf <- AddStep(data, step)
```
If the step involved more than one data source, a list of data sources could be provided as first parameter. You can find examples using more than one data source further in this guide.
It is not possible for now to define workflows with more than one step, but this is not a crucial gap since a step function can contain more than one statistical analysis procedure. Furthermore, it is usually enough to perform only the first or two first steps of the analysis workflow on the HPCs because, after these steps, the volume of data involved is reduced substantially and the analysis can go on with conventional methods.
Once the data sources are declared and the workflow is defined, you can proceed to specify the execution parameters (including which platform to run on) and trigger the execution with the `Compute()` function.
Next, a few examples are shown with `Compute()` calls to trigger the processing of a dataset locally (only on the machine where the R session is running) and on two different HPCs (the Earth Sciences fat nodes and CTE-Power9). However, let's first define a `Start()` call that involves a smaller subset of data in order not to make the examples too heavy.
```r
library(startR)
repos <- '/esarchive/exp/ecmwf/system5_m1/6hourly/$var$/$var$_$sdate$.nc'
data <- Start(dat = repos,
var = 'tas',
sdate = '19930101',
ensemble = 'all',
time = indices(1:2),
latitude = 'all',
longitude = 'all')
fun <- function(x) {
# Expected inputs:
# x: array with dimensions ('ensemble', 'time')
apply(x + 1, 2, mean)
}
step <- Step(fun,
target_dims = c('ensemble', 'time'),
output_dims = c('time'))
wf <- AddStep(data, step)
```
The output of the `Start()` call follows, where the size of the selected data is reported.
```r
* Exploring files... This will take a variable amount of time depending
* on the issued request and the performance of the file server...
* Detected dimension sizes:
* dat: 1
* var: 1
* sdate: 1
* ensemble: 25
* time: 2
* latitude: 640
* longitude: 1296
* Total size of involved data:
* 1 x 1 x 1 x 25 x 2 x 640 x 1296 x 8 bytes = 316.4 Mb
* Successfully discovered data dimensions.
Warning messages:
1: ! Warning: Parameter 'pattern_dims' not specified. Taking the first dimension,
! 'dat' as 'pattern_dims'.
2: ! Warning: Could not find any pattern dim with explicit data set descriptions (in
! the form of list of lists). Taking the first pattern dim, 'dat', as
! dimension with pattern specifications.
```
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
When only your own workstation is available, startR can still be useful to process a very large dataset by chunks, thus avoiding a RAM memory overload and consequent crash of the workstation. startR will simply load and process the dataset by chunks, one after the other. The operations defined in the workflow will be applied to each chunk, and the results will be stored on a temporary file. `Compute()` will finally gather and merge the results of each chunk and return a single data object, including one or multiple multidimensional data arrays, and additional metadata.
A list of the dimensions which to split the data along, and the number of slices (or "chunks") to make for each, is the only piece of information required for `Compute()` to run locally. It will only be possible to request chunks for those dimensions not required by any of the steps in the workflow built by `Step()` and `AddStep()`.
Following the worklfow of steps defined in the example, where the step uses 'time' and 'ensemble' as target dimensions, the dimensions remaining for chunking would be 'dat', 'var', 'sdate', 'latitude' and 'longitude'. Note that defining a step which has many target dimensions should be avoided as it will reduce the chunking options.
As an example, we could request for computation performing two chunks along the 'latitude' dimension, and two chunks along the 'longitude' dimension. This would result in the data being processed in 4 chunks of about 80MB (the size of the involved data, 316MB, divided by 4).
Calculate the size of the chunks before executing the computation, and make sure they fit in the RAM of your machine. You should have as much free RAM as 2x or 3x the expected size of one chunk. Read more on adjusting the number of chunks in the section "How to choose the number of chunks, jobs and cores".
```r
res <- Compute(wf,
chunks = list(latitude = 2,
longitude = 2))
```
Once the `Compute()` call is executed, the R session will wait for it to return the result. Progress messages will be shown, with a remaining time estimate after the first chunk has been processed.
```r
* Processing chunks... remaining time estimate soon...
* Loading chunk 1 out of 4 ...
* Processing...
* Remaining time estimate (at 2019-01-27 23:31:07) (neglecting merge
* time): 1.662909 mins
* Loading chunk 2 out of 4 ...
* Processing...
* Loading chunk 3 out of 4 ...
* Processing...
* Loading chunk 4 out of 4 ...
* Processing...
```
After all chunks have been processed, a summary of the execution and its performance is reported.
```r
* Computation ended successfully.
* Number of chunks: 4
* Max. number of concurrent chunks (jobs): 1
* Requested cores per job: NA
* Load threads per chunk: 1
* Compute threads per chunk: 1
* Total time (s): 123.636801242828
* Chunking setup: 0.000516414642333984
* Data upload to cluster: 0
* All chunks: 123.416877985001
* Transfer results from cluster: 0
* Merge: 0.219406843185425
* Each chunk:
* queue:
* mean: 0
* min: 0
* max: 0
* job setup:
* mean: 0
* min: 0
* max: 0
* load:
* mean: 3.07621091604233
* min: 2.05077886581421
* max: 5.05771064758301
* compute:
* mean: 27.7766432762146
* min: 27.4815557003021
* max: 28.199684381485
```
Also, some warning messages will be displayed corresponding to the execution of the `Start()` function to retrieve the data for each chunk.
`Compute()` will return a list of data arrays in your R session, one data array for each output returned by the last step in the specified workflow. In the example here, ony one array is returned.
```r
str(res)
```
```r
List of 1
$ output1: num [1:2, 1, 1, 1, 1:640, 1:1296] 250 250 251 250 251 ...
- attr(*, "startR_compute_profiling")=List of 14
..$ nchunks : num 4
..$ concurrent_chunks: num 1
..$ cores_per_job : logi NA
..$ threads_load : num 1
..$ threads_compute : num 1
..$ bychunks_setup : num 0.000516
..$ transfer : num 0
..$ queue : num 0
..$ job_setup : num 0
..$ load : num [1:4] 5.06 3.14 2.06 2.05
..$ compute : num [1:4] 28.2 27.8 27.5 27.7
..$ transfer_back : num 0
..$ merge : num 0.219
..$ total : num 124
```
The configuration details and profiling information are attached as attributes to the returned list of arrays.
If you check the dimensions of one of the ouput arrays, you will see that it preserves named dimensions, and that the output dimensions of the workflow steps appear on the first prositions (left-most).
```r
dim(res$output1)
```
```r
time dat var sdate latitude longitude
2 1 1 1 640 1296
```
In addition to performing the computation in chunks, you can adjust the number of execution threads to use for the data retrieval stage (with `threads_load`) and for the computation (with `threads_compute`). Using more than 2 threads for the retrieval will usually be perjudicial, since two will already be able to make full use of the bandwidth between the workstation and the data repository. The optimal number of threads for the computation will depend on the number of processors in your machine, the number of cores they have, and the number of threads supported by each of them.
```r
res <- Compute(wf,
chunks = list(latitude = 2,
longitude = 2),
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
threads_load = 2,
threads_compute = 4)
```
```r
* Computation ended successfully.
* Number of chunks: 4
* Max. number of concurrent chunks (jobs): 1
* Requested cores per job: NA
* Load threads per chunk: 2
* Compute threads per chunk: 4
* Total time (s): 44.6467976570129
* Chunking setup: 0.000483036041259766
* Data upload to cluster: 0
* All chunks: 44.4387269020081
* Transfer results from cluster: 0
* Merge: 0.207587718963623
* Each chunk:
* queue:
* mean: 0
* min: 0
* max: 0
* job setup:
* mean: 0
* min: 0
* max: 0
* load:
* mean: 3.08622789382935
* min: 2.77512788772583
* max: 3.93441939353943
* compute:
* mean: 8.02220791578293
* min: 8.01178908348083
* max: 8.03660178184509
In order to run the computation on a HPC, such as the BSC CTE-Power 9, you will need to make sure the passwordless connection with the login node of that HPC is configured, as shown at the beginning of this guide. If possible, in both directions. Also, you will need to know whether there is a shared file system between your workstation and that HPC, and will need information on the number of nodes, cores per node, threads per core, RAM memory per node, and type of workload used by that HPC (Slurm, PBS and LSF supported).
You will need to add two parameters to your `Compute()` call: `cluster` and `ecflow_suite_dir`.
The parameter `ecflow_suite_dir` expects a path to a folder in the workstation where to store temporary files generated for the automatic management of the workflow. As you will see later, the EC-Flow workflow manager is used transparently for this purpose.
The parameter `cluster` expects a list with a number of components that will have to be provided a bit differently depending on the HPC you want to run on. You can see next an example cluster configuration that will execute the previously defined workflow on CTE-Power 9.
```r
res <- Compute(wf,
chunks = list(latitude = 2,
longitude = 2),
threads_load = 2,
threads_compute = 4,
cluster = list(queue_host = 'p9login1.bsc.es',
temp_dir = '/gpfs/scratch/bsc32/bsc32473/startR_hpc/',
lib_dir = '/gpfs/projects/bsc32/share/R_libs/3.5/',
r_module = 'R/3.5.0-foss-2018b',
cores_per_job = 4,
extra_queue_params = list('#SBATCH --mem-per-cpu=3000'),
bidirectional = FALSE,
polling_period = 10
ecflow_suite_dir = '/home/Earth/nmanuben/startR_local/'
)
The cluster components and options are explained next:
- `queue_host`: 'queue_host' must match the 'short_name_of_the_host' you associated to the login node of the selected HPC in your .ssh/config file.
- `queue_type`:
- `temp_dir`:
- `lib_dir`:
- `r_module`:
- `cores_per_job`:
- `job_wallclock`:
- `max_jobs`:
- `extra_queue_params`: better qos, memory, special release: '#SBATCH --reservation=test-rhel-7.5'
- `bidirectional`:
- `polling_period`:
EEEEEEEEEEEEEEEEEEEEEEEEE
After the `Compute()` call is executed, an EC-Flow server is automatically started on your workstation, which will orchestrate the work and dispatch jobs onto the HPC. Thanks to the use of EC-Flow, you will also be able to monitor visually the progress of the execution. See the "Collect and the EC-Flow GUI" section.
The following messages will be displayed upon execution:
* ATTENTION: Dispatching chunks on a remote cluster. Make sure
* passwordless access is properly set in both directions.
ping server(bscearth329:5678) succeeded in 00:00:00.001342 ~1 milliseconds
server is already started
* Processing chunks...
* Remaining time estimate soon...
```
At this point, you may want to check the jobs are being dispatched and executed properly onto the HPC. For that, you can either use the EC-Flow GUI (covered in the next section), or you can `ssh` to the login node of the HPC and check the status of the queue with `squeue` or `qstat`, as shown below.
```
[bsc32473@p9login1 ~]$ squeue
JOBID PARTITION NAME USER ST TIME NODES NODELIST(REASON)
1142418 main /STARTR_ bsc32473 R 0:12 1 p9r3n08
1142419 main /STARTR_ bsc32473 R 0:12 1 p9r3n08
1142420 main /STARTR_ bsc32473 R 0:12 1 p9r3n08
1142421 main /STARTR_ bsc32473 R 0:12 1 p9r3n08
```
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
Here the output of the execution on CTE-Power 9 after waiting for about a minute:
```r
* Remaining time estimate (neglecting queue and merge time) (at
* 2019-01-28 01:16:59): 0 mins (46.22883 secs per chunk)
* Computation ended successfully.
* Number of chunks: 4
* Max. number of concurrent chunks (jobs): 4
* Requested cores per job: 4
* Load threads per chunk: 2
* Compute threads per chunk: 4
* Total time (s): 58.3834805488586
* Chunking setup: 1.13428020477295
* Data upload to cluster: 0
* All chunks: 52.9015641212463
* Transfer results from cluster: 4.05326294898987
* Merge: 0.294373273849487
* Each chunk:
* queue:
* mean: 10.5
* min: 7
* max: 14
* job setup:
* mean: 3.3917236328125
* min: 1.16518902778625
* max: 5.61825823783875
* load:
* mean: 5.77260076999664
* min: 5.27595114707947
* max: 6.26925039291382
* compute:
* mean: 23.7037765979767
* min: 23.4937765598297
* max: 23.9137766361237
* Computation ended successfully.
Warning messages:
1: ! Warning: Parameter 'ecflow_server' has not been specified but execution on
!cluster has been requested. An ecFlow server instance will
!be created on localhost:5678.
2: ! Warning: ATTENTION: The source chunks will be removed from the
!system. Store the result after Collect() ends if needed.
```
As you have probably realized, this execution has been slower than the local execution, even if 4 simultaneous jobs have been executed on CTE-Power. This is due to the small size of the data being processed here. The overhead of queuing and starting jobs at CTE-Power is large compared to the required computation time for this amount of data. The benefit would be obvious in use cases with larger inputs.
Usually, in use cases with larger data inputs, it will be preferrable to add the parameter `wait = FALSE` to your `Compute()` call. With this parameter, `Compute()` will return an object with all the information on your startR execution which you will be able to store in your disk. After doing that, you will be able to close your R session and collect the results later on with the `Collect()` function. This is discussed in the next section.
As mentioned above in the definition of the `cluster` parameters, it is strongly recommended to check the section on "How to choose the number of chunks, jobs and cores".
#### Compute() on the fat nodes and other HPCs
The `Compute()` call with the parameters to run the example in this section on the BSC ES fat nodes is provided below (you will need to adjust some of the parameters before using it). As you can see, the only thing that needs to be changed to execute startR on a different HPC is the definition of the `cluster` parameters.
The `cluster` configuration for the fat nodes, CTE-Power 9, Marenostrum 4, Nord III, Minotauro and ECMWF cca/ccb are all provided at the very end of this guide.
```r
res <- Compute(wf,
chunks = list(latitude = 2,
longitude = 2),
threads_load = 2,
threads_compute = 4,
cluster = list(queue_host = 'bsceslogin01.bsc.es',
temp_dir = '/home/Earth/nmanuben/startR_hpc/',
cores_per_job = 2,
job_wallclock = '00:10:00',
max_jobs = 4,
ecflow_suite_dir = '/home/Earth/nmanuben/startR_local/')
```
### Collect() and the EC-Flow GUI
EEEEEEEEEEEEEE
```r
result <- Collect(res, wait = TRUE)
#### How to choose the number of chunks, jobs and cores
EEEEEEEEEEEEE
#### How to clean a failed execution
EEEEEEEEEEEE
#### Visualizing the profiling of the execution
EEEEEEEEEEEEEE
Discuss PlotProfiling()
#### What to do if your function has too many target dimensions
Computation of weekly means with startR is still pending future work. By now, it is not possible to do that because the metadata associated to each chunk, such as the dates, is not being sent to the `Compute()` function.
### Other examples
#### Using experimental and (date-corresponding) observational data
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
```r
repos <- paste0('/esnas/exp/ecmwf/system4_m1/6hourly/',
'$var$/$var$_$sdate$.nc')
system4 <- Start(dat = repos,
var = 'sfcWind',
#sdate = paste0(1981:2015, '1101'),
sdate = paste0(1981:1984, '1101'),
#time = indices((30*4+1):(120*4)),
time = indices((30*4+1):(30*4+4)),
ensemble = 'all',
#ensemble = indices(1:6),
#latitude = 'all',
latitude = indices(1:10),
#longitude = 'all',
longitude = indices(1:10),
return_vars = list(latitude = NULL,
longitude = NULL,
time = c('sdate')))
repos <- paste0('/esnas/recon/ecmwf/erainterim/6hourly/',
'$var$/$var$_$file_date$.nc')
dates <- attr(system4, 'Variables')$common$time
dates_file <- sort(unique(gsub('-', '', sapply(as.character(dates),
substr, 1, 7))))
erai <- Start(dat = repos,
var = 'sfcWind',
file_date = dates_file,
time = values(dates),
#latitude = 'all',
latitude = indices(1:10),
#longitude = 'all',
longitude = indices(1:10),
time_var = 'time',
time_tolerance = as.difftime(1, units = 'hours'),
time_across = 'file_date',
return_vars = list(latitude = NULL,
longitude = NULL,
time = 'file_date'),
merge_across_dims = TRUE,
split_multiselected_dims = TRUE)
step <- Step(eqmcv_atomic,
list(a = c('ensemble', 'sdate'),
b = c('sdate')),
list(c = c('ensemble', 'sdate')))
res <- Compute(step, list(system4, erai),
chunks = list(latitude = 5,
longitude = 5,
time = 2),
cluster = list(queue_host = 'bsceslogin01.bsc.es',
max_jobs = 4,
cores_per_job = 2),
shared_dir = '/esnas/scratch/nmanuben/test_bychunk',
wait = FALSE)
```
#### Example of computation of weekly means
#### Example with data on an irregular grid with selection of a region
#### Example on MareNostrum 4
```r
library(startR)
repos <- paste0('/esarchive/exp/ecmwf/system5_m1/6hourly/',
'$var$-longitudeS1latitudeS1all/$var$_$sdate$.nc')
# Slower alternative, using files with a less efficient NetCDF
# compression configuration
#repos <- '/esarchive/exp/ecmwf/system5_m1/6hourly/$var$/$var$_$sdate$.nc'
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
data <- Start(dat = repos,
var = 'tas',
sdate = indices(1),
ensemble = 'all',
time = 'all',
latitude = indices(1:40),
longitude = indices(1:40),
retrieve = FALSE)
lons <- attr(data, 'Variables')$common$longitude
lats <- attr(data, 'Variables')$common$latitude
fun <- function(x) apply(x + 1, 2, mean)
step <- Step(fun, c('ensemble', 'time'), c('time'))
wf <- AddStep(data, step)
res <- Compute(wf,
chunks = list(latitude = 2,
longitude = 2),
threads_load = 1,
threads_compute = 2,
cluster = list(queue_host = 'mn2.bsc.es',
queue_type = 'slurm',
data_dir = '/gpfs/projects/bsc32/share/startR_data_repos/',
temp_dir = '/gpfs/scratch/pr1efe00/pr1efe03/startR_tests/',
lib_dir = '/gpfs/projects/bsc32/share/R_libs/3.4/',
r_module = 'R/3.4.0',
cores_per_job = 2,
job_wallclock = '00:10:00',
max_jobs = 4,
extra_queue_params = list('#SBATCH --qos=prace'),
bidirectional = FALSE,
polling_period = 10,
special_setup = 'marenostrum4'
),
ecflow_suite_dir = '/home/Earth/nmanuben/test_remove/',
ecflow_server = NULL,
silent = FALSE,
debug = FALSE,
wait = TRUE)
```
#### Example on CTE-Power using GPUs
#### Seasonal forecast verification example on cca
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
```r
crps <- function(x, y) {
mean(SpecsVerification::EnsCrps(x, y, R.new = Inf))
}
library(startR)
repos <- '/perm/ms/spesiccf/c3ah/qa4seas/data/seasonal/g1x1/ecmf-system4/msmm/atmos/seas/tprate/12/ecmf-system4_msmm_atmos_seas_sfc_$date$_tprate_g1x1_init12.nc'
data <- Start(dat = repos,
var = 'tprate',
date = 'all',
time = 'all',
number = 'all',
latitude = 'all',
longitude = 'all',
return_vars = list(time = 'date'))
dates <- attr(data, 'Variables')$common$time
repos <- '/perm/ms/spesiccf/c3ah/qa4seas/data/ecmf-ei_msmm_atmos_seas_sfc_19910101-20161201_t2m_g1x1_init02.nc'
obs <- Start(dat = repos,
var = 't2m',
time = values(dates),
latitude = 'all',
longitude = 'all',
split_multiselected_dims = TRUE)
s <- Step(crps, target_dims = list(c('date', 'number'), c('date')),
output_dims = NULL)
wf <- AddStep(list(data, obs), s)
r <- Compute(wf,
chunks = list(latitude = 10,
longitude = 3),
cluster = list(queue_host = 'cca',
queue_type = 'pbs',
max_jobs = 10,
init_commands = list('module load ecflow'),
r_module = 'R/3.3.1',
extra_queue_params = list('#PBS -l EC_billing_account=spesiccf')),
ecflow_output_dir = '/perm/ms/spesiccf/c3ah/startR_test/',
is_ecflow_output_dir_shared = FALSE
)
```
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
### Compute() cluster templates
#### CTE-Power9
```r
cluster = list(queue_host = 'p9login1.bsc.es',
queue_type = 'slurm',
temp_dir = '/gpfs/scratch/bsc32/bsc32473/startR_hpc/',
lib_dir = '/gpfs/projects/bsc32/share/R_libs/3.5/',
r_module = 'R/3.5.0-foss-2018b',
cores_per_job = 4,
job_wallclock = '00:10:00',
max_jobs = 4,
bidirectional = FALSE,
polling_period = 10
)
```
#### BSC ES fat nodes
```r
cluster = list(queue_host = 'bsceslogin01.bsc.es',
queue_type = 'slurm',
temp_dir = '/home/Earth/nmanuben/startR_hpc/',
cores_per_job = 2,
job_wallclock = '00:10:00',
max_jobs = 4,
bidirectional = TRUE
)
```
#### Marenostrum 4
```r
cluster = list(queue_host = 'mn2.bsc.es',
queue_type = 'slurm',
data_dir = '/gpfs/projects/bsc32/share/startR_data_repos/',
temp_dir = '/gpfs/scratch/pr1efe00/pr1efe03/startR_hpc/',
lib_dir = '/gpfs/projects/bsc32/share/R_libs/3.4/',
r_module = 'R/3.4.0',
cores_per_job = 2,
job_wallclock = '00:10:00',
max_jobs = 4,
extra_queue_params = list('#SBATCH --qos=prace'),
bidirectional = FALSE,
polling_period = 10,
special_setup = 'marenostrum4'
)
```
#### Nord III
```r
cluster = list(queue_host = 'nord1.bsc.es',
queue_type = 'lsf',
data_dir = '/gpfs/projects/bsc32/share/startR_data_repos/',
temp_dir = '/gpfs/scratch/bsc32/bsc32473/startR_hpc/',
lib_dir = '/gpfs/projects/bsc32/share/R_libs/3.3/',
init_commands = list('module load intel/16.0.1'),
r_module = 'R/3.3.0',
cores_per_job = 2,
job_wallclock = '00:10',
max_jobs = 4,
extra_queue_params = list('#BSUB -q bsc_es'),
bidirectional = FALSE,
polling_period = 10,
special_setup = 'marenostrum4'
)
```
```r
cluster = list(queue_host = 'mt1.bsc.es',
queue_type = 'slurm',
data_dir = '/gpfs/projects/bsc32/share/startR_data_repos/',
temp_dir = '/gpfs/scratch/bsc32/bsc32473/startR_hpc/',
lib_dir = '/gpfs/projects/bsc32/share/R_libs/3.3/',
r_module = 'R/3.3.3',
cores_per_job = 2,
job_wallclock = '00:10:00',
max_jobs = 4,
extra_queue_params = list('#SBATCH --qos=bsc_es'),
bidirectional = FALSE,
polling_period = 10,
special_setup = 'marenostrum4'
)
```
#### ECMWF ecgate (workstation) + cca/ccb (HPC)
```r
cluster = list(queue_host = 'cca',
queue_type = 'pbs',
max_jobs = 10,
init_commands = list('module load ecflow'),
r_module = 'R/3.3.1',
extra_queue_params = list('#PBS -l EC_billing_account=spesiccf')
)
```