Skip to content

Commit

Permalink
explict single threaded coverage
Browse files Browse the repository at this point in the history
  • Loading branch information
Rob Patro committed Nov 18, 2024
1 parent 68f4df8 commit 9d2a9da
Show file tree
Hide file tree
Showing 3 changed files with 35 additions and 23 deletions.
2 changes: 1 addition & 1 deletion Cargo.toml
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
[package]
name = "oarfish"
version = "0.6.3"
version = "0.6.4"
edition = "2021"
authors = [
"Zahra Zare Jousheghani <zzare@umd.edu>",
Expand Down
41 changes: 21 additions & 20 deletions src/util/binomial_probability.rs
Original file line number Diff line number Diff line change
Expand Up @@ -183,26 +183,9 @@ pub fn binomial_continuous_prob(txps: &mut [TranscriptInfo], bins: &u32, threads

let _log_span = info_span!("binomial_continuous_prob").entered();
info!("computing coverage probabilities");
rayon::ThreadPoolBuilder::new()
.num_threads(threads)
.build()
.unwrap();

txps.par_iter_mut().enumerate().for_each(|(_i, t)| {
let compute_txp_coverage_probs = |_i: usize, t: &mut TranscriptInfo| {
let temp_prob: Vec<f64> = if *bins != 0 {
/*
let bin_counts: Vec<f32>;
let bin_lengths: Vec<f32>;
let _num_discarded_read_temp: usize;
let _bin_coverage: Vec<f64>;
(
bin_counts,
bin_lengths,
_num_discarded_read_temp,
_bin_coverage,
) = bin_transcript_normalize_counts(t, bins); //binning the transcript length and obtain the counts and length vectors
//==============================================================================================
*/
let min_cov = t.total_weight / 100.;
t.coverage_bins.iter_mut().for_each(|elem| *elem += min_cov);
let (bin_counts, bin_lengths) = t.get_normalized_counts_and_lengths();
Expand All @@ -216,8 +199,26 @@ pub fn binomial_continuous_prob(txps: &mut [TranscriptInfo], bins: &u32, threads
} else {
std::unimplemented!("coverage model with 0 bins is not currently implemented");
};

t.coverage_prob = temp_prob;
});
};

// if we are requesting only a single thread, then don't bother with
// the overhead of e.g. creating a thread pool and doing parallel
// iteration, etc.
if threads > 1 {
let pool = rayon::ThreadPoolBuilder::new()
.num_threads(threads)
.build()
.unwrap();
pool.install(|| {
txps.par_iter_mut()
.enumerate()
.for_each(|(_i, t)| compute_txp_coverage_probs(_i, t));
});
} else {
txps.iter_mut()
.enumerate()
.for_each(|(_i, t)| compute_txp_coverage_probs(_i, t));
}
info!("done");
}
15 changes: 13 additions & 2 deletions src/util/write_function.rs
Original file line number Diff line number Diff line change
Expand Up @@ -276,12 +276,23 @@ pub fn write_out_prob(
txps.clear();
txp_probs.clear();

const DISPLAY_THRESH: f64 = 0.001;
let mut denom2 = 0.0_f64;

for (a, p, cp) in izip!(alns, probs, coverage_probs) {
let target_id = a.ref_id as usize;
let prob = *p as f64;
let cov_prob = if model_coverage { *cp } else { 1.0 };
txps.push(target_id);
txp_probs.push(((prob * cov_prob) / denom).clamp(0.0, 1.0));
let nprob = ((prob * cov_prob) / denom).clamp(0.0, 1.0);
if nprob >= DISPLAY_THRESH {
txps.push(target_id);
txp_probs.push(nprob);
denom2 += nprob;
}
}

for p in txp_probs.iter_mut() {
*p /= denom2;
}

let txp_ids = txps
Expand Down

0 comments on commit 9d2a9da

Please sign in to comment.