website/archive/binaries/mac/src/tcoffee/t_coffee_source/util_dp_gotoh_nw.c

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