Autocorrelation (#2)

* 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>

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 

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<f64> {
+    let n = timeseries.len();
+    let mean = statistics::mean(timeseries, 1, timeseries.len());
+    let d_mean = timeseries.iter().map(|x| x-mean).collect::<Vec<f64>>();
+    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::<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.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(&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

@@ -66,7 +66,7 @@ pub fn run_bootstrap(cfg: &Config, ds: Dataset, num_runs: usize) -> (Vec<f64>,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);

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<Histogram> {
-	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<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();
+            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::<f64>()
+                    .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<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
     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<Histogram> {
             (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<Vec<f64>> = read_timeseries(window_file, cfg)?;
+    
+    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];
         }
 
-        {
-
-            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))?;
-            }
-
-            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;
-            }
+        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();
     }
 
     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() {

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)
     }
 }
 

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