src/me.c

   1 /* me.c    2008-01-14 */
   2
   3 /* Copyright 2007-2008 Olivier Gascuel, Rick Desper,
   4    R port by Vincent Lefort, me_*() 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 #include "me.h"
  10
  11 //functions from me_balanced.c
  12 tree *BMEaddSpecies(tree *T, node *v, double **D, double **A);
  13 void assignBMEWeights(tree *T, double **A);
  14 void makeBMEAveragesTable(tree *T, double **D, double **A);
  15 //functions from me_ols.c
  16 tree *GMEaddSpecies(tree *T, node *v, double **D, double **A);
  17 void assignOLSWeights(tree *T, double **A);
  18 void makeOLSAveragesTable(tree *T, double **D, double **A);
  19 //functions from bNNI.c
  20 void bNNI(tree *T, double **avgDistArray, int *count, double **D, int numSpecies);
  21 //functions from NNI.c
  22 void NNI(tree *T, double **avgDistArray, int *count, double **D, int numSpecies);
  23
  24
  25 void me_b(double *X, int *N, char **labels, int *nni,
  26           int *edge1, int *edge2, double *el, char **tl)
  27 {
  28   double **D, **A;
  29   set *species, *slooper;
  30   node *addNode;
  31   tree *T;
  32   int n, nniCount;
  33
  34   n = *N;
  35   T = NULL;
  36   nniCount = 0;
  37   species = (set *) malloc(sizeof(set));
  38   species->firstNode = NULL;
  39   species->secondNode = NULL;
  40   D = loadMatrix(X, labels, n, species);
  41   A = initDoubleMatrix(2*n - 2);
  42
  43   for(slooper = species; NULL != slooper; slooper = slooper->secondNode)
  44   {
  45     addNode = copyNode(slooper->firstNode);
  46     T = BMEaddSpecies(T, addNode, D, A);
  47   }
  48   // Compute bNNI
  49   if (*nni == 1)
  50     bNNI(T, A, &nniCount, D, n);
  51   assignBMEWeights(T,A);
  52
  53   tree2phylo(T, edge1, edge2, el, tl, n);
  54
  55   freeMatrix(D,n);
  56   freeMatrix(A,2*n - 2);
  57   freeSet(species);
  58   freeTree(T);
  59   T = NULL;
  60 }
  61
  62 void me_o(double *X, int *N, char **labels, int *nni,
  63           int *edge1, int *edge2, double *el, char **tl)
  64 {
  65   double **D, **A;
  66   set *species, *slooper;
  67   node *addNode;
  68   tree *T;
  69   int n, nniCount;
  70
  71   n = *N;
  72   T = NULL;
  73   nniCount = 0;
  74   species = (set *) malloc(sizeof(set));
  75   species->firstNode = NULL;
  76   species->secondNode = NULL;
  77
  78   D = loadMatrix (X, labels, n, species);
  79   A = initDoubleMatrix(2 * n - 2);
  80
  81   for(slooper = species; NULL != slooper; slooper = slooper->secondNode)
  82   {
  83     addNode = copyNode(slooper->firstNode);
  84     T = GMEaddSpecies(T,addNode,D,A);
  85   }
  86   makeOLSAveragesTable(T,D,A);
  87   // Compute NNI
  88   if (*nni == 1)
  89     NNI(T,A,&nniCount,D,n);
  90   assignOLSWeights(T,A);
  91
  92   tree2phylo(T, edge1, edge2, el, tl, n);
  93
  94   freeMatrix(D,n);
  95   freeMatrix(A,2*n - 2);
  96   freeSet(species);
  97   freeTree(T);
  98   T = NULL;
  99 }
 100
 101 /*************************************************************************
 102
 103                            MATRIX FUNCTIONS
 104
 105 *************************************************************************/
 106
 107 double **initDoubleMatrix(int d)
 108 {
 109   int i,j;
 110   double **A;
 111   A = (double **) malloc(d*sizeof(double *));
 112   for(i=0;i<d;i++)
 113     {
 114       A[i] = (double *) malloc(d*sizeof(double));
 115       for(j=0;j<d;j++)
 116         A[i][j] = 0.0;
 117     }
 118   return(A);
 119 }
 120
 121 double **loadMatrix (double *X, char **labels, int n, set *S)
 122 {
 123   char nextString[MAX_LABEL_LENGTH];
 124   node *v;
 125   double **table;
 126   int i, j, a, b;
 127
 128   table = (double **) calloc(n,sizeof(double *));
 129   for(i=0; i<n; i++)
 130     table[i] = (double *) calloc(n,sizeof(double));
 131
 132   for(i=0; i<n; i++)
 133     {
 134       strncpy (nextString, labels[i], MAX_LABEL_LENGTH);
 135 //      ReplaceForbiddenChars (nextString, '_');
 136       v = makeNewNode(nextString,-1);
 137       v->index2 = i;
 138       S = addToSet(v,S);
 139       for (j=i; j<n; j++) {
 140         a=i+1;
 141         b=j+1;
 142         table[j][i] = X[XINDEX(a,b)];
 143         table[i][j] = X[XINDEX(a,b)];
 144         if (i==j)
 145           table[i][j] = 0;
 146       }
 147     }
 148   return (table);
 149 }
 150
 151 /*************************************************************************
 152
 153                            GRAPH FUNCTIONS
 154
 155 *************************************************************************/
 156
 157 set *addToSet(node *v, set *X)
 158 {
 159   if (NULL == X)
 160     {
 161       X = (set *) malloc(sizeof(set));
 162       X->firstNode = v;
 163       X->secondNode = NULL;
 164     }
 165   else if (NULL == X->firstNode)
 166     X->firstNode = v;
 167   else
 168     X->secondNode = addToSet(v,X->secondNode);
 169   return(X);
 170 }
 171
 172 node *makeNewNode(char *label, int i)
 173 {
 174   return(makeNode(label,NULL,i));
 175 }
 176
 177 node *makeNode(char *label, edge *parentEdge, int index)
 178 {
 179   node *newNode;  /*points to new node added to the graph*/
 180   newNode = (node *) malloc(sizeof(node));
 181   strncpy(newNode->label,label,NODE_LABEL_LENGTH);
 182   newNode->index = index;
 183   newNode->index2 = -1;
 184   newNode->parentEdge = parentEdge;
 185   newNode->leftEdge = NULL;
 186   newNode->middleEdge = NULL;
 187   newNode->rightEdge = NULL;
 188   /*all fields have been initialized*/
 189   return(newNode);
 190 }
 191
 192 /*copyNode returns a copy of v which has all of the fields identical to those
 193 of v, except the node pointer fields*/
 194 node *copyNode(node *v)
 195 {
 196   node *w;
 197   w = makeNode(v->label,NULL,v->index);
 198   w->index2 = v->index2;
 199   return(w);
 200 }
 201
 202 edge *siblingEdge(edge *e)
 203 {
 204   if(e == e->tail->leftEdge)
 205     return(e->tail->rightEdge);
 206   else
 207     return(e->tail->leftEdge);
 208 }
 209
 210 edge *makeEdge(char *label, node *tail, node *head, double weight)
 211 {
 212   edge *newEdge;
 213   newEdge = (edge *) malloc(sizeof(edge));
 214   strncpy(newEdge->label,label,EDGE_LABEL_LENGTH);
 215   newEdge->tail = tail;
 216   newEdge->head = head;
 217   newEdge->distance = weight;
 218   newEdge->totalweight = 0.0;
 219   return(newEdge);
 220 }
 221
 222 tree *newTree()
 223 {
 224   tree *T;
 225   T = (tree *) malloc(sizeof(tree));
 226   T->root = NULL;
 227   T->size = 0;
 228   T->weight = -1;
 229   return(T);
 230 }
 231
 232 void updateSizes(edge *e, int direction)
 233 {
 234   edge *f;
 235   switch(direction)
 236     {
 237     case UP:
 238       f = e->head->leftEdge;
 239       if (NULL != f)
 240         updateSizes(f,UP);
 241       f = e->head->rightEdge;
 242       if (NULL != f)
 243         updateSizes(f,UP);
 244       e->topsize++;
 245       break;
 246     case DOWN:
 247       f = siblingEdge(e);
 248       if (NULL != f)
 249         updateSizes(f,UP);
 250       f = e->tail->parentEdge;
 251       if (NULL != f)
 252         updateSizes(f,DOWN);
 253       e->bottomsize++;
 254       break;
 255     }
 256 }
 257
 258 /*detrifurcate takes the (possibly trifurcated) input tree
 259   and reroots the tree to a leaf*/
 260 /*assumes tree is only trifurcated at root*/
 261 tree *detrifurcate(tree *T)
 262 {
 263   node *v, *w;
 264   edge *e, *f;
 265   v = T->root;
 266   if(leaf(v))
 267     return(T);
 268   if (NULL != v->parentEdge)
 269     {
 270       Rprintf ("Error: root %s is poorly rooted.\n",v->label);
 271       exit(0);
 272     }
 273   for(e = v->middleEdge, v->middleEdge = NULL; NULL != e; e = f )
 274     {
 275       w = e->head;
 276       v = e->tail;
 277       e->tail = w;
 278       e->head = v;
 279       f = w->leftEdge;
 280       v->parentEdge = e;
 281       w->leftEdge = e;
 282       w->parentEdge = NULL;
 283     }
 284   T->root = w;
 285   return(T);
 286 }
 287
 288 void compareSets(tree *T, set *S)
 289 {
 290   edge *e;
 291   node *v,*w;
 292   set *X;
 293   e = depthFirstTraverse(T,NULL);
 294   while (NULL != e)
 295     {
 296       v = e->head;
 297       for(X = S; NULL != X; X = X->secondNode)
 298         {
 299           w = X->firstNode;
 300           if (0 == strcmp(v->label,w->label))
 301             {
 302               v->index2 = w->index2;
 303             w->index2 = -1;
 304             break;
 305             }
 306         }
 307       e = depthFirstTraverse(T,e);
 308     }
 309   v = T->root;
 310   for(X = S; NULL != X; X = X->secondNode)
 311     {
 312       w = X->firstNode;
 313       if (0 == strcmp(v->label,w->label))
 314         {
 315           v->index2 = w->index2;
 316           w->index2 = -1;
 317           break;
 318         }
 319     }
 320   if (-1 == v->index2)
 321     {
 322       Rprintf("Error leaf %s in tree not in distance matrix.\n",v->label);
 323       exit(0);
 324     }
 325   e = depthFirstTraverse(T,NULL);
 326   while (NULL != e)
 327     {
 328       v = e->head;
 329       if ((leaf(v)) && (-1 == v->index2))
 330         {
 331           Rprintf("Error leaf %s in tree not in distance matrix.\n",v->label);
 332           exit(0);
 333         }
 334       e = depthFirstTraverse(T,e);
 335       }
 336   for(X = S; NULL != X; X = X->secondNode)
 337     if (X->firstNode->index2 > -1)
 338       {
 339         Rprintf("Error node %s in matrix but not a leaf in tree.\n",X->firstNode->label);
 340         exit(0);
 341       }
 342   return;
 343 }
 344
 345 void partitionSizes(tree *T)
 346 {
 347   edge *e;
 348   e = depthFirstTraverse(T,NULL);
 349   while (NULL != e)
 350     {
 351       if (leaf(e->head))
 352         e->bottomsize = 1;
 353       else
 354         e->bottomsize = e->head->leftEdge->bottomsize
 355           + e->head->rightEdge->bottomsize;
 356       e->topsize = (T->size + 2)/2 - e->bottomsize;
 357       e = depthFirstTraverse(T,e);
 358     }
 359 }
 360
 361 /*************************************************************************
 362
 363                            TRAVERSE FUNCTIONS
 364
 365 *************************************************************************/
 366
 367 edge *depthFirstTraverse(tree *T, edge *e)
 368      /*depthFirstTraverse returns the edge f which is least in T according
 369        to the depth-first order, but which is later than e in the search
 370        pattern.  If e is null, f is the least edge of T*/
 371 {
 372   edge *f;
 373   if (NULL == e)
 374     {
 375       f = T->root->leftEdge;
 376       if (NULL != f)
 377         f = findBottomLeft(f);
 378       return(f);  /*this is the first edge of this search pattern*/
 379     }
 380   else /*e is non-null*/
 381     {
 382       if (e->tail->leftEdge == e)
 383         /*if e is a left-oriented edge, we skip the entire
 384           tree cut below e, and find least edge*/
 385         f = moveRight(e);
 386       else  /*if e is a right-oriented edge, we have already looked at its
 387               sibling and everything below e, so we move up*/
 388         f = e->tail->parentEdge;
 389     }
 390   return(f);
 391 }
 392
 393 edge *findBottomLeft(edge *e)
 394      /*findBottomLeft searches by gottom down in the tree and to the left.*/
 395 {
 396   edge *f;
 397   f = e;
 398   while (NULL != f->head->leftEdge)
 399     f = f->head->leftEdge;
 400   return(f);
 401 }
 402
 403 edge *moveRight(edge *e)
 404 {
 405   edge *f;
 406   f = e->tail->rightEdge; /*this step moves from a left-oriented edge
 407                             to a right-oriented edge*/
 408   if (NULL != f)
 409     f = findBottomLeft(f);
 410   return(f);
 411 }
 412
 413 edge *topFirstTraverse(tree *T, edge *e)
 414      /*topFirstTraverse starts from the top of T, and from there moves stepwise
 415        down, left before right*/
 416      /*assumes tree has been detrifurcated*/
 417 {
 418   edge *f;
 419   if (NULL == e)
 420     return(T->root->leftEdge); /*first Edge searched*/
 421   else if (!(leaf(e->head)))
 422     return(e->head->leftEdge); /*down and to the left is preferred*/
 423   else /*e->head is a leaf*/
 424     {
 425       f = moveUpRight(e);
 426       return(f);
 427     }
 428 }
 429
 430 edge *moveUpRight(edge *e)
 431 {
 432   edge *f;
 433   f = e;
 434   while ((NULL != f) && ( f->tail->leftEdge != f))
 435     f = f->tail->parentEdge;
 436   /*go up the tree until f is a leftEdge*/
 437   if (NULL == f)
 438     return(f); /*triggered at end of search*/
 439   else
 440     return(f->tail->rightEdge);
 441   /*and then go right*/
 442 }
 443
 444 /*************************************************************************
 445
 446                            FREE FUNCTIONS
 447
 448 *************************************************************************/
 449
 450 void freeMatrix(double **D, int size)
 451 {
 452   int i;
 453   for(i=0;i<size;i++)
 454     free(D[i]);
 455   free(D);
 456 }
 457
 458 void freeSet(set *S)
 459 {
 460   if (NULL != S)
 461     freeSet(S->secondNode);
 462   free(S);
 463 }
 464
 465 void freeTree(tree *T)
 466 {
 467   node *v;
 468   v = T->root;
 469   if (NULL != v->leftEdge)
 470     freeSubTree(v->leftEdge);
 471   free(T->root);
 472   free(T);
 473 }
 474
 475 void freeSubTree(edge *e)
 476 {
 477   node *v;
 478   edge *e1, *e2;
 479   v = e->head;
 480   e1 = v->leftEdge;
 481   if (NULL != e1)
 482     freeSubTree(e1);
 483   e2 = v->rightEdge;
 484   if (NULL != e2)
 485     freeSubTree(e2);
 486   free(v);
 487   e->tail = NULL;
 488   e->head = NULL;
 489   free(e);
 490 }