### Setup: Parameters and Logging

Source: https://github.com/stephenslab/susier/blob/master/inst/notebooks/stochastic_ld_benchmark.ipynb

Initializes SuSiE-R, sets up simulation parameters, and configures logging to a file for progress monitoring.

```R
# ============================================================ 
# Setup: parameters and logging
# ============================================================ 
library(susieR)

if (requireNamespace("Rfast", quietly = TRUE)) {
  library(Rfast)
  fast_cor <- Rfast::cora
  cat("Using Rfast::cora for correlation matrices (much faster)\n")
} else {
  fast_cor <- cor
  cat("Rfast not available, using base cor() (consider installing Rfast)\n")
}

# --- Simulation parameters ---
p         <- 5000
n         <- 100000
n_causal  <- 4
pve_per   <- 0.0012      # PVE per causal
B_values  <- c(2000, 5000, 10000)
n_rep     <- 50         # set to 250 for full run
L         <- 10
seed      <- 999

# --- Log file for monitoring progress from the terminal ---
# When running via jupyter nbconvert, cat() output is captured into the
# notebook and NOT shown on screen. To monitor progress, open a second
# terminal and run:   tail -f stochastic_ld_benchmark.log
LOG_FILE <- "stochastic_ld_benchmark.log"
log_con  <- file(LOG_FILE, open = "w")

log_msg <- function(fmt, ...) {
  msg <- sprintf(paste0("[%s] ", fmt, "\n"), format(Sys.time(), "%H:%M:%S"), ...)
  cat(msg)                           # notebook cell output
  cat(msg, file = log_con)           # log file for tail -f
  flush(log_con)
}

set.seed(seed)
log_msg("Stochastic LD Benchmark started")
log_msg("  susieR version: %s", packageVersion("susieR"))
log_msg("  Parameters: p=%d, n=%d, n_causal=%d, pve_per=%.4f, n_rep=%d",
        p, n, n_causal, pve_per, n_rep)
log_msg("  B values: %s", paste(B_values, collapse = ", "))
log_msg("  Log file: %s  (run 'tail -f %s' to monitor)", LOG_FILE, LOG_FILE)

```

--------------------------------

### Output Directory Setup

Source: https://github.com/stephenslab/susier/blob/master/inst/notebooks/small_sample_benchmark.ipynb

Sets up the output directory for storing simulation results. All results are consolidated into a single RDS file.

```R
outdir  <- "benchmark_results"
all_rds <- file.path(outdir, "all_results.rds")
```

--------------------------------

### Build Package Website with pkgdown

Source: https://github.com/stephenslab/susier/blob/master/README.md

Run this command to build the package website using pkgdown. Ensure you have pkgdown installed.

```bash
make pkgdown
```

--------------------------------

### Install susieR from GitHub

Source: https://github.com/stephenslab/susier/blob/master/README.md

Install the latest development version of the susieR package directly from its GitHub repository. Ensure the 'remotes' package is installed first.

```R
# install.packages("remotes")
remotes::install_github("stephenslab/susieR")
```

--------------------------------

### Initialize SusieR Benchmark Configuration

Source: https://github.com/stephenslab/susier/blob/master/inst/notebooks/small_sample_benchmark.ipynb

Sets up parallel processing and defines simulation parameters for the SusieR benchmark. Ensure 'parallelly' is installed for core availability.

```R
library(susieR)
library(glmnet)
library(digest)
library(future)
library(future.apply)

# --- Configuration ---
ncores    <- max(1, parallelly::availableCores() - 2)
n_rep     <- 200       # replicates per setting
L         <- 10
N_vals    <- c(30, 50, 70, 100)
h2_sparse <- c(0.25, 0.50, 0.75)
L_causal  <- c(1, 2, 3, 4, 5)

# Use multisession (PSOCK) to avoid fork + BLAS crashes
plan(multisession, workers = ncores)

cat(sprintf("susieR version : %s\n", packageVersion("susieR")))
cat(sprintf("Workers        : %d  (multisession / PSOCK)\n", ncores))
cat(sprintf("N values       : %s\n", paste(N_vals, collapse = ", ")))
cat(sprintf("h2_sparse      : %s\n", paste(h2_sparse, collapse = ", ")))
cat(sprintf("n_causal       : %s\n", paste(L_causal, collapse = ", ")))
cat(sprintf("Settings       : %d\n", length(N_vals) * length(h2_sparse) * length(L_causal)))
cat(sprintf("Total reps     : %d\n",
    length(N_vals) * length(h2_sparse) * length(L_causal) * n_rep))
```

--------------------------------

### Simulation Setup and Checkpointing

Source: https://github.com/stephenslab/susier/blob/master/inst/notebooks/small_sample_benchmark.ipynb

Defines functions for managing simulation seeds and checkpointing results. Ensures reproducibility and efficient continuation of simulation runs.

```R
cat("Functions defined:\n")
cat("  Seeds:       make_seed(rep_i, N, purpose, h2, nc)\n")
cat("  Checkpoint:  compute_config_md5, checkpoint_init, checkpoint_completed_reps,\n")
cat("               checkpoint_load, checkpoint_save, checkpoint_save_meta\n")
cat("  Simulation:  get_lasso_residual, gen_signflip_noise, run_one_rep\n")
```

--------------------------------

### Install susieR from CRAN

Source: https://github.com/stephenslab/susier/blob/master/README.md

Use this command to install the stable version of the susieR package from the Comprehensive R Archive Network (CRAN).

```R
install.packages("susieR")
```

--------------------------------

### checkpoint_init

