Multithreading

Cargo.lock

@@ -16,6 +16,14 @@ dependencies = [
  "winapi 0.3.6 (registry+https://github.com/rust-lang/crates.io-index)",
 ]
 
+[[package]]
+name = "arrayvec"
+version = "0.4.7"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+dependencies = [
+ "nodrop 0.1.13 (registry+https://github.com/rust-lang/crates.io-index)",
+]
+
 [[package]]
 name = "atty"
 version = "0.2.11"
@@ -90,6 +98,42 @@ dependencies = [
  "bitflags 1.0.4 (registry+https://github.com/rust-lang/crates.io-index)",
 ]
 
+[[package]]
+name = "crossbeam-deque"
+version = "0.2.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+dependencies = [
+ "crossbeam-epoch 0.3.1 (registry+https://github.com/rust-lang/crates.io-index)",
+ "crossbeam-utils 0.2.2 (registry+https://github.com/rust-lang/crates.io-index)",
+]
+
+[[package]]
+name = "crossbeam-epoch"
+version = "0.3.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+dependencies = [
+ "arrayvec 0.4.7 (registry+https://github.com/rust-lang/crates.io-index)",
+ "cfg-if 0.1.6 (registry+https://github.com/rust-lang/crates.io-index)",
+ "crossbeam-utils 0.2.2 (registry+https://github.com/rust-lang/crates.io-index)",
+ "lazy_static 1.2.0 (registry+https://github.com/rust-lang/crates.io-index)",
+ "memoffset 0.2.1 (registry+https://github.com/rust-lang/crates.io-index)",
+ "nodrop 0.1.13 (registry+https://github.com/rust-lang/crates.io-index)",
+ "scopeguard 0.3.3 (registry+https://github.com/rust-lang/crates.io-index)",
+]
+
+[[package]]
+name = "crossbeam-utils"
+version = "0.2.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+dependencies = [
+ "cfg-if 0.1.6 (registry+https://github.com/rust-lang/crates.io-index)",
+]
+
+[[package]]
+name = "either"
+version = "1.5.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+
 [[package]]
 name = "error-chain"
 version = "0.12.0"
@@ -112,11 +156,34 @@ name = "fuchsia-zircon-sys"
 version = "0.3.3"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 
+[[package]]
+name = "lazy_static"
+version = "1.2.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+
 [[package]]
 name = "libc"
 version = "0.2.43"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 
+[[package]]
+name = "memoffset"
+version = "0.2.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+
+[[package]]
+name = "nodrop"
+version = "0.1.13"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+
+[[package]]
+name = "num_cpus"
+version = "1.8.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+dependencies = [
+ "libc 0.2.43 (registry+https://github.com/rust-lang/crates.io-index)",
+]
+
 [[package]]
 name = "pkg-config"
 version = "0.3.14"
@@ -147,6 +214,27 @@ name = "rand_core"
 version = "0.3.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 
+[[package]]
+name = "rayon"
+version = "1.0.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+dependencies = [
+ "crossbeam-deque 0.2.0 (registry+https://github.com/rust-lang/crates.io-index)",
+ "either 1.5.0 (registry+https://github.com/rust-lang/crates.io-index)",
+ "rayon-core 1.4.1 (registry+https://github.com/rust-lang/crates.io-index)",
+]
+
+[[package]]
+name = "rayon-core"
+version = "1.4.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+dependencies = [
+ "crossbeam-deque 0.2.0 (registry+https://github.com/rust-lang/crates.io-index)",
+ "lazy_static 1.2.0 (registry+https://github.com/rust-lang/crates.io-index)",
+ "libc 0.2.43 (registry+https://github.com/rust-lang/crates.io-index)",
+ "num_cpus 1.8.0 (registry+https://github.com/rust-lang/crates.io-index)",
+]
+
 [[package]]
 name = "redox_syscall"
 version = "0.1.40"
@@ -165,6 +253,11 @@ name = "rustc-demangle"
 version = "0.1.9"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 
