version up

* calculate autocorrelation stats_ineff and tau
* uncorr flag
* read timeseries into vector
* uncorrelate data
* README
* Update README.md
* Update README.md

Co-authored-by: Daniel Bauer <bauer@cbs.tu-darmstadt.de>

Cargo.lock

@@ -1,16 +1,5 @@
 # This file is automatically @generated by Cargo.
 # It is not intended for manual editing.
-[[package]]
-name = "GSL"
-version = "1.1.0"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "7830156ea389bcbbdc8f01bf140b609b892bf7cbd0ec6ccf9957ea2be6f25ad3"
-dependencies = [
- "c_vec",
- "libc",
- "pkg-config",
-]
-
 [[package]]
 name = "addr2line"
 version = "0.13.0"
@@ -78,12 +67,6 @@ version = "1.2.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "cf1de2fe8c75bc145a2f577add951f8134889b4795d47466a54a5c846d691693"
 
-[[package]]
-name = "c_vec"
-version = "1.0.12"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "aa9e1d9f7d49e289f36f19effbf3d5a5e30163ecf9c7a3c9be94d5374dec5b9a"
-
 [[package]]
 name = "cfg-if"
 version = "0.1.10"
@@ -215,9 +198,9 @@ checksum = "e2abad23fbc42b3700f2f279844dc832adb2b2eb069b2df918f455c4e18cc646"
 
 [[package]]
 name = "libc"
-version = "0.2.79"
+version = "0.2.80"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "2448f6066e80e3bfc792e9c98bf705b4b0fc6e8ef5b43e5889aff0eaa9c58743"
+checksum = "4d58d1b70b004888f764dfbf6a26a3b0342a1632d33968e4a179d8011c760614"
 
 [[package]]
 name = "memoffset"
@@ -254,12 +237,6 @@ version = "0.21.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "37fd5004feb2ce328a52b0b3d01dbf4ffff72583493900ed15f22d4111c51693"
 
-[[package]]
-name = "pkg-config"
-version = "0.3.19"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "3831453b3449ceb48b6d9c7ad7c96d5ea673e9b470a1dc578c2ce6521230884c"
-
 [[package]]
 name = "ppv-lite86"
 version = "0.2.9"
@@ -385,9 +362,8 @@ checksum = "cccddf32554fecc6acb585f82a32a72e28b48f8c4c1883ddfeeeaa96f7d8e519"
 
 [[package]]
 name = "wham"