Source: https://github.com/stephenslab/susier/blob/master/inst/notebooks/small_sample_benchmark.ipynb

Initializes the checkpoint directory. It verifies if existing results match the current configuration's MD5 hash. If they mismatch, existing results are cleared.

```APIDOC
## checkpoint_init()

### Description
Initializes the checkpoint directory. It verifies existing results against the provided configuration signature. If a mismatch is detected, all existing .rds files in the directory are deleted.

### Parameters
- **outdir** (string) - The directory path for storing checkpoints.
- **config_sig** (list) - A list containing the configuration and its MD5 hash, typically generated by `compute_config_md5()`.

### Returns
A list with the following elements:
- **status**: A string, either "match" if existing results are valid, or "fresh" if new results need to be generated.
- **meta**: The metadata from the previous run if status is "match", otherwise NULL.
```

--------------------------------

### Monitor Progress

Source: https://github.com/stephenslab/susier/blob/master/inst/notebooks/stochastic_ld_benchmark.ipynb

Provides instructions on how to monitor the simulation progress in real-time by tailing the log file.

```bash
tail -f stochastic_ld_benchmark.log
```

--------------------------------

### Get LASSO Residuals

Source: https://github.com/stephenslab/susier/blob/master/inst/notebooks/small_sample_benchmark.ipynb

Calculates residuals using LASSO regression. Useful for estimating noise components in simulation studies.

```R
get_lasso_residual <- function(X, y, seed = 42) {
  set.seed(seed)
  n <- nrow(X)
  cv_fit    <- cv.glmnet(X, y, alpha = 1, nfolds = min(10, n))
  lasso_fit <- glmnet(X, y, alpha = 1, lambda = cv_fit$lambda.1se)
  as.vector(y - predict(lasso_fit, X))
}
```

--------------------------------

### Initialize Checkpoint Directory

Source: https://github.com/stephenslab/susier/blob/master/inst/notebooks/small_sample_benchmark.ipynb

Initializes the checkpoint directory, verifying existing results against the new configuration's MD5 hash. If a mismatch occurs, old results are cleared.

```R
checkpoint_init <- function(outdir, config_sig) {
  dir.create(outdir, showWarnings = FALSE, recursive = TRUE)
  meta_path <- file.path(outdir, "_meta.rds")

  if (file.exists(meta_path)) {
    old_meta <- tryCatch(readRDS(meta_path), error = function(e) NULL)
    if (!is.null(old_meta) && !is.null(old_meta$md5) &&
        identical(old_meta$md5, config_sig$md5)) {
      cat(sprintf("Config MD5 verified: %s\n", config_sig$md5))
      cat(sprintf("  Previous run: n_rep=%d, timestamp=%s\n",
          old_meta$n_rep, old_meta$timestamp))
      return(list(status = "match", meta = old_meta))
    }
    old_md5 <- if (!is.null(old_meta$md5)) old_meta$md5 else "(missing/corrupt)"
    cat(sprintf("Config MD5 MISMATCH — clearing old results.\n"))
    cat(sprintf("  Old: %s\n  New: %s\n", old_md5, config_sig$md5))
    old_files <- list.files(outdir, pattern = "[.]rds$", full.names = TRUE)
    if (length(old_files) > 0) {
      file.remove(old_files)
      cat(sprintf("  Removed %d old file(s).\n", length(old_files)))
    }
    return(list(status = "fresh", meta = NULL))
  }

  cat(sprintf("No existing checkpoint. Config MD5: %s\n", config_sig$md5))
  list(status = "fresh", meta = NULL)
}
```

--------------------------------

### Get Completed Replication IDs

Source: https://github.com/stephenslab/susier/blob/master/inst/notebooks/small_sample_benchmark.ipynb

Retrieves a sorted vector of completed replication IDs for a given tag. Returns an empty integer vector if the file does not exist or is corrupt.

```R
checkpoint_completed_reps <- function(outdir, tag) {
  path <- file.path(outdir, paste0(tag, ".rds"))
  if (!file.exists(path)) return(integer(0))
  tryCatch({
    x <- readRDS(path)
    if (is.data.frame(x) && "rep" %in% names(x) && nrow(x) > 0)
      sort(unique(as.integer(x$rep)))
    else
      integer(0)
  }, error = function(e) integer(0))
}
```

--------------------------------

### Initialize SuSiE with susie_init_coef

Source: https://context7.com/stephenslab/susier/llms.txt

Creates a susie initialization object from known or estimated non-zero coefficient locations and values. Useful for warm-starting model fitting. Requires the 'susieR' library.

```r
library(susieR)
set.seed(1)
n <- 1000; p <- 1000
beta <- rep(0, p)
beta[sample(1:p, 4)] <- 1
X <- matrix(rnorm(n * p), nrow = n, ncol = p)
X <- scale(X, center = TRUE, scale = TRUE)
y <- drop(X %*% beta + rnorm(n))

# Initialize susie at true coefficient locations
true_idx   <- which(beta != 0)
true_vals  <- beta[true_idx]
s_init <- susie_init_coef(
  coef_index = true_idx,
  coef_value = true_vals,
  p          = p
)

# Warm-start susie from ground truth
fit_warmstart <- susie(X, y, L = 10, model_init = s_init)

# Compare to cold-start
fit_cold <- susie(X, y, L = 10)
cat("Warm-start iterations:", fit_warmstart$niter, "\n")
cat("Cold-start iterations:", fit_cold$niter, "\n")
```

--------------------------------

### Run Single Method for Benchmark

