+ created necessary interfaces
[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
14 len = 0;
15 limits = 0;
16 penalties = 0;
17 max_len = 0;
18 min_len = 0;
19 cache = 0;
20 enum ETransformType transform = T_LINEAR;
21 id = 0;
22 name = 0;
23 use_svm = 0;
24 }
25
26 Alignment::~Alignment() {}
27
28 void Alignment::myalign(int nr_paths, char* dna, int dna_len, char* est,
29 int est_len, struct penalty_struct h, double* matchmatrix, int mm_len,
30 double* donor, int d_len, double* acceptor, int a_len,
31 bool remove_duplicate_scores, bool print_matrix) {}
32
33 void Alignment::getSpliceAlign(){}
34 void Alignment::getEstAlign(){}
35 void Alignment::getWeightMatch(){}
36 void Alignment::getTotalScores(){}
37 void Alignment::getDNAEST(){}
38
39 #if 0
40 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
41 {
42
43
44
45 /***************************************************************************/
46 // right number of arguments?
47 /***************************************************************************/
48 if(nrhs != 9)
49 {
50 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");
51 return ;
52 }
53
54 if(nlhs != 5)
55 {
56 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");
57 return ;
58 }
59
60
61 /***************************************************************************/
62 // variables
63 /***************************************************************************/
64 const int mlen = 6; // score matrix: length of 6 for "- A C G T N"
65 int arg = 0 ; // counter variable
66 int failure ; //
67
68
69 /***************************************************************************/
70 // 0. nr_paths
71 /***************************************************************************/
72 int nr_paths = (int)*mxGetPr(prhs[arg]); // 0. input (double)
73 arg++; // arg -> 1
74
75
76 /***************************************************************************/
77 // 1. dna
78 /***************************************************************************/
79 int dna_len = mxGetN(prhs[arg]) + 1; // mxGetN: number of columns (int)
80 char* dna = (char*)mxCalloc(dna_len, sizeof(char));
81
82 failure = mxGetString(prhs[arg], dna, dna_len); // NULL at end of string
83 if (failure) { // mxGetString returns 0 if success, 1 on failure
84 mexPrintf("couldn't open dna");
85 return;
86 }
87 arg++; // arg -> 2
88
89
90 /***************************************************************************/
91 // 2. est
92 /***************************************************************************/
93 int est_len = mxGetN(prhs[arg]) + 1;
94 char* est = (char*)mxCalloc(est_len, sizeof(char));
95
96 failure = mxGetString(prhs[arg], est, est_len);
97 if (failure) {
98 mexPrintf("couldn't open est");
99 return;
100 }
101 arg++; // arg -> 3
102
103
104 /***************************************************************************/
105 // 3. h (cell array)
106 /***************************************************************************/
107 penalty_struct* functions;
108 mwSize anz_func = 1;
109
110 functions = read_penalty_struct_from_cell(prhs[arg], anz_func);
111 //std::cout << " bla " << functions->max_len << " \n";
112 //std::cout << "lookup_penalty(5) = " << lookup_penalty(functions, 5, NULL, 0)<<std::endl ;
113 arg++; // arg -> 4
114
115
116 /***************************************************************************/
117 // 4. match matrix (columnwise)
118 /***************************************************************************/
119 int nr_entries = mxGetN(prhs[arg]) * mxGetM(prhs[arg]);
120 //mexPrintf("\nNumber of entries in matchmatrix: %d\n", nr_entries) ;
121 double *matchmatrix;
122 matchmatrix = mxGetPr(prhs[arg]);
123 //mexPrintf("matchmatrix\n") ;
124 //for (int i=0; i< nr_entries; i++) {
125 // mexPrintf("entry %d: <%f>\n", i, matchmatrix[i]);
126 //}
127 arg++; // arg -> 5
128
129
130 /***************************************************************************/
131 // 5. donor
132 /***************************************************************************/
133 int donor_length = mxGetN(prhs[arg]);
134 double *donor;
135 donor = mxGetPr(prhs[arg]);
136 //mexPrintf("Donor eingelesen: <%f>\n", donor[0]);
137 arg++; // arg -> 6
138
139
140 /***************************************************************************/
141 // 6. acceptor
142 /***************************************************************************/
143 int acceptor_length = mxGetN(prhs[arg]);
144 double *acceptor;
145 acceptor = mxGetPr(prhs[arg]);
146 //mexPrintf("Acceptor eingelesen: <%f>\n", acceptor[0]);
147 arg++; // arg -> 7
148
149
150 /***************************************************************************/
151 // 7.remove_duplicate_scores
152 /***************************************************************************/
153 bool remove_duplicate_scores = (bool)*mxGetPr(prhs[arg]);
154 arg++; // arg -> 8
155
156
157 /***************************************************************************/
158 // 8. print_matrix
159 /***************************************************************************/
160 int print_matrix = (int)*mxGetPr(prhs[arg]);
161 arg++; // arg -> 9
162
163
164 /***************************************************************************/
165 // initialize alignment matrices and call fill_matrix()
166 /***************************************************************************/
167
168 //possible donor positions
169 int nr_donor_sites = 0 ;
170 for (int ii=0; ii<donor_length; ii++) {
171 if(isfinite(donor[ii])) {
172 nr_donor_sites++ ;
173 }
174 }
175 //mwSize donor_sites[nr_donor_sites] ;
176 mwSize* donor_sites = (mwSize*)mxCalloc(nr_donor_sites,sizeof(mwSize));
177 int donor_idx = 0 ;
178 for (int ii=0; ii<donor_length; ii++) {
179 if(isfinite(donor[ii])) {
180 donor_sites[donor_idx] = ii+1 ;
181 donor_idx++ ;
182 }
183 }
184
185
186 //int max_score_positions[nr_paths*2]; //local alignment: where is the highest alignment score [i1,j1,i2,j2...]
187 mwSize* max_score_positions = (mwSize*)mxCalloc(nr_paths*2,sizeof(mwSize));
188
189 //Pre_score* matrices[nr_paths];
190 Pre_score** matrices = (Pre_score**)mxCalloc(nr_paths,sizeof(Pre_score*));
191 for (int z=0; z<nr_paths; z++) {
192 //matrices[z] = new Pre_score[dna_len * est_len];
193 //mexPrintf("size of one matrice: %d\n", est_len*dna_len*sizeof(Pre_score)) ;
194 matrices[z] = (Pre_score*)mxCalloc(dna_len*est_len,sizeof(Pre_score));
195 }
196
197 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);
198
199 //just some info
200 //for (int z=0; z<nr_paths; z++) {
201 // mexPrintf("(%d, %d)\n", max_score_positions[2*z], max_score_positions[2*z +1]) ;
202 //}
203
204 /***************************************************************************/
205 // return arguments etc.
206 /***************************************************************************/
207 int result_length; //Eine Variable fuer alle Matrizen
208
209 //int splice_align[(dna_len-1)*nr_paths]; //SpliceAlign for DNA - 1 line per matrix
210 //int est_align[(est_len-1)*nr_paths]; //SpliceAlign for EST - 1 line per matrix
211 //int mmatrix_param[(mlen*mlen)*nr_paths];//Matrix fuer Scorematrixparameter - 1 Zeile pro Matrix
212 int* splice_align = (int*)mxCalloc((dna_len-1)*nr_paths,sizeof(int));
213 int* est_align = (int*)mxCalloc((est_len-1)*nr_paths,sizeof(int));
214 int* mmatrix_param = (int*)mxCalloc((mlen*mlen)*nr_paths,sizeof(int));
215 double alignmentscores[nr_paths]; //alignment score for each path/matrix
216
217
218 memset((char*)splice_align, -1, (dna_len-1)*nr_paths*sizeof(int)); // fills splice_align with zeros
219 memset((char*)est_align, -1, (est_len-1)*nr_paths*sizeof(int)); // fills est_align with zeros
220 memset((char*)mmatrix_param, 0, mlen*mlen*nr_paths*sizeof(int)); //fills mmatrix_param with zeros
221 memset(alignmentscores, -1, nr_paths*sizeof(double)); //fills alignmentscores with zeros
222
223 // dnaest
224 plhs[nlhs-1] = mxCreateCellMatrix(2, nr_paths);
225
226
227
228 //mexPrintf("DEBUG: decoding %d paths\n",nr_paths);
229 for (int z=0; z<nr_paths; z++) {
230 result_length = 0 ;
231
232 int* s_align = splice_align + (dna_len-1)*z; //pointer
233 int* e_align = est_align + (est_len-1)*z ; //pointer
234 int* mparam = mmatrix_param + (mlen*mlen)*z; //pointer
235
236 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);
237
238 // dnaest
239 mxArray* mxdna = mxCreateDoubleMatrix(result_length, 1, mxREAL);
240 mxArray* mxest = mxCreateDoubleMatrix(result_length, 1, mxREAL);
241 double *DNA = mxGetPr(mxdna) ;
242 double *EST = mxGetPr(mxest) ;
243 mxSetCell(plhs[nlhs-1], z*2, mxdna) ;
244 mxSetCell(plhs[nlhs-1], z*2+1, mxest) ;
245
246 //backtracking
247 mwSize i = max_score_positions[2*z] ; //i (est)
248 mwSize j = max_score_positions[2*z +1] ; //j (dna)
249 mwSize prev_i = ((Pre_score*)matrices[z] + i*dna_len +j)->prev_i;
250 mwSize prev_j = ((Pre_score*)matrices[z] + i*dna_len +j)->prev_j;
251 mwSize prev_m_no = ((Pre_score*)matrices[z] + i*dna_len +j)->prev_matrix_no;
252
253 for (int ii=result_length; ii>0; ii--) {
254 if ((prev_i == (i-1)) && (prev_j == (j-1))) { // match or mismatch
255 DNA[ii-1] = check_char(dna[j-1]) ;
256 EST[ii-1] = check_char(est[i-1]) ;
257 }
258 else if ((prev_i == (i)) && (prev_j == (j-1))) {// gap on EST
259 DNA[ii-1] = check_char(dna[j-1]) ;
260 EST[ii-1] = 0 ; //gap
261 }
262 else if ((prev_i == (i-1)) && (prev_j == (j))) { // gap on DNA
263 DNA[ii-1] = 0 ; //gap
264 EST[ii-1] = check_char(est[i-1]) ;
265 }
266 else {// splice site
267 for (j; j > prev_j; j--) {
268 DNA[ii-1] = check_char(dna[j-1]) ;
269 EST[ii-1] = 6 ; //intron
270 ii-- ;
271 }
272 ii++ ; // last ii-- too much (because done in for-loop)
273 }
274
275 i = prev_i;
276 j = prev_j;
277 prev_i = ((Pre_score*)matrices[prev_m_no] + i*dna_len + j)->prev_i;
278 prev_j = ((Pre_score*)matrices[prev_m_no] + i*dna_len + j)->prev_j;
279 prev_m_no = ((Pre_score*)matrices[prev_m_no] + i*dna_len + j)->prev_matrix_no;
280 }
281
282 } //end of z
283
284
285
286 //mexPrintf("DEBUG:Returning arguments\n");
287 int index = 0 ;
288
289 //splice_align (mxArray):
290 const mwSize dims0[] = {(dna_len-1), nr_paths};
291 //const int dims0[] = {(dna_len-1), nr_paths};
292 plhs[index] = mxCreateNumericArray(2, dims0, mxINT32_CLASS, mxREAL);
293 mwSize* start0 = (mwSize*)(mxGetData(plhs[index]));
294 memcpy(start0, splice_align, mxGetM(plhs[index])*mxGetN(plhs[index])*mxGetElementSize(plhs[index]));
295 index++ ; // ->1
296
297 //est_align (mxArray):
298 const mwSize dims1[] = {(est_len-1), nr_paths};
299 //const int dims1[] = {(est_len-1), nr_paths};
300 plhs[index] = mxCreateNumericArray(2, dims1, mxINT32_CLASS, mxREAL);
301 mwSize* start1 = (mwSize*)(mxGetData(plhs[index]));
302 memcpy(start1, est_align, mxGetM(plhs[index])*mxGetN(plhs[index])*mxGetElementSize(plhs[index]));
303 index++ ; // ->2
304
305 //weightmatch (mxArray)
306 const mwSize dims2[] = {mlen*mlen, nr_paths};
307 //const int dims2[] = {mlen*mlen, nr_paths};
308 plhs[index] = mxCreateNumericArray(2, dims2, mxINT32_CLASS, mxREAL);
309 mwSize* start2= (mwSize*)(mxGetData(plhs[index]));
310 memcpy(start2, mmatrix_param, mxGetM(plhs[index])*mxGetN(plhs[index])*mxGetElementSize(plhs[index]));
311 index++ ; // ->3
312
313 //alignmentscore
314 plhs[index] = mxCreateDoubleMatrix(1,nr_paths,mxREAL);
315 if (plhs[index] == NULL)
316 {
317 printf("%s : Out of memory on line %d\n", __FILE__, __LINE__);
318 printf("Unable to create mxArray.\n");
319 exit(1);
320 }
321 double* start3 = (double*)(mxGetData(plhs[index]));
322 memcpy(start3, alignmentscores, mxGetM(plhs[index])*mxGetN(plhs[index])*mxGetElementSize(plhs[index]));
323 index++ ; // ->4
324
325 assert(nlhs==index+1) ;
326
327 //Free Memory
328 //mexPrintf("DEBUG:Free Memory\n");
329 mxFree(donor_sites);
330 mxFree(max_score_positions);
331 for (int i=nr_paths-1;i>=0; i--)
332 {
333 //delete[] matrices[i];
334 mxFree(matrices[i]);
335 }
336 mxFree(matrices);
337 mxFree(splice_align);
338 mxFree(est_align);
339 mxFree(mmatrix_param);
340 mxFree(dna);
341 mxFree(est);
342
343 return;
344
345 }
346 #endif