version up

commit eaebf0dcbb259decbc0d8f5bbffa62244303f6c7
Author: Daniel Bauer <bauer@cbs.tu-darmstadt.de>
Date:   Sat Jul 17 14:12:30 2021 +0200

    error for empty timeseries

commit 1fe5383c5079ca6c0ad102d4e2cb32d5b1227e80
Author: Daniel Bauer <bauer@cbs.tu-darmstadt.de>
Date:   Sat Jul 17 14:09:13 2021 +0200

    refractored convdt slices calculation

commit 0b1fb7fb6a72edc50b013b59623551b2ccab6913
Author: Daniel Bauer <bauer@cbs.tu-darmstadt.de>
Date:   Sat Jul 17 13:27:59 2021 +0200

    fix histogram building without convdt

commit f9882ca4cece57451cd2971d0993a2d3621b0cdf
Author: Daniel Bauer <bauer@cbs.tu-darmstadt.de>
Date:   Sat Jul 17 12:53:06 2021 +0200

    fix tests not compiling

commit e6550e20bde3824f8199432b08777be41fc79fa8
Author: Daniel Bauer <bauer@cbs.tu-darmstadt.de>
Date:   Sat Jul 17 12:43:36 2021 +0200

    refractoring

commit 18c77a9b6694d491ef23cebb3918ccda8fcc44a5
Author: Daniel Bauer <bauer@cbs.tu-darmstadt.de>
Date:   Fri Jul 16 08:16:34 2021 +0200

    run and output for multiple datasets

commit 069f318f72207c416387007435eeff9867494e7a
Author: Daniel Bauer <bauer@cbs.tu-darmstadt.de>
Date:   Fri Jul 16 07:58:01 2021 +0200

    cleanup io.rs

commit 10efa428a0c6bf490c9f2d7a4c1df185de402a11
Author: Daniel Bauer <bauer@cbs.tu-darmstadt.de>
Date:   Thu Jul 15 19:27:22 2021 +0200

    parse multiple datasets with convdt

Cargo.lock

@@ -362,7 +362,7 @@ checksum = "cccddf32554fecc6acb585f82a32a72e28b48f8c4c1883ddfeeeaa96f7d8e519"
 
 [[package]]
 name = "wham"
