+ rewrote C interface for SWIG/Python
authorfabio <fabio@e1793c9e-67f9-0310-80fc-b846ff1f7b36>
Mon, 3 Dec 2007 12:19:09 +0000 (12:19 +0000)
committerfabio <fabio@e1793c9e-67f9-0310-80fc-b846ff1f7b36>
Mon, 3 Dec 2007 12:19:09 +0000 (12:19 +0000)
+ extended qpalma algorithm
+ added evalutation script for alignments
TODO

- added preprocessing scripts
- create an evaluation score

git-svn-id: http://svn.tuebingen.mpg.de/ag-raetsch/projects/QPalma@6915 e1793c9e-67f9-0310-80fc-b846ff1f7b36

QPalmaDP/QPalmaDP.i
QPalmaDP/qpalma_dp.cpp
QPalmaDP/qpalma_dp.h
python/.qpalma.py.swp
python/computeSpliceWeights.py
python/qpalma.py
tools/calculateAlignmentQuality.m [new file with mode: 0644]

index a592bb0..b0e8647 100644 (file)
@@ -17,7 +17,7 @@
 %include "qpalma_dp.h"
 
 %array_functions(int, intArray) ;
-%array_function(double, doubleArray) ;
+%array_functions(double, doubleArray) ;
 %array_class(Pre_score, Pre_scoreArray) ;
 
 %pythoncode
