bootstrapping

src/cli.yml

@@ -61,6 +61,7 @@ args:
         required: false
         help: Stop WHAM after this many iterations without convergence (defaults to 100,000).
     - cyclic:
+        short: c
         long: cyclic
         help: For periodic reaction coordinates. If this is set, the first and last coordinate bin in each dimension are treated as neighbors for the bias calculation.
     - verbose:
@@ -71,12 +72,17 @@ args:
     - temperature:
         short: T
         long: temperature
-        help: WHAM temperature in Kelvin
+        help: WHAM temperature in Kelvin.
         takes_value: true
         required: true
     - output:
         short: o
         long: output
-        help: Free energy output file (defaults to wham.out)
+        help: Free energy output file (defaults to wham.out).
         takes_value: true
-        required: false
+        required: false
+    - bootstrap:
+        long: bt
+        help: Number of bayesian bootstrapping runs for error analysis by assigning random weights (defaults to 0).
+        takes_value: true
+        required: false

src/error_analysis.rs

@@ -0,0 +1,95 @@
+use rand;
+use super::histogram::{Dataset,Histogram};
+use super::k_B;
+use super::perform_wham;
+use super::Config;
+use rgsl::statistics;
+
+fn generate_random_weights(num_windows: usize) -> Vec<f64> {
+    // TODO this is ugly.
+    // 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>>();
+    tmp.sort_by(|a,b| { a.partial_cmp(b).unwrap() });
+    let mut rnds = vec![0.0];
+    rnds.append(&mut tmp);
+    rnds.append(&mut vec![1.0]);
+
+    // weights of window i is the difference between rnd[i+1] and rnd[i]
+    let mut weights = vec![0.0; num_windows];
+    for i in 0..num_windows {
+        weights[i] = rnds[i+1] - rnds[i]
+    }
+    return weights
+}
+
+
+fn generate_random_weighted_dataset(ds: Dataset) -> Dataset {
+    let weights = generate_random_weights(ds.num_windows);
+    Dataset::new_weighted(ds, weights )
+}
+
+pub fn run_bootstrap(cfg: &Config, ds: Dataset, P: &[f64], num_runs: usize) -> (Vec<f64>,Vec<f64>) {
+    let 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());
+        perform_wham(cfg, &rnd_weighted_dataset).unwrap().0
+    }).collect();
+    let mut P_std = vec![0.0; ds.num_bins];
+    for bin in 0..ds.num_bins {
+        let Ps = Ps.iter().map(|window| window[bin]).collect::<Vec<f64>>();
+        P_std[bin] = statistics::sd(&Ps, 1, num_runs);
+    }
+    let A_std = P_std.iter().zip(P.iter()).map(|(std,P)| ds.kT*1.0/P*std).collect();
+    (P_std, A_std)
+}
+
+mod tests {
+    use super::*;
+
+    fn build_hist() -> Histogram {
+		Histogram::new(
+			22, // num_points
+			vec![1.0, 1.0, 3.0, 5.0, 12.0] // bins
+		)
+	}
+
+	fn build_hist_set() -> Dataset {
+		let h1 = build_hist();
+		let h2 = build_hist();
+		let h3 = build_hist();
+		Dataset::new(
+			5, // num bins
+			vec![3],
+			vec![1.0], // bin width
+			vec![0.0], // hist min
+			vec![9.0], // hist max
+			vec![4.5, 4.5, 4.5], // x0
+			vec![10.0, 10.0, 10.0], // fc
+			300.0*k_B, // kT
+			vec![h1, h2, h3], // hists
+			false // cyclic
+		)
+	}
+
+    #[test]
+    fn random_weights() {
+        let num_windows = 5;
+        let weights = generate_random_weights(num_windows);
+        assert_eq!(num_windows, weights.len());
+        for w in weights {
+            assert!(0.0 < w && w < 1.0);
+        }
+    }
+
+    #[test]
+    fn random_weighted_dataset() {
+        let ds = build_hist_set();
+        let rnd_weights_ds = generate_random_weighted_dataset(ds);
+        println!("{:?}", rnd_weights_ds.weights);
+        for w in rnd_weights_ds.weights {
+            assert!(w > 0.0);
+            assert!(w < 1.0);
+        }
+    }
+
+}

src/histogram.rs

@@ -2,7 +2,7 @@ use std::fmt;
 use std::cell::RefCell;
 
 // One histogram
