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/correlation_analysis.rs b/src/correlation_analysis.rs index 054c838..a26bb25 100644 --- a/src/correlation_analysis.rs +++ b/src/correlation_analysis.rs @@ -6,7 +6,7 @@ use rgsl::statistics; // 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" -fn statistical_ineff(timeseries: &[f64]) -> f64 { +pub fn statistical_ineff(timeseries: &[f64]) -> f64 { let n = timeseries.len(); let autocorr = autocorrelation(timeseries); @@ -43,7 +43,7 @@ fn autocorrelation(timeseries: &[f64]) -> Vec { // The autocorrelation time of a timeseries can be deduced from the // `statistical_ineff` by (g-1)/2.0 -fn autocorrelation_time(g: f64) -> f64 { +pub fn autocorrelation_time(g: f64) -> f64 { (g - 1.0) / 2.0 } diff --git a/src/error_analysis.rs b/src/error_analysis.rs index b64b36e..45219e8 100644 --- a/src/error_analysis.rs +++ b/src/error_analysis.rs @@ -71,7 +71,6 @@ mod tests { use super::*; use super::super::k_B; use super::super::histogram::Histogram; - use rand::prelude::*; fn build_hist() -> Histogram { Histogram::new( diff --git a/src/io.rs b/src/io.rs index 843cbb2..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}; @@ -155,6 +156,38 @@ fn read_timeseries(window_file: &str, cfg: &Config) -> Result>> { 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 @@ -166,9 +199,11 @@ 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(); - let timeseries: Vec> = read_timeseries(window_file, cfg)?; + let mut timeseries: Vec> = read_timeseries(window_file, cfg)?; - // TODO decorrelate + if cfg.uncorr { + timeseries = uncorrelate(timeseries, cfg); + } for i in 0..timeseries[0].len() { let mut values: Vec = vec![f64::NAN; cfg.dimens+1]; 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() {