A Simulation Study using Piecemeal

Pavel N. Krivitsky

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.

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:

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:

outdir <- file.path(tempdir(), "piecemeal_demo")

Initialise the Piecemeal object.

sim <- piecemeal::init(outdir) # short for `piecemeal::Piecemeal$new(outdir)`

The following will clear all simulation results.

sim$reset()
#> Finding individual runs

The cluster has 2 nodes. (We could also pass a preexisting cluster object.)

sim$cluster(2)

Set a factorial design with \(x = 1, 2\) and \(y = 1, 3, 9, 27\):

sim$factorial(x = 2^(0:1), y = 3^(0:3))

Use 3 replications per combination of \(x\) and \(y\):

sim$nrep(3)

Set the function to be called for each treatment combination:

sim$worker(f)

Here is our set up so far:

sim
#> A Piecemeal simulation
#> Output directory: /tmp/RtmpGEUHQA/piecemeal_demo 
#> 
#> Design: 8 treatment configurations by 3 seeds = 24 runs
#> 
#> No cluster set up.
#> 
#> Call for each configuration and seed:
#>   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)
#>       dbl(p = p, u = u)
#>   } 
#> 
#> Options:
#>   directory split: 1 level(s) by treatment and 1 level(s) by seed
#>   errored runs: auto 
#> 
#> Ready to execute? Yes.

Note that we could have chained all these together:

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

head(sim$todo(), 2)
#> Finding individual runs
#> [[1]]
#> [[1]]$seed
#> [1] 1
#> 
#> [[1]]$treatment
#> [[1]]$treatment$x
#> [1] 1
#> 
#> [[1]]$treatment$y
#> [1] 1
#> 
#> attr(,"hash")
#> [1] "d5b52728060d2c841d2a84fb1e9aaada"
#> 
#> [[1]]$fn
#> [1] "d5b52728060d2c841d2a84fb1e9aaada.1.rds"
#> 
#> [[1]]$subdirs
#> [1] "d5" "1" 
#> 
#> 
#> [[2]]
#> [[2]]$seed
#> [1] 2
#> 
#> [[2]]$treatment
#> [[2]]$treatment$x
#> [1] 1
#> 
#> [[2]]$treatment$y
#> [1] 1
#> 
#> attr(,"hash")
#> [1] "d5b52728060d2c841d2a84fb1e9aaada"
#> 
#> [[2]]$fn
#> [1] "d5b52728060d2c841d2a84fb1e9aaada.2.rds"
#> 
#> [[2]]$subdirs
#> [1] "d5" "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.

sim$test()
#> Finding individual runs
#> ── Running configuration "cd548e90fc70067d87a8e7b66a920ba1" ────────── seed 1 ──
#> Error in `map()`:
#> ℹ In index: 1.
#> Caused by error in `<worker>`:
#> ! object 'a' not found

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

sim$export_vars("a")
sim$test()
#> Finding individual runs
#> ── Running configuration "214c32b04c7319f2a0a5e53ec022240e" ────────── seed 1 ──
#> Error in `map()`:
#> ℹ In index: 1.
#> Caused by error in `dbl()`:
#> ! could not find function "dbl"

We now have a complaint that we haven’t loaded the rlang package. Let’s fix that:

sim$setup({library(rlang)})
sim$test()
#> Finding individual runs
#> ── Running configuration "d5b52728060d2c841d2a84fb1e9aaada" ────────── seed 1 ──
#> [[1]]
#> [[1]]$seed
#> [1] 1
#> 
#> [[1]]$treatment
#> [[1]]$treatment$x
#> [1] 1
#> 
#> [[1]]$treatment$y
#> [1] 1
#> 
#> attr(,"hash")
#> [1] "d5b52728060d2c841d2a84fb1e9aaada"
#> 
#> [[1]]$fn
#> [1] "d5b52728060d2c841d2a84fb1e9aaada.1.rds"
#> 
#> [[1]]$subdirs
#> [1] "d5" "1" 
#> 
#> [[1]]$output
#>         p         u 
#> 1.0000000 0.2655087

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.

