src/newick.c

   1 /*#include <stdio.h>
   2 #include <stdlib.h>
   3 #include <math.h>
   4 #include <string.h>
   5 #include "graph.h"
   6 #include "newick.h"
   7 #include "main.h"
   8 */
   9 #include "me.h"
  10
  11 int nodeCount;
  12 int edgeCount;
  13
  14 int whiteSpace(char c)
  15 {
  16   if ((' ' == c) || ('\t' == c) || ('\n' == c))
  17     return(1);
  18   else return(0);
  19 }
  20
  21 int leaf(node *v)
  22 {
  23   int count = 0;
  24   if (NULL != v->parentEdge)
  25     count++;
  26   if (NULL != v->leftEdge)
  27     count++;
  28   if (NULL != v->rightEdge)
  29     count++;
  30   if (NULL != v->middleEdge)
  31     count++;
  32   if (count > 1)
  33     return(0);
  34   return(1);
  35 }
  36
  37 /*decodeNewickSubtree is used to turn a string of the form
  38   "(v1:d1,v2:d2,(subtree) v3:d3....vk:dk) subroot:d," into a subtree
  39   rooted at subrooted, with corresponding subtrees and leaves at v1
  40   through vk.  It is called recursively on subtrees*/
  41
  42 node *decodeNewickSubtree(char *treeString, tree *T, int *uCount)
  43 {
  44   node *thisNode;
  45   node *centerNode;
  46   double thisWeight;
  47   edge *thisEdge;
  48 //  char label[MAX_LABEL_LENGTH];
  49   char stringWeight[MAX_LABEL_LENGTH];
  50   int state;
  51   int i = 0;
  52   int j;
  53   int left,right;
  54   int parcount;
  55 //  snprintf(label,14,"Default_Label\0");
  56   left = right = 0;
  57   parcount = 0;
  58   state = ReadOpenParenthesis;
  59   if('(' == treeString[0])
  60     parcount++;
  61   //centerNode = makeNode(label,NULL,nodeCount++);
  62   centerNode = makeNode("",NULL,nodeCount++);
  63   T->size++;
  64   while(parcount > 0)
  65     {
  66       while(whiteSpace(treeString[i]))
  67         i++;
  68       switch(state)
  69         {
  70         case(ReadOpenParenthesis):
  71           if('(' != treeString[0])
  72             {
  73               Rprintf("Error reading subtree.\n");
  74               exit(0);
  75             }
  76           i++;
  77           state = ReadSubTree;
  78           break;
  79         case(ReadSubTree):
  80           if('(' == treeString[i])  /*if treeString[i] is a left parenthesis,
  81                                       we scan down the string until we find its partner.
  82                                       the relative number of '('s vs. ')'s is counted
  83                                       by the variable parcount*/
  84             {
  85               left = i++;
  86               parcount++;
  87               while(parcount > 1)
  88                 {
  89                   while (('(' != treeString[i]) && (')' != treeString[i]))
  90                     i++;  /*skip over characters which are not parentheses*/
  91                   if('(' == treeString[i])
  92                     parcount++;
  93                   else if (')' == treeString[i])
  94                     parcount--;
  95                   i++;
  96                 }  /*end while */
  97               right = i;  /*at this point, the subtree string goes from
  98                             treeString[left] to treeString[right - 1]*/
  99               thisNode = decodeNewickSubtree(treeString + left,T,uCount);  /*note that this
 100                                                                       step will put
 101                                                                       thisNode in T*/
 102               i = right;  /*having created the node for the subtree, we move
 103                             to assigning the label for the new node.
 104                             treeString[right] will be the start of this label */
 105             } /* end if ('(' == treeString[i]) condition */
 106           else
 107             {
 108               //thisNode = makeNode(label,NULL,nodeCount++);
 109               thisNode = makeNode("",NULL,nodeCount++);
 110               T->size++;
 111             }
 112           state = ReadLabel;
 113           break;
 114         case(ReadLabel):
 115           left = i;  /*recall "left" is the left marker for the substring, "right" the right*/
 116           if (':' == treeString[i])   /*this means an internal node?*/
 117             {
 118               //sprintf(thisNode->label,"I%d",(*uCount)++);
 119               //snprintf(thisNode->label,MAX_LABEL_LENGTH,"I%d",(*uCount)++);
 120               (*uCount)++;
 121               right = i;
 122             }
 123           else
 124             {
 125               while((':' != treeString[i]) && (',' != treeString[i]) && (')' != treeString[i]))
 126                 i++;
 127               right = i;
 128               j = 0;
 129               for(i = left; i < right;i++)
 130                 if(!(whiteSpace(treeString[i])))
 131                   thisNode->label[j++] = treeString[i];
 132               thisNode->label[j] = '\0';
 133             }
 134           if(':' == treeString[right])
 135             state = ReadWeight;
 136           else
 137             {
 138               state = AddEdge;
 139               thisWeight = 0.0;
 140             }
 141           i = right + 1;
 142           break;
 143         case(ReadWeight):
 144           left = i;
 145           while
 146             ((',' != treeString[i]) && (')' != treeString[i]))
 147             i++;
 148           right = i;
 149           j = 0;
 150           for(i = left; i < right; i++)
 151             stringWeight[j++] = treeString[i];
 152           stringWeight[j] = '\0';
 153           thisWeight = atof(stringWeight);
 154           state=AddEdge;
 155           break;
 156         case(AddEdge):
 157           //thisEdge = makeEdge(label,centerNode,thisNode,thisWeight);
 158           thisEdge = makeEdge("",centerNode,thisNode,thisWeight);
 159           thisNode->parentEdge = thisEdge;
 160           if (NULL == centerNode->leftEdge)
 161             centerNode->leftEdge = thisEdge;
 162           else if (NULL == centerNode->rightEdge)
 163             centerNode->rightEdge = thisEdge;
 164           else if (NULL == centerNode->middleEdge)
 165             centerNode->middleEdge = thisEdge;
 166           else
 167             {
 168               Rprintf("Error: node %s has too many (>3) children.\n",centerNode->label);
 169               exit(0);
 170             }
 171           //sprintf(thisEdge->label,"E%d",edgeCount++);
 172           //snprintf(thisEdge->label,MAX_LABEL_LENGTH,"E%d",edgeCount++);
 173           edgeCount++;
 174           i = right + 1;
 175           if (',' == treeString[right])
 176             state = ReadSubTree;
 177           else
 178             parcount--;
 179           break;
 180         }
 181     }
 182   return(centerNode);
 183 }
 184
 185 tree *readNewickString (char *str, int numLeaves)
 186 {
 187   tree *T;
 188   node *centerNode;
 189   int i = 0;
 190   int j = 0;
 191   int inputLength;
 192   int uCount = 0;
 193   int parCount = 0;
 194   char rootLabel[MAX_LABEL_LENGTH];
 195   nodeCount = edgeCount = 0;
 196
 197   T = newTree();
 198
 199   if ('(' != str[0])
 200     {
 201       Rprintf("Error reading generated tree - does not start with '('.\n");
 202       exit(0);
 203     }
 204   inputLength = strlen (str)+1;
 205   for(i = 0; i < inputLength; i++)
 206     {
 207       if ('(' == str[i])
 208         parCount++;
 209       else if (')' == str[i])
 210         parCount--;
 211       if (parCount > 0)
 212         ;
 213       else if (0 == parCount)
 214         {
 215           i++;
 216 /*        if(';' == str[i])
 217             sprintf(rootLabel,"URoot");
 218           else
 219             {*/
 220               while(';' != str[i])
 221                 if (!(whiteSpace (str[i++])))
 222                   rootLabel[j++] = str[i-1];  /*be careful here */
 223                 rootLabel[j] = '\0';
 224 //          }
 225           i = inputLength;
 226         }
 227       else if (parCount < 0)
 228         {
 229           Rprintf("Error reading generated tree. Too many right parentheses.\n");
 230           exit(0);
 231         }
 232     }
 233   centerNode = decodeNewickSubtree (str, T, &uCount);
 234   snprintf (centerNode->label, MAX_LABEL_LENGTH, rootLabel);
 235   T->root = centerNode;
 236   return (T);
 237 }
 238
 239 tree *loadNewickTree(FILE *ifile, int numLeaves)
 240 {
 241 //  char label[] = "EmptyEdge";
 242   tree *T;
 243   node *centerNode;
 244   int i = 0;
 245   int j = 0;
 246   int inputLength;
 247   int uCount = 0;
 248   int parCount = 0;
 249   char c;
 250   int Comment;
 251   char *nextString;
 252   char rootLabel[MAX_LABEL_LENGTH];
 253   nodeCount = edgeCount = 0;
 254   T = newTree();
 255   nextString = (char *) malloc(numLeaves*INPUT_SIZE*sizeof(char));
 256   if (NULL == nextString)
 257     nextString = (char *) malloc(MAX_INPUT_SIZE*sizeof(char));
 258   Comment = 0;
 259   while(1 == fscanf(ifile,"%c",&c))
 260     {
 261       if('[' == c)
 262         Comment = 1;
 263       else if (']' == c)
 264         Comment = 0;
 265       else if (!(Comment))
 266         {
 267           if(whiteSpace(c))
 268             {
 269               if (i > 0)
 270                 nextString[i++] = ' ';
 271             }
 272           else  /*note that this else goes with if(whiteSpace(c))*/
 273             nextString[i++] = c;
 274           if (';' == c)
 275             break;
 276         }
 277     }
 278   if ('(' != nextString[0])
 279     {
 280       fprintf(stderr,"Error reading input file - does not start with '('.\n");
 281       exit(EXIT_FAILURE);
 282     }
 283   inputLength = i;
 284   for(i = 0; i < inputLength;i++)
 285     {
 286       if ('(' == nextString[i])
 287         parCount++;
 288       else if (')' == nextString[i])
 289         parCount--;
 290       if (parCount > 0)
 291         ;
 292       else if (0 == parCount)
 293         {
 294           i++;
 295 /*        if(';' == nextString[i])
 296             sprintf(rootLabel,"URoot");
 297           else
 298             {*/
 299               while(';' != nextString[i])
 300                 if(!(whiteSpace(nextString[i++])))
 301                   rootLabel[j++] = nextString[i-1];  /*be careful here */
 302               rootLabel[j] = '\0';
 303 //          }
 304           i = inputLength;
 305         }
 306       else if (parCount < 0)
 307         {
 308           fprintf(stderr,"Error reading tree input file.  Too many right parentheses.\n");
 309           exit(EXIT_FAILURE);
 310         }
 311     }
 312   centerNode = decodeNewickSubtree(nextString,T,&uCount);
 313   snprintf(centerNode->label, MAX_LABEL_LENGTH, rootLabel);
 314   T->root = centerNode;
 315   free(nextString);
 316   return(T);
 317 }
 318
 319 double GetSubTreeLength (tree *T, edge *e)
 320 {
 321   double ret = 0;
 322
 323   if ( (NULL != e) && (! leaf(e->head) )) {
 324     ret += GetSubTreeLength (T, e->head->leftEdge);
 325     ret += GetSubTreeLength (T, e->head->rightEdge);
 326   }
 327   ret += e->distance;
 328   return ret;
 329 }
 330
 331 void NewickPrintSubtree(tree *T, edge *e, char *str)
 332 {
 333   char *tmp;
 334   if (NULL == e)
 335     {
 336       Rprintf("Error with Newick Printing routine.\n");
 337       exit(0);
 338     }
 339   if(!(leaf(e->head)))
 340     {
 341       if (strlen (str) < MAX_INPUT_SIZE -2)
 342         strncat (str, "(", 1);
 343       NewickPrintSubtree(T,e->head->leftEdge,str);
 344       if (strlen (str) < MAX_INPUT_SIZE -2)
 345         strncat (str, ",", 1);
 346       NewickPrintSubtree(T,e->head->rightEdge,str);
 347       if (strlen (str) < MAX_INPUT_SIZE -2)
 348         strncat (str, ")", 1);
 349     }
 350   if (strlen (str) < MAX_INPUT_SIZE - strlen (e->head->label) -1)
 351     strncat (str, e->head->label, strlen (e->head->label));
 352
 353   if (strlen (str) < MAX_INPUT_SIZE - 2)
 354     strncat (str, ":", 1);
 355
 356   tmp = (char *)R_alloc(INPUT_SIZE, sizeof(char));
 357   /* added by EP */
 358   if (strlen(tmp))
 359     strncpy(tmp, "", strlen(tmp));
 360   /* end */
 361   snprintf (tmp, INPUT_SIZE, "%lf", e->distance);
 362   if (strlen (str) < MAX_INPUT_SIZE - strlen (tmp) -1)
 363     strncat (str, tmp, strlen (tmp));
 364
 365   /* free (tmp); */
 366   return;
 367 }
 368
 369 double GetBinaryTreeLength (tree *T)
 370 {
 371   double ret = 0;
 372   edge *e, *f;
 373   node *rootchild;
 374   e = T->root->leftEdge;
 375   rootchild = e->head;
 376
 377   f = rootchild->leftEdge;
 378   if (NULL != f)
 379     ret += GetSubTreeLength (T, f);
 380   f = rootchild->rightEdge;
 381   if (NULL != f)
 382     ret += GetSubTreeLength (T, f);
 383   ret += e->distance;
 384   return ret;
 385 }
 386
 387 void NewickPrintBinaryTree(tree *T, char *str)
 388 {
 389   edge *e, *f;
 390   node *rootchild;
 391   char *tmp;
 392   e = T->root->leftEdge;
 393   rootchild = e->head;
 394   if (strlen (str) < MAX_INPUT_SIZE -2)
 395     strncat (str, "(", 1);
 396   f = rootchild->leftEdge;
 397   if (NULL != f)
 398     {
 399       NewickPrintSubtree(T,f,str);
 400       if (strlen (str) < MAX_INPUT_SIZE -2)
 401         strncat (str, ",", 1);
 402     }
 403   f = rootchild->rightEdge;
 404   if (NULL != f)
 405     {
 406       NewickPrintSubtree(T,f,str);
 407       if (strlen (str) < MAX_INPUT_SIZE -2)
 408         strncat (str, ",", 1);
 409     }
 410   if (strlen (str) < MAX_INPUT_SIZE - strlen (T->root->label) -1)
 411     strncat (str, T->root->label, strlen (T->root->label));
 412
 413   if (strlen (str) < MAX_INPUT_SIZE - 2)
 414     strncat (str, ":", 1);
 415
 416   tmp = (char *)R_alloc(INPUT_SIZE, sizeof(char));
 417   /* added by EP */
 418   if (strlen(tmp))
 419     strncpy(tmp, "", strlen(tmp));
 420   /* end */
 421   snprintf (tmp, INPUT_SIZE, "%lf", e->distance);
 422   if (strlen (str) < MAX_INPUT_SIZE - strlen (tmp) -1)
 423     strncat (str, tmp, strlen (tmp));
 424
 425   if (strlen (str) < MAX_INPUT_SIZE - 2)
 426     strncat (str, ")", 1);
 427
 428   if (NULL != rootchild->label)
 429     if (strlen (str) < MAX_INPUT_SIZE - strlen (rootchild->label) -1)
 430       strncat (str, T->root->label, strlen (rootchild->label));
 431
 432   if (strlen (str) < MAX_INPUT_SIZE - 3)
 433     strncat (str, ";\n", 2);
 434
 435   /* free (tmp); */
 436   return;
 437 }
 438
 439 double GetTrinaryTreeLength (tree *T)
 440 {
 441   double ret = 0;
 442   edge *f;
 443   f = T->root->leftEdge;
 444   if (NULL != f)
 445     ret += GetSubTreeLength (T, f);
 446   f = T->root->rightEdge;
 447   if (NULL != f)
 448     ret += GetSubTreeLength (T, f);
 449   f = T->root->middleEdge;
 450   if (NULL != f)
 451     ret += GetSubTreeLength (T, f);
 452
 453   return ret;
 454 }
 455
 456 void NewickPrintTrinaryTree(tree *T, char *str)
 457 {
 458   edge *f;
 459   f = T->root->leftEdge;
 460   if (strlen (str) < MAX_INPUT_SIZE -2)
 461     strncat (str, "(", 1);
 462   if (NULL != f)
 463     {
 464       NewickPrintSubtree(T,f,str);
 465       if (strlen (str) < MAX_INPUT_SIZE -2)
 466         strncat (str, ",", 1);
 467     }
 468   f = T->root->rightEdge;
 469   if (NULL != f)
 470     {
 471       NewickPrintSubtree(T,f,str);
 472       if (strlen (str) < MAX_INPUT_SIZE -2)
 473         strncat (str, ",", 1);
 474     }
 475   f = T->root->middleEdge;
 476   if (NULL != f)
 477     {
 478       NewickPrintSubtree(T,f,str);
 479       if (strlen (str) < MAX_INPUT_SIZE -2)
 480         strncat (str, ")", 1);
 481     }
 482   if (NULL != T->root->label)
 483     if (strlen (str) < MAX_INPUT_SIZE - strlen (T->root->label) -1)
 484       strncat (str, T->root->label, strlen (T->root->label));
 485   if (strlen (str) < MAX_INPUT_SIZE - 3)
 486     strncat (str, ";\n", 2);
 487   return;
 488 }
 489
 490 void NewickPrintTreeStr(tree *T, char *str)
 491 {
 492   if (leaf(T->root))
 493     NewickPrintBinaryTree(T,str);
 494   else
 495     NewickPrintTrinaryTree(T,str);
 496 }
 497
 498 double GetTreeLength (tree *T)
 499 {
 500   double ret = 0;
 501   if (leaf(T->root))
 502     ret = GetBinaryTreeLength (T);
 503   else
 504     ret = GetTrinaryTreeLength (T);
 505   return ret;
 506 }
 507 /*
 508 void NewickPrintTree(tree *T, FILE *ofile)
 509 {
 510   if (leaf(T->root))
 511     NewickPrintBinaryTree(T,ofile);
 512   else
 513     NewickPrintTrinaryTree(T,ofile);
 514 }
 515 */
 516 //edge *depthFirstTraverse(tree *T, edge *e);
 517