Source: https://github.com/stephenslab/susier/blob/master/inst/notebooks/stochastic_ld_benchmark.ipynb

Executes a single fine-mapping method (e.g., susie_rss) for a given replicate. It handles potential errors during fitting and logs performance metrics like elapsed time, power, FDR, number of credible sets (CS), and mean CS size. Returns a data frame with results for the replicate.

```R
run_one_method <- function(method_name, z, causal_idx, n, L, ...) {
  t0 <- proc.time()
  fit <- tryCatch(
    susie_rss(z = z, n = n, L = L, max_iter = 200, verbose = FALSE, ...),
    error = function(e) { log_msg("    %s ERROR: %s", method_name, e$message); NULL }
  )
  if (is.null(fit)) return(NULL)
  ev <- evaluate_fit(fit, causal_idx)
  elapsed <- (proc.time() - t0)[3]
  # Compute per-replicate rates for the log line
  rep_power <- ev$n_found / ev$n_causal
  rep_fdr   <- if (ev$n_cs > 0) ev$n_false_cs / ev$n_cs else 0
  log_msg("    %-25s %5.1fs | power=%.2f  fdr=%.2f  n_cs=%d  size=%.0f",
          method_name, elapsed, rep_power, rep_fdr, ev$n_cs, ev$mean_size)
  data.frame(method = method_name,
             n_cs = ev$n_cs, n_true_cs = ev$n_true_cs, n_false_cs = ev$n_false_cs,
             n_found = ev$n_found, n_causal = ev$n_causal,
             mean_size = ev$mean_size,
             stringsAsFactors = FALSE)
}
```

--------------------------------

### Build LD Structure

Source: https://github.com/stephenslab/susier/blob/master/inst/notebooks/stochastic_ld_benchmark.ipynb

Initializes the LD structure, which is a prerequisite for genotype and phenotype simulation. This step is computationally fast.

```R
# ============================================================ 
# Phase 1: Build LD structure
# ============================================================ 
log_msg("Phase 1: Building LD structure (p=%d)...", p)
ld_struct <- build_ld_structure(p, seed = seed)
log_msg("  %d blocks, sizes %d-%d, rho %.2f-%.2f",
        ld_struct$n_blocks,
        min(sapply(ld_struct$blocks, function(b) b$size)),
        max(sapply(ld_struct$blocks, function(b) b$size)),
        min(sapply(ld_struct$blocks, function(b) b$rho)),
        max(sapply(ld_struct$blocks, function(b) b$rho)))
```

--------------------------------

### Execute Notebook and Monitor Progress

Source: https://github.com/stephenslab/susier/blob/master/inst/notebooks/stochastic_ld_benchmark.ipynb

Provides commands to execute the Jupyter notebook and monitor its log file. Adjust the timeout and ensure the notebook path is correct.

```bash
cd /path/to/notebooks

# Terminal 1: run the notebook
jupyter nbconvert --execute --to notebook \
    --output stochastic_ld_benchmark_executed.ipynb \
    --ExecutePreprocessor.timeout=7200 \
    stochastic_ld_benchmark.ipynb

# Terminal 2: monitor progress
tail -f stochastic_ld_benchmark.log
```

--------------------------------

### susie_init_coef

Source: https://context7.com/stephenslab/susier/llms.txt

Creates a `susie` initialization object from known or estimated non-zero coefficient locations and values. Useful for warm-starting model fitting.

```APIDOC
## susie_init_coef — Initialize SuSiE from Known Coefficients

Creates a `susie` initialization object from known or estimated non-zero coefficient locations and values. Useful for warm-starting model fitting from a previous result or from an external method (e.g., LASSO).

### Method Signature
```r
susie_init_coef(coef_index, coef_value, p)
```

### Parameters
- **coef_index**: A vector of indices where coefficients are non-zero.
- **coef_value**: A vector of the corresponding non-zero coefficient values.
- **p**: The total number of variables.

### Returns
A `susie_init` object that can be passed to the `model_init` argument of `susie()` or `susie_ss()`.

### Example
```r
library(susieR)
set.seed(1)
n <- 1000; p <- 1000
beta <- rep(0, p)
beta[sample(1:p, 4)] <- 1
X <- matrix(rnorm(n * p), nrow = n, ncol = p)
X <- scale(X, center = TRUE, scale = TRUE)
y <- drop(X %*% beta + rnorm(n))

# Initialize susie at true coefficient locations
true_idx   <- which(beta != 0)
true_vals  <- beta[true_idx]
s_init <- susie_init_coef(
  coef_index = true_idx,
  coef_value = true_vals,
  p          = p
)

# Warm-start susie from ground truth
fit_warmstart <- susie(X, y, L = 10, model_init = s_init)

# Compare to cold-start
fit_cold <- susie(X, y, L = 10)
cat("Warm-start iterations:", fit_warmstart$niter, "\n")
catalogue("Cold-start iterations:", fit_cold$niter, "\n")
```
```

--------------------------------

### Fit SuSiE RSS with summary statistics

Source: https://github.com/stephenslab/susier/blob/master/inst/code/python_example/N3finemapping_python.ipynb

Performs univariate regression to obtain summary statistics (bhat, shat) and then fits the SuSiE-RSS model using these statistics, LD matrix (R), sample size (n), and variance of Y. Residual variance is estimated.

