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

calculate total bin count during i/o

Browse Source
This commit is contained in:
Daniel Bauer
2018-10-22 10:01:34 +02:00
parent 5842a33f44
commit b40ad9b516
2 changed files with 45 additions and 20 deletions
Download Patch File Download Diff File Expand all files Collapse all files

58
src/histogram.rs
View File

@@ -54,7 +54,10 @@ pub struct Dataset {
bias_fc: Vec<f64>,
// bias value cache
bias: RefCell<Vec<Option<f64>>>
bias: RefCell<Vec<Option<f64>>>,
// sum of bin count for windows
pub bin_count: Vec<f64>
}
impl Dataset {
@@ -62,6 +65,9 @@ impl Dataset {
pub fn new(num_bins: usize, dimens_lengths: Vec<usize>, bin_width: Vec<f64>, hist_min: Vec<f64>, hist_max: Vec<f64>, bias_pos: Vec<f64>, bias_fc: Vec<f64>, kT: f64, histograms: Vec<Histogram>, cyclic: bool) -> Dataset {
let num_windows = histograms.len();
let bias: RefCell<Vec<Option<f64>>> = RefCell::new(vec![None; num_bins*num_windows]);
let bin_count = (0..num_bins).map(|bin| {
histograms.iter().map(|h| h.bins[bin]).sum()
}).collect();
Dataset{
num_windows,
num_bins,
@@ -75,6 +81,7 @@ impl Dataset {
bias_pos,
bias_fc,
bias,
bin_count,
}
}
@@ -166,12 +173,12 @@ mod tests {
fn build_hist_set() -> Dataset {
let h = build_hist();
Dataset::new(
9, // num bins
5, // num bins
vec![1],
vec![1.0], // bin width
vec![0.0], // hist min
vec![9.0], // hist max
vec![7.5], // x0
vec![4.5], // x0
vec![10.0], // fc
300.0*k_B, // kT
vec![h], // hists
@@ -183,13 +190,13 @@ mod tests {
fn calc_bias() {
let ds = build_hist_set(); // k = 10
// 7th element -> x=7.5, x0=7.5
assert_delta!(1.0, ds.calc_bias(7, 0), 0.00000001);
// 3th element -> x=3.5, x0=3.5
assert_delta!(0.134722337796, ds.calc_bias(3, 0), 0.00000001);
// 8th element -> x=8.5, x0=7.5
assert_delta!(0.13472233779, ds.calc_bias(8,0), 0.00000001);
// 8th element -> x=8.5, x0=3.5
assert_delta!(1.0, ds.calc_bias(4,0), 0.00000001);
// 1st element -> x=0.5, x0=7.5. non-cyclic!
// 1st element -> x=0.5, x0=3.5. non-cyclic!
assert_delta!(0.0, ds.calc_bias(0,0), 0.0000001);
}
@@ -198,18 +205,18 @@ mod tests {
let mut ds = build_hist_set();
ds.cyclic = true;
// 7th element -> x=7.5, x0=7.5
assert_delta!(1.0, ds.calc_bias(7, 0), 0.00000001);
// 7th element -> x=3.5, x0=3.5
assert_delta!(0.134722337796, ds.calc_bias(3, 0), 0.00000001);
// 8th element -> x=8.5, x0=7.5
assert_delta!(0.13472233779, ds.calc_bias(8, 0), 0.00000001);
// 8th element -> x=4.5, x0=3.5
assert_delta!(1.0, ds.calc_bias(4, 0), 0.00000001);
// 1th element -> x=0.5, x0=7.5
// 1th element -> x=0.5, x0=3.5
// cyclic flag makes bin 0 neighboring bin 9, so the distance is actually 2
assert_delta!(0.00032942643, ds.calc_bias(0, 0), 0.00000001);
assert_delta!(0.0000000000000117769, ds.calc_bias(0, 0), 0.00000001);
// 2nd element -> x=1.5, x0=7.5
// 2nd element -> x=1.5, x0=3.5
assert_delta!(0.00000001, ds.calc_bias(1, 0), 0.00000001);
}
@@ -222,4 +229,25 @@ mod tests {
assert_eq!(expected[i], ds.get_coords_for_bin(i)[0]);
}
}
#[test]
fn get_bin_count() {
let ds = Dataset::new(
5, // num bins
vec![1],
vec![1.0], // bin width
vec![0.0], // hist min
vec![9.0], // hist max
vec![7.5], // x0
vec![10.0], // fc
300.0*k_B, // kT
vec![build_hist(), build_hist()], // hists
false // cyclic
);
assert_delta!(2.0, ds.bin_count[0], 0.0000000001);
assert_delta!(2.0, ds.bin_count[1], 0.0000000001);
assert_delta!(6.0, ds.bin_count[2], 0.0000000001);
assert_delta!(10.0, ds.bin_count[3], 0.0000000001);
assert_delta!(24.0, ds.bin_count[4], 0.0000000001);
}
}

5
src/lib.rs
View File

@@ -49,14 +49,11 @@ fn is_converged(old_F: &[f64], new_F: &[f64], tolerance: f64) -> bool {
// This evaluates the first WHAM equation for each bin.
fn calc_bin_probability(bin: usize, ds: &Dataset, F: &[f64]) -> f64 {
let mut denom_sum: f64 = 0.0;
let mut bin_count: f64 = 0.0;
// TODO calculate bin_count before wham iterations for performance
for (window, h) in ds.histograms.iter().enumerate() {
bin_count += h.bins[bin];
let bias = ds.calc_bias(bin, window);
denom_sum += (h.num_points as f64) * bias * F[window];
}
bin_count / denom_sum
ds.bin_count[bin] / denom_sum
}
// estimate the bias offset F of the histogram based on given probabilities.