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