sim$run()
#> Finding individual runs
#> Finding individual runs
#> ℹ Starting 24 runs (0 already done).
#> Run summary:
#>   Error : condition 0: 4
#>   Error : condition 4: 4
#>   OK: 16

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.

sim$status() # or summary(sim)
#> Finding individual runs
#> Finding running workers
#> Scanning individual runs
#> A Piecemeal simulation
#> Output directory: /tmp/RtmpGEUHQA/piecemeal_demo 
#> 
#>                Result Freq
#> 1                Done   16
#> 2 Error : condition 0    4
#> 3 Error : condition 4    4
#> Last successful completion: 2026-10-08 16:27:46 AEDT (0.03 secs ago) 
#> Last consolidation: never 
#> 
#> A Piecemeal simulation ETA calculation
#> Output directory: /tmp/RtmpGEUHQA/piecemeal_demo 
#> Based on 23 completions in 0.006 secs 
#> 
#> Time per completion: 3e-04 secs 
#> Completion rate: 3563 per sec 
#> Estimated time left: 0 secs 
#> Estimated completion time: 2026-10-08 16:27:46 
#> 
#> In progress (approximate): 0

Some runs have succeeded. Here’s what the individual run output files look like:

list.files(outdir, recursive = TRUE)
#>  [1] "1d/1/1d4603447b8408049a1013a0acfcdfaf.1.rds"
#>  [2] "1d/2/1d4603447b8408049a1013a0acfcdfaf.2.rds"
#>  [3] "1d/3/1d4603447b8408049a1013a0acfcdfaf.3.rds"
#>  [4] "21/1/214c32b04c7319f2a0a5e53ec022240e.1.rds"
#>  [5] "21/2/214c32b04c7319f2a0a5e53ec022240e.2.rds"
#>  [6] "21/3/214c32b04c7319f2a0a5e53ec022240e.3.rds"
#>  [7] "34/1/3416c4f9b19a5d4617d9e48c67b63f58.1.rds"
#>  [8] "34/2/3416c4f9b19a5d4617d9e48c67b63f58.2.rds"
#>  [9] "34/3/3416c4f9b19a5d4617d9e48c67b63f58.3.rds"
#> [10] "81/1/815a82aa86afb73a77ac055a7f7bed69.1.rds"
#> [11] "81/2/815a82aa86afb73a77ac055a7f7bed69.2.rds"
#> [12] "81/3/815a82aa86afb73a77ac055a7f7bed69.3.rds"
#> [13] "cc/1/ccca271c0f476c532ef2e95afca7af78.1.rds"
#> [14] "cc/2/ccca271c0f476c532ef2e95afca7af78.2.rds"
#> [15] "cc/3/ccca271c0f476c532ef2e95afca7af78.3.rds"
#> [16] "cd/1/cd548e90fc70067d87a8e7b66a920ba1.1.rds"
#> [17] "cd/2/cd548e90fc70067d87a8e7b66a920ba1.2.rds"
#> [18] "cd/3/cd548e90fc70067d87a8e7b66a920ba1.3.rds"
#> [19] "d5/1/d5b52728060d2c841d2a84fb1e9aaada.1.rds"
#> [20] "d5/2/d5b52728060d2c841d2a84fb1e9aaada.2.rds"
#> [21] "d5/3/d5b52728060d2c841d2a84fb1e9aaada.3.rds"
#> [22] "fc/1/fc8e1832ac3ef120e1eb3126cae3bdb8.1.rds"
#> [23] "fc/2/fc8e1832ac3ef120e1eb3126cae3bdb8.2.rds"
#> [24] "fc/3/fc8e1832ac3ef120e1eb3126cae3bdb8.3.rds"
#> [25] "last_OK"

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:

sim$run()
#> Finding individual runs
#> ℹ Starting 0 runs (24 already done).
#> Run summary:
#>   SKIPPED: 24

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.

