LANDISDisp: optional OpenMP threads for the spiral scan

NEWS.md

@@ -1,5 +1,12 @@
 # LandR (development version)
 
+* `LANDISDisp()`'s C++ spiral scan can use OpenMP threads: `options(LandR.LANDISDisp.threads = n)`
+  (default `1`). Only the search for receivers with a source at each spiral offset is
+  parallel, with a static schedule and per-thread buffers joined in thread order; the uniform
+  draws stay serial on the main thread, so output is bit-identical for any thread count and to
+  `useCpp = FALSE`. `src/Makevars` and `src/Makevars.win` add `$(SHLIB_OPENMP_CXXFLAGS)`; without
+  OpenMP the package builds and runs serially.
+
 * `prepInputsFireYear()` has a new argument `latestYear` (default `Inf`, no change for existing calls) that drops fires after that year before rasterizing. A pixel that burned twice holds only its most recent fire, so for an earlier data year `replaceAgeInFires()` skipped it (the later fire was after `startTime`) and the pixel kept its satellite age. `prepInputsStandAgeMap()` now passes `dataYear` when it builds the fire raster, so a pixel burned in 1980 and 2015 gets the age since 1980 for `dataYear` 1990. A `firePerimeters` raster supplied by the caller is used as given.
 
 * `updateCohortData()` is faster when it re-groups all pixels (after seed dispersal): it builds the community string of each `pixelGroup` once and gives each pixel its group's string, instead of expanding `cohortData` to one row per pixel and cohort (20M rows at 10M cohorts) and pasting a string for each. The `pixelGroup` numbers, the rows kept and their order, the columns and the RNG state are those of the expanded version (a dispersal event took 46 s instead of 37 s with 10.6M cohorts and 1.3M treed pixels, and 18.5 s instead of 8.7 s with 0.4M cohorts and 1.3M treed pixels, with about 2 GB less peak memory in both). `calcSiteShade()` computes its two sums per `pixelGroup` in one pass and the shade cut once per ecoregion and `bAM` instead of once per `pixelGroup` (the first cohort of a `pixelGroup` is still the one read, as before), and `rmMissingCohorts()` works on the map values as a vector instead of building a table with a row per cell. Results are identical.

R/RcppExports.R

@@ -9,8 +9,12 @@ cellsInPgsCpp <- function(pgv, pgs) {
     .Call(`_LandR_cellsInPgsCpp`, pgv, pgs)
 }
 
-spiralLoopCpp <- function(pixelIndex_in, speciesCode_in, rowOrig_in, colOrig_in, seeddist_max_perRow, spiralRow, spiralCol, spiralCurDist, pgmRows, pgmCols, pgv, srcPg, srcSpeciesCode, numSp, wardProbByDist, activeSpMaxDist, activeSpMax, cellSize, successionTimestep, verbose, wardAlreadyExp, debug = FALSE) {
-    .Call(`_LandR_spiralLoopCpp`, pixelIndex_in, speciesCode_in, rowOrig_in, colOrig_in, seeddist_max_perRow, spiralRow, spiralCol, spiralCurDist, pgmRows, pgmCols, pgv, srcPg, srcSpeciesCode, numSp, wardProbByDist, activeSpMaxDist, activeSpMax, cellSize, successionTimestep, verbose, wardAlreadyExp, debug)
+landisDispHasOpenMP <- function() {
+    .Call(`_LandR_landisDispHasOpenMP`)
+}
+
+spiralLoopCpp <- function(pixelIndex_in, speciesCode_in, rowOrig_in, colOrig_in, seeddist_max_perRow, spiralRow, spiralCol, spiralCurDist, pgmRows, pgmCols, pgv, srcPg, srcSpeciesCode, numSp, wardProbByDist, activeSpMaxDist, activeSpMax, cellSize, successionTimestep, verbose, wardAlreadyExp, debug = FALSE, nThreads = 1L) {
+    .Call(`_LandR_spiralLoopCpp`, pixelIndex_in, speciesCode_in, rowOrig_in, colOrig_in, seeddist_max_perRow, spiralRow, spiralCol, spiralCurDist, pgmRows, pgmCols, pgv, srcPg, srcSpeciesCode, numSp, wardProbByDist, activeSpMaxDist, activeSpMax, cellSize, successionTimestep, verbose, wardAlreadyExp, debug, nThreads)
 }
 
 windowCountsByClassCpp <- function(lccVals, candClasses, nrow, ncol, cells0, kx, ky) {

R/seedDispersalLANDIS.R

@@ -34,6 +34,19 @@ utils::globalVariables(c(
 #' To fall back to the R reference (e.g., for debugging), pass
 #' `useCpp = FALSE` or set `options(LandR.LANDISDisp.useCpp = FALSE)`.
 #'
+#' The scan for receivers that have a source at each spiral offset can use several
+#' OpenMP threads: set `options(LandR.LANDISDisp.threads = n)` (default `1`).
+#' Random draws stay on the main thread in the same order, so the output is identical
+#' for any number of threads. Without OpenMP support the option has no effect.
+#' Two cautions when setting it above `1`. First, every R process that calls
+#' `LANDISDisp()` starts its own threads, so with several parallel workers keep
+#' workers x threads at or below the free cores, or leave the option at `1`.
+#' Second, GCC's OpenMP runtime is not fork-safe: a process that has already run
+#' OpenMP threads and then forks (e.g., `parallel::mclapply()` or a fork cluster)
+#' can hang when the child uses OpenMP. Set the option inside each worker after it
+#' starts, or only in single-process runs. With the default `1` no OpenMP threads
+#' are started, so neither applies.
+#'
 #' @param dtSrc data.table
 #'
 #' @param dtRcv data.table
@@ -942,7 +955,8 @@ spiralSeedDispersalCpp <- function(speciesTable, pixelGroupMap, dtRcvLong,
     successionTimestep  = as.integer(successionTimestep),
     verbose             = as.integer(verbose),
     wardAlreadyExp      = TRUE,
-    debug               = isTRUE(getOption("LandR.LANDISDisp.debug", FALSE))
+    debug               = isTRUE(getOption("LandR.LANDISDisp.debug", FALSE)),
+    nThreads            = as.integer(getOption("LandR.LANDISDisp.threads", 1L))
   )
 
   whSuccess <- which(res$Success)

src/Makevars

@@ -0,0 +1,2 @@
+PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS)
+PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS)

src/Makevars.win

@@ -0,0 +1,2 @@
+PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS)
+PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS)

