uncorrelate data

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/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<f64> {
 
 // 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
 }
 

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(

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<Vec<Vec<f64>>> {
     Ok(timeseries)
 }
 
+// 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
+}
+
 // parse a time series file into a histogram
 fn read_window_file(window_file: &str, cfg: &Config) -> Result<Histogram> {
     // 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<Histogram> {
             (cfg.hist_max[idx] - cfg.hist_min[idx])/(cfg.num_bins[idx] as f64)
         }).collect();
 
-    let timeseries: Vec<Vec<f64>> = read_timeseries(window_file, cfg)?;
+    let mut timeseries: Vec<Vec<f64>> = 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<f64> = vec![f64::NAN; cfg.dimens+1];

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() {