sim$result_df()
#> Finding individual runs
#> ! 8/24 runs returned an error.
#>    x  y  p         u .seed
#> 1  1  9  9 0.2655087     1
#> 2  1  9  9 0.1848823     2
#> 3  1  9  9 0.1680415     3
#> 4  1  3  3 0.2655087     1
#> 5  1  3  3 0.1848823     2
#> 6  1  3  3 0.1680415     3
#> 7  2 27 54 0.1848823     2
#> 8  2  9 18 0.1848823     2
#> 9  2  3  6 0.1848823     2
#> 10 2  1  2 0.1848823     2
#> 11 1  1  1 0.2655087     1
#> 12 1  1  1 0.1848823     2
#> 13 1  1  1 0.1680415     3
#> 14 1 27 27 0.2655087     1
#> 15 1 27 27 0.1848823     2
#> 16 1 27 27 0.1680415     3

The remaining errors are stochastic. We can see which treatments and seeds led to unsuccessful runs:

head(sim$erred(), 2)
#> Finding individual runs
#> [[1]]
#> [[1]]$seed
#> [1] 1
#> 
#> [[1]]$treatment
#> [[1]]$treatment$x
#> [1] 2
#> 
#> [[1]]$treatment$y
#> [1] 27
#> 
#> attr(,"hash")
#> [1] "3416c4f9b19a5d4617d9e48c67b63f58"
#> 
#> [[1]]$fn
#> [1] "3416c4f9b19a5d4617d9e48c67b63f58.1.rds"
#> 
#> [[1]]$subdirs
#> [1] "34" "1" 
#> 
#> [[1]]$output
#> [1] "Error : condition 4\n"
#> attr(,"class")
#> [1] "try-error"
#> attr(,"condition")
#> <simpleError: condition 4>
#> 
#> [[1]]$OK
#> [1] FALSE
#> 
#> 
#> [[2]]
#> [[2]]$seed
#> [1] 3
#> 
#> [[2]]$treatment
#> [[2]]$treatment$x
#> [1] 2
#> 
#> [[2]]$treatment$y
#> [1] 27
#> 
#> attr(,"hash")
#> [1] "3416c4f9b19a5d4617d9e48c67b63f58"
#> 
#> [[2]]$fn
#> [1] "3416c4f9b19a5d4617d9e48c67b63f58.3.rds"
#> 
#> [[2]]$subdirs
#> [1] "34" "3" 
#> 
#> [[2]]$output
#> [1] "Error : condition 4\n"
#> attr(,"class")
#> [1] "try-error"
#> attr(,"condition")
#> <simpleError: condition 4>
#> 
#> [[2]]$OK
#> [1] FALSE

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:

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:

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?

sim$worker(f)
sim
#> A Piecemeal simulation
#> Output directory: /tmp/RtmpGEUHQA/piecemeal_demo 
#> 
#> Design: 8 treatment configurations by 3 seeds = 24 runs
#> 
#> No cluster set up.
#> 
#> Setup:
#>   initialise each node with:
#>     {
#>         library(rlang)
#>     } 
#>   variables to copy from each environment:
#>     R_GlobalEnv: 'a'
#> 
#> Call for each configuration and seed:
#>   function (x, y) 
#>   {
#>       p <- x * y
#>       u <- runif(1)
#>       dbl(p = p, u = u)
#>   } 
#> 
#> Options:
#>   directory split: 1 level(s) by treatment and 1 level(s) by seed
#>   errored runs: auto 
#> 
#> Ready to execute? Yes.

Run again:

sim$run()
#> Finding individual runs
#> ℹ Deleting...
#> ✔ 8 failed runs deleted.
#> Finding individual runs
#> ℹ Starting 8 runs (16 already done).
#> Run summary:
#>   OK: 8
#>   SKIPPED: 16

Success!

We have all 24 combinations!

