binaries/src/tcoffee/t_coffee_source/dp_lib/util_dp_gotoh_nw.c

   1 /******************************COPYRIGHT NOTICE*******************************/
   2 /*  (c) Centro de Regulacio Genomica                                                        */
   3 /*  and                                                                                     */
   4 /*  Cedric Notredame                                                                        */
   5 /*  12 Aug 2014 - 22:07.                                                                    */
   6 /*All rights reserved.                                                                      */
   7 /*This file is part of T-COFFEE.                                                            */
   8 /*                                                                                          */
   9 /*    T-COFFEE is free software; you can redistribute it and/or modify                      */
  10 /*    it under the terms of the GNU General Public License as published by                  */
  11 /*    the Free Software Foundation; either version 2 of the License, or                     */
  12 /*    (at your option) any later version.                                                   */
  13 /*                                                                                          */
  14 /*    T-COFFEE is distributed in the hope that it will be useful,                           */
  15 /*    but WITHOUT ANY WARRANTY; without even the implied warranty of                        */
  16 /*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the                         */
  17 /*    GNU General Public License for more details.                                          */
  18 /*                                                                                          */
  19 /*    You should have received a copy of the GNU General Public License                     */
  20 /*    along with Foobar; if not, write to the Free Software                                 */
  21 /*    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA             */
  22 /*...............................................                                           */
  23 /*  If you need some more information                                                       */
  24 /*  cedric.notredame@europe.com                                                             */
  25 /*...............................................                                           */
  26 /******************************COPYRIGHT NOTICE*******************************/
  27 #include <stdio.h>
  28 #include <stdlib.h>
  29 #include <math.h>
  30 #include <stdarg.h>
  31 #include <string.h>
  32 #include <limits.h>
  33
  34 #include "io_lib_header.h"
  35 #include "util_lib_header.h"
  36 #include "define_header.h"
  37 #include "dp_lib_header.h"
  38
  39 #define addE(i,j,d,s,list,n)\
  40   if (list)\
  41     {\
  42       int max_n;                                \
  43       Memcontrol *ppp;                          \
  44       if (!list[0])\
  45         {\
  46           list[0]=(int**)vcalloc ( 1000, sizeof (int*));\
  47         }\
  48                                                 \
  49       ppp=(Memcontrol*)list[0];                 \
  50       ppp-=2;                                           \
  51       max_n=ppp[0].size/sizeof(int*);\
  52       if (n[0]>=max_n){max_n+=1000;list[0]=(int**)vrealloc (list[0], max_n*sizeof (int*));} \
  53       if (!list[0][n[0]])list[0][n[0]]=(int*)vcalloc (7, sizeof (int)); \
  54       list[0][n[0]][0]=i;                                               \
  55       list[0][n[0]][1]=j;                                               \
  56       list[0][n[0]][3]=d;                                               \
  57       list[0][n[0]][2]=s;                                               \
  58       n[0]++;                                                           \
  59     }                                                                   \
  60
  61
  62
  63 void debug_list (char *tag, int n, int **list, int ex);
  64 void debug_list (char *tag, int n, int **list, int ex)
  65 {
  66   int a;
  67   for (a=0; a<n; a++)
  68     fprintf (stdout, "\n_%s_: %d %d %d %d\n", tag,list[a][0], list[a][1],list[a][2],list[a][4]);
  69   if (ex) exit(0);
  70 }
  71
  72 /*******************************************************************************/
  73 /*                idscore_pairseq: measure the % id without delivering thze aln*/
  74 /*                                                                             */
  75 /*      makes DP between the the ns[0] sequences and the ns[1]                 */
  76 /*                                                                             */
  77 /*      for MODE, see the function  get_dp_cost                                */
  78 /*******************************************************************************/
  79 int idscore_pairseq (char *s1, char *s2, int gop, int gep, int **m, char *comp_mode)
  80 {
  81   int **I, **D, **M, *P;
  82   int i, j,l1, l2, l,score, id, igop,match;
  83
  84
  85   l1=strlen (s1); l2=strlen (s2);
  86   lower_string (s1); lower_string (s2);
  87
  88   I=declare_int (6,l2+1);D=declare_int (6,l2+1);M=declare_int (6,l2+1);
  89   for (j=0; j<=l2; j++)
  90     {
  91       D[0][j]=gep*j;M[0][j]=2*gep*j;D[4][j]=0;
  92     }
  93
  94   for (i=1; i<=l1; i++)
  95     {
  96
  97       I[1][0]=i*gep;
  98       M[1][0]=2*i*gep;
  99
 100       for (j=1; j<=l2; j++)
 101         {
 102           score=m[s1[i-1]-'a'][s2[j-1]-'a'];
 103           id=(s1[i-1]==s2[j-1])?1:0;
 104
 105           igop=(i==l1 || j==l2)?0:gop;
 106
 107           if   ((D[0][j]+gep)>(M[0][j]+igop+gep))   {D[1][j]=D[0][j]+gep;      D[3][j]=D[2][j];        D[5][j]=D[4][j];}
 108           else                                      {D[1][j]=M[0][j]+igop+gep; D[3][j]=M[2][j];        D[5][j]=M[4][j];}
 109
 110           if ( (I[1][j-1]+gep)>(M[1][j-1]+igop+gep)){I[1][j]=I[1][j-1]+gep;      I[3][j]=I[3][j-1];    I[5][j]=I[5][j-1];}
 111           else                                      {I[1][j]=M[1][j-1]+igop+gep; I[3][j]=M[3][j-1];    I[5][j]=M[5][j-1];}
 112
 113           match=M[0][j-1]+score;
 114           if (I[1][j]>match && I[1][j]>D[1][j])     {M[1][j]=I[1][j]           ; M[3][j]=I[3][j];      M[5][j]=I[5][j];}
 115           else if (D[1][j]>match)                   {M[1][j]=D[1][j]           ; M[3][j]=D[3][j];      M[5][j]=D[5][j];}
 116           else                                      {M[1][j]=match             ; M[3][j]=M[2][j-1]+id; M[5][j]=M[4][j-1]+1;}
 117         }
 118       P=I[0]; I[0]=I[1]; I[1]=P;
 119       P=I[2]; I[2]=I[3]; I[3]=P;
 120       P=I[4]; I[4]=I[5]; I[5]=P;
 121
 122       P=D[0]; D[0]=D[1]; D[1]=P;
 123       P=D[2]; D[2]=D[3]; D[3]=P;
 124       P=D[4]; D[4]=D[5]; D[5]=P;
 125
 126       P=M[0]; M[0]=M[1]; M[1]=P;
 127       P=M[2]; M[2]=M[3]; M[3]=P;
 128       P=M[4]; M[4]=M[5]; M[5]=P;
 129     }
 130
 131
 132
 133
 134   if ( strstr (comp_mode, "sim2"))
 135     {
 136        l=MIN(l1,l2);
 137        score=(l==0)?0:(M[2][l2]*100)/l;
 138     }
 139   else if ( strstr (comp_mode, "sim3"))
 140     {
 141        l=MAX(l1,l2);
 142        score=(l==0)?0:(M[2][l2]*100)/l;
 143     }
 144   else if ( strstr (comp_mode, "cov"))
 145     {
 146       l=MAX(l1,l2);
 147       score=(l==0)?0:((M[4][l2]*100)/l);
 148     }
 149   else
 150     {
 151       //default: simple sim
 152       l=M[4][l2];
 153       score=(l==0)?0:(M[2][l2]*100)/l;
 154     }
 155
 156   free_int (I, -1);
 157   free_int (D, -1);
 158   free_int (M, -1);
 159
 160   return score;
 161 }
 162
 163 int test_pair_wise (Alignment *A, int *ns, int **l_s, Constraint_list *CL)
 164 {
 165   int a,l0, l1, n;
 166   char buf[VERY_LONG_STRING];
 167   char *gap, *seq;
 168
 169   l0=strlen (A->seq_al[l_s[0][0]]);
 170   l1=strlen (A->seq_al[l_s[1][0]]);
 171
 172   n=(l0<5)?l0/2:5;
 173   gap=generate_null(l1-n);
 174   for (a=0;a<ns[0]; a++)
 175     {
 176       seq=A->seq_al[l_s[0][a]];
 177       sprintf (buf, "%s%s",seq, gap);
 178       sprintf (seq, "%s", buf);
 179     }
 180   vfree (gap);
 181   gap=generate_null(l0-n);
 182
 183   for (a=0;a<ns[1]; a++)
 184     {
 185       seq=A->seq_al[l_s[1][a]];
 186       sprintf (buf, "%s%s",seq, gap);
 187       sprintf (seq, "%s", buf);
 188     }
 189   vfree(gap);
 190
 191
 192   A->len_aln=strlen (A->seq_al[l_s[0][0]]);
 193   A->score=A->score_aln=100;
 194   return 100;
 195 }
 196
 197 int idscore_pair_wise (Alignment *A,int*ns, int **l_s,Constraint_list *CL)
 198 {
 199
 200   A->score_aln=A->score=idscore_pairseq (A->seq_al[l_s[0][0]], A->seq_al[l_s[1][0]], CL->gop, CL->gep,CL->M, "sim3");
 201   return A->score_aln;
 202 }
 203 int dp_max (int *trace, int n, ...);
 204 int dp_max (int *trace, int n, ...)
 205 {
 206   va_list ap;
 207   int a, v, t, best_v=0;
 208
 209   va_start (ap, n);
 210   for (a=0; a< n; a++)
 211     {
 212       t=va_arg (ap, int);
 213       v=va_arg (ap, int);
 214
 215       if (a==0)
 216         {
 217           best_v=v;
 218           trace[0]=t;
 219         }
 220       else
 221         {
 222           if (v>best_v)
 223             {
 224               best_v=v;
 225               trace[0]=t;
 226             }
 227         }
 228     }
 229
 230   return best_v;
 231 }
 232 int is_tied (int *trace, int n, ...);
 233 int is_tied(int *trace, int n, ...)
 234 {
 235   va_list ap;
 236   int a, v, t, best_v=0;
 237   int nties=0;
 238
 239   va_start (ap, n);
 240   for (a=0; a< n; a++)
 241     {
 242       t=va_arg (ap, int);
 243       v=va_arg (ap, int);
 244
 245       if (a==0)
 246         {
 247           best_v=v;
 248           trace[0]=t;
 249         }
 250       else
 251         {
 252           if (v>best_v)
 253             {
 254               best_v=v;
 255               trace[0]=t;
 256             }
 257         }
 258     }
 259   va_end(ap);
 260   va_start (ap,n);
 261   for (a=0; a<n; a++)
 262     {
 263       t=va_arg (ap, int);
 264       v=va_arg (ap, int);
 265       if (v==best_v && trace[0]!=t)
 266         nties++;
 267     }
 268   va_end (ap);
 269   return nties;
 270 }
 271
 272 void display_mat (int **M, int l1, int l2, char *title);
 273 void display_mat (int **M, int l1, int l2, char *title)
 274 {
 275   int a, b;
 276
 277   fprintf ( stdout, "\n\nTitle %s\n", title);
 278   for ( a=0; a<=l1; a++)
 279     {
 280       fprintf ( stdout, "\n");
 281       for ( b=0; b<=l2; b++)
 282         fprintf ( stdout, "%3d ", M[a][b]);
 283     }
 284 }
 285 int glocal_pair_wise ( Alignment *A, int *ns, int **l_s, Constraint_list *CL)
 286 {
 287   int ***t, ***m;
 288   int i,j, l1, l2, n, sub, trace,ntrace, a, b, c, score;
 289   int gop,rgop,tgop, gep, unmatch;
 290   int M1, M2, I1, D1, LEN;
 291   char **al, *char_buf, **aln;
 292   int **pos0;
 293
 294
 295   l1=strlen (A->seq_al[l_s[0][0]]);
 296   l2=strlen (A->seq_al[l_s[1][0]]);
 297
 298   n=1;
 299   M1=n++;D1=n++;I1=n++;M2=n++;
 300   t=(int***)declare_arrayN(3, sizeof (int),n, l1+1, l2+1);
 301   m=(int***)declare_arrayN(3, sizeof (int),n, l1+1, l2+1);
 302
 303
 304   gop=CL->gop*SCORE_K;
 305   gep=CL->gep*SCORE_K;
 306   tgop=gop;
 307   unmatch=gep;
 308
 309   pos0=aln2pos_simple ( A,-1, ns, l_s);
 310
 311
 312   for (j=1; j<=l2; j++)
 313     {
 314       m[D1][0][j]=gep*j;
 315
 316       m[M1][0][j]=2*gep*j;
 317       m[M2][0][j]=4*gep*j;
 318     }
 319
 320
 321   for (i=1; i<=l1; i++)
 322     {
 323       m[I1][i][0]=i*gep;
 324       m[M2][i][0]=4*i*gep;
 325       m[M1][i][0]=2*i*gep;
 326
 327       for ( j=1; j<=l2; j++)
 328         {
 329           rgop=(i==l1 || j==1)?0:gop;
 330           rgop=gop;
 331           sub=(CL->get_dp_cost) (A, pos0, ns[0], l_s[0], i-1, pos0, ns[1], l_s[1],j-1,CL);
 332           m[M1][i][j]=dp_max (&trace,4,M1,m[M1][i-1][j-1],I1, m[I1][i-1][j-1],D1,m[D1][i-1][j-1],M2,m[M2][i-1][j-1])+sub;
 333           t[M1][i][j]=trace;
 334
 335           m[D1][i][j]=dp_max (&trace,3, M1,m[M1][i][j-1]+rgop,D1, m[D1][i][j-1]+gep, M2, m[M2][i][j-1]);
 336           t[D1][i][j]=trace;
 337
 338           m[I1][i][j]=dp_max (&trace,3, M1,m[M1][i-1][j]+rgop, I1, m[I1][i-1][j]+gep, M2, m[M2][i-1][j]);
 339           t[I1][i][j]=trace;
 340
 341           m[M2][i][j]=dp_max (&trace,4,M1,m[M1][i-1][j-1]+tgop,I1, m[I1][i-1][j-1]+tgop,D1,m[D1][i-1][j-1]+tgop,M2,m[M2][i-1][j-1])+unmatch;
 342           t[M2][i][j]=trace;
 343
 344         }
 345
 346     }
 347   score=dp_max (&trace,4, M1,m[M1][l1][l2],D1,m[D1][l1][l2],I1, m[I1][l1][l2],M2,m[M2][l1][l2]);
 348   LEN=0;i=l1;j=l2;
 349   al=declare_char (2, l1+l2+1);
 350
 351
 352   trace=t[trace][i][j];
 353   while (!(i==0 &&j==0))
 354     {
 355
 356       ntrace=t[trace][i][j];
 357       if (i==0)
 358         {
 359           al[0][LEN]=0;
 360           al[1][LEN]=1;
 361           j--;
 362           LEN++;
 363         }
 364       else if ( j==0)
 365         {
 366           al[0][LEN]=1;
 367           al[1][LEN]=0;
 368           i--;
 369           LEN++;
 370         }
 371       else if ( trace==M1)
 372         {
 373           al[0][LEN]=1;
 374           al[1][LEN]=1;
 375           i--; j--;
 376           LEN++;
 377         }
 378       else if ( trace==M2)
 379         {
 380           al[0][LEN]=1;
 381           al[1][LEN]=0;
 382           LEN++;
 383
 384           al[0][LEN]=0;
 385           al[1][LEN]=1;
 386           LEN++;
 387
 388           i--; j--;
 389
 390         }
 391       else if ( trace==D1)
 392         {
 393           al[0][LEN]=0;
 394           al[1][LEN]=1;
 395           j--;
 396           LEN++;
 397         }
 398       else if ( trace == I1)
 399         {
 400           al[0][LEN]=1;
 401           al[1][LEN]=0;
 402           i--;
 403           LEN++;
 404         }
 405       trace=ntrace;
 406
 407     }
 408
 409   invert_list_char ( al[0], LEN);
 410   invert_list_char ( al[1], LEN);
 411   if ( A->declared_len<=LEN)A=realloc_aln2  ( A,A->max_n_seq, 2*LEN);
 412
 413   aln=A->seq_al;
 414   char_buf=(char*)vcalloc (LEN+1, sizeof (char));
 415   for ( c=0; c< 2; c++)
 416     {
 417       for ( a=0; a< ns[c]; a++)
 418         {
 419           int ch=0;
 420           for ( b=0; b< LEN; b++)
 421             {
 422               if (al[c][b]==1)
 423                 char_buf[b]=aln[l_s[c][a]][ch++];
 424               else
 425                 char_buf[b]='-';
 426             }
 427           char_buf[b]='\0';
 428           sprintf (aln[l_s[c][a]],"%s", char_buf);
 429         }
 430     }
 431
 432
 433   A->len_aln=LEN;
 434   A->nseq=ns[0]+ns[1];
 435   free_arrayN((void *)m, 3);
 436   free_arrayN((void *)t, 3);
 437   vfree (char_buf);
 438   free_char (al, -1);
 439   return score;
 440 }
 441
 442 int ** aln2local_penalties (Alignment *A, int n, int *ls, Constraint_list *CL, int **lg);
 443 int ** aln2local_penalties (Alignment *A, int n, int *ls, Constraint_list *CL, int **lg)
 444 {
 445   //adapted from gap_count in MAFFT V 5.5
 446   int p,s,l, c1, c2;
 447   int gep,gop;
 448   int open=3, close=4, gap=5;
 449
 450   gop=CL->gop*SCORE_K;
 451   gep=CL->gep*SCORE_K;
 452
 453   l=strlen (A->seq_al[ls[0]]);
 454
 455   if (!lg)
 456     {
 457       lg=declare_int (6, l);
 458     }
 459
 460   if ( read_array_size_new (lg[0])<l)
 461     {
 462       free_int (lg, -1);
 463       lg=declare_int (6, l);
 464     }
 465
 466   for( s=0; s<n; s++ )
 467         {
 468           c1='x';
 469           for (p=0; p<l; p++)
 470             {
 471               c2=A->seq_al[ls[s]][p];
 472
 473               if (c1!='-' && c2=='-')lg[open][p]++;
 474               if (c1=='-' && c2!='-')lg[close][p]++;
 475               if ( c1=='-')lg[gap][p]++;
 476               c1=c2;
 477             }
 478         }
 479
 480   for (p=0; p<l; p++)
 481     {
 482       float go, gc, nn;
 483       nn=n;
 484       go=lg[open ][p];
 485       gc=lg[close][p];
 486
 487
 488       lg[GOP][p]=0.1*(1-(go/nn))*gop;
 489       lg[GCP][p]=0.1*(1-(gc/nn))*gop;
 490
 491       lg[GEP][p]=gep;//NO LOCAL GEP
 492       //(1-((float)lg[gap][p]/(float)n))*gep;//LOCAL_GEP
 493
 494       lg[open][p]=lg[close][p]=lg[gap][p]=0;
 495
 496     }
 497
 498   return lg;
 499 }
 500 int free_gotoh_pair_wise_lgp()
 501 {
 502   return gotoh_pair_wise_lgp (NULL, NULL, NULL, NULL);
 503 }
 504 int gotoh_pair_wise_lgp ( Alignment *A, int *ns, int **l_s, Constraint_list *CL)
 505 {
 506   int i,j, li, lj, n, sub, trace,ntrace, a, b, c, score;
 507   int I, J;
 508   int M1, I1, D1, LEN;
 509   char **al, *char_buf, **aln;
 510   int **pos0, **pos;
 511   Alignment *Aln;
 512
 513   int gop[2], gcp[2], gep[2];
 514   static int ***gpl, ***t, ***m;
 515   static int max_li, max_lj;
 516
 517
 518
 519   //gotoh_pair_wise ( A, ns, l_s,CL);
 520   //ungap_sub_aln (A, ns[0], l_s[0]);
 521   //ungap_sub_aln (A, ns[1], l_s[1]);
 522
 523   if (!A)
 524     {
 525       free_arrayN((void**)gpl, 3);
 526       free_arrayN((void**)t, 3);
 527       free_arrayN((void**)m, 3);
 528       max_li=max_lj=0;
 529       return 0;
 530     }
 531
 532   I=0;J=1;
 533
 534
 535   li=strlen (A->seq_al[l_s[I][0]]);
 536   lj=strlen (A->seq_al[l_s[J][0]]);
 537
 538   if ( !gpl)gpl=(int***)vcalloc ( 2, sizeof (int**));
 539   gpl[I]=aln2local_penalties (A,ns[I], l_s[I], CL,gpl[I]);
 540   gpl[J]=aln2local_penalties (A,ns[J], l_s[J], CL,gpl[J]);
 541
 542
 543   n=1;
 544   M1=n++;D1=n++;I1=n++;
 545
 546   if ( li>max_li ||lj>max_lj )
 547     {
 548       free_arrayN((void**)t, 3);
 549       free_arrayN((void**)m, 3);
 550
 551
 552       max_li=li;
 553       max_lj=lj;
 554       t=(int***)declare_arrayN(3, sizeof (int),n, max_li+1, max_lj+1);
 555       m=(int***)declare_arrayN(3, sizeof (int),n, max_li+1, max_lj+1);
 556
 557     }
 558   pos0=aln2pos_simple ( A,-1, ns, l_s);
 559
 560   //Compatibility with Macro
 561   Aln=A;
 562   pos=pos0;
 563
 564   for (j=1; j<=lj; j++)
 565     {
 566       gep[J]=gpl[J][GEP][j-1];
 567       m[D1][0][j]=gep[J]*j;
 568       m[I1][0][j]=m[D1][0][j]-1;
 569       m[M1][0][j]=m[D1][0][j]-1;
 570     }
 571
 572   //D1: gap in sequence I
 573   //I1: gap in sequence J
 574
 575
 576   for (i=1; i<=li; i++)
 577     {
 578       gep[I]=gpl[I][GEP][i-1];
 579       gop[I]=gpl[I][GOP][i-1];
 580       gcp[I]=gpl[I][GCP][i-1];
 581
 582       m[I1][i][0]=i*gep[I];
 583       m[D1][i][0]= m[I1][i][0]-1;
 584       m[M1][i][0]= m[I1][i][0]-1;
 585
 586
 587
 588       gop[I]=(i==1 || i==li )?0:gop[I];
 589       gcp[I]=(i==1 || i==li )?0:gcp[I];
 590
 591
 592       for ( j=1; j<=lj; j++)
 593         {
 594
 595           gep[J]=gpl[J][GEP][j-1];
 596           gop[J]=gpl[J][GOP][j-1];
 597           gcp[J]=gpl[J][GCP][j-1];
 598
 599           //gep[J]=gep[I]=(gep[J]+gep[I])/2;
 600           //gop[J]=gop[I]=(gop[J]+gop[I])/2;
 601           //gcp[J]=gcp[I]=(gcp[J]+gcp[I])/2;
 602
 603
 604           gop[J]=(j==1 || j==lj )?0:gop[J];
 605           gcp[J]=(j==1 || j==lj )?0:gcp[J];
 606
 607
 608           //sub=(CL->get_dp_cost) (A, pos0, ns[0], l_s[0], i-1, pos0, ns[1], l_s[1],j-1,CL);
 609           sub=TC_SCORE((i-1), (j-1));
 610
 611           m[M1][i][j]=dp_max (&trace,3,M1,m[M1][i-1][j-1],I1, m[I1][i-1][j-1]+gcp[I],D1,m[D1][i-1][j-1]+gcp[J])+sub;
 612           t[M1][i][j]=trace;
 613
 614
 615           m[D1][i][j]=dp_max (&trace,2, M1,m[M1][i][j-1]+gop[J]+gep[J],D1, m[D1][i][j-1]+gep[J]);
 616           t[D1][i][j]=trace;
 617
 618
 619           m[I1][i][j]=dp_max (&trace,2, M1,m[M1][i-1][j]+gop[I]+gep[I],I1, m[I1][i-1][j]+gep[I]);
 620           t[I1][i][j]=trace;
 621
 622         }
 623
 624     }
 625   score=dp_max (&trace,3, M1,m[M1][li][lj],D1,m[D1][li][lj],I1, m[I1][li][lj]);
 626
 627   LEN=0;i=li;j=lj;
 628   al=declare_char (2, li+lj);
 629
 630
 631   trace=t[trace][i][j];
 632   while (!(i==0 &&j==0))
 633     {
 634
 635       ntrace=t[trace][i][j];
 636
 637
 638       if (i==0)
 639         {
 640           al[0][LEN]=0;
 641           al[1][LEN]=1;
 642           j--;
 643           LEN++;
 644         }
 645       else if ( j==0)
 646         {
 647           al[0][LEN]=1;
 648           al[1][LEN]=0;
 649           i--;
 650           LEN++;
 651         }
 652       else if ( trace==M1)
 653         {
 654           al[0][LEN]=1;
 655           al[1][LEN]=1;
 656           i--; j--;
 657           LEN++;
 658         }
 659       else if ( trace==D1)
 660         {
 661           al[0][LEN]=0;
 662           al[1][LEN]=1;
 663           j--;
 664           LEN++;
 665         }
 666       else if ( trace == I1)
 667         {
 668           al[0][LEN]=1;
 669           al[1][LEN]=0;
 670           i--;
 671           LEN++;
 672         }
 673       trace=ntrace;
 674
 675     }
 676
 677   invert_list_char ( al[0], LEN);
 678   invert_list_char ( al[1], LEN);
 679   if ( A->declared_len<=LEN)A=realloc_aln  ( A,2*LEN+1);
 680
 681   aln=A->seq_al;
 682   char_buf=(char*) vcalloc (LEN+1, sizeof (char));
 683   for ( c=0; c< 2; c++)
 684     {
 685       for ( a=0; a< ns[c]; a++)
 686         {
 687           int ch=0;
 688           for ( b=0; b< LEN; b++)
 689             {
 690               if (al[c][b]==1)
 691                 char_buf[b]=aln[l_s[c][a]][ch++];
 692               else
 693                 char_buf[b]='-';
 694             }
 695           char_buf[b]='\0';
 696           sprintf (aln[l_s[c][a]],"%s", char_buf);
 697         }
 698     }
 699
 700
 701   A->len_aln=LEN;
 702   A->nseq=ns[0]+ns[1];
 703   vfree (char_buf);
 704   free_char (al, -1);
 705   free_int (pos0, -1);
 706   return score;
 707 }
 708 /*******************************************************************************/
 709 /*                GLOCAL 2                                                     */
 710 /*                                                                             */
 711 /*      makes DP between the the ns[0] sequences and the ns[1]                 */
 712 /*                                                                             */
 713 /*      for MODE, see the function  get_dp_cost                                */
 714 /*******************************************************************************/
 715 int glocal2_pair_wise (Alignment *IN,int*ns, int **ls,Constraint_list *CL)
 716 {
 717   int a, b, s=0;
 718   Alignment *A, *R,*L;
 719   char *seq, *buf;
 720
 721   buf=(char*)vcalloc (1000, sizeof (char));
 722   seq=(char*)vcalloc (1000, sizeof (char));
 723
 724   A=copy_aln (IN,NULL);
 725   L=copy_aln (IN,NULL);
 726   R=copy_aln (IN,NULL);
 727
 728   gotoh_pair_wise_sw (A, ns, ls, CL);
 729
 730
 731   for (a=0; a<2; a++)
 732     {
 733       for (b=0; b<ns[a]; b++)
 734         {
 735           s=ls[a][b];
 736           sprintf ( seq,"%s", IN->seq_al[s]);
 737
 738           seq[A->order[s][2]]='\0';
 739           sprintf (L->seq_al[s], "%s", seq);
 740           sprintf (R->seq_al[s], "%s", seq+A->order[s][3]+1);
 741         }
 742     }
 743   HERE ("2");
 744   print_sub_aln (A, ns, ls);
 745   gotoh_pair_wise(L, ns, ls, CL);
 746   print_sub_aln (L, ns, ls);
 747   gotoh_pair_wise(R, ns, ls, CL);
 748   print_sub_aln (R, ns, ls);
 749
 750   IN=realloc_aln (IN, A->len_aln+L->len_aln+R->len_aln+1);
 751   for (a=0; a<2; a++)
 752     {
 753       for (b=0; b<ns[a]; b++)
 754         {
 755           s=ls[a][b];
 756           sprintf (IN->seq_al[s], "%s%s%s",L->seq_al[s], A->seq_al[s], R->seq_al[s]);
 757         }
 758     }
 759   IN->len_aln=strlen (IN->seq_al[s]);
 760
 761   print_sub_aln (IN, ns, ls);
 762   vfree (seq); vfree (buf);
 763   free_aln (A); free_aln (L);free_aln (R);
 764   return IN->score_aln;
 765 }
 766
 767
 768 int gotoh_pair_wise (Alignment *A,int*ns, int **l_s,Constraint_list *CL)
 769         {
 770 /*******************************************************************************/
 771 /*                NEEDLEMAN AND WUNSCH (GOTOH)                                 */
 772 /*                                                                             */
 773 /*      makes DP between the the ns[0] sequences and the ns[1]                 */
 774 /*                                                                             */
 775 /*      for MODE, see the function  get_dp_cost                                */
 776 /*******************************************************************************/
 777
 778
 779 /*TREATMENT OF THE TERMINAL GAP PENALTIES*/
 780 /*TG_MODE=0---> gop and gep*/
 781 /*TG_MODE=1---> ---     gep*/
 782 /*TG_MODE=2---> ---     ---*/
 783
 784
 785             int TG_MODE;
 786             int l_gop, l_gep;
 787             int gop, gep;
 788             int maximise;
 789 /*VARIANLES FOR THE MULTIPLE SEQUENCE ALIGNMENT*/
 790         int a, b, i, j;
 791
 792         int *cc;
 793         int *dd,*ddg;
 794         int   e, eg;
 795
 796         int lenal[2], len;
 797         int t, c=0,s, ch;
 798         int sub;
 799         int fop;
 800         int score=0;
 801         int **pos0;
 802         static char **al;
 803         char **aln;
 804         int ala, alb,LEN;
 805         char *buffer;
 806         char *char_buf;
 807 /*trace back variables       */
 808         static int **trace;
 809         static int **bit;
 810         int *bi;
 811         int *tr;
 812         int dim;
 813         int ibit=0;
 814         int k;
 815         static int sample=0;//road==0, random tie; road=1: upper road; road=2 lower road;
 816         static int set_sample;
 817
 818         /********Prepare penalties*******/
 819         gop=CL->gop*SCORE_K;
 820         gep=CL->gep*SCORE_K;
 821         TG_MODE=CL->TG_MODE;
 822         maximise=CL->maximise;
 823
 824
 825 /********************************/
 826 /*CLEAN UP AFTER USE*/
 827         if ( A==NULL)
 828            {
 829            free_int (trace,-1);
 830            free_int (bit, -1);
 831            trace=NULL;
 832            bit=NULL;
 833
 834            free_char (al,-1);
 835            al=NULL;
 836            return 0;
 837            }
 838
 839 /*DO MEMORY ALLOCATION FOR DP*/
 840
 841         if (!set_sample)
 842           {
 843             set_sample=1;
 844             sample=atoigetenv ("SAMPLE_DP_4_TCOFFEE");
 845           }
 846
 847         lenal[0]=strlen (A->seq_al[l_s[0][0]]);
 848         lenal[1]=strlen (A->seq_al[l_s[1][0]]);
 849         len= MAX(lenal[0],lenal[1])+1;
 850
 851         buffer=(char*)vcalloc ( 2*len, sizeof (char));
 852         al=declare_char (2, 2*len);
 853
 854         char_buf=(char*) vcalloc (2*len, sizeof (char));
 855
 856
 857         dd =(int*)vcalloc (len, sizeof (int));
 858
 859
 860         cc =(int*) vcalloc (len, sizeof (int));
 861         ddg=(int*)vcalloc (len, sizeof (int));
 862
 863
 864
 865
 866
 867         dim=(trace==NULL)?0:read_size_int ( trace,sizeof (int*));
 868         trace    =realloc_int ( trace,dim,dim,MAX(0,len-dim), MAX(0,len-dim));
 869         bit      =realloc_int ( bit,dim,dim,MAX(0,len-dim), MAX(0,len-dim));
 870
 871 /*END OF MEMORY ALLOCATION*/
 872
 873
 874                 /*
 875                 0(s)   +(dd)
 876                   \      |
 877                    \     |
 878                     \    |
 879                      \   |
 880                       \  |
 881                        \ |
 882                         \|
 883                 -(e)----O
 884                 */
 885
 886         pos0=aln2pos_simple ( A,-1, ns, l_s);
 887
 888
 889         cc[0]=0;
 890         tr=trace[0];
 891         bi=bit[0];
 892         tr[0]=1;
 893         for ( j=1; j<=lenal[1]; j++)tr[j]=-1;
 894
 895         t=(TG_MODE==0)?gop:0;
 896
 897
 898         for (cc[0]=0,j=1; j<=lenal[1]; j++)
 899             {
 900
 901             l_gop=(TG_MODE==0)?gop:0;
 902             l_gep=(TG_MODE==2)?0:gep;
 903
 904             cc[j]=t=t+l_gep;
 905             dd[j]=  t+  gop;
 906             }
 907
 908         t=(TG_MODE==0)?gop:0;
 909
 910         for (i=1; i<=lenal[0];i++)
 911                         {
 912                           tr=trace[i];
 913                           bi=bit[i];
 914                           s=cc[0];
 915
 916                           l_gop=(TG_MODE==0)?gop:0;
 917                           l_gep=(TG_MODE==2)?0:gep;
 918
 919
 920
 921                           cc[0]=c=t=t+l_gep;
 922                           e=t+  gop;
 923                           tr[0]=1;
 924
 925
 926
 927                           for (eg=0,j=1; j<=lenal[1];j++)
 928                             {
 929
 930                               sub=(CL->get_dp_cost) (A, pos0, ns[0], l_s[0], i-1, pos0, ns[1], l_s[1],j-1,CL);
 931
 932                               /*get the best Insertion*/
 933                               l_gop=(i==lenal[0] || i==1 )?((TG_MODE==0)?gop:0):gop;
 934                               l_gep=(i==lenal[0] || i==1)?((TG_MODE==2)?0:gep):gep;
 935
 936
 937                               if ( a_better_than_b ( e,c+l_gop, maximise))eg++;
 938                               else eg=1;
 939                               e=best_of_a_b (e, c+l_gop, maximise)+l_gep;
 940
 941                               /*Get the best deletion*/
 942                               l_gop=(j==lenal[1] || j==1)?((TG_MODE==0)?gop:0):gop;
 943                               l_gep=(j==lenal[1] || j==1)?((TG_MODE==2)?0:gep):gep;
 944
 945
 946                               if ( a_better_than_b ( dd[j], cc[j]+l_gop, maximise))ddg[j]++;
 947                               else ddg[j]=1;
 948                               dd[j]=best_of_a_b( dd[j], cc[j]+l_gop,maximise)+l_gep;
 949
 950
 951
 952                               c=best_int(3,maximise,&fop, e, s+sub,dd[j]);
 953
 954
 955                               if (sample==1)
 956                                 {
 957                                   int rr[3];
 958                                   int nn=0;
 959                                   int fop2;
 960                                   int ind;
 961                                   if (c==e)rr[nn++]=0;
 962                                   if (c==(s+sub))rr[nn++]=1;
 963                                   if (c==dd[j])rr[nn++]=2;
 964                                   ind=rand()%(nn);
 965                                   fop=rr[ind];
 966                                   if (nn>1)
 967                                     {
 968                                       // HERE ("NN=%d index=%d",nn, ind);
 969                                       //HERE ("%d ->%d", fop, fop2);
 970                                       //HERE ("%d %d %d",  e, s+sub,dd[j]);
 971                                       ;
 972                                     }
 973                                 }
 974                               else if (sample==0)
 975                                 {
 976                                   /*Chose Substitution for tie breaking*/
 977                                   if ( fop==0 && (s+sub)==e)fop=1;
 978                                   else if ( fop==2 && (s+sub)==dd[j])fop=1;
 979                                   /*Chose Deletion for tie breaking*/
 980                                   else if ( fop==2 && e==dd[j])fop=2;
 981                                 }
 982                               else if (sample==-1)
 983                                 {
 984
 985                                   if ( fop==0 && (s+sub)==e)fop=1;
 986                                   else if ( fop==1 && (s+sub)==dd[j])fop=2;
 987                                   /*Chose Deletion for tie breaking*/
 988                                   else if ( fop==2 && e==dd[j])fop=1;
 989                                 }
 990                               bi[j]=0;
 991                               if (c==e){bi[j]++;}
 992                               if (c==(s+sub)){bi[j]++;}
 993                               if (c==(dd[j])){bi[j]++;}
 994                               //bi[j]--;
 995
 996
 997                               fop-=1;
 998                               s=cc[j];
 999                               cc[j]=c;
1000
1001
1002
1003
1004
1005                               if ( fop<0)
1006                                 {tr[j]=(TRACE_TYPE)fop*eg;
1007                                 }
1008                               else if ( fop>0)
1009                                 {tr[j]=(TRACE_TYPE)fop*ddg[j];
1010                                 }
1011                               else if (fop==0)
1012                                 {tr[j]=(TRACE_TYPE)0;
1013                                 }
1014                               fop= -2;
1015                             }
1016
1017                         }
1018
1019         score=c;
1020
1021         i=lenal[0];
1022         j=lenal[1];
1023         ala=alb=0;
1024
1025         if (!A->ibit)A->ibit=1; //set the bit counter on
1026         while (i>=0 && j>=0 && ((i+j)!=0))
1027                         {
1028                           if ( i==0)
1029                                 k=-1;
1030                         else if ( j==0)
1031                                 k=1;
1032                         else if ( j==0 && i==0)
1033                                 k=1;
1034                         else
1035                           {
1036                             k=trace[i][j];
1037                             A->ibit*=bit[i][j];
1038                           }
1039
1040
1041                         if (k==0)
1042                           {
1043
1044                             al[0][ala++]=1;
1045                             al[1][alb++]=1;
1046                             i--;
1047                             j--;
1048                           }
1049                         else if (k>0)
1050                           {
1051
1052                             for ( a=0; a< k; a++)
1053                               {
1054                                 al[0][ala++]=1;
1055                                 al[1][alb++]=0;
1056                                 i--;
1057                               }
1058                           }
1059                         else if (k<0)
1060                           {
1061
1062                             for ( a=0; a>k; a--)
1063                               {
1064                                 al[0][ala++]=0;
1065                                 al[1][alb++]=1;
1066                                 j--;
1067                               }
1068                           }
1069                         }
1070
1071         LEN=ala;
1072         c=LEN-1;
1073
1074
1075
1076         invert_list_char ( al[0], LEN);
1077         invert_list_char ( al[1], LEN);
1078         if ( A->declared_len<=LEN)A=realloc_aln2  ( A,A->max_n_seq, 2*LEN);
1079         aln=A->seq_al;
1080
1081         for ( c=0; c< 2; c++)
1082             {
1083             for ( a=0; a< ns[c]; a++)
1084                 {
1085                 ch=0;
1086                 for ( b=0; b< LEN; b++)
1087                     {
1088                     if (al[c][b]==1)
1089                         char_buf[b]=aln[l_s[c][a]][ch++];
1090                     else
1091                         char_buf[b]='-';
1092                    }
1093                 char_buf[b]='\0';
1094                 sprintf (aln[l_s[c][a]],"%s", char_buf);
1095                 }
1096              }
1097
1098
1099         A->len_aln=LEN;
1100         A->nseq=ns[0]+ns[1];
1101
1102
1103         vfree ( cc);
1104         vfree (dd);
1105         vfree (ddg);
1106         vfree (buffer);
1107         vfree (char_buf);
1108
1109         free_char ( al, -1);
1110         free_int (pos0, -1);
1111
1112
1113
1114
1115         return score;
1116         }
1117
1118
1119 int get_transition_cost (Alignment *A, int **posi, int ni, int *li, int i, int **posj, int nj, int *lj, int j,Constraint_list *CL);
1120 int gotoh_pair_wise_lgp_sticky ( Alignment *A, int *ns, int **l_s, Constraint_list *CL)
1121 {
1122   int i,j, li, lj, n, sub, trace,ntrace, a, b, c, score;
1123   int I, J;
1124   int M1, I1, D1, LEN;
1125   char **al, *char_buf, **aln;
1126   int **pos0;
1127
1128   int gop[2], gcp[2], gep[2];
1129   static int ***gpl, ***t, ***m;
1130   static int max_li, max_lj;
1131
1132
1133
1134   //gotoh_pair_wise ( A, ns, l_s,CL);
1135   //ungap_sub_aln (A, ns[0], l_s[0]);
1136   //ungap_sub_aln (A, ns[1], l_s[1]);
1137
1138   I=0;J=1;
1139
1140
1141   li=strlen (A->seq_al[l_s[I][0]]);
1142   lj=strlen (A->seq_al[l_s[J][0]]);
1143
1144   if ( !gpl)gpl=(int***)vcalloc ( 2, sizeof (int**));
1145   gpl[I]=aln2local_penalties (A,ns[I], l_s[I], CL,gpl[I]);
1146   gpl[J]=aln2local_penalties (A,ns[J], l_s[J], CL,gpl[J]);
1147
1148
1149   n=1;
1150   M1=n++;D1=n++;I1=n++;
1151
1152   if ( li>max_li ||lj>max_lj )
1153     {
1154       free_arrayN((void**)t, 3);
1155       free_arrayN((void**)m, 3);
1156
1157
1158       max_li=li;
1159       max_lj=lj;
1160       t=(int***)declare_arrayN(3, sizeof (int),n, max_li+1, max_lj+1);
1161       m=(int***)declare_arrayN(3, sizeof (int),n, max_li+1, max_lj+1);
1162
1163     }
1164   pos0=aln2pos_simple ( A,-1, ns, l_s);
1165
1166
1167   for (j=1; j<=lj; j++)
1168     {
1169       gep[J]=gpl[J][GEP][j-1];
1170       m[D1][0][j]=gep[J]*j;
1171       m[I1][0][j]=m[D1][0][j]-1;
1172       m[M1][0][j]=m[D1][0][j]-1;
1173     }
1174
1175   //D1: gap in sequence I
1176   //I1: gap in sequence J
1177
1178
1179   for (i=1; i<=li; i++)
1180     {
1181       gep[I]=gpl[I][GEP][i-1];
1182       gop[I]=gpl[I][GOP][i-1];
1183       gcp[I]=gpl[I][GCP][i-1];
1184
1185       m[I1][i][0]=i*gep[I];
1186       m[D1][i][0]= m[I1][i][0]-1;
1187       m[M1][i][0]= m[I1][i][0]-1;
1188
1189
1190
1191       gop[I]=(i==1 || i==li )?0:gop[I];
1192       gcp[I]=(i==1 || i==li )?0:gcp[I];
1193
1194
1195       for ( j=1; j<=lj; j++)
1196         {
1197           int transition;
1198
1199           gep[J]=gpl[J][GEP][j-1];
1200           gop[J]=gpl[J][GOP][j-1];
1201           gcp[J]=gpl[J][GCP][j-1];
1202
1203           //gep[J]=gep[I]=(gep[J]+gep[I])/2;
1204           //gop[J]=gop[I]=(gop[J]+gop[I])/2;
1205           //gcp[J]=gcp[I]=(gcp[J]+gcp[I])/2;
1206
1207
1208           gop[J]=(j==1 || j==lj )?0:gop[J];
1209           gcp[J]=(j==1 || j==lj )?0:gcp[J];
1210
1211
1212           sub=(CL->get_dp_cost) (A, pos0, ns[0], l_s[0], i-1, pos0, ns[1], l_s[1],j-1,CL);
1213           transition=get_transition_cost (A, pos0, ns[0], l_s[0], i-1, pos0, ns[1], l_s[1],j-1,CL);
1214
1215           m[M1][i][j]=dp_max (&trace,3,M1,m[M1][i-1][j-1]+transition,I1, m[I1][i-1][j-1]+gcp[I],D1,m[D1][i-1][j-1]+gcp[J])+sub;
1216           t[M1][i][j]=trace;
1217
1218
1219           m[D1][i][j]=dp_max (&trace,2, M1,m[M1][i][j-1]+gop[J]+gep[J],D1, m[D1][i][j-1]+gep[J]);
1220           t[D1][i][j]=trace;
1221
1222
1223           m[I1][i][j]=dp_max (&trace,2, M1,m[M1][i-1][j]+gop[I]+gep[I],I1, m[I1][i-1][j]+gep[I]);
1224           t[I1][i][j]=trace;
1225
1226         }
1227
1228     }
1229   score=dp_max (&trace,3, M1,m[M1][li][lj],D1,m[D1][li][lj],I1, m[I1][li][lj]);
1230
1231   LEN=0;i=li;j=lj;
1232   al=declare_char (2, li+lj);
1233
1234
1235   trace=t[trace][i][j];
1236   while (!(i==0 &&j==0))
1237     {
1238
1239       ntrace=t[trace][i][j];
1240
1241
1242       if (i==0)
1243         {
1244           al[0][LEN]=0;
1245           al[1][LEN]=1;
1246           j--;
1247           LEN++;
1248         }
1249       else if ( j==0)
1250         {
1251           al[0][LEN]=1;
1252           al[1][LEN]=0;
1253           i--;
1254           LEN++;
1255         }
1256       else if ( trace==M1)
1257         {
1258           al[0][LEN]=1;
1259           al[1][LEN]=1;
1260           i--; j--;
1261           LEN++;
1262         }
1263       else if ( trace==D1)
1264         {
1265           al[0][LEN]=0;
1266           al[1][LEN]=1;
1267           j--;
1268           LEN++;
1269         }
1270       else if ( trace == I1)
1271         {
1272           al[0][LEN]=1;
1273           al[1][LEN]=0;
1274           i--;
1275           LEN++;
1276         }
1277       trace=ntrace;
1278
1279     }
1280
1281   invert_list_char ( al[0], LEN);
1282   invert_list_char ( al[1], LEN);
1283   if ( A->declared_len<=LEN)A=realloc_aln  ( A,2*LEN+1);
1284
1285   aln=A->seq_al;
1286   char_buf=(char*) vcalloc (LEN+1, sizeof (char));
1287   for ( c=0; c< 2; c++)
1288     {
1289       for ( a=0; a< ns[c]; a++)
1290         {
1291           int ch=0;
1292           for ( b=0; b< LEN; b++)
1293             {
1294               if (al[c][b]==1)
1295                 char_buf[b]=aln[l_s[c][a]][ch++];
1296               else
1297                 char_buf[b]='-';
1298             }
1299           char_buf[b]='\0';
1300           sprintf (aln[l_s[c][a]],"%s", char_buf);
1301         }
1302     }
1303
1304
1305   A->len_aln=LEN;
1306   A->nseq=ns[0]+ns[1];
1307   vfree (char_buf);
1308   free_char (al, -1);
1309   free_int (pos0, -1);
1310   return score;
1311 }
1312 int get_transition_cost (Alignment *A, int **posi, int ni, int *li, int i, int **posj, int nj, int *lj, int j,Constraint_list *CL)
1313 {
1314   /*counts the number of identical transitions between position i-1, i and j-1..j*/
1315   float t=0;
1316   int a,s;
1317   Sequence *S;
1318
1319   if (i==0 || j==0)return 0;
1320
1321   for (a=0; a<ni; a++)
1322     {
1323       s=li[a];
1324       if (posi[s][i]<0 || posi[s][i-1]<0)continue;
1325       if (S->seq[li[a]][i-1]==S->seq[li[a]][i-1])t++;
1326     }
1327
1328   for (a=0; a<nj; a++)
1329     {
1330       s=lj[a];
1331       if (posj[s][j]<0 || posj[s][j-1]<0)continue;
1332       if (S->seq[li[a]][j-1]==S->seq[li[a]][j-1])t++;
1333     }
1334
1335   t=(t*10)/(float)(ni+nj);
1336   return t;
1337 }
1338 /*******************************************************************************/
1339 /*                idscore_pairseq: measure the % id without delivering thze aln*/
1340 /*                                                                             */
1341 /*      makes DP between the the ns[0] sequences and the ns[1]                 */
1342 /*                                                                             */
1343 /*      for MODE, see the function  get_dp_cost                                */
1344 /*******************************************************************************/
1345
1346 int cl2pair_list ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in, int mode, int ndiag);
1347 int cl2pair_list_ref ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
1348 int cl2pair_list_ecf ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
1349 int cl2pair_list_diag ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in, int add);
1350 int cl2list_borders   (Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
1351 int cl2diag_cap (Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in); //add one element at the end of each segment so that they can be joined
1352 int** cl2sorted_diagonals   ( Alignment *A, int *ns, int **ls, Constraint_list *CL);
1353 int** cl2sorted_diagonals_mat  ( Alignment *A, int *ns, int **ls, Constraint_list *CL);
1354 int** cl2sorted_diagonals_cs   ( Alignment *A, int *ns, int **ls, Constraint_list *CL);
1355 int list2nodup_list (Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
1356 int fill_matrix ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
1357
1358 int list2nodup_list ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in)
1359 {
1360   int **list;
1361   int n, a, b, c;
1362
1363   list=list_in[0];
1364   n=n_in[0];
1365
1366   if ( !A)return 0;
1367
1368
1369   sort_list_int (list,7, 1, 0, n-1);
1370   for (b=a=1; a<n; a++)
1371     {
1372       if (list[a][0]==list[b-1][0] && list[a][1]==list[b-1][1])
1373         {
1374           //HERE ("Duplicate");
1375            list[b-1][2]=MAX(list[b-1][2],list[a][2]);
1376         }
1377       else
1378         {
1379           for (c=0; c<4; c++)list[b][c]=list[a][c];
1380           b++;
1381         }
1382
1383         }
1384   n_in[0]=b;
1385   list_in[0]=list;
1386   return b;
1387 }
1388
1389 int cl2list_borders  (Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in)
1390 {
1391   int  a,n, p1, p2, l1, l2;
1392   int **list;
1393   int **pos;
1394   if (!A)return 0;
1395
1396
1397
1398
1399   l1=strlen (A->seq_al[ls[0][0]]);
1400   l2=strlen (A->seq_al[ls[1][0]]);
1401
1402   for (p1=0; p1<=l1; p1++)
1403     {
1404       if (p1==0 || p1==l1)
1405         {
1406           for (p2=0; p2<=l2; p2++)
1407             {
1408               addE(p1,p2,((l1-(p1))+(p2)),((CL->gep)*SCORE_K*p2), list_in, n_in);
1409             }
1410         }
1411       else
1412         {
1413           for (a=0; a<2; a++)
1414             {
1415               p2=(a==0)?0:l2;
1416
1417               addE(p1,p2,((l1-(p1))+(p2)),((CL->gep)*SCORE_K*p1), list_in, n_in);
1418             }
1419         }
1420     }
1421
1422   return read_array_size (list_in[0], sizeof (int*));
1423 }
1424
1425 int cl2diag_cap_r390 (Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in)
1426 {
1427   int **list;
1428   int n, in, a, b, al1, al2;
1429   int max_n;
1430   int cap=0;
1431
1432   if (!A) return 0;
1433
1434   al1=strlen (A->seq_al[ls[0][0]]);
1435   al2=strlen (A->seq_al[ls[1][0]]);
1436
1437   list=list_in[0];
1438   n=n_in[0];
1439   max_n=read_array_size (list, sizeof (int*));
1440
1441
1442
1443
1444   for (a=0; a< n; a++)
1445     {
1446       b=list[a][3];
1447       list[a][3]=list[a][0];
1448       list[a][0]=b;
1449
1450     }
1451   sort_list_int (list, 4, 1, 0, n-1);
1452   for (a=0; a< n; a++)
1453     {
1454       b=list[a][3];
1455       list[a][3]=list[a][0];
1456       list[a][0]=b;
1457     }
1458
1459
1460   in=n;
1461
1462   for (a=0; a<in; a++)
1463     {
1464       int i, j, pi, pj, ni, nj;
1465       if (list[a][2]==0)continue;
1466       i=list[a][0];
1467       j=list[a][1];
1468
1469       if (a==0){pi=-10;pj=-10;}
1470       else {pi=list[a-1][0];pj=list[a-1][1];}
1471
1472       if (a==in-1){ni=-10; nj=-10;}
1473       else {ni=list[a+1][0]; nj=list[a+1][1];}
1474
1475
1476       if (i==0 || j==0);
1477       else if ( i==pi || j==pj);
1478       else if ( i-pi!=1 || j-pj!=1)
1479         {
1480
1481
1482           int x;
1483           int delta=MAX((i-pi),(j-pj));
1484           for (x=1; x<=delta; x++)
1485             {
1486               if (n==max_n){max_n+=1000;list=(int**)vrealloc (list, max_n*sizeof (int*));}
1487               if (!list[n])list[n]=(int*)vcalloc (7, sizeof (int));
1488               if ((i-x)<=1)continue;
1489               if ((j-x)<=1)continue;
1490               list[n][0]=i-x;
1491               list[n][1]=j-x;
1492               list[n][3]=list[a][3];
1493               list[n][2]=cap;
1494               n++;
1495             }
1496         }
1497
1498
1499       if (i==al1 || j==al2);
1500       else if ( i==ni || j==nj);
1501       else if ( ni-i!=1 || nj-j!=1)
1502         {
1503           int x;
1504           int delta=MAX((ni-i),(nj-j));
1505           for (x=1; x<=delta; x++)
1506             {
1507                 if (n==max_n){max_n+=1000;list=(int**)vrealloc (list, max_n*sizeof (int*));}
1508                 if (!list[n])list[n]=(int*)vcalloc (7, sizeof (int));
1509
1510                 if (i+x>=al1)continue;
1511                 if (j+x>=al2) continue;
1512
1513
1514                 list[n][0]=i+x;
1515                 list[n][1]=j+x;
1516                 list[n][3]=list[a][3];
1517                 list[n][2]=cap;
1518                 n++;
1519               }
1520
1521         }
1522
1523     }
1524   list_in[0]=list;
1525   n_in[0]=n;
1526   return n;
1527 }
1528
1529 int cl2diag_cap (Alignment *A, int *nns, int **ls, Constraint_list *CL, int ***list, int *n)
1530 {
1531   int *sortseq;
1532
1533   int in, a, b, al1, al2;
1534   int max_n;
1535   int cap=0;
1536   int k=0;
1537
1538   static int **ll;
1539   static int max_ll;
1540   int nll=0;
1541
1542   int ns=0;
1543   int nt=0;
1544   int i,j,si,sj,ti,tj;
1545
1546   if ( !A)vfree (ll);max_ll=0;
1547
1548   al1=strlen (A->seq_al[ls[0][0]]);
1549   al2=strlen (A->seq_al[ls[1][0]]);
1550
1551   sortseq=(int*)vcalloc (7, sizeof (int));
1552   sortseq[0]=3;sortseq[1]=0;sortseq[2]=-1;
1553   sort_list_int2 (list[0], sortseq,4, 0, n[0]-1);
1554   vfree(sortseq);
1555   in=n[0];
1556
1557
1558   if (!ll){max_ll=100;ll=(int**)vcalloc(max_ll,sizeof(int*));}
1559
1560   for (a=0; a<in; a++)
1561     {
1562       int i, j, pi, pj, ni, nj;
1563       if (list[0][a][2]==0)continue;//this is where borders are avoided
1564       i=list[0][a][0];
1565       j=list[0][a][1];
1566
1567       if (a==0){pi=-10;pj=-10;}
1568       else {pi=list[0][a-1][0];pj=list[0][a-1][1];}
1569
1570       if (a==in-1){ni=-10; nj=-10;}
1571       else {ni=list[0][a+1][0]; nj=list[0][a+1][1];}
1572
1573
1574       if ((i==0 || j==0));
1575       else if ( i==pi || j==pj);
1576       else if ( i-pi!=1 || j-pj!=1)
1577         {
1578           if (nll>=max_ll){max_ll+=1000;ll=(int**)vrealloc (ll, max_ll*sizeof (int*));}
1579           ll[nll++]=list[0][a];
1580           list[0][a][6]=_START;
1581         }
1582
1583       if (i==al1 || j==al2);
1584       else if ( i==ni || j==nj);
1585       else if ( ni-i!=1 || nj-j!=1)
1586         {
1587           if (nll>=max_ll){max_ll+=1000;ll=(int**)vrealloc (ll, max_ll*sizeof (int*));}
1588           ll[nll++]=list[0][a];
1589           list[0][a][6]=_TERM;
1590         }
1591     }
1592
1593   sortseq=(int*)vcalloc (7, sizeof (int));
1594   sortseq[0]=0;sortseq[1]=1;sortseq[2]=-1;
1595   sort_list_int2 (ll, sortseq,4, 0,nll-1);
1596   vfree (sortseq);
1597
1598   for (a=0; a<nll; a++)
1599     {
1600       int ci, nl,max_nl,best_d,d,best_s;
1601       max_nl=100;
1602       if (ll[a][6]!=_TERM)continue;
1603
1604       ti=ll[a][0];
1605       tj=ll[a][1];
1606       ci=ti;
1607
1608       for (nl=0,best_d=-1,b=a+1;b<nll && nl<max_nl; b++)
1609         {
1610           if (ll[b][6]!=_START)continue;
1611
1612           si=ll[b][0];
1613           sj=ll[b][1];
1614
1615           if (si>ci){nl++;ci=si;}
1616           d=MIN((si-ti), (sj-tj));
1617           if (d<=0);
1618           else if (best_d==-1 || best_d>d){best_d=d; best_s=b;}
1619         }
1620       if (best_d==-1)continue;
1621
1622       si=ll[best_s][0];
1623       sj=ll[best_s][1];
1624
1625       for (i=ti, j=tj; (i<=si && j<=sj); i++, j++)//extend the top diagonal
1626         {
1627           addE(i,j,(al1-i+j),cap, list,n);
1628         }
1629
1630       for (i=si, j=sj; (i>=ti && j>=tj); i--, j--)//extend the bottom diagonal
1631         {
1632           addE(i,j,(al1-i+j),cap, list,n);
1633         }
1634     }
1635
1636   for (a=0; a<nll; a++)ll[a][6]=0;
1637
1638   return n[0];
1639 }
1640 int cl2pair_list_ext ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
1641
1642 /**
1643  * Calculates scores for diagonal segments.
1644  *
1645  * \param Alignment The sequences.
1646  * \param ns Number of sequences in each group
1647  * \param ls sequences in in groups (ls[0][x] sequences in group 1, ls[1][x] squences in group 2).
1648  * \param CL the constraint list
1649  * \param list_in the diagonals
1650  * \param n_in number of sequences?
1651  */
1652
1653 int fork_cl2pair_list_ext ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in, int njobs);
1654 int cl2pair_list_ext ( Alignment *A, int *ins, int **ils, Constraint_list *CL, int ***list_in, int *n_in)
1655 {
1656   int *ns;
1657   int **ls;
1658   int ret;
1659   int njobs;
1660
1661   if (!CL || !CL->S || !CL->residue_index) return 0;
1662
1663   if (read_size_int(ins,sizeof (int))==3 && ins[2]!=-1)
1664     {
1665       ns=(int*)vcalloc (2, sizeof (int));
1666       ls=declare_int (2,1);
1667       ns[0]=ns[1]=1;
1668       ls[0][0]=ils[0][ins[0]];
1669       ls[1][0]=ils[1][ins[1]];
1670     }
1671   else
1672     {
1673       ns=ins;
1674       ls=ils;
1675     }
1676
1677   if (strstr ( CL->multi_thread, "pairwise"))njobs=get_nproc();
1678   else njobs=1;
1679
1680   ret=fork_cl2pair_list_ext(A,ns,ls,CL,list_in,n_in,njobs);
1681
1682   if (read_size_int (ins, sizeof(int))==3)
1683     {
1684       vfree (ns);
1685       free_int (ls,-1);
1686     }
1687   return ret;
1688 }
1689
1690
1691 int fork_cl2pair_list_ext ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in,int njobs)
1692 {
1693   int p1, p2,diag, si, s, r, t_s, t_r,t_w, t_s2, t_r2, t_w2;
1694   int a, b,g, l1, l2;
1695   int **pos;
1696   int ll[3];
1697
1698   int nused;
1699   int *used_list;
1700   int *sl2,*sl1, **inv_pos;
1701
1702   int normalisation_mode=0;
1703
1704   float nscore, score, tot, filter, avg=0;
1705   float **used;
1706   float *norm;
1707   //Maria added this
1708   int norm2=0;  //
1709   static int **nr;
1710   static int max_nr;
1711   //variables for fork
1712   FILE *fp;
1713   char **pid_tmpfile;
1714   int sjobs,j;
1715   int **sl;
1716
1717   if (!A) return 0;
1718
1719
1720
1721   pos=aln2pos_simple ( A,-1, ns, ls);
1722   inv_pos=(int**)vcalloc ((CL->S)->nseq, sizeof (int*));
1723   for (a=0; a<ns[1]; a++)inv_pos[ls[1][a]] =seq2inv_pos(A->seq_al[ls[1][a]]);
1724
1725   ll[0]=l1=strlen (A->seq_al[ls[0][0]]);
1726   ll[1]=l2=strlen (A->seq_al[ls[1][0]]);
1727   ll[2]=MAX(l1,l2);
1728   sl1=(int*)vcalloc ((CL->S)->nseq, sizeof (int));
1729   sl2=(int*)vcalloc ((CL->S)->nseq, sizeof (int));
1730
1731   for (a=0;a<ns[0]; a++)sl1[ls[0][a]]=1;
1732   for (a=0;a<ns[1]; a++)sl2[ls[1][a]]=1;
1733   norm= (float*)vcalloc ( l1+1, sizeof (float));
1734
1735   //data structure used for norm 2
1736
1737
1738   if (!nr || ll[2]>max_nr)
1739     {
1740       if (nr)free_int (nr, -1);
1741       max_nr=ll[2];
1742       nr=declare_int (2, max_nr+1);
1743     }
1744
1745   for (g=0; g<2; g++)
1746     {
1747       for (a=0; a<ll[g]; a++)
1748         {
1749           nr[g][a+1]=0;
1750           for (b=0; b<ns[g]; b++)
1751             if (A->seq_al[ls[g][b]][a]!='-')nr[g][a+1]++;
1752         }
1753     }
1754
1755   sl=n2splits (njobs,l1+1);
1756   pid_tmpfile=(char**)vcalloc (njobs, sizeof (char*));
1757
1758   used=declare_float (l2+1,2);
1759   used_list=(int*)vcalloc (l2+1, sizeof (int));
1760   nused=0;
1761
1762   for (sjobs=0, j=0; j<njobs; j++)
1763     {
1764       pid_tmpfile[j]=vtmpnam(NULL);
1765       if (vvfork (NULL)==0)
1766         {
1767           initiate_vtmpnam(NULL);
1768           fp=vfopen (pid_tmpfile[j], "w");
1769           for (p1=sl[j][0]; p1<sl[j][1]; p1++)
1770             {
1771               for (tot=0,nused=0,si=0;p1>0 && si<ns[0]; si++)
1772                 {
1773                   s=ls [0][si];r=pos[s][p1-1];
1774                   for (a=1; r>0 && a<CL->residue_index[s][r][0];a+=ICHUNK)
1775                     {
1776                       t_s=CL->residue_index[s][r][a+SEQ2];
1777                       t_r=CL->residue_index[s][r][a+R2];
1778                       t_w=CL->residue_index[s][r][a+WE];
1779                       if (sl1[t_s])continue;//do not extend within a profile
1780
1781                       norm[p1]++;
1782                       norm2=0;
1783                       for (b=0; b<CL->residue_index[t_s][t_r][0];)
1784                         {
1785                           if (b==0){t_s2=t_s;t_r2=t_r;t_w2=t_w;b++;}
1786                           else
1787                             {
1788                               t_s2=CL->residue_index[t_s][t_r][b+SEQ2];
1789                               t_r2=CL->residue_index[t_s][t_r][b+R2];
1790                               t_w2=CL->residue_index[t_s][t_r][b+WE];
1791                               b+=ICHUNK;
1792                               norm2++;
1793                             }
1794                           if (sl2[t_s2])
1795                             {
1796                               p2=inv_pos[t_s2][t_r2];
1797                               score=MIN(((float)t_w/(float)NORM_F),((float)t_w2/(float)NORM_F));
1798
1799                               if (!used[p2][1] && score>0)
1800                                 {
1801                                   used_list[nused++]=p2;
1802                                 }
1803
1804                               tot+=score;
1805                               used[p2][0]+=score;
1806                               used[p2][1]++;
1807                             }
1808                         }
1809                     }
1810                 }
1811               filter=0.01;
1812               for (a=0; a<nused; a++)
1813                 {
1814
1815                   p2=used_list[a];
1816                   nscore=used[p2][0]/tot; //Normalized score used for filtering
1817                   score =used[p2][0];
1818                   used[p2][0]=used[p2][1]=0;
1819
1820                   if (nscore>filter && p1!=0 && p2!=0 && p1!=l1 && p2!=l2)
1821                     {
1822                       if (normalisation_mode==0)
1823                         {
1824                           score=((norm[p1]>0)?score/norm[p1]:0);
1825                         }
1826                       else if (normalisation_mode==1)
1827                         {
1828                           score/=(float)((CL->S)->nseq*nr[0][p1]*nr[1][p2]);
1829                         }
1830
1831                       score*=NORM_F;
1832                       fprintf (fp, "%d %d %d %f ", p1, p2, ((l1-(p1))+(p2)), score);
1833                     }
1834                 }
1835             }
1836           vfclose (fp);
1837           myexit (EXIT_SUCCESS);
1838         }
1839       else
1840         {
1841           sjobs++;
1842         }
1843     }
1844   while (sjobs>=0){vwait(NULL); sjobs--;}//wait for all jobs to complete
1845   for (j=0; j<njobs; j++)
1846     {
1847       fp=vfopen (pid_tmpfile[j], "r");
1848       while ((fscanf(fp, "%d %d %d %f ", &p1,&p2, &diag, &score))==4)
1849         addE (p1,p2,((l1-(p1))+(p2)),score,list_in, n_in);
1850       vfclose (fp);
1851       remove (pid_tmpfile[j]);
1852     }
1853
1854   free_float (used, -1);
1855   vfree (used_list);
1856   free_int (inv_pos, -1);
1857   free_int (pos, -1);
1858   vfree (sl2);vfree (sl1);
1859   vfree(norm);
1860   return n_in[0];
1861 }
1862
1863
1864
1865 int list2linked_pair_wise ( Alignment *A, int *ns, int **l_s, Constraint_list *CL, int **list, int n, char ***al, int *len);
1866 int linked_pair_wise ( Alignment *A, int *ns, int **ls, Constraint_list *CL)
1867 {
1868   int n=0;
1869   static int **list=NULL;
1870   int score, a;
1871   char **al;
1872   int len=0;
1873   int invert=0;
1874   int tr0,tr1;
1875
1876   if ( !A)free_int (list, -1);
1877   if ( !CL->residue_index)return myers_miller_pair_wise (A, ns,ls,CL);
1878
1879
1880   tr0=ns[0]*strlen (A->seq_al[ls[0][0]]);
1881   tr1=ns[1]*strlen (A->seq_al[ls[1][0]]);
1882
1883   if (tr0>tr1)
1884     {
1885       int *ins;
1886       int **ils;
1887       int a,b,c;
1888       invert=1;
1889       ins=(int*)vcalloc (2, sizeof(int));
1890       ils=declare_int (2, (CL->S)->nseq);
1891
1892       for ( a=0; a<2; a++)
1893         {
1894           ins[a]=ns[a];
1895           for (b=0; b<ns[a]; b++)ils[a][b]=ls[a][b];
1896         }
1897
1898       for (c=1,a=0; a<2; a++,c--)
1899         {
1900           ns[c]=ins[a];
1901           for (b=0; b<ins[a]; b++)
1902             ls[c][b]=ils[a][b];
1903         }
1904
1905       vfree (ins);
1906       free_int (ils, -1);
1907     }
1908
1909
1910
1911
1912
1913   /*Prepare the list*/
1914
1915
1916   cl2pair_list_ext    (A, ns, ls, CL, &list, &n);
1917   cl2diag_cap         (A, ns, ls, CL, &list, &n);
1918   cl2list_borders     (A, ns, ls, CL, &list, &n);
1919   list2nodup_list     (A, ns, ls, CL, &list, &n);
1920   /*Do the DP*/
1921   score=list2linked_pair_wise (A, ns, ls, CL, list, n, &al,&len);
1922   free_char (al, -1);
1923
1924   if (invert)
1925     {
1926       int *ins;
1927       int **ils;
1928       int a,b,c;
1929
1930       ins=(int*)vcalloc (2, sizeof(int));
1931       ils=declare_int (2, (CL->S)->nseq);
1932
1933       for ( a=0; a<2; a++)
1934         {
1935           ins[a]=ns[a];
1936           for (b=0; b<ns[a]; b++)ils[a][b]=ls[a][b];
1937         }
1938
1939       for (c=1,a=0; a<2; a++,c--)
1940         {
1941           ns[c]=ins[a];
1942           for (b=0; b<ins[a]; b++)
1943             ls[c][b]=ils[a][b];
1944         }
1945
1946       vfree (ins);
1947       free_int (ils, -1);
1948     }
1949
1950   /*Free the list*/
1951   return score;
1952 }
1953
1954 #define LIN(a,b,c) a[b*5+c]
1955 int list2linked_pair_wise( Alignment *A, int *ns, int **l_s, Constraint_list *CL, int **list, int n, char ***tb, int *len)
1956 {
1957   int a,b,c, i, j, LEN=0, start_trace;
1958   int pi, pj,ij, delta_i, delta_j, prev_i, prev_j;
1959   static int **slist;
1960   static long *MI, *MJ, *MM,*MT2;
1961   static int *sortseq;
1962   static int max_size;
1963   int gop, gep, igop, igep;
1964   int l1, l2, l, ls;
1965   char **al;
1966   char **aln,*char_buf;
1967   int ni=0, nj=0;
1968   long score;
1969   int nomatch;
1970
1971
1972
1973   l1=strlen (A->seq_al[l_s[0][0]]);
1974   l2=strlen (A->seq_al[l_s[1][0]]);
1975   al=declare_char (2,l1+l2+1);
1976   tb[0]=al;
1977
1978
1979   //Penalties: max score is NORM_F
1980   //Penalties must be negative
1981   igop=CL->gop;
1982   gep=igep=CL->gep;
1983
1984   if (n>max_size)
1985     {
1986       max_size=n;
1987
1988       vfree (MI);vfree (MJ); vfree (MM);
1989       free_int (slist, -1);
1990
1991       slist=declare_int (n,3);
1992
1993       MI=(long int*)vcalloc (5*n, sizeof (long));
1994       MJ=(long int*)vcalloc (5*n, sizeof (long));
1995       MM=(long int*)vcalloc (5*n, sizeof (long));
1996
1997     }
1998   else
1999     {
2000       for (a=0; a<n; a++)
2001         for (b=0; b<5; b++)LIN(MI,a,b)=LIN(MJ,a,b)=LIN(MJ,a,b)=-1000000;
2002     }
2003
2004   /*New Bit: Start*/
2005   if (!sortseq) sortseq=(int*)vcalloc( 7, sizeof (int));
2006   sortseq[0]=0; sortseq[1]=1;sortseq[2]=-1;
2007   sort_list_int2 (list, sortseq,7, 0, n-1);
2008
2009   for (a=0; a<n; a++)
2010     {
2011
2012       slist[a][0]=a;
2013       list[a][4]=a;
2014     }
2015
2016   sortseq[0]=1; sortseq[1]=0;sortseq[2]=-1;
2017   sort_list_int2 (list, sortseq,7, 0, n-1);
2018   for (a=0; a<n; a++)
2019   {
2020         slist[a][1]=list[a][4];
2021         list[a][5]=a;
2022   }
2023
2024   sortseq[0]=3; sortseq[1]=0;sortseq[2]=1;sortseq[3]=-1;
2025   sort_list_int2 (list, sortseq,7, 0, n-1);
2026   for (a=0; a<n; a++)
2027     {
2028       slist[a][2]=list[a][4];
2029       list[a][6]=a;
2030     }
2031
2032   sortseq[0]=0; sortseq[1]=1;sortseq[2]=-1;
2033   sort_list_int2 (list, sortseq,7, 0, n-1);
2034
2035   /*New Bit: EnD*/
2036
2037
2038
2039
2040
2041
2042   for (a=0; a<n; a++)
2043     {
2044
2045
2046       i=list[a][0];
2047       j=list[a][1];
2048
2049
2050       if (i==l1 || j==l2)gop=0;
2051       else gop=igop;
2052
2053       if (i==l1 && j==l2)start_trace=a;
2054       else if ( i==0 || j==0)
2055         {
2056           LIN(MM,a,0)=-1000000;
2057           if (j==0)
2058             {
2059
2060               LIN(MJ,a,0)=-10000000;
2061               LIN(MI,a,0)=gep*i;
2062
2063             }
2064           else if (i==0)
2065             {
2066
2067               LIN(MI,a,0)=-10000000;
2068               LIN(MJ,a,0)=gep*j;
2069
2070             }
2071
2072           LIN(MI,a,1)=LIN(MJ,a,1)=-1;
2073           LIN(MI,a,2)=LIN(MJ,a,2)=i;
2074           LIN(MI,a,3)=LIN(MJ,a,3)=j;
2075           continue;
2076         }
2077
2078       pi=list[a][5];
2079       pi=slist[pi-1][1];
2080
2081       pj=list[a][4];
2082       pj=slist[pj-1][0];
2083
2084       ij=list[a][6];
2085       ij=slist[ij-1][2];
2086
2087
2088       ij=list[a][6];
2089       ij=slist[ij-1][2];
2090
2091
2092
2093
2094
2095       prev_i=list[pi][0];
2096       prev_j=list[pj][1];
2097
2098       delta_i=list[a][0]-list[pi][0];
2099       delta_j=list[a][1]-list[pj][1];
2100
2101       /*Linear Notation*/
2102       LIN(MI,a,0)=MAX(LIN(MI,pi,0),(LIN(MM,pi,0)+gop))+delta_i*gep;
2103       LIN(MI,a,1)=pi;
2104       LIN(MI,a,2)=delta_i;
2105       LIN(MI,a,3)=0;
2106       LIN(MI,a,4)=(LIN(MI,pi,0)>=(LIN(MM,pi,0)+gop))?'i':'m';
2107
2108
2109       LIN(MJ,a,0)=MAX(LIN(MJ,pj,0),(LIN(MM,pj,0)+gop))+delta_j*gep;
2110       LIN(MJ,a,1)=pj;
2111       LIN(MJ,a,2)=0;
2112       LIN(MJ,a,3)=delta_j;
2113
2114       LIN(MJ,a,4)=(LIN(MJ,pj,0)>=LIN(MM,pj,0)+gop)?'j':'m';
2115
2116
2117
2118       if (a>1 && (ls=list[a][0]-list[ij][0])==(list[a][1]-list[ij][1]))
2119         {
2120           LIN(MM,a,0)=MAX3(LIN(MM,ij,0),LIN(MI,ij,0),LIN(MJ,ij,0))+list[a][2]-(ls*CL->nomatch);
2121
2122           LIN(MM,a,1)=ij;
2123           LIN(MM,a,2)=ls;
2124           LIN(MM,a,3)=ls;
2125           if ( LIN(MM,ij,0)>=LIN(MI,ij,0) && LIN(MM,ij,0)>=LIN(MJ,ij,0))LIN(MM,a,4)='m';
2126           else if ( LIN(MI,ij,0) >= LIN(MJ,ij,0))LIN(MM,a,4)='i';
2127           else LIN(MM,a,4)='j';
2128
2129         }
2130       else
2131         {
2132           LIN(MM,a,0)=UNDEFINED;
2133           LIN(MM,a,1)=-1;
2134         }
2135     }
2136
2137   a=start_trace;
2138   if (LIN(MM,a,0)>=LIN(MI,a,0) && LIN(MM,a,0) >=LIN(MJ,a,0))MT2=MM;
2139   else if ( LIN(MI,a,0)>=LIN(MJ,a,0))MT2=MI;
2140   else MT2=MJ;
2141
2142   score=MAX3(LIN(MM,a,0), LIN(MI,a,0), LIN(MJ,a,0));
2143
2144   i=l1;
2145   j=l2;
2146
2147
2148   while (!(i==0 &&j==0))
2149     {
2150       int next_a;
2151       l=MAX(LIN(MT2,a,2),LIN(MT2,a,3));
2152       // HERE ("%c from %c %d %d SCORE=%d [%d %d] [%2d %2d]", T2[a][5],T2[a][4], T2[a][2], T2[a][3], T2[a][0], gop, gep, i, j);
2153       if (i==0)
2154         {
2155           while ( j>0)
2156             {
2157               al[0][LEN]=0;
2158               al[1][LEN]=1;
2159               j--; LEN++;
2160             }
2161         }
2162       else if (j==0)
2163         {
2164           while ( i>0)
2165             {
2166               al[0][LEN]=1;
2167               al[1][LEN]=0;
2168               i--; LEN++;
2169             }
2170         }
2171
2172       else if (l==0) {HERE ("L=0 i=%d j=%d",l, i, j);exit (0);}
2173       else
2174         {
2175           for (b=0; b<l; b++, LEN++)
2176             {
2177               if (LIN(MT2,a,2)){al[0][LEN]=1;i--;ni++;}
2178               else al[0][LEN]=0;
2179
2180               if (LIN(MT2,a,3)){al[1][LEN]=1;j--;nj++;}
2181               else al[1][LEN]=0;
2182             }
2183
2184           next_a=LIN(MT2,a,1);
2185           if (LIN(MT2,a,4)=='m')MT2=MM;
2186           else if (LIN(MT2,a,4)=='i')MT2=MI;
2187           else if (LIN(MT2,a,4)=='j')MT2=MJ;
2188           a=next_a;
2189         }
2190     }
2191
2192
2193
2194   invert_list_char ( al[0], LEN);
2195   invert_list_char ( al[1], LEN);
2196
2197
2198   if ( A->declared_len<=LEN)A=realloc_aln  ( A,2*LEN+1);
2199   aln=A->seq_al;
2200   char_buf=(char*)vcalloc (LEN+1, sizeof (char));
2201
2202   for ( c=0; c< 2; c++)
2203     {
2204       for ( a=0; a< ns[c]; a++)
2205         {
2206           int ch=0;
2207           for ( b=0; b< LEN; b++)
2208             {
2209               if (al[c][b]==1)
2210                 char_buf[b]=aln[l_s[c][a]][ch++];
2211               else
2212                 char_buf[b]='-';
2213             }
2214           char_buf[b]='\0';
2215           sprintf (aln[l_s[c][a]],"%s", char_buf);
2216         }
2217     }
2218
2219   A->len_aln=LEN;
2220   A->nseq=ns[0]+ns[1];
2221
2222   vfree (char_buf);
2223   len[0]=LEN;
2224   return score;
2225 }
2226
2227
2228
2229
2230 //linked_pair_wise_collapse
2231 //Collapses the CL as it proceeds during the progressive alignment
2232 //Cannot be parralelized
2233
2234
2235 void display_ns (Alignment *A,Constraint_list *CL, int *ns, int **ls, char *txt);
2236 int cl2pair_list_collapse ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
2237 Constraint_list* collapse_list (Alignment *A,int *ns, int **ls, char**al, int len, Constraint_list *CL);
2238 int ns2s (int *ns, int **ls, int *is1, int *is2, int *is);
2239 int linked_pair_wise_collapse ( Alignment *A, int *ns, int **ls, Constraint_list *CL)
2240 {
2241   int n=0;
2242   static int **list=NULL;
2243   int score, a;
2244   char **al;
2245   int len=0;
2246
2247
2248   if ( !A)free_int (list, -1);
2249   if ( !CL->residue_index)return myers_miller_pair_wise (A, ns,ls,CL);
2250
2251   /*Prepare the list*/
2252
2253   cl2pair_list_collapse (A, ns, ls, CL, &list, &n);
2254   cl2diag_cap     (A, ns, ls, CL, &list, &n);
2255   cl2list_borders (A, ns, ls, CL, &list, &n);
2256   list2nodup_list (A, ns, ls, CL, &list, &n);
2257
2258   /*Do the DP*/
2259
2260   score=list2linked_pair_wise (A, ns, ls, CL, list, n, &al,&len);
2261   CL=collapse_list (A,ns, ls, al, len, CL);
2262   free_char (al, -1);
2263
2264   /*Free the list*/
2265   return score;
2266 }
2267
2268
2269 Constraint_list* collapse_list (Alignment *A,int *ns, int **ls, char **al, int len, Constraint_list *CL)
2270 {
2271   int s1, s2,s, cs1, cs2, cr1, cr2,l,ll;
2272   int **lu;
2273   int a,b,c,d;
2274   static char *add;
2275   static int  *p;
2276   FILE *fp;
2277
2278   if (!add)
2279     {
2280       add=vtmpnam (NULL);
2281       p=(int*)vcalloc ( 100, sizeof (int));
2282     }
2283
2284   lu=declare_int (2, len+1);
2285   for (a=0; a<2; a++)
2286     for (c=0,b=0; b<len; b++)if (al[a][b]){lu[a][++c]=b+1;}
2287
2288
2289   ns2s (ns, ls, &s1, &s2, &s);
2290
2291
2292   s1=name_is_in_list (A->name[s1], (CL->S)->name, (CL->S)->nseq, 100);
2293   s2=name_is_in_list (A->name[s2], (CL->S)->name, (CL->S)->nseq, 100);
2294   s =name_is_in_list (A->name[s ], (CL->S)->name, (CL->S)->nseq, 100);
2295
2296
2297   CL->residue_index[s]=(int**)vrealloc (CL->residue_index[s], (len+2)*sizeof (int*));
2298   for (a=0; a<=len; a++)
2299     {
2300       if (!CL->residue_index[s][a])
2301         {
2302           CL->residue_index[s][a]=(int*)vcalloc (1, sizeof (int));
2303           CL->residue_index[s][a][0]=1;
2304         }
2305     }
2306
2307   fp=vfopen (add, "w");
2308   CL->ne=0;
2309   for (cs1=0; cs1<(CL->S)->nseq; cs1++)
2310     {
2311       cr1=1;
2312       while (CL->residue_index[cs1][cr1])
2313         {
2314           for (ll=l=1; l<CL->residue_index[cs1][cr1][0]; l+=ICHUNK)
2315             {
2316               cs2=CL->residue_index[cs1][cr1][l+SEQ2];
2317
2318               if (cs1==s1 || cs1==s2 || cs2==s1 || cs2==s2)
2319                 {
2320                   p[SEQ1]=cs1;
2321                   p[SEQ2]=CL->residue_index[cs1][cr1][l+SEQ2];
2322                   p[R1]  =cr1;
2323                   p[R2]  =CL->residue_index[cs1][cr1][l+R2];
2324                   p[CONS]=CL->residue_index[cs1][cr1][l+CONS];
2325                   p[MISC]=CL->residue_index[cs1][cr1][l+MISC];
2326                   p[WE]=CL->residue_index[cs1][cr1][l+WE];
2327                   if (cs1==s1)
2328                     {
2329                       p[SEQ1]=s;
2330                       p[R1]=lu[0][p[R1]];
2331                     }
2332                   else if (cs1==s2)
2333                     {
2334                       p[SEQ1]=s;
2335                       p[R1]=lu[1][p[R1]];
2336                     }
2337
2338                   if (cs2==s1)
2339                     {
2340                       p[SEQ2]=s;
2341                       p[R2]=lu[0][p[R2]];
2342                     }
2343                   else if (cs2==s2)
2344                     {
2345
2346                       p[SEQ2]=s;
2347                       p[R2]=lu[1][p[R2]];
2348                     }
2349
2350                   if (p[SEQ1]==p[SEQ2]);
2351                   else for (d=0; d<CL->entry_len; d++)fprintf (fp, "%d ", p[d]);
2352                 }
2353               else
2354                 {
2355                   for (d=0; d<ICHUNK; d++) CL->residue_index[cs1][cr1][ll++]=CL->residue_index[cs1][cr1][d+l];
2356                   CL->ne++;
2357                 }
2358             }
2359           CL->residue_index[cs1][cr1][0]=ll;
2360           cr1++;
2361         }
2362     }
2363   vfclose (fp);
2364   CL=undump_constraint_list (CL,add);
2365   return CL;
2366 }
2367
2368
2369
2370
2371 int cl2pair_list_collapse ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in)
2372 {
2373   int si, r1,r2,t_s, t_r,t_w, t_s2, t_r2, t_w2, s1, s2;
2374   int a, b, l1, l2;
2375
2376   int nused;
2377   int *used_list;
2378
2379   float nscore, score, tot, filter, avg=0;
2380   float **used;
2381   int *norm;
2382
2383
2384
2385   if ( !A) return 0;
2386
2387   ns2s (ns, ls, &s1, &s2,NULL);
2388
2389   l1=strlen (A->seq_al[s1]);
2390   l2=strlen (A->seq_al[s2]);
2391   used=declare_float (l2+1,2);  used_list=(int*)vcalloc (l2+1, sizeof (int));
2392   nused=0;
2393   norm=(int*)vcalloc (l1+2, sizeof(int));
2394
2395   s1=name_is_in_list (A->name[s1], (CL->S)->name, (CL->S)->nseq, 100);
2396   s2=name_is_in_list (A->name[s2], (CL->S)->name, (CL->S)->nseq, 100);
2397
2398   for (r1=1; r1<=l1; r1++)
2399     {
2400       tot=0; nused=0;
2401       for (a=1; r1>0 && a<CL->residue_index[s1][r1][0];a+=ICHUNK)
2402         {
2403           t_s=CL->residue_index[s1][r1][a+SEQ2];
2404           t_r=CL->residue_index[s1][r1][a+R2];
2405           t_w=CL->residue_index[s1][r1][a+WE];
2406           norm[r1]++;
2407           for (b=0; b<CL->residue_index[t_s][t_r][0];)
2408             {
2409               if (b==0){t_s2=t_s;t_r2=t_r;t_w2=t_w;b++;}
2410               else
2411                 {
2412                   t_s2=CL->residue_index[t_s][t_r][b+SEQ2];
2413                   t_r2=CL->residue_index[t_s][t_r][b+R2];
2414                   t_w2=CL->residue_index[t_s][t_r][b+WE];
2415                   b+=ICHUNK;
2416                 }
2417
2418               if (t_s2==s2)
2419                 {
2420                   score=MIN(((float)t_w/(float)NORM_F),((float)t_w2/(float)NORM_F));
2421
2422                   if (!used[t_r2][1] && score>0)
2423                     {
2424                       used_list[nused++]=t_r2;
2425                     }
2426
2427                   tot+=score;
2428                   used[t_r2][0]+=score;
2429                   used[t_r2][1]++;
2430                 }
2431             }
2432         }
2433
2434       //FILTER: Keep in the graph the edges where (p1->p2/(Sum (P1->x))>0.01
2435       filter=0.01;
2436
2437       for (a=0; a<nused; a++)
2438         {
2439
2440           r2=used_list[a];
2441           nscore=used[r2][0]/tot; //Normalized score used for filtering
2442           score =used[r2][0];
2443
2444           used[r2][0]=used[r2][1]=0;
2445
2446           if (nscore>filter && r1!=0 && r2!=0 && r1!=l1 && r2!=l2)
2447             {
2448               score=((norm[r1]>0)?score/norm[r1]:0)*NORM_F;
2449               addE (r1,r2,((l1-(r1))+(r2)),score,list_in, n_in);
2450             }
2451         }
2452     }
2453
2454   free_float (used, -1);
2455   vfree (used_list);
2456   vfree (norm);
2457   return n_in[0];
2458 }
2459 int ns2s (int *ns, int **ls, int *is1, int *is2, int *is)
2460 {
2461   int a, b;
2462   int s1, s2, s;
2463
2464   s1=s2=s=-1;
2465
2466   for (a=0; a< 2; a++)
2467     for (b=0; b<ns[a]; b++)
2468       {
2469         if (a==0)s1=MAX(s1,(ls[a][b]));
2470         if (a==1)s2=MAX(s2,(ls[a][b]));
2471       }
2472   s=MAX((s1),(s2));
2473   if (is1)is1[0]=s1;
2474   if (is2)is2[0]=s2;
2475   if (is) is [0]=s;
2476
2477   return s;
2478 }
2479
2480 /////////////////////////////////////////////////////////////////////////////////////////////////////////
2481 ///////////////                                                                           ///////////////
2482 ///////////////                                                                           ///////////////
2483 ///////////////                       OLD linked Pw                                       ///////////////
2484 ///////////////                                                                           ///////////////
2485 ///////////////                                                                           ///////////////
2486 /////////////////////////////////////////////////////////////////////////////////////////////////////////
2487
2488
2489 int cl2pair_list390 ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
2490 int cl2list_borders390   (Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
2491 int cl2diag_cap390 (Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
2492 int list2nodup_list390 ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
2493 int list2linked_pair_wise390( Alignment *A, int *ns, int **l_s, Constraint_list *CL, int **list, int n);
2494
2495 void sort_list_int2_390 ( int **V,int *list,int N_F, int left, int right);
2496 int cmp_list_int2_390 (const int**a, const int**b);
2497 int cmp_list_int390 (const int**a, const int**b);
2498 void sort_list_int390 ( int **V,int N_F, int F, int left, int right);
2499 void sort_int390 ( int **V,int N_F, int F, int left, int right);
2500 void qsort390(void *base, size_t nmemb, size_t size, int (*compar)(const void *, const void *));
2501
2502 int procoffee_pair_wise ( Alignment *A, int *nsi, int **lsi, Constraint_list *CL)
2503 {
2504   int n=0;
2505   int **list=NULL;
2506   int score, a;
2507   int *ns, **ls;
2508   char **al;
2509   int len=0;
2510   if ( !CL->residue_index)return myers_miller_pair_wise (A, nsi,lsi,CL);
2511
2512
2513   ns=(int*)vcalloc (2, sizeof (int));
2514   ns[0]=nsi[1]; ns[1]=nsi[0];
2515
2516
2517   ls=declare_int (2, ns[0]+ns[1]);
2518   for (a=0; a<ns[1]; a++)
2519     ls[1][a]=lsi[0][a];
2520   for (a=0; a<ns[0]; a++)
2521     ls[0][a]=lsi[1][a];
2522
2523   cl2list_borders(A, ns, ls, CL, &list, &n);
2524   cl2pair_list_ext    (A, ns, ls, CL, &list, &n);
2525
2526   cl2diag_cap390    (A, ns, ls, CL, &list, &n);
2527
2528   sort_list_int     (list,7, 1, 0, n-1);
2529   list2nodup_list (A, ns, ls, CL, &list, &n);
2530   score=list2linked_pair_wise (A, ns, ls, CL, list, n, &al,&len);
2531   free_int (list,-1);
2532   return score;
2533 }
2534
2535 int cl2diag_cap390 (Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in)
2536 {
2537   int **list;
2538   int n, in, a, b, al1, al2;
2539   int max_n;
2540   int cap=0;
2541
2542   if (!A) return 0;
2543
2544   al1=strlen (A->seq_al[ls[0][0]]);
2545   al2=strlen (A->seq_al[ls[1][0]]);
2546
2547   list=list_in[0];
2548   n=n_in[0];
2549   max_n=read_array_size (list, sizeof (int*));
2550
2551
2552
2553
2554   for (a=0; a< n; a++)
2555     {
2556       b=list[a][3];
2557       list[a][3]=list[a][0];
2558       list[a][0]=b;
2559
2560     }
2561   sort_list_int (list, 4, 1, 0, n-1);
2562   for (a=0; a< n; a++)
2563     {
2564       b=list[a][3];
2565       list[a][3]=list[a][0];
2566       list[a][0]=b;
2567     }
2568
2569
2570   in=n;
2571
2572   for (a=0; a<in; a++)
2573     {
2574       int i, j, pi, pj, ni, nj;
2575       if (list[a][2]==0)continue;
2576       i=list[a][0];
2577       j=list[a][1];
2578
2579       if (a==0){pi=-10;pj=-10;}
2580       else {pi=list[a-1][0];pj=list[a-1][1];}
2581
2582       if (a==in-1){ni=-10; nj=-10;}
2583       else {ni=list[a+1][0]; nj=list[a+1][1];}
2584
2585
2586       if (i==0 || j==0);
2587       else if ( i==pi || j==pj);
2588       else if ( i-pi!=1 || j-pj!=1)
2589         {
2590
2591
2592           int x;
2593           int delta=MAX((i-pi),(j-pj));
2594           for (x=1; x<=delta; x++)
2595             {
2596               if (n==max_n){max_n+=1000;list=(int**)vrealloc (list, max_n*sizeof (int*));}
2597               if (!list[n])list[n]=(int*)vcalloc (7, sizeof (int));
2598               if ((i-x)<=1)continue;
2599               if ((j-x)<=1)continue;
2600               list[n][0]=i-x;
2601               list[n][1]=j-x;
2602               list[n][3]=list[a][3];
2603               list[n][2]=cap;
2604               n++;
2605             }
2606         }
2607
2608
2609       if (i==al1 || j==al2);
2610       else if ( i==ni || j==nj);
2611       else if ( ni-i!=1 || nj-j!=1)
2612         {
2613           int x;
2614           int delta=MAX((ni-i),(nj-j));
2615           for (x=1; x<=delta; x++)
2616             {
2617                 if (n==max_n){max_n+=1000;list=(int**)vrealloc (list, max_n*sizeof (int*));}
2618                 if (!list[n])list[n]=(int*)vcalloc (7, sizeof (int));
2619
2620                 if (i+x>=al1)continue;
2621                 if (j+x>=al2) continue;
2622
2623
2624                 list[n][0]=i+x;
2625                 list[n][1]=j+x;
2626                 list[n][3]=list[a][3];
2627                 list[n][2]=cap;
2628                 n++;
2629               }
2630
2631         }
2632
2633     }
2634   list_in[0]=list;
2635   n_in[0]=n;
2636   return n;
2637 }
2638
2639 /////////////////////////////////////////////////////////////////////////////////////////////////////////
2640 ///////////////                                                                           ///////////////
2641 ///////////////                                                                           ///////////////
2642 ///////////////                       Profile_PW                                          ///////////////
2643 ///////////////                                                                           ///////////////
2644 ///////////////                                                                           ///////////////
2645 /////////////////////////////////////////////////////////////////////////////////////////////////////////
2646
2647 int **aln2prf (Alignment *A, char *matrix)
2648 {
2649   int **mat=read_matrice (matrix);
2650   int *vector;
2651   int **prf;
2652   int r,a,b,c;
2653   vector=(int*)vcalloc (27, sizeof (int));
2654   prf=declare_int (A->len_aln, 26);
2655
2656   for (a=0; a<A->len_aln; a++)
2657     {
2658       for (b=0; b<26; b++)vector[b]=0;
2659       for ( b=0; b<A->nseq; b++)
2660         {
2661           r=tolower(A->seq_al[b][a]);
2662           if (r=='-')continue;
2663           r-='a';
2664           vector[r]++;
2665         }
2666
2667       for (b=0; b<26; b++)
2668         {
2669           double score=0;
2670           double tot=0;
2671           for (c=0; c<26; c++)
2672             {
2673
2674               score+=mat[b][c]*vector[c];
2675               tot  +=vector[c];
2676             }
2677           prf[a][b]=(int)((score/tot)*100);
2678         }
2679     }
2680   vfree (vector);
2681
2682   return prf;
2683 }
2684 char *add_sequence2prf (char *seq,char *al, int **prf, int lj, int gop, int gep)
2685 {
2686   static int **match, **score,**indel,**tb;
2687   int m,g;
2688   int li, i, j,l;
2689   static int max_i;
2690   static int max_j;
2691   gop*=100;
2692   gep*=100;
2693   if (seq)li=strlen (seq);
2694
2695   if (li>max_i ||lj>max_j|| !seq)
2696     {
2697       if (match)free_int (match, -1);
2698       if (score)free_int (score, -1);
2699       if (indel)free_int (indel, -1);
2700       if (!seq) return NULL;
2701
2702       match=declare_int( li+1, lj+1);
2703       score=declare_int( li+1, lj+1);
2704       indel=declare_int( li+1, lj+1);
2705       tb=declare_int( li+1, lj+1);
2706
2707       max_i=li;
2708       max_j=lj;
2709     }
2710
2711   if (!al)al=(char*)vcalloc ( li+lj+1, sizeof (char));
2712   lower_string (seq);
2713   for (i=1; i<=li; i++)//seq
2714     {
2715       match[i][0]=2*gop;
2716       indel[i][0]=gop;
2717       for (j=1; j<=lj; j++)//prf
2718         {
2719           g=indel[i][j]=MAX((match[i][j-1]+gop),(indel[i][j-1]+gep));
2720
2721           m=match[i][j]=score[i-1][j-1]+prf[j-1][seq[i-1]-'a'];
2722           if (m>g){score[i][j]=m;tb[i][j]=1;}
2723           else    {score[i][j]=g;tb[i][j]=0;}
2724         }
2725     }
2726
2727   l=0;
2728   i=li;
2729   j=lj;
2730   while (l<lj)
2731     {
2732
2733       //HERE ("%d %d -- %d",i,j, tb[i][j]);
2734       if (tb[i][j]==1){al[l]=seq[i-1];i--;j--;}
2735       else {j--;al[l]='-';}
2736       l++;
2737     }
2738   al[l]='\0';
2739
2740   invert_string2 (al);
2741   return al;
2742 }