fixed bug in SNV trans mat, added test/data

.gitignore

@@ -1,9 +1,10 @@
 *.[oa]
 build/*
 *~
-test/
 release/
 release*
 *.bk
 bk/
 recompile*
+.DS_Store
+test/results/*

src/clone-predict.cpp

@@ -65,7 +65,7 @@ void Clone::set_TransMat_cna( gsl_matrix * Trans, int chr){
 
 void Clone::set_TransMat_snv(){
   if (TransMat_snv==NULL){
-    TransMat_snv = new gsl_matrix * [cnaEmit->nSamples];
+    TransMat_snv = new gsl_matrix * [snvEmit->nSamples];
     for (int s=0; s<snvEmit->nSamples;s++) TransMat_snv[s] = NULL;
   }
   for (int s=0; s<snvEmit->nSamples;s++){

src/clone-prior.cpp

@@ -50,7 +50,7 @@ void Clone::initialize_snv_prior_param(){// SNV prior, conditional on max-tcn
   if (initial_snv_prior_param != NULL) gsl_matrix_free(initial_snv_prior_param);
   initial_snv_prior_param = gsl_matrix_calloc( maxtcn+1, maxtcn+1);
   gsl_matrix_set( initial_snv_prior_param, 0, 0, snv_fpr);
-  double p = snv_pen_high;// initial penalty for higher genotypes
+  double p = snv_pen_high;// penalty for higher genotypes
   for (int cn=1; cn <= maxtcn; cn++){
     if ( all_maxtcn.count(cn) == 0 ) continue;
     gsl_vector_view subrow = gsl_matrix_subrow( initial_snv_prior_param, cn, 0, cn+1);
@@ -68,9 +68,10 @@ void Clone::set_snv_prior( gsl_matrix * snv_prior_param){
   double fpr = gsl_matrix_get( snv_prior_param, 0, 0);
   snv_prior.clear();
   std::map<int, vector<int> >::iterator it;
+  gsl_vector * mem = gsl_vector_calloc(nLevels);    
   for (it=maxtcn_per_clone.begin(); it != maxtcn_per_clone.end(); it++){
     int chr = it->first;
-    snv_prior[chr] = gsl_vector_calloc(nLevels);
+    mem->data[0] = 0;
     for (int i=1; i<nLevels; i++){
       double p=1.0;
       for (int j=0; j<nClones; j++){
@@ -83,20 +84,24 @@ void Clone::set_snv_prior( gsl_matrix * snv_prior_param){
 	  break;
 	}
       }
-      gsl_vector_set( snv_prior[chr], i, p);
+      //gsl_vector_set( snv_prior[chr], i, p);
+      mem->data[i] = p;
     }
     //normalize and logify...
-    double norm = gsl_blas_dasum(snv_prior[chr]);
+    double norm = gsl_blas_dasum(mem);
     if (norm<=0) abort();
-    gsl_vector_scale( snv_prior[chr], (1.0-fpr) / norm);
-    gsl_vector_set( snv_prior[chr], 0, fpr);
+    gsl_vector_scale( mem, (1.0-fpr) / norm);
+    gsl_vector_set( mem, 0, fpr);
     if (snvEmit->log_space){
       for (int l=0; l<nLevels; l++){
-	double val = gsl_vector_get( snv_prior[chr], l);
-	gsl_vector_set( snv_prior[chr], l, val>0 ? log(val) : logzero);
+	double val = gsl_vector_get( mem, l);
+	gsl_vector_set( mem, l, val>0 ? log(val) : logzero);
       }
     }
+    snv_prior[chr] = gsl_vector_calloc(nLevels);
+    gsl_vector_memcpy(snv_prior[chr],mem);
   }
+  gsl_vector_free(mem);
 }
 
 //CNA + BAF (+SNV) mode...

src/cloneHD-inference.cpp

@@ -100,7 +100,7 @@ int infer_clones( gsl_matrix * Clones, gsl_vector * Mass, Clone * myClone, cmdl_
     best_priors[n] = NULL;
   }
   //exit(0);
-  double bic=0, llh=0, max_llh=0,cna_llh=0,baf_llh=0,snv_llh=0, max_bic=0; 
+  double bic=0, llh=0, max_llh=0, cna_llh=0, baf_llh=0, snv_llh=0, max_bic=0; 
   int steps=0,btrial=0,bestn=0;
   int nT = myClone->nTimes;
   // *** NO CLONE SCENARIO (n==0) ***

src/emission.cpp

@@ -908,7 +908,7 @@ void Emission::coarse_grain_jumps( int sample, double plow, int range){
     while(it != remaining.end()){
       it++;
       p = it->second;
-      if( p>10.0*pmin || it->first != lidx+1 || p<10.0*gmin) break;
+      if ( p>10.0*pmin || it->first != lidx+1 || p<10.0*gmin) break;
       if (it->first < idx+range) ps *= 1.0 - p;
       pmin = min(pmin,p);
       lidx = it->first;