+ changed interface array_function to array_class
[qpalma.git] / QPalmaDP / qpalma_dp.cpp
index 1e86e49..cb6cb8c 100644 (file)
@@ -1,4 +1,6 @@
 #include "qpalma_dp.h"
+#include <cstring>
+using namespace std;
 
 /*[splice_align, est_align, weightMatch, alignmentscores, dnaest] = ...
   myalign([nr_paths], dna, est, {h}, matchmatrix, donor, acceptor, remove_duplicate_scores, ...
    nrhs - Number of inputs (Right Hand Side)
    prhs - Array of pointers to input data. */
 
-  myalign([nr_paths], dna, est, {h}, matchmatrix, donor, acceptor, remove_duplicate_scores, ...
-  print_matrix) */
-
-void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
-{
-  
-  
-
-  /***************************************************************************/
-  // 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 ;
-  }
-
+Alignment::Alignment() {
+      len = 0;
+      limits = 0;
+      penalties = 0;
+      max_len = 0;
+      min_len = 0;
+      cache = 0;
+      enum ETransformType transform = T_LINEAR;
+      id = 0;
+      name = 0;
+      use_svm = 0;
+}
 
-  /***************************************************************************/
-  // variables 
-  /***************************************************************************/
-  const int mlen = 6; // score matrix: length of 6 for "- A C G T N"
-  int arg = 0 ; // counter variable
-  int failure ; // 
+Alignment::~Alignment() {}
 
+void Alignment::setQualityMatrix(double* qMat, int length){
+   qualityMatrix = new double[length];
+   for(int i=0; i<length; i++)
+      qualityMatrix[i] = qMat[i];
+}
+void Alignment::getSpliceAlign(){}
+void Alignment::getEstAlign(){}
+void Alignment::getWeightMatch(){}
+void Alignment::getTotalScores(){}
+void Alignment::getDNAEST(){}
 
-  /***************************************************************************/
-  // 0. nr_paths                              
-  /***************************************************************************/
-  int nr_paths = (int)*mxGetPr(prhs[arg]); // 0. input (double)
-  arg++; // arg -> 1
+void Alignment::myalign(int nr_paths_p, char* dna, int dna_len_p, char* est,
+      int est_len_p, 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) {
 
-  
-  /***************************************************************************/
-  // 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
-  
+   printf("Entering myalign...\n");
 
-  /***************************************************************************/
-  // 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
-  
+   mlen = 6; // score matrix: length of 6 for "- A C G T N"
+   dna_len =  dna_len_p;
+   est_len =  est_len_p;
+   nr_paths = nr_paths_p;
 
+#if 0
   /***************************************************************************/
   // 3. h (cell array)
   /***************************************************************************/
-  penalty_struct* functions;
   mwSize anz_func = 1;
   
   functions = read_penalty_struct_from_cell(prhs[arg], anz_func);
@@ -100,41 +73,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()  
@@ -142,65 +81,55 @@ 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);
+  printf("calling fill_matrix...\n");
+
+  fill_matrix(nr_paths, matrices, est_len, dna_len, est, dna, &h, 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]) ;
-  //}
+  printf("after call to fill_matrix...\n");
   /***************************************************************************/ 
   // 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));
-  double alignmentscores[nr_paths]; //alignment score for each path/matrix
-  
-  
+  printf("before init of splice_align and rest...\n");
+
+  splice_align = new int[(dna_len-1)*nr_paths];
+  est_align = new int[(est_len-1)*nr_paths];
+  mmatrix_param = new int[(mlen*mlen)*nr_paths];
+  alignmentscores = new double[nr_paths]; //alignment score for each path/matrix
+
+  printf("Lengths are %d %d %d...\n",dna_len,est_len,mlen);
+  printf("before memset...\n");
+
   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);
-
+  printf("after memset...\n");
 
-  //mexPrintf("DEBUG: decoding %d paths\n",nr_paths);
   for (int z=0; z<nr_paths; z++) {      
     result_length = 0 ;
     
@@ -211,19 +140,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
@@ -255,67 +180,31 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
     }
     
   } //end of z
-  
-  
 
-  //mexPrintf("DEBUG:Returning arguments\n");
-  int index = 0 ;
+  printf("Leaving myalign...\n");
+}
 
-  //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 
+void Alignment::getAlignmentResults(int* s_align, int* e_align,
+      int* mmatrix_p, double* alignscores) {
 
-  //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
+   printf("Entering getAlignmentResults...\n");
+   uint splice_align_size = (dna_len-1)*nr_paths;
+   uint est_align_size = (est_len-1)*nr_paths;
+   uint mmatrix_param_size = (mlen*mlen)*nr_paths;
+   uint alignmentscores_size = nr_paths; //alignment score for each path/matrix
 
-  assert(nlhs==index+1) ;
+   for(int idx=0; idx<splice_align_size; idx++)
+      s_align[idx] = splice_align[idx];
 
-  //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);
+   for(int idx=0; idx<est_align_size; idx++)
+      e_align[idx] =  est_align[idx];
 
-  return;    
-  
+   for(int idx=0; idx<mmatrix_param_size; idx++)
+      mmatrix_p[idx] = mmatrix_param[idx];
+
+   for(int idx=0; idx<alignmentscores_size; idx++)
+      alignscores[idx] = alignmentscores[idx];
+   
+   printf("Leaving getAlignmentResults...\n");
 }