+ rewrote C interface for SWIG/Python
[qpalma.git] / QPalmaDP / qpalma_dp.cpp
1 #include "qpalma_dp.h"
2
3 /*[splice_align, est_align, weightMatch, alignmentscores, dnaest] = ...
4 myalign([nr_paths], dna, est, {h}, matchmatrix, donor, acceptor, remove_duplicate_scores, ...
5 print_matrix) */
6
7 /* nlhs - Number of expected mxArrays (Left Hand Side)
8 plhs - Array of pointers to expected outputs
9 nrhs - Number of inputs (Right Hand Side)
10 prhs - Array of pointers to input data. */
11
12 Alignment::Alignment() {
13 len = 0;
14 limits = 0;
15 penalties = 0;
16 max_len = 0;
17 min_len = 0;
18 cache = 0;
19 enum ETransformType transform = T_LINEAR;
20 id = 0;
21 name = 0;
22 use_svm = 0;
23 }
24
25 Alignment::~Alignment() {}
26
27 void Alignment::setQualityMatrix(double* qualityMat, int length){}
28 void Alignment::getSpliceAlign(){}
29 void Alignment::getEstAlign(){}
30 void Alignment::getWeightMatch(){}
31 void Alignment::getTotalScores(){}
32 void Alignment::getDNAEST(){}
33
34 void Alignment::myalign(int nr_paths, char* dna, int dna_len, char* est,
35 int est_len, struct penalty_struct h, double* matchmatrix, int mm_len,
36 double* donor, int d_len, double* acceptor, int a_len,
37 bool remove_duplicate_scores, bool print_matrix) {
38
39
40 /***************************************************************************/
41 // variables
42 /***************************************************************************/
43 const int mlen = 6; // score matrix: length of 6 for "- A C G T N"
44 int arg = 0 ; // counter variable
45 int failure ; //
46
47 penalty_struct* functions;
48 #if 0
49 /***************************************************************************/
50 // 3. h (cell array)
51 /***************************************************************************/
52 mwSize anz_func = 1;
53
54 functions = read_penalty_struct_from_cell(prhs[arg], anz_func);
55 //std::cout << " bla " << functions->max_len << " \n";
56 //std::cout << "lookup_penalty(5) = " << lookup_penalty(functions, 5, NULL, 0)<<std::endl ;
57 arg++; // arg -> 4
58
59
60 /***************************************************************************/
61 // 4. match matrix (columnwise)
62 /***************************************************************************/
63 int nr_entries = mxGetN(prhs[arg]) * mxGetM(prhs[arg]);
64 //mexPrintf("\nNumber of entries in matchmatrix: %d\n", nr_entries) ;
65 double *matchmatrix;
66 matchmatrix = mxGetPr(prhs[arg]);
67 //mexPrintf("matchmatrix\n") ;
68 //for (int i=0; i< nr_entries; i++) {
69 // mexPrintf("entry %d: <%f>\n", i, matchmatrix[i]);
70 //}
71 arg++; // arg -> 5
72 #endif
73
74 /***************************************************************************/
75 // initialize alignment matrices and call fill_matrix()
76 /***************************************************************************/
77
78 //possible donor positions
79 int nr_donor_sites = 0 ;
80 for (int ii=0; ii<d_len; ii++) {
81 if(isfinite(donor[ii])) {
82 nr_donor_sites++ ;
83 }
84 }
85
86 int* donor_sites = new int[nr_donor_sites];
87 int donor_idx = 0 ;
88 for (int ii=0; ii<d_len; ii++) {
89 if(isfinite(donor[ii])) {
90 donor_sites[donor_idx] = ii+1 ;
91 donor_idx++ ;
92 }
93 }
94
95 int* max_score_positions = new int[nr_paths*2];
96
97 Pre_score** matrices = new Pre_score*[nr_paths];
98 for (int z=0; z<nr_paths; z++) {
99 matrices[z] = new Pre_score[dna_len * est_len];
100 }
101
102 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);
103
104 /***************************************************************************/
105 // return arguments etc.
106 /***************************************************************************/
107 int result_length; //Eine Variable fuer alle Matrizen
108
109 int* splice_align = new int[(dna_len-1)*nr_paths];
110 int* est_align = new int[(est_len-1)*nr_paths];
111 int* mmatrix_param = new int[(mlen*mlen)*nr_paths];
112 double alignmentscores[nr_paths]; //alignment score for each path/matrix
113
114 //memset((char*)splice_align, -1, (dna_len-1)*nr_paths*sizeof(int)); // fills splice_align with zeros
115 //memset((char*)est_align, -1, (est_len-1)*nr_paths*sizeof(int)); // fills est_align with zeros
116 //memset((char*)mmatrix_param, 0, mlen*mlen*nr_paths*sizeof(int)); //fills mmatrix_param with zeros
117 //memset(alignmentscores, -1, nr_paths*sizeof(double)); //fills alignmentscores with zeros
118
119 for (int z=0; z<nr_paths; z++) {
120 result_length = 0 ;
121
122 int* s_align = splice_align + (dna_len-1)*z; //pointer
123 int* e_align = est_align + (est_len-1)*z ; //pointer
124 int* mparam = mmatrix_param + (mlen*mlen)*z; //pointer
125
126 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);
127
128 // dnaest
129 double *DNA = new double[result_length];
130 double *EST = new double[result_length];
131
132 //backtracking
133 int i = max_score_positions[2*z] ; //i (est)
134 int j = max_score_positions[2*z +1] ; //j (dna)
135 int prev_i = ((Pre_score*)matrices[z] + i*dna_len +j)->prev_i;
136 int prev_j = ((Pre_score*)matrices[z] + i*dna_len +j)->prev_j;
137 int prev_m_no = ((Pre_score*)matrices[z] + i*dna_len +j)->prev_matrix_no;
138
139 for (int ii=result_length; ii>0; ii--) {
140 if ((prev_i == (i-1)) && (prev_j == (j-1))) { // match or mismatch
141 DNA[ii-1] = check_char(dna[j-1]) ;
142 EST[ii-1] = check_char(est[i-1]) ;
143 }
144 else if ((prev_i == (i)) && (prev_j == (j-1))) {// gap on EST
145 DNA[ii-1] = check_char(dna[j-1]) ;
146 EST[ii-1] = 0 ; //gap
147 }
148 else if ((prev_i == (i-1)) && (prev_j == (j))) { // gap on DNA
149 DNA[ii-1] = 0 ; //gap
150 EST[ii-1] = check_char(est[i-1]) ;
151 }
152 else {// splice site
153 for (j; j > prev_j; j--) {
154 DNA[ii-1] = check_char(dna[j-1]) ;
155 EST[ii-1] = 6 ; //intron
156 ii-- ;
157 }
158 ii++ ; // last ii-- too much (because done in for-loop)
159 }
160
161 i = prev_i;
162 j = prev_j;
163 prev_i = ((Pre_score*)matrices[prev_m_no] + i*dna_len + j)->prev_i;
164 prev_j = ((Pre_score*)matrices[prev_m_no] + i*dna_len + j)->prev_j;
165 prev_m_no = ((Pre_score*)matrices[prev_m_no] + i*dna_len + j)->prev_matrix_no;
166 }
167
168 } //end of z
169 }