From 08e7c9332732f9babf1dea23cdafb6516517c038 Mon Sep 17 00:00:00 2001
From: Daniel Bauer <bauer@t480s>
Date: Sun, 8 Mar 2020 12:29:52 +0100
Subject: [PATCH] seedable rng

---
 src/cli.yml           |  5 +++++
 src/error_analysis.rs | 15 +++++++++------
 src/lib.rs            |  7 +++++--
 src/main.rs           | 10 +++++++++-
 4 files changed, 28 insertions(+), 9 deletions(-)

diff --git a/src/cli.yml b/src/cli.yml
index 35add1e..545e738 100644
--- a/src/cli.yml
+++ b/src/cli.yml
@@ -86,6 +86,11 @@ args:
         help: Number of bayesian bootstrapping runs for error analysis by assigning random weights (defaults to 0).
         takes_value: true
         required: false
+    - bootstrap_seed:
+        long: seed
+        help: Random seed for bootstrapping runs.
+        takes_value: true
+        required: false
     - start:
         long: start
         help: Skip rows in timeseries with an index smaller than this value (defaults to 0)
diff --git a/src/error_analysis.rs b/src/error_analysis.rs
index 32cf472..5b7d3af 100644
--- a/src/error_analysis.rs
+++ b/src/error_analysis.rs
@@ -1,4 +1,4 @@
-use rand;
+use rand::{SeedableRng, StdRng, Rng};
 use super::histogram::{Dataset};
 use super::perform_wham;
 use super::Config;
@@ -7,9 +7,9 @@ use rgsl::statistics;
 // returns a set of num_windows continious weights by
 // a) generate num_windows-1 random variables and sort them
 // b) each weight n is the difference between n+1 and n, where n0=0 and nN+1=1
-fn generate_random_weights(num_windows: usize) -> Vec<f64> {
+fn generate_random_weights(num_windows: usize, rng: &mut StdRng) -> Vec<f64> {
     // create a list of num_windows - 1 sorted random numbers and append/prepend 0 and 1
-    let mut tmp = (0..num_windows-1).map(|_| rand::random::<f64>()).collect::<Vec<f64>>();
+    let mut tmp = (0..num_windows-1).map(|_| rng.gen::<f64>()).collect::<Vec<f64>>();
     tmp.sort_by(|a,b| { a.partial_cmp(b).unwrap() });
     let mut rnds = vec![0.0];
     rnds.append(&mut tmp);
@@ -24,8 +24,8 @@ fn generate_random_weights(num_windows: usize) -> Vec<f64> {
 }
 
 // Generate a random weighted dataset from the given dataset by changing the weights
-fn generate_random_weighted_dataset(ds: Dataset) -> Dataset {
-    let weights = generate_random_weights(ds.num_windows);
+fn generate_random_weighted_dataset(ds: Dataset, rng: &mut StdRng) -> Dataset {
+    let weights = generate_random_weights(ds.num_windows, rng);
     Dataset::new_weighted(ds, weights)
 }
 
@@ -33,10 +33,13 @@ fn generate_random_weighted_dataset(ds: Dataset) -> Dataset {
 // datasets. The standard deviation is calculated on the bootstrapped probabilities of each bin. The
 // standard deviation of the free eneergy is then deduced by error propagation (A_std = kT*1/P*P_std)
 pub fn run_bootstrap(cfg: &Config, ds: Dataset, P: &[f64], num_runs: usize) -> (Vec<f64>,Vec<f64>) {
+    // seed the rng
+    let mut rng: StdRng = SeedableRng::seed_from_u64(cfg.bootstrap_seed);
+
     // Calculate bootstrapped probabilities
     let bootstrapped_Ps: Vec<Vec<f64>> = (0..num_runs).map(|x| {
         println!("Bootstrap run {}/{}", x, num_runs);
-        let rnd_weighted_dataset = generate_random_weighted_dataset(ds.clone());
+        let rnd_weighted_dataset = generate_random_weighted_dataset(ds.clone(), &mut rng);
         perform_wham(cfg, &rnd_weighted_dataset).unwrap().0
     }).collect();
 
diff --git a/src/lib.rs b/src/lib.rs
index ffe8bd5..4d5bf2e 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -38,14 +38,17 @@ pub struct Config {
 	pub cyclic: bool,
 	pub output: String,
 	pub bootstrap: usize,
+    pub bootstrap_seed: u64,
     pub start: f64,
     pub end: f64,
 }
 
 impl fmt::Display for Config {
 	 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
-         write!(f, "Metadata={}, hist_min={:?}, hist_max={:?}, bins={:?} verbose={}, tolerance={}, iterations={}, temperature={}, cyclic={:?}, bootstrap={:?}", self.metadata_file, self.hist_min, self.hist_max, self.num_bins,
-                self.verbose, self.tolerance, self.max_iterations, self.temperature, self.cyclic, self.bootstrap)
+         write!(f, "Metadata={}, hist_min={:?}, hist_max={:?}, bins={:?} verbose={}, tolerance={}, iterations={}, temperature={}, cyclic={:?}, bootstrap={:?}, seed={:?}",
+            self.metadata_file, self.hist_min, self.hist_max, self.num_bins,
+            self.verbose, self.tolerance, self.max_iterations, self.temperature,
+            self.cyclic, self.bootstrap, self.bootstrap_seed)
     }
 }
 
diff --git a/src/main.rs b/src/main.rs
index 02bd99d..913fa05 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -1,6 +1,7 @@
 extern crate wham;
 #[macro_use]
 extern crate clap;
+extern crate rand;
 
 use clap::App;
 use wham::Config;
@@ -46,6 +47,13 @@ fn cli() -> Result<Config> {
         .split(',').map(|x| { x.parse().unwrap() }).collect();
 	let bootstrap: usize = matches.value_of("bootstrap").unwrap_or("0").parse()
 		.chain_err(|| "Cannot parse bootstrap iteration.")?;
+    let bootstrap_seed: u64 = matches.value_of("bootstrap_seed")
+        .unwrap_or({
+            use rand::Rng;
+            let mut rng = rand::thread_rng();
+            &rng.gen::<u32>().to_string()
+        }).parse()
+        .chain_err(|| "Cannot parse bootstrap iteration.")?;
     let start: f64 = matches.value_of("start").unwrap_or("0").parse()
         .chain_err(|| "Cannot parse start time.")?;
     let end: f64 = matches.value_of("end").unwrap_or("1e+20").parse()
@@ -61,7 +69,7 @@ fn cli() -> Result<Config> {
 
 	Ok(wham::Config{metadata_file, hist_min, hist_max, num_bins, dimens,
 		verbose, tolerance, max_iterations, temperature, cyclic, output,
-		bootstrap, start, end})
+		bootstrap, bootstrap_seed, start, end})
 }
 
 fn main() {