diff --git a/src/correlation_analysis.rs b/src/correlation_analysis.rs new file mode 100644 index 0000000..7928118 --- /dev/null +++ b/src/correlation_analysis.rs @@ -0,0 +1,75 @@ +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" +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::>(); + let cov = statistics::covariance(timeseries, 1, timeseries, 1, n); + + let mut g = 1.0; + for t in 1..(n-1) { + + // normalized autocorr C(t) = ( - ^2) / ( - ^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::() / (2.0 * (n as f64-t as f64)*cov); + + if c <= 0.0 { // terminate at first 0 (autocorr gets noisy from here) + break; + } + g = g + (2.0*c*(1.0-t as f64/n as f64)) + } + if g < 1.0 { + 1.0 + } else { + g + } +} + +// The autocorrelation time of a timeseries can be deduced from the +// `statistical_ineff` by (g-1)/2.0 +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 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..b64b36e 100644 --- a/src/error_analysis.rs +++ b/src/error_analysis.rs @@ -66,11 +66,12 @@ 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; use super::super::histogram::Histogram; + use rand::prelude::*; fn build_hist() -> Histogram { Histogram::new( @@ -99,8 +100,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 +111,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/lib.rs b/src/lib.rs index 90ec06f..497a172 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;