From 07fc8344fe170669b680d6f3b2c09fa6feaa58e6 Mon Sep 17 00:00:00 2001
From: Daniel Bauer <bauer@cbs.tu-darmstadt.de>
Date: Sun, 25 Oct 2020 19:42:33 +0100
Subject: [PATCH] read timeseries into vector

---
 src/io.rs | 111 +++++++++++++++++++++++++++++++++++++-----------------
 1 file changed, 76 insertions(+), 35 deletions(-)

diff --git a/src/io.rs b/src/io.rs
index 4eba21f..843cbb2 100644
--- a/src/io.rs
+++ b/src/io.rs
@@ -115,12 +115,48 @@ fn is_in_time_boundaries(time: f64, cfg: &Config) -> bool {
     false
 }
 
-// parse a time series file into a histogram
-fn read_window_file(window_file: &str, cfg: &Config) -> Result<Histogram> {
-	let f = File::open(window_file)
+// Read a multidimensional timeseries
+// The resulting vector contains one vector per dimension
+fn read_timeseries(window_file: &str, cfg: &Config) -> Result<Vec<Vec<f64>>> {
+    let f = File::open(window_file)
         .chain_err(|| format!("Failed to open sample data file {}.", window_file))?;
     let mut buf = BufReader::new(&f);
 
+    let mut timeseries = vec![Vec::new(); cfg.dimens+1];
+
+    // read and parse each timeseries line
+    let mut line = String::new();
+    let mut linecount = 0;
+    while buf.read_line(&mut line).chain_err(|| "Failed to read line")? > 0 {
+        linecount += 1;
+
+        // skip comments and empty lines
+        if line.starts_with('#') || line.starts_with('@') || line.is_empty() {
+            line.clear();
+            continue;
+        }
+
+        {
+            let split: Vec<&str> = line.split_whitespace().collect();
+            if split.len() < cfg.dimens+1 {
+                bail!(format!("Wrong number of columns in line {} of window file {}. Empty Line?.", linecount, window_file));
+            }
+
+            for i in 0..cfg.dimens+1 {
+                timeseries[i].push(split[i].parse::<f64>()
+                    .chain_err(|| format!("Failed to parse line {} of window file {}.", linecount, window_file))?
+
+                );
+            }
+        }
+        
+        line.clear();
+    }
+    Ok(timeseries)
+}
+
+// parse a time series file into a histogram
+fn read_window_file(window_file: &str, cfg: &Config) -> Result<Histogram> {
     // total number of bins is the product of all dimensions length
     let total_bins = cfg.num_bins.iter().product();
     let mut hist = vec![0.0; total_bins];
@@ -130,40 +166,24 @@ fn read_window_file(window_file: &str, cfg: &Config) -> Result<Histogram> {
             (cfg.hist_max[idx] - cfg.hist_min[idx])/(cfg.num_bins[idx] as f64)
         }).collect();
 
-    // read and parse each timeseries line
-    let mut line = String::new();
-    let mut linecount = 0;
-    while buf.read_line(&mut line).chain_err(|| "Failed to read line")? > 0 {
-        linecount += 1;
-        // skip comments and empty lines
-        if line.starts_with('#') || line.starts_with('@') || line.is_empty() {
-            line.clear();
-            continue;
+    let timeseries: Vec<Vec<f64>> = read_timeseries(window_file, cfg)?;
+    
+    // TODO decorrelate
+
+    for i in 0..timeseries[0].len() {
+        let mut values: Vec<f64> = vec![f64::NAN; cfg.dimens+1];
+        for j in 0..values.len() {
+            values[j] = timeseries[j][i];
         }
 
-        {
-
-            let split: Vec<&str> = line.split_whitespace().collect();
-            if split.len() < cfg.dimens+1 {
-                bail!(format!("Wrong number of columns in line {} of window file {}. Empty Line?.", linecount, window_file));
-            }
-
-            let mut values: Vec<f64> = vec![f64::NAN; cfg.dimens+1];
-            for i in 0..values.len() {
-                values[i] = split[i].parse::<f64>()
-                    .chain_err(|| format!("Failed to parse line {} of window file {}.", linecount, window_file))?;
-            }
-
-            if is_in_hist_boundaries(&values[1..], cfg) && is_in_time_boundaries(values[0], cfg) {
-                let bin_indeces: Vec<usize> = (0..cfg.dimens).map(|dimen: usize| {
-                    let val = values[dimen+1];
-                    ((val - cfg.hist_min[dimen]) / bin_width[dimen]) as usize
-                }).collect();
-                let index = flat_index(&bin_indeces, &cfg.num_bins);
-                hist[index] += 1.0;
-            }
+        if is_in_hist_boundaries(&values[1..], cfg) && is_in_time_boundaries(values[0], cfg) {
+            let bin_indeces: Vec<usize> = (0..cfg.dimens).map(|dimen: usize| {
+                let val = values[dimen+1];
+                ((val - cfg.hist_min[dimen]) / bin_width[dimen]) as usize
+            }).collect();
+            let index = flat_index(&bin_indeces, &cfg.num_bins);
+            hist[index] += 1.0;
         }
-        line.clear();
     }
 
     let num_points: f64 = hist.iter().sum();
@@ -214,7 +234,8 @@ mod tests {
             bootstrap: 0,
             bootstrap_seed: 1234,
             start: 0.0,
-            end: 1e+20
+            end: 1e+20,
+            uncorr: false,
         }
     }
 
@@ -233,6 +254,26 @@ mod tests {
         assert_approx_eq!(0.0, h.bins[7]);
     }
 
+    #[test]
+    fn read_timeseries() {
+        let f = "example/1d_cyclic/COLVAR+0.0.xvg";
+        let cfg = cfg();
+        let ts = super::read_timeseries(&f, &cfg).unwrap();
+        let expected = [
+            -0.153_145,
+            -0.377_860,
+            0.010_992,
+            0.123_074,
+            0.108_291,
+            0.261_607,
+        ];
+        assert!(ts.len() == 2);
+        assert!(ts[0].len() == 5000);
+        println!("{:?}", ts);
+        for (actual, expected) in ts[1].iter().zip(expected.iter()) {
+            assert!((actual-expected).abs() < 0.001, format!("{:?} != {:?}", actual, expected));
+        }
+    }
 
     #[test]
     fn read_data() {