+[[package]]
+name = "scopeguard"
+version = "0.3.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+
 [[package]]
 name = "strsim"
 version = "0.7.0"
@@ -200,12 +293,13 @@ source = "registry+https://github.com/rust-lang/crates.io-index"
 
 [[package]]
 name = "wham"
-version = "0.9.0"
+version = "0.9.1"
 dependencies = [
  "GSL 1.0.0 (registry+https://github.com/rust-lang/crates.io-index)",
  "clap 2.32.0 (registry+https://github.com/rust-lang/crates.io-index)",
  "error-chain 0.12.0 (registry+https://github.com/rust-lang/crates.io-index)",
  "rand 0.5.5 (registry+https://github.com/rust-lang/crates.io-index)",
+ "rayon 1.0.3 (registry+https://github.com/rust-lang/crates.io-index)",
 ]
 
 [[package]]
@@ -235,6 +329,7 @@ source = "registry+https://github.com/rust-lang/crates.io-index"
 [metadata]
 "checksum GSL 1.0.0 (registry+https://github.com/rust-lang/crates.io-index)" = "2e70167131757876cda6ab95869b57ad6b55c627ff28421f5935f55f9b97d85c"
 "checksum ansi_term 0.11.0 (registry+https://github.com/rust-lang/crates.io-index)" = "ee49baf6cb617b853aa8d93bf420db2383fab46d314482ca2803b40d5fde979b"
+"checksum arrayvec 0.4.7 (registry+https://github.com/rust-lang/crates.io-index)" = "a1e964f9e24d588183fcb43503abda40d288c8657dfc27311516ce2f05675aef"
 "checksum atty 0.2.11 (registry+https://github.com/rust-lang/crates.io-index)" = "9a7d5b8723950951411ee34d271d99dddcc2035a16ab25310ea2c8cfd4369652"
 "checksum backtrace 0.3.9 (registry+https://github.com/rust-lang/crates.io-index)" = "89a47830402e9981c5c41223151efcced65a0510c13097c769cede7efb34782a"
 "checksum backtrace-sys 0.1.24 (registry+https://github.com/rust-lang/crates.io-index)" = "c66d56ac8dabd07f6aacdaf633f4b8262f5b3601a810a0dcddffd5c22c69daa0"
@@ -244,17 +339,28 @@ source = "registry+https://github.com/rust-lang/crates.io-index"
 "checksum cfg-if 0.1.6 (registry+https://github.com/rust-lang/crates.io-index)" = "082bb9b28e00d3c9d39cc03e64ce4cea0f1bb9b3fde493f0cbc008472d22bdf4"
 "checksum clap 2.32.0 (registry+https://github.com/rust-lang/crates.io-index)" = "b957d88f4b6a63b9d70d5f454ac8011819c6efa7727858f458ab71c756ce2d3e"
 "checksum cloudabi 0.0.3 (registry+https://github.com/rust-lang/crates.io-index)" = "ddfc5b9aa5d4507acaf872de71051dfd0e309860e88966e1051e462a077aac4f"
+"checksum crossbeam-deque 0.2.0 (registry+https://github.com/rust-lang/crates.io-index)" = "f739f8c5363aca78cfb059edf753d8f0d36908c348f3d8d1503f03d8b75d9cf3"
+"checksum crossbeam-epoch 0.3.1 (registry+https://github.com/rust-lang/crates.io-index)" = "927121f5407de9956180ff5e936fe3cf4324279280001cd56b669d28ee7e9150"
+"checksum crossbeam-utils 0.2.2 (registry+https://github.com/rust-lang/crates.io-index)" = "2760899e32a1d58d5abb31129f8fae5de75220bc2176e77ff7c627ae45c918d9"
+"checksum either 1.5.0 (registry+https://github.com/rust-lang/crates.io-index)" = "3be565ca5c557d7f59e7cfcf1844f9e3033650c929c6566f511e8005f205c1d0"
 "checksum error-chain 0.12.0 (registry+https://github.com/rust-lang/crates.io-index)" = "07e791d3be96241c77c43846b665ef1384606da2cd2a48730abe606a12906e02"
 "checksum fuchsia-zircon 0.3.3 (registry+https://github.com/rust-lang/crates.io-index)" = "2e9763c69ebaae630ba35f74888db465e49e259ba1bc0eda7d06f4a067615d82"
 "checksum fuchsia-zircon-sys 0.3.3 (registry+https://github.com/rust-lang/crates.io-index)" = "3dcaa9ae7725d12cdb85b3ad99a434db70b468c09ded17e012d86b5c1010f7a7"
