binaries/src/tcoffee/t_coffee_source/util_dp_sim.c

   1 #include <stdio.h>
   2 #include <math.h>
   3 #include <stdarg.h>
   4 #include <string.h>
   5
   6 #include "io_lib_header.h"
   7 #include "util_lib_header.h"
   8 #include "define_header.h"
   9 #include "dp_lib_header.h"
  10 /* extern char name0[], name1[]; */
  11 /* extern int match, mismh; */
  12
  13
  14
  15 static Constraint_list *CL;
  16 static int * ns;
  17 static int **l_s;
  18 static Alignment *Aln;
  19 static int **pos;
  20 static int *seqc0, *seqc1;
  21 static int min0,min1,max0,max1,mins;
  22
  23 static void* sim_vcalloc( size_t nobj, size_t size);
  24 static void sim_free_all ();
  25 static int sim_reset_static_variable ();
  26 static int big_pass(int *A,int *B,int M,int N,int K, int nseq) ;
  27 static int locate(int *A,int *B,int nseq);
  28 static int small_pass(int *A,int *B,int count,int nseq);
  29 static int no_cross ();
  30 static int diff_sim( int *A,int *B,int M,int N,int tb,int te);
  31 int calcons(int *aa0,int n0,int *aa1,int n1,int *res,int *nc,int *nident, Alignment *A, int *ns, int **l_s, Constraint_list *CL);
  32
  33
  34
  35 #define SIM_GAP -1
  36 #define min(x,y) ((x)<=(y) ? (x) : (y))
  37 //#define TC_SCORE_SIM(x,y) TC_SCORE (x,y)
  38
  39 static int q, r;                        /* gap penalties */
  40 static int qr;                          /* qr = q + r */
  41
  42
  43
  44 typedef struct ONE { int COL ;  struct ONE  *NEXT ;} pair, *pairptr;
  45 pairptr *row, z, z1;                    /* for saving used aligned pairs */
  46
  47
  48 #define PAIRNULL (pairptr)NULL
  49 static int tt;
  50
  51 typedef struct SIM_NODE
  52         { int  SIM_SCORE;
  53           int  SIM_STARI;
  54           int  SIM_STARJ;
  55           int  SIM_ENDI;
  56           int  SIM_ENDJ;
  57           int  SIM_TOP;
  58           int  SIM_BOT;
  59           int  SIM_LEFT;
  60           int  SIM_RIGHT; }  vertex,
  61 #ifdef FAR_PTR
  62  far *vertexptr;
  63 #else
  64      *vertexptr;
  65 #endif
  66
  67 vertexptr  *LIST;                       /* an array for saving k best scores */
  68 vertexptr  low = 0;                     /* lowest score node in LIST */
  69 vertexptr  most = 0;                    /* latestly accessed node in LIST */
  70 static int numnode;                     /* the number of nodes in LIST */
  71
  72 static int *CC, *DD;                    /* saving matrix scores */
  73 static int *RR, *SS, *EE, *FF;          /* saving start-points */
  74 static int *HH, *WW;                    /* saving matrix scores */
  75 static int *II, *JJ, *XX, *YY;          /* saving start-points */
  76 static int  m1, mm, n1, nn;             /* boundaries of recomputed area */
  77 static int  rl, cl;                     /* left and top boundaries */
  78 static int  lmin;                       /* minimum score in LIST */
  79 static int flag;                        /* indicate if recomputation necessary*/
  80
  81 /* DIAG() assigns value to x if (ii,jj) is never used before */
  82 #define DIAG(ii, jj, x, value)                          \
  83 { for ( tt = 1, z = row[(ii)]; z != PAIRNULL; z = z->NEXT )     \
  84     if ( z->COL == (jj) )                               \
  85       { tt = 0; break; }                                \
  86   if ( tt )                                             \
  87     x = ( value );                                      \
  88 }
  89
  90 /* replace (ss1, xx1, yy1) by (ss2, xx2, yy2) if the latter is large */
  91 #define ORDER(ss1, xx1, yy1, ss2, xx2, yy2)             \
  92 { if ( ss1 < ss2 )                                      \
  93     { ss1 = ss2; xx1 = xx2; yy1 = yy2; }                \
  94   else                                                  \
  95     if ( ss1 == ss2 )                                   \
  96       { if ( xx1 < xx2 )                                \
  97           { xx1 = xx2; yy1 = yy2; }                     \
  98         else                                            \
  99           if ( xx1 == xx2 && yy1 < yy2 )                \
 100             yy1 = yy2;                                  \
 101       }                                                 \
 102 }
 103
 104 /* The following definitions are for function diff() */
 105
 106 static int  zero = 0;                           /* int type zero        */
 107 #define gap(k)  ((k) <= 0 ? 0 : q+r*(k))        /* k-symbol indel score */
 108
 109 static int *sapp;                               /* Current script append ptr */
 110 static int  last;                               /* Last script op appended */
 111
 112 static int I, J;                                /* current positions of A ,B */
 113 static int no_mat;                              /* number of matches */
 114 static int no_mis;                              /* number of mismatches */
 115 static int al_len;                              /* length of alignment */
 116                                                 /* Append "Delete k" op */
 117 #define DEL(k) \
 118 { I += k;\
 119   al_len += k;\
 120   if (last < 0)\
 121     last = sapp[-1] -= (k);\
 122   else\
 123     last = *sapp++ = -(k);\
 124 }
 125                                                 /* Append "Insert k" op */
 126 #define INS(k) \
 127 { J += k;\
 128   al_len += k;\
 129   if (last < 0)\
 130     { sapp[-1] = (k); *sapp++ = last; } \
 131   else\
 132     last = *sapp++ = (k);\
 133 }
 134
 135                                                 /* Append "Replace" op */
 136 #define REP \
 137 { last = *sapp++ = 0;\
 138   al_len += 1;\
 139 }
 140
 141
 142 /*
 143 int sim_pair_wise_lalign (Alignment *in_A, int *in_ns, int **in_l_s,Constraint_list *in_CL)
 144 {
 145   if ( in_ns[0]==1 && in_ns[1]==1)
 146     return sim_pair_wise_lalign (in_A, in_ns, in_l_s,in_CL);
 147   else
 148   */
 149
 150
 151
 152
 153
 154 int sim_pair_wise_lalign (Alignment *in_A, int *in_ns, int **in_l_s,Constraint_list *in_CL)
 155 /* SIM(A,B,M,N,K,V,Q,R) reports K best non-intersecting alignments of
 156    the segments of A and B in order of similarity scores, where
 157    V[a][b] is the score of aligning a and b, and -(Q+R*i) is the score
 158    of an i-symbol indel.
 159 */
 160 {
 161   int endi, endj, stari, starj; /* endpoint and startpoint */
 162   int  score;                           /* the max score in LIST */
 163   int count;                            /* maximum size of list */
 164   int i;
 165   int  *S;                              /* saving operations for diff */
 166   int nc, nident;               /* for display */
 167   vertexptr cur;                        /* temporary pointer */
 168   vertexptr findmax();                  /* return the largest score node */
 169   double percent;
 170   int t1, t2, g1, g2, r1, r2;
 171   int a, b, c, d, e;
 172 /*cedric was here 11/2/99*/
 173   int CEDRIC_MAX_N_ALN=999;
 174   int CEDRIC_THRESHOLD=50;
 175   int *A, *B;
 176   int M, N, K, maxl;
 177   int nseq;
 178   int R, Q;
 179   Alignment *DA;
 180
 181
 182   DA=in_A;
 183
 184   Aln=copy_aln (in_A, NULL);
 185
 186
 187
 188   l_s=in_l_s;
 189   ns=in_ns;
 190   CL=in_CL;
 191   K=CL->lalign_n_top;
 192
 193   M=strlen (Aln->seq_al[l_s[0][0]]);
 194   N=strlen (Aln->seq_al[l_s[1][0]]);
 195   maxl=M+N+1;
 196
 197   pos=aln2pos_simple (Aln,-1, ns, l_s);
 198
 199   seqc0=(int*)sim_vcalloc (maxl,sizeof (int));
 200   A=(int*)sim_vcalloc (maxl,sizeof (int));
 201   for ( a=0; a<maxl; a++){seqc0[a]=A[a]=a;}
 202   A[M+1]='\0';
 203
 204   seqc1=(int*)sim_vcalloc (maxl,sizeof (int));
 205   B=(int*)sim_vcalloc (maxl,sizeof (int));
 206   for ( a=0; a<maxl; a++){seqc1[a]=B[a]=a;}
 207   B[N+1]='\0';
 208
 209   nseq=(l_s[0][0]!=l_s[1][0])?2:1;
 210
 211
 212   Q=MAX(CL->gop, -CL->gop)*SCORE_K;
 213   R=MAX(CL->gep, -CL->gep)*SCORE_K;
 214
 215
 216
 217   if ( K==CEDRIC_MAX_N_ALN)K--;
 218   else if ( K<0)
 219       {
 220
 221        CEDRIC_THRESHOLD=-K;
 222        K=CEDRIC_MAX_N_ALN;
 223       }
 224
 225   /* allocate space for all vectors */
 226
 227   CC = ( int * ) sim_vcalloc(N+1, sizeof(int));
 228   DD = ( int * ) sim_vcalloc(N+1, sizeof(int));
 229   RR = ( int * ) sim_vcalloc(N+1, sizeof(int));
 230   SS = ( int * ) sim_vcalloc(N+1, sizeof(int));
 231   EE = ( int * ) sim_vcalloc(N+1, sizeof(int));
 232   FF = ( int * ) sim_vcalloc(N+1, sizeof(int));
 233
 234   HH = ( int * ) sim_vcalloc(M + 1, sizeof(int));
 235   WW = ( int * ) sim_vcalloc(M + 1, sizeof(int));
 236   II = ( int * ) sim_vcalloc(M + 1, sizeof(int));
 237   JJ = ( int * ) sim_vcalloc(M + 1, sizeof(int));
 238   XX = ( int * ) sim_vcalloc(M + 1, sizeof(int));
 239   YY = ( int * ) sim_vcalloc(M + 1, sizeof(int));
 240   S = ( int * )  sim_vcalloc(min(M,N)*5/4+1, sizeof (int));
 241   row = ( pairptr * ) sim_vcalloc( (M + 1), sizeof(pairptr));
 242
 243
 244   /* set up list for each row */
 245   if (nseq == 2) for ( i = 1; i <= M; i++ ) row[i]= PAIRNULL;
 246   else {
 247           z = ( pairptr )sim_vcalloc (M,(int)sizeof(pair));
 248           for ( i = 1; i <= M; i++,z++) {
 249                   row[i] = z;
 250                   z->COL = i;
 251                   z->NEXT = PAIRNULL;
 252           }
 253   }
 254
 255
 256   q = Q;
 257   r = R;
 258   qr = q + r;
 259
 260   LIST = ( vertexptr * ) sim_vcalloc( K, sizeof(vertexptr));
 261   for ( i = 0; i < K ; i++ )
 262     LIST[i] = ( vertexptr )sim_vcalloc( 1, sizeof(vertex));
 263
 264
 265   numnode = lmin = 0;
 266   big_pass(A,B,M,N,K,nseq);
 267
 268
 269
 270   /* Report the K best alignments one by one. After each alignment is
 271      output, recompute part of the matrix. First determine the size
 272      of the area to be recomputed, then do the recomputation         */
 273
 274
 275   for ( count = K - 1; count >= 0; count-- )
 276     { if ( numnode == 0 )
 277         {
 278
 279           padd_aln (in_A);
 280           /*fatal("The number of alignments computed is too large");*/
 281           sim_free_all();
 282           return 1;
 283         }
 284
 285       cur = findmax();  /* Return a pointer to a node with max score*/
 286       score = cur->SIM_SCORE;
 287       if ( K==CEDRIC_MAX_N_ALN && score<CEDRIC_THRESHOLD)break;
 288       stari = ++cur->SIM_STARI;
 289       starj = ++cur->SIM_STARJ;
 290       endi = cur->SIM_ENDI;
 291       endj = cur->SIM_ENDJ;
 292       m1 = cur->SIM_TOP;
 293       mm = cur->SIM_BOT;
 294       n1 = cur->SIM_LEFT;
 295       nn = cur->SIM_RIGHT;
 296       rl = endi - stari + 1;
 297       cl = endj - starj + 1;
 298       I = stari - 1;
 299       J = starj - 1;
 300       sapp = S;
 301       last = 0;
 302       al_len = 0;
 303       no_mat = 0;
 304       no_mis = 0;
 305       diff_sim(&A[stari]-1, &B[starj]-1,rl,cl,q,q);
 306
 307
 308       min0 = stari;
 309       min1 = starj;
 310       max0 = stari+rl-1;
 311       max1 = starj+cl-1;
 312       calcons(A+1,M,B+1,N,S,&nc,&nident, Aln,ns, l_s, CL);
 313       percent = (double)nident*100.0/(double)nc;
 314
 315
 316
 317       /*Min0: index of the last residue before the first in a 1..N+1 numerotation*/
 318
 319
 320
 321       if (!DA->A)DA->A=copy_aln(Aln, DA->A);
 322       DA->A=realloc_alignment (DA->A,nc+1);
 323
 324
 325       DA=DA->A;
 326       DA->A=NULL;
 327
 328       for ( c=0; c< 2; c++)
 329             {
 330             for ( a=0; a< ns[c]; a++)
 331                 {
 332                   e=(c==0)?min0:min1;
 333                   for ( d=0; d<e; d++)
 334                     {
 335                       DA->order[l_s[c][a]][1]+=1-is_gap(Aln->seq_al[l_s[c][a]][d]);
 336                     }
 337                 }
 338             }
 339
 340
 341       for ( t1=min0,t2=min1,a=0; a<nc; a++)
 342         {
 343           r1=seqc0[a];
 344           r2=seqc1[a];
 345
 346           g1=(r1==SIM_GAP || r1>M)?0:1;
 347           g2=(r2==SIM_GAP || r2>N)?0:1;
 348           t1+=g1;
 349           t2+=g2;
 350           for (b=0; b<ns[0]; b++)DA->seq_al[l_s[0][b]][a]=(g1)?Aln->seq_al[l_s[0][b]][A[t1-1]]:'-';
 351           for (b=0; b<ns[1]; b++)DA->seq_al[l_s[1][b]][a]=(g2)?Aln->seq_al[l_s[1][b]][B[t2-1]]:'-';
 352         }
 353       for (b=0; b<ns[0]; b++){DA->seq_al[l_s[0][b]][a]='\0';}
 354       for (b=0; b<ns[1]; b++){DA->seq_al[l_s[1][b]][a]='\0';}
 355
 356       DA->nseq=ns[0]+ns[1];
 357       DA->len_aln=nc;
 358       DA->score=percent;
 359       DA->score_aln=score;
 360       fflush(stdout);
 361
 362
 363       if ( count )
 364         { flag = 0;
 365           locate(A,B,nseq);
 366           if ( flag )
 367             small_pass(A,B,count,nseq);
 368         }
 369     }
 370   padd_aln (in_A);
 371
 372   sim_free_all();
 373   free_int (pos, -1);
 374   free_aln (Aln);
 375
 376
 377
 378   return 1;
 379
 380 }
 381 int sim_reset_static_variable ()
 382 {
 383   CC=DD=RR=SS=EE=FF=HH=WW=II=JJ=XX=YY=sapp=NULL;
 384   min0=min1=max0=max1=mins=q=r=qr=tt=numnode=m1=n1=nn=rl=cl=lmin=flag=zero=last=I=J=no_mat=no_mis=al_len=0;
 385   most=low=NULL;/*Very important: cause a bug if not reset*/
 386   LIST=NULL;    /*Very important: cause a bug if not reset*/
 387   return 0;
 388 }
 389 /* A big pass to compute K best classes */
 390
 391
 392 int big_pass(int *A,int *B,int M,int N,int K, int nseq)
 393 { register  int  i, j;                  /* row and column indices */
 394   register  int  c;                     /* best score at current point */
 395   register  int  f;                     /* best score ending with insertion */
 396   register  int  d;                     /* best score ending with deletion */
 397   register  int  p;                     /* best score at (i-1, j-1) */
 398   register  int  ci, cj;                /* end-point associated with c */
 399   register  int  di, dj;                /* end-point associated with d */
 400   register  int  fi, fj;                /* end-point associated with f */
 401   register  int  pi, pj;                /* end-point associated with p */
 402
 403   int   addnode();                      /* function for inserting a node */
 404
 405
 406         /* Compute the matrix and save the top K best scores in LIST
 407            CC : the scores of the current row
 408            RR and EE : the starting point that leads to score CC
 409            DD : the scores of the current row, ending with deletion
 410            SS and FF : the starting point that leads to score DD        */
 411         /* Initialize the 0 th row */
 412         for ( j = 1; j <= N ; j++ )
 413           {  CC[j] = 0;
 414              RR[j] = 0;
 415              EE[j] = j;
 416              DD[j] = - (q);
 417              SS[j] = 0;
 418              FF[j] = j;
 419           }
 420         for ( i = 1; i <= M; i++)
 421           {  c = 0;                             /* Initialize column 0 */
 422              f = - (q);
 423              ci = fi = i;
 424              if ( nseq == 2 )
 425                { p = 0;
 426                  pi = i - 1;
 427                  cj = fj = pj = 0;
 428                }
 429              else
 430                { p = CC[i];
 431                  pi = RR[i];
 432                  pj = EE[i];
 433                  cj = fj = i;
 434                }
 435              for ( j = (nseq == 2 ? 1 : (i+1)) ; j <= N ; j++ )
 436                {  f = f - r;
 437                   c = c - qr;
 438                   ORDER(f, fi, fj, c, ci, cj)
 439                   c = CC[j] - qr;
 440                   ci = RR[j];
 441                   cj = EE[j];
 442                   d = DD[j] - r;
 443                   di = SS[j];
 444                   dj = FF[j];
 445                   ORDER(d, di, dj, c, ci, cj)
 446                   c = 0;
 447
 448                   DIAG(i, j, c, p+TC_SCORE(A[i-1],B[j-1]))              /* diagonal */
 449
 450                   if ( c <= 0 )
 451                     { c = 0; ci = i; cj = j; }
 452                   else
 453                     { ci = pi; cj = pj; }
 454                   ORDER(c, ci, cj, d, di, dj)
 455                   ORDER(c, ci, cj, f, fi, fj)
 456                   p = CC[j];
 457                   CC[j] = c;
 458                   pi = RR[j];
 459                   pj = EE[j];
 460                   RR[j] = ci;
 461                   EE[j] = cj;
 462                   DD[j] = d;
 463                   SS[j] = di;
 464                   FF[j] = dj;
 465                   if ( c > lmin )       /* add the score into list */
 466                     lmin = addnode(c, ci, cj, i, j, K, lmin);
 467                 }
 468           }
 469 return 1;
 470 }
 471
 472 /* Determine the left and top boundaries of the recomputed area */
 473
 474 int locate(int *A,int *B,int nseq)
 475 { register  int  i, j;                  /* row and column indices */
 476   register  int  c;                     /* best score at current point */
 477   register  int  f;                     /* best score ending with insertion */
 478   register  int  d;                     /* best score ending with deletion */
 479   register  int  p;                     /* best score at (i-1, j-1) */
 480   register  int  ci, cj;                /* end-point associated with c */
 481   register  int  di=0, dj=0;            /* end-point associated with d */
 482   register  int  fi, fj;                /* end-point associated with f */
 483   register  int  pi, pj;                /* end-point associated with p */
 484   int  cflag, rflag;                    /* for recomputation */
 485   int   addnode();                      /* function for inserting a node */
 486   int  limit;                           /* the bound on j */
 487
 488         /* Reverse pass
 489            rows
 490            CC : the scores on the current row
 491            RR and EE : the endpoints that lead to CC
 492            DD : the deletion scores
 493            SS and FF : the endpoints that lead to DD
 494
 495            columns
 496            HH : the scores on the current columns
 497            II and JJ : the endpoints that lead to HH
 498            WW : the deletion scores
 499            XX and YY : the endpoints that lead to WW
 500         */
 501         for ( j = nn; j >= n1 ; j-- )
 502           {  CC[j] = 0;
 503              EE[j] = j;
 504              DD[j] = - (q);
 505              FF[j] = j;
 506              if ( nseq == 2 || j > mm )
 507                 RR[j] = SS[j] = mm + 1;
 508              else
 509                 RR[j] = SS[j] = j;
 510           }
 511
 512         for ( i = mm; i >= m1; i-- )
 513           {  c = p = 0;
 514              f = - (q);
 515              ci = fi = i;
 516              pi = i + 1;
 517              cj = fj = pj = nn + 1;
 518
 519              if ( nseq == 2 || n1 > i )
 520                 limit = n1;
 521              else
 522                 limit = i + 1;
 523              for ( j = nn; j >= limit ; j-- )
 524                {  f = f - r;
 525                   c = c - qr;
 526                   ORDER(f, fi, fj, c, ci, cj)
 527                   c = CC[j] - qr;
 528                   ci = RR[j];
 529                   cj = EE[j];
 530                   d = DD[j] - r;
 531                   di = SS[j];
 532                   dj = FF[j];
 533                   ORDER(d, di, dj, c, ci, cj)
 534                   c = 0;
 535                   DIAG(i, j, c, p+TC_SCORE(A[i-1],B[j-1]))              /* diagonal */
 536
 537                   if ( c <= 0 )
 538                     { c = 0; ci = i; cj = j; }
 539                   else
 540                     { ci = pi; cj = pj; }
 541                   ORDER(c, ci, cj, d, di, dj)
 542                   ORDER(c, ci, cj, f, fi, fj)
 543                   p = CC[j];
 544                   CC[j] = c;
 545                   pi = RR[j];
 546                   pj = EE[j];
 547                   RR[j] = ci;
 548                   EE[j] = cj;
 549                   DD[j] = d;
 550                   SS[j] = di;
 551                   FF[j] = dj;
 552                   if ( c > lmin )
 553                     flag = 1;
 554                 }
 555              if ( nseq == 2 || i < n1 )
 556                { HH[i] = CC[n1];
 557                  II[i] = RR[n1];
 558                  JJ[i] = EE[n1];
 559                  WW[i] = DD[n1];
 560                  XX[i] = SS[n1];
 561                  YY[i] = FF[n1];
 562                }
 563           }
 564
 565   for ( rl = m1, cl = n1; ; )
 566     { for ( rflag = cflag = 1; ( rflag && m1 > 1 ) || ( cflag && n1 > 1 ) ;  )
 567         { if ( rflag && m1 > 1 )        /* Compute one row */
 568             { rflag = 0;
 569               m1--;
 570               c = p = 0;
 571               f = - (q);
 572               ci = fi = m1;
 573               pi = m1 + 1;
 574               cj = fj = pj = nn + 1;
 575
 576               for ( j = nn; j >= n1 ; j-- )
 577                 { f = f - r;
 578                   c = c - qr;
 579                   ORDER(f, fi, fj, c, ci, cj)
 580                   c = CC[j] - qr;
 581                   ci = RR[j];
 582                   cj = EE[j];
 583                   d = DD[j] - r;
 584                   di = SS[j];
 585                   dj = FF[j];
 586                   ORDER(d, di, dj, c, ci, cj)
 587                   c = 0;
 588                   DIAG(m1, j, c, TC_SCORE(A[m1-1],B[j-1]))              /* diagonal */
 589
 590                   if ( c <= 0 )
 591                     { c = 0; ci = m1; cj = j; }
 592                   else
 593                     { ci = pi; cj = pj; }
 594                   ORDER(c, ci, cj, d, di, dj)
 595                   ORDER(c, ci, cj, f, fi, fj)
 596                   p = CC[j];
 597                   CC[j] = c;
 598                   pi = RR[j];
 599                   pj = EE[j];
 600                   RR[j] = ci;
 601                   EE[j] = cj;
 602                   DD[j] = d;
 603                   SS[j] = di;
 604                   FF[j] = dj;
 605                   if ( c > lmin )
 606                      flag = 1;
 607                   if ( ! rflag && ( (ci > rl && cj > cl) || (di > rl && dj > cl)
 608                                     || (fi > rl && fj > cl) ) )
 609                       rflag = 1;
 610                 }
 611               HH[m1] = CC[n1];
 612               II[m1] = RR[n1];
 613               JJ[m1] = EE[n1];
 614               WW[m1] = DD[n1];
 615               XX[m1] = SS[n1];
 616               YY[m1] = FF[n1];
 617               if ( ! cflag && ( (ci > rl && cj > cl) || (di > rl && dj > cl)
 618                                 || (fi > rl && fj > cl )) )
 619                  cflag = 1;
 620             }
 621
 622           if ( nseq == 1 && n1 == (m1 + 1) && ! rflag )
 623              cflag = 0;
 624           if ( cflag && n1 > 1 )        /* Compute one column */
 625             { cflag = 0;
 626               n1--;
 627               c = 0;
 628               f = - (q);
 629               cj = fj = n1;
 630               if ( nseq == 2 || mm < n1 )
 631                 { p = 0;
 632                   ci = fi = pi = mm + 1;
 633                   pj = n1 + 1;
 634                   limit = mm;
 635                 }
 636               else
 637                 { p = HH[n1];
 638                   pi = II[n1];
 639                   pj = JJ[n1];
 640                   ci = fi = n1;
 641                   limit = n1 - 1;
 642                 }
 643               for ( i = limit; i >= m1 ; i-- )
 644                 { f = f - r;
 645                   c = c - qr;
 646                   ORDER(f, fi, fj, c, ci, cj)
 647                   c = HH[i] - qr;
 648                   ci = II[i];
 649                   cj = JJ[i];
 650                   d = WW[i] - r;
 651                   di = XX[i];
 652                   dj = YY[i];
 653                   ORDER(d, di, dj, c, ci, cj)
 654                   c = 0;
 655                   DIAG(i, n1, c, p+TC_SCORE(A[i-1], B[n1-1]))
 656
 657
 658
 659                   if ( c <= 0 )
 660                     { c = 0; ci = i; cj = n1; }
 661                   else
 662                     { ci = pi; cj = pj; }
 663                   ORDER(c, ci, cj, d, di, dj)
 664                   ORDER(c, ci, cj, f, fi, fj)
 665                   p = HH[i];
 666                   HH[i] = c;
 667                   pi = II[i];
 668                   pj = JJ[i];
 669                   II[i] = ci;
 670                   JJ[i] = cj;
 671                   WW[i] = d;
 672                   XX[i] = di;
 673                   YY[i] = dj;
 674                   if ( c > lmin )
 675                      flag = 1;
 676                   if ( ! cflag && ( (ci > rl && cj > cl) || (di > rl && dj > cl)
 677                                     || (fi > rl && fj > cl )) )
 678                      cflag = 1;
 679                 }
 680               CC[n1] = HH[m1];
 681               RR[n1] = II[m1];
 682               EE[n1] = JJ[m1];
 683               DD[n1] = WW[m1];
 684               SS[n1] = XX[m1];
 685               FF[n1] = YY[m1];
 686               if ( ! rflag && ( (ci > rl && cj > cl) || (di > rl && dj > cl)
 687                                 || (fi > rl && fj > cl) ) )
 688                  rflag = 1;
 689             }
 690         }
 691       if ( (m1 == 1 && n1 == 1) || no_cross() )
 692          break;
 693    }
 694   m1--;
 695   n1--;
 696 return 1;
 697 }
 698
 699 /* recompute the area on forward pass */
 700 int small_pass(int *A,int *B,int count,int nseq)
 701 { register  int  i, j;                  /* row and column indices */
 702   register  int  c;                     /* best score at current point */
 703   register  int  f;                     /* best score ending with insertion */
 704   register  int  d;                     /* best score ending with deletion */
 705   register  int  p;                     /* best score at (i-1, j-1) */
 706   register  int  ci, cj;                /* end-point associated with c */
 707   register  int  di, dj;                /* end-point associated with d */
 708   register  int  fi, fj;                /* end-point associated with f */
 709   register  int  pi, pj;                /* end-point associated with p */
 710   int   addnode();                      /* function for inserting a node */
 711   int  limit;                           /* lower bound on j */
 712
 713         for ( j = n1 + 1; j <= nn ; j++ )
 714           {  CC[j] = 0;
 715              RR[j] = m1;
 716              EE[j] = j;
 717              DD[j] = - (q);
 718              SS[j] = m1;
 719              FF[j] = j;
 720           }
 721         for ( i = m1 + 1; i <= mm; i++)
 722           {  c = 0;                             /* Initialize column 0 */
 723              f = - (q);
 724              ci = fi = i;
 725
 726              if ( nseq == 2 || i <= n1 )
 727                { p = 0;
 728                  pi = i - 1;
 729                  cj = fj = pj = n1;
 730                  limit = n1 + 1;
 731                }
 732              else
 733                { p = CC[i];
 734                  pi = RR[i];
 735                  pj = EE[i];
 736                  cj = fj = i;
 737                  limit = i + 1;
 738                }
 739              for ( j = limit ; j <= nn ; j++ )
 740                {  f = f - r;
 741                   c = c - qr;
 742                   ORDER(f, fi, fj, c, ci, cj)
 743                   c = CC[j] - qr;
 744                   ci = RR[j];
 745                   cj = EE[j];
 746                   d = DD[j] - r;
 747                   di = SS[j];
 748                   dj = FF[j];
 749                   ORDER(d, di, dj, c, ci, cj)
 750                   c = 0;
 751                   DIAG(i, j, c, p+TC_SCORE(A[i-1], B[j-1]))             /* diagonal */
 752                     //checked
 753
 754                   if ( c <= 0 )
 755                     { c = 0; ci = i; cj = j; }
 756                   else
 757                     { ci = pi; cj = pj; }
 758                   ORDER(c, ci, cj, d, di, dj)
 759                   ORDER(c, ci, cj, f, fi, fj)
 760                   p = CC[j];
 761                   CC[j] = c;
 762                   pi = RR[j];
 763                   pj = EE[j];
 764                   RR[j] = ci;
 765                   EE[j] = cj;
 766                   DD[j] = d;
 767                   SS[j] = di;
 768                   FF[j] = dj;
 769                   if ( c > lmin )       /* add the score into list */
 770                     lmin = addnode(c, ci, cj, i, j, count, lmin);
 771                 }
 772           }
 773 return 1;
 774 }
 775
 776 /* Add a new node into list.  */
 777
 778 int addnode(c, ci, cj, i, j, K, cost)  int c, ci, cj, i, j, K, cost;
 779 { int found;                            /* 1 if the node is in LIST */
 780   register int d;
 781
 782   found = 0;
 783   if ( most != 0 && most->SIM_STARI == ci && most->SIM_STARJ == cj )
 784     found = 1;
 785   else
 786      for ( d = 0; d < numnode ; d++ )
 787         { most = LIST[d];
 788           if ( most->SIM_STARI == ci && most->SIM_STARJ == cj )
 789             { found = 1;
 790               break;
 791             }
 792         }
 793   if ( found )
 794     { if ( most->SIM_SCORE < c )
 795         { most->SIM_SCORE = c;
 796           most->SIM_ENDI = i;
 797           most->SIM_ENDJ = j;
 798         }
 799       if ( most->SIM_TOP > i ) most->SIM_TOP = i;
 800       if ( most->SIM_BOT < i ) most->SIM_BOT = i;
 801       if ( most->SIM_LEFT > j ) most->SIM_LEFT = j;
 802       if ( most->SIM_RIGHT < j ) most->SIM_RIGHT = j;
 803     }
 804   else
 805     { if ( numnode == K )       /* list full */
 806          most = low;
 807       else
 808          most = LIST[numnode++];
 809       most->SIM_SCORE = c;
 810       most->SIM_STARI = ci;
 811       most->SIM_STARJ = cj;
 812       most->SIM_ENDI = i;
 813       most->SIM_ENDJ = j;
 814       most->SIM_TOP = most->SIM_BOT = i;
 815       most->SIM_LEFT = most->SIM_RIGHT = j;
 816     }
 817   if ( numnode == K )
 818     { if ( low == most || ! low )
 819         { for ( low = LIST[0], d = 1; d < numnode ; d++ )
 820             if ( LIST[d]->SIM_SCORE < low->SIM_SCORE )
 821               low = LIST[d];
 822         }
 823       return ( low->SIM_SCORE ) ;
 824     }
 825   else
 826     return cost;
 827 }
 828
 829 /* Find and remove the largest score in list */
 830
 831 vertexptr findmax()
 832 { vertexptr  cur;
 833   register int i, j;
 834
 835   for ( j = 0, i = 1; i < numnode ; i++ )
 836     if ( LIST[i]->SIM_SCORE > LIST[j]->SIM_SCORE )
 837        j = i;
 838   cur = LIST[j];
 839   if ( j != --numnode )
 840     { LIST[j] = LIST[numnode];
 841       LIST[numnode] =  cur;
 842     }
 843   most = LIST[0];
 844   if ( low == cur ) low = LIST[0];
 845   return ( cur );
 846 }
 847
 848 /* return 1 if no node in LIST share vertices with the area */
 849
 850 int no_cross()
 851 { vertexptr  cur;
 852   register int i;
 853
 854       for ( i = 0; i < numnode; i++ )
 855         { cur = LIST[i];
 856           if ( cur->SIM_STARI <= mm && cur->SIM_STARJ <= nn && cur->SIM_BOT >= m1-1 &&
 857                cur->SIM_RIGHT >= n1-1 && ( cur->SIM_STARI < rl || cur->SIM_STARJ < cl ))
 858              { if ( cur->SIM_STARI < rl ) rl = cur->SIM_STARI;
 859                if ( cur->SIM_STARJ < cl ) cl = cur->SIM_STARJ;
 860                flag = 1;
 861                break;
 862              }
 863         }
 864       if ( i == numnode )
 865         return 1;
 866       else
 867         return 0;
 868 }
 869
 870 /* diff(A,B,M,N,tb,te) returns the score of an optimum conversion between
 871    A[1..M] and B[1..N] that begins(ends) with a delete if tb(te) is zero
 872    and appends such a conversion to the current script.                   */
 873
 874 int diff_sim( int *A,int *B,int M,int N,int tb,int te)
 875
 876 { int   midi, midj, type;       /* Midpoint, type, and cost */
 877   int midc;
 878
 879   {
 880     register int   i, j;
 881     register int c, e, d, s;
 882     int t;
 883
 884
 885     /* Boundary cases: M <= 1 or N == 0 */
 886
 887     if (N <= 0)
 888       { if (M > 0) DEL(M)
 889           return - gap(M);
 890       }
 891     if (M <= 1)
 892       { if (M <= 0)
 893           { INS(N);
 894             return - gap(N);
 895           }
 896         if (tb > te) tb = te;
 897         midc = - (tb + r + gap(N) );
 898         midj = 0;
 899
 900         for (j = 1; j <= N; j++)
 901           {  for ( tt = 1, z = row[I+1]; z != PAIRNULL; z = z->NEXT )
 902               if ( z->COL == j+J )
 903                 { tt = 0; break; }
 904             if ( tt )
 905               { c = TC_SCORE (A[0],B[j-1]) - ( gap(j-1) + gap(N-j) );
 906                 //checked
 907
 908                 if (c > midc)
 909                   { midc = c;
 910                     midj = j;
 911                   }
 912               }
 913           }
 914         if (midj == 0)
 915           { INS(N) DEL(1) }
 916         else
 917           { if (midj > 1) INS(midj-1)
 918               REP
 919               if ( A[1] == B[midj] )
 920                 no_mat += 1;
 921               else
 922                 no_mis += 1;
 923             /* mark (A[I],B[J]) as used: put J into list row[I] */
 924             I++; J++;
 925
 926
 927             z = ( pairptr )sim_vcalloc(1,sizeof(pair));
 928             z->COL = J;
 929             z->NEXT = row[I];
 930             row[I] = z;
 931             if (midj < N) INS(N-midj)
 932               }
 933         return midc;
 934       }
 935
 936     /* Divide: Find optimum midpoint (midi,midj) of cost midc */
 937
 938     midi = M/2;                 /* Forward phase:                          */
 939     CC[0] = 0;                  /*   Compute C(M/2,k) & D(M/2,k) for all k */
 940     t = -q;
 941     for (j = 1; j <= N; j++)
 942       { CC[j] = t = t-r;
 943         DD[j] = t-q;
 944       }
 945     t = -tb;
 946     for (i = 1; i <= midi; i++)
 947       { s = CC[0];
 948         CC[0] = c = t = t-r;
 949         e = t-q;
 950
 951         for (j = 1; j <= N; j++)
 952           { if ((c = c - qr) > (e = e - r)) e = c;
 953             if ((c = CC[j] - qr) > (d = DD[j] - r)) d = c;
 954             DIAG(i+I, j+J, c, s+TC_SCORE(A[i-1], B[j-1]))
 955               //checked
 956
 957               if (c < d) c = d;
 958             if (c < e) c = e;
 959             s = CC[j];
 960             CC[j] = c;
 961             DD[j] = d;
 962           }
 963       }
 964     DD[0] = CC[0];
 965
 966     RR[N] = 0;                  /* Reverse phase:                          */
 967     t = -q;                     /*   Compute R(M/2,k) & S(M/2,k) for all k */
 968     for (j = N-1; j >= 0; j--)
 969       { RR[j] = t = t-r;
 970         SS[j] = t-q;
 971       }
 972     t = -te;
 973     for (i = M-1; i >= midi; i--)
 974       { s = RR[N];
 975         RR[N] = c = t = t-r;
 976         e = t-q;
 977
 978         for (j = N-1; j >= 0; j--)
 979           { if ((c = c - qr) > (e = e - r)) e = c;
 980             if ((c = RR[j] - qr) > (d = SS[j] - r)) d = c;
 981             DIAG(i+1+I, j+1+J, c, s+TC_SCORE (A[i],B[j])) /*not -1 on purpose*/
 982
 983               if (c < d) c = d;
 984             if (c < e) c = e;
 985             s = RR[j];
 986             RR[j] = c;
 987             SS[j] = d;
 988           }
 989       }
 990     SS[N] = RR[N];
 991
 992     midc = CC[0]+RR[0];         /* Find optimal midpoint */
 993     midj = 0;
 994     type = 1;
 995     for (j = 0; j <= N; j++)
 996       if ((c = CC[j] + RR[j]) >= midc)
 997         if (c > midc || (CC[j] != DD[j] && RR[j] == SS[j]))
 998           { midc = c;
 999             midj = j;
1000           }
1001     for (j = N; j >= 0; j--)
1002       if ((c = DD[j] + SS[j] + q) > midc)
1003         { midc = c;
1004           midj = j;
1005           type = 2;
1006         }
1007   }
1008
1009   /* Conquer: recursively around midpoint */
1010
1011   if (type == 1)
1012     { diff_sim(A,B,midi,midj,tb,q);
1013       diff_sim(A+midi,B+midj,M-midi,N-midj,q,te);
1014     }
1015   else
1016     { diff_sim(A,B,midi-1,midj,tb,zero);
1017       DEL(2);
1018       diff_sim(A+midi+1,B+midj,M-midi-1,N-midj,zero,te);
1019     }
1020   return midc;
1021 }
1022
1023
1024
1025
1026 int calcons(int *aa0,int n0,int *aa1,int n1,int *res,int *nc,int *nident, Alignment *A, int *ns, int **l_s, Constraint_list *CL)
1027 {
1028   int i0, i1;
1029   int op, nid, lenc, nd;
1030   int *sp0, *sp1;
1031   int *rp;
1032   int a, b, id_col, tot_col, r0, r1;
1033
1034   min0--; min1--;
1035
1036   sp0 = seqc0+mins;
1037   sp1 = seqc1+mins;
1038   rp = res;
1039   lenc = nid = op = 0;
1040   i0 = min0;
1041   i1 = min1;
1042
1043   while (i0 < max0 || i1 < max1) {
1044     if (op == 0 && *rp == 0) {
1045       op = *rp++;
1046       *sp0 = aa0[i0++];
1047       *sp1 = aa1[i1++];
1048
1049
1050       for (id_col=tot_col=0,a=0; a< ns[0]; a++)
1051         for ( b=0; b< ns[1]; b++)
1052           {
1053             r0=Aln->seq_al[l_s[0][a]][*sp0-1];
1054             r1=Aln->seq_al[l_s[1][a]][*sp1-1];
1055
1056             if ( !is_gap(r0) && r1==r0)id_col++;
1057             if ( !is_gap(r0) && !is_gap(r1))tot_col++;
1058           }
1059       nid+=(tot_col)?(id_col/tot_col):0;
1060       lenc++;
1061       sp0++; sp1++;
1062     }
1063     else {
1064       if (op==0) op = *rp++;
1065       if (op>0) {
1066         *sp0++ = SIM_GAP;
1067         *sp1++ = aa1[i1++];
1068         op--;
1069         lenc++;
1070       }
1071       else {
1072         *sp0++ = aa0[i0++];
1073         *sp1++ = SIM_GAP;
1074         op++;
1075         lenc++;
1076       }
1077     }
1078   }
1079
1080   *nident = nid;
1081   *nc = lenc;
1082
1083   nd = 0;
1084   return mins+lenc+nd;
1085 }
1086
1087 /*Memory management */
1088 struct Mem
1089     {
1090     void   *p;
1091     struct Mem *next;
1092     };
1093
1094 typedef struct Mem Mem;
1095
1096 Mem *first_mem;
1097 Mem *last_mem;
1098
1099 void *sim_vcalloc ( size_t nobj, size_t size)
1100 {
1101   void *p;
1102   Mem *new_mem;
1103
1104   p=vcalloc (nobj, size);
1105
1106
1107   new_mem=vcalloc (1, sizeof (Mem));
1108   if ( last_mem==NULL)first_mem=last_mem=new_mem;
1109   else
1110     {
1111       last_mem->next=new_mem;
1112       last_mem=new_mem;
1113     }
1114   last_mem->p=p;
1115   return p;
1116 }
1117
1118 void sim_free_all()
1119 {
1120   Mem *p1, *p2;
1121   p1=first_mem;
1122
1123
1124   while (p1)
1125     {
1126       p2=p1->next;
1127       vfree(p1->p);
1128       vfree(p1);
1129       p1=p2;
1130     }
1131   first_mem=last_mem=NULL;
1132   sim_reset_static_variable();
1133 }
1134
1135 /*********************************COPYRIGHT NOTICE**********************************/
1136 /*© Centro de Regulacio Genomica */
1137 /*and */
1138 /*Cedric Notredame */
1139 /*Tue Oct 27 10:12:26 WEST 2009. */
1140 /*All rights reserved.*/
1141 /*This file is part of T-COFFEE.*/
1142 /**/
1143 /*    T-COFFEE is free software; you can redistribute it and/or modify*/
1144 /*    it under the terms of the GNU General Public License as published by*/
1145 /*    the Free Software Foundation; either version 2 of the License, or*/
1146 /*    (at your option) any later version.*/
1147 /**/
1148 /*    T-COFFEE is distributed in the hope that it will be useful,*/
1149 /*    but WITHOUT ANY WARRANTY; without even the implied warranty of*/
1150 /*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the*/
1151 /*    GNU General Public License for more details.*/
1152 /**/
1153 /*    You should have received a copy of the GNU General Public License*/
1154 /*    along with Foobar; if not, write to the Free Software*/
1155 /*    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA*/
1156 /*...............................................                                                                                      |*/
1157 /*  If you need some more information*/
1158 /*  cedric.notredame@europe.com*/
1159 /*...............................................                                                                                                                                     |*/
1160 /**/
1161 /**/
1162 /*      */
1163 /*********************************COPYRIGHT NOTICE**********************************/