diff --git a/README.md b/README.md index 6553ccd..3793af9 100644 --- a/README.md +++ b/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 @@ -67,6 +68,8 @@ 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,7 +88,6 @@ OPTIONS: -T, --temperature WHAM temperature in Kelvin. -t, --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): @@ -134,6 +136,18 @@ To perform bayesian bootstrapping in WHAM, use the ```-bt ``` 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. +Autocorrelation analysis +--- +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. + +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* + + Examples --- The example folder contains input and output files for two simple test systems: @@ -144,7 +158,6 @@ The example folder contains input and output files for two simple test systems: TODO --- -- Autocorrelation - Replica exchange License & Citing diff --git a/example/1d_cyclic/wham_uncorrelated.out b/example/1d_cyclic/wham_uncorrelated.out new file mode 100644 index 0000000..8ca4fdc --- /dev/null +++ b/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 diff --git a/src/cli.yml b/src/cli.yml index 633adfd..7c3f714 100644 --- a/src/cli.yml +++ b/src/cli.yml @@ -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 diff --git a/src/correlation_analysis.rs b/src/correlation_analysis.rs new file mode 100644 index 0000000..a26bb25 --- /dev/null +++ b/src/correlation_analysis.rs @@ -0,0 +1,102 @@ +use rgsl::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 { + let n = timeseries.len(); + let mean = statistics::mean(timeseries, 1, timeseries.len()); + let d_mean = timeseries.iter().map(|x| x-mean).collect::>(); + let cov = statistics::covariance(timeseries, 1, timeseries, 1, n); + + 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::() / (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 { + let mut timeseries: Vec = 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::>()[1].parse::().unwrap(); + timeseries.push(val); + } + timeseries + + } + + #[test] + fn autocorrelation() { + let timeseries = read_timeseries("example/1d_cyclic/COLVAR-2.5.xvg"); + let autocorr = super::autocorrelation(×eries); + let expected = [ + 0.6919008655979143, 0.5331719399671355, + 0.20472620956463589, -0.002850876458920514, + -0.13842850077938146, -0.2652923552973232, + -0.31427198272235385, -0.2617505151557693, + -0.20594864730290338, -0.13310019091811812, + -0.1887568901193426, -0.1944936625424933, + -0.1963599673189461, -0.11391587833838003]; + 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(×eries); + 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(×eries); + let tau = super::autocorrelation_time(g); + println!("{:?}", tau); + assert!((tau - 1.430).abs() < 0.001) + } + +} \ No newline at end of file diff --git a/src/error_analysis.rs b/src/error_analysis.rs index ce8bb9e..45219e8 100644 --- a/src/error_analysis.rs +++ b/src/error_analysis.rs @@ -66,7 +66,7 @@ pub fn run_bootstrap(cfg: &Config, ds: Dataset, num_runs: usize) -> (Vec,Ve (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); diff --git a/src/io.rs b/src/io.rs index 4eba21f..6c2b904 100644 --- a/src/io.rs +++ b/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}; @@ -115,12 +116,80 @@ 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 { - let f = File::open(window_file) +// Read a multidimensional timeseries +// The resulting vector contains one vector per dimension +fn read_timeseries(window_file: &str, cfg: &Config) -> Result>> { + 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(); + continue; + } + + { + 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)); + } + + for i in 0..cfg.dimens+1 { + timeseries[i].push(split[i].parse::() + .chain_err(|| format!("Failed to parse line {} of window file {}.", linecount, window_file))? + + ); + } + } + + line.clear(); + } + Ok(timeseries) +} + +// calculates the inefficiency for every collective variable +// filters the timeseries based on the highest inefficiency +fn uncorrelate(timeseries: Vec>, cfg: &Config) -> Vec> { + // calculate inefficiencies and find the highest one + let gs: Vec = 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::>() + }).collect::>>(); + + 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 +} + +// parse a time series file into a histogram +fn read_window_file(window_file: &str, cfg: &Config) -> Result { // 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,40 +199,26 @@ fn read_window_file(window_file: &str, cfg: &Config) -> Result { (cfg.hist_max[idx] - cfg.hist_min[idx])/(cfg.num_bins[idx] as f64) }).collect(); - // 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(); - continue; + let mut timeseries: Vec> = read_timeseries(window_file, cfg)?; + + if cfg.uncorr { + timeseries = uncorrelate(timeseries, cfg); + } + + for i in 0..timeseries[0].len() { + let mut values: Vec = vec![f64::NAN; cfg.dimens+1]; + for j in 0..values.len() { + values[j] = timeseries[j][i]; } - { - - 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 = vec![f64::NAN; cfg.dimens+1]; - for i in 0..values.len() { - values[i] = split[i].parse::() - .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 = (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; - } + if is_in_hist_boundaries(&values[1..], cfg) && is_in_time_boundaries(values[0], cfg) { + let bin_indeces: Vec = (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(); } let num_points: f64 = hist.iter().sum(); @@ -214,7 +269,8 @@ mod tests { bootstrap: 0, bootstrap_seed: 1234, start: 0.0, - end: 1e+20 + end: 1e+20, + uncorr: false, } } @@ -233,6 +289,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() { diff --git a/src/lib.rs b/src/lib.rs index 90ec06f..c0fbea3 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -13,6 +13,7 @@ extern crate assert_approx_eq; pub mod io; pub mod histogram; pub mod error_analysis; +pub mod correlation_analysis; 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) } } diff --git a/src/main.rs b/src/main.rs index 6764e7a..6aadf80 100644 --- a/src/main.rs +++ b/src/main.rs @@ -58,6 +58,8 @@ fn cli() -> Result { .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 { 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() { diff --git a/tests/examples.rs b/tests/examples.rs index e2643f2..783b453 100644 --- a/tests/examples.rs +++ b/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() {