-version = "1.0.0"
+version = "1.1.0"
 dependencies = [
  "assert_approx_eq",
  "clap",

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "wham"
-version = "1.0.0"
+version = "1.1.0"
 authors = ["Daniel Bauer <bauer@cbs.tu-darmstadt.de>"]
 description = "An implementation of the weighted histogram analysis method"
 license = "GPL-3.0"
@@ -8,6 +8,9 @@ repository = "https://github.com/danijoo/WHAM"
 readme = "README.md"
 categories = ["science", "command-line-utilities", "algorithms"]
 keywords = ["math", "statistics", "histogram", "bioinformatics", "molecular-dynamics"]
+exclude = [
+    "example/*"
+]
 
 [dependencies]
 clap = {version="2.32.0", features=['yaml']}

README.md

@@ -20,12 +20,6 @@ Features
 
 Installation
 ---
-WHAM requires the GSL library to be installed: 
-```bash
-# on debian/ubuntu:
-sudo apt-get install libgsl0-dev
-```
-
 Installation from source via cargo:
 ```bash
 # cargo installation
@@ -40,8 +34,8 @@ wham has a convenient command line interface. You can see all options with
 ```wham -h```:
 
 ```
-wham 1.0.0
+wham 1.1.0
-D. Bauer <bauer@bio.tu-darmstadt.de>
+D. Bauer <bauer@cbs.tu-darmstadt.de>
 wham is a fast implementation of the weighted histogram analysis method (WHAM) written in Rust. It currently supports
 potential of mean force (PMF) calculations in multiple dimensions at constant temperature.
 
@@ -62,7 +56,7 @@ The first column will be ignored and is followed by N reaction coordinates x.
 Shipped under the GPLv3 license.
 
 USAGE:
-   wham [FLAGS] [OPTIONS] --bins <BINS> --max <HIST_MAX> --file <METADATA> --min <HIST_MIN> --temperature <temperature>
+    wham [FLAGS] [OPTIONS] --bins <BINS> --max <HIST_MAX> --file <METADATA> --min <HIST_MIN> --temperature <temperature>
 
 FLAGS:
     -c, --cyclic     For periodic reaction coordinates. If this is set, the first and last coordinate bin in each
@@ -78,6 +72,10 @@ OPTIONS:
         --bt <bootstrap>               Number of bayesian bootstrapping runs for error analysis by assigning random
                                        weights (defaults to 0).
         --seed <bootstrap_seed>        Random seed for bootstrapping runs.
+        --convdt <convdt>              Performs WHAM for slices with the given delta in time and returns an output file
+                                       for each slice. THis is useful to check the result for convergence. Example: with
+                                       --convdt 100 and a timeseries ranging from 0-300, free energy surfaces for slices
+                                       0-100, 0-200 and 0-300 will be given returned.
         --end <end>                    Skip rows in timeseries with an index larger than this value (defaults to 1e+20)
     -i, --iterations <ITERATIONS>      Stop WHAM after this many iterations without convergence (defaults to 100,000).
         --max <HIST_MAX>               Histogram maxima (comma separated). Also accepts "pi".
@@ -158,6 +156,7 @@ tempering simulations, JCTC 3(1):26-41*
 
 TODO
 ---
+- Option to output histograms
 - Replica exchange
 
 License & Citing

src/cli.yml

@@ -1,5 +1,5 @@
 name: wham
-version: "1.0.0"
+version: "1.1.0"
 author: D. Bauer <bauer@cbs.tu-darmstadt.de>
 about: |
   wham is a fast implementation of the weighted histogram analysis method (WHAM) written in Rust. It currently supports potential of mean force (PMF) calculations in multiple dimensions at constant temperature.
@@ -106,4 +106,9 @@ args:
         long: uncorr
         help: Estimates statistical inefficiency of each timeseries via autocorrelation and removes correlated samples (default is off).
         takes_value: false
-        required: false 
+        required: false
+    - convdt:
+        long: convdt
+        help: "Performs WHAM for slices with the given delta in time and returns an output file for each slice. THis is useful to check the result for convergence. Example: with --convdt 100 and a timeseries ranging from 0-300, free energy surfaces for slices 0-100, 0-200 and 0-300 will be given returned."
+        takes_value: true
+        required: false

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::<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,29 +52,18 @@ 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");
         let g = super::statistical_ineff(&timeseries);
         println!("{:?}", g);
-        assert!((g - 3.859).abs() < 0.001)
+        assert!((g - 3.859).abs() < 0.001);
+
+        // a "random" timeseries with g < 1.0
+        let timeseries = [1_f64, 4_f64, 921_f64, 121213_f64, 23192_f64,
+            8913_f64, 1232_f64, 2_f64, 151_f64, 123091_f64];
+        let g = super::statistical_ineff(&timeseries);
+        assert_approx_eq!(g, 1.0);
     }
 
     #[test]
@@ -99,4 +75,4 @@ mod tests {
         assert!((tau - 1.430).abs() < 0.001)
     }
 
-}
+}

src/io.rs

@@ -2,6 +2,7 @@ use super::histogram::Dataset;
 use super::histogram::Histogram;
 use super::Config;
 use super::correlation_analysis::{statistical_ineff, autocorrelation_time};
+use std::fs::OpenOptions;
 use std::fs::File;
 use std::io::prelude::*;
 use std::io::{BufReader,BufWriter};
@@ -26,23 +27,26 @@ pub fn vprintln(s: String, verbose: bool) {
 }
 
 // Read input data into a histogram set by iterating over input files
