From b6881a6db1a59cb499530c1baec74b25c3b2807b Mon Sep 17 00:00:00 2001 From: Daniel Bauer Date: Thu, 22 Jul 2021 15:48:14 +0200 Subject: [PATCH] code formatting --- src/error_analysis.rs | 44 ++--- src/histogram.rs | 396 +++++++++++++++++++++--------------------- src/io.rs | 14 +- src/lib.rs | 186 ++++++++++---------- src/main.rs | 54 +++--- 5 files changed, 347 insertions(+), 347 deletions(-) diff --git a/src/error_analysis.rs b/src/error_analysis.rs index 0c2a3fd..9393830 100644 --- a/src/error_analysis.rs +++ b/src/error_analysis.rs @@ -73,29 +73,29 @@ mod tests { use super::super::histogram::Histogram; fn build_hist() -> Histogram { - Histogram::new( - 22, // num_points - vec![1.0, 1.0, 3.0, 5.0, 12.0] // bins - ) - } + Histogram::new( + 22, // num_points + vec![1.0, 1.0, 3.0, 5.0, 12.0] // bins + ) + } - fn build_hist_set() -> Dataset { - let h1 = build_hist(); - let h2 = build_hist(); - let h3 = build_hist(); - Dataset::new( - 5, // num bins - vec![3], - vec![1.0], // bin width - vec![0.0], // hist min - vec![9.0], // hist max - vec![4.5, 4.5, 4.5], // x0 - vec![10.0, 10.0, 10.0], // fc - 300.0*k_B, // kT - vec![h1, h2, h3], // hists - false // cyclic - ) - } + fn build_hist_set() -> Dataset { + let h1 = build_hist(); + let h2 = build_hist(); + let h3 = build_hist(); + Dataset::new( + 5, // num bins + vec![3], + vec![1.0], // bin width + vec![0.0], // hist min + vec![9.0], // hist max + vec![4.5, 4.5, 4.5], // x0 + vec![10.0, 10.0, 10.0], // fc + 300.0*k_B, // kT + vec![h1, h2, h3], // hists + false // cyclic + ) + } #[test] fn random_weights() { diff --git a/src/histogram.rs b/src/histogram.rs index cf32eb5..6839205 100644 --- a/src/histogram.rs +++ b/src/histogram.rs @@ -3,266 +3,266 @@ use std::fmt; // One histogram #[derive(Debug,Clone)] pub struct Histogram { - // total number of data points stored in the histogram - pub num_points: u32, + // total number of data points stored in the histogram + pub num_points: u32, - // histogram bins - pub bins: Vec + // histogram bins + pub bins: Vec } impl Histogram { - pub fn new(num_points: u32, bins: Vec) -> Histogram { - Histogram {num_points, bins} - } + pub fn new(num_points: u32, bins: Vec) -> Histogram { + Histogram {num_points, bins} + } } // a set of histograms #[derive(Debug,Clone)] pub struct Dataset { - // number of histogram windows (number of simulations) - pub num_windows: usize, + // number of histogram windows (number of simulations) + pub num_windows: usize, - // total number of bins - pub num_bins: usize, + // total number of bins + pub num_bins: usize, - // number of bins in each dimension - pub dimens_lengths: Vec, + // number of bins in each dimension + pub dimens_lengths: Vec, - // min values of the histogram in each dimension - hist_min: Vec, + // min values of the histogram in each dimension + hist_min: Vec, - // max values of the histogram in each dimension - hist_max: Vec, + // max values of the histogram in each dimension + hist_max: Vec, - // width of a bin in unit of its dimension - bin_width: Vec, + // width of a bin in unit of its dimension + bin_width: Vec, - // value of kT - pub kT: f64, + // value of kT + pub kT: f64, - // histogram for each window - pub histograms: Vec, + // histogram for each window + pub histograms: Vec, - // flag for cyclic reaction coordinates - pub cyclic: bool, + // flag for cyclic reaction coordinates + pub cyclic: bool, - // locations of biases - bias_pos: Vec, + // locations of biases + bias_pos: Vec, - // force constants of biases - bias_fc: Vec, + // force constants of biases + bias_fc: Vec, - // bias value cache - bias: Vec, + // bias value cache + bias: Vec, - // histogram weight - pub weights: Vec, + // histogram weight + pub weights: Vec, } impl Dataset { - pub fn new(num_bins: usize, dimens_lengths: Vec, bin_width: Vec, + 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: Vec = vec![0.0; num_bins*num_windows]; - let weights = vec![1.0; num_windows]; - let mut ds = Dataset{ - num_windows, - num_bins, - dimens_lengths, - bin_width, - hist_min, - hist_max, - kT, - histograms, - cyclic, - bias_pos, - bias_fc, - bias, - weights - }; - for window in 0..num_windows { - for bin in 0..num_bins { - let ndx = window * num_bins + bin; - ds.bias[ndx] = ds.calc_bias(bin, window); - } - } - ds + let num_windows = histograms.len(); + let bias: Vec = vec![0.0; num_bins*num_windows]; + let weights = vec![1.0; num_windows]; + let mut ds = Dataset{ + num_windows, + num_bins, + dimens_lengths, + bin_width, + hist_min, + hist_max, + kT, + histograms, + cyclic, + bias_pos, + bias_fc, + bias, + weights + }; + for window in 0..num_windows { + for bin in 0..num_bins { + let ndx = window * num_bins + bin; + ds.bias[ndx] = ds.calc_bias(bin, window); + } + } + ds - } + } - pub fn new_weighted(ds: Dataset, weights: Vec) -> Dataset { - Dataset { - weights, - ..ds - } - } + pub fn new_weighted(ds: Dataset, weights: Vec) -> Dataset { + Dataset { + weights, + ..ds + } + } - pub fn get_weighted_bin_count(&self, bin: usize) -> f64 { - self.histograms.iter().enumerate().map(|(idx,h)| self.weights[idx]*h.bins[bin]).sum() - } + pub fn get_weighted_bin_count(&self, bin: usize) -> f64 { + self.histograms.iter().enumerate().map(|(idx,h)| self.weights[idx]*h.bins[bin]).sum() + } - fn expand_index(&self, bin: usize, lengths: &[usize]) -> Vec { - let mut tmp = bin; - let mut idx = vec![0; lengths.len()]; - for dimen in (1..lengths.len()).rev() { - let denom: usize = lengths.iter().take(dimen).product(); - idx[dimen] = tmp / denom; - tmp %= denom; - } - idx[0] = tmp; - idx - } + fn expand_index(&self, bin: usize, lengths: &[usize]) -> Vec { + let mut tmp = bin; + let mut idx = vec![0; lengths.len()]; + for dimen in (1..lengths.len()).rev() { + let denom: usize = lengths.iter().take(dimen).product(); + idx[dimen] = tmp / denom; + tmp %= denom; + } + idx[0] = tmp; + idx + } - // get center x value for a bin - pub fn get_coords_for_bin(&self, bin: usize) -> Vec { - self.expand_index(bin, &self.dimens_lengths).iter().enumerate().map(|(i, dimen_bin)| { - self.hist_min[i] + self.bin_width[i]*(*dimen_bin as f64 + 0.5) - }).collect() - } + // get center x value for a bin + pub fn get_coords_for_bin(&self, bin: usize) -> Vec { + self.expand_index(bin, &self.dimens_lengths).iter().enumerate().map(|(i, dimen_bin)| { + self.hist_min[i] + self.bin_width[i]*(*dimen_bin as f64 + 0.5) + }).collect() + } - pub fn get_bias(&self, bin: usize, window: usize) -> f64 { - let ndx = window * self.num_bins + bin; - self.bias[ndx] - } + pub fn get_bias(&self, bin: usize, window: usize) -> f64 { + let ndx = window * self.num_bins + bin; + self.bias[ndx] + } - // Harmonic bias calculation: bias = 0.5*k(dx)^2 - // if cyclic is true, lowest and highest bins are assumed to be - // neighbors. This returns exp(U/kT) instead of U for better performance. - fn calc_bias(&self, bin: usize, window: usize) -> f64 { - let dimens = self.dimens_lengths.len(); - // index of the bias value depends on the window und dimension - let bias_ndx: Vec = (0..dimens) - .map(|dimen| { window * dimens + dimen }).collect(); + // Harmonic bias calculation: bias = 0.5*k(dx)^2 + // if cyclic is true, lowest and highest bins are assumed to be + // neighbors. This returns exp(U/kT) instead of U for better performance. + fn calc_bias(&self, bin: usize, window: usize) -> f64 { + let dimens = self.dimens_lengths.len(); + // index of the bias value depends on the window und dimension + let bias_ndx: Vec = (0..dimens) + .map(|dimen| { window * dimens + dimen }).collect(); - // find the N coords, force constants and bias coords - let coord = self.get_coords_for_bin(bin); - let bias_fc: Vec = bias_ndx.iter().map(|ndx| { self.bias_fc[*ndx] }).collect(); - let bias_pos: Vec = bias_ndx.iter().map(|ndx| { self.bias_pos[*ndx] }).collect(); + // find the N coords, force constants and bias coords + let coord = self.get_coords_for_bin(bin); + let bias_fc: Vec = bias_ndx.iter().map(|ndx| { self.bias_fc[*ndx] }).collect(); + let bias_pos: Vec = bias_ndx.iter().map(|ndx| { self.bias_pos[*ndx] }).collect(); - let mut bias_sum = 0.0; - for i in 0..dimens { - let mut dist = (coord[i] - bias_pos[i]).abs(); - if self.cyclic { // periodic conditions - let hist_len = self.hist_max[i] - self.hist_min[i]; - if dist > 0.5 * hist_len { - dist -= hist_len; - } - } - // store exp(U/kT) for better performance - bias_sum += 0.5 * bias_fc[i] * dist * dist - } - (-bias_sum/self.kT).exp() - } + let mut bias_sum = 0.0; + for i in 0..dimens { + let mut dist = (coord[i] - bias_pos[i]).abs(); + if self.cyclic { // periodic conditions + let hist_len = self.hist_max[i] - self.hist_min[i]; + if dist > 0.5 * hist_len { + dist -= hist_len; + } + } + // store exp(U/kT) for better performance + bias_sum += 0.5 * bias_fc[i] * dist * dist + } + (-bias_sum/self.kT).exp() + } } impl fmt::Display for Dataset { - fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { - let mut datapoints: u32 = 0; - for h in &self.histograms { - datapoints += h.num_points; - } - write!(f, "{} windows, {} datapoints", self.num_windows, datapoints) + fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { + let mut datapoints: u32 = 0; + for h in &self.histograms { + datapoints += h.num_points; + } + write!(f, "{} windows, {} datapoints", self.num_windows, datapoints) } } #[cfg(test)] mod tests { - use super::*; - use super::super::k_B; + use super::*; + use super::super::k_B; - macro_rules! assert_delta { + 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 - vec![1.0, 1.0, 3.0, 5.0, 12.0] // bins - ) - } + fn build_hist() -> Histogram { + Histogram::new( + 22, // num_points + vec![1.0, 1.0, 3.0, 5.0, 12.0] // bins + ) + } - fn build_hist_set() -> Dataset { - let h = build_hist(); - Dataset::new( - 5, // num bins - vec![1], - vec![1.0], // bin width - vec![0.0], // hist min - vec![9.0], // hist max - vec![4.5], // x0 - vec![10.0], // fc - 300.0*k_B, // kT - vec![h], // hists - false // cyclic - ) - } + fn build_hist_set() -> Dataset { + let h = build_hist(); + Dataset::new( + 5, // num bins + vec![1], + vec![1.0], // bin width + vec![0.0], // hist min + vec![9.0], // hist max + vec![4.5], // x0 + vec![10.0], // fc + 300.0*k_B, // kT + vec![h], // hists + false // cyclic + ) + } - #[test] - fn calc_bias() { - let ds = build_hist_set(); // k = 10 + #[test] + fn calc_bias() { + let ds = build_hist_set(); // k = 10 - // 3th element -> x=3.5, x0=3.5 - assert_delta!(0.134_722_337_796, ds.calc_bias(3, 0), 0.000_000_01); + // 3th element -> x=3.5, x0=3.5 + assert_delta!(0.134_722_337_796, ds.calc_bias(3, 0), 0.000_000_01); - // 8th element -> x=8.5, x0=3.5 - assert_delta!(1.0, ds.calc_bias(4,0), 0.000_000_01); - - // 1st element -> x=0.5, x0=3.5. non-cyclic! - assert_delta!(0.0, ds.calc_bias(0,0), 0.000_000_1); - } + // 8th element -> x=8.5, x0=3.5 + assert_delta!(1.0, ds.calc_bias(4,0), 0.000_000_01); + + // 1st element -> x=0.5, x0=3.5. non-cyclic! + assert_delta!(0.0, ds.calc_bias(0,0), 0.000_000_1); + } - #[test] - fn calc_biascyclic() { - let mut ds = build_hist_set(); - ds.cyclic = true; + #[test] + fn calc_biascyclic() { + let mut ds = build_hist_set(); + ds.cyclic = true; - // 7th element -> x=3.5, x0=3.5 - assert_delta!(0.134_722_337_796, ds.calc_bias(3, 0), 0.000_000_01); + // 7th element -> x=3.5, x0=3.5 + assert_delta!(0.134_722_337_796, ds.calc_bias(3, 0), 0.000_000_01); - // 8th element -> x=4.5, x0=3.5 - assert_delta!(1.0, ds.calc_bias(4, 0), 0.000_000_01); - + // 8th element -> x=4.5, x0=3.5 + assert_delta!(1.0, ds.calc_bias(4, 0), 0.000_000_01); + - // 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.000_000_000_000_011_776_9, ds.calc_bias(0, 0), 0.000_000_01); + // 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.000_000_000_000_011_776_9, ds.calc_bias(0, 0), 0.000_000_01); - // 2nd element -> x=1.5, x0=3.5 - assert_delta!(0.000_000_01, ds.calc_bias(1, 0), 0.000_000_01); - } + // 2nd element -> x=1.5, x0=3.5 + assert_delta!(0.000_000_01, ds.calc_bias(1, 0), 0.000_000_01); + } - #[test] - fn get_x_for_bin() { - let ds = build_hist_set(); - let expected: Vec = vec![0,1,2,3,4,5,6,7,8].iter() - .map(|x| *x as f64 + 0.5).collect(); + #[test] + fn get_x_for_bin() { + let ds = build_hist_set(); + let expected: Vec = vec![0,1,2,3,4,5,6,7,8].iter() + .map(|x| *x as f64 + 0.5).collect(); expected.iter().enumerate().for_each(|(i, exp)| { assert_approx_eq!(exp, &ds.get_coords_for_bin(i)[0]); }) - } + } - #[test] - fn get_bin_count() { - let ds = Dataset::new( - 5, // num bins - vec![1], - vec![1.0, 1.0], // bin width - vec![0.0, 0.0], // hist min - vec![5.0, 5.0], // hist max - vec![7.5, 7.5], // x0 - vec![10.0, 10.0], // fc - 300.0*k_B, // kT - vec![build_hist(), build_hist()], // hists - false // cyclic - ); - assert_delta!(2.0, ds.get_weighted_bin_count(0), 0.000_000_000_1); - assert_delta!(2.0, ds.get_weighted_bin_count(1), 0.000_000_000_1); - assert_delta!(6.0, ds.get_weighted_bin_count(2), 0.000_000_000_1); - assert_delta!(10.0, ds.get_weighted_bin_count(3), 0.000_000_000_1); - assert_delta!(24.0, ds.get_weighted_bin_count(4), 0.000_000_000_1); - } + #[test] + fn get_bin_count() { + let ds = Dataset::new( + 5, // num bins + vec![1], + vec![1.0, 1.0], // bin width + vec![0.0, 0.0], // hist min + vec![5.0, 5.0], // hist max + vec![7.5, 7.5], // x0 + vec![10.0, 10.0], // fc + 300.0*k_B, // kT + vec![build_hist(), build_hist()], // hists + false // cyclic + ); + assert_delta!(2.0, ds.get_weighted_bin_count(0), 0.000_000_000_1); + assert_delta!(2.0, ds.get_weighted_bin_count(1), 0.000_000_000_1); + assert_delta!(6.0, ds.get_weighted_bin_count(2), 0.000_000_000_1); + assert_delta!(10.0, ds.get_weighted_bin_count(3), 0.000_000_000_1); + assert_delta!(24.0, ds.get_weighted_bin_count(4), 0.000_000_000_1); + } } \ No newline at end of file diff --git a/src/io.rs b/src/io.rs index e46d20e..d356bb5 100644 --- a/src/io.rs +++ b/src/io.rs @@ -30,8 +30,8 @@ pub fn vprintln(s: String, verbose: bool) { // given in the metadata file. This generates at least one Dataset, // or multiple Datasets if convdt is set in the config pub fn read_data(cfg: &Config) -> Result> { - let mut bias_pos: Vec = Vec::new(); - let mut bias_fc: Vec = Vec::new(); + let mut bias_pos: Vec = Vec::new(); + let mut bias_fc: Vec = Vec::new(); let mut timeseries_lengths: Vec = Vec::new(); let mut paths = Vec::new(); @@ -43,7 +43,7 @@ pub fn read_data(cfg: &Config) -> Result> { // start..convdt, start..2*convdt, ... let mut histograms = vec![Vec::new(); dataset_boundaries.len()]; - let kT = cfg.temperature * k_B; + let kT = cfg.temperature * k_B; let bin_width: Vec = (0..cfg.dimens).map(|idx| { (cfg.hist_max[idx] - cfg.hist_min[idx])/(cfg.num_bins[idx] as f64) }).collect(); @@ -55,14 +55,14 @@ pub fn read_data(cfg: &Config) -> Result> { // read each metadata file line and parse it for (line_num,l) in buf.lines().enumerate() { - let line = l.chain_err(|| "Failed to read line")?; + let line = l.chain_err(|| "Failed to read line")?; // skip comments and empty lines if line.starts_with('#') || line.is_empty() { - continue; - } + continue; + } - let split: Vec<&str> = line.split_whitespace().collect(); + let split: Vec<&str> = line.split_whitespace().collect(); if split.len() < 1 + cfg.dimens * 2 { bail!(format!("Wrong number of columns in line {} of metadata file. Empty Line?", line_num+1)); } diff --git a/src/lib.rs b/src/lib.rs index 31b6881..2919676 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -31,18 +31,18 @@ static k_B: f64 = 0.008_314_462_1; // kJ/mol*K // Application config #[derive(Debug)] pub struct Config { - pub metadata_file: String, - pub hist_min: Vec, - pub hist_max: Vec, - pub num_bins: Vec, - pub dimens: usize, - pub verbose: bool, - pub tolerance: f64, - pub max_iterations: usize, - pub temperature: f64, - pub cyclic: bool, - pub output: String, - pub bootstrap: usize, + pub metadata_file: String, + pub hist_min: Vec, + pub hist_max: Vec, + pub num_bins: Vec, + pub dimens: usize, + pub verbose: bool, + pub tolerance: f64, + pub max_iterations: usize, + pub temperature: f64, + pub cyclic: bool, + pub output: String, + pub bootstrap: usize, pub bootstrap_seed: u64, pub start: f64, pub end: f64, @@ -52,7 +52,7 @@ pub struct Config { } impl fmt::Display for Config { - fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { + fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { write!(f, "Metadata={}, hist_min={:?}, hist_max={:?}, bins={:?}, verbose={}, tolerance={}, iterations={}, temperature={}, cyclic={:?}, uncorr={:?}, bootstrap={:?}, seed={:?}, @@ -70,7 +70,7 @@ impl fmt::Display for Config { fn is_converged(old_F: &[f64], new_F: &[f64], tolerance: f64) -> bool { // calculates abs diff between every old and new F and checks if any // is larger than tolerance - !new_F.iter() + !new_F.iter() .zip(old_F.iter()) .map(|x| { (x.0-x.1).abs() }) .any(|diff| { diff > tolerance }) @@ -81,13 +81,13 @@ fn is_converged(old_F: &[f64], new_F: &[f64], tolerance: f64) -> bool { // P(x) = \frac {\sum_{i=1}^N{n_i(x)}} // {\sum_{i=1}^N{ N_i exp(\beta [F_i - U_{bias,i}(x)])}} fn calc_bin_probability(bin: usize, dataset: &Dataset, F: &[f64]) -> f64 { - let mut denom_sum: f64 = 0.0; - let bin_count: f64 = dataset.get_weighted_bin_count(bin); + let mut denom_sum: f64 = 0.0; + let bin_count: f64 = dataset.get_weighted_bin_count(bin); for (window, h) in dataset.histograms.iter().enumerate() { - let bias = dataset.get_bias(bin, window); + let bias = dataset.get_bias(bin, window); denom_sum += (dataset.weights[window] * h.num_points as f64) * bias * F[window]; - } + } bin_count / denom_sum } @@ -109,26 +109,26 @@ fn calc_window_F(window: usize, dataset: &Dataset, P: &[f64]) -> f64 { // offsets F based on previous bias offsets F_prev. This updates the values in // vectors F and P. fn perform_wham_iteration(dataset: &Dataset, F_prev: &[f64], F: &mut Vec, P: &mut Vec) { - // Update P + // Update P // evaluate first WHAM equation for each bin to - // estimate probabilities based on previous offsets (F_prev)) + // estimate probabilities based on previous offsets (F_prev)) (0..dataset.num_bins).into_par_iter() - .map(|bin| { calc_bin_probability(bin, dataset, F_prev) }) - .collect_into_vec(P); + .map(|bin| { calc_bin_probability(bin, dataset, F_prev) }) + .collect_into_vec(P); // Update F - // evaluate second WHAM equation for each window to - // estimate new bias offsets from propabilities - (0..dataset.num_windows).into_par_iter() - .map(|window| {calc_window_F(window, dataset, P)} ) - .collect_into_vec(F); + // evaluate second WHAM equation for each window to + // estimate new bias offsets from propabilities + (0..dataset.num_windows).into_par_iter() + .map(|window| {calc_window_F(window, dataset, P)} ) + .collect_into_vec(F); } // Full WHAM calculation. Calls `perform_wham_iteration` until convergence // criteria are met or max iterations reached. pub fn perform_wham(cfg: &Config, dataset: &Dataset) -> Result<(Vec, Vec, Vec)> { - // allocate required vectors. + // allocate required vectors. // bin probability let mut P: Vec = vec![f64::NAN; dataset.num_bins]; @@ -177,10 +177,10 @@ pub fn perform_wham(cfg: &Config, dataset: &Dataset) } if iteration == cfg.max_iterations { - bail!("WHAM not converged! (max iterations reached)"); + bail!("WHAM not converged! (max iterations reached)"); } - Ok((P, F, F_prev)) + Ok((P, F, F_prev)) } pub fn run(cfg: &Config) -> Result<()>{ @@ -227,17 +227,17 @@ pub fn run(cfg: &Config) -> Result<()>{ // get average difference between two bias offset sets fn diff_avg(F: &[f64], F_prev: &[f64]) -> f64 { - let mut F_sum: f64 = 0.0; - for i in 0..F.len() { - F_sum += (F[i]-F_prev[i]).abs() - } - F_sum / F.len() as f64 + let mut F_sum: f64 = 0.0; + for i in 0..F.len() { + F_sum += (F[i]-F_prev[i]).abs() + } + F_sum / F.len() as f64 } // calculate the normalized free energy from probability values fn calc_free_energy(dataset: &Dataset, P: &[f64]) -> Vec { let mut minimum = f64::MAX; - let mut free_energy: Vec = P.iter() + let mut free_energy: Vec = P.iter() .map(|p| { -dataset.kT * p.ln() }) @@ -257,28 +257,28 @@ fn calc_free_energy(dataset: &Dataset, P: &[f64]) -> Vec { // Print the current WHAM iteration state. Dumps the PMF and associated vectors fn dump_state(dataset: &Dataset, F: &[f64], F_prev: &[f64], P: &[f64], P_std: &[f64], A: &[f64], A_std: &[f64]) { - // TODO fix output of F/F_prev - let out = std::io::stdout(); + // TODO fix output of F/F_prev + let out = std::io::stdout(); let mut lock = out.lock(); - writeln!(lock, "# PMF").unwrap(); - writeln!(lock, "#bin\t\tFree Energy\t\t+/-\t\tP(x)\t\t+/-").unwrap(); - for bin in 0..dataset.num_bins { - writeln!(lock, "{:9.5}\t{:9.5}\t{:9.5}\t{:9.5}\t{:9.5}", + writeln!(lock, "# PMF").unwrap(); + writeln!(lock, "#bin\t\tFree Energy\t\t+/-\t\tP(x)\t\t+/-").unwrap(); + for bin in 0..dataset.num_bins { + writeln!(lock, "{:9.5}\t{:9.5}\t{:9.5}\t{:9.5}\t{:9.5}", bin, A[bin], A_std[bin], P[bin], P_std[bin]).unwrap(); - } - writeln!(lock, "# Bias offsets").unwrap(); - writeln!(lock, "#Window\t\tF\t\tF_prev").unwrap(); - for window in 0..dataset.num_windows { - writeln!(lock, "{}\t{:9.5}\t{:8.8}", + } + writeln!(lock, "# Bias offsets").unwrap(); + writeln!(lock, "#Window\t\tF\t\tF_prev").unwrap(); + for window in 0..dataset.num_windows { + writeln!(lock, "{}\t{:9.5}\t{:8.8}", window, F[window], (F[window]-F_prev[window]).abs()).unwrap(); - } + } } #[cfg(test)] mod tests { - use super::histogram::{Dataset,Histogram}; - use std::f64; + use super::histogram::{Dataset,Histogram}; + use std::f64; use super::k_B; macro_rules! assert_delta { @@ -288,65 +288,65 @@ mod tests { } - fn create_test_dataset() -> Dataset { - let h1 = Histogram::new(10, vec![0.0, 1.0, 1.0, 8.0, 0.0]); - let h2 = Histogram::new(10, vec![0.0, 0.0, 8.0, 1.0, 1.0]); - Dataset::new(5, vec![5], vec![1.0], vec![0.0], vec![4.0], + fn create_test_dataset() -> Dataset { + let h1 = Histogram::new(10, vec![0.0, 1.0, 1.0, 8.0, 0.0]); + let h2 = Histogram::new(10, vec![0.0, 0.0, 8.0, 1.0, 1.0]); + Dataset::new(5, vec![5], vec![1.0], vec![0.0], vec![4.0], vec![1.0, 1.0], vec![10.0, 10.0], 300.0*k_B, vec![h1, h2], false) - } - - #[test] - fn is_converged() { - let new = vec![1.0,1.0]; - let old = vec![0.95, 1.0]; - let tolerance = 0.1; - let converged = super::is_converged(&old, &new, tolerance); - assert!(converged); - - let old = vec![0.8, 1.0]; - let converged = super::is_converged(&old, &new, tolerance); - assert!(!converged); - } + } #[test] - fn calc_bin_probability() { - let dataset = create_test_dataset(); - let F = vec![1.0; dataset.num_bins] ; + fn is_converged() { + let new = vec![1.0,1.0]; + let old = vec![0.95, 1.0]; + let tolerance = 0.1; + let converged = super::is_converged(&old, &new, tolerance); + assert!(converged); + + let old = vec![0.8, 1.0]; + let converged = super::is_converged(&old, &new, tolerance); + assert!(!converged); + } + + #[test] + fn calc_bin_probability() { + let dataset = create_test_dataset(); + let F = vec![1.0; dataset.num_bins] ; let expected = vec!(0.0, 0.082_529_668_703_131_6, 40.923_558_470_974_93, 124_226.700_033_77, 2_308_526_035.528_374_7); expected.iter().enumerate().for_each(|(i, exp)| { - let p = super::calc_bin_probability(i, &dataset, &F); - assert_delta!(exp, p, 0.000_000_1); + let p = super::calc_bin_probability(i, &dataset, &F); + assert_delta!(exp, p, 0.000_000_1); }) } #[test] - fn calc_bias_offset() { - let dataset = create_test_dataset(); - let probability = vec!(0.0, 0.1, 0.2, 0.3, 0.4); + fn calc_bias_offset() { + let dataset = create_test_dataset(); + let probability = vec!(0.0, 0.1, 0.2, 0.3, 0.4); let expected = vec!(15.927_477_169_990_633, 15.927_477_169_990_633); expected.iter().enumerate().for_each(|(i, exp)| { let F = super::calc_window_F(i, &dataset, &probability); assert_delta!(exp, F, 0.000_000_1); }) - } + } - #[test] - fn perform_wham_iteration() { - let dataset = create_test_dataset(); - let prev_F = vec![1.0; dataset.num_windows]; - let mut F = vec![f64::NAN; dataset.num_windows]; - let mut P = vec![f64::NAN; dataset.num_bins]; - super::perform_wham_iteration(&dataset, &prev_F, &mut F, &mut P); + #[test] + fn perform_wham_iteration() { + let dataset = create_test_dataset(); + let prev_F = vec![1.0; dataset.num_windows]; + let mut F = vec![f64::NAN; dataset.num_windows]; + let mut P = vec![f64::NAN; dataset.num_bins]; + super::perform_wham_iteration(&dataset, &prev_F, &mut F, &mut P); let expected_F = vec!(1.0, 1.0); - let expected_P = vec!(0.0, 0.082_529_668_703_131_6, 40.923_558_470_974_93, + let expected_P = vec!(0.0, 0.082_529_668_703_131_6, 40.923_558_470_974_93, 124_226.700_033_77, 2_308_526_035.528_374_7); - for bin in 0..dataset.num_bins { - assert_delta!(expected_P[bin], P[bin], 0.01) - } - for window in 0..dataset.num_windows { - assert_delta!(expected_F[window], F[window], 0.01) - } - - } + for bin in 0..dataset.num_bins { + assert_delta!(expected_P[bin], P[bin], 0.01) + } + for window in 0..dataset.num_windows { + assert_delta!(expected_F[window], F[window], 0.01) + } + + } } \ No newline at end of file diff --git a/src/main.rs b/src/main.rs index 0e81f7a..adacd6d 100644 --- a/src/main.rs +++ b/src/main.rs @@ -13,20 +13,20 @@ use std::process; // Parse command line arguments into a Config struct fn cli() -> Result { - let yaml = load_yaml!("cli.yml"); - let matches = App::from_yaml(yaml).get_matches(); - let metadata_file = matches.value_of("metadata").unwrap().to_string(); - let verbose: bool = matches.is_present("verbose"); - let temperature: f64 = matches.value_of("temperature").unwrap().parse() - .chain_err(|| "Cannot read temperature.")?; - let tolerance: f64 = matches.value_of("tolerance").unwrap_or("0.000001").parse() - .chain_err(|| "Cannot read tolerance.")?; - let max_iterations: usize = matches.value_of("iterations").unwrap_or("100000").parse() - .chain_err(|| "Cannot parse iterations.")?; - let output = matches.value_of("output").unwrap_or("wham.out").to_string(); + let yaml = load_yaml!("cli.yml"); + let matches = App::from_yaml(yaml).get_matches(); + let metadata_file = matches.value_of("metadata").unwrap().to_string(); + let verbose: bool = matches.is_present("verbose"); + let temperature: f64 = matches.value_of("temperature").unwrap().parse() + .chain_err(|| "Cannot read temperature.")?; + let tolerance: f64 = matches.value_of("tolerance").unwrap_or("0.000001").parse() + .chain_err(|| "Cannot read tolerance.")?; + let max_iterations: usize = matches.value_of("iterations").unwrap_or("100000").parse() + .chain_err(|| "Cannot parse iterations.")?; + let output = matches.value_of("output").unwrap_or("wham.out").to_string(); let cyclic: bool = matches.is_present("cyclic"); - let hist_min: Vec = matches.value_of("min_hist").unwrap() + let hist_min: Vec = matches.value_of("min_hist").unwrap() .split(',').map(|x| { if x.to_ascii_lowercase() == "pi" { std::f64::consts::PI @@ -36,7 +36,7 @@ fn cli() -> Result { x.parse().unwrap() } }).collect(); - let hist_max: Vec = matches.value_of("max_hist").unwrap() + let hist_max: Vec = matches.value_of("max_hist").unwrap() .split(',').map(|x| { if x.to_ascii_lowercase() == "pi" { std::f64::consts::PI @@ -46,10 +46,10 @@ fn cli() -> Result { x.parse().unwrap() } }).collect(); - let num_bins: Vec = matches.value_of("bins").unwrap() + let num_bins: Vec = matches.value_of("bins").unwrap() .split(',').map(|x| { x.parse().unwrap() }).collect(); - let bootstrap: usize = matches.value_of("bootstrap").unwrap_or("0").parse() - .chain_err(|| "Cannot parse bootstrap iteration.")?; + let bootstrap: usize = matches.value_of("bootstrap").unwrap_or("0").parse() + .chain_err(|| "Cannot parse bootstrap iteration.")?; let bootstrap_seed: u64 = matches.value_of("bootstrap_seed") .unwrap_or({ let mut rng = rand::thread_rng(); @@ -78,20 +78,20 @@ fn cli() -> Result { let ignore_empty: bool = matches.is_present("ignore_empty"); - Ok(wham::Config{metadata_file, hist_min, hist_max, num_bins, dimens, - verbose, tolerance, max_iterations, temperature, cyclic, output, - bootstrap, bootstrap_seed, start, end, uncorr, convdt, ignore_empty}) + Ok(wham::Config{metadata_file, hist_min, hist_max, num_bins, dimens, + verbose, tolerance, max_iterations, temperature, cyclic, output, + bootstrap, bootstrap_seed, start, end, uncorr, convdt, ignore_empty}) } fn main() { - let cfg = cli().expect("Failed to parse CLI."); - if let Err(error) = wham::run(&cfg) { - eprintln!("Error: {}", error); + let cfg = cli().expect("Failed to parse CLI."); + if let Err(error) = wham::run(&cfg) { + eprintln!("Error: {}", error); - for e in error.iter().skip(1) { - eprintln!("Reason: {}", e) - } - process::exit(1); - } + for e in error.iter().skip(1) { + eprintln!("Reason: {}", e) + } + process::exit(1); + } }