index 5aeaa80..171aa05 100644 (file)
@@ -10,7 +10,6 @@
    prhs - Array of pointers to input data. */
 
 Alignment::Alignment() {
-
       len = 0;
       limits = 0;
       penalties = 0;
@@ -25,39 +24,19 @@ Alignment::Alignment() {
 
 Alignment::~Alignment() {}
 
-void Alignment::myalign(int nr_paths, char* dna, int dna_len, char* est,
-      int est_len, struct penalty_struct h, double* matchmatrix, int mm_len,
-      double* donor, int d_len, double* acceptor, int a_len,
-      bool remove_duplicate_scores, bool print_matrix) {}
-
+void Alignment::setQualityMatrix(double* qualityMat, int length){}
 void Alignment::getSpliceAlign(){}
 void Alignment::getEstAlign(){}
 void Alignment::getWeightMatch(){}
 void Alignment::getTotalScores(){}
 void Alignment::getDNAEST(){}
 
-#if 0
-void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
-{
-  
-  
+void Alignment::myalign(int nr_paths, char* dna, int dna_len, char* est,
+      int est_len, struct penalty_struct h, double* matchmatrix, int mm_len,
+      double* donor, int d_len, double* acceptor, int a_len,
+      bool remove_duplicate_scores, bool print_matrix) {
 
-  /***************************************************************************/
-  // right number of arguments? 
-  /***************************************************************************/
-  if(nrhs != 9)
-  {
-    mexPrintf("10 input arguments please: \n [splice_align, est_align, weightMatch, alignmentscores, dnaest] = myalign(nr_paths, dna, est, cellarray, matchmatrix, donor, acceptor, remove_duplicate_scores, print_matrix)\n");
-    return ;
-  }
   
-  if(nlhs != 5)
-  {
-    mexPrintf("5 output arguments please\n [splice_align, est_align, weightMatch, alignmentscores, dnaest] = myalign(nr_paths, dna, est, cellarray, matchmatrix, donor, acceptor, remove_duplicate_scores, print_matrix)\n");
-    return ;
-  }
-
-
   /***************************************************************************/
   // variables 
   /***************************************************************************/
@@ -65,46 +44,11 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
   int arg = 0 ; // counter variable
   int failure ; // 
 
-
-  /***************************************************************************/
-  // 0. nr_paths                              
-  /***************************************************************************/
-  int nr_paths = (int)*mxGetPr(prhs[arg]); // 0. input (double)
-  arg++; // arg -> 1
-
-  
-  /***************************************************************************/
-  // 1. dna
-  /***************************************************************************/
-  int dna_len = mxGetN(prhs[arg]) + 1; // mxGetN: number of columns (int)
-  char* dna = (char*)mxCalloc(dna_len, sizeof(char)); 
-  
-  failure = mxGetString(prhs[arg], dna, dna_len); // NULL at end of string
-  if (failure) { // mxGetString returns 0 if success, 1 on failure
-    mexPrintf("couldn't open dna");
-    return;
-  }
-  arg++; // arg -> 2
-  
-
-  /***************************************************************************/
-  // 2. est
-  /***************************************************************************/
-  int est_len = mxGetN(prhs[arg]) + 1;
-  char* est = (char*)mxCalloc(est_len, sizeof(char));
-  
-  failure = mxGetString(prhs[arg], est, est_len);
-  if (failure) {
-    mexPrintf("couldn't open est");
-    return;
-  }
-  arg++; // arg -> 3
-  
-
+  penalty_struct* functions;
+#if 0
   /***************************************************************************/
   // 3. h (cell array)
   /***************************************************************************/
-  penalty_struct* functions;
   mwSize anz_func = 1;
   
   functions = read_penalty_struct_from_cell(prhs[arg], anz_func);
@@ -125,41 +69,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
   //  mexPrintf("entry %d: <%f>\n", i, matchmatrix[i]);
   //}
   arg++; // arg -> 5
-  
-
-  /***************************************************************************/
-  // 5. donor
-  /***************************************************************************/
-  int donor_length = mxGetN(prhs[arg]);
-  double *donor;
-  donor = mxGetPr(prhs[arg]);
-  //mexPrintf("Donor eingelesen: <%f>\n", donor[0]);
-  arg++; // arg -> 6
-  
-
-  /***************************************************************************/
-  // 6. acceptor
-  /***************************************************************************/
-  int acceptor_length = mxGetN(prhs[arg]);
-  double *acceptor;
-  acceptor = mxGetPr(prhs[arg]);
-  //mexPrintf("Acceptor eingelesen: <%f>\n", acceptor[0]);
-  arg++; // arg -> 7
-
-  
-  /***************************************************************************/
-  // 7.remove_duplicate_scores
-  /***************************************************************************/
-  bool remove_duplicate_scores = (bool)*mxGetPr(prhs[arg]);
-  arg++; // arg -> 8
-  
-
-  /***************************************************************************/
-  // 8. print_matrix
-  /***************************************************************************/
-  int print_matrix = (int)*mxGetPr(prhs[arg]);
-  arg++; // arg -> 9
-
+#endif
 
   /***************************************************************************/ 
   // initialize alignment matrices and call fill_matrix()  
@@ -167,65 +77,45 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
   
   //possible donor positions
   int nr_donor_sites = 0 ;
-  for (int ii=0; ii<donor_length; ii++) {
+  for (int ii=0; ii<d_len; ii++) {
     if(isfinite(donor[ii])) {
       nr_donor_sites++ ;
     }
   }
-  //mwSize donor_sites[nr_donor_sites] ;
-  mwSize* donor_sites = (mwSize*)mxCalloc(nr_donor_sites,sizeof(mwSize));
+
+  int* donor_sites = new int[nr_donor_sites];
   int donor_idx = 0 ;
-  for (int ii=0; ii<donor_length; ii++) {
+  for (int ii=0; ii<d_len; ii++) {
     if(isfinite(donor[ii])) {
       donor_sites[donor_idx] = ii+1 ;
       donor_idx++ ;
     }  
   }
 
+  int* max_score_positions = new int[nr_paths*2];
 
-  //int max_score_positions[nr_paths*2]; //local alignment: where is the highest alignment score [i1,j1,i2,j2...]
-  mwSize* max_score_positions = (mwSize*)mxCalloc(nr_paths*2,sizeof(mwSize));
-
-  //Pre_score* matrices[nr_paths];
-  Pre_score** matrices = (Pre_score**)mxCalloc(nr_paths,sizeof(Pre_score*));
+  Pre_score** matrices = new Pre_score*[nr_paths];
   for (int z=0; z<nr_paths; z++) {
-      //matrices[z] = new Pre_score[dna_len * est_len];
-      //mexPrintf("size of one matrice: %d\n", est_len*dna_len*sizeof(Pre_score)) ;
-      matrices[z] = (Pre_score*)mxCalloc(dna_len*est_len,sizeof(Pre_score));
+      matrices[z] = new Pre_score[dna_len * est_len];
   }
   
   fill_matrix(nr_paths, matrices, est_len, dna_len, est, dna, functions, matchmatrix, donor, acceptor, remove_duplicate_scores, nr_donor_sites, donor_sites, max_score_positions);
   
-  //just some info
-  //for (int z=0; z<nr_paths; z++) {
-  //  mexPrintf("(%d, %d)\n", max_score_positions[2*z], max_score_positions[2*z +1]) ;
-  //}
   /***************************************************************************/ 
   // return arguments etc.
   /***************************************************************************/
   int result_length; //Eine Variable fuer alle Matrizen
 
-  //int splice_align[(dna_len-1)*nr_paths]; //SpliceAlign for DNA - 1 line per matrix
-  //int est_align[(est_len-1)*nr_paths]; //SpliceAlign for EST - 1 line per matrix
-  //int mmatrix_param[(mlen*mlen)*nr_paths];//Matrix fuer Scorematrixparameter - 1 Zeile pro Matrix
-  int* splice_align = (int*)mxCalloc((dna_len-1)*nr_paths,sizeof(int));
-  int* est_align = (int*)mxCalloc((est_len-1)*nr_paths,sizeof(int));
-  int* mmatrix_param = (int*)mxCalloc((mlen*mlen)*nr_paths,sizeof(int));
+  int* splice_align = new int[(dna_len-1)*nr_paths];
+  int* est_align = new int[(est_len-1)*nr_paths];
+  int* mmatrix_param = new int[(mlen*mlen)*nr_paths];
   double alignmentscores[nr_paths]; //alignment score for each path/matrix
-  
-  
-  memset((char*)splice_align, -1, (dna_len-1)*nr_paths*sizeof(int)); // fills splice_align with zeros
-  memset((char*)est_align, -1, (est_len-1)*nr_paths*sizeof(int)); // fills est_align with zeros
-  memset((char*)mmatrix_param, 0, mlen*mlen*nr_paths*sizeof(int)); //fills mmatrix_param with zeros
-  memset(alignmentscores, -1, nr_paths*sizeof(double)); //fills alignmentscores with zeros
-  
-  // dnaest
-  plhs[nlhs-1] = mxCreateCellMatrix(2, nr_paths);
-
 
-  //mexPrintf("DEBUG: decoding %d paths\n",nr_paths);
+  //memset((char*)splice_align, -1, (dna_len-1)*nr_paths*sizeof(int)); // fills splice_align with zeros
+  //memset((char*)est_align, -1, (est_len-1)*nr_paths*sizeof(int)); // fills est_align with zeros
+  //memset((char*)mmatrix_param, 0, mlen*mlen*nr_paths*sizeof(int)); //fills mmatrix_param with zeros
+  //memset(alignmentscores, -1, nr_paths*sizeof(double)); //fills alignmentscores with zeros
+  
   for (int z=0; z<nr_paths; z++) {      
     result_length = 0 ;
     
@@ -236,19 +126,15 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
     bool no_more_path = result_align(matrices, z, est_len, dna_len, &result_length, est, dna, s_align, e_align, mparam, alignmentscores, max_score_positions);
       
     // dnaest
-    mxArray* mxdna = mxCreateDoubleMatrix(result_length, 1, mxREAL);
-    mxArray* mxest = mxCreateDoubleMatrix(result_length, 1, mxREAL);
-    double *DNA = mxGetPr(mxdna) ;
-    double *EST = mxGetPr(mxest) ;
-    mxSetCell(plhs[nlhs-1], z*2, mxdna) ;
-    mxSetCell(plhs[nlhs-1], z*2+1, mxest) ;
+    double *DNA = new double[result_length];
+    double *EST = new double[result_length];
     
     //backtracking
-    mwSize i = max_score_positions[2*z] ; //i (est)
-    mwSize j = max_score_positions[2*z +1] ; //j (dna)
-    mwSize prev_i = ((Pre_score*)matrices[z] + i*dna_len +j)->prev_i; 
-    mwSize prev_j = ((Pre_score*)matrices[z] + i*dna_len +j)->prev_j;
-    mwSize prev_m_no = ((Pre_score*)matrices[z] + i*dna_len +j)->prev_matrix_no;
+    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
@@ -280,67 +166,4 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
     }
     
   } //end of z
-  
-  
-
-  //mexPrintf("DEBUG:Returning arguments\n");
-  int index = 0 ;
-
-  //splice_align (mxArray):
-  const mwSize dims0[] = {(dna_len-1), nr_paths};
-  //const int dims0[] = {(dna_len-1), nr_paths};
-  plhs[index] = mxCreateNumericArray(2, dims0, mxINT32_CLASS, mxREAL);
-  mwSize* start0 = (mwSize*)(mxGetData(plhs[index]));
-  memcpy(start0, splice_align, mxGetM(plhs[index])*mxGetN(plhs[index])*mxGetElementSize(plhs[index]));
-  index++ ; // ->1
-  //est_align (mxArray):
-  const mwSize dims1[] = {(est_len-1), nr_paths};
-  //const int dims1[] = {(est_len-1), nr_paths};
-  plhs[index] = mxCreateNumericArray(2, dims1, mxINT32_CLASS, mxREAL);
-  mwSize* start1 = (mwSize*)(mxGetData(plhs[index]));
-  memcpy(start1, est_align, mxGetM(plhs[index])*mxGetN(plhs[index])*mxGetElementSize(plhs[index]));
-  index++ ; // ->2 
-
-  //weightmatch (mxArray)
-  const mwSize dims2[] = {mlen*mlen, nr_paths};
-  //const int dims2[] = {mlen*mlen, nr_paths};
-  plhs[index] = mxCreateNumericArray(2, dims2, mxINT32_CLASS, mxREAL);
-  mwSize* start2= (mwSize*)(mxGetData(plhs[index]));
-  memcpy(start2, mmatrix_param, mxGetM(plhs[index])*mxGetN(plhs[index])*mxGetElementSize(plhs[index]));
-  index++ ; // ->3
-  
-  //alignmentscore 
-  plhs[index] = mxCreateDoubleMatrix(1,nr_paths,mxREAL);
-  if (plhs[index] == NULL)
-  {
-    printf("%s : Out of memory on line %d\n", __FILE__, __LINE__); 
-    printf("Unable to create mxArray.\n");
-    exit(1);
-  }
-  double* start3 = (double*)(mxGetData(plhs[index]));
-  memcpy(start3, alignmentscores, mxGetM(plhs[index])*mxGetN(plhs[index])*mxGetElementSize(plhs[index]));
-  index++ ; // ->4
-
-  assert(nlhs==index+1) ;
-
-  //Free Memory
-  //mexPrintf("DEBUG:Free Memory\n");
-  mxFree(donor_sites);
-  mxFree(max_score_positions);
-  for (int i=nr_paths-1;i>=0; i--)
-  {
-      //delete[] matrices[i];
-      mxFree(matrices[i]);
-  }
-  mxFree(matrices);
-  mxFree(splice_align);
-  mxFree(est_align);
-  mxFree(mmatrix_param);
-  mxFree(dna);
-  mxFree(est);
-
-  return;    
-  
 }
-#endif
index 5b3a4f7..0825904 100644 (file)
@@ -91,7 +91,7 @@ class Alignment {
             limits[i] = lts[i];
       }
       
-
+      void setQualityMatrix(double* qualityMat, int length);
       void getSpliceAlign();
       void getEstAlign();
       void getWeightMatch();
index b87f9e8..88f30fa 100644 (file)
Binary files a/python/.qpalma.py.swp and b/python/.qpalma.py.swp differ
index 51c6e9f..762b003 100644 (file)
@@ -72,7 +72,7 @@ def computeSpliceWeights(d, a, h, SpliceAlign, don_supp, acc_supp):
    #assert(all(AcceptorScores~=-Inf)) ;
 
    #Den Vektor Acceptorstellen durchgehen und die Gewichtsvektoren belasten:
-   AcceptorScores = [elem for pos,elem in enumerate(acc_supp) if SpliceAlign[pos] == 2]
+   AcceptorScores = [elem for pos,elem in enumerate(acc_supp) if SpliceAlign[pos+1] == 2]
    #assert not ( -inf in AcceptorScores )
 
    weightAcc = calculateWeights(a,AcceptorScores)
index 2c70418..ef82edd 100644 (file)
@@ -12,6 +12,7 @@ from computeSpliceAlign import *
 
 from computeSpliceWeights import *
 import QPalmaDP
+from SIQP_CPX import SIQPSolver
 
 from penalty_lookup import *
 from compute_donacc import *
@@ -19,21 +20,14 @@ from compute_donacc import *
 """
 A training method for the QPalma project
 
-Overall :
+Overall procedure:
  
-   1. Load the data 
-   via paths_load_data
-
+   1. Load the data -> via paths_load_data
    2. Create a QP using SIQP (paths_create_qp)
-   
    3. Set initial params using set_param_palma -> palma_utils
-
    4. computeSpliceAlign
-
    5. compute_SpliceAlignLocal
-
    6. computeSpliceWeights
-
    7. myalign_local
 
 """
@@ -94,7 +88,6 @@ def run():
    random_N = 100 ; # number of constraints to generate per iteration
    iteration_steps = 200 ; #upper bound on iteration steps
 
-   # myalign
    remove_duplicate_scores = 0
    anzpath = 2
    print_matrix = 0
@@ -114,54 +107,33 @@ def run():
    #############################################################################################
    # Training
    #############################################################################################
-   print 'Starting training...'
-
-      #The alignment matrix M looks like 
-
-      #        DNA
-
-      #      A C G T -
-      #   ------------
-      #   A
-      # 
-      #E  C
-      #S
-      #T  G
-
-      #   T
-
-      #   - 
-
-
-      #where 
-
-      #   - is a gap, 
+   plog(logfh,'Starting training...')
 
+   numExamples = N
+   C=1.0
 
    numDonorSupportPoints         = 30
    numAcceptorSupportPoints      = 30
    numIntronLengthSupportPoints  = 30 
-   numQualityScoreSupportPoints  =  0
+   sizeMMatrix                   = 36
+   numQualityScoreSupportPoints  = 16*0
 
    numFeatures = numDonorSupportPoints\
                + numAcceptorSupportPoints\
                + numIntronLengthSupportPoints\
+               + sizeMMatrix\
                + numQualityScoreSupportPoints
 
-   numExamples = N
-
-   from SIQP_CPX import SIQPSolver
-   C=1.0
+   qualityMatrix = [0.0]*numQualityScoreSupportPoints
    
-   logfile = open('qpalma.log','w+')
+   plog(logfh,'Initializing problem...\n')
    #problem = SIQPSolver(numFeatures,numExamples,C,logfile)
-
-   #xis = zeros(1,N) ; #ninitialize slack variables
    #num_path = anzpath*ones(1,N) ; # nr of alignments done (best path, second-best path etc.)
    #gap = zeros(1,N) ; %? sum of differences between true alignment score and 'false' alignment scores
+   #for step_nr=1:iteration_steps,
 
+   plog(logfh,'Starting training...\n')
 
-   #for step_nr=1:iteration_steps,
    iteration_nr = 1
 
    while True:
@@ -170,6 +142,7 @@ def run():
       for exampleId in range(numExamples):
          if exampleId % 1000 == 0:
             print 'Current example nr %d' % exampleId
+
          dna = Sequences[exampleId] 
          est = Ests[exampleId] 
 
@@ -181,12 +154,25 @@ def run():
          # Berechne die Parameter des wirklichen Alignments (but with untrained d,a,h ...)    
          trueSpliceAlign, trueWeightMatch = computeSpliceAlign(dna, exons)
 
+         # Calculate
          trueWeightDon, trueWeightAcc, trueWeightIntron = computeSpliceWeights(d, a, h, trueSpliceAlign, don_supp, acc_supp)
 
+         # Reshape currentW param 
+         currentW = zeros((numFeatures,1))
+         currentW[0:numDonorSupportPoints,1] = trueWeigthDon[:]
+         currentW[0:numAcceptorSupportPoints,1] = trueWeigthAcc[:]
+         currentW[0:numIntronLengthSupportPoints,1] = trueWeigthIntron[:]
+         currentW[0:numDonorSupportPoints,1] = trueWeigthMatch[:]
+
+         #currentPhi = d.penalties a h mmatrix
+
+         # Calculate w'phi(x,y) the total score of the alignment
+         #alignmentScore = currentW.T * currentPhi
+
+         #
+         # Calculate wrong alignments
          #
-         # compute wrong alignments
-         #  
-         nr_paths = 1
+         nr_paths = 2
          dna = 'acgtagct'
          dna_len = len(dna)
 
@@ -203,6 +189,9 @@ def run():
          print_matrix = False
 
          currentAlignment = QPalmaDP.Alignment()
+
+         qualityMat = QPalmaDP.createDoubleArrayFromList(qualityMatrix)
+         currentAlignment.setQualityMatrix(qualityMat,numQualitySupportPoints)
          ps = h.convert2SWIG()
 
          currentAlignment.myalign( nr_paths, dna, dna_len,\
@@ -232,29 +221,6 @@ def run():
    logfh.close()
 
    """
-     for id = 1:N  
-       %fprintf('%i\n', id) ;
-       if (mod(id,50)==0), fprintf(1,'.'); end
-      
-       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-       % True Alignment and Comparison with wrong ones
-       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-       
-       %Berechne die Parameter des wirklichen Alignments (but with untrained d,a,h ...)    
-       [true_SpliceAlign, true_weightMatch] = compute_SpliceAlign_local(dna, exons) ;
-       [true_weightDon, true_weightAcc, true_weightIntron] = ...
-           compute_SpliceWeights(d, a, h, true_SpliceAlign, don_supp, acc_supp) ;
-       
-       double(true_weightMatch) ;
-       
-       % Berechne Gewichtsmatrix fuer aktuelle id (matrix anlegen)
-       Weights = zeros(num_path(id)+1, 126) ; %first line: true alignment, rest: wrong alignments
-       Weights(1,:) = [true_weightIntron, true_weightDon, true_weightAcc, true_weightMatch] ;
-       
-       %Score of whole alignment
-       AlignmentScores = zeros(1, num_path(id)+1) ; %first entry: true alignment
-       AlignmentScores(1) = Weights(1,:) * [h.penalties' ; d.penalties' ; a.penalties' ; mmatrix(:)] ;
-      
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Wrong Alignments
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
@@ -269,14 +235,10 @@ def run():
            myalign_local(num_path(id), dna, est, {h}, mmatrix, donor, acceptor, ...
                          remove_duplicate_scores, print_matrix);
        
-       %[SpliceAlign, EstAlign, weightMatch, Gesamtscores, dnaest] = myalign_local(1, dna, est, {h}, mmatrix, donor,
-       %acceptor, remove_duplicate_scores, print_matrix);
-       
        %return values (are int32, have to be double later cause we compute scores
        SpliceAlign = double(SpliceAlign') ; %column
        weightMatch = double(weightMatch') ;
 
-
        %test
        for pfadNo = 1:num_path(id)
          assert(sum(weightMatch(pfadNo,7:end)) == sum(SpliceAlign(pfadNo,:)==0)) ;
@@ -288,18 +250,12 @@ def run():
          end
          assert(all(new_weightMatch == weightMatch(pfadNo,:))) ;
          assert(sum(new_weightMatch(7:end)) == sum(SpliceAlign(pfadNo,:)==0)) ;
-         %reshape(new_weightMatch,6,6)
-         %reshape(weightMatch(pfadNo,:),6,6)
        end
        
        %assert that it is the right file 
        assert(all(dna(find(SpliceAlign(1,:)==1)) == 'g')) ;
        
        %just some info
-       %print_align(dnaest,1) ;
-       %exons(3,2)-exons(1,1)
-       %length(dnaest{1,1})
-       %keyboard
        
        %Berechne Gewichte der durch Alignment berechneten Pfade    
        true_map = zeros(1,1+num_path(id)) ;
@@ -363,45 +319,18 @@ def run():
      end       
      fprintf('\n') ;
       #################################################################################
-      const_added = solver.addConstraint(deltas, loss, pos, marginRescaling)
 
+      const_added = solver.addConstraint(deltas, loss, pos, marginRescaling)
+      
       objValue,w,self.slacks = solver.solve()
-     
-     # [param,xis]  = paths_qp_solve(lpenv,Q,f,A,b,LB,UB);
-     # sum_xis = sum(xis)
-     # maxgap=max(gap) 
-     [mean(xis>=1) mean(gap>=1) mean(xis>=1e-3) mean(gap>=1e-3)]
-
-      if max(gap)<=column_eps:
-         break
 
-   assert(all(maxgap<=column_eps))
-   assert(length(take_idx)==N)
+      sum_xis = 0
+      for elem in self.slacks:
+         sum_xis +=  elem
+     
    """
+
    print 'Training completed'
 
 if __name__ == '__main__':
    run()
-
-   """
-      def __init__(self):
-         numAcceptorParams = 0
-         numDonorParams = 0
-         numIntronLengthParams = 30
-         
-         # how many supporting points do we have for the 
-         # quality function per base 
-         numQualityParamsPerBase = 5
-
-         # we have quality params for each possible 
-         # for all possible combinations of 
-         # {A,C,G,T} x {A,C,G,T,-}
-         # 4 x 5
-         numQualityParams = numQualityParamsPerBase * 4 * 5
-
-         # for all possible combinations of 
-         # {A,C,G,T} x {A,C,G,T,-,*} plus
-         # {-} x {A,C,G,T,-}
-         # 1 x 5
-         numParams = numQualityParams + 1*5
-   """
diff --git a/tools/calculateAlignmentQuality.m b/tools/calculateAlignmentQuality.m
new file mode 100644 (file)
index 0000000..a26989b
--- /dev/null
@@ -0,0 +1,53 @@
+load confirmed_sequences.mat;
+testrun = genes;
+load /fml/ag-raetsch/share/projects/altsplicedata/zebrafish/confirmed_sequences.mat;
+ground_truth = genes;
+clear genes;
+
+for i = 1:length(ground_truth)
+
+   currentTru = ground_truth(i);
+   currentTruExons = currentTru.exons{1};
+
+   for j = 1:length(testrun)
+      currentPred = testrun(j);
+
+      if strcmp(currentTru.chr,currentPred.chr) && strcmp(currentTru.strands,currentPred.strands)
+         continue;
+      end
+
+      truStart = currentTru.start;
+      truStop  = currentTru.stop;
+      predStart = currentPred.start;
+      predStop = currentPred.stop;
+
+      if predStop < truStart || predStart > truStop
+         continue;
+      end
+
+      currentPredExons = currentPred.exons{1};
+
+      rows = size(currentTruExons,1);
+      for ii = 1:rows
+         exonStart = currentTruExons(ii,1);
+         exonStop = currentTruExons(ii,2);
+         for jj = 1:size(currentPredExons,1)
+            predExonStart = currentPredExons(jj,1);
+            predExonStop = currentPredExons(jj,2);
+
+            if exonStart >= predExonStart && exonStart <= predExonStop
+               disp('Overlapping');
+            % end of exonth is nested inside predExoniction
+            elseif exonStop >= predExonStart && exonStop <= predExonStop
+               disp('Overlapping');
+            % predExoniction is nested inside exonth
+            elseif exonStart <= predExonStart && predExonStop <= exonStop
+               disp('Overlapping');
+            else
+               d=1;
+            end
+            
+         end
+      end
+   end
+end