From c9e18ebda9b02470acb5c47b8b85fdab230b5249 Mon Sep 17 00:00:00 2001 From: Daniel Bauer Date: Sat, 17 Jul 2021 14:14:25 +0200 Subject: [PATCH] Squashed commit of the following: commit eaebf0dcbb259decbc0d8f5bbffa62244303f6c7 Author: Daniel Bauer Date: Sat Jul 17 14:12:30 2021 +0200 error for empty timeseries commit 1fe5383c5079ca6c0ad102d4e2cb32d5b1227e80 Author: Daniel Bauer Date: Sat Jul 17 14:09:13 2021 +0200 refractored convdt slices calculation commit 0b1fb7fb6a72edc50b013b59623551b2ccab6913 Author: Daniel Bauer Date: Sat Jul 17 13:27:59 2021 +0200 fix histogram building without convdt commit f9882ca4cece57451cd2971d0993a2d3621b0cdf Author: Daniel Bauer Date: Sat Jul 17 12:53:06 2021 +0200 fix tests not compiling commit e6550e20bde3824f8199432b08777be41fc79fa8 Author: Daniel Bauer Date: Sat Jul 17 12:43:36 2021 +0200 refractoring commit 18c77a9b6694d491ef23cebb3918ccda8fcc44a5 Author: Daniel Bauer Date: Fri Jul 16 08:16:34 2021 +0200 run and output for multiple datasets commit 069f318f72207c416387007435eeff9867494e7a Author: Daniel Bauer Date: Fri Jul 16 07:58:01 2021 +0200 cleanup io.rs commit 10efa428a0c6bf490c9f2d7a4c1df185de402a11 Author: Daniel Bauer Date: Thu Jul 15 19:27:22 2021 +0200 parse multiple datasets with convdt --- src/cli.yml | 7 +- src/io.rs | 306 +++++++++++++++++++++++++++++++++++++++++----------- src/lib.rs | 54 ++++++---- src/main.rs | 4 +- 4 files changed, 285 insertions(+), 86 deletions(-) diff --git a/src/cli.yml b/src/cli.yml index 9b1f067..407f416 100644 --- a/src/cli.yml +++ b/src/cli.yml @@ -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 diff --git a/src/io.rs b/src/io.rs index 226dcdb..6d4355a 100644 --- a/src/io.rs +++ b/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 { +// 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> { let mut bias_pos: Vec = Vec::new(); let mut bias_fc: Vec = Vec::new(); - let mut histograms: Vec = Vec::new(); + let mut histograms: Vec> = Vec::new(); let mut timeseries_lengths: Vec = Vec::new(); + let mut paths = Vec::new(); let kT = cfg.temperature * k_B; let bin_width: Vec = (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 { 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 { .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(×eries[0], &cfg); + for (idx, interval) in convdt_stops.iter().enumerate() { + // build histogram for slice start.._stop + let (start, stop) = interval; + let timeseries_mask: Vec = (0..timeseries[0].len()).map(|i| { + is_in_time_boundaries(timeseries[0][i], *start, *stop) + }).collect(); + let hist = build_histogram_from_timeseries(×eries, ×eries_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 = (0..num_datasets).map(|idx| { + let mut dataset_histograms: Vec = 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::>>().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 +156,69 @@ pub fn read_data(cfg: &Config) -> Result { let total_n = timeseries_lengths.iter().sum::() as f64; let total_h = histograms.iter().map(|h| h.num_points).sum::() 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; + } + 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], 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 = (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 = 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 = (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 { - 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 = (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>, usize)> { let mut timeseries: Vec> = 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 = 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::>>(); + + 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 = 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 = (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 +351,24 @@ fn uncorrelate(timeseries: Vec>, cfg: &Config) -> Vec> { } // 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) -> 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::>().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(×eries, &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()); @@ -355,13 +474,70 @@ mod tests { #[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 = (0..31).map(|i| i as f64).collect(); + println!("{:?}", timeseries); + cfg.start = 10.0; cfg.end = 20.0; - assert!(super::is_in_time_boundaries(15.0, &cfg)); - assert!(super::is_in_time_boundaries(10.0, &cfg)); - assert!(super::is_in_time_boundaries(20.0, &cfg)); - assert!(!super::is_in_time_boundaries(9.9999999, &cfg)); - assert!(!super::is_in_time_boundaries(20.000001, &cfg)); + cfg.convdt = 10.0; + let test = super::get_convdt_boundaries(×eries, &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(×eries, &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 = (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(×eries, &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(×eries, &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(×eries, &cfg); + println!("{:?}", test); + assert!(test.len() == 1); + assert_approx_eq!(test[0].0, 10.0); + assert_approx_eq!(test[0].1, 20.0); } } \ No newline at end of file diff --git a/src/lib.rs b/src/lib.rs index a604ba4..14d0129 100644 --- a/src/lib.rs +++ b/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.")?; - 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(()) } diff --git a/src/main.rs b/src/main.rs index 6aadf80..94e7a56 100644 --- a/src/main.rs +++ b/src/main.rs @@ -68,10 +68,12 @@ fn cli() -> Result { } 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() {