src/BIONJ.c

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