Introduction

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.

Monte Carlo methodology for uncertainty analysis

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:

1. Characterise uncertain model inputs/parameters with pdfs.
2. Repeatedly sample from the pdfs of uncertain inputs/parameters.
3. Run model with sampled inputs and store model outputs.
4. Compute summary statistics of model outputs.

Uncertainty propagation analysis with ‘spup’

# 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:

1. Modify ASCII input files.
2. Run the external model.

Modify ASCII input files - rendering

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:

1. The additional extension .template.
2. Input that needs to be modified is replaced by mustache-style tags.

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")

[1] "3 4"

The corresponding template file should be named as input.txt.template. It contains:

read_lines("examples/input.txt.template")

[1] "{{b0}} {{b1}}"

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")

[1] "Hello Winnie the Pooh. How are you doing?"

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

  X Y
1 1 a
2 2 b
3 3 c
4 4 d
5 5 e

my_template %>% 
    render(MY_TABLE = unname(rowSplit(my_table))) %>%
    cat

| x | y |
|---|---|
| 1 | a |
| 2 | b |
| 3 | c |
| 4 | d |
| 5 | e |

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

[1] "{{b0}} {{b1}}"

Rendering will create a new file, called input.txt.

my_template %>% 
    render(b0 = 3, b1 = 4)

[1] "examples/input.txt"

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

[1] "3 4"

Running external models

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:

#include <stdio.h>

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:

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:

1. Render the template.
2. Run the model.
3. Read the results.
4. Process the results.

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)