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

dump state periodically

Browse Source
This commit is contained in:
Daniel Bauer
2018-10-05 21:09:50 +02:00
parent 530ae7b7be
commit 1ec2f3ac45
3 changed files with 36 additions and 45 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
src/histogram.rs
View File

@@ -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)
}
}

10
src/io.rs
View File

@@ -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<HistogramSet> {
@@ -43,9 +49,7 @@ pub fn read_data(cfg: &Config) -> Option<HistogramSet> {
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);

67
src/lib.rs
View File

@@ -153,13 +153,14 @@ fn free_energy(hs: &HistogramSet, P: &mut Vec<f32>, A: &mut Vec<f32>) {
}
pub fn run(cfg: &Config) -> Result<(), Box<Error>>{
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<Error>>{
// 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<Error>>{
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<f32>, F_prev: &Vec<f32>, P: &Vec<f32>, A: &Vec<f32>) {
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};