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