---
title: "A Simulation Study using Piecemeal"
author: Pavel N. Krivitsky
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{A Simulation Study using Piecemeal}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
  )
options(cli.ansi = FALSE, cli.dynamic = FALSE)
```

Quick note: This vignette is intended to demonstrate one possible way
debugging of this code might have gone. Because some parts of this
process are stochastic (even outside the worker function) we set a
starting seed that will provide the most informative series of
results.
```{r}
set.seed(1)
```

Suppose that we want to simulate a function of two variables, $x$ and
$y$, for every combination of $x = 1, 2$ and $y = 1, 3, 9, 27$, in a
factorial design.

The function ultimately returns $x\times y$ and a random number
between 0 and 1, but to get there, it relies on an external constant
`a` equalling to 8 that is defined in its environment and on the
package `rlang`.

Unfortunately, our initial implementation also has some bugs that
cause it to crash if the sum of the product and the hundredths digit
of the random number is divisible by 4, with an error message that
depends on the remainder of dividing it by 8.

Here it is:
```{r}
a <- 8
f <- function(x, y) {
  p <- x*y
  u <- runif(1)

  errcond <- p + floor(u * 100) %% 10 + a
  if(errcond %% 4 == 0) stop("condition ", errcond %% 8, call. = FALSE)

  # rlang::dbl
  dbl(p = p, u = u)
}
```

We begin by setting the directory for run results. We will use R's
temporary directory for this demonstration:
```{r}
outdir <- file.path(tempdir(), "piecemeal_demo")
```

Initialise the `Piecemeal` object.
```{r}
sim <- piecemeal::init(outdir) # short for `piecemeal::Piecemeal$new(outdir)`
```

The following will clear all simulation results.
```{r}
sim$reset()
```

The cluster has 2 nodes. (We could also pass a preexisting `cluster` object.)
```{r, eval = FALSE}
sim$cluster(2)
```
Set a factorial design with $x = 1, 2$ and $y = 1, 3, 9, 27$:
```{r}
sim$factorial(x = 2^(0:1), y = 3^(0:3))
```
Use 3 replications per combination of $x$ and $y$:
```{r}
sim$nrep(3)
```
Set the function to be called for each treatment combination:
```{r}
sim$worker(f)
```

Here is our set up so far:
```{r}
sim
```

Note that we could have chained all these together:
```{r, eval = FALSE}
sim$
  cluster(2)$
  factorial(x = 2^(0:1), y = 3^(0:3))$
  nrep(3)$
  worker(f)
```

We can obtain a list of configurations (combinations of treatment and
seed) to be run. (Here, first two are shown.)
```{r}
head(sim$todo(), 2)
```

Let's try a quick test run. The following will run the worker for a
random configuration on the local system (ignoring the cluster
settings). It can be useful to test for basic problems before
submitting the batch job, though it is likely to work best when the
individual runs are relatively short.
```{r, error = TRUE}
sim$test()
```
(Each treatment is identified by a hash which can be used to reference
it for future evaluations of `$test()` if desired.)

Notice that despite running locally, the worker function can't see the
variable `a`. This is because we want the results of `$test()` to be
as representative as possible of what would happen if we were to run
it on a cluster. (Note that certain aspects of system setup, such as
which libraries are loaded, cannot be replicated exactly while
allowing interactive debugging.)

We need to export `a` to the worker nodes and test again:
```{r, error = TRUE}
sim$export_vars("a")
sim$test()
```

We now have a complaint that we haven't loaded the `rlang` package. Let's fix that:
```{r, error = TRUE}
sim$setup({library(rlang)})
sim$test()
```
It works (assuming you are using the same random seed as this vignette)!

We can now execute this setup. It will print a summary of the results of each run.
```{r}
sim$run()
```

We can also view a summary of completed, pending, and ongoing runs any
time. If any runs have been completed, an estimated time to finish the
rest at the current rate will also be estimated.
```{r}
sim$status() # or summary(sim)
```

Some runs have succeeded. Here's what the individual run output files look like:
```{r}
list.files(outdir, recursive = TRUE)
```

Notice that they are split up into subdirectories. This improves
performance on some file systems. It can be controlled by the
`$options(split=)` setting.

If we run again, completed runs (both successful and erred) will be skipped:
```{r}
sim$run()
```

At any time, we can obtain the data frame of successful runs. If your
configuration or output data structure is more complex, you may need
to use custom functions for `$result_df()` or use `$result_list()` instead.
```{r}
sim$result_df()
```

The remaining errors are stochastic. We can see which treatments and
seeds led to unsuccessful runs:
```{r}
head(sim$erred(), 2)
```

We can run an erred configuration (by default the first one) on the
local system (ignoring the cluster settings) with `options(error =
recover)` set. Alternatively, we can debug it from the beginning with
a special argument value:
```{r, eval = FALSE}
sim$debug()
sim$debug(error = "debug")
```

Let's say that based on the above, we were able to fix the bug. Our
new functions is as follows:
```{r}
f <- function(x, y) {
  p <- x*y
  u <- runif(1)
  dbl(p = p, u = u)
}
```

Replace the worker function. What does our simulation look like now?
```{r}
sim$worker(f)
sim
```

Run again:
```{r}
sim$run()
```
Success!

We have all 24 combinations!
```{r}
sim$result_df()
```

In a large simulation, we could end up with tens of thousands of small
files, which can be slow and inefficient on some file systems and
cause [inode
exhaustion](https://en.wikipedia.org/wiki/Inode#Potential_for_inode_exhaustion_and_solutions),
particularly on a shared system where inode counts might be limited by
a quota. To avoid this, a transparent mechanism is provided for
consolidating the result files into an SQLite database. Simply run:
```{r}
sim$consolidate()
list.files(outdir, recursive = TRUE)
```
The individual result files (only the successful ones) have been
replaced by a single database. All of `Piecemeal`'s methods access
these transparently.

Lastly, suppose that we want to run additional 2 replications of each
treatment combination. We can add more replications, and the
simulation will pick up where it left off.
```{r}
sim$nrep(5)
```

Incidentally, we can estimate how much longer the simulation will take
based on the past runs. (This can also be run while the simulation is
running.)
```{r}
sim$eta()
```
`sim$status()` will also calculate and print the estimate, but it is typically much slower.

In any case, we resume our run:
```{r}
sim$run()
sim$result_df()
```