bauer/WHAM
1
0
Fork 0
mirror of https://github.com/dnlbauer/WHAM.git synced 2026-09-10 22:25:31 +00:00
Code Issues Packages Projects Releases Wiki Activity

improved I/O error messages

Browse Source
This commit is contained in:
Daniel Bauer
2018-11-03 18:27:48 +01:00
parent 344066cdb8
commit 78c077a6f2
9 changed files with 217 additions and 78 deletions
Download Patch File Download Diff File Expand all files Collapse all files

51
Cargo.lock generated
View File

@@ -16,11 +16,42 @@ dependencies = [
"winapi 0.3.6 (registry+https://github.com/rust-lang/crates.io-index)",
]
[[package]]
name = "backtrace"
version = "0.3.9"
source = "registry+https://github.com/rust-lang/crates.io-index"
dependencies = [
"backtrace-sys 0.1.24 (registry+https://github.com/rust-lang/crates.io-index)",
"cfg-if 0.1.6 (registry+https://github.com/rust-lang/crates.io-index)",
"libc 0.2.43 (registry+https://github.com/rust-lang/crates.io-index)",
"rustc-demangle 0.1.9 (registry+https://github.com/rust-lang/crates.io-index)",
"winapi 0.3.6 (registry+https://github.com/rust-lang/crates.io-index)",
]
[[package]]
name = "backtrace-sys"
version = "0.1.24"
source = "registry+https://github.com/rust-lang/crates.io-index"
dependencies = [
"cc 1.0.25 (registry+https://github.com/rust-lang/crates.io-index)",
"libc 0.2.43 (registry+https://github.com/rust-lang/crates.io-index)",
]
[[package]]
name = "bitflags"
version = "1.0.4"
source = "registry+https://github.com/rust-lang/crates.io-index"
[[package]]
name = "cc"
version = "1.0.25"
source = "registry+https://github.com/rust-lang/crates.io-index"
[[package]]
name = "cfg-if"
version = "0.1.6"
source = "registry+https://github.com/rust-lang/crates.io-index"
[[package]]
name = "clap"
version = "2.32.0"
@@ -36,6 +67,14 @@ dependencies = [
"yaml-rust 0.3.5 (registry+https://github.com/rust-lang/crates.io-index)",
]
[[package]]
name = "error-chain"
version = "0.12.0"
source = "registry+https://github.com/rust-lang/crates.io-index"
dependencies = [
"backtrace 0.3.9 (registry+https://github.com/rust-lang/crates.io-index)",
]
[[package]]
name = "libc"
version = "0.2.43"
@@ -54,6 +93,11 @@ dependencies = [
"redox_syscall 0.1.40 (registry+https://github.com/rust-lang/crates.io-index)",
]
[[package]]
name = "rustc-demangle"
version = "0.1.9"
source = "registry+https://github.com/rust-lang/crates.io-index"
[[package]]
name = "strsim"
version = "0.7.0"
@@ -92,6 +136,7 @@ name = "wham"
version = "0.9.0"
dependencies = [
"clap 2.32.0 (registry+https://github.com/rust-lang/crates.io-index)",
"error-chain 0.12.0 (registry+https://github.com/rust-lang/crates.io-index)",
]
[[package]]
@@ -121,11 +166,17 @@ source = "registry+https://github.com/rust-lang/crates.io-index"
[metadata]
"checksum ansi_term 0.11.0 (registry+https://github.com/rust-lang/crates.io-index)" = "ee49baf6cb617b853aa8d93bf420db2383fab46d314482ca2803b40d5fde979b"
"checksum atty 0.2.11 (registry+https://github.com/rust-lang/crates.io-index)" = "9a7d5b8723950951411ee34d271d99dddcc2035a16ab25310ea2c8cfd4369652"
"checksum backtrace 0.3.9 (registry+https://github.com/rust-lang/crates.io-index)" = "89a47830402e9981c5c41223151efcced65a0510c13097c769cede7efb34782a"
"checksum backtrace-sys 0.1.24 (registry+https://github.com/rust-lang/crates.io-index)" = "c66d56ac8dabd07f6aacdaf633f4b8262f5b3601a810a0dcddffd5c22c69daa0"
"checksum bitflags 1.0.4 (registry+https://github.com/rust-lang/crates.io-index)" = "228047a76f468627ca71776ecdebd732a3423081fcf5125585bcd7c49886ce12"
"checksum cc 1.0.25 (registry+https://github.com/rust-lang/crates.io-index)" = "f159dfd43363c4d08055a07703eb7a3406b0dac4d0584d96965a3262db3c9d16"
"checksum cfg-if 0.1.6 (registry+https://github.com/rust-lang/crates.io-index)" = "082bb9b28e00d3c9d39cc03e64ce4cea0f1bb9b3fde493f0cbc008472d22bdf4"
"checksum clap 2.32.0 (registry+https://github.com/rust-lang/crates.io-index)" = "b957d88f4b6a63b9d70d5f454ac8011819c6efa7727858f458ab71c756ce2d3e"
"checksum error-chain 0.12.0 (registry+https://github.com/rust-lang/crates.io-index)" = "07e791d3be96241c77c43846b665ef1384606da2cd2a48730abe606a12906e02"
"checksum libc 0.2.43 (registry+https://github.com/rust-lang/crates.io-index)" = "76e3a3ef172f1a0b9a9ff0dd1491ae5e6c948b94479a3021819ba7d860c8645d"
"checksum redox_syscall 0.1.40 (registry+https://github.com/rust-lang/crates.io-index)" = "c214e91d3ecf43e9a4e41e578973adeb14b474f2bee858742d127af75a0112b1"
"checksum redox_termios 0.1.1 (registry+https://github.com/rust-lang/crates.io-index)" = "7e891cfe48e9100a70a3b6eb652fef28920c117d366339687bd5576160db0f76"
"checksum rustc-demangle 0.1.9 (registry+https://github.com/rust-lang/crates.io-index)" = "bcfe5b13211b4d78e5c2cadfebd7769197d95c639c35a50057eb4c05de811395"
"checksum strsim 0.7.0 (registry+https://github.com/rust-lang/crates.io-index)" = "bb4f380125926a99e52bc279241539c018323fab05ad6368b56f93d9369ff550"
"checksum termion 1.5.1 (registry+https://github.com/rust-lang/crates.io-index)" = "689a3bdfaab439fd92bc87df5c4c78417d3cbe537487274e9b0b2dce76e92096"
"checksum textwrap 0.10.0 (registry+https://github.com/rust-lang/crates.io-index)" = "307686869c93e71f94da64286f9a9524c0f308a9e1c87a583de8e9c9039ad3f6"

