binaries/src/tcoffee/t_coffee_source/util_dp_generic_fasta_nw.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <math.h>
   4 #include <stdarg.h>
   5 #include <string.h>
   6
   7
   8 #include "io_lib_header.h"
   9 #include "util_lib_header.h"
  10 #include "define_header.h"
  11 #include "dp_lib_header.h"
  12
  13
  14
  15 /*********************************************************************************/
  16 /*                                                                               */
  17 /*                                                                               */
  18 /*                       Generic DP                                              */
  19 /*                                                                               */
  20 /*                                                                               */
  21 /*********************************************************************************/
  22
  23
  24 Dp_Result * make_fast_generic_dp_pair_wise (Alignment *A, int*ns, int **l_s,Dp_Model *M)
  25         {
  26
  27           /*SIZE VARIABLES*/
  28
  29           int ndiag;
  30           int l0, l1, len_al,len_diag;
  31           static int max_len_al, max_len_diag;
  32           static int mI, mJ;
  33           /*Evaluation*/
  34           int **pos0;
  35
  36
  37
  38           /*DP VARIABLES*/
  39           static int *Mat, *LMat, *trace;
  40           int a, i, j,l;
  41           int state, cur_state, prev_state;
  42           int pos_i=0,  pos_j=0;
  43           int last_i=0, last_j=0;
  44           int prev_i, prev_j;
  45           int len_i, len_j, len;
  46           int t, e, em;
  47
  48           int prev_score;
  49           int pc, best_pc;
  50
  51           int *prev;
  52           int model_index;
  53           /*TRACEBACK*/
  54           Dp_Result *DPR;
  55           int k=0, next_k;
  56           int new_i, new_j;
  57
  58
  59           /*Cleqanning CALL*/
  60           if ( A==NULL)
  61             {
  62               max_len_al=0; max_len_diag=0;mI=0;mJ=0;
  63               vfree (Mat); vfree(LMat);vfree(trace);
  64               Mat=trace=LMat=NULL;
  65               return NULL;
  66             }
  67
  68           ndiag=M->diag[0];
  69
  70           l0=strlen (A->seq_al[l_s[0][0]]);
  71           l1=strlen (A->seq_al[l_s[1][0]]);
  72           len_al =l0+l1+1;
  73           len_diag=ndiag+4;
  74
  75
  76
  77           if ( (len_al>max_len_al || len_diag>max_len_diag))
  78             {
  79
  80               vfree (Mat);
  81               vfree (LMat);
  82               vfree(trace);
  83               max_len_diag=max_len_al=0;
  84             }
  85
  86           if (max_len_al==0)
  87             {
  88               max_len_al=len_al;
  89               max_len_diag=len_diag;
  90               mI=max_len_al*max_len_diag;
  91               mJ=max_len_diag;
  92
  93
  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));
  97
  98             }
  99
 100           prev=vcalloc ( M->nstate, sizeof (int));
 101           DPR=vcalloc ( 1, sizeof ( Dp_Result));
 102           DPR->traceback=vcalloc (max_len_al, sizeof (int));
 103
 104 /*PREPARE THE EVALUATION*/
 105
 106
 107           pos0=aln2pos_simple ( A,-1, ns, l_s);
 108
 109 /*INITIALIZATION OF THE DP MATRICES*/
 110
 111         for (i=0; i<=l0;i++)
 112           {
 113             for (j=0; j<=ndiag+1;j++)
 114               {
 115                 for ( state=0; state<M->nstate; state++)
 116                   {
 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;
 120                   }
 121               }
 122           }
 123
 124         M->diag[0]=1;
 125         M->diag[ndiag+1]=M->diag[ndiag];
 126
 127         for (i=0; i<=l0; i++)
 128           for ( j=0; j<=ndiag+1; j++)
 129             {
 130               pos_j=M->diag[j]-l0+i;
 131               pos_i=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)
 135                   {
 136                   for ( a=0; a< M->nstate; a++)
 137                     {
 138                      Mat  [a*mI+i*mJ+j]=0;
 139                      LMat [a*mI+i*mJ+j]=0;
 140                      trace[a*mI+i*mJ+j]=M->START;
 141                     }
 142                 }
 143               else
 144                 {
 145                   l=MAX(pos_i,pos_j);
 146                   for ( state=0; state<M->START; state++)
 147                     {
 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;
 150
 151
 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);*/
 155
 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;
 159                     }
 160                 }
 161             }
 162
 163 /*DYNAMIC PROGRAMMING: Forward Pass*/
 164
 165         /*Diagonals:
 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
 170         */
 171
 172         for (i=1; i<=l0;i++)
 173           {
 174             for (j=1; j<=ndiag;j++)
 175               {
 176                 pos_j=M->diag[j]-l0+i;
 177                 pos_i=i;
 178
 179                 if (pos_j<=0 || pos_j>l1 )continue;
 180                 last_i=i;
 181                 last_j=j;
 182
 183                 for (cur_state=0; cur_state<M->START; cur_state++)
 184                   {
 185                     if (M->model_properties[cur_state][M->DELTA_J])
 186                       {
 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]));
 189
 190                       }
 191                     else
 192                       {
 193                         prev_j=j;
 194                         prev_i=i+M->model_properties[cur_state][M->DELTA_I];
 195                       }
 196
 197
 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);
 201
 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);*/
 204
 205                     for (pc=best_pc=UNDEFINED, model_index=1; model_index<=M->bounded_model[cur_state][0]; model_index++)
 206                       {
 207                         prev_state=M->bounded_model[cur_state][model_index];
 208
 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];
 212                         e=em;
 213
 214                         if   (prev_score==UNDEFINED || len==UNDEFINED)e=UNDEFINED;
 215                         else if (len==0|| e==UNDEFINED)e=UNDEFINED;
 216                         else e=e*len;
 217
 218                         if (is_defined_int(3,prev_score,e, t))
 219                           {
 220                             pc=prev_score+t+e;
 221                           }
 222                         else  pc=UNDEFINED;
 223
 224                         /*Identify the best previous score*/
 225                         if (best_pc==UNDEFINED || (pc>best_pc && pc!=UNDEFINED))
 226                           {
 227                             prev[cur_state]=prev_state;
 228                             best_pc=pc;
 229
 230                           }
 231                       }
 232
 233                     Mat[cur_state*mI+i*mJ+j]=best_pc;
 234
 235
 236
 237                     if ( Mat[cur_state*mI+i*mJ+j]==UNDEFINED)
 238                       {
 239                         LMat[cur_state*mI+i*mJ+j]=UNDEFINED;
 240                         trace[cur_state*mI+i*mJ+j]=UNDEFINED;
 241                         continue;
 242                       }
 243
 244                     else if ( prev[cur_state]==cur_state)
 245                       {
 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];
 248                       }
 249                     else
 250                       {
 251                         LMat[cur_state*mI+i*mJ+j]=len;
 252                         trace[cur_state*mI+i*mJ+j]=prev[cur_state];
 253                       }
 254                   }
 255               }
 256           }
 257
 258
 259         i=last_i;
 260         j=last_j;
 261         for (pc=best_pc=UNDEFINED, state=0; state<M->START; state++)
 262           {
 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);
 265
 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);*/
 267
 268             l=LMat[state*mI+i*mJ+j];
 269
 270
 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];
 274
 275
 276             if (best_pc==UNDEFINED || (pc>best_pc && pc!=UNDEFINED))
 277               {
 278                 k=state;
 279                 best_pc=pc;
 280               }
 281           }
 282          DPR->score=best_pc;
 283
 284 /*TRACEBACK*/
 285
 286
 287         e=0;
 288         len=0;
 289
 290
 291         while (k!=M->START)
 292           {
 293             next_k=trace[k*mI+i*mJ+j];
 294
 295             new_i=i;
 296             new_j=j;
 297             l=LMat[k*mI+i*mJ+j];
 298             for (a=0; a< l; a++)
 299               {
 300                 DPR->traceback[len++]=k;
 301               }
 302            new_i+=M->model_properties[k][M->DELTA_I]*l;
 303
 304
 305            if ( M->model_properties[k][M->DELTA_J])
 306              {
 307                while ( next_k!=M->START && FABS((M->diag[j]-M->diag[new_j]))!=l)new_j+=M->model_properties[k][M->DELTA_J];
 308              }
 309
 310            i=new_i;
 311            j=new_j;
 312            k=next_k;
 313           }
 314         DPR->len=len;
 315         DPR->traceback[DPR->len++]=M->START;
 316         invert_list_int  (DPR->traceback,DPR->len);
 317         DPR->traceback[DPR->len]=M->END;
 318
 319         vfree (prev);
 320         free_int (pos0, -1);
 321         return DPR;
 322
 323
 324         }
 325
 326
 327 Constraint_list* free_dp_model  (Dp_Model *D)
 328        {
 329          Constraint_list *CL;
 330
 331          if ( !D)return NULL;
 332          CL=D->CL;
 333          vfree (D->diag);
 334          free_int (D->model, -1);
 335          free_int (D->model_properties, -1);
 336          free_int (D->bounded_model, -1);
 337
 338          vfree (D);
 339          return CL;
 340        }
 341
 342 Dp_Result * free_dp_result (Dp_Result *D )
 343         {
 344           if (!D) return NULL;
 345           vfree ( D->traceback);
 346           vfree (D);
 347           return NULL;
 348         }
 349
 350 /*********************************COPYRIGHT NOTICE**********************************/
 351 /*© Centro de Regulacio Genomica */
 352 /*and */
 353 /*Cedric Notredame */
 354 /*Tue Oct 27 10:12:26 WEST 2009. */
 355 /*All rights reserved.*/
 356 /*This file is part of T-COFFEE.*/
 357 /**/
 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.*/
 362 /**/
 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.*/
 367 /**/
 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 /*...............................................                                                                                                                                     |*/
 375 /**/
 376 /**/
 377 /*      */
 378 /*********************************COPYRIGHT NOTICE**********************************/