+ fixded index bug in feature count
[qpalma.git] / QPalmaDP / result_align.cpp
index 1bfe6b1..85778f8 100644 (file)
@@ -5,29 +5,64 @@ using namespace std;
 
 void increaseFeatureCount(penalty_struct* qparam, int dnanum, int estnum, double estprb) {
 
+   //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<qparam->len;idx++) {
-      if (qparam->limits[idx] <= value)
+
+   for (idx=0;idx<currentStruct.len;idx++) {
+      if (currentStruct.limits[idx] <= value)
          Lower++;
    }
 
-   if (Lower == 0)
-         qparam->penalties[0] += 1;
+   if (Lower == 0) {
+         currentStruct.penalties[0] += 1;
+         return;
+   }
 
-   if (Lower == qparam->len)
-         qparam->penalties[qparam->len-1] += 1;
+   if (Lower == (currentStruct.len-1)) {
+         currentStruct.penalties[currentStruct.len-1] += 1;
+         return;
+   }
 
    Lower -= 1;
    int Upper = Lower+1; // x-werte bleiben fest
-   double weightup  = 1.0*(value - qparam->limits[Lower]) / (qparam->limits[Upper] - qparam->limits[Lower]);
-   double weightlow = 1.0*(qparam->limits[Upper] - value) / (qparam->limits[Upper] - qparam->limits[Lower]);
-   qparam->penalties[Upper] = qparam->penalties[Upper] + weightup;
-   qparam->penalties[Lower] = qparam->penalties[Lower] + weightlow;
+
+   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[currentPos] = currentStruct;
+
+   //printf("after\n");
+   //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");
 }
 
-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
 
 {  
@@ -41,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 ;
@@ -82,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) {
@@ -104,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;
@@ -129,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);
@@ -189,5 +231,3 @@ bool result_align(Pre_score* matrices[], int z, int est_len, int dna_len, int* r
 
   return 0;
 }
-
-