bauer/WHAM
1
0
Fork 0
mirror of https://github.com/dnlbauer/WHAM.git synced 2026-09-10 14:15:31 +00:00
Code Issues Packages Projects Releases Wiki Activity

inlined autocorrelation calculation for better performance

Browse Source
This commit is contained in:
daniel
2020-10-26 11:24:48 +01:00
parent cc4448f1d8
commit 034586d08e
1 changed files with 7 additions and 37 deletions
Download Patch File Download Diff File Expand all files Collapse all files

44
src/correlation_analysis.rs
View File

@@ -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::<Vec<f64>>();
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::<f64>() / (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<f64> {
let n = timeseries.len();
let mean = statistics::mean(timeseries);
let d_mean = timeseries.iter().map(|x| x-mean).collect::<Vec<f64>>();
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::<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 {
@@ -65,23 +52,6 @@ mod tests {
}
#[test]
fn autocorrelation() {
let timeseries = read_timeseries("example/1d_cyclic/COLVAR-2.5.xvg");
let autocorr = super::autocorrelation(&timeseries);
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)
}
}
}