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

commenting

Browse Source
This commit is contained in:
Daniel Bauer
2018-11-13 09:32:20 +01:00
parent cffa9158d0
commit fe5ff0a9af
1 changed files with 15 additions and 5 deletions
Download Patch File Download Diff File Expand all files Collapse all files

20
src/error_analysis.rs
View File

@@ -4,8 +4,10 @@ use super::perform_wham;
use super::Config; use super::Config;
use rgsl::statistics; use rgsl::statistics;
// returns a set of num_windows continious weights by
// a) generate num_windows-1 random variables and sort them
// b) each weight n is the difference between n+1 and n, where n0=0 and nN+1=1
fn generate_random_weights(num_windows: usize) -> Vec<f64> { fn generate_random_weights(num_windows: usize) -> Vec<f64> {
// TODO this is ugly.
// create a list of num_windows - 1 sorted random numbers and append/prepend 0 and 1 // create a list of num_windows - 1 sorted random numbers and append/prepend 0 and 1
let mut tmp = (0..num_windows-1).map(|_| rand::random::<f64>()).collect::<Vec<f64>>(); let mut tmp = (0..num_windows-1).map(|_| rand::random::<f64>()).collect::<Vec<f64>>();
tmp.sort_by(|a,b| { a.partial_cmp(b).unwrap() }); tmp.sort_by(|a,b| { a.partial_cmp(b).unwrap() });
@@ -21,23 +23,31 @@ fn generate_random_weights(num_windows: usize) -> Vec<f64> {
return weights return weights
} }
// Generate a random weighted dataset from the given dataset by changing the weights
fn generate_random_weighted_dataset(ds: Dataset) -> Dataset { fn generate_random_weighted_dataset(ds: Dataset) -> Dataset {
let weights = generate_random_weights(ds.num_windows); let weights = generate_random_weights(ds.num_windows);
Dataset::new_weighted(ds, weights ) Dataset::new_weighted(ds, weights)
} }
// Perform bootstrap error analysis. This runs the WHAM analysis num_runs times on random weighted
// datasets. The standard deviation is calculated on the bootstrapped probabilities of each bin. The
// standard deviation of the free eneergy is then deduced by error propagation (A_std = kT*1/P*P_std)
pub fn run_bootstrap(cfg: &Config, ds: Dataset, P: &[f64], num_runs: usize) -> (Vec<f64>,Vec<f64>) { pub fn run_bootstrap(cfg: &Config, ds: Dataset, P: &[f64], num_runs: usize) -> (Vec<f64>,Vec<f64>) {
let Ps: Vec<Vec<f64>> = (0..num_runs).map(|x| { // Calculate bootstrapped probabilities
let bootstrapped_Ps: Vec<Vec<f64>> = (0..num_runs).map(|x| {
println!("Bootstrap run {}/{}", x, num_runs); println!("Bootstrap run {}/{}", x, num_runs);
let rnd_weighted_dataset = generate_random_weighted_dataset(ds.clone()); let rnd_weighted_dataset = generate_random_weighted_dataset(ds.clone());
perform_wham(cfg, &rnd_weighted_dataset).unwrap().0 perform_wham(cfg, &rnd_weighted_dataset).unwrap().0
}).collect(); }).collect();
// Evaulate standard deviation of P per bin
let mut P_std = vec![0.0; ds.num_bins]; let mut P_std = vec![0.0; ds.num_bins];
for bin in 0..ds.num_bins { for bin in 0..ds.num_bins {
let Ps = Ps.iter().map(|window| window[bin]).collect::<Vec<f64>>(); let Ps = bootstrapped_Ps.iter().map(|window| window[bin]).collect::<Vec<f64>>();
P_std[bin] = statistics::sd(&Ps, 1, num_runs); P_std[bin] = statistics::sd(&Ps, 1, num_runs);
} }
// A_std by error propagation
let A_std = P_std.iter().zip(P.iter()).map(|(std,P)| ds.kT*1.0/P*std).collect(); let A_std = P_std.iter().zip(P.iter()).map(|(std,P)| ds.kT*1.0/P*std).collect();
(P_std, A_std) (P_std, A_std)
} }