is_check <- ("CheckExEnv" %in% search()) || any(c("_R_CHECK_TIMINGS_", "_R_CHECK_LICENSE_") %in% names(Sys.getenv())) knitr::opts_chunk$set(eval = !is_check)
Sys.sleep(100)
knitr::opts_chunk$set( comment = NA, quiet = TRUE, progress = FALSE, tidy = FALSE, cache = FALSE, message = FALSE, error = FALSE, # FALSE: always stop execution. warning = TRUE, dpi = 100 )
knitr::opts_knit$set(global.par = TRUE)
par(mar = c(3, 3, 2, 2), mgp = c(1.7, 0.5, 0), las = 1, cex.main = 1, tcl = -0.2, cex.axis = 0.8, cex.lab = 0.8)
Often, environmental models are developed in other languages than R, for example C or FORTRAN. It can significantly speed up processing. In this simple example, it is shown how to perform uncertainty analysis with a model developed in a different language than R. We use an example with a basic model written in C.
The adapted uncertainty propagation analysis approach is based on the Monte Carlo method that computes the output of the model repeatedly, with input values that are randomly sampled from their pdfs. The set of model outputs forms a random sample from the output pdf. The method thus consists of the following steps:
# the 'spup' library contains functions described below # and it loads all the dependencies library(spup) library(dplyr) # a grammar of data manipulation library(readr) # fast I/O library(whisker) # rendering methods library(purrr) # set seed set.seed(12345)
Spatial (or other) inputs to the models are often stored in ASCII files. In that case, when using external models in R we need additional code to:
For rendering ASCII input files, the mustache templating framework is implemented (https://mustache.github.io). In R this is implemented in the package whisker
.
Function template()
allow user to define a 'container' class to store all templates with model inputs. The aim of this class is to organise model input files and perform necessary checks. A print()
method is also provided for the class "template".
A template is simply a model input file with:
.template
. For example, suppose we have a model that needs the input file: input.txt
. This input file contains two parameters "b0" and "b1". The contents of the original file may look like:
read_lines("examples/input.txt")
The corresponding template file should be named as input.txt.template
. It contains:
read_lines("examples/input.txt.template")
We can see that the original numbers are replaced by symbols b0 and b1 placed in moustaches {{
and }}
.
Rendering is the process of replacing the tags in moustaches by text. For this, render-methods for class "character" and "template" are provided. For example:
my_template <- "Hello {{name}}. How are you doing?" my_template %>% render(name = "Winnie the Pooh")
The example above calls method render()
for the class "character". It is also possible to fill an entire table:
my_template <- c( "| x | y |", "|---|---|", "{{#MY_TABLE}}", "| {{X}} | {{Y}} |", "{{/MY_TABLE}}" ) my_table <- data.frame(X = 1:5, Y = letters[1:5]) my_table my_template %>% render(MY_TABLE = unname(rowSplit(my_table))) %>% cat
A template stored as a file will always be rendered on disk. Let's return to our template on disk:
my_template <- template("examples/input.txt.template")
with contents:
my_template %>%
read_lines
Rendering will create a new file, called input.txt
.
my_template %>% render(b0 = 3, b1 = 4)
As can be seen above, the path of this file is also the return value of the render
method. This facilitates further processing by means of the pipe-operator:
my_template %>% render(b0 = 3, b1 = 4) %>% read_lines
An external model can be called from R by means of the system
or system2
function. To facilitate this, spup includes the wrapper function executable()
.
Below is an example of an external model written in the C language:
```{c, eval = FALSE}
int main() {
FILE *fp; double x[9] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0}; double y; double b0; double b1; int i;
fp = fopen("input.txt", "r"); if (fp == NULL) { printf("Can't read input.txt\n"); return 1; } fscanf(fp, "%lf %lf\n", &b0, &b1); fclose(fp);
fp = fopen("output.txt", "w"); if (fp == NULL) { printf("Can't create output.txt\n"); return 1; } else { for (i=0; i<9; i++) { y = b0 + b1 * x[i]; fprintf(fp, "%10.2f\n", y); } fclose(fp); }
return 0; }
You can copy this code into a text file and save with the extension ".C". For example, `dummy_model.C`. This model calculates a value of a simple simple linear regression model. To compile this code to a MS-Windows executable you can use a free C compiler GCC running command `gcc dummy_model.c -o dummy_model`. This will create a file `dummy_model.exe`. We can now wrap this model in R as follows: ```r dummy_model <- executable("examples/dummy_model.exe")
Running the rendering procedure allows to pass any values for b0 ad b1 and the model gives:
# create template my_template <- template("examples/input.txt.template") # render the template render(my_template, b0 = 3.1, b1 = 4.2) # run external model dummy_model() # read output (output file of dummy_model is "output.txt") scan(file = "examples/output.txt", quiet = TRUE)
To perform the uncertainty propagation analysis we need to derive multiple realizations of the model output in steps as follows:
For example:
# number of Monte Carlo runs n_realizations <- 100 n_realizations %>% purrr::rerun({ # render template render(my_template, b0 = rnorm(n = 1), b1 = runif(n = 1)) # run model dummy_model() # read output scan("examples/output.txt", quiet = TRUE) }) %>% set_names(paste0("r", 1:n_realizations)) %>% as_data_frame %>% apply(MARGIN = 1, FUN = quantile)
Thanks to Dennis Walvoort for his valuable contribution to the development of the 'spup' package.
This project has received funding from the European Union’s Seventh Framework Programme for research, technological development and demonstration under grant agreement no 607000.
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.