cleanup + version

v0.9.2
better error messages when reading timeseries
2026-09-10 22:25:31 +00:00 · 2018-11-30 13:06:58 +01:00 · 2018-11-30 12:58:53 +01:00 · 2018-11-30 12:58:24 +01:00 · 2018-11-18 12:09:03 +01:00 · 2018-11-18 12:04:58 +01:00
12 changed files with 241 additions and 64 deletions
--- a/.gitignore
+++ b/.gitignore
@@ -1,4 +1,4 @@
 target
 reference
 .idea
-bench.sh
+*bench.*
--- 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.2"
 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"
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -1,6 +1,6 @@
 [package]
 name = "wham"
-version = "0.9.0"
+version = "0.9.2"
 authors = ["D. Bauer <bauer@bio.tu-darmstadt.de>"]
 description = "An implementation of the weighted histogram analysis method"
 license = "GPL-3.0"
@@ -11,9 +11,10 @@ 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

 [features]
-default = ["GSL/v2"]
+default = ["GSL/v2"]
--- a/README.md
+++ b/README.md
@@ -1,4 +1,6 @@
-[![Build Status](https://travis-ci.com/danijoo/WHAM.svg?branch=master)](https://travis-ci.com/danijoo/WHAM) [![crates.io](https://img.shields.io/badge/crates.io-orange.svg?longCache=true)](https://www.crates.io/crates/wham)
+[![Build Status](https://travis-ci.com/danijoo/WHAM.svg?branch=master)](https://travis-ci.com/danijoo/WHAM) [![crates.io](https://img.shields.io/badge/crates.io-orange.svg?longCache=true)](https://www.crates.io/crates/wham) [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.1488598.svg)](https://doi.org/10.5281/zenodo.1488598)
+
+

 Weighted Histogram Analysis Method (WHAM)
 ===
@@ -10,6 +12,7 @@ from umbrella sampling simulations. For more details on the method, I suggest *R
 Features
 ---
 - Fast, especially for small systems
+- Multithreaded
 - Multidimensional
 - Error analysis
 - Unit tested
@@ -91,14 +94,15 @@ The example folder contains input and output files for two simple test systems:

 TODO
 ---
- Multithreading (?)
 - Autocorrelation
 - Replica exchange

-License 
+License & Citing
 ---
 WHAM is licensed under the GPL-3.0 license. Please read the LICENSE file in this
 repository for more information.

+There's no publication for this WHAM implementation. However, there is a citeabe DOI. If you use this software for your work, please consider citing it: *Bauer, D, WHAM - An efficient weighted histogram analysis implementation written in Rust, Zenodo.  https://doi.org/10.5281/zenodo.1488597*
+
 Parts of this work, especially some perfomance optimizations and the I/O format, are inspired by the
 implementation of A. Grossfield (*Grossfield, A, WHAM: the weighted histogram analysis method, http://membrane.urmc.rochester.edu/content/wham*).
--- 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<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,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<f64>) -> 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<usize> = (0..dimens)
+			.map(|dimen| { window * dimens + dimen }).collect();

-				// index of the bias value depends on the window und dimension
-				let bias_ndx: Vec<usize> = (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<f64> = bias_ndx.iter().map(|ndx| { self.bias_fc[*ndx] }).collect();
+		let bias_pos: Vec<f64> = 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<f64> = bias_ndx.iter().map(|ndx| { self.bias_fc[*ndx] }).collect();
-				let bias_pos: Vec<f64> = 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
--- a/src/io.rs
+++ b/src/io.rs
@@ -7,6 +7,7 @@ use std::io::{BufReader,BufWriter};
 use k_B;
 use std::path::Path;
 use super::errors::*;
+use f64;

 // Returns the path to path2 relative to path1
 // path1: "path/to/file.dat"
@@ -112,7 +113,7 @@ fn is_in_hist_boundaries(values: &Vec<f64>, cfg: &Config) -> bool {
 // parse a timeseries file into a histogram
 fn read_window_file(window_file: &str, cfg: &Config) -> Result<Histogram> {
 	let f = File::open(window_file)
-        .chain_err(|| format!("Failed to open sample data file {}", window_file))?;
+        .chain_err(|| format!("Failed to open sample data file {}.", window_file))?;
    let mut buf = BufReader::new(&f);

    // total number of bins is the product of all dimensions length
@@ -126,7 +127,9 @@ fn read_window_file(window_file: &str, cfg: &Config) -> Result<Histogram> {

    // read and parse each timeseries line
    let mut line = String::new();
+    let mut linecount = 0;
    while buf.read_line(&mut line).chain_err(|| "Failed to read line")? > 0 {
+        linecount += 1;
        // skip comments and empty lines
        if line.starts_with("#") || line.starts_with("@") || line.len() == 0 {
            line.clear();
@@ -134,12 +137,17 @@ fn read_window_file(window_file: &str, cfg: &Config) -> Result<Histogram> {
        }

        {
-            let mut split = line.split_whitespace();

-            split.next(); // skip time/step column
-            let values: Vec<f64> = (0..cfg.dimens).collect::<Vec<usize>>().iter().map(|_| {
-                split.next().unwrap().parse::<f64>().unwrap()
-            }).collect();
+            let split: Vec<&str> = line.split_whitespace().collect();
+            if split.len() < cfg.dimens+1 {
+                bail!(format!("Wrong number of columns in line {} of window file {}. Empty Line?.", linecount, window_file));
+            }
+
+            let mut values: Vec<f64> = vec![f64::NAN; cfg.dimens];
+            for i in 0..values.len() {
+                values[i] = split[i+1].parse::<f64>()
+                    .chain_err(|| format!("Failed to parse line {} of window file {}.", linecount, window_file))?;
+            }

            if is_in_hist_boundaries(&values, cfg) {
                let bin_indeces = (0..cfg.dimens).map(|dimen: usize| {
--- 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<f64>, P: &mut Vec<f64>) {
 	// 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<f64>, Vec<f64>, Vec<f64>)> {
--- a/tests/COLVAR_unparseable_1.xvg
+++ b/tests/COLVAR_unparseable_1.xvg
@@ -0,0 +1,11 @@
+84.840004  -2.878293
+84.860004  -2.827750
+84.880004  -2.962375
+84.900004  
+84.920004  -2.748589
+84.940004  -2.726197
+84.960004  -2.802982
+84.980004  -2.771980
+85.000004  -2.706482
+85.020004  -2.776983
+85.040004  -2.883180
--- a/tests/COLVAR_unparseable_2.xvg
+++ b/tests/COLVAR_unparseable_2.xvg
@@ -0,0 +1,11 @@
+84.840004  -2.878293
+84.860004  -2.827750
+84.880004  -2.962375
+84.900004   qwee
+84.920004  -2.748589
+84.940004  -2.726197
+84.960004  -2.802982
+84.980004  -2.771980
+85.000004  -2.706482
+85.020004  -2.776983
+85.040004  -2.883180
--- a/tests/integration.rs
+++ b/tests/integration.rs
@@ -77,6 +77,36 @@ mod integration {
        ));
    }

+    #[test]
+    fn unparseable_timeseries() {
+        let output = Command::new("./target/debug/wham")
+            .args(&["--bins", "100", "--min", "-3.0", "--max", "3.0", "-T", "300"])
+            .args(&["-f", "tests/metadata_unparseable3.dat"])
+            .args(&["-o", "/dev/null"])
+            .output()
+            .expect("failed to execute process");
+        let output = String::from_utf8_lossy(&output.stderr);
+        println!("{}", output);
+        assert!(output.to_string().contains(
+            "Failed to parse line"
+        ));
+    }
+
+    #[test]
+    fn unparseable_timeseries_empty() {
+        let output = Command::new("./target/debug/wham")
+            .args(&["--bins", "100", "--min", "-3.0", "--max", "3.0", "-T", "300"])
+            .args(&["-f", "tests/metadata_unparseable4.dat"])
+            .args(&["-o", "/dev/null"])
+            .output()
+            .expect("failed to execute process");
+        let output = String::from_utf8_lossy(&output.stderr);
+        println!("{}", output);
+        assert!(output.to_string().contains(
+            "Wrong number of columns in line"
+        ));
+    }
+
    #[test]
    fn wham_1d() {;
        Command::new("./target/debug/wham")
@@ -101,7 +131,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");
--- a/tests/metadata_unparseable3.dat
+++ b/tests/metadata_unparseable3.dat
@@ -0,0 +1 @@
+COLVAR_unparseable_2.xvg     -3.0       100
--- a/tests/metadata_unparseable4.dat
+++ b/tests/metadata_unparseable4.dat
@@ -0,0 +1 @@
+COLVAR_unparseable_1.xvg     -3.0       100
Author	SHA1	Message	Date
Daniel Bauer	99f2c9069d	cleanup + version	2018-11-30 13:06:58 +01:00
Daniel Bauer	7109d27543	v0.9.2	2018-11-30 12:58:53 +01:00
Daniel Bauer	c03be1723e	better error messages when reading timeseries	2018-11-30 12:58:24 +01:00
Daniel Bauer	4aa3552f1d	gitignore benchmark files	2018-11-18 12:09:03 +01:00
Daniel Bauer	5650d2a65a	Merge branch 'master' of github.com:danijoo/WHAM	2018-11-18 12:04:58 +01:00
Daniel Bauer	6ed0e4c49c	Multithreading	2018-11-18 12:04:27 +01:00
Daniel Bauer	5caa2b5f20	DOI	2018-11-15 10:19:39 +01:00