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