-#[derive(Debug)]
+#[derive(Debug,Clone)]
 pub struct Histogram {
 	// total number of data points stored in the histogram
 	pub num_points: u32,
@@ -18,7 +18,7 @@ impl Histogram {
 }
 
 // a set of histograms
-#[derive(Debug)]
+#[derive(Debug,Clone)]
 pub struct Dataset {
 	// number of histogram windows (number of simulations)
 	pub num_windows: usize,
@@ -56,8 +56,8 @@ pub struct Dataset {
 	// bias value cache
 	bias: RefCell<Vec<Option<f64>>>,
 
-	// sum of bin count for windows
-	pub bin_count: Vec<f64>
+	// histogram weight
+	pub weights: Vec<f64>,
 }
 
 impl Dataset {
@@ -65,9 +65,7 @@ impl Dataset {
 	pub fn new(num_bins: usize, dimens_lengths: Vec<usize>, bin_width: Vec<f64>, hist_min: Vec<f64>, hist_max: Vec<f64>, bias_pos: Vec<f64>, bias_fc: Vec<f64>, kT: f64, histograms: Vec<Histogram>, cyclic: bool) -> Dataset {
 		let num_windows = histograms.len();
 		let bias: RefCell<Vec<Option<f64>>> = RefCell::new(vec![None; num_bins*num_windows]);
-		let bin_count = (0..num_bins).map(|bin| {
-			histograms.iter().map(|h| h.bins[bin]).sum()
-		}).collect();
+		let weights = vec![1.0; num_windows];
 		Dataset{
 			num_windows,
 			num_bins,
@@ -81,10 +79,21 @@ impl Dataset {
 			bias_pos,
 			bias_fc,
 			bias,
-			bin_count,
+			weights
 		}
 	}
 
+	pub fn new_weighted(ds: Dataset, weights: Vec<f64>) -> Dataset {
+		Dataset {
+			weights: weights,
+			..ds
+		}
+	}
+
+	pub fn get_weighted_bin_count(&self, bin: usize) -> f64 {
+		self.histograms.iter().enumerate().map(|(idx,h)| self.weights[idx]*h.bins[bin]).sum()
+	}
+
 	fn expand_index(&self, bin: usize, lengths: &[usize]) -> Vec<usize> {
     	let mut tmp = bin;
     	let mut idx = vec![0; lengths.len()];
@@ -247,10 +256,10 @@ mod tests {
 			vec![build_hist(), build_hist()], // hists
 			false // cyclic
 		);
-		assert_delta!(2.0, ds.bin_count[0], 0.0000000001);
-		assert_delta!(2.0, ds.bin_count[1], 0.0000000001);
-		assert_delta!(6.0, ds.bin_count[2], 0.0000000001);
-		assert_delta!(10.0, ds.bin_count[3], 0.0000000001);
-		assert_delta!(24.0, ds.bin_count[4], 0.0000000001);
+		assert_delta!(2.0, ds.get_weighted_bin_count(0), 0.0000000001);
+		assert_delta!(2.0, ds.get_weighted_bin_count(1), 0.0000000001);
+		assert_delta!(6.0, ds.get_weighted_bin_count(2), 0.0000000001);
+		assert_delta!(10.0, ds.get_weighted_bin_count(3), 0.0000000001);
+		assert_delta!(24.0, ds.get_weighted_bin_count(4), 0.0000000001);
 	}
 }

src/io.rs

@@ -158,20 +158,20 @@ fn read_window_file(window_file: &str, cfg: &Config) -> Result<Histogram> {
 }
 
 // Write WHAM calculation results to out_file.
-pub fn write_results(out_file: &str, ds: &Dataset, free: &Vec<f64>, prob: &Vec<f64>) -> Result<()> {
+pub fn write_results(out_file: &str, ds: &Dataset, free: &Vec<f64>, free_std: &Vec<f64>, prob: &Vec<f64>, prob_std: &Vec<f64>) -> Result<()> {
     let output = File::create(out_file)
         .chain_err(|| format!("Failed to create file with path {}", out_file))?;
     let mut buf = BufWriter::new(output);
 
     let header: String = (0..ds.dimens_lengths.len()).map(|d| format!("coord{}", d+1))
         .collect::<Vec<String>>().join("    ");
-    writeln!(buf, "#{}    {}    {}", header, "Free Energy", "Probability");
+    writeln!(buf, "#{}    {}    {}    {}    {}", header, "Free Energy", "+/-", "Probability", "+/-");
 
     for bin in 0..free.len() {
         let coords = ds.get_coords_for_bin(bin);
         let coords_str: String = coords.iter().map(|c| {format!("{:8.6}    ", c)})
             .collect::<Vec<String>>().join("\t");
-        writeln!(buf, "{}{:8.6}    {:8.6}", coords_str, free[bin], prob[bin])
+        writeln!(buf, "{}{:8.6}    {:8.6}    {:8.6}    {:8.6}", coords_str, free[bin], free_std[bin], prob[bin], prob_std[bin])
             .chain_err(|| "Failed to write to file.")?;
     }
     Ok(())
@@ -194,6 +194,7 @@ mod tests {
             temperature: 300.0,
             cyclic: false,
             output: "qwert".to_string(),
+            bootstrap: 0,
         }
     }
 

src/lib.rs

@@ -2,9 +2,12 @@
 
 #[macro_use]
 extern crate error_chain;
+extern crate rand;
+extern crate rgsl;
 
 pub mod io;
 pub mod histogram;
+pub mod error_analysis;
 
 use histogram::Dataset;
 use std::f64;
@@ -32,12 +35,13 @@ pub struct Config {
 	pub temperature: f64,
 	pub cyclic: bool,
 	pub output: String,
+	pub bootstrap: usize,
 }
 
 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={:?}", self.metadata_file, self.hist_min, self.hist_max, self.num_bins,
-                self.verbose, self.tolerance, self.max_iterations, self.temperature, self.cyclic)
+         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)
     }
 }
 
@@ -54,11 +58,12 @@ fn is_converged(old_F: &[f64], new_F: &[f64], tolerance: f64) -> bool {
 // This evaluates the first WHAM equation for each bin.
 fn calc_bin_probability(bin: usize, dataset: &Dataset, F: &[f64]) -> f64 {
     let mut denom_sum: f64 = 0.0;
+	let bin_count: f64 = dataset.get_weighted_bin_count(bin);
     for (window, h) in dataset.histograms.iter().enumerate() {
 		let bias = dataset.calc_bias(bin, window);
-        denom_sum += (h.num_points as f64) * bias * F[window];
+        denom_sum += (dataset.weights[window] * h.num_points as f64) * bias * F[window];
 	}
-    dataset.bin_count[bin] / denom_sum
+    bin_count / denom_sum
 }
 
 // estimate the bias offset F of the histogram based on given probabilities.
@@ -150,12 +155,23 @@ pub fn run(cfg: &Config) -> Result<()>{
 
     let (P, F, F_prev) = perform_wham(&cfg, &dataset)?;
 
+	let P_std: Vec<f64>;
+	let free_energy_std: Vec<f64>;
+	if cfg.bootstrap > 0 {
+		let error_est = error_analysis::run_bootstrap(&cfg, dataset.clone(), &P, cfg.bootstrap);
+		P_std = error_est.0;
+		free_energy_std = error_est.1;
+	} else {
+		P_std = vec![0.0; P.len()];
+		free_energy_std = vec![0.0; P.len()];
+	}
+
     // calculate free energy and dump state
     println!("Finished. Dumping final PMF");
 	let free_energy = calc_free_energy(&dataset, &P);
-    dump_state(&dataset, &F, &F_prev, &P, &free_energy);
+    dump_state(&dataset, &F, &F_prev, &P, &P_std, &free_energy, &free_energy_std);
 
-    io::write_results(&cfg.output, &dataset, &free_energy, &P)
+    io::write_results(&cfg.output, &dataset, &free_energy, &free_energy_std, &P, &P_std)
 		.chain_err(|| "Could not write results to output file")?;
 
     Ok(())
@@ -191,16 +207,17 @@ fn calc_free_energy(dataset: &Dataset, P: &[f64]) -> Vec<f64> {
     free_energy
 }
 
-fn dump_state(dataset: &Dataset, F: &[f64], F_prev: &[f64], P: &[f64], A: &[f64]) {
+fn dump_state(dataset: &Dataset, F: &[f64], F_prev: &[f64], P: &[f64], P_std: &[f64], A: &[f64], A_std: &[f64]) {
+	// TODO fix output of F/F_prev
 	let out = std::io::stdout();
     let mut lock = out.lock();
 	writeln!(lock, "# PMF");
-	writeln!(lock, "#bin\t\tFree Energy\t\tP(x)");
+	writeln!(lock, "#bin\t\tFree Energy\t\t+/-\t\tP(x)\t\t+/-");
 	for bin in 0..dataset.num_bins {
-		writeln!(lock, "{:9.5}\t{:9.5}\t{:9.5}", bin, A[bin], P[bin]);
+		writeln!(lock, "{:9.5}\t{:9.5}\t{:9.5}\t{:9.5}\t{:9.5}", bin, A[bin], A_std[bin], P[bin], P_std[bin]);
 	}
 	writeln!(lock, "# Bias offsets");
-	writeln!(lock, "#Window\t\tF\t\tdF");
+	writeln!(lock, "#Window\t\tF\t\tF_prev");
 	for window in 0..dataset.num_windows {
 		writeln!(lock, "{}\t{:9.5}\t{:8.8}", window, F[window], (F[window]-F_prev[window]).abs());
 	}

src/main.rs

@@ -28,6 +28,8 @@ fn cli() -> Result<Config> {
         .split(',').map(|x| { x.parse().unwrap() }).collect();
 	let num_bins: Vec<usize> = matches.value_of("bins").unwrap()
         .split(',').map(|x| { x.parse().unwrap() }).collect();
+	let bootstrap: usize = matches.value_of("bootstrap").unwrap_or("0").parse()
+		.chain_err(|| "Cannot parse bootstrap iteration.")?;
 
     if num_bins.len() != hist_max.len() || num_bins.len() != hist_max.len() {
         eprintln!("Input dimensions do not match (min: {}, max: {}, bins: {})",
@@ -38,7 +40,8 @@ fn cli() -> Result<Config> {
     let dimens = num_bins.len();
 
 	Ok(wham::Config{metadata_file, hist_min, hist_max, num_bins, dimens,
-		verbose, tolerance, max_iterations, temperature, cyclic, output})
+		verbose, tolerance, max_iterations, temperature, cyclic, output,
+		bootstrap})
 }
 
 fn main() {