06. Large-Scale and Parallel RTM Simulation

library(ToolsRTM)
library(doParallel)
library(foreach)

Every simulation so far used a plain sapply()/loop over LUT rows – fine for hundreds of rows, slow for thousands. This page shows the doParallel/foreach pattern this package’s own course scripts use (Scripts/R/ForPROSAIL/1-GetSimulationsLUTs.R) to distribute simulations across CPU cores.

1. Build a LUT

n.samples <- 300
LUT <- as.data.frame(getLUT(inputs = ToolsRTM::inputsPROSAIL, nLUT = n.samples, setseed = 1234))
rsoil <- rep(0.15, 2101)

2. Sequential baseline

t_seq <- system.time({
  sims_seq <- lapply(seq_len(n.samples), function(i) {
    sail_i <- foursail(inputLUT = LUT[i, ], rsoil = rsoil, LeafModel = "PROSPECT-D")
    sail_i$rsot
  })
})
cat("Sequential:", round(t_seq[["elapsed"]], 2), "s for", n.samples, "simulations\n")
#> Sequential: 1.49 s for 300 simulations

3. Parallel with `doParallel`/`foreach`

no_cores <- max(1, parallel::detectCores() - 2)
# CRAN check machines cap parallel workers at 2; respect that limit so this
# vignette rebuilds cleanly under R CMD check --as-cran.
if (nzchar(Sys.getenv("_R_CHECK_LIMIT_CORES_"))) no_cores <- min(no_cores, 2L)
cl <- makeCluster(no_cores)
registerDoParallel(cl)

t_par <- system.time({
  sims_par <- foreach(i = 1:n.samples, .packages = "ToolsRTM") %dopar% {
    sail_i <- foursail(inputLUT = LUT[i, ], rsoil = rsoil, LeafModel = "PROSPECT-D")
    sail_i$rsot
  }
})
stopCluster(cl)

cat("Parallel (", no_cores, "cores):", round(t_par[["elapsed"]], 2), "s for", n.samples, "simulations\n")
#> Parallel ( 10 cores): 3.33 s for 300 simulations
cat("Speedup:", round(t_seq[["elapsed"]] / t_par[["elapsed"]], 2), "x\n")
#> Speedup: 0.45 x

Three things matter for this pattern to actually work, not just look parallel:

makeCluster()/registerDoParallel()/stopCluster() always paired – an unclosed cluster leaks worker processes across the whole R session, compounding if this chunk runs more than once (e.g. re-knitting interactively).
.packages = "ToolsRTM" inside foreach() – each worker process starts fresh, with no packages attached; without this, foursail() isn’t found on the workers and every iteration errors.
no_cores <- detectCores() - 2, not detectCores() itself – leaves headroom for the OS and this R session’s own main process, a convention this package’s own course scripts (ForPROSAIL, Sensibility) use throughout.

identical(sims_seq[[1]], sims_par[[1]])  # same result, sequential or parallel
#> [1] TRUE

4. When parallel is (and isn’t) worth it

cat("Per-simulation cost: sequential", round(1000 * t_seq[["elapsed"]] / n.samples, 2),
    "ms, parallel", round(1000 * t_par[["elapsed"]] / n.samples, 2), "ms\n")
#> Per-simulation cost: sequential 4.97 ms, parallel 11.1 ms

foursail() itself is fast (milliseconds per call) – the parallel speedup above has to overcome the fixed overhead of starting a cluster and shipping ToolsRTM/data to every worker, so it only pays off once n.samples is large enough (thousands, not hundreds) or each simulation is individually expensive (e.g. SPART()’s atmosphere calculations, Tutorial 03; the full SCOPE energy balance, covered in the SCOPEinR package’s own tutorials). For a few hundred rows, the sequential lapply()/sapply() used throughout Tutorials 01-05 is often simpler and not meaningfully slower once cluster startup is accounted for.

5. Filtering non-finite rows

Large parallel runs sometimes include a few unstable simulations (e.g. Liberty’s leaf solver, Tutorial 02, near the edges of its own parameter range) that come back as NaN rather than erroring outright. Always check before using the results downstream:

spectra_mat <- do.call(rbind, sims_par)
ok <- apply(spectra_mat, 1, function(r) all(is.finite(r)))
cat(sum(!ok), "of", n.samples, "simulations were non-finite and would be dropped.\n")
#> 0 of 300 simulations were non-finite and would be dropped.
spectra_mat <- spectra_mat[ok, , drop = FALSE]
LUT_ok <- LUT[ok, ]

What’s next