+"checksum lazy_static 1.2.0 (registry+https://github.com/rust-lang/crates.io-index)" = "a374c89b9db55895453a74c1e38861d9deec0b01b405a82516e9d5de4820dea1"
 "checksum libc 0.2.43 (registry+https://github.com/rust-lang/crates.io-index)" = "76e3a3ef172f1a0b9a9ff0dd1491ae5e6c948b94479a3021819ba7d860c8645d"
+"checksum memoffset 0.2.1 (registry+https://github.com/rust-lang/crates.io-index)" = "0f9dc261e2b62d7a622bf416ea3c5245cdd5d9a7fcc428c0d06804dfce1775b3"
+"checksum nodrop 0.1.13 (registry+https://github.com/rust-lang/crates.io-index)" = "2f9667ddcc6cc8a43afc9b7917599d7216aa09c463919ea32c59ed6cac8bc945"
+"checksum num_cpus 1.8.0 (registry+https://github.com/rust-lang/crates.io-index)" = "c51a3322e4bca9d212ad9a158a02abc6934d005490c054a2778df73a70aa0a30"
 "checksum pkg-config 0.3.14 (registry+https://github.com/rust-lang/crates.io-index)" = "676e8eb2b1b4c9043511a9b7bea0915320d7e502b0a079fb03f9635a5252b18c"
 "checksum rand 0.5.5 (registry+https://github.com/rust-lang/crates.io-index)" = "e464cd887e869cddcae8792a4ee31d23c7edd516700695608f5b98c67ee0131c"
 "checksum rand_core 0.2.2 (registry+https://github.com/rust-lang/crates.io-index)" = "1961a422c4d189dfb50ffa9320bf1f2a9bd54ecb92792fb9477f99a1045f3372"
 "checksum rand_core 0.3.0 (registry+https://github.com/rust-lang/crates.io-index)" = "0905b6b7079ec73b314d4c748701f6931eb79fd97c668caa3f1899b22b32c6db"
+"checksum rayon 1.0.3 (registry+https://github.com/rust-lang/crates.io-index)" = "373814f27745b2686b350dd261bfd24576a6fb0e2c5919b3a2b6005f820b0473"
+"checksum rayon-core 1.4.1 (registry+https://github.com/rust-lang/crates.io-index)" = "b055d1e92aba6877574d8fe604a63c8b5df60f60e5982bf7ccbb1338ea527356"
 "checksum redox_syscall 0.1.40 (registry+https://github.com/rust-lang/crates.io-index)" = "c214e91d3ecf43e9a4e41e578973adeb14b474f2bee858742d127af75a0112b1"
 "checksum redox_termios 0.1.1 (registry+https://github.com/rust-lang/crates.io-index)" = "7e891cfe48e9100a70a3b6eb652fef28920c117d366339687bd5576160db0f76"
 "checksum rustc-demangle 0.1.9 (registry+https://github.com/rust-lang/crates.io-index)" = "bcfe5b13211b4d78e5c2cadfebd7769197d95c639c35a50057eb4c05de811395"