```python
conversion_rules = numpy2ri.converter + ro.default_converter
with (conversion_rules).context():   
    sumstats = susie.univariate_regression(
        N3[N3_names['X']],
        np.array(N3[N3_names['Y']])[:,0]
    )
    z_scores = np.array(sumstats[0])/np.array(sumstats[1])
    R = np.corrcoef(N3[N3_names['X']], rowvar=False)
    fitted_rss1 = susie.susie_rss(
        bhat=sumstats[0], 
        shat=sumstats[1], 
        R=R,
        n=n,
        var_y=np.var(np.array(N3[N3_names['Y']])[:,0]),
        L=10,
        estimate_residual_variance=True
    )
    print(ro.r.summary(fitted_rss1)[1])
```

--------------------------------

### Fit SuSiE with Summary Statistics (RSS)

Source: https://context7.com/stephenslab/susier/llms.txt

Use `susie_rss` to fit SuSiE models when only summary statistics (betahat, sebetahat) and LD information (R) are available. Can estimate residual variance and handle different LD inputs.

```r
sumstats <- univariate_regression(X, Y[, 1])
z_scores <- sumstats$betahat / sumstats$sebetahat
R <- cor(X)
fit_rss1 <- susie_rss(
  bhat    = sumstats$betahat,
  shat    = sumstats$sebetahat,
  n       = n,
  R       = R,
  var_y   = var(Y[, 1]),
  L       = 10,
  estimate_residual_variance = TRUE
)
summary(fit_rss1)$cs  # 3 credible sets matching true causal variants
```

```r
fit_rss2 <- susie_rss(z = z_scores, R = R, n = n, L = 10,
                      estimate_residual_variance = TRUE)
all.equal(fit_rss1$pip, fit_rss2$pip)  # TRUE
```

```r
set.seed(42)
X_ref <- matrix(rnorm(500 * ncol(X)), 500, ncol(X))
R_ref <- cor(X_ref)
fit_rss3 <- susie_rss(z_scores, R_ref, n = n, L = 10)
susie_plot(fit_rss3, y = "PIP", b = true_coef[, 1])
```

```r
fit_rss_finite <- susie_rss(
  z      = z_scores,
  R      = R_ref,
  n      = n,
  L      = 10,
  R_finite = 500  # reference panel size
)
# Diagnostics available in fit_rss_finite$R_finite_diagnostics
```

```r
fit_rss_mm <- susie_rss(
  z         = z_scores,
  R         = R_ref,
  n         = n,
  L         = 10,
  R_mismatch = "eb"
)
```

--------------------------------

### Compute config signature and initialize checkpoint

Source: https://github.com/stephenslab/susier/blob/master/inst/notebooks/small_sample_benchmark.ipynb

Computes the MD5 hash of the simulation configuration and initializes the checkpoint directory. This is used to ensure reproducibility and manage simulation states.

```R
config_sig <- compute_config_md5(X_full, y_full, N_vals, h2_sparse,
                                 L_causal, L, noise_scale_factor)
init <- checkpoint_init(outdir, config_sig)
```

--------------------------------

### Update Package Documentation

Source: https://github.com/stephenslab/susier/blob/master/README.md

Execute this command after making changes to roxygen2 markup to update the package NAMESPACE and documentation files.

```bash
make document
```

--------------------------------

### Order and Display Simulation Summary

Source: https://github.com/stephenslab/susier/blob/master/inst/notebooks/stochastic_ld_benchmark.ipynb

Orders the simulation summary data frame by method and prints key performance metrics. This is useful for reviewing the accumulated results over all replicates.

```R
# ============================================================ 
# Load and display summary (from accumulated raw counts)
# ============================================================ 
# To reload from saved file:
# dat <- readRDS(rds_file)
# results <- dat$results
# summary_df <- dat$summary

# Order methods logically
method_order <- c("insample_R")
for (B in B_values) {
  method_order <- c(method_order,
                    paste0("subsample_B", B),
                    paste0("stoch_B", B, "_nocorr"),
                    paste0("stoch_B", B, "_corr"),
                    paste0("stoch_B", B, "_NIG_corr"))
}
summary_df$method <- factor(summary_df$method, levels = method_order)
summary_df <- summary_df[order(summary_df$method), ]

cat("\n=== Summary (accumulated over all replicates) ===\n")
print(summary_df[, c("method", "FDR", "power", "total_cs", "total_false_cs",
                      "total_found", "total_causal", "mean_size")], row.names = FALSE)
```

--------------------------------

### Load and prepare N3finemapping data

Source: https://github.com/stephenslab/susier/blob/master/inst/code/python_example/N3finemapping_python.ipynb

Loads the N3finemapping dataset from the susieR package and sets the random seed. It also extracts dimensions and names for further processing.

```python
conversion_rules = numpy2ri.converter + ro.default_converter
with (conversion_rules).context():   
    ro.r['set.seed'](1)
    N3 = data(susie).fetch('N3finemapping')['N3finemapping']
    N3_names = {v:i for i,v in enumerate(N3.names)}
    n = N3[N3_names['X']].dim[0]
```

--------------------------------

### Run Simulation Replicates and Save Results

Source: https://github.com/stephenslab/susier/blob/master/inst/notebooks/stochastic_ld_benchmark.ipynb

Runs simulation replicates using the generated genotypes and LD approximations. Results are then saved to an RDS file and a summary is printed.

