diff --git a/src/histogram.rs b/src/histogram.rs index cda7847..c0d81a4 100644 --- a/src/histogram.rs +++ b/src/histogram.rs @@ -54,7 +54,10 @@ pub struct Dataset { bias_fc: Vec, // bias value cache - bias: RefCell>> + bias: RefCell>>, + + // sum of bin count for windows + pub bin_count: Vec } impl Dataset { @@ -62,6 +65,9 @@ impl Dataset { pub fn new(num_bins: usize, dimens_lengths: Vec, bin_width: Vec, hist_min: Vec, hist_max: Vec, bias_pos: Vec, bias_fc: Vec, kT: f64, histograms: Vec, cyclic: bool) -> Dataset { let num_windows = histograms.len(); let bias: RefCell>> = 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); + } } \ No newline at end of file diff --git a/src/lib.rs b/src/lib.rs index e800d40..51b3bfb 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -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]; + for (window, h) in ds.histograms.iter().enumerate() { 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.