src/me.c

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