1
Cargo.toml
View File

@@ -5,6 +5,7 @@ authors = ["D. Bauer <bauer@bio.tu-darmstadt.de>"]
[dependencies]
clap = {version="2.32.0", features=['yaml']}
error-chain = "0.12.0"
[profile.release]
opt-level = 2

2
README.md
View File

@@ -61,10 +61,8 @@ TODO
---
- Multithreading (?)
- Error analysis / bootstrapping
- Better error messages during file I/O
- Autocorrelation
- Replica exchange
- Unit tests for 2d/Nd WHAM
License
---

100
src/io.rs
View File

@@ -5,11 +5,8 @@ use std::fs::File;
use std::io::prelude::*;
use std::io::{BufReader,BufWriter};
use k_B;
use std::process;
use std::option::Option;
use std::path::Path;
use std::error::Error;
use std::result::Result;
use super::errors::*;
// Returns the path to path2 relative to path1
// path1: "path/to/file.dat"
@@ -28,7 +25,7 @@ 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) -> Option<Dataset> {
pub fn read_data(cfg: &Config) -> Result<Dataset> {
let mut bias_pos: Vec<f64> = Vec::new();
let mut bias_fc: Vec<f64> = Vec::new();
let mut histograms: Vec<Histogram> = Vec::new();
@@ -40,62 +37,51 @@ pub fn read_data(cfg: &Config) -> Option<Dataset> {
let num_bins = cfg.num_bins.iter().fold(1, |state, &bins| state*bins);
let dimens_length = cfg.num_bins.clone();
let f = File::open(&cfg.metadata_file).unwrap_or_else(|x| {
eprintln!("Failed to read metadata from {}. {}", &cfg.metadata_file, x);
process::exit(1)
});
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 l in buf.lines() {
let line = l.unwrap();
for (line_num,l) in buf.lines().enumerate() {
let line = l.chain_err(|| "Failed to read line")?;
// skip comments and empty lines
if line.starts_with("#") || line.len() == 0 {
continue;
}
let mut split = line.split_whitespace();
let split: Vec<&str> = line.split_whitespace().collect();
if split.len() < 1 + cfg.dimens * 2 {
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.next()?);
match read_window_file(&path, cfg) {
Some(h) => {
histograms.push(h);
vprintln(format!("{}, {} data points added.", &path, histograms.last().unwrap().num_points), cfg.verbose);
},
None => {
eprintln!("No data points inside histogram boundaries: {}", &path);
process::exit(1)
}
let path = get_relative_path(&cfg.metadata_file, split[0]);
let h = 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);
vprintln(format!("{}, {} data points added.", &path,
histograms.last().unwrap().num_points), cfg.verbose);
// parse bias force constants and positions
for _ in 0..cfg.dimens {
match split.next()?.parse() {
Ok(x) => bias_pos.push(x),
_ => {
eprintln!("Failed to read bias coordinate.");
process::exit(1);
}
}
for i in 1..cfg.dimens+1 {
let pos = split[i].parse()
.chain_err(|| format!("Failed to read bias position in line {} of metadata file", line_num+1))?;
bias_pos.push(pos);
}
for _ in 0..cfg.dimens {
match split.next()?.parse() {
Ok(x) => bias_fc.push(x),
_ => {
eprintln!("Failed to read bias force constant.");
process::exit(1);
}
}
for i in (1+cfg.dimens)..(1+2*cfg.dimens) {
let fc = split[i].parse()
.chain_err(|| format!("Failed to read bias fc in line {} of metadata file", line_num+1))?;
bias_fc.push(fc);
}
}
if histograms.len() > 0 {
Some(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(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 {
None
bail!("Histogram has no datapoints.")
}
}
@@ -124,11 +110,9 @@ fn is_in_hist_boundaries(values: &Vec<f64>, cfg: &Config) -> bool {
}
// parse a timeseries file into a histogram
fn read_window_file(window_file: &str, cfg: &Config) -> Option<Histogram> {
let f = File::open(window_file).unwrap_or_else(|x| {
eprintln!("Failed to read sample data from {}. {}", window_file, x);
process::exit(1)
});
fn read_window_file(window_file: &str, cfg: &Config) -> Result<Histogram> {
let f = File::open(window_file)
.chain_err(|| format!("Failed to open sample data file {}", window_file))?;
let mut buf = BufReader::new(&f);
// total number of bins is the product of all dimensions length
@@ -142,7 +126,7 @@ fn read_window_file(window_file: &str, cfg: &Config) -> Option<Histogram> {
// read and parse each timeseries line
let mut line = String::new();
while buf.read_line(&mut line).unwrap() > 0 {
while buf.read_line(&mut line).chain_err(|| "Failed to read line")? > 0 {
// skip comments and empty lines
if line.starts_with("#") || line.starts_with("@") || line.len() == 0 {
line.clear();
@@ -170,24 +154,24 @@ fn read_window_file(window_file: &str, cfg: &Config) -> Option<Histogram> {
}
let num_points: f64 = hist.iter().sum();
if num_points == 0.0 {
return None
}
Some(Histogram::new(num_points as u32, hist))
Ok(Histogram::new(num_points as u32, hist))
}
pub fn write_results(out_file: &str, ds: &Dataset, free: &Vec<f64>, prob: &Vec<f64>) -> Result<(), Box<Error>> {
let output = File::create(out_file)?;
pub fn write_results(out_file: &str, ds: &Dataset, free: &Vec<f64>, prob: &Vec<f64>) -> Result<()> {
let output = File::create(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(" ");
let header: String = (0..ds.dimens_lengths.len()).map(|d| format!("coord{}", d+1))
.collect::<Vec<String>>().join(" ");
writeln!(buf, "#{} {} {}", header, "Free Energy", "Probability");
for bin in 0..free.len() {
let coords = ds.get_coords_for_bin(bin);
let coords_str: String = coords.iter().map(|c| {format!("{:8.6} ", c)})
.collect::<Vec<String>>().join("\t");
writeln!(buf, "{}{:8.6} {:8.6}", coords_str, free[bin], prob[bin])?;
writeln!(buf, "{}{:8.6} {:8.6}", coords_str, free[bin], prob[bin])
.chain_err(|| "Failed to write to file.")?;
}
Ok(())
}
@@ -231,9 +215,7 @@ mod tests {
#[test]
fn read_data() {
let cfg = cfg();
let ds = super::read_data(&cfg);
assert!(ds.is_some());
let ds = ds.unwrap();
let ds = super::read_data(&cfg).unwrap();
println!("{:?}", ds);
assert_eq!(25, ds.num_windows);
assert_eq!(cfg.num_bins.len(), ds.dimens_lengths.len());

17
src/lib.rs
View File

@@ -1,15 +1,20 @@
#![allow(non_snake_case)]
#[macro_use]
extern crate error_chain;
pub mod io;
pub mod histogram;
use std::error::Error;
use std::result::Result;
use histogram::Dataset;
use std::f64;
use std::fmt;
use std::io::prelude::*;
// init error chain
pub mod errors { error_chain!{} }
use errors::*;
#[allow(non_upper_case_globals)]
static k_B: f64 = 0.0083144621; // kJ/mol*K
@@ -84,13 +89,12 @@ fn perform_wham_iteration(ds: &Dataset, F_prev: &[f64], F: &mut [f64], P: &mut [
}
}
pub fn run(cfg: &Config) -> Result<(), Box<Error>>{
pub fn run(cfg: &Config) -> Result<()>{
println!("Supplied WHAM options: {}", &cfg);
println!("Reading input files.");
// TODO Better error handling with nice error messages instead of a panic!
let histograms = io::read_data(&cfg)
.expect("No datapoints in histogram boundaries.");
let histograms = io::read_data(&cfg).chain_err(|| "Failed to create histogram.")?;
println!("{}",&histograms);
// allocate required vectors.
@@ -146,7 +150,8 @@ pub fn run(cfg: &Config) -> Result<(), Box<Error>>{
println!("!!!!! WHAM not converged! (max iterations reached) !!!!!");
}
io::write_results(&cfg.output, &histograms, &free_energy, &P)?;
io::write_results(&cfg.output, &histograms, &free_energy, &P)
.chain_err(|| "Could not write results to output file")?;
Ok(())
}

25
src/main.rs
View File

@@ -4,19 +4,21 @@ extern crate clap;
use clap::App;
use wham::Config;
use std::error::Error;
use std::result::Result;
use wham::errors::*;
use std::process;
// Parse command line arguments into a Config struct
fn cli() -> Result<Config, Box<Error>> {
fn cli() -> Result<Config> {
let yaml = load_yaml!("cli.yml");
let matches = App::from_yaml(yaml).get_matches();
let metadata_file = matches.value_of("metadata").unwrap().to_string();
let verbose: bool = matches.is_present("verbose");
let temperature: f64 = matches.value_of("temperature").unwrap().parse()?;
let tolerance: f64 = matches.value_of("tolerance").unwrap_or("0.000001").parse()?;
let max_iterations: usize = matches.value_of("iterations").unwrap_or("100000").parse()?;
let temperature: f64 = matches.value_of("temperature").unwrap().parse()
.chain_err(|| "Cannot read temperature.")?;
let tolerance: f64 = matches.value_of("tolerance").unwrap_or("0.000001").parse()
.chain_err(|| "Cannot read tolerance.")?;
let max_iterations: usize = matches.value_of("iterations").unwrap_or("100000").parse()
.chain_err(|| "Cannot parse iterations.")?;
let output = matches.value_of("output").unwrap_or("wham.out").to_string();
let cyclic: bool = matches.is_present("cyclic");
@@ -40,9 +42,14 @@ fn cli() -> Result<Config, Box<Error>> {
}
fn main() {
let cfg = cli().expect("Failed to parse CLI.");
match wham::run(&cfg) {
Err(_) => process::exit(1),
_ => {}
if let Err(error) = wham::run(&cfg) {
eprintln!("Error: {}", error);
for e in error.iter().skip(1) {
eprintln!("Reason: {}", e)
}
process::exit(1);
}
}

49
tests/integration.rs
View File

@@ -16,11 +16,56 @@ mod integration {
assert!(out.contains(expected_error));
}
#[test]
fn calling_wham_with_data_out_of_bonds() {
let output = Command::new("./target/debug/wham")
.args(&["--bins", "100", "--min", "2.0", "--max", "3.0", "-T", "300"])
.args(&["-f", "example/1d/metadata.dat"])
.args(&["-o", "/dev/null"])
.output()
.expect("failed to execute process");
let output = String::from_utf8_lossy(&output.stderr);
println!("{}", output);
assert!(output.to_string().contains(
"No data points in histogram boundaries"
));
}
#[test]
fn calling_wham_with_unparseable_bias_pos() {
let output = Command::new("./target/debug/wham")
.args(&["--bins", "100", "--min", "-3.0", "--max", "3.0", "-T", "300"])
.args(&["-f", "tests/metadata_unparseable1.dat"])
.args(&["-o", "/dev/null"])
.output()
.expect("failed to execute process");
let output = String::from_utf8_lossy(&output.stderr);
println!("{}", output);
assert!(output.to_string().contains(
"Failed to read bias position in line 1"
));
}
#[test]
fn calling_wham_with_unparseable_bias_fc() {
let output = Command::new("./target/debug/wham")
.args(&["--bins", "100", "--min", "-3.0", "--max", "3.0", "-T", "300"])
.args(&["-f", "tests/metadata_unparseable2.dat"])
.args(&["-o", "/dev/null"])
.output()
.expect("failed to execute process");
let output = String::from_utf8_lossy(&output.stderr);
println!("{}", output);
assert!(output.to_string().contains(
"Failed to read bias fc in line 1"
));
}
#[test]
fn calling_wham_1d() {;
Command::new("./target/debug/wham")
.args(&["--bins", "100", "--max", "3.14", "--min", "-3.14", "-T", "300", "--cyclic"])
.args(&["-f", "example/1d/metadata.dat"])
.args(&["-f", "example/1d/metadata_unparseable1.dat"])
.args(&["-o", "/tmp/wham_test_1d.out"])
.output()
.expect("failed to execute process");
@@ -40,7 +85,7 @@ mod integration {
fn calling_wham_2d() {;
Command::new("./target/debug/wham")
.args(&["--bins", "100,100", "--max", "3.14,3.14", "--min", "-3.14,-3.14", "-T", "300", "--cyclic"])
.args(&["-f", "example/2d/metadata.dat"])
.args(&["-f", "example/2d/metadata_unparseable1.dat"])
.args(&["-o", "/tmp/wham_test_2d.out"])
.output()
.expect("failed to execute process");

25
tests/metadata_unparseable1.dat Normal file
View File

@@ -0,0 +1,25 @@
../example/1d/COLVAR+0.0.xvg qwert 100
../example/1d/COLVAR+0.25.xvg 0.25 100
../example/1d/COLVAR+0.5.xvg 0.5 100
../example/1d/COLVAR+0.75.xvg 0.75 100
../example/1d/COLVAR+1.0.xvg 1.0 100
../example/1d/COLVAR+1.25.xvg 1.25 100
../example/1d/COLVAR+1.5.xvg 1.5 100
../example/1d/COLVAR+1.75.xvg 1.75 100
../example/1d/COLVAR+2.0.xvg 2.0 100
../example/1d/COLVAR+2.25.xvg 2.25 100
../example/1d/COLVAR+2.5.xvg 2.5 100
../example/1d/COLVAR+2.75.xvg 2.75 100
../example/1d/COLVAR+3.0.xvg 3.0 100
../example/1d/COLVAR-0.25.xvg -0.25 100
../example/1d/COLVAR-0.5.xvg -0.5 100
../example/1d/COLVAR-0.75.xvg -0.75 100
../example/1d/COLVAR-1.0.xvg -1.0 100
../example/1d/COLVAR-1.25.xvg -1.25 100
../example/1d/COLVAR-1.5.xvg -1.5 100
../example/1d/COLVAR-1.75.xvg -1.75 100
../example/1d/COLVAR-2.0.xvg -2.0 100
../example/1d/COLVAR-2.25.xvg -2.25 100
../example/1d/COLVAR-2.5.xvg -2.5 100
../example/1d/COLVAR-2.75.xvg -2.75 100
../example/1d/COLVAR-3.0.xvg -3.0 100

25
tests/metadata_unparseable2.dat Normal file
View File

@@ -0,0 +1,25 @@
../example/1d/COLVAR+0.0.xvg 0.0 qwert
../example/1d/COLVAR+0.25.xvg 0.25 100
../example/1d/COLVAR+0.5.xvg 0.5 100
../example/1d/COLVAR+0.75.xvg 0.75 100
../example/1d/COLVAR+1.0.xvg 1.0 100
../example/1d/COLVAR+1.25.xvg 1.25 100
../example/1d/COLVAR+1.5.xvg 1.5 100
../example/1d/COLVAR+1.75.xvg 1.75 100
../example/1d/COLVAR+2.0.xvg 2.0 100
../example/1d/COLVAR+2.25.xvg 2.25 100
../example/1d/COLVAR+2.5.xvg 2.5 100
../example/1d/COLVAR+2.75.xvg 2.75 100
../example/1d/COLVAR+3.0.xvg 3.0 100
../example/1d/COLVAR-0.25.xvg -0.25 100
../example/1d/COLVAR-0.5.xvg -0.5 100
../example/1d/COLVAR-0.75.xvg -0.75 100
../example/1d/COLVAR-1.0.xvg -1.0 100
../example/1d/COLVAR-1.25.xvg -1.25 100
../example/1d/COLVAR-1.5.xvg -1.5 100
../example/1d/COLVAR-1.75.xvg -1.75 100
../example/1d/COLVAR-2.0.xvg -2.0 100
../example/1d/COLVAR-2.25.xvg -2.25 100
../example/1d/COLVAR-2.5.xvg -2.5 100
../example/1d/COLVAR-2.75.xvg -2.75 100
../example/1d/COLVAR-3.0.xvg -3.0 100