Skip to content

Commit 78b6d1c

Browse files
pmd: many places the vector indicating array_status had been passed around despite only being used in a few of those places. These were deleted
1 parent 7ef3bd3 commit 78b6d1c

File tree

3 files changed

+134
-170
lines changed

3 files changed

+134
-170
lines changed

src/analysis/pmd.cpp

Lines changed: 6 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -533,24 +533,19 @@ get_domain_scores(const vector<bool> &classes,
533533

534534

535535
static void
536-
build_domains(const bool VERBOSE,
537-
const vector<SimpleGenomicRegion> &bins,
538-
const vector<double> &post_scores,
536+
build_domains(const vector<SimpleGenomicRegion> &bins,
539537
const vector<size_t> &reset_points,
540538
const vector<bool> &classes,
541-
const vector<size_t> &dists_btwn_bins,
542539
vector<GenomicRegion> &domains) {
543540

544541
size_t n_bins = 0, reset_idx = 1, prev_end = 0;
545542
bool in_domain = false;
546-
double score = 0;
547543
for (size_t i = 0; i < classes.size(); ++i) {
548544
if (reset_points[reset_idx] == i) {
549545
if (in_domain) {
550546
domains.back().set_end(prev_end);
551547
domains.back().set_score(n_bins);
552548
n_bins = 0;
553-
score = 0;
554549
in_domain = false;
555550
}
556551
++reset_idx;
@@ -561,13 +556,11 @@ build_domains(const bool VERBOSE,
561556
in_domain = true;
562557
}
563558
++n_bins;
564-
score += post_scores[i];
565559
}
566560
else if (in_domain) {
567561
domains.back().set_end(prev_end);
568562
domains.back().set_score(n_bins);
569563
n_bins = 0;
570-
score = 0;
571564
in_domain = false;
572565
}
573566
prev_end = bins[i].get_end();
@@ -1029,10 +1022,9 @@ get_min_reads_for_confidence(const double conf_level) {
10291022
// ADS: this function will return num_lim<size_t>::max() if the
10301023
// fraction of "good" bins is zero for all attempted bin sizes.
10311024
static size_t
1032-
binsize_selection(const bool &VERBOSE, const size_t resolution,
1033-
const size_t min_bin_sz, const size_t max_bin_sz,
1034-
const double conf_level, const double min_frac_passed,
1035-
const string &cpgs_file) {
1025+
binsize_selection(const size_t resolution, const size_t min_bin_sz,
1026+
const size_t max_bin_sz, const double conf_level,
1027+
const double min_frac_passed, const string &cpgs_file) {
10361028

10371029
const size_t min_cov_to_pass = get_min_reads_for_confidence(conf_level);
10381030

@@ -1239,8 +1231,7 @@ main_pmd(int argc, const char **argv) {
12391231
for (size_t i = 0; i < n_replicates && !insufficient_data; ++i) {
12401232
const bool arrayData = check_if_array_data(cpgs_file[i]);
12411233
if (!arrayData) {
1242-
bin_size = binsize_selection(VERBOSE, resolution,
1243-
min_bin_size, max_bin_size,
1234+
bin_size = binsize_selection(resolution, min_bin_size, max_bin_size,
12441235
confidence_interval, prop_accept,
12451236
cpgs_file[i]);
12461237
if (bin_size == num_lim<size_t>::max())
@@ -1414,8 +1405,7 @@ main_pmd(int argc, const char **argv) {
14141405
}
14151406

14161407
vector<GenomicRegion> domains;
1417-
build_domains(VERBOSE, bins[0], scores, reset_points, classes,
1418-
dists_btwn_bins, domains);
1408+
build_domains(bins[0], reset_points, classes, domains);
14191409

14201410
size_t good_pmd_count = 0;
14211411
vector<GenomicRegion> good_domains;

src/common/TwoStateHMM_PMD.cpp

Lines changed: 70 additions & 85 deletions
Original file line numberDiff line numberDiff line change
@@ -18,8 +18,6 @@
1818

1919
#include "TwoStateHMM_PMD.hpp"
2020

21-
// #pragma omp <rest of pragma>
22-
2321
using std::vector;
2422
using std::pair;
2523
using std::setw;
@@ -115,7 +113,7 @@ TwoStateHMM::TransitionPosteriors_rep(
115113
start_trans[0], start_trans[1],
116114
trans[0][0], trans[0][1], end_trans[0],
117115
trans[1][0], trans[1][1], end_trans[1],
118-
fg_distro, bg_distro, transition, array_status, llr_scores);
116+
fg_distro, bg_distro, transition, llr_scores);
119117
}
120118

121119
////////////////////////////////////////////////////////////////////////////////
@@ -133,8 +131,7 @@ TwoStateHMM::forward_algorithm_rep(
133131
const double lp_bb, const double lp_bt,
134132
const vector<unique_ptr<EmissionDistribution> > &fg_distro,
135133
const vector<unique_ptr<EmissionDistribution> > &bg_distro,
136-
vector<pair<double, double> > &f,
137-
const vector<bool> &array_status) const {
134+
vector<pair<double, double> > &f) const {
138135
size_t NREP = vals.size();
139136
f[start].first = lp_sf;
140137
f[start].second = lp_sb;
@@ -163,17 +160,15 @@ TwoStateHMM::forward_algorithm_rep(
163160

164161

165162
double
166-
TwoStateHMM::backward_algorithm_rep(
167-
const vector<vector<pair<double, double> > > &vals,
168-
const size_t start, const size_t end,
169-
const double lp_sf, const double lp_sb,
170-
const double lp_ff, const double lp_fb,
171-
const double lp_ft, const double lp_bf,
172-
const double lp_bb, const double lp_bt,
173-
const vector<unique_ptr<EmissionDistribution> > &fg_distro,
174-
const vector<unique_ptr<EmissionDistribution> > &bg_distro,
175-
vector<pair<double, double> > &b,
176-
const vector<bool> &array_status) const {
163+
TwoStateHMM::backward_algorithm_rep(const vector<vector<pair<double, double> > > &vals,
164+
const size_t start, const size_t end,
165+
const double lp_sf, const double lp_sb,
166+
const double lp_ff, const double lp_fb,
167+
const double lp_ft, const double lp_bf,
168+
const double lp_bb, const double lp_bt,
169+
const vector<unique_ptr<EmissionDistribution> > &fg_distro,
170+
const vector<unique_ptr<EmissionDistribution> > &bg_distro,
171+
vector<pair<double, double> > &b) const {
177172

178173
size_t NREP = vals.size();
179174
b[end - 1].first = lp_ft;
@@ -220,12 +215,10 @@ TwoStateHMM::estimate_transitions_rep(
220215
const vector<unique_ptr<EmissionDistribution> > &bg_distro,
221216
const double lp_ff, const double lp_fb,
222217
const double lp_bf, const double lp_bb,
223-
const double lp_ft, const double lp_bt,
224218
vector<double> &ff_vals,
225219
vector<double> &fb_vals,
226220
vector<double> &bf_vals,
227-
vector<double> &bb_vals,
228-
const vector<bool> &array_status) const {
221+
vector<double> &bb_vals) const {
229222
size_t NREP = vals.size();
230223
for (size_t i = start + 1; i < end; ++i) {
231224
const size_t k = i - 1;
@@ -261,8 +254,7 @@ TwoStateHMM::single_iteration_rep(const vector<vector<pair<double, double> > > &
261254
double &p_ff, double &p_fb, double &p_ft,
262255
double &p_bf, double &p_bb, double &p_bt,
263256
vector<unique_ptr<EmissionDistribution> > &fg_distro,
264-
vector<unique_ptr<EmissionDistribution> > &bg_distro,
265-
const vector<bool> &array_status) const {
257+
vector<unique_ptr<EmissionDistribution> > &bg_distro) const {
266258

267259
size_t NREP= values.size();
268260
vector<double> log_fg_expected;
@@ -299,7 +291,7 @@ TwoStateHMM::single_iteration_rep(const vector<vector<pair<double, double> > > &
299291
lp_ff, lp_fb, lp_ft,
300292
lp_bf, lp_bb, lp_bt,
301293
fg_distro, bg_distro,
302-
forward, array_status);
294+
forward);
303295

304296
const double backward_score =
305297
backward_algorithm_rep(values,
@@ -309,7 +301,7 @@ TwoStateHMM::single_iteration_rep(const vector<vector<pair<double, double> > > &
309301
lp_ff, lp_fb, lp_ft,
310302
lp_bf, lp_bb, lp_bt,
311303
fg_distro, bg_distro,
312-
backward, array_status);
304+
backward);
313305

314306
if (DEBUG && (fabs(score - backward_score)/
315307
max(score, backward_score)) > 1e-10) {
@@ -319,8 +311,8 @@ TwoStateHMM::single_iteration_rep(const vector<vector<pair<double, double> > > &
319311

320312
estimate_transitions_rep(values, reset_points[i], reset_points[i + 1],
321313
forward, backward, score, fg_distro, bg_distro,
322-
lp_ff, lp_fb, lp_bf, lp_bb, lp_ft, lp_bt,
323-
ff_vals, fb_vals, bf_vals, bb_vals, array_status);
314+
lp_ff, lp_fb, lp_bf, lp_bb,
315+
ff_vals, fb_vals, bf_vals, bb_vals);
324316

325317
total_score += score;
326318
}
@@ -530,7 +522,7 @@ TwoStateHMM::BaumWelchTraining_rep(
530522
p_sf_est, p_sb_est,
531523
p_ff_est, p_fb_est, p_ft_est,
532524
p_bf_est, p_bb_est, p_bt_est,
533-
fg_distro, bg_distro,array_status);
525+
fg_distro, bg_distro);
534526

535527
if (DEBUG) {
536528
cerr << "S_F\t" << p_sf_est << endl
@@ -582,17 +574,18 @@ TwoStateHMM::BaumWelchTraining_rep(
582574

583575

584576
void
585-
TwoStateHMM::PosteriorScores_rep(
586-
const vector<vector<pair<double, double> > > &values,
587-
const vector<size_t> &reset_points,
588-
const vector<double> &start_trans,
589-
const vector<vector<double> > &trans,
590-
const vector<double> &end_trans,
591-
const vector<double> fg_alpha, const vector<double> fg_beta,
592-
const vector<double> bg_alpha, const vector<double> bg_beta,
593-
const vector<bool> &classes,
594-
vector<double> &llr_scores,
595-
const vector<bool> &array_status) const {
577+
TwoStateHMM::PosteriorScores_rep(const vector<vector<pair<double, double> > > &values,
578+
const vector<size_t> &reset_points,
579+
const vector<double> &start_trans,
580+
const vector<vector<double> > &trans,
581+
const vector<double> &end_trans,
582+
const vector<double> fg_alpha,
583+
const vector<double> fg_beta,
584+
const vector<double> bg_alpha,
585+
const vector<double> bg_beta,
586+
const vector<bool> &classes,
587+
vector<double> &llr_scores,
588+
const vector<bool> &array_status) const {
596589

597590
vector<unique_ptr<EmissionDistribution> > fg_distro;
598591
vector<unique_ptr<EmissionDistribution> > bg_distro;
@@ -619,23 +612,19 @@ TwoStateHMM::PosteriorScores_rep(
619612
start_trans[0], start_trans[1],
620613
trans[0][0], trans[0][1], end_trans[0],
621614
trans[1][0], trans[1][1], end_trans[1],
622-
fg_distro, bg_distro, classes, llr_scores, array_status);
615+
fg_distro, bg_distro, classes, llr_scores);
623616
}
624617

625-
626618
void
627619
TwoStateHMM::PosteriorScores_rep(
628-
const vector<vector<pair<double, double> > > &values,
629-
const vector<size_t> &reset_points,
630-
double p_sf, double p_sb,
631-
double p_ff, double p_fb, double p_ft,
632-
double p_bf, double p_bb, double p_bt,
633-
const vector<unique_ptr<EmissionDistribution> > &fg_distro,
634-
const vector<unique_ptr<EmissionDistribution> > &bg_distro,
635-
const vector<bool> &classes,
636-
vector<double> &llr_scores,
637-
const vector<bool> &array_status) const {
638-
620+
const vector<vector<pair<double, double>>> &values,
621+
const vector<size_t> &reset_points,
622+
double p_sf, double p_sb,
623+
double p_ff, double p_fb,
624+
double p_ft, double p_bf, double p_bb, double p_bt,
625+
const vector<unique_ptr<EmissionDistribution>> &fg_distro,
626+
const vector<unique_ptr<EmissionDistribution>> &bg_distro,
627+
const vector<bool> &classes, vector<double> &llr_scores) const {
639628
double total_score = 0;
640629

641630
const double lp_sf = log(p_sf);
@@ -664,7 +653,7 @@ TwoStateHMM::PosteriorScores_rep(
664653
lp_ff, lp_fb, lp_ft,
665654
lp_bf, lp_bb, lp_bt,
666655
fg_distro, bg_distro,
667-
forward, array_status);
656+
forward);
668657

669658
const double backward_score =
670659
backward_algorithm_rep(values,
@@ -674,7 +663,7 @@ TwoStateHMM::PosteriorScores_rep(
674663
lp_ff, lp_fb, lp_ft,
675664
lp_bf, lp_bb, lp_bt,
676665
fg_distro, bg_distro,
677-
backward,array_status);
666+
backward);
678667

679668
if (DEBUG && (fabs(score - backward_score)/
680669
max(score, backward_score)) > 1e-10)
@@ -695,26 +684,24 @@ TwoStateHMM::PosteriorScores_rep(
695684
}
696685
}
697686

698-
699687
double
700-
TwoStateHMM::PosteriorDecoding_rep(
701-
const vector<vector<pair<double, double> > > &values,
702-
const vector<size_t> &reset_points,
703-
const vector<double> &start_trans,
704-
const vector<vector<double> > &trans,
705-
const vector<double> &end_trans,
706-
const vector<double> fg_alpha, const vector<double> fg_beta,
707-
const vector<double> bg_alpha, const vector<double> bg_beta,
708-
vector<bool> &classes,
709-
vector<double> &llr_scores,
710-
const vector<bool> &array_status) const {
688+
TwoStateHMM::PosteriorDecoding_rep(const vector<vector<pair<double, double> > > &values,
689+
const vector<size_t> &reset_points,
690+
const vector<double> &start_trans,
691+
const vector<vector<double> > &trans,
692+
const vector<double> &end_trans,
693+
const vector<double> fg_alpha, const vector<double> fg_beta,
694+
const vector<double> bg_alpha, const vector<double> bg_beta,
695+
vector<bool> &classes,
696+
vector<double> &llr_scores,
697+
const vector<bool> &array_status) const {
711698

712699
vector<unique_ptr<EmissionDistribution> > fg_distro;
713700
vector<unique_ptr<EmissionDistribution> > bg_distro;
714701

715702
size_t NREP = values.size();
716703
for (size_t i = 0; i < NREP; ++i) {
717-
if(array_status[i]) {
704+
if (array_status[i]) {
718705
fg_distro.emplace_back(new Beta(fg_alpha[i], fg_beta[i]));
719706
bg_distro.emplace_back(new Beta(bg_alpha[i], bg_beta[i]));
720707
}
@@ -734,20 +721,20 @@ TwoStateHMM::PosteriorDecoding_rep(
734721
start_trans[0], start_trans[1],
735722
trans[0][0], trans[0][1], end_trans[0],
736723
trans[1][0], trans[1][1], end_trans[1],
737-
fg_distro, bg_distro, classes, llr_scores,array_status);
724+
fg_distro, bg_distro, classes, llr_scores);
738725
}
739726

740727

741728
void
742729
TwoStateHMM::TransitionPosteriors_rep(
743-
const vector<vector<pair<double, double> > > &values,
730+
const vector<vector<pair<double, double>>> &values,
744731
const vector<size_t> &reset_points,
745732
double p_sf, double p_sb,
746733
double p_ff, double p_fb, double p_ft,
747734
double p_bf, double p_bb, double p_bt,
748-
const vector<std::unique_ptr<EmissionDistribution> > &fg_distro,
749-
const vector<std::unique_ptr<EmissionDistribution> > &bg_distro,
750-
const size_t transition, const vector<bool> &array_status,
735+
const vector<unique_ptr<EmissionDistribution>> &fg_distro,
736+
const vector<unique_ptr<EmissionDistribution>> &bg_distro,
737+
const size_t transition,
751738
vector<double> &scores) const {
752739

753740
double total_score = 0;
@@ -776,12 +763,12 @@ TwoStateHMM::TransitionPosteriors_rep(
776763
lp_sf, lp_sb,
777764
lp_ff, lp_fb, lp_ft,
778765
lp_bf, lp_bb, lp_bt,
779-
fg_distro, bg_distro, forward, array_status);
766+
fg_distro, bg_distro, forward);
780767
const double backward_score =
781768
backward_algorithm_rep(values, reset_points[i], reset_points[i + 1],
782769
lp_sf, lp_sb, lp_ff, lp_fb, lp_ft,
783770
lp_bf, lp_bb, lp_bt,
784-
fg_distro, bg_distro, backward, array_status);
771+
fg_distro, bg_distro, backward);
785772

786773
if (DEBUG && (fabs(score - backward_score)/
787774
max(score, backward_score)) > 1e-10)
@@ -832,17 +819,15 @@ TwoStateHMM::TransitionPosteriors_rep(
832819

833820

834821
double
835-
TwoStateHMM::PosteriorDecoding_rep(
836-
const vector<vector<pair<double, double> > > &values,
837-
const vector<size_t> &reset_points,
838-
double p_sf, double p_sb,
839-
double p_ff, double p_fb, double p_ft,
840-
double p_bf, double p_bb, double p_bt,
841-
const vector<unique_ptr<EmissionDistribution> > &fg_distro,
842-
const vector<unique_ptr<EmissionDistribution> > &bg_distro,
843-
vector<bool> &classes,
844-
vector<double> &llr_scores,
845-
const vector<bool> &array_status) const {
822+
TwoStateHMM::PosteriorDecoding_rep(const vector<vector<pair<double, double> > > &values,
823+
const vector<size_t> &reset_points,
824+
double p_sf, double p_sb,
825+
double p_ff, double p_fb, double p_ft,
826+
double p_bf, double p_bb, double p_bt,
827+
const vector<unique_ptr<EmissionDistribution> > &fg_distro,
828+
const vector<unique_ptr<EmissionDistribution> > &bg_distro,
829+
vector<bool> &classes,
830+
vector<double> &llr_scores) const {
846831
double total_score = 0;
847832

848833
const double lp_sf = log(p_sf);
@@ -870,7 +855,7 @@ TwoStateHMM::PosteriorDecoding_rep(
870855
lp_ff, lp_fb, lp_ft,
871856
lp_bf, lp_bb, lp_bt,
872857
fg_distro, bg_distro,
873-
forward,array_status);
858+
forward);
874859

875860
const double backward_score =
876861
backward_algorithm_rep(values,
@@ -880,7 +865,7 @@ TwoStateHMM::PosteriorDecoding_rep(
880865
lp_ff, lp_fb, lp_ft,
881866
lp_bf, lp_bb, lp_bt,
882867
fg_distro, bg_distro,
883-
backward,array_status);
868+
backward);
884869

885870
if (DEBUG && (fabs(score - backward_score)/
886871
max(score, backward_score)) > 1e-10)

0 commit comments

Comments
 (0)