faster tests

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.2"
 dependencies = [
  "assert_approx_eq",
  "clap",

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "wham"
-version = "1.0.0"
+version = "1.1.2"
 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
-D. Bauer <bauer@bio.tu-darmstadt.de>
+wham 1.1.0
+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

example/1d_cyclic/wham.out

@@ -1,101 +1,101 @@
 #coord1    Free Energy    +/-    Probability    +/-
--3.110177    7.158102    0.071438    0.003494    0.000068
--3.047345    5.365727    0.069449    0.007168    0.000135
--2.984513    3.873190    0.067495    0.013039    0.000232
--2.921681    2.953162    0.067589    0.018855    0.000343
--2.858849    1.949554    0.064640    0.028195    0.000480
--2.796017    1.391747    0.063570    0.035261    0.000584
--2.733186    1.128270    0.061710    0.039189    0.000620
--2.670354    0.839970    0.060445    0.043991    0.000667
--2.607522    0.624769    0.060762    0.047955    0.000739
--2.544690    0.663757    0.060786    0.047211    0.000731
--2.481858    1.051932    0.059774    0.040407    0.000617
--2.419026    1.463048    0.060422    0.034267    0.000527
--2.356194    1.990616    0.055242    0.027734    0.000360
--2.293363    2.190692    0.044458    0.025597    0.000210
--2.230531    2.553036    0.042144    0.022136    0.000159
--2.167699    2.572522    0.043313    0.021964    0.000165
--2.104867    2.472360    0.039799    0.022863    0.000157
--2.042035    2.517562    0.036792    0.022453    0.000203
--1.979203    2.469778    0.036324    0.022887    0.000237
--1.916372    2.223125    0.049496    0.025266    0.000565
--1.853540    2.080157    0.039205    0.026756    0.000453
--1.790708    1.793841    0.041488    0.030011    0.000558
--1.727876    1.458784    0.042362    0.034326    0.000565
--1.665044    0.914441    0.037561    0.042697    0.000682
--1.602212    0.303813    0.028701    0.054540    0.000761
--1.539380    0.268335    0.025765    0.055321    0.000771
--1.476549    0.000000    0.023371    0.061604    0.000866
--1.413717    0.537061    0.026370    0.049671    0.000715
--1.350885    1.439940    0.025620    0.034586    0.000547
--1.288053    2.391838    0.027653    0.023614    0.000393
--1.225221    3.779470    0.027093    0.013538    0.000229
--1.162389    5.685555    0.027762    0.006305    0.000111
--1.099557    7.661896    0.028996    0.002855    0.000051
--1.036726    9.946594    0.032118    0.001142    0.000021
--0.973894    12.359408    0.042801    0.000434    0.000009
--0.911062    14.954655    0.052443    0.000153    0.000004
--0.848230    17.742242    0.063261    0.000050    0.000002
--0.785398    20.555785    0.068133    0.000016    0.000001
--0.722566    22.811160    0.075518    0.000007    0.000000
--0.659734    25.178094    0.086226    0.000003    0.000000
--0.596903    26.442288    0.091783    0.000002    0.000000
--0.534071    27.897565    0.093175    0.000001    0.000000
--0.471239    29.062473    0.096027    0.000001    0.000000
--0.408407    30.384521    0.096855    0.000000    0.000000
--0.345575    31.638454    0.099956    0.000000    0.000000
--0.282743    32.817727    0.107689    0.000000    0.000000
--0.219911    33.770369    0.108946    0.000000    0.000000
--0.157080    34.505503    0.112902    0.000000    0.000000
--0.094248    35.431659    0.120856    0.000000    0.000000
--0.031416    35.615810    0.122046    0.000000    0.000000
-0.031416    35.561946    0.118225    0.000000    0.000000
-0.094248    35.382089    0.108652    0.000000    0.000000
-0.157080    34.934827    0.118384    0.000000    0.000000
-0.219911    33.673460    0.119195    0.000000    0.000000
-0.282743    32.731563    0.121592    0.000000    0.000000
-0.345575    31.261855    0.129698    0.000000    0.000000
-0.408407    29.717377    0.139789    0.000000    0.000000
-0.471239    28.076620    0.137289    0.000001    0.000000
-0.534071    26.481097    0.135712    0.000002    0.000000
-0.596903    24.487358    0.135156    0.000003    0.000000
-0.659734    22.344251    0.131027    0.000008    0.000000
-0.722566    20.241543    0.133145    0.000018    0.000001
-0.785398    18.341869    0.131952    0.000039    0.000002
-0.848230    16.261582    0.134676    0.000091    0.000006
-0.911062    14.301801    0.134417    0.000199    0.000013
-0.973894    12.603788    0.131018    0.000394    0.000024
-1.036726    11.249601    0.131276    0.000678    0.000043
-1.099557    10.087886    0.132498    0.001079    0.000070
-1.162389    9.443303    0.132990    0.001398    0.000092
-1.225221    9.152799    0.132250    0.001570    0.000104
-1.288053    9.331937    0.133178    0.001462    0.000099
-1.350885    9.905546    0.133357    0.001161    0.000078
-1.413717    11.042050    0.133807    0.000736    0.000051
-1.476549    12.598167    0.132722    0.000395    0.000027
-1.539380    14.520167    0.131816    0.000183    0.000012
-1.602212    16.569783    0.131773    0.000080    0.000005
-1.665044    18.687390    0.132601    0.000034    0.000002
-1.727876    20.775408    0.133867    0.000015    0.000001
-1.790708    22.905200    0.129265    0.000006    0.000000
-1.853540    24.643852    0.128894    0.000003    0.000000
-1.916372    26.301740    0.130266    0.000002    0.000000
-1.979203    27.372071    0.128620    0.000001    0.000000
-2.042035    28.697726    0.133263    0.000001    0.000000
-2.104867    29.417901    0.133513    0.000000    0.000000
-2.167699    30.008351    0.130925    0.000000    0.000000
-2.230531    30.406016    0.124398    0.000000    0.000000
-2.293363    30.171275    0.122493    0.000000    0.000000
-2.356194    29.884646    0.130278    0.000000    0.000000
-2.419026    29.428153    0.130844    0.000000    0.000000
-2.481858    28.546982    0.148114    0.000001    0.000000
-2.544690    27.757520    0.133045    0.000001    0.000000
-2.607522    26.505787    0.134364    0.000001    0.000000
-2.670354    24.491866    0.115061    0.000003    0.000000
-2.733186    22.320664    0.110036    0.000008    0.000000
-2.796017    20.052723    0.107894    0.000020    0.000001
-2.858849    17.655650    0.105964    0.000052    0.000002
-2.921681    15.471590    0.107005    0.000125    0.000005
-2.984513    13.138167    0.099378    0.000318    0.000012
-3.047345    11.092386    0.087770    0.000722    0.000022
-3.110177    9.065722    0.077092    0.001626    0.000038
+-3.110177    7.158102    0.000000    0.003494    0.000000
+-3.047345    5.365727    0.000000    0.007168    0.000000
+-2.984513    3.873190    0.000000    0.013039    0.000000
+-2.921681    2.953162    0.000000    0.018855    0.000000
+-2.858849    1.949554    0.000000    0.028195    0.000000
+-2.796017    1.391747    0.000000    0.035261    0.000000
+-2.733186    1.128270    0.000000    0.039189    0.000000
+-2.670354    0.839970    0.000000    0.043991    0.000000
+-2.607522    0.624769    0.000000    0.047955    0.000000
+-2.544690    0.663757    0.000000    0.047211    0.000000
+-2.481858    1.051932    0.000000    0.040407    0.000000
+-2.419026    1.463048    0.000000    0.034267    0.000000
+-2.356194    1.990616    0.000000    0.027734    0.000000
+-2.293363    2.190692    0.000000    0.025597    0.000000
+-2.230531    2.553036    0.000000    0.022136    0.000000
+-2.167699    2.572522    0.000000    0.021964    0.000000
+-2.104867    2.472360    0.000000    0.022863    0.000000
+-2.042035    2.517562    0.000000    0.022453    0.000000
+-1.979203    2.469778    0.000000    0.022887    0.000000
+-1.916372    2.223125    0.000000    0.025266    0.000000
+-1.853540    2.080157    0.000000    0.026756    0.000000
+-1.790708    1.793841    0.000000    0.030011    0.000000
+-1.727876    1.458784    0.000000    0.034326    0.000000
+-1.665044    0.914441    0.000000    0.042697    0.000000
+-1.602212    0.303813    0.000000    0.054540    0.000000
+-1.539380    0.268335    0.000000    0.055321    0.000000
+-1.476549    0.000000    0.000000    0.061604    0.000000
+-1.413717    0.537061    0.000000    0.049671    0.000000
+-1.350885    1.439940    0.000000    0.034586    0.000000
+-1.288053    2.391838    0.000000    0.023614    0.000000
+-1.225221    3.779470    0.000000    0.013538    0.000000
+-1.162389    5.685555    0.000000    0.006305    0.000000
+-1.099557    7.661896    0.000000    0.002855    0.000000
+-1.036726    9.946594    0.000000    0.001142    0.000000
+-0.973894    12.359408    0.000000    0.000434    0.000000
+-0.911062    14.954655    0.000000    0.000153    0.000000
+-0.848230    17.742242    0.000000    0.000050    0.000000
+-0.785398    20.555785    0.000000    0.000016    0.000000
+-0.722566    22.811160    0.000000    0.000007    0.000000
+-0.659734    25.178094    0.000000    0.000003    0.000000
+-0.596903    26.442288    0.000000    0.000002    0.000000
+-0.534071    27.897565    0.000000    0.000001    0.000000
+-0.471239    29.062473    0.000000    0.000001    0.000000
+-0.408407    30.384521    0.000000    0.000000    0.000000
+-0.345575    31.638454    0.000000    0.000000    0.000000
+-0.282743    32.817727    0.000000    0.000000    0.000000
+-0.219911    33.770369    0.000000    0.000000    0.000000
+-0.157080    34.505503    0.000000    0.000000    0.000000
+-0.094248    35.431659    0.000000    0.000000    0.000000
+-0.031416    35.615810    0.000000    0.000000    0.000000
+0.031416    35.561946    0.000000    0.000000    0.000000
+0.094248    35.382089    0.000000    0.000000    0.000000
+0.157080    34.934827    0.000000    0.000000    0.000000
+0.219911    33.673460    0.000000    0.000000    0.000000
+0.282743    32.731563    0.000000    0.000000    0.000000
+0.345575    31.261855    0.000000    0.000000    0.000000
+0.408407    29.717377    0.000000    0.000000    0.000000
+0.471239    28.076620    0.000000    0.000001    0.000000
+0.534071    26.481097    0.000000    0.000002    0.000000
+0.596903    24.487358    0.000000    0.000003    0.000000
+0.659734    22.344251    0.000000    0.000008    0.000000
+0.722566    20.241543    0.000000    0.000018    0.000000
+0.785398    18.341869    0.000000    0.000039    0.000000
+0.848230    16.261582    0.000000    0.000091    0.000000
+0.911062    14.301801    0.000000    0.000199    0.000000
+0.973894    12.603788    0.000000    0.000394    0.000000
+1.036726    11.249601    0.000000    0.000678    0.000000
+1.099557    10.087886    0.000000    0.001079    0.000000
+1.162389    9.443303    0.000000    0.001398    0.000000
+1.225221    9.152799    0.000000    0.001570    0.000000
+1.288053    9.331937    0.000000    0.001462    0.000000
+1.350885    9.905546    0.000000    0.001161    0.000000
+1.413717    11.042050    0.000000    0.000736    0.000000
+1.476549    12.598167    0.000000    0.000395    0.000000
+1.539380    14.520167    0.000000    0.000183    0.000000
+1.602212    16.569783    0.000000    0.000080    0.000000
+1.665044    18.687390    0.000000    0.000034    0.000000
+1.727876    20.775408    0.000000    0.000015    0.000000
+1.790708    22.905200    0.000000    0.000006    0.000000
+1.853540    24.643852    0.000000    0.000003    0.000000
+1.916372    26.301740    0.000000    0.000002    0.000000
+1.979203    27.372071    0.000000    0.000001    0.000000
+2.042035    28.697726    0.000000    0.000001    0.000000
+2.104867    29.417901    0.000000    0.000000    0.000000
+2.167699    30.008351    0.000000    0.000000    0.000000
+2.230531    30.406016    0.000000    0.000000    0.000000
+2.293363    30.171275    0.000000    0.000000    0.000000
+2.356194    29.884646    0.000000    0.000000    0.000000
+2.419026    29.428153    0.000000    0.000000    0.000000
+2.481858    28.546982    0.000000    0.000001    0.000000
+2.544690    27.757520    0.000000    0.000001    0.000000
+2.607522    26.505787    0.000000    0.000001    0.000000
+2.670354    24.491866    0.000000    0.000003    0.000000
+2.733186    22.320664    0.000000    0.000008    0.000000
+2.796017    20.052723    0.000000    0.000020    0.000000
+2.858849    17.655650    0.000000    0.000052    0.000000
+2.921681    15.471590    0.000000    0.000125    0.000000
+2.984513    13.138167    0.000000    0.000318    0.000000
+3.047345    11.092386    0.000000    0.000722    0.000000
+3.110177    9.065722    0.000000    0.001626    0.000000

example/1d_cyclic/wham_bt.out

@@ -0,0 +1,101 @@
+#coord1    Free Energy    +/-    Probability    +/-
+-3.110177    7.158102    0.071438    0.003494    0.000068
+-3.047345    5.365727    0.069449    0.007168    0.000135
+-2.984513    3.873190    0.067495    0.013039    0.000232
+-2.921681    2.953162    0.067589    0.018855    0.000343
+-2.858849    1.949554    0.064640    0.028195    0.000480
+-2.796017    1.391747    0.063570    0.035261    0.000584
+-2.733186    1.128270    0.061710    0.039189    0.000620
+-2.670354    0.839970    0.060445    0.043991    0.000667
+-2.607522    0.624769    0.060762    0.047955    0.000739
+-2.544690    0.663757    0.060786    0.047211    0.000731
+-2.481858    1.051932    0.059774    0.040407    0.000617
+-2.419026    1.463048    0.060422    0.034267    0.000527
+-2.356194    1.990616    0.055242    0.027734    0.000360
+-2.293363    2.190692    0.044458    0.025597    0.000210
+-2.230531    2.553036    0.042144    0.022136    0.000159
+-2.167699    2.572522    0.043313    0.021964    0.000165
+-2.104867    2.472360    0.039799    0.022863    0.000157
+-2.042035    2.517562    0.036792    0.022453    0.000203
+-1.979203    2.469778    0.036324    0.022887    0.000237
+-1.916372    2.223125    0.049496    0.025266    0.000565
+-1.853540    2.080157    0.039205    0.026756    0.000453
+-1.790708    1.793841    0.041488    0.030011    0.000558
+-1.727876    1.458784    0.042362    0.034326    0.000565
+-1.665044    0.914441    0.037561    0.042697    0.000682
+-1.602212    0.303813    0.028701    0.054540    0.000761
+-1.539380    0.268335    0.025765    0.055321    0.000771
+-1.476549    0.000000    0.023371    0.061604    0.000866
+-1.413717    0.537061    0.026370    0.049671    0.000715
+-1.350885    1.439940    0.025620    0.034586    0.000547
+-1.288053    2.391838    0.027653    0.023614    0.000393
+-1.225221    3.779470    0.027093    0.013538    0.000229
+-1.162389    5.685555    0.027762    0.006305    0.000111
+-1.099557    7.661896    0.028996    0.002855    0.000051
+-1.036726    9.946594    0.032118    0.001142    0.000021
+-0.973894    12.359408    0.042801    0.000434    0.000009
+-0.911062    14.954655    0.052443    0.000153    0.000004
+-0.848230    17.742242    0.063261    0.000050    0.000002
+-0.785398    20.555785    0.068133    0.000016    0.000001
+-0.722566    22.811160    0.075518    0.000007    0.000000
+-0.659734    25.178094    0.086226    0.000003    0.000000
+-0.596903    26.442288    0.091783    0.000002    0.000000
+-0.534071    27.897565    0.093175    0.000001    0.000000
+-0.471239    29.062473    0.096027    0.000001    0.000000
+-0.408407    30.384521    0.096855    0.000000    0.000000
+-0.345575    31.638454    0.099956    0.000000    0.000000
+-0.282743    32.817727    0.107689    0.000000    0.000000
+-0.219911    33.770369    0.108946    0.000000    0.000000
+-0.157080    34.505503    0.112902    0.000000    0.000000
+-0.094248    35.431659    0.120856    0.000000    0.000000
+-0.031416    35.615810    0.122046    0.000000    0.000000
+0.031416    35.561946    0.118225    0.000000    0.000000
+0.094248    35.382089    0.108652    0.000000    0.000000
+0.157080    34.934827    0.118384    0.000000    0.000000
+0.219911    33.673460    0.119195    0.000000    0.000000
+0.282743    32.731563    0.121592    0.000000    0.000000
+0.345575    31.261855    0.129698    0.000000    0.000000
+0.408407    29.717377    0.139789    0.000000    0.000000
+0.471239    28.076620    0.137289    0.000001    0.000000
+0.534071    26.481097    0.135712    0.000002    0.000000
+0.596903    24.487358    0.135156    0.000003    0.000000
+0.659734    22.344251    0.131027    0.000008    0.000000
+0.722566    20.241543    0.133145    0.000018    0.000001
+0.785398    18.341869    0.131952    0.000039    0.000002
+0.848230    16.261582    0.134676    0.000091    0.000006
+0.911062    14.301801    0.134417    0.000199    0.000013
+0.973894    12.603788    0.131018    0.000394    0.000024
+1.036726    11.249601    0.131276    0.000678    0.000043
+1.099557    10.087886    0.132498    0.001079    0.000070
+1.162389    9.443303    0.132990    0.001398    0.000092
+1.225221    9.152799    0.132250    0.001570    0.000104
+1.288053    9.331937    0.133178    0.001462    0.000099
+1.350885    9.905546    0.133357    0.001161    0.000078
+1.413717    11.042050    0.133807    0.000736    0.000051
+1.476549    12.598167    0.132722    0.000395    0.000027
+1.539380    14.520167    0.131816    0.000183    0.000012
+1.602212    16.569783    0.131773    0.000080    0.000005
+1.665044    18.687390    0.132601    0.000034    0.000002
+1.727876    20.775408    0.133867    0.000015    0.000001
+1.790708    22.905200    0.129265    0.000006    0.000000
+1.853540    24.643852    0.128894    0.000003    0.000000
+1.916372    26.301740    0.130266    0.000002    0.000000
+1.979203    27.372071    0.128620    0.000001    0.000000
+2.042035    28.697726    0.133263    0.000001    0.000000
+2.104867    29.417901    0.133513    0.000000    0.000000
+2.167699    30.008351    0.130925    0.000000    0.000000
+2.230531    30.406016    0.124398    0.000000    0.000000
+2.293363    30.171275    0.122493    0.000000    0.000000
+2.356194    29.884646    0.130278    0.000000    0.000000
+2.419026    29.428153    0.130844    0.000000    0.000000
+2.481858    28.546982    0.148114    0.000001    0.000000
+2.544690    27.757520    0.133045    0.000001    0.000000
+2.607522    26.505787    0.134364    0.000001    0.000000
+2.670354    24.491866    0.115061    0.000003    0.000000
+2.733186    22.320664    0.110036    0.000008    0.000000
+2.796017    20.052723    0.107894    0.000020    0.000001
+2.858849    17.655650    0.105964    0.000052    0.000002
+2.921681    15.471590    0.107005    0.000125    0.000005
+2.984513    13.138167    0.099378    0.000318    0.000012
+3.047345    11.092386    0.087770    0.000722    0.000022
+3.110177    9.065722    0.077092    0.001626    0.000038

src/cli.yml

@@ -1,5 +1,5 @@
 name: wham
-version: "1.0.0"
+version: "1.1.2"
 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,14 @@ 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
+    - ignore_empty:
+        long: ignore_empty
+        help: If this is set, do not fail if a histogram is empty.
+        takes_value: false
+        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
-pub fn read_data(cfg: &Config) -> Result<Dataset> {
+// given in the metadata file. This generates at least one 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,87 @@ 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() {
+
+    // Datasets are created from histograms.
+    // Empty histograms result in an error when its the final dataset,
+    // and a warning otherwise.
+    let num_datasets: usize = histograms.iter().map(|h| h.len()).max().unwrap();
+    let dataset_boundaries: Vec<(f64, f64)> = (0..num_datasets).map(|idx| {
+        (cfg.start, cfg.start+(idx as f64 + 1.0)*cfg.convdt) }
+    ).collect();
+    vprintln(format!("Generating {} datasets from histograms.", num_datasets), cfg.verbose);
+    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 for interval {}-{} in histogram boundaries: {}.",
+                    dataset_boundaries[idx].0, dataset_boundaries[idx].1 ,&path);
+                if !cfg.ignore_empty && idx+1 == num_datasets {
+                    bail!(warning + " This is the final dataset.");
+                } 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!();
+            println!("Timeseries Correlation:");
             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 +161,67 @@ 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))
