From 36e0c31b1e912820974636849843e34a83571536 Mon Sep 17 00:00:00 2001
From: Matthew Mah <mmah@broadinstitute.org>
Date: Tue, 8 Jul 2025 17:21:53 -0400
Subject: [PATCH] option to set RNG seed from Philip Johnson July 2017,
 formerly called version 1.0-12

---
 R/output.R      | 1 +
 exec/estimate.R | 9 +++++++++
 2 files changed, 10 insertions(+)
 mode change 100755 => 100644 exec/estimate.R

diff --git a/R/output.R b/R/output.R
index 1454895..75695cf 100644
--- a/R/output.R
+++ b/R/output.R
@@ -80,6 +80,7 @@ print.contamMix <- function(x, tabDelim = FALSE, ...) {
   if (tabDelim) {
     cat(x$e, anal$MAP, anal$ci,
         if (length(x$chains) > 1) coda::gelman.diag(anal$mcmcObj)$psrf else "",
+        x$seed,
         sep="\t");
     cat("\n")
   } else {
diff --git a/exec/estimate.R b/exec/estimate.R
old mode 100755
new mode 100644
index 05f2c7b..d9b8f93
--- a/exec/estimate.R
+++ b/exec/estimate.R
@@ -18,6 +18,7 @@ opt = getopt(rbind(
   c('saveMN', NA, 0, "logical"),    #save MN intermediate for debugging
   c('recordAll', NA, 0, "logical"), #record *all* proportions from MCMC chain
   c('tabOutput', NA, 0, "logical"), #simple tab-delimited output
+  c('seed', NA, 0, "integer"), #set rng seed
   c('nrThreads', NA, 1, "integer")) # Number of threads to use. Defaults to 1. 
   )
 if (is.null(opt$nIter))  opt$nIter = 50000
@@ -61,6 +62,14 @@ cat(paste0("Using ",used_cores, " threads.\n"), file=stderr())
 
 ## do the work
 data = loadSAM(samFn=opt$samFn, malnFn=opt$malnFn, endogId=opt$consId, baseqThreshold=opt$baseq, transverOnly=opt$transverOnly, trimBases=opt$trimBases, saveMN=opt$saveMN)
+## set/record RNG seed
+if (is.null(opt$seed)) {
+    data$seed = round(runif(1,max=1e5))
+} else {
+    data$seed = opt$seed
+}
+RNGkind("L'Ecuyer-CMRG")
+set.seed(data$seed)
 res = runChains(nChains=opt$nChains, nIter=opt$nIter, data=data, alpha=opt$alpha)
 
 ## results in text format