```R
# ============================================================ 
# Phase 4: Run simulation replicates
# ============================================================ 
n_methods <- 1 + length(B_values) * 4  # insample + 4 per B
log_msg("Phase 4: Running %d replicates x %d methods...", n_rep, n_methods)
sim_output <- run_simulation(X, R_true, sketches, R_subs, ld_struct,
                             B_values, n_rep, n_causal, pve_per, L,
                             seed = seed + 100)
results <- sim_output$results
summary_df <- sim_output$summary
log_msg("Simulation complete. %d result rows.", nrow(results))

# ============================================================ 
# Save results to RDS (both per-replicate and summary)
# ============================================================ 
rds_file <- sprintf("stochastic_ld_benchmark_n%d_p%d_nrep%d.rds", n, p, n_rep)
saveRDS(list(results = results, summary = summary_df), file = rds_file)
log_msg("Results saved to %s", rds_file)
close(log_con)  # close log file

cat("\n=== FINAL SUMMARY ===\n")
print(summary_df[, c("method", "FDR", "power", "total_cs", "total_false_cs", "mean_size")],
      row.names = FALSE)
```

--------------------------------

### checkpoint_load

Source: https://github.com/stephenslab/susier/blob/master/inst/notebooks/small_sample_benchmark.ipynb

Loads checkpoint data for a specific simulation setting tag. Returns NULL if the data is absent or corrupt.

```APIDOC
## checkpoint_load()

### Description
Loads checkpoint data for a given setting tag. Returns NULL if the checkpoint file does not exist or is found to be corrupt.

### Parameters
- **outdir** (string) - The directory path for storing checkpoints.
- **tag** (string) - The tag identifying the specific simulation setting.

### Returns
A data.frame containing the checkpoint data if found and valid, otherwise NULL.
```

--------------------------------

### Architecture Diagram

Source: https://github.com/stephenslab/susier/blob/master/inst/misc/README_susie_v2.md

Illustrates the flow of data and control within the susieR 2.0 architecture, from user interface to backend methods.

```text
Interface → Constructor → Workhorse → IBSS Core → Backend Methods
               ↓            ↓
          (data, params) → model
```

--------------------------------

### Run One Simulation Replication

Source: https://github.com/stephenslab/susier/blob/master/inst/notebooks/small_sample_benchmark.ipynb

Runs a single replication of the simulation study, generating data, applying SuSiE, and extracting results. Handles both Gaussian and Non-Gaussian noise models.

```R
run_one_rep <- function(rep_i, n_causal, h2_sp, X, y_real,
                        U_pre = NULL, resid_pre = NULL, L = 10,
                        noise_scale = noise_scale_factor) {
  suppressWarnings({
    n <- nrow(X); p <- ncol(X)

    if (!is.null(U_pre) && !is.null(resid_pre)) {
      U     <- U_pre
      resid <- resid_pre
    } else {
      seed_sub <- make_seed(rep_i, n, "subsample")
      resid    <- get_lasso_residual(X, y_real, seed = seed_sub)
      Xs       <- scale(X, center = TRUE, scale = FALSE)
      svd_X    <- svd(Xs, nu = n, nv = 0)
      U        <- svd_X$u
    }

    # Sign-flip noise, scaled to match real data's noise level
    target_noise_sd <- sd(y_real) * noise_scale
    seed_flip <- make_seed(rep_i, n, "flip", h2 = h2_sp, nc = n_causal)
    noise <- gen_signflip_noise(U, resid, target_noise_sd, seed = seed_flip)

    # Causal signal: h2 = var(signal) / var(y)
    seed_causal <- make_seed(rep_i, n, "causal", h2 = h2_sp, nc = n_causal)
    set.seed(seed_causal)
    causal_idx <- sample(p, n_causal)
    y <- noise
    for (j in causal_idx) {
      bj <- sqrt(h2_sp * var(noise) / ((1 - h2_sp) * n_causal * var(X[, j])))
      y  <- y + X[, j] * bj
    }

    # Fit both methods
    fit_gaus <- tryCatch(
      susie(X, y, L = L, verbose = FALSE),
      error = function(e) NULL)
    fit_ss <- tryCatch(
      susie(X, y, L = L, estimate_residual_method = "NIG",
            verbose = FALSE),
      error = function(e) NULL)

    extract <- function(fit, tag) {
      na_row <- data.frame(
        method = tag, rep = rep_i,
        discovered = NA_real_, n_true_cs = NA_real_, n_cs = NA_real_,
        mean_size = NA_real_, sigma2 = NA_real_,
        mean_V = NA_real_, max_V = NA_real_, sum_V = NA_real_,
        stringsAsFactors = FALSE)
      if (is.null(fit)) return(na_row)

      cs_obj <- susie_get_cs(fit, X = X, min_abs_corr = 0.5)
      cs     <- cs_obj$cs
      ncs    <- length(cs)

      discovered <- 0; n_true_cs <- 0; avg_size <- NA_real_
      if (ncs > 0) {
        discovered <- length(intersect(unique(unlist(cs)), causal_idx))
        n_true_cs  <- sum(sapply(cs, function(s) any(causal_idx %in% s)))
        avg_size   <- mean(sapply(cs, length))
      }

      V_vec <- fit$V
      if (is.null(V_vec)) V_vec <- rep(NA_real_, L)
      if (length(V_vec) == 1) V_vec <- rep(V_vec, L)

      data.frame(
        method     = tag,
        rep        = rep_i,
        discovered = discovered,
        n_true_cs  = n_true_cs,
        n_cs       = ncs,
        mean_size  = avg_size,
        sigma2     = fit$sigma2,
        mean_V     = mean(V_vec, na.rm = TRUE),
        max_V      = max(V_vec, na.rm = TRUE),
        sum_V      = sum(V_vec, na.rm = TRUE),
        stringsAsFactors = FALSE)
    }

    rbind(extract(fit_gaus, "Gaussian"), extract(fit_ss, "SS"))
  })
}
```

--------------------------------

### Main Simulation Loop

