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