diff --git a/Cargo.lock b/Cargo.lock index a605004..76baf0b 100644 --- a/Cargo.lock +++ b/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" diff --git a/Cargo.toml b/Cargo.toml index 97f4bb1..5dfdf4a 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "wham" -version = "0.9.0" +version = "0.9.1" authors = ["D. Bauer "] 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 diff --git a/README.md b/README.md index c15eb84..6e307e0 100644 --- a/README.md +++ b/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 diff --git a/src/histogram.rs b/src/histogram.rs index 8ca9f9d..6c27d85 100644 --- a/src/histogram.rs +++ b/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, // bias value cache - bias: RefCell>>, + bias: Vec, // histogram weight pub weights: Vec, @@ -64,9 +63,9 @@ impl Dataset { pub fn new(num_bins: usize, dimens_lengths: Vec, bin_width: Vec, hist_min: Vec, hist_max: Vec, bias_pos: Vec, bias_fc: Vec, kT: f64, histograms: Vec, cyclic: bool) -> Dataset { let num_windows = histograms.len(); - let bias: RefCell>> = RefCell::new(vec![None; num_bins*num_windows]); + let bias: Vec = 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,7 +79,15 @@ 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) -> Dataset { @@ -113,45 +120,40 @@ 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 { - 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(); + fn calc_bias(&self, bin: usize, window: usize) -> f64 { + let dimens = self.dimens_lengths.len(); + // index of the bias value depends on the window und dimension + let bias_ndx: Vec = (0..dimens) + .map(|dimen| { window * dimens + dimen }).collect(); - // index of the bias value depends on the window und dimension - let bias_ndx: Vec = (0..dimens) - .map(|dimen| { window * dimens + dimen }).collect(); + // find the N coords, force constants and bias coords + let coord = self.get_coords_for_bin(bin); + let bias_fc: Vec = bias_ndx.iter().map(|ndx| { self.bias_fc[*ndx] }).collect(); + let bias_pos: Vec = bias_ndx.iter().map(|ndx| { self.bias_pos[*ndx] }).collect(); - // find the N coords, force constants and bias coords - let coord = self.get_coords_for_bin(bin); - let bias_fc: Vec = bias_ndx.iter().map(|ndx| { self.bias_fc[*ndx] }).collect(); - let bias_pos: Vec = bias_ndx.iter().map(|ndx| { self.bias_pos[*ndx] }).collect(); - - let mut bias_sum = 0.0; - for i in 0..dimens { - let mut dist = (coord[i] - bias_pos[i]).abs(); - if self.cyclic { // periodic conditions - let hist_len = self.hist_max[i] - self.hist_min[i]; - if dist > 0.5 * hist_len { - dist -= hist_len; - } - } - // store exp(U/kT) for better performance - bias_sum += 0.5 * bias_fc[i] * dist * dist + let mut bias_sum = 0.0; + for i in 0..dimens { + let mut dist = (coord[i] - bias_pos[i]).abs(); + if self.cyclic { // periodic conditions + let hist_len = self.hist_max[i] - self.hist_min[i]; + if dist > 0.5 * hist_len { + dist -= hist_len; } - let bias_sum = (-bias_sum/self.kT).exp(); - cache[ndx] = Some(bias_sum); - bias_sum } + // store exp(U/kT) for better performance + bias_sum += 0.5 * bias_fc[i] * dist * dist } + let bias_sum = (-bias_sum/self.kT).exp(); + 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![0.0], // hist min - vec![9.0], // hist max - vec![7.5], // x0 - vec![10.0], // fc + vec![1.0, 1.0], // bin width + vec![0.0, 0.0], // hist min + vec![5.0, 5.0], // hist max + vec![7.5, 7.5], // x0 + vec![10.0, 10.0], // fc 300.0*k_B, // kT vec![build_hist(), build_hist()], // hists false // cyclic diff --git a/src/lib.rs b/src/lib.rs index f022de6..9d9cccc 100644 --- a/src/lib.rs +++ b/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!{} } @@ -57,10 +59,10 @@ fn is_converged(old_F: &[f64], new_F: &[f64], tolerance: f64) -> bool { // estimate the probability of a bin of the histogram set based on given bias offsets (F) // 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 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, P: &mut Vec) { // evaluate first WHAM equation for each bin to - // estimage probabilities based on previous offsets (F_prev) - for bin in 0..dataset.num_bins { - P[bin] = calc_bin_probability(bin, dataset, F_prev); - } + // estimage probabilities based on previous offsets (F_prev)) + (0..dataset.num_bins).into_par_iter() + .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 { - F[window] = calc_window_F(window, dataset, P); - } + (0..dataset.num_windows).into_par_iter() + .map(|window| {calc_window_F(window, dataset, P)} ) + .collect_into_vec(F); } pub fn perform_wham(cfg: &Config, dataset: &Dataset) -> Result<(Vec, Vec, Vec)> { diff --git a/tests/integration.rs b/tests/integration.rs index d7cd8c7..1be8480 100644 --- a/tests/integration.rs +++ b/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");