Source: https://github.com/stephenslab/susier/blob/master/inst/notebooks/stochastic_ld_benchmark.ipynb

The main function orchestrates the simulation process. For each replicate, it simulates a phenotype, computes summary statistics (z-scores), runs all specified fine-mapping methods, and collects the results. Finally, it prints a summary of the aggregated results.

```R
# Main loop: for each replicate, simulate phenotype, compute z,
# run all methods, collect results. Print final summary.

```

--------------------------------

### checkpoint_save_meta

Source: https://github.com/stephenslab/susier/blob/master/inst/notebooks/small_sample_benchmark.ipynb

Writes the checkpoint metadata, including configuration, MD5 hash, number of repetitions, and timestamp, after successful completion.

```APIDOC
## checkpoint_save_meta()

### Description
Writes the checkpoint metadata to the `_meta.rds` file after a simulation run has successfully completed. This includes the configuration details, MD5 hash, the total number of repetitions, and the completion timestamp.

### Parameters
- **outdir** (string) - The directory path for storing checkpoints.
- **config_sig** (list) - A list containing the simulation configuration and its MD5 hash.
- **n_rep** (integer) - The total number of simulation repetitions completed.

### Returns
This function saves a metadata file and does not return any value.
```

--------------------------------

### Generate Genotypes and Compute LD Approximations

Source: https://github.com/stephenslab/susier/blob/master/inst/notebooks/stochastic_ld_benchmark.ipynb

Generates genotypes and then computes LD approximations. This involves pre-computing sketches and subset correlations for different block sizes (B_values).

```R
log_msg("Phase 2: Generating genotypes (%d x %d)...", n, p)
t0 <- proc.time()
geno <- generate_genotypes(n, ld_struct, seed = seed + 1)
X      <- geno$X
R_true <- geno$R_true
rm(geno)
log_msg("Phase 2 complete in %.1f min", (proc.time() - t0)[3] / 60)

# ============================================================ 
# Phase 3: Pre-compute LD approximations
# ============================================================ 
log_msg("Phase 3: Computing LD approximations for B = {%s}...",
        paste(B_values, collapse = ", "))
t0 <- proc.time()
ld_approx <- compute_ld_approximations(X, B_values, seed = seed + 2)
sketches <- ld_approx$sketches
R_subs   <- ld_approx$R_subs
rm(ld_approx)
log_msg("Phase 3 complete in %.1f min", (proc.time() - t0)[3] / 60)
```

--------------------------------

### Set up plotting theme and colors

Source: https://github.com/stephenslab/susier/blob/master/inst/notebooks/small_sample_benchmark.ipynb

Defines a custom theme for plots using cowplot and ggplot2, including color palettes for 'SS' and 'Gaussian' methods. This theme is applied to ensure consistent visualization of benchmark results.

```R
library(ggplot2)
library(cowplot)
library(gridExtra)
library(grid)

figdir <- "benchmark_results"

methods_colors <- c("SS" = "#D41159", "Gaussian" = "#1A85FF")

perf_theme <- theme_cowplot(font_size = 16) +
  theme(
    legend.position  = "none",
    panel.grid.major.y = element_line(color = "gray80"),
    panel.grid.major.x = element_blank(),
    panel.grid.minor   = element_blank(),
    axis.line   = element_line(linewidth = 1, color = "black"),
    axis.ticks  = element_line(linewidth = 1, color = "black"),
    axis.ticks.length = unit(0.25, "cm"),
    plot.margin = margin(t = 2, r = 2, b = 2, l = 2, unit = "mm"),
    axis.text   = element_text(size = 14, face = "bold"),
    axis.title  = element_text(size = 16, face = "bold"),
    plot.title  = element_text(size = 16, face = "bold")
  )
dot_size <- 4
cat("Plot theme set.\n")
```

--------------------------------

### Compute MD5 Hash of Simulation Parameters

Source: https://github.com/stephenslab/susier/blob/master/inst/notebooks/small_sample_benchmark.ipynb

Calculates an MD5 hash for simulation parameters to detect changes. This is useful for automatically invalidating previous results if configuration changes.

```R
compute_config_md5 <- function(X, y, N_vals, h2_sparse, L_causal,
                               L, noise_scale) {
  nr <- min(50, nrow(X)); nc_dim <- min(50, ncol(X))
  config <- list(
    data_nrow   = nrow(X),
    data_ncol   = ncol(X),
    data_corner = sum(X[1:nr, 1:nc_dim]),
    var_y       = round(var(as.vector(y)), 8),
    noise_scale = round(noise_scale, 8),
    L           = as.integer(L),
    N_vals      = sort(as.integer(N_vals)),
    h2_sparse   = sort(round(h2_sparse, 6)),
    L_causal    = sort(as.integer(L_causal))
  )
  md5 <- digest(config, algo = "md5")
  list(config = config, md5 = md5)
}
```

--------------------------------

### Print R session information

Source: https://github.com/stephenslab/susier/blob/master/inst/code/python_example/N3finemapping_python.ipynb

Prints the session information for the R environment, including R version, platform, loaded packages, and locale settings. This is useful for reproducibility.

```python
print(ro.r.sessionInfo())
```

--------------------------------

### Fit SuSiE with Regularized RSS Likelihood

Source: https://context7.com/stephenslab/susier/llms.txt

Use `susie_rss_lambda` for numerically unstable LD matrices. It requires a regularization parameter `lambda` and can estimate it or use a provided factor matrix.