-// given in the metadata file
+// given in the metadata file. This generates at least one Dataset,
-pub fn read_data(cfg: &Config) -> Result<Dataset> {
+// or multiple Datasets if convdt is set in the config
+pub fn read_data(cfg: &Config) -> Result<Vec<Dataset>> {
 	let mut bias_pos: Vec<f64> = Vec::new();
 	let mut bias_fc: Vec<f64> = Vec::new();
-    let mut histograms: Vec<Histogram> = Vec::new();
+    let mut histograms: Vec<Vec<Histogram>> = Vec::new();
     let mut timeseries_lengths: Vec<usize> = Vec::new();
+    let mut paths = Vec::new();
 
 	let kT = cfg.temperature * k_B;
     let bin_width: Vec<f64> = (0..cfg.dimens).map(|idx| {
             (cfg.hist_max[idx] - cfg.hist_min[idx])/(cfg.num_bins[idx] as f64)
         }).collect();
-    let num_bins = cfg.num_bins.iter().product();
+    let num_bins: usize = cfg.num_bins.iter().product();
     let dimens_length = cfg.num_bins.clone();
 
     let f = File::open(&cfg.metadata_file).chain_err(|| "Failed to open metadata file")?;
     let buf = BufReader::new(&f);
 
+
     // read each metadata file line and parse it
     for (line_num,l) in buf.lines().enumerate() {
     	let line = l.chain_err(|| "Failed to read line")?;
@@ -57,18 +61,6 @@ pub fn read_data(cfg: &Config) -> Result<Dataset> {
             bail!(format!("Wrong number of columns in line {} of metadata file. Empty Line?", line_num+1));
         }
 
-        // parse histogram data
-        let path = get_relative_path(&cfg.metadata_file, split[0]);
-        let (h, timeseries_inital_length) = read_window_file(&path, cfg)
-            .chain_err(|| format!("Failed to parse process data file {}", &path))?;
-        if h.num_points == 0 {
-            bail!(format!("No data points in histogram boundaries: {}", &path))
-        }
-        histograms.push(h);
-        timeseries_lengths.push(timeseries_inital_length);
-        vprintln(format!("{}, {} data points added.", &path,
-                         histograms.last().unwrap().num_points), cfg.verbose);
-
         // parse bias force constants and positions
         for val in split.iter().skip(1).take(cfg.dimens) {
             let pos = val.parse()
@@ -80,12 +72,82 @@ pub fn read_data(cfg: &Config) -> Result<Dataset> {
                 .chain_err(|| format!("Failed to read bias fc in line {} of metadata file", line_num+1))?;
             bias_fc.push(fc);
         }
+
+        // parse histogram data
+        let path = get_relative_path(&cfg.metadata_file, split[0]);
+        paths.push(path.clone());
+        let (timeseries, timeseries_initial_lengths) = read_window_file(&path, cfg)
+            .chain_err(|| format!("Failed to read time series from {}", &path))?;
+        timeseries_lengths.push(timeseries_initial_lengths);
+
+
+        // for each timeseries, histograms are build for slices according to
+        // start..convdt, start..2*convdt, ...
+        histograms.push(Vec::new());
+        let h_idx = histograms.len()-1;
+        let convdt_stops = get_convdt_boundaries(&timeseries[0], &cfg);
+        for (idx, interval) in convdt_stops.iter().enumerate() {
+            // build histogram for slice start.._stop
+            let (start, stop) = interval;
+            let timeseries_mask: Vec<bool> = (0..timeseries[0].len()).map(|i| {
+                is_in_time_boundaries(timeseries[0][i], *start, *stop)
+            }).collect();
+            let hist = build_histogram_from_timeseries(&timeseries, &timeseries_mask, cfg);
+            histograms[h_idx].push(hist);
+
+            if (cfg.convdt == 0.00) || idx+1 == convdt_stops.len() {
+                vprintln(format!("{}, {} data points added.", &path,
+                    histograms[h_idx].last().unwrap().num_points), cfg.verbose);
+                break
+            }
+        }
     }
-    
+
-    if !histograms.is_empty() {
+    // Histograms are stored as timeseries x convdt right now,
+    // but we need convdt x timeseries to create Datasets
+    // this transposes the data
+    let num_datasets: usize = histograms.iter().map(|h| h.len()).max().unwrap();
+    let datasets: Vec<Dataset> = (0..num_datasets).map(|idx| {
+        let mut dataset_histograms: Vec<Histogram> = Vec::with_capacity(histograms.len());
+        for (hs, path) in histograms.iter().zip(&paths) {
+            if hs.len() > idx {
+                dataset_histograms.push(hs[idx].clone())
+            } else {
+                let warning = format!("No data points in histogram boundaries: {}", &path);
+                if idx+1 == num_datasets {
+                    bail!(warning);
+                } else {
+                    eprintln!("{}", warning);
+                }
+            }
+        }
+
+        Ok(Dataset::new(num_bins, dimens_length.clone(), bin_width.clone(),
+            cfg.hist_min.clone(), cfg.hist_max.clone(), bias_pos.clone(),
+            bias_fc.clone(), kT, dataset_histograms, cfg.cyclic))
+    }).collect::<Result<Vec<Dataset>>>().chain_err(|| "Failed to create datasets.")?;
+
+    if datasets.is_empty() {
+        bail!("No datasets created.")
+    } else if datasets[0].histograms.is_empty() {
+        bail!("Dataset has no associated data points.")
+    } else {
+        if datasets.len() > 1 {
+            println!("Datasets:");
+            println!("Dataset\t\tTime interval\t\tWindows\t\tN_total");
+            for (idx, dataset) in datasets.iter().enumerate() {
+                let n: u32 = dataset.histograms.iter().map(|h| h.num_points).sum();
+                let mut stop = cfg.start+cfg.convdt*(idx+1) as f64;
+                if stop > cfg.end {
+                    stop = cfg.end;
+                }
+                println!("{:?}\t\t{:?}-{:?}\t\t{:?}\t\t{:?}", idx+1, cfg.start, stop, dataset.histograms.len(), n);
+            }
+        }
+
+        let histograms = &datasets.last().unwrap().histograms;
         if cfg.uncorr {
-            println!("Timeseries Correlation");
+            println!("Timeseries Correlation:");
-            println!();
             println!("Window\t\tN\t\tN_uncorr\tN/N_uncorr");
             for (idx, (n, h)) in timeseries_lengths.iter().zip(histograms.iter()).enumerate() {
                 println!("{:?}\t\t{:?}\t\t{:?}\t\t{:.2}",
@@ -94,15 +156,69 @@ pub fn read_data(cfg: &Config) -> Result<Dataset> {
             let total_n = timeseries_lengths.iter().sum::<usize>() as f64;
             let total_h = histograms.iter().map(|h| h.num_points).sum::<u32>() as f64;
             println!("\t\t\t\t\tTotal:\t{:.2}", total_h/total_n);
-               
         }
 
-        Ok(Dataset::new(num_bins, dimens_length, bin_width, cfg.hist_min.clone(), cfg.hist_max.clone(), bias_pos, bias_fc, kT, histograms, cfg.cyclic))
+        Ok(datasets)
-    } else {
-        bail!("Histogram has no datapoints.")
     }
 }
 
+// builds a time boundaries for datasets from convdt, timeseries start and end
+fn get_convdt_boundaries(timeseries: &[f64], cfg: &Config) -> Vec<(f64, f64)> {
+    let mut last_timestep = *timeseries.last().unwrap();
+    if last_timestep > cfg.end {
+        last_timestep = cfg.end;
+    }
+    let mut first_timestep = *timeseries.first().unwrap();
+    if first_timestep < cfg.start {
+        first_timestep = cfg.start;
+    }
+    println!("{} to {} with dt={}", first_timestep, last_timestep, cfg.convdt);
+    if cfg.convdt == 0.0 {
+        vec![(0.0, last_timestep)]
+    } else {
+        let intervals: usize = ((last_timestep - first_timestep) / cfg.convdt).ceil() as usize;
+        println!("{:?}", intervals);
+        (1..intervals+1).map(|i| {
+            i as f64 * cfg.convdt + first_timestep
+        }).map(|end| { (first_timestep, end) }).collect()
+    }
+}
+
+// build a histogram from a timeseries
+// mask is used to filter the timeseries for selected frames
+fn build_histogram_from_timeseries(timeseries: &[Vec<f64>], mask: &[bool],
+    cfg: &Config) -> Histogram {
+
+    // total number of bins is the product of all dimensions length
+    let total_bins = cfg.num_bins.iter().product();
+
+    // bin width for each dimension: (max-min)/bins
+    let bin_width: Vec<f64> = (0..cfg.dimens).map(|idx| {
+        (cfg.hist_max[idx] - cfg.hist_min[idx])/(cfg.num_bins[idx] as f64)
+    }).collect();
+
+    // build histogram for slice start..convdt_stop
+    let mut hist = vec![0.0; total_bins];
+    for i in (0..timeseries[0].len()).filter(|i| mask[*i]) {
+        let mut values: Vec<f64> = vec![f64::NAN; cfg.dimens+1];
+        for j in 0..values.len() {
+            values[j] = timeseries[j][i];
+        }
+
+        if is_in_hist_boundaries(&values[1..], 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;
+        }
+    }
+
+    let num_points: f64 = hist.iter().sum();
+    Histogram::new(num_points as u32, hist)    
+}
+
 // transforms a multidimensional index into a one dimensional index
 // indeces: multidimensional indeces
 // lengths: length of the matrix in each dimension
@@ -125,50 +241,40 @@ fn is_in_hist_boundaries(values: &[f64], cfg: &Config) -> bool {
 }
 
 // returns true given time in inside the time boundaries defined by cfg
-fn is_in_time_boundaries(time: f64, cfg: &Config) -> bool {
+fn is_in_time_boundaries(time: f64, start: f64, end: f64) -> bool {
-    if cfg.start <= time && time <= cfg.end {
+    if start <= time && time <= end {
         return true
     }
     false
 }
 
-
+// parse a time series file
-// parse a time series file into a histogram
+fn read_window_file(window_file: &str, cfg: &Config) -> Result<(Vec<Vec<f64>>, usize)> {
-fn read_window_file(window_file: &str, cfg: &Config) -> Result<(Histogram, usize)> {
-    // 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];
-
-    // bin width for each dimension: (max-min)/bins
-    let bin_width: Vec<f64> = (0..cfg.dimens).map(|idx| {
-            (cfg.hist_max[idx] - cfg.hist_min[idx])/(cfg.num_bins[idx] as f64)
-        }).collect();
-
     let mut timeseries: Vec<Vec<f64>> = read_timeseries(window_file, cfg)?;
-    let timeseries_inital_length = timeseries[0].len();
 
+    // filter the timeseries based on start/end parameters
+    let time_series_mask: Vec<bool> = timeseries[0].iter()
+        .map(|t| is_in_time_boundaries(*t, cfg.start, cfg.end)).collect();
+    timeseries = timeseries.into_iter().map(|ts| {
+        ts.into_iter().zip(time_series_mask.iter()).filter_map(|(val, mask)| {
+            if *mask {
+                Some(val)
+            } else {
+                None
+            }
+        }).collect()
+    }).collect::<Vec<Vec<f64>>>();
+
+    let timeseries_inital_length = timeseries[0].len();
     if cfg.uncorr {
         timeseries = uncorrelate(timeseries, cfg);
     }
 
-    for i in 0..timeseries[0].len() {
+    if timeseries[0].is_empty() {
-        let mut values: Vec<f64> = vec![f64::NAN; cfg.dimens+1];
+        bail!("Time series is empty")
-        for j in 0..values.len() {
-            values[j] = timeseries[j][i];
-        }
-
-        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;
-        }
     }
 
-    let num_points: f64 = hist.iter().sum();
+    Ok((timeseries, timeseries_inital_length))
-    Ok((Histogram::new(num_points as u32, hist), timeseries_inital_length))
 }
 
 // Read a multidimensional timeseries
@@ -245,14 +351,24 @@ fn uncorrelate(timeseries: Vec<Vec<f64>>, cfg: &Config) -> Vec<Vec<f64>> {
 }
 
 // Write WHAM calculation results to out_file.
-pub fn write_results(out_file: &str, ds: &Dataset, free: &[f64],
+pub fn write_results(out_file: &str, append: bool, ds: &Dataset, free: &[f64],
-    free_std: &[f64], prob: &[f64], prob_std: &[f64]) -> Result<()> {
+    free_std: &[f64], prob: &[f64], prob_std: &[f64], index: Option<usize>) -> Result<()> {
-    let output = File::create(out_file)
+
+    if !append && Path::new(out_file).exists() {
+        std::fs::remove_file(out_file).chain_err(|| "Failed to delete file.")?;
+    }
+    let output = OpenOptions::new().write(true)
+        .append(true)
+        .create(true)
+        .open(out_file)
         .chain_err(|| format!("Failed to create file with path {}", out_file))?;
     let mut buf = BufWriter::new(output);
 
     let header: String = (0..ds.dimens_lengths.len()).map(|d| format!("coord{}", d+1))
         .collect::<Vec<String>>().join("    ");
+    if let Some(index) = index {
+        writeln!(buf, "#Dataset {}", index).unwrap();
+    }
     writeln!(buf, "#{}    Free Energy    +/-    Probability    +/-", header).unwrap();
 
     for bin in 0..free.len() {
@@ -290,6 +406,7 @@ mod tests {
             start: 0.0,
             end: 1e+20,
             uncorr: false,
+            convdt: 0.0,
         }
     }
 
@@ -297,7 +414,9 @@ mod tests {
     fn read_window_file() {
         let f = "example/1d_cyclic/COLVAR+0.0.xvg";
         let cfg = cfg();
-        let (h, timeseries_inital_length) = super::read_window_file(&f, &cfg).unwrap();
+        let (timeseries, timeseries_inital_length) = super::read_window_file(&f, &cfg).unwrap();
+        let mask = vec![true; timeseries[0].len()];
+        let h = build_histogram_from_timeseries(&timeseries, &mask, &cfg);
         println!("{:?}", h);
         assert_eq!(5000, timeseries_inital_length);
         assert_eq!(5000, h.num_points);
@@ -333,7 +452,7 @@ mod tests {
     #[test]
     fn read_data() {
         let cfg = cfg();
-        let ds = super::read_data(&cfg).unwrap();
+        let ds = &super::read_data(&cfg).unwrap()[0];
         println!("{:?}", ds);
         assert_eq!(25, ds.num_windows);
         assert_eq!(cfg.num_bins.len(), ds.dimens_lengths.len());
@@ -352,4 +471,73 @@ mod tests {
         let relative3 = super::get_relative_path(&path1, &path3);
         assert_eq!("path/to/subfolder/another_file.dat" ,relative3);
     }
+
+    #[test]
+    fn is_in_time_boundaries() {
+        let start = 10.0;
+        let end = 20.0;
+        assert!(super::is_in_time_boundaries(15.0, start, end));
+        assert!(super::is_in_time_boundaries(10.0, start, end));
+        assert!(super::is_in_time_boundaries(20.0, start, end));
+        assert!(!super::is_in_time_boundaries(9.9999999, start, end));
+        assert!(!super::is_in_time_boundaries(20.000001, start, end));   
+    }
+
+    #[test]
+    fn get_convdt_boundaries() {
+        let mut cfg = cfg();
+        
+        let timeseries: Vec<f64> = (0..31).map(|i| i as f64).collect();
+        println!("{:?}", timeseries);
+
+        cfg.start = 10.0;
+        cfg.end = 20.0;
+        cfg.convdt = 10.0;
+        let test = super::get_convdt_boundaries(&timeseries, &cfg);
+        println!("{:?}", test);
+        assert!(test.len() == 1);
+        assert_approx_eq!(test[0].0, 10.0);
+        assert_approx_eq!(test[0].1, 20.0);
+
+        cfg.start = 10.0;
+        cfg.end = 20.0;
+        cfg.convdt = 5.0;
+        let test = super::get_convdt_boundaries(&timeseries, &cfg);
+        println!("{:?}", test);
+        assert!(test.len() == 2);
+        assert_approx_eq!(test[0].0, 10.0);
+        assert_approx_eq!(test[0].1, 15.0);
+        assert_approx_eq!(test[1].0, 10.0);
+        assert_approx_eq!(test[1].1, 20.0);
+
+        let timeseries: Vec<f64> = (10..21).map(|i| i as f64).collect();
+        println!("{:?}", timeseries);
+
+        cfg.start = 10.0;
+        cfg.end = 20.0;
+        cfg.convdt = 10.0;
+        let test = super::get_convdt_boundaries(&timeseries, &cfg);
+        println!("{:?}", test);
+        assert!(test.len() == 1);
+        assert_approx_eq!(test[0].0, 10.0);
+        assert_approx_eq!(test[0].1, 20.0);
+
+        cfg.start = 5.0;
+        cfg.end = 20.0;
+        cfg.convdt = 10.0;
+        let test = super::get_convdt_boundaries(&timeseries, &cfg);
+        println!("{:?}", test);
+        assert!(test.len() == 1);
+        assert_approx_eq!(test[0].0, 10.0);
+        assert_approx_eq!(test[0].1, 20.0);
+
+        cfg.start = 5.0;
+        cfg.end = 30.0;
+        cfg.convdt = 10.0;
+        let test = super::get_convdt_boundaries(&timeseries, &cfg);
+        println!("{:?}", test);
+        assert!(test.len() == 1);
+        assert_approx_eq!(test[0].0, 10.0);
+        assert_approx_eq!(test[0].1, 20.0);
+    }
 }

src/lib.rs

@@ -47,16 +47,19 @@ pub struct Config {
     pub start: f64,
     pub end: f64,
     pub uncorr: bool,
+    pub convdt: f64,
 }
 
 impl fmt::Display for Config {
 	 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
          write!(f, "Metadata={}, hist_min={:?}, hist_max={:?}, bins={:?}, 
             verbose={}, tolerance={}, iterations={}, temperature={},
-            cyclic={:?}, uncorr={:?}, bootstrap={:?}, seed={:?}",
+            cyclic={:?}, uncorr={:?}, bootstrap={:?}, seed={:?},
+            uncorr={:?}, start={:?}, end={:?}, convdt={:?}",
             self.metadata_file, self.hist_min, self.hist_max, self.num_bins,
             self.verbose, self.tolerance, self.max_iterations, self.temperature,
-            self.cyclic, self.uncorr, self.bootstrap, self.bootstrap_seed)
+            self.cyclic, self.uncorr, self.bootstrap, self.bootstrap_seed,
+            self.uncorr, self.start, self.end, self.convdt)
     }
 }
 
@@ -183,27 +186,40 @@ pub fn run(cfg: &Config) -> Result<()>{
     println!("Supplied WHAM options: {}", &cfg);
 
     println!("Reading input files.");
-    let dataset = io::read_data(&cfg).chain_err(|| "Failed to create histogram.")?;
+    let datasets = io::read_data(&cfg).chain_err(|| "Failed to read data.")?;
-    println!("{}", &dataset);
 
-    let (P, F, F_prev) = perform_wham(&cfg, &dataset)?;
+    for (idx, dataset) in datasets.iter().enumerate() {
-    println!("WHAM converged.");
+        if datasets.len() > 1 {
+            println!("Dataset {}/{}: {}", idx+1, datasets.len(), &dataset);
+        }
+        else {
+            println!("{}", &dataset);
+        }
+        let (P, F, F_prev) = perform_wham(&cfg, &dataset)?;
+        println!("WHAM converged.");
 
-	let (P_std, free_energy_std) = if cfg.bootstrap > 0 {
+        let (P_std, free_energy_std) = if cfg.bootstrap > 0 {
-        println!("Bootstrapping..");      
+            println!("Bootstrapping..");      
-		error_analysis::run_bootstrap(&cfg, dataset.clone(), cfg.bootstrap)
+            error_analysis::run_bootstrap(&cfg, dataset.clone(), cfg.bootstrap)
-	} else {
+        } else {
-		(vec![0.0; P.len()], vec![0.0; P.len()])
+            (vec![0.0; P.len()], vec![0.0; P.len()])
-	};
+        };
 
-    // calculate free energy and dump state
+        // calculate free energy and dump state
-    println!("Finished. Dumping final PMF");
+        println!("Finished. Dumping PMF");
-	let free_energy = calc_free_energy(&dataset, &P);
+        let free_energy = calc_free_energy(&dataset, &P);
-    dump_state(&dataset, &F, &F_prev, &P, &P_std, &free_energy, &free_energy_std);
-
-    io::write_results(&cfg.output, &dataset, &free_energy, &free_energy_std, &P, &P_std)
-		.chain_err(|| "Could not write results to output file")?;
 
+        dump_state(&dataset, &F, &F_prev, &P, &P_std, &free_energy, &free_energy_std);
+        let append = idx > 0 && datasets.len() > 1;
+        let index = if datasets.len() > 1 {
+            Some(idx)
+        } else {
+            None
+        };
+        io::write_results(&cfg.output, append, &dataset, &free_energy, &free_energy_std, &P, &P_std, index)
+            .chain_err(|| "Could not write results to output file")?;
+    }
+    
     Ok(())
 }
 

src/main.rs

@@ -68,10 +68,12 @@ fn cli() -> Result<Config> {
     }
 
     let dimens = num_bins.len();
+    let convdt: f64 = matches.value_of("convdt").unwrap_or("0").parse()
+        .chain_err(|| "Cannot parse convdt.")?;
 
 	Ok(wham::Config{metadata_file, hist_min, hist_max, num_bins, dimens,
 		verbose, tolerance, max_iterations, temperature, cyclic, output,
-		bootstrap, bootstrap_seed, start, end, uncorr})
+		bootstrap, bootstrap_seed, start, end, uncorr, convdt})
 }
 
 fn main() {