Merge pull request #447 from UCD-SERG/sim-pop-multi-no-skip

try to not skip

Author: d-morrison

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: serocalculator
 Title: Estimating Infection Rates from Serological Data
-Version: 1.4.0.9015
+Version: 1.4.0.9016
 Authors@R: c(
     person("Kristina", "Lai", , "kwlai@ucdavis.edu", role = c("aut", "cre")),
     person("Chris", "Orwa", role = "aut"),

NEWS.md

@@ -16,6 +16,13 @@
 
 ## Bug fixes
 
+* `sim_pop_data()` and `sim_pop_data_multi()` now produce identical results
+  across operating systems. Simulated inter-infection times are now rounded to
+  whole days, so the number of random draws consumed no longer depends on
+  platform-specific floating-point results of `log()` (which previously
+  shifted the random-number stream out of sync and made simulated values, and
+  their snapshots, differ between macOS, Windows, and Linux). Simulated
+  values change slightly as a result of this fix. (#447)
 * Removed lingering terminology discrepancies
 * `load_noise_params()` and `load_sr_params()` now fail gracefully with informative messages when internet resources are unavailable, complying with CRAN policy (#505)
 * Added Version Crosswalk article to pkgdown website to help users migrate code from v1.3.0 to v1.4.0

R/quantize_t_next.R

@@ -0,0 +1,20 @@
+#' Quantize an inter-infection time to whole days
+#'
+#' Rounds a continuous inter-infection time (in days) to the nearest whole
+#' day, with a floor of one day. The antibody-kinetics simulation in
+#' `simresp.tinf()` already operates on a daily grid, but the number of
+#' `while()`-loop iterations (and hence the number of random draws consumed)
+#' otherwise depends on the platform-specific floating-point result of
+#' `log()`. Quantizing to whole days keeps random-number consumption
+#' identical across operating systems, so simulated data and its snapshots
+#' are reproducible on macOS, Windows, and Linux.
+#'
+#' The floor of one day prevents a zero-length step, which would otherwise
+#' stall the `while()` loop in `simresp.tinf()`.
+#'
+#' @param t_next a positive [numeric] scalar inter-infection time, in days
+#' @returns a whole-number-valued [numeric] scalar, at least 1
+#' @noRd
+quantize_t_next <- function(t_next) {
+  max(1, round(t_next))
+}

R/simresp.tinf.R

@@ -42,7 +42,7 @@ simresp.tinf <- function(# nolint: object_name_linter
   mcsize <- dim(predpar)[3]
   nmc <- n_mcmc_samples
 
-  day2yr <- 365.25
+  days_per_year <- 365.25
 
   if (n_mcmc_samples == 0) {
     nmc <- sample.int(n = mcsize, size = 1)
@@ -59,8 +59,19 @@ simresp.tinf <- function(# nolint: object_name_linter
   t_inf <- c()
 
   t_next <- -log(runif(1, 0, 1)) / lambda # time to first infection...
-
-  age <- if_else(!is.na(age_fixed), age_fixed, t_next / day2yr)
+  # Quantize inter-infection times to whole days. The simulation already
+  # operates on a daily grid (`seq(0, t_next, by = 1)`, and ages are sampled
+  # over integer days), but the `while()` loop boundary and the per-episode
+  # grid lengths otherwise depend on the platform-specific floating-point
+  # result of `log()`. Those ULP-level differences flip the number of loop
+  # iterations on some operating systems, which changes how many random
+  # draws are consumed and desynchronizes the (otherwise platform-independent)
+  # L'Ecuyer-CMRG stream -- making snapshot output irreproducible across
+  # macOS/Windows/Linux. Rounding to whole days keeps RNG consumption
+  # identical across platforms. See `quantize_t_next()`.
+  t_next <- quantize_t_next(t_next)
+
+  age <- if_else(!is.na(age_fixed), age_fixed, t_next / days_per_year)
 
   mcpar <- ldpar(
     age = age,
@@ -116,7 +127,7 @@ simresp.tinf <- function(# nolint: object_name_linter
 
     if (!renew_params) {
       par_now <- ldpar(
-        if (!is.na(age_fixed)) age_fixed else t0 / day2yr,
+        if (!is.na(age_fixed)) age_fixed else t0 / days_per_year,
         antigen_isos,
         nmc,
         predpar = predpar, ...
@@ -129,6 +140,10 @@ simresp.tinf <- function(# nolint: object_name_linter
     # note: the renew_params == TRUE case is handled many lines below
 
     t_next <- -log(runif(1, 0, 1)) / lambda
+    # Quantize to whole days for cross-platform reproducibility (see the
+    # comment on the first-infection draw above). `t_end` and `t0` are whole
+    # days, so the boundary clamp below stays integer-valued too.
+    t_next <- quantize_t_next(t_next)
     if (t0 <= t_end && t0 + t_next > t_end) {
       t_next <- t_end - t0
     }
@@ -150,7 +165,7 @@ simresp.tinf <- function(# nolint: object_name_linter
     b <- rbind(b, b_now)
     y_mat <- rbind(y_mat, y_now)
 
-    y_end <- y_mat %>% tail(1)
+    y_end <- y_mat |> tail(1)
 
     if (renew_params) {
       if (n_mcmc_samples == 0) {
@@ -163,7 +178,7 @@ simresp.tinf <- function(# nolint: object_name_linter
     age <- if_else(
       !is.na(age_fixed),
       age_fixed,
-      (t0 + t_next) / day2yr
+      (t0 + t_next) / days_per_year
     )
 
     par_now <- ldpar(