binaries/src/tcoffee/t_coffee_source/fsa_dp.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <math.h>
   4 #include <stdarg.h>
   5 #include <string.h>
   6
   7 #include "io_lib_header.h"
   8 #include "util_lib_header.h"
   9 #include "define_header.h"
  10 #include "dp_lib_header.h"
  11
  12 #define hmm_add(x,y) ((x==UNDEFINED || y==UNDEFINED)?UNDEFINED:(x+y))
  13 #define MAX_EMISSION 256
  14
  15 /*********************************************************************************/
  16 /*                                                                               */
  17 /*                                                                               */
  18 /*                       Procons dp                                              */
  19 /*                                                                               */
  20 /*                                                                               */
  21 /*********************************************************************************/
  22 char alphabetDefault[] = "ARNDCQEGHILKMFPSTWYV";
  23 double emitPairsDefault[20][20] = {
  24   {0.02373072f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
  25   {0.00244502f, 0.01775118f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
  26   {0.00210228f, 0.00207782f, 0.01281864f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
  27   {0.00223549f, 0.00161657f, 0.00353540f, 0.01911178f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
  28   {0.00145515f, 0.00044701f, 0.00042479f, 0.00036798f, 0.01013470f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
  29   {0.00219102f, 0.00253532f, 0.00158223f, 0.00176784f, 0.00032102f, 0.00756604f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
  30   {0.00332218f, 0.00268865f, 0.00224738f, 0.00496800f, 0.00037956f, 0.00345128f, 0.01676565f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
  31   {0.00597898f, 0.00194865f, 0.00288882f, 0.00235249f, 0.00071206f, 0.00142432f, 0.00214860f, 0.04062876f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
  32   {0.00114353f, 0.00132105f, 0.00141205f, 0.00097077f, 0.00026421f, 0.00113901f, 0.00131767f, 0.00103704f, 0.00867996f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
  33   {0.00318853f, 0.00138145f, 0.00104273f, 0.00105355f, 0.00094040f, 0.00100883f, 0.00124207f, 0.00142520f, 0.00059716f, 0.01778263f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
  34   {0.00449576f, 0.00246811f, 0.00160275f, 0.00161966f, 0.00138494f, 0.00180553f, 0.00222063f, 0.00212853f, 0.00111754f, 0.01071834f, 0.03583921f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
  35   {0.00331693f, 0.00595650f, 0.00257310f, 0.00252518f, 0.00046951f, 0.00312308f, 0.00428420f, 0.00259311f, 0.00121376f, 0.00157852f, 0.00259626f, 0.01612228f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
  36   {0.00148878f, 0.00076734f, 0.00063401f, 0.00047808f, 0.00037421f, 0.00075546f, 0.00076105f, 0.00066504f, 0.00042237f, 0.00224097f, 0.00461939f, 0.00096120f, 0.00409522f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
  37   {0.00165004f, 0.00090768f, 0.00084658f, 0.00069041f, 0.00052274f, 0.00059248f, 0.00078814f, 0.00115204f, 0.00072545f, 0.00279948f, 0.00533369f, 0.00087222f, 0.00116111f, 0.01661038f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
  38   {0.00230618f, 0.00106268f, 0.00100282f, 0.00125381f, 0.00034766f, 0.00090111f, 0.00151550f, 0.00155601f, 0.00049078f, 0.00103767f, 0.00157310f, 0.00154836f, 0.00046718f, 0.00060701f, 0.01846071f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
  39   {0.00631752f, 0.00224540f, 0.00301397f, 0.00285226f, 0.00094867f, 0.00191155f, 0.00293898f, 0.00381962f, 0.00116422f, 0.00173565f, 0.00250962f, 0.00312633f, 0.00087787f, 0.00119036f, 0.00180037f, 0.01346609f, 0.0f, 0.0f, 0.0f, 0.0f},
  40   {0.00389995f, 0.00186053f, 0.00220144f, 0.00180488f, 0.00073798f, 0.00154526f, 0.00216760f, 0.00214841f, 0.00077747f, 0.00248968f, 0.00302273f, 0.00250862f, 0.00093371f, 0.00107595f, 0.00147982f, 0.00487295f, 0.01299436f, 0.0f, 0.0f, 0.0f},
  41   {0.00039119f, 0.00029139f, 0.00021006f, 0.00016015f, 0.00010666f, 0.00020592f, 0.00023815f, 0.00038786f, 0.00019097f, 0.00039549f, 0.00076736f, 0.00028448f, 0.00016253f, 0.00085751f, 0.00015674f, 0.00026525f, 0.00024961f, 0.00563625f, 0.0f, 0.0f},
  42   {0.00131840f, 0.00099430f, 0.00074960f, 0.00066005f, 0.00036626f, 0.00070192f, 0.00092548f, 0.00089301f, 0.00131038f, 0.00127857f, 0.00219713f, 0.00100817f, 0.00054105f, 0.00368739f, 0.00047608f, 0.00102648f, 0.00094759f, 0.00069226f, 0.00999315f, 0.0f},
  43   {0.00533241f, 0.00169359f, 0.00136609f, 0.00127915f, 0.00119152f, 0.00132844f, 0.00178697f, 0.00194579f, 0.00071553f, 0.01117956f, 0.00914460f, 0.00210897f, 0.00197461f, 0.00256159f, 0.00135781f, 0.00241601f, 0.00343452f, 0.00038538f, 0.00148001f, 0.02075171f}
  44 };
  45
  46 static void DisplayMatState ( MatState *S, char *s);
  47
  48
  49
  50
  51
  52
  53 void check_viterbiL ( Alignment *A,int *ns, int **ls, Constraint_list *CL);
  54 MatState *viterbi2path2 ( double ***Sc, int ***St, Hmm *H, MatState *S, MatState *E);
  55 void testfunc ( MatState *S, char *s);
  56
  57 #ifdef IN_PGROGRESS
  58 /*********************************************************************************/
  59 /*                                                                               */
  60 /*                                                                               */
  61 /*                     MSA Analyzer                                              */
  62 /*                                                                               */
  63 /*                                                                               */
  64 /*********************************************************************************/
  65 Alignment * analyze_alignment ( Alignment *A)
  66 {
  67   evaluate_alignment (A);
  68   H=define_msa_model (-100);
  69   M=seq_viterbi_hmm (A->seq_al[0], H);
  70   path=seq_viterbi2path ( seq, H, M);
  71 }
  72
  73
  74 Hmm* define_msa_model(double penalty)
  75 {
  76   Hmm *H;
  77   double freeT=0;
  78   int n=0;
  79   HmmState *S;
  80
  81
  82   H=declare_hmm(2);
  83   H->freeT=freeT=0;
  84
  85   H->forbiden=FORBIDEN;
  86   H->start=START_STATE;
  87   H->end=END_STATE;
  88
  89   /*define START*/
  90   S=H->S[n];
  91   sprintf (S->name, "START"); S->state=n;
  92
  93   S->DI=0;
  94   S->DJ=0;
  95   S->em=freeT;
  96
  97   sprintf ( (S->T[S->nT])->name, "C") ;(S->T[S->nT])->tr=freeT  ;S->nT++;
  98   sprintf ( (S->T[S->nT])->name, "W");(S->T[S->nT])->tr=freeT  ;S->nT++;
  99   n++;
 100   /*define END*/
 101   S=H->S[n];
 102   sprintf (S->name, "END"); S->state=n;
 103   S->DI=0;
 104   S->DJ=0;
 105   S->em=freeT;
 106   n++;
 107
 108   /*define Correct*/
 109   S=H->S[n];
 110   sprintf (S->name, "C"); S->state=n;
 111   S->DI=1;
 112   S->DJ=0;
 113   S->em=H->forbiden;
 114   S->em_func=em_correct_msa;
 115
 116   sprintf ( (S->T[S->nT])->name, "C") ;(S->T[S->nT])->tr=freeT;S->nT++;
 117   sprintf ( (S->T[S->nT])->name, "W");(S->T[S->nT])->tr=penalty  ;S->nT++;
 118   sprintf ( (S->T[S->nT])->name, "END");   (S->T[S->nT])->tr=freeT;S->nT++;
 119   n++;
 120
 121   /*define Wrong*/
 122   S=H->S[n];
 123   sprintf (S->name, "INSERT"); S->state=n;
 124   S->DI=1;
 125   S->DJ=0;
 126   S->em=H->forbiden;
 127   S->em_func=em_wrong_msa;
 128   sprintf ( (S->T[S->nT])->name, "C") ; (S->T[S->nT])->tr=penalty;S->nT++;
 129   sprintf ( (S->T[S->nT])->name, "W"); (S->T[S->nT])->tr=freeT;S->nT++;
 130   sprintf ( (S->T[S->nT])->name, "END");    (S->T[S->nT])->tr=-gop;S->nT++;
 131   n++;
 132
 133   /*define LInsert*/
 134   S=H->S[n];
 135   sprintf (S->name, "LINSERT"); S->state=n;
 136   S->DI=1;
 137   S->DJ=0;
 138   S->em=lgep;
 139
 140   sprintf ( (S->T[S->nT])->name, "INSERT")  ;(S->T[S->nT])->tr=freeT;S->nT++;
 141   sprintf ( (S->T[S->nT])->name, "LINSERT");(S->T[S->nT])->tr=freeT;S->nT++;
 142   n++;
 143
 144   H=bound_hmm ( H);
 145   return H;
 146 }
 147 #endif
 148 /*********************************************************************************/
 149 /*                                                                               */
 150 /*                                                                               */
 151 /*                     simple HMM: Viterbi                                       */
 152 /*                                                                               */
 153 /*                                                                               */
 154 /*********************************************************************************/
 155 double pavie_em_func (Hmm*H, HmmState *S, int v);
 156 Hmm* define_full_model(int nstate,char **state_list, char *model_name,Generic_em_func evaluation_func );
 157 char **produce_state_name (int nstate,char **list, char *model_name, Hmm* H);
 158 double** seq_viterbi_hmm (char *seq, Hmm *H);
 159 int * seq_viterbi2path (char *s, Hmm *H, double **M);
 160 double analyze_sequence ( char *seq, Hmm*H);
 161
 162 double pavie_emission (Hmm*H, HmmState *S, int v)
 163 {
 164   char *n;
 165
 166
 167   n=S->name;
 168
 169   if ( v==n[0] || ( v=='*' && n[0]=='E')) return H->freeT;
 170   return H->forbiden;
 171 }
 172 Hmm* define_full_model(int nstate, char **list, char *model_name, Generic_em_func emission_function)
 173 {
 174   /*list: a list of the state names: does not include START or END*/
 175   /*model_name: a string that will be appended to the names in list*/
 176
 177   Hmm *H;
 178   int a,n;
 179   HmmState *S;
 180
 181
 182   H=declare_hmm(nstate+2);
 183   H->freeT=0;
 184
 185   H->forbiden=FORBIDEN;
 186   H->start=START_STATE;
 187   H->end=END_STATE;
 188
 189   list=produce_state_name (nstate,list, model_name, H);
 190   nstate+=2;
 191
 192   for (n=0; n<nstate; n++)
 193     {
 194
 195       S=H->S[n];
 196       S->state=n;
 197       sprintf ( S->name, "%s", list[n]);
 198       S->em_func2=emission_function;
 199       if (n==H->end || n==H->start){S->DI=0;S->DJ=0;}
 200       else S->DI=1;S->DJ=0;
 201
 202       /*Emmissions*/
 203       S->em_func2=emission_function;
 204       for (a=0; a< MAX_EMISSION; a++)S->em2[a]=H->freeT;
 205
 206       for (a=0; a<nstate && n!=H->end; a++)
 207         {
 208           if (a!=H->start && !(n==H->start && a==H->end) )
 209           {
 210             sprintf ( (S->T[S->nT])->name, "%s", list[a]);
 211             (S->T[S->nT])->tr=H->freeT;
 212             S->nT++;
 213           }
 214         }
 215     }
 216   return H;
 217 }
 218
 219 char **produce_state_name (int nstate,char **list, char *model_name, Hmm* H)
 220 {
 221   int a,b,c;
 222   char **new_list;
 223   nstate+=2;
 224
 225   new_list=declare_char ( nstate, 100);
 226   for ( a=0, b=0, c=0; a< nstate; a++)
 227     {
 228       if ( a==H->start)sprintf ( new_list[a], "START");
 229       else if ( a==H->end)sprintf ( new_list[a], "END");
 230       else if ( list==NULL){sprintf ( new_list[a], "%c%s", 'a'+b, (model_name)?model_name:"");b++;}
 231       else {sprintf ( new_list[a], "%c%s", list[a][c], (model_name)?model_name:"");c++;}
 232     }
 233   return new_list;
 234 }
 235
 236 int seq_viterbi_pair_wise (Alignment *A,int*ns, int **ls,Constraint_list *CL)
 237 {
 238   ungap(A->seq_al[0]);
 239   analyze_sequence (A->seq_al[0], NULL);
 240   myexit (EXIT_FAILURE);
 241   return 1;
 242 }
 243 double analyze_sequence ( char *seq, Hmm *H)
 244 {
 245
 246   double **M;
 247   int *path;
 248
 249   if ( H==NULL)
 250     {
 251       H=define_full_model(5, NULL,"_first", pavie_emission);
 252       H=bound_hmm(H);
 253       DisplayHmm (H);
 254     }
 255   M=seq_viterbi_hmm (seq, H);
 256   path=seq_viterbi2path (seq, H, M);
 257   return M[H->end][strlen (seq)];
 258 }
 259
 260
 261 double** seq_viterbi_hmm (char *seq, Hmm *H)
 262 {
 263   /*Given a model H and a sequence seq*/
 264   double **M;
 265   double e, v, max;
 266   int i,pi, bestk, s, k, l1;
 267   HmmState *S1, *S2;
 268
 269
 270   l1=strlen (seq);
 271   M=declare_double (H->nS*2,l1+2);
 272
 273   /*Handle the start*/
 274   M[H->start][0]=0;
 275   for ( i=0; i<=l1; i++)
 276     {
 277       for ( s=0; s< H->nS; s++)
 278         {
 279           S1=H->S[s];
 280           pi=i-S1->DI;
 281           max=H->forbiden;
 282           bestk=H->forbiden;
 283           if ( pi<0){M[s][i]=H->forbiden;}/*Boundary*/
 284           else
 285             {
 286               if (pi==0) {max=H->T[(int)H->start][s];bestk=H->start;}/*Start*/
 287               else
 288                 {
 289                   for (k=1; k<=H->fromM[S1->state][0]; k++)
 290                     {
 291                       S2=H->S[H->fromM[s][k]];
 292                       if ( S2->state==H->start || S2->state==H->end)continue;
 293                       v=hmm_add((M[S2->state][pi]),(H->T[S2->state][S1->state]));
 294                       if ( v!=H->forbiden && (max==H->forbiden || v>max)){max=v;bestk=S2->state;}
 295                     }
 296                 }
 297               if (S1->em2)e=S1->em2[(int)seq[pi]];
 298               else e=S1->em_func2(H,S1, (int)seq[pi]);
 299
 300               e=hmm_add (e,max);
 301
 302               M[s][i]=e;
 303               M[s+H->nS][i]=bestk;
 304               }
 305           }
 306       }
 307   /*Terminate viterbi: connect the path to the END state*/
 308   max=UNDEFINED;
 309   bestk=UNDEFINED;
 310   for (k=0; k< H->nS; k++)
 311     {
 312       if (k==H->start || k==H->end);
 313       else
 314         {
 315           v=(M[k][l1]==H->forbiden || H->T[k][H->end]==H->forbiden)?H->forbiden:M[k][l1]+H->T[k][H->end];
 316           if ( max==H->forbiden || v>max){bestk=k;max=v;}
 317         }
 318     }
 319   M[H->end][l1]=max;
 320   M[H->nS+H->end][l1]=bestk;
 321   return M;
 322 }
 323
 324 int * seq_viterbi2path (char *s, Hmm *H, double **M)
 325 {
 326   int i,l,l1;
 327   int *path;
 328   HmmState *S1;
 329   int cs;
 330
 331   l1=strlen (s);
 332   path=vcalloc (l1+1, sizeof (int));
 333   i=l1;
 334   l=0;
 335   cs=M[H->nS+H->end][i];
 336
 337   while (i>0)
 338     {
 339
 340       S1=H->S[cs];
 341       path[l++]=cs;
 342
 343       cs=M[H->nS+cs][i];
 344       i-=S1->DI;
 345       /*fprintf ( stderr, "%d", cs);*/
 346     }
 347   invert_list_int (path, l);
 348   path[l++]=H->forbiden;
 349
 350   return path;
 351 }
 352
 353 /*********************************************************************************/
 354 /*                                                                               */
 355 /*                                                                               */
 356 /*                     pairHMM: Viterbi                                              */
 357 /*                                                                               */
 358 /*                                                                               */
 359 /*********************************************************************************/
 360 Hmm* define_mnm_model(Constraint_list *CL);
 361 int viterbi_pair_wise_OLD (Alignment *A,int*ns, int **ls,Constraint_list *CL)
 362 {
 363   int l1, l2, a;
 364   double ***M;
 365   int *path;
 366   Hmm * H;
 367
 368   A->pos=aln2pos_simple( A, -1, ns, ls);
 369
 370   //  H=define_mnm_model (CL);
 371   H=define_two_mat_model (CL);
 372
 373   l1=strlen (A->seq_al[ls[0][0]]);
 374   l2=strlen (A->seq_al[ls[1][0]]);
 375   M=viterbi_hmm (A, ns, ls, H, CL);
 376   path=viterbi2path (l1,l2, H,M);
 377   A=viterbipath2aln (A,ns,ls,path, H);
 378   A->score=A->score_aln=M[H->end][l1][l2];
 379   for ( a=0; a< H->nS*2; a++)free_double (M[a], -1);
 380   vfree (M);
 381   free_int (A->pos, -1);
 382   A->pos=NULL;
 383
 384   free_Hmm (H);
 385   vfree (path);
 386
 387   return A->score_aln;
 388 }
 389
 390 Alignment * viterbipath2aln (Alignment *A, int *ns,int **ls,int *tb, Hmm *H)
 391 {
 392   char **aln;
 393   char *char_buf;
 394   int a, b, c, len, ch;
 395   HmmState *S;
 396   int l[2];
 397
 398   len=0;while (tb[len]!=H->forbiden)len++;
 399
 400   if ( A->declared_len<=len)A=realloc_aln2  ( A,A->max_n_seq,2*len);
 401   aln=A->seq_al;
 402
 403   char_buf=vcalloc (len+1, sizeof (char));
 404   l[0]=strlen ( A->seq_al[ls[0][0]]);
 405   l[1]=strlen ( A->seq_al[ls[1][0]]);
 406
 407   for ( c=0; c< 2; c++)
 408     for ( a=0; a< ns[c]; a++)
 409       {
 410         for (ch=0, b=0; b<len; b++)
 411           {
 412             S=H->S[tb[b]];
 413             if ( (c==0 && S->DI)|| (c==1 && S->DJ) )
 414               char_buf[b]=aln[ls[c][a]][ch++];
 415             else
 416               char_buf[b]='-';
 417           }
 418         char_buf[b]='\0';
 419         sprintf (aln[ls[c][a]],"%s", char_buf);
 420         if ( l[c]!=ch){fprintf (stderr, "\nERROR: Wrong Size Of Alignmnent (Real %d, Observed %d)[FATAL:%s]",l[c], ch, PROGRAM);}
 421       }
 422   A->len_aln=len;
 423   A->nseq=ns[0]+ns[1];
 424
 425   vfree(char_buf);
 426   return A;
 427 }
 428
 429 double*** viterbi_hmm (Alignment *A,int *ns, int **ls, Hmm *H, Constraint_list *CL)
 430 {
 431   double ***M;
 432   double e, v, max;
 433   int a, i,pi, bestk,j,pj, s, k, l1, l2;
 434   HmmState *S1, *S2;
 435
 436   l1=strlen (A->seq_al[ls[0][0]]);
 437   l2=strlen (A->seq_al[ls[1][0]]);
 438
 439   M=vcalloc (H->nS*2, sizeof (double**));
 440   for ( a=0; a<H->nS*2; a++)M[a]=declare_double (l1+2, l2+2);
 441
 442   /*Handle the start*/
 443
 444   M[H->start][0][0]=0;
 445   for ( i=0; i<=l1; i++)
 446     for ( j=0; j<=l2; j++)
 447       {
 448         for ( s=0; s< H->nS; s++)
 449           {
 450             S1=H->S[s];
 451             pi=i-S1->DI;
 452             pj=j-S1->DJ;
 453             max=H->forbiden;
 454             bestk=H->forbiden;
 455             if ( pi<0 ||pj<0){M[s][i][j]=H->forbiden;}/*Boundary*/
 456             else
 457               {
 458                 if (pi+pj==0) {max=H->T[H->start][s];bestk=H->start;}/*Start*/
 459                 else
 460                   {
 461                     for (k=1; k<=H->fromM[S1->state][0]; k++)
 462                       {
 463                         S2=H->S[H->fromM[s][k]];
 464                         if ( S2->state==H->start || S2->state==H->end)continue;
 465                         v=(M[S2->state][pi][pj]==H->forbiden)?H->forbiden:(M[S2->state][pi][pj]+H->T[S2->state][S1->state]);
 466                         if ( v!=H->forbiden && (max==H->forbiden || v>max)){max=v;bestk=S2->state;}
 467                       }
 468                   }
 469
 470                 e=(S1->em==H->forbiden)?S1->em_func (A, A->pos, ns[0], ls[0],i-1, A->pos,ns[1], ls[1], j-1, CL):S1->em;
 471                 e=(max==H->forbiden || e==H->forbiden)?H->forbiden:e+max;
 472
 473                 M[s][i][j]=e;
 474                 M[s+H->nS][i][j]=bestk;
 475               }
 476           }
 477       }
 478
 479   /*Terminate viterbi: connect the path to the END state*/
 480   max=UNDEFINED;
 481   bestk=UNDEFINED;
 482   for (k=0; k< H->nS; k++)
 483     {
 484       if (k==H->start || k==H->end);
 485       else
 486         {
 487           v=(M[k][l1][l2]==H->forbiden || H->T[k][H->end]==H->forbiden)?H->forbiden:M[k][l1][l2]+H->T[k][H->end];
 488           if ( max==H->forbiden || v>max){bestk=k;max=v;}
 489         }
 490     }
 491   M[H->end][l1][l2]=max;
 492   M[H->nS+H->end][l1][l2]=bestk;
 493
 494   return M;
 495 }
 496
 497 /*********************************************************************************/
 498 /*                                                                               */
 499 /*                                                                               */
 500 /*                      HMM: Decode/Traceback                                       */
 501 /*                                                                               */
 502 /*                                                                               */
 503 /*********************************************************************************/
 504 int * traceback (Alignment *A,int *ns, int **ls, Hmm *H, Constraint_list *CL,MatState *S, MatState *E, int **seg_list)
 505
 506 {
 507   int *path;
 508   int l=0;
 509   MatState *N;
 510   int l1, l2;
 511
 512   l1=strlen (A->seq_al[ls[0][0]]);
 513   l2=strlen (A->seq_al[ls[1][0]]);
 514   path=vcalloc ( l1+l2+1, sizeof (int));
 515
 516   while ( S->st!=H->end)
 517     {
 518       DisplayMatState (S, "\n\tTraceback");
 519       N=S->n;
 520       if ( N && S && (((N->i-S->i)>1) ||((N->j-S->j)>1)))
 521         {
 522           RviterbiD_hmm (A,ns,ls,H,CL,S,N,seg_list);
 523           N=S->n;
 524         }
 525
 526       path[l++]=S->st;
 527       ManageMatState (FREE,S);
 528       S=N;
 529     }
 530
 531   path[l]=H->forbiden;
 532   return path;
 533 }
 534
 535 int * viterbi2path (int l1,int l2, Hmm *H, double ***M)
 536 {
 537   int i, j,l;
 538   int *path;
 539   HmmState *S1;
 540   int cs;
 541
 542   l=0;
 543   path=vcalloc (l1+l2+1, sizeof (int));
 544   i=l1;j=l2;
 545   l=0;
 546   cs=M[H->nS+H->end][i][j];
 547
 548   while (i>0|| j>0)
 549     {
 550
 551       S1=H->S[cs];
 552       path[l++]=cs;
 553
 554       cs=M[H->nS+cs][i][j];
 555       i-=S1->DI;
 556       j-=S1->DJ;
 557       /*fprintf ( stderr, "%d", cs);*/
 558     }
 559   invert_list_int (path, l);
 560   path[l++]=H->forbiden;
 561
 562   return path;
 563 }
 564
 565 /*********************************************************************************/
 566 /*                                                                               */
 567 /*                                                                               */
 568 /*                      HMM Viterbi Linear                                       */
 569 /*                                                                               */
 570 /*                                                                               */
 571 /*********************************************************************************/
 572
 573
 574 int viterbiL_pair_wise (Alignment *A,int*ns, int **ls,Constraint_list *CL)
 575 {
 576   int l1, l2;
 577   int *path;
 578   Hmm * H;
 579   MatState *Start;
 580   MatState *End;
 581
 582   A->pos=aln2pos_simple( A, -1, ns, ls);
 583   Start=ManageMatState ( DECLARE, NULL);
 584   End=ManageMatState ( DECLARE, NULL);
 585   H=define_simple_model (CL);
 586   l1=strlen (A->seq_al[ls[0][0]]);
 587   l2=strlen (A->seq_al[ls[1][0]]);
 588
 589
 590
 591   Start->i=0 ;Start->j=0 ; Start->st=H->start;Start->sc=0;
 592   End->i  =l1;  End->j=l2; End  ->st=H->end;
 593   Start=RviterbiL_hmm (A, ns, ls, H, CL, Start,End);
 594   path=traceback (A, ns, ls, H, CL, Start,NULL, NULL);
 595
 596   A=viterbipath2aln (A,ns,ls,path, H);
 597
 598   free_Hmm (H);
 599   free_int (A->pos, -1);
 600   A->pos=NULL;
 601
 602   return A->score_aln;
 603 }
 604
 605
 606 MatState* RviterbiL_hmm (Alignment *A,int *ns, int **ls, Hmm *H, Constraint_list *CL,MatState *S, MatState *E)
 607 {
 608   MatState *MS, *ME;
 609   MS=S;
 610   ME=E;
 611
 612   viterbiL_hmm (A,ns, ls,H, CL, S, E);
 613
 614
 615   if ( S->n==E)return S;
 616   if ( E->sc==H->forbiden)
 617         {
 618           DisplayHmm (H);
 619           fprintf ( stderr, "\nERROR: The Requested Model (Cf Above) Cannot Produce the Pair-Alignment\nYou must allow extra possible transitions\n[FATAL:%s]", PROGRAM);
 620           myexit ( EXIT_FAILURE);
 621         }
 622   E=S->n;
 623
 624   while (S!=ME)
 625     {
 626       int d1, d2, align;
 627       d1=MinDeltaMatState(S,E);
 628       d2=MaxDeltaMatState(S,E);
 629       align=((d1==1 && d2==1) || ( d1==0))?0:1;
 630       if (align)RviterbiL_hmm (A,ns, ls,H, CL,S,E);
 631       S=E;
 632       E=S->n;
 633     }
 634   return MS;
 635 }
 636
 637 #define Dim(i,j) (i*LenJ+j)
 638 MatState* viterbiL_hmm (Alignment *A,int *ns, int **ls, Hmm *H, Constraint_list *CL, MatState *S,MatState *E)
 639 {
 640   int current, previous,row, prow;
 641   double v;
 642   int a,i,j,pi,pj, s, k;
 643   int start_i, start_j, end_i, end_j, l1, l2;
 644   HmmState *S1, *S2;
 645   MatState *CC, *PCC,*tS, *tE, *mark=NULL;
 646   int midpoint;
 647
 648
 649   static MatState ***M;
 650
 651
 652   static int LenJ, LenI;
 653   int MaxDelta=50, DeltaI, DeltaJ;
 654
 655
 656
 657   DisplayMatState (S, "\n\tS");
 658   DisplayMatState (E, "\n\tE");
 659
 660
 661   if ( A==NULL)
 662     {
 663       for ( a=0; a<2; a++)memset(M[a],0,LenJ*LenI*sizeof (MatState*));
 664       free_arrayN((void **)M, 3);M=NULL;
 665       ManageMatState ( FREE_ALL, NULL);
 666       return NULL;
 667     }
 668
 669
 670   if ( MatStateAreIdentical ( S, E))return NULL;
 671   l1=strlen (A->seq_al[ls[0][0]]);l2=strlen (A->seq_al[ls[1][0]]);
 672
 673   midpoint=S->i+((E->i-S->i)/2);
 674   DeltaI=E->i-S->i;
 675   DeltaJ=E->j-S->j;
 676
 677   start_i=S->i;end_i=E->i;start_j=S->j;end_j=E->j;
 678   current=0;previous=1;
 679
 680
 681   if ( !M)
 682     {
 683       LenI=l2+1;
 684       LenJ=H->nS;
 685       M=declare_arrayN(3, sizeof ( MatState),2,LenI*LenJ,0);
 686     }
 687
 688
 689   /*MAKE THE VITERBI FROM S(tart) to E(nd)*/
 690   mark=ManageMatState ( MARK, mark);
 691   for (i=start_i; i<=end_i; i++)
 692     {
 693       row=current;
 694       for ( j=start_j; j<=end_j; j++)
 695         {
 696           DeltaJ=((FABS(j-i))<MaxDelta)?1:MaxDelta+1; /*Keep Pointers on a DealtaMax Band around the main diagonal*/
 697           for ( s=H->nS-1;s>=0; s--)
 698             {
 699               S1=H->S[s];pi=i-S1->DI;prow=S1->DI;pj=j-S1->DJ;
 700
 701               CC=M[row][Dim(j,s)]=CopyMatState(NULL, M[row][Dim(j,s)]);
 702               CC->i=i; CC->j=j; CC->st=s;PCC=NULL;
 703
 704               if (i==start_i && j==start_j && s==S->st){CC=CopyMatState(S,CC);}
 705               else if ( i==end_i && j==end_j && E->st!=H->end && s!=E->st)CC->sc=H->forbiden;
 706               else if ( pi<start_i || pj<start_j)        {CC->sc=H->forbiden;}
 707               else
 708                 {
 709                   for (k=1; k<=H->fromM[S1->state][0]; k++)
 710                     {
 711                       S2=H->S[H->fromM[s][k]];
 712                       PCC=M[prow][Dim((j-S1->DJ),(S2->state))];
 713
 714                       if ( !PCC)PCC=NULL;
 715                       else if      ( pi+pj!=0 && S2->state==H->start);
 716                       else if ( !(pi==l1 && pj==l2) && s==H->end);
 717                       else
 718                         {
 719
 720                           v=hmm_add(CC->sc,H->T[PCC->st][CC->st]);
 721
 722                           v=lu_RviterbiD_hmm(A,ns, ls, H, CL,PCC,CC, NULL);
 723                           if ( v!=H->forbiden && (CC->sc==H->forbiden || v> CC->sc)){CC->sc=v; CC->pst=S2->state;CC->p=PCC;}
 724                         }
 725                     }
 726                 }
 727               if (CC->sc==H->forbiden);
 728               else if (i==midpoint || DeltaI<=MaxDelta||DeltaJ<=MaxDelta ||(i==start_i && j==start_j && s==S->st) )
 729                 {
 730                   CC->m=(CC->p)?(CC->p)->m:NULL;
 731                   PCC=CopyMatState(CC,NULL);
 732                   PCC->m=CC->m;CC->m=PCC;
 733                 }
 734               else CC->m=(CC->p)?(CC->p)->m:NULL;
 735             }
 736         }
 737       prow=previous;
 738       for ( j=start_j; j<=end_j && i!=end_i; j++)
 739         {
 740           for ( s=H->nS-1;s>=0; s--)
 741             {
 742
 743               CC=(M[prow][Dim(j,s)]);M[prow][Dim(j,s)]=M[row][Dim(j,s)];M[row][Dim(j,s)]=CC;
 744               if (M[prow][Dim(j,s)]) M[row ][Dim(j,s)]=CopyMatState ( M[prow][Dim(j,s)], M[row][Dim(j,s)]);
 745
 746             }
 747         }
 748
 749     }
 750
 751   mark=ManageMatState ( MARK,mark);
 752   row=current;
 753
 754
 755   if ( E->st==H->end || E->st==H->forbiden){E=CopyMatState ((M[row][Dim(end_j,E->st)]),E);}
 756
 757
 758
 759
 760   PCC=CopyMatState (M[row][Dim(end_j,E->st)], NULL);
 761
 762   if ( MatStateAreIdentical(PCC,PCC->m))PCC=PCC->m;
 763     tS=tE=PCC;
 764   while (PCC->m)
 765     {
 766       tS=CopyMatState (PCC->m,NULL); tS->n=PCC; PCC->p=tS;PCC=tS;
 767     }
 768
 769   if (tS==tE);
 770   else
 771     {
 772       S->n=tS->n; (S->n)->p=S;
 773       E->p=tE->p; (E->p)->n=E;
 774     }
 775   for ( a=0; a<2; a++)memset(M[a],0,LenJ*LenI*sizeof (MatState*));
 776   ManageMatState ( FREE_MARK,mark);
 777
 778
 779   while (S && S->p!=E){S->m=NULL;S=S->n;}/*Clean the memory of the rturned Cells*/
 780   return NULL;
 781 }
 782
 783 /*********************************************************************************/
 784 /*                                                                               */
 785 /*                                                                               */
 786 /*                      HMM Viterbi Diagonals                                    */
 787 /*                                                                               */
 788 /*                                                                               */
 789 /*********************************************************************************/
 790 int viterbiD_pair_wise (Alignment *A,int*ns, int **ls,Constraint_list *CL)
 791 {
 792   int l1, l2;
 793   int *path;
 794   Hmm * H;
 795   MatState *Start;
 796   MatState *End;
 797   int **seg_list;
 798   int a, b, c;
 799   int main_i;
 800   int main_j;
 801
 802
 803   A->pos=aln2pos_simple( A, -1, ns, ls);
 804
 805   Start=ManageMatState ( DECLARE, NULL);
 806   End=ManageMatState ( DECLARE, NULL);
 807   H=define_simple_model (CL);
 808   l1=strlen (A->seq_al[ls[0][0]]);
 809   l2=strlen (A->seq_al[ls[1][0]]);
 810
 811   main_i=MAX(1,(l2-l1)+1);
 812   main_j=MAX(1,(l1-l2)+1);
 813
 814   seg_list=declare_arrayN(2, sizeof (int), l1+l2+3, 3);
 815   seg_list[0][0]=DIAGONALS;
 816
 817
 818   c=1;
 819   for ( b=1,a=l1; a>= 1; a--)
 820     {
 821       if (a<50 || (b==main_i && a==main_j))
 822         {
 823         seg_list[c][0]=a;
 824         seg_list[c][1]=b;
 825         seg_list[c][2]=MIN((l1-a), (l2-b));
 826         c++;
 827         }
 828       }
 829
 830
 831   for ( b=2,a=1; b<= l2; b++, c++)
 832     {
 833       if (b<50 || (b==main_i && a==main_j))
 834         {
 835           seg_list[c][0]=a;
 836           seg_list[c][1]=b;
 837           seg_list[c][2]=MIN((l1-a), (l2-b));
 838         }
 839       }
 840
 841
 842   seg_list[c][0]=FORBIDEN;
 843
 844   Start->i=0 ;Start->j=0 ; Start->st=H->start;Start->sc=0;
 845   End->i  =l1;  End->j=l2; End  ->st=H->end;
 846   Start=RviterbiD_hmm (A, ns, ls, H, CL, Start,End,seg_list);
 847
 848
 849   path=traceback (A, ns, ls, H, CL, Start,NULL, NULL);
 850
 851
 852
 853   A=viterbipath2aln (A,ns,ls,path, H);
 854
 855   viterbiD_hmm (NULL, ns, ls, H, CL, Start,End, seg_list);
 856   free_Hmm (H);
 857   free_int (A->pos, -1);
 858   free_arrayN((void **)seg_list, 2);
 859
 860   A->pos=NULL;
 861   return A->score_aln;
 862 }
 863
 864
 865 double lu_RviterbiD_hmm (Alignment *A,int *ns, int **ls, Hmm *H, Constraint_list *CL,MatState *S, MatState *E, int **seg_list)
 866 {
 867   HmmState *S1;
 868   double sc, sc2,e, t;
 869   static MatState *cS=NULL, *cE=NULL;
 870   double min, max;
 871   max=MAX((E->i-S->i), (E->j-S->j));
 872   min=MIN((E->i-S->i), (E->j-S->j));
 873
 874
 875   if ( S->sc==H->forbiden) return H->forbiden;
 876   else if (min==0)
 877     {
 878       e=hmm_add(S->sc,H->T[S->st][E->st]);
 879       if (  H->T[E->st][E->st]!=H->forbiden)e=hmm_add(e, (max-1)*H->T[E->st][E->st]);
 880       if ( (H->S[E->st])->em!=H->forbiden)  e=hmm_add(e, max    *(H->S[E->st])->em );
 881       return e;
 882     }
 883   else if ( min>0 && max>1)
 884     {
 885
 886       fprintf ( stderr, "\nWarning: Disjoined Diagonals");
 887       DisplayMatState (S, "\n\tS");
 888       DisplayMatState (E, "\n\tE");
 889
 890
 891       cS=CopyMatState ( S,cS);
 892       cE=CopyMatState ( E,cE);
 893       cE->sc=H->forbiden;
 894       viterbiD_hmm (A,ns,ls, H,CL,cS, cE, NULL);
 895       sc2=cE->sc;
 896
 897       return sc2;
 898     }
 899   else
 900     {
 901       S1=H->S[E->st];
 902       t=H->T[S->st][E->st];
 903       e=(S1->em==H->forbiden)?S1->em_func (A, A->pos, ns[0], ls[0],E->i-1, A->pos,ns[1], ls[1], E->j-1, CL):S1->em;
 904       sc=hmm_add(S->sc,t);
 905       sc=hmm_add(sc,e);
 906       return sc;
 907     }
 908   return H->forbiden;
 909 }
 910
 911
 912 MatState* RviterbiD_hmm (Alignment *A,int *ns, int **ls, Hmm *H, Constraint_list *CL,MatState *S, MatState *E, int **seg_list)
 913 {
 914   MatState *MS, *ME;
 915   MS=S;
 916   ME=E;
 917
 918   viterbiD_hmm (A,ns, ls,H, CL, S, E, seg_list);
 919
 920
 921   if ( S->n==E)return S;
 922   if ( E->sc==H->forbiden)
 923         {
 924           DisplayHmm (H);
 925           fprintf ( stderr, "\nERROR: The Requested Model (Cf Above) Cannot Produce the Pair-Alignment\nYou must allow extra possible transitions\n[FATAL:%s]", PROGRAM);
 926           myexit ( EXIT_FAILURE);
 927         }
 928   E=S->n;
 929
 930   while (S!=ME)
 931     {
 932       int d1, d2, align;
 933       d1=MinDeltaMatState(S,E);
 934       d2=MaxDeltaMatState(S,E);
 935       align=((d1==1 && d2==1) || ( d1==0))?0:1;
 936       if (align)RviterbiD_hmm (A,ns, ls,H, CL,S,E, seg_list);
 937       S=E;
 938       E=S->n;
 939     }
 940   return MS;
 941 }
 942
 943 #define Dim(i,j) (i*LenJ+j)
 944 MatState* viterbiD_hmm (Alignment *A,int *ns, int **ls, Hmm *H, Constraint_list *CL, MatState *S,MatState *E, int **seg_list)
 945 {
 946   int current, previous,row, prow;
 947   double v;
 948   int a,b,i,j,pi,pj, s, k;
 949   int start_i, start_j, end_i, end_j, l1, l2;
 950   HmmState *S1, *S2;
 951   MatState *CC, *PCC,*tS, *tE, *mark=NULL;
 952   int midpoint;
 953
 954   int dj;
 955   int dc;
 956   int *jlist=NULL;
 957   static int **main_jlist;
 958   static MatState ***M;
 959   static int *toclean;
 960   int ntoclean;
 961   static int LenJ, LenI;
 962   int MaxDelta=50, DeltaI, DeltaJ;
 963   int mode;
 964
 965   DisplayMatState (S, "\n\tS");
 966   DisplayMatState (E, "\n\tE");
 967
 968
 969   if ( A==NULL)
 970     {
 971       free_arrayN((void **)main_jlist, 2);main_jlist=NULL;
 972
 973       for ( a=0; a<2; a++)memset(M[a],0,LenJ*LenI*sizeof (MatState*));
 974       free_arrayN((void **)M, 3);M=NULL;
 975       vfree (toclean);
 976       ManageMatState ( FREE_ALL, NULL);
 977       return NULL;
 978     }
 979
 980
 981   if ( MatStateAreIdentical ( S, E))return NULL;
 982   l1=strlen (A->seq_al[ls[0][0]]);l2=strlen (A->seq_al[ls[1][0]]);
 983
 984   midpoint=S->i+((E->i-S->i)/2);
 985   DeltaI=E->i-S->i;
 986
 987
 988   start_i=S->i;end_i=E->i;start_j=S->j;end_j=E->j;
 989   current=0;previous=1;
 990
 991
 992   if ( !M)
 993     {
 994       LenI=l2+1;
 995       LenJ=H->nS;
 996       M=declare_arrayN(3, sizeof ( MatState),2,LenI*LenJ,0);
 997       toclean=vcalloc ( LenI*LenJ, sizeof (int));
 998     }
 999
1000   if ( !main_jlist)main_jlist= seglist2table(seg_list, l1, l2);
1001
1002
1003   /*MAKE THE VITERBI FROM S(tart) to E(nd)*/
1004   mark=ManageMatState ( MARK, mark);
1005   mode=(!seg_list)?ALL:seg_list[0][0];
1006
1007   for (ntoclean=0,i=start_i; i<=end_i; i++)
1008     {
1009       row=current;
1010
1011       if ( mode==ALL)jlist=main_jlist[0];
1012       else if ( mode==DIAGONALS)jlist=(i==0)?main_jlist[0]:main_jlist[1];
1013       else if ( mode==SEGMENTS)  jlist=main_jlist[i+2];
1014
1015
1016       for ( dj=1; dj<=jlist[0]; dj++)
1017         {
1018           DeltaJ=((FABS(dj-i))<MaxDelta)?1:MaxDelta+1; /*Keep Pointers on a 2*DeltaMax Band around the main diagonal*/
1019
1020           dc=(mode==DIAGONALS && dj!=1)?i:0;/*Make sure the diagonal Mode uses the sides of the array*/
1021           j=jlist[dj]+dc;
1022
1023
1024           if ( j<start_j || j>end_j)continue;
1025           for ( s=H->nS-1;s>=0; s--)
1026             {
1027               S1=H->S[s];pi=i-S1->DI;prow=S1->DI;
1028
1029               if ( S1->DI && S1->DJ){pj=j-S1->DJ;}
1030               else if ( !S1->DJ)pj=j;
1031               else if ( dj>1)pj=jlist[dj-S1->DJ]+dc;
1032               else pj=-1;
1033
1034               if (!M[row][Dim(j,s)])toclean[ntoclean]=Dim(j,s);
1035
1036               CC=M[row][Dim(j,s)]=CopyMatState(NULL, M[row][Dim(j,s)]);
1037               CC->i=i; CC->j=j; CC->st=s;PCC=NULL;
1038
1039               if (i==start_i && j==start_j && s==S->st){CC=CopyMatState(S,CC);}
1040               else if ( i==end_i && j==end_j && E->st!=H->end && s!=E->st)CC->sc=H->forbiden;
1041               else if ( pi<start_i || pj<start_j)        {CC->sc=H->forbiden;}
1042               else
1043                 {
1044                   for (k=1; k<=H->fromM[S1->state][0]; k++)
1045                     {
1046                       S2=H->S[H->fromM[s][k]];
1047
1048                       if ( S1->DI && S1->DJ)PCC=M[prow][Dim((j-S1->DJ),(S2->state))];
1049                       else PCC=M[prow][Dim((jlist[dj-S1->DJ]+dc),(S2->state))];
1050
1051                       if ( !PCC)PCC=NULL;
1052                       else if      ( pi+pj!=0 && S2->state==H->start);
1053                       else if ( !(pi==l1 && pj==l2) && s==H->end);
1054                       else
1055                         {
1056                           v=lu_RviterbiD_hmm(A,ns, ls, H, CL,PCC,CC, NULL);
1057                           if ( v!=H->forbiden && (CC->sc==H->forbiden || v> CC->sc)){CC->sc=v; CC->pst=S2->state;CC->p=PCC;}
1058                         }
1059                     }
1060                 }
1061               if (CC->sc==H->forbiden);
1062               else if (i==midpoint || DeltaI<=MaxDelta||DeltaJ<=MaxDelta ||(i==start_i && j==start_j && s==S->st) )
1063                 {
1064                   CC->m=(CC->p)?(CC->p)->m:NULL;
1065                   PCC=CopyMatState(CC,NULL);
1066                   PCC->m=CC->m;CC->m=PCC;
1067                 }
1068               else CC->m=(CC->p)?(CC->p)->m:NULL;
1069             }
1070         }
1071       prow=previous;
1072       for ( dj=1; dj<=jlist[0] && i!=end_i; dj++)
1073         {
1074           dc=(mode==DIAGONALS && dj!=1)?i:0;
1075           j=jlist[dj]+dc;
1076           if ( j<start_j || j>end_j)continue;
1077
1078           for ( s=H->nS-1;s>=0; s--)
1079             {
1080
1081               CC=(M[prow][Dim(j,s)]);M[prow][Dim(j,s)]=M[row][Dim(j,s)];M[row][Dim(j,s)]=CC;
1082               if (!M[row][Dim(j,s)])toclean[ntoclean++]=Dim(j,s);
1083               if (M[prow][Dim(j,s)]) M[row ][Dim(j,s)]=CopyMatState ( M[prow][Dim(j,s)], M[row][Dim(j,s)]);
1084
1085             }
1086         }
1087
1088     }
1089
1090   mark=ManageMatState ( MARK,mark);
1091   row=current;
1092
1093
1094   if ( E->st==H->end || E->st==H->forbiden){E=CopyMatState ((M[row][Dim(end_j,E->st)]),E);}
1095
1096
1097
1098
1099   PCC=CopyMatState (M[row][Dim(end_j,E->st)], NULL);
1100
1101   if ( MatStateAreIdentical(PCC,PCC->m))PCC=PCC->m;
1102     tS=tE=PCC;
1103   while (PCC->m)
1104     {
1105       tS=CopyMatState (PCC->m,NULL); tS->n=PCC; PCC->p=tS;PCC=tS;
1106     }
1107
1108   if (tS==tE);
1109   else
1110     {
1111       S->n=tS->n; (S->n)->p=S;
1112       E->p=tE->p; (E->p)->n=E;
1113     }
1114
1115   ManageMatState ( FREE_MARK,mark);
1116
1117
1118   for ( a=0; a<ntoclean; a++)
1119     {
1120       for ( b=0; b< 2; b++){M[b][toclean[a]]=NULL;}
1121       toclean[a]=0;
1122     }
1123
1124   while (S && S->p!=E){S->m=NULL;S=S->n;}/*Clean the memory of the rturned Cells*/
1125   return NULL;
1126 }
1127
1128 /*********************************************************************************/
1129 /*                                                                               */
1130 /*                                                                               */
1131 /*                      HMM Viterbi Diagonals  GLOBAL/LOCAL                                  */
1132 /*                                                                               */
1133 /*                                                                               */
1134 /*********************************************************************************/
1135
1136 int viterbiDGL_pair_wise (Alignment *A,int*ns, int **ls,Constraint_list *CL)
1137 {
1138   int l1, l2;
1139   int *path;
1140   Hmm * H;
1141   MatState *Start;
1142   MatState *End;
1143   int **seg_list;
1144   int a, b, c;
1145   int main_i;
1146   int main_j;
1147
1148
1149   A->pos=aln2pos_simple( A, -1, ns, ls);
1150
1151   Start=ManageMatState ( DECLARE, NULL);
1152   End=ManageMatState ( DECLARE, NULL);
1153   H=define_simple_model (CL);
1154   l1=strlen (A->seq_al[ls[0][0]]);
1155   l2=strlen (A->seq_al[ls[1][0]]);
1156
1157   main_i=MAX(1,(l2-l1)+1);
1158   main_j=MAX(1,(l1-l2)+1);
1159
1160   seg_list=declare_arrayN(2, sizeof (int), l1+l2+3, 3);
1161   seg_list[0][0]=DIAGONALS;
1162
1163
1164   c=1;
1165   for ( b=1,a=l1; a>= 1; a--)
1166     {
1167       if (a<50 || (b==main_i && a==main_j))
1168         {
1169         seg_list[c][0]=a;
1170         seg_list[c][1]=b;
1171         seg_list[c][2]=MIN((l1-a), (l2-b));
1172         c++;
1173         }
1174       }
1175
1176
1177   for ( b=2,a=1; b<= l2; b++, c++)
1178     {
1179       if (b<50 || (b==main_i && a==main_j))
1180         {
1181           seg_list[c][0]=a;
1182           seg_list[c][1]=b;
1183           seg_list[c][2]=MIN((l1-a), (l2-b));
1184         }
1185       }
1186
1187
1188   seg_list[c][0]=FORBIDEN;
1189
1190   Start->i=0 ;Start->j=0 ; Start->st=H->start;Start->sc=0;
1191   End->i  =l1;  End->j=l2; End  ->st=H->end;
1192   Start=RviterbiDGL_hmm (A, ns, ls, H, CL, Start,End,seg_list);
1193
1194
1195   path=traceback (A, ns, ls, H, CL, Start,NULL, NULL);
1196
1197
1198
1199   A=viterbipath2aln (A,ns,ls,path, H);
1200
1201   viterbiD_hmm (NULL, ns, ls, H, CL, Start,End, seg_list);
1202   free_Hmm (H);
1203   free_int (A->pos, -1);
1204   free_arrayN((void **)seg_list, 2);
1205
1206   A->pos=NULL;
1207   return A->score_aln;
1208 }
1209
1210
1211 double lu_RviterbiDGL_hmm (Alignment *A,int *ns, int **ls, Hmm *H, Constraint_list *CL,MatState *S, MatState *E, int **seg_list)
1212 {
1213   HmmState *S1;
1214   double sc, sc2,e, t;
1215   static MatState *cS=NULL, *cE=NULL;
1216   double min, max;
1217   max=MAX((E->i-S->i), (E->j-S->j));
1218   min=MIN((E->i-S->i), (E->j-S->j));
1219
1220
1221
1222
1223   if ( S==NULL || E==NULL || S->sc==H->forbiden) return H->forbiden;
1224   else if ( S->st==H->start) return 0;
1225   else if ( E->st==H->end) return S->sc;
1226   else if (min==0)
1227     {
1228       e=hmm_add(S->sc,H->T[S->st][E->st]);
1229       if (  H->T[E->st][E->st]!=H->forbiden)e=hmm_add(e, (max-1)*H->T[E->st][E->st]);
1230       if ( (H->S[E->st])->em!=H->forbiden)  e=hmm_add(e, max    *(H->S[E->st])->em );
1231       return e;
1232     }
1233   else if ( min>0 && max>1)
1234     {
1235
1236       fprintf ( stderr, "\nWarning: Disjoined Diagonals");
1237       DisplayMatState (S, "\n\tS");
1238       DisplayMatState (E, "\n\tE");
1239
1240
1241       cS=CopyMatState ( S,cS);
1242       cE=CopyMatState ( E,cE);
1243       cE->sc=H->forbiden;
1244       viterbiD_hmm (A,ns,ls, H,CL,cS, cE, NULL);
1245       sc2=cE->sc;
1246
1247       return sc2;
1248     }
1249   else
1250     {
1251       S1=H->S[E->st];
1252       t=H->T[S->st][E->st];
1253       e=(S1->em==H->forbiden)?S1->em_func (A, A->pos, ns[0], ls[0],E->i-1, A->pos,ns[1], ls[1], E->j-1, CL):S1->em;
1254       sc=hmm_add(S->sc,t);
1255       sc=hmm_add(sc,e);
1256       return sc;
1257     }
1258   return H->forbiden;
1259 }
1260
1261
1262 MatState* RviterbiDGL_hmm (Alignment *A,int *ns, int **ls, Hmm *H, Constraint_list *CL,MatState *S, MatState *E, int **seg_list)
1263 {
1264   MatState *MS, *ME;
1265   MS=S;
1266   ME=E;
1267
1268   viterbiDGL_hmm (A,ns, ls,H, CL, S, E, seg_list);
1269
1270
1271   if ( S->n==E)return S;
1272   if ( E->sc==H->forbiden)
1273         {
1274           DisplayHmm (H);
1275           fprintf ( stderr, "\nERROR: The Requested Model (Cf Above) Cannot Produce the Pair-Alignment\nYou must allow extra possible transitions\n[FATAL:%s]", PROGRAM);
1276           myexit ( EXIT_FAILURE);
1277         }
1278   E=S->n;
1279
1280   while (S!=ME)
1281     {
1282       int d1, d2, align;
1283       d1=MinDeltaMatState(S,E);
1284       d2=MaxDeltaMatState(S,E);
1285       align=((d1==1 && d2==1) || ( d1==0))?0:1;
1286       if (align)RviterbiDGL_hmm (A,ns, ls,H, CL,S,E, seg_list);
1287       S=E;
1288       E=S->n;
1289     }
1290   return MS;
1291 }
1292
1293
1294 #define Dim(i,j) (i*LenJ+j)
1295 MatState* viterbiDGL_hmm (Alignment *A,int *ns, int **ls, Hmm *H, Constraint_list *CL, MatState *S,MatState *E, int **seg_list)
1296 {
1297   int current, previous,row, prow;
1298   double v;
1299   int a,i,j,pi,pj, s, k;
1300   int start_i, start_j, end_i, end_j, l1, l2;
1301   HmmState *S1, *S2;
1302   MatState *CC, *PCC,*tS, *tE,*bestE,*bestS, *mark=NULL;
1303   int midpoint;
1304
1305   int dj;
1306   int dc;
1307   int *jlist=NULL;
1308   static int **main_jlist;
1309   static MatState ***M;
1310   static int *toclean;
1311   int ntoclean;
1312   static int LenJ, LenI;
1313   int MaxDelta=50, DeltaI, DeltaJ;
1314   int mode;
1315
1316
1317
1318   DisplayMatState (S, "\n\tS");
1319   DisplayMatState (E, "\n\tE");
1320
1321
1322   if ( A==NULL)
1323     {
1324       free_arrayN((void **)main_jlist, 2);main_jlist=NULL;
1325
1326       for ( a=0; a<2; a++)memset(M[a],0,LenJ*LenI*sizeof (MatState*));
1327       free_arrayN((void **)M, 3);M=NULL;
1328       vfree (toclean);
1329       ManageMatState ( FREE_ALL, NULL);
1330       return NULL;
1331     }
1332
1333
1334   if ( MatStateAreIdentical ( S, E))return NULL;
1335   l1=strlen (A->seq_al[ls[0][0]]);l2=strlen (A->seq_al[ls[1][0]]);
1336
1337   midpoint=S->i+((E->i-S->i)/2);
1338   DeltaI=E->i-S->i;
1339
1340
1341   start_i=S->i;end_i=E->i;start_j=S->j;end_j=E->j;
1342   current=0;previous=1;
1343
1344
1345   if ( !M)
1346     {
1347       LenI=l2+1;
1348       LenJ=H->nS;
1349       M=declare_arrayN(3, sizeof ( MatState),2,LenI*LenJ,0);
1350       toclean=vcalloc ( LenI*LenJ, sizeof (int));
1351     }
1352
1353   if ( !main_jlist)main_jlist= seglist2table(seg_list, l1, l2);
1354
1355
1356   /*MAKE THE VITERBI FROM S(tart) to E(nd)*/
1357   mark=ManageMatState ( MARK, mark);
1358   mode=(!seg_list)?ALL:seg_list[0][0];
1359   bestE=CopyMatState (E, NULL);
1360   bestS=CopyMatState (NULL, NULL);
1361   for (ntoclean=0,i=start_i; i<=end_i; i++)
1362     {
1363       row=current;
1364
1365       if ( mode==ALL)jlist=main_jlist[0];
1366       else if ( mode==DIAGONALS)jlist=(i==0)?main_jlist[0]:main_jlist[1];
1367       else if ( mode==SEGMENTS)  jlist=main_jlist[i+2];
1368
1369
1370       for ( dj=1; dj<=jlist[0]; dj++)
1371         {
1372           DeltaJ=(FABS(dj-i)<MaxDelta)?1:MaxDelta+1; /*Keep Pointers on a 2*DeltaMax Band around the main diagonal*/
1373
1374           dc=(mode==DIAGONALS && dj!=1)?i:0;/*Make sure the diagonal Mode uses the sides of the array*/
1375           j=jlist[dj]+dc;
1376
1377
1378           if ( j<start_j || j>end_j)continue;
1379           for ( s=H->nS-1;s>=0; s--)
1380             {
1381               if ( s==S->st)continue;
1382               S1=H->S[s];pi=i-S1->DI;prow=S1->DI;
1383
1384               if ( S1->DI && S1->DJ){pj=j-S1->DJ;}
1385               else if ( !S1->DJ)pj=j;
1386               else if ( dj>1)pj=jlist[dj-S1->DJ]+dc;
1387               else pj=-1;
1388
1389               if (!M[row][Dim(j,s)])toclean[ntoclean]=Dim(j,s);
1390
1391               CC=M[row][Dim(j,s)]=CopyMatState(NULL, M[row][Dim(j,s)]);
1392               CC->i=i; CC->j=j; CC->st=s;PCC=NULL;
1393
1394               if (i==start_i && j==start_j && s==S->st){CC=CopyMatState(S,CC);}
1395               else if ( s==S->st);
1396               else if ( i==end_i && j==end_j && E->st!=H->end && s!=E->st)CC->sc=H->forbiden;
1397
1398               else if ( pi<start_i || pj<start_j)        {CC->sc=H->forbiden;}
1399               else
1400                 {
1401                   for (k=1; k<=H->fromM[S1->state][0]; k++)
1402                     {
1403                       S2=H->S[H->fromM[s][k]];
1404
1405                       if ( S1->DI && S1->DJ)PCC=M[prow][Dim((j-S1->DJ),(S2->state))];
1406                       else PCC=M[prow][Dim((jlist[dj-S1->DJ]+dc),(S2->state))];
1407
1408                       if ( S2->state==H->start){PCC=bestS;PCC->st=0;PCC->sc=0;PCC->m=PCC->n=PCC->p=NULL;}
1409
1410                       v=lu_RviterbiDGL_hmm(A,ns, ls, H, CL,PCC,CC, NULL);
1411                       if ( v!=H->forbiden && (CC->sc==H->forbiden || v> CC->sc)){CC->sc=v; CC->pst=S2->state;CC->p=PCC;}
1412                     }
1413                 }
1414               if ( CC->sc==H->forbiden);
1415               else if ( bestE->sc==H->forbiden || bestE->sc>CC->sc)
1416                 {
1417                   bestE=CopyMatState(CC, bestE);
1418                   bestE->m=(CC->p)->m;
1419                 }
1420               else if (CC->p && (CC->p)->st==H->start)
1421                 {
1422                   CC->m=CopyMatState (CC->p, NULL);
1423                 }
1424               else if (i==midpoint || DeltaI<=MaxDelta||DeltaJ<=MaxDelta ||(i==start_i && j==start_j && s==S->st) )
1425                 {
1426                   CC->m=(CC->p)?(CC->p)->m:NULL;
1427                   PCC=CopyMatState(CC,NULL);
1428                   PCC->m=CC->m;CC->m=PCC;
1429                 }
1430               else CC->m=(CC->p)?(CC->p)->m:NULL;
1431             }
1432         }
1433       prow=previous;
1434       for ( dj=1; dj<=jlist[0] && i!=end_i; dj++)
1435         {
1436           dc=(mode==DIAGONALS && dj!=1)?i:0;
1437           j=jlist[dj]+dc;
1438           if ( j<start_j || j>end_j)continue;
1439
1440           for ( s=H->nS-1;s>=0; s--)
1441             {
1442
1443               CC=(M[prow][Dim(j,s)]);M[prow][Dim(j,s)]=M[row][Dim(j,s)];M[row][Dim(j,s)]=CC;
1444               /*if (!M[row][Dim(j,s)])toclean[ntoclean++]=Dim(j,s);*/
1445               if (M[prow][Dim(j,s)]) M[row ][Dim(j,s)]=CopyMatState ( M[prow][Dim(j,s)], M[row][Dim(j,s)]);
1446
1447             }
1448         }
1449
1450     }
1451
1452   mark=ManageMatState ( MARK,mark);
1453   row=current;
1454
1455
1456   if ( E->st==H->end || E->st==H->forbiden){E=CopyMatState ((M[row][Dim(end_j,E->st)]),E);}
1457   PCC=CopyMatState (bestE, NULL);
1458
1459   if ( MatStateAreIdentical(PCC,PCC->m))PCC=PCC->m;
1460     tS=tE=PCC;
1461   while (PCC->m)
1462     {
1463       tS=CopyMatState (PCC->m,NULL); tS->n=PCC; PCC->p=tS;PCC=tS;
1464     }
1465
1466   if (tS==tE);
1467   else
1468     {
1469       CopyMatState ( tS, S);
1470       CopyMatState ( tE, E);
1471     }
1472   ManageMatState ( FREE_MARK,mark);
1473
1474   for ( a=0; a<2; a++)memset(M[a],0,LenJ*LenI*sizeof (MatState*));
1475
1476   while (S && S->p!=E){S->m=NULL;S=S->n;}/*Clean the memory of the rturned Cells*/
1477   return NULL;
1478 }
1479
1480
1481 /*********************************************************************************/
1482 /*                                                                               */
1483 /*                                                                               */
1484 /*                      HMM Viterbi Diagonals  PROCESSING                        */
1485 /*                                                                               */
1486 /*                                                                               */
1487 /*********************************************************************************/
1488 int **seglist2table ( int **seglist,int l1, int l2)
1489   {
1490     int **valuesT;
1491     int *bvalues;
1492     int line, a,si, sj, ei, j, c;
1493
1494     /*All: 0*/
1495     valuesT=vcalloc ((l1+2)+3, sizeof (int*));
1496     valuesT[0]=vcalloc (l2+2, sizeof (int));
1497     for (a=0; a<=l2; a++)valuesT[0][++valuesT[0][0]]=a;
1498
1499     if ( !seglist) return valuesT;
1500     /*Diagonals: 1*/
1501     valuesT[1]=vcalloc (l1+l2+2, sizeof (int));
1502     bvalues=vcalloc (l1+l2+2, sizeof (int));
1503     c=1;
1504     while (seglist[c][0]!=FORBIDEN)
1505       {
1506
1507         si=seglist[c][0];
1508         sj=seglist[c][1];
1509
1510         bvalues[(sj-si)+l1]=1;
1511         c++;
1512       }
1513     valuesT[1][++valuesT[1][0]]=0;
1514     for (a=0; a<=(l1+l2); a++)
1515       {
1516         if (bvalues[a])
1517           {
1518             valuesT[1][++valuesT[1][0]]=a-l1;
1519           }
1520
1521       }
1522     vfree (bvalues);
1523
1524     /*Segments: 2*/
1525     valuesT[2]=vcalloc (l2+2, sizeof (int));
1526     for (a=0; a<=l2; a++)valuesT[2][++valuesT[2][0]]=a;
1527
1528     bvalues=vcalloc (l2+2, sizeof (int));
1529     for ( line=1; line<=l1; line++)
1530       {
1531         bvalues[0]=c=0;
1532         bvalues[++bvalues[0]]=0;
1533         while (seglist[c][0]!=FORBIDEN)
1534           {
1535                 si=seglist[c][0];
1536                 ei=si+seglist[c][2];
1537                 sj=seglist[c][1];
1538                 j=sj+(line-si);
1539                 if ( line<si || line>ei);
1540                 else if (j>=0 && j<=l2 && seglist[c][2])
1541                   {
1542                     bvalues[++bvalues[0]]=j;
1543                   }
1544                 c++;
1545           }
1546         valuesT[line+2]=vcalloc (bvalues[0]+1, sizeof (int));
1547         for ( a=0; a<=bvalues[0]; a++) valuesT[line+2][a]=bvalues[a];
1548       }
1549     vfree (bvalues);
1550     return valuesT;
1551   }
1552
1553
1554
1555 /*********************************************************************************/
1556 /*                                                                               */
1557 /*                                                                               */
1558 /*                      HMM modeling                                             */
1559 /*                                                                               */
1560 /*                                                                               */
1561 /*********************************************************************************/
1562
1563 Hmm* declare_hmm(int n)
1564   {
1565     Hmm *H;
1566     int a, b;
1567
1568     H=vcalloc (1, sizeof (Hmm));
1569     H->nS=n;
1570     H->S=vcalloc (H->nS, sizeof (HmmState*));
1571     for (a=0; a<H->nS; a++)
1572       {
1573         H->S[a]=vcalloc (1, sizeof (HmmState));
1574         (H->S[a])->em2=vcalloc (MAX_EMISSION, sizeof (double));
1575
1576         (H->S[a])->T=vcalloc ( H->nS, sizeof (StateTrans*));
1577         for ( b=0; b< H->nS; b++)
1578           (H->S[a])->T[b]=vcalloc (1, sizeof (StateTrans));
1579       }
1580     return H;
1581   }
1582
1583 Hmm* free_Hmm(Hmm*H)
1584   {
1585     int a, b;
1586
1587     H=vcalloc (1, sizeof (Hmm));
1588     free_double (H->T, -1);
1589     free_int ( H->fromM, -1);
1590     free_int ( H->toM, -1);
1591
1592     for (a=0; a< H->nS; a++)
1593       {
1594
1595         for ( b=0; b< H->nS; b++)
1596           {
1597             vfree ((H->S[a])->em2);
1598             vfree((H->S[a])->T[b]);
1599           }
1600         vfree((H->S[a])->T);
1601         vfree(H->S[a]);
1602       }
1603     vfree (H->S);
1604     vfree (H);
1605     return NULL;
1606   }
1607
1608 void DisplayHmm ( Hmm *H)
1609 {
1610   int a, b;
1611   HmmState *S1, *S2;
1612
1613   for ( a=0; a< H->nS; a++)
1614     {
1615       S1=H->S[a];
1616       fprintf ( stderr, "\nSTATE %d: %s\n",S1->state,S1->name);
1617       fprintf ( stderr, "\n\tDI %d", S1->DI);
1618       fprintf ( stderr, "\n\tDJ %d", S1->DJ);
1619       fprintf ( stderr, "\n\tE  %f", (float)S1->em);
1620
1621       fprintf ( stderr, "\nReached FROM: ");
1622       for ( b=1; b<=H->fromM[a][0]; b++)
1623         {
1624           S2=H->S[H->fromM[a][b]];
1625           fprintf ( stderr, "[ %s %f ] ", S2->name, H->T[S2->state][S1->state]);
1626         }
1627       fprintf ( stderr, "\nGoes TO: ");
1628       for ( b=1; b<=H->toM[a][0]; b++)
1629         {
1630           S2=H->S[H->toM[a][b]];
1631           fprintf ( stderr, "[ %s %f ] ", S2->name, H->T[S1->state][S2->state]);
1632         }
1633     }
1634   return;
1635 }
1636 Hmm * bound_hmm ( Hmm *H)
1637 {
1638   int a, b, c;
1639   char **name;
1640   HmmState *S;
1641
1642   name=declare_char(H->nS, 100);
1643   H->T=declare_double ( H->nS, H->nS);
1644
1645   for ( a=0; a< H->nS; a++)
1646     {
1647       sprintf ( name[a], "%s", (H->S[a])->name);
1648       H->order=MAX(H->order, (H->S[a])->DI);
1649       H->order=MAX(H->order, (H->S[a])->DJ);
1650     }
1651
1652   for ( a=0; a< H->nS; a++)for (b=0; b<H->nS; b++)H->T[a][b]=H->forbiden;
1653   for (a=0; a< H->nS; a++)
1654     {
1655       S=H->S[a];
1656       for ( b=0; b< S->nT; b++)
1657         {
1658           c=name_is_in_list ((S->T[b])->name, name, H->nS, 100);
1659           if ( c!=-1)H->T[a][c]=(S->T[b])->tr;
1660         }
1661     }
1662
1663   /*Bound the model:
1664     bM[state][0]=n_allowed transitions
1665     bM[state][1]=first allowed transition
1666   */
1667
1668   H->toM=declare_int ( H->nS, H->nS);
1669   H->fromM=declare_int ( H->nS, H->nS);
1670
1671   for ( a=0; a< H->nS; a++)
1672     for ( b=0; b< H->nS; b++)
1673       {
1674         if ( H->T[a][b]!=H->forbiden )
1675           {
1676             {H->fromM[b][0]++; H->fromM[b][H->fromM[b][0]]=a;}
1677             {H->toM[a][0]++; H->toM[a][H->toM[a][0]]=b;}
1678           }
1679       }
1680   for ( a=0; a< H->nS; a++)
1681     {
1682       if (( H->S[a])->em!=H->forbiden)( H->S[a])->em*=SCORE_K;
1683       for ( b=0; b< H->nS; b++)
1684         if ( H->T[a][b]!=H->forbiden)H->T[a][b]*=SCORE_K;
1685     }
1686   free_arrayN((void**)name, 2);
1687   return H;
1688 }
1689
1690 /*********************************************************************************/
1691 /*                                                                               */
1692 /*                                                                               */
1693 /*                      Memory Management                                        */
1694 /*                                                                               */
1695 /*                                                                               */
1696 /*********************************************************************************/
1697
1698 MatState * ManageMatState(int Mode, MatState *C)
1699 {
1700   static MatState *Fheap;
1701   static MatState *Aheap;
1702   MatState *Cmark, *Pmark;
1703   static int alloc, free;
1704   if (!Fheap || Fheap->Hp==NULL)
1705     {
1706       int c=0;
1707       int extension=1000;
1708       if (!Fheap){Fheap=vcalloc (1, sizeof (MatState));Fheap->free=1;free++;}
1709       if (!Aheap)Aheap=vcalloc (1, sizeof (MatState));
1710       while ( c!=extension)
1711         {
1712           C=vcalloc ( 1, sizeof (MatState));
1713           C->free=1;Fheap->Hn=C;C->Hp=Fheap;
1714           Fheap=C;
1715           c++;
1716           free++;
1717         }
1718     }
1719
1720   if ( Mode==DECLARE)
1721     {
1722
1723       C=Fheap;
1724       Fheap=Fheap->Hp;
1725       C->Hn=C->Hp=NULL;
1726       if ( Aheap){Aheap->Hn=C;C->Hp=Aheap;Aheap=C;}
1727       else Aheap=C;
1728       alloc++;
1729       free--;
1730       C->free=0;
1731       C=CopyMatState(NULL, C);
1732       return C;
1733     }
1734   else if ( Mode==FREE)
1735     {
1736       if ( !C || C->free==1);
1737       else
1738         {
1739           C=CopyMatState(NULL, C);
1740           C->free=1;
1741           if (C->Hp==NULL && C==Aheap)crash ("");
1742           if (C==Aheap)Aheap=C->Hp;
1743           if (C->Hn){(C->Hn)->Hp=C->Hp;}
1744           if (C->Hp){(C->Hp)->Hn=C->Hn;}
1745           C->Hp=C->Hn=NULL;
1746           Fheap->Hn=C;C->Hp=Fheap;
1747           Fheap=C;
1748           alloc--;
1749           free++;
1750         }
1751       return NULL;
1752     }
1753   else if ( Mode==FREE_ALL)
1754     {
1755       while ( Aheap)
1756         {
1757           C=Aheap->Hp;
1758           vfree (Aheap);
1759           Aheap=C;
1760         }
1761       while ( Fheap)
1762         {
1763           C=Fheap->Hp;
1764           vfree (Fheap);
1765           Fheap=C;
1766         }
1767     }
1768   else if ( Mode==INFO)
1769     {
1770       fprintf ( stderr, "\nAllocated: %d Free %d", alloc, free);
1771     }
1772   else if ( Mode==MARK)
1773     {
1774
1775       if (C==NULL);
1776       else {C->Mn=Aheap;Aheap->Mp=C;}
1777
1778       return Aheap;
1779     }
1780   else if ( Mode==UNMARK)
1781     {
1782       Pmark=Cmark=NULL;
1783     }
1784   else if ( Mode == FREE_MARK)
1785     {
1786       Cmark=C;
1787       Pmark=C->Mp;
1788
1789       if ( Cmark==Pmark)return NULL;
1790       else if ( Cmark==Aheap)
1791         {Aheap=Pmark;C=Pmark->Hn;Pmark->Hn=NULL;}
1792       else
1793         {
1794           (Cmark->Hn)->Hp=Pmark;
1795           C=Pmark->Hn;
1796           Pmark->Hn=Cmark->Hn;
1797         }
1798
1799       Fheap->Hn=C;
1800       C->Hp=Fheap;
1801       Fheap=Cmark;
1802       Fheap->Hn=NULL;
1803
1804       C=Fheap;
1805       while (C && !C->free)
1806         {
1807           free++;alloc--;
1808           C->free=1;
1809           C=C->Hp;
1810         }
1811
1812     }
1813   return NULL;
1814 }
1815
1816
1817 MatState* CopyMatState ( MatState*I, MatState*O)
1818 {
1819   if (O==NULL || O->free==1) O=ManageMatState(DECLARE, NULL);
1820   if (I==NULL || I->free==1)I=NULL;
1821   O->i  =(I)?I->i:0;
1822   O->j  =(I)?I->j:0;
1823   O->st =(I)?I->st:FORBIDEN;
1824   O->pst=(I)?I->pst:FORBIDEN;
1825   O->sc =(I)?I->sc:FORBIDEN;
1826   O->n  =(I)?I->n:NULL;
1827   O->p  =(I)?I->p:NULL;
1828   O->m  =(I)?I->m:NULL;
1829   O->s  =(I)?I->m:NULL;
1830
1831   return O;
1832 }
1833
1834 /*********************************************************************************/
1835 /*                                                                               */
1836 /*                                                                               */
1837 /*                      Comparisons                                              */
1838 /*                                                                               */
1839 /*                                                                               */
1840 /*********************************************************************************/
1841 int MaxDeltaMatState (MatState*S, MatState*E)
1842 {
1843   if ( !S || !E) return -1;
1844   else return MAX((E->i-S->i),(E->j-S->j));
1845 }
1846 int MinDeltaMatState (MatState*S, MatState*E)
1847 {
1848   if ( !S || !E) return -1;
1849   else return MIN((E->i-S->i),(E->j-S->j));
1850 }
1851 int MatStateAreIdentical (MatState*I, MatState*O)
1852 {
1853   if ( !I || !O)return 0;
1854
1855   if ( I->i!=O->i)return 0;
1856   if ( I->j!=O->j)return 0;
1857   if ( I->st!=O->st)return 0;
1858   return 1;
1859 }
1860
1861
1862
1863 Hmm* define_probcons_model(Constraint_list *CL)
1864 {
1865   Hmm *H;
1866   double gop=-10;
1867   double gep=-1;
1868   double lgop=-100;
1869   double lgep=-100;
1870   double freeT=0;
1871   int n=0;
1872   HmmState *S;
1873
1874
1875   H=declare_hmm(7);
1876   H->freeT=freeT=0;
1877
1878   H->forbiden=FORBIDEN;
1879   H->start=START_STATE;
1880   H->end=END_STATE;
1881
1882   /*define START*/
1883   S=H->S[n];
1884   sprintf (S->name, "START"); S->state=n;
1885
1886   S->DI=0;
1887   S->DJ=0;
1888   S->em=freeT;
1889
1890   sprintf ( (S->T[S->nT])->name, "MATCH") ;(S->T[S->nT])->tr=freeT  ;S->nT++;
1891   sprintf ( (S->T[S->nT])->name, "INSERT");(S->T[S->nT])->tr=freeT  ;S->nT++;
1892   sprintf ( (S->T[S->nT])->name, "DELETE");(S->T[S->nT])->tr=freeT  ;S->nT++;
1893   sprintf ( (S->T[S->nT])->name, "END")   ;(S->T[S->nT])->tr=freeT  ;S->nT++;
1894
1895   n++;
1896   /*define END*/
1897   S=H->S[n];
1898   sprintf (S->name, "END"); S->state=n;
1899   S->DI=0;
1900   S->DJ=0;
1901   S->em=freeT;
1902   n++;
1903
1904   /*define Match*/
1905   S=H->S[n];
1906   sprintf (S->name, "MATCH"); S->state=n;
1907   S->DI=1;
1908   S->DJ=1;
1909   S->em=H->forbiden;
1910   S->em_func=CL->get_dp_cost;
1911
1912   sprintf ( (S->T[S->nT])->name, "MATCH") ;(S->T[S->nT])->tr=freeT;S->nT++;
1913   sprintf ( (S->T[S->nT])->name, "INSERT");(S->T[S->nT])->tr=gop  ;S->nT++;
1914   sprintf ( (S->T[S->nT])->name, "DELETE");(S->T[S->nT])->tr=gop  ;S->nT++;
1915   sprintf ( (S->T[S->nT])->name, "END");   (S->T[S->nT])->tr=freeT;S->nT++;
1916
1917   n++;
1918
1919   /*define Insert*/
1920   S=H->S[n];
1921   sprintf (S->name, "INSERT"); S->state=n;
1922   S->DI=1;
1923   S->DJ=0;
1924   S->em=gep;
1925   sprintf ( (S->T[S->nT])->name, "MATCH") ; (S->T[S->nT])->tr=freeT;S->nT++;
1926   sprintf ( (S->T[S->nT])->name, "INSERT"); (S->T[S->nT])->tr=freeT;S->nT++;
1927   sprintf ( (S->T[S->nT])->name, "LINSERT");(S->T[S->nT])->tr=lgop ;S->nT++;
1928   sprintf ( (S->T[S->nT])->name, "END");    (S->T[S->nT])->tr=-gop;S->nT++;
1929
1930   n++;
1931
1932   /*define LInsert*/
1933   S=H->S[n];
1934   sprintf (S->name, "LINSERT"); S->state=n;
1935   S->DI=1;
1936   S->DJ=0;
1937   S->em=lgep;
1938
1939   sprintf ( (S->T[S->nT])->name, "INSERT")  ;(S->T[S->nT])->tr=freeT;S->nT++;
1940   sprintf ( (S->T[S->nT])->name, "LINSERT");(S->T[S->nT])->tr=freeT;S->nT++;
1941
1942   n++;
1943
1944
1945   /*define Delete*/
1946   S=H->S[n];
1947   sprintf (S->name, "DELETE"); S->state=n;
1948   S->DI=0;
1949   S->DJ=1;
1950   S->em=gep;
1951
1952   sprintf ( (S->T[S->nT])->name, "MATCH")   ;(S->T[S->nT])->tr=freeT;S->nT++;
1953    sprintf ( (S->T[S->nT])->name, "DELETE")  ;(S->T[S->nT])->tr=freeT;S->nT++;
1954   sprintf ( (S->T[S->nT])->name, "LDELETE") ;(S->T[S->nT])->tr=lgop ;S->nT++;
1955   sprintf ( (S->T[S->nT])->name, "END")     ;(S->T[S->nT])->tr=-gop;S->nT++;
1956
1957   n++;
1958
1959   /*define LDelete*/
1960   S=H->S[n];
1961   sprintf (S->name, "LDELETE"); S->state=n;
1962   S->DI=0;
1963   S->DJ=1;
1964   S->em=lgep;
1965   sprintf ( (S->T[S->nT])->name, "DELETE")  ;(S->T[S->nT])->tr=freeT;S->nT++;
1966   sprintf ( (S->T[S->nT])->name, "LDELETE");(S->T[S->nT])->tr=freeT;S->nT++;
1967
1968   n++;
1969
1970
1971   H=bound_hmm ( H);
1972   return H;
1973 }
1974
1975 Hmm* define_mnm_model(Constraint_list *CL)
1976 {
1977   Hmm *H;
1978   double gop=20;
1979
1980
1981
1982   double freeT=0;
1983   int n=0;
1984   HmmState *S;
1985
1986
1987   H=declare_hmm(6);
1988   H->freeT=freeT=0;
1989
1990   H->forbiden=FORBIDEN;
1991   H->start=START_STATE;
1992   H->end=END_STATE;
1993
1994   /*define START*/
1995   S=H->S[n];
1996   sprintf (S->name, "START"); S->state=n;
1997
1998   S->DI=0;
1999   S->DJ=0;
2000   S->em=freeT;
2001
2002   sprintf ( (S->T[S->nT])->name, "MATCH") ;(S->T[S->nT])->tr=freeT  ;S->nT++;
2003   sprintf ( (S->T[S->nT])->name, "INSERT");(S->T[S->nT])->tr=freeT  ;S->nT++;
2004   sprintf ( (S->T[S->nT])->name, "DELETE");(S->T[S->nT])->tr=freeT  ;S->nT++;
2005   sprintf ( (S->T[S->nT])->name, "NOMATCH");(S->T[S->nT])->tr=freeT  ;S->nT++;
2006   sprintf ( (S->T[S->nT])->name, "END")   ;(S->T[S->nT])->tr=freeT  ;S->nT++;
2007
2008   n++;
2009   /*define END*/
2010   S=H->S[n];
2011   sprintf (S->name, "END"); S->state=n;
2012   S->DI=0;
2013   S->DJ=0;
2014   S->em=freeT;
2015   n++;
2016
2017   /*define Match*/
2018   S=H->S[n];
2019   sprintf (S->name, "MATCH"); S->state=n;
2020   S->DI=1;
2021   S->DJ=1;
2022   S->em=H->forbiden;
2023   S->em_func=CL->get_dp_cost;
2024
2025   sprintf ( (S->T[S->nT])->name, "MATCH")  ;(S->T[S->nT])->tr=freeT;S->nT++;
2026   sprintf ( (S->T[S->nT])->name, "NOMATCH");(S->T[S->nT])->tr=gop  ;S->nT++;
2027   sprintf ( (S->T[S->nT])->name, "END");   (S->T[S->nT])->tr=freeT;S->nT++;
2028
2029   n++;
2030
2031   /*define NOMatch*/
2032   S=H->S[n];
2033   sprintf (S->name, "NOMATCH"); S->state=n;
2034   S->DI=1;
2035   S->DJ=1;
2036   S->em=freeT;
2037   S->em_func=NULL;
2038
2039   sprintf ( (S->T[S->nT])->name, "NOMATCH")  ;(S->T[S->nT])->tr=freeT;S->nT++;
2040   sprintf ( (S->T[S->nT])->name, "MATCH")  ;(S->T[S->nT])->tr=freeT;S->nT++;
2041   sprintf ( (S->T[S->nT])->name, "INSERT");(S->T[S->nT])->tr=freeT  ;S->nT++;
2042   sprintf ( (S->T[S->nT])->name, "DELETE");(S->T[S->nT])->tr=freeT  ;S->nT++;
2043   sprintf ( (S->T[S->nT])->name, "END");   (S->T[S->nT])->tr=freeT;S->nT++;
2044
2045   n++;
2046   /*define Insert*/
2047   S=H->S[n];
2048   sprintf (S->name, "INSERT"); S->state=n;
2049   S->DI=1;
2050   S->DJ=0;
2051   S->em=freeT;
2052   sprintf ( (S->T[S->nT])->name, "NOMATCH") ; (S->T[S->nT])->tr=freeT;S->nT++;
2053   sprintf ( (S->T[S->nT])->name, "INSERT")  ; (S->T[S->nT])->tr=freeT;S->nT++;
2054   sprintf ( (S->T[S->nT])->name, "END");    (S->T[S->nT])->tr=freeT;S->nT++;
2055
2056   n++;
2057
2058   /*define Delete*/
2059   S=H->S[n];
2060   sprintf (S->name, "DELETE"); S->state=n;
2061   S->DI=0;
2062   S->DJ=1;
2063   S->em=freeT;
2064
2065   sprintf ( (S->T[S->nT])->name, "NOMATCH")   ;(S->T[S->nT])->tr=freeT;S->nT++;
2066   sprintf ( (S->T[S->nT])->name, "DELETE")  ;(S->T[S->nT])->tr=freeT;S->nT++;
2067   sprintf ( (S->T[S->nT])->name, "END")     ;(S->T[S->nT])->tr=freeT;S->nT++;
2068
2069   n++;
2070
2071
2072   H=bound_hmm ( H);
2073   return H;
2074 }
2075
2076 Hmm* define_simple_model(Constraint_list *CL)
2077 {
2078   Hmm *H;
2079   double gop=-10;
2080   double gep=-1;
2081   double freeT=0;
2082   int n=0;
2083   HmmState *S;
2084
2085
2086   H=declare_hmm(5);
2087   H->freeT=freeT=0;
2088
2089   H->forbiden=FORBIDEN;
2090   H->start=START_STATE;
2091   H->end=END_STATE;
2092
2093   /*define START*/
2094   S=H->S[n];
2095   sprintf (S->name, "START"); S->state=n;
2096
2097   S->DI=0;
2098   S->DJ=0;
2099   S->em=freeT;
2100
2101   sprintf ( (S->T[S->nT])->name, "MATCH") ;(S->T[S->nT])->tr=freeT  ;S->nT++;
2102   sprintf ( (S->T[S->nT])->name, "INSERT");(S->T[S->nT])->tr=freeT  ;S->nT++;
2103   sprintf ( (S->T[S->nT])->name, "DELETE");(S->T[S->nT])->tr=freeT  ;S->nT++;
2104   sprintf ( (S->T[S->nT])->name, "END")   ;(S->T[S->nT])->tr=freeT  ;S->nT++;
2105
2106   n++;
2107   /*define END*/
2108   S=H->S[n];
2109   sprintf (S->name, "END"); S->state=n;
2110   S->DI=0;
2111   S->DJ=0;
2112   S->em=freeT;
2113   n++;
2114
2115   /*define Match*/
2116   S=H->S[n];
2117   sprintf (S->name, "MATCH"); S->state=n;
2118   S->DI=1;
2119   S->DJ=1;
2120   S->em=H->forbiden;
2121   S->em_func=CL->get_dp_cost;
2122
2123   sprintf ( (S->T[S->nT])->name, "MATCH") ;(S->T[S->nT])->tr=freeT;S->nT++;
2124   sprintf ( (S->T[S->nT])->name, "INSERT");(S->T[S->nT])->tr=gop  ;S->nT++;
2125   sprintf ( (S->T[S->nT])->name, "DELETE");(S->T[S->nT])->tr=gop  ;S->nT++;
2126   sprintf ( (S->T[S->nT])->name, "END");   (S->T[S->nT])->tr=freeT;S->nT++;
2127
2128   n++;
2129
2130   /*define Insert*/
2131   S=H->S[n];
2132   sprintf (S->name, "INSERT"); S->state=n;
2133   S->DI=1;
2134   S->DJ=0;
2135   S->em=gep;
2136   sprintf ( (S->T[S->nT])->name, "MATCH") ; (S->T[S->nT])->tr=freeT;S->nT++;
2137   sprintf ( (S->T[S->nT])->name, "INSERT"); (S->T[S->nT])->tr=freeT;S->nT++;
2138   sprintf ( (S->T[S->nT])->name, "DELETE"); (S->T[S->nT])->tr=freeT;S->nT++;
2139
2140   sprintf ( (S->T[S->nT])->name, "END");    (S->T[S->nT])->tr=-gop;S->nT++;
2141
2142   n++;
2143
2144
2145   /*define Delete*/
2146   S=H->S[n];
2147   sprintf (S->name, "DELETE"); S->state=n;
2148   S->DI=0;
2149   S->DJ=1;
2150   S->em=gep;
2151
2152   sprintf ( (S->T[S->nT])->name, "MATCH")   ;(S->T[S->nT])->tr=freeT;S->nT++;
2153   sprintf ( (S->T[S->nT])->name, "DELETE")  ;(S->T[S->nT])->tr=freeT;S->nT++;
2154   sprintf ( (S->T[S->nT])->name, "INSERT"); (S->T[S->nT])->tr=freeT;S->nT++;
2155   sprintf ( (S->T[S->nT])->name, "END")     ;(S->T[S->nT])->tr=-gop;S->nT++;
2156
2157   n++;
2158
2159   H=bound_hmm ( H);
2160   return H;
2161 }
2162
2163 Hmm* define_two_mat_model(Constraint_list *CL)
2164 {
2165   Hmm *H;
2166   double gop=-15;
2167   double gep=-2;
2168   double lgop=-6;
2169   double lgep=-1;
2170   double freeT=0;
2171   int n=0;
2172   HmmState *S;
2173
2174
2175   H=declare_hmm(8);
2176   H->freeT=freeT=0;
2177
2178   H->forbiden=FORBIDEN;
2179   H->start=START_STATE;
2180   H->end=END_STATE;
2181
2182   /*define START*/
2183   S=H->S[n];
2184   sprintf (S->name, "START"); S->state=n;
2185
2186   S->DI=0;
2187   S->DJ=0;
2188   S->em=freeT;
2189
2190   sprintf ( (S->T[S->nT])->name, "MATCH1") ;(S->T[S->nT])->tr=freeT  ;S->nT++;
2191   sprintf ( (S->T[S->nT])->name, "MATCH2") ;(S->T[S->nT])->tr=freeT  ;S->nT++;
2192
2193   sprintf ( (S->T[S->nT])->name, "INSERT");(S->T[S->nT])->tr=freeT  ;S->nT++;
2194   sprintf ( (S->T[S->nT])->name, "DELETE");(S->T[S->nT])->tr=freeT  ;S->nT++;
2195   sprintf ( (S->T[S->nT])->name, "END")   ;(S->T[S->nT])->tr=freeT  ;S->nT++;
2196
2197   n++;
2198   /*define END*/
2199   S=H->S[n];
2200   sprintf (S->name, "END"); S->state=n;
2201   S->DI=0;
2202   S->DJ=0;
2203   S->em=freeT;
2204   n++;
2205
2206   /*define Match*/
2207   S=H->S[n];
2208   sprintf (S->name, "MATCH1"); S->state=n;
2209   S->DI=1;
2210   S->DJ=1;
2211   S->em=H->forbiden;
2212   S->em_func=get_dp_cost_pam_matrix;
2213
2214   sprintf ( (S->T[S->nT])->name, "MATCH1") ;(S->T[S->nT])->tr=freeT;S->nT++;
2215   sprintf ( (S->T[S->nT])->name, "INSERT");(S->T[S->nT])->tr=gop  ;S->nT++;
2216   sprintf ( (S->T[S->nT])->name, "DELETE");(S->T[S->nT])->tr=gop  ;S->nT++;
2217   sprintf ( (S->T[S->nT])->name, "END");   (S->T[S->nT])->tr=freeT;S->nT++;
2218
2219   n++;
2220
2221   /*define Match*/
2222   S=H->S[n];
2223   sprintf (S->name, "MATCH2"); S->state=n;
2224   S->DI=1;
2225   S->DJ=1;
2226   S->em=H->forbiden;
2227   S->em_func=get_dp_cost_blosum_matrix;
2228
2229   sprintf ( (S->T[S->nT])->name, "MATCH2") ;(S->T[S->nT])->tr=freeT;S->nT++;
2230   sprintf ( (S->T[S->nT])->name, "INSERT");(S->T[S->nT])->tr=gop  ;S->nT++;
2231   sprintf ( (S->T[S->nT])->name, "DELETE");(S->T[S->nT])->tr=gop  ;S->nT++;
2232   sprintf ( (S->T[S->nT])->name, "END");   (S->T[S->nT])->tr=freeT;S->nT++;
2233
2234   n++;
2235
2236   /*define Insert*/
2237   S=H->S[n];
2238   sprintf (S->name, "INSERT"); S->state=n;
2239   S->DI=1;
2240   S->DJ=0;
2241   S->em=gep;
2242   sprintf ( (S->T[S->nT])->name, "MATCH2") ; (S->T[S->nT])->tr=freeT;S->nT++;
2243   sprintf ( (S->T[S->nT])->name, "MATCH1") ; (S->T[S->nT])->tr=freeT;S->nT++;
2244
2245   sprintf ( (S->T[S->nT])->name, "INSERT"); (S->T[S->nT])->tr=freeT;S->nT++;
2246   sprintf ( (S->T[S->nT])->name, "LINSERT");(S->T[S->nT])->tr=lgop ;S->nT++;
2247   sprintf ( (S->T[S->nT])->name, "END");    (S->T[S->nT])->tr=-gop;S->nT++;
2248
2249   n++;
2250
2251   /*define LInsert*/
2252   S=H->S[n];
2253   sprintf (S->name, "LINSERT"); S->state=n;
2254   S->DI=1;
2255   S->DJ=0;
2256   S->em=lgep;
2257
2258   sprintf ( (S->T[S->nT])->name, "INSERT")  ;(S->T[S->nT])->tr=freeT;S->nT++;
2259   sprintf ( (S->T[S->nT])->name, "LINSERT");(S->T[S->nT])->tr=freeT;S->nT++;
2260
2261   n++;
2262
2263
2264   /*define Delete*/
2265   S=H->S[n];
2266   sprintf (S->name, "DELETE"); S->state=n;
2267   S->DI=0;
2268   S->DJ=1;
2269   S->em=gep;
2270
2271   sprintf ( (S->T[S->nT])->name, "MATCH2")   ;(S->T[S->nT])->tr=freeT;S->nT++;
2272   sprintf ( (S->T[S->nT])->name, "MATCH1")   ;(S->T[S->nT])->tr=freeT;S->nT++;
2273
2274   sprintf ( (S->T[S->nT])->name, "DELETE")  ;(S->T[S->nT])->tr=freeT;S->nT++;
2275   sprintf ( (S->T[S->nT])->name, "LDELETE") ;(S->T[S->nT])->tr=lgop ;S->nT++;
2276   sprintf ( (S->T[S->nT])->name, "END")     ;(S->T[S->nT])->tr=-gop;S->nT++;
2277   n++;
2278
2279   /*define LDelete*/
2280   S=H->S[n];
2281   sprintf (S->name, "LDELETE"); S->state=n;
2282   S->DI=0;
2283   S->DJ=1;
2284   S->em=lgep;
2285   sprintf ( (S->T[S->nT])->name, "DELETE")  ;(S->T[S->nT])->tr=freeT;S->nT++;
2286   sprintf ( (S->T[S->nT])->name, "LDELETE");(S->T[S->nT])->tr=freeT;S->nT++;
2287   n++;
2288
2289   if ( n!=H->nS)
2290     {
2291       fprintf ( stderr, "\nERROR in HMM definition [FATAL:%s]", PROGRAM);
2292       myexit (EXIT_FAILURE);
2293     }
2294
2295   H=bound_hmm ( H);
2296   return H;
2297 }
2298 void DisplayMatState ( MatState *S, char *s)
2299 {
2300   if ( S==NULL)fprintf ( stderr, "%s: Cell is undefined", s);
2301   else  fprintf ( stderr, "%s: i=%d j=%d st=%d pst=%d sc=%d Free %d", s, S->i, S->j, S->st, S->pst, (int)S->sc, S->free);
2302 }
2303 void testfunc ( MatState *S, char *s)
2304 {
2305   if ( S==NULL)return;
2306   fprintf ( stderr, "\n#### %s ", s);
2307   while ( S){DisplayMatState ( S,"\n\t");S=S->n;}
2308   fprintf ( stderr, "\n");
2309 }
2310
2311
2312
2313
2314
2315
2316
2317
2318
2319
2320
2321
2322
2323
2324
2325
2326
2327
2328
2329
2330
2331
2332
2333
2334
2335
2336
2337
2338
2339
2340
2341
2342
2343 #ifdef BACKHERE
2344
2345           if ( i>0 && j>0)
2346           m=emit_pair_default[alphabetDefault[seq1[a]]][alphabetDefault[seq2[a]]];
2347           /*Match*/
2348           F[M][i][j]=F[M][i-1][j-1];
2349
2350
2351
2352         M[Match][i][j]=m+log_add3( M[Match][i-step_i][j-step_j],M[I][i-step_i][j],M[D][i][j-step_j]);
2353         M[D    ][i][j]=log_add3(gep,M[Match][i       ][j-step_j]+gop,M[D][i       ][j-step_j]);
2354         M[I    ][i][j]=log_add3(gep,M[Match][i-step_i][j       ]+gop,M[I][i-step_i][j       ]);
2355
2356         /*Long gaps
2357         M[Match][i][j]=log_add3(M[Match][i][j], M[LI][i-step_i][j],M[LD][i][j-step_j]);
2358         M[LI    ][i][j]=log_add3(lgep, M[I][i-step_i][j       ]+lgop,M[LI][i-step_i][j      ]);
2359         M[LD    ][i][j]=log_add3(lgep, M[D][i       ][j-step_j]+lgop,M[LD][i       ][j-step_j]);
2360         */
2361
2362         }
2363     }
2364   retun M;
2365 MatState* RviterbiL_hmm (Alignment *A,int *ns, int **ls, Hmm *H, Constraint_list *CL,MatState *S, MatState *E)
2366 {
2367   MatState *Mid;
2368   Mid=viterbiL_hmm (A,ns, ls,H, CL, S, E);
2369
2370   if (!Mid)
2371     {
2372       return S;
2373     }
2374   else if ( Mid->n)
2375     {
2376       return Mid;
2377     }
2378
2379   else
2380     {
2381       Mid->p=S;S->n=Mid;
2382       Mid->n=E;E->p=Mid;
2383       RviterbiL_hmm (A,ns, ls,H, CL,S,   Mid);
2384       RviterbiL_hmm (A,ns, ls,H, CL,Mid, E);
2385       return S;
2386     }
2387 }
2388
2389 MatState* viterbiL_hmm (Alignment *A,int *ns, int **ls, Hmm *H, Constraint_list *CL, MatState *S,MatState *E)
2390 {
2391   int current,memory, dim;
2392   double e, v,t;
2393   int i,j,pi,pj, s, k;
2394   int start_i, start_j, end_i, end_j, l1, l2;
2395   HmmState *S1, *S2;
2396   static MatState ****M;
2397   static int maxl;
2398   MatState *Mid=NULL;
2399
2400   MatState *CC, *PCC;
2401   int midpoint;
2402   int Delta;
2403
2404
2405   if ( MatStateAreIdentical ( S, E))return NULL;
2406
2407   l1=strlen (A->seq_al[ls[0][0]]);
2408   l2=strlen (A->seq_al[ls[1][0]]);
2409
2410   midpoint=S->i+(E->i-S->i)/2;
2411   Delta=E->i-S->i;
2412
2413   start_i=S->i;end_i=E->i;
2414   start_j=S->j;end_j=E->j;
2415
2416   dim=H->order+2;current=0;memory=H->order+1;
2417   if (!M || (l2+1)>maxl)
2418     { free_arrayN((void **)M, 4);
2419       M=declare_arrayN(4, sizeof ( MatState), dim, maxl=(l2+1), H->nS,1);
2420     }
2421
2422   /*MAKE THE VITERBI FROM S(tart) to E(nd)*/
2423   for ( i=start_i; i<=end_i; i++)
2424     {
2425       M= (MatState****)recycle  ( (void **)M,H->order+1,1);
2426       for ( j=start_j; j<=end_j; j++)
2427         {
2428           for ( s=H->nS-1;s>=0; s--)
2429             {
2430
2431               S1=H->S[s];
2432               pi=i-S1->DI;pj=j-S1->DJ;
2433               CC=M[current][j][s];
2434               CC->i=i; CC->j=j; CC->st=s;CC->sc=H->forbiden;CC->p=CC->n=CC->m=NULL;CC->sc=H->forbiden;
2435               if (i==start_i && j==start_j && s==S->st)  {CopyMatState(S,CC);}
2436               else if ( i==end_i && j==end_j && s==E->st && s!=H->end)
2437                 {
2438                   S2=H->S[E->pst];
2439                   CopyMatState(E,CC);
2440                   CC->p=M[S1->DI][j-S1->DJ][S2->state];
2441                 }
2442               else if ( pi<start_i || pj<start_j)        {CC->sc=H->forbiden;}
2443               else
2444                 {
2445                   for (k=1; k<=H->fromM[S1->state][0]; k++)
2446                     {
2447                       S2=H->S[H->fromM[s][k]];
2448                       PCC=M[S1->DI][j-S1->DJ][S2->state];
2449
2450                       if      ( pi+pj!=0 && S2->state==H->start) {t=H->forbiden;}
2451                       else if ( !(pi==l1 && pj==l2) && s==H->end){t=H->forbiden;}
2452                       else   t=H->T[S2->state][S1->state];
2453
2454                       v=hmm_add(t,PCC->sc);
2455                       if ( v!=H->forbiden && (CC->sc==H->forbiden || v> CC->sc)){CC->sc=v; CC->pst=S2->state;CC->p=PCC;}
2456                     }
2457
2458                   e=(S1->em==H->forbiden)?S1->em_func (A, A->pos, ns[0], ls[0],i-1, A->pos,ns[1], ls[1], j-1, CL):S1->em;
2459                   CC->sc=hmm_add(CC->sc,e);
2460                 }
2461
2462               if (i==midpoint)CC->m=CopyMatState(CC, M[memory][j][s]);
2463               else if (i>midpoint && CC->sc!=H->forbiden) CC->m=(M[S1->DI][j-S1->DJ][CC->pst])->m;
2464             }
2465         }
2466     }
2467
2468   if ( E->st==H->end)CopyMatState ((M[current][end_j][E->st]),E);
2469
2470   if ( Delta>1)
2471     {
2472       Mid=CopyMatState ((M[current][end_j][E->st])->m,NULL);
2473     }
2474   else if ( Delta==1)
2475     {
2476       CC=M[current][E->j][E->st];
2477       Mid=E;
2478       while (!MatStateAreIdentical (CC->p, S) )
2479         {
2480           Mid->p=CopyMatState(CC->p,NULL);
2481           (Mid->p)->n=Mid;
2482           Mid=Mid->p;CC=CC->p;
2483         }
2484       Mid->p=S;
2485       S->n=Mid;
2486       Mid=S;
2487     }
2488
2489   return Mid;
2490 }
2491 #endif
2492 /*********************************COPYRIGHT NOTICE**********************************/
2493 /*© Centre National de la Recherche Scientifique (CNRS) */
2494 /*and */
2495 /*Please Cite: Notredame*/
2496 /*Mon May 17 20:15:35 MDT 2004. */
2497 /*All rights reserved.*/
2498 /*NOTICE:                                                                                                                                     |*/
2499 /*  This file is an integral part of the */
2500 /*  ALIGN_TWO_SEQ Software. */
2501 /*  Its content is protected and all */
2502 /*  the conditions mentioned in the licensing */
2503 /*  agreement of the software apply to this file.*/
2504 /*...............................................                                                                                      |*/
2505 /*  If you need some more information, or if you */
2506 /*  wish to obtain a full license, please contact: */
2507 /*  cedric.notredame@europe.com*/
2508 /*...............................................                                                                                                                                     |*/
2509 /**/
2510 /**/
2511 /*      */
2512 /*********************************COPYRIGHT NOTICE**********************************/
2513 /*********************************COPYRIGHT NOTICE**********************************/
2514 /*© Centro de Regulacio Genomica */
2515 /*and */
2516 /*Cedric Notredame */
2517 /*Tue Oct 27 10:12:26 WEST 2009. */
2518 /*All rights reserved.*/
2519 /*This file is part of T-COFFEE.*/
2520 /**/
2521 /*    T-COFFEE is free software; you can redistribute it and/or modify*/
2522 /*    it under the terms of the GNU General Public License as published by*/
2523 /*    the Free Software Foundation; either version 2 of the License, or*/
2524 /*    (at your option) any later version.*/
2525 /**/
2526 /*    T-COFFEE is distributed in the hope that it will be useful,*/
2527 /*    but WITHOUT ANY WARRANTY; without even the implied warranty of*/
2528 /*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the*/
2529 /*    GNU General Public License for more details.*/
2530 /**/
2531 /*    You should have received a copy of the GNU General Public License*/
2532 /*    along with Foobar; if not, write to the Free Software*/
2533 /*    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA*/
2534 /*...............................................                                                                                      |*/
2535 /*  If you need some more information*/
2536 /*  cedric.notredame@europe.com*/
2537 /*...............................................                                                                                                                                     |*/
2538 /**/
2539 /**/
2540 /*      */
2541 /*********************************COPYRIGHT NOTICE**********************************/