diff --git a/src/histogram.rs b/src/histogram.rs index 80d1fe7..128e89f 100644 --- a/src/histogram.rs +++ b/src/histogram.rs @@ -34,9 +34,8 @@ impl Histogram { } // a set of histograms - #[derive(Debug)] -pub struct HistogramSet { +pub struct Dataset { // number of histogram windows (number of simulations) pub num_windows: usize, @@ -68,11 +67,11 @@ pub struct HistogramSet { pub cyclic: bool, } -impl HistogramSet { +impl Dataset { - pub fn new(num_bins: usize, bin_width: f32, hist_min: f32, hist_max: f32, bias_x0: Vec, bias_fc: Vec, kT: f32, histograms: Vec, cyclic: bool) -> HistogramSet { + pub fn new(num_bins: usize, bin_width: f32, hist_min: f32, hist_max: f32, bias_x0: Vec, bias_fc: Vec, kT: f32, histograms: Vec, cyclic: bool) -> Dataset { let num_windows = histograms.len(); - HistogramSet{num_windows, num_bins, bin_width, hist_min, hist_max, bias_x0, bias_fc, kT, histograms, cyclic} + Dataset{num_windows, num_bins, bin_width, hist_min, hist_max, bias_x0, bias_fc, kT, histograms, cyclic} } @@ -98,7 +97,7 @@ impl HistogramSet { } -impl fmt::Display for HistogramSet { +impl fmt::Display for Dataset { fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { let mut datapoints: u32 = 0; for h in &self.histograms { @@ -122,9 +121,9 @@ mod tests { ) } - fn build_hist_set() -> HistogramSet { + fn build_hist_set() -> Dataset { let h = build_hist(); - HistogramSet::new( + Dataset::new( 7, // num bins 1.0, // bin width 0.0, // hist min @@ -152,51 +151,51 @@ mod tests { #[test] fn calc_bias() { - let hs = build_hist_set(); + let ds = build_hist_set(); // 7th element -> x=7.5, x0=7.5 - assert_eq!(0.0, hs.calc_bias(7, 0)); + assert_eq!(0.0, ds.calc_bias(7, 0)); // 8th element -> x=8.5, x0=7.5 - assert_eq!(5.0, hs.calc_bias(8, 0)); + assert_eq!(5.0, ds.calc_bias(8, 0)); // 9th element -> x=9.5, x0=7.5 - assert_eq!(20.0, hs.calc_bias(9, 0)); + assert_eq!(20.0, ds.calc_bias(9, 0)); // 1st element -> x=0.5, x0=7.5. non-cyclic! - assert_eq!(245.0, hs.calc_bias(0, 0)); + assert_eq!(245.0, ds.calc_bias(0, 0)); } #[test] fn calc_bias_offset_cyclic() { - let mut hs = build_hist_set(); - hs.cyclic = true; + let mut ds = build_hist_set(); + ds.cyclic = true; // 7th element -> x=7.5, x0=7.5 - assert_eq!(0.0, hs.calc_bias(7, 0)); + assert_eq!(0.0, ds.calc_bias(7, 0)); // 8th element -> x=8.5, x0=7.5 - assert_eq!(5.0, hs.calc_bias(8, 0)); + assert_eq!(5.0, ds.calc_bias(8, 0)); // 9th element -> x=9.5, x0=7.5 - assert_eq!(20.0, hs.calc_bias(9, 0)); + assert_eq!(20.0, ds.calc_bias(9, 0)); // 1st element -> x=0.5, x0=7.5 // cyclic flag makes bin 0 neighboring bin 9, so the distance is actually 2 - assert_eq!(20.0, hs.calc_bias(0, 0)); + assert_eq!(20.0, ds.calc_bias(0, 0)); } #[test] fn get_x_for_bin() { - let hs = build_hist_set(); + let ds = build_hist_set(); let expected: Vec = vec![0,1,2,3,4,5,6,7,8].iter() .map(|x| *x as f32 + 0.5).collect(); for i in 0..9 { - assert_eq!(expected[i], hs.get_x_for_bin(i)); + assert_eq!(expected[i], ds.get_x_for_bin(i)); } } } \ No newline at end of file diff --git a/src/io.rs b/src/io.rs index 4e6eb67..83d756f 100644 --- a/src/io.rs +++ b/src/io.rs @@ -1,4 +1,4 @@ -use super::histogram::HistogramSet; +use super::histogram::Dataset; use super::histogram::Histogram; use super::Config; use std::fs::File; @@ -26,7 +26,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 { +pub fn read_data(cfg: &Config) -> Option { let mut bias_x0: Vec = Vec::new(); let mut bias_fc: Vec = Vec::new(); let mut histograms: Vec = Vec::new(); @@ -71,7 +71,7 @@ pub fn read_data(cfg: &Config) -> Option { if histograms.len() > 0 { let bin_width = (cfg.hist_max - cfg.hist_min)/(cfg.num_bins as f32); - Some(HistogramSet::new(cfg.num_bins, bin_width, cfg.hist_min, cfg.hist_max, bias_x0, bias_fc, kT, histograms, cfg.cyclic)) + Some(Dataset::new(cfg.num_bins, bin_width, cfg.hist_min, cfg.hist_max, bias_x0, bias_fc, kT, histograms, cfg.cyclic)) } else { None } @@ -170,20 +170,20 @@ mod tests { #[test] fn read_data() { let cfg = cfg(); - let hs = super::read_data(&cfg); - assert!(hs.is_some()); - let hs = hs.unwrap(); - println!("{:?}", hs); - assert_eq!(2, hs.num_windows); - assert_eq!(cfg.num_bins, hs.num_bins); - assert_eq!(cfg.hist_min, hs.hist_min); - assert_eq!(cfg.hist_max, hs.hist_max); + let ds = super::read_data(&cfg); + assert!(ds.is_some()); + let ds = ds.unwrap(); + println!("{:?}", ds); + assert_eq!(2, ds.num_windows); + assert_eq!(cfg.num_bins, ds.num_bins); + assert_eq!(cfg.hist_min, ds.hist_min); + assert_eq!(cfg.hist_max, ds.hist_max); let expected_bin_width = (cfg.hist_max - cfg.hist_min)/cfg.num_bins as f32; - assert_eq!(expected_bin_width, hs.bin_width); - assert_eq!(vec![0.0, 1.0], hs.bias_x0); - assert_eq!(vec![100.0, 200.0], hs.bias_fc); - assert_eq!(cfg.temperature * k_B, hs.kT); - assert_eq!(2, hs.histograms.len()) + assert_eq!(expected_bin_width, ds.bin_width); + assert_eq!(vec![0.0, 1.0], ds.bias_x0); + assert_eq!(vec![100.0, 200.0], ds.bias_fc); + assert_eq!(cfg.temperature * k_B, ds.kT); + assert_eq!(2, ds.histograms.len()) } #[test] diff --git a/src/lib.rs b/src/lib.rs index 2f17862..c33bde1 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -8,7 +8,7 @@ pub mod histogram; use std::error::Error; use std::result::Result; -use histogram::{HistogramSet,Histogram}; +use histogram::{Dataset,Histogram}; use std::f32; use std::fmt; @@ -52,16 +52,16 @@ fn is_converged(old_F: &Vec, new_F: &Vec, tolerance: f32) -> bool { // estimate the probability of a bin of the histogram set based on F values // This evaluates the first WHAM equation for each bin -fn calc_bin_probability(bin: usize, hs: &HistogramSet, F: &Vec) -> f32 { +fn calc_bin_probability(bin: usize, ds: &Dataset, F: &Vec) -> f32 { let mut denom_sum = 0.0; let mut bin_count = 0.0; - for window in 0..hs.num_windows { - let h: &Histogram = &hs.histograms[window]; + for window in 0..ds.num_windows { + let h: &Histogram = &ds.histograms[window]; if let Some(count) = h.get_bin_count(bin) { bin_count += count; } - let bias = hs.calc_bias(bin, window); - let bias_offset = ((F[window] - bias) / hs.kT).exp(); + let bias = ds.calc_bias(bin, window); + let bias_offset = ((F[window] - bias) / ds.kT).exp(); denom_sum += (h.num_points as f32) * bias_offset; } bin_count / denom_sum @@ -69,77 +69,77 @@ fn calc_bin_probability(bin: usize, hs: &HistogramSet, F: &Vec) -> f32 { // estimate the bias offset F of the histogram based on given probabilities // This evaluates the second WHAM equation for each window -fn calc_window_F(window: usize, hs: &HistogramSet, P: &Vec) -> f32 { +fn calc_window_F(window: usize, ds: &Dataset, P: &Vec) -> f32 { let mut ln_sum = 0.0; - for bin in 0..hs.num_bins { - let bias = hs.calc_bias(bin, window); - ln_sum += P[bin] * (-bias/hs.kT).exp() + for bin in 0..ds.num_bins { + let bias = ds.calc_bias(bin, window); + ln_sum += P[bin] * (-bias/ds.kT).exp() } - -hs.kT * ln_sum.ln() + -ds.kT * ln_sum.ln() } // One full WHAM iteration includes calculation of new probabilities P and // new bias offsets F based on previous bias offsets F_prev. This updates // the values in vectors F and P -fn perform_wham_iteration(hs: &HistogramSet, F_prev: &Vec,F: &mut Vec, P: &mut Vec) { +fn perform_wham_iteration(ds: &Dataset, F_prev: &Vec,F: &mut Vec, P: &mut Vec) { // reset bias offsets - for window in 0..hs.num_windows { + for window in 0..ds.num_windows { F[window] = 0.0; } - // for bin in 0..hs.num_bins { - // let x = get_x_for_bin(bin, hs.hist_min, hs.bin_width); + // for bin in 0..ds.num_bins { + // let x = get_x_for_bin(bin, ds.hist_min, ds.bin_width); // let mut num = 0.0; // let mut denom = 0.0; - // for window in 0..hs.num_windows { - // match hs.histograms[window].get_bin_count(bin) { + // for window in 0..ds.num_windows { + // match ds.histograms[window].get_bin_count(bin) { // Some(c) => num += c, // _ => {} // } // let bias = calc_bias( - // hs.bias_fc[window], - // hs.bias_x0[window], + // ds.bias_fc[window], + // ds.bias_x0[window], // x); - // let bf = ((F_prev[window]-bias) / hs.kT).exp(); - // denom += hs.histograms[window].num_points as f32* bf + // let bf = ((F_prev[window]-bias) / ds.kT).exp(); + // denom += ds.histograms[window].num_points as f32* bf // } // P[bin] = num / denom; - // for window in 0..hs.num_windows { + // for window in 0..ds.num_windows { // let bias = calc_bias( - // hs.bias_fc[window], - // hs.bias_x0[window], + // ds.bias_fc[window], + // ds.bias_x0[window], // x); - // let bf = (-bias/hs.kT).exp() * P[bin]; + // let bf = (-bias/ds.kT).exp() * P[bin]; // F[window] += bf; // } // } - // for window in 0..hs.num_windows { - // F[window] = -hs.kT * F[window].ln(); + // for window in 0..ds.num_windows { + // F[window] = -ds.kT * F[window].ln(); // } // let norm = F[0]; - // for window in 0..hs.num_windows { + // for window in 0..ds.num_windows { // F[window] = F[window] - norm; // } // evaluate first WHAM equation for each bin to // estimage probabilities based on previous offsets (F_prev) - for bin in 0..hs.num_bins { - P[bin] = calc_bin_probability(bin, hs, F_prev); + for bin in 0..ds.num_bins { + P[bin] = calc_bin_probability(bin, ds, F_prev); } // evaluate second WHAM equation for each window to // estimate new bias offsets from propabilities - for window in 0..hs.num_windows { - F[window] = calc_window_F(window, hs, P); + for window in 0..ds.num_windows { + F[window] = calc_window_F(window, ds, P); } // normalize F let norm = F[0]; - for window in 0..hs.num_windows { + for window in 0..ds.num_windows { F[window] = F[window] - norm; } } @@ -155,27 +155,27 @@ fn diff_avg(F: &Vec, F_prev: &Vec) -> f32 { // calculate the normalized free energy from normalized probability values -fn free_energy(hs: &HistogramSet, P: &mut Vec, A: &mut Vec) { +fn free_energy(ds: &Dataset, P: &mut Vec, A: &mut Vec) { let mut bin_min = f32::MAX; // Free energy calculation - for bin in 0..hs.num_bins { - A[bin] = -hs.kT*P[bin].ln(); + for bin in 0..ds.num_bins { + A[bin] = -ds.kT*P[bin].ln(); if A[bin] < bin_min { bin_min = A[bin]; } } // Make A relative to minimum - for bin in 0..hs.num_bins { + for bin in 0..ds.num_bins { A[bin] -= bin_min; } // Normalize P // let mut P_sum = 0.0; - // for bin in 0..hs.num_bins { + // for bin in 0..ds.num_bins { // P_sum += P[bin]; // } - // for bin in 0..hs.num_bins { + // for bin in 0..ds.num_bins { // P[bin] /= P_sum; // } @@ -231,23 +231,23 @@ pub fn run(cfg: &Config) -> Result<(), Box>{ Ok(()) } -fn dump_state(hs: &HistogramSet, F: &Vec, F_prev: &Vec, P: &Vec, A: &Vec) { +fn dump_state(ds: &Dataset, 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 = hs.get_x_for_bin(bin); + for bin in 0..ds.num_bins { + let x = ds.get_x_for_bin(bin); 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 { + for window in 0..ds.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}; + use super::histogram::{Dataset,Histogram}; use std::f32; #[test] @@ -263,10 +263,10 @@ mod tests { assert!(!converged); } - fn create_test_hs() -> HistogramSet { + fn create_test_ds() -> Dataset { let h1 = Histogram::new(0, 2, 10, vec![3.0, 4.0, 3.0]); let h2 = Histogram::new(0, 3, 20, vec![3.0, 2.0, 5.0, 10.0]); - HistogramSet::new(4, 1.0, 0.0, 4.0, vec![1.0, 2.0], vec![10.0, 10.0], 2.479, vec![h1, h2], false) + Dataset::new(4, 1.0, 0.0, 4.0, vec![1.0, 2.0], vec![10.0, 10.0], 2.479, vec![h1, h2], false) } fn assert_near(a: f32, b: f32, tolerance: f32) { @@ -276,46 +276,46 @@ mod tests { #[test] fn calc_bias_offset() { - let hs = create_test_hs(); + let ds = create_test_ds(); let probability = vec!(0.959, 0.331, 0.656, 46.750); let expected = vec!(0.596, -0.250); - for window in 0..hs.num_windows { - let F = super::calc_window_F(window, &hs, &probability); + for window in 0..ds.num_windows { + let F = super::calc_window_F(window, &ds, &probability); assert_near(expected[window], F, 0.001); } } #[test] fn calc_bin_probability() { - let hs = create_test_hs(); + let ds = create_test_ds(); let F = vec!(0.0, 0.0); let expected = vec!(0.959, 0.331, 0.656, 46.750); for b in 0..4 { - let p = super::calc_bin_probability(b, &hs, &F); + let p = super::calc_bin_probability(b, &ds, &F); assert_near(expected[b], p, 0.001); } let F = vec!(1.0, 1.0); let expected = vec!(0.641, 0.221, 0.439, 31.232); for b in 0..4 { - let p = super::calc_bin_probability(b, &hs, &F); + let p = super::calc_bin_probability(b, &ds, &F); assert_near(expected[b], p, 0.001); } } #[test] fn perform_wham_iteration() { - let hs = create_test_hs(); - let prev_F = vec![0.0; hs.num_windows]; - let mut F = vec![0.0; hs.num_windows]; - let mut P = vec![f32::NAN; hs.num_bins]; - super::perform_wham_iteration(&hs, &prev_F, &mut F, &mut P); + let ds = create_test_ds(); + let prev_F = vec![0.0; ds.num_windows]; + let mut F = vec![0.0; ds.num_windows]; + let mut P = vec![f32::NAN; ds.num_bins]; + super::perform_wham_iteration(&ds, &prev_F, &mut F, &mut P); let expected_F = vec!(0.0, -0.846); let expected_P = vec!(0.959, 0.331, 0.656, 46.750); - for bin in 0..hs.num_bins { + for bin in 0..ds.num_bins { assert_near(expected_P[bin], P[bin], 0.01) } - for window in 0..hs.num_windows { + for window in 0..ds.num_windows { assert_near(expected_F[window], F[window], 0.01) }