From b064ecea6d37578c3aa8aacf6d8ed3466adc5963 Mon Sep 17 00:00:00 2001
From: Daniel Bauer <bauer@cbs.tu-darmstadt.de>
Date: Wed, 17 Oct 2018 10:21:19 +0200
Subject: [PATCH] calc exp(F/kT) instead of F for each window

---
 src/lib.rs | 25 +++++++++++++++----------
 1 file changed, 15 insertions(+), 10 deletions(-)

diff --git a/src/lib.rs b/src/lib.rs
index ac2e63f..8400bfa 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -64,15 +64,14 @@ fn calc_bin_probability(bin: usize, ds: &Dataset, F: &[f64]) -> f64 {
 // 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, ds: &Dataset, P: &[f64]) -> f64 {
-    let bf_sum: f64 = (0..ds.num_bins).zip(P.iter()) // zip bins and P
-		.filter_map(|count_and_prob: (usize, &f64)| {
-            if count_and_prob.1 == &0.0 { // skip zeros for speed
+    (0..ds.num_bins).zip(P.iter()) // zip bins and P
+		.filter_map(|bin_and_prob: (usize, &f64)| {
+            if bin_and_prob.1 == &0.0 { // skip zeros for speed
                 None
             } else {
-                Some(count_and_prob.1 * (-(-ds.kT*ds.calc_bias(count_and_prob.0, window).ln()) / ds.kT).exp())
+                Some(bin_and_prob.1 * ds.calc_bias(bin_and_prob.0, window))
             }
-		}).sum();
-	-ds.kT * bf_sum.ln()
+		}).sum()
 }
 
 // One full WHAM iteration includes calculation of new probabilities P and
@@ -231,11 +230,18 @@ fn dump_state(ds: &Dataset, F: &[f64], F_prev: &[f64], P: &[f64], A: &[f64]) {
 	}
 }
 
+
 #[cfg(test)]
 mod tests {
 	use super::histogram::{Dataset,Histogram};
 	use std::f64;
 
+    macro_rules! assert_delta {
+        ($x:expr, $y:expr, $d:expr) => {
+            assert!(($x-$y).abs() < $d, "{} != {}", $x, $y)
+        }
+    }
+
 	#[test]
 	fn is_converged() {
 		let new = vec![1.0,1.0];
@@ -262,15 +268,14 @@ mod tests {
 	}
 
 	#[test]
-    #[ignore]
 	fn calc_bias_offset() {
 		let ds = create_test_ds();
 		let probability = vec!(0.959, 0.331, 0.656, 46.750);
-		let expected = vec!(0.596, -0.250);
+        let expected = vec!(0.786289183, 1.10629119);
 		for window in 0..ds.num_windows {
 			let F = super::calc_window_F(window, &ds, &probability);
-			assert_near(expected[window], F, 0.001);
-		}
+            assert_delta!(expected[window], F, 0.0000001);
+        }
 	}
 
 	#[test]