src/me.c

   1 /* me.c    2012-02-09 */
   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       error("root %s is poorly rooted.", v->label);
 278     }
 279   for(e = v->middleEdge, v->middleEdge = NULL; NULL != e; e = f )
 280     {
 281       w = e->head;
 282       v = e->tail;
 283       e->tail = w;
 284       e->head = v;
 285       f = w->leftEdge;
 286       v->parentEdge = e;
 287       w->leftEdge = e;
 288       w->parentEdge = NULL;
 289     }
 290   T->root = w;
 291   return(T);
 292 }
 293
 294 void compareSets(tree *T, set *S)
 295 {
 296   edge *e;
 297   node *v,*w;
 298   set *X;
 299   e = depthFirstTraverse(T,NULL);
 300   while (NULL != e)
 301     {
 302       v = e->head;
 303       for(X = S; NULL != X; X = X->secondNode)
 304         {
 305           w = X->firstNode;
 306           if (0 == strcmp(v->label,w->label))
 307             {
 308               v->index2 = w->index2;
 309             w->index2 = -1;
 310             break;
 311             }
 312         }
 313       e = depthFirstTraverse(T,e);
 314     }
 315   v = T->root;
 316   for(X = S; NULL != X; X = X->secondNode)
 317     {
 318       w = X->firstNode;
 319       if (0 == strcmp(v->label,w->label))
 320         {
 321           v->index2 = w->index2;
 322           w->index2 = -1;
 323           break;
 324         }
 325     }
 326   if (-1 == v->index2)
 327     {
 328       error("leaf %s in tree not in distance matrix.", v->label);
 329     }
 330   e = depthFirstTraverse(T,NULL);
 331   while (NULL != e)
 332     {
 333       v = e->head;
 334       if ((leaf(v)) && (-1 == v->index2))
 335         {
 336           error("leaf %s in tree not in distance matrix.", v->label);
 337         }
 338       e = depthFirstTraverse(T,e);
 339       }
 340   for(X = S; NULL != X; X = X->secondNode)
 341     if (X->firstNode->index2 > -1)
 342       {
 343         error("node %s in matrix but not a leaf in tree.", X->firstNode->label);
 344       }
 345   return;
 346 }
 347
 348 void partitionSizes(tree *T)
 349 {
 350   edge *e;
 351   e = depthFirstTraverse(T,NULL);
 352   while (NULL != e)
 353     {
 354       if (leaf(e->head))
 355         e->bottomsize = 1;
 356       else
 357         e->bottomsize = e->head->leftEdge->bottomsize
 358           + e->head->rightEdge->bottomsize;
 359       e->topsize = (T->size + 2)/2 - e->bottomsize;
 360       e = depthFirstTraverse(T,e);
 361     }
 362 }
 363
 364 /*************************************************************************
 365
 366                            TRAVERSE FUNCTIONS
 367
 368 *************************************************************************/
 369
 370 edge *depthFirstTraverse(tree *T, edge *e)
 371      /*depthFirstTraverse returns the edge f which is least in T according
 372        to the depth-first order, but which is later than e in the search
 373        pattern.  If e is null, f is the least edge of T*/
 374 {
 375   edge *f;
 376   if (NULL == e)
 377     {
 378       f = T->root->leftEdge;
 379       if (NULL != f)
 380         f = findBottomLeft(f);
 381       return(f);  /*this is the first edge of this search pattern*/
 382     }
 383   else /*e is non-null*/
 384     {
 385       if (e->tail->leftEdge == e)
 386         /*if e is a left-oriented edge, we skip the entire
 387           tree cut below e, and find least edge*/
 388         f = moveRight(e);
 389       else  /*if e is a right-oriented edge, we have already looked at its
 390               sibling and everything below e, so we move up*/
 391         f = e->tail->parentEdge;
 392     }
 393   return(f);
 394 }
 395
 396 edge *findBottomLeft(edge *e)
 397      /*findBottomLeft searches by gottom down in the tree and to the left.*/
 398 {
 399   edge *f;
 400   f = e;
 401   while (NULL != f->head->leftEdge)
 402     f = f->head->leftEdge;
 403   return(f);
 404 }
 405
 406 edge *moveRight(edge *e)
 407 {
 408   edge *f;
 409   f = e->tail->rightEdge; /*this step moves from a left-oriented edge
 410                             to a right-oriented edge*/
 411   if (NULL != f)
 412     f = findBottomLeft(f);
 413   return(f);
 414 }
 415
 416 edge *topFirstTraverse(tree *T, edge *e)
 417      /*topFirstTraverse starts from the top of T, and from there moves stepwise
 418        down, left before right*/
 419      /*assumes tree has been detrifurcated*/
 420 {
 421   edge *f;
 422   if (NULL == e)
 423     return(T->root->leftEdge); /*first Edge searched*/
 424   else if (!(leaf(e->head)))
 425     return(e->head->leftEdge); /*down and to the left is preferred*/
 426   else /*e->head is a leaf*/
 427     {
 428       f = moveUpRight(e);
 429       return(f);
 430     }
 431 }
 432
 433 edge *moveUpRight(edge *e)
 434 {
 435   edge *f;
 436   f = e;
 437   while ((NULL != f) && ( f->tail->leftEdge != f))
 438     f = f->tail->parentEdge;
 439   /*go up the tree until f is a leftEdge*/
 440   if (NULL == f)
 441     return(f); /*triggered at end of search*/
 442   else
 443     return(f->tail->rightEdge);
 444   /*and then go right*/
 445 }
 446
 447 /*************************************************************************
 448
 449                            FREE FUNCTIONS
 450
 451 *************************************************************************/
 452
 453 void freeMatrix(double **D, int size)
 454 {
 455   int i;
 456   for(i=0;i<size;i++)
 457     free(D[i]);
 458   free(D);
 459 }
 460
 461 void freeSet(set *S)
 462 {
 463   if (NULL != S)
 464     freeSet(S->secondNode);
 465   free(S);
 466 }
 467
 468 void freeTree(tree *T)
 469 {
 470   node *v;
 471   v = T->root;
 472   if (NULL != v->leftEdge)
 473     freeSubTree(v->leftEdge);
 474   free(T->root);
 475   free(T);
 476 }
 477
 478 void freeSubTree(edge *e)
 479 {
 480   node *v;
 481   edge *e1, *e2;
 482   v = e->head;
 483   e1 = v->leftEdge;
 484   if (NULL != e1)
 485     freeSubTree(e1);
 486   e2 = v->rightEdge;
 487   if (NULL != e2)
 488     freeSubTree(e2);
 489   free(v);
 490   e->tail = NULL;
 491   e->head = NULL;
 492   free(e);
 493 }