+ minor fixes
authorfabio <fabio@e1793c9e-67f9-0310-80fc-b846ff1f7b36>
Thu, 17 Jan 2008 16:30:44 +0000 (16:30 +0000)
committerfabio <fabio@e1793c9e-67f9-0310-80fc-b846ff1f7b36>
Thu, 17 Jan 2008 16:30:44 +0000 (16:30 +0000)
git-svn-id: http://svn.tuebingen.mpg.de/ag-raetsch/projects/QPalma@7499 e1793c9e-67f9-0310-80fc-b846ff1f7b36

QPalmaDP/fill_matrix.cpp
QPalmaDP/qpalma_dp.cpp
QPalmaDP/qpalma_dp.h
QPalmaDP/result_align.cpp
python/qpalma.py

index e226266..81f46a3 100644 (file)
@@ -74,8 +74,8 @@ int check_char(char base)
 
 double getScore(struct penalty_struct* qualityScores, int mlen, int dnaChar, int estChar, double baseScore) {
    double score = 0.0;
-   FA(estChar > 0);
-   FA(dnaChar > -1);
+   FA(estChar > 0 && estChar < 6);
+   FA(dnaChar > -1 && dnaChar < 6);
 
    //printf("Starting scoring / scheme is USE_QUALITY_SCORES\n");
    //printf("estChar,dnaChar are: %d,%d,%d",estChar,dnaChar);
@@ -88,7 +88,7 @@ double getScore(struct penalty_struct* qualityScores, int mlen, int dnaChar, int
 
        INT idx = 0 ;
        for (INT i=0; i<currentPen.len; i++)
-               if (currentPen.limits[i]<=baseScore)
+               if (currentPen.limits[i] <= baseScore)
                        idx++ ;
        
        if (idx==0)
@@ -171,7 +171,6 @@ void fill_matrix(int nr_paths, Pre_score* matrices[], int est_len /*counter var:
     }
   }
 
-       
   /*********************************************************************************************/      
   // Element (0,0) for first/best matrice (z=0)
   /*********************************************************************************************/       
index 23513cf..931b46c 100644 (file)
@@ -28,11 +28,13 @@ Alignment::Alignment(int numQPlifs, int numq, bool use_qscores) {
       mlen = 6; // score matrix: length of 6 for "- A C G T N"
 
       numQualSuppPoints = numq;
-      totalNumPlifs = numQPlifs;
+      numPlifs = numQPlifs;
       use_quality_scores = use_qscores;
 
+      printf("number of support points: %d\n",numQualSuppPoints);
+      printf("number of plifs: %d\n",numPlifs );
       FA( numQualSuppPoints >= 0 );
-      FA( totalNumPlifs >= 0 );
+      FA( numPlifs >= 0 );
 }
 
 void Alignment::getDNAEST(){}
@@ -84,10 +86,20 @@ void Alignment::myalign(int nr_paths_p, char* dna, int dna_len_p, char* est,
   for (int z=0; z<nr_paths; z++) {
       matrices[z] = new Pre_score[dna_len * est_len];
   }
+
+   //printf("qualityScores\n");
+   //for(int qidx=0;qidx<numPlifs;qidx++) {
+   //   penalty_struct p; 
+   //   p = qualityScores[qidx];
+   //   printf("%d ",qidx);
+   //   for(int pidx=0;pidx<p.len;pidx++)
+   //      printf("%f ",pidx,p.penalties[pidx]);
+   //   printf("\n");
+   //}
   
-  printf("calling fill_matrix...\n");
+  //printf("calling fill_matrix...\n");
   fill_matrix(nr_paths, matrices, est_len, dna_len, est, dna, prb, &h, matchmatrix, qualityScores, donor, acceptor, remove_duplicate_scores, nr_donor_sites, donor_sites, max_score_positions,currentMode);
-  printf("after call to fill_matrix...\n");
+  //printf("after call to fill_matrix...\n");
 
   /***************************************************************************/ 
   // return arguments etc.
@@ -110,26 +122,13 @@ void Alignment::myalign(int nr_paths_p, char* dna, int dna_len_p, char* est,
    }
 
    alignmentscores_size = nr_paths; //alignment score for each path/matrix
-   qScores_size = totalNumPlifs*nr_paths; //alignment score for each path/matrix
-
-   splice_align = new int[splice_align_size];
-   est_align = new int[est_align_size];
-   mmatrix_param = new int[mmatrix_param_size];
-   alignmentscores = new double[nr_paths]; //alignment score for each path/matrix
-   qualityScoresAllPaths = new penalty_struct[qScores_size];
-
-  int qidx, pidx;
-  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%totalNumPlifs].limits[pidx];
-
-   p.penalties = (REAL*) calloc(p.len,sizeof(REAL));
-   qualityScoresAllPaths[qidx] = p;
-  }
+   numPathsPlifs = numPlifs*nr_paths; //alignment score for each path/matrix
+
+   splice_align         = new int[splice_align_size];
+   est_align            = new int[est_align_size];
+   mmatrix_param        = new int[mmatrix_param_size];
+   alignmentscores      = new double[nr_paths]; //alignment score for each path/matrix
+
 
   // printf("before memset...\n");
   memset((char*)splice_align, -1, (dna_len-1)*nr_paths*sizeof(int)); // fills splice_align with zeros
@@ -137,6 +136,9 @@ void Alignment::myalign(int nr_paths_p, char* dna, int dna_len_p, char* est,
   memset((char*)mmatrix_param, 0, mmatrix_size*nr_paths*sizeof(int)); //fills mmatrix_param with zeros
   memset(alignmentscores, -1, nr_paths*sizeof(double)); //fills alignmentscores with zeros
   //printf("after memset...\n");
+  
+
+  qualityScoresAllPaths= new penalty_struct*[nr_paths];
 
   for (int z=0; z<nr_paths; z++) {
     result_length = 0 ;
@@ -144,17 +146,35 @@ 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 + mmatrix_size*z; //pointer
-    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, currentMode);
-    printf("after call to result_align...\n");
+   qualityScoresAllPaths[z] = new penalty_struct[numPlifs];
+
+   int qidx, pidx;
+   for(qidx=0;qidx<numPlifs;qidx++) {
+      penalty_struct p;
+      p.len = numQualSuppPoints;
+      p.limits = (REAL*) calloc(p.len,sizeof(REAL));
 
-   //for(qidx=0;qidx<totalNumPlifs*nr_paths;qidx++) {
+      for(pidx=0;pidx<p.len;pidx++)
+         p.limits[pidx] = qualityScores[qidx%numPlifs].limits[pidx];
+
+      p.penalties = (REAL*) calloc(p.len,sizeof(REAL));
+      qualityScoresAllPaths[z][qidx] = p;
+   }
+
+   //penalty_struct* qparam = qualityScoresAllPaths + (numPlifs*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, qualityScoresAllPaths[z] , currentMode);
+    //printf("after call to result_align...\n");
+
+   //printf("z is %d\n",z);
+   //for(qidx=0;qidx<numPlifs;qidx++) {
    //   penalty_struct p;
-   //   p = qualityScoresAllPaths[qidx];
+   //   p = qualityScoresAllPaths[qidx+numPlifs*z];
+   //   printf("%d ",qidx);
    //   for(pidx=0;pidx<p.len;pidx++)
-   //      printf("%f ",p.penalties[pidx]);
+   //      printf("%f ",pidx,p.penalties[pidx]);
    //   printf("\n");
    //}
 
@@ -230,16 +250,26 @@ void Alignment::getAlignmentResults(int* s_align, int* e_align,
       alignscores[idx] = alignmentscores[idx];
    
    penalty_struct currentPlif;
-   int pidx;
    int ctr=0;
-   for(idx=0; idx<qScores_size; idx++) {
-      currentPlif = qualityScoresAllPaths[idx];
-      //printf("Size is %dx%d\n",qScores_size,currentPlif.len);
-      for(pidx=0; pidx<currentPlif.len; pidx++) {
-         qScores[ctr] = currentPlif.penalties[pidx];
-         ctr++;
-      }
-   }
+   for (int z=0; z<nr_paths; z++) {
+      for(int estChar=1;estChar<6;estChar++) {
+         for(int dnaChar=0;dnaChar<6;dnaChar++) {
+            int currentPos = (estChar-1)*6+dnaChar;
+            currentPlif = qualityScoresAllPaths[z][currentPos];
+            for(int pidx=0; pidx<currentPlif.len; pidx++) {
+               qScores[ctr] = currentPlif.penalties[pidx];
+               ctr++;
+            }
+   }}}
+
+   //for(idx=0; idx<numPathsPlifs; idx++) {
+   //   currentPlif = qualityScoresAllPaths[idx];
+   //   //printf("Size is %dx%d\n",numPathsPlifs,currentPlif.len);
+   //   for(pidx=0; pidx<currentPlif.len; pidx++) {
+   //      qScores[ctr] = currentPlif.penalties[pidx];
+   //      ctr++;
+   //   }
+   //}
 
    //printf("\nctr is %d\n",ctr);
    //printf("Leaving getAlignmentResults...\n");
index 070b6c9..4049221 100644 (file)
@@ -71,7 +71,7 @@ class Alignment {
       int* est_align;
       int* mmatrix_param;
       double* alignmentscores;
-      struct penalty_struct* qualityScoresAllPaths;
+      struct penalty_struct** qualityScoresAllPaths;
 
       int dna_len;
       int est_len;
@@ -79,14 +79,14 @@ class Alignment {
       int nr_paths;
 
       int numQualSuppPoints;
-      int totalNumPlifs;
+      int numPlifs;
       bool use_quality_scores;
 
       uint splice_align_size ;
       uint est_align_size ;
       uint mmatrix_param_size ;
       uint alignmentscores_size ;
-      uint qScores_size ;
+      uint numPathsPlifs ;
 
       INT len;
       REAL *limits;
index e1f92f7..7994419 100644 (file)
@@ -6,13 +6,14 @@ using namespace std;
 
 void increaseFeatureCount(penalty_struct* qparam, int dnaChar, int estChar, double estprb) {
 
-   FA(estChar > 0);
-   FA(dnaChar > -1);
+   FA(estChar > 0 && estChar < 6);
+   FA(dnaChar > -1 && dnaChar < 6);
 
    int currentPos = (estChar-1)*6+dnaChar;
+
    penalty_struct currentStruct = qparam[currentPos];
 
-   printf("Current index %d est/dna: %d %d\n",currentPos,estChar,dnaChar);
+   printf("Current index %d est/dna: %d %d score %f\n",currentPos,estChar,dnaChar,estprb);
 
    //printf("before\n");
    //int p_idx;
@@ -41,7 +42,7 @@ void increaseFeatureCount(penalty_struct* qparam, int dnaChar, int estChar, doub
          return;
    }
 
-   if (Lower == (currentStruct.len)) {
+   if (Lower == currentStruct.len) {
          currentStruct.penalties[currentStruct.len-1] += 1;
          qparam[currentPos] = currentStruct;
          return;
index 1be8ea1..22c623c 100644 (file)
@@ -204,10 +204,9 @@ class QPalma:
             dna_len = len(dna)
             est_len = len(est)
 
-            dna ='t'*dna_len
-            est ='a'*est_len
-            #quality = [40]*len(quality)
-
+            #dna ='t'*dna_len
+            #est ='t'*est_len
+            #quality = [math.fabs(elem) for elem in quality]
             #mmatrix = zeros((6,1))
             #for pos,elem in enumerate(qualityPlifs):
             #   elem.penalties = [0.0]*len(elem.penalties)
@@ -228,6 +227,10 @@ class QPalma:
             # Create the alignment object representing the interface to the C/C++ code.
             currentAlignment = QPalmaDP.Alignment(Configuration.numQualPlifs,Configuration.numQualSuppPoints, use_quality_scores)
 
+            #for pos,elem in enumerate(qualityPlifs):
+            #   elem.penalties = [1.0]*len(elem.penalties)
+            #   qualityPlifs[pos] = elem
+
             c_qualityPlifs = QPalmaDP.createPenaltyArrayFromList([elem.convert2SWIG() for elem in qualityPlifs])
 
             # calculates SpliceAlign, EstAlign, weightMatch, Gesamtscores, dnaest
@@ -321,18 +324,20 @@ class QPalma:
                dpen = mat(d.penalties).reshape(len(d.penalties),1)
                apen = mat(a.penalties).reshape(len(a.penalties),1)
 
-               totalQualityPenalties = param[-totalQualSP:]
+               #totalQualityPenalties = ones((Configuration.totalQualSuppPoints,1))
                features = numpy.vstack([hpen, dpen, apen, mmatrix[:], totalQualityPenalties])
 
                AlignmentScores[pathNr+1] = (allWeights[:,pathNr+1].T * features)[0,0]
 
-               pdb.set_trace()
+               #pdb.set_trace()
 
                # Check wether scalar product + loss equals viterbi score
                #assert math.fabs(AlignmentScores[pathNr] - AlignmentScores[pathNr+1]) < 1e-6,\
                #'Scalar prod + loss is not equal Viterbi score. Respective values are %f, %f' % \
                #(AlignmentScores[pathNr],AlignmentScores[pathNr+1])
                print '%f vs %f' % (newDPScores[pathNr], AlignmentScores[pathNr+1] + path_loss[pathNr])
+               assert math.fabs(newDPScores[pathNr] - (AlignmentScores[pathNr+1]\
+               + path_loss[pathNr])) <= 1e-5
 
                #  # if the pathNr-best alignment is very close to the true alignment consider it as true
                if norm( allWeights[:,0] - allWeights[:,pathNr+1] ) < 1e-5: