fixed bug in snv prior map; added openMP conditional compilation

ivazquez · Mar 1, 2014 · 8544816 · 8544816
1 parent a98996a
commit 8544816
Show file tree

Hide file tree

Showing 8 changed files with 75 additions and 23 deletions.
diff --git 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
@@ -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
@@ -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
@@ -601,7 +601,9 @@ void Clone::get_cnaEmitLog(){//only ever used for BAF data
     for (int s=0; s<cnaEmit->nSamples; 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; evt<cnaEmit->nEvents[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; s<bafEmit->nSamples; 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; evt<bafEmit->nEvents[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<<flush;
       int evt;
+#ifdef _OPENMP
 #pragma omp parallel for schedule( dynamic, 1) default(shared)
+#endif
       for (evt=0; evt<snvEmit->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
@@ -12,7 +12,11 @@
 #include <set>
 #include <vector>
 #include <algorithm>
+
+#ifdef _OPENMP
 #include <omp.h>
+#endif
+
 //#include <unordered_map>
 
 // GSL headers...

diff --git 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; s<snvEmit->nSamples; 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; s<snvEmit->nSamples; 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; s<snvEmit->nSamples; 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
@@ -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<nSamples; sample++){
     double llh = JumpDiffusion::do_Fwd(sample);
+#ifdef _OPENMP
 #pragma omp critical
+#endif
     {
       total_llh += llh;
     }
@@ -187,7 +192,9 @@ int JumpDiffusion::set_DiffProp( gsl_matrix * propagator, double sd){
 void JumpDiffusion::set_pstay(){
   int s;
   // SAMPLES:
+#ifdef _OPENMP
 #pragma omp parallel for schedule( dynamic, 1) default(shared)
+#endif
   for (s=0; s<nSamples; s++){
     pstay[s][0] = 1.0;
     for (int l=1; l< nSites[s]; l++){

diff --git a/src/jump-diffusion.h b/src/jump-diffusion.h
@@ -10,6 +10,9 @@
 #include <string>
 #include <map>
 #include <vector>
+#ifdef _OPENMP
+#include <omp.h>
+#endif
 
 // GSL headers...
 #include "gsl/gsl_vector.h"