binaries/src/tcoffee/t_coffee_source/util_dp_suboptimal_nw.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <math.h>
   4 #include <stdarg.h>
   5 #include <string.h>
   6 #include <ctype.h>
   7 #include "io_lib_header.h"
   8 #include "util_lib_header.h"
   9 #include "define_header.h"
  10
  11
  12
  13
  14 //Values as provided in Probcons V1.1
  15 static float EXP_UNDERFLOW_THRESHOLD = -4.60f;
  16 static float LOG_UNDERFLOW_THRESHOLD = 7.50f;
  17 //static float LOG_ZERO = -FLT_MAX;
  18 static float LOG_ZERO=-200000004008175468544.000000;
  19 static float LOG_ONE = 0.0f;
  20
  21 //Probabilistic Alignment DNA
  22
  23
  24
  25
  26 static float DNAinitDistrib2Default[] = { 0.9615409374f, 0.0000004538f, 0.0000004538f, 0.0192291681f, 0.0192291681f };
  27 static float DNAgapOpen2Default[] = { 0.0082473317f, 0.0082473317f, 0.0107844425f, 0.0107844425f };
  28 static float DNAgapExtend2Default[] = { 0.3210460842f, 0.3210460842f, 0.3298229277f, 0.3298229277f };
  29
  30 static char DNAalphabetDefault[] = "ACGUTN";
  31 static float DNAemitSingleDefault[6] = {0.2270790040f, 0.2422080040f, 0.2839320004f, 0.2464679927f, 0.2464679927f, 0.0003124650f};
  32
  33 static float DNAemitPairsDefault[6][6] = {
  34   { 0.1487240046f, 0.0184142999f, 0.0361397006f, 0.0238473993f, 0.0238473993f, 0.0000375308f },
  35   { 0.0184142999f, 0.1583919972f, 0.0275536999f, 0.0389291011f, 0.0389291011f, 0.0000815823f },
  36   { 0.0361397006f, 0.0275536999f, 0.1979320049f, 0.0244289003f, 0.0244289003f, 0.0000824765f },
  37   { 0.0238473993f, 0.0389291011f, 0.0244289003f, 0.1557479948f, 0.1557479948f, 0.0000743985f },
  38   { 0.0238473993f, 0.0389291011f, 0.0244289003f, 0.1557479948f, 0.1557479948f, 0.0000743985f },
  39   { 0.0000375308f, 0.0000815823f, 0.0000824765f, 0.0000743985f, 0.0000743985f, 0.0000263252f }
  40 };
  41
  42 //Probabilistic ALignment Blosum62mt
  43
  44
  45
  46
  47 static float initDistrib2Default[] = { 0.6814756989f, 8.615339902e-05f, 8.615339902e-05f, 0.1591759622f, 0.1591759622 };
  48 static float gapOpen2Default[] = { 0.0119511066f, 0.0119511066f, 0.008008334786f, 0.008008334786 };
  49 static float gapExtend2Default[] = { 0.3965826333f, 0.3965826333f, 0.8988758326f, 0.8988758326 };
  50
  51 static char alphabetDefault[] = "ARNDCQEGHILKMFPSTWYV";
  52 static float emitSingleDefault[20] = {
  53   0.07831005f, 0.05246024f, 0.04433257f, 0.05130349f, 0.02189704f,
  54   0.03585766f, 0.05615771f, 0.07783433f, 0.02601093f, 0.06511648f,
  55   0.09716489f, 0.05877077f, 0.02438117f, 0.04463228f, 0.03940142f,
  56   0.05849916f, 0.05115306f, 0.01203523f, 0.03124726f, 0.07343426f
  57 };
  58
  59 static float emitPairsDefault[20][20] = {
  60   {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},
  61   {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},
  62   {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},
  63   {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},
  64   {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},
  65   {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},
  66   {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},
  67   {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},
  68   {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},
  69   {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},
  70   {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},
  71   {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},
  72   {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},
  73   {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},
  74   {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},
  75   {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},
  76   {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},
  77   {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},
  78   {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},
  79   {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}
  80 };
  81
  82
  83 static int suboptimal_pair_wise ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int mode);
  84 static int *** forward_so_dp ( Alignment *A, int *ns, int **ls, int **pos,int I, int J, int gop, int gep,int gop2, int gep2,Constraint_list *CL);
  85 static int *** backward_so_dp ( Alignment *A, int *ns, int **ls,int **pos,int I, int J, int gop, int gep,int gop2, int gep2,Constraint_list *CL);
  86 static int *** forward_so_dp_biphasic ( Alignment *A, int *ns, int **ls, int **pos,int I, int J, int gop, int gep,int gop2, int gep2,Constraint_list *CL);
  87 static int *** backward_so_dp_biphasic ( Alignment *A, int *ns, int **ls,int **pos,int I, int J, int gop, int gep,int gop2, int gep2,Constraint_list *CL);
  88 static int *** forward_so_dp_glocal ( Alignment *A, int *ns, int **ls, int **pos,int I, int J, int gop, int gep,int gop2, int gep2,Constraint_list *CL);
  89 static int *** backward_so_dp_glocal ( Alignment *A, int *ns, int **ls,int **pos,int I, int J, int gop, int gep,int gop2, int gep2,Constraint_list *CL);
  90
  91 static int match=0;
  92 static int ins=1;
  93 static int del=2;
  94 static int umatch=3;
  95 static int ins2=3;
  96 static int del2=4;
  97 float **    get_emitPairs (char *mat, char *alp, float **p, float *s);
  98 int subop1_pair_wise ( Alignment *A, int *ns, int **ls, Constraint_list *CL)
  99 {
 100   return suboptimal_pair_wise ( A, ns, ls, CL, 1);
 101 }
 102
 103 int subop2_pair_wise ( Alignment *A, int *ns, int **ls, Constraint_list *CL)
 104 {
 105   return suboptimal_pair_wise ( A, ns, ls, CL, 3);
 106 }
 107
 108
 109
 110 int suboptimal_pair_wise ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int mode)
 111 {
 112   int ***F=NULL;
 113   int ***B=NULL;
 114   int **pos0;
 115   int gop, gep,gop2, gep2;
 116   int i, I, j, J, n, s1, s2;
 117   char *seqI, *seqJ;
 118   int id;
 119   int *entry;
 120   float opt, min, score, nscore, thres;
 121   int l1, l2, set;
 122
 123
 124   gop=CL->gop*SCORE_K;
 125   gep=CL->gep*SCORE_K;
 126
 127   /*gop2=CL->gop*10*SCORE_K;*/
 128   gop2=CL->gop*2*SCORE_K;
 129   gep2=0;
 130
 131   //Values Adapted from Probcons 1.1
 132   gop=-132;
 133   gep=-27;
 134
 135   gop2=-144;
 136   gep2=-3;
 137
 138   ungap(A->seq_al[ls[0][0]]);
 139   ungap(A->seq_al[ls[1][0]]);
 140
 141   seqI=A->seq_al[ls[0][0]];
 142   seqJ=A->seq_al[ls[1][0]];
 143
 144   I=strlen (seqI); J=strlen (seqJ);
 145   pos0=aln2pos_simple ( A,-1, ns, ls);
 146   l1=strlen (A->seq_al[ls[0][0]]);
 147   l2=strlen (A->seq_al[ls[1][0]]);
 148
 149   if ( mode==1)
 150     {
 151       F=forward_so_dp (A, ns, ls, pos0,I, J,gop, gep,gop2, gep2,CL);
 152       B=backward_so_dp (A, ns, ls, pos0,I, J,gop, gep,gop2, gep2, CL);
 153     }
 154   else if ( mode ==2)
 155     {
 156       F=forward_so_dp_glocal (A, ns, ls, pos0,I, J,gop, gep,gop2, gep2,CL);
 157       B=backward_so_dp_glocal (A, ns, ls, pos0,I, J,gop, gep,gop2, gep2, CL);
 158     }
 159   else if ( mode ==3)
 160     {
 161       F=forward_so_dp_biphasic (A, ns, ls, pos0,I, J,gop, gep,gop2, gep2,CL);
 162       B=backward_so_dp_biphasic (A, ns, ls, pos0,I, J,gop, gep,gop2, gep2, CL);
 163     }
 164   if ( MAX5(F[match][l1][l2], F[ins][l1][l2], F[del][l1][l2],F[ins2][l1][l2], F[del2][l1][l2] )!=MAX5( B[match][1][1], B[ins][1][1], B[del][1][1], B[ins2][1][1], B[del2][1][1]))
 165     {
 166       HERE ("ERROR in subop_pair");
 167       fprintf ( stdout, "\nForward:  %d", MAX3(F[match][l1][l2], F[ins][l1][l2], F[del][l1][l2]));
 168       fprintf ( stdout, "\nBackWard: %d \n\n",MAX3( B[match][1][1], B[ins][1][1], B[del][1][1]));
 169     }
 170
 171
 172   for (opt=0,min=0, set=0, i=1; i<=I; i++)
 173     for (j=1; j<=J; j++)
 174       {
 175         if ( F[match][i][j]==UNDEFINED)continue;
 176         F[match][i][j]+=B[match][i][j]-(CL->get_dp_cost) (A, pos0, ns[0], ls[0], i-1, pos0, ns[1], ls[1],j-1,CL);
 177         if (set==0)
 178           {set=1; opt=F[match][i][j];min=F[match][i][j];}
 179         opt=MAX(F[match][i][j],opt);
 180         min=MIN(F[match][i][j],min);
 181       }
 182
 183
 184   s1=name_is_in_list (A->name[ls[0][0]], (CL->S)->name, (CL->S)->nseq, 100);
 185   s2=name_is_in_list (A->name[ls[1][0]], (CL->S)->name, (CL->S)->nseq, 100);
 186
 187   id=idscore_pairseq(seqI,seqJ,-12, -1, CL->M, "idmat");
 188
 189   entry=vcalloc ( CL->entry_len, CL->el_size);
 190   entry[SEQ1]=s1;entry[SEQ2]=s2;
 191
 192   thres=opt;
 193   for ( n=0,i=1; i<=I; i++)
 194     {
 195       for (j=1; j<=J; j++)
 196         {
 197           score=F[0][i][j];
 198           nscore=((score-min))/(opt-min);
 199
 200           if (score==opt)
 201             {
 202               n++;
 203               entry[R1]=i;entry[R2]=j;
 204               entry[WE]=id;
 205               entry[CONS]=1;
 206               add_entry2list (entry,A->CL);
 207             }
 208         }
 209     }
 210
 211   vfree (entry);
 212   free_int (pos0, -1);
 213   free_arrayN (F, 3);
 214   free_arrayN (B, 3);
 215
 216   return A->score_aln;
 217 }
 218 /************************************************************************************************************************/
 219 /*                                                                                                                      */
 220 /*                                                                                                                      */
 221 /*                                                     GLOCAL                                                           */
 222 /*                                                                                                                      */
 223 /*                                                                                                                      */
 224 /************************************************************************************************************************/
 225 int *** forward_so_dp_glocal ( Alignment *A, int *ns, int **ls, int **pos0,int I, int J,int gop, int gep,int gop2, int gep2,Constraint_list *CL)
 226 {
 227   int i,j;
 228   int c;
 229   int sub;
 230   int ***M;
 231   int match=0, del=1, ins=2;
 232
 233   M=declare_arrayN (3, sizeof (int), 5, I+1, J+1);
 234
 235   for ( i=0; i<=I; i++)for (j=0; j<=J; j++)for (c=0; c<5; c++)M[c][i][j]=-999999;
 236
 237   M[match][0][0]=0;
 238
 239   for (i=1; i<=I; i++){M[del]  [i][0]=i*gep;M[umatch][i][0]=i*gep2+gop2;}
 240   for (j=1; j<=J; j++){M[ins]  [0][j]=j*gep;M[umatch][0][j]=j*gep2+gop2;}
 241
 242
 243   for (i=1; i<=I; i++)
 244     {
 245       for ( j=1; j<=J; j++)
 246         {
 247         sub=(CL->get_dp_cost) (A, pos0, ns[0], ls[0], i-1, pos0, ns[1], ls[1],j-1,CL);
 248
 249         M[match][i][j] =MAX4  (M[match][i-1][j-1],M[del][i-1][j-1], M[ins][i-1][j-1],M[umatch][i-1][j-1])+sub;
 250         M[del][i][j]   =MAX2       ((M[match][i-1][j]+gop), M[del][i-1][j])+gep;
 251         M[ins][i][j]   =MAX2       ((M[match][i][j-1]+gop), M[ins][i][j-1])+gep;
 252         M[umatch][i][j]=MAX6 (M[match][i-1][j-1]+gop2, M[match][i][j-1]+gop2, M[match][i-1][j]+gop2,M[umatch][i-1][j-1], M[umatch][i-1][j], M[umatch][i][j-1])+gep2;
 253         }
 254     }
 255   return M;
 256 }
 257 int *** backward_so_dp_glocal ( Alignment *A, int *ns, int **ls, int **pos0, int I, int J, int gop, int gep,int gop2, int gep2,Constraint_list *CL)
 258 {
 259   int i,j;
 260   int c;
 261   int sub;
 262   int ***M;
 263
 264
 265
 266   M=declare_arrayN (3, sizeof (int), 5, I+2, J+2);
 267   for ( i=I+1; i>=0; i--)for (j=J+1; j>=0; j--)for (c=0; c<5; c++)M[c][i][j]=-999999;
 268   M[match][I+1][J+1]=0;
 269
 270   for (i=I; i>0; i--){M[ins]  [i][J+1]=i*gep;M[umatch]  [i][J+1]=i*gep2+gop2;}
 271   for (j=J; j>0; j--){M[del]  [I+1][j]=j*gep;M[umatch]  [I+1][j]=j*gep2+gop2;}
 272
 273   for (i=I; i>0; i--)
 274     {
 275       for ( j=J; j>0; j--)
 276         {
 277         sub=(CL->get_dp_cost) (A, pos0, ns[0], ls[0], i-1, pos0, ns[1], ls[1],j-1,CL);
 278
 279         M[match ][i][j]  =MAX4 ((M[del][i+1][j+1]+gop), (M[ins][i+1][j+1]+gop), M[match][i+1][j+1], M[umatch][i+1][j+1]+gop2)+sub;
 280         M[del   ][i][j]  =MAX2 (M[match][i+1][j], M[del][i+1][j])+gep;
 281         M[ins   ][i][j]  =MAX2 (M[match][i][j+1], M[ins][i][j+1])+gep;
 282         M[umatch][i][j]  =MAX6 (M[match][i+1][j+1], M[match][i+1][j],M[match][i][j+1], M[umatch][i+1][j+1], M[umatch][i+1][j], M[umatch][i][j+1])+gep2;
 283
 284         }
 285     }
 286   return M;
 287 }
 288
 289
 290
 291
 292 /************************************************************************************************************************/
 293 /*                                                                                                                      */
 294 /*                                                                                                                      */
 295 /*                                                     SIMPLE                                                           */
 296 /*                                                                                                                      */
 297 /*                                                                                                                      */
 298 /************************************************************************************************************************/
 299
 300 int *** forward_so_dp ( Alignment *A, int *ns, int **ls, int **pos0,int I, int J,int gop, int gep,int gop2, int gep2,Constraint_list *CL)
 301 {
 302   int i,j;
 303   int c;
 304   int sub;
 305   int ***M;
 306   int lgop;
 307
 308
 309
 310   M=declare_arrayN (3, sizeof (int), 5, I+1, J+1);
 311   for ( i=0; i<=I; i++)for (j=0; j<=J; j++)for (c=0; c<3; c++)M[c][i][j]=-999999;
 312
 313   M[match][0][0]=0;
 314   for (i=1; i<=I; i++){M[del]  [i][0]=i*gep;}
 315   for (j=1; j<=J; j++){M[ins]  [0][j]=j*gep;}
 316
 317
 318
 319   for (i=1; i<=I; i++)
 320     {
 321       for ( j=1; j<=J; j++)
 322         {
 323           lgop=(i==I || j==J)?0:gop;
 324           sub=(CL->get_dp_cost) (A, pos0, ns[0], ls[0], i-1, pos0, ns[1], ls[1],j-1,CL);
 325
 326           M[match][i][j]=MAX3  (M[del][i-1][j-1], M[ins][i-1][j-1], M[match][i-1][j-1])+sub;
 327           M[del][i][j]  =MAX ((M[match][i-1][j]+lgop),M[del][i-1][j])+gep;
 328           M[ins][i][j]  =MAX ((M[match][i][j-1]+lgop),  M[ins][i][j-1])+gep;
 329         }
 330
 331     }
 332
 333   return M;
 334   }
 335 int *** backward_so_dp ( Alignment *A, int *ns, int **ls, int **pos0, int I, int J, int gop, int gep,int gop2, int gep2,Constraint_list *CL)
 336 {
 337   int i,j, a, b;
 338
 339
 340   int ***M, ***T;
 341
 342
 343   for (a=0; a<2; a++)
 344     for (b=0; b<ns[a]; b++)
 345       {
 346         invert_string2(A->seq_al[ls[a][b]]);
 347         invert_string2((CL->S)->seq[A->order[ls[a][b]][0]]);
 348       }
 349   T=forward_so_dp(A,ns,ls,pos0, I, J, gop, gep, gop2, gep2, CL);
 350   for (a=0; a<2; a++)
 351     for (b=0; b<ns[a]; b++)
 352       {
 353         invert_string2(A->seq_al[ls[a][b]]);
 354         invert_string2((CL->S)->seq[A->order[ls[a][b]][0]]);
 355       }
 356
 357   M=declare_arrayN (3, sizeof (int), 5, I+2, J+2);
 358
 359
 360   for (i=0; i<=I; i++)
 361     for (j=0; j<=J; j++)
 362       {
 363         M[match][i+1][j+1]=T[match][I-i][J-j];
 364         M[ins][i+1][j+1]=T[ins][I-i][J-j];
 365         M[del][i+1][j+1]=T[del][I-i][J-j];
 366       }
 367   return M;
 368 }
 369
 370 /************************************************************************************************************************/
 371 /*                                                                                                                      */
 372 /*                                                                                                                      */
 373 /*                                                     BI-PHASIC                                                        */
 374 /*                                                                                                                      */
 375 /*                                                                                                                      */
 376 /************************************************************************************************************************/
 377 int biphasic_pair_wise ( Alignment *A, int *ns, int **ls, Constraint_list *CL)
 378 {
 379   int i,j,a,b;
 380   int c;
 381   int sub;
 382   int ***m, ***t;
 383   int M1, D1, D2, I1, I2, LEN;
 384   int I, J;
 385   int n=1;
 386   char **al, **aln, *char_buf;
 387   int gop1, gop2, gep1, gep2;
 388   int **pos0;
 389   int score, trace, ntrace;
 390   M1=n++; D1=n++; D2=n++; I1=n++, I2=n++;
 391
 392   I=strlen (A->seq_al[ls[0][0]]);
 393   J=strlen (A->seq_al[ls[1][0]]);
 394   m=declare_arrayN (3, sizeof (int),n, I+1, J+1);
 395   t=declare_arrayN (3, sizeof (int),n, I+1, J+1);
 396   pos0=aln2pos_simple ( A,-1, ns, ls);
 397   al=declare_char (2, I+J+1);
 398   for ( i=0; i<=I; i++)for (j=0; j<=J; j++)for (c=0; c<n; c++)m[c][i][j]=-999999;
 399
 400   gop1=CL->gop*SCORE_K*2;
 401   gep1=CL->gep*SCORE_K/2;
 402
 403   gop2=CL->gop*SCORE_K/2;
 404   gep2=CL->gep*SCORE_K*2;
 405
 406   m[M1][0][0]=0;
 407   for (i=1; i<=I; i++){m[I1][i][0]=gep1*i;}
 408   for (j=1; j<=J; j++){m[D1][0][j]=gep1*j;}
 409
 410   for (i=1; i<=I; i++){m[I2]  [i][0]=gep2*i;}
 411   for (j=1; j<=J; j++){m[D2]  [0][j]=gep2*j;}
 412
 413   for (i=1; i<=I; i++)
 414     {
 415       for ( j=1; j<=J; j++)
 416         {
 417           sub=(CL->get_dp_cost) (A, pos0, ns[0], ls[0], i-1, pos0, ns[1], ls[1],j-1,CL);
 418           m[M1][i][j]=max_int  (&t[M1][i][j],D1,m[D1][i-1][j-1],I1,m[I1][i-1][j-1], M1, m[M1][i-1][j-1],D2,m[D2][i-1][j-1],I2,m[I2][i-1][j-1], -1)+sub;
 419
 420           m[D1][i][j]=max_int  (&t[D1][i][j],M1,(m[M1][i][j-1]+gop1),D1,m[D1][i][j-1], -1)+gep1;
 421           m[I1][i][j]=max_int  (&t[I1][i][j],M1,(m[M1][i-1][j]+gop1),I1,m[I1][i-1][j], -1)+gep1;
 422
 423           m[D2][i][j]=max_int  (&t[D2][i][j],M1,(m[M1][i][j-1]+gop2),D2,m[D2][i][j-1], -1)+gep2;
 424           m[I2][i][j]=max_int  (&t[I2][i][j],M1,(m[M1][i-1][j]+gop2),I2,m[I2][i-1][j], -1)+gep2;
 425         }
 426     }
 427
 428   score=max_int (&trace,M1,m[M1][I][J],D1,m[D1][I][J],I1, m[I1][I][J],D2,m[D2][I][J],I2,m[I2][I][J], -1);
 429   LEN=0;i=I;j=J;
 430
 431
 432   trace=t[trace][i][j];
 433   while (!(i==0 &&j==0))
 434     {
 435
 436       ntrace=t[trace][i][j];
 437       if (i==0)
 438         {
 439           al[0][LEN]=0;
 440           al[1][LEN]=1;
 441           j--;
 442           LEN++;
 443         }
 444       else if ( j==0)
 445         {
 446           al[0][LEN]=1;
 447           al[1][LEN]=0;
 448           i--;
 449           LEN++;
 450         }
 451       else if ( trace==M1)
 452         {
 453           al[0][LEN]=1;
 454           al[1][LEN]=1;
 455           i--; j--;
 456           LEN++;
 457         }
 458
 459       else if ( trace==D1 || trace==D2)
 460         {
 461           al[0][LEN]=0;
 462           al[1][LEN]=1;
 463           j--;
 464           LEN++;
 465         }
 466       else if ( trace == I1 || trace==I2)
 467         {
 468           al[0][LEN]=1;
 469           al[1][LEN]=0;
 470           i--;
 471           LEN++;
 472         }
 473       trace=ntrace;
 474
 475     }
 476
 477   invert_list_char ( al[0], LEN);
 478   invert_list_char ( al[1], LEN);
 479   if ( A->declared_len<=LEN)A=realloc_aln2  ( A,A->max_n_seq, 2*LEN);
 480
 481   aln=A->seq_al;
 482   char_buf= vcalloc (LEN+1, sizeof (char));
 483   for ( c=0; c< 2; c++)
 484     {
 485       for ( a=0; a< ns[c]; a++)
 486         {
 487           int ch=0;
 488           for ( b=0; b< LEN; b++)
 489             {
 490               if (al[c][b]==1)
 491                 char_buf[b]=aln[ls[c][a]][ch++];
 492               else
 493                 char_buf[b]='-';
 494             }
 495           char_buf[b]='\0';
 496           sprintf (aln[ls[c][a]],"%s", char_buf);
 497         }
 498     }
 499
 500
 501   A->len_aln=LEN;
 502   A->nseq=ns[0]+ns[1];
 503   free_arrayN((void *)m, 3);
 504   free_arrayN((void *)t, 3);
 505   vfree (char_buf);
 506   free_char (al, -1);
 507   return score;
 508   }
 509 int *** forward_so_dp_biphasic ( Alignment *A, int *ns, int **ls, int **pos0,int I, int J,int gop1, int gep1,int gop2, int gep2,Constraint_list *CL)
 510 {
 511   int i,j;
 512   int c;
 513   int sub;
 514   int ***M;
 515   int match=0, del=1, ins=2;
 516   int lgop1, lgop2, lgep1, lgep2;
 517
 518   M=declare_arrayN (3, sizeof (int), 5, I+1, J+1);
 519
 520   for ( i=0; i<=I; i++)for (j=0; j<=J; j++)for (c=0; c<5; c++)M[c][i][j]=-999999;
 521
 522   M[match][0][0]=0;
 523
 524   for (i=1; i<=I; i++){M[del]  [i][0]=gep1*i+gop1;}
 525   for (j=1; j<=J; j++){M[ins]  [0][j]=gep1*j+gop1;}
 526
 527   for (i=1; i<=I; i++){M[del2]  [i][0]=gep2*i+gop2;}
 528   for (j=1; j<=J; j++){M[ins2]  [0][j]=gep2*j+gop2;}
 529
 530   for (i=1; i<=I; i++)
 531     {
 532       for ( j=1; j<=J; j++)
 533         {
 534           lgop1=(i==I || j==J)?gop1:gop1;
 535           lgop2=(i==I || j==J)?gop2:gop2;
 536           lgep1=gep1;
 537           lgep2=gep2;
 538
 539           sub=(CL->get_dp_cost) (A, pos0, ns[0], ls[0], i-1, pos0, ns[1], ls[1],j-1,CL);
 540           M[match][i][j]=MAX5  (M[del][i-1][j-1], M[ins][i-1][j-1], M[match][i-1][j-1], M[ins2][i-1][j-1], M[del2][i-1][j-1])+sub;
 541
 542           M[del ][i][j] =MAX2 ((M[match][i-1][j]+lgop1), M[del ][i-1][j])+lgep1;
 543           M[del2][i][j] =MAX2 ((M[match][i-1][j]+lgop2), M[del2][i-1][j])+lgep2;
 544
 545           M[ins ][i][j] =MAX2 ((M[match][i][j-1]+lgop1), M[ins ][i][j-1] )+lgep1;
 546           M[ins2][i][j] =MAX2 ((M[match][i][j-1]+lgop2), M[ins2][i][j-1] )+lgep2;
 547         }
 548     }
 549   return M;
 550   }
 551 int *** backward_so_dp_biphasic ( Alignment *A, int *ns, int **ls, int **pos0, int I, int J, int gop, int gep,int gop2, int gep2,Constraint_list *CL)
 552 {
 553   int i,j, a, b;
 554
 555
 556   int ***M, ***T;
 557
 558
 559   for (a=0; a<2; a++)
 560     for (b=0; b<ns[a]; b++)
 561       {
 562         invert_string2(A->seq_al[ls[a][b]]);
 563         invert_string2((CL->S)->seq[A->order[ls[a][b]][0]]);
 564       }
 565   T=forward_so_dp_biphasic(A,ns,ls,pos0, I, J, gop, gep, gop2, gep2, CL);
 566   for (a=0; a<2; a++)
 567     for (b=0; b<ns[a]; b++)
 568       {
 569         invert_string2(A->seq_al[ls[a][b]]);
 570         invert_string2((CL->S)->seq[A->order[ls[a][b]][0]]);
 571       }
 572
 573   M=declare_arrayN (3, sizeof (int), 5, I+2, J+2);
 574
 575
 576   for (i=0; i<=I; i++)
 577     for (j=0; j<=J; j++)
 578       {
 579         M[match][i+1][j+1]=T[match][I-i][J-j];
 580         M[ins][i+1][j+1]=T[ins][I-i][J-j];
 581         M[del][i+1][j+1]=T[del][I-i][J-j];
 582         M[ins2][i+1][j+1]=T[ins2][I-i][J-j];
 583         M[del2][i+1][j+1]=T[del2][I-i][J-j];
 584       }
 585   free_arrayN(T,3);
 586   return M;
 587 }
 588
 589
 590 int get_tot_prob (Alignment *A1,Alignment *A2, int *ns, int **ls, int nstates, float **matchProb, float **insProb, float *TmatchProb, float ***TinsProb, Constraint_list *CL);
 591
 592 float * forward_proba_pair_wise  ( char *seq1, char *seq2, int NumMatrixTypes, int NumInsertStates, float **transMat, float *initialDistribution,float *TmatchProb, float ***TinsProb, float **transProb);
 593 float * backward_proba_pair_wise ( char *seq1, char *seq2, int NumMatrixTypes, int NumInsertStates, float **transMat, float *initialDistribution,float *TmatchProb, float ***TinsProb,float **transProb);
 594 float ComputeTotalProbability (int seq1Length, int seq2Length,int NumMatrixTypes, int NumInsertStates,float *forward, float *backward) ;
 595 int ProbabilisticModel (int NumMatrixTypes, int NumInsertStates,float *initDistribMat,float *emitSingle,  float** emitPairs, float *gapOpen, float *gapExtend, float **transMat, float *initialDistribution, float **matchProb, float **insProb, float **transProb);
 596
 597 Constraint_list *ProbaMatrix2CL (Alignment *A, int *ns, int **ls, int NumMatrixTypes, int NumInsertStates, float *forward, float *backward, float thr, Constraint_list *CL);
 598
 599
 600 void free_proba_pair_wise ()
 601 {
 602   proba_pair_wise (NULL, NULL, NULL, NULL);
 603 }
 604 int proba_pair_wise ( Alignment *A, int *ns, int **ls, Constraint_list *CL)
 605 {
 606    int NumMatrixTypes=5;
 607    int NumInsertStates=2;
 608    static float **transMat, **insProb, **matchProb, *initialDistribution, **transProb, **emitPairs, *emitSingle, ***TinsProb, *TmatchProb;
 609    static int TinsProb_ml, TmatchProb_ml;
 610    int i, j,I, J;
 611    float *F, *B, *P;
 612    float tot;
 613    int l, s1, s2;
 614    float thr=0.01;//ProbCons Default
 615    char *alphabet;
 616
 617
 618
 619    //Free all the memory
 620    if (A==NULL)
 621      {
 622        free_float (transMat, -1);transMat=NULL;
 623        free_float (insProb, -1);insProb=NULL;
 624        free_float (matchProb, -1);matchProb=NULL;
 625        vfree (initialDistribution); initialDistribution=NULL;
 626        free_float (transProb, -1);transProb=NULL;
 627        free_float (emitPairs, -1);emitPairs=NULL;
 628        vfree (emitSingle);emitSingle=NULL;
 629
 630
 631        free_arrayN((void***)TinsProb, 3);TinsProb=NULL;
 632        vfree (TmatchProb);TmatchProb=NULL;
 633        TinsProb_ml=0; TmatchProb_ml=0;
 634
 635        forward_proba_pair_wise (NULL, NULL, 0,0,NULL,NULL,NULL,NULL,NULL);
 636        backward_proba_pair_wise (NULL, NULL, 0,0,NULL,NULL,NULL,NULL,NULL);
 637        ProbaMatrix2CL(NULL, NULL, NULL, 0, 0, NULL, NULL, 0, NULL);
 638        return 0;
 639      }
 640
 641    if (!transMat && (strm (retrieve_seq_type(), "DNA") ||strm (retrieve_seq_type(), "RNA")) )
 642      {
 643        static float **p;
 644        static float *s;
 645        if (!p)
 646          {
 647            int l,a,b;
 648            l=strlen (DNAalphabetDefault);
 649            p=declare_float (l,l);
 650            s=vcalloc (l, sizeof (float));
 651            for (a=0; a<l; a++)
 652              {
 653                s[a]=DNAemitSingleDefault[a];
 654                for (b=0; b<l; b++)
 655                  p[a][b]=DNAemitPairsDefault[a][b];
 656              }
 657          }
 658        p=get_emitPairs (CL->method_matrix, DNAalphabetDefault,p,s);
 659        alphabet=DNAalphabetDefault;
 660        emitPairs=declare_float (256, 256);
 661        emitSingle=vcalloc (256, sizeof (float));
 662        for (i=0; i<256; i++)
 663          {
 664            emitSingle[i]=1e-5;
 665            for (j=0; j<256; j++)
 666              emitPairs[i][j]=1e-10;
 667          }
 668        l=strlen (alphabet);
 669
 670        for (i=0; i<l; i++)
 671          {
 672            char C1,c1, C2,c2;
 673            c1=tolower(alphabet[i]);
 674            C1=toupper(alphabet[i]);
 675            emitSingle[c1]=s[i];
 676            emitSingle[C1]=s[i];
 677            for (j=0; j<=i; j++)
 678              {
 679                c2=tolower(alphabet[j]);
 680                C2=toupper(alphabet[j]);
 681
 682                emitPairs[c1][c2]=p[i][j];
 683                emitPairs[C1][c2]=p[i][j];
 684                emitPairs[C1][C2]=p[i][j];
 685                emitPairs[c1][C2]=p[i][j];
 686                emitPairs[c2][c1]=p[i][j];
 687                emitPairs[C2][c1]=p[i][j];
 688                emitPairs[C2][C1]=p[i][j];
 689                emitPairs[c2][C1]=p[i][j];
 690              }
 691          }
 692
 693
 694        transMat=declare_float (2*NumInsertStates+1, 2*NumInsertStates+1);
 695        transProb=declare_float (2*NumInsertStates+1,2* NumInsertStates+1);
 696        insProb=declare_float (256,NumMatrixTypes);
 697        matchProb=declare_float (256, 256);
 698        initialDistribution=vcalloc (2*NumMatrixTypes+1, sizeof (float));
 699
 700        ProbabilisticModel (NumMatrixTypes,NumInsertStates,initDistrib2Default, emitSingle,emitPairs,DNAgapOpen2Default,DNAgapExtend2Default, transMat,initialDistribution,matchProb, insProb,transProb);
 701      }
 702    else if ( !transMat && strm (retrieve_seq_type(), "PROTEIN"))
 703      {
 704        static float **p;
 705        static float *s;
 706        if (!p)
 707          {
 708            int l,a,b;
 709            l=strlen (alphabetDefault);
 710            p=declare_float (l,l);
 711            s=vcalloc (l, sizeof (float));
 712            for (a=0; a<l; a++)
 713              {
 714                s[a]=emitSingleDefault[a];
 715                for (b=0; b<l; b++)
 716                  p[a][b]=emitPairsDefault[a][b];
 717              }
 718          }
 719        p=get_emitPairs (CL->method_matrix, alphabetDefault,p,s);
 720        alphabet=alphabetDefault;
 721        emitPairs=declare_float (256, 256);
 722        emitSingle=vcalloc (256, sizeof (float));
 723        for (i=0; i<256; i++)
 724          {
 725            //emitSingle[i]=1e-5;
 726            emitSingle[i]=1;
 727            for (j=0; j<256; j++)
 728              //emitPairs[i][j]=1e-10;
 729              emitPairs[i][j]=1;
 730
 731          }
 732        l=strlen (alphabet);
 733
 734        for (i=0; i<l; i++)
 735          {
 736            char C1,c1, C2,c2;
 737            c1=tolower(alphabet[i]);
 738            C1=toupper(alphabet[i]);
 739            emitSingle[c1]=s[i];
 740            emitSingle[C1]=s[i];
 741            for (j=0; j<=i; j++)
 742              {
 743                c2=tolower(alphabet[j]);
 744                C2=toupper(alphabet[j]);
 745
 746                emitPairs[c1][c2]=p[i][j];
 747                emitPairs[C1][c2]=p[i][j];
 748                emitPairs[C1][C2]=p[i][j];
 749                emitPairs[c1][C2]=p[i][j];
 750                emitPairs[c2][c1]=p[i][j];
 751                emitPairs[C2][c1]=p[i][j];
 752                emitPairs[C2][C1]=p[i][j];
 753                emitPairs[c2][C1]=p[i][j];
 754
 755              }
 756          }
 757
 758
 759        transMat=declare_float (2*NumInsertStates+1, 2*NumInsertStates+1);
 760        transProb=declare_float (2*NumInsertStates+1,2* NumInsertStates+1);
 761        insProb=declare_float (256,NumMatrixTypes);
 762        matchProb=declare_float (256, 256);
 763        initialDistribution=vcalloc (2*NumMatrixTypes+1, sizeof (float));
 764
 765        ProbabilisticModel (NumMatrixTypes,NumInsertStates,initDistrib2Default, emitSingle,emitPairs,gapOpen2Default,gapExtend2Default, transMat,initialDistribution,matchProb, insProb,transProb);
 766      }
 767
 768    I=strlen (A->seq_al[ls[0][0]]);
 769    J=strlen (A->seq_al[ls[1][0]]);
 770    //TmatchProb=vcalloc ((I+1)*(J+1), sizeof (float));
 771    //TinsProb=declare_arrayN (3, sizeof (float),2,NumMatrixTypes,MAX(I,J)+1);
 772
 773    l=(I+1)*(J+1);
 774    if (l>TmatchProb_ml)
 775      {
 776        TmatchProb_ml=l;
 777        if (TmatchProb)TmatchProb=vrealloc(TmatchProb,TmatchProb_ml*sizeof (float));
 778        else TmatchProb=vcalloc ( l, sizeof (float));
 779      }
 780    l=MAX(I,J)+1;
 781    if ( l>TinsProb_ml)
 782      {
 783        TinsProb_ml=l;
 784        if (TinsProb)free_arrayN (TinsProb, 3);
 785        TinsProb=declare_arrayN (3, sizeof (float),2,NumMatrixTypes,TinsProb_ml);
 786      }
 787
 788    get_tot_prob (A,A, ns,ls,NumMatrixTypes, matchProb, insProb,TmatchProb,TinsProb, CL);
 789
 790    F=forward_proba_pair_wise (A->seq_al[ls[0][0]], A->seq_al[ls[1][0]], NumMatrixTypes,NumInsertStates,transMat, initialDistribution,TmatchProb,TinsProb, transProb);
 791    B=backward_proba_pair_wise (A->seq_al[ls[0][0]], A->seq_al[ls[1][0]], NumMatrixTypes,NumInsertStates,transMat, initialDistribution,TmatchProb,TinsProb, transProb);
 792    A->CL=ProbaMatrix2CL(A,ns, ls,NumMatrixTypes,NumInsertStates, F, B, thr,CL);
 793
 794    //free_proba_pair_wise();
 795    return 1;
 796    }
 797
 798 int get_tot_prob (Alignment *A1,Alignment *A2, int *ns, int **ls, int nstates, float **matchProb, float **insProb, float *TmatchProb, float ***TinsProb, Constraint_list *CL)
 799 {
 800   int i, j, a, b, c,d, k, n,n1,n2, ij;
 801   char c1, c2;
 802   int I, J;
 803   int ***VA1,***VA2, *observed, index;
 804   //Pre-computation of the pairwise scores in order to use potential profiles
 805   //The profiles are vectorized AND Compressed so that the actual alphabet size (proteins/DNA) does not need to be considered
 806
 807
 808   if (ns[0]==1 && ns[1]==1 )
 809     {
 810       int s1, s2;
 811       int *nns, **nls;
 812       Alignment *NA1, *NA2;
 813
 814       nns=vcalloc ( 2, sizeof (int));
 815       nls=vcalloc (2, sizeof (int*));
 816
 817       s1=A1->order[ls[0][0]][0];
 818       s2=A2->order[ls[1][0]][0];
 819       NA1=seq2R_template_profile (CL->S,s1);
 820       NA2=seq2R_template_profile (CL->S,s2);
 821
 822       if (NA1 || NA2)
 823         {
 824           if (NA1)
 825             {
 826               nns[0]=NA1->nseq;
 827               nls[0]=vcalloc (NA1->nseq, sizeof (int));
 828               for (a=0; a<NA1->nseq; a++)
 829                 nls[0][a]=a;
 830             }
 831           else
 832             {
 833             NA1=A1;
 834             nns[0]=ns[0];
 835             nls[0]=vcalloc (ns[0], sizeof (int));
 836             for (a=0; a<ns[0]; a++)
 837               nls[0][a]=ls[0][a];
 838             }
 839
 840           if (NA2)
 841             {
 842               nns[1]=NA2->nseq;
 843               nls[1]=vcalloc (NA2->nseq, sizeof (int));
 844               for (a=0; a<NA2->nseq; a++)
 845                 nls[1][a]=a;
 846             }
 847           else
 848             {
 849             NA2=A2;
 850             nns[1]=ns[1];
 851             nls[1]=vcalloc (ns[1], sizeof (int));
 852             for (a=0; a<ns[1]; a++)
 853               nls[1][a]=ls[1][a];
 854             }
 855
 856           get_tot_prob (NA1, NA2, nns, nls, nstates, matchProb, insProb, TmatchProb, TinsProb, CL);
 857           vfree (nns); free_int (nls,-1);
 858           return 1;
 859         }
 860     }
 861
 862   I=strlen (A1->seq_al[ls[0][0]]);
 863   J=strlen (A2->seq_al[ls[1][0]]);
 864
 865
 866
 867   //get Ins for I
 868   for (i=1; i<=I; i++)
 869     {
 870       for (k=0; k<nstates; k++)
 871         {
 872           TinsProb[0][k][i]=0;
 873           for (n=0,b=0; b<ns[0]; b++)
 874             {
 875               c1=A1->seq_al[ls[0][b]][i-1];
 876               if (c1!='-')
 877                 {
 878                   TinsProb[0][k][i]+=insProb[c1][k];
 879
 880                   n++;
 881                 }
 882             }
 883           if (n)TinsProb[0][k][i]/=n;
 884         }
 885     }
 886   //Get Ins for J
 887   for (j=1; j<=J; j++)
 888     {
 889       for (k=0; k<nstates; k++)
 890         {
 891           TinsProb[1][k][j]=0;
 892           for (n=0,b=0; b<ns[1]; b++)
 893             {
 894             c2=A2->seq_al[ls[1][b]][j-1];
 895             if (c2!='-')
 896               {
 897               TinsProb[1][k][j]+=insProb[c2][k];
 898
 899               n++;
 900               }
 901             }
 902           if (n)TinsProb[1][k][j]/=n;
 903         }
 904     }
 905
 906   observed=vcalloc ( 26, sizeof (int));
 907   VA1=declare_arrayN (3, sizeof (int),2,26,I);
 908   for (i=0; i<I; i++)
 909     {
 910       for (index=0, b=0; b<ns[0]; b++)
 911         {
 912           int in;
 913           c1=tolower(A1->seq_al[ls[0][b]][i]);
 914           if ( c1=='-')continue;
 915           c1-='a';
 916
 917           if (!(in=observed[c1])){in=observed[c1]=++index;}
 918
 919           VA1[0][in-1][i]=c1;
 920           VA1[1][in-1][i]++;
 921         }
 922
 923       VA1[0][index][i]=-1;
 924       for (b=0; b<26; b++)observed[b]=0;
 925     }
 926
 927   VA2=declare_arrayN (3, sizeof (int),2,26,J);
 928   for (i=0; i<J; i++)
 929     {
 930       for (index=0, b=0; b<ns[1]; b++)
 931         {
 932           int in;
 933
 934           c1=tolower(A2->seq_al[ls[1][b]][i]);
 935           if ( c1=='-')continue;
 936           c1-='a';
 937
 938           if (!(in=observed[c1])){in=observed[c1]=++index;}
 939
 940           VA2[0][in-1][i]=c1;
 941           VA2[1][in-1][i]++;
 942         }
 943       VA2[0][index][i]=-1;
 944       for (b=0; b<26; b++)observed[b]=0;
 945     }
 946   vfree (observed);
 947
 948   for ( ij=0,i=0; i<=I; i++)
 949     {
 950       for ( j=0; j<=J ; j++, ij++)
 951         {
 952           n=0;
 953           TmatchProb[ij]=0;
 954           if (i==0 || j==0);
 955           else
 956             {
 957               c=0;
 958               while (VA1[0][c][i-1]!=-1)
 959                 {
 960                   c1=VA1[0][c][i-1]+'a';
 961                   n1=VA1[1][c][i-1];
 962                   d=0;
 963                   while (VA2[0][d][j-1]!=-1)
 964                     {
 965                       c2=VA2[0][d][j-1]+'a';
 966                       n2=VA2[1][d][j-1];
 967                       TmatchProb[ij]+=matchProb[c1][c2]*(double)n1*(double)n2;
 968                       n+=n1*n2;
 969                       d++;
 970                     }
 971                   c++;
 972                 }
 973             }
 974           if (n)TmatchProb[ij]/=n;
 975         }
 976     }
 977
 978   free_arrayN ((void **)VA1, 3);
 979   free_arrayN ((void **)VA2, 3);
 980   return 1;
 981 }
 982
 983
 984
 985
 986 Constraint_list *ProbaMatrix2CL (Alignment *A, int *ns, int **ls, int NumMatrixTypes, int NumInsertStates, float *forward, float *backward, float thr, Constraint_list *CL)
 987 {
 988   float totalProb;
 989   int ij, i, j,k, I, J, s1, s2;
 990   static int *entry;
 991   static int **list;
 992   static int list_max;
 993   int sim;
 994   int list_size;
 995   int list_n;
 996   int old_n=0;
 997   double v;
 998   static float F=4; //potential number of full suboptimal alignmnents incorporated in the library
 999   static int tot_old, tot_new;
1000
1001   if (!A)
1002     {
1003       free_int (list, -1);list=NULL;
1004       list_max=0;
1005
1006       vfree(entry); entry=NULL;
1007       return NULL;
1008     }
1009
1010   I=strlen (A->seq_al[ls[0][0]]);
1011   J=strlen (A->seq_al[ls[1][0]]);
1012   s1=name_is_in_list (A->name[ls[0][0]], (CL->S)->name, (CL->S)->nseq, 100);
1013   s2=name_is_in_list (A->name[ls[1][0]], (CL->S)->name, (CL->S)->nseq, 100);
1014
1015   list_size=I*J;
1016
1017   if ( list_max<list_size)
1018     {
1019       free_int (list, -1);
1020       list_max=list_size;
1021       list=declare_int (list_max, 3);
1022     }
1023
1024
1025   totalProb = ComputeTotalProbability (I,J,NumMatrixTypes, NumInsertStates,forward, backward);
1026
1027   ij = 0;
1028   for (list_n=0,ij=0,i =0; i <= I; i++)
1029     {
1030       for (j =0; j <= J; j++, ij+=NumMatrixTypes)
1031         {
1032           v= EXP (MIN(LOG_ONE,(forward[ij] + backward[ij] - totalProb)));
1033           if (v>thr)//Conservative reduction of the list size to speed up the sorting
1034             {
1035               list[list_n][0]=i;
1036               list[list_n][1]=j;
1037               list[list_n][2]=(int)((float)v*(float)NORM_F);
1038               list_n++;
1039             }
1040           if (v>0.01)old_n++;
1041         }
1042     }
1043
1044   sort_int_inv (list, 3, 2, 0, list_n-1);
1045   if (!entry)entry=vcalloc ( CL->entry_len, CL->el_size);
1046   list_n=MIN(list_n,(F*MIN(I,J)));
1047   for (i=0; i<list_n; i++)
1048     {
1049        entry[SEQ1]=s1;
1050        entry[SEQ2]=s2;
1051        entry[R1]  =list[i][0];
1052        entry[R2]  =list[i][1];
1053        entry[WE]  =list[i][2];
1054        entry[CONS]=1;
1055        add_entry2list (entry,A->CL);
1056     }
1057   tot_new+=list_n;
1058   tot_old+=old_n;
1059   // HERE ("LIB_SIZE NEW: %d (new) %d (old) [%.2f]", list_n, old_n, (float)tot_new/(float)tot_old);
1060   return A->CL;
1061 }
1062
1063 Constraint_list *ProbaMatrix2CL_old (Alignment *A, int *ns, int **ls, int NumMatrixTypes, int NumInsertStates, float *forward, float *backward, float thr, Constraint_list *CL)
1064 {
1065   float totalProb;
1066   int ij, i, j,k, I, J, s1, s2;
1067   static int *entry;
1068   int lib_size=0;
1069   double v;
1070
1071   I=strlen (A->seq_al[ls[0][0]]);
1072   J=strlen (A->seq_al[ls[1][0]]);
1073   s1=name_is_in_list (A->name[ls[0][0]], (CL->S)->name, (CL->S)->nseq, 100);
1074   s2=name_is_in_list (A->name[ls[1][0]], (CL->S)->name, (CL->S)->nseq, 100);
1075
1076   totalProb = ComputeTotalProbability (I,J,NumMatrixTypes, NumInsertStates,forward, backward);
1077   if (!entry)entry=vcalloc ( CL->entry_len, CL->el_size);
1078   ij = 0;
1079   thr=0.01;
1080
1081
1082   for (ij=0,i =0; i <= I; i++)
1083     {
1084     for (j =0; j <= J; j++)
1085       {
1086         v= EXP (MIN(LOG_ONE,(forward[ij] + backward[ij] - totalProb)));
1087         if (i && j && v>=thr)
1088
1089           {
1090             entry[SEQ1]=s1;entry[SEQ2]=s2;
1091             entry[R1]=i;entry[R2]=j;
1092             entry[WE]=(int)((float)v*(float)NORM_F);
1093             entry[CONS]=1;
1094             add_entry2list (entry,A->CL);
1095             lib_size++;
1096           }
1097         ij += NumMatrixTypes;
1098       }
1099     }
1100   HERE ("LIB_SIZE_OLD: %d", lib_size);
1101   return A->CL;
1102 }
1103
1104 float ComputeTotalProbability (int seq1Length, int seq2Length,int NumMatrixTypes, int NumInsertStates,float *forward, float *backward)
1105 {
1106
1107     float totalForwardProb = LOG_ZERO;
1108     float totalBackwardProb = LOG_ZERO;
1109     int k;
1110
1111     for (k = 0; k < NumMatrixTypes; k++)
1112       {
1113       LOG_PLUS_EQUALS (&totalForwardProb,forward[k + NumMatrixTypes * ((seq1Length+1) * (seq2Length+1) - 1)] + backward[k + NumMatrixTypes * ((seq1Length+1) * (seq2Length+1) - 1)]);
1114       }
1115
1116     totalBackwardProb =forward[0 + NumMatrixTypes * (1 * (seq2Length+1) + 1)] +backward[0 + NumMatrixTypes * (1 * (seq2Length+1) + 1)];
1117
1118     for (k = 0; k < NumInsertStates; k++)
1119       {
1120       LOG_PLUS_EQUALS (&totalBackwardProb,forward[2*k+1 + NumMatrixTypes * (1 * (seq2Length+1) + 0)] +backward[2*k+1 + NumMatrixTypes * (1 * (seq2Length+1) + 0)]);
1121       LOG_PLUS_EQUALS (&totalBackwardProb,forward[2*k+2 + NumMatrixTypes * (0 * (seq2Length+1) + 1)] +backward[2*k+2 + NumMatrixTypes * (0 * (seq2Length+1) + 1)]);
1122       }
1123     return (totalForwardProb + totalBackwardProb) / 2;
1124   }
1125
1126
1127 float * backward_proba_pair_wise ( char *seq1, char *seq2, int NumMatrixTypes, int NumInsertStates, float **transMat, float *initialDistribution,float *matchProb, float ***insProb, float **transProb)
1128 {
1129   static float *backward;
1130   static int max_l;
1131
1132
1133   int k, i, j,ij, i1j1, i1j, ij1,a, l, seq1Length, seq2Length, m;
1134   char c1, c2;
1135   char *iter1, *iter2;
1136
1137   if (!seq1)
1138     {
1139       vfree (backward);
1140       backward=NULL; max_l=0;
1141       return NULL;
1142     }
1143
1144   iter1=seq1-1;
1145   iter2=seq2-1;
1146   seq1Length=strlen (seq1);
1147   seq2Length=strlen (seq2);
1148   l=(seq1Length+1)*(seq2Length+1)*NumMatrixTypes;
1149
1150   if (!backward)
1151     {
1152       backward=vcalloc (l, sizeof (float));
1153       max_l=l;
1154     }
1155   else if (max_l<l)
1156     {
1157       backward=vrealloc (backward, l*sizeof(float));
1158       max_l=l;
1159     }
1160
1161   for (a=0; a<l; a++)backward[a]=LOG_ZERO;
1162
1163   for (k = 0; k < NumMatrixTypes; k++)
1164     backward[NumMatrixTypes * ((seq1Length+1) * (seq2Length+1) - 1) + k] = initialDistribution[k];
1165
1166   // remember offset for each index combination
1167   ij = (seq1Length+1) * (seq2Length+1) - 1;
1168
1169   i1j = ij + seq2Length + 1;
1170   ij1 = ij + 1;
1171   i1j1 = ij + seq2Length + 2;
1172   ij *= NumMatrixTypes;
1173   i1j *= NumMatrixTypes;
1174   ij1 *= NumMatrixTypes;
1175   i1j1 *= NumMatrixTypes;
1176
1177   // compute backward scores
1178   for (i = seq1Length; i >= 0; i--)
1179     {
1180       c1 = (i == seq1Length) ? '~' : (unsigned char) iter1[i+1];
1181       for (j = seq2Length; j >= 0; j--)
1182         {
1183           c2 = (j == seq2Length) ? '~' : (unsigned char) iter2[j+1];
1184
1185           if (i < seq1Length && j < seq2Length)
1186             {
1187               m=((i+1)*(seq2Length+1))+j+1;//The backward and the forward are offset by 1
1188               float ProbXY = backward[0 + i1j1] + matchProb[m];
1189
1190
1191               for (k = 0; k < NumMatrixTypes; k++)
1192                 {
1193                   LOG_PLUS_EQUALS (&backward[k + ij], ProbXY + transProb[k][0]);
1194                 }
1195             }
1196           if (i < seq1Length)
1197             {
1198               for (k = 0; k < NumInsertStates; k++)
1199                 {
1200                 LOG_PLUS_EQUALS (&backward[0 + ij], backward[2*k+1 + i1j] + insProb[0][k][i+1] + transProb[0][2*k+1]);
1201                 LOG_PLUS_EQUALS (&backward[2*k+1 + ij], backward[2*k+1 + i1j] + insProb[0][k][i+1] + transProb[2*k+1][2*k+1]);
1202                 }
1203             }
1204         if (j < seq2Length)
1205           {
1206             for (k = 0; k < NumInsertStates; k++)
1207               {
1208                 //+1 because the backward and the forward are offset by 1
1209                 LOG_PLUS_EQUALS (&backward[0 + ij], backward[2*k+2 + ij1] + insProb[1][k][j+1] + transProb[0][2*k+2]);
1210                 LOG_PLUS_EQUALS (&backward[2*k+2 + ij], backward[2*k+2 + ij1] + insProb[1][k][j+1] + transProb[2*k+2][2*k+2]);
1211               }
1212           }
1213
1214         ij -= NumMatrixTypes;
1215         i1j -= NumMatrixTypes;
1216         ij1 -= NumMatrixTypes;
1217         i1j1 -= NumMatrixTypes;
1218         }
1219     }
1220
1221   return backward;
1222 }
1223 float * forward_proba_pair_wise ( char *seq1, char *seq2, int NumMatrixTypes, int NumInsertStates, float **transMat, float *initialDistribution,float *matchProb, float ***insProb, float **transProb)
1224 {
1225   static float *forward;
1226   static int max_l;
1227   int k, i, j,ij, i1j1, i1j, ij1, seq1Length, seq2Length, m;
1228   char *iter1, *iter2;
1229   int l,a;
1230
1231   if (!seq1)
1232     {
1233       vfree (forward);
1234       forward=NULL; max_l=0;
1235       return NULL;
1236     }
1237   iter1=seq1-1;
1238   iter2=seq2-1;
1239   seq1Length=strlen (seq1);
1240   seq2Length=strlen (seq2);
1241   l=(seq1Length+1)*(seq2Length+1)*NumMatrixTypes;
1242
1243   if (!forward)
1244     {
1245       forward=vcalloc (l, sizeof (float));
1246       max_l=l;
1247     }
1248   else if (max_l<l)
1249     {
1250       forward=vrealloc (forward, l*sizeof(float));
1251       max_l=l;
1252     }
1253   for (a=0; a<l; a++)forward[a]=LOG_ZERO;
1254
1255
1256   forward[0 + NumMatrixTypes * (1 * (seq2Length+1) + 1)] = initialDistribution[0] + matchProb[seq2Length+2];
1257
1258   for (k = 0; k < NumInsertStates; k++)
1259     {
1260       forward[2*k+1 + NumMatrixTypes * (1 * (seq2Length+1) + 0)] = initialDistribution[2*k+1] + insProb[0][k][1];
1261       forward[2*k+2 + NumMatrixTypes * (0 * (seq2Length+1) + 1)] = initialDistribution[2*k+2] + insProb[1][k][1];
1262     }
1263
1264   // remember offset for each index combination
1265     ij = 0;
1266     i1j = -seq2Length - 1;
1267     ij1 = -1;
1268     i1j1 = -seq2Length - 2;
1269
1270     ij *= NumMatrixTypes;
1271     i1j *= NumMatrixTypes;
1272     ij1 *= NumMatrixTypes;
1273     i1j1 *= NumMatrixTypes;
1274
1275
1276     // compute forward scores
1277     for (m=0,i = 0; i <= seq1Length; i++)
1278       {
1279         for (j = 0; j <= seq2Length; j++, m++)
1280         {
1281           if (i > 1 || j > 1)
1282           {
1283             if (i > 0 && j > 0)
1284               {
1285                 //Sum over all possible alignments
1286                 forward[0 + ij] = forward[0 + i1j1] + transProb[0][0];
1287                 for (k = 1; k < NumMatrixTypes; k++)
1288                   {
1289                     LOG_PLUS_EQUALS (&forward[0 + ij], forward[k + i1j1] + transProb[k][0]);
1290                   }
1291                 forward[0 + ij] += matchProb[m];
1292               }
1293             if ( i > 0)
1294               {
1295               for (k = 0; k < NumInsertStates; k++)
1296                 {
1297                   forward[2*k+1 + ij] = insProb[0][k][i] + LOG_ADD (forward[0 + i1j] + transProb[0][2*k+1],forward[2*k+1 + i1j] + transProb[2*k+1][2*k+1]);
1298                 }
1299               }
1300           if (j > 0)
1301             {
1302             for (k = 0; k < NumInsertStates; k++)
1303               {
1304                 forward[2*k+2 + ij] = insProb[1][k][j] +LOG_ADD (forward[0 + ij1] + transProb[0][2*k+2],forward[2*k+2 + ij1] + transProb[2*k+2][2*k+2]);
1305               }
1306             }
1307           }
1308
1309         ij += NumMatrixTypes;
1310         i1j += NumMatrixTypes;
1311         ij1 += NumMatrixTypes;
1312         i1j1 += NumMatrixTypes;
1313       }
1314
1315     }
1316     return forward;
1317   }
1318 int ProbabilisticModel (int NumMatrixTypes, int NumInsertStates,float *initDistribMat,float *emitSingle,  float **emitPairs, float *gapOpen, float *gapExtend, float **transMat, float *initialDistribution, float **matchProb, float **insProb, float **transProb)
1319 {
1320
1321
1322     // build transition matrix
1323   int i, j;
1324
1325
1326   transMat[0][0] = 1;
1327   for (i = 0; i < NumInsertStates; i++)
1328     {
1329     transMat[0][2*i+1] = gapOpen[2*i];
1330     transMat[0][2*i+2] = gapOpen[2*i+1];
1331     transMat[0][0] -= (gapOpen[2*i] + gapOpen[2*i+1]);
1332
1333     transMat[2*i+1][2*i+1] = gapExtend[2*i];
1334     transMat[2*i+2][2*i+2] = gapExtend[2*i+1];
1335     transMat[2*i+1][2*i+2] = 0;
1336     transMat[2*i+2][2*i+1] = 0;
1337     transMat[2*i+1][0] = 1 - gapExtend[2*i];
1338     transMat[2*i+2][0] = 1 - gapExtend[2*i+1];
1339     }
1340
1341
1342
1343   // create initial and transition probability matrices
1344   for (i = 0; i < NumMatrixTypes; i++){
1345     initialDistribution[i] = (float)log ((float)initDistribMat[i]);
1346     for (j = 0; j < NumMatrixTypes; j++)
1347       transProb[i][j] = (float)log ((float)transMat[i][j]);
1348   }
1349
1350   // create insertion and match probability matrices
1351   for (i = 0; i < 256; i++)
1352     {
1353       for (j = 0; j < NumMatrixTypes; j++)
1354         {
1355           insProb[i][j] = (float)log((float)emitSingle[i]);
1356         }
1357       for (j = 0; j < 256; j++)
1358         {
1359           matchProb[i][j] = (float)log((float)emitPairs[i][j]);
1360         }
1361     }
1362   return 1;
1363 }
1364
1365
1366 int viterbi_pair_wise ( Alignment *A, int *ns, int **ls, Constraint_list *CL)
1367 {
1368   char C1,c1, C2,c2;
1369   char *alphabet, *char_buf;
1370   char **al, **aln;
1371   int seq1Length, seq2Length, I, J;
1372   int i, j,ij, i1j1, i1j, ij1, k, a, b,l, LEN, r, c, m, state;
1373   int NumMatrixTypes=5;
1374   int NumInsertStates=2;
1375   int *traceback;
1376   float bestProb;
1377   static float **transMat, **insProb, **matchProb, *initialDistribution, **transProb, **emitPairs, *emitSingle, *TmatchProb, ***TinsProb;
1378   float *viterbi;
1379
1380   ungap_sub_aln (A, ns[0],ls[0]);
1381   ungap_sub_aln (A, ns[1],ls[1]);
1382
1383   seq1Length=I=strlen (A->seq_al[ls[0][0]]);
1384   seq2Length=J=strlen (A->seq_al[ls[1][0]]);
1385
1386   if (!transMat)
1387     {
1388        alphabet=alphabetDefault;
1389        emitPairs=declare_float (256, 256);
1390        emitSingle=vcalloc (256, sizeof (float));
1391        for (i=0; i<256; i++)
1392          {
1393            emitSingle[i]=1e-5;
1394            for (j=0; j<256; j++)
1395              emitPairs[i][j]=1e-10;
1396          }
1397        l=strlen (alphabet);
1398
1399        for (i=0; i<l; i++)
1400          {
1401
1402            c1=tolower(alphabet[i]);
1403            C1=toupper(alphabet[i]);
1404            emitSingle[c1]=emitSingleDefault[i];
1405            emitSingle[C1]=emitSingleDefault[i];
1406            for (j=0; j<=i; j++)
1407              {
1408                c2=tolower(alphabet[j]);
1409                C2=toupper(alphabet[j]);
1410
1411                emitPairs[c1][c2]=emitPairsDefault[i][j];
1412                emitPairs[C1][c2]=emitPairsDefault[i][j];
1413                emitPairs[C1][C2]=emitPairsDefault[i][j];
1414                emitPairs[c1][C2]=emitPairsDefault[i][j];
1415                emitPairs[c2][c1]=emitPairsDefault[i][j];
1416                emitPairs[C2][c1]=emitPairsDefault[i][j];
1417                emitPairs[C2][C1]=emitPairsDefault[i][j];
1418                emitPairs[c2][C1]=emitPairsDefault[i][j];
1419              }
1420          }
1421
1422
1423        transMat=declare_float (2*NumInsertStates+1, 2*NumInsertStates+1);
1424        transProb=declare_float (2*NumInsertStates+1,2* NumInsertStates+1);
1425        insProb=declare_float (256,NumMatrixTypes);
1426        matchProb=declare_float (256, 256);
1427        initialDistribution=vcalloc (2*NumMatrixTypes+1, sizeof (float));
1428
1429        ProbabilisticModel (NumMatrixTypes,NumInsertStates,initDistrib2Default, emitSingle,emitPairs,gapOpen2Default,gapExtend2Default, transMat,initialDistribution,matchProb, insProb,transProb);
1430      }
1431
1432
1433    TmatchProb=vcalloc ((I+1)*(J+1), sizeof (float));
1434    TinsProb=declare_arrayN (3, sizeof (float),2,NumMatrixTypes,MAX(I,J)+1);
1435    get_tot_prob (A,A, ns,ls,NumMatrixTypes, matchProb, insProb,TmatchProb,TinsProb, CL);
1436
1437    // create viterbi matrix
1438    l=NumMatrixTypes * (seq1Length+1) * (seq2Length+1);
1439    viterbi =vcalloc (l, sizeof (float));
1440    for (a=0; a<l; a++)viterbi[a]=LOG_ZERO;
1441    traceback=vcalloc (l, sizeof (int));
1442    for (a=0; a<l; a++)traceback[a]=-1;
1443
1444    // initialization condition
1445    for (k = 0; k < NumMatrixTypes; k++)
1446      viterbi[k] = initialDistribution[k];
1447
1448    // remember offset for each index combination
1449    ij = 0;
1450    i1j = -seq2Length - 1;
1451    ij1 = -1;
1452    i1j1 = -seq2Length - 2;
1453
1454    ij *= NumMatrixTypes;
1455    i1j *= NumMatrixTypes;
1456    ij1 *= NumMatrixTypes;
1457    i1j1 *= NumMatrixTypes;
1458
1459    // compute viterbi scores
1460    for (m=0,i = 0; i <= seq1Length; i++)
1461      {
1462        for ( j = 0; j <= seq2Length; j++, m++)
1463          {
1464            if (i > 0 && j > 0)
1465              {
1466                for (k = 0; k < NumMatrixTypes; k++)
1467                  {
1468                    float newVal = viterbi[k + i1j1] + transProb[k][0] + TmatchProb[m];
1469                    if (viterbi[0 + ij] < newVal)
1470                      {
1471                        viterbi[0 + ij] = newVal;
1472                        traceback[0 + ij] = k;
1473                      }
1474                  }
1475              }
1476            if (i > 0)
1477              {
1478                for (k = 0; k < NumInsertStates; k++)
1479                  {
1480                    float valFromMatch = TinsProb[0][k][i] + viterbi[0 + i1j] + transProb[0][2*k+1];
1481                    float valFromIns = TinsProb[0][k][i] + viterbi[2*k+1 + i1j] + transProb[2*k+1][2*k+1];
1482                    if (valFromMatch >= valFromIns){
1483                      viterbi[2*k+1 + ij] = valFromMatch;
1484                      traceback[2*k+1 + ij] = 0;
1485                    }
1486                    else {
1487                      viterbi[2*k+1 + ij] = valFromIns;
1488                      traceback[2*k+1 + ij] = 2*k+1;
1489                    }
1490                  }
1491              }
1492            if (j > 0)
1493              {
1494                for (k = 0; k < NumInsertStates; k++){
1495                  float valFromMatch = TinsProb[1][k][j] + viterbi[0 + ij1] + transProb[0][2*k+2];
1496                  float valFromIns = TinsProb[1][k][j] + viterbi[2*k+2 + ij1] + transProb[2*k+2][2*k+2];
1497                  if (valFromMatch >= valFromIns){
1498                    viterbi[2*k+2 + ij] = valFromMatch;
1499                    traceback[2*k+2 + ij] = 0;
1500                  }
1501                  else
1502                    {
1503                      viterbi[2*k+2 + ij] = valFromIns;
1504                      traceback[2*k+2 + ij] = 2*k+2;
1505                    }
1506                }
1507              }
1508
1509            ij += NumMatrixTypes;
1510            i1j += NumMatrixTypes;
1511            ij1 += NumMatrixTypes;
1512            i1j1 += NumMatrixTypes;
1513          }
1514      }
1515
1516    // figure out best terminating cell
1517    bestProb = LOG_ZERO;
1518    state = -1;
1519    for (k = 0; k < NumMatrixTypes; k++)
1520      {
1521        float thisProb = viterbi[k + NumMatrixTypes * ((seq1Length+1)*(seq2Length+1) - 1)] + initialDistribution[k];
1522        if (bestProb < thisProb)
1523          {
1524            bestProb = thisProb;
1525            state = k;
1526          }
1527      }
1528
1529
1530
1531    // compute traceback
1532    al=declare_char(2,seq1Length+seq2Length);
1533    LEN=0;
1534    r = seq1Length, c = seq2Length;
1535    while (r != 0 || c != 0)
1536      {
1537        int newState = traceback[state + NumMatrixTypes * (r * (seq2Length+1) + c)];
1538
1539        if (state == 0){ c--; r--; al[0][LEN]=1;al[1][LEN]=1;}
1540        else if (state % 2 == 1) {r--; al[0][LEN]=1;al[1][LEN]=0;}
1541        else { c--; al[0][LEN]=0;al[1][LEN]=1;}
1542        LEN++;
1543        state = newState;
1544      }
1545
1546
1547    invert_list_char ( al[0], LEN);
1548    invert_list_char ( al[1], LEN);
1549    if ( A->declared_len<=LEN)A=realloc_aln2  ( A,A->max_n_seq, 2*LEN);
1550    aln=A->seq_al;
1551    char_buf= vcalloc (LEN+1, sizeof (char));
1552    for ( c=0; c< 2; c++)
1553      {
1554        for ( a=0; a< ns[c]; a++)
1555          {
1556            int ch=0;
1557            for ( b=0; b< LEN; b++)
1558              {
1559                if (al[c][b]==1)
1560                  char_buf[b]=aln[ls[c][a]][ch++];
1561                else
1562                  char_buf[b]='-';
1563              }
1564            char_buf[b]='\0';
1565            sprintf (aln[ls[c][a]],"%s", char_buf);
1566          }
1567      }
1568
1569
1570    A->len_aln=LEN;
1571    A->nseq=ns[0]+ns[1];
1572    vfree (char_buf);
1573    free_char (al, -1);
1574
1575
1576
1577
1578
1579    return (int)(bestProb*(float)1000);
1580 }
1581
1582 float ** get_emitPairs (char *mat, char *alp, float **p, float *s)
1583   {
1584     static char *rmat;
1585     float k=0, t=0;
1586     int a, b, c, l;
1587     int **M;
1588
1589     if (!rmat)rmat=vcalloc (100, sizeof (char));
1590
1591     if (!mat || !mat[0] || strm (mat, "default"))return p;
1592     else if (strm (rmat, mat))return p;
1593
1594     sprintf (rmat,"%s", mat);
1595
1596     M=read_matrice (mat);
1597     l=strlen (alp);
1598
1599     k=log (2)/2;
1600     for (a=0; a<l; a++)
1601       for (b=0; b<l; b++)
1602         {
1603           int sc;
1604           float e;
1605           e=s[a]*s[b];
1606           sc=M[alp[a]-'A'][alp[b]-'A'];
1607           p[a][b]=e*exp ((double)sc*k);
1608         }
1609
1610     for (a=0; a<l; a++)
1611       for (b=0; b<l; b++)
1612         t+=p[a][b];
1613
1614     for (a=0; a<l; a++)
1615       for (b=0; b<l; b++)
1616         p[a][b]=p[a][b]/t;
1617
1618     t=0;
1619
1620     for (a=0; a<l; a++)
1621       for (b=0; b<l; b++)
1622         t+=p[a][b];
1623
1624     return p;
1625   }
1626
1627 /*********************************COPYRIGHT NOTICE**********************************/
1628 /*© Centro de Regulacio Genomica */
1629 /*and */
1630 /*Cedric Notredame */
1631 /*Tue Oct 27 10:12:26 WEST 2009. */
1632 /*All rights reserved.*/
1633 /*This file is part of T-COFFEE.*/
1634 /**/
1635 /*    T-COFFEE is free software; you can redistribute it and/or modify*/
1636 /*    it under the terms of the GNU General Public License as published by*/
1637 /*    the Free Software Foundation; either version 2 of the License, or*/
1638 /*    (at your option) any later version.*/
1639 /**/
1640 /*    T-COFFEE is distributed in the hope that it will be useful,*/
1641 /*    but WITHOUT ANY WARRANTY; without even the implied warranty of*/
1642 /*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the*/
1643 /*    GNU General Public License for more details.*/
1644 /**/
1645 /*    You should have received a copy of the GNU General Public License*/
1646 /*    along with Foobar; if not, write to the Free Software*/
1647 /*    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA*/
1648 /*...............................................                                                                                      |*/
1649 /*  If you need some more information*/
1650 /*  cedric.notredame@europe.com*/
1651 /*...............................................                                                                                                                                     |*/
1652 /**/
1653 /**/
1654 /*      */
1655 /*********************************COPYRIGHT NOTICE**********************************/