-    } else {
-        bail!("Histogram has no datapoints.")
+        Ok(datasets)
     }
 }
 
+// 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;
+    }
+    if cfg.convdt == 0.0 {
+        vec![(0.0, last_timestep)]
+    } else {
+        let intervals: usize = ((last_timestep - first_timestep) / cfg.convdt).ceil() as usize;
+        (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 +244,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 {
-    if cfg.start <= time && time <= cfg.end {
+fn is_in_time_boundaries(time: f64, start: f64, end: f64) -> bool {
+    if start <= time && time <= end {
         return true
     }
     false
 }
 
-
-// parse a time series file into a histogram
-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();
-
+// parse a time series file
+fn read_window_file(window_file: &str, cfg: &Config) -> Result<(Vec<Vec<f64>>, usize)> {
     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() {
-        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) && 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 timeseries[0].is_empty() {
+        bail!("Time series is empty")
     }
 
-    let num_points: f64 = hist.iter().sum();
-    Ok((Histogram::new(num_points as u32, hist), timeseries_inital_length))
+    Ok((timeseries, timeseries_inital_length))
 }
 
 // Read a multidimensional timeseries
@@ -245,14 +354,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],
-    free_std: &[f64], prob: &[f64], prob_std: &[f64]) -> Result<()> {
-    let output = File::create(out_file)
+pub fn write_results(out_file: &str, append: bool, ds: &Dataset, free: &[f64],
+    free_std: &[f64], prob: &[f64], prob_std: &[f64], index: Option<usize>) -> Result<()> {
+
+    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 +409,8 @@ mod tests {
             start: 0.0,
             end: 1e+20,
             uncorr: false,
+            convdt: 0.0,
+            ignore_empty: false,
         }
     }
 
@@ -297,7 +418,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 +456,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 +475,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,20 @@ pub struct Config {
     pub start: f64,
     pub end: f64,
     pub uncorr: bool,
+    pub convdt: f64,
+    pub ignore_empty: bool
 }
 
 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={:?}, ignore_empty={:?}",
             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, self.ignore_empty)
     }
 }
 
