set.seed() en reproduceerbaarheid

# Check current RNG settings
RNGkind()
# [1] 'Mersenne-Twister' 'Inversion' 'Rejection'
# [1] kind, normal.kind, sample.kind

# Set explicitly for maximum reproducibility
RNGkind('Mersenne-Twister', 'Inversion', 'Rejection')
set.seed(42)
runif(3)

# Alternative RNG algorithms (rarely needed)
# RNGkind('L'Ecuyer-CMRG') # for parallel RNG
# RNGkind('Super-Duper')   # older algorithm

cat('Default Mersenne-Twister has period 2^19937 - 1')

set.seed(42)
# .Random.seed is created in .GlobalEnv after first use
length(.Random.seed)  # 626 integers for Mersenne-Twister
.Random.seed[1]       # encodes RNG type

# Generate one number (state advances)
x1 <- runif(1)
state_after <- .Random.seed

# Generate another
x2 <- runif(1)

# Restore state to after x1
.Random.seed <<- state_after
x2_replay <- runif(1)

identical(x2, x2_replay)  # TRUE - we replayed the state
cat('x2:', x2, 'x2_replay:', x2_replay)

set.seed(100)
# Generate some random numbers first
runif(10)

# Save current state
saved_state <- .Random.seed

# Generate block of interest
block1 <- rnorm(5)
print(block1)

# ... later, restore and replay
.Random.seed <<- saved_state
block1_replay <- rnorm(5)

identical(block1, block1_replay)  # TRUE!
cat('RNG state saved and restored successfully')

# Set up L'Ecuyer-CMRG for parallel use
RNGkind('L\'Ecuyer-CMRG')
set.seed(42)

# Each parallel worker gets its own stream
# (parallel package handles this automatically)
# library(parallel)
# cl <- makeCluster(4)
# clusterSetRNGStream(cl, iseed = 42)

# Generate numbers with this generator
samples <- rnorm(5)
print(samples)

# Reset to default Mersenne-Twister
RNGkind('Mersenne-Twister', 'Inversion', 'Rejection')
cat('L\'Ecuyer-CMRG: safe for parallel simulation')

# Function that optionally accepts a seed
simulate_data <- function(n, seed = NULL) {
  if (!is.null(seed)) set.seed(seed)
  list(
    x = rnorm(n),
    y = rnorm(n)
  )
}

# Reproducible call
result1 <- simulate_data(5, seed = 42)
result2 <- simulate_data(5, seed = 42)
identical(result1$x, result2$x)  # TRUE

# Without seed: different each time
result3 <- simulate_data(5)
result4 <- simulate_data(5)
identical(result3$x, result4$x)  # FALSE

# R 3.6+ default sample kind
RNGkind(sample.kind = 'Rejection')  # default
set.seed(42)
sample(1:10, 5)
# [1] 1 5 10 8 2  (with Rejection)

# Old behavior (R < 3.6) for legacy code
RNGkind(sample.kind = 'Rounding')
set.seed(42)
sample(1:10, 5)
# Different result!

# Best practice: reset to modern defaults
RNGkind('Mersenne-Twister', 'Inversion', 'Rejection')
set.seed(42)
sample(1:10, 5)

# Reproducibility header for analysis scripts:
# R version: R.version$version.string
# Seed: 42
# RNGkind: Mersenne-Twister / Inversion / Rejection

# Capture R version
R.version$version.string

# Capture full session info
# sessionInfo()  # shows all packages and versions

# In RMarkdown, set seed in setup chunk:
# ```{r setup}
# knitr::opts_chunk$set(echo = TRUE)
# set.seed(42)
# ```

# This ensures all chunks share the same RNG state
set.seed(42)
cat('R version:', R.version$version.string, '\n')
cat('Seed: 42, RNGkind: Mersenne-Twister')

# Generate per-replicate seeds from a master seed
set.seed(2024)
n_replicates <- 5
replicate_seeds <- sample.int(1e6, n_replicates)
print(replicate_seeds)

# Each replicate uses its own seed
run_replicate <- function(rep_id) {
  set.seed(replicate_seeds[rep_id])
  mean(rnorm(1000))  # estimate with this seed
}

results <- sapply(seq_len(n_replicates), run_replicate)
print(round(results, 4))

# Any single replicate can be re-run independently:
set.seed(replicate_seeds[3])
mean(rnorm(1000))  # exactly matches results[3]

# withr::with_seed - restore RNG state automatically
# library(withr)
# result <- with_seed(42, rnorm(5))

# Simulation without withr (manual state save)
old_seed <- if (exists('.Random.seed')) .Random.seed
old_kind <- RNGkind()

set.seed(42)
result <- rnorm(5)

# Restore
if (!is.null(old_seed)) .Random.seed <<- old_seed
cat('Result:', result, '\n')
cat('Caller\'s RNG state preserved')

# set.seed in tests ensures deterministic assertions:
set.seed(1)
stopifnot(round(rnorm(1), 6) == round(-0.6264538, 6))

# Reproducibility checklist:
# 1. Set seed at script top
set.seed(42)

# 2. Document RNGkind
RNGkind()  # Mersenne-Twister / Inversion / Rejection

# 3. Save RNG state if needed
saved <- .Random.seed

# 4. For functions, accept seed argument
my_sim <- function(n, seed = NULL) {
  if (!is.null(seed)) set.seed(seed)
  rnorm(n)
}

# 5. For parallel: use L'Ecuyer-CMRG
# RNGkind('L\'Ecuyer-CMRG'); set.seed(42)

cat('Reproducibility = trust in your results')