src/RcppExports.cpp

@@ -36,9 +36,19 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// landisDispHasOpenMP
+bool landisDispHasOpenMP();
+RcppExport SEXP _LandR_landisDispHasOpenMP() {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    rcpp_result_gen = Rcpp::wrap(landisDispHasOpenMP());
+    return rcpp_result_gen;
+END_RCPP
+}
 // spiralLoopCpp
-List spiralLoopCpp(IntegerVector pixelIndex_in, IntegerVector speciesCode_in, IntegerVector rowOrig_in, IntegerVector colOrig_in, IntegerVector seeddist_max_perRow, IntegerVector spiralRow, IntegerVector spiralCol, NumericVector spiralCurDist, int pgmRows, int pgmCols, IntegerVector pgv, IntegerVector srcPg, IntegerVector srcSpeciesCode, int numSp, NumericMatrix wardProbByDist, NumericVector activeSpMaxDist, NumericVector activeSpMax, double cellSize, int successionTimestep, int verbose, bool wardAlreadyExp, bool debug);
-RcppExport SEXP _LandR_spiralLoopCpp(SEXP pixelIndex_inSEXP, SEXP speciesCode_inSEXP, SEXP rowOrig_inSEXP, SEXP colOrig_inSEXP, SEXP seeddist_max_perRowSEXP, SEXP spiralRowSEXP, SEXP spiralColSEXP, SEXP spiralCurDistSEXP, SEXP pgmRowsSEXP, SEXP pgmColsSEXP, SEXP pgvSEXP, SEXP srcPgSEXP, SEXP srcSpeciesCodeSEXP, SEXP numSpSEXP, SEXP wardProbByDistSEXP, SEXP activeSpMaxDistSEXP, SEXP activeSpMaxSEXP, SEXP cellSizeSEXP, SEXP successionTimestepSEXP, SEXP verboseSEXP, SEXP wardAlreadyExpSEXP, SEXP debugSEXP) {
+List spiralLoopCpp(IntegerVector pixelIndex_in, IntegerVector speciesCode_in, IntegerVector rowOrig_in, IntegerVector colOrig_in, IntegerVector seeddist_max_perRow, IntegerVector spiralRow, IntegerVector spiralCol, NumericVector spiralCurDist, int pgmRows, int pgmCols, IntegerVector pgv, IntegerVector srcPg, IntegerVector srcSpeciesCode, int numSp, NumericMatrix wardProbByDist, NumericVector activeSpMaxDist, NumericVector activeSpMax, double cellSize, int successionTimestep, int verbose, bool wardAlreadyExp, bool debug, int nThreads);
+RcppExport SEXP _LandR_spiralLoopCpp(SEXP pixelIndex_inSEXP, SEXP speciesCode_inSEXP, SEXP rowOrig_inSEXP, SEXP colOrig_inSEXP, SEXP seeddist_max_perRowSEXP, SEXP spiralRowSEXP, SEXP spiralColSEXP, SEXP spiralCurDistSEXP, SEXP pgmRowsSEXP, SEXP pgmColsSEXP, SEXP pgvSEXP, SEXP srcPgSEXP, SEXP srcSpeciesCodeSEXP, SEXP numSpSEXP, SEXP wardProbByDistSEXP, SEXP activeSpMaxDistSEXP, SEXP activeSpMaxSEXP, SEXP cellSizeSEXP, SEXP successionTimestepSEXP, SEXP verboseSEXP, SEXP wardAlreadyExpSEXP, SEXP debugSEXP, SEXP nThreadsSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
@@ -64,7 +74,8 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< int >::type verbose(verboseSEXP);
     Rcpp::traits::input_parameter< bool >::type wardAlreadyExp(wardAlreadyExpSEXP);
     Rcpp::traits::input_parameter< bool >::type debug(debugSEXP);
-    rcpp_result_gen = Rcpp::wrap(spiralLoopCpp(pixelIndex_in, speciesCode_in, rowOrig_in, colOrig_in, seeddist_max_perRow, spiralRow, spiralCol, spiralCurDist, pgmRows, pgmCols, pgv, srcPg, srcSpeciesCode, numSp, wardProbByDist, activeSpMaxDist, activeSpMax, cellSize, successionTimestep, verbose, wardAlreadyExp, debug));
+    Rcpp::traits::input_parameter< int >::type nThreads(nThreadsSEXP);
+    rcpp_result_gen = Rcpp::wrap(spiralLoopCpp(pixelIndex_in, speciesCode_in, rowOrig_in, colOrig_in, seeddist_max_perRow, spiralRow, spiralCol, spiralCurDist, pgmRows, pgmCols, pgv, srcPg, srcSpeciesCode, numSp, wardProbByDist, activeSpMaxDist, activeSpMax, cellSize, successionTimestep, verbose, wardAlreadyExp, debug, nThreads));
     return rcpp_result_gen;
 END_RCPP
 }
@@ -89,7 +100,8 @@ END_RCPP
 static const R_CallMethodDef CallEntries[] = {
     {"_LandR_pixelUnifCpp", (DL_FUNC) &_LandR_pixelUnifCpp, 4},
     {"_LandR_cellsInPgsCpp", (DL_FUNC) &_LandR_cellsInPgsCpp, 2},
-    {"_LandR_spiralLoopCpp", (DL_FUNC) &_LandR_spiralLoopCpp, 22},
+    {"_LandR_landisDispHasOpenMP", (DL_FUNC) &_LandR_landisDispHasOpenMP, 0},
+    {"_LandR_spiralLoopCpp", (DL_FUNC) &_LandR_spiralLoopCpp, 23},
     {"_LandR_windowCountsByClassCpp", (DL_FUNC) &_LandR_windowCountsByClassCpp, 7},
     {NULL, NULL, 0}
 };

src/spiralSeedDispersal.cpp

@@ -27,6 +27,9 @@
 #include <Rcpp.h>
 #include <cstdint>
 #include <vector>
+#ifdef _OPENMP
+#include <omp.h>
+#endif
 
 using namespace Rcpp;
 
@@ -92,6 +95,20 @@ IntegerVector cellsInPgsCpp(IntegerVector pgv, IntegerVector pgs) {
 //   Success        : LogicalVector length numRcv
 //   DistOfSuccess  : NumericVector length numRcv (NA_real_ if no success)
 
+// TRUE if the package was compiled with OpenMP.
+// [[Rcpp::export]]
+bool landisDispHasOpenMP() {
+#ifdef _OPENMP
+  return true;
+#else
+  return false;
+#endif
+}
+
+// Steps with fewer active receivers than this are scanned serially: opening a
+// parallel region costs more than it saves.
+static const size_t kOmpMinActive = 20000;
+
 // [[Rcpp::export]]
 List spiralLoopCpp(IntegerVector pixelIndex_in,
                    IntegerVector speciesCode_in,
@@ -114,7 +131,8 @@ List spiralLoopCpp(IntegerVector pixelIndex_in,
                    int successionTimestep,
                    int verbose,
                    bool wardAlreadyExp,
-                   bool debug = false) {
+                   bool debug = false,
+                   int nThreads = 1) {
 
   const int numRcv = pixelIndex_in.size();
   const int numSpiral = spiralRow.size();
@@ -189,6 +207,24 @@ List spiralLoopCpp(IntegerVector pixelIndex_in,
   // Number of currently-live receivers (for early termination)
   int numLive = numRcv;
 
+  if (nThreads < 1) nThreads = 1;
+#ifdef _OPENMP
+  const int maxThreads = omp_get_max_threads();
+  if (nThreads > maxThreads) nThreads = maxThreads;
+#else
+  nThreads = 1;
+#endif
+  // One hit buffer per thread; concatenated in thread order so hasSpRows has
+  // the same order as the serial scan.
+  std::vector< std::vector<int> > threadHits(nThreads);
+
+  // Raw pointers: no R API is touched inside the parallel scan.
+  const int* pgvP = pgv.begin();
+  const std::uint64_t* maskP = srcPgBitmask.data();
+  const int* rowActP = rowAct.data();
+  const int* colActP = colAct.data();
+  const int* spActP  = spAct.data();
+
   GetRNGstate();
 
   // Workspace arrays reused across iterations to avoid reallocation
@@ -271,26 +307,44 @@ List spiralLoopCpp(IntegerVector pixelIndex_in,
     // — important for RNG-stream parity with the R reference.
     hasSpRows.clear();
 
-    for (size_t idx = 0; idx < active.size(); ++idx) {
-      const int j = active[idx];
-      const int sp = spAct[j];
+    // The scan uses no RNG. Above the size threshold it runs in a parallel
+    // region with a static schedule (contiguous chunks in thread order).
+    const int nAct = (int) active.size();
+    const int* activeP = active.data();
+    // Scan one receiver; append to `out` if it has a source at this offset.
+#define LANDIS_SCAN(IDX, OUT)                                                 \
+    {                                                                         \
+      const int j = activeP[IDX];                                             \
+      const int newRow = rowActP[j] + sRow;                                   \
+      const int newCol = colActP[j] + sCol;                                   \
+      if (newRow >= 1 && newRow <= pgmRows && newCol >= 1 && newCol <= pgmCols) { \
+        const int pg = pgvP[(newRow - 1) * pgmCols + (newCol - 1)];           \
+        if (pg != NA_INTEGER && pg >= 0 && pg <= maxPg &&                     \
+            ((maskP[(size_t) pg] >> (spActP[j] - 1)) & 1ULL) != 0ULL) {       \
+          (OUT).push_back(j);                                                 \
+        }                                                                     \
+      }                                                                       \
+    }
 
-      const int newRow = rowAct[j] + sRow;
-      const int newCol = colAct[j] + sCol;
-      if (newRow < 1 || newRow > pgmRows ||
-          newCol < 1 || newCol > pgmCols) {
-        continue; // out of bounds → no source
+    hasSpRows.clear();
+#ifdef _OPENMP
+    if (nThreads > 1 && active.size() >= kOmpMinActive) {
+      for (int t = 0; t < nThreads; ++t) threadHits[t].clear();
+#pragma omp parallel num_threads(nThreads)
+      {
+        std::vector<int>& mine = threadHits[omp_get_thread_num()];
+#pragma omp for schedule(static)
+        for (int idx = 0; idx < nAct; ++idx) LANDIS_SCAN(idx, mine)
+      }
+      for (int t = 0; t < nThreads; ++t) {
+        hasSpRows.insert(hasSpRows.end(), threadHits[t].begin(), threadHits[t].end());
       }
-      // Cell index using terra's row-major numbering, 0-based
-      const int cell0 = (newRow - 1) * pgmCols + (newCol - 1);
-      const int pg = pgv[cell0];
-      if (pg == NA_INTEGER) continue;          // masked cell
-      if (pg < 0 || pg > maxPg) continue;      // pg outside dtSrc range
-      const std::uint64_t mask = srcPgBitmask[(size_t) pg];
-      if (((mask >> (sp - 1)) & 1ULL) == 0ULL) continue;
-
-      hasSpRows.push_back(j);
+    } else
+#endif
+    {
+      for (int idx = 0; idx < nAct; ++idx) LANDIS_SCAN(idx, hasSpRows)
     }
+#undef LANDIS_SCAN
 
     const int sumHasSp = (int) hasSpRows.size();
     if (sumHasSp == 0) {

tests/testthat/test-LANDISDisp-threads.R

@@ -0,0 +1,20 @@
+## Thread count must not change LANDISDisp results.
+
+test_that("LANDISDisp output is identical for threads = 1, 2, 4 and the R reference", {
+  skip_if_not_installed("terra")
+  ## the OpenMP path only engages above ~20000 active receivers per step
+  fix <- makeLANDISDispFixture(size = "large", fixtureSeed = 11L, successionTimestep = 10L)
+  ompAvail <- LandR:::landisDispHasOpenMP()
+
+  run <- function(nThr, useCpp = TRUE) {
+    withr::with_options(list(LandR.LANDISDisp.threads = nThr),
+                        runLANDISDispOnFixture(fix, runSeed = 42L, useCpp = useCpp))
+  }
+  ref <- run(1L)
+  expect_gt(nrow(ref), 0L)
+  expect_identical(run(1L, useCpp = FALSE), ref)
+  for (n in c(2L, 4L)) {
+    skip_if(!ompAvail, "OpenMP not available")
+    expect_identical(run(n), ref, info = paste("threads =", n))
+  }
+})