Reproducibility: RNG and Seed Conventions

The default case: seed sets R’s own RNG state once

Every Design constructor accepts a seed argument. Internally this does not call set.seed() immediately; it stores private$seed and calls private$maybe_set_seed() — if (!is.null(private$seed)) set.seed(private$seed) — once, immediately before each call that consumes randomness (draw_ws_raw()/assign_w_to_all_subjects()). This means:

• Two designs constructed with the same seed and then drawn from once each produce identical allocations.
• Calling a draw method a second time on the same object re-seeds again (since maybe_set_seed() runs before every draw), so it reproduces the first draw again rather than advancing to a fresh one — draws are not incremented automatically across repeated calls on one object. If you want r independent replicate columns, request them in one call (draw_ws_according_to_design(r)), not via r separate single-draw calls.
• seed = NULL (the default) leaves R’s ambient RNG state untouched — draws consume whatever state R’s global stream happens to be in, exactly like any other call to sample()/runif().

This is the mechanism behind the overwhelming majority of Design subclasses (DesignFixedBernoulli, DesignFixedFactorial, DesignFixedBlocking, the DesignSeqOneByOne* family, etc.). Three design families implement “seed means R’s set.seed() governs the draw” via a different, local-generator mechanism, documented in detail below: DesignFixedGreedy (parallel-safe per-thread generators), DesignFixedGreedyDOptimal (a single local generator seeded from R’s stream), and DesignFixedOptimal (per-chain generators for its annealing solver plus an R-level mirror coin; see its own subsection below). Every concrete design is seed-reproducible; the package currently has no exceptions.

DesignFixedOptimal: deterministic solves, seeded coin and chains

DesignFixedOptimal computes one optimal allocation rather than drawing from a randomization distribution, so most of its “draw” is not random at all:

• The exact "ompr" MILP path is fully deterministic — same data, same arguments, same allocation, no RNG consumed by the solve itself.
• The mirror coin is R-seeded. With mirror_coin = TRUE (the default) and a verified co-optimal mirror at prob_T = 0.5, one runif(1) draw from R’s live stream decides between w* and 1 - w*. Under the constructor’s seed, maybe_set_seed() runs before the solve, so the flip is reproducible; the MILP path is therefore “deterministic up to the seeded label flip.”
• The annealing solver seeds one edi_rng::RRng per chain (not per thread), each from one R::unif_rand() draw before the parallel region — so a given set.seed() reproduces the identical search under any OpenMP thread count, a deliberately stronger guarantee than DesignFixedGreedy’s per-thread seeding. When initial_temp is auto-calibrated, the calibration probe also consumes R’s stream (and is therefore covered by the same seed). Annealing carries certificate "annealing_converged", never "global": Hajek (1988) guarantees convergence in probability only under a logarithmic cooling schedule, and the practical geometric schedule used here is asymptotically motivated, not a finite-time proof.
• BRT replicates replay the coin-inclusive mechanism: each replicate’s re-optimization (reduced annealing by default, solver_args$brt_*) and its own mirror flip consume the worker’s seeded stream, so bootstrap randomization p-values are reproducible under set.seed() like every other BRT in the package.

The portable cross-language RNG: edi_rng::RRng

Most C++ backends that need bulk random draws (Pocock-Simon minimization’s pocock_simon_assign_cpp, the bootstrap-index generators bootstrap_indices.cpp/bootstrap_match_indices.cpp, weighted-distance sampling, and others) do not call back into R’s unif_rand() for every individual draw — that has real per-call overhead. Instead they use a two-step pattern:

1. Draw exactly one value from R’s live RNG stream via R::unif_rand(), and convert it to a 32-bit integer seed (edi_rng::seed_from_unif01()).
2. Seed a fresh edi_rng::RRng instance (RNG.h) from that one integer, and do every subsequent draw against that instance instead of R’s own generator.

edi_rng::RRng is a portable, from-scratch reimplementation of R’s own Mersenne-Twister + Inversion generator (not a wrapper around it) — the same generator, byte-for-byte, callable from C++ or Python without linking against R’s runtime. This is what makes a given seed produce identical draws in R and in the Python bindings (edi_kernels) using the same core, and is why individual function docstrings describe this as “seeded from one R::unif_rand() draw into edi_rng::RRng.”

Consequence worth internalizing: because only one unif_rand() value is consumed from R’s stream per call, R’s own .Random.seed afterward has advanced by exactly one draw, not by however many draws happened inside the C++ function — so R-level code interleaving several such calls still gets a well-defined, reproducible (if not obviously predictable) sequence of R-level draws in between them.

Seed-reproducible via a local generator: DesignFixedGreedyDOptimal

