diff --git a/changelog.md b/changelog.md index 6746ac8..5a0dba1 100644 --- a/changelog.md +++ b/changelog.md @@ -1,10 +1,12 @@ # changelog for cloneHD/filterHD -## v1.17.1 / to come +## v1.17.1 / 01.03.2014 * BAF posterior symmetrized for output * CNA transition matrix penalizes clones with zero copies of a segment - +* fixed bug in SNV prior computation +* added pre-processor directives for conditional openMP compilation + ## v1.17.0 / 25.02.2014 major release ### changed the way SNV priors are computed: diff --git a/src/clone-llh.cpp b/src/clone-llh.cpp index 5fb68fc..c1dd1a9 100644 --- a/src/clone-llh.cpp +++ b/src/clone-llh.cpp @@ -40,12 +40,16 @@ double Clone::get_all_total_llh(){ baf_total_llh = 0.0; snv_total_llh = 0.0; int sample; +#ifdef _OPENMP int nt = min( cnaEmit->nSamples, omp_get_max_threads()); #pragma omp parallel for schedule( dynamic, 1) default(shared) num_threads(nt) +#endif for ( sample=0; sample < cnaEmit->nSamples; sample++){ double llh,ent; Clone::do_cna_Fwd( sample, llh); +#ifdef _OPENMP #pragma omp critical +#endif { cna_total_llh += llh; } @@ -70,12 +74,16 @@ double Clone::get_all_total_llh(){ // // BAF if ( bafEmit->is_set ){ +#ifdef _OPENMP int nt = min( bafEmit->nSamples, omp_get_max_threads()); #pragma omp parallel for schedule( dynamic, 1) default(shared) num_threads(nt) +#endif for ( sample=0; sample < bafEmit->nSamples; sample++){//START PARALLEL FOR double llh=0,ent=0; Clone::do_baf_Fwd( sample, llh); +#ifdef _OPENMP #pragma omp critical +#endif { baf_total_llh += llh; } @@ -89,12 +97,16 @@ double Clone::get_all_total_llh(){ // // SNV if( snvEmit->is_set ){ +#ifdef _OPENMP int nt = min( snvEmit->nSamples, omp_get_max_threads()); #pragma omp parallel for schedule( dynamic, 1) default(shared) num_threads(nt) +#endif for ( sample=0; sample < snvEmit->nSamples; sample++){//START PARALLEL FOR double llh = 0.0; Clone::do_snv_Fwd(sample, llh); +#ifdef _OPENMP #pragma omp critical +#endif { snv_total_llh += llh; } @@ -127,12 +139,16 @@ double Clone::get_cna_total_llh(){ int sample; save_cna_alpha = 0; cna_total_llh = 0.0; +#ifdef _OPENMP int nt = min( cnaEmit->nSamples, omp_get_max_threads()); #pragma omp parallel for schedule( dynamic, 1) default(shared) num_threads(nt) +#endif for ( sample=0; sample < cnaEmit->nSamples; sample++){ double llh; Clone::do_cna_Fwd( sample, llh); +#ifdef _OPENMP #pragma omp critical +#endif { cna_total_llh += llh; } @@ -146,12 +162,16 @@ double Clone::get_baf_total_llh(){ save_baf_alpha = 0; baf_total_llh = 0.0; int sample; +#ifdef _OPENMP int nt = min( bafEmit->nSamples, omp_get_max_threads()); #pragma omp parallel for schedule( dynamic, 1) default(shared) num_threads(nt) +#endif for ( sample=0; sample< bafEmit->nSamples; sample++){ double llh; Clone::do_baf_Fwd( sample, llh); +#ifdef _OPENMP #pragma omp critical +#endif { baf_total_llh += llh; } @@ -165,12 +185,16 @@ double Clone::get_snv_total_llh(){ int sample; save_snv_alpha = 0; snv_total_llh = 0.0; +#ifdef _OPENMP int nt = min( snvEmit->nSamples, omp_get_max_threads()); #pragma omp parallel for schedule( dynamic, 1) default(shared) num_threads(nt) +#endif for ( sample=0; sample< snvEmit->nSamples; sample++){ double llh; Clone::do_snv_Fwd( sample, llh); +#ifdef _OPENMP #pragma omp critical +#endif { snv_total_llh += llh; } diff --git a/src/clone-prior.cpp b/src/clone-prior.cpp index ba739c5..808c6d3 100644 --- a/src/clone-prior.cpp +++ b/src/clone-prior.cpp @@ -139,25 +139,21 @@ void Clone::set_snv_prior_map(){//either via BAF or else via CNA } } } - else if (cnaEmit->is_set){//via CNA posterior only... - if ( snv_prior_from_cna_map == NULL){//allocate - snv_prior_from_cna_map = gsl_matrix_alloc( maxtcn+1, maxtcn+1); - } - gsl_matrix_set_zero( snv_prior_from_cna_map); - double p = (snvEmit->connect) ? 1.0 : snv_pen;// penalty for high genotypes - for (int cn=0; cn <= maxtcn; cn++){ - for (int i=0; i <= cn; i++){ - gsl_matrix_set( snv_prior_from_cna_map, i, cn, pow(p,i)); - } - //normalize... - gsl_vector_view col = gsl_matrix_column( snv_prior_from_cna_map, cn); - double norm = gsl_blas_dasum(&col.vector); - if (norm <=0.0) abort(); - gsl_vector_scale( &col.vector, 1.0 / norm); - } + //via CNA posterior only... + if ( snv_prior_from_cna_map == NULL){//allocate + snv_prior_from_cna_map = gsl_matrix_alloc( maxtcn+1, maxtcn+1); } - else{ - abort(); + gsl_matrix_set_zero( snv_prior_from_cna_map); + double p = (snvEmit->connect) ? 1.0 : snv_pen;// penalty for high genotypes + for (int cn=0; cn <= maxtcn; cn++){ + for (int i=0; i <= cn; i++){ + gsl_matrix_set( snv_prior_from_cna_map, i, cn, pow(p,i)); + } + //normalize... + gsl_vector_view col = gsl_matrix_column( snv_prior_from_cna_map, cn); + double norm = gsl_blas_dasum(&col.vector); + if (norm <=0.0) abort(); + gsl_vector_scale( &col.vector, 1.0 / norm); } } diff --git a/src/clone-update.cpp b/src/clone-update.cpp index f9dbc0c..0f9472f 100644 --- a/src/clone-update.cpp +++ b/src/clone-update.cpp @@ -601,7 +601,9 @@ void Clone::get_cnaEmitLog(){//only ever used for BAF data for (int s=0; snSamples; s++){ cnaEmitLog[t][s] = gsl_matrix_calloc( cnaEmit->nEvents[s], cnaEmit->gridSize+1); int evt; -#pragma omp parallel for schedule( dynamic, 1) default(shared) +#ifdef _OPENMP +#pragma omp parallel for schedule( dynamic, 1) default(shared) +#endif for ( evt=0; evtnEvents[s]; evt++){ gsl_vector * mem = gsl_vector_calloc(cnaEmit->gridSize+1); gsl_vector * rnd = gsl_vector_calloc(cnaEmit->gridSize+1); @@ -648,7 +650,9 @@ void Clone::get_bafEmitLog(){//only ever used for BAF data for (int s=0; snSamples; s++){ bafEmitLog[t][s] = gsl_matrix_calloc( bafEmit->nEvents[s], bafEmit->gridSize+1); int evt; +#ifdef _OPENMP #pragma omp parallel for schedule( dynamic, 1) default(shared) +#endif for (evt=0; evtnEvents[s]; evt++){ unsigned int n,N; gsl_vector * mem = gsl_vector_calloc(bafEmit->gridSize+1); @@ -703,7 +707,9 @@ void Clone::get_snvEmitLog(){//only ever used for SNV data printf("\r%2i/%2i %2i/%2i...",t+1,nTimes,s+1,snvEmit->nSamples); cout<nEvents[s]; evt++){ gsl_matrix * mem = gsl_matrix_calloc(s1,s2); gsl_matrix * rnd_vec = gsl_matrix_calloc(s1,s2); diff --git a/src/clone.h b/src/clone.h index fe4f02e..f1b984e 100755 --- a/src/clone.h +++ b/src/clone.h @@ -12,7 +12,11 @@ #include #include #include + +#ifdef _OPENMP #include +#endif + //#include // GSL headers... diff --git a/src/cloneHD-inference.cpp b/src/cloneHD-inference.cpp index e097fab..06ca262 100644 --- a/src/cloneHD-inference.cpp +++ b/src/cloneHD-inference.cpp @@ -458,7 +458,9 @@ double get_clones_cna( gsl_matrix *& clones, myClone->gamma_cna = new gsl_matrix * [cnaEmit->nSamples]; int s; snv_llh=0; +#ifdef _OPENMP #pragma omp parallel for schedule( dynamic, 1) default(shared) +#endif for ( s=0; snSamples; s++){ int cna_sample = cnaEmit->idx_of[snvEmit->chr[s]]; myClone->alpha_cna[cna_sample]=NULL; @@ -468,7 +470,9 @@ double get_clones_cna( gsl_matrix *& clones, myClone->map_mean_tcn( cnaEmit, cna_sample, snvEmit); double l; myClone->do_snv_Fwd( s, l); - #pragma omp critical +#ifdef _OPENMP +#pragma omp critical +#endif { snv_llh += l; } @@ -534,7 +538,9 @@ double get_clones_cna( gsl_matrix *& clones, int s; snv_llh = 0.0; myClone->save_snv_alpha = 0; +#ifdef _OPENMP #pragma omp parallel for schedule( dynamic, 1) default(shared) +#endif for ( s=0; snSamples; s++){ int cna_sample = cnaEmit->idx_of[snvEmit->chr[s]]; myClone->alpha_cna[cna_sample]=NULL; @@ -544,7 +550,9 @@ double get_clones_cna( gsl_matrix *& clones, myClone->map_mean_tcn( cnaEmit, cna_sample, snvEmit); double l=0; myClone->do_snv_Fwd(s,l); +#ifdef _OPENMP #pragma omp critical +#endif { snv_llh += l; } @@ -1190,7 +1198,9 @@ void snv_bulk_update(Clone * myClone){ } int cna_sample=0; int s; +#ifdef _OPENMP #pragma omp parallel for schedule( dynamic, 1) default(shared) +#endif for ( s=0; snSamples; s++){ int snvChr = snvEmit->chr[s]; if (cnaEmit->is_set && cnaEmit->chrs.count(snvChr) == 1){ diff --git a/src/jump-diffusion.cpp b/src/jump-diffusion.cpp index 6dac905..e3b9640 100755 --- a/src/jump-diffusion.cpp +++ b/src/jump-diffusion.cpp @@ -73,10 +73,15 @@ double JumpDiffusion::get_total_llh(){ int sample; total_llh = 0.0; // SAMPLES: -#pragma omp parallel for schedule( dynamic, 1) default(shared) +#ifdef _OPENMP + int nt = min( nSamples, omp_get_max_threads()); +#pragma omp parallel for schedule( dynamic, 1) default(shared) num_threads(nt) +#endif for (sample=0; sample #include #include +#ifdef _OPENMP +#include +#endif // GSL headers... #include "gsl/gsl_vector.h"