remove gsl

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<f64> {
     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::<Vec<f64>>();
-    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) {

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<f64>,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::<Vec<f64>>();
-        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<f64>,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::<Vec<f64>>();
-        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)

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;

src/statistics.rs

@@ -0,0 +1,59 @@
+pub fn mean(x: &[f64]) -> f64 {
+    x.iter().sum::<f64>() / 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::<f64>() / 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::<f64>();
+    (1.0/(n-1.0) * sum).sqrt()
+}
+
+#[cfg(test)]
+mod tests {
+    use assert_approx_eq::assert_approx_eq;
+
+    // a sine wave
+    fn dataset() -> Vec<f64> {
+        (0..100).map(|i| (i as f64 / 100.0 * std::f64::consts::PI).sin())
+            .collect::<Vec<f64>>()
+    }
+
+    #[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);
+    }
+}