diff --git a/src/histogram.rs b/src/histogram.rs index e07b92e..a070d0d 100644 --- a/src/histogram.rs +++ b/src/histogram.rs @@ -78,9 +78,7 @@ impl fmt::Display for HistogramSet { for h in &self.histograms { datapoints += h.num_points; } - - write!(f, "{} windows and {} bins ranging from {} to {} (bin width: {}).\nDatapoints: {}", - self.num_windows, self.num_bins, self.hist_min, self.hist_max, self.bin_width, datapoints) + write!(f, "{} windows, {} datapoints", self.num_windows, datapoints) } } diff --git a/src/io.rs b/src/io.rs index a819b9f..83e236a 100644 --- a/src/io.rs +++ b/src/io.rs @@ -18,6 +18,12 @@ fn get_relative_path(path1: &str, path2: &str) -> String { path1.parent().unwrap().join(path2).to_str().unwrap().to_string() } +pub fn vprintln(s: String, verbose: bool) { + if verbose { + println!("{}", s); + } +} + // Read input data into a histogram set by iterating over input files // given in the metadata file pub fn read_data(cfg: &Config) -> Option { @@ -43,9 +49,7 @@ pub fn read_data(cfg: &Config) -> Option { match read_window_file(&path, cfg) { Some(h) => { histograms.push(h); - if cfg.verbose { - println!("File: {}, {} Data points added.", &path, histograms.last().unwrap().num_points); - } + vprintln(format!("{}, {} data points added.", &path, histograms.last().unwrap().num_points), cfg.verbose); }, None => { eprintln!("No data points inside histogram boundaries: {}", &path); diff --git a/src/lib.rs b/src/lib.rs index 6588bbf..995e2be 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -153,13 +153,14 @@ fn free_energy(hs: &HistogramSet, P: &mut Vec, A: &mut Vec) { } pub fn run(cfg: &Config) -> Result<(), Box>{ - println!("{}", &cfg); + println!("Supplied WHAM options: {}", &cfg); // read input data into the histograms object + println!("Reading input files."); let histograms = io::read_data(&cfg).unwrap(); - println!("Input consists of {}",&histograms); + println!("{}",&histograms); - // allocate data arrays for performance + // allocate only once for better performance let mut F_prev = vec![f32::INFINITY; histograms.num_windows]; let mut F = vec![0.0; histograms.num_windows]; let mut P = vec![f32::NAN; histograms.num_bins]; @@ -168,6 +169,9 @@ pub fn run(cfg: &Config) -> Result<(), Box>{ // perform WHAM until convergence let mut iteration = 0; while !is_converged(&F_prev, &F, cfg.tolerance) && iteration < cfg.max_iterations { + use std::thread; + use std::time; + thread::sleep(time::Duration::from_millis(10)); iteration += 1; // store F values before the next iteration F_prev.copy_from_slice(&F[..]); @@ -176,57 +180,42 @@ pub fn run(cfg: &Config) -> Result<(), Box>{ perform_wham_iteration(&histograms, &F_prev, &mut F, &mut P); // output some stats during calculation - if iteration % 100 == 0 { + if iteration % 10 == 0 { println!("Iteration {}: dF={}", &iteration, &diff_avg(&F_prev, &F)); } // Dump free energy and bias offsets - if iteration % 1000 == 0 { + if iteration % 100 == 0 { free_energy(&histograms, &mut P, &mut A); - println!("#PMF"); - println!("#x\t\tFree Energy\t\tP(x)"); - for bin in 0..histograms.num_bins { - let x = get_x_for_bin(bin, histograms.hist_min, histograms.bin_width); - println!("{:9.5}\t{:9.5}\t{:9.5}", x, A[bin], P[bin]); - } - println!("#Bias offsets"); - println!("#Window\t\tF\t\tdF"); - for window in 0..histograms.num_windows { - println!("{}\t{:9.5}\t{:8.8}", window, F[window], (F[window]-F_prev[window]).abs()); - } + dump_state(&histograms, &F, &F_prev, &P, &A); } } + // final free energy calculation and state dump + free_energy(&histograms, &mut P, &mut A); + dump_state(&histograms, &F, &F_prev, &P, &A); + if iteration == cfg.max_iterations { println!("!!!!! WHAM not converged! (max iterations reached) !!!!!"); } - - // println!("#Window\t\t F"); - // for window in 0..histograms.num_windows { - // println!("#{}\t\t{}", window, F[window]); - // } - // let mut free_energy = vec![0.0; histograms.num_bins]; - // for bin in 0..histograms.num_bins { - // let P = calc_bin_probability(bin, &histograms, &F); - // free_energy[bin] = -histograms.kT * P.ln(); - // } - // let mut min = &f32::MAX; - // for bin in 0..histograms.num_bins { - // if &free_energy[bin] < min { - // min = &free_energy[bin] - // } - // } - // println!("#x\t\tFree energy"); - // for bin in 0..histograms.num_bins { - // let x = get_x_for_bin(bin, histograms.hist_min, histograms.bin_width); - // let free_normalized = free_energy[bin] - min; - // println!("{}\t\t{}", x, free_normalized); - // } - Ok(()) } +fn dump_state(hs: &HistogramSet, F: &Vec, F_prev: &Vec, P: &Vec, A: &Vec) { + println!("# PMF"); + println!("#x\t\tFree Energy\t\tP(x)"); + for bin in 0..hs.num_bins { + let x = get_x_for_bin(bin, hs.hist_min, hs.bin_width); + println!("{:9.5}\t{:9.5}\t{:9.5}", x, A[bin], P[bin]); + } + println!("# Bias offsets"); + println!("#Window\t\tF\t\tdF"); + for window in 0..hs.num_windows { + println!("{}\t{:9.5}\t{:8.8}", window, F[window], (F[window]-F_prev[window]).abs()); + } +} + #[cfg(test)] mod tests { use super::histogram::{HistogramSet,Histogram};