sim$result_df()
#> Finding individual runs
#>    x  y  p         u .seed
#> 1  1  9  9 0.2655087     1
#> 2  1  9  9 0.1848823     2
#> 3  1  9  9 0.1680415     3
#> 4  1  3  3 0.2655087     1
#> 5  1  3  3 0.1848823     2
#> 6  1  3  3 0.1680415     3
#> 7  2 27 54 0.2655087     1
#> 8  2 27 54 0.1848823     2
#> 9  2 27 54 0.1680415     3
#> 10 2  9 18 0.2655087     1
#> 11 2  9 18 0.1848823     2
#> 12 2  9 18 0.1680415     3
#> 13 2  3  6 0.2655087     1
#> 14 2  3  6 0.1848823     2
#> 15 2  3  6 0.1680415     3
#> 16 2  1  2 0.2655087     1
#> 17 2  1  2 0.1848823     2
#> 18 2  1  2 0.1680415     3
#> 19 1  1  1 0.2655087     1
#> 20 1  1  1 0.1848823     2
#> 21 1  1  1 0.1680415     3
#> 22 1 27 27 0.2655087     1
#> 23 1 27 27 0.1848823     2
#> 24 1 27 27 0.1680415     3

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

sim$consolidate()
#> Finding unconsolidated runs
#> ✔ 24 files consolidated.
list.files(outdir, recursive = TRUE)
#> [1] "consolidated.db" "last_OK"

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.

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

sim$eta()
#> Finding running workers
#> Finding individual runs
#> Finding consolidated runs
#> Scanning consolidated runs
#> A Piecemeal simulation ETA calculation
#> Output directory: /tmp/RtmpGEUHQA/piecemeal_demo 
#> Based on 23 completions in 0.2 secs 
#> 
#> Time per completion: 0.009 secs 
#> Completion rate: 111 per sec 
#> Estimated time left: 0.1 secs 
#> Estimated completion time: 2026-10-08 16:27:47 
#> 
#> In progress (approximate): 0

sim$status() will also calculate and print the estimate, but it is typically much slower.

In any case, we resume our run:

sim$run()
#> Finding individual runs
#> Finding consolidated runs
#> ℹ Starting 16 runs (24 already done).
#> Run summary:
#>   OK: 16
#>   SKIPPED: 24
sim$result_df()
#> Finding individual runs
#> Finding consolidated runs
#>    x  y  p         u .seed
#> 1  1  9  9 0.5858003     4
#> 2  1  9  9 0.2002145     5
#> 3  1  3  3 0.5858003     4
#> 4  1  3  3 0.2002145     5
#> 5  2 27 54 0.5858003     4
#> 6  2 27 54 0.2002145     5
#> 7  2  9 18 0.5858003     4
#> 8  2  9 18 0.2002145     5
#> 9  2  3  6 0.5858003     4
#> 10 2  3  6 0.2002145     5
#> 11 2  1  2 0.5858003     4
#> 12 2  1  2 0.2002145     5
#> 13 1  1  1 0.5858003     4
#> 14 1  1  1 0.2002145     5
#> 15 1 27 27 0.5858003     4
#> 16 1 27 27 0.2002145     5
#> 17 1  9  9 0.2655087     1
#> 18 1  9  9 0.1848823     2
#> 19 1  9  9 0.1680415     3
#> 20 1  3  3 0.2655087     1
#> 21 1  3  3 0.1848823     2
#> 22 1  3  3 0.1680415     3
#> 23 2 27 54 0.2655087     1
#> 24 2 27 54 0.1848823     2
#> 25 2 27 54 0.1680415     3
#> 26 2  9 18 0.2655087     1
#> 27 2  9 18 0.1848823     2
#> 28 2  9 18 0.1680415     3
#> 29 2  3  6 0.2655087     1
#> 30 2  3  6 0.1848823     2
#> 31 2  3  6 0.1680415     3
#> 32 2  1  2 0.2655087     1
#> 33 2  1  2 0.1848823     2
#> 34 2  1  2 0.1680415     3
#> 35 1  1  1 0.2655087     1
#> 36 1  1  1 0.1848823     2
#> 37 1  1  1 0.1680415     3
#> 38 1 27 27 0.2655087     1
#> 39 1 27 27 0.1848823     2
#> 40 1 27 27 0.1680415     3