diff --git a/src/histogram.rs b/src/histogram.rs index cd4f02f..20c7c13 100644 --- a/src/histogram.rs +++ b/src/histogram.rs @@ -4,33 +4,16 @@ use std::cell::RefCell; // One histogram #[derive(Debug)] pub struct Histogram { - // offset of this histogram bins from the global histogram - first: usize, - - // offset of the last element of the histogram. TODO required? - last: usize, - // total number of data points stored in the histogram pub num_points: u32, // histogram bins - bins: Vec + pub bins: Vec } impl Histogram { - pub fn new(first: usize, last: usize, num_points: u32, bins: Vec) -> Histogram { - assert_eq!(last-first+1, bins.len(), "histogram length does not match first/last."); - Histogram {first, last, num_points, bins} - } - - // Returns the value of a bin if the bin is present in this - // histogram - pub fn get_bin_count(&self, bin: usize) -> Option { - if bin < self.first || bin > self.last { - None - } else { - Some(self.bins[bin-self.first]) - } + pub fn new(num_points: u32, bins: Vec) -> Histogram { + Histogram {num_points, bins} } } @@ -40,17 +23,20 @@ pub struct Dataset { // number of histogram windows (number of simulations) pub num_windows: usize, - // number of global histogram bins + // total number of bins pub num_bins: usize, - // min value of the histogram - hist_min: f64, + // number of bins in each dimension + pub dimens_lengths: Vec, - // max value of the histogram - hist_max: f64, + // min values of the histogram in each dimension + hist_min: Vec, - // width of a bin in unit of x - bin_width: f64, + // max values of the histogram in each dimension + hist_max: Vec, + + // width of a bin in unit of its dimension + bin_width: Vec, // value of kT pub kT: f64, @@ -72,16 +58,17 @@ pub struct Dataset { } impl Dataset { - - pub fn new(num_bins: usize, bin_width: f64, hist_min: f64, hist_max: f64, bias_pos: Vec, bias_fc: Vec, kT: f64, histograms: Vec, cyclic: bool) -> 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]); Dataset{ num_windows, num_bins, + dimens_lengths, bin_width, hist_min, - hist_max, + hist_max, kT, histograms, cyclic, @@ -91,36 +78,60 @@ impl Dataset { } } - + fn expand_index(&self, bin: usize, lengths: &Vec) -> Vec { + let mut tmp = bin; + let mut idx = vec![0; lengths.len()]; + for dimen in (1..lengths.len()).rev() { + let denom = lengths.iter().take(dimen).fold(1, |s,&x| s*x); + idx[dimen] = tmp / denom; + tmp = 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() + } + // Harmonic bias calculation: bias = 0.5*k(dx)^2 // if cyclic is true, lowest and highest bins are assumed to be // neighbors pub fn calc_bias(&self, bin: usize, window: usize) -> f64 { - let ndx = bin + (self.num_bins*window); + let ndx = window * self.num_bins + bin; let mut cache = self.bias.borrow_mut(); match cache[ndx] { Some(val) => val, None => { - let x = self.get_x_for_bin(bin); - let mut dx = (x-self.bias_pos[window]).abs(); - if self.cyclic { - let hist_len = self.hist_max-self.hist_min; - if dx > 0.5*hist_len { - dx -= hist_len; + // TODO optimize this part! + let dimens = self.hist_min.len(); + + let bias_ndx: Vec = (0..dimens) + .map(|dimen| { window * dimens + dimen }).collect(); + 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 { + let hist_len = self.hist_max[i] - self.hist_min[i]; + if dist > 0.5 * hist_len { + dist -= hist_len; + } } + bias_sum += 0.5 * bias_fc[i] * dist * dist } - let bias = 0.5*self.bias_fc[window]*dx*dx; - cache[ndx] = Some(bias); - bias + cache[ndx] = Some(bias_sum); + bias_sum } } } - // get center x value for a bin - pub fn get_x_for_bin(&self, bin: usize) -> f64 { - self.hist_min + self.bin_width * ((bin as f64) + 0.5) - } - } impl fmt::Display for Dataset { @@ -140,8 +151,6 @@ mod tests { fn build_hist() -> Histogram { Histogram::new( - 5, // first - 9, // last 22, // num_points vec![1.0, 1.0, 3.0, 5.0, 12.0] // bins ) diff --git a/src/io.rs b/src/io.rs index bf01a57..a6b89f0 100644 --- a/src/io.rs +++ b/src/io.rs @@ -32,15 +32,24 @@ pub fn read_data(cfg: &Config) -> Option { let mut bias_pos: Vec = Vec::new(); let mut bias_fc: Vec = Vec::new(); let mut histograms: Vec = Vec::new(); + 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(); + let num_bins = cfg.num_bins.iter().fold(1, |state, &bins| state*bins); + let dimens_length = cfg.num_bins.clone(); + let f = File::open(&cfg.metadata_file).unwrap_or_else(|x| { eprintln!("Failed to read metadata from {}. {}", &cfg.metadata_file, x); process::exit(1) - }); + }); let buf = BufReader::new(&f); + // read each metadata file line and parse it for l in buf.lines() { let line = l.unwrap(); + // skip comments and empty lines if line.starts_with("#") || line.len() == 0 { continue; @@ -84,15 +93,36 @@ pub fn read_data(cfg: &Config) -> Option { } if histograms.len() > 0 { - let bin_width: Vec = (0..cfg.dimens).map(|idx| { - (cfg.hist_max[idx] - cfg.hist_min[idx])/(cfg.num_bins[idx] as f64) - }).collect(); - Some(Dataset::new(cfg.num_bins[0], bin_width[0], cfg.hist_min[0], cfg.hist_max[0], bias_pos, bias_fc, kT, histograms, cfg.cyclic)) + Some(Dataset::new(num_bins, dimens_length, bin_width, cfg.hist_min.clone(), cfg.hist_max.clone(), bias_pos, bias_fc, kT, histograms, cfg.cyclic)) } else { None } } +// transforms a multidimensional index into a one dimensional index +// indeces: multidimensional indeces +// lengths: length of the matrix in each dimension +// returns an index if the matrix is flattened to a one dimensional vector +// example for 3 dimensions N,M,O: idx = i_O + l_O*l_M*i_M + l_O*l_M*l_N*i_N +fn flat_index(indeces: &Vec, lengths: &Vec) -> usize { + let mut idx = 0; + for i in 0..indeces.len() { + idx += indeces[i]*lengths[0..i].iter() + .fold(1, |state, &l| { state * l }); + } + idx +} + +// returns true if the values are inside the histogram boundaries defined by cfg +fn is_in_hist_boundaries(values: &Vec, cfg: &Config) -> bool { + for dimen in 0..cfg.dimens { + if values[dimen] < cfg.hist_min[dimen] || values[dimen] > cfg.hist_max[dimen] { + return false + } + } + true +} + // parse a timeseries file into a histogram fn read_window_file(window_file: &str, cfg: &Config) -> Option { let f = File::open(window_file).unwrap_or_else(|x| { @@ -100,63 +130,61 @@ fn read_window_file(window_file: &str, cfg: &Config) -> Option { process::exit(1) }); let buf = BufReader::new(&f); - - let mut global_hist = vec![0.0; cfg.num_bins[0]]; - let bin_width = (cfg.hist_max[0] - cfg.hist_min[0])/(cfg.num_bins[0] as f64); + + // total number of bins is the product of all dimensions length + let total_bins = cfg.num_bins.iter().fold(1, |s, &x| { s*x }); + let mut hist = vec![0.0; total_bins]; + + // bin width for each dimension: (max-min)/bins + let bin_width: Vec = (0..cfg.dimens).map(|idx| { + (cfg.hist_max[idx] - cfg.hist_min[idx])/(cfg.num_bins[idx] as f64) + }).collect(); + + // read and parse each timeseries line for l in buf.lines() { let line = l.unwrap(); + // skip comments and empty lines if line.starts_with("#") || line.starts_with("@") || line.len() == 0 { continue; } - let (_, x) = scan_fmt!(&line, "{} {}", f64, f64); + let mut split = line.split_whitespace(); + split.next(); // skip time/step column - match x { - Some(x) => { - if x > cfg.hist_min[0] && x < cfg.hist_max[0] { - let bin_ndx = ((x-cfg.hist_min[0]) / bin_width) as usize; - global_hist[bin_ndx] += 1.0; - } - } - None => { - eprintln!("{}, Failed to read datapoint from line: {}", &window_file, &line); - process::exit(1); - } - } - } - let mut max_bin: usize = 0; - let mut min_bin: usize =(cfg.num_bins[0]-1) as usize; - for bin in 0..global_hist.len() { - if global_hist[bin] != 0.0 && bin > max_bin { - max_bin = bin; - } - if global_hist[bin] != 0.0 && bin < min_bin { - min_bin = bin; + let values: Vec = (0..cfg.dimens).collect::>().iter().map(|_| { + split.next().unwrap().parse::().unwrap() + }).collect(); + + if is_in_hist_boundaries(&values, cfg) { + let bin_indeces = (0..cfg.dimens).map(|dimen: usize| { + let val = values[dimen]; + ((val-cfg.hist_min[dimen]) / bin_width[dimen]) as usize + }).collect(); + let index = flat_index(&bin_indeces, &cfg.num_bins); + hist[index] += 1.0; } } - if (max_bin == min_bin && global_hist[max_bin] == 0.0) || max_bin < min_bin { - None // zero length histogram - } else { - // trim global hist to save memory - global_hist.truncate(max_bin+1); - global_hist.drain(..min_bin); - - let num_points: f64 = global_hist.iter().sum(); - Some(Histogram::new(min_bin, max_bin, num_points as u32, global_hist)) + let num_points: f64 = hist.iter().sum(); + if num_points == 0.0 { + return None } + Some(Histogram::new(num_points as u32, hist)) } +// TODO multidimensional output pub fn write_results(out_file: &str, ds: &Dataset, free: &Vec, prob: &Vec) -> Result<(), Box> { - let mut output = File::create(out_file)?; - writeln!(output, "#{:8}\t{:8}\t{:8}", "x", "Free Energy", "Probability"); - for bin in 0..free.len() { - let x = ds.get_x_for_bin(bin); - writeln!(output, "{:8.6}\t{:8.6}\t{:8.6}", x, free[bin], prob[bin])?; - } - Ok(()) + let mut output = File::create(out_file)?; + writeln!(output, "#{}\t{}\t{}", "x", "Free Energy", "Probability"); // TODO better format (coord1, coord2..) + for bin in 0..free.len() { + let coords = ds.get_coords_for_bin(bin); + let coords_str: String = coords.iter().map(|c| {format!("{:8.6}", c)}) + .collect::>().join("\t"); + writeln!(output, "{}\t{:8.6}\t{:8.6}", coords_str, free[bin], prob[bin])?; + } + Ok(()) } #[cfg(test)] @@ -188,9 +216,9 @@ mod tests { // assert_eq!(1, h.first); // assert_eq!(6, h.last); assert_eq!(11, h.num_points); - assert_eq!(2.0, h.get_bin_count(1).unwrap()); - assert_eq!(2.0, h.get_bin_count(2).unwrap()); - assert_eq!(1.0, h.get_bin_count(6).unwrap()); + assert_eq!(2.0, h.bins[1]); + assert_eq!(2.0, h.bins[2]); + assert_eq!(1.0, h.bins[6]); } diff --git a/src/lib.rs b/src/lib.rs index d52d792..ec418e7 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,8 +1,5 @@ #![allow(non_snake_case)] -#[macro_use] -extern crate scan_fmt; - pub mod io; pub mod histogram; @@ -55,9 +52,7 @@ fn calc_bin_probability(bin: usize, ds: &Dataset, F: &Vec) -> f64 { let mut bin_count: f64 = 0.0; for window in 0..ds.num_windows { let h: &Histogram = &ds.histograms[window]; - if let Some(count) = h.get_bin_count(bin) { - bin_count += count; - } + bin_count += h.bins[bin]; let bias = ds.calc_bias(bin, window); let bias_offset = ((F[window] - bias) / ds.kT).exp(); denom_sum += (h.num_points as f64) * bias_offset; @@ -245,7 +240,7 @@ fn dump_state(ds: &Dataset, F: &Vec, F_prev: &Vec, P: &Vec, A: &V println!("# PMF"); println!("#x\t\tFree Energy\t\tP(x)"); for bin in 0..ds.num_bins { - let x = ds.get_x_for_bin(bin); + let x = ds.get_coords_for_bin(bin)[0]; // TODO println!("{:9.5}\t{:9.5}\t{:9.5}", x, A[bin], P[bin]); } println!("# Bias offsets"); @@ -274,8 +269,8 @@ mod tests { } fn create_test_ds() -> Dataset { - let h1 = Histogram::new(0, 2, 10, vec![3.0, 4.0, 3.0]); - let h2 = Histogram::new(0, 3, 20, vec![3.0, 2.0, 5.0, 10.0]); + let h1 = Histogram::new(10, vec![0.0, 0.0, 3.0, 4.0, 3.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]); + let h2 = Histogram::new(20, vec![0.0, 0.0, 0.0, 3.0, 2.0, 5.0, 10.0, 0.0, 0.0, 0.0, 0.0]); Dataset::new(4, 1.0, 0.0, 4.0, vec![1.0, 2.0], vec![10.0, 10.0], 2.479, vec![h1, h2], false) }