8 #include "io_lib_header.h"
9 #include "util_lib_header.h"
10 #include "define_header.h"
11 #include "dp_lib_header.h"
15 /*********************************************************************************/
21 /*********************************************************************************/
24 Dp_Result * make_fast_generic_dp_pair_wise (Alignment *A, int*ns, int **l_s,Dp_Model *M)
30 int l0, l1, len_al,len_diag;
31 static int max_len_al, max_len_diag;
39 static int *Mat, *LMat, *trace;
41 int state, cur_state, prev_state;
43 int last_i=0, last_j=0;
45 int len_i, len_j, len;
62 max_len_al=0; max_len_diag=0;mI=0;mJ=0;
63 vfree (Mat); vfree(LMat);vfree(trace);
70 l0=strlen (A->seq_al[l_s[0][0]]);
71 l1=strlen (A->seq_al[l_s[1][0]]);
77 if ( (len_al>max_len_al || len_diag>max_len_diag))
83 max_len_diag=max_len_al=0;
89 max_len_diag=len_diag;
90 mI=max_len_al*max_len_diag;
94 Mat =vcalloc ( M->nstate*max_len_al*max_len_diag, sizeof (int));
95 LMat =vcalloc ( M->nstate*max_len_al*max_len_diag, sizeof (int));
96 trace=vcalloc ( M->nstate*max_len_al*max_len_diag, sizeof (int));
100 prev=vcalloc ( M->nstate, sizeof (int));
101 DPR=vcalloc ( 1, sizeof ( Dp_Result));
102 DPR->traceback=vcalloc (max_len_al, sizeof (int));
104 /*PREPARE THE EVALUATION*/
107 pos0=aln2pos_simple ( A,-1, ns, l_s);
109 /*INITIALIZATION OF THE DP MATRICES*/
113 for (j=0; j<=ndiag+1;j++)
115 for ( state=0; state<M->nstate; state++)
117 Mat [state*mI+i*mJ+j]=UNDEFINED;
118 LMat [state*mI+i*mJ+j]=UNDEFINED;
119 trace [state*mI+i*mJ+j]=M->START;
125 M->diag[ndiag+1]=M->diag[ndiag];
127 for (i=0; i<=l0; i++)
128 for ( j=0; j<=ndiag+1; j++)
130 pos_j=M->diag[j]-l0+i;
132 if (!(pos_j==0 || pos_i==0))continue;
133 if ( pos_j<0 || pos_i<0)continue;
134 if ( pos_i==0 && pos_j==0)
136 for ( a=0; a< M->nstate; a++)
139 LMat [a*mI+i*mJ+j]=0;
140 trace[a*mI+i*mJ+j]=M->START;
146 for ( state=0; state<M->START; state++)
148 if (pos_j==0 && M->model_properties[state][M->LEN_J])continue;
149 if (pos_i==0 && M->model_properties[state][M->LEN_I])continue;
152 t=M->model[M->START][state];
153 e=((M->model_emission_function)[state][M->START_EMISSION])(A, pos0, ns[0], l_s[0], pos_i-1, pos0, ns[1], l_s[1],pos_j-1,M->CL);
154 /*e=((M->get_dp_cost_list)[M->model_properties[state][M->START_EMISSION]])(A, pos0, ns[0], l_s[0], pos_i-1, pos0, ns[1], l_s[1],pos_j-1,M->CL);*/
156 Mat [state*mI+i*mJ+j]=t+e*l;
157 LMat [state*mI+i*mJ+j]=l;
158 trace [state*mI+i*mJ+j]=M->START;
163 /*DYNAMIC PROGRAMMING: Forward Pass*/
166 M->diag[0]=Number of diagonals being considered
167 M->diag[1]=First diagonal being considered
168 Diagonals are numbered 1...L0+l1-1
169 1 is the bottom-left diag
174 for (j=1; j<=ndiag;j++)
176 pos_j=M->diag[j]-l0+i;
179 if (pos_j<=0 || pos_j>l1 )continue;
183 for (cur_state=0; cur_state<M->START; cur_state++)
185 if (M->model_properties[cur_state][M->DELTA_J])
187 prev_j=j+M->model_properties[cur_state][M->DELTA_J];
188 prev_i=i+M->model_properties[cur_state][M->DELTA_I]*FABS((M->diag[j]-M->diag[prev_j]));
194 prev_i=i+M->model_properties[cur_state][M->DELTA_I];
198 len_i=FABS((i-prev_i));
199 len_j=FABS((M->diag[prev_j]-M->diag[j]));
200 len=MAX(len_i, len_j);
202 em=((M->model_emission_function[cur_state][M->EMISSION]))(A, pos0, ns[0], l_s[0], pos_i-1, pos0, ns[1], l_s[1],pos_j-1,M->CL);
203 /*em=((M->get_dp_cost_list)[M->model_properties[cur_state][M->EMISSION]])(A, pos0, ns[0], l_s[0], pos_i-1, pos0, ns[1], l_s[1],pos_j-1,M->CL);*/
205 for (pc=best_pc=UNDEFINED, model_index=1; model_index<=M->bounded_model[cur_state][0]; model_index++)
207 prev_state=M->bounded_model[cur_state][model_index];
209 if(prev_i<0 || prev_j<0 ||prev_i>l0 || prev_j>ndiag || len==UNDEFINED)prev_score=UNDEFINED;
210 else prev_score=Mat[prev_state*mI+prev_i*mJ+prev_j];
211 t=M->model[prev_state][cur_state];
214 if (prev_score==UNDEFINED || len==UNDEFINED)e=UNDEFINED;
215 else if (len==0|| e==UNDEFINED)e=UNDEFINED;
218 if (is_defined_int(3,prev_score,e, t))
224 /*Identify the best previous score*/
225 if (best_pc==UNDEFINED || (pc>best_pc && pc!=UNDEFINED))
227 prev[cur_state]=prev_state;
233 Mat[cur_state*mI+i*mJ+j]=best_pc;
237 if ( Mat[cur_state*mI+i*mJ+j]==UNDEFINED)
239 LMat[cur_state*mI+i*mJ+j]=UNDEFINED;
240 trace[cur_state*mI+i*mJ+j]=UNDEFINED;
244 else if ( prev[cur_state]==cur_state)
246 LMat [cur_state*mI+i*mJ+j]= LMat [cur_state*mI+prev_i*mJ+prev_j]+len;
247 trace[cur_state*mI+i*mJ+j]= trace[cur_state*mI+prev_i*mJ+prev_j];
251 LMat[cur_state*mI+i*mJ+j]=len;
252 trace[cur_state*mI+i*mJ+j]=prev[cur_state];
261 for (pc=best_pc=UNDEFINED, state=0; state<M->START; state++)
263 t=M->model[state][M->END];
264 e=( M->model_emission_function[state][M->TERM_EMISSION])(A, pos0, ns[0], l_s[0], pos_i-1, pos0, ns[1], l_s[1],pos_j-1,M->CL);
266 /*e=((M->get_dp_cost_list)[M->model_properties[state][M->TERM_EMISSION]])(A, pos0, ns[0], l_s[0], pos_i-1, pos0, ns[1], l_s[1],pos_j-1,M->CL);*/
268 l=LMat[state*mI+i*mJ+j];
271 if (!is_defined_int(4,t,e,Mat[state*mI+i*mJ+j],l))Mat[state*mI+i*mJ+j]=UNDEFINED;
272 else Mat[state*mI+i*mJ+j]+=t+e*(l);
273 pc=Mat[state*mI+i*mJ+j];
276 if (best_pc==UNDEFINED || (pc>best_pc && pc!=UNDEFINED))
293 next_k=trace[k*mI+i*mJ+j];
300 DPR->traceback[len++]=k;
302 new_i+=M->model_properties[k][M->DELTA_I]*l;
305 if ( M->model_properties[k][M->DELTA_J])
307 while ( next_k!=M->START && FABS((M->diag[j]-M->diag[new_j]))!=l)new_j+=M->model_properties[k][M->DELTA_J];
315 DPR->traceback[DPR->len++]=M->START;
316 invert_list_int (DPR->traceback,DPR->len);
317 DPR->traceback[DPR->len]=M->END;
327 Constraint_list* free_dp_model (Dp_Model *D)
334 free_int (D->model, -1);
335 free_int (D->model_properties, -1);
336 free_int (D->bounded_model, -1);
342 Dp_Result * free_dp_result (Dp_Result *D )
345 vfree ( D->traceback);
350 /*********************************COPYRIGHT NOTICE**********************************/
351 /*© Centro de Regulacio Genomica */
353 /*Cedric Notredame */
354 /*Tue Oct 27 10:12:26 WEST 2009. */
355 /*All rights reserved.*/
356 /*This file is part of T-COFFEE.*/
358 /* T-COFFEE is free software; you can redistribute it and/or modify*/
359 /* it under the terms of the GNU General Public License as published by*/
360 /* the Free Software Foundation; either version 2 of the License, or*/
361 /* (at your option) any later version.*/
363 /* T-COFFEE is distributed in the hope that it will be useful,*/
364 /* but WITHOUT ANY WARRANTY; without even the implied warranty of*/
365 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the*/
366 /* GNU General Public License for more details.*/
368 /* You should have received a copy of the GNU General Public License*/
369 /* along with Foobar; if not, write to the Free Software*/
370 /* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA*/
371 /*............................................... |*/
372 /* If you need some more information*/
373 /* cedric.notredame@europe.com*/
374 /*............................................... |*/
378 /*********************************COPYRIGHT NOTICE**********************************/