-version = "0.9.9"
+version = "1.0.0"
 dependencies = [
- "GSL",
  "assert_approx_eq",
  "clap",
  "error-chain",

Cargo.toml

@@ -1,19 +1,18 @@
 [package]
 name = "wham"
-version = "0.9.9"
+version = "1.0.0"
 authors = ["Daniel Bauer <bauer@cbs.tu-darmstadt.de>"]
 description = "An implementation of the weighted histogram analysis method"
 license = "GPL-3.0"
 repository = "https://github.com/danijoo/WHAM"
 readme = "README.md"
-categories = ["science", "command-line-utilities", "molecular-dynamics", "algorithms"]
+categories = ["science", "command-line-utilities", "algorithms"]
 keywords = ["math", "statistics", "histogram", "bioinformatics", "molecular-dynamics"]
 
 [dependencies]
 clap = {version="2.32.0", features=['yaml']}
 error-chain = "0.12.0"
 rand = "0.7.*"
-GSL = "1.1"
 rayon = "1.0.3"
 
 [dev-dependencies]

README.md

@@ -12,9 +12,10 @@ from umbrella sampling simulations. For more details on the method, I suggest *R
 Features
 ---
 - Fast, especially for small systems
-- Multithreaded
-- Multidimensional
-- Error analysis
+- Multithreaded (automatically runs on all available cores) 
+- Multidimensional (any number of collective variables are possible)
+- Autocorrelation to remove correlated samples
+- Error analysis via bootstrapping
 - Unit tested
 
 Installation
@@ -39,7 +40,7 @@ wham has a convenient command line interface. You can see all options with
 ```wham -h```:
 
 ```
-wham 0.9.9
+wham 1.0.0
 D. Bauer <bauer@bio.tu-darmstadt.de>
 wham is a fast implementation of the weighted histogram analysis method (WHAM) written in Rust. It currently supports
 potential of mean force (PMF) calculations in multiple dimensions at constant temperature.
@@ -61,12 +62,14 @@ The first column will be ignored and is followed by N reaction coordinates x.
 Shipped under the GPLv3 license.
 
 USAGE:
-    wham [FLAGS] [OPTIONS] --bins <BINS> --max <HIST_MAX> --file <METADATA> --min <HIST_MIN> --temperature <temperature>
+   wham [FLAGS] [OPTIONS] --bins <BINS> --max <HIST_MAX> --file <METADATA> --min <HIST_MIN> --temperature <temperature>
 
 FLAGS:
     -c, --cyclic     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.
     -h, --help       Prints help information
+    -g, --uncorr     Estimates statistical inefficiency of each timeseries via autocorrelation and removes correlated
+                     samples (default is off).
     -V, --version    Prints version information
     -v, --verbose    Enables verbose output.
 
@@ -85,10 +88,18 @@ OPTIONS:
     -T, --temperature <temperature>    WHAM temperature in Kelvin.
     -t, --tolerance <TOLERANCE>        Abortion criteria for WHAM calculation. WHAM stops if abs(F_new - F_old) <
                                        tolerance (defaults to 0.000001).
-
 ```
 
-To run the two dimensional example (simulation of dialanine phi and psi angle):
+Examples
+---
+The example folder contains input and output files for two simple test systems:
+
+- 1d_cyclic: Phi torsion angle of dialanine in vaccum
+- 2d_cyclic: Phi and psi torsion angles of the same system
+
+The command below will run the two dimensional example (simulation of dialanine phi and psi angle) and calculate the free energy based on the two collective variables
+in the range of -3.14 to 3.14, with 100 bins in each dimension and periodic collective variables:
+ 
 ```bash
 wham --max 3.14,3.14 --min -3.14,-3.14 -T 300 --bins 100,100 --cyclic -f example/2d/metadata.dat       
 > Supplied WHAM options: Metadata=example/2d/metadata.dat, hist_min=[-3.14, -3.14], hist_max=[3.14, 3.14], bins=[100, 100] verbose=false, tolerance=0.000001, iterations=100000, temperature=300, cyclic=true
@@ -105,7 +116,6 @@ wham --max 3.14,3.14 --min -3.14,-3.14 -T 300 --bins 100,100 --cyclic -f example
 ```
 After convergence, final bias offsets (F) and the free energy will be dumped to stdout and the output file is written.
 
-
 The output file contains the free energy and probability for each bin. Probabilities are normalized to sum to P=1.0 and
 the smallest free energy is set to 0 (with other free energies based on that).
 ```
@@ -125,7 +135,7 @@ the smallest free energy is set to 0 (with other free energies based on that).
 Error analysis
 ---
 WHAM can perform error analysis using the bayesian bootstrapping method. Every simulation window is assumed to be an
-individual set of data point. By calculating probabilities N times with randomly assigned weights for each window,
+individual set of data points. By calculating probabilities N times with randomly assigned weights for each window,
 one can estimate the error as standard deviation between the N bootstrapping runs. For more details see
 *Van der Spoel, D. et al. (2010). g_wham—A Free Weighted Histogram Analysis Implementation Including Robust Error and
 Autocorrelation Estimates, JCTC, 6(12), 3713-3720*.
@@ -134,17 +144,20 @@ To perform bayesian bootstrapping in WHAM, use the ```-bt <RUNS>``` flag to perf
 runs. The error estimates of bin probabilities and free energy will be given as standard error (SE) in a 
 separate column (+/-) in the output file. If no error analysis is performed, these columns are set to 0.0.
 
-Examples
+Autocorrelation analysis
 ---
-The example folder contains input and output files for two simple test systems:
+With the ```--uncorr``` flag, WHAM calculates the autocorrelation time ```tau``` for all timeseries and all collective
+variables. Timeseries are then filtered based on their highest autocorrelation time to remove correlated samples from
+the dataset. This reduces the number of data points but can improve the accuracy of the result.
 
-- 1d_cyclic: Phi torsion angle of dialanine in vaccum
-- 2d_cyclic: Phi and psi torsion angles of the same system
+For filtering, the statistical inefficiency `g` is calculated: ```g = 1 + 2*tau```, and only every `g`th element of the
+timeseries is used for unbiasing. A more detailed description of the method can be found in
+*Chodera, J.D. et al. (2007). Use of the weighted histogram analysis method for the analysis of simulated and parallel
+tempering simulations, JCTC 3(1):26-41*
 
 
 TODO
 ---
-- Autocorrelation
 - Replica exchange
 
 License & Citing
@@ -152,7 +165,7 @@ 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*
+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*).

example/1d_cyclic/wham_uncorrelated.out

@@ -0,0 +1,101 @@
+#coord1    Free Energy    +/-    Probability    +/-
+-3.110177    7.531315    0.000000    0.003080    0.000000
+-3.047345    5.690157    0.000000    0.006443    0.000000
+-2.984513    4.243063    0.000000    0.011509    0.000000
+-2.921681    3.334686    0.000000    0.016564    0.000000
+-2.858849    2.277349    0.000000    0.025309    0.000000
+-2.796017    1.723296    0.000000    0.031604    0.000000
+-2.733186    1.246264    0.000000    0.038265    0.000000
+-2.670354    1.099867    0.000000    0.040578    0.000000
+-2.607522    0.771910    0.000000    0.046279    0.000000
+-2.544690    0.770616    0.000000    0.046303    0.000000
+-2.481858    1.265507    0.000000    0.037971    0.000000
+-2.419026    1.562335    0.000000    0.033711    0.000000
+-2.356194    1.891577    0.000000    0.029542    0.000000
+-2.293363    2.227858    0.000000    0.025816    0.000000
+-2.230531    2.488355    0.000000    0.023256    0.000000
+-2.167699    2.502265    0.000000    0.023127    0.000000
+-2.104867    2.358037    0.000000    0.024503    0.000000
+-2.042035    2.278147    0.000000    0.025301    0.000000
+-1.979203    2.974067    0.000000    0.019141    0.000000
+-1.916372    2.696600    0.000000    0.021393    0.000000
+-1.853540    2.361827    0.000000    0.024466    0.000000
+-1.790708    1.516746    0.000000    0.034332    0.000000
+-1.727876    1.526829    0.000000    0.034194    0.000000
+-1.665044    0.884114    0.000000    0.044244    0.000000
+-1.602212    0.323912    0.000000    0.055385    0.000000
+-1.539380    0.197985    0.000000    0.058253    0.000000
+-1.476549    0.000000    0.000000    0.063065    0.000000
+-1.413717    0.458247    0.000000    0.052481    0.000000
+-1.350885    1.389410    0.000000    0.036131    0.000000
+-1.288053    2.386522    0.000000    0.024225    0.000000
+-1.225221    3.743253    0.000000    0.014062    0.000000
+-1.162389    5.566654    0.000000    0.006770    0.000000
+-1.099557    7.822800    0.000000    0.002740    0.000000
+-1.036726    10.128719    0.000000    0.001087    0.000000
+-0.973894    12.199246    0.000000    0.000474    0.000000
+-0.911062    14.488129    0.000000    0.000189    0.000000
+-0.848230    16.902310    0.000000    0.000072    0.000000
+-0.785398    18.910200    0.000000    0.000032    0.000000
+-0.722566    21.241681    0.000000    0.000013    0.000000
+-0.659734    22.706373    0.000000    0.000007    0.000000
+-0.596903    24.531129    0.000000    0.000003    0.000000
+-0.534071    25.936227    0.000000    0.000002    0.000000
+-0.471239    27.000262    0.000000    0.000001    0.000000
+-0.408407    28.673293    0.000000    0.000001    0.000000
+-0.345575    29.335203    0.000000    0.000000    0.000000
+-0.282743    30.841118    0.000000    0.000000    0.000000
+-0.219911    31.983859    0.000000    0.000000    0.000000
+-0.157080    32.144015    0.000000    0.000000    0.000000
+-0.094248    33.885395    0.000000    0.000000    0.000000
+-0.031416    33.783105    0.000000    0.000000    0.000000
+0.031416    34.243727    0.000000    0.000000    0.000000
+0.094248    33.975567    0.000000    0.000000    0.000000
+0.157080    32.994787    0.000000    0.000000    0.000000
+0.219911    32.607398    0.000000    0.000000    0.000000
+0.282743    31.401902    0.000000    0.000000    0.000000
+0.345575    29.911670    0.000000    0.000000    0.000000
+0.408407    28.603574    0.000000    0.000001    0.000000
+0.471239    26.925443    0.000000    0.000001    0.000000
+0.534071    25.298070    0.000000    0.000002    0.000000
+0.596903    23.638560    0.000000    0.000005    0.000000
+0.659734    21.156231    0.000000    0.000013    0.000000
+0.722566    19.126480    0.000000    0.000029    0.000000
+0.785398    17.351953    0.000000    0.000060    0.000000
+0.848230    15.135525    0.000000    0.000146    0.000000
+0.911062    13.188112    0.000000    0.000319    0.000000
+0.973894    11.536983    0.000000    0.000618    0.000000
+1.036726    10.158328    0.000000    0.001074    0.000000
+1.099557    9.109036    0.000000    0.001636    0.000000
+1.162389    8.282343    0.000000    0.002279    0.000000
+1.225221    8.022102    0.000000    0.002530    0.000000
+1.288053    8.162415    0.000000    0.002391    0.000000
+1.350885    8.600135    0.000000    0.002006    0.000000
+1.413717    9.837348    0.000000    0.001222    0.000000
+1.476549    11.363156    0.000000    0.000663    0.000000
+1.539380    13.077849    0.000000    0.000333    0.000000
+1.602212    15.353594    0.000000    0.000134    0.000000
+1.665044    17.565051    0.000000    0.000055    0.000000
+1.727876    19.710884    0.000000    0.000023    0.000000
+1.790708    21.721260    0.000000    0.000010    0.000000
+1.853540    23.567649    0.000000    0.000005    0.000000
+1.916372    25.008817    0.000000    0.000003    0.000000
+1.979203    26.405367    0.000000    0.000002    0.000000
+2.042035    28.070821    0.000000    0.000001    0.000000
+2.104867    28.877213    0.000000    0.000001    0.000000
+2.167699    29.378146    0.000000    0.000000    0.000000
+2.230531    31.093267    0.000000    0.000000    0.000000
+2.293363    30.704994    0.000000    0.000000    0.000000
+2.356194    30.563093    0.000000    0.000000    0.000000
+2.419026    31.215952    0.000000    0.000000    0.000000
+2.481858    30.331416    0.000000    0.000000    0.000000
+2.544690    29.005123    0.000000    0.000001    0.000000
+2.607522    27.674618    0.000000    0.000001    0.000000
+2.670354    24.911788    0.000000    0.000003    0.000000
+2.733186    22.746417    0.000000    0.000007    0.000000
+2.796017    20.608288    0.000000    0.000016    0.000000
+2.858849    18.018280    0.000000    0.000046    0.000000
+2.921681    15.949215    0.000000    0.000105    0.000000
+2.984513    13.617809    0.000000    0.000268    0.000000
+3.047345    11.432193    0.000000    0.000645    0.000000
+3.110177    9.461494    0.000000    0.001420    0.000000

src/cli.yml

@@ -1,6 +1,6 @@
 name: wham
-version: "0.9.9"
-author: D. Bauer <bauer@bio.tu-darmstadt.de>
+version: "1.0.0"
+author: D. Bauer <bauer@cbs.tu-darmstadt.de>
 about: |
   wham is a fast implementation of the weighted histogram analysis method (WHAM) written in Rust. It currently supports potential of mean force (PMF) calculations in multiple dimensions at constant temperature.
 
@@ -101,3 +101,9 @@ args:
         help: Skip rows in timeseries with an index larger than this value (defaults to 1e+20)
         takes_value: true
         required: false
+    - uncorr:
+        short: g
+        long: uncorr
+        help: Estimates statistical inefficiency of each timeseries via autocorrelation and removes correlated samples (default is off).
+        takes_value: false
+        required: false 

src/correlation_analysis.rs

@@ -0,0 +1,102 @@
+use super::statistics;
+
+// calculates the statistical inefficiency g of the given timeseries
+// the quantity g can be thought of: N/g is the number of uncorrelated
+// configurations in the timeseries, where samples are separated by
+// the a multiple of g 
+// For details, see "Chodera et al. (2007). Use of a Weighted Histogram Analysis
+// Method for the Analysis of Simulated and Parallel Tempering Simulations, JCTC"
+pub fn statistical_ineff(timeseries: &[f64]) -> f64 {
+    let n = timeseries.len();
+    let autocorr = autocorrelation(timeseries);
+
+    let mut g = 1.0;
+    for t in 1..(n-1) {
+        let c = autocorr[t-1];
+        if c <= 0.0 {
+            break;
+        }
+        g = g + (2.0*c*(1.0-t as f64/n as f64))
+    }
+    if g < 1.0 {
+        1.0
+    } else {
+        g
+    }
+}
+
+// calculates the autocorrelation of a simeseries
+fn autocorrelation(timeseries: &[f64]) -> Vec<f64> {
+    let n = timeseries.len();
+    let mean = statistics::mean(timeseries);
+    let d_mean = timeseries.iter().map(|x| x-mean).collect::<Vec<f64>>();
+    let cov = statistics::autocov(timeseries);
+
+    let mut autocorr = Vec::new();
+    for t in 1..(n-1) {
+        let tmp = d_mean[0..n-t].iter().zip(d_mean[t..n].iter()).map(|(x, y)| x*y);
+        let c: f64 = tmp.map(|x| x+x).sum::<f64>() / (2.0 * (n as f64-t as f64)*cov);
+        autocorr.push(c);
+    }
+    autocorr
+}
+
+// The autocorrelation time of a timeseries can be deduced from the
+// `statistical_ineff` by (g-1)/2.0 
+pub fn autocorrelation_time(g: f64) -> f64 {
+    (g - 1.0) / 2.0
+}
+
+#[cfg(test)]
+mod tests {
+    use std::io::{BufRead, BufReader};
+    use std::fs::File;
+    
+    fn read_timeseries(filename: &str) -> Vec<f64> {
+        let mut timeseries: Vec<f64> = Vec::new();
+        let file = File::open(filename).unwrap();
+        let reader = BufReader::new(&file);
+        for line in reader.lines() {
+            let val = line.unwrap().split_whitespace()
+                .collect::<Vec<&str>>()[1].parse::<f64>().unwrap();
+            timeseries.push(val);
+        }
+        timeseries
+
+    }
+
+    #[test]
+    fn autocorrelation() {
+        let timeseries = read_timeseries("example/1d_cyclic/COLVAR-2.5.xvg");
+        let autocorr = super::autocorrelation(&timeseries);
+        let expected = [
+            0.691_900_865_597_914_3, 0.533_171_939_967_135_5,
+            0.204_726_209_564_635_89, -0.002_850_876_458_920_514,
+            -0.138_428_500_779_381_46, -0.265_292_355_297_323_2,
+            -0.314_271_982_722_353_85, -0.261_750_515_155_769_3,
+            -0.205_948_647_302_903_38, -0.133_100_190_918_118_12,
+            -0.188_756_890_119_342_6, -0.194_493_662_542_493_3,
+            -0.196_359_967_318_946_1, -0.113_915_878_338_380_03];
+        for (actual, expected) in autocorr.iter().zip(expected.iter()) {
+            assert!((actual-expected).abs() < 0.001);
+        }
+    }
+
+    #[test]
+    fn statistical_ineff() {
+        let timeseries = read_timeseries("example/1d_cyclic/COLVAR-2.5.xvg");
+        let g = super::statistical_ineff(&timeseries);
+        println!("{:?}", g);
+        assert!((g - 3.859).abs() < 0.001)
+    }
+
+    #[test]
+    fn autocorrelation_time() {
+        let timeseries = read_timeseries("example/1d_cyclic/COLVAR-2.5.xvg");
+        let g = super::statistical_ineff(&timeseries);
+        let tau = super::autocorrelation_time(g);
+        println!("{:?}", tau);
+        assert!((tau - 1.430).abs() < 0.001)
+    }
+
+}

src/error_analysis.rs

@@ -2,7 +2,7 @@ use rand::prelude::*;
 use super::histogram::{Dataset};
 use super::perform_wham;
 use super::{Config,calc_free_energy};
-use rgsl::statistics;
+use super::statistics;
 
 // returns a set of num_windows continious weights by
 // a) generate num_windows-1 random variables and sort them
@@ -48,7 +48,7 @@ pub fn run_bootstrap(cfg: &Config, ds: Dataset, num_runs: usize) -> (Vec<f64>,Ve
     let mut P_se = vec![0.0; ds.num_bins];
     for bin in 0..ds.num_bins {
         let Ps = bootstrapped_Ps.iter().map(|window| window[bin]).collect::<Vec<f64>>();
-        P_se[bin] = statistics::sd(&Ps, 1, num_runs)/(num_runs as f64).sqrt();
+        P_se[bin] = statistics::sd(&Ps)/(num_runs as f64).sqrt();
     }
 
     // SE of A
@@ -60,13 +60,13 @@ pub fn run_bootstrap(cfg: &Config, ds: Dataset, num_runs: usize) -> (Vec<f64>,Ve
     let mut A_se = vec![0.0; ds.num_bins];
     for bin in 0..ds.num_bins {
         let As = bootstrapped_As.iter().map(|window| window[bin]).collect::<Vec<f64>>();
-        A_se[bin] = statistics::sd(&As, 1, num_runs)/(num_runs as f64).sqrt();
+        A_se[bin] = statistics::sd(&As)/(num_runs as f64).sqrt();
     }
 
     (P_se, A_se)
 }
 
-#[cfg(tests)]
+#[cfg(test)]
 mod tests {
     use super::*;
     use super::super::k_B;
@@ -99,8 +99,9 @@ mod tests {
 
     #[test]
     fn random_weights() {
+        let mut rng = StdRng::from_entropy();
         let num_windows = 5;
-        let weights = generate_random_weights(num_windows);
+        let weights = generate_random_weights(num_windows, &mut rng);
         assert_eq!(num_windows, weights.len());
         for w in weights {
             assert!(0.0 < w && w < 1.0);
@@ -109,8 +110,9 @@ mod tests {
 
     #[test]
     fn random_weighted_dataset() {
+        let mut rng = StdRng::from_entropy();
         let ds = build_hist_set();
-        let rnd_weights_ds = generate_random_weighted_dataset(ds);
+        let rnd_weights_ds = generate_random_weighted_dataset(ds, &mut rng);
         println!("{:?}", rnd_weights_ds.weights);
         for w in rnd_weights_ds.weights {
             assert!(w > 0.0);

src/io.rs

@@ -1,6 +1,7 @@
 use super::histogram::Dataset;
 use super::histogram::Histogram;
 use super::Config;
+use super::correlation_analysis::{statistical_ineff, autocorrelation_time};
 use std::fs::File;
 use std::io::prelude::*;
 use std::io::{BufReader,BufWriter};
@@ -30,6 +31,7 @@ pub fn read_data(cfg: &Config) -> Result<Dataset> {
 	let mut bias_pos: Vec<f64> = Vec::new();
 	let mut bias_fc: Vec<f64> = Vec::new();
     let mut histograms: Vec<Histogram> = Vec::new();
+    let mut timeseries_lengths: Vec<usize> = Vec::new();
 
 	let kT = cfg.temperature * k_B;
     let bin_width: Vec<f64> = (0..cfg.dimens).map(|idx| {
@@ -57,12 +59,13 @@ pub fn read_data(cfg: &Config) -> Result<Dataset> {
 
         // parse histogram data
         let path = get_relative_path(&cfg.metadata_file, split[0]);
-        let h = read_window_file(&path, cfg)
+        let (h, timeseries_inital_length) = read_window_file(&path, cfg)
             .chain_err(|| format!("Failed to parse process data file {}", &path))?;
         if h.num_points == 0 {
             bail!(format!("No data points in histogram boundaries: {}", &path))
         }
         histograms.push(h);
+        timeseries_lengths.push(timeseries_inital_length);
         vprintln(format!("{}, {} data points added.", &path,
                          histograms.last().unwrap().num_points), cfg.verbose);
 
@@ -80,6 +83,20 @@ pub fn read_data(cfg: &Config) -> Result<Dataset> {
     }
     
     if !histograms.is_empty() {
+        if cfg.uncorr {
+            println!("Timeseries Correlation");
+            println!();
+            println!("Window\t\tN\t\tN_uncorr\tN/N_uncorr");
+            for (idx, (n, h)) in timeseries_lengths.iter().zip(histograms.iter()).enumerate() {
+                println!("{:?}\t\t{:?}\t\t{:?}\t\t{:.2}",
+                    idx+1, n, h.num_points, h.num_points as f64 / *n as f64);
+            }
+            let total_n = timeseries_lengths.iter().sum::<usize>() as f64;
+            let total_h = histograms.iter().map(|h| h.num_points).sum::<u32>() as f64;
+            println!("\t\t\t\t\tTotal:\t{:.2}", total_h/total_n);
+               
+        }
+
         Ok(Dataset::new(num_bins, dimens_length, bin_width, cfg.hist_min.clone(), cfg.hist_max.clone(), bias_pos, bias_fc, kT, histograms, cfg.cyclic))
     } else {
         bail!("Histogram has no datapoints.")
@@ -115,12 +132,9 @@ fn is_in_time_boundaries(time: f64, cfg: &Config) -> bool {
     false
 }
 
-// parse a time series 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))?;
-    let mut buf = BufReader::new(&f);
 
+// parse a time series file into a histogram
+fn read_window_file(window_file: &str, cfg: &Config) -> Result<(Histogram, usize)> {
     // total number of bins is the product of all dimensions length
     let total_bins = cfg.num_bins.iter().product();
     let mut hist = vec![0.0; total_bins];
@@ -130,11 +144,48 @@ fn read_window_file(window_file: &str, cfg: &Config) -> Result<Histogram> {
             (cfg.hist_max[idx] - cfg.hist_min[idx])/(cfg.num_bins[idx] as f64)
         }).collect();
 
+    let mut timeseries: Vec<Vec<f64>> = read_timeseries(window_file, cfg)?;
+    let timeseries_inital_length = timeseries[0].len();
+
+    if cfg.uncorr {
+        timeseries = uncorrelate(timeseries, cfg);
+    }
+
+    for i in 0..timeseries[0].len() {
+        let mut values: Vec<f64> = vec![f64::NAN; cfg.dimens+1];
+        for j in 0..values.len() {
+            values[j] = timeseries[j][i];
+        }
+
+        if is_in_hist_boundaries(&values[1..], cfg) && is_in_time_boundaries(values[0], cfg) {
+            let bin_indeces: Vec<usize> = (0..cfg.dimens).map(|dimen: usize| {
+                let val = values[dimen+1];
+                ((val - cfg.hist_min[dimen]) / bin_width[dimen]) as usize
+            }).collect();
+            let index = flat_index(&bin_indeces, &cfg.num_bins);
+            hist[index] += 1.0;
+        }
+    }
+
+    let num_points: f64 = hist.iter().sum();
+    Ok((Histogram::new(num_points as u32, hist), timeseries_inital_length))
+}
+
+// Read a multidimensional timeseries
+// The resulting vector contains one vector per dimension
+fn read_timeseries(window_file: &str, cfg: &Config) -> Result<Vec<Vec<f64>>> {
+    let f = File::open(window_file)
+        .chain_err(|| format!("Failed to open sample data file {}.", window_file))?;
+    let mut buf = BufReader::new(&f);
+
+    let mut timeseries = vec![Vec::new(); cfg.dimens+1];
+
     // 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.is_empty() {
             line.clear();
@@ -142,32 +193,55 @@ fn read_window_file(window_file: &str, cfg: &Config) -> Result<Histogram> {
         }
 
         {
-
             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+1];
-            for i in 0..values.len() {
-                values[i] = split[i].parse::<f64>()
-                    .chain_err(|| format!("Failed to parse line {} of window file {}.", linecount, window_file))?;
-            }
+            for i in 0..cfg.dimens+1 {
+                timeseries[i].push(split[i].parse::<f64>()
+                    .chain_err(|| format!("Failed to parse line {} of window file {}.", linecount, window_file))?
 
-            if is_in_hist_boundaries(&values[1..], cfg) && is_in_time_boundaries(values[0], cfg) {
-                let bin_indeces: Vec<usize> = (0..cfg.dimens).map(|dimen: usize| {
-                    let val = values[dimen+1];
-                    ((val - cfg.hist_min[dimen]) / bin_width[dimen]) as usize
-                }).collect();
-                let index = flat_index(&bin_indeces, &cfg.num_bins);
-                hist[index] += 1.0;
+                );
             }
         }
+        
         line.clear();
     }
+    Ok(timeseries)
+}
 
-    let num_points: f64 = hist.iter().sum();
-    Ok(Histogram::new(num_points as u32, hist))
+
+// calculates the inefficiency for every collective variable
+// filters the timeseries based on the highest inefficiency
+fn uncorrelate(timeseries: Vec<Vec<f64>>, cfg: &Config) -> Vec<Vec<f64>> {
+    // calculate inefficiencies and find the highest one
+    let gs: Vec<f64> = timeseries[1..].iter().map(|ts| statistical_ineff(ts)).collect();
+    let mut max_g = 1.0;
+    for g in gs {
+        if g > max_g {
+            max_g = g;
+        }
+    }
+
+    // round g up
+    let mut trunc_g = max_g.trunc() as usize;
+    if (trunc_g as f64 - max_g).abs() > 0.000_000_000_1 {
+        trunc_g += 1;
+    }
+
+    // filter correlated samples from timeseries
+    let prev_len = timeseries[0].len();
+    let timeseries = timeseries.into_iter().map(|ts| {
+        ts.into_iter().step_by(trunc_g).collect::<Vec<f64>>()
+    }).collect::<Vec<Vec<f64>>>();
+
+    let new_len = timeseries[0].len();
+    if cfg.verbose {
+        let tau = autocorrelation_time(max_g)* (timeseries[0][1]-timeseries[0][0]);
+        vprintln(format!("{:?}/{:?} samples are uncorrelated. {:?} samples removed from timeseries (tau={:.5})", new_len, prev_len, prev_len-new_len, tau), true);
+    }
+    timeseries
 }
 
 // Write WHAM calculation results to out_file.
@@ -214,7 +288,8 @@ mod tests {
             bootstrap: 0,
             bootstrap_seed: 1234,
             start: 0.0,
-            end: 1e+20
+            end: 1e+20,
+            uncorr: false,
         }
     }
 
@@ -222,8 +297,9 @@ mod tests {
     fn read_window_file() {
         let f = "example/1d_cyclic/COLVAR+0.0.xvg";
         let cfg = cfg();
-        let h = super::read_window_file(&f, &cfg).unwrap();
+        let (h, timeseries_inital_length) = super::read_window_file(&f, &cfg).unwrap();
         println!("{:?}", h);
+        assert_eq!(5000, timeseries_inital_length);
         assert_eq!(5000, h.num_points);
         assert_approx_eq!(0.0, h.bins[2]);
         assert_approx_eq!(11.0, h.bins[3]);
@@ -233,6 +309,26 @@ mod tests {
         assert_approx_eq!(0.0, h.bins[7]);
     }
 
+    #[test]
+    fn read_timeseries() {
+        let f = "example/1d_cyclic/COLVAR+0.0.xvg";
+        let cfg = cfg();
+        let ts = super::read_timeseries(&f, &cfg).unwrap();
+        let expected = [
+            -0.153_145,
+            -0.377_860,
+            0.010_992,
+            0.123_074,
+            0.108_291,
+            0.261_607,
+        ];
+        assert!(ts.len() == 2);
+        assert!(ts[0].len() == 5000);
+        println!("{:?}", ts);
+        for (actual, expected) in ts[1].iter().zip(expected.iter()) {
+            assert!((actual-expected).abs() < 0.001, format!("{:?} != {:?}", actual, expected));
+        }
+    }
 
     #[test]
     fn read_data() {

src/lib.rs

@@ -3,7 +3,6 @@
 #[macro_use]
 extern crate error_chain;
 extern crate rand;
-extern crate rgsl;
 extern crate rayon;
 #[cfg(test)]
 #[macro_use]
@@ -13,6 +12,8 @@ extern crate assert_approx_eq;
 pub mod io;
 pub mod histogram;
 pub mod error_analysis;
+pub mod correlation_analysis;
+pub mod statistics;
 
 use histogram::Dataset;
 use std::f64;
@@ -45,16 +46,17 @@ pub struct Config {
     pub bootstrap_seed: u64,
     pub start: f64,
     pub end: f64,
+    pub uncorr: bool,
 }
 
 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={:?}, seed={:?}",
+            cyclic={:?}, uncorr={:?}, 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)
+            self.cyclic, self.uncorr, self.bootstrap, self.bootstrap_seed)
     }
 }
 

src/main.rs

@@ -58,6 +58,8 @@ fn cli() -> Result<Config> {
         .chain_err(|| "Cannot parse start time.")?;
     let end: f64 = matches.value_of("end").unwrap_or("1e+20").parse()
         .chain_err(|| "Cannot parse end time.")?;
+
+    let uncorr: bool = matches.is_present("uncorr");
      
     if num_bins.len() != hist_max.len() || num_bins.len() != hist_max.len() {
         eprintln!("Input dimensions do not match (min: {}, max: {}, bins: {})",
@@ -69,7 +71,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, bootstrap_seed, start, end})
+		bootstrap, bootstrap_seed, start, end, uncorr})
 }
 
 fn main() {

src/statistics.rs

@@ -0,0 +1,59 @@
+pub fn mean(x: &[f64]) -> f64 {
+    x.iter().sum::<f64>() / x.len() as f64
+}
+
+pub fn autocov(x: &[f64]) -> f64 {
+    let x_mean = mean(x);
+
+    x.iter().map(|xi| {
+        (xi-x_mean).powi(2)
+    }).sum::<f64>() / x.len() as f64
+}
+
+pub fn sd(x: &[f64]) -> f64 {
+    let x_mean = mean(x);
+    let n = x.len() as f64;
+
+    let sum = x.iter().map(|xi| {
+        (xi-x_mean).powi(2)
+    }).sum::<f64>();
+    (1.0/(n-1.0) * sum).sqrt()
+}
+
+#[cfg(test)]
+mod tests {
+    use assert_approx_eq::assert_approx_eq;
+
+    // a sine wave
+    fn dataset() -> Vec<f64> {
+        (0..100).map(|i| (i as f64 / 100.0 * std::f64::consts::PI).sin())
+            .collect::<Vec<f64>>()
+    }
+
+    #[test]
+    fn mean() {
+        let ds = dataset();
+        let expected =0.6366;
+        
+        let m = super::mean(&ds);
+        assert_approx_eq!(m, expected, 0.0001);
+    }
+
+    #[test]
+    fn autocorr() {
+        let ds = dataset();
+        let expected = 0.094_782;
+        
+        let m = super::autocov(&ds);
+        assert_approx_eq!(m, expected, 0.000_001);
+    }
+
+    #[test]
+    fn sd() {
+        let ds = dataset();
+        let expected = 0.309_418;
+        
+        let m = super::sd(&ds);
+        assert_approx_eq!(m, expected, 0.000_001);
+    }
+}

tests/examples.rs

@@ -27,6 +27,26 @@ mod integration {
         assert_eq!(output_len, 0);
     }
 
+    #[test]
+    fn wham_1d_cyclic_uncorrelated() {
+        get_command()
+            .args(&["--bins", "100", "--max", "pi", "--min", "-pi", "-T", "300", "--cyclic", "--uncorr"])
+            .args(&["--seed", "1234"])
+            .args(&["-f", "example/1d_cyclic/metadata.dat"])
+            .args(&["-o", "/tmp/wham_test_1d_cyclic.out"])
+            .output()
+            .expect("failed to execute process");
+
+        assert!(fs::metadata("/tmp/wham_test_1d_cyclic.out").is_ok());
+        let output = Command::new("diff")
+            .arg("/tmp/wham_test_1d_cyclic.out")
+            .arg("example/1d_cyclic/wham_uncorrelated.out")
+            .output()
+            .expect("failed to run diff");
+        let output_len = String::from_utf8_lossy(&output.stdout).len();
+        assert_eq!(output_len, 0);
+    }
+
     #[test]
     #[ignore] // expensive
     fn wham_2d_cyclic() {