+ fixded index bug in feature count
[qpalma.git] / QPalmaDP / result_align.cpp
index 3b95634..85778f8 100644 (file)
@@ -5,59 +5,64 @@ using namespace std;
 
 void increaseFeatureCount(penalty_struct* qparam, int dnanum, int estnum, double estprb) {
 
-   if(dnanum == 0) {
-      //printf("Zero entry\n");
-      return;
-   }
-
-   penalty_struct currentStruct = qparam[(estnum-1)*5+(dnanum-1)];
-
-   //   printf("before\n");
-   //   int p_idx;
-   //   for(p_idx=0;p_idx<10;p_idx++) {
-   //      printf("%f ",currentStruct.limits[p_idx]);
-   //   }
-   //   printf("\n");
-   //
-   //   for(p_idx=0;p_idx<10;p_idx++) {
-   //      printf("%f ",currentStruct.penalties[p_idx]);
-   //   }
-   //   printf("\n");
+   //printf("Current index %d dna/est: %d %d\n",(estnum-1)*6+(dnanum),dnanum,estnum);
+   int currentPos = (estnum-1)*6+(dnanum);
+   penalty_struct currentStruct = qparam[currentPos];
+
+      //printf("before\n");
+      //int p_idx;
+      //for(p_idx=0;p_idx<currentStruct.len;p_idx++) {
+      //   printf("%f ",currentStruct.limits[p_idx]);
+      //}
+      //printf("\n");
+   
+      //for(p_idx=0;p_idx<currentStruct.len;p_idx++) {
+      //   printf("%f ",currentStruct.penalties[p_idx]);
+      //}
+      //printf("\n");
 
    double value = estprb;
    int Lower = 0;
    int idx;
+
    for (idx=0;idx<currentStruct.len;idx++) {
       if (currentStruct.limits[idx] <= value)
          Lower++;
    }
 
-   if (Lower == 0)
+   if (Lower == 0) {
          currentStruct.penalties[0] += 1;
+         return;
+   }
 
-   if (Lower == currentStruct.len)
+   if (Lower == (currentStruct.len-1)) {
          currentStruct.penalties[currentStruct.len-1] += 1;
+         return;
+   }
 
    Lower -= 1;
    int Upper = Lower+1; // x-werte bleiben fest
-   currentStruct.penalties[Upper] += 1; 
-   currentStruct.penalties[Lower] += 1; 
+
+   double weightup  = 1.0*(value - currentStruct.limits[Lower]) / (currentStruct.limits[Upper] - currentStruct.limits[Lower]);
+   double weightlow = 1.0*(currentStruct.limits[Upper] - value) / (currentStruct.limits[Upper] - currentStruct.limits[Lower]);
+   currentStruct.penalties[Upper] += weightup;
+   currentStruct.penalties[Lower] += weightlow; 
 
    //printf("estprb/Lower/Upper %f %d %d\n",estprb,Lower,Upper);
 
-   qparam[(estnum-1)*5+(dnanum-1)] = currentStruct;
+   qparam[currentPos] = currentStruct;
 
    //printf("after\n");
-   //for(p_idx=0;p_idx<10;p_idx++) {
+   //for(p_idx=0;p_idx<currentStruct.len;p_idx++) {
    //   printf("%f ",currentStruct.limits[p_idx]);
    //}
    //printf("\n");
-   //for(p_idx=0;p_idx<10;p_idx++)
+   //for(p_idx=0;p_idx<currentStruct.len;p_idx++)
    //   printf("%f ",currentStruct.penalties[p_idx]);
    //printf("\n");
 }
 
-bool result_align(Pre_score* matrices[], int z, int est_len, int dna_len, int* result_length_ptr, char* est, char* dna, double* prb, double* chastity, int* s_align, int* e_align, int* mparam, double* alignmentscores, int* max_score_positions, penalty_struct* qparam)
+bool result_align(Pre_score* matrices[], int z, int est_len, int dna_len, int* result_length_ptr, char* est, char* dna, double* prb, double* chastity, int* s_align, int* e_align, int* mparam, double* alignmentscores, int* max_score_positions, penalty_struct* qparam, mode currentMode)
 //Gibt 1 zurueck, wenn kein weiterer Pfad existiert, sonst 0
 
 {  
@@ -71,6 +76,7 @@ bool result_align(Pre_score* matrices[], int z, int est_len, int dna_len, int* r
 
   double prbnum ;
   double chastitynum ;
+
   
   int dna_pos ;
   int est_pos ;
@@ -112,6 +118,7 @@ bool result_align(Pre_score* matrices[], int z, int est_len, int dna_len, int* r
   splice_state = 0 ; //exon
   est_state = 1 ; //[112222111 ... for exons of length 2,4,3 ...
 
+
   //compute length of alignment (between max(m,n) and m+n)     
   //local_alignment: backtracking to first Zero: look at actual value
   while(!(((Pre_score*)matrices[act_m_no] + i*dna_len +j)->value)==0.0) {
@@ -134,9 +141,12 @@ bool result_align(Pre_score* matrices[], int z, int est_len, int dna_len, int* r
 
       prbnum = prb[i-1];
       chastitynum = chastity[i-1];
-      increaseFeatureCount(qparam,dnanum,estnum,prbnum);
 
-      mparam[mlen*dnanum +estnum] ++ ;
+      if(currentMode == USE_QUALITY_SCORES)
+         increaseFeatureCount(qparam,dnanum,estnum,prbnum);
+      else
+         mparam[mlen*dnanum+estnum] ++ ;
+
     }
     else if ((prev_i == (i)) && (prev_j == (j-1))) {// gap on EST
       (*result_length_ptr)= (*result_length_ptr) + 1;
@@ -159,9 +169,11 @@ bool result_align(Pre_score* matrices[], int z, int est_len, int dna_len, int* r
 
       prbnum = prb[i-1];
       chastitynum = chastity[i-1];
-      increaseFeatureCount(qparam,dnanum,estnum,prbnum);
 
-      mparam[mlen*dnanum +estnum] ++ ;
+      if(currentMode == USE_QUALITY_SCORES)
+         increaseFeatureCount(qparam,dnanum,estnum,prbnum);
+      else
+         mparam[mlen*dnanum +estnum] ++ ;
     }
     else {// splice site
       (*result_length_ptr) =  (*result_length_ptr) + (j - prev_j);
@@ -219,5 +231,3 @@ bool result_align(Pre_score* matrices[], int z, int est_len, int dna_len, int* r
 
   return 0;
 }
-
-