Historical note: the two classes this design merges (DesignFixedAOptimal/DesignFixedDOptimal) originally seeded their exchange-search kernels’ initial random shuffles from std::random_device (hardware entropy), which made their draws genuinely not reproducible via seed. The RNG migration (see the SEXP-removal spec’s RNG section) replaced that with a local edi_rng::RRng generator seeded from R’s own stream (R::unif_rand()) inside d_optimal_search_cpp()/a_optimal_search_cpp(), and the merged DesignFixedGreedyDOptimal class inherits that behavior: repeated calls with the same seed return identical allocations (verified empirically in test-greedy-d-optimal-merged.R, and reflected in the seed_reproducible_draw registry field). Because randomization-based inference (compute_rand_two_sided_pval(), compute_rand_confidence_interval()) generates its reference distribution by calling the design’s own draw_ws_according_to_design() (see below), randomization p-values/CIs against this design are exactly reproducible with seed set, like every other design’s.

Seed-reproducible despite being parallel: DesignFixedGreedy

DesignFixedGreedy‘s search (greedy_design_search_cpp()) runs r independent searches in parallel via OpenMP, each with its own std::mt19937 generator — but unlike the A-/D-optimal kernels, these per-thread generators are seeded from R’s own RNG state (GetRNGstate()/unif_rand()) before the parallel region begins, so private$maybe_set_seed() does govern the resulting allocation, and the result is identical regardless of how many OpenMP threads (RhpcBLASctl/ set_num_cores()) are actually used at draw time. This is the template other designs’ parallel kernels should follow if they need both speed and seed-reproducibility simultaneously — draw all per-worker seeds from R’s stream up front, before fanning out.

Randomization inference reuses the design’s own draw mechanism

InferenceRand’s generate_permutations(r) does not implement its own permutation-drawing logic; it duplicates the design object (des_obj$duplicate()) and calls that duplicate’s draw_ws_according_to_design(r) — the exact same entry point assign_w_to_all_subjects() uses. Consequences:

• RNG/seed-reproducibility of a randomization p-value or CI is exactly whatever the underlying Design subclass’s is — and every concrete design is currently seed-reproducible (see above).
• Draw reuse / caching: the generated permutation matrix is cached, keyed on r and a stable signature of the design’s structural parameters (class, n, prob_T, m, strata_cols). A second call for the same r against a structurally-identical design reuses the cached matrix rather than drawing again — so, for example, computing both a randomization p-value and a randomization confidence interval (which internally makes several p-value evaluations at different delta) against the same object draws the reference permutation set once, not once per evaluation.

Bootstrap resampling

• Non-parametric bootstrap: the default fallback is sample.int(n, n, replace = TRUE) — plain R-level sampling, governed by R’s ambient RNG state / set.seed() exactly like any base-R code. Block/pair/cluster-structured designs instead call des_obj$draw_bootstrap_indices(bootstrap_type), which for most designs routes to the one-draw-seeded edi_rng::RRng C++ backends described above (draw_matching_bootstrap_sample_cpp, stratified_bootstrap_indices_cpp, resample_group_rows_cpp) — reproducible via set.seed() acting on R’s stream at the point the one seeding draw is taken.
• Bayesian bootstrap: weights are drawn as stats::rgamma(length(idx), shape = 1, rate = 1), one Gamma(1,1) draw per exchangeable resampling unit (subject or block), which is the standard construction of Dirichlet\((1,\dots,1)\) weights (a Dirichlet draw is a Gamma\((1,1)\) vector normalized to sum to the unit count — see vignette("notation-glossary")’s “Resampling and randomization” section). This is plain R-level rgamma(), governed by R’s ambient RNG state.
• Parametric bootstrap / warm-start / factorization reuse: several inference classes cache a factorization or warm-start state across bootstrap replicates purely for speed (see get_warm_start_dispatch_policy()/set_warm_start_dispatch_policy()); this caching does not change which random draws are made, only how fast each replicate’s model fit converges.

Machine-dependent performance defaults: tune_EDI_for_this_machine()

Every performance-policy default mentioned above — whether an inference class’s C++ backend uses a smart_cold_start OLS warm-up or a plain zero start, whether resampling reuses a previous replicate’s warm start (and at what sample size that stops paying off), which optimizer algorithm a family uses by default, and at what sample size parallel bootstrapping starts to beat serial execution — was measured empirically on the maintainer’s machine. These are speed judgments, not statistical ones: core count, cache sizes, and BLAS backend all affect which setting wins, so a default that is net-positive on the maintainer’s machine can be net-negative on yours, and vice versa.

tune_EDI_for_this_machine() re-runs those same benchmarks on your own hardware, keeps only the settings that win by a real margin (median improvement past a noise threshold, not any transient win), and persists the result to a per-user config file that every subsequent library(EDI) re-applies automatically. This does not weaken anything this vignette documents: tuning never changes which random draws are made or which estimate/CI a fit produces — only how fast it gets there. Concretely, every accepted change is re-fit once under both settings before being kept, and any disagreement in the result discards that change rather than applying it (the same “measure, don’t assume” discipline this vignette applies to RNG behavior, applied here to timing behavior). The one axis with a partial exception is the parallel/core-count benchmark itself, since forked workers draw from an independent RNG stream by construction — that axis compares the (core-count-invariant) point estimate instead of the resampling distribution, for exactly the reason a bootstrap CI is expected to differ across independent Monte Carlo draws.

