src/me.c

   1 /* me.c    2012-04-30 */
   2
   3 /* Copyright 2007-2008 Olivier Gascuel, Rick Desper,
   4    R port by Vincent Lefort and 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, int *labels,
  30           int *nni, int *spr, int *tbr,
  31           int *edge1, int *edge2, double *el)
  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, labels, 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, int *labels, int *nni,
  70           int *edge1, int *edge2, double *el)
  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, labels, 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 double **loadMatrix (double *X, int *labels, int n, set *S)
 130 {
 131 //  char nextString[MAX_LABEL_LENGTH];
 132   node *v;
 133   double **table;
 134   int i, j, a, b;
 135
 136   table = (double **) calloc(n,sizeof(double *));
 137   for(i=0; i<n; i++)
 138     table[i] = (double *) calloc(n,sizeof(double));
 139
 140   for(i=0; i<n; i++)
 141     {
 142 //      strncpy (nextString, labels[i], MAX_LABEL_LENGTH);
 143 //      ReplaceForbiddenChars (nextString, '_');
 144 //      v = makeNewNode(nextString,-1);
 145       v = makeNewNode(labels[i], -1);
 146       v->index2 = i;
 147       S = addToSet(v,S);
 148       for (j=i; j<n; j++) {
 149         a=i+1;
 150         b=j+1;
 151         table[j][i] = X[XINDEX(a,b)];
 152         table[i][j] = X[XINDEX(a,b)];
 153         if (i==j)
 154           table[i][j] = 0;
 155       }
 156     }
 157   return (table);
 158 }
 159
 160 /*
 161
 162   -- GRAPH FUNCTIONS --
 163
 164 */
 165
 166 set *addToSet(node *v, set *X)
 167 {
 168   if (NULL == X)
 169     {
 170       X = (set *) malloc(sizeof(set));
 171       X->firstNode = v;
 172       X->secondNode = NULL;
 173     }
 174   else if (NULL == X->firstNode)
 175     X->firstNode = v;
 176   else
 177     X->secondNode = addToSet(v,X->secondNode);
 178   return(X);
 179 }
 180
 181 //node *makeNewNode(char *label, int i)
 182 node *makeNewNode(int label, int i)
 183 {
 184   return(makeNode(label,NULL,i));
 185 }
 186
 187 //node *makeNode(char *label, edge *parentEdge, int index)
 188 node *makeNode(int label, edge *parentEdge, int index)
 189 {
 190   node *newNode;  /*points to new node added to the graph*/
 191   newNode = (node *) malloc(sizeof(node));
 192 //  strncpy(newNode->label,label,NODE_LABEL_LENGTH);
 193   newNode->label = label;
 194   newNode->index = index;
 195   newNode->index2 = -1;
 196   newNode->parentEdge = parentEdge;
 197   newNode->leftEdge = NULL;
 198   newNode->middleEdge = NULL;
 199   newNode->rightEdge = NULL;
 200   /*all fields have been initialized*/
 201   return(newNode);
 202 }
 203
 204 /*copyNode returns a copy of v which has all of the fields identical to those
 205 of v, except the node pointer fields*/
 206 node *copyNode(node *v)
 207 {
 208   node *w;
 209   w = makeNode(v->label,NULL,v->index);
 210   w->index2 = v->index2;
 211   return(w);
 212 }
 213
 214 edge *siblingEdge(edge *e)
 215 {
 216   if(e == e->tail->leftEdge)
 217     return(e->tail->rightEdge);
 218   else
 219     return(e->tail->leftEdge);
 220 }
 221
 222 edge *makeEdge(char *label, node *tail, node *head, double weight)
 223 {
 224   edge *newEdge;
 225   newEdge = (edge *) malloc(sizeof(edge));
 226   strncpy(newEdge->label,label,EDGE_LABEL_LENGTH);
 227   newEdge->tail = tail;
 228   newEdge->head = head;
 229   newEdge->distance = weight;
 230   newEdge->totalweight = 0.0;
 231   return(newEdge);
 232 }
 233
 234 tree *newTree()
 235 {
 236   tree *T;
 237   T = (tree *) malloc(sizeof(tree));
 238   T->root = NULL;
 239   T->size = 0;
 240   T->weight = -1;
 241   return(T);
 242 }
 243
 244 void updateSizes(edge *e, int direction)
 245 {
 246   edge *f;
 247   switch(direction)
 248     {
 249     case UP:
 250       f = e->head->leftEdge;
 251       if (NULL != f)
 252         updateSizes(f,UP);
 253       f = e->head->rightEdge;
 254       if (NULL != f)
 255         updateSizes(f,UP);
 256       e->topsize++;
 257       break;
 258     case DOWN:
 259       f = siblingEdge(e);
 260       if (NULL != f)
 261         updateSizes(f,UP);
 262       f = e->tail->parentEdge;
 263       if (NULL != f)
 264         updateSizes(f,DOWN);
 265       e->bottomsize++;
 266       break;
 267     }
 268 }
 269
 270 /*detrifurcate takes the (possibly trifurcated) input tree
 271   and reroots the tree to a leaf*/
 272 /*assumes tree is only trifurcated at root*/
 273 tree *detrifurcate(tree *T)
 274 {
 275   node *v, *w;
 276   edge *e, *f;
 277   v = T->root;
 278   if(leaf(v))
 279     return(T);
 280   if (NULL != v->parentEdge)
 281     {
 282       error("root %d is poorly rooted.", v->label);
 283     }
 284   for(e = v->middleEdge, v->middleEdge = NULL; NULL != e; e = f )
 285     {
 286       w = e->head;
 287       v = e->tail;
 288       e->tail = w;
 289       e->head = v;
 290       f = w->leftEdge;
 291       v->parentEdge = e;
 292       w->leftEdge = e;
 293       w->parentEdge = NULL;
 294     }
 295   T->root = w;
 296   return(T);
 297 }
 298
 299 void compareSets(tree *T, set *S)
 300 {
 301   edge *e;
 302   node *v,*w;
 303   set *X;
 304   e = depthFirstTraverse(T,NULL);
 305   while (NULL != e)
 306     {
 307       v = e->head;
 308       for(X = S; NULL != X; X = X->secondNode)
 309         {
 310           w = X->firstNode;
 311 //        if (0 == strcmp(v->label,w->label))
 312           if (v->label == w->label)
 313             {
 314               v->index2 = w->index2;
 315             w->index2 = -1;
 316             break;
 317             }
 318         }
 319       e = depthFirstTraverse(T,e);
 320     }
 321   v = T->root;
 322   for(X = S; NULL != X; X = X->secondNode)
 323     {
 324       w = X->firstNode;
 325 //      if (0 == strcmp(v->label,w->label))
 326       if (v->label == w->label)
 327         {
 328           v->index2 = w->index2;
 329           w->index2 = -1;
 330           break;
 331         }
 332     }
 333   if (-1 == v->index2)
 334     {
 335       error("leaf %d in tree not in distance matrix.", v->label);
 336     }
 337   e = depthFirstTraverse(T,NULL);
 338   while (NULL != e)
 339     {
 340       v = e->head;
 341       if ((leaf(v)) && (-1 == v->index2))
 342         {
 343           error("leaf %d in tree not in distance matrix.", v->label);
 344         }
 345       e = depthFirstTraverse(T,e);
 346       }
 347   for(X = S; NULL != X; X = X->secondNode)
 348     if (X->firstNode->index2 > -1)
 349       {
 350         error("node %d in matrix but not a leaf in tree.", X->firstNode->label);
 351       }
 352   return;
 353 }
 354
 355 void partitionSizes(tree *T)
 356 {
 357   edge *e;
 358   e = depthFirstTraverse(T,NULL);
 359   while (NULL != e)
 360     {
 361       if (leaf(e->head))
 362         e->bottomsize = 1;
 363       else
 364         e->bottomsize = e->head->leftEdge->bottomsize
 365           + e->head->rightEdge->bottomsize;
 366       e->topsize = (T->size + 2)/2 - e->bottomsize;
 367       e = depthFirstTraverse(T,e);
 368     }
 369 }
 370
 371 /*************************************************************************
 372
 373                            TRAVERSE FUNCTIONS
 374
 375 *************************************************************************/
 376
 377 edge *depthFirstTraverse(tree *T, edge *e)
 378      /*depthFirstTraverse returns the edge f which is least in T according
 379        to the depth-first order, but which is later than e in the search
 380        pattern.  If e is null, f is the least edge of T*/
 381 {
 382   edge *f;
 383   if (NULL == e)
 384     {
 385       f = T->root->leftEdge;
 386       if (NULL != f)
 387         f = findBottomLeft(f);
 388       return(f);  /*this is the first edge of this search pattern*/
 389     }
 390   else /*e is non-null*/
 391     {
 392       if (e->tail->leftEdge == e)
 393         /*if e is a left-oriented edge, we skip the entire
 394           tree cut below e, and find least edge*/
 395         f = moveRight(e);
 396       else  /*if e is a right-oriented edge, we have already looked at its
 397               sibling and everything below e, so we move up*/
 398         f = e->tail->parentEdge;
 399     }
 400   return(f);
 401 }
 402
 403 edge *findBottomLeft(edge *e)
 404      /*findBottomLeft searches by gottom down in the tree and to the left.*/
 405 {
 406   edge *f;
 407   f = e;
 408   while (NULL != f->head->leftEdge)
 409     f = f->head->leftEdge;
 410   return(f);
 411 }
 412
 413 edge *moveRight(edge *e)
 414 {
 415   edge *f;
 416   f = e->tail->rightEdge; /*this step moves from a left-oriented edge
 417                             to a right-oriented edge*/
 418   if (NULL != f)
 419     f = findBottomLeft(f);
 420   return(f);
 421 }
 422
 423 edge *topFirstTraverse(tree *T, edge *e)
 424      /*topFirstTraverse starts from the top of T, and from there moves stepwise
 425        down, left before right*/
 426      /*assumes tree has been detrifurcated*/
 427 {
 428   edge *f;
 429   if (NULL == e)
 430     return(T->root->leftEdge); /*first Edge searched*/
 431   else if (!(leaf(e->head)))
 432     return(e->head->leftEdge); /*down and to the left is preferred*/
 433   else /*e->head is a leaf*/
 434     {
 435       f = moveUpRight(e);
 436       return(f);
 437     }
 438 }
 439
 440 edge *moveUpRight(edge *e)
 441 {
 442   edge *f;
 443   f = e;
 444   while ((NULL != f) && ( f->tail->leftEdge != f))
 445     f = f->tail->parentEdge;
 446   /*go up the tree until f is a leftEdge*/
 447   if (NULL == f)
 448     return(f); /*triggered at end of search*/
 449   else
 450     return(f->tail->rightEdge);
 451   /*and then go right*/
 452 }
 453
 454 /*************************************************************************
 455
 456                            FREE FUNCTIONS
 457
 458 *************************************************************************/
 459
 460 void freeMatrix(double **D, int size)
 461 {
 462   int i;
 463   for(i=0;i<size;i++)
 464     free(D[i]);
 465   free(D);
 466 }
 467
 468 void freeSet(set *S)
 469 {
 470   if (NULL != S)
 471     freeSet(S->secondNode);
 472   free(S);
 473 }
 474
 475 void freeTree(tree *T)
 476 {
 477   node *v;
 478   v = T->root;
 479   if (NULL != v->leftEdge)
 480     freeSubTree(v->leftEdge);
 481   free(T->root);
 482   free(T);
 483 }
 484
 485 void freeSubTree(edge *e)
 486 {
 487   node *v;
 488   edge *e1, *e2;
 489   v = e->head;
 490   e1 = v->leftEdge;
 491   if (NULL != e1)
 492     freeSubTree(e1);
 493   e2 = v->rightEdge;
 494   if (NULL != e2)
 495     freeSubTree(e2);
 496   free(v);
 497   e->tail = NULL;
 498   e->head = NULL;
 499   free(e);
 500 }