src/BIONJ.c

   1 /*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   2 ;                                                                           ;
   3 ;                         BIONJ program                                     ;
   4 ;                                                                           ;
   5 ;                         Olivier Gascuel                                   ;
   6 ;                                                                           ;
   7 ;                         GERAD - Montreal- Canada                          ;
   8 ;                         olivierg@crt.umontreal.ca                         ;
   9 ;                                                                           ;
  10 ;                         LIRMM - Montpellier- France                       ;
  11 ;                         gascuel@lirmm.fr                                  ;
  12 ;                                                                           ;
  13 ;                         UNIX version, written in C                        ;
  14 ;                         by Hoa Sien Cuong (Univ. Montreal)                ;
  15 ;                                                                           ;
  16 \*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*/
  17
  18 //#include "BIONJ.h"
  19 #include "me.h"
  20
  21 void Initialize(float **delta, double *X, char **labels, int n, POINTERS *trees);
  22 void Print_outputChar(int i, POINTERS *trees, char *output);
  23 void bionj (double *X, int *N, char **labels, char **treeStr);
  24 int Symmetrize(float **delta, int n);
  25 void Concatenate(char chain1[MAX_LABEL_LENGTH], int ind, POINTERS *trees, int post);
  26 float Distance(int i, int j, float **delta);
  27 float Variance(int i, int j, float **delta);
  28 int Emptied(int i, float **delta);
  29 float Sum_S(int i, float **delta);
  30 void Compute_sums_Sx(float **delta, int n);
  31 void Best_pair(float **delta, int r, int *a, int *b, int n);
  32 float Finish_branch_length(int i, int j, int k, float **delta);
  33 void FinishStr (float **delta, int n, POINTERS *trees, char *StrTree);
  34 float Agglomerative_criterion(int i, int j, float **delta, int r);
  35 float Branch_length(int a, int b, float **delta, int r);
  36 float Reduction4(int a, float la, int b, float lb, int i, float lamda, float **delta);
  37 float Reduction10(int a, int b, int i, float lamda, float vab, float **delta);
  38 float Lamda(int a, int b, float vab, float **delta, int n, int r);
  39
  40 /* void printMat(float **delta, int n); */
  41
  42 /*;;;;;;;;;;;  INPUT, OUTPUT, INITIALIZATION ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  43 ;                                                                           ;
  44 ;                                                                           ;
  45 ;              The delta matrix is read from the input-file.                ;
  46 ;              It is recommended to put it and the executable in            ;
  47 ;              a special directory. The input-file and output-file          ;
  48 ;              can be given as arguments to the executable by               ;
  49 ;              typing them after the executable (Bionj input-file           ;
  50 ;              output-file) or by typing them when asked by the             ;
  51 ;              program. The input-file has to be formated according         ;
  52 ;              the PHYLIP standard. The output file is formated             ;
  53 ;              according to the NEWWICK standard.                           ;
  54 ;                                                                           ;
  55 ;              The lower-half of the delta matrix is occupied by            ;
  56 ;              dissimilarities. The upper-half of the matrix is             ;
  57 ;              occupied by variances. The first column                      ;
  58 ;              is initialized as 0; during the algorithm some               ;
  59 ;              indices are no more used, and the corresponding              ;
  60 ;              positions in the first column are set to 1.                  ;
  61 ;                                                                           ;
  62 ;              This delta matix is made symmetrical using the rule:         ;
  63 ;              Dij = Dji <- (Dij + Dji)/2. The diagonal is set to 0;        ;
  64 ;              during the further steps of the algorithm, it is used        ;
  65 ;              to store the sums Sx.                                        ;
  66 ;                                                                           ;
  67 ;              A second array, trees, is used to store taxon names.         ;
  68 ;              During the further steps of the algoritm, some               ;
  69 ;              positions in this array are emptied while the others         ;
  70 ;              are used to store subtrees.                                  ;
  71 ;                                                                           ;
  72 \*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*/
  73
  74
  75 /*;;;;;;;;;;;;;;;;;;;;;;;;;; Initialize        ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*\
  76 ;                                                                           ;
  77 ; Description : This function reads an input file and return the            ;
  78 ;               delta matrix and trees: the list of taxa.                   ;
  79 ;                                                                           ;
  80 ; input       :                                                             ;
  81 ;              float **delta : delta matrix                                 ;
  82 ;              FILE *input    : pointer to input file                       ;
  83 ;              int n          : number of taxa                              ;
  84 ;              char **trees   : list of taxa                                ;
  85 ;                                                                           ;
  86 ; return value:                                                             ;
  87 ;              float **delta : delta matrix                                 ;
  88 ;              char *trees    : list of taxa                                ;
  89 ;                                                                           ;
  90 \*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*/
  91
  92 //void Initialize(float **delta, FILE *input, int n, POINTERS *trees)
  93 void Initialize(float **delta, double *X, char **labels, int n, POINTERS *trees)
  94 {
  95   int lig;                                          /* matrix line       */
  96   int col;                                          /* matrix column     */
  97 //  float distance;
  98   //char name_taxon[LEN];                             /* taxon name        */
  99   char name_taxon[MAX_LABEL_LENGTH];
 100   char format[MAX_DIGITS];
 101   WORD *name;
 102
 103   snprintf (format, MAX_DIGITS, "%%%ds", MAX_LABEL_LENGTH);
 104
 105   for(lig=1; lig <= n; lig++)
 106     {
 107       //fscanf(input,"%s",name_taxon);                  /* read taxon name   */
 108       //fscanf (input, format, name_taxon);             /* read taxon name   */
 109       strncpy (name_taxon, labels[lig-1], MAX_LABEL_LENGTH);
 110       name=(WORD *)calloc(1,sizeof(WORD));            /* taxon name is     */
 111       if(name == NULL)                                /* put in trees      */
 112         {
 113           Rprintf("Out of memories !!");
 114           exit(0);
 115         }
 116       else
 117         {
 118           strncpy (name->name, name_taxon, MAX_LABEL_LENGTH);
 119           name->suiv=NULL;
 120           trees[lig].head=name;
 121           trees[lig].tail=name;
 122           for(col=lig; col <= n; col++)
 123             {
 124               //fscanf(input,"%f",&distance);             /* read the distance  */
 125 //            &distance = X[XINDEX(lig,col)];
 126               delta[col][lig]=X[XINDEX(lig,col)];
 127               delta[lig][col]=X[XINDEX(lig,col)];
 128               if (lig==col)
 129                 delta[lig][col]=0;
 130             }
 131         }
 132     }
 133   return;
 134 }
 135
 136 /*;;;;;;;;;;;;;;;;;;;;;;;;;;; Print_output;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*\
 137 ;                                                                           ;
 138 ;                                                                           ;
 139 ; Description : This function prints out the tree in the output file.       ;
 140 ;                                                                           ;
 141 ; input       :                                                             ;
 142 ;              POINTERS *trees : pointer to the subtrees.                   ;
 143 ;              int i          : indicate the subtree i to be printed.       ;
 144 :              FILE *output   : pointer to the output file.                 ;
 145 ;                                                                           ;
 146 ; return value: The phylogenetic tree in the output file.                   ;
 147 ;                                                                           ;
 148 \*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*/
 149
 150 void Print_outputChar(int i, POINTERS *trees, char *output)
 151 {
 152   WORD *parcour;
 153   parcour=trees[i].head;
 154   while (parcour != NULL && (strlen (output) + strlen (parcour->name) < MAX_INPUT_SIZE))
 155     {
 156       output = strncat (output, parcour->name, strlen (parcour->name));
 157       parcour=parcour->suiv;
 158     }
 159   return;
 160 }
 161
 162 /*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 163 ;                                                                           ;
 164 ;                         Main program                                      ;
 165 ;                                                                           ;
 166 ;                         argc is the number of arguments                   ;
 167 ;                         **argv contains the arguments:                    ;
 168 ;                         the first argument has to be BIONJ;               ;
 169 ;                         the second is the inptu-file;                     ;
 170 ;                         the third is the output-file.                     ;
 171 ;                         When the input and output files are               ;
 172 ;                         not given, the user is asked for them.            ;
 173 ;                                                                           ;
 174 \*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*/
 175
 176
 177 //tree *bionj (FILE *input, boolean isNJ)
 178 void bionj (double *X, int *N, char **labels, char **treeStr)
 179 {
 180   POINTERS *trees;                          /* list of subtrees            */
 181   tree *T;                                  /* the returned tree           */
 182   char *chain1;                             /* stringized branch-length    */
 183 //  char *chain2;                             /* idem                        */
 184   char *str;                                /* the string containing final tree */
 185   int *a, *b;                               /* pair to be agglomerated     */
 186   float **delta;                            /* delta matrix                */
 187   float la;                                 /* first taxon branch-length   */
 188   float lb;                                 /* second taxon branch-length  */
 189   float vab;                                /* variance of Dab             */
 190   float lamda = 0.5;
 191   int i;
 192 //  int ok;
 193   int r;                                    /* number of subtrees          */
 194   int n;                                    /* number of taxa              */
 195   int x, y;
 196 //  double t;
 197   a=(int*)calloc(1,sizeof(int));
 198   b=(int*)calloc(1,sizeof(int));
 199   chain1=(char *)calloc(MAX_LABEL_LENGTH,sizeof(char));
 200   str = (char *)R_alloc(MAX_INPUT_SIZE, sizeof(char));
 201   /* added by EP */
 202   if (strlen(str))
 203     strncpy(str, "", strlen(str));
 204   /* end */
 205
 206 //  fscanf(input,"%d",&n);
 207   n = *N;
 208   /*      Create the delta matrix     */
 209   delta=(float **)calloc(n+1,sizeof(float*));
 210   for(i=1; i<= n; i++)
 211     {
 212       delta[i]=(float *)calloc(n+1, sizeof(float));
 213       if(delta[i] == NULL)
 214         {
 215           Rprintf("Out of memories!!");
 216           exit(0);
 217         }
 218     }
 219   trees=(POINTERS *)calloc(n+1,sizeof(POINTERS));
 220   if(trees == NULL)
 221     {
 222       Rprintf("Out of memories!!");
 223       exit(0);
 224     }
 225   /*   initialise and symmetrize the running delta matrix    */
 226     r=n;
 227     *a=0;
 228     *b=0;
 229     Initialize(delta, X, labels, n, trees);
 230 //    ok=Symmetrize(delta, n);
 231
 232 //    if(!ok)
 233 //     Rprintf("\n The matrix is not symmetric.\n ");
 234     while (r > 3)                             /* until r=3                 */
 235       {
 236         Compute_sums_Sx(delta, n);             /* compute the sum Sx       */
 237         Best_pair(delta, r, a, b, n);          /* find the best pair by    */
 238         vab=Variance(*a, *b, delta);           /* minimizing (1)           */
 239         la=Branch_length(*a, *b, delta, r);    /* compute branch-lengths   */
 240         lb=Branch_length(*b, *a, delta, r);    /* using formula (2)        */
 241 //      if (!isNJ)
 242           lamda=Lamda(*a, *b, vab, delta, n, r); /* compute lambda* using (9)*/
 243         for(i=1; i <= n; i++)
 244           {
 245             if(!Emptied(i,delta) && (i != *a) && (i != *b))
 246               {
 247                 if(*a > i)
 248                   {
 249                     x=*a;
 250                     y=i;
 251                   }
 252                 else
 253                   {
 254                     x=i;
 255                     y=*a;                           /* apply reduction formulae */
 256                   }                                 /* 4 and 10 to delta        */
 257                 delta[x][y]=Reduction4(*a, la, *b, lb, i, lamda, delta);
 258                 delta[y][x]=Reduction10(*a, *b, i, lamda, vab, delta);
 259               }
 260           }
 261         strncpy(chain1,"",1);                  /* agglomerate the subtrees */
 262         strncat(chain1,"(",1);                 /* a and b together with the*/
 263         Concatenate(chain1, *a, trees, 0);     /* branch-lengths according */
 264         strncpy(chain1,"",1);                  /* to the NEWWICK format    */
 265         strncat(chain1,":",1);
 266         snprintf(chain1+strlen(chain1),MAX_LABEL_LENGTH,"%f",la);
 267
 268         strncat(chain1,",",1);
 269         Concatenate(chain1,*a, trees, 1);
 270         trees[*a].tail->suiv=trees[*b].head;
 271         trees[*a].tail=trees[*b].tail;
 272         strncpy(chain1,"",1);
 273         strncat(chain1,":",1);
 274         snprintf(chain1+strlen(chain1),MAX_LABEL_LENGTH,"%f",lb);
 275
 276         strncat(chain1,")",1);
 277         Concatenate(chain1, *a, trees, 1);
 278         delta[*b][0]=1.0;                     /* make the b line empty     */
 279         trees[*b].head=NULL;
 280         trees[*b].tail=NULL;
 281         r=r-1;
 282       }
 283
 284     FinishStr (delta, n, trees, str);       /* compute the branch-lengths*/
 285                                             /* of the last three subtrees*/
 286                                             /* and print the tree in the */
 287                                             /* output-file               */
 288     T = readNewickString (str, n);
 289     T = detrifurcate(T);
 290 //    compareSets(T,species);
 291     partitionSizes(T);
 292     *treeStr = (char *)R_alloc(MAX_INPUT_SIZE, sizeof(char));
 293     /* added by EP */
 294     if (strlen(*treeStr))
 295       strncpy(*treeStr, "", strlen(*treeStr));
 296     /* end */
 297     NewickPrintTreeStr (T, *treeStr);
 298
 299     for(i=1; i<=n; i++)
 300       {
 301         delta[i][0]=0.0;
 302         trees[i].head=NULL;
 303         trees[i].tail=NULL;
 304       }
 305   free(delta);
 306   free(trees);
 307   /* free (str); */
 308   free (chain1);
 309   free (a);
 310   free (b);
 311   freeTree(T);
 312   T = NULL;
 313 //  return (ret);
 314   return;
 315 }
 316
 317 /*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*\
 318 ;                                                                           ;
 319 ;                             Utilities                                     ;
 320 ;                                                                           ;
 321 \*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*/
 322
 323
 324
 325 /*;;;;;;;;;;;;;;;;;;;;;;;;;;; Symmetrize  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*\
 326 ;                                                                           ;
 327 ; Description : This function verifies if the delta matrix is symmetric;    ;
 328 ;               if not the matrix is made symmetric.                        ;
 329 ;                                                                           ;
 330 ; input       :                                                             ;
 331 ;              float **delta : delta matrix                                 ;
 332 ;              int n          : number of taxa                              ;
 333 ;                                                                           ;
 334 ; return value:                                                             ;
 335 ;              int symmetric  : indicate if the matrix has been made        ;
 336 ;                               symmetric or not                            ;
 337 ;                                                                           ;
 338 \*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*/
 339
 340 int Symmetrize(float **delta, int n)
 341 {
 342   int lig;                                         /* matrix line        */
 343   int col;                                         /* matrix column      */
 344   float value;                                     /* symmetrized value  */
 345   int symmetric;
 346
 347   symmetric=1;
 348   for(lig=1; lig  <=  n; lig++)
 349     {
 350       for(col=1; col < lig; col++)
 351         {
 352           if(delta[lig][col] != delta[col][lig])
 353             {
 354               value= (delta[lig][col]+delta[col][lig])/2;
 355               delta[lig][col]=value;
 356               delta[col][lig]=value;
 357               symmetric=0;
 358             }
 359         }
 360     }
 361   return(symmetric);
 362 }
 363
 364
 365
 366
 367 /*;;;;;;;;;;;;;;;;;;;;;;;;;;; Concatenate ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*\
 368 ;                                                                           ;
 369 ;                                                                           ;
 370 ; Description : This function concatenates a string to another.             ;
 371 ;                                                                           ;
 372 ; input       :                                                             ;
 373 ;      char *chain1    : the string to be concatenated.                     ;
 374 ;      int ind         : indicate the subtree to which concatenate the      ;
 375 ;                        string                                             ;
 376 ;      POINTERS *trees  : pointer to subtrees.                              ;
 377 ;      int post        : position to which concatenate (front (0) or        ;
 378 ;                        end (1))                                           ;
 379 ;                                                                           ;
 380 \*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*/
 381
 382 //void Concatenate(char chain1[LEN], int ind, POINTERS *trees, int post)
 383 void Concatenate(char chain1[MAX_LABEL_LENGTH], int ind, POINTERS *trees, int post)
 384 {
 385   WORD *bran;
 386
 387   bran=(WORD *)calloc(1,sizeof(WORD));
 388   if(bran == NULL)
 389     {
 390       Rprintf("Out of memories");
 391       exit(0);
 392     }
 393   else
 394     {
 395       strncpy(bran->name,chain1,MAX_LABEL_LENGTH);
 396       bran->suiv=NULL;
 397     }
 398   if(post == 0)
 399     {
 400       bran->suiv=trees[ind].head;
 401       trees[ind].head=bran;
 402     }
 403   else
 404     {
 405       trees[ind].tail->suiv=bran;
 406       trees[ind].tail=trees[ind].tail->suiv;
 407     }
 408 }
 409
 410
 411 /*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Distance;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*\
 412 ;                                                                           ;
 413 ; Description : This function retrieve ant return de distance between taxa  ;
 414 ;               i and j from the delta matrix.                              ;
 415 ;                                                                           ;
 416 ; input       :                                                             ;
 417 ;               int i          : taxon i                                    ;
 418 ;               int j          : taxon j                                    ;
 419 ;               float **delta : the delta matrix                            ;
 420 ;                                                                           ;
 421 ; return value:                                                             ;
 422 ;               float distance : dissimilarity between the two taxa         ;
 423 ;                                                                           ;
 424 \*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*/
 425
 426 float Distance(int i, int j, float **delta)
 427 {
 428   if(i > j)
 429     return(delta[i][j]);
 430   else
 431     return(delta[j][i]);
 432 }
 433
 434
 435 /*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Variance;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*\
 436 ;                                                                           ;
 437 ; Description : This function retrieve and return the variance of the       ;
 438 ;               distance between i and j, from the delta matrix.            ;
 439 ;                                                                           ;
 440 ; input       :                                                             ;
 441 ;               int i           : taxon i                                   ;
 442 ;               int j           : taxon j                                   ;
 443 ;               float **delta  : the delta matrix                           ;
 444 ;                                                                           ;
 445 ; return value:                                                             ;
 446 ;               float distance : the variance of  Dij                       ;
 447 ;                                                                           ;
 448 \*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*/
 449
 450 float Variance(int i, int j, float **delta)
 451 {
 452   if(i > j)
 453     return(delta[j][i]);
 454   else
 455     return(delta[i][j]);
 456 }
 457
 458
 459 /*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Emptied ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*\
 460 ;                                                                           ;
 461 ; Description : This function verifie if a line is emptied or not.          ;
 462 ;                                                                           ;
 463 ; input       :                                                             ;
 464 ;               int i          : subtree (or line) i                        ;
 465 ;               float **delta : the delta matrix                            ;
 466 ;                                                                           ;
 467 ; return value:                                                             ;
 468 ;               0              : if not emptied.                            ;
 469 ;               1              : if emptied.                                ;
 470 ;                                                                           ;
 471 \*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*/
 472
 473 int Emptied(int i, float **delta)      /* test if the ith line is emptied */
 474 {
 475   return((int)delta[i][0]);
 476 }
 477
 478
 479 /*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Sum_S;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*\
 480 ;                                                                           ;
 481 ;  Description : This function retrieves the sum Sx from the diagonal       ;
 482 ;                of the delta matrix.                                       ;
 483 ;                                                                           ;
 484 ;  input       :                                                            ;
 485 ;               int i          : subtree i                                  ;
 486 ;               float **delta : the delta matrix                            ;
 487 ;                                                                           ;
 488 ;  return value:                                                            ;
 489 ;                float delta[i][i] : sum Si                                 ;
 490 ;                                                                           ;
 491 \*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*/
 492
 493 float Sum_S(int i, float **delta)          /* get sum Si form the diagonal */
 494 {
 495   return(delta[i][i]);
 496 }
 497
 498
 499 /*;;;;;;;;;;;;;;;;;;;;;;;Compute_sums_Sx;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*\
 500 ;                                                                           ;
 501 ; Description : This function computes the sums Sx and store them in the    ;
 502 ;               diagonal the delta matrix.                                  ;
 503 ;                                                                           ;
 504 ; input       :                                                             ;
 505 ;                float **delta : the delta matrix.                      ;
 506 ;                int n          : the number of taxa                    ;
 507 ;                                                                           ;
 508 \*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*/
 509
 510 void Compute_sums_Sx(float **delta, int n)
 511 {
 512   float sum;
 513   int i;
 514   int j;
 515
 516   for(i= 1; i <= n ; i++)
 517     {
 518       if(!Emptied(i,delta))
 519         {
 520           sum=0;
 521           for(j=1; j <=n; j++)
 522             {
 523               if(i != j && !Emptied(j,delta))           /* compute the sum Si */
 524                 sum=sum + Distance(i,j,delta);
 525             }
 526         }
 527       delta[i][i]=sum;                           /* store the sum Si in */
 528     }                                               /* delta\92s diagonal    */
 529 }
 530
 531
 532 /*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Best_pair;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*\
 533 ;                                                                           ;
 534 ;  Description : This function finds the best pair to be agglomerated by    ;
 535 ;                minimizing the agglomerative criterion (1).                ;
 536 ;                                                                           ;
 537 ;  input       :                                                            ;
 538 ;                float **delta : the delta matrix                           ;
 539 ;                int r          : number of subtrees                        ;
 540 ;                int *a         : contain the first taxon of the pair       ;
 541 ;                int *b         : contain the second taxon of the pair      ;
 542 ;                int n          : number of taxa                            ;
 543 ;                                                                           ;
 544 ;  return value:                                                            ;
 545 ;                int *a         : the first taxon of the pair               ;
 546 ;                int *b         : the second taxon of the pair              ;
 547 \*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*/
 548
 549 void Best_pair(float **delta, int r, int *a, int *b, int n)
 550 {
 551   float Qxy;                         /* value of the criterion calculated*/
 552   int x,y;                           /* the pair which is tested         */
 553   float Qmin;                        /* current minimun of the criterion */
 554
 555   Qmin=1.0e300;
 556   for(x=1; x <= n; x++)
 557     {
 558       if(!Emptied(x,delta))
 559         {
 560           for(y=1; y < x; y++)
 561             {
 562               if(!Emptied(y,delta))
 563                 {
 564                   Qxy=Agglomerative_criterion(x,y,delta,r);
 565                   if(Qxy < Qmin-0.000001)
 566                     {
 567                       Qmin=Qxy;
 568                       *a=x;
 569                       *b=y;
 570                     }
 571                 }
 572             }
 573         }
 574     }
 575 }
 576
 577
 578 /*;;;;;;;;;;;;;;;;;;;;;;Finish_branch_length;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*\
 579 ;                                                                           ;
 580 ;  Description :  Compute the length of the branch attached                 ;
 581 ;                 to the subtree i, during the final step                   ;
 582 ;                                                                           ;
 583 ;  input       :                                                            ;
 584 ;                int i          : position of subtree i                     ;
 585 ;                int j          : position of subtree j                     ;
 586 ;                int k          : position of subtree k                     ;
 587 ;                float **delta :                                            ;
 588 ;                                                                           ;
 589 ;  return value:                                                            ;
 590 ;                float length  : The length of the branch                   ;
 591 ;                                                                           ;
 592 \*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*/
 593
 594 float Finish_branch_length(int i, int j, int k, float **delta)
 595 {
 596   float length;
 597   length=0.5*(Distance(i,j,delta) + Distance(i,k,delta)
 598               -Distance(j,k,delta));
 599   return(length);
 600 }
 601
 602
 603 /*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Finish;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*\
 604 ;                                                                           ;
 605 ;  Description : This function compute the length of the lasts three        ;
 606 ;                subtrees and write the tree in the output file.            ;
 607 ;                                                                           ;
 608 ;  input       :                                                            ;
 609 ;                float **delta  : the delta matrix                          ;
 610 ;                int n           : the number of taxa                       ;
 611 ;                WORD *trees   : list of subtrees                           ;
 612 ;                                                                           ;
 613 ;                                                                           ;
 614 \*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*/
 615 /*
 616 void Finish(float **delta, int n, POINTERS *trees, FILE *output)
 617 {
 618   int l=1;
 619   int i=0;
 620   float length;
 621   WORD *bidon;
 622   WORD *ele;
 623   int last[3];                            // the last three subtrees
 624
 625   while(l <= n)
 626     {                                     // find the last tree subtree
 627       if(!Emptied(l, delta))
 628         {
 629           last[i]=l;
 630           i++;
 631         }
 632       l++;
 633     }
 634
 635   length=Finish_branch_length(last[0],last[1],last[2],delta);
 636   fprintf(output,"(");
 637   Print_output(last[0],trees,output);
 638   fprintf(output,":");
 639 //   gcvt(length,PREC, str);
 640 //   fprintf(output,"%s,",str);
 641   fprintf(output,"%f,",length);
 642
 643   length=Finish_branch_length(last[1],last[0],last[2],delta);
 644   Print_output(last[1],trees,output);
 645   fprintf(output,":");
 646 //   gcvt(length,PREC, str);
 647 //   fprintf(output,"%s,",str);
 648   fprintf(output,"%f,",length);
 649
 650   length=Finish_branch_length(last[2],last[1],last[0],delta);
 651   Print_output(last[2],trees,output);
 652   fprintf(output,":");
 653 //   gcvt(length,PREC,str);
 654 //   fprintf(output,"%s",str);
 655   fprintf(output,"%f",length);
 656   fprintf(output,");");
 657   fprintf(output,"\n");
 658
 659   for(i=0; i < 3; i++)
 660     {
 661       bidon=trees[last[i]].head;
 662       ele=bidon;
 663       while(bidon!=NULL)
 664         {
 665           ele=ele->suiv;
 666           free(bidon);
 667           bidon=ele;
 668         }
 669     }
 670 }
 671 */
 672
 673 void FinishStr (float **delta, int n, POINTERS *trees, char *StrTree)
 674 {
 675   int l=1;
 676   int i=0;
 677   float length;
 678   char *tmp;
 679   WORD *bidon;
 680   WORD *ele;
 681   int last[3];                            /* the last three subtrees     */
 682
 683   while(l <= n)
 684     {                                     /* find the last tree subtree  */
 685       if(!Emptied(l, delta))
 686         {
 687           last[i]=l;
 688           i++;
 689         }
 690       l++;
 691     }
 692   tmp = (char*) calloc (12, sizeof(char));
 693   StrTree[0]='(';
 694
 695   length=Finish_branch_length(last[0],last[1],last[2],delta);
 696   Print_outputChar (last[0], trees, StrTree);
 697   snprintf (tmp, 12, "%f,", length);
 698   if (strlen (StrTree) + strlen (tmp) < MAX_INPUT_SIZE-1) {
 699     strncat (StrTree, ":", 1);
 700     strncat (StrTree, tmp, strlen (tmp));
 701   }
 702
 703   length=Finish_branch_length(last[1],last[0],last[2],delta);
 704   Print_outputChar (last[1], trees, StrTree);
 705   snprintf (tmp, 12, "%f,", length);
 706   if (strlen (StrTree) + strlen (tmp) < MAX_INPUT_SIZE-1) {
 707     strncat (StrTree, ":", 1);
 708     strncat (StrTree, tmp, strlen (tmp));
 709   }
 710
 711   length=Finish_branch_length(last[2],last[1],last[0],delta);
 712   Print_outputChar (last[2], trees, StrTree);
 713   snprintf (tmp, 12, "%f", length);
 714   if (strlen (StrTree) + strlen (tmp) < MAX_INPUT_SIZE-1) {
 715     strncat (StrTree, ":", 1);
 716     strncat (StrTree, tmp, strlen (tmp));
 717   }
 718
 719   if (strlen (StrTree) < MAX_INPUT_SIZE-3)
 720     strncat (StrTree, ");\n", 3);
 721
 722   free (tmp);
 723   for(i=0; i < 3; i++)
 724     {
 725       bidon=trees[last[i]].head;
 726       ele=bidon;
 727       while(bidon!=NULL)
 728         {
 729           ele=ele->suiv;
 730           free(bidon);
 731           bidon=ele;
 732         }
 733     }
 734   return;
 735 }
 736
 737 /*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*\
 738 ;                                                                           ;
 739 ;                          Formulae                                         ;
 740 ;                                                                           ;
 741 \*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*/
 742
 743
 744 float Agglomerative_criterion(int i, int j, float **delta, int r)
 745 {
 746   float Qij;
 747   Qij=(r-2)*Distance(i,j,delta)                           /* Formula (1) */
 748     -Sum_S(i,delta)
 749     -Sum_S(j,delta);
 750
 751   return(Qij);
 752 }
 753
 754
 755 float Branch_length(int a, int b, float **delta, int r)
 756 {
 757   float length;
 758   length=0.5*(Distance(a,b,delta)                         /* Formula (2) */
 759               +(Sum_S(a,delta)
 760                 -Sum_S(b,delta))/(r-2));
 761   return(length);
 762 }
 763
 764
 765 float Reduction4(int a, float la, int b, float lb, int i, float lamda, float **delta)
 766 {
 767   float Dui;
 768   Dui=lamda*(Distance(a,i,delta)-la)
 769     +(1-lamda)*(Distance(b,i,delta)-lb);                /* Formula (4) */
 770   return(Dui);
 771 }
 772
 773
 774 float Reduction10(int a, int b, int i, float lamda, float vab, float **delta)
 775 {
 776   float Vci;
 777   Vci=lamda*Variance(a,i,delta)+(1-lamda)*Variance(b,i,delta)
 778     -lamda*(1-lamda)*vab;                              /*Formula (10)  */
 779   return(Vci);
 780 }
 781
 782 float Lamda(int a, int b, float vab, float **delta, int n, int r)
 783 {
 784   float lamda=0.0;
 785   int i;
 786
 787   if(vab==0.0)
 788     lamda=0.5;
 789   else
 790     {
 791       for(i=1; i <= n ; i++)
 792         {
 793           if(a != i && b != i && !Emptied(i,delta))
 794             lamda=lamda + (Variance(b,i,delta) - Variance(a,i,delta));
 795         }
 796       lamda=0.5 + lamda/(2*(r-2)*vab);
 797     }                                              /* Formula (9) and the  */
 798   if(lamda > 1.0)                                  /* constraint that lamda*/
 799     lamda = 1.0;                                   /* belongs to [0,1]     */
 800   if(lamda < 0.0)
 801     lamda=0.0;
 802   return(lamda);
 803 }
 804
 805
 806 /* void printMat(float **delta, int n) */
 807 /* { */
 808 /*   int i, j; */
 809 /*   for (i=1; i<=n; i++) { */
 810 /*     for (j=1; j<=n; j++) */
 811 /*       Rprintf ("%f ", delta[i][j]); */
 812 /*     Rprintf("\n"); */
 813 /*   } */
 814 /*   Rprintf("\n"); */
 815 /*   return; */
 816 /* } */