```r
library(susieR)
data(N3finemapping)
attach(N3finemapping)

sumstats <- univariate_regression(X, Y[, 1])
z_scores <- sumstats$betahat / sumstats$sebetahat
R        <- cor(X)

# Fit with fixed lambda = 0 (standard path, no regularization)
fit_lambda0 <- susie_rss_lambda(
  z      = z_scores,
  R      = R,
  n      = nrow(X),
  L      = 10,
  lambda = 0
)
```

```r
# Estimate lambda from null-space residual
fit_lambda_est <- susie_rss_lambda(
  z      = z_scores,
  R      = R,
  n      = nrow(X),
  L      = 10,
  lambda = "estimate"
)
cat("Estimated lambda:", fit_lambda_est$lambda, "\n")
```

```r
# Provide factor matrix X_ref instead of R (avoids forming p x p matrix)
set.seed(1)
X_ref <- matrix(rnorm(300 * ncol(X)), 300, ncol(X))
fit_factor <- susie_rss_lambda(
  z      = z_scores,
  X      = X_ref,
  n      = nrow(X),
  L      = 10,
  lambda = 1e-3,
  prior_variance = 50
)
```

--------------------------------

### Incremental simulation loop

Source: https://github.com/stephenslab/susier/blob/master/inst/notebooks/small_sample_benchmark.ipynb

Iterates through different sample sizes (N) and simulation settings (h2, L_causal). It checks for completed simulations, precomputes necessary components for new simulations, and runs or extends existing ones, saving results incrementally.

```R
all_results <- list()
t_total     <- proc.time()
n_full      <- nrow(X_full)
ns          <- noise_scale_factor
any_new     <- FALSE

for (N in N_vals) {
  # ── Which rep_ids does ANY setting at this N still need? ──
  # We precompute LASSO+SVD only for these, saving time when extending.
  needed_reps <- integer(0)
  for (h2 in h2_sparse) {
    for (nc in L_causal) {
      tag  <- sprintf("N%d_h2%03d_nc%d", N, round(h2 * 100), nc)
      done <- checkpoint_completed_reps(outdir, tag)
      todo <- setdiff(seq_len(n_rep), done)
      needed_reps <- union(needed_reps, todo)
    }
  }
  needed_reps <- sort(needed_reps)

  if (length(needed_reps) == 0) {
    # Everything cached — just load
    cat(sprintf("\n=== N = %d: all settings complete, loading ===\n", N))
    for (h2 in h2_sparse) {
      for (nc in L_causal) {
        tag <- sprintf("N%d_h2%03d_nc%d", N, round(h2 * 100), nc)
        res <- checkpoint_load(outdir, tag)
        # Trim to exactly n_rep reps
        done_ids <- sort(unique(res$rep))
        if (length(done_ids) > n_rep)
          res <- res[res$rep %in% done_ids[1:n_rep], ]
        all_results[[length(all_results) + 1]] <- res
      }
    }
    next
  }

  # ── Precompute LASSO + SVD only for needed reps ──
  cat(sprintf("\n=== N = %d: precomputing %d / %d subsamples ===\n",
      N, length(needed_reps), n_rep))
  t_pre <- proc.time()

  keep_list  <- vector("list", n_rep)
  U_list     <- vector("list", n_rep)
  resid_list <- vector("list", n_rep)

  for (i in needed_reps) {
    seed_i <- make_seed(i, N, "subsample")
    set.seed(seed_i)
    keep <- sample(n_full, N)
    keep_list[[i]] <- keep
    Xi <- X_full[keep, ]
    yi <- y_full[keep]
    resid_list[[i]] <- get_lasso_residual(Xi, yi, seed = seed_i)
    Xs <- scale(Xi, center = TRUE, scale = FALSE)
    svd_i <- svd(Xs, nu = N, nv = 0)
    U_list[[i]] <- svd_i$u
  }
  cat(sprintf("  Precompute: %.0f sec\n", (proc.time() - t_pre)[3]))

  # ── Run or extend each setting ──
  for (h2 in h2_sparse) {
    for (nc in L_causal) {
      tag       <- sprintf("N%d_h2%03d_nc%d", N, round(h2 * 100), nc)
      done_reps <- checkpoint_completed_reps(outdir, tag)
      todo_reps <- sort(setdiff(seq_len(n_rep), done_reps))

      if (length(todo_reps) == 0) {
        # Already complete — load
        cat(sprintf("[DONE] %s (%d reps)\n", tag, length(done_reps)))
        res <- checkpoint_load(outdir, tag)
        done_ids <- sort(unique(res$rep))
        if (length(done_ids) > n_rep)
          res <- res[res$rep %in% done_ids[1:n_rep], ]
        all_results[[length(all_results) + 1]] <- res
        next
      }

      if (length(done_reps) > 0) {
        cat(sprintf("[EXT]  %s: %d -> %d reps ... ",
            tag, length(done_reps), length(done_reps) + length(todo_reps)))
      } else {
        cat(sprintf("[RUN]  %s (%d reps) ... ", tag, length(todo_reps)))
      }
      t0 <- proc.time()
      any_new <- TRUE

      new_res <- do.call(rbind, future_lapply(todo_reps, function(i) {
        run_one_rep(i, nc, h2,
                    X_full[keep_list[[i]], ], y_full[keep_list[[i]]],
                    U_pre = U_list[[i]], resid_pre = resid_list[[i]],
                    L = L, noise_scale = ns)
      }, future.seed = TRUE))
      new_res$N         <- N
      new_res$h2_sparse <- h2
      new_res$n_causal  <- nc

      # Save with dedup, then collect
      res <- checkpoint_save(outdir, tag, new_res)
      all_results[[length(all_results) + 1]] <- res
      cat(sprintf("done (%.0f sec)\n", (proc.time() - t0)[3]))
    }
  }
}

# --- Finalize ---
results <- do.call(rbind, all_results)
saveRDS(results, all_rds)
checkpoint_save_meta(outdir, config_sig, n_rep)

total_settings <- length(N_vals) * length(h2_sparse) * length(L_causal)
if (any_new) {
  cat(sprintf("\nCompleted: %d rows across %d settings, %.1f min\n",
      nrow(results), total_settings, (proc.time() - t_total)[3] / 60))
} else {
  cat(sprintf("\nAll %d settings x %d reps loaded from cache (%d rows)\n",
      total_settings, n_rep, nrow(results)))
}
cat(sprintf("Config MD5: %s\nSaved: %s\n", config_sig$md5, all_rds))
```