@@ -183,27 +187,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.")?;
-    println!("{}", &dataset);
+    let datasets = io::read_data(&cfg).chain_err(|| "Failed to read data.")?;
 
-    let (P, F, F_prev) = perform_wham(&cfg, &dataset)?;
-    println!("WHAM converged.");
+    for (idx, dataset) in datasets.iter().enumerate() {
+        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 {
-        println!("Bootstrapping..");      
-		error_analysis::run_bootstrap(&cfg, dataset.clone(), cfg.bootstrap)
-	} else {
-		(vec![0.0; P.len()], vec![0.0; P.len()])
-	};
+        let (P_std, free_energy_std) = if cfg.bootstrap > 0 {
+            println!("Bootstrapping..");      
+            error_analysis::run_bootstrap(&cfg, dataset.clone(), cfg.bootstrap)
+        } else {
+            (vec![0.0; P.len()], vec![0.0; P.len()])
+        };
 
-    // calculate free energy and dump state
-    println!("Finished. Dumping final PMF");
-	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")?;
+        // calculate free energy and dump state
+        println!("Finished. Dumping PMF");
+        let free_energy = calc_free_energy(&dataset, &P);
 
+        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,14 @@ 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.")?;
+
+    let ignore_empty: bool = matches.is_present("ignore_empty");
 
 	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, ignore_empty})
 }
 
 fn main() {

tests/examples.rs

@@ -6,64 +6,154 @@ mod integration {
     use std::process::Command;
     use std::fs;
     use super::command::get_command;
+    use std::fs::OpenOptions;
+    use std::io::prelude::*;
 
     #[test]
     fn wham_1d_cyclic() {
+        let output_file = "/tmp/wham_test_1d_cyclic.out";
         get_command()
             .args(&["--bins", "100", "--max", "pi", "--min", "-pi", "-T", "300", "--cyclic"])
-            .args(&["--bt", "100", "--seed", "1234"])
+            .args(&["--seed", "1234"])
             .args(&["-f", "example/1d_cyclic/metadata.dat"])
-            .args(&["-o", "/tmp/wham_test_1d_cyclic.out"])
+            .args(&["-o", output_file])
             .output()
             .expect("failed to execute process");
 
-        assert!(fs::metadata("/tmp/wham_test_1d_cyclic.out").is_ok());
+        assert!(fs::metadata(output_file).is_ok());
         let output = Command::new("diff")
-            .arg("/tmp/wham_test_1d_cyclic.out")
+            .arg(output_file)
             .arg("example/1d_cyclic/wham.out")
             .output()
             .expect("failed to run diff");
         let output_len = String::from_utf8_lossy(&output.stdout).len();
         assert_eq!(output_len, 0);
+        std::fs::remove_file(output_file).unwrap();
     }
 
     #[test]
     fn wham_1d_cyclic_uncorrelated() {
+        let output_file = "/tmp/wham_test_1d_cyclic.out";
         get_command()
             .args(&["--bins", "100", "--max", "pi", "--min", "-pi", "-T", "300", "--cyclic", "--uncorr"])
             .args(&["--seed", "1234"])
             .args(&["-f", "example/1d_cyclic/metadata.dat"])
-            .args(&["-o", "/tmp/wham_test_1d_cyclic.out"])
+            .args(&["-o", output_file])
             .output()
             .expect("failed to execute process");
 
         assert!(fs::metadata("/tmp/wham_test_1d_cyclic.out").is_ok());
         let output = Command::new("diff")
-            .arg("/tmp/wham_test_1d_cyclic.out")
+            .arg(output_file)
             .arg("example/1d_cyclic/wham_uncorrelated.out")
             .output()
             .expect("failed to run diff");
         let output_len = String::from_utf8_lossy(&output.stdout).len();
         assert_eq!(output_len, 0);
+        std::fs::remove_file(output_file).unwrap();
+    }
+
+    #[test]
+    // Test if convdt runs have the same result as normal runs
+    fn wham_convdt() {
+        // run wham with convdt
+        let output_file = "/tmp/wham_test_convdt.out";
+        get_command()
+            .args(&["--bins", "10", "--max", "pi", "--min", "-pi", "-T", "300", "--cyclic"])
+            .args(&["--seed", "1234", "--tolerance", "0.001"])
+            .args(&["--start", "0", "--end", "10"])
+            .args(&["--convdt", "1"])
+            .args(&["-f", "example/1d_cyclic/metadata.dat"])
+            .args(&["-o", output_file])
+            .output()
+            .expect("failed to execute process");
+        assert!(fs::metadata(output_file).is_ok());
+
+        // run wham for individual sets
+        for i in 1..11 {
+            let output_file_single = format!("/tmp/wham_test_convdt_{}.out", i);
+            get_command()
+                .args(&["--bins", "10", "--max", "pi", "--min", "-pi", "-T", "300", "--cyclic"])
+                .args(&["--seed", "1234", "--tolerance", "0.001"])
+                .args(&["--start", "0", "--end", &i.to_string()])
+                .args(&["-f", "example/1d_cyclic/metadata.dat"])
+                .args(&["-o", &output_file_single])
+                .output()
+                .expect("failed to execute process");
+            assert!(fs::metadata(output_file_single).is_ok());
+        }
+
+        // combine individual runs
+        let output_combined = "/tmp/wham_test_convdt_combined.out";
+        let mut file = OpenOptions::new()
+            .create(true)
+            .write(true)
+            .open(output_combined)
+            .unwrap();
+        for i in 1..11 {
+            let output_file_single = format!("/tmp/wham_test_convdt_{}.out", i);
+            println!("{}", output_file_single);
+            file.write_all(format!("#Dataset {}\n", i-1).as_bytes()).unwrap();
+            file.write_all(fs::read_to_string(output_file_single.clone()).unwrap().as_bytes()).unwrap();
+            std::fs::remove_file(output_file_single).unwrap();
+        }
+
+        // compare combined runs with single run
+        let output = Command::new("diff")
+            .arg(output_file)
+            .arg(output_combined)
+            .output()
+            .expect("failed to run diff");
+        let output_len = String::from_utf8_lossy(&output.stdout).len();
+        assert_eq!(output_len, 0);
+
+        std::fs::remove_file(output_combined).unwrap();
+        std::fs::remove_file(output_file).unwrap();
+        
+    }
+
+    #[test]
+    #[ignore]
+    fn wham_1d_cyclic_bootstrap() {
+        let output_file = "/tmp/wham_test_1d_cyclic_bt.out";
+        get_command()
+            .args(&["--bins", "100", "--max", "pi", "--min", "-pi", "-T", "300", "--cyclic"])
+            .args(&["--seed", "1234", "--bt", "100"])
+            .args(&["-f", "example/1d_cyclic/metadata.dat"])
+            .args(&["-o", output_file])
+            .output()
+            .expect("failed to execute process");
+
+        assert!(fs::metadata(output_file).is_ok());
+        let output = Command::new("diff")
+            .arg(output_file)
+            .arg("example/1d_cyclic/wham_bt.out")
+            .output()
+            .expect("failed to run diff");
+        let output_len = String::from_utf8_lossy(&output.stdout).len();
+        assert_eq!(output_len, 0);
+        std::fs::remove_file(output_file).unwrap();
     }
 
     #[test]
     #[ignore] // expensive
     fn wham_2d_cyclic() {
+        let output_file = "/tmp/wham_test_2d_cyclic.out";
         get_command()
             .args(&["--bins", "100,100", "--max", "pi,pi", "--min", "-pi,-pi", "-T", "300", "--cyclic"])
             .args(&["-f", "example/2d_cyclic/metadata.dat"])
-            .args(&["-o", "/tmp/wham_test_2d_cyclic.out"])
+            .args(&["-o", output_file])
             .output()
             .expect("failed to execute process");
 
-        assert!(fs::metadata("/tmp/wham_test_2d_cyclic.out").is_ok());
+        assert!(fs::metadata(output_file).is_ok());
         let output = Command::new("diff")
-            .arg("/tmp/wham_test_2d_cyclic.out")
+            .arg(output_file)
             .arg("example/2d_cyclic/wham.out")
             .output()
             .expect("failed to run diff");
         let output_len = String::from_utf8_lossy(&output.stdout).len();
         assert_eq!(output_len, 0);
+        std::fs::remove_file(output_file).unwrap();
     }
-}
+}