diff --git a/Cargo.lock b/Cargo.lock index 4ffbe8e..b760c8b 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -1,16 +1,5 @@ # This file is automatically @generated by Cargo. # It is not intended for manual editing. -[[package]] -name = "GSL" -version = "1.1.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "7830156ea389bcbbdc8f01bf140b609b892bf7cbd0ec6ccf9957ea2be6f25ad3" -dependencies = [ - "c_vec", - "libc", - "pkg-config", -] - [[package]] name = "addr2line" version = "0.13.0" @@ -78,12 +67,6 @@ version = "1.2.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "cf1de2fe8c75bc145a2f577add951f8134889b4795d47466a54a5c846d691693" -[[package]] -name = "c_vec" -version = "1.0.12" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "aa9e1d9f7d49e289f36f19effbf3d5a5e30163ecf9c7a3c9be94d5374dec5b9a" - [[package]] name = "cfg-if" version = "0.1.10" @@ -215,9 +198,9 @@ checksum = "e2abad23fbc42b3700f2f279844dc832adb2b2eb069b2df918f455c4e18cc646" [[package]] name = "libc" -version = "0.2.79" +version = "0.2.80" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "2448f6066e80e3bfc792e9c98bf705b4b0fc6e8ef5b43e5889aff0eaa9c58743" +checksum = "4d58d1b70b004888f764dfbf6a26a3b0342a1632d33968e4a179d8011c760614" [[package]] name = "memoffset" @@ -254,12 +237,6 @@ version = "0.21.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "37fd5004feb2ce328a52b0b3d01dbf4ffff72583493900ed15f22d4111c51693" -[[package]] -name = "pkg-config" -version = "0.3.19" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "3831453b3449ceb48b6d9c7ad7c96d5ea673e9b470a1dc578c2ce6521230884c" - [[package]] name = "ppv-lite86" version = "0.2.9" @@ -387,7 +364,6 @@ checksum = "cccddf32554fecc6acb585f82a32a72e28b48f8c4c1883ddfeeeaa96f7d8e519" name = "wham" version = "0.9.9" dependencies = [ - "GSL", "assert_approx_eq", "clap", "error-chain", diff --git a/Cargo.toml b/Cargo.toml index 637f06a..5ddb8cb 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -13,7 +13,6 @@ keywords = ["math", "statistics", "histogram", "bioinformatics", "molecular-dyna clap = {version="2.32.0", features=['yaml']} error-chain = "0.12.0" rand = "0.7.*" -GSL = "1.1" rayon = "1.0.3" [dev-dependencies] diff --git a/src/correlation_analysis.rs b/src/correlation_analysis.rs index a26bb25..6ff24d0 100644 --- a/src/correlation_analysis.rs +++ b/src/correlation_analysis.rs @@ -1,4 +1,4 @@ -use rgsl::statistics; +use super::statistics; // calculates the statistical inefficiency g of the given timeseries // the quantity g can be thought of: N/g is the number of uncorrelated @@ -28,9 +28,9 @@ 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, 1, timeseries.len()); + let mean = statistics::mean(timeseries); let d_mean = timeseries.iter().map(|x| x-mean).collect::>(); - let cov = statistics::covariance(timeseries, 1, timeseries, 1, n); + let cov = statistics::autocov(timeseries); let mut autocorr = Vec::new(); for t in 1..(n-1) { diff --git a/src/error_analysis.rs b/src/error_analysis.rs index 45219e8..0c2a3fd 100644 --- a/src/error_analysis.rs +++ b/src/error_analysis.rs @@ -2,7 +2,7 @@ use rand::prelude::*; use super::histogram::{Dataset}; use super::perform_wham; use super::{Config,calc_free_energy}; -use rgsl::statistics; +use super::statistics; // returns a set of num_windows continious weights by // a) generate num_windows-1 random variables and sort them @@ -48,7 +48,7 @@ pub fn run_bootstrap(cfg: &Config, ds: Dataset, num_runs: usize) -> (Vec,Ve let mut P_se = vec![0.0; ds.num_bins]; for bin in 0..ds.num_bins { let Ps = bootstrapped_Ps.iter().map(|window| window[bin]).collect::>(); - P_se[bin] = statistics::sd(&Ps, 1, num_runs)/(num_runs as f64).sqrt(); + P_se[bin] = statistics::sd(&Ps)/(num_runs as f64).sqrt(); } // SE of A @@ -60,7 +60,7 @@ pub fn run_bootstrap(cfg: &Config, ds: Dataset, num_runs: usize) -> (Vec,Ve let mut A_se = vec![0.0; ds.num_bins]; for bin in 0..ds.num_bins { let As = bootstrapped_As.iter().map(|window| window[bin]).collect::>(); - A_se[bin] = statistics::sd(&As, 1, num_runs)/(num_runs as f64).sqrt(); + A_se[bin] = statistics::sd(&As)/(num_runs as f64).sqrt(); } (P_se, A_se) diff --git a/src/lib.rs b/src/lib.rs index c0fbea3..a604ba4 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -3,7 +3,6 @@ #[macro_use] extern crate error_chain; extern crate rand; -extern crate rgsl; extern crate rayon; #[cfg(test)] #[macro_use] @@ -14,6 +13,7 @@ pub mod io; pub mod histogram; pub mod error_analysis; pub mod correlation_analysis; +pub mod statistics; use histogram::Dataset; use std::f64; diff --git a/src/statistics.rs b/src/statistics.rs new file mode 100644 index 0000000..e69c1bc --- /dev/null +++ b/src/statistics.rs @@ -0,0 +1,59 @@ +pub fn mean(x: &[f64]) -> f64 { + x.iter().sum::() / x.len() as f64 +} + +pub fn autocov(x: &[f64]) -> f64 { + let x_mean = mean(x); + + x.iter().map(|xi| { + (xi-x_mean).powi(2) + }).sum::() / x.len() as f64 +} + +pub fn sd(x: &[f64]) -> f64 { + let x_mean = mean(x); + let n = x.len() as f64; + + let sum = x.iter().map(|xi| { + (xi-x_mean).powi(2) + }).sum::(); + (1.0/(n-1.0) * sum).sqrt() +} + +#[cfg(test)] +mod tests { + use assert_approx_eq::assert_approx_eq; + + // a sine wave + fn dataset() -> Vec { + (0..100).map(|i| (i as f64 / 100.0 * std::f64::consts::PI).sin()) + .collect::>() + } + + #[test] + fn mean() { + let ds = dataset(); + let expected =0.6366; + + let m = super::mean(&ds); + assert_approx_eq!(m, expected, 0.0001); + } + + #[test] + fn autocorr() { + let ds = dataset(); + let expected = 0.094_782; + + let m = super::autocov(&ds); + assert_approx_eq!(m, expected, 0.000_001); + } + + #[test] + fn sd() { + let ds = dataset(); + let expected = 0.309_418; + + let m = super::sd(&ds); + assert_approx_eq!(m, expected, 0.000_001); + } +} \ No newline at end of file