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

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <math.h>
   4 #include <stdarg.h>
   5 #include <string.h>
   6 #include <ctype.h>
   7 #include "io_lib_header.h"
   8 #include "util_lib_header.h"
   9 #include "define_header.h"
  10
  11 #include "dp_lib_header.h"
  12 float compute_lambda (int **matrix,char *alphabet);
  13 /*********************************************************************************************/
  14 /*                                                                                           */
  15 /*         FUNCTIONS FOR EVALUATING THE CONSISTENCY BETWEEN ALN AND CL                       */
  16 /*                                                                                           */
  17 /*********************************************************************************************/
  18
  19 /*Fast:         score= extended_res/max_extended_residue_for the whole aln
  20   slow:         score= extended_res/sum all extended score for that residue
  21   non_extended  score= non_ext     /sum all non extended  score for that residue
  22   heuristic     score= extended    /sum of extended score of all pairs in the library
  23                                     (i.e. Not ALL the possible pairs)
  24 */
  25 Alignment * main_coffee_evaluate_output2 ( Alignment *IN,Constraint_list *CL, const char *mode );
  26 int sub_aln2ecl_raw_score (Alignment *A, Constraint_list *CL, int ns, int *ls)
  27 {
  28   int **pos;
  29   int p1,r1, r2, s1, s2;
  30   int score=0;
  31
  32   if ( !A) return 0;
  33   A=reorder_aln (A, (CL->S)->name,(CL->S)->nseq);
  34   pos=aln2pos_simple ( A,A->nseq);
  35
  36
  37
  38   for (p1=0; p1<A->len_aln; p1++)
  39     {
  40       for (s1=0; s1<ns-1; s1++)
  41         {
  42           for (s2=s1+1; s2<ns; s2++)
  43             {
  44
  45               r1=pos[ls[s1]][p1];
  46               r2=pos[ls[s2]][p1];
  47               if (r1>0 && r2>0)
  48                 {
  49                   score+= residue_pair_extended_list_pc (CL,ls[s1], r1, ls[s2], r2)*SCORE_K;
  50                 }
  51             }
  52         }
  53     }
  54   free_int (pos, -1);
  55   return score;
  56   return (score/(((ns*ns)-ns)/2))/A->len_aln;
  57 }
  58 int aln2ecl_raw_score (Alignment *A, Constraint_list *CL)
  59 {
  60   int **pos;
  61   int p1,r1, r2, s1, s2;
  62   int score=0;
  63   if ( !A) return 0;
  64   A=reorder_aln (A, (CL->S)->name,(CL->S)->nseq);
  65   pos=aln2pos_simple ( A,A->nseq);
  66
  67
  68
  69   for (p1=0; p1<A->len_aln; p1++)
  70     {
  71       for (s1=0; s1<A->nseq-1; s1++)
  72         {
  73           for (s2=s1+1; s2<A->nseq; s2++)
  74             {
  75
  76               r1=pos[s1][p1];
  77               r2=pos[s2][p1];
  78               if (r1>0 && r2>0)score+= residue_pair_extended_list_pc (CL,s1, r1, s2, r2);
  79             }
  80         }
  81     }
  82   free_int (pos, -1);
  83   return score;
  84   return (score/(((A->nseq*A->nseq)-A->nseq)/2))/A->len_aln;
  85 }
  86 int node2sub_aln_score    (Alignment *A,Constraint_list *CL, char *mode, NT_node T)
  87 {
  88   if ( !T || !T->right ||!T->left)return 0;
  89   else
  90     {
  91       int *ns;
  92       int **ls;
  93       ns=vcalloc (2, sizeof (int));
  94       ls=vcalloc (2, sizeof (int*));
  95       ns[0]= (T->left)->nseq;
  96       ns[1]=(T->right)->nseq;
  97       ls[0]= (T->left)->lseq;
  98       ls[1]=(T->right)->lseq;
  99
 100       return sub_aln2sub_aln_score (A, CL, mode, ns, ls);
 101     }
 102   return -1;
 103 }
 104 int sub_aln2sub_aln_score ( Alignment *A,Constraint_list *CL, const char *mode, int *ns, int **ls)
 105 {
 106   /*Warning: Does Not Take Gaps into account*/
 107
 108   int **pos;
 109   int a;
 110   float score=0;
 111
 112   /*Make sure evaluation functions update their cache if needed*/
 113   A=update_aln_random_tag (A);
 114   pos=aln2pos_simple ( A, -1, ns, ls);
 115   for (a=0; a< A->len_aln; a++)
 116     score+=CL->get_dp_cost (A, pos, ns[0], ls[0], a, pos, ns[1],ls[1], a, CL);
 117   free_int (pos, -1);
 118
 119   score=(int)(((float)score)/(A->len_aln*SCORE_K));
 120   score=(int)(CL->residue_index && CL->normalise)?((score*MAXID)/(CL->normalise)):(score);
 121   return (int)score;
 122 }
 123 int sub_aln2sub_aln_raw_score ( Alignment *A,Constraint_list *CL, const char *mode, int *ns, int **ls)
 124 {
 125   /*Warning: Does Not Take Gaps into account*/
 126
 127   int **pos;
 128   int a;
 129   float score=0;
 130
 131   /*Make sure evaluation functions update their cache if needed*/
 132   A=update_aln_random_tag (A);
 133   pos=aln2pos_simple ( A, -1, ns, ls);
 134   for (a=0; a< A->len_aln; a++)
 135     score+=CL->get_dp_cost (A, pos, ns[0], ls[0], a, pos, ns[1],ls[1], a, CL);
 136   free_int (pos, -1);
 137   return (int) score;
 138 }
 139
 140 Alignment* main_coffee_evaluate_output_sub_aln ( Alignment *A,Constraint_list *CL, const char *mode, int *n_s, int **l_s)
 141 {
 142   Alignment *SUB1, *SUB2, *SUB3;
 143   int a, b, c,*list_seq;
 144
 145
 146   if (strm ( CL->evaluate_mode, "no"))return NULL;
 147   else
 148     {
 149         list_seq=vcalloc (n_s[0]+n_s[1], sizeof (int));
 150         for (b=0, a=0; a< 2; a++){for (c=0;c< n_s[a]; c++)list_seq[b++]=l_s[a][c];}
 151
 152
 153         SUB1=copy_aln (A, NULL);
 154         SUB2=extract_sub_aln (SUB1,n_s[0]+n_s[1],list_seq);
 155         SUB3=main_coffee_evaluate_output (SUB2,CL,CL->evaluate_mode);
 156         free_aln (SUB1);
 157         free_aln (SUB2);
 158         vfree (list_seq);
 159
 160         return SUB3;
 161     }
 162 }
 163 Alignment * overlay_alignment_evaluation     ( Alignment *I, Alignment *O)
 164 {
 165   int a, b, c, r, i;
 166   int *buf;
 167
 168   if ( !I || !O) return O;
 169   if ( I->len_aln!=O->len_aln)printf_exit (EXIT_FAILURE, stderr, "ERROR: Incompatible alignments in overlay_alignment_evaluation");
 170
 171   buf=vcalloc ( MAX(I->len_aln, O->len_aln), sizeof (int));
 172
 173   for (a=0; a<O->nseq; a++)
 174     {
 175       if (!strm (I->name[a], O->name[a]))printf_exit (EXIT_FAILURE, stderr, "ERROR: Incompatible alignments in overlay_alignment_evaluation");
 176       for (b=0; b<O->len_aln; b++)
 177         {
 178           r=I->seq_al[a][b];
 179           if ( islower(r))O->seq_al[a][b]=0;
 180           else if (r<=9 || (r>='0' && r<='9'))O->seq_al[a][b]=I->seq_al[a][b];
 181         }
 182     }
 183   return O;
 184 }
 185
 186 Alignment * main_coffee_evaluate_output ( Alignment *IN,Constraint_list *CL, const char *mode )
 187 {
 188
 189   Alignment *TopS=NULL, *LastS=NULL, *CurrentS=NULL;
 190
 191
 192   if ( IN->A)IN=IN->A;
 193   while (IN)
 194     {
 195
 196       CurrentS= main_coffee_evaluate_output2(IN, CL, mode);
 197       if (!TopS)LastS=TopS=CurrentS;
 198       else
 199         {
 200           LastS->A=CurrentS;
 201           LastS=CurrentS;
 202         }
 203       IN=IN->A;
 204     }
 205   return TopS;
 206 }
 207
 208 Alignment * main_coffee_evaluate_output2 ( Alignment *IN,Constraint_list *CL, const char *mode )
 209 {
 210
 211   /*Make sure evaluation functions update their cache if needed*/
 212   IN=update_aln_random_tag (IN);
 213
 214   if ( CL->evaluate_residue_pair==evaluate_matrix_score || CL->ne==0 ||strm ( mode , "categories") || strm ( mode , "matrix")|| strm(mode, "sar")|| strstr (mode, "boxshade") )
 215      {
 216
 217        if ( strm ( mode , "categories")) return categories_evaluate_output (IN, CL);
 218        else if ( strm ( mode , "matrix"))return matrix_evaluate_output (IN, CL);
 219        else if ( strm ( mode, "sar"))return sar_evaluate_output (IN, CL);
 220        else if ( strstr ( mode, "boxshade"))return boxshade_evaluate_output (IN, CL, atoi (strstr(mode, "_")+1));
 221
 222        else if ( CL->evaluate_residue_pair==evaluate_matrix_score) return matrix_evaluate_output (IN, CL);
 223        else if ( CL->ne==0) return matrix_evaluate_output (IN, CL);
 224      }
 225    else if ( strm (mode, "no"))return NULL;
 226    else if ( strstr( mode, "triplet"))
 227      {
 228        return triplet_coffee_evaluate_output ( IN,CL);
 229      }
 230    else if ( strstr ( mode, "fast"))
 231      {
 232        return fast_coffee_evaluate_output ( IN,CL);
 233      }
 234    else if ( strstr ( mode, "slow"))
 235      {
 236        return slow_coffee_evaluate_output ( IN,CL);
 237      }
 238
 239    else if ( strstr ( mode, "non_extended"))
 240      {
 241        return non_extended_t_coffee_evaluate_output ( IN,CL);
 242      }
 243
 244    else if ( strm (mode, "sequences"))
 245      {
 246        return coffee_seq_evaluate_output ( IN,CL);
 247      }
 248    else
 249      {
 250        fprintf ( stderr, "\nUNKNOWN MODE FOR ALIGNMENT EVALUATION: *%s* [FATAL:%s]",mode, PROGRAM);
 251        crash ("");
 252        return NULL;
 253      }
 254   return IN;
 255 }
 256
 257
 258
 259 Alignment * coffee_evaluate_output ( Alignment *IN,Constraint_list *CL)
 260     {
 261     fprintf ( stderr, "\n[WARNING:%s]THE FUNCTION coffee_evaluate_output IS NOT ANYMORE SUPPORTED\n", PROGRAM);
 262     fprintf ( stderr, "\n[WARNING]fast_coffee_evaluate_output WILL BE USED INSTEAD\n");
 263
 264     return fast_coffee_evaluate_output (IN,CL);
 265     }
 266 Alignment * matrix_evaluate_output ( Alignment *IN,Constraint_list *CL)
 267     {
 268     int a,b, c,r, s, r1, r2;
 269     Alignment *OUT=NULL;
 270
 271     double **tot_res;
 272     double **max_res;
 273
 274     double **tot_seq;
 275     double **max_seq;
 276
 277     double **tot_col;
 278     double **max_col;
 279
 280     double max_aln=0;
 281     double tot_aln=0;
 282
 283
 284     /*
 285       Residue x: sum of observed extended X.. /sum of possible X..
 286     */
 287
 288
 289     if ( !CL->M)CL->M=read_matrice ("blosum62mt");
 290
 291     OUT=copy_aln (IN, OUT);
 292
 293
 294     tot_res=declare_double ( IN->nseq, IN->len_aln);
 295     max_res=declare_double ( IN->nseq, IN->len_aln);
 296
 297     tot_seq=declare_double ( IN->nseq, 1);
 298     max_seq=declare_double ( IN->nseq, 1);
 299     tot_col=declare_double ( IN->len_aln,1);
 300     max_col=declare_double ( IN->len_aln,1);
 301
 302     max_aln=tot_aln=0;
 303
 304     for (a=0; a< IN->len_aln; a++)
 305         {
 306             for ( b=0; b< IN->nseq; b++)
 307                 {
 308                   r1=tolower(IN->seq_al[b][a]);
 309                   if ( is_gap(r1))continue;
 310                   r= CL->M[r1-'A'][r1-'A'];
 311                   r= 1;
 312                   for ( c=0; c<IN->nseq; c++)
 313                     {
 314                       r2=tolower(IN->seq_al[c][a]);
 315                       if (b==c || is_gap (r2))continue;
 316
 317                       s=CL->M[r2-'A'][r1-'A'];
 318                       s=(s<=0)?0:1;
 319
 320                       tot_res[b][a]+=s;
 321                       max_res[b][a]+=r;
 322
 323                       tot_col[a][0]+=s;
 324                       max_col[a][0]+=r;
 325
 326                       tot_seq[b][0]+=s;
 327                       max_seq[b][0]+=r;
 328
 329                       tot_aln+=s;
 330                       max_aln+=r;
 331                     }
 332                 }
 333         }
 334
 335
 336     for ( a=0; a< IN->nseq; a++)
 337       {
 338         if ( !max_seq[a][0])continue;
 339
 340         OUT->score_seq[a]=(tot_seq[a][0]*100)/max_seq[a][0];
 341         for (b=0; b< IN->len_aln; b++)
 342           {
 343             r1=IN->seq_al[a][b];
 344             if ( is_gap(r1) || !max_res[a][b])continue;
 345
 346             r1=(tot_res[a][b]*10)/max_res[a][b];
 347             r1=(r1>=10)?9:r1;
 348             r1=r1<0?0:r1;
 349             (OUT)->seq_al[a][b]=r1+'0';
 350           }
 351       }
 352
 353     for ( a=0; a< IN->len_aln; a++)
 354       {
 355         r1=(max_col[a][0]==0)?0:((tot_col[a][0]*10)/max_col[a][0]);
 356         r1=(r1>=10)?9:r1;
 357         (OUT)->seq_al[OUT->nseq][a]=r1+'0';
 358       }
 359     sprintf ( OUT->name[IN->nseq], "cons");
 360     if (max_aln)OUT->score_seq[OUT->nseq]=OUT->score_aln=(100*tot_aln)/max_aln;
 361
 362
 363     free_double (tot_res,-1);
 364     free_double (max_res,-1);
 365
 366     free_double (tot_seq,-1);
 367     free_double (max_seq,-1);
 368
 369     return OUT;
 370     }
 371
 372 Alignment * sar_evaluate_output ( Alignment *IN,Constraint_list *CL)
 373     {
 374       int a,b, c,r, s, r1, r2;
 375       Alignment *OUT=NULL;
 376
 377       double **tot_res;
 378       double **max_res;
 379
 380       double **tot_seq;
 381       double **max_seq;
 382
 383       double **tot_col;
 384       double **max_col;
 385
 386       double max_aln=0;
 387       double tot_aln=0;
 388
 389
 390     /*
 391       Residue x: sum of observed extended X.. /sum of possible X..
 392     */
 393
 394
 395     if ( !CL->M)CL->M=read_matrice ("blosum62mt");
 396
 397     OUT=copy_aln (IN, OUT);
 398
 399
 400     tot_res=declare_double ( IN->nseq, IN->len_aln);
 401     max_res=declare_double ( IN->nseq, IN->len_aln);
 402
 403     tot_seq=declare_double ( IN->nseq, 1);
 404     max_seq=declare_double ( IN->nseq, 1);
 405     tot_col=declare_double ( IN->len_aln,1);
 406     max_col=declare_double ( IN->len_aln,1);
 407
 408     max_aln=tot_aln=0;
 409
 410     for (a=0; a< IN->len_aln; a++)
 411         {
 412           for (b=0; b< IN->nseq; b++)
 413                 {
 414                   r1=tolower(IN->seq_al[b][a]);
 415                   for (c=0; c<IN->nseq; c++)
 416                     {
 417                       r2=tolower(IN->seq_al[c][a]);
 418                       if (b==c)continue;
 419
 420                       if ( is_gap(r1) && is_gap(r2))s=0;
 421                       else s=(r1==r2)?1:0;
 422
 423                       r=1;
 424
 425
 426                       tot_res[b][a]+=s;
 427                       max_res[b][a]+=r;
 428
 429                       tot_col[a][0]+=s;
 430                       max_col[a][0]+=r;
 431
 432                       tot_seq[b][0]+=s;
 433                       max_seq[b][0]+=r;
 434
 435                       tot_aln+=s;
 436                       max_aln+=r;
 437                     }
 438                 }
 439         }
 440
 441     for ( a=0; a< IN->nseq; a++)
 442       {
 443         if ( !max_seq[a][0])continue;
 444         OUT->score_seq[a]=(max_seq[a][0]*100)/max_seq[a][0];
 445         for (b=0; b< IN->len_aln; b++)
 446           {
 447             r1=IN->seq_al[a][b];
 448             if ( is_gap(r1) || !max_res[a][b])continue;
 449             r1=(tot_res[a][b]*10)/max_res[a][b];
 450             r1=(r1>=10)?9:r1;
 451             r1=r1<0?0:r1;
 452             (OUT)->seq_al[a][b]=r1+'0';
 453           }
 454       }
 455
 456     for ( a=0; a< IN->len_aln; a++)
 457       {
 458         r1=(max_col[a][0]==0)?0:((tot_col[a][0]*10)/max_col[a][0]);
 459         r1=(r1>=10)?9:r1;
 460         (OUT)->seq_al[OUT->nseq][a]=r1+'0';
 461       }
 462     sprintf ( OUT->name[IN->nseq], "cons");
 463     if (max_aln)OUT->score_aln=(100*tot_aln)/max_aln;
 464
 465
 466     free_double (tot_res,-1);
 467     free_double (max_res,-1);
 468
 469     free_double (tot_seq,-1);
 470     free_double (max_seq,-1);
 471
 472     return OUT;
 473     }
 474 Alignment * boxshade_evaluate_output ( Alignment *IN,Constraint_list *CL, int T)
 475     {
 476       Alignment *OUT=NULL;
 477       int **aa;
 478       int r,br, bs, a, b;
 479       float f;
 480
 481
 482     /*
 483       Residue x: sum of observed extended X.. /sum of possible X..
 484     */
 485
 486
 487     OUT=copy_aln (IN, OUT);
 488     aa=declare_int (26, 2);
 489
 490     for ( a=0; a< OUT->len_aln; a++)
 491       {
 492         for ( b=0; b< 26; b++){aa[b][1]=0;aa[b][0]='a'+b;}
 493         for ( b=0; b< OUT->nseq; b++)
 494           {
 495             r=tolower(OUT->seq_al[b][a]);
 496             if ( !is_gap(r))aa[r-'a'][1]++;
 497           }
 498         sort_int ( aa, 2, 1, 0,25);
 499         f=(aa[25][1]*100)/OUT->nseq;
 500
 501         if (f<T);
 502         else
 503           {
 504             bs=((f-T)*10)/(100-T);
 505             br=aa[25][0];
 506
 507             if (bs==10)bs--;bs+='0';
 508             for ( b=0; b< OUT->nseq; b++)
 509               {
 510                 r=tolower(OUT->seq_al[b][a]);
 511                 if (r==br && bs>'1')OUT->seq_al[b][a]=bs;
 512               }
 513             OUT->seq_al[b][a]=bs;
 514           }
 515       }
 516     sprintf ( OUT->name[IN->nseq], "cons");
 517
 518     return OUT;
 519     }
 520
 521 Alignment * categories_evaluate_output ( Alignment *IN,Constraint_list *CL)
 522     {
 523
 524     Alignment *OUT=NULL;
 525     int a, b, r;
 526     int *aa;
 527     float score, nseq2, tot_aln;
 528     float n;
 529     /*
 530       Residue x: sum of observed extended X.. /sum of possible X..
 531     */
 532     OUT=copy_aln (IN, OUT);
 533     aa=vcalloc ( 26, sizeof (int));
 534     nseq2=IN->nseq*IN->nseq;
 535
 536     for (tot_aln=0, a=0; a< IN->len_aln; a++)
 537       {
 538         for (n=0,b=0; b< IN->nseq; b++)
 539           {
 540             r=IN->seq_al[b][a];
 541
 542             if ( is_gap(r))n++;
 543             else
 544               {
 545                 aa[tolower(r)-'a']++;
 546                 n++;
 547               }
 548           }
 549         n=n*n;
 550         for ( score=0,b=0; b< 26; b++){score+=aa[b]*aa[b];aa[b]=0;}
 551         /*score/=nseq2;*/
 552         score=(n==0)?0:score/n;
 553         tot_aln+=score;
 554         r=score*10;
 555         r=(r>=10)?9:r;
 556         (OUT)->seq_al[OUT->nseq][a]='0'+r;
 557       }
 558     OUT->score_aln=(tot_aln/OUT->len_aln)*100;
 559     sprintf ( OUT->name[IN->nseq], "cons");
 560     vfree(aa);
 561     return OUT;
 562     }
 563
 564 Alignment * categories_evaluate_output_old ( Alignment *IN,Constraint_list *CL)
 565     {
 566
 567     Alignment *OUT=NULL;
 568     int nc,a, b, r;
 569     int *aa, ng;
 570     float score, nseq2, tot_aln, min=0;
 571
 572     /*
 573       Residue x: sum of observed extended X.. /sum of possible X..
 574     */
 575     OUT=copy_aln (IN, OUT);
 576     aa=vcalloc ( 26, sizeof (int));
 577     nseq2=IN->nseq*IN->nseq;
 578
 579     for (tot_aln=0, a=0; a< IN->len_aln; a++)
 580       {
 581         for (ng=0,b=0; b< IN->nseq; b++)
 582           {
 583             r=IN->seq_al[b][a];
 584
 585             if ( is_gap(r))ng++;
 586             else
 587               {
 588                 aa[tolower(r)-'a']++;
 589               }
 590           }
 591         for (nc=0, b=0; b<26; b++)
 592           {
 593             if ( aa[b])nc++;
 594             aa[b]=0;
 595           }
 596         if (nc>9)score=0;
 597         else score=9-nc;
 598
 599         score=(2*min)/IN->nseq;
 600
 601         tot_aln+=score;
 602         r=score*10;
 603         r=(r>=10)?9:r;
 604         (OUT)->seq_al[OUT->nseq][a]='0'+r;
 605       }
 606
 607     OUT->score_aln=(tot_aln/OUT->len_aln)*100;
 608     sprintf ( OUT->name[IN->nseq], "cons");
 609     vfree(aa);
 610     return OUT;
 611     }
 612 Alignment * fork_triplet_coffee_evaluate_output ( Alignment *IN,Constraint_list *CL);
 613 Alignment * nfork_triplet_coffee_evaluate_output ( Alignment *IN,Constraint_list *CL);
 614 Alignment * triplet_coffee_evaluate_output ( Alignment *IN,Constraint_list *CL)
 615 {
 616
 617   if (!IN || !CL || !CL->residue_index) return IN;
 618
 619
 620   if ( get_nproc()==1)return  nfork_triplet_coffee_evaluate_output (IN,CL);
 621   else if (strstr ( CL->multi_thread, "evaluate"))return  fork_triplet_coffee_evaluate_output (IN,CL);
 622   else return nfork_triplet_coffee_evaluate_output (IN,CL);
 623 }
 624 Alignment * fork_triplet_coffee_evaluate_output ( Alignment *IN,Constraint_list *CL)
 625     {
 626       Alignment *OUT=NULL;
 627       int **pos;
 628
 629       double score_col=0, score_aln=0, score_res=0, score=0;
 630       double max_col=0, max_aln=0, max_res=0;
 631       double *max_seq, *score_seq;
 632
 633       int a,b, x, y,res;
 634       int s1,r1,s2,r2,w2,s3,r3,w3;
 635       int **lu;
 636
 637       //multi-threading
 638       FILE *fp;
 639       char **pid_tmpfile;
 640       int sjobs, njobs;
 641       int **sl;
 642       int j;
 643
 644       OUT=copy_aln (IN, OUT);
 645       pos=aln2pos_simple(IN, IN->nseq);
 646       sprintf ( OUT->name[IN->nseq], "cons");
 647
 648       max_seq=vcalloc ( IN->nseq, sizeof (double));
 649       score_seq=vcalloc ( IN->nseq, sizeof (double));
 650       lu=declare_int (IN->nseq, IN->len_aln+1);
 651
 652       //multi Threading stuff
 653       njobs=get_nproc();
 654       sl=n2splits (njobs,IN->len_aln);
 655       pid_tmpfile=vcalloc (njobs, sizeof (char*));
 656
 657       for (sjobs=0,j=0; sjobs<njobs; j++)
 658         {
 659           pid_tmpfile[j]=vtmpnam(NULL);
 660           if (vvfork (NULL)==0)
 661             {
 662               initiate_vtmpnam(NULL);
 663               fp=vfopen (pid_tmpfile[j], "w");
 664               score_aln=max_aln=0;
 665               for (a=sl[j][0]; a<sl[j][1]; a++)
 666                 {
 667                   if (j==0)output_completion (stderr,a,sl[0][1],1, "Final Evaluation");
 668                   score_col=max_col=0;
 669                   for (b=0; b<IN->nseq; b++)
 670                     {
 671                       s1=IN->order[b][0];
 672                       r1=pos[b][a];
 673                       if (r1>=0)lu[s1][r1]=1;
 674                     }
 675                   for (b=0; b<IN->nseq; b++)
 676                     {
 677                       score_res=max_res=NORM_F;
 678                       res=NO_COLOR_RESIDUE;
 679                       s1=IN->order[b][0];
 680                       r1=pos[b][a];
 681                       if (r1<=0)
 682                         {
 683                           fprintf (fp, "%d ", res);
 684                           continue;
 685                         }
 686
 687                       for (x=1;x<CL->residue_index[s1][r1][0];x+=ICHUNK)
 688                         {
 689                           s2=CL->residue_index[s1][r1][x+SEQ2];
 690                           r2=CL->residue_index[s1][r1][x+R2];
 691                           w2=CL->residue_index[s1][r1][x+WE];
 692                           for (y=1; y<CL->residue_index[s2][r2][0];y+=ICHUNK)
 693                             {
 694                               s3=CL->residue_index[s2][r2][y+SEQ2];
 695                               r3=CL->residue_index[s2][r2][y+R2];
 696                               w3=CL->residue_index[s2][r2][y+WE];
 697
 698                               score=MIN(w2,w3);
 699
 700                               max_res+=score;
 701                               max_col+=score;
 702                               max_seq[b]+=score;
 703                               max_aln+=score;
 704
 705                               if (lu[s3][r3])
 706                                 {
 707                                   score_res+=score;
 708                                   score_col+=score;
 709                                   score_seq[b]+=score;
 710                                   score_aln+=score;
 711                                 }
 712                             }
 713                         }
 714                       res=(max_res==0)?NO_COLOR_RESIDUE:((score_res*10)/max_res);
 715                       res=(res==NO_COLOR_RESIDUE)?res:(MIN(res,9));
 716                       fprintf ( fp, "%d ", res);
 717                     }
 718                   for (b=0; b<IN->nseq; b++)
 719                     {
 720                       s1=IN->order[b][0];
 721                       r1=pos[b][a];
 722                       if (r1>0)lu[s1][r1]=0;
 723                     }
 724
 725                   res=(max_col==0)?NO_COLOR_RESIDUE:((score_col*10)/max_col);
 726                   res=(res==NO_COLOR_RESIDUE)?res:(MIN(res,9));
 727                   fprintf (fp, "%d ", res);
 728                   for (b=0; b<IN->nseq; b++)fprintf (fp, "%f %f ", score_seq[b], max_seq[b]);
 729                 }
 730               fprintf (fp, "%f %f ", score_aln, max_aln);
 731               vfclose (fp);
 732               myexit (EXIT_SUCCESS);
 733             }
 734           else
 735             {
 736               sjobs++;
 737             }
 738         }
 739
 740       while (sjobs>=0){vwait(NULL); sjobs--;}//wait for all jobs to complete
 741
 742       for (j=0; j<njobs; j++)
 743         {
 744           float sseq, mseq;
 745           fp=vfopen (pid_tmpfile[j], "r");
 746           for (a=sl[j][0];a<sl[j][1]; a++)
 747             {
 748               for (b=0; b<=IN->nseq; b++)//don't forget the consensus
 749                 {
 750                   fscanf (fp, "%d ", &res);
 751                   OUT->seq_al[b][a]=res;
 752                 }
 753               for (b=0; b<IN->nseq; b++)
 754                 {
 755                   fscanf (fp, "%f %f", &sseq, &mseq);
 756                   score_seq[b]+=sseq;
 757                   max_seq[b]+=mseq;
 758                 }
 759             }
 760           fscanf (fp, "%f %f", &sseq, &mseq);
 761           score_aln+=sseq;
 762           max_aln+=mseq;
 763           vfclose (fp);
 764           remove (pid_tmpfile[j]);
 765         }
 766       fprintf ( stderr, "\n");
 767
 768       IN->score_aln=OUT->score_aln=(max_aln==0)?0:((score_aln*100)/max_aln);
 769       for ( a=0; a< OUT->nseq; a++)
 770         {
 771           OUT->score_seq[a]=(max_seq[a]==0)?0:((score_seq[a]*100)/max_seq[a]);
 772         }
 773
 774       free_int (lu,-1);
 775       free_int (pos , -1);
 776       vfree (pid_tmpfile);
 777       free_int (sl, -1);
 778       vfree ( score_seq);
 779       vfree ( max_seq);
 780       return OUT;
 781     }
 782
 783 Alignment * nfork_triplet_coffee_evaluate_output ( Alignment *IN,Constraint_list *CL)
 784     {
 785       Alignment *OUT=NULL;
 786       int **pos;
 787
 788       double score_col=0, score_aln=0, score_res=0, score=0;
 789       double max_col=0, max_aln=0, max_res=0;
 790       double *max_seq, *score_seq;
 791
 792       int a,b, x, y,res;
 793       int s1,r1,s2,r2,w2,s3,r3,w3;
 794       int **lu;
 795
 796
 797
 798       OUT=copy_aln (IN, OUT);
 799       pos=aln2pos_simple(IN, IN->nseq);
 800       sprintf ( OUT->name[IN->nseq], "cons");
 801
 802       max_seq=vcalloc ( IN->nseq, sizeof (double));
 803       score_seq=vcalloc ( IN->nseq, sizeof (double));
 804       lu=declare_int (IN->nseq, IN->len_aln+1);
 805
 806
 807
 808
 809       score_aln=max_aln=0;
 810       for (a=0; a<IN->len_aln; a++)
 811         {
 812           output_completion (stderr,a,IN->len_aln,1, "Final Evaluation");
 813           score_col=max_col=0;
 814           for (b=0; b<IN->nseq; b++)
 815             {
 816               s1=IN->order[b][0];
 817               r1=pos[b][a];
 818               if (r1>=0)lu[s1][r1]=1;
 819             }
 820           for (b=0; b<IN->nseq; b++)
 821             {
 822               score_res=max_res=NORM_F;
 823               OUT->seq_al[b][a]=NO_COLOR_RESIDUE;
 824               s1=IN->order[b][0];
 825               r1=pos[b][a];
 826               if (r1<=0)continue;
 827               for (x=1;x<CL->residue_index[s1][r1][0];x+=ICHUNK)
 828                 {
 829                   s2=CL->residue_index[s1][r1][x+SEQ2];
 830                   r2=CL->residue_index[s1][r1][x+R2];
 831                   w2=CL->residue_index[s1][r1][x+WE];
 832                   for (y=1; y<CL->residue_index[s2][r2][0];y+=ICHUNK)
 833                     {
 834                       s3=CL->residue_index[s2][r2][y+SEQ2];
 835                       r3=CL->residue_index[s2][r2][y+R2];
 836                       w3=CL->residue_index[s2][r2][y+WE];
 837
 838                       score=MIN(w2,w3);
 839
 840                       max_res+=score;
 841                       max_col+=score;
 842                       max_seq[b]+=score;
 843                       max_aln+=score;
 844
 845                       if (lu[s3][r3])
 846                         {
 847                           score_res+=score;
 848                           score_col+=score;
 849                           score_seq[b]+=score;
 850                           score_aln+=score;
 851                         }
 852                     }
 853                 }
 854               res=(max_res==0)?NO_COLOR_RESIDUE:((score_res*10)/max_res);
 855               res=(res==NO_COLOR_RESIDUE)?res:(MIN(res,9));
 856               OUT->seq_al[b][a]=res;
 857             }
 858           for (b=0; b<IN->nseq; b++)
 859             {
 860               s1=IN->order[b][0];
 861               r1=pos[b][a];
 862               if (r1>0)lu[s1][r1]=0;
 863             }
 864
 865           res=(max_col==0)?NO_COLOR_RESIDUE:((score_col*10)/max_col);
 866           res=(res==NO_COLOR_RESIDUE)?res:(MIN(res,9));
 867           OUT->seq_al[IN->nseq][a]=res;
 868         }
 869       fprintf ( stderr, "\n");
 870
 871       IN->score_aln=OUT->score_aln=(max_aln==0)?0:((score_aln*100)/max_aln);
 872       for ( a=0; a< OUT->nseq; a++)
 873         {
 874           OUT->score_seq[a]=(max_seq[a]==0)?0:((score_seq[a]*100)/max_seq[a]);
 875         }
 876
 877       free_int (lu,-1);
 878       free_int (pos , -1);
 879       vfree ( score_seq);
 880       vfree ( max_seq);
 881       return OUT;
 882     }
 883 Alignment * fast_coffee_evaluate_output ( Alignment *IN,Constraint_list *CL)
 884     {
 885     int a,b, c, m,res, s, s1, s2, r1, r2;
 886     Alignment *OUT=NULL;
 887     int **pos, **pos2;
 888
 889     double score_col=0, score_aln=0, score_res=0;
 890     double max_col, max_aln;
 891     double *max_seq, *score_seq;
 892     int local_m;
 893     int local_nseq;
 894
 895
 896     /*NORMALIZE: with the highest scoring pair found in the multiple alignment*/
 897
 898
 899     if ( !CL->evaluate_residue_pair){fprintf ( stderr, "\nWARNING: CL->evaluate_residue_pair Not set\nSet to: extend_residue_pair\n");CL->evaluate_residue_pair= extend_residue_pair; }
 900
 901     OUT=copy_aln (IN, OUT);
 902     pos=aln2pos_simple(IN, IN->nseq);
 903     pos2=aln2defined_residues (IN, CL);
 904
 905     max_seq=vcalloc ( IN->nseq, sizeof (double));
 906     score_seq=vcalloc ( IN->nseq, sizeof (double));
 907
 908
 909
 910     /*1: Identify the highest scoring pair within the alignment*/
 911
 912     for ( m=0, a=0; a< IN->len_aln; a++)
 913         {
 914             for ( b=0; b< IN->nseq; b++)
 915                 {
 916                 s1=IN->order[b][0];
 917                 r1=pos[b][a];
 918
 919
 920                 for ( c=0; c< IN->nseq; c++)
 921                     {
 922                     s2=IN->order[c][0];
 923                     r2=pos[c][a];
 924                     if ( s1==s2 && !CL->do_self)continue;
 925
 926                     if ( s1< s2)s=(CL->evaluate_residue_pair)( CL, s1, r1, s2, r2);
 927                     else        s=(CL->evaluate_residue_pair)( CL, s2, r2, s1, r1);
 928
 929                     s=(s!=UNDEFINED)?s:0;
 930                     m=MAX(m, s);
 931                     }
 932                 }
 933         }
 934
 935     local_m=m;
 936
 937     sprintf ( OUT->name[IN->nseq], "cons");
 938     for ( max_aln=0,score_aln=0,a=0; a< IN->len_aln; a++)
 939         {
 940         OUT->seq_al[IN->nseq][a]=NO_COLOR_RESIDUE;
 941         for ( local_nseq=0,b=0; b<IN->nseq; b++){local_nseq+=(pos[b][a]>0 && pos2[b][a])?1:0;}
 942         local_m=m*(local_nseq-1);
 943
 944         for ( max_col=0, score_col=0,b=0; b< IN->nseq; b++)
 945             {
 946             OUT->seq_al[b][a]=NO_COLOR_RESIDUE;
 947             s1=IN->order[b][0];
 948             r1=pos[b][a];
 949
 950             if (r1<=0 || !pos2[b][a])
 951               {
 952                 continue;
 953               }
 954
 955             for ( score_res=0,c=0; c< IN->nseq; c++)
 956                 {
 957                     s2=IN->order[c][0];
 958                     r2=pos[c][a];
 959
 960                     if ((s1==s2 && !CL->do_self) || r2<=0 || !pos2[c][a]){continue;}
 961                     max_col   +=m;
 962                     max_seq[b]+=m;
 963                     max_aln   +=m;
 964
 965                     if ( s1< s2)s=(CL->evaluate_residue_pair)( CL, s1, r1, s2, r2);
 966                     else        s=(CL->evaluate_residue_pair)( CL, s2, r2, s1, r1);
 967                     s=(s!=UNDEFINED)?s:0;
 968
 969                     score_res+=s;
 970                     score_col+=s;
 971                     score_seq[b]+=s;
 972                     score_aln+=s;
 973                 }
 974
 975             res=(local_m==0)?NO_COLOR_RESIDUE:((score_res*10)/local_m);
 976             (OUT)->seq_al[b][a]=(res==NO_COLOR_RESIDUE)?res:(MIN(res,9));
 977
 978
 979             }
 980
 981         res=(max_col==0)?NO_COLOR_RESIDUE:((score_col*10)/max_col);
 982         OUT->seq_al[IN->nseq][a]=(res==NO_COLOR_RESIDUE)?res:(MIN(res,9));
 983
 984         }
 985
 986     IN->score_aln=OUT->score_aln=(max_aln==0)?0:((score_aln*100)/max_aln);
 987     for ( a=0; a< OUT->nseq; a++)
 988         {
 989         OUT->score_seq[a]=(max_seq[a]==0)?0:((score_seq[a]*100)/max_seq[a]);
 990         }
 991
 992     free_int (pos , -1);
 993     free_int (pos2, -1);
 994
 995     vfree ( score_seq);
 996     vfree ( max_seq);
 997     return OUT;
 998     }
 999
