uncorr flag

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

@@ -8,18 +8,12 @@ use rgsl::statistics;
 // Method for the Analysis of Simulated and Parallel Tempering Simulations, JCTC"
 fn statistical_ineff(timeseries: &[f64]) -> 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 autocorr = autocorrelation(timeseries);
 
     let mut g = 1.0;
     for t in 1..(n-1) {
-
-        // normalized autocorr C(t) = (<x_n*x_(n+t)> - <x_n>^2) / (<x_n^2> - <x_n>^2)
-        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);
-
-        if c <= 0.0 {  // terminate at first 0 (autocorr gets noisy from here)
+        let c = autocorr[t-1];
+        if c <= 0.0 {
             break;
         }
         g = g + (2.0*c*(1.0-t as f64/n as f64))
@@ -31,6 +25,22 @@ fn statistical_ineff(timeseries: &[f64]) -> f64 {
     }
 }
 
+// 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 
 fn autocorrelation_time(g: f64) -> f64 {
@@ -55,6 +65,23 @@ mod tests {
 
     }
 
+    #[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");

src/lib.rs

@@ -46,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() {