diff --git a/src/bulk.rs b/src/bulk.rs index 57ff1734..8b2f1442 100644 --- a/src/bulk.rs +++ b/src/bulk.rs @@ -174,7 +174,7 @@ fn perform_inference_and_write_output( // if the user requested bootstrap replicates, // compute and write those out now. if args.num_bootstraps > 0 { - let breps = em::bootstrap(&emi, args.num_bootstraps, args.threads); + let breps = em::bootstrap(&emi, args.num_bootstraps, args.threads, args.seed); let mut new_arrays = vec![]; let mut bs_fields = vec![]; diff --git a/src/em.rs b/src/em.rs index 42cb022b..a4e97afe 100644 --- a/src/em.rs +++ b/src/em.rs @@ -5,7 +5,6 @@ use crate::util::oarfish_types::{AlnInfo, EMInfo, TranscriptInfo}; use atomic_float::AtomicF64; use itertools::izip; use num_format::{Locale, ToFormattedString}; -use rand::rng as trng; use rayon::iter::{IntoParallelIterator, IntoParallelRefIterator, ParallelIterator}; use tracing::{info, span, trace}; @@ -270,8 +269,20 @@ pub fn em(em_info: &EMInfo, _nthreads: usize) -> Vec { do_em(em_info, make_iter, true) } -pub fn do_bootstrap(em_info: &EMInfo) -> Vec { - let mut rng = trng(); +pub fn do_bootstrap(em_info: &EMInfo, seed: Option, rep_idx: u64) -> Vec { + use rand::SeedableRng; + use rand::rngs::StdRng; + + let mut rng: StdRng = match seed { + Some(s) => { + // Derive a deterministic per-replicate seed so that parallel execution + // yields the same results regardless of scheduling. + let derived = s.wrapping_add(0x9E37_79B9_7F4A_7C15u64.wrapping_mul(rep_idx + 1)); + StdRng::seed_from_u64(derived) + } + None => StdRng::from_os_rng(), + }; + let n = em_info.eq_map.len(); let inds = bootstrap::get_sample_inds(n, &mut rng); @@ -289,7 +300,12 @@ pub fn do_bootstrap(em_info: &EMInfo) -> Vec { do_em(em_info, make_iter, false) } -pub fn bootstrap(em_info: &EMInfo, num_boot: u32, nthreads: usize) -> Vec> { +pub fn bootstrap( + em_info: &EMInfo, + num_boot: u32, + nthreads: usize, + seed: Option, +) -> Vec> { let span = span!(tracing::Level::INFO, "bootstrap"); let _guard = span.enter(); @@ -307,7 +323,7 @@ pub fn bootstrap(em_info: &EMInfo, num_boot: u32, nthreads: usize) -> Vec anyhow::Result<()> { let mut args = Args::parse(); + // If a seed is provided force single-threaded execution + if let Some(seed) = args.seed { + info!("reproducibility: using seed {seed}; forcing single-threaded execution"); + args.threads = 1; + } + // change the logging filter if the user specified quiet or // verbose. if args.quiet { diff --git a/src/prog_opts.rs b/src/prog_opts.rs index e4ec3c53..4fbfcfa6 100644 --- a/src/prog_opts.rs +++ b/src/prog_opts.rs @@ -418,4 +418,9 @@ pub struct Args { /// use a KDE model of the observed fragment length distribution #[arg(short, long, hide = true)] pub use_kde: bool, + + /// RNG seed to make stochastic components (e.g. bootstrapping) deterministic. + /// If provided, oarfish will also force the EM implementation to be single-threaded + #[arg(long, help_heading = "reproducibility")] + pub seed: Option, }