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

cache exp(-U/kT) instead of U

Browse Source
This commit is contained in:
Daniel Bauer
2018-10-17 09:51:11 +02:00
parent 7fd78fc2a5
commit 40bb99a612
2 changed files with 22 additions and 12 deletions
Download Patch File Download Diff File Expand all files Collapse all files

29
src/histogram.rs
View File

@@ -126,6 +126,7 @@ impl Dataset {
}
bias_sum += 0.5 * bias_fc[i] * dist * dist
}
let bias_sum = (-bias_sum/self.kT).exp();
cache[ndx] = Some(bias_sum);
bias_sum
}
@@ -149,6 +150,12 @@ mod tests {
use super::*;
use super::super::k_B;
macro_rules! assert_delta {
($x:expr, $y:expr, $d:expr) => {
assert!(($x-$y).abs() < $d, "{} != {}", $x, $y)
}
}
fn build_hist() -> Histogram {
Histogram::new(
22, // num_points
@@ -158,7 +165,7 @@ mod tests {
fn build_hist_set() -> Dataset {
let h = build_hist();
Dataset::new(
Dataset::new(
9, // num bins
vec![1],
vec![1.0], // bin width
@@ -186,38 +193,40 @@ mod tests {
}
}
#[test]
fn calc_bias() {
let ds = build_hist_set();
let ds = build_hist_set(); // k = 10
// 7th element -> x=7.5, x0=7.5
assert_eq!(0.0, ds.calc_bias(7, 0));
assert_delta!(1.0, ds.calc_bias(7, 0), 0.00000001);
// 8th element -> x=8.5, x0=7.5
assert_eq!(5.0, ds.calc_bias(8, 0));
assert_delta!(0.13472233779, ds.calc_bias(8,0), 0.00000001);
// 1st element -> x=0.5, x0=7.5. non-cyclic!
assert_eq!(245.0, ds.calc_bias(0, 0));
assert_delta!(0.0, ds.calc_bias(0,0), 0.0000001);
}
#[test]
fn calc_bias_offset_cyclic() {
fn calc_biascyclic() {
let mut ds = build_hist_set();
ds.cyclic = true;
// 7th element -> x=7.5, x0=7.5
assert_eq!(0.0, ds.calc_bias(7, 0));
assert_delta!(1.0, ds.calc_bias(7, 0), 0.00000001);
// 8th element -> x=8.5, x0=7.5
assert_eq!(5.0, ds.calc_bias(8, 0));
assert_delta!(0.13472233779, ds.calc_bias(8, 0), 0.00000001);
// 1th 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, ds.calc_bias(0, 0));
assert_delta!(0.00032942643, ds.calc_bias(0, 0), 0.00000001);
// 2nd element -> x=1.5, x0=7.5
assert_eq!(45.0, ds.calc_bias(1, 0));
assert_delta!(0.00000001, ds.calc_bias(1, 0), 0.00000001);
}
#[test]

5
src/lib.rs
View File

@@ -53,7 +53,7 @@ fn calc_bin_probability(bin: usize, ds: &Dataset, F: &[f64]) -> f64 {
let mut bin_count: f64 = 0.0;
for (window, h) in ds.histograms.iter().enumerate() {
bin_count += h.bins[bin];
let bias = ds.calc_bias(bin, window);
let bias = -ds.kT*ds.calc_bias(bin, window).ln();
let bias_offset = ((F[window] - bias) / ds.kT).exp();
denom_sum += (h.num_points as f64) * bias_offset;
}
@@ -69,7 +69,7 @@ fn calc_window_F(window: usize, ds: &Dataset, P: &[f64]) -> f64 {
if count_and_prob.1 == &0.0 { // skip zeros for speed
None
} else {
Some(count_and_prob.1 * (-ds.calc_bias(count_and_prob.0, window) / ds.kT).exp())
Some(count_and_prob.1 * (-(-ds.kT*ds.calc_bias(count_and_prob.0, window).ln()) / ds.kT).exp())
}
}).sum();
-ds.kT * bf_sum.ln()
@@ -262,6 +262,7 @@ mod tests {
}
#[test]
#[ignore]
fn calc_bias_offset() {
let ds = create_test_ds();
let probability = vec!(0.959, 0.331, 0.656, 46.750);