diff --git a/src/correlation_analysis.rs b/src/correlation_analysis.rs index fea772c..f7cff7d 100644 --- a/src/correlation_analysis.rs +++ b/src/correlation_analysis.rs @@ -8,11 +8,14 @@ use super::statistics; // 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 mean = statistics::mean(timeseries); + let d_mean = timeseries.iter().map(|x| x-mean).collect::>(); + let cov = statistics::autocov(timeseries); + let mut g = 1.0; for t in 1..(n-1) { - let c = autocorr[t-1]; + let tmp = d_mean[0..n-t].iter().zip(d_mean[t..n].iter()).map(|(x,y)| x*y); + let c = tmp.map(|x| x+x).sum::() / (2.0 * (n as f64 - t as f64)*cov); if c <= 0.0 { break; } @@ -25,22 +28,6 @@ pub fn statistical_ineff(timeseries: &[f64]) -> f64 { } } -// calculates the autocorrelation of a simeseries -fn autocorrelation(timeseries: &[f64]) -> Vec { - let n = timeseries.len(); - let mean = statistics::mean(timeseries); - let d_mean = timeseries.iter().map(|x| x-mean).collect::>(); - let cov = statistics::autocov(timeseries); - - 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::() / (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 { @@ -65,23 +52,6 @@ mod tests { } - #[test] - fn autocorrelation() { - let timeseries = read_timeseries("example/1d_cyclic/COLVAR-2.5.xvg"); - let autocorr = super::autocorrelation(×eries); - let expected = [ - 0.691_900_865_597_914_3, 0.533_171_939_967_135_5, - 0.204_726_209_564_635_89, -0.002_850_876_458_920_514, - -0.138_428_500_779_381_46, -0.265_292_355_297_323_2, - -0.314_271_982_722_353_85, -0.261_750_515_155_769_3, - -0.205_948_647_302_903_38, -0.133_100_190_918_118_12, - -0.188_756_890_119_342_6, -0.194_493_662_542_493_3, - -0.196_359_967_318_946_1, -0.113_915_878_338_380_03]; - 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"); @@ -99,4 +69,4 @@ mod tests { assert!((tau - 1.430).abs() < 0.001) } -} \ No newline at end of file +}