+"checksum scopeguard 0.3.3 (registry+https://github.com/rust-lang/crates.io-index)" = "94258f53601af11e6a49f722422f6e3425c52b06245a5cf9bc09908b174f5e27"
 "checksum strsim 0.7.0 (registry+https://github.com/rust-lang/crates.io-index)" = "bb4f380125926a99e52bc279241539c018323fab05ad6368b56f93d9369ff550"
 "checksum termion 1.5.1 (registry+https://github.com/rust-lang/crates.io-index)" = "689a3bdfaab439fd92bc87df5c4c78417d3cbe537487274e9b0b2dce76e92096"
 "checksum textwrap 0.10.0 (registry+https://github.com/rust-lang/crates.io-index)" = "307686869c93e71f94da64286f9a9524c0f308a9e1c87a583de8e9c9039ad3f6"

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "wham"
-version = "0.9.0"
+version = "0.9.1"
 authors = ["D. Bauer <bauer@bio.tu-darmstadt.de>"]
 description = "An implementation of the weighted histogram analysis method"
 license = "GPL-3.0"
@@ -11,6 +11,7 @@ clap = {version="2.32.0", features=['yaml']}
 error-chain = "0.12.0"
 rand = "0.5.5"
 GSL = "1.0.0"
+rayon = "1.0.3"
 
 [profile.release]
 opt-level = 2

README.md

@@ -10,6 +10,7 @@ from umbrella sampling simulations. For more details on the method, I suggest Ro
 Features
 ---
 - Fast, especially for small systems
+- Multithreaded
 - Multidimensional
 - Error analysis
 - Unit tested
@@ -91,7 +92,6 @@ The example folder contains input and output files for two simple test systems:
 
 TODO
 ---
-- Multithreading (?)
 - Autocorrelation
 - Replica exchange
 

src/histogram.rs

@@ -1,5 +1,4 @@
 use std::fmt;
-use std::cell::RefCell;
 
 // One histogram
 #[derive(Debug,Clone)]