--------------------------------

### Execute Notebook

Source: https://github.com/stephenslab/susier/blob/master/inst/notebooks/small_sample_benchmark.ipynb

Command to execute the Jupyter notebook, disabling cell timeouts for long-running simulations. Use this to run the benchmark.

```bash
jupyter nbconvert --to notebook --execute \
  --ExecutePreprocessor.timeout=0 \
  --output small_sample_benchmark_executed.ipynb \
  small_sample_benchmark.ipynb
```

--------------------------------

### Draw Posterior Samples with susie_get_posterior_samples

Source: https://context7.com/stephenslab/susier/llms.txt

Draws Monte Carlo samples from the posterior distribution of regression coefficients. Useful for propagating uncertainty into downstream analyses. Requires the 'susieR' library.

```r
library(susieR)
set.seed(1)
n <- 500; p <- 300
beta <- rep(0, p); beta[c(20, 150)] <- c(1, -1)
X <- matrix(rnorm(n * p), n, p)
y <- X %*% beta + rnorm(n)
fit <- susie(X, y, L = 10)

# Draw 1000 posterior samples
samples <- susie_get_posterior_samples(fit, num_samples = 1000)
# samples$b:     300 x 1000 matrix of effect sizes
# samples$gamma: 300 x 1000 matrix of causal inclusion indicators (0/1)

# Posterior mean of b from samples (approximates fit$pip * posterior_mean)
b_mean_sample <- rowMeans(samples$b)

# Inclusion probability from samples
pip_sample <- rowMeans(samples$gamma)

# Compare to analytic PIPs
plot(fit$pip, pip_sample,
     xlab = "Analytic PIP", ylab = "Monte Carlo PIP")
abline(0, 1, col = "blue")

# Credible intervals from samples (e.g., 95% CI for variable 20)
ci_20 <- quantile(samples$b[20, ], c(0.025, 0.975))
cat("95% CI for variable 20:", ci_20, "\n")
```

--------------------------------

### susie_ss() — SuSiE with Sufficient Statistics

Source: https://context7.com/stephenslab/susier/llms.txt

Fits the SuSiE model using precomputed sufficient statistics (XtX, Xty, yty, n). This is efficient when the number of samples (n) is much larger than the number of variables (p), or when fitting the model for multiple response vectors with the same design matrix.

```APIDOC
## susie_ss(XtX, Xty, yty, n, X_colmeans, y_mean, L, ...)

### Description
Fits SuSiE from precomputed sufficient statistics `(XtX, Xty, yty, n)` instead of raw data. This is preferred when `n >> p` and multiple response vectors `y` share the same design matrix `X` — `XtX` is computed once and reused. Results are numerically equivalent to `susie()` when inputs are derived from the same data.

### Parameters
- **XtX** (matrix) - The matrix product X'X.
- **Xty** (vector) - The matrix-vector product X'y.
- **yty** (scalar) - The scalar product y'y.
- **n** (integer) - The number of samples.
- **X_colmeans** (vector) - The column means of X.
- **y_mean** (scalar) - The mean of y.
- **L** (integer) - Maximum number of single effects to assume.
- **...** - Additional arguments passed to the fitting function.

### Request Example
```r
library(susieR)
data(N3finemapping)
attach(N3finemapping)

# Compute sufficient statistics from (X, y)
ss <- compute_suff_stat(X, Y[, 1])

# Fit using sufficient statistics
fit_ss <- susie_ss(
  XtX        = ss$XtX,
  Xty        = ss$Xty,
  yty        = ss$yty,
  n          = ss$n,
  X_colmeans = ss$X_colmeans,
  y_mean     = ss$y_mean,
  L          = 10,
  estimate_residual_variance = TRUE
)
```

### Response
#### Success Response
A `susie` fit object containing:
- **pip** (vector) - Per-variable posterior inclusion probabilities.
- **sets** (list) - Credible sets.

### Response Example
```r
# Verify equivalence with susie()
fit_raw <- susie(X, Y[, 1], L = 10)
all.equal(fit_raw$pip, fit_ss$pip)  # Should return TRUE
```
```

--------------------------------

### Fit SuSiE with raw data

Source: https://github.com/stephenslab/susier/blob/master/inst/code/python_example/N3finemapping_python.ipynb

Fits the standard SuSiE model directly using the raw genotype (X) and phenotype (Y) data. The number of credible sets (L) is set to 10.

```python
conversion_rules = numpy2ri.converter + ro.default_converter
with (conversion_rules).context():   
    fitted = susie.susie(N3[N3_names['X']], 
                         np.array(N3[N3_names['Y']])[:,0], 
                         L=10
                        )
    print(ro.r['all.equal'](fitted[-2], fitted_rss1[-1]))
```