forester/archive/RIO/others/phylip_mod/src/cons.c

   1 #include "phylip.h"
   2 #include "cons.h"
   3
   4 int tree_pairing;
   5
   6 Char outfilename[FNMLNGTH], intreename[FNMLNGTH], intree2name[FNMLNGTH], outtreename[FNMLNGTH];
   7 node *root;
   8
   9 long numopts, outgrno, col, setsz;
  10 long maxgrp;               /* max. no. of groups in all trees found  */
  11
  12 boolean trout, firsttree, noroot, outgropt, didreroot, prntsets,
  13                progress, treeprint, goteof, strict, mr=false, mre=false,
  14                ml=false; /* initialized all false for Treedist */
  15 pointarray nodep;
  16 pointarray treenode;
  17 group_type **grouping, **grping2, **group2;/* to store groups found  */
  18 double *lengths, *lengths2;
  19 long **order, **order2, lasti;
  20 group_type *fullset;
  21 node *grbg;
  22 long tipy;
  23
  24 double **timesseen, **tmseen2, **times2 ;
  25 double trweight, ntrees, mlfrac;
  26
  27 /* prototypes */
  28 void censor(void);
  29 boolean compatible(long, long);
  30 void elimboth(long);
  31 void enternohash(group_type*, long*);
  32 void enterpartition (group_type*, long*);
  33 void reorient(node* n);
  34
  35 /* begin hash table code */
  36
  37 #define NUM_BUCKETS 100
  38
  39 typedef struct namenode {
  40   struct namenode *next;
  41   plotstring naym;
  42   int hitCount;
  43 } namenode;
  44
  45 typedef namenode **hashtype;
  46
  47 hashtype hashp;
  48
  49 long namesGetBucket(plotstring);
  50 void namesAdd(plotstring);
  51 boolean namesSearch(plotstring);
  52 void namesDelete(plotstring);
  53 void namesClearTable(void);
  54 void namesCheckTable(void);
  55 void missingnameRecurs(node *p);
  56
  57 /**
  58  * namesGetBucket - return the bucket for a given name
  59  */
  60 long namesGetBucket(plotstring searchname) {
  61   long i;
  62   long sum = 0;
  63
  64   for (i = 0; (i < MAXNCH) && (searchname[i] != '\0'); i++) {
  65     sum += searchname[i];
  66   }
  67   return (sum % NUM_BUCKETS);
  68 }
  69
  70
  71 /**
  72  * namesAdd - add a name to the hash table
  73  *
  74  * The argument is added at the head of the appropriate linked list.  No
  75  * checking is done for duplicates.  The caller can call
  76  * namesSearch to check for an existing name prior to calling
  77  * namesAdd.
  78  */
  79 void namesAdd(plotstring addname) {
  80   long bucket = namesGetBucket(addname);
  81   namenode *hp, *temp;
  82
  83   temp = hashp[bucket];
  84   hashp[bucket] = (namenode *)Malloc(sizeof(namenode));
  85   hp = hashp[bucket];
  86   strcpy(hp->naym, addname);
  87   hp->next = temp;
  88   hp->hitCount = 0;
  89 }
  90
  91 /**
  92  * namesSearch - search for a name in the hash table
  93  *
  94  * Return true if the name is found, else false.
  95  */
  96 boolean namesSearch(plotstring searchname) {
  97   long i = namesGetBucket(searchname);
  98   namenode *p;
  99
 100   p = hashp[i];
 101   if (p == NULL) {
 102     return false;
 103   }
 104   do {
 105     if (strcmp(searchname,p->naym) == 0) {
 106       p->hitCount++;
 107       return true;
 108     }
 109     p = p->next;
 110   } while (p != NULL);
 111
 112   return false;
 113 }
 114
 115 /**
 116  * Go through hash table and check that the hit count on all entries is one.
 117  * If it is zero, then a species was missed, if it is two, then there is a
 118  * duplicate species.
 119  */
 120
 121 void namesCheckTable(void) {
 122   namenode *p;
 123   long i;
 124
 125   for (i=0; i< NUM_BUCKETS; i++) {
 126     p = hashp[i];
 127     while (p != NULL){
 128       if(p->hitCount >1){
 129         printf("\n\nERROR in user tree: duplicate name found: ");
 130         puts(p->naym);
 131         printf("\n\n");
 132         exxit(-1);
 133       } else if(p->hitCount == 0){
 134         printf("\n\nERROR in user tree: name %s not found\n\n\n",
 135                p->naym);
 136         exxit(-1);
 137       }
 138       p->hitCount = 0;
 139       p = p->next;
 140     }
 141   }
 142 }
 143
 144 /**
 145  * namesClearTable - empty names out of the table and
 146  *                   return allocated memory
 147  */
 148 void namesClearTable(void) {
 149   long i;
 150   namenode *p, *temp;
 151
 152   for (i=0; i< NUM_BUCKETS; i++) {
 153     p = hashp[i];
 154     if (p != NULL) {
 155       do {
 156         temp = p;
 157         p = p->next;
 158         free(temp);
 159       } while (p != NULL);
 160     hashp[i] = NULL;
 161     }
 162   }
 163 }
 164 /* end hash table code */
 165
 166 void initconsnode(node **p, node **grbg, node *q, long len, long nodei,
 167                         long *ntips, long *parens, initops whichinit,
 168                         pointarray treenode, pointarray nodep, Char *str,
 169                         Char *ch, FILE *intree)
 170 {
 171   /* initializes a node */
 172   long i;
 173   char c;
 174   boolean minusread;
 175   double valyew, divisor, fracchange;
 176
 177   switch (whichinit) {
 178   case bottom:
 179     gnu(grbg, p);
 180     (*p)->index = nodei;
 181     (*p)->tip = false;
 182     for (i=0; i<MAXNCH; i++)
 183       (*p)->nayme[i] = '\0';
 184     nodep[(*p)->index - 1] = (*p);
 185     (*p)->v = 0;
 186     break;
 187   case nonbottom:
 188     gnu(grbg, p);
 189     (*p)->index = nodei;
 190     (*p)->v = 0;
 191     break;
 192   case tip:
 193     (*ntips)++;
 194     gnu(grbg, p);
 195     nodep[(*ntips) - 1] = *p;
 196     setupnode(*p, *ntips);
 197     (*p)->tip = true;
 198     strncpy ((*p)->nayme, str, MAXNCH);
 199     if (firsttree && prntsets) {
 200       fprintf(outfile, "  %ld. ", *ntips);
 201       for (i = 0; i < len; i++)
 202         putc(str[i], outfile);
 203       putc('\n', outfile);
 204       if ((*ntips > 0) && (((*ntips) % 10) == 0))
 205         putc('\n', outfile);
 206     }
 207     (*p)->v = 0;
 208     break;
 209   case length:
 210     processlength(&valyew, &divisor, ch, &minusread, intree, parens);
 211     fracchange = 1.0;
 212     (*p)->v = valyew / divisor / fracchange;
 213     break;
 214   case treewt:
 215     if (!eoln(intree)) {
 216       fscanf(intree, "%lf", &trweight);
 217       getch(ch, parens, intree);
 218       if (*ch != ']') {
 219         printf("\n\nERROR: Missing right square bracket\n\n");
 220         exxit(-1);
 221       } else {
 222         getch(ch, parens, intree);
 223         if (*ch != ';') {
 224           printf("\n\nERROR: Missing semicolon after square brackets\n\n");
 225           exxit(-1);
 226         }
 227       }
 228     }
 229     break;
 230   case unittrwt:
 231     /* This comes not only when setting trweight but also at the end of
 232      * any tree. The following code saves the current position in a
 233      * file and reads to a new line. If there is a new line then we're
 234      * at the end of tree, otherwise warn the user. This function should
 235      * really leave the file alone, so once we're done with 'intree'
 236      * we seek the position back so that it doesn't look like we did
 237      * anything */
 238     trweight = 1.0 ;
 239     i = ftell (intree);
 240     c = ' ';
 241     while (c == ' ')  {
 242       if (eoff(intree)) {
 243         fseek(intree,i,SEEK_SET);
 244         return;
 245       }
 246       c = gettc(intree);
 247     }
 248     fseek(intree,i,SEEK_SET);
 249     if ( c != '\n' && c!= '\r')
 250       printf("WARNING: Tree weight set to 1.0\n");
 251     if ( c == '\r' )
 252       if ( (c == gettc(intree)) != '\n')
 253         ungetc(c, intree);
 254     break;
 255   case hsnolength:
 256     (*p)->v = -1;         /* signal value that a length is missing */
 257     break;
 258   default:                /* cases hslength, iter, hsnolength      */
 259     break;                /* should there be an error message here?*/
 260   }
 261 } /* initconsnode */
 262
 263
 264 void censor(void)
 265 {
 266   /* delete groups that are too rare to be in the consensus tree */
 267   long i;
 268
 269   i = 1;
 270   do {
 271     if (timesseen[i-1])
 272       if (!(mre || (mr && (2*(*timesseen[i-1]) > ntrees))
 273                 || (ml && ((*timesseen[i-1]) > mlfrac*ntrees))
 274                 || (strict && ((*timesseen[i-1]) == ntrees)))) {
 275         free(grouping[i - 1]);
 276         free(timesseen[i - 1]);
 277         grouping[i - 1] = NULL;
 278         timesseen[i - 1] = NULL;
 279     }
 280     i++;
 281   } while (i < maxgrp);
 282 } /* censor */
 283
 284
 285 void compress(long *n)
 286 {
 287   /* push all the nonempty subsets to the front end of their array */
 288   long i, j;
 289
 290   i = 1;
 291   j = 1;
 292   do {
 293     while (grouping[i - 1] != NULL)
 294       i++;
 295     if (j <= i)
 296       j = i + 1;
 297     while ((grouping[j - 1] == NULL) && (j < maxgrp))
 298       j++;
 299     if (j < maxgrp) {
 300       grouping[i - 1] = (group_type *)Malloc(setsz * sizeof(group_type));
 301       timesseen[i - 1] = (double *)Malloc(sizeof(double));
 302       memcpy(grouping[i - 1], grouping[j - 1], setsz * sizeof(group_type));
 303       *timesseen[i - 1] = *timesseen[j - 1];
 304       free(grouping[j - 1]);
 305       free(timesseen[j - 1]);
 306       grouping[j - 1] = NULL;
 307       timesseen[j - 1] = NULL;
 308     }
 309   } while (j != maxgrp);
 310   (*n) = i - 1;
 311 }  /* compress */
 312
 313
 314 void sort(long n)
 315 {
 316   /* Shell sort keeping grouping, timesseen in same order */
 317   long gap, i, j;
 318   group_type *stemp;
 319   double rtemp;
 320
 321   gap = n / 2;
 322   stemp = (group_type *)Malloc(setsz * sizeof(group_type));
 323   while (gap > 0) {
 324     for (i = gap + 1; i <= n; i++) {
 325       j = i - gap;
 326       while (j > 0) {
 327         if (*timesseen[j - 1] < *timesseen[j + gap - 1]) {
 328           memcpy(stemp, grouping[j - 1], setsz * sizeof(group_type));
 329           memcpy(grouping[j - 1], grouping[j + gap - 1], setsz * sizeof(group_type));
 330           memcpy(grouping[j + gap - 1], stemp, setsz * sizeof(group_type));
 331           rtemp = *timesseen[j - 1];
 332           *timesseen[j - 1] = *timesseen[j + gap - 1];
 333           *timesseen[j + gap - 1] = rtemp;
 334         }
 335         j -= gap;
 336       }
 337     }
 338     gap /= 2;
 339   }
 340   free(stemp);
 341 }  /* sort */
 342
 343
 344 boolean compatible(long i, long j)
 345 {
 346   /* are groups i and j compatible? */
 347   boolean comp;
 348   long k;
 349
 350   comp = true;
 351   for (k = 0; k < setsz; k++)
 352     if ((grouping[i][k] & grouping[j][k]) != 0)
 353       comp = false;
 354   if (!comp) {
 355     comp = true;
 356     for (k = 0; k < setsz; k++)
 357       if ((grouping[i][k] & ~grouping[j][k]) != 0)
 358         comp = false;
 359     if (!comp) {
 360       comp = true;
 361       for (k = 0; k < setsz; k++)
 362         if ((grouping[j][k] & ~grouping[i][k]) != 0)
 363           comp = false;
 364       if (!comp) {
 365         comp = noroot;
 366         if (comp) {
 367           for (k = 0; k < setsz; k++)
 368             if ((fullset[k] & ~grouping[i][k] & ~grouping[j][k]) != 0)
 369               comp = false;
 370         }
 371       }
 372     }
 373   }
 374   return comp;
 375 } /* compatible */
 376
 377
 378 void eliminate(long *n, long *n2)
 379 {
 380   /* eliminate groups incompatible with preceding ones */
 381   long i, j, k;
 382   boolean comp;
 383
 384   for (i = 2; i <= (*n); i++) {
 385     comp = true;
 386     for (j = 0; comp && (j <= i - 2); j++) {
 387       if ((timesseen[j] != NULL) && *timesseen[j] > 0) {
 388         comp = compatible(i-1,j);
 389         if (!comp) {
 390           (*n2)++;
 391           times2[(*n2) - 1] = (double *)Malloc(sizeof(double));
 392           group2[(*n2) - 1] = (group_type *)Malloc(setsz * sizeof(group_type));
 393           *times2[(*n2) - 1] = *timesseen[i - 1];
 394           memcpy(group2[(*n2) - 1], grouping[i - 1], setsz * sizeof(group_type));
 395           *timesseen[i - 1] = 0.0;
 396           for (k = 0; k < setsz; k++)
 397             grouping[i - 1][k] = 0;
 398         }
 399       }
 400     }
 401     if (*timesseen[i - 1] == 0.0) {
 402       free(grouping[i - 1]);
 403       free(timesseen[i -  1]);
 404       timesseen[i - 1] = NULL;
 405       grouping[i - 1] = NULL;
 406     }
 407   }
 408 }  /* eliminate */
 409
 410
 411 void printset(long n)
 412 {
 413   /* print out the n sets of species */
 414   long i, j, k, size;
 415   boolean noneprinted;
 416
 417   fprintf(outfile, "\nSet (species in order)   ");
 418   for (i = 1; i <= spp - 25; i++)
 419     putc(' ', outfile);
 420   fprintf(outfile, "  How many times out of %7.2f\n\n", ntrees);
 421   noneprinted = true;
 422   for (i = 0; i < n; i++) {
 423     if ((timesseen[i] != NULL) && (*timesseen[i] > 0)) {
 424       size = 0;
 425       k = 0;
 426       for (j = 1; j <= spp; j++) {
 427         if (j == ((k+1)*SETBITS+1)) k++;
 428         if (((1L << (j - 1 - k*SETBITS)) & grouping[i][k]) != 0)
 429           size++;
 430       }
 431       if (size != 1 && !(noroot && size >= (spp-1))) {
 432         noneprinted = false;
 433         k = 0;
 434         for (j = 1; j <= spp; j++) {
 435           if (j == ((k+1)*SETBITS+1)) k++;
 436           if (((1L << (j - 1 - k*SETBITS)) & grouping[i][k]) != 0)
 437             putc('*', outfile);
 438           else
 439             putc('.', outfile);
 440           if (j % 10 == 0)
 441             putc(' ', outfile);
 442         }
 443         for (j = 1; j <= 23 - spp; j++)
 444           putc(' ', outfile);
 445         fprintf(outfile, "    %5.2f\n", *timesseen[i]);
 446       }
 447     }
 448   }
 449   if (noneprinted)
 450     fprintf(outfile, " NONE\n");
 451 }  /* printset */
 452
 453
 454 void bigsubset(group_type *st, long n)
 455 {
 456   /* Find a maximal subset of st among the n groupings,
 457      to be the set at the base of the tree.  */
 458   long i, j;
 459   group_type *su;
 460   boolean max, same;
 461
 462   su = (group_type *)Malloc(setsz * sizeof(group_type));
 463   for (i = 0; i < setsz; i++)
 464     su[i] = 0;
 465   for (i = 0; i < n; i++) {
 466     max = true;
 467     for (j = 0; j < setsz; j++)
 468       if ((grouping[i][j] & ~st[j]) != 0)
 469         max = false;
 470     if (max) {
 471       same = true;
 472       for (j = 0; j < setsz; j++)
 473         if (grouping[i][j] != st[j])
 474           same = false;
 475       max = !same;
 476     }
 477     if (max) {
 478       for (j = 0; j < setsz; j ++)
 479         if ((su[j] & ~grouping[i][j]) != 0)
 480           max = false;
 481       if (max) {
 482         same = true;
 483         for (j = 0; j < setsz; j ++)
 484           if (su[j] != grouping[i][j])
 485             same = false;
 486         max = !same;
 487       }
 488       if (max)
 489         memcpy(su, grouping[i], setsz * sizeof(group_type));
 490     }
 491   }
 492   memcpy(st, su, setsz * sizeof(group_type));
 493   free(su);
 494 }  /* bigsubset */
 495
 496
 497 void recontraverse(node **p, group_type *st, long n, long *nextnode)
 498 {
 499   /* traverse to add next node to consensus tree */
 500   long i, j = 0, k = 0, l = 0;
 501
 502   boolean found, same = 0, zero, zero2;
 503   group_type *tempset, *st2;
 504   node *q, *r;
 505
 506   for (i = 1; i <= spp; i++) {  /* count species in set */
 507     if (i == ((l+1)*SETBITS+1)) l++;
 508     if (((1L << (i - 1 - l*SETBITS)) & st[l]) != 0) {
 509       k++;               /* k  is the number of species in the set */
 510       j = i;             /* j  is set to last species in the set */
 511     }
 512   }
 513   if (k == 1) {           /* if only 1, set up that tip */
 514     *p = nodep[j - 1];
 515     (*p)->tip = true;
 516     (*p)->index = j;
 517     return;
 518   }
 519   gnu(&grbg, p);          /* otherwise make interior node */
 520   (*p)->tip = false;
 521   (*p)->index = *nextnode;
 522   nodep[*nextnode - 1] = *p;
 523   (*nextnode)++;
 524   (*p)->deltav = 0.0;
 525   for (i = 0; i < n; i++) { /* go through all sets */
 526     same = true;            /* to find one which is this one */
 527     for (j = 0; j < setsz; j++)
 528       if (grouping[i][j] != st[j])
 529         same = false;
 530     if (same)
 531       (*p)->deltav = *timesseen[i];
 532   }
 533   tempset = (group_type *)Malloc(setsz * sizeof(group_type));
 534   memcpy(tempset, st, setsz * sizeof(group_type));
 535   q = *p;
 536   st2 = (group_type *)Malloc(setsz * sizeof(group_type));
 537   memcpy(st2, st, setsz * sizeof(group_type));
 538   zero = true;      /* having made two copies of the set ... */
 539   for (j = 0; j < setsz; j++)       /* see if they are empty */
 540     if (tempset[j] != 0)
 541       zero = false;
 542   if (!zero)
 543     bigsubset(tempset, n);        /* find biggest set within it */
 544   zero = zero2 = false;           /* ... tempset is that subset */
 545   while (!zero && !zero2) {
 546     zero = zero2 = true;
 547     for (j = 0; j < setsz; j++) {
 548       if (st2[j] != 0)
 549         zero = false;
 550       if (tempset[j] != 0)
 551         zero2 = false;
 552     }
 553     if (!zero && !zero2) {
 554       gnu(&grbg, &q->next);
 555       q->next->index = q->index;
 556       q = q->next;
 557       q->tip = false;
 558       r = *p;
 559       recontraverse(&q->back, tempset, n, nextnode); /* put it on tree */
 560       *p = r;
 561       q->back->back = q;
 562       for (j = 0; j < setsz; j++)
 563         st2[j] &= ~tempset[j];     /* remove that subset from the set */
 564       memcpy(tempset, st2, setsz * sizeof(group_type));  /* that becomes set */
 565       found = false;
 566       i = 1;
 567       while (!found && i <= n) {
 568         if (grouping[i - 1] != 0) {
 569           same = true;
 570           for (j = 0; j < setsz; j++)
 571             if (grouping[i - 1][j] != tempset[j])
 572               same = false;
 573         }
 574         if ((grouping[i - 1] != 0) && same)
 575           found = true;
 576         else
 577           i++;
 578       }
 579       zero = true;
 580       for (j = 0; j < setsz; j++)
 581         if (tempset[j] != 0)
 582           zero = false;
 583       if (!zero && !found)
 584         bigsubset(tempset, n);
 585     }
 586   }
 587   q->next = *p;
 588   free(tempset);
 589   free(st2);
 590 }  /* recontraverse */
 591
 592
 593 void reconstruct(long n)
 594 {
 595   /* reconstruct tree from the subsets */
 596   long nextnode;
 597   group_type *s;
 598
 599   nextnode = spp + 1;
 600   s = (group_type *)Malloc(setsz * sizeof(group_type));
 601   memcpy(s, fullset, setsz * sizeof(group_type));
 602   recontraverse(&root, s, n, &nextnode);
 603   free(s);
 604 }  /* reconstruct */
 605
 606
 607 void coordinates(node *p, long *tipy)
 608 {
 609   /* establishes coordinates of nodes */
 610   node *q, *first, *last;
 611   long maxx;
 612
 613   if (p->tip) {
 614     p->xcoord = 0;
 615     p->ycoord = *tipy;
 616     p->ymin = *tipy;
 617     p->ymax = *tipy;
 618     (*tipy) += down;
 619     return;
 620   }
 621   q = p->next;
 622   maxx = 0;
 623   while (q != p) {
 624     coordinates(q->back, tipy);
 625     if (!q->back->tip) {
 626       if (q->back->xcoord > maxx)
 627         maxx = q->back->xcoord;
 628     }
 629     q = q->next;
 630   }
 631   first = p->next->back;
 632   q = p;
 633   while (q->next != p)
 634     q = q->next;
 635   last = q->back;
 636   p->xcoord = maxx + OVER;
 637   p->ycoord = (long)((first->ycoord + last->ycoord) / 2);
 638   p->ymin = first->ymin;
 639   p->ymax = last->ymax;
 640 }  /* coordinates */
 641
 642
 643 void drawline(long i)
 644 {
 645   /* draws one row of the tree diagram by moving up tree */
 646   node *p, *q;
 647   long n, j;
 648   boolean extra, done, trif;
 649   node *r,  *first = NULL, *last = NULL;
 650   boolean found;
 651
 652   p = root;
 653   q = root;
 654   fprintf(outfile, "  ");
 655   extra = false;
 656   trif = false;
 657   do {
 658     if (!p->tip) {
 659       found = false;
 660       r = p->next;
 661       while (r != p && !found) {
 662         if (i >= r->back->ymin && i <= r->back->ymax) {
 663           q = r->back;
 664           found = true;
 665         } else
 666           r = r->next;
 667       }
 668       first = p->next->back;
 669       r = p;
 670       while (r->next != p)
 671         r = r->next;
 672       last = r->back;
 673     }
 674     done = (p->tip || p == q);
 675     n = p->xcoord - q->xcoord;
 676     if (extra) {
 677       n--;
 678       extra = false;
 679     }
 680     if (q->ycoord == i && !done) {
 681       if (trif)
 682         putc('-', outfile);
 683       else
 684         putc('+', outfile);
 685       trif = false;
 686       if (!q->tip) {
 687         for (j = 1; j <= n - 7; j++)
 688           putc('-', outfile);
 689         if (noroot && (root->next->next->next == root) &&
 690             (((root->next->back == q) && root->next->next->back->tip)
 691              || ((root->next->next->back == q) && root->next->back->tip)))
 692           fprintf(outfile, "------|");
 693         else {
 694           if (!strict) {   /* write number of times seen */
 695             if (q->deltav >= 100)
 696               fprintf(outfile, "%5.1f-|", (double)q->deltav);
 697             else if (q->deltav >= 10)
 698               fprintf(outfile, "-%4.1f-|", (double)q->deltav);
 699             else
 700               fprintf(outfile, "--%3.1f-|", (double)q->deltav);
 701           } else
 702             fprintf(outfile, "------|");
 703         }
 704         extra = true;
 705         trif = true;
 706       } else {
 707         for (j = 1; j < n; j++)
 708           putc('-', outfile);
 709       }
 710     } else if (!p->tip && last->ycoord > i && first->ycoord < i &&
 711                (i != p->ycoord || p == root)) {
 712       putc('|', outfile);
 713       for (j = 1; j < n; j++)
 714         putc(' ', outfile);
 715     } else {
 716       for (j = 1; j <= n; j++)
 717         putc(' ', outfile);
 718       if (trif)
 719         trif = false;
 720     }
 721     if (q != p)
 722       p = q;
 723   } while (!done);
 724   if (p->ycoord == i && p->tip) {
 725     for (j = 0; (j < MAXNCH) && (p->nayme[j] != '\0'); j++)
 726       putc(p->nayme[j], outfile);
 727   }
 728   putc('\n', outfile);
 729 }  /* drawline */
 730
 731
 732 void printree()
 733 {
 734   /* prints out diagram of the tree */
 735   long i;
 736   long tipy;
 737
 738   if (treeprint) {
 739     fprintf(outfile, "\nCONSENSUS TREE:\n");
 740     if (mr || mre || ml) {
 741       if (noroot) {
 742         fprintf(outfile, "the numbers on the branches indicate the number\n");
 743  fprintf(outfile, "of times the partition of the species into the two sets\n");
 744         fprintf(outfile, "which are separated by that branch occurred\n");
 745       } else {
 746         fprintf(outfile, "the numbers forks indicate the number\n");
 747         fprintf(outfile, "of times the group consisting of the species\n");
 748         fprintf(outfile, "which are to the right of that fork occurred\n");
 749       }
 750       fprintf(outfile, "among the trees, out of %6.2f trees\n", ntrees);
 751       if (ntrees <= 1.001)
 752         fprintf(outfile, "(trees had fractional weights)\n");
 753     }
 754     tipy = 1;
 755     coordinates(root, &tipy);
 756     putc('\n', outfile);
 757     for (i = 1; i <= tipy - down; i++)
 758       drawline(i);
 759     putc('\n', outfile);
 760   }
 761   if (noroot) {
 762     fprintf(outfile, "\n  remember:");
 763     if (didreroot)
 764       fprintf(outfile, " (though rerooted by outgroup)");
 765     fprintf(outfile, " this is an unrooted tree!\n");
 766   }
 767   putc('\n', outfile);
 768 }  /* printree */
 769
 770
 771 void enternohash(group_type *s, long *n)
 772 {
 773   /* if set s is already there, enter it into groupings in the next slot
 774      (without hash-coding).  n is number of sets stored there and is updated */
 775   long i, j;
 776   boolean found;
 777
 778   found = false;
 779   for (i = 0; i < (*n); i++) {  /* go through looking whether it is there */
 780     found = true;
 781     for (j = 0; j < setsz; j++) {  /* check both parts of partition */
 782       found = found && (grouping[i][j] == s[j]);
 783       found = found && (group2[i][j] == (fullset[j] & (~s[j])));
 784     }
 785     if (found)
 786       break;
 787   }
 788   if (!found) {    /* if not, add it to the slot after the end,
 789                       which must be empty */
 790     grouping[i] = (group_type *)Malloc(setsz * sizeof(group_type));
 791     timesseen[i] = (double *)Malloc(sizeof(double));
 792     group2[i] = (group_type *)Malloc(setsz * sizeof(group_type));
 793     for (j = 0; j < setsz; j++)
 794       grouping[i][j] = s[j];
 795     *timesseen[i] = 1;
 796     (*n)++;
 797   }
 798 } /* enternohash */
 799
 800
 801 void enterpartition (group_type *s1, long *n)
 802 {
 803   /* try to put this partition in list of partitions.  If implied by others,
 804      don't bother.  If others implied by it, replace them.  If this one
 805      vacuous because only one element in s1, forget it */
 806   long i, j;
 807   boolean found;
 808
 809 /* this stuff all to be rewritten but left here so pieces can be used */
 810   found = false;
 811   for (i = 0; i < (*n); i++) {  /* go through looking whether it is there */
 812     found = true;
 813     for (j = 0; j < setsz; j++) {  /* check both parts of partition */
 814       found = found && (grouping[i][j] == s1[j]);
 815       found = found && (group2[i][j] == (fullset[j] & (~s1[j])));
 816     }
 817     if (found)
 818       break;
 819   }
 820   if (!found) {    /* if not, add it to the slot after the end,
 821                       which must be empty */
 822     grouping[i] = (group_type *)Malloc(setsz * sizeof(group_type));
 823     timesseen[i] = (double *)Malloc(sizeof(double));
 824     group2[i] = (group_type *)Malloc(setsz * sizeof(group_type));
 825     for (j = 0; j < setsz; j++)
 826       grouping[i][j] = s1[j];
 827     *timesseen[i] = 1;
 828     (*n)++;
 829   }
 830 } /* enterpartition */
 831
 832
 833 void elimboth(long n)
 834 {
 835   /* for Adams case: eliminate pairs of groups incompatible with each other */
 836   long i, j;
 837   boolean comp;
 838
 839   for (i = 0; i < n-1; i++) {
 840     for (j = i+1; j < n; j++) {
 841       comp = compatible(i,j);
 842       if (!comp) {
 843         *timesseen[i] = 0.0;
 844         *timesseen[j] = 0.0;
 845       }
 846     }
 847     if (*timesseen[i] == 0.0) {
 848       free(grouping[i]);
 849       free(timesseen[i]);
 850       timesseen[i] = NULL;
 851       grouping[i] = NULL;
 852     }
 853   }
 854   if (*timesseen[n-1] == 0.0) {
 855     free(grouping[n-1]);
 856     free(timesseen[n-1]);
 857     timesseen[n-1] = NULL;
 858     grouping[n-1] = NULL;
 859   }
 860 }  /* elimboth */
 861
 862
 863 void consensus(pattern_elm ***pattern_array, long trees_in)
 864 {
 865   long i, n, n2, tipy;
 866
 867   group2 = (group_type **)  Malloc(maxgrp*sizeof(group_type *));
 868   for (i = 0; i < maxgrp; i++)
 869     group2[i] = NULL;
 870   times2 = (double **)Malloc(maxgrp*sizeof(double *));
 871   for (i = 0; i < maxgrp; i++)
 872     times2[i] = NULL;
 873   n2 = 0;
 874   censor();                /* drop groups that are too rare */
 875   compress(&n);            /* push everybody to front of array */
 876   if (!strict) {           /* drop those incompatible, if any */
 877     sort(n);
 878     eliminate(&n, &n2);
 879     compress(&n);
 880     }
 881   reconstruct(n);
 882   tipy = 1;
 883   coordinates(root, &tipy);
 884   if (prntsets) {
 885     fprintf(outfile, "\nSets included in the consensus tree\n");
 886     printset(n);
 887     for (i = 0; i < n2; i++) {
 888       if (!grouping[i]) {
 889         grouping[i] = (group_type *)Malloc(setsz * sizeof(group_type));
 890         timesseen[i] = (double *)Malloc(sizeof(double));
 891         }
 892       memcpy(grouping[i], group2[i], setsz * sizeof(group_type));
 893       *timesseen[i] = *times2[i];
 894     }
 895     n = n2;
 896     fprintf(outfile, "\n\nSets NOT included in consensus tree:");
 897     if (n2 == 0)
 898       fprintf(outfile, " NONE\n");
 899     else {
 900       putc('\n', outfile);
 901       printset(n);
 902     }
 903   }
 904   putc('\n', outfile);
 905   if (strict)
 906     fprintf(outfile, "\nStrict consensus tree\n");
 907   if (mre)
 908     fprintf(outfile, "\nExtended majority rule consensus tree\n");
 909   if (ml) {
 910     fprintf(outfile, "\nM  consensus tree (l = %4.2f)\n", mlfrac);
 911     fprintf(outfile, " l\n");
 912   }
 913   if (mr)
 914     fprintf(outfile, "\nMajority rule consensus tree\n");
 915   printree();
 916   free(nayme);
 917   for (i = 0; i < maxgrp; i++)
 918     free(grouping[i]);
 919   free(grouping);
 920   for (i = 0; i < maxgrp; i++)
 921     free(order[i]);
 922   free(order);
 923   for (i = 0; i < maxgrp; i++)
 924     if (timesseen[i] != NULL)
 925       free(timesseen[i]);
 926   free(timesseen);
 927 }  /* consensus */
 928
 929
 930 void rehash()
 931 {
 932   group_type *s;
 933   long i, j, k;
 934   double temp, ss, smult;
 935   boolean done;
 936
 937   smult = (sqrt(5.0) - 1) / 2;
 938   s = (group_type *)Malloc(setsz * sizeof(group_type));
 939   for (i = 0; i < maxgrp/2; i++) {
 940     k = *order[i];
 941     memcpy(s, grouping[k], setsz * sizeof(group_type));
 942     ss = 0.0;
 943     for (j = 0; j < setsz; j++)
 944       ss += s[j] /* pow(2, SETBITS*j)*/;
 945     temp = ss * smult;
 946     j = (long)(maxgrp * (temp - floor(temp)));
 947     done = false;
 948     while (!done) {
 949       if (!grping2[j]) {
 950         grping2[j] = (group_type *)Malloc(setsz * sizeof(group_type));
 951         order2[i] = (long *)Malloc(sizeof(long));
 952         tmseen2[j] = (double *)Malloc(sizeof(double));
 953         memcpy(grping2[j], grouping[k], setsz * sizeof(group_type));
 954         *tmseen2[j] = *timesseen[k];
 955         *order2[i] = j;
 956         grouping[k] = NULL;
 957         timesseen[k] = NULL;
 958         order[i] = NULL;
 959         done = true;
 960       } else {
 961         j++;
 962         if (j >= maxgrp) j -= maxgrp;
 963       }
 964     }
 965   }
 966   free(s);
 967 }  /* rehash */
 968
 969
 970 void enternodeset(node* r)
 971 { /* enter a set of species into the hash table */
 972   long i, j, start;
 973   double ss, n;
 974   boolean done, same;
 975   double times ;
 976   group_type *s;
 977
 978   s = r->nodeset;
 979   same = true;
 980   for (i = 0; i < setsz; i++)
 981     if (s[i] != fullset[i])
 982       same = false;
 983   if (same)
 984     return;
 985   times = trweight;
 986   ss = 0.0;                        /* compute the hashcode for the set */
 987   n = ((sqrt(5.0) - 1.0) / 2.0);   /* use an irrational multiplier */
 988   for (i = 0; i < setsz; i++)
 989     ss += s[i] * n;
 990   i = (long)(maxgrp * (ss - floor(ss))) + 1; /* use fractional part of code */
 991   start = i;
 992   done = false;                   /* go through seeing if it is there */
 993   while (!done) {
 994     if (grouping[i - 1]) {        /* ... i.e. if group is absent, or  */
 995       same = false;               /* (will be false if timesseen = 0) */
 996       if (!(timesseen[i-1] == 0)) {  /* ... if timesseen = 0 */
 997         same = true;
 998         for (j = 0; j < setsz; j++) {
 999           if (s[j] != grouping[i - 1][j])
1000             same = false;
1001         }
1002       }
1003     }
1004     if (grouping[i - 1] && same) {  /* if it is there, increment timesseen */
1005       *timesseen[i - 1] += times;
1006       lengths[i - 1] = nodep[r->index - 1]->v;
1007       done = true;
1008     } else if (!grouping[i - 1]) {  /* if not there and slot empty ... */
1009       grouping[i - 1] = (group_type *)Malloc(setsz * sizeof(group_type));
1010       lasti++;
1011       order[lasti] = (long *)Malloc(sizeof(long));
1012       timesseen[i - 1] = (double *)Malloc(sizeof(double));
1013       memcpy(grouping[i - 1], s, setsz * sizeof(group_type));
1014       *timesseen[i - 1] = times;
1015       *order[lasti] = i - 1;
1016       done = true;
1017       lengths[i - 1] = nodep[r->index -1]->v;
1018     } else {  /* otherwise look to put it in next slot ... */
1019       i++;
1020       if (i > maxgrp) i -= maxgrp;
1021     }
1022     if (!done && i == start) {  /* if no place to put it, expand hash table */
1023       maxgrp = maxgrp*2;
1024       tmseen2 = (double **)Malloc(maxgrp*sizeof(double *));
1025       for (j = 0; j < maxgrp; j++)
1026         tmseen2[j] = NULL;
1027       grping2 = (group_type **)Malloc(maxgrp*sizeof(group_type *));
1028       for (j = 0; j < maxgrp; j++)
1029         grping2[j] = NULL;
1030       order2 = (long **)Malloc(maxgrp*sizeof(long *));
1031       for (j = 0; j < maxgrp; j++)
1032         order2[j] = NULL;
1033       lengths2 = (double *)Malloc(maxgrp*sizeof(double));
1034       for (j = 0; j < maxgrp; j++)
1035         lengths2[j] = 0.0;
1036       memcpy(lengths2,lengths,maxgrp*sizeof(double) / 2);
1037       rehash();
1038       free(lengths);
1039       free(timesseen);
1040       free(grouping);
1041       free(order);
1042       timesseen = tmseen2;
1043       grouping = grping2;
1044       lengths = lengths2;
1045       order = order2;
1046       done = true;
1047       lasti = maxgrp/2 - 1;
1048       enternodeset(r);
1049     }
1050   }
1051 }  /* enternodeset */
1052
1053
1054 void accumulate(node *r)
1055 {
1056   node *q;
1057   long i;
1058
1059   if (r->tip) {
1060     if (!r->nodeset)
1061       r->nodeset = (group_type *)Malloc(setsz * sizeof(group_type));
1062     for (i = 0; i < setsz; i++)
1063       r->nodeset[i] = 0L;
1064     i = (r->index-1) / (long)SETBITS;
1065     r->nodeset[i] = 1L << (r->index - 1 - i*SETBITS);
1066   }
1067   else {
1068     q = r->next;
1069     while (q != r) {
1070       accumulate(q->back);
1071       q = q->next;
1072     }
1073     q = r->next;
1074     if (!r->nodeset)
1075       r->nodeset = (group_type *)Malloc(setsz * sizeof(group_type));
1076     for (i = 0; i < setsz; i++)
1077       r->nodeset[i] = 0;
1078     while (q != r) {
1079       for (i = 0; i < setsz; i++)
1080         r->nodeset[i] |= q->back->nodeset[i];
1081       q = q->next;
1082     }
1083   }
1084   if ((!r->tip && (r->next->next != r)) || r->tip)
1085     enternodeset(r);
1086 }  /* accumulate */
1087
1088
1089 void dupname2(Char *name, node *p, node *this)
1090 {
1091   /* search for a duplicate name recursively */
1092   node *q;
1093
1094   if (p->tip) {
1095     if (p != this) {
1096       if (namesSearch(p->nayme)) {
1097         printf("\n\nERROR in user tree: duplicate name found: ");
1098         puts(p->nayme);
1099         printf("\n\n");
1100         exxit(-1);
1101       } else {
1102         namesAdd(p->nayme);
1103       }
1104     }
1105   } else {
1106     q = p;
1107     while (p->next != q) {
1108       dupname2(name, p->next->back, this);
1109       p = p->next;
1110     }
1111   }
1112 }  /* dupname2 */
1113
1114
1115 void dupname(node *p)
1116 {
1117   /* search for a duplicate name in tree */
1118   node *q;
1119
1120   if (p->tip) {
1121     if (namesSearch(p->nayme)) {
1122       printf("\n\nERROR in user tree: duplicate name found: ");
1123       puts(p->nayme);
1124       printf("\n\n");
1125       exxit(-1);
1126     } else {
1127       namesAdd(p->nayme);
1128     }
1129   } else {
1130     q = p;
1131     while (p->next != q) {
1132       dupname(p->next->back);
1133       p = p->next;
1134     }
1135   }
1136 }  /* dupname */
1137
1138
1139 void missingnameRecurs(node *p)
1140 {
1141   /* search for missing names in first tree */
1142   node *q;
1143
1144   if (p->tip) {
1145     if (!namesSearch(p->nayme)) {
1146       printf("\n\nERROR in user tree: name %s not found in first tree\n\n\n",
1147              p->nayme);
1148       exxit(-1);
1149     }
1150   } else {
1151     q = p;
1152     while (p->next != q) {
1153       missingnameRecurs(p->next->back);
1154       p = p->next;
1155     }
1156   }
1157 }  /* missingnameRecurs */
1158
1159 /**
1160  * wrapper for recursive missingname function
1161  */
1162 void missingname(node *p){
1163   missingnameRecurs(p);
1164   namesCheckTable();
1165 } /* missingname */
1166
1167 void gdispose(node *p)
1168 {
1169   /* go through tree throwing away nodes */
1170   node *q, *r;
1171
1172   if (p->tip) {
1173     chuck(&grbg, p);
1174     return;
1175   }
1176   q = p->next;
1177   while (q != p) {
1178     gdispose(q->back);
1179     r = q;
1180     q = q->next;
1181     chuck(&grbg, r);
1182   }
1183   chuck(&grbg, q);
1184 }  /* gdispose */
1185
1186
1187 void initreenode(node *p)
1188 {
1189   /* traverse tree and assign species names to tip nodes */
1190   node *q;
1191
1192   if (p->tip) {
1193     memcpy(nayme[p->index - 1], p->nayme, MAXNCH);
1194   } else {
1195     q = p->next;
1196     while (q && q != p) {
1197       initreenode(q->back);
1198       q = q->next;
1199     }
1200   }
1201 } /* initreenode */
1202
1203
1204 void reroot(node *outgroup, long *nextnode)
1205 {
1206   /* reorients tree, putting outgroup in desired position. */
1207   long i;
1208   boolean nroot;
1209   node *p, *q;
1210   double newv;
1211
1212   nroot = false;
1213   p = root->next;
1214   while (p != root) {
1215     if ((outgroup->back == p) && (root->next->next->next == root)) {
1216       nroot = true;
1217       p = root;
1218     } else
1219       p = p->next;
1220   }
1221   if (nroot && root->next->next->next == root) {
1222     root->next->next->back->v += root->next->back->v;
1223     root->next->back->v = 0;
1224   }
1225   if (nroot) return;
1226
1227   p = root;
1228   i = 0;
1229   while (p->next != root) {
1230     p = p->next;
1231     i++;
1232   }
1233   if (i == 2) {
1234     newv = root->next->back->v + root->next->next->back->v;
1235     root->next->back->back = p->back;
1236     p->back->back = root->next->back;
1237     q = root->next;
1238     p->back->v = newv;
1239     q->back->v = newv;
1240   } else {
1241     p->next = root->next;
1242     nodep[root->index-1] = root->next;
1243     gnu(&grbg, &root->next);
1244     q = root->next;
1245     gnu(&grbg, &q->next);
1246     p = q->next;
1247     p->next = root;
1248     q->tip = false;
1249     p->tip = false;
1250     nodep[*nextnode] = root;
1251     (*nextnode)++;
1252     root->index = *nextnode;
1253     root->next->index = root->index;
1254     root->next->next->index = root->index;
1255   }
1256   newv = outgroup->v;
1257   q->back = outgroup;
1258   p->back = outgroup->back;
1259   outgroup->back->back = p;
1260   outgroup->back = q;
1261   outgroup->v = 0;
1262   outgroup->back->v = 0;
1263   root->v = 0;
1264   p->v = newv;
1265   p->back->v = newv;
1266   reorient(root);
1267 }  /* reroot */
1268
1269
1270 void reorient(node* n) {
1271   node* p;
1272
1273   if ( n->tip ) return;
1274   if ( nodep[n->index - 1] != n )  {
1275     nodep[n->index - 1] = n;
1276     if ( n->back )
1277       n->v = n->back->v;
1278   }
1279
1280   for ( p = n->next ; p != n ; p = p->next)
1281     reorient(p->back);
1282 }
1283
1284
1285 void store_pattern (pattern_elm ***pattern_array,
1286                         double *timesseen_changes, int trees_in_file)
1287 {
1288   /* put a tree's groups into a pattern array.
1289      Don't forget that when not Adams, grouping[] is not compressed. . . */
1290   long i, total_groups=0, j=0, k;
1291
1292   /* First, find out how many groups exist in the given tree. */
1293   for (i = 0 ; i < maxgrp ; i++)
1294     if ((grouping[i] != NULL) &&
1295        (*timesseen[i] > timesseen_changes[i]))
1296       /* If this is group exists and is present in the current tree, */
1297       total_groups++ ;
1298
1299   /* Then allocate a space to store the bit patterns. . . */
1300   for (i = 0 ; i < setsz ; i++) {
1301     pattern_array[i][trees_in_file]
1302       = (pattern_elm *) Malloc(sizeof(pattern_elm)) ;
1303     pattern_array[i][trees_in_file]->apattern =
1304       (group_type *) Malloc (total_groups * sizeof (group_type)) ;
1305     pattern_array[i][trees_in_file]->length =
1306       (double *) Malloc (maxgrp * sizeof (double)) ;
1307       for ( j = 0 ; j < maxgrp ; j++ ) {
1308         pattern_array[i][trees_in_file]->length[j] = -1;
1309       }
1310     pattern_array[i][trees_in_file]->patternsize = (long *)Malloc(sizeof(long));
1311   }
1312   j = 0;
1313   /* Then go through groupings again, and copy in each element
1314      appropriately. */
1315   for (i = 0 ; i < maxgrp ; i++)
1316     if (grouping[i] != NULL) {
1317       if (*timesseen[i] > timesseen_changes[i]) {
1318         for (k = 0 ; k < setsz ; k++)
1319           pattern_array[k][trees_in_file]->apattern[j] = grouping[i][k] ;
1320         pattern_array[0][trees_in_file]->length[j] = lengths[i];
1321         j++ ;
1322         timesseen_changes[i] = *timesseen[i] ;
1323       }
1324     }
1325   *pattern_array[0][trees_in_file]->patternsize = total_groups;
1326 }  /* store_pattern */
1327
1328
1329 boolean samename(naym name1, plotstring name2)
1330 {
1331   return !(strncmp(name1, name2, MAXNCH));
1332 }  /* samename */
1333
1334
1335 void reordertips()
1336 {
1337   /* matchs tip nodes to species names first read in */
1338   long i, j;
1339   boolean done;
1340   node *p, *q, *r;
1341   for (i = 0;  i < spp; i++) {
1342     j = 0;
1343     done = false;
1344     do {
1345       if (samename(nayme[i], nodep[j]->nayme)) {
1346         done =  true;
1347         if (i != j) {
1348           p = nodep[i];
1349           q = nodep[j];
1350           r = p->back;
1351           p->back->back = q;
1352           q->back->back = p;
1353           p->back =  q->back;
1354           q->back = r;
1355           memcpy(q->nayme, p->nayme, MAXNCH);
1356           memcpy(p->nayme, nayme[i], MAXNCH);
1357         }
1358       }
1359       j++;
1360     } while (j < spp && !done);
1361   }
1362 }  /* reordertips */
1363
1364
1365 void read_groups (pattern_elm ****pattern_array,double *timesseen_changes,
1366                         long trees_in_1, long total_trees, FILE *intree)
1367 {
1368   /* read the trees.  Accumulate sets. */
1369   int i, j, k;
1370   boolean haslengths, initial;
1371   long nextnode, trees_read = 0;
1372
1373   /* set up the groupings array and the timesseen array */
1374   grouping  = (group_type **)  Malloc(maxgrp*sizeof(group_type *));
1375   lengths  = (double *)  Malloc(maxgrp*sizeof(double));
1376   for (i = 0; i < maxgrp; i++)
1377     grouping[i] = NULL;
1378   order     = (long **) Malloc(maxgrp*sizeof(long *));
1379   for (i = 0; i < maxgrp; i++)
1380     order[i] = NULL;
1381   timesseen = (double **)Malloc(maxgrp*sizeof(double *));
1382   for (i = 0; i < maxgrp; i++)
1383     timesseen[i] = NULL;
1384
1385   firsttree = true;
1386   grbg = NULL;
1387   initial = true;
1388   while (!eoff(intree)) {          /* go till end of input tree file */
1389     for (i = 0; i < maxgrp; i++) {
1390       lengths[i] = -1;
1391     }
1392     goteof = false;
1393     nextnode = 0;
1394     haslengths = true;
1395     allocate_nodep(&nodep, &intree, &spp);
1396     if (firsttree)
1397       nayme = (naym *)Malloc(spp*sizeof(naym));
1398     treeread(intree, &root, treenode, &goteof, &firsttree, nodep,
1399               &nextnode, &haslengths, &grbg, initconsnode,true,-1);
1400     if (!initial) {
1401       reordertips();
1402       missingname(root);
1403     } else {
1404       initial = false;
1405       hashp = (hashtype)Malloc(sizeof(namenode) * NUM_BUCKETS);
1406       for (i=0;i<NUM_BUCKETS;i++) {
1407         hashp[i] = NULL;
1408       }
1409       dupname(root);
1410       initreenode(root);
1411       setsz = (long)ceil((double)spp/(double)SETBITS);
1412       if (tree_pairing != NO_PAIRING)
1413         {
1414           /* Now that we know setsz, we can malloc pattern_array
1415           and pattern_array[n] accordingly. */
1416           (*pattern_array) =
1417             (pattern_elm ***)Malloc(setsz * sizeof(pattern_elm **));
1418
1419           /* For this assignment, let's assume that there will be no
1420              more than maxtrees. */
1421           for (j = 0 ; j < setsz ; j++)
1422               (*pattern_array)[j] =
1423                 (pattern_elm **)Malloc(total_trees * sizeof(pattern_elm *));
1424         }
1425
1426       fullset = (group_type *)Malloc(setsz * sizeof(group_type));
1427       for (j = 0; j < setsz; j++)
1428         fullset[j] = 0L;
1429       k = 0;
1430       for (j = 1; j <= spp; j++) {
1431         if (j == ((k+1)*SETBITS+1)) k++;
1432         fullset[k] |= 1L << (j - k*SETBITS - 1);
1433       }
1434     }
1435     if (goteof)
1436       continue;
1437     ntrees += trweight;
1438     if (noroot) {
1439       reroot(nodep[outgrno - 1], &nextnode);
1440       didreroot = outgropt;
1441     }
1442     accumulate(root);
1443     gdispose(root);
1444     for (j = 0; j < 2*(1+spp); j++)
1445       nodep[j] = NULL;
1446     free(nodep);
1447     /* Added by Dan F. */
1448     if (tree_pairing != NO_PAIRING) {
1449         /* If we're computing pairing or need separate tree sets, store the
1450            current pattern as an element of it's trees array. */
1451       store_pattern ((*pattern_array), timesseen_changes, trees_read) ;
1452       trees_read++ ;
1453     }
1454   }
1455 } /* read_groups */
1456
1457