+ changed interface for decoded plif features -> directly access double array
[qpalma.git] / QPalmaDP / qpalma_dp.cpp
index b9cd6da..c2af278 100644 (file)
@@ -6,7 +6,7 @@ using namespace std;
   myalign([nr_paths], dna, est, {h}, matchmatrix, donor, acceptor, remove_duplicate_scores, ...
   print_matrix) */
 
-Alignment::Alignment(int numq) {
+Alignment::Alignment(int numQPlifs, int numq) {
       len = 0;
       limits = 0;
       penalties = 0;
@@ -23,15 +23,9 @@ Alignment::Alignment(int numq) {
       mmatrix_param = 0;
       alignmentscores = 0;
       qualityScoresAllPaths = 0;
-      numQualSuppPoints = numq;
-}
 
-void Alignment::setQualityMatrix(double* qMat, int length){
-   if(qualityMatrix == 0)
-      qualityMatrix = new double[length];
-
-   for(int i=0; i<length; i++)
-      qualityMatrix[i] = qMat[i];
+      numQualSuppPoints = numq;
+      totalNumPlifs = numQPlifs;
 }
 
 void Alignment::getDNAEST(){}
@@ -50,6 +44,10 @@ void Alignment::myalign(int nr_paths_p, char* dna, int dna_len_p, char* est,
   /***************************************************************************/ 
   // initialize alignment matrices and call fill_matrix()  
   /***************************************************************************/
+
+    // dnaest
+   double *DNA_ARRAY = 0;
+   double *EST_ARRAY = 0;
   
   //possible donor positions
   int nr_donor_sites = 0 ;
@@ -94,7 +92,7 @@ void Alignment::myalign(int nr_paths_p, char* dna, int dna_len_p, char* est,
    est_align_size = (est_len-1)*nr_paths;
    mmatrix_param_size = (mlen*mlen)*nr_paths;
    alignmentscores_size = nr_paths; //alignment score for each path/matrix
-   qScores_size = 25*nr_paths; //alignment score for each path/matrix
+   qScores_size = totalNumPlifs*nr_paths; //alignment score for each path/matrix
 
    //printf("dna_len is %d est_len is %d mmatrix_len is %d",splice_align_size, est_align_size, mmatrix_param_size);
 
@@ -104,13 +102,14 @@ void Alignment::myalign(int nr_paths_p, char* dna, int dna_len_p, char* est,
   alignmentscores = new double[nr_paths]; //alignment score for each path/matrix
   qualityScoresAllPaths = new penalty_struct[qScores_size];
   int qidx, pidx;
-  for(qidx=0;qidx<25*nr_paths;qidx++) {
+
+  for(qidx=0;qidx<totalNumPlifs*nr_paths;qidx++) {
    penalty_struct p;
    init_penalty_struct(p);
    p.len = numQualSuppPoints;
    p.limits = (REAL*) calloc(p.len,sizeof(REAL));
    for(pidx=0;pidx<p.len;pidx++)
-      p.limits[pidx] = qualityScores[qidx%25].limits[pidx];
+      p.limits[pidx] = qualityScores[qidx%totalNumPlifs].limits[pidx];
 
    p.penalties = (REAL*) calloc(p.len,sizeof(REAL));
 
@@ -132,9 +131,6 @@ void Alignment::myalign(int nr_paths_p, char* dna, int dna_len_p, char* est,
   memset(alignmentscores, -1, nr_paths*sizeof(double)); //fills alignmentscores with zeros
   // printf("after memset...\n");
   
-    // dnaest
-   double *DNA_ARRAY = 0;
-   double *EST_ARRAY = 0;
 
   for (int z=0; z<nr_paths; z++) {
     result_length = 0 ;
@@ -142,90 +138,81 @@ void Alignment::myalign(int nr_paths_p, char* dna, int dna_len_p, char* est,
     int* s_align = splice_align + (dna_len-1)*z;  //pointer
     int* e_align = est_align + (est_len-1)*z ;    //pointer
     int* mparam  = mmatrix_param + (mlen*mlen)*z; //pointer
-    penalty_struct* qparam = qualityScoresAllPaths + (25*z);
+    penalty_struct* qparam = qualityScoresAllPaths + (totalNumPlifs*z);
 
     printf("before call to result_align...\n");
-    bool no_more_path = result_align(matrices, z, est_len, dna_len, &result_length, est, dna, prb, chastity, s_align, e_align, mparam, alignmentscores, max_score_positions, qparam);
+    bool no_more_path = result_align(matrices, z, est_len, dna_len, &result_length, est, dna, prb, chastity, s_align, e_align, mparam, alignmentscores, max_score_positions, qparam, currentMode);
     printf("after call to result_align...\n");
 
-  int qidx;
-  for(qidx=0;qidx<25*nr_paths;qidx++) {
-      penalty_struct p;
-      p = qualityScoresAllPaths[qidx];
-      int pidx;
-     for(pidx=0;pidx<p.len;pidx++)
-         printf("%f ",p.penalties[pidx]);
-      printf("\n");
+  for(qidx=0;qidx<totalNumPlifs*nr_paths;qidx++) {
+   penalty_struct p;
+   p = qualityScoresAllPaths[qidx];
+   for(pidx=0;pidx<p.len;pidx++)
+      printf("%f ",p.penalties[pidx]);
+   printf("\n");
   }
 
+ printf("after printf...\n");
 
-   if(DNA_ARRAY != 0) {
-       delete[] DNA_ARRAY;
-       delete[] EST_ARRAY;
-    }
-
-   DNA_ARRAY = new double[result_length];
-   EST_ARRAY = new double[result_length];
+   //if(DNA_ARRAY != 0) {
+   //    delete[] DNA_ARRAY;
+   //    delete[] EST_ARRAY;
+   // }
 
-    //backtracking
-    int i = max_score_positions[2*z] ; //i (est)
-    int j = max_score_positions[2*z +1] ; //j (dna)
-    int prev_i = ((Pre_score*)matrices[z] + i*dna_len +j)->prev_i; 
-    int prev_j = ((Pre_score*)matrices[z] + i*dna_len +j)->prev_j;
-    int prev_m_no = ((Pre_score*)matrices[z] + i*dna_len +j)->prev_matrix_no;
-    
-    for (int ii=result_length; ii>0; ii--) {
-      if ((prev_i == (i-1)) && (prev_j == (j-1))) { // match or mismatch
-       DNA_ARRAY[ii-1] = check_char(dna[j-1]) ;
-       EST_ARRAY[ii-1] = check_char(est[i-1]) ;
-      }
-      else if ((prev_i == (i)) && (prev_j == (j-1))) {// gap on EST_ARRAY
-       DNA_ARRAY[ii-1] = check_char(dna[j-1]) ;
-       EST_ARRAY[ii-1] = 0 ; //gap
-      }
-      else if ((prev_i == (i-1)) && (prev_j == (j))) { // gap on DNA_ARRAY
-       DNA_ARRAY[ii-1] = 0 ; //gap
-       EST_ARRAY[ii-1] = check_char(est[i-1]) ;
-      }
-      else {// splice site
-       for (j; j > prev_j; j--) {
-         DNA_ARRAY[ii-1] = check_char(dna[j-1]) ;
-         EST_ARRAY[ii-1] = 6 ; //intron
-         ii-- ;
-       }
-       ii++ ; // last ii-- too much (because done in for-loop) 
-      }
-      
-      i = prev_i;
-      j = prev_j;    
-      prev_i = ((Pre_score*)matrices[prev_m_no] + i*dna_len + j)->prev_i; 
-      prev_j = ((Pre_score*)matrices[prev_m_no] + i*dna_len + j)->prev_j; 
-      prev_m_no = ((Pre_score*)matrices[prev_m_no] + i*dna_len + j)->prev_matrix_no;
-    }
-    
-  } //end of z
-
-   // printf("Content of qualityScoresAllPaths\n");
-   // for(int iii=0;iii<qScores_size;iii++) {
-   //    for(int jjj=0;jjj<qualityScoresAllPaths[iii].len;jjj++) {
-   //       printf("%f ",qualityScoresAllPaths[iii].penalties[jjj]);
-   //    }
-   //    printf("\n");
+   //DNA_ARRAY = new double[result_length];
+   //EST_ARRAY = new double[result_length];
+
+   // //backtracking
+   // int i = max_score_positions[2*z] ; //i (est)
+   // int j = max_score_positions[2*z +1] ; //j (dna)
+   // int prev_i = ((Pre_score*)matrices[z] + i*dna_len +j)->prev_i; 
+   // int prev_j = ((Pre_score*)matrices[z] + i*dna_len +j)->prev_j;
+   // int prev_m_no = ((Pre_score*)matrices[z] + i*dna_len +j)->prev_matrix_no;
+   // 
+   // for (int ii=result_length; ii>0; ii--) {
+   //   if ((prev_i == (i-1)) && (prev_j == (j-1))) { // match or mismatch
+       //DNA_ARRAY[ii-1] = check_char(dna[j-1]) ;
+       //EST_ARRAY[ii-1] = check_char(est[i-1]) ;
+   //   }
+   //   else if ((prev_i == (i)) && (prev_j == (j-1))) {// gap on EST_ARRAY
+       //DNA_ARRAY[ii-1] = check_char(dna[j-1]) ;
+       //EST_ARRAY[ii-1] = 0 ; //gap
+   //   }
+   //   else if ((prev_i == (i-1)) && (prev_j == (j))) { // gap on DNA_ARRAY
+       //DNA_ARRAY[ii-1] = 0 ; //gap
+       //EST_ARRAY[ii-1] = check_char(est[i-1]) ;
+   //   }
+   //   else {// splice site
+       //for (j; j > prev_j; j--) {
+       //  DNA_ARRAY[ii-1] = check_char(dna[j-1]) ;
+       //  EST_ARRAY[ii-1] = 6 ; //intron
+       //  ii-- ;
+       //}
+       //ii++ ; // last ii-- too much (because done in for-loop)       
+   //   }
+   //   
+   //   i = prev_i;
+   //   j = prev_j;    
+   //   prev_i = ((Pre_score*)matrices[prev_m_no] + i*dna_len + j)->prev_i; 
+   //   prev_j = ((Pre_score*)matrices[prev_m_no] + i*dna_len + j)->prev_j; 
+   //   prev_m_no = ((Pre_score*)matrices[prev_m_no] + i*dna_len + j)->prev_matrix_no;
    // }
+   // 
+   } //end of z
 
-   if(DNA_ARRAY != 0) {
-    delete[] DNA_ARRAY;
-    delete[] EST_ARRAY;
-    }
+   //if(DNA_ARRAY != 0) {
+   // delete[] DNA_ARRAY;
+   // delete[] EST_ARRAY;
+   // }
 
    for (int i=nr_paths-1;i>=0; i--)
       delete[] matrices[i];
 
-   //printf("Leaving myalign...\n");
+   printf("Leaving myalign...\n");
 }
 
 void Alignment::getAlignmentResults(int* s_align, int* e_align,
-      int* mmatrix_p, double* alignscores, penalty_struct* qScores) {
+      int* mmatrix_p, double* alignscores, double* qScores) {
 
    //printf("Entering getAlignmentResults...\n");
 
@@ -242,10 +229,17 @@ void Alignment::getAlignmentResults(int* s_align, int* e_align,
    for(idx=0; idx<alignmentscores_size; idx++)
       alignscores[idx] = alignmentscores[idx];
    
+   penalty_struct currentPlif;
+   int pidx;
+   int ctr=0;
    for(idx=0; idx<qScores_size; idx++) {
-      init_penalty_struct(qScores[idx]);
-      qScores[idx] = qualityScoresAllPaths[idx];
-      copy_penalty_struct(&qualityScoresAllPaths[idx],&qScores[idx]);
+      currentPlif = qualityScoresAllPaths[idx];
+      for(pidx=0; pidx<currentPlif.len; pidx++) {
+         qScores[ctr] = currentPlif.penalties[pidx];
+         ctr++;
+      }
    }
+
+   //printf("ctr is %d\n",ctr);
    //printf("Leaving getAlignmentResults...\n");
 }