@@ -54,7 +53,7 @@ pub struct Dataset {
 	bias_fc: Vec<f64>,
 
 	// bias value cache
-	bias: RefCell<Vec<Option<f64>>>,
+	bias: Vec<f64>,
 
 	// histogram weight
 	pub weights: Vec<f64>,
@@ -64,9 +63,9 @@ 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 bias: Vec<f64> = vec![0.0; num_bins*num_windows];
 		let weights = vec![1.0; num_windows];
-		Dataset{
+		let mut ds = Dataset{
 			num_windows,
 			num_bins,
 			dimens_lengths,
@@ -80,8 +79,16 @@ impl Dataset {
 			bias_fc,
 			bias,
 			weights
+		};
+		for window in 0..num_windows {
+			for bin in 0..num_bins {
+				let ndx = window * num_bins + bin;
+				ds.bias[ndx] = ds.calc_bias(bin, window);
 			}
 		}
+		ds
+
+	}
 
 	pub fn new_weighted(ds: Dataset, weights: Vec<f64>) -> Dataset {
 		Dataset {
@@ -113,17 +120,16 @@ impl Dataset {
 		}).collect()
 	}
 
+	pub fn get_bias(&self, bin: usize, window: usize) -> f64 {
+		let ndx = window * self.num_bins + bin;
+		self.bias[ndx]
+	}
+
 	// Harmonic bias calculation: bias = 0.5*k(dx)^2
 	// if cyclic is true, lowest and highest bins are assumed to be
 	// neighbors. This returns exp(U/kT) instead of U for better performance.
-	pub fn calc_bias(&self, bin: usize, window: usize) -> f64 {
+	fn calc_bias(&self, bin: usize, window: usize) -> f64 {
-		let ndx = window * self.num_bins + bin;
-		let mut cache = self.bias.borrow_mut();
-		match cache[ndx] {
-			Some(val) => val,
-			None => {
 		let dimens = self.dimens_lengths.len();
-
 		// index of the bias value depends on the window und dimension
 		let bias_ndx: Vec<usize> = (0..dimens)
 			.map(|dimen| { window * dimens + dimen }).collect();
@@ -146,12 +152,8 @@ impl Dataset {
 			bias_sum += 0.5 * bias_fc[i] * dist * dist
 		}
 		let bias_sum = (-bias_sum/self.kT).exp();
-				cache[ndx] = Some(bias_sum);
 		bias_sum
 	}
-		}
-	}
-
 }
 
 impl fmt::Display for Dataset {
@@ -247,11 +249,11 @@ mod tests {
 		let ds = Dataset::new(
 			5, // num bins
 			vec![1],
-			vec![1.0], // bin width
+			vec![1.0, 1.0], // bin width
-			vec![0.0], // hist min
+			vec![0.0, 0.0], // hist min
-			vec![9.0], // hist max
+			vec![5.0, 5.0], // hist max
-			vec![7.5], // x0
+			vec![7.5, 7.5], // x0
-			vec![10.0], // fc
+			vec![10.0, 10.0], // fc
 			300.0*k_B, // kT
 			vec![build_hist(), build_hist()], // hists
 			false // cyclic

src/lib.rs

@@ -4,6 +4,7 @@
 extern crate error_chain;
 extern crate rand;
 extern crate rgsl;
+extern crate rayon;
 
 pub mod io;
 pub mod histogram;
@@ -13,6 +14,7 @@ use histogram::Dataset;
 use std::f64;
 use std::fmt;
 use std::io::prelude::*;
+use rayon::prelude::*;
 
 // init error chain
 pub mod errors { error_chain!{} }
@@ -60,7 +62,7 @@ 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);
+		let bias = dataset.get_bias(bin, window);
         denom_sum += (dataset.weights[window] * h.num_points as f64) * bias * F[window];
 	}
     bin_count / denom_sum
@@ -71,7 +73,7 @@ fn calc_bin_probability(bin: usize, dataset: &Dataset, F: &[f64]) -> f64 {
 fn calc_window_F(window: usize, dataset: &Dataset, P: &[f64]) -> f64 {
     let f: f64 = (0..dataset.num_bins).zip(P.iter()) // zip bins and P
         .map(|bin_and_prob: (usize, &f64)| {
-            let bias = dataset.calc_bias(bin_and_prob.0, window);
+            let bias = dataset.get_bias(bin_and_prob.0, window);
             bin_and_prob.1 * bias
         }).sum();
     1.0/f
@@ -80,18 +82,18 @@ fn calc_window_F(window: usize, dataset: &Dataset, P: &[f64]) -> f64 {
 // One full WHAM iteration includes calculation of new probabilities P and
 // new bias offsets F based on previous bias offsets F_prev. This updates
 // the values in vectors F and P
-fn perform_wham_iteration(dataset: &Dataset, F_prev: &[f64], F: &mut [f64], P: &mut [f64]) {
+fn perform_wham_iteration(dataset: &Dataset, F_prev: &[f64], F: &mut Vec<f64>, P: &mut Vec<f64>) {
 	// evaluate first WHAM equation for each bin to
-	// estimage probabilities based on previous offsets (F_prev)
+	// estimage probabilities based on previous offsets (F_prev))
-	for bin in 0..dataset.num_bins {
+	(0..dataset.num_bins).into_par_iter()
-		P[bin] = calc_bin_probability(bin, dataset, F_prev);
+		.map(|bin| { calc_bin_probability(bin, dataset, F_prev) })
-	}
+		.collect_into_vec(P);
 
 	// evaluate second WHAM equation for each window to
 	// estimate new bias offsets from propabilities
-	for window in 0..dataset.num_windows {
+	(0..dataset.num_windows).into_par_iter()
-		F[window] = calc_window_F(window, dataset, P);
+		.map(|window| {calc_window_F(window, dataset, P)} )
-	}
+		.collect_into_vec(F);
 }
 
 pub fn perform_wham(cfg: &Config, dataset: &Dataset) -> Result<(Vec<f64>, Vec<f64>, Vec<f64>)> {

tests/integration.rs

@@ -101,7 +101,7 @@ mod integration {
     fn wham_2d() {;
         Command::new("./target/debug/wham")
             .args(&["--bins", "100,100", "--max", "3.14,3.14", "--min", "-3.14,-3.14", "-T", "300", "--cyclic"])
-            .args(&["-f", "example/2d/metadata_unparseable1.dat"])
+            .args(&["-f", "example/2d/metadata.dat"])
             .args(&["-o", "/tmp/wham_test_2d.out"])
             .output()
             .expect("failed to execute process");