To see what has been tuned on your machine, call get_local_EDI_optimization(); to discard it and return to the shipped defaults, call clear_local_EDI_optimization(). Setting EDI_SKIP_LOCAL_TUNING=1 (e.g. before library(EDI)) skips the automatic import for one session without deleting the saved file — useful when isolating whether a saved tuning is responsible for an observed timing difference.

Simulation (SimulationFramework): per-replication and per-cache-job seeds

SimulationFramework$new(seed = ...) does not rely on a single global set.seed() call covering the entire run. At the top of run(), if seed is non-NULL, R’s .Random.seed is saved (to be restored via on.exit() when run() returns, so a simulation run never leaks RNG state into the caller’s session) and set.seed(private$seed) is called once — but the more important mechanism is per-unit deterministic seed derivation:

• Each replication i (within a w-rep loop) is dispatched with rep_seed = seed + i, and the worker executing that replication calls set.seed(rep_seed) itself before drawing anything.
• Each cache-building job (pre-generating the design/SE caches used across replications) is dispatched with cache_seed = seed + 1000003L + job_idx — a large additive offset specifically chosen so the cache-job seed range and the replication seed range do not collide for any realistic Nrep_W/ cell count.

Why this matters for parallel execution: because every unit of work carries its own explicit seed and calls set.seed() itself, the ambient RNG state of whichever worker process executes it is irrelevant — this is what makes a SimulationFramework run reproducible regardless of num_cores, regardless of whether the fork-cluster or mirai-daemon backend is used, and regardless of the order in which the scheduler happens to dispatch replications/cache jobs across workers. Each saved on-disk cache record additionally stores the RNG state present right after that cache object was built (rng_after); restore_rng (default FALSE on cache loads) controls whether a cache hit replays that saved state into the current session. Since a cache hit skips the computation that would have consumed that randomness anyway, leaving restore_rng = FALSE is the correct default — cache hits are RNG-inert (they neither consume nor need to replay randomness), while a cache miss (an actual fresh .run_simulation_cache_job() call) still explicitly seeds itself via the same seed + 1000003L + job_idx derivation as any other cache job.

Consequence: a SimulationFramework run’s results are reproducible across separate run() invocations with the same seed and the same (design_classes_and_params, inference_classes_and_params, n, p, betaT, ...) configuration. (Historically this carried an exception for the two pre-merge optimal-design classes, whose kernels were not seeded from R’s stream; since the RNG migration and the DesignFixedGreedyDOptimal merge, no such exception exists.)

Monte Carlo error

None of B_boot, r_rand, or Nrep_W/Nrep_Y_w has a closed-form “this value is large enough” answer baked into the package — larger values reduce simulation noise at the cost of runtime, and the right value is estimand/design-specific. Rules of thumb used elsewhere in statistics apply directly here:

• A randomization or bootstrap p-value built from r/B draws has Monte Carlo standard error on the order of \(\sqrt{\hat p (1-\hat p) / r}\) (treating “did this draw’s statistic exceed the observed one” as a Bernoulli(\(\hat p\)) indicator) — e.g. r = 999 gives a Monte Carlo SE of roughly \(0.016\) at \(\hat p \approx 0.5\), tighter near the tails that usually matter for a decision at \(\alpha = 0.05\).
• A bootstrap confidence interval’s endpoints (percentile or BCa) are themselves noisy quantile estimates from B draws; their Monte Carlo error shrinks roughly like \(O(1/\sqrt{B})\), but unlike the p-value case there is no single clean formula — the standard practical guidance is to re-run with a different seed and confirm the interval doesn’t move appreciably before trusting a B_boot choice for a final reported result.
• A SimulationFramework operating characteristic (MSE, coverage, power/size in SimulationFrameworkReport$summarize()) is itself a Monte Carlo estimate over Nrep_W * Nrep_Y_w replications; coverage_pval/size_pval (exact two-sided binomial test p-values against the nominal \(1-\alpha\)/\(\alpha\) target) are provided specifically so a large enough Nrep_W can be chosen to distinguish “genuinely miscalibrated” from “within Monte Carlo noise of nominal” rather than eyeballing a point estimate.

Reproducing a documented example

To exactly reproduce a design allocation, bootstrap replicate, or simulation run shown in this package’s own examples/vignettes/published comparisons:

1. Use the same R version and the default RNGkind() (c("Mersenne-Twister", "Inversion", "Rejection")) — the portable edi_rng::RRng reimplementation and pocock_simon_redraw_w_cpp’s live-stream continuation both assume the Mersenne-Twister + Inversion normal-sampling kind specifically.
2. Pass an explicit seed to the Design/SimulationFramework constructor rather than relying on ambient RNG state, and do not call any other RNG-consuming code between construction and the draw you want to reproduce (anything that advances R’s global stream in between — including, per above, a second call to a draw method on the same object — changes what gets drawn next).
3. For parallel SimulationFramework runs, num_cores and the fork/mirai backend choice do not need to match the original run for reproducibility (per-unit seed derivation makes them irrelevant) — only seed and the simulation configuration do.