1000 Alignment * slow_coffee_evaluate_output ( Alignment *IN,Constraint_list *CL)
1001     {
1002     int a,b, c,res, s, s1, s2, r1, r2;
1003     Alignment *OUT=NULL;
1004     int **pos, **pos2;
1005     double max_score_r, score_r, max;
1006     double score_col=0, score_aln=0;
1007     double max_score_col, max_score_aln;
1008     double *max_score_seq, *score_seq;
1009     int ***res_extended_weight;
1010     int n_res_in_col;
1011
1012
1013     /*
1014       Residue x: sum of observed extended X.. /sum of possible X..
1015     */
1016
1017
1018
1019
1020     if ( !CL->evaluate_residue_pair){fprintf ( stderr, "\nWARNING: CL->evaluate_residue_pair Not set\nSet to: extend_residue_pair\n");CL->evaluate_residue_pair= extend_residue_pair; }
1021
1022
1023     OUT=copy_aln (IN, OUT);
1024     pos=aln2pos_simple(IN, IN->nseq);
1025     pos2=aln2defined_residues (IN, CL);
1026
1027     max_score_seq=vcalloc ( IN->nseq, sizeof (double));
1028     score_seq=vcalloc ( IN->nseq, sizeof (double));
1029     res_extended_weight=declare_arrayN(3,sizeof(int), (CL->S)->nseq, (CL->S)->max_len+1, 2);
1030     max=(CL->normalise)?(100*CL->normalise)*SCORE_K:100;
1031
1032     for (a=0; a< IN->len_aln; a++)
1033         {
1034           for ( b=0; b< IN->nseq-1; b++)
1035             {
1036               s1=IN->order[b][0];
1037               r1=pos[b][a];
1038               for ( c=b+1; c< IN->nseq; c++)
1039                 {
1040                   s2=IN->order[c][0];
1041                   r2=pos[c][a];
1042                   if ( s1==s2 && !CL->do_self)continue;
1043                   else if ( r1<=0 || r2<=0)   continue;
1044                   else
1045                     {
1046                       s=(CL->evaluate_residue_pair)( CL, s1, r1, s2, r2);
1047                       res_extended_weight[s1][r1][0]+=s*100;
1048                       res_extended_weight[s2][r2][0]+=s*100;
1049                       res_extended_weight[s1][r1][1]+=max;
1050                       res_extended_weight[s2][r2][1]+=max;
1051                     }
1052                 }
1053             }
1054         }
1055
1056
1057     sprintf ( OUT->name[IN->nseq], "cons");
1058     for ( max_score_aln=0,score_aln=0,a=0; a< IN->len_aln; a++)
1059       {
1060         OUT->seq_al[IN->nseq][a]=NO_COLOR_RESIDUE;
1061         for ( n_res_in_col=0,b=0; b<IN->nseq; b++){n_res_in_col+=(pos[b][a]>0 && pos2[b][a]>0)?1:0;}
1062         for ( max_score_col=0, score_col=0,b=0; b< IN->nseq; b++)
1063             {
1064             OUT->seq_al[b][a]=NO_COLOR_RESIDUE;
1065             s1=IN->order[b][0];
1066             r1=pos[b][a];
1067             if (r1<=0 || pos2[b][a]<1)continue;
1068             else
1069               {
1070                 max_score_r  =res_extended_weight[s1][r1][1];
1071                 score_r      =res_extended_weight[s1][r1][0];
1072                 if      ( max_score_r==0 && n_res_in_col>1)res=0;
1073                 else if ( n_res_in_col==1)res=NO_COLOR_RESIDUE;
1074                 else res=((score_r*10)/max_score_r);
1075
1076
1077                 (OUT)->seq_al[b][a]=(res==NO_COLOR_RESIDUE)?res:(MIN(res, 9));
1078                 max_score_col+=max_score_r;
1079                     score_col+=score_r;
1080                 max_score_seq[b]+=max_score_r;
1081                     score_seq[b]+=score_r;
1082                 max_score_aln+=max_score_r;
1083                     score_aln+=score_r;
1084               }
1085             if      ( max_score_col==0 && n_res_in_col>1)res=0;
1086             else if ( n_res_in_col<2)res=NO_COLOR_RESIDUE;
1087             else res=((score_col*10)/max_score_col);
1088
1089             OUT->seq_al[IN->nseq][a]=(res==NO_COLOR_RESIDUE)?res:(MIN(res,9));
1090             }
1091         }
1092     IN->score_aln=OUT->score_aln=(max_score_aln==0)?0:((score_aln*100)/max_score_aln);
1093     for ( a=0; a< OUT->nseq; a++)
1094       {
1095         OUT->score_seq[a]=(max_score_seq[a]==0)?0:((score_seq[a]*100)/max_score_seq[a]);
1096       }
1097
1098
1099     vfree ( score_seq);
1100     vfree ( max_score_seq);
1101     free_arrayN((void*)res_extended_weight, 3);
1102
1103
1104     free_int (pos, -1);
1105     free_int (pos2, -1);
1106     return OUT;
1107     }
1108
1109 double genepred2sd    (Sequence *S);
1110 double genepred2avg   (Sequence *S);
1111 double genepred2zsum2 (Sequence *S);
1112 int *  genepred2orf_len (Sequence *S);
1113 double genepred2acc (Sequence *S, Sequence *PS);
1114 double seq_genepred2acc (Sequence *S, Sequence *PS, char *name);
1115 Alignment *coffee_seq_evaluate_output ( Alignment *IN, Constraint_list *CL)
1116 {
1117   int **w, a,b,c,l,i;
1118   int avg, min_avg, best_cycle;
1119   char **pred_list,**weight_list;
1120   FILE *fp;
1121   Sequence *S, *PS;
1122   int ncycle=20;
1123   float nsd, best_score,score, acc;
1124   int min_ne;
1125   Alignment *A;
1126   float *count;
1127   min_ne=CL->ne/10; //avoid side effects due to very relaxed libraries
1128   w=NULL;
1129   pred_list=declare_char (ncycle, 100);
1130   weight_list=declare_char(ncycle,100);
1131   S=CL->S;
1132   fprintf ( stderr, "\nSCORE_MODE: %s\n", CL->genepred_score);
1133   for (a=0; a<ncycle && CL->ne>min_ne; a++)
1134     {
1135       if (w);free_int (w, -1);
1136       w=list2residue_total_weight (CL);
1137       PS=dnaseq2geneseq (S, w);
1138       nsd=(float)(genepred2sd(PS)/genepred2avg(PS));
1139       if   ( strm (CL->genepred_score, "nsd"))score=nsd;
1140       else score=1-seq_genepred2acc (S, PS, CL->genepred_score);
1141       fprintf ( stderr, "Cycle: %2d LIB: %6d AVG LENGTH: %6.2f NSD: %8.2f SCORE: %.2f ",a+1,CL->ne,(float) genepred2avg(PS),nsd,score);
1142
1143       vfree (display_accuracy (genepred_seq2accuracy_counts (S, PS, NULL),stderr));
1144       pred_list[a]=vtmpnam(NULL);
1145       weight_list[a]=vtmpnam(NULL);
1146       output_fasta_seqS(pred_list[a],PS);
1147       output_array_int (weight_list[a],w);
1148       free_sequence (PS, -1);
1149
1150       if (a==0 || score<best_score)
1151         {
1152           best_score=score;
1153           best_cycle=a;
1154         }
1155       CL=relax_constraint_list_4gp(CL);
1156     }
1157   if (w)free_int(w,-1);
1158   S=get_fasta_sequence (pred_list[best_cycle], NULL);
1159
1160   fprintf ( stderr, "\nSelected Cycle: %d (SCORE=%.4f)----> ", best_cycle+1, best_score);
1161   vfree (display_accuracy (genepred_seq2accuracy_counts ((CL->S), S, NULL),stderr));
1162
1163   genepred_seq2accuracy_counts4all((CL->S), S);
1164
1165   IN=seq2aln (S,NULL, 1);
1166   IN->score_res=input_array_int (weight_list[best_cycle]);
1167   return IN;
1168 }
1169
1170 double seq_genepred2acc (Sequence *S, Sequence *PS, char *name)
1171 {
1172   float *count;
1173   float *acc;
1174   float r;
1175   int ref, target;
1176
1177   if ( strm (name, "best"))return genepred2acc (S, PS);
1178
1179   ref=name_is_in_list (name, S->name, S->nseq, 100);
1180   target=name_is_in_list (name, PS->name, PS->nseq, 100);
1181
1182   if ( target==-1 || ref==-1)
1183     {
1184       printf_exit (EXIT_FAILURE,stderr, "\nERROR: %s is not a valid sequence", name);
1185     }
1186   count=genepred2accuracy_counts      (S->seq[ref],PS->seq[target],NULL);
1187   acc=counts2accuracy (count);
1188
1189   r=acc[3];
1190   vfree (acc);
1191   vfree (count);
1192   return r;
1193 }
1194
1195
1196
1197 double genepred2acc (Sequence *S, Sequence *PS)
1198 {
1199   float *count;
1200   float *acc;
1201   float r;
1202   count=genepred_seq2accuracy_counts (S, PS, NULL);
1203   acc=counts2accuracy (count);
1204
1205   r=acc[3];
1206   vfree (acc);
1207   vfree (count);
1208   return r;
1209 }
1210
1211 double genepred2sd (Sequence *S)
1212 {
1213   double sum=0, sum2=0;
1214   int a, b;
1215   int *len;
1216   len=genepred2orf_len (S);
1217
1218   for (a=0; a<S->nseq; a++)
1219     {
1220       sum+=(double)len[a];
1221       sum2+=(double)len[a]*(double)len[a];
1222     }
1223   sum/=(double)S->nseq;
1224   sum2/=S->nseq;
1225   vfree (len);
1226   return sqrt (sum2-(sum*sum));
1227 }
1228 double genepred2avg (Sequence *S)
1229 {
1230   int a, b;
1231   double avg=0;
1232   int *len;
1233   len=genepred2orf_len (S);
1234
1235   for (a=0; a<S->nseq; a++)avg+=len[a];
1236   vfree (len);
1237   return avg/(double)S->nseq;
1238 }
1239 double genepred2zsum2 (Sequence *S)
1240 {
1241   double sum=0, sum2=0, zscore=0, zsum2=0;
1242   int a, b;
1243   int *len;
1244   double sd, avg;
1245   sd=genepred2sd(S);
1246   avg=genepred2avg (S);
1247
1248   len=genepred2orf_len (S);
1249   for (a=0; a<S->nseq; a++)
1250     {
1251       zscore=((double)len[a]-avg)/sd;
1252       zsum2+=FABS(zscore);
1253     }
1254   zsum2/=(float)S->nseq;
1255   vfree (len);
1256   return zsum2;
1257 }
1258 int *genepred2orf_len (Sequence *S)
1259 {
1260   int a,b, *len;
1261   len=vcalloc (S->nseq, sizeof (int));
1262    for (a=0; a<S->nseq; a++)
1263     for (b=0; b<S->len[a]; b++)
1264       len[a]+=(isupper(S->seq[a][b]));
1265    return len;
1266 }
1267
1268 Alignment *coffee_seq_evaluate_output_old2 ( Alignment *IN, Constraint_list *CL)
1269 {
1270   int **w, a,b,c,l,i;
1271   int avg, min_avg, best_cycle;
1272   char **pred_list; FILE *fp;
1273   Sequence *S, *PS;
1274   int ncycle=100;
1275   w=NULL;
1276   pred_list=declare_char (ncycle, 100);
1277   S=CL->S;
1278
1279   //CL=expand_constraint_list_4gp(CL);
1280   min_avg=constraint_list2avg(CL);
1281
1282   for (a=1; a<ncycle && CL->ne>0; a++)
1283     {
1284       int t;
1285       if (w);free_int (w, -1);
1286       w=list2residue_total_weight (CL);
1287       CL=relax_constraint_list_4gp (CL);
1288       fprintf (stderr,"\nRELAX CYCLE: %2d AVG: %5d [%10d] ", a, avg=constraint_list2avg(CL), CL->ne);
1289       for (b=0; b<S->nseq; b++)for (c=0; c<S->len[b]; c++)w[b][c]-=1;
1290
1291       //rescale nuclotide coding weights
1292
1293       PS=dnaseq2geneseq (S, w);
1294       vfree (display_accuracy (genepred_seq2accuracy_counts (S, PS, NULL),stderr));
1295
1296
1297       free_sequence (PS, PS->nseq);
1298       if ( avg<min_avg)
1299         {
1300           min_avg=avg;
1301           best_cycle=a;
1302         }
1303     }
1304
1305   S=get_fasta_sequence (pred_list[best_cycle], NULL);
1306   vfree (display_accuracy (genepred_seq2accuracy_counts ((CL->S), S, NULL),stderr));myexit (0);
1307   for (a=0; a<S->nseq; a++)
1308     HERE (">%s\n%s", S->name[a], S->seq[a]);
1309   myexit (0);
1310   return seq2aln (S,NULL, 1);
1311 }
1312
1313
1314
1315 Alignment * non_extended_t_coffee_evaluate_output ( Alignment *IN,Constraint_list *CL)
1316     {
1317     int a,b, c,res, s1, s2, r1, r2;
1318     Alignment *OUT=NULL;
1319     int **pos;
1320     int max_score_r, score_r;
1321     double score_col=0, score_aln=0;
1322     int max_score_col, max_score_aln;
1323     double *max_score_seq, *score_seq;
1324     int local_nseq;
1325     int **tot_non_extended_weight;
1326     int **res_non_extended_weight;
1327     int *l;
1328     CLIST_TYPE  *entry=NULL;
1329     int p;
1330     int max_score=0;
1331
1332     entry=vcalloc (CL->entry_len+1, CL->el_size);
1333     if ( !CL->evaluate_residue_pair){fprintf ( stderr, "\nWARNING: CL->evaluate_residue_pair Not set\nSet to: extend_residue_pair\n");CL->evaluate_residue_pair= extend_residue_pair; }
1334
1335     OUT=copy_aln (IN, OUT);
1336     pos=aln2pos_simple(IN, IN->nseq);
1337
1338
1339     max_score_seq=vcalloc ( IN->nseq, sizeof (double));
1340     score_seq=vcalloc ( IN->nseq, sizeof (double));
1341
1342     tot_non_extended_weight=list2residue_total_weight(CL);
1343     res_non_extended_weight=declare_int ((CL->S)->nseq, (CL->S)->max_len+1);
1344
1345     for (a=0; a< IN->len_aln; a++)
1346         {
1347             for ( b=0; b< IN->nseq-1; b++)
1348                 {
1349                 s1=IN->order[b][0];
1350                 r1=pos[b][a];
1351                 for ( c=b+1; c< IN->nseq; c++)
1352                     {
1353                     s2=IN->order[c][0];
1354                     r2=pos[c][a];
1355                     if ( s1==s2 && !CL->do_self)continue;
1356                     else if ( r1<=0 || r2<=0)   continue;
1357                     else
1358                       {
1359                         entry[SEQ1]=s1;
1360                         entry[SEQ2]=s2;
1361                         entry[R1]=r1;
1362                         entry[R2]=r2;
1363                         if ((l=main_search_in_list_constraint (entry,&p,4,CL))!=NULL)
1364                           {
1365                             res_non_extended_weight[s1][r1]+=l[WE];
1366                             res_non_extended_weight[s2][r2]+=l[WE];
1367                           }
1368                         entry[SEQ1]=s2;
1369                         entry[SEQ2]=s1;
1370                         entry[R1]=r2;
1371                         entry[R2]=r1;
1372                         if ((l=main_search_in_list_constraint (entry,&p,4,CL))!=NULL)
1373                             {
1374                               res_non_extended_weight[s1][r1]+=l[WE];
1375                               res_non_extended_weight[s2][r2]+=l[WE];
1376                             }
1377                         max_score=MAX(max_score,res_non_extended_weight[s1][r1]);
1378                         max_score=MAX(max_score,res_non_extended_weight[s2][r2]);
1379
1380                       }
1381                     }
1382                 }
1383         }
1384
1385     sprintf ( OUT->name[IN->nseq], "cons");
1386     for ( max_score_aln=0,score_aln=0,a=0; a< IN->len_aln; a++)
1387         {
1388         OUT->seq_al[IN->nseq][a]=NO_COLOR_RESIDUE;
1389         for ( local_nseq=0,b=0; b<IN->nseq; b++){local_nseq+=(pos[b][a]>0)?1:0;}
1390
1391         for ( max_score_col=0, score_col=0,b=0; b< IN->nseq; b++)
1392             {
1393             OUT->seq_al[b][a]=NO_COLOR_RESIDUE;
1394             s1=IN->order[b][0];
1395             r1=pos[b][a];
1396             if (r1<=0)continue;
1397             else
1398               {
1399                 max_score_r  =max_score;/*tot_non_extended_weight[s1][r1];*/
1400                 score_r=res_non_extended_weight[s1][r1];
1401                 res=(max_score_r==0 || local_nseq<2 )?NO_COLOR_RESIDUE:((score_r*10)/max_score_r);
1402
1403                 (OUT)->seq_al[b][a]=(res==NO_COLOR_RESIDUE)?res:(MIN(res, 9));
1404                 max_score_col+=max_score_r;
1405                     score_col+=score_r;
1406                 max_score_seq[b]+=max_score_r;
1407                     score_seq[b]+=score_r;
1408                 max_score_aln+=max_score_r;
1409                     score_aln+=score_r;
1410               }
1411             res=(max_score_col==0 || local_nseq<2)?NO_COLOR_RESIDUE:((score_col*10)/max_score_col);
1412             OUT->seq_al[IN->nseq][a]=(res==NO_COLOR_RESIDUE)?res:(MIN(res,9));
1413             }
1414         }
1415     IN->score_aln=OUT->score_aln=(max_score_aln==0)?0:((score_aln*100)/max_score_aln);
1416     for ( a=0; a< OUT->nseq; a++)
1417       {
1418         OUT->score_seq[a]=(max_score_seq[a]==0)?0:((score_seq[a]*100)/max_score_seq[a]);
1419         OUT->score_seq[a]=(OUT->score_seq[a]>100)?100:OUT->score_seq[a];
1420       }
1421     OUT->score_aln=(OUT->score_aln>100)?100:OUT->score_aln;
1422
1423     vfree ( score_seq);
1424     vfree ( max_score_seq);
1425
1426     free_int (tot_non_extended_weight, -1);
1427     free_int (res_non_extended_weight, -1);
1428     vfree(entry);
1429     free_int (pos, -1);
1430
1431     return OUT;
1432     }
1433
1434
1435 /*********************************************************************************************/
1436 /*                                                                                           */
1437 /*        PROFILE/PRofile Functions                                                          */
1438 /*                                                                                           */
1439 /*********************************************************************************************/
1440 int channel_profile_profile (int *prf1, int *prf2, Constraint_list *CL);
1441
1442 Profile_cost_func get_profile_mode_function (char *name, Profile_cost_func func)
1443 {
1444   int a;
1445   static int nfunc;
1446   static Profile_cost_func *flist;
1447   static char **nlist;
1448
1449
1450
1451   /*The first time: initialize the list of pairwse functions*/
1452   /*If func==NULL:REturns a pointer to the function associated with a name*/
1453   /*If name is empty:Prints the name of the function  associated with name*/
1454
1455     if ( nfunc==0)
1456       {
1457         flist=vcalloc ( 100, sizeof (Pwfunc));
1458         nlist=declare_char (100, 100);
1459
1460         flist[nfunc]=cw_profile_profile;
1461         sprintf (nlist[nfunc], "cw_profile_profile");
1462         nfunc++;
1463
1464         flist[nfunc]=muscle_profile_profile;
1465         sprintf (nlist[nfunc], "muscle_profile_profile");
1466         nfunc++;
1467
1468         flist[nfunc]=channel_profile_profile;
1469         sprintf (nlist[nfunc], "channel_profile_profile");
1470         nfunc++;
1471       }
1472
1473    for ( a=0; a<nfunc; a++)
1474       {
1475         if ( (name && strm (nlist[a],name)) || flist[a]==func)
1476              {
1477                if (name)sprintf (name,"%s", nlist[a]);
1478                return flist[a];
1479              }
1480       }
1481     fprintf ( stderr, "\n[%s] is an unknown profile_profile function[FATAL:%s]\n",name, PROGRAM);
1482     crash ( "\n");
1483     return NULL;
1484   }
1485
1486 int generic_evaluate_profile_score     (Constraint_list *CL,Alignment *Profile1,int s1, int r1, Alignment *Profile2,int s2, int r2, Profile_cost_func prf_prf)
1487     {
1488       int *prf1, *prf2;
1489       static int *dummy;
1490       int score;
1491
1492
1493       /*Generic profile function*/
1494       if( !dummy)
1495         {
1496          dummy=vcalloc (10, sizeof(int));
1497          dummy[0]=1;/*Number of Amino acid types on colum*/
1498          dummy[1]=5;/*Length of Dummy*/
1499          dummy[3]='\0';/*Amino acid*/
1500          dummy[4]=1; /*Number of occurences*/
1501          dummy[5]=100; /*Frequency in the MSA column*/
1502
1503         }
1504
1505       if (r1>0 && r2>0)
1506             {
1507             r1--;
1508             r2--;
1509
1510             prf1=(Profile1)?(Profile1->P)->count2[r1]:NULL;
1511             prf2=(Profile2)?(Profile2->P)->count2[r2]:NULL;
1512
1513             if (!prf1)     {prf1=dummy; prf1[3]=(CL->S)->seq[s1][r1];}
1514             else if (!prf2){prf2=dummy; prf2[3]=(CL->S)->seq[s2][r2];}
1515
1516             score=((prf_prf==NULL)?cw_profile_profile:prf_prf) (prf1, prf2, CL);
1517             return score;
1518             }
1519       else
1520         return 0;
1521     }
1522
1523 int cw_profile_profile_count    (int *prf1, int *prf2, Constraint_list *CL)
1524     {
1525       /*see function aln2count2 for prf structure*/
1526       int a, b, n;
1527       int res1, res2;
1528       double score=0;
1529
1530
1531       for ( n=0,a=3; a<prf1[1]; a+=3)
1532         for ( b=3; b<prf2[1]; b+=3)
1533           {
1534
1535             res1=prf1[a];
1536             res2=prf2[b];
1537
1538             score+=prf1[a+1]*prf2[b+1]*CL->M[res1-'A'][res2-'A'];
1539             n+=prf1[a+1]*prf2[b+1];
1540           }
1541
1542
1543       score=(score*SCORE_K)/n;
1544       return score;
1545     }
1546 int muscle_profile_profile    (int *prf1, int *prf2, Constraint_list *CL)
1547     {
1548       /*see function aln2count2 for prf structure*/
1549       int a, b;
1550       int res1, res2;
1551       double score=0, fg1, fg2, fi, fj, m;
1552       static double *exp_lu;
1553       if (exp_lu==NULL)
1554         {
1555           exp_lu=vcalloc ( 10000, sizeof (double));
1556           exp_lu+=2000;
1557           for ( a=-1000; a<1000; a++)
1558             exp_lu[a]=exp((double)a);
1559         }
1560
1561
1562
1563       for (a=3; a<prf1[1]; a+=3)
1564         {
1565           res1=prf1[a];
1566           fi=(double)prf1[a+2]/100;
1567
1568           for ( b=3; b<prf2[1]; b+=3)
1569             {
1570               res2=prf2[b];
1571               fj=(double)prf2[b+2]/100;
1572               /*m=exp((double)CL->M[res1-'A'][res2-'A']);*/
1573               m=exp_lu[CL->M[res1-'A'][res2-'A']];
1574               score+=m*fi*fj;
1575             }
1576         }
1577
1578       fg1=(double)prf1[2]/100;
1579       fg2=(double)prf2[2]/100;
1580       score=(score==0)?0:log(score)*(1-fg1)*(1-fg2);
1581       score=(score*SCORE_K);
1582       /*if ( score<-100)fprintf ( stderr, "\nSCORE %d %d", (int)score, cw_profile_profile(prf1, prf2, CL));*/
1583
1584       return (int)score;
1585     }
1586 int cw_profile_profile    (int *prf1, int *prf2, Constraint_list *CL)
1587     {
1588       /*see function aln2count2 for prf structure*/
1589       int a, b, n,p;
1590       int res1, res2;
1591       double score=0;
1592
1593
1594       for ( n=0,a=3; a<prf1[1]; a+=3)
1595         for ( b=3; b<prf2[1]; b+=3)
1596           {
1597
1598             res1=prf1[a];
1599             res2=prf2[b];
1600             p=prf1[a+1]*prf2[b+1];
1601
1602             n+=p;
1603             score+=p*CL->M[res1-'A'][res2-'A'];
1604           }
1605       score=(score*SCORE_K)/((double)(n==0)?1:n);
1606       return score;
1607     }
1608 int cw_profile_profile_old    (int *prf1, int *prf2, Constraint_list *CL)
1609     {
1610       /*see function aln2count2 for prf structure*/
1611       int a, b, n,p;
1612       int res1, res2;
1613       double score=0;
1614
1615
1616
1617       for ( n=0,a=3; a<prf1[1]; a+=3)
1618         for ( b=3; b<prf2[1]; b+=3)
1619           {
1620
1621             res1=prf1[a];
1622             res2=prf2[b];
1623             p=prf1[a+1]*prf2[b+1];
1624
1625             n+=p;
1626             score+=p*CL->M[res1-'A'][res2-'A'];
1627           }
1628       score=(score*SCORE_K)/((double)(n==0)?1:n);
1629       return score;
1630     }
1631 int channel_profile_profile ( int *prf1, int *prf2, Constraint_list *CL)
1632 {
1633
1634   int score=0;
1635
1636   prf1+=prf1[1];
1637   prf2+=prf2[1];
1638
1639
1640   if (prf1[0]!=prf1[0]){fprintf ( stderr, "\nERROR: Inconsistent number of channels [channel_profile_profile::FATAL%s]", PROGRAM);}
1641   else
1642     {
1643       int a, n;
1644       for (a=1, n=0; a<=prf1[0]; a++)
1645         {
1646           if (prf1[a]>0 && prf2[a]>0)
1647             {
1648               n++;score+=CL->M[prf1[a]-'A'][prf2[a]-'A'];
1649
1650             }
1651         }
1652
1653       if ( n==0)return 0;
1654
1655       score=(n==0)?0:(score*SCORE_K)/n;
1656
1657     }
1658   return score;
1659 }
1660
1661 /*********************************************************************************************/
1662 /*                                                                                           */
1663 /*         FUNCTIONS FOR GETING THE COST : (Sequences) ->evaluate_residue_pair               */
1664 /*                                                                                           */
1665 /*********************************************************************************************/
1666 int initialize_scoring_scheme (Constraint_list *CL)
1667 {
1668   if (!CL) return 0;
1669   else if (!CL->evaluate_residue_pair)return 0;
1670   else if ( !CL->S) return 0;
1671   else if ( (CL->S)->nseq<2) return 0;
1672   else if ( strlen ((CL->S)->seq[0])==0)return 0;
1673   else if ( strlen ((CL->S)->seq[1])==0)return 0;
1674   else
1675     {
1676       //CL->evaluate_residue_pair (CL,13,1,16,1);
1677       CL->evaluate_residue_pair (CL,0,1,1,1);
1678
1679     }
1680   return 1;
1681 }
1682
1683 int evaluate_blast_profile_score (Constraint_list *CL, int s1, int r1, int s2, int r2)
1684 {
1685   Alignment *PRF1;
1686   Alignment *PRF2;
1687
1688
1689   PRF1=(Alignment*)atop(seq2T_value (CL->S, s1, "A", "_RB_"));
1690   PRF2=(Alignment*)atop(seq2T_value (CL->S, s2, "A", "_RB_"));
1691
1692   return generic_evaluate_profile_score     (CL,PRF1,s1, r1, PRF2,s2, r2, CL->profile_mode);
1693 }
1694
1695 int evaluate_aln_profile_score (Constraint_list *CL, int s1, int r1, int s2, int r2)
1696 {
1697
1698   return generic_evaluate_profile_score   (CL,seq2R_template_profile((CL->S),s1),s1, r1, seq2R_template_profile(CL->S,s2),s2, r2, CL->profile_mode);
1699 }
1700
1701
1702 int evaluate_profile_score     (Constraint_list *CL,Alignment *Prf1, int s1, int r1, Alignment *Prf2,int s2, int r2)
1703 {
1704
1705   return generic_evaluate_profile_score (CL, Prf1, s1,r1,Prf2, s2,r2,CL->profile_mode);
1706 }
1707
1708 int evaluate_cdna_matrix_score (Constraint_list *CL, int s1, int r1, int s2, int r2)
1709     {
1710       char a1, a2;
1711
1712       if (r1>0 && r2>0)
1713        {
1714          r1--;
1715          r2--;
1716
1717          a1=translate_dna_codon((CL->S)->seq[s1]+r1,'x');
1718          a2=translate_dna_codon((CL->S)->seq[s2]+r2,'x');
1719
1720
1721
1722          if (a1=='x' || a2=='x')return 0;
1723          else return CL->M[a1-'A'][a2-'A']*SCORE_K;
1724        }
1725       else
1726         {
1727           return 0;
1728         }
1729     }
1730 int evaluate_physico_score ( Constraint_list *CL, int s1, int r1, int s2, int r2)
1731 {
1732   int a, b, p;
1733   double tot;
1734   static float **prop_table;
1735   static int n_prop;
1736   static double max;
1737   if (r1<0 || r2<0)return 0;
1738   if ( !prop_table)
1739     {
1740       prop_table= initialise_aa_physico_chemical_property_table(&n_prop);
1741       for (p=0; p< n_prop; p++)max+=100;
1742       max=sqrt(max);
1743     }
1744   a=tolower (( CL->S)->seq[s1][r1]);
1745   b=tolower (( CL->S)->seq[s2][r2]);
1746
1747   for (tot=0,p=0; p< n_prop; p++)
1748     {
1749       tot+=(double)(prop_table[p][a]-prop_table[p][b])*(prop_table[p][a]-prop_table[p][b]);
1750     }
1751
1752   tot=(sqrt(tot)/max)*10;
1753
1754   return (int) tot*SCORE_K;
1755 }
1756
1757
1758
1759 int evaluate_diaa_matrix_score ( Constraint_list *CL, int s1, int r1, int s2, int r2)
1760     {
1761
1762       static int ****m;
1763       static int *alp;
1764
1765       if (m==NULL)
1766         {
1767           FILE *fp;
1768           char k1[2], k2[2];
1769           int v1, v2, c;
1770           char *buf=NULL;
1771           int a;
1772
1773           m=declare_arrayN(4, sizeof (int), 26, 26, 26, 26);
1774           fp=vfopen ("diaa_mat.mat", "r");
1775           while ((c=fgetc (fp))!=EOF)
1776             {
1777
1778               ungetc (c, fp);
1779               buf=vfgets(buf, fp);
1780
1781               if (c=='#');
1782               else
1783                 {
1784                   sscanf (buf, "%s %s %d %d", k1, k2, &v1, &v2);
1785
1786                   m[k1[0]-'a'][k1[1]-'a'][k2[0]-'a'][k2[1]-'a']=v1;
1787                   m[k2[0]-'a'][k2[1]-'a'][k1[0]-'a'][k1[1]-'a']=v1;
1788                 }
1789             }
1790           vfclose (fp);
1791           alp=vcalloc (256, sizeof (int));
1792           for (a=0; a<26; a++)alp[a+'a']=1;
1793           alp['b']=0;
1794           alp['j']=0;
1795           alp['o']=0;
1796           alp['u']=0;
1797           alp['x']=0;
1798           alp['z']=0;
1799         }
1800
1801
1802       if (r1>0 && r2>0)
1803           {
1804             int s=0, n=0;
1805             char aa1, aa2, aa3, aa4, u;
1806
1807             r1--;
1808             r2--;
1809
1810             if (r1>0 && r2>0)
1811               {
1812                 aa1=tolower((CL->S)->seq[s1][r1-1]);
1813                 aa2=tolower((CL->S)->seq[s1][r1]);
1814                 aa3=tolower((CL->S)->seq[s2][r2-1]);
1815                 aa4=tolower((CL->S)->seq[s2][r2]);
1816                 u=alp[(int)aa1];u+=alp[(int)aa2];u+=alp[(int)aa3];u+=alp[(int)aa4];
1817                 if (u==4)
1818                   {
1819                     s+=m[aa1-'a'][aa2-'a'][aa3-'a'][aa4-'a'];
1820                     n++;
1821                   }
1822               }
1823
1824             aa1=tolower((CL->S)->seq[s1][r1]);
1825             aa2=tolower((CL->S)->seq[s1][r1+1]);
1826             aa3=tolower((CL->S)->seq[s2][r2]);
1827             aa4=tolower((CL->S)->seq[s2][r2+1]);
1828             u=alp[(int)aa1];u+=alp[(int)aa2];u+=alp[(int)aa3];u+=alp[(int)aa4];
1829             if (u==4)
1830               {
1831                 s+=m[aa1-'a'][aa2-'a'][aa3-'a'][aa4-'a'];
1832                 n++;
1833               }
1834             if (n)return (s*SCORE_K)/n;
1835             else return 0;
1836             }
1837       return 0;}
1838 int evaluate_monoaa_matrix_score ( Constraint_list *CL, int s1, int r1, int s2, int r2)
1839     {
1840
1841       static int **m;
1842       static int *alp;
1843
1844       if (m==NULL)
1845         {
1846           FILE *fp;
1847           char k1[2], k2[2];
1848           int v1, v2, c;
1849           char *buf=NULL;
1850           int a;
1851
1852           m=declare_arrayN(2, sizeof (int), 26, 26);
1853           fp=vfopen ("monoaa_mat.mat", "r");
1854           while ((c=fgetc (fp))!=EOF)
1855             {
1856
1857               ungetc (c, fp);
1858               buf=vfgets(buf, fp);
1859
1860               if (c=='#');
1861               else
1862                 {
1863                   sscanf (buf, "%s %s %d %d", k1, k2, &v1, &v2);
1864
1865                   m[k1[0]-'a'][k2[0]-'a']=v1;
1866                   m[k2[0]-'a'][k1[0]-'a']=v1;
1867                 }
1868             }
1869           vfclose (fp);
1870           alp=vcalloc (256, sizeof (int));
1871           for (a=0; a<26; a++)alp[a+'a']=1;
1872           alp['b']=0;
1873           alp['j']=0;
1874           alp['o']=0;
1875           alp['u']=0;
1876           alp['x']=0;
1877           alp['z']=0;
1878         }
1879
1880
1881       if (r1>0 && r2>0)
1882           {
1883             int s=0, n=0;
1884             char aa1, aa3, u;
1885
1886             r1--;
1887             r2--;
1888
1889             if (r1>0 && r2>0)
1890               {
1891                 aa1=tolower((CL->S)->seq[s1][r1]);
1892                 aa3=tolower((CL->S)->seq[s2][r2]);
1893                 u=alp[(int)aa1];u+=alp[(int)aa3];
1894                 if (u==2)
1895                   {
1896                     s+=m[aa1-'a'][aa3-'a'];
1897                     n++;
1898                   }
1899               }
1900
1901             if (n)return (s*SCORE_K)/n;
1902             else return 0;
1903             }
1904       return 0;}
1905 int evaluate_matrix_score ( Constraint_list *CL, int s1, int r1, int s2, int r2)
1906     {
1907
1908       if ( seq2R_template_profile (CL->S,s1) ||seq2R_template_profile (CL->S,s2))
1909         {
1910           return evaluate_aln_profile_score ( CL, s1,r1, s2, r2);
1911         }
1912
1913       if (r1>0 && r2>0)
1914             {
1915             r1--;
1916             r2--;
1917
1918             return CL->M[(CL->S)->seq[s1][r1]-'A'][(CL->S)->seq[s2][r2]-'A']*SCORE_K;
1919             }
1920         else
1921             return 0;
1922     }
1923 int *get_curvature ( int s1, Constraint_list *CL);
1924 int evaluate_curvature_score( Constraint_list *CL, int s1, int r1, int s2, int r2)
1925 {
1926   static int **st;
1927   int score;
1928   CL->gop=0;
1929   CL->gep=0;
1930
1931   if (!st) st= vcalloc ((CL->S)->nseq, sizeof (char*));
1932   if (!st[s1])
1933     {
1934       st[s1]=get_curvature (s1, CL);
1935     }
1936   if (!st[s2])
1937     {
1938       st[s2]=get_curvature (s2, CL);
1939     }
1940
1941   if (r1>0 && r2>0)
1942     {
1943       char p1, p2;
1944
1945       r1--;
1946       r2--;
1947
1948       p1=st[s1][r1];
1949       p2=st[s2][r2];
1950
1951       score=p1-p2;
1952       score=FABS(score);
1953       score=20-score;
1954       //HERE ("%d", score);
1955       //if (score<0)HERE ("%d", score);
1956       return score;
1957     }
1958   else
1959     {
1960       return 0;
1961     }
1962
1963 }
1964 int *get_curvature ( int s1, Constraint_list *CL)
1965 {
1966   int *array, n=0, a;
1967   char c;
1968   FILE *fp;
1969   char name [1000], b1[100], b2[100];
1970   float f;
1971   sprintf ( name, "%s.curvature", (CL->S)->name[s1]);
1972   array=vcalloc (strlen ((CL->S)->seq[s1]), sizeof (int));
1973   fp=vfopen ( name, "r");
1974   while ( fscanf (fp, "%s %d %c %f\n",b1, &a, &c,&f )==4)
1975     {
1976       if ( c!=(CL->S)->seq[s1][n]){HERE ("ERROR: %c %c", c,(CL->S)->seq[s1][n] );myexit (0);}
1977       else array[n++]=(int)(float)100*(float)f;
1978     }
1979   vfclose (fp);
1980   return array;
1981 }
1982 int evaluate_tm_matrix_score ( Constraint_list *CL, int s1, int r1, int s2, int r2)
1983 {
1984   static char **st;
1985   float F, RF;
1986
1987   if (!st)
1988     {
1989       st= vcalloc ((CL->S)->nseq, sizeof (char*));
1990       RF=atoigetenv ("TM_FACTOR_4_TCOFFEE");
1991       if ( !RF)RF=10;
1992     }
1993
1994   if (!st[s1])st[s1]=seq2T_template_string((CL->S),s1);
1995   if (!st[s2])st[s2]=seq2T_template_string((CL->S),s2);
1996
1997   if (r1<=0 || r2<=0)return 0;
1998   else if (st[s1] && st[s2] && (st[s1][r1-1]==st[s2][r2-1]))F=RF;
1999   else if (st[s1] && st[s2] && (st[s1][r1-1]!=st[s2][r2-1]))F=0;
2000   else F=0;
2001
2002   return (int)(evaluate_matrix_score (CL, s1, r1, s2, r2))+F;
2003 }
2004 int evaluate_ssp_matrix_score ( Constraint_list *CL, int s1, int r1, int s2, int r2)
2005 {
2006   static char **st;
2007   float F, RF;
2008
2009   if (!st)
2010     {
2011       st= vcalloc ((CL->S)->nseq, sizeof (char*));
2012       RF=atoigetenv ("SSP_FACTOR_4_TCOFFEE");
2013       if ( !RF)RF=10;
2014     }
2015
2016   if (!st[s1])st[s1]=seq2T_template_string((CL->S),s1);
2017   if (!st[s2])st[s2]=seq2T_template_string((CL->S),s2);
2018
2019   if (r1<=0 || r2<=0)return 0;
2020   else if (st[s1] && st[s2] && (st[s1][r1-1]==st[s2][r2-1]))F=RF;
2021   else if (st[s1] && st[s2] && (st[s1][r1-1]!=st[s2][r2-1]))F=0;
2022   else F=0;
2023
2024   return (int)(evaluate_matrix_score (CL, s1, r1, s2, r2))+F;
2025 }
2026
2027
2028 int evaluate_combined_matrix_score ( Constraint_list *CL, int s1, int r1, int s2, int r2)
2029     {
2030       /*
2031         function documentation: start
2032
2033         int evaluate_matrix_score ( Constraint_list *CL, int s1, int s2, int r1, int r2)
2034
2035         this function evaluates the score for matching residue S1(r1) wit residue S2(r2)
2036         using Matrix CL->M;
2037
2038         function documentation: end
2039       */
2040
2041       if ( seq2R_template_profile (CL->S,s1) ||seq2R_template_profile (CL->S,s2))
2042         {
2043           return evaluate_aln_profile_score ( CL, s1,r1, s2, r2);
2044         }
2045
2046       if (r1>0 && r2>0)
2047             {
2048             r1--;
2049             r2--;
2050             if (r1==0 || r2==0)return CL->M[(CL->S)->seq[s1][r1]-'A'][(CL->S)->seq[s2][r2]-'A']*SCORE_K;
2051             else
2052               {
2053                 int A1, A2, B1, B2, a2, b2;
2054                 int score;
2055
2056                 A1=toupper((CL->S)->seq[s1][r1-1]);
2057                 A2=toupper((CL->S)->seq[s1][r1]);
2058                 B1=toupper((CL->S)->seq[s2][r2-1]);
2059                 B2=toupper((CL->S)->seq[s2][r2]);
2060
2061                 a2=tolower(A2);
2062                 b2=tolower(B2);
2063                 A1-='A';A2-='A';B1-='A'; B2-='A';a2-='A';b2-='A';
2064
2065                 score=CL->M[a2][b2]-FABS((CL->M[A1][A2])-(CL->M[B1][B2]));
2066                 score*=SCORE_K;
2067                 return score;
2068               }
2069             }
2070       else
2071         return 0;
2072     }
2073
2074 int residue_pair_non_extended_list ( Constraint_list *CL, int s1, int r1, int s2, int r2 )
2075         {
2076
2077         /*
2078           This is the generic Function->works with everything
2079
2080           int residue_pair_non_extended_list ( Constraint_list *CL, int s1, int r1, int s2, int r2, int field );
2081
2082           Computes the non extended score for aligning residue seq1(r1) Vs seq2(r2)
2083
2084           This function can compare a sequence with itself.
2085
2086           Associated functions: See util constraint list, list extention functions.
2087           function documentation: end
2088         */
2089
2090         int p;
2091
2092
2093         static int *entry;
2094         int *r;
2095         int field;
2096
2097         field = CL->weight_field;
2098
2099
2100         if ( r1<=0 || r2<=0)return 0;
2101         else if ( !CL->extend_jit)
2102            {
2103             if ( !entry) entry=vcalloc (LIST_N_FIELDS , sizeof (int));
2104             entry[SEQ1]=s1;
2105             entry[SEQ2]=s2;
2106             entry[R1]=r1;
2107             entry[R2]=r2;
2108             if ( r1==r2 && s1==s2) return UNDEFINED;
2109             r=main_search_in_list_constraint( entry,&p,4,CL);
2110             if (r==NULL)return 0;
2111             else return r[field]*SCORE_K;
2112            }
2113         else
2114           return UNDEFINED;/*ERROR*/
2115
2116
2117         }
2118
2119
2120
2121 int residue_pair_extended_list_mixt (Constraint_list *CL, int s1, int r1, int s2, int r2 )
2122         {
2123           int score=0;
2124
2125           score+= residue_pair_extended_list_quadruplet(CL, s1, r1, s2, r2);
2126           score+= residue_pair_extended_list (CL, s1, r1, s2, r2);
2127
2128           return score*SCORE_K;
2129         }
2130
2131 int residue_pair_extended_list_quadruplet (Constraint_list *CL, int s1, int r1, int s2, int r2 )
2132         {
2133           double score=0;
2134
2135           int t_s, t_r, t_w, q_s, q_r, q_w;
2136           int a, b;
2137           static int **hasch;
2138
2139           int field;
2140           /* This measure the quadruplets cost on a pair of residues*/
2141
2142
2143
2144           field=CL->weight_field;
2145
2146           if ( r1<=0 || r2<=0)return 0;
2147           if ( !hasch)
2148             {
2149               hasch=vcalloc ( (CL->S)->nseq, sizeof (int*));
2150               for ( a=0; a< (CL->S)->nseq; a++)hasch[a]=vcalloc ( (CL->S)->len[a]+1, sizeof (int));
2151             }
2152
2153
2154           hasch[s1][r1]=100000;
2155           for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
2156             {
2157               t_s=CL->residue_index[s1][r1][a+SEQ2];
2158               t_r=CL->residue_index[s1][r1][a+R2];
2159               t_w=CL->residue_index[s1][r1][a+WE];
2160               if ( CL->seq_for_quadruplet[t_s])
2161                 {
2162                   for ( b=1; b<CL->residue_index[t_s][t_r][0]; b+=ICHUNK)
2163                     {
2164                       q_s=CL->residue_index[t_s][t_r][b+SEQ2];
2165                       q_r=CL->residue_index[t_s][t_r][b+R2];
2166                       q_w=CL->residue_index[t_s][t_r][b+WE];
2167                       if (CL-> seq_for_quadruplet[q_s])
2168                         hasch[q_s][q_r]=MIN(q_w,t_w);
2169
2170                     }
2171                 }
2172             }
2173
2174
2175           for (a=1; a< CL->residue_index[s2][r2][0]; a+=ICHUNK)
2176             {
2177               t_s=CL->residue_index[s2][r2][a+SEQ2];
2178               t_r=CL->residue_index[s2][r2][a+R2];
2179               t_w=CL->residue_index[s2][r2][a+WE];
2180               if ( CL->seq_for_quadruplet[t_s])
2181                 {
2182                   for ( b=1; b<CL->residue_index[t_s][t_r][0]; b+=ICHUNK)
2183                     {
2184                       q_s=CL->residue_index[t_s][t_r][b+SEQ2];
2185                       q_r=CL->residue_index[t_s][t_r][b+R2];
2186                       q_w=CL->residue_index[t_s][t_r][b+WE];
2187                       if (hasch[q_s][q_r] && CL->seq_for_quadruplet[q_s])
2188                         score+=MIN(hasch[q_s][q_r],MIN(q_w,t_w));
2189                     }
2190                 }
2191             }
2192
2193           score=(CL->normalise)?((score*CL->normalise)/CL->max_ext_value):score;
2194
2195           for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
2196             {
2197               t_s=CL->residue_index[s1][r1][a+SEQ2];
2198               t_r=CL->residue_index[s1][r1][a+R2];
2199               t_w=CL->residue_index[s1][r1][a+WE];
2200               if ( CL->seq_for_quadruplet[t_s])
2201                 {
2202                   for ( b=1; b<CL->residue_index[t_s][t_r][0]; b+=ICHUNK)
2203                     {
2204                       q_s=CL->residue_index[t_s][t_r][b+SEQ2];
2205                       q_r=CL->residue_index[t_s][t_r][b+R2];
2206                       hasch[q_s][q_r]=0;
2207                     }
2208                 }
2209             }
2210
2211           return (int)(score*SCORE_K);
2212         }
2213
2214
2215 Constraint_list * R_extension ( Constraint_list *CL, Constraint_list *R);
2216 int residue_pair_extended_list4rna4 (Constraint_list *CL,int s1, int r1, int s2, int r2 )
2217 {
2218   static int rna_lib;
2219
2220   if (!rna_lib)
2221     {
2222       sprintf ( CL->rna_lib, "%s", seq2rna_lib (CL->S, NULL));
2223       rna_lib=1;
2224     }
2225   return residue_pair_extended_list4rna2 (CL, s1, r1, s2,r2);
2226 }
2227 int residue_pair_extended_list4rna1 (Constraint_list *CL, int s1, int r1, int s2, int r2 )
2228 {
2229   static Constraint_list *R;
2230   if (!R)R=read_rna_lib (CL->S, CL->rna_lib);
2231   return residue_pair_extended_list4rna (CL, R, s1, r1, s2, r2);
2232 }
2233
2234 int residue_pair_extended_list4rna3 (Constraint_list *CL,int s1, int r1, int s2, int r2 )
2235 {
2236   static Constraint_list *R;
2237   if (!R)
2238     {
2239       R=read_rna_lib (CL->S, CL->rna_lib);
2240       rna_lib_extension (CL,R);
2241     }
2242   return residue_pair_extended_list (CL, s1,r1, s2,r2);
2243 }
2244
2245 int residue_pair_extended_list4rna2 (Constraint_list *CL,int s1, int r1, int s2, int r2 )
2246 {
2247   static Constraint_list *R;
2248
2249
2250   if (!R)
2251     {
2252
2253       R=read_rna_lib (CL->S, CL->rna_lib);
2254       rna_lib_extension (CL,R);
2255
2256     }
2257
2258   return residue_pair_extended_list4rna (CL, R, s1, r1, s2, r2);
2259 }
2260 int residue_pair_extended_list4rna ( Constraint_list *CL,Constraint_list *R, int s1, int r1, int s2, int r2 )
2261 {
2262
2263   int a, b, n1, n2;
2264   int list1[100];
2265   int list2[100];
2266   int score=0, score2;
2267
2268
2269
2270   if ( r1<0 || r2<0)return 0;
2271   n1=n2=0;
2272
2273   list1[n1++]=r1;
2274   for (a=1; a<R->residue_index[s1][r1][0]; a+=ICHUNK)
2275     {
2276       list1[n1++]=R->residue_index[s1][r1][a+R2];
2277     }
2278
2279
2280   list2[n2++]=r2;
2281   for (a=1; a<R->residue_index[s2][r2][0]; a+=ICHUNK)
2282     {
2283       list2[n2++]=R->residue_index[s2][r2][a+R2];
2284     }
2285
2286
2287   score=residue_pair_extended_list ( CL, s1,list1[0], s2,list2[0]);
2288
2289   for (score2=0,a=1; a<n1; a++)
2290     for ( b=1; b<n2; b++)
2291       score2=MAX(residue_pair_extended_list ( CL, s1,list1[a], s2,list2[b]), score2);
2292
2293   if ( n1==1 || n2==1);
2294   else score=MAX(score,score2);
2295
2296   return score;
2297 }
2298
2299
2300 int residue_pair_extended_list4rna_ref ( Constraint_list *CL, int s1, int r1, int s2, int r2 )
2301 {
2302   static Constraint_list *R;
2303   int a, b, n1, n2;
2304   int list1[100];
2305   int list2[100];
2306   int score=0, score2;
2307
2308   if (R==NULL)
2309     {
2310       char **list;
2311       int n,a;
2312       list=read_lib_list ( CL->rna_lib, &n);
2313
2314       R=declare_constraint_list ( CL->S,NULL, NULL, 0,NULL, NULL);
2315
2316       for (a=0; a< n; a++)
2317         {
2318           R=read_constraint_list_file (R, list[a]);
2319         }
2320
2321     }
2322
2323   if ( r1<0 || r2<0)return 0;
2324   n1=n2=0;
2325
2326   list1[n1++]=r1;
2327   for (a=1; a<R->residue_index[s1][r1][0]; a+=ICHUNK)
2328     {
2329       list1[n1++]=R->residue_index[s1][r1][a+R2];
2330     }
2331
2332
2333   list2[n2++]=r2;
2334   for (a=1; a<R->residue_index[s2][r2][0]; a+=ICHUNK)
2335     {
2336       list2[n2++]=R->residue_index[s2][r2][a+R2];
2337     }
2338
2339
2340   score=residue_pair_extended_list ( CL, s1,list1[0], s2,list2[0]);
2341
2342   for (score2=0,a=1; a<n1; a++)
2343     for ( b=1; b<n2; b++)
2344       score2=MAX(residue_pair_extended_list ( CL, s1,list1[a], s2,list2[b]), score2);
2345
2346   if ( n1==1 || n2==1);
2347   else score=MAX(score,score2);
2348
2349   return score;
2350 }
2351
2352 static int ** clean_residue_pair_hasch (int s1, int r1, int s2, int r2,int **hasch, Constraint_list *CL);
2353 int residue_pair_extended_list_raw ( Constraint_list *CL, int s1, int r1, int s2, int r2 )
2354         {
2355         double score=0;
2356
2357
2358
2359         int a, t_s, t_r;
2360         static int **hasch;
2361         static int max_len;
2362         int field;
2363         double delta;
2364         /*
2365
2366           function documentation: start
2367
2368           int residue_pair_extended_list ( Constraint_list *CL, int s1, int r1, int s2, int r2);
2369
2370           Computes the extended score for aligning residue seq1(r1) Vs seq2(r2)
2371           Computes: matrix_score
2372                     non extended score
2373                     extended score
2374
2375           The extended score depends on the function index_res_constraint_list.
2376           This function can compare a sequence with itself.
2377
2378           Associated functions: See util constraint list, list extention functions.
2379
2380           function documentation: end
2381         */
2382
2383         field=CL->weight_field;
2384         field=WE;
2385
2386         if ( r1<=0 || r2<=0)return 0;
2387         if ( !hasch || max_len!=(CL->S)->max_len)
2388                {
2389                max_len=(CL->S)->max_len;
2390                if ( hasch) free_int ( hasch, -1);
2391                hasch=declare_int ( (CL->S)->nseq, (CL->S)->max_len+1);
2392                }
2393
2394
2395
2396         /* Check matches for R1 in the indexed lib*/
2397         hasch[s1][r1]=FORBIDEN;
2398         for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
2399           {
2400             t_s=CL->residue_index[s1][r1][a+SEQ2];
2401             t_r=CL->residue_index[s1][r1][a+R2];
2402             hasch[t_s][t_r]=CL->residue_index[s1][r1][a+field];
2403           }
2404
2405         /*Check Matches for r1 <-> r2 in the indexed lib */
2406         for (a=1; a< CL->residue_index[s2][r2][0]; a+=ICHUNK)
2407           {
2408             t_s=CL->residue_index[s2][r2][a+SEQ2];
2409             t_r=CL->residue_index[s2][r2][a+R2];
2410
2411
2412             if (hasch[t_s][t_r])
2413               {
2414                 if (hasch[t_s][t_r]==FORBIDEN)
2415                   {
2416                     score+=CL->residue_index[s2][r2][a+field];
2417                   }
2418                 else
2419                   {
2420                     delta=MIN(hasch[t_s][t_r],CL->residue_index[s2][r2][a+field]);
2421                     score+=delta;
2422                   }
2423               }
2424           }
2425
2426         clean_residue_pair_hasch ( s1, r1,s2, r2, hasch, CL);
2427         return score;
2428         }
2429 int residue_pair_extended_list_4gp ( Constraint_list *CL, int s1, int r1, int s2, int r2 )
2430         {
2431           double score=0;
2432
2433
2434
2435         int a, t_s, t_r;
2436         static int **hasch;
2437         static int max_len;
2438         int field;
2439         double delta;
2440
2441         /*
2442
2443           function documentation: start
2444
2445           int residue_pair_extended_list ( Constraint_list *CL, int s1, int r1, int s2, int r2);
2446
2447           Computes the extended score for aligning residue seq1(r1) Vs seq2(r2)
2448           Computes: matrix_score
2449                     non extended score
2450                     extended score
2451
2452           The extended score depends on the function index_res_constraint_list.
2453           This function can compare a sequence with itself.
2454
2455           Associated functions: See util constraint list, list extention functions.
2456
2457           function documentation: end
2458         */
2459
2460         field=CL->weight_field;
2461         field=WE;
2462
2463         if ( r1<=0 || r2<=0)return 0;
2464         if ( !hasch || max_len!=(CL->S)->max_len)
2465                {
2466                max_len=(CL->S)->max_len;
2467                if ( hasch) free_int ( hasch, -1);
2468                hasch=declare_int ( (CL->S)->nseq, (CL->S)->max_len+1);
2469                }
2470
2471
2472
2473         /* Check matches for R1 in the indexed lib*/
2474         hasch[s1][r1]=FORBIDEN;
2475         for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
2476           {
2477             t_s=CL->residue_index[s1][r1][a+SEQ2];
2478             t_r=CL->residue_index[s1][r1][a+R2];
2479             hasch[t_s][t_r]=CL->residue_index[s1][r1][a+field];
2480           }
2481
2482         /*Check Matches for r1 <-> r2 in the indexed lib */
2483         for (a=1; a< CL->residue_index[s2][r2][0]; a+=ICHUNK)
2484           {
2485             t_s=CL->residue_index[s2][r2][a+SEQ2];
2486             t_r=CL->residue_index[s2][r2][a+R2];
2487
2488
2489             if (hasch[t_s][t_r])
2490               {
2491                 if (hasch[t_s][t_r]==FORBIDEN)
2492                   {
2493                     score+=((float)CL->residue_index[s2][r2][a+2]/NORM_F);
2494                   }
2495                 else
2496                   {
2497                     delta=MAX((((float)hasch[t_s][t_r]/NORM_F)),(((float)CL->residue_index[s2][r2][a+field]/NORM_F)));
2498                     score+=delta;
2499                   }
2500               }
2501           }
2502
2503         clean_residue_pair_hasch ( s1, r1,s2, r2, hasch, CL);
2504         score/=(CL->S)->nseq;
2505         score *=NORM_F;
2506
2507         return score;
2508         }
2509
2510 int residue_pair_extended_list_pc ( Constraint_list *CL, int s1, int r1, int s2, int r2 )
2511         {
2512           //old norm does not take into account the number of effective intermediate sequences
2513
2514           double score=0;
2515           int a, t_s, t_r;
2516           static int **hasch;
2517           int nclean;
2518           static int max_len;
2519           int field;
2520           double delta;
2521           float norm1=0;
2522           float norm2=0;
2523
2524           /*
2525
2526           function documentation: start
2527
2528           int residue_pair_extended_list ( Constraint_list *CL, int s1, int r1, int s2, int r2);
2529
2530           Computes the extended score for aligning residue seq1(r1) Vs seq2(r2)
2531           Computes: matrix_score
2532                     non extended score
2533                     extended score
2534
2535           The extended score depends on the function index_res_constraint_list.
2536           This function can compare a sequence with itself.
2537
2538           Associated functions: See util constraint list, list extention functions.
2539
2540           function documentation: end
2541         */
2542
2543         field=CL->weight_field;
2544         field=WE;
2545
2546         if ( r1<=0 || r2<=0)return 0;
2547         if ( !hasch || max_len!=(CL->S)->max_len)
2548                {
2549                max_len=(CL->S)->max_len;
2550                if ( hasch) free_int ( hasch, -1);
2551                hasch=declare_int ( (CL->S)->nseq, (CL->S)->max_len+1);
2552
2553                }
2554
2555
2556         /* Check matches for R1 in the indexed lib*/
2557         hasch[s1][r1]=FORBIDEN;
2558         for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
2559           {
2560             t_s=CL->residue_index[s1][r1][a+SEQ2];
2561             t_r=CL->residue_index[s1][r1][a+R2];
2562             hasch[t_s][t_r]=CL->residue_index[s1][r1][a+field];
2563             norm1++;
2564
2565           }
2566
2567
2568         /*Check Matches for r1 <-> r2 in the indexed lib */
2569         for (a=1; a< CL->residue_index[s2][r2][0]; a+=ICHUNK)
2570           {
2571             t_s=CL->residue_index[s2][r2][a+SEQ2];
2572             t_r=CL->residue_index[s2][r2][a+R2];
2573             norm2++;
2574
2575             if (hasch[t_s][t_r])
2576               {
2577                 if (hasch[t_s][t_r]==FORBIDEN)
2578                   {
2579                     score+=((float)CL->residue_index[s2][r2][a+field]/NORM_F);
2580                   }
2581                 else
2582                   {
2583                     delta=MIN((((float)hasch[t_s][t_r]/NORM_F)),(((float)CL->residue_index[s2][r2][a+field]/NORM_F)));
2584                     score+=delta;
2585                   }
2586               }
2587           }
2588
2589         clean_residue_pair_hasch ( s1, r1,s2, r2, hasch, CL);
2590
2591         //New Normalization
2592         norm1=MIN(norm1,norm2);
2593
2594         //Old Normailzation: on the number of sequences, useless when not doiang an all against all
2595         //norm1=(CL->S)->nseq;
2596
2597         score=(norm1)?score/norm1:0;
2598
2599
2600         return score*NORM_F;
2601         }
2602 int residue_pair_extended_list_pc_old ( Constraint_list *CL, int s1, int r1, int s2, int r2 )
2603         {
2604           //old norm does not take into account the number of effective intermediate sequences
2605
2606           double score=0;
2607           int a, t_s, t_r;
2608           static int **hasch;
2609
2610           static int max_len;
2611           int field;
2612           double delta;
2613
2614           /*
2615
2616           function documentation: start
2617
2618           int residue_pair_extended_list ( Constraint_list *CL, int s1, int r1, int s2, int r2);
2619
2620           Computes the extended score for aligning residue seq1(r1) Vs seq2(r2)
2621           Computes: matrix_score
2622                     non extended score
2623                     extended score
2624
2625           The extended score depends on the function index_res_constraint_list.
2626           This function can compare a sequence with itself.
2627
2628           Associated functions: See util constraint list, list extention functions.
2629
2630           function documentation: end
2631         */
2632
2633         field=CL->weight_field;
2634         field=WE;
2635
2636         if ( r1<=0 || r2<=0)return 0;
2637         if ( !hasch || max_len!=(CL->S)->max_len)
2638                {
2639                max_len=(CL->S)->max_len;
2640                if ( hasch) free_int ( hasch, -1);
2641                hasch=declare_int ( (CL->S)->nseq, (CL->S)->max_len+1);
2642                }
2643
2644
2645         /* Check matches for R1 in the indexed lib*/
2646         hasch[s1][r1]=FORBIDEN;
2647         for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
2648           {
2649             t_s=CL->residue_index[s1][r1][a+SEQ2];
2650             t_r=CL->residue_index[s1][r1][a+R2];
2651             hasch[t_s][t_r]=CL->residue_index[s1][r1][a+field];
2652           }
2653
2654
2655         /*Check Matches for r1 <-> r2 in the indexed lib */
2656         for (a=1; a< CL->residue_index[s2][r2][0]; a+=ICHUNK)
2657           {
2658             t_s=CL->residue_index[s2][r2][a+SEQ2];
2659             t_r=CL->residue_index[s2][r2][a+R2];
2660
2661
2662             if (hasch[t_s][t_r])
2663               {
2664                 if (hasch[t_s][t_r]==FORBIDEN)
2665                   {
2666                     score+=((float)CL->residue_index[s2][r2][a+2]/NORM_F);
2667                   }
2668                 else
2669                   {
2670                     delta=MIN((((float)hasch[t_s][t_r]/NORM_F)),(((float)CL->residue_index[s2][r2][a+field]/NORM_F)));
2671                     score+=delta;
2672                   }
2673               }
2674           }
2675
2676         clean_residue_pair_hasch ( s1, r1,s2, r2, hasch, CL);
2677         score/=(CL->S)->nseq;
2678         return score*NORM_F;
2679         }
2680
2681
2682 int residue_pair_extended_list ( Constraint_list *CL, int s1, int r1, int s2, int r2 )
2683         {
2684         double score=0;
2685         double max_score=0;
2686         double max_val=0;
2687
2688         int a, t_s, t_r;
2689         static int **hasch;
2690         static int max_len;
2691         int field;
2692
2693         /*
2694           new function: self normalized
2695           function documentation: start
2696
2697           int residue_pair_extended_list ( Constraint_list *CL, int s1, int r1, int s2, int r2);
2698
2699           Computes the extended score for aligning residue seq1(r1) Vs seq2(r2)
2700           Computes: matrix_score
2701                     non extended score
2702                     extended score
2703
2704           The extended score depends on the function index_res_constraint_list.
2705           This function can compare a sequence with itself.
2706
2707           Associated functions: See util constraint list, list extention functions.
2708
2709           function documentation: end
2710         */
2711
2712         field=CL->weight_field;
2713         field=WE;
2714
2715         if ( r1<=0 || r2<=0)return 0;
2716         if ( !hasch || max_len!=(CL->S)->max_len)
2717                {
2718                max_len=(CL->S)->max_len;
2719                if ( hasch) free_int ( hasch, -1);
2720                hasch=declare_int ( (CL->S)->nseq, (CL->S)->max_len+1);
2721                }
2722
2723
2724
2725         /* Check matches for R1 in the indexed lib*/
2726         hasch[s1][r1]=FORBIDEN;
2727         for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
2728           {
2729             t_s=CL->residue_index[s1][r1][a+SEQ2];
2730             t_r=CL->residue_index[s1][r1][a+R2];
2731             hasch[t_s][t_r]=CL->residue_index[s1][r1][a+field];
2732             max_score+=CL->residue_index[s1][r1][a+field];
2733           }
2734
2735         /*Check Matches for r1 <-> r2 in the indexed lib */
2736         for (a=1; a< CL->residue_index[s2][r2][0]; a+=ICHUNK)
2737           {
2738             t_s=CL->residue_index[s2][r2][a+SEQ2];
2739             t_r=CL->residue_index[s2][r2][a+R2];
2740
2741
2742             if (hasch[t_s][t_r])
2743               {
2744                 if (hasch[t_s][t_r]==FORBIDEN)
2745                   {
2746                     score+=CL->residue_index[s2][r2][a+field];
2747                     max_score+=CL->residue_index[s2][r2][a+field];
2748                   }
2749                 else
2750                   {
2751                     double delta;
2752                     delta=MIN(hasch[t_s][t_r],CL->residue_index[s2][r2][a+field]);
2753
2754                     score+=delta;
2755                     max_score-=hasch[t_s][t_r];
2756                     max_score+=delta;
2757                     max_val=MAX(max_val,delta);
2758                   }
2759               }
2760             else
2761               {
2762                 max_score+=CL->residue_index[s2][r2][a+field];
2763               }
2764           }
2765
2766         max_score-=hasch[s2][r2];
2767         clean_residue_pair_hasch ( s1, r1,s2, r2, hasch, CL);
2768
2769
2770         if ( max_score==0)score=0;
2771         else if ( CL->normalise)
2772           {
2773             score=((score*CL->normalise)/max_score)*SCORE_K;
2774             if (max_val> CL->normalise)
2775               {
2776                 score*=max_val/(double)CL->normalise;
2777               }
2778           }
2779         return (int) score;
2780         }
2781 int ** clean_residue_pair_hasch (int s1, int r1, int s2, int r2,int **hasch, Constraint_list *CL)
2782   {
2783     int a, t_s, t_r;
2784     if ( !hasch) return hasch;
2785
2786     for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
2787       {
2788         t_s=CL->residue_index[s1][r1][a+SEQ2];
2789         t_r=CL->residue_index[s1][r1][a+R2];
2790         hasch[t_s][t_r]=0;
2791       }
2792     hasch[s1][r1]=hasch[s2][r2]=0;
2793     return hasch;
2794   }
2795
2796 int residue_pair_test_function ( Constraint_list *CL, int s1, int r1, int s2, int r2 )
2797         {
2798           double score=0;
2799
2800           int a, t_s, t_r;
2801           static int **hasch;
2802           static int max_len;
2803           int cons1;
2804           int cons2;
2805
2806           int field;
2807         /*
2808           function documentation: start
2809
2810           int residue_pair_extended_list ( Constraint_list *CL, int s1, int r1, int s2, int r2);
2811
2812           Computes the extended score for aligning residue seq1(r1) Vs seq2(r2)
2813           Computes: matrix_score
2814                     non extended score
2815                     extended score
2816
2817           The extended score depends on the function index_res_constraint_list.
2818           This function can compare a sequence with itself.
2819
2820           Associated functions: See util constraint list, list extention functions.
2821
2822           function documentation: end
2823         */
2824
2825
2826         CL->weight_field=WE;
2827         field=CL->weight_field;
2828
2829
2830         if ( r1<=0 || r2<=0)return 0;
2831         if ( !hasch || max_len!=(CL->S)->max_len)
2832                {
2833                max_len=(CL->S)->max_len;
2834                if ( hasch) free_int ( hasch, -1);
2835                hasch=declare_int ( (CL->S)->nseq, (CL->S)->max_len+1);
2836                }
2837
2838
2839
2840         hasch[s1][r1]=1000;
2841         for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
2842           {
2843             t_s=CL->residue_index[s1][r1][a+SEQ2];
2844             t_r=CL->residue_index[s1][r1][a+R2];
2845             hasch[t_s][t_r]=CL->residue_index[s1][r1][a+field];
2846           }
2847
2848
2849         for (a=1; a< CL->residue_index[s2][r2][0]; a+=ICHUNK)
2850           {
2851             t_s=CL->residue_index[s2][r2][a+SEQ2];
2852             t_r=CL->residue_index[s2][r2][a+R2];
2853             if (hasch[t_s][t_r])
2854               {
2855                 cons1=hasch[t_s][t_r];
2856                 cons2=CL->residue_index[s2][r2][a+field];
2857                 score +=MIN(cons1,cons2);
2858               }
2859           }
2860
2861
2862         score=(CL->normalise)?((score*CL->normalise)/CL->max_ext_value):score;
2863
2864         for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
2865           {
2866             t_s=CL->residue_index[s1][r1][a+SEQ2];
2867             t_r=CL->residue_index[s1][r1][a+R2];
2868             hasch[t_s][t_r]=0;
2869           }
2870         hasch[s1][r1]=hasch[s2][r2]=0;
2871
2872
2873         return (int)(score*SCORE_K);
2874         }
2875
2876 int residue_pair_relative_extended_list ( Constraint_list *CL, int s1, int r1, int s2, int r2 )
2877         {
2878         int a, t_s, t_r;
2879         static int **hasch;
2880         static int max_len;
2881         int score=0;
2882         int total_score=0;
2883         int field;
2884         /*
2885           function documentation: start
2886
2887           int residue_pair_extended_list ( Constraint_list *CL, int s1, int r1, int s2, int r2);
2888
2889           Computes the extended score for aligning residue seq1(r1) Vs seq2(r2)
2890           Computes: matrix_score
2891                     non extended score
2892                     extended score
2893
2894           The extended score depends on the function index_res_constraint_list.
2895           This function can compare a sequence with itself.
2896
2897           Associated functions: See util constraint list, list extention functions.
2898
2899           function documentation: end
2900         */
2901
2902
2903
2904         field=CL->weight_field;
2905         field=WE;
2906
2907         if ( r1<=0 || r2<=0)return 0;
2908         if ( !hasch || max_len!=(CL->S)->max_len)
2909                {
2910                max_len=(CL->S)->max_len;
2911                if ( hasch) free_int ( hasch, -1);
2912                hasch=declare_int ( (CL->S)->nseq, (CL->S)->max_len+1);
2913                }
2914
2915
2916
2917         hasch[s1][r1]=100000;
2918         for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
2919           {
2920             t_s=CL->residue_index[s1][r1][a+SEQ2];
2921             t_r=CL->residue_index[s1][r1][a+R2];
2922             hasch[t_s][t_r]=CL->residue_index[s1][r1][a+field];
2923             total_score+=CL->residue_index[s1][r1][a+field];
2924           }
2925
2926
2927         for (a=1; a< CL->residue_index[s2][r2][0]; a+=ICHUNK)
2928           {
2929             t_s=CL->residue_index[s2][r2][a+SEQ2];
2930             t_r=CL->residue_index[s2][r2][a+R2];
2931             total_score+=CL->residue_index[s1][r1][a+field];
2932             if (hasch[t_s][t_r])
2933               {
2934                 if (field==WE){score+=2*MIN(hasch[t_s][t_r],CL->residue_index[s2][r2][a+field]);}
2935               }
2936           }
2937
2938         score=((CL->normalise*score)/total_score);
2939
2940
2941         for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
2942           {
2943             t_s=CL->residue_index[s1][r1][a+SEQ2];
2944             t_r=CL->residue_index[s1][r1][a+R2];
2945             hasch[t_s][t_r]=0;
2946           }
2947         hasch[s1][r1]=hasch[s2][r2]=0;
2948
2949         return score*SCORE_K;
2950         }
2951 int residue_pair_extended_list_g_coffee_quadruplet ( Constraint_list *CL, int s1, int r1, int s2, int r2 )
2952 {
2953    int t_s, t_r, t_w, q_s, q_r, q_w;
2954           int a, b;
2955           static int **hasch;
2956           int score=0, s=0;
2957
2958           int field;
2959           /* This measure the quadruplets cost on a pair of residues*/
2960
2961           field=CL->weight_field;
2962
2963           if ( r1<=0 || r2<=0)return 0;
2964           if ( !hasch)
2965             {
2966               hasch=vcalloc ( (CL->S)->nseq, sizeof (int*));
2967               for ( a=0; a< (CL->S)->nseq; a++)hasch[a]=vcalloc ( (CL->S)->len[a]+1, sizeof (int));
2968             }
2969
2970
2971           hasch[s1][r1]=100000;
2972           for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
2973             {
2974               t_s=CL->residue_index[s1][r1][a+SEQ2];
2975               t_r=CL->residue_index[s1][r1][a+R2];
2976               t_w=CL->residue_index[s1][r1][a+field];
2977               if ( CL->seq_for_quadruplet[t_s])
2978                 {
2979                   for ( b=1; b<CL->residue_index[t_s][t_r][0]; b+=ICHUNK)
2980                     {
2981                       q_s=CL->residue_index[t_s][t_r][b+SEQ2];
2982                       q_r=CL->residue_index[t_s][t_r][b+R2];
2983                       if (CL-> seq_for_quadruplet[q_s])
2984                         {
2985
2986                           hasch[q_s][q_r]=MIN(CL->residue_index[t_s][t_r][b+2],t_w);
2987                         }
2988                     }
2989                 }
2990             }
2991
2992
2993           for (s=0,score=0,a=1; a< CL->residue_index[s2][r2][0]; a+=ICHUNK)
2994             {
2995               t_s=CL->residue_index[s2][r2][a+SEQ2];
2996               t_r=CL->residue_index[s2][r2][a+R2];
2997               t_w=CL->residue_index[s2][r2][a+field];
2998               if ( CL->seq_for_quadruplet[t_s])
2999                 {
3000                   for ( b=1; b<CL->residue_index[t_s][t_r][0]; b+=ICHUNK)
3001                     {
3002                       q_s=CL->residue_index[t_s][t_r][b+SEQ2];
3003                       q_r=CL->residue_index[t_s][t_r][b+R2];
3004                       q_w=CL->residue_index[t_s][t_r][b+field];
3005                       if (hasch[q_s][q_r] && CL->seq_for_quadruplet[q_s])
3006                         s=MIN(hasch[q_s][q_r],MIN(CL->residue_index[t_s][t_r][b+2],q_w));
3007                       score=MAX(score, s);
3008                     }
3009                 }
3010             }
3011
3012           for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
3013             {
3014               t_s=CL->residue_index[s1][r1][a+SEQ2];
3015               t_r=CL->residue_index[s1][r1][a+R2];
3016               t_w=CL->residue_index[s1][r1][a+field];
3017               if ( CL->seq_for_quadruplet[t_s])
3018                 {
3019                   for ( b=1; b<CL->residue_index[t_s][t_r][0]; b+=ICHUNK)
3020                     {
3021                       q_s=CL->residue_index[t_s][t_r][b+SEQ2];
3022                       q_r=CL->residue_index[t_s][t_r][b+R2];
3023                       hasch[q_s][q_r]=0;
3024                     }
3025                 }
3026             }
3027
3028           return score*SCORE_K;
3029         }
3030 int residue_pair_extended_list_g_coffee ( Constraint_list *CL, int s1, int r1, int s2, int r2 )
3031         {
3032         int a, t_s, t_r;
3033         static int **hasch;
3034         int score=0,s;
3035
3036         int field;
3037         /*
3038           function documentation: start
3039
3040           int residue_pair_extended_list ( Constraint_list *CL, int s1, int r1, int s2, int r2);
3041
3042           Computes the extended score for aligning residue seq1(r1) Vs seq2(r2)
3043           Computes: matrix_score
3044                     non extended score
3045                     extended score
3046
3047           The extended score depends on the function index_res_constraint_list.
3048           This function can compare a sequence with itself.
3049
3050           Associated functions: See util constraint list, list extention functions.
3051
3052           function documentation: end
3053         */
3054
3055         field=CL->weight_field;
3056
3057         if ( r1<=0 || r2<=0)return 0;
3058         if ( !hasch)
3059                {
3060                hasch=vcalloc ( (CL->S)->nseq, sizeof (int*));
3061                for ( a=0; a< (CL->S)->nseq; a++)hasch[a]=vcalloc ( (CL->S)->len[a]+1, sizeof (int));
3062                }
3063
3064
3065
3066         hasch[s1][r1]=100000;
3067         for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
3068           {
3069             t_s=CL->residue_index[s1][r1][a+SEQ2];
3070             t_r=CL->residue_index[s1][r1][a+R2];
3071             hasch[t_s][t_r]=CL->residue_index[s1][r1][a+field];
3072           }
3073
3074
3075         for (s=0, score=0,a=1; a< CL->residue_index[s2][r2][0]; a+=ICHUNK)
3076           {
3077             t_s=CL->residue_index[s2][r2][a+SEQ2];
3078             t_r=CL->residue_index[s2][r2][a+R2];
3079
3080             if (hasch[t_s][t_r])
3081               {
3082                 if (field==WE)
3083                   {s=MIN(hasch[t_s][t_r],CL->residue_index[s2][r2][a+WE]);
3084                    score=MAX(s,score);
3085                   }
3086               }
3087           }
3088
3089         for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
3090           {
3091             t_s=CL->residue_index[s1][r1][a+SEQ2];
3092             t_r=CL->residue_index[s1][r1][a+R2];
3093             hasch[t_s][t_r]=0;
3094           }
3095         hasch[s1][r1]=hasch[s2][r2]=0;
3096
3097         return score*SCORE_K;
3098         }
3099
3100 int extend_residue_pair ( Constraint_list *CL, int s1, int r1, int s2, int r2)
3101         {
3102         double score=0;
3103
3104         int a, t_s, t_r, p;
3105         static int **hasch;
3106
3107         static int *entry;
3108         int *r;
3109         int field;
3110
3111
3112
3113         /*
3114           This is the generic Function->works with everything
3115           should be gradually phased out
3116
3117
3118           int extend_residue_pair ( Constraint_list *CL, int s1, int r1, int s2, int r2, int field )
3119
3120           Computes the extended score for aligning residue seq1(r1) Vs seq2(r2)
3121           Computes: matrix_score
3122                     non extended score
3123                     extended score
3124
3125           The extended score depends on the function index_res_constraint_list.
3126           This function can compare a sequence with itself.
3127
3128           Associated functions: See util constraint list, list extention functions.
3129           function documentation: end
3130         */
3131
3132
3133         field=CL->weight_field;
3134
3135         if ( r1<=0 || r2<=0)return 0;
3136         else if ( !CL->residue_index && CL->M)
3137            {
3138             return evaluate_matrix_score (CL, s1,r1, s2, r2);
3139            }
3140
3141         else if ( !CL->extend_jit)
3142            {
3143             if ( !entry) entry=vcalloc (LIST_N_FIELDS , sizeof (int));
3144             entry[SEQ1]=s1;
3145             entry[SEQ2]=s2;
3146             entry[R1]=r1;
3147             entry[R2]=r2;
3148             r=main_search_in_list_constraint( entry,&p,4,CL);
3149             if (r==NULL)return 0;
3150             else return r[field];
3151            }
3152         else
3153            {
3154            if ( !hasch)
3155                {
3156                hasch=vcalloc ( (CL->S)->nseq, sizeof (int*));
3157                for ( a=0; a< (CL->S)->nseq; a++)hasch[a]=vcalloc ( (CL->S)->len[a]+1, sizeof (int));
3158                }
3159
3160
3161
3162            hasch[s1][r1]=100000;
3163            for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
3164                 {
3165                 t_s=CL->residue_index[s1][r1][a+SEQ2];
3166                 t_r=CL->residue_index[s1][r1][a+R2];
3167                 hasch[t_s][t_r]=CL->residue_index[s1][r1][a+WE];
3168                 }
3169
3170
3171            for (a=1; a< CL->residue_index[s2][r2][0]; a+=ICHUNK)
3172                {
3173                t_s=CL->residue_index[s2][r2][a+SEQ2];
3174                t_r=CL->residue_index[s2][r2][a+R2];
3175                if (hasch[t_s][t_r])
3176                    {
3177                    if (field==WE)score+=MIN(hasch[t_s][t_r],CL->residue_index[s2][r2][a+WE] );
3178
3179                    }
3180                }
3181            score=(CL->normalise)?((score*CL->normalise)/CL->max_ext_value):score;
3182            for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
3183                {
3184                t_s=CL->residue_index[s1][r1][a+SEQ2];
3185                t_r=CL->residue_index[s1][r1][a+R2];
3186                hasch[t_s][t_r]=0;
3187                }
3188            hasch[s1][r1]=hasch[s2][r2]=0;
3189
3190            return (int)(score*SCORE_K);
3191            }
3192         }
3193 /*********************************************************************************************/
3194 /*                                                                                           */
3195 /*         FUNCTIONS FOR GETTING THE PW COST :  CL->get_dp_cost                              */
3196 /*                                                                                           */
3197 /*********************************************************************************************/
3198 int get_dp_cost_blosum_matrix (Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
3199 {
3200   int s1, r1, s2, r2;
3201   static int **matrix;
3202
3203   if (!matrix) matrix=read_matrice ("blosum62mt");
3204   s1=A->order[list1[0]][0];
3205   s2=A->order[list2[0]][0];
3206   r1=pos1[list1[0]][col1];
3207   r2=pos2[list2[0]][col2];
3208
3209   /*dp cost function: works only with two sequences*/
3210
3211   if ( seq2R_template_profile (CL->S,s1) ||seq2R_template_profile (CL->S,s2))
3212     return evaluate_aln_profile_score ( CL, s1,r1, s2, r2) -CL->nomatch*SCORE_K;
3213   else if (r1>0 && r2>0)
3214             {
3215             r1--;
3216             r2--;
3217
3218
3219             return matrix [(CL->S)->seq[s1][r1]-'A'][(CL->S)->seq[s2][r2]-'A']*SCORE_K -CL->nomatch*SCORE_K;
3220
3221             }
3222   else
3223     return -CL->nomatch*SCORE_K ;
3224 }
3225 int get_dp_cost_pam_matrix (Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
3226 {
3227   int s1, r1, s2, r2;
3228   static int **matrix;
3229
3230   if (!matrix) matrix=read_matrice ("pam250mt");
3231   s1=A->order[list1[0]][0];
3232   s2=A->order[list2[0]][0];
3233   r1=pos1[list1[0]][col1];
3234   r2=pos2[list2[0]][col2];
3235
3236   /*dp cost function: works only with two sequences*/
3237
3238
3239   if ( seq2R_template_profile (CL->S,s1) ||seq2R_template_profile (CL->S,s2))
3240     return evaluate_aln_profile_score ( CL, s1,r1, s2, r2) -CL->nomatch*SCORE_K;
3241   else if (r1>0 && r2>0)
3242             {
3243             r1--;
3244             r2--;
3245
3246
3247             return matrix [(CL->S)->seq[s1][r1]-'A'][(CL->S)->seq[s2][r2]-'A']*SCORE_K -CL->nomatch*SCORE_K;
3248
3249             }
3250   else
3251     return -CL->nomatch*SCORE_K ;
3252 }
3253
3254 int get_dp_cost_pw_matrix (Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
3255 {
3256   int s1, r1, s2, r2;
3257
3258   s1=A->order[list1[0]][0];
3259   s2=A->order[list2[0]][0];
3260   r1=pos1[list1[0]][col1];
3261   r2=pos2[list2[0]][col2];
3262
3263   /*dp cost function: works only with two sequences*/
3264   if ( seq2R_template_profile (CL->S,s1) ||seq2R_template_profile (CL->S,s2))
3265     return evaluate_aln_profile_score ( CL, s1,r1, s2, r2) -CL->nomatch*SCORE_K;
3266   else if (r1>0 && r2>0)
3267             {
3268             r1--;
3269             r2--;
3270
3271
3272             return CL->M[(CL->S)->seq[s1][r1]-'A'][(CL->S)->seq[s2][r2]-'A']*SCORE_K -CL->nomatch*SCORE_K;
3273
3274             }
3275   else
3276     return -CL->nomatch*SCORE_K ;
3277 }
3278
3279 /*********************************************************************************************/
3280 /*                                                                                           */
3281 /*         FUNCTIONS FOR GETTING THE COST :  CL->get_dp_cost                                     */
3282 /*                                                                                           */
3283 /*********************************************************************************************/
3284
3285
3286
3287 int get_cdna_best_frame_dp_cost (Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
3288     {
3289         int a, b;
3290         int n=4;
3291         int s;
3292         char a1, a2;
3293         static int l1, l2;
3294         static Alignment *B;
3295         static int **score;
3296
3297         if ( !score)score=declare_int(3, 2);
3298
3299         if (!A)
3300            {
3301                free_aln(B);
3302                B=NULL;
3303                return UNDEFINED;
3304            }
3305         if (!B)
3306            {
3307                if (ns1+ns2>2){fprintf ( stderr, "\nERROR: get_cdna_dp_cost mode is only for pair-wise ALN [FATAL]\n");crash("");}
3308                free_aln (B);
3309                B=copy_aln (A, NULL);
3310
3311                l1=(int)strlen ( A->seq_al[list1[0]]);
3312                for ( b=0; b<l1; b++)
3313                        B->seq_al[list1[0]][b]=translate_dna_codon (A->seq_al[list1[0]]+b, 'x');
3314                l2=(int)strlen ( A->seq_al[list2[0]]);
3315                for ( b=0; b<l2; b++)
3316                        B->seq_al[list2[0]][b]=translate_dna_codon (A->seq_al[list2[0]]+b, 'x');
3317            }
3318
3319 /*Set the frame*/
3320
3321         for ( a=0; a< 3; a++)score[a][0]=score[a][1]=0;
3322         for ( a=col1-(n*3),b=col2-(n*3); a<col1+(n*3) ; a++, b++)
3323                 {
3324                     if ( a<0 || b<0 || a>=l1 || b>=l2)continue;
3325
3326                     a1=tolower(B->seq_al[list1[0]][a]);
3327                     a2=tolower(B->seq_al[list2[0]][b]);
3328
3329                     score[a%3][0]+=(a1=='x' || a2=='x')?0:CL->M[a1-'A'][a2-'A'];
3330                     score[a%3][1]++;
3331                  }
3332
3333         for ( a=0; a< 3; a++)score[a][0]=(score[a][1]>0)?(score[a][0]/score[a][1]):0;
3334         if ( score[0][0]>score[1][0] &&  score[0][0]>score[2][0])
3335             s=score[0][0];
3336         else if ( score[1][0]>score[0][0] &&  score[1][0]>score[2][0])
3337             s=score[1][0];
3338         else s=score[2][0];
3339
3340         return s*SCORE_K;
3341
3342         }
3343
3344 int get_dp_cost_quadruplet ( Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
3345         {
3346           int score;
3347
3348
3349           if ( ns1==1 || ns2==1)
3350              score=slow_get_dp_cost ( A, pos1, ns1, list1,col1, pos2, ns2, list2, col2, CL);
3351           else
3352             score=fast_get_dp_cost_quadruplet ( A, pos1, ns1, list1,col1, pos2, ns2, list2, col2, CL);
3353
3354           return (score==UNDEFINED)?UNDEFINED:(score-SCORE_K*CL->nomatch);
3355         }
3356 int get_dp_cost ( Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
3357         {
3358         int MODE=0;
3359         int score;
3360
3361
3362
3363         if (A==NULL)return 0;
3364
3365         if (MODE!=2 || MODE==0  || (!CL->residue_index && CL->M) || (!CL->residue_index && CL->T)|| ns1==1 || ns2==1)
3366           score=slow_get_dp_cost ( A, pos1, ns1, list1,col1, pos2, ns2, list2, col2, CL);
3367         else if (MODE==1 || MODE==2)
3368           score=fast_get_dp_cost ( A, pos1, ns1, list1,col1, pos2, ns2, list2, col2, CL);
3369         else
3370           score=0;
3371
3372
3373
3374         return (score==UNDEFINED)?UNDEFINED:(score-SCORE_K*CL->nomatch);
3375         }
3376 int ***make_cw_lu (int **cons, int l, Constraint_list *CL);
3377 int ***make_cw_lu (int **cons, int l, Constraint_list *CL)
3378 {
3379   int ***lu;
3380   int p, a,r;
3381
3382   lu=declare_arrayN(3, sizeof (int),l,26, 2);
3383   for ( p=0; p<l ; p++)
3384     {
3385       for (r=0; r<26; r++)
3386         {
3387           for (a=3; a<cons[p][1]; a+=3)
3388             {
3389               lu[p][r][0]+=cons[p][a+1]*CL->M[r][cons[p][a]-'A'];
3390               lu[p][r][1]+=cons[p][a+1];
3391             }
3392         }
3393     }
3394   return lu;
3395 }
3396
3397 int cw_profile_get_dp_cost ( Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
3398 {
3399           static int last_tag;
3400           static int *pr, ***lu;
3401           int score;
3402           static int *list[2], ns[2], **cons[2], ref;
3403           int  eva_col,ref_col, a, p, r;
3404           float t1, t2;
3405
3406
3407
3408
3409           if (last_tag!=A->random_tag)
3410             {
3411               int n1, n2;
3412
3413               last_tag=A->random_tag;
3414               list[0]=list1;ns[0]=ns1;
3415               list[1]=list2;ns[1]=ns2;
3416               free_int (cons[0],-1);free_int (cons[1],-1);free_arrayN((void*)lu,3);
3417               cons[0]=NULL;cons[1]=NULL;lu=NULL;
3418
3419               n1=sub_aln2nseq_prf (A, ns[0], list[0]);
3420               n2=sub_aln2nseq_prf (A, ns[1], list[1]);
3421               if ( n1>1 || n2>1)
3422                 {
3423                   cons[0]=sub_aln2count_mat2 (A, ns[0], list[0]);
3424                   cons[1]=sub_aln2count_mat2 (A, ns[1], list[1]);
3425                   ref=(ns[0]>ns[1])?0:1;
3426                   lu=make_cw_lu(cons[ref],(int)strlen(A->seq_al[list[ref][0]]), CL);
3427                 }
3428             }
3429
3430           if (!lu)
3431             {
3432               char r1, r2;
3433               r1=A->seq_al[list1[0]][col1];
3434               r2=A->seq_al[list2[0]][col2];
3435               if ( r1!='-' && r2!='-')
3436                 return CL->M[r1-'A'][r2-'A']*SCORE_K -SCORE_K*CL->nomatch;
3437               else
3438                 return -SCORE_K*CL->nomatch;
3439             }
3440           else
3441             {
3442               eva_col= (ref==0)?col2:col1;
3443               ref_col= (ref==0)?col1:col2;
3444               pr=cons[1-ref][eva_col];
3445               t1=t2=0;
3446               for (a=3; a< pr[1]; a+=3)
3447                 {
3448                   r=tolower(pr[a]);
3449                   p=  pr[a+1];
3450
3451                   t1+=lu[ref_col][r-'a'][0]*p;
3452                   t2+=lu[ref_col][r-'a'][1]*p;
3453                 }
3454               score=(t2==0)?0:(t1*SCORE_K)/t2;
3455               score -=SCORE_K*CL->nomatch;
3456               return score;
3457             }
3458 }
3459 int cw_profile_get_dp_cost_old ( Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
3460 {
3461           static int last_tag;
3462           static int **cons1, **cons2;
3463           int score;
3464
3465
3466           if (last_tag!=A->random_tag)
3467             {
3468               last_tag=A->random_tag;
3469               free_int (cons1,-1);free_int (cons2,-1);
3470               cons1=sub_aln2count_mat2 (A, ns1, list1);
3471               cons2=sub_aln2count_mat2 (A, ns2, list2);
3472             }
3473           score=cw_profile_profile (cons1[col1], cons2[col2], CL)-SCORE_K*CL->nomatch;
3474           return score;
3475 }
3476
3477 int cw_profile_get_dp_cost_window ( Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
3478 {
3479           static int last_tag;
3480           static int **cons1, **cons2;
3481           int a, score, window_size=5, n, p1, p2;
3482
3483
3484         if (last_tag!=A->random_tag)
3485           {
3486             last_tag=A->random_tag;
3487             free_int (cons1,-1);free_int (cons2,-1);
3488             cons1=sub_aln2count_mat2 (A, ns1, list1);
3489             cons2=sub_aln2count_mat2 (A, ns2, list2);
3490           }
3491
3492         for (n=0,score=0,a=0; a<window_size; a++)
3493           {
3494             p1=col1-a;
3495             p2=col2-a;
3496             if ( p1<0 || cons1[p1][0]==END_ARRAY)continue;
3497             if ( p2<0 || cons2[p2][0]==END_ARRAY)continue;
3498             score+=cw_profile_profile (cons1[col1], cons2[col2], CL)-SCORE_K*CL->nomatch;
3499             n++;
3500           }
3501         if (n>0)score/=n;
3502
3503         return score;
3504         }
3505
3506 int consensus_get_dp_cost ( Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
3507         {
3508           static int last_tag;
3509           static char *seq1, *seq2;
3510
3511
3512         /*Works only with matrix*/
3513           if (last_tag !=A->random_tag)
3514             {
3515               int ls1, ls2, lenls1, lenls2;
3516
3517               last_tag=A->random_tag;
3518               vfree (seq1);vfree (seq2);
3519               seq1=sub_aln2cons_seq_mat (A, ns1, list1, "blosum62mt");
3520               seq2=sub_aln2cons_seq_mat (A, ns2, list2, "blosum62mt");
3521               ls1=list1[ns1-1];ls2=list2[ns2-1];
3522               lenls1=(int)strlen (A->seq_al[ls1]); lenls2=(int)strlen (A->seq_al[ls2]);
3523           }
3524
3525         return (CL->M[seq1[col1]-'A'][seq2[col2]-'A']*SCORE_K)-SCORE_K*CL->nomatch;
3526         }
3527
3528 int fast_get_dp_cost_quadruplet ( Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
3529         {
3530         /*WARNING: WORKS ONLY WITH List to Extend*/
3531           /*That function does a quruple extension beween two columns by pooling the residues together*/
3532
3533           double score=0;
3534
3535           int a,b, c;
3536           int n_gap1=0;
3537           int n_gap2=0;
3538
3539           int s1, rs1, r1, t_r, t_s,t_w, q_r, q_s, q_w, s2, rs2, r2;
3540           int **buf_pos, buf_ns, *buf_list, buf_col;
3541
3542           static int **hasch1;
3543           static int **hasch2;
3544
3545           static int **n_hasch1;
3546           static int **n_hasch2;
3547
3548           static int **is_in_col1;
3549           static int **is_in_col2;
3550
3551
3552           if (ns2>ns1)
3553             {
3554               buf_pos=pos1;
3555               buf_ns=ns1;
3556               buf_list=list1;
3557               buf_col=col1;
3558
3559               pos1=pos2;
3560               ns1=ns2;
3561               list1=list2;
3562               col1=col2;
3563
3564               pos2=buf_pos;
3565               ns2=buf_ns;
3566               list2=buf_list;
3567               col2=buf_col;
3568             }
3569
3570
3571         if ( !hasch1)
3572             {
3573
3574             hasch1=declare_int( (CL->S)->nseq, (CL->S)->max_len+1);
3575             hasch2=declare_int( (CL->S)->nseq, (CL->S)->max_len+1);
3576             n_hasch1=declare_int ( (CL->S)->nseq, (CL->S)->max_len+1);
3577             n_hasch2=declare_int( (CL->S)->nseq, (CL->S)->max_len+1);
3578             is_in_col1=declare_int( (CL->S)->nseq, (CL->S)->max_len+1);
3579             is_in_col2=declare_int( (CL->S)->nseq, (CL->S)->max_len+1);
3580             }
3581
3582         for ( a=0; a< ns1; a++)
3583             {
3584                 rs1= list1[a];
3585                 s1=A->order[rs1][0];
3586                 r1=pos1[rs1][col1];
3587
3588                 if (r1<0)n_gap1++;
3589                 else
3590                    {
3591                    is_in_col1[s1][r1]=1;
3592                    for (b=1; b< CL->residue_index[s1][r1][0]; b+=ICHUNK)
3593                            {
3594                            t_s=CL->residue_index[s1][r1][b+SEQ2];
3595                            t_r=CL->residue_index[s1][r1][b+R2];
3596                            t_w=CL->residue_index[s1][r1][b+WE];
3597                            for ( c=1; c< CL->residue_index[t_s][t_r][0]; c+=ICHUNK)
3598                              {
3599                                q_s=CL->residue_index[t_s][t_r][c+SEQ2];
3600                                q_r=CL->residue_index[t_s][t_r][c+R2];
3601                                q_w=CL->residue_index[t_s][t_r][c+WE];
3602                                hasch1[q_s][q_r]+=MIN(q_w, t_w);
3603                                n_hasch1[q_s][q_r]++;
3604                              }
3605                            }
3606                    }
3607             }
3608
3609         for ( a=0; a< ns2; a++)
3610             {
3611                 rs2=list2[a];
3612                 s2=A->order[rs2][0];
3613                 r2=pos2[rs2][col2];
3614
3615                 if (r2<0)n_gap2++;
3616                 else
3617                    {
3618                    is_in_col2[s2][r2]=1;
3619                    for (b=1; b< CL->residue_index[s2][r2][0]; b+=ICHUNK)
3620                            {
3621                            t_s=CL->residue_index[s2][r2][b+SEQ2];
3622                            t_r=CL->residue_index[s2][r2][b+R2];
3623                            t_w=CL->residue_index[s2][r2][b+WE];
3624                            for ( c=1; c< CL->residue_index[t_s][t_r][0]; c+=ICHUNK)
3625                              {
3626                                q_s=CL->residue_index[t_s][t_r][c+SEQ2];
3627                                q_r=CL->residue_index[t_s][t_r][c+R2];
3628                                q_w=CL->residue_index[t_s][t_r][c+WE];
3629                                hasch2[q_s][q_r]+=MIN(t_w, q_w);
3630                                n_hasch2[q_s][q_r]++;
3631                              }
3632                            }
3633                    }
3634             }
3635
3636
3637         for ( a=0; a< ns2; a++)
3638           {
3639             rs2=list2[a];
3640             s2=A->order[rs2][0];
3641             r2=pos1[rs2][col2];
3642
3643             if (r2<0);
3644             else
3645               {
3646                 for (b=1; b< CL->residue_index[s2][r2][0]; b+=ICHUNK)
3647                   {
3648                     t_s=CL->residue_index[s2][r2][b+SEQ2];
3649                     t_r=CL->residue_index[s2][r2][b+R2];
3650
3651                     for ( c=1; c< CL->residue_index[t_s][t_r][0]; c+=ICHUNK)
3652                       {
3653                         q_s=CL->residue_index[t_s][t_r][c+SEQ2];
3654                         q_r=CL->residue_index[t_s][t_r][c+R2];
3655                         if ( hasch2[q_s][q_r] && hasch1[q_s][q_r]&& !(is_in_col1[q_s][q_r] || is_in_col2[q_s][q_r]))
3656                           {
3657                             score+=MIN(hasch2[q_s][q_r]*(n_hasch1[q_s][q_r]),hasch1[q_s][q_r]*(n_hasch2[q_s][q_r]));
3658                           }
3659                         else if ( hasch2[q_s][q_r] && is_in_col1[q_s][q_r])
3660                           {
3661                             score+=hasch2[q_s][q_r]*(n_hasch1[q_s][q_r]+1);
3662                           }
3663                         else if (hasch1[q_s][q_r] && is_in_col2[q_s][q_r])
3664                           {
3665                             score+=hasch1[q_s][q_r]*(n_hasch2[q_s][q_r]+1);
3666                           }
3667                         hasch2[q_s][q_r]=0;
3668                         n_hasch2[q_s][q_r]=0;
3669                       }
3670                   }
3671                 hasch2[s2][r2]=0;
3672                 is_in_col2[s2][r2]=0;
3673               }
3674           }
3675
3676
3677         for ( a=0; a< ns1; a++)
3678           {
3679             rs1= list1[a];
3680             s1=A->order[rs1][0];
3681             r1=pos1[rs1][col1];
3682
3683             if (r1<0);
3684             else
3685               {
3686                 is_in_col1[s1][r1]=0;
3687                 hasch1[s1][r1]=0;
3688                 for (b=1; b< CL->residue_index[s1][r1][0]; b+=ICHUNK)
3689                   {
3690                     t_s=CL->residue_index[s1][r1][b+SEQ2];
3691                     t_r=CL->residue_index[s1][r1][b+R2];
3692                     for ( c=1; c< CL->residue_index[t_s][t_r][0]; c+=ICHUNK)
3693                       {
3694                         q_s=CL->residue_index[t_s][t_r][c+SEQ2];
3695                         q_r=CL->residue_index[t_s][t_r][c+R2];
3696                         hasch1[q_s][q_r]=0;
3697                         n_hasch1[q_s][q_r]=0;
3698                       }
3699                   }
3700               }
3701           }
3702
3703
3704         score=(score*SCORE_K)/((ns1-n_gap1)*(ns2-n_gap2));
3705         score=(CL->normalise)?((score*CL->normalise)/CL->max_ext_value):score;
3706
3707         return (int)(score-SCORE_K*CL->nomatch);
3708         }
3709
3710
3711 int fast_get_dp_cost ( Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
3712         {
3713         /*WARNING: WORKS ONLY WITH List to Extend*/
3714
3715           double score=0;
3716
3717           int a,b;
3718           int n_gap1=0;
3719           int n_gap2=0;
3720
3721           int s1, rs1, r1, t_r, t_s, s2, rs2, r2;
3722           static int **hasch1;
3723           static int **hasch2;
3724
3725           static int **n_hasch1;
3726           static int **n_hasch2;
3727
3728           static int **is_in_col1;
3729           static int **is_in_col2;
3730
3731
3732
3733
3734         if ( !hasch1)
3735             {
3736
3737             hasch1=declare_int( (CL->S)->nseq, (CL->S)->max_len+1);
3738             hasch2=declare_int( (CL->S)->nseq, (CL->S)->max_len+1);
3739             n_hasch1=declare_int ( (CL->S)->nseq, (CL->S)->max_len+1);
3740             n_hasch2=declare_int( (CL->S)->nseq, (CL->S)->max_len+1);
3741             is_in_col1=declare_int( (CL->S)->nseq, (CL->S)->max_len+1);
3742             is_in_col2=declare_int( (CL->S)->nseq, (CL->S)->max_len+1);
3743             }
3744
3745         for ( a=0; a< ns1; a++)
3746             {
3747                 rs1= list1[a];
3748                 s1=A->order[rs1][0];
3749                 r1=pos1[rs1][col1];
3750
3751                 if (r1<0)n_gap1++;
3752                 else
3753                    {
3754                    is_in_col1[s1][r1]=1;
3755                    for (b=1; b< CL->residue_index[s1][r1][0]; b+=ICHUNK)
3756                            {
3757                            t_s=CL->residue_index[s1][r1][b+SEQ2];
3758                            t_r=CL->residue_index[s1][r1][b+R2];
3759                            hasch1[t_s][t_r]+=CL->residue_index[s1][r1][b+WE];
3760                            n_hasch1[t_s][t_r]++;
3761                            }
3762                    }
3763             }
3764
3765
3766         for ( a=0; a< ns2; a++)
3767             {
3768                 rs2=list2[a];
3769                 s2=A->order[rs2][0];
3770                 r2=pos2[rs2][col2];
3771
3772                 if (r2<0)n_gap2++;
3773                 else
3774                    {
3775                    is_in_col2[s2][r2]=1;
3776                    for (b=1; b< CL->residue_index[s2][r2][0]; b+=ICHUNK)
3777                            {
3778                            t_s=CL->residue_index[s2][r2][b+SEQ2];
3779                            t_r=CL->residue_index[s2][r2][b+R2];
3780
3781                            hasch2[t_s][t_r]+=CL->residue_index[s2][r2][b+WE];
3782                            n_hasch2[t_s][t_r]++;
3783                            }
3784                    }
3785             }
3786         /*return 2;*/
3787
3788         if ( ns2<ns1)
3789             {
3790                 for ( a=0; a< ns2; a++)
3791                     {
3792                     rs2=list2[a];
3793                     s2=A->order[rs2][0];
3794                     r2=pos1[rs2][col2];
3795
3796                     if (r2<0);
3797                     else
3798                         {
3799                         for (b=1; b< CL->residue_index[s2][r2][0]; b+=ICHUNK)
3800                             {
3801                             t_s=CL->residue_index[s2][r2][b+SEQ2];
3802                             t_r=CL->residue_index[s2][r2][b+R2];
3803
3804                             if ( hasch2[t_s][t_r] && hasch1[t_s][t_r]&& !(is_in_col1[t_s][t_r] || is_in_col2[t_s][t_r]))
3805                                 {
3806                                 score+=MIN(hasch2[t_s][t_r]*(n_hasch1[t_s][t_r]),hasch1[t_s][t_r]*(n_hasch2[t_s][t_r]));
3807                                 }
3808                             else if ( hasch2[t_s][t_r] && is_in_col1[t_s][t_r])
3809                                 {
3810                                 score+=hasch2[t_s][t_r]*(n_hasch1[t_s][t_r]+1);
3811                                 }
3812                             else if (hasch1[t_s][t_r] && is_in_col2[t_s][t_r])
3813                                 {
3814                                 score+=hasch1[t_s][t_r]*(n_hasch2[t_s][t_r]+1);
3815                                 }
3816                             hasch2[t_s][t_r]=0;
3817                             n_hasch2[t_s][t_r]=0;
3818                             }
3819                         hasch2[s2][r2]=0;
3820                         is_in_col2[s2][r2]=0;
3821                         }
3822                     }
3823
3824
3825                 for ( a=0; a< ns1; a++)
3826                     {
3827                     rs1= list1[a];
3828                     s1=A->order[rs1][0];
3829                     r1=pos1[rs1][col1];
3830
3831                     if (r1<0);
3832                     else
3833                         {
3834                         is_in_col1[s1][r1]=0;
3835                         hasch1[s1][r1]=0;
3836                         for (b=1; b< CL->residue_index[s1][r1][0]; b+=ICHUNK)
3837                             {
3838                             t_s=CL->residue_index[s1][r1][b+SEQ2];
3839                             t_r=CL->residue_index[s1][r1][b+R2];
3840
3841                             hasch1[t_s][t_r]=0;
3842                             n_hasch1[t_s][t_r]=0;
3843                             }
3844                         }
3845                     }
3846             }
3847         else
3848            {
3849                 for ( a=0; a< ns1; a++)
3850                     {
3851                     rs1=list1[a];
3852                     s1=A->order[rs1][0];
3853                     r1=pos1[rs1][col1];
3854
3855                     if (r1<0);
3856                     else
3857                         {
3858                         for (b=1; b< CL->residue_index[s1][r1][0]; b+=ICHUNK)
3859                             {
3860                             t_s=CL->residue_index[s1][r1][b+SEQ2];
3861                             t_r=CL->residue_index[s1][r1][b+R2];
3862
3863                             if ( hasch1[t_s][t_r] && hasch2[t_s][t_r]&& !(is_in_col2[t_s][t_r] || is_in_col1[t_s][t_r]))
3864                                 {
3865                                 score+=MIN(hasch1[t_s][t_r]*(n_hasch2[t_s][t_r]),hasch2[t_s][t_r]*(n_hasch1[t_s][t_r]));
3866                                 }
3867                             else if ( hasch1[t_s][t_r] && is_in_col2[t_s][t_r])
3868                                 {
3869                                 score+=hasch1[t_s][t_r]*(n_hasch2[t_s][t_r]+1);
3870                                 }
3871                             else if (hasch2[t_s][t_r] && is_in_col1[t_s][t_r])
3872                                 {
3873                                 score+=hasch2[t_s][t_r]*(n_hasch1[t_s][t_r]+1);
3874                                 }
3875                             hasch1[t_s][t_r]=0;
3876                             n_hasch1[t_s][t_r]=0;
3877                             }
3878                         hasch1[s1][r1]=0;
3879                         is_in_col1[s1][r1]=0;
3880                         }
3881                     }
3882
3883
3884                 for ( a=0; a< ns2; a++)
3885                     {
3886                     rs2= list2[a];
3887                     s2=A->order[rs2][0];
3888                     r2=pos1[rs2][col2];
3889
3890                     if (r2<0);
3891                     else
3892                         {
3893                         is_in_col2[s2][r2]=0;
3894                         hasch1[s2][r2]=0;
3895                         for (b=1; b< CL->residue_index[s2][r2][0]; b+=ICHUNK)
3896                             {
3897                             t_s=CL->residue_index[s2][r2][b+SEQ2];
3898                             t_r=CL->residue_index[s2][r2][b+R2];
3899
3900                             hasch2[t_s][t_r]=0;
3901                             n_hasch2[t_s][t_r]=0;
3902                             }
3903                         }
3904                     }
3905             }
3906         score=(score*SCORE_K)/((ns1-n_gap1)*(ns2-n_gap2));
3907         score=(CL->normalise)?((score*CL->normalise)/CL->max_ext_value):score;
3908
3909         return (int)(score-SCORE_K*CL->nomatch);
3910         }
3911
3912 int fast_get_dp_cost_2 ( Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
3913         {
3914         double score=0;
3915
3916         int a, b, s1, s2,r1, r2;
3917         static int n_pair;
3918
3919         int s;
3920         int res_res=0;
3921         int rs1, rs2;
3922         static int last_ns1;
3923         static int last_ns2;
3924         static int last_top1;
3925         static int last_top2;
3926         static int **check_list;
3927
3928
3929         /*New heuristic get dp_cost on a limited number of pairs*/
3930         /*This is the current default*/
3931
3932
3933         if ( last_ns1==ns1 && last_top1==list1[0] && last_ns2==ns2 && last_top2==list2[0]);
3934         else
3935           {
3936
3937
3938             last_ns1=ns1;
3939             last_ns2=ns2;
3940             last_top1=list1[0];
3941             last_top2=list2[0];
3942             if ( check_list) free_int (check_list, -1);
3943             check_list=declare_int ( (CL->S)->nseq*(CL->S)->nseq, 3);
3944
3945             for ( n_pair=0,a=0; a< ns1; a++)
3946               {
3947                 s1 =list1[a];
3948                 rs1=A->order[s1][0];
3949                 for ( b=0; b< ns2; b++, n_pair++)
3950                   {
3951                     s2 =list2[b];
3952                     rs2=A->order[s2][0];
3953                     check_list[n_pair][0]=s1;
3954                     check_list[n_pair][1]=s2;
3955                     check_list[n_pair][2]=(!CL->DM)?0:(CL->DM)->similarity_matrix[rs1][rs2];
3956                   }
3957                 sort_int ( check_list, 3, 2, 0, n_pair-1);
3958               }
3959           }
3960
3961         if ( !CL->evaluate_residue_pair){fprintf ( stderr, "\nWARNING: CL->evaluate_residue_pair Not set\nSet to: extend_residue_pair\n");CL->evaluate_residue_pair= extend_residue_pair;}
3962
3963
3964         for ( a=n_pair-1; a>=0; a--)
3965           {
3966             s1= check_list[a][0];
3967             rs1=A->order[s1][0];
3968             r1 =pos1[s1][col1];
3969
3970             s2= check_list[a][1];
3971             rs2=A->order[s2][0];
3972             r2 =pos2[s2][col2];
3973
3974             if ( r1>0 && r2 >0)
3975               {res_res++;}
3976             if ( rs1>rs2)
3977               {
3978                 SWAP (rs1, rs2);
3979                 SWAP (r1, r2);
3980               }
3981
3982             if ((s=(CL->evaluate_residue_pair)(CL, rs1, r1, rs2, r2))!=UNDEFINED) score+=s;
3983             else
3984               {
3985
3986                 return UNDEFINED;
3987               }
3988             if ( res_res>=CL->max_n_pair && CL->max_n_pair!=0)a=0;
3989           }
3990
3991         score=(res_res==0)?0:( (score)/res_res);
3992         score=score-SCORE_K*CL->nomatch;
3993
3994         return (int)score;
3995         }
3996
3997
3998
3999
4000 int slow_get_dp_cost ( Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
4001         {
4002         double score=0;
4003
4004         int a, b, s1, s2,r1, r2;
4005         int s;
4006         int gap_gap=0;
4007         int gap_res=0;
4008         int res_res=0;
4009         int rs1, rs2;
4010         static int last_tag;
4011         static int *dummy;
4012
4013         if (col1<0 || col2<0) return 0;
4014         if ( last_tag !=A->random_tag)
4015           {
4016             last_tag=A->random_tag;
4017             if (!dummy)
4018               {
4019                 dummy=vcalloc (10, sizeof(int));
4020                 dummy[0]=1;/*Number of Amino acid types on colum*/
4021                 dummy[1]=5;/*Length of Dummy*/
4022                 dummy[3]='\0';/*Amino acid*/
4023                 dummy[4]=1; /*Number of occurences*/
4024                 dummy[5]=100; /*Frequency in the MSA column*/
4025               }
4026           }
4027
4028         if ( !CL->evaluate_residue_pair){fprintf ( stderr, "\nWARNING: CL->evaluate_residue_pair Not set\nSet to: extend_residue_pair\n");CL->evaluate_residue_pair= extend_residue_pair;}
4029
4030         for ( a=0; a< ns1; a++)
4031           {
4032             for ( b=0; b<ns2; b++)
4033               {
4034
4035                 s1 =list1[a];
4036                 rs1=A->order[s1][0];
4037                 r1 =pos1[s1][col1];
4038
4039                 s2 =list2[b];
4040                 rs2=A->order[s2][0];
4041                 r2 =pos2[s2][col2];
4042
4043                 if ( rs1>rs2)
4044                   {
4045                     SWAP (rs1, rs2);
4046                     SWAP (r1, r2);
4047                   }
4048
4049                 if (r1==0 && r2==0)gap_gap++;
4050                 else if ( r1<0 || r2<0) gap_res++;
4051                 else
4052                   {
4053                     res_res++;
4054                     if ((s=(CL->evaluate_residue_pair)(CL, rs1, r1, rs2, r2))!=UNDEFINED) score+=s;
4055                     else
4056                       {
4057
4058                         return UNDEFINED;
4059                       }
4060                   }
4061
4062               }
4063           }
4064
4065
4066         score=(res_res==0)?0:( (score)/res_res);
4067
4068         return score-SCORE_K*CL->nomatch;
4069
4070         }
4071 int slow_get_dp_cost_pc ( Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
4072         {
4073         double score=0;
4074
4075         int a, b, s1, s2,r1, r2;
4076         int s;
4077         int gap_gap=0;
4078         int gap_res=0;
4079         int res_res=0;
4080         int rs1, rs2;
4081         static int last_tag;
4082         static int *dummy;
4083
4084         if (col1<0 || col2<0) return 0;
4085         if ( last_tag !=A->random_tag)
4086           {
4087             last_tag=A->random_tag;
4088             if (!dummy)
4089               {
4090                 dummy=vcalloc (10, sizeof(int));
4091                 dummy[0]=1;/*Number of Amino acid types on colum*/
4092                 dummy[1]=5;/*Length of Dummy*/
4093                 dummy[3]='\0';/*Amino acid*/
4094                 dummy[4]=1; /*Number of occurences*/
4095                 dummy[5]=100; /*Frequency in the MSA column*/
4096               }
4097           }
4098
4099         if ( !CL->evaluate_residue_pair){fprintf ( stderr, "\nWARNING: CL->evaluate_residue_pair Not set\nSet to: extend_residue_pair\n");CL->evaluate_residue_pair= extend_residue_pair;}
4100
4101         for ( a=0; a< ns1; a++)
4102           {
4103             for ( b=0; b<ns2; b++)
4104               {
4105
4106                 s1 =list1[a];
4107                 rs1=A->order[s1][0];
4108                 r1 =pos1[s1][col1];
4109
4110                 s2 =list2[b];
4111                 rs2=A->order[s2][0];
4112                 r2 =pos2[s2][col2];
4113
4114                 if ( rs1>rs2)
4115                   {
4116                     SWAP (rs1, rs2);
4117                     SWAP (r1, r2);
4118                   }
4119
4120                 if (r1==0 && r2==0)gap_gap++;
4121                 else if ( r1<0 || r2<0) gap_res++;
4122                 else
4123                   {
4124                     res_res++;
4125                     if ((s=residue_pair_extended_list_pc(CL, rs1, r1, rs2, r2))!=UNDEFINED) score+=s;
4126                     else
4127                       {
4128
4129                         return UNDEFINED;
4130                       }
4131                   }
4132
4133               }
4134           }
4135
4136
4137         score=(res_res==0)?0:( (score)/res_res);
4138
4139         return score;
4140
4141         }
4142 int slow_get_dp_cost_test ( Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
4143         {
4144         double score=0;
4145
4146         int a, b, s1, s2,r1, r2;
4147         int gap_gap=0, gap_res=0, res_res=0, rs1, rs2;
4148
4149         for ( a=0; a< ns1; a++)
4150                 {
4151                 for ( b=0; b<ns2; b++)
4152                         {
4153
4154                         s1 =list1[a];
4155                         rs1=A->order[s1][0];
4156                         r1 =pos1[s1][col1];
4157
4158                         s2 =list2[b];
4159                         rs2=A->order[s2][0];
4160                         r2 =pos2[s2][col2];
4161
4162                         if ( rs1>rs2)
4163                            {
4164                            SWAP (rs1, rs2);
4165                            SWAP (r1, r2);
4166                            }
4167
4168                         if (r1==0 && r2==0)gap_gap++;
4169                         else if ( r1<0 || r2<0) gap_res++;
4170                         else
4171                             {
4172                             res_res++;
4173                             score+=residue_pair_extended_list_raw (CL, rs1, r1, rs2, r2);
4174                             }
4175                         }
4176                 }
4177
4178         return (int)(score*10)/(ns1*ns2);
4179         }
4180
4181 int sw_get_dp_cost     ( Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
4182         {
4183           int a, b;
4184           int s1,r1,rs1;
4185           int s2,r2,rs2;
4186
4187
4188
4189
4190           for ( a=0; a< ns1; a++)
4191             {
4192               s1 =list1[a];
4193               rs1=A->order[s1][0];
4194               r1 =pos1[s1][col1];
4195               if ( r1<=0)continue;
4196               for ( b=0; b< ns2; b++)
4197                 {
4198
4199
4200                   s2 =list2[b];
4201                   rs2=A->order[s2][0];
4202                   r2 =pos2[s2][col2];
4203
4204                   if (r2<=0)continue;
4205
4206
4207                   if (sw_pair_is_defined (CL, rs1, r1, rs2, r2)==UNDEFINED)return UNDEFINED;
4208                 }
4209             }
4210
4211           return slow_get_dp_cost ( A, pos1, ns1, list1, col1, pos2, ns2, list2,col2, CL);
4212
4213         }
4214
4215
4216
4217
4218
4219
4220
4221
4222 int get_domain_dp_cost ( Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2,Constraint_list *CL , int scale , int gop, int gep)
4223
4224         {
4225         int a, b, s1, s2,r1, r2;
4226         static int *entry;
4227         int *r;
4228         int score=0;
4229         int gap_gap=0;
4230         int gap_res=0;
4231         int res_res=0;
4232         int rs1, rs2;
4233         int flag_list_is_aa_sub_mat=0;
4234         int p;
4235
4236 /*Needs to be cleanned After Usage*/
4237
4238
4239
4240         if ( entry==NULL) entry=vcalloc (LIST_N_FIELDS , sizeof (int));
4241
4242         for (a=0; a< ns1; a++)
4243                 {
4244                 s1=list1[a];
4245                 rs1=A->order[s1][0];
4246                 for ( b=0; b<ns2; b++)
4247                         {
4248                         s2 =list2[b];
4249                         rs2=A->order[s2][0];
4250
4251                         entry[SEQ1]=rs1;
4252                         entry[SEQ2]=rs2;
4253                         r1=entry[R1]=pos1[s1][col1];
4254                         r2=entry[R2]=pos2[s2][col2];
4255
4256                         if ( !flag_list_is_aa_sub_mat)
4257                             {
4258                             if ( r1==r2 && rs1==rs2)
4259                                {
4260
4261                                  return UNDEFINED;
4262                                }
4263                             else if (r1==0 && r2==0)
4264                                {
4265                                gap_gap++;
4266                                }
4267                             else if ( r1<=0 || r2<=0)
4268                                {
4269                                gap_res++;
4270                                }
4271                             else if ((r=main_search_in_list_constraint ( entry,&p,4,CL))!=NULL)
4272                                 {
4273                                 res_res++;
4274
4275                                 if (r[WE]!=UNDEFINED)
4276                                     {
4277                                     score+=(r[WE]*SCORE_K)+scale;
4278                                     }
4279                                 else
4280                                   {
4281                                     fprintf ( stderr, "**");
4282                                     return UNDEFINED;
4283                                   }
4284                                 }
4285                             }
4286                         }
4287                 }
4288         return score;
4289         score=((res_res+gap_res)==0)?0:score/(res_res+gap_res);
4290         return score;
4291         }
4292
4293 /*********************************************************************************************/
4294 /*                                                                                           */
4295 /*         FUNCTIONS FOR ANALYSING AL OR MATRIX                                              */
4296 /*                                                                                           */
4297 /*********************************************************************************************/
4298
4299 int aln2n_res ( Alignment *A, int start, int end)
4300     {
4301     int a, b;
4302     int score=0;
4303
4304     for ( a=start; a<end; a++)for ( b=0; b< A->nseq; b++)score+=!is_gap(A->seq_al[b][a]);
4305     return score;
4306     }
4307
4308 float get_gop_scaling_factor ( int **matrix,float id, int l1, int l2)
4309     {
4310         return id* get_avg_matrix_mm(matrix, AA_ALPHABET);
4311     }
4312
4313 float get_avg_matrix_mm ( int **matrix, char *alphabet)
4314     {
4315         int a, b;
4316         float naa;
4317         float gop;
4318         int l;
4319
4320
4321         l=MIN(20,(int)strlen (alphabet));
4322         for (naa=0, gop=0,a=0; a<l; a++)
4323             for ( b=0; b<l; b++)
4324                 {
4325                     if ( a!=b)
4326                         {
4327                         gop+=matrix[alphabet[a]-'A'][alphabet[b]-'A'];
4328                         naa++;
4329                         }
4330                 }
4331
4332         gop=gop/naa;
4333         return gop;
4334     }
4335 float get_avg_matrix_match ( int **matrix, char *alphabet)
4336     {
4337         int a;
4338         float gop;
4339         int l;
4340
4341
4342
4343
4344         l=MIN(20,(int)strlen (alphabet));
4345         for (gop=0,a=0; a<l; a++)
4346           gop+=matrix[alphabet[a]-'A'][alphabet[a]-'A'];
4347
4348         gop=gop/l;
4349         return gop;
4350     }
4351
4352 float get_avg_matrix_diff ( int **matrix1,int **matrix2, char *alphabet)
4353     {
4354         int a, b;
4355         float naa;
4356         float gop;
4357         int l,v1;
4358
4359
4360
4361
4362         l=MIN(20,(int)strlen (alphabet));
4363         for (naa=0, gop=0,a=0; a<l; a++)
4364           for (b=0; b<l; b++)
4365             {
4366               v1=matrix1[alphabet[a]-'A'][alphabet[a]-'A']-matrix2[alphabet[b]-'A'][alphabet[b]-'A'];
4367               gop+=v1*v1;
4368               naa++;
4369             }
4370
4371         gop=gop/l;
4372         return gop;
4373     }
4374
4375 float* set_aa_frequencies ()
4376    {
4377
4378      float *frequency;
4379      /*frequencies tqken from psw*/
4380
4381      frequency=vcalloc (100, sizeof (float));
4382      frequency ['x'-'A']=0.0013;
4383      frequency ['a'-'A']=0.0076;
4384      frequency ['c'-'A']=0.0176;
4385      frequency ['d'-'A']=0.0529;
4386      frequency ['e'-'A']=0.0628;
4387      frequency ['f'-'A']=0.0401;
4388      frequency ['g'-'A']=0.0695;
4389      frequency ['h'-'A']=0.0224;
4390      frequency ['i'-'A']=0.0561;
4391      frequency ['k'-'A']=0.0584;
4392      frequency ['l'-'A']=0.0922;
4393      frequency ['m'-'A']=0.0236;
4394      frequency ['n'-'A']=0.0448;
4395      frequency ['p'-'A']=0.0500;
4396      frequency ['q'-'A']=0.0403;
4397      frequency ['r'-'A']=0.0523;
4398      frequency ['s'-'A']=0.0715;
4399      frequency ['t'-'A']=0.0581;
4400      frequency ['v'-'A']=0.0652;
4401      frequency ['w'-'A']=0.0128;
4402      frequency ['y'-'A']=0.0321;
4403      frequency ['X'-'A']=0.0013;
4404      frequency ['A'-'A']=0.0076;
4405      frequency ['C'-'A']=0.0176;
4406      frequency ['D'-'A']=0.0529;
4407      frequency ['E'-'A']=0.0628;
4408      frequency ['F'-'A']=0.0401;
4409      frequency ['G'-'A']=0.0695;
4410      frequency ['H'-'A']=0.0224;
4411      frequency ['I'-'A']=0.0561;
4412      frequency ['J'-'A']=0.0584;
4413      frequency ['L'-'A']=0.0922;
4414      frequency ['M'-'A']=0.0236;
4415      frequency ['N'-'A']=0.0448;
4416      frequency ['P'-'A']=0.0500;
4417      frequency ['Q'-'A']=0.0403;
4418      frequency ['R'-'A']=0.0523;
4419      frequency ['S'-'A']=0.0715;
4420      frequency ['T'-'A']=0.0581;
4421      frequency ['V'-'A']=0.0652;
4422      frequency ['W'-'A']=0.0128;
4423      frequency ['Y'-'A']=0.0321;
4424      return frequency;
4425    }
4426
4427 float measure_matrix_pos_avg (int **matrix,char *alphabet)
4428 {
4429         float naa=0, tot=0;
4430         int a, b;
4431
4432
4433         for ( tot=0,a=0; a< 20; a++)
4434            for ( b=0; b<20; b++)
4435              {
4436                if (matrix[alphabet[a]-'A'][alphabet[b]-'A']>=0){naa++;tot+=matrix[alphabet[a]-'A'][alphabet[b]-'A'];}
4437
4438              }
4439         return tot/naa;
4440 }
4441
4442 float measure_matrix_enthropy (int **matrix,char *alphabet)
4443    {
4444
4445      int a, b;
4446      double s, p, q, h=0, tq=0;
4447      float lambda;
4448      float *frequency;
4449      /*frequencies tqken from psw*/
4450
4451      frequency=set_aa_frequencies ();
4452
4453
4454      lambda=compute_lambda(matrix,alphabet);
4455      fprintf ( stderr, "\nLambda=%f", (float)lambda);
4456
4457      for ( a=0; a< 20; a++)
4458        for ( b=0; b<=a; b++)
4459          {
4460            s=matrix[alphabet[a]-'A'][alphabet[b]-'A'];
4461
4462
4463            p=frequency[alphabet[a]-'A']*frequency[alphabet[b]-'A'];
4464
4465            if ( p==0)continue;
4466
4467            q=exp(lambda*s+log(p));
4468
4469            tq+=q;
4470            h+=q*log(q/p)*log(2);
4471
4472          }
4473
4474      fprintf ( stderr,"\ntq=%f\n", (float)tq);
4475
4476      return (float) h;
4477     }
4478 float compute_lambda (int **matrix,char *alphabet)
4479 {
4480
4481   int a, b;
4482   double lambda, best_lambda=0, delta, best_delta=0, p, tq,s;
4483   static float *frequency;
4484
4485   if ( !frequency)frequency=set_aa_frequencies ();
4486
4487   for ( lambda=0.001; lambda<1; lambda+=0.005)
4488      {
4489        tq=0;
4490        for ( a=0; a< 20; a++)
4491          for ( b=0; b<20; b++)
4492            {
4493              p=frequency[alphabet[a]-'A']*frequency[alphabet[b]-'A'];
4494              s=matrix[alphabet[a]-'A'][alphabet[b]-'A'];
4495              tq+=exp(lambda*s+log(p));
4496            }
4497        delta=fabs(1-tq);
4498        if (lambda==0.001)
4499          {
4500            best_delta=delta;
4501            best_lambda=lambda;
4502          }
4503        else
4504          {
4505            if (delta<best_delta)
4506              {
4507                best_delta=delta;
4508                best_lambda=lambda;
4509              }
4510          }
4511
4512        fprintf ( stderr, "\n%f %f ", lambda, tq);
4513        if ( tq>1)break;
4514      }
4515    fprintf ( stderr, "\nRESULT: %f %f ", best_lambda, best_delta);
4516    return (float) best_lambda;
4517 }
4518
4519
4520
4521 float evaluate_random_match (char  *mat, int n, int len,char *alp)
4522 {
4523   int **matrix;
4524   matrix=read_matrice ( mat);
4525   fprintf ( stderr, "Matrix=%15s ", mat);
4526   return evaluate_random_match2 (matrix, n,len,alp);
4527
4528 }
4529
4530 float evaluate_random_match2 (int **matrix, int n, int len,char *alp)
4531 {
4532   int a, b, c, d, c1, c2, tot;
4533   static int *list;
4534   static float *freq;
4535   float score_random=0;
4536   float score_id=0;
4537   float score_good=0;
4538   float tot_len=0;
4539   float tot_try=0;
4540
4541
4542   if ( !list)
4543     {
4544       vsrand(0);
4545       freq=set_aa_frequencies ();
4546       list=vcalloc ( 10000, sizeof (char));
4547     }
4548
4549   for (tot=0,c=0,a=0;a<20; a++)
4550     {
4551       b=freq[alp[a]-'A']*1000;
4552       tot+=b;
4553       for (d=0; d<b; d++, c++)
4554         {
4555           list[c]=alp[a];
4556         }
4557     }
4558
4559
4560   for (a=0; a< len*n; a++)
4561     {
4562       c1=rand()%tot;
4563       c2=rand()%tot;
4564       score_random+=matrix[list[c1]-'A'][list[c2]-'A'];
4565       score_id+=matrix[list[c1]-'A'][list[c1]-'A'];
4566     }
4567   while (tot_len< len*n)
4568     {
4569       tot_try++;
4570       c1=rand()%tot;
4571       c2=rand()%tot;
4572       if ( matrix[list[c1]-'A'][list[c2]-'A']>=0){score_good+=matrix[list[c1]-'A'][list[c2]-'A']; tot_len++;}
4573     }
4574
4575
4576   score_random=score_random/tot_len;
4577   score_id=score_id/tot_len;
4578   score_good=score_good/tot_len;
4579
4580   fprintf ( stderr, "Random=%8.3f Id=%8.3f Good=%8.3f [%7.2f]\n",score_random, score_id, score_good, tot_len/tot_try);
4581
4582   return score_random;
4583 }
4584 float compare_two_mat (char  *mat1,char*mat2, int n, int len,char *alp)
4585 {
4586   int **matrix1, **matrix2;
4587
4588  evaluate_random_match (mat1, n, len,alp);
4589  evaluate_random_match (mat2, n, len,alp);
4590  matrix1=read_matrice ( mat1);
4591  matrix2=read_matrice ( mat2);
4592  matrix1=rescale_matrix(matrix1, 10, alp);
4593  matrix2=rescale_matrix(matrix2, 10, alp);
4594  compare_two_mat_array(matrix1,matrix2, n, len,alp);
4595  return 0;
4596 }
4597
4598
4599 int ** rescale_two_mat (char  *mat1,char*mat2, int n, int len,char *alp)
4600 {
4601   float lambda;
4602   int **matrix1, **matrix2;
4603
4604   lambda=measure_lambda2 (mat1, mat2, n, len, alp)*10;
4605
4606   fprintf ( stderr, "\nLambda=%.2f", lambda);
4607   matrix2=read_matrice(mat2);
4608   matrix2=neg_matrix2pos_matrix(matrix2);
4609   matrix2=rescale_matrix( matrix2, lambda,"abcdefghiklmnpqrstvwxyz");
4610
4611   matrix1=read_matrice(mat1);
4612   matrix1=neg_matrix2pos_matrix(matrix1);
4613   matrix1=rescale_matrix( matrix1,10,"abcdefghiklmnpqrstvwxyz");
4614
4615   output_matrix_header ( "stdout", matrix2, alp);
4616   evaluate_random_match2(matrix1, 1000, 100, alp);
4617   evaluate_random_match2(matrix2, 1000, 100, alp);
4618   compare_two_mat_array(matrix1,matrix2, n, len,alp);
4619
4620   return matrix2;
4621 }
4622 float measure_lambda2(char  *mat1,char*mat2, int n, int len,char *alp)
4623 {
4624   int **m1, **m2;
4625   float f1, f2;
4626
4627   m1=read_matrice (mat1);
4628   m2=read_matrice (mat2);
4629
4630   m1=neg_matrix2pos_matrix(m1);
4631   m2=neg_matrix2pos_matrix(m2);
4632
4633   f1=measure_matrix_pos_avg( m1, alp);
4634   f2=measure_matrix_pos_avg( m2, alp);
4635
4636   return f1/f2;
4637 }
4638
4639
4640 float measure_lambda (char  *mat1,char*mat2, int n, int len,char *alp)
4641 {
4642   int c;
4643   int **matrix1, **matrix2, **mat;
4644   float a;
4645   float best_quality=0, quality=0, best_lambda=0;
4646
4647   matrix1=read_matrice ( mat1);
4648   matrix2=read_matrice ( mat2);
4649   matrix1=rescale_matrix(matrix1, 10, alp);
4650   matrix2=rescale_matrix(matrix2, 10, alp);
4651
4652   for (c=0, a=0.1; a< 2; a+=0.05)
4653     {
4654       fprintf ( stderr, "Lambda=%.2f\n", a);
4655       mat=duplicate_int (matrix2,-1,-1);
4656       mat=rescale_matrix(mat, a, alp);
4657       quality=compare_two_mat_array(matrix1,mat, n, len,alp);
4658       quality=MAX((-quality),quality);
4659
4660       if (c==0 || (best_quality>quality))
4661         {
4662           c=1;
4663           fprintf ( stderr, "*");
4664           best_quality=quality;
4665           best_lambda=a;
4666         }
4667
4668
4669       evaluate_random_match2(mat, 1000, 100, alp);
4670       evaluate_random_match2(matrix1, 1000, 100, alp);
4671       free_int (mat, -1);
4672     }
4673
4674   return best_lambda;
4675
4676 }
4677
4678 float compare_two_mat_array (int  **matrix1,int **matrix2, int n, int len,char *alp)
4679 {
4680   int a, b, c, d, c1, c2, tot;
4681   static int *list;
4682   static float *freq;
4683   float delta_random=0;
4684   float delta2_random=0;
4685
4686   float delta_id=0;
4687   float delta2_id=0;
4688
4689   float delta_good=0;
4690   float delta2_good=0;
4691
4692   float delta;
4693
4694   float tot_len=0;
4695   float tot_try=0;
4696
4697
4698
4699   if ( !list)
4700     {
4701       vsrand(0);
4702       freq=set_aa_frequencies ();
4703       list=vcalloc ( 10000, sizeof (char));
4704     }
4705
4706   for (tot=0,c=0,a=0;a<20; a++)
4707     {
4708       b=freq[alp[a]-'A']*1000;
4709       tot+=b;
4710       for (d=0; d<b; d++, c++)
4711         {
4712           list[c]=alp[a];
4713         }
4714     }
4715
4716
4717   for (a=0; a< len*n; a++)
4718     {
4719       c1=rand()%tot;
4720       c2=rand()%tot;
4721       delta=matrix1[list[c1]-'A'][list[c2]-'A']-matrix2[list[c1]-'A'][list[c2]-'A'];
4722       delta_random+=delta;
4723       delta2_random+=MAX(delta,(-delta));
4724
4725       delta=matrix1[list[c1]-'A'][list[c1]-'A']-matrix2[list[c1]-'A'][list[c1]-'A'];
4726       delta_id+=delta;
4727       delta2_id+=MAX(delta,(-delta));
4728     }
4729   while (tot_len< len*n)
4730     {
4731       tot_try++;
4732       c1=rand()%tot;
4733       c2=rand()%tot;
4734       if ( matrix1[list[c1]-'A'][list[c2]-'A']>=0 || matrix2[list[c1]-'A'][list[c2]-'A'] )
4735         {
4736           delta=matrix1[list[c1]-'A'][list[c2]-'A']-matrix2[list[c1]-'A'][list[c2]-'A'];
4737           delta_good+=delta;
4738           delta2_good+=MAX(delta,(-delta));
4739           tot_len++;
4740         }
4741     }
4742
4743
4744   delta_random=delta_random/tot_len;
4745   delta2_random=delta2_random/tot_len;
4746
4747
4748   delta_id=delta_id/tot_len;
4749   delta2_id=delta2_id/tot_len;
4750
4751   delta_good=delta_good/tot_len;
4752   delta2_good=delta2_good/tot_len;
4753
4754
4755   fprintf ( stderr, "\tRand=%8.3f %8.3f\n\tId  =%8.3f %8.3f\n\tGood=%8.3f %8.3f\n",delta_random, delta2_random, delta_id,delta2_id, delta_good,delta2_good);
4756
4757   return delta_good;
4758 }
4759
4760
4761
4762 int ** rescale_matrix ( int **matrix, float lambda, char *alp)
4763 {
4764   int a, b;
4765
4766
4767   for ( a=0; a< 20; a++)
4768     for ( b=0; b< 20; b++)
4769       {
4770         matrix[alp[a]-'A'][alp[b]-'A']= matrix[alp[a]-'A'][alp[b]-'A']*lambda;
4771       }
4772   return matrix;
4773 }
4774 int **mat2inverted_mat (int **matrix, char *alp)
4775 {
4776   int a, b, min, max, v,l;
4777   int c1,c2, C1, C2;
4778
4779   l=(int)strlen (alp);
4780   min=max=matrix[alp[0]-'A'][alp[0]-'A'];
4781   for ( a=0; a<l; a++)
4782     for ( b=0; b< l; b++)
4783       {
4784         v=matrix[alp[a]-'A'][alp[b]-'A'];
4785         min=MIN(min,v);
4786         max=MAX(max,v);
4787       }
4788   for ( a=0; a<l; a++)
4789     for ( b=0; b< l; b++)
4790       {
4791         v=matrix[alp[a]-'A'][alp[b]-'A'];
4792         v=max-v;
4793
4794         c1=tolower(alp[a])-'A';
4795         c2=tolower(alp[b])-'A';
4796
4797         C1=toupper(alp[a])-'A';
4798         C2=toupper(alp[b])-'A';
4799         matrix[C1][C2]=matrix[c1][c2]=matrix[C1][c2]=matrix[c1][C2]=v;
4800       }
4801   return matrix;
4802 }
4803 void output_matrix_header ( char *name, int **matrix, char *alp)
4804 {
4805   int a, b;
4806   FILE *fp;
4807   char *nalp;
4808   int l;
4809   nalp=vcalloc ( 1000, sizeof (char));
4810
4811
4812
4813   sprintf ( nalp, "ABCDEFGHIKLMNPQRSTVWXYZ");
4814   l=(int)strlen (nalp);
4815
4816   fp=vfopen ( name, "w");
4817   fprintf (fp, "\nint []={\n");
4818
4819   for (a=0; a<l; a++)
4820     {
4821       for ( b=0; b<=a; b++)
4822         fprintf ( fp, "%3d, ",matrix[nalp[a]-'A'][nalp[b]-'A']);
4823       fprintf ( fp, "\n");
4824     }
4825   fprintf (fp, "};\n");
4826
4827   vfclose (fp);
4828 }
4829
4830
4831 float ** initialise_aa_physico_chemical_property_table (int *n)
4832 {
4833   FILE *fp;
4834   int a, b,c;
4835   float **p;
4836   char c1;
4837   float v1, max, min;
4838   int in=0;
4839
4840   n[0]=0;
4841   fp=vfopen ("properties.txt", "r");
4842   while ((c=fgetc(fp))!=EOF)n[0]+=(c=='#');
4843   vfclose (fp);
4844
4845
4846
4847   p=declare_float ( n[0], 'z'+1);
4848   fp=vfopen ("properties.txt", "r");
4849   n[0]=0;
4850   while ((c=fgetc(fp))!=EOF)
4851          {
4852            if (c=='#'){in=0;while (fgetc(fp)!='\n');}
4853            else
4854              {
4855                if ( !in){in=1; ++n[0];}
4856                fscanf (fp, "%f %c %*s",&v1, &c1 );
4857                p[n[0]-1][tolower(c1)]=v1;
4858                while ( (c=fgetc(fp))!='\n');
4859              }
4860          }
4861   vfclose (fp);
4862
4863   /*rescale so that Delta max=10*/
4864   for (a=0; a<n[0]; a++)
4865     {
4866       min=max=p[a]['a'];
4867       for (b='a'; b<='z'; b++)
4868         {
4869           min=(p[a][b]<min)?p[a][b]:min;
4870           max=(p[a][b]>max)?p[a][b]:max;
4871         }
4872       for (b='a'; b<='z'; b++)
4873         {
4874           p[a][b]=((p[a][b]-min)/(max-min))*10;
4875
4876         }
4877     }
4878
4879   return p;
4880 }
4881 Constraint_list * choose_extension_mode ( char *extend_mode, Constraint_list *CL)
4882 {
4883   //evaluation_functions: residues start at 1, sequences at 0;
4884
4885   if ( !CL)
4886     {
4887       fprintf ( stderr, "\nWarning: CL was not set");
4888       return CL;
4889     }
4890   else if ( strm ( extend_mode, "rna0"))
4891     {
4892       CL->evaluate_residue_pair=residue_pair_extended_list;
4893       CL->get_dp_cost      =slow_get_dp_cost;
4894     }
4895   else if ( strm ( extend_mode, "rna1") || strm (extend_mode, "rna"))
4896     {
4897       CL->evaluate_residue_pair=residue_pair_extended_list4rna1;
4898       CL->get_dp_cost      =slow_get_dp_cost;
4899     }
4900   else if ( strm ( extend_mode, "rna2"))
4901     {
4902       CL->evaluate_residue_pair=residue_pair_extended_list4rna2;
4903       CL->get_dp_cost      =slow_get_dp_cost;
4904     }
4905   else if ( strm ( extend_mode, "rna3"))
4906     {
4907       CL->evaluate_residue_pair=residue_pair_extended_list4rna3;
4908       CL->get_dp_cost      =slow_get_dp_cost;
4909     }
4910   else if ( strm ( extend_mode, "rna4"))
4911     {
4912       CL->evaluate_residue_pair=residue_pair_extended_list4rna4;
4913       CL->get_dp_cost      =slow_get_dp_cost;
4914     }
4915   else if ( strm ( extend_mode, "pc") && !CL->M)
4916     {
4917       CL->evaluate_residue_pair=residue_pair_extended_list_pc;
4918       CL->get_dp_cost      =slow_get_dp_cost;
4919     }
4920   else if ( strm ( extend_mode, "triplet") && !CL->M)
4921     {
4922       CL->evaluate_residue_pair=residue_pair_extended_list;
4923       CL->get_dp_cost      =get_dp_cost;
4924     }
4925   else if ( strm ( extend_mode, "relative_triplet") && !CL->M)
4926     {
4927       CL->evaluate_residue_pair=residue_pair_relative_extended_list;
4928       CL->get_dp_cost      =fast_get_dp_cost_2;
4929     }
4930   else if ( strm ( extend_mode, "g_coffee") && !CL->M)
4931     {
4932       CL->evaluate_residue_pair=residue_pair_extended_list_g_coffee;
4933       CL->get_dp_cost      =slow_get_dp_cost;
4934     }
4935   else if ( strm ( extend_mode, "g_coffee_quadruplets") && !CL->M)
4936     {
4937       CL->evaluate_residue_pair=residue_pair_extended_list_g_coffee_quadruplet;
4938       CL->get_dp_cost      =slow_get_dp_cost;
4939     }
4940   else if ( strm ( extend_mode, "fast_triplet") && !CL->M)
4941     {
4942       CL->evaluate_residue_pair=residue_pair_extended_list;
4943       CL->get_dp_cost      =fast_get_dp_cost;
4944     }
4945
4946   else if ( strm ( extend_mode, "very_fast_triplet") && !CL->M)
4947     {
4948       CL->evaluate_residue_pair=residue_pair_extended_list;
4949       CL->get_dp_cost      =fast_get_dp_cost_2;
4950     }
4951   else if ( strm ( extend_mode, "slow_triplet") && !CL->M)
4952     {
4953       CL->evaluate_residue_pair=residue_pair_extended_list;
4954       CL->get_dp_cost      =slow_get_dp_cost;
4955     }
4956   else if (  strm ( extend_mode, "mixt") && !CL->M)
4957     {
4958       CL->evaluate_residue_pair=residue_pair_extended_list_mixt;
4959       CL->get_dp_cost=slow_get_dp_cost;
4960     }
4961   else if (  strm ( extend_mode, "quadruplet") && !CL->M)
4962     {
4963       CL->evaluate_residue_pair=residue_pair_extended_list_quadruplet;
4964       CL->get_dp_cost      =get_dp_cost_quadruplet;
4965     }
4966   else if (  strm ( extend_mode, "test") && !CL->M)
4967     {
4968       CL->evaluate_residue_pair=residue_pair_test_function;
4969       CL->get_dp_cost      =slow_get_dp_cost_test;
4970     }
4971   else if (  strm ( extend_mode, "ssp"))
4972     {
4973
4974       CL->evaluate_residue_pair=evaluate_ssp_matrix_score;
4975       CL->get_dp_cost=slow_get_dp_cost;
4976       CL->normalise=1;
4977     }
4978   else if (  strm ( extend_mode, "tm"))
4979     {
4980
4981       CL->evaluate_residue_pair=evaluate_tm_matrix_score;
4982       CL->get_dp_cost=slow_get_dp_cost;
4983       CL->normalise=1;
4984     }
4985   else if (  strm ( extend_mode, "matrix"))
4986     {
4987
4988       CL->evaluate_residue_pair=evaluate_matrix_score;
4989       CL->get_dp_cost=cw_profile_get_dp_cost;
4990       CL->normalise=1;
4991     }
4992   else if ( strm ( extend_mode, "curvature"))
4993     {
4994        CL->evaluate_residue_pair=evaluate_curvature_score;
4995        CL->get_dp_cost=slow_get_dp_cost;
4996        CL->normalise=1;
4997     }
4998   else if ( CL->M)
4999     {
5000       CL->evaluate_residue_pair=evaluate_matrix_score;
5001       CL->get_dp_cost=cw_profile_get_dp_cost;
5002       CL->normalise=1;
5003     }
5004   else
5005     {
5006       fprintf ( stderr, "\nERROR: %s is an unknown extend_mode[FATAL:%s]\n", extend_mode, PROGRAM);
5007       myexit (EXIT_FAILURE);
5008     }
5009
5010   return CL;
5011 }
5012 int ** combine_two_matrices ( int **mat1, int **mat2)
5013 {
5014   int naa, re1, re2, Re1, Re2, a, b, u, l;
5015
5016   naa=(int)strlen (BLAST_AA_ALPHABET);
5017   for ( a=0; a< naa; a++)
5018     for ( b=0; b< naa; b++)
5019       {
5020         re1=BLAST_AA_ALPHABET[a];
5021         re2=BLAST_AA_ALPHABET[b];
5022         if (re1=='*' || re2=='*');
5023         else
5024           {
5025
5026             Re1=toupper(re1);Re2=toupper(re2);
5027             re1-='A';re2-='A';Re1-='A';Re2-='A';
5028
5029             l=mat1[re1][re2];
5030             u=mat2[re1][re2];
5031             mat1[re1][re2]=mat2[re1][re2]=l;
5032             mat2[Re1][Re2]=mat2[Re1][Re2]=u;
5033           }
5034       }
5035   return mat1;
5036 }
5037
5038 /* Off the shelves evaluations */
5039 /*********************************************************************************************/
5040 /*                                                                                           */
5041 /*         OFF THE SHELVES EVALUATION                                              */
5042 /*                                                                                           */
5043 /*********************************************************************************************/
5044
5045
5046 int lat_sum_pair (Alignment *A, char *mat)
5047 {
5048   int a,b,c, tot=0, v1, v2, score;
5049   int **matrix;
5050
5051   matrix=read_matrice (mat);
5052
5053   for (a=0; a<A->nseq; a++)
5054     for ( b=0; b<A->nseq; b++)
5055       {
5056         for (c=1; c<A->len_aln; c++)
5057           {
5058             char r11, r12;
5059
5060             r11=A->seq_al[a][c-1];
5061             r12=A->seq_al[a][c];
5062             if (is_gap(r11) || is_gap(r12))continue;
5063             else v1=matrix[r11-'A'][r12-'A'];
5064
5065             r11=A->seq_al[b][c-1];
5066             r12=A->seq_al[b][c];
5067             if (is_gap(r11) || is_gap(r12))continue;
5068             else v2=matrix[r11-'A'][r12-'A'];
5069
5070             score+=(v1-v2)*(v1-v2);
5071             tot++;
5072           }
5073       }
5074   score=(100*score)/tot;
5075   return (float)score;
5076 }
5077
5078
5079
5080 /* Off the shelves evaluations */
5081 /*********************************************************************************************/
5082 /*                                                                                           */
5083 /*         OFF THE SHELVES EVALUATION                                              */
5084 /*                                                                                           */
5085 /*********************************************************************************************/
5086
5087 int comp_pair ( int len,char *sa, char *sb, int seqA, int seqB,int *tgp_a, int *tgp_b,int gap_op,int gap_ex, int start, int end,int **matrix,int MODE);
5088 int score_gap ( int len, char *sa, char *sb,int seqA, int seqB,int *tgp_a, int *tgp_b,  int op, int ex, int start, int end, int MODE);
5089 void evaluate_tgp_decoded_chromosome ( Alignment *A,int **TGP,int start, int end,int MODE);
5090 int gap_type ( char a, char b);
5091
5092
5093
5094 float  sum_pair ( Alignment*A,char *mat_name, int gap_op, int gap_ex)
5095     {
5096         int a,b;
5097         float pscore=0;
5098
5099         int start, end;
5100         static int **tgp;
5101         double score=0;
5102         int MODE=1;
5103         int **matrix;
5104
5105         matrix=read_matrice (mat_name);
5106         matrix=mat2inverted_mat (matrix, "acdefghiklmnpqrstvwy");
5107
5108         start=0;
5109         end=A->len_aln-1;
5110
5111         if ( tgp==NULL)
5112             tgp= declare_int (A->nseq,2);
5113
5114         evaluate_tgp_decoded_chromosome ( A,tgp,start, end,MODE);
5115
5116         for ( a=0; a< A->nseq-1; a++)
5117                  for (b=a+1; b<A->nseq; b++)
5118                         {
5119                           pscore= comp_pair (A->len_aln,A->seq_al[a], A->seq_al[b],a, b,tgp[a], tgp[b],gap_op,gap_ex, start, end,matrix, MODE);
5120                           score+=pscore*100;
5121                           /*score+=(double)pscore*(int)(PARAM->OFP)->weight[A->order[a][0]][A->order[b][0]];*//*NO WEIGHTS*/
5122                         }
5123
5124         score=score/(A->nseq*A->nseq);
5125
5126         return (float)score;
5127         }
5128
5129 int comp_pair ( int len,char *sa, char *sb, int seqA, int seqB,int *tgp_a, int *tgp_b,int gap_op,int gap_ex, int start, int end,int **matrix,int MODE)
5130         {
5131           int score=0, a, ex;
5132
5133
5134
5135         if ( end-start>=0)
5136             score+= score_gap (len, sa,sb, seqA, seqB,tgp_a, tgp_b, gap_op,gap_ex, start, end,MODE);
5137
5138         ex=gap_ex;
5139
5140
5141         for (a=start; a<=end; a++)
5142                 {
5143                   if ( is_gap(sa[a]) || is_gap(sb[a]))
5144                     {
5145                       if (is_gap(sa[a]) && is_gap(sb[a]));
5146                       else
5147                         {
5148
5149                           score +=ex;
5150                         }
5151                     }
5152                   else
5153                         {
5154                         score += matrix [sa[a]-'A'][sb[a]-'A'];
5155
5156                         }
5157                 }
5158         return score;
5159         }
5160 int score_gap ( int len, char *sa, char *sb,int seqA, int seqB,int *tgp_a, int *tgp_b,  int op, int ex, int start, int end, int MODE)
5161         {
5162           int a,b;
5163         int ga=0,gb=0;
5164         int score=0;
5165
5166
5167         int right_gap, left_gap;
5168
5169
5170
5171
5172
5173         int type;
5174         int flag1=0;
5175         int flag2=0;
5176         int continue_loop;
5177         int sequence_pattern[2][3];
5178         int null_gap;
5179         int natural_gap=1;
5180
5181         /*op= gor_gap_op ( 0,seqA, seqB, PARAM);
5182         ex= gor_gap_ext ( 0, seqA, seqB, PARAM);*/
5183
5184
5185
5186         for (a=start; a<=end; ++a)
5187                 {
5188
5189                 type= gap_type ( sa[a], sb[a]);
5190
5191                 if ( type==2 && ga<=gb)
5192                         {++ga;
5193                          gb=0;
5194                          score += op;
5195                         }
5196                 else if (type==1 && ga >=gb)
5197                         {
5198                         ++gb;
5199                         ga=0;
5200                         score +=op;
5201                         }
5202                 else if (type==0)
5203                         {
5204                         ga++;
5205                         gb++;
5206                         }
5207
5208                 else if (type== -1)
5209                         ga=gb=0;
5210
5211
5212                 if (natural_gap==0)
5213                     {
5214                     if ( type== -1)
5215                         flag1=flag2=0;
5216                     else if ( type==0)
5217                         flag2=1;
5218                     else if ( (type==flag1) && flag2==1)
5219                         {
5220                         score+=op;
5221                         flag2=0;
5222                         }
5223                     else if ( (type!=flag1) && flag2==1)
5224                         {
5225                         flag1=type;
5226                         flag2=0;
5227                         }
5228                     else if ( flag2==0)
5229                         flag1=type;
5230                     }
5231                  }
5232           /*gap_type -/-:0, X/X:-1 X/-:1, -/X:2*/
5233 /*evaluate the pattern of gaps*/
5234
5235         continue_loop=1;
5236         sequence_pattern[0][0]=sequence_pattern[1][0]=0;
5237         for ( a=start; a<=end && continue_loop==1; a++)
5238             {
5239             left_gap= gap_type ( sa[a], sb[a]);
5240             if ( left_gap!= 0)
5241                 {
5242                 if ( left_gap==-1)
5243                     {
5244                     sequence_pattern[0][0]=sequence_pattern[1][0]=0;
5245                     continue_loop=0;
5246                     }
5247                 else
5248                     {
5249                     null_gap=0;
5250                     for (b=a; b<=end && continue_loop==1; b++)
5251                         {type=gap_type( sa[b], sb[b]);
5252                         if (type==0)
5253                             null_gap++;
5254                         if ( type!=left_gap && type !=0)
5255                             {
5256                             continue_loop=0;
5257                             sequence_pattern[2-left_gap][0]= b-a-null_gap;
5258                             sequence_pattern [1-(2-left_gap)][0]=0;
5259                             }
5260                          }
5261                      if ( continue_loop==1)
5262                         {
5263                         continue_loop=0;
5264                         sequence_pattern[2-left_gap][0]= b-a-null_gap;
5265                         sequence_pattern [1-(2-left_gap)][0]=0;
5266                         }
5267                      }
5268                   }
5269                }
5270
5271            sequence_pattern[0][2]=sequence_pattern[1][2]=1;
5272            for ( a=start; a<=end; a++)
5273                 {
5274                   if ( !is_gap(sa[a]))
5275                     sequence_pattern[0][2]=0;
5276                   if ( !is_gap(sb[a]))
5277                     sequence_pattern[1][2]=0;
5278
5279                 }
5280            continue_loop=1;
5281            sequence_pattern[0][1]=sequence_pattern[1][1]=0;
5282            for ( a=end; a>=start && continue_loop==1; a--)
5283             {
5284             right_gap= gap_type ( sa[a], sb[a]);
5285             if ( right_gap!= 0)
5286                 {
5287                 if ( right_gap==-1)
5288                     {
5289                     sequence_pattern[0][1]=sequence_pattern[1][1]=0;
5290                     continue_loop=0;
5291                     }
5292                 else
5293                     {
5294                     null_gap=0;
5295                     for (b=a; b>=start && continue_loop==1; b--)
5296                         {type=gap_type( sa[b], sb[b]);
5297                         if ( type==0)
5298                             null_gap++;
5299                         if ( type!=right_gap && type !=0)
5300                             {
5301                             continue_loop=0;
5302                             sequence_pattern[2-right_gap][1]= a-b-null_gap;
5303                             sequence_pattern [1-(2-right_gap)][1]=0;
5304                             }
5305                         }
5306                      if ( continue_loop==1)
5307                         {
5308                         continue_loop=0;
5309                         sequence_pattern[2-right_gap][1]= a-b-null_gap;
5310                         sequence_pattern [1-(2-right_gap)][1]=0;
5311                         }
5312                      }
5313                   }
5314                }
5315
5316 /*
5317 printf ( "\n*****************************************************");
5318 printf ( "\n%c\n%c", sa[start],sb[start]);
5319 printf ( "\n%d %d %d",sequence_pattern[0][0] ,sequence_pattern[0][1], sequence_pattern[0][2]);
5320 printf ( "\n%d %d %d",sequence_pattern[1][0] ,sequence_pattern[1][1], sequence_pattern[1][2]);
5321 printf ( "\n*****************************************************");
5322 */
5323
5324 /*correct the scoring*/
5325
5326
5327         if ( MODE==0)
5328             {
5329             if ( FABS(tgp_a[0])>1 && (FABS(tgp_a[0])>FABS( tgp_b[0])))
5330                 score-= (sequence_pattern[0][0]>0)?op:0;
5331             if ( FABS(tgp_b[0])>1 && (FABS(tgp_b[0])> FABS(tgp_a[0])))
5332                 score-= (sequence_pattern[1][0]>0)?op:0;
5333             }
5334         else if ( MODE ==1 || MODE ==2)
5335             {
5336             if ( FABS(tgp_a[0])>1 && (FABS(tgp_a[0])>FABS( tgp_b[0])) && (tgp_a[1]!=1 || sequence_pattern[0][2]==0))
5337                 score-= (sequence_pattern[0][0]>0)?op:0;
5338             if ( FABS(tgp_b[0])>1 && (FABS(tgp_b[0])> FABS(tgp_a[0])) && (tgp_b[1]!=1 || sequence_pattern[1][2]==0))
5339                 score-= (sequence_pattern[1][0]>0)?op:0;
5340
5341
5342             if ( tgp_a[0]>=1 && tgp_a[0]==tgp_b[0])
5343                 score -=(sequence_pattern[0][0]>0)?op:0;
5344             if ( tgp_b[0]>=1 && tgp_a[0]==tgp_b[0])
5345                 score-= (sequence_pattern[1][0]>0)?op:0;
5346
5347
5348             if ( tgp_a[1]==1 && sequence_pattern[0][2]==0)
5349                 score -= ( sequence_pattern[0][1]>0)?op:0;
5350             else if ( tgp_a[1]==1 && sequence_pattern[0][2]==1 && tgp_a[0]<=0)
5351                 score -= ( sequence_pattern[0][1]>0)?op:0;
5352
5353
5354             if ( tgp_b[1]==1 && sequence_pattern[1][2]==0)
5355                 score -= ( sequence_pattern[1][1]>0)?op:0;
5356             else if ( tgp_b[1]==1 && sequence_pattern[1][2]==1 && tgp_b[0]<=0)
5357                 score -= ( sequence_pattern[1][1]>0)?op:0;
5358
5359             if ( MODE==2)
5360                 {
5361                 if ( tgp_a[0]>0)
5362                     score -=sequence_pattern[0][0]*ex;
5363                 if ( tgp_b[0]>0)
5364                     score -= sequence_pattern[1][0]*ex;
5365                 if ( tgp_a[1]>0)
5366                     score-=sequence_pattern[0][1]*ex;
5367                 if ( tgp_b[1]>0)
5368                     score-=sequence_pattern[1][1]*ex;
5369                 }
5370             }
5371
5372
5373         return score;
5374
5375
5376
5377         }
5378 void evaluate_tgp_decoded_chromosome ( Alignment *A,int **TGP,int start, int end,int MODE)
5379     {
5380     int a,b;
5381     int continue_loop;
5382
5383
5384
5385     if (MODE==11 || MODE==13|| MODE==14)
5386         {
5387         if ( start==0)for ( a=0; a<A->nseq; a++)TGP[a][0]=-1;
5388         else for ( a=0; a<A->nseq; a++)TGP[a][0]=(is_gap(A->seq_al[a][start-1])==1)?0:1;
5389
5390         if ( end==A->len_aln-1)for ( a=0; a<A->nseq; a++)TGP[a][1]=-1;
5391         else for ( a=0; a<A->nseq; a++)TGP[a][1]=(is_gap(A->seq_al[a][start-1])==1)?0:1;
5392         }
5393     else
5394         {
5395         /* 0: in the middle of the alignement
5396         1: natural end
5397         2: q left gap is the continuation of another gap that was open outside the bloc ( don't open it)
5398         */
5399
5400         for ( a=0; a< A->nseq; a++)
5401           {
5402             TGP[a][0]=1;
5403             TGP[a][1]=1;
5404             for ( b=0; b< start; b++)
5405               if ( !is_gap(A->seq_al[a][b]))
5406                 TGP[a][0]=0;
5407             if ( start>0 )
5408               {
5409                 if (is_gap(A->seq_al[a][start-1]) && TGP[a][0]!=1)
5410                   {TGP[a][0]=-1;
5411                     continue_loop=1;
5412                     for ( b=(start-1); b>=0 && continue_loop==1; b--)
5413                       {TGP[a][0]-= ( is_gap(A->seq_al[a][b])==1)?1:0;
5414                         continue_loop= (is_gap(A->seq_al[a][b])==1)?continue_loop:0;
5415                       }
5416                   }
5417               }
5418               else if (is_gap(A->seq_al[a][start-1]) && TGP[a][0]==1)
5419                 {
5420                   TGP[a][0]=1;
5421                   continue_loop=1;
5422                   for ( b=(start-1); b>=0 && continue_loop==1; b--)
5423                     {TGP[a][0]+= ( is_gap(A->seq_al[a][b])==1)?1:0;
5424                       continue_loop= (is_gap(A->seq_al[a][b])==1)?continue_loop:0;
5425                     }
5426                 }
5427             for ( b=(A->len_aln-1); b>end; b--)
5428               if ( !is_gap(A->seq_al[a][b]))
5429                 TGP[a][1]=0;
5430           }
5431         }
5432     }
5433 int gap_type ( char a, char b)
5434     {
5435     /*gap_type -/-:0, X/X:-1 X/-:1, -/STAR:2*/
5436
5437       if ( is_gap(a) && is_gap(b))
5438         return 0;
5439       else if ( !is_gap(a) && !is_gap(b))
5440         return -1;
5441       else if ( !is_gap(a))
5442         return 1;
5443       else if ( !is_gap(b))
5444         return 2;
5445       else
5446         return -1;
5447     }
5448
5449 /******************************COPYRIGHT NOTICE*******************************/
5450 /*© Centro de Regulacio Genomica */
5451 /*and */
5452 /*Cedric Notredame */
5453 /*Fri Feb 18 08:27:45 CET 2011 - Revision 596. */
5454 /*All rights reserved.*/
5455 /*This file is part of T-COFFEE.*/
5456 /**/
5457 /*    T-COFFEE is free software; you can redistribute it and/or modify*/
5458 /*    it under the terms of the GNU General Public License as published by*/
5459 /*    the Free Software Foundation; either version 2 of the License, or*/
5460 /*    (at your option) any later version.*/
5461 /**/
5462 /*    T-COFFEE is distributed in the hope that it will be useful,*/
5463 /*    but WITHOUT ANY WARRANTY; without even the implied warranty of*/
5464 /*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the*/
5465 /*    GNU General Public License for more details.*/
5466 /**/
5467 /*    You should have received a copy of the GNU General Public License*/
5468 /*    along with Foobar; if not, write to the Free Software*/
5469 /*    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA*/
5470 /*...............................................                                                                                      |*/
5471 /*  If you need some more information*/
5472 /*  cedric.notredame@europe.com*/
5473 /*...............................................                                                                                                                                     |*/
5474 /**/
5475 /**/
5476 /*      */
5477 /******************************COPYRIGHT NOTICE*******************************/