source: tags/initial/GDE/PHYLIP/consense.c

Last change on this file was 2, checked in by oldcode, 24 years ago

Initial revision
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 29.7 KB
Line 
1#include "phylip.h"
2
3/*
4 * version 3.56c. (c) Copyright 1993 by Joseph Felsenstein. Written by Joseph
5 * Felsenstein, Hisashi Horino, Akiko Fuseki, Sean Lamont, and Andrew Keeffe.
6 * Permission is granted to copy and use this program provided no fee is
7 * charged for it and provided that this copyright notice is not removed.
8 */
9
10#define nmlngth         20      /* max. no. of characters in species name  */
11
12#define ibmpc0          false
13#define ansi0           true
14#define vt520           false
15
16/* nodes will form a binary tree */
17
18typedef Char    naym[nmlngth];
19
20typedef struct node {
21        /* describes a tip species or an ancestor */
22        struct node    *next, *back;    /* pointers to nodes                */
23        double          index;  /* number of the node               */
24        boolean         tip;    /* present species are tips of tree */
25        naym            naime;
26        long           *nodeset;/* used by accumulate               */
27        short           xcoord, ycoord, ymin;   /* used by printree                 */
28        short           ymax;
29}               node;
30
31typedef node  **pointarray;
32
33Static node    *root;
34Static double   trweight, ntrees;
35Static FILE    *infile, *outfile, *treefile;
36Static short    spp, numopts, outgrno, i, j, col, setsz;
37Static long     maxgrp;         /* max. no. of groups in all trees found  */
38Static boolean  anerror, trout, first, noroot, outgropt, didreroot, prntsets, progress,
39                treeprint, ibmpc, vt52, ansi, goteof;
40Static pointarray treenode;     /* pointers to all nodes in tree */
41Static naym    *nayme;          /* names of species       */
42Static long   **grouping, **grping2, **group2;  /* to store groups found  */
43Static long   **order, **order2, lasti;
44Static double **timesseen, **tmseen2, **times2; /* how often they're seen */
45Static long    *fullset;
46Static node    *garbage;
47
48Static short    tipy;
49
50openfile(fp, filename, mode, application, perm)
51        FILE          **fp;
52        char           *filename;
53        char           *mode;
54        char           *application;
55        char           *perm;
56{
57        FILE           *of;
58        char            file[100];
59        strcpy(file, filename);
60        while (1) {
61                of = fopen(file, mode);
62                if (of)
63                        break;
64                else {
65                        switch (*mode) {
66                        case 'r':
67                                printf("%s:  can't read %s\n", application, file);
68                                printf("Please enter a new filename>");
69                                gets(file);
70                                break;
71                        case 'w':
72                                printf("%s: can't write %s\n", application, file);
73                                printf("Please enter a new filename>");
74                                gets(file);
75                                break;
76                        }
77                }
78        }
79        *fp = of;
80        if (perm != NULL)
81                strcpy(perm, file);
82}
83
84
85Static Void
86gnu(p)
87        node          **p;
88{
89        /*
90         * this and the following are do-it-yourself garbage collectors. Make
91         * a new node or pull one off the garbage list
92         */
93        if (garbage != NULL) {
94                *p = garbage;
95                garbage = garbage->next;
96        } else
97                *p = (node *) Malloc(sizeof(node));
98        (*p)->next = NULL;
99        (*p)->tip = false;
100}                               /* gnu */
101
102
103Static Void
104chuck(p)
105        node           *p;
106{
107        /* collect garbage on p -- put it on front of garbage list */
108        p->next = garbage;
109        garbage = p;
110}                               /* chuck */
111
112
113Static Void
114gdispose(p)
115        node           *p;
116{
117        /* go through tree throwing away nodes */
118        node           *q, *r;
119
120        if (p->tip)
121                return;
122        q = p->next;
123        while (q != p) {
124                gdispose(q->back);
125                r = q;
126                q = q->next;
127                chuck(r);
128        }
129        chuck(q);
130}                               /* gdispose */
131
132
133Static Void
134uppercase(ch)
135        Char           *ch;
136{                               /* convert ch to upper case -- either ASCII
137                                 * or EBCDIC */
138        *ch = (islower(*ch) ? toupper(*ch) : (*ch));
139}                               /* uppercase */
140
141
142Static Void
143getoptions()
144{
145        /* interactively set options */
146        Char            ch;
147        boolean         done, done1;
148
149        fprintf(outfile, "\nMajority-rule and strict consensus tree");
150        fprintf(outfile, " program, version %s\n\n", VERSION);
151        putchar('\n');
152        anerror = false;
153        noroot = true;
154        numopts = 0;
155        outgrno = 1;
156        outgropt = false;
157        trout = true;
158        prntsets = true;
159        progress = true;
160        treeprint = true;
161        do {
162                if (ansi)
163                        printf("\033[2J\033[H");
164                else if (vt52)
165                        printf("\033E\033H");
166                else
167                        putchar('\n');
168                printf("\nMajority-rule and strict consensus tree");
169                printf(" program, version %s\n\n", VERSION);
170                printf("Settings for this run:\n");
171                if (noroot) {
172                        printf("  O                        Outgroup root?");
173                        if (outgropt)
174                                printf("  Yes, at species number%3hd\n", outgrno);
175                        else
176                                printf("  No, use as outgroup species%3hd\n", outgrno);
177                }
178                printf("  R        Trees to be treated as Rooted?");
179                if (noroot)
180                        printf("  No\n");
181                else
182                        printf("  Yes\n");
183                printf("  0   Terminal type (IBM PC, VT52, ANSI)?");
184                if (ibmpc)
185                        printf("  IBM PC\n");
186                if (ansi)
187                        printf("  ANSI\n");
188                if (vt52)
189                        printf("  VT52\n");
190                if (!(ibmpc || vt52 || ansi))
191                        printf("  (none)\n");
192                printf("  1         Print out the sets of species");
193                if (prntsets)
194                        printf("  Yes\n");
195                else
196                        printf("  No\n");
197                printf("  2  Print indications of progress of run  %s\n",
198                       (progress ? "Yes" : "No"));
199                printf("  3                        Print out tree");
200                if (treeprint)
201                        printf("  Yes\n");
202                else
203                        printf("  No\n");
204                printf("  4       Write out trees onto tree file?");
205                if (trout)
206                        printf("  Yes\n");
207                else
208                        printf("  No\n");
209                printf("\nAre these settings correct? (type Y or the letter for one to change)\n");
210                scanf("%c%*[^\n]", &ch);
211                getchar();
212                uppercase(&ch);
213                done = (ch == 'Y');
214                if (!done) {
215                        if ((noroot && (ch == 'O')) || ch == 'R' || ch == '0' || ch == '1' ||
216                            ch == '2' || ch == '3' || ch == '4') {
217                                switch (ch) {
218
219                                case 'O':
220                                        outgropt = !outgropt;
221                                        if (outgropt) {
222                                                numopts++;
223                                                done1 = true;
224                                                do {
225                                                        printf("Type number of the outgroup:\n");
226                                                        scanf("%hd%*[^\n]", &outgrno);
227                                                        getchar();
228                                                        done1 = (outgrno >= 1);
229                                                        if (!done1) {
230                                                                printf("BAD OUTGROUP NUMBER: %4hd\n", outgrno);
231                                                                printf("  Must be greater than zero\n");
232                                                        }
233                                                } while (done1 != true);
234                                        }
235                                        break;
236
237                                case 'R':
238                                        noroot = !noroot;
239                                        break;
240
241                                case '0':
242                                        if (ibmpc) {
243                                                ibmpc = false;
244                                                vt52 = true;
245                                        } else {
246                                                if (vt52) {
247                                                        vt52 = false;
248                                                        ansi = true;
249                                                } else if (ansi)
250                                                        ansi = false;
251                                                else
252                                                        ibmpc = true;
253                                        }
254                                        break;
255
256                                case '1':
257                                        prntsets = !prntsets;
258                                        break;
259
260                                case '2':
261                                        progress = !progress;
262                                        break;
263
264                                case '3':
265                                        treeprint = !treeprint;
266                                        break;
267
268                                case '4':
269                                        trout = !trout;
270                                        break;
271                                }
272                        } else
273                                printf("Not a possible option!\n");
274                }
275        } while (!done);
276}                               /* getoptions */
277
278void 
279getch(c, parens)
280        Char           *c;
281        short          *parens;
282{
283        /* get next nonblank character */
284        do {
285                if (eoln(infile)) {
286                        fscanf(infile, "%*[^\n]");
287                        getc(infile);
288                }
289                *c = getc(infile);
290        } while (!(*c != ' ' || eof(infile)));
291        if (*c == '(')
292                (*parens)++;
293        if (*c == ')')
294                (*parens)--;
295}                               /* getch */
296
297Static Void
298setupnode(p, i)
299        node           *p;
300        short           i;
301{
302        /* initialization of node pointers, variables */
303
304        p->next = NULL;
305        p->back = NULL;
306        p->index = i * 1.0;
307        p->tip = false;
308}                               /* setupnode */
309
310void 
311dupname(name, p, found)
312        Char            name[nmlngth];
313        node           *p;
314        boolean        *found;
315{
316        /* search for a duplicate name recursively */
317        if (p) {
318                if (p->next) {
319                        dupname(name, p->next->back, found);
320                        if (p->next->next)
321                                dupname(name, p->next->next->back, found);
322                }
323                if (p->tip) {
324                        for (i = 0; i < nmlngth; i++) {
325                                if (name[i] != p->naime[i])
326                                        *found = true;
327                        }
328                }
329        }
330}
331
332void 
333addelement(p, q, ch, parens, nextnode)
334        node          **p, *q;
335        Char           *ch;
336        short          *nextnode, *parens;
337{
338        /* recursive procedure adds nodes to user-defined tree */
339        node           *r, *pfirst;
340        short           i, n, len;
341        boolean         found, notlast;
342        Char            str[nmlngth];
343
344        if ((*ch) == '(') {
345                (*nextnode)++;
346                gnu(p);
347                pfirst = *p;
348                notlast = true;
349                while (notlast && !anerror) {
350                        gnu(&(*p)->next);
351                        r = (*p)->next;
352                        getch(ch, parens);
353                        addelement(&r->back, r, ch, parens, nextnode);
354                        if (anerror)
355                                break;
356                        if (r->back != NULL)
357                                *p = r;
358                        if ((*ch) == ')') {
359                                notlast = false;
360                                do {
361                                        getch(ch, parens);
362                                } while ((*ch) != ',' && (*ch) != ')' && (*ch) != '[' && (*ch) != ';');
363                        }
364                }
365                if (!anerror) {
366                        (*p)->next = pfirst;
367                        *p = pfirst;
368                }
369        } else if ((*ch) != ')') {
370                if (!anerror) {
371                        for (i = 0; i < nmlngth; i++)
372                                str[i] = ' ';
373                        n = 1;
374                        do {
375                                if ((*ch) == '_')
376                                        (*ch) = ' ';
377                                str[n - 1] = (*ch);
378                                if (eoln(infile)) {
379                                        fscanf(infile, "%*[^\n]");
380                                        getc(infile);
381                                }
382                                (*ch) = getc(infile);
383                                n++;
384                        } while ((*ch) != ':' && (*ch) != ',' && (*ch) != ')' &&
385                              (*ch) != '[' && (*ch) != ';' && n <= nmlngth);
386                        if ((*ch) == ')')
387                                (*parens)--;
388                        len = n - 1;
389                        if (first && !anerror) {
390                                spp++;
391                                if (!anerror) {
392                                        found = false;
393                                        dupname(str, root, &found);
394                                        if (found) {
395                                                printf("\nERROR IN USER TREE: DUPLICATE NAME FOUND: ");
396                                                for (i = 0; i < nmlngth; i++)
397                                                        putchar(str[i]);
398                                                putchar('\n');
399                                                anerror = true;
400                                        }
401                                }
402                                if (!anerror) {
403                                        gnu(p);
404                                        setupnode(*p, spp);
405                                        (*p)->tip = true;
406                                        for (i = 0; i < nmlngth; i++)
407                                                (*p)->naime[i] = ' ';
408                                        for (i = 0; i < len; i++)
409                                                (*p)->naime[i] = str[i];
410                                        if (prntsets) {
411                                                fprintf(outfile, "  ");
412                                                for (i = 0; i < len; i++)
413                                                        putc(str[i], outfile);
414                                                putc('\n', outfile);
415                                        }
416                                }
417                        } else {
418                                n = 1;
419                                do {
420                                        found = true;
421                                        for (i = 0; i < nmlngth; i++)
422                                                found = (found && str[i] == nayme[n - 1][i]);
423                                        if (found)
424                                                *p = treenode[n - 1];
425                                        else
426                                                n++;
427                                } while (!(n > spp || found));
428                                if (n > spp) {
429                                        anerror = true;
430                                        printf("ERROR: CANNOT FIND SPECIES: ");
431                                        for (i = 0; i < nmlngth; i++)
432                                                putchar(str[i]);
433                                        putchar('\n');
434                                }
435                        }
436                }
437        } else
438                getch(ch, parens);
439        if ((*ch) == ':') {
440                do {
441                        getch(ch, parens);
442                } while ((*ch) != ',' && (*ch) != ')' && (*ch) != '[' && (*ch) != ';');
443        }
444        if ((*ch) == '[' && !anerror) {
445                if (!eoln(infile)) {
446                        fscanf(infile, "%lf", &trweight);
447                        getch(ch, parens);
448                        if (*ch != ']') {
449                                printf("ERROR: MISSING RIGHT SQUARE BRACKET\n");
450                                anerror = true;
451                        } else {
452                                getch(ch, parens);
453                                if (*ch != ';') {
454                                        printf("ERROR: MISSING SEMICOLON AFTER SQUARE BRACKETS\n");
455                                        anerror = true;
456                                }
457                        }
458                }
459        } else if (*ch == ';' && !anerror) {
460                trweight = 1.0;
461                if (!eoln(infile))
462                        printf("WARNING: TREE WEIGHT SET TO 1.0\n");
463        }
464        if (*p != NULL && !anerror)
465                (*p)->back = q;
466}                               /* addelement */
467
468
469void 
470initreenode(p)
471        node           *p;
472{
473        /* traverse tree and assign tips to treenode */
474        node           *q;
475
476        if (p) {
477                q = p->next;
478                while (q && q != p) {
479                        initreenode(q->back);
480                        q = q->next;
481                }
482                if (p->tip) {
483                        treenode[(int) p->index - 1] = p;
484                        memcpy(nayme[(int) p->index - 1], p->naime, nmlngth);
485                }
486        }
487}                               /* initreenode */
488
489
490Static Void
491treeread()
492{
493        /* read in user-defined tree and set it up */
494        char            ch;
495        short           parens, nextnode;
496
497        goteof = false;
498        parens = 0;
499        getch(&ch, &parens);
500        while (eoln(infile) && !eof(infile)) {
501                fscanf(infile, "%*[^\n]");
502                getc(infile);
503        }
504        if (eof(infile)) {
505                goteof = true;
506                return;
507        }
508        addelement(&root, NULL, &ch, &parens, &nextnode);
509        fscanf(infile, "%*[^\n]");
510        getc(infile);
511        if (parens != 0 && !anerror) {
512                printf("\nERROR IN TREE FILE:  UNMATCHED PARENTHESES\n");
513                anerror = true;
514        }
515        if (outgrno > spp && !anerror) {
516                anerror = true;
517                printf("ERROR IN OUTGROUP OPTION: SPECIES NUMBER%3hd DOES NOT EXIST\n",
518                       outgrno);
519        }
520}                               /* treeread */
521
522
523Static Void
524reroot(outgroup)
525        node           *outgroup;
526{
527        /* reorients tree, putting outgroup in desired position. */
528        short           i;
529        boolean         nroot;
530        node           *p, *q;
531
532        nroot = false;
533        p = root->next;
534        while (p != root) {
535                if (outgroup->back == p) {
536                        nroot = true;
537                        p = root;
538                } else
539                        p = p->next;
540        }
541        if (nroot)
542                return;
543        p = root;
544        i = 0;
545        while (p->next != root) {
546                p = p->next;
547                i++;
548        }
549        if (i == 2) {
550                root->next->back->back = p->back;
551                p->back->back = root->next->back;
552                q = root->next;
553        } else {
554                p->next = root->next;
555                gnu(&root->next);
556                q = root->next;
557                gnu(&q->next);
558                p = q->next;
559                p->next = root;
560                q->tip = false;
561                p->tip = false;
562        }
563        q->back = outgroup;
564        p->back = outgroup->back;
565        outgroup->back->back = p;
566        outgroup->back = q;
567}                               /* reroot */
568
569void 
570rehash()
571{
572        long           *s, k;
573        short           i, j;
574        double          temp, ss;
575        boolean         done;
576
577        s = (long *) Malloc(setsz * sizeof(long));
578        for (i = 0; i < maxgrp / 2; i++) {
579                k = *order[i];
580                memcpy(s, grouping[k], setsz * sizeof(long));
581                ss = 0.0;
582                for (j = 0; j < setsz; j++)
583                        ss += s[j] /* pow(2, SETBITS*j) */ ;
584                temp = ss * ((sqrt(5.0) - 1) / 2);
585                j = (long) (maxgrp * (temp - floor(temp)));
586                done = false;
587                while (!done) {
588                        if (!grping2[j]) {
589                                grping2[j] = (long *) Malloc(setsz * sizeof(long));
590                                order2[i] = (long *) Malloc(sizeof(long));
591                                tmseen2[j] = (double *) Malloc(sizeof(double));
592                                memcpy(grping2[j], grouping[k], setsz * sizeof(long));
593                                *tmseen2[j] = *timesseen[k];
594                                *order2[i] = j;
595                                grouping[k] = NULL;
596                                timesseen[k] = NULL;
597                                order[i] = NULL;
598                                done = true;
599                        } else {
600                                j++;
601                                if (j >= maxgrp)
602                                        j -= maxgrp;
603                        }
604                }
605        }
606        free(s);
607}                               /* rehash */
608
609
610void 
611enterset(r)
612        node           *r;
613{
614        long           *s, i, j, start;
615        double          times, ss, n;
616        boolean         done, same;
617
618        s = (long *) Malloc(setsz * sizeof(long));
619        memcpy(s, r->nodeset, setsz * sizeof(long));
620        same = true;
621        for (i = 0; i < setsz; i++)
622                if (s[i] != fullset[i])
623                        same = false;
624        if (same)
625                return;
626        times = trweight;
627        ss = 0.0;
628        n = ((sqrt(5.0) - 1) / 2);
629        for (i = 0; i < setsz; i++)
630                ss += s[i] * n;
631        i = (long) (maxgrp * (ss - floor(ss))) + 1;
632        start = i;
633        done = false;
634        while (!done) {
635                if (grouping[i - 1]) {
636                        same = true;
637                        for (j = 0; j < setsz; j++) {
638                                if (s[j] != grouping[i - 1][j])
639                                        same = false;
640                        }
641                }
642                if (grouping[i - 1] && same) {
643                        *timesseen[i - 1] += times;
644                        done = true;
645                } else if (!grouping[i - 1]) {
646                        grouping[i - 1] = (long *) Malloc(setsz * sizeof(long));
647                        lasti++;
648                        order[lasti] = (long *) Malloc(sizeof(long));
649                        timesseen[i - 1] = (double *) Malloc(sizeof(double));
650                        memcpy(grouping[i - 1], s, setsz * sizeof(long));
651                        *timesseen[i - 1] = times;
652                        *order[lasti] = i - 1;
653                        done = true;
654                } else {
655                        i++;
656                        if (i > maxgrp)
657                                i -= maxgrp;
658                }
659                if (!done && i == start) {
660                        maxgrp = maxgrp * 2;
661                        tmseen2 = (double **) Malloc(maxgrp * sizeof(double *));
662                        grping2 = (long **) Malloc(maxgrp * sizeof(long *));
663                        order2 = (long **) Malloc(maxgrp * sizeof(long *));
664                        rehash();
665                        free(timesseen);
666                        free(grouping);
667                        free(order);
668                        timesseen = tmseen2;
669                        grouping = grping2;
670                        order = order2;
671                        done = true;
672                        lasti = maxgrp / 2 - 1;
673                        enterset(r);
674                }
675        }
676        free(s);
677}                               /* enterset */
678
679
680Static Void
681accumulate(r_)
682        node           *r_;
683{
684        node           *r;
685        node           *q;
686        short           i;
687        boolean         done;
688
689        r = r_;
690        if (r->tip) {
691                if (!r->nodeset)
692                        r->nodeset = (long *) Malloc(setsz * sizeof(long));
693                for (i = 0; i < setsz; i++)
694                        r->nodeset[i] = 0L;
695                done = false;
696                i = 0;
697                while (!done && i < setsz) {
698                        if ((int) r->index < (i + 1) * SETBITS) {
699                                r->nodeset[i] = 1L << ((int) r->index - i * SETBITS);
700                                done = true;
701                        } else
702                                i++;
703                }
704        } else {
705                q = r->next;
706                while (q != r) {
707                        accumulate(q->back);
708                        q = q->next;
709                }
710                q = r->next;
711                if (!r->nodeset)
712                        r->nodeset = (long *) Malloc(setsz * sizeof(long));
713                for (i = 0; i < setsz; i++)
714                        r->nodeset[i] = 0;
715                while (q != r) {
716                        for (i = 0; i < setsz; i++)
717                                r->nodeset[i] |= q->back->nodeset[i];
718                        q = q->next;
719                }
720        }
721        if (!r->tip) {
722                if (r->next->next != r)
723                        enterset(r);
724        } else
725                enterset(r);
726}                               /* accumulate */
727
728#define down            2
729#define over            5
730
731
732void 
733compress(n)
734        short          *n;
735{
736        /* push all the nonempty subsets to the front end of their array */
737        short           i, j;
738
739        i = 1;
740        j = 1;
741        do {
742                while (grouping[i - 1])
743                        i++;
744                if (j <= i)
745                        j = i + 1;
746                while (!grouping[j - 1] && j < maxgrp)
747                        j++;
748                if (j < maxgrp) {
749                        grouping[i - 1] = (long *) Malloc(setsz * sizeof(long));
750                        timesseen[i - 1] = (double *) Malloc(sizeof(double));
751                        memcpy(grouping[i - 1], grouping[j - 1], setsz * sizeof(long));
752                        *timesseen[i - 1] = *timesseen[j - 1];
753                        grouping[j - 1] = NULL;
754                        timesseen[j - 1] = NULL;
755                }
756        } while (j != maxgrp);
757        (*n) = i - 1;
758}                               /* compress */
759
760
761void 
762sort(n)
763        short           n;
764{
765        /* Shell sort keeping grouping, timesseen in same order */
766        short           gap, i, j;
767        double          rtemp;
768        long           *stemp;
769
770        gap = n / 2;
771        stemp = (long *) Malloc(setsz * sizeof(long));
772        while (gap > 0) {
773                for (i = gap + 1; i <= n; i++) {
774                        j = i - gap;
775                        while (j > 0) {
776                                if (*timesseen[j - 1] < *timesseen[j + gap - 1]) {
777                                        memcpy(stemp, grouping[j - 1], setsz * sizeof(long));
778                                        memcpy(grouping[j - 1], grouping[j + gap - 1], setsz * sizeof(long));
779                                        memcpy(grouping[j + gap - 1], stemp, setsz * sizeof(long));
780                                        rtemp = *timesseen[j - 1];
781                                        *timesseen[j - 1] = *timesseen[j + gap - 1];
782                                        *timesseen[j + gap - 1] = rtemp;
783                                }
784                                j -= gap;
785                        }
786                }
787                gap /= 2;
788        }
789        free(stemp);
790}                               /* sort */
791
792
793void 
794eliminate(n, n2, fullset)
795        short          *n, *n2;
796        long           *fullset;
797{
798        /* eliminate groups incompatible with preceding ones */
799        short           i, j, k;
800        boolean         comp;
801
802        for (i = 2; i <= (*n); i++) {
803                for (j = 0; j <= i - 2; j++) {
804                        if (timesseen[j] && *timesseen[j] > 0) {
805                                comp = true;
806                                for (k = 0; k < setsz; k++)
807                                        if ((grouping[i - 1][k] & grouping[j][k]) != 0)
808                                                comp = false;
809                                if (!comp) {
810                                        comp = true;
811                                        for (k = 0; k < setsz; k++)
812                                                if ((grouping[i - 1][k] & ~grouping[j][k]) != 0)
813                                                        comp = false;
814                                        if (!comp) {
815                                                comp = true;
816                                                for (k = 0; k < setsz; k++)
817                                                        if ((grouping[j][k] & ~grouping[i - 1][k]) != 0)
818                                                                comp = false;
819                                                if (!comp) {
820                                                        comp = noroot;
821                                                        if (comp) {
822                                                                for (k = 0; k < setsz; k++)
823                                                                        if ((fullset[k] & ~grouping[i - 1][k] & ~grouping[j][k]) != 0)
824                                                                                comp = false;
825                                                        }
826                                                }
827                                        }
828                                }
829                                if (!comp) {
830                                        (*n2)++;
831                                        times2[(*n2) - 1] = (double *) Malloc(sizeof(double));
832                                        group2[(*n2) - 1] = (long *) Malloc(setsz * sizeof(long));
833                                        *times2[(*n2) - 1] = *timesseen[i - 1];
834                                        memcpy(group2[(*n2) - 1], grouping[i - 1], setsz * sizeof(long));
835                                        *timesseen[i - 1] = 0.0;
836                                        for (k = 0; k < setsz; k++)
837                                                grouping[i - 1][k] = 0;
838                                }
839                        }
840                }
841                if (*timesseen[i - 1] == 0.0) {
842                        timesseen[i - 1] = NULL;
843                        grouping[i - 1] = NULL;
844                }
845        }
846}                               /* eliminate */
847
848
849void 
850printset(n)
851        short           n;
852{
853        /* print out a set of species */
854        short           i, j, k, size;
855
856        fprintf(outfile, "\nSet (species in order)   ");
857        for (i = 1; i <= spp - 25; i++)
858                putc(' ', outfile);
859        fprintf(outfile, "  How many times out of%7.2f\n\n", ntrees);
860        for (i = 0; i < n; i++) {
861                if (*timesseen[i] > 0) {
862                        size = 0;
863                        k = 0;
864                        for (j = 1; j <= spp; j++) {
865                                if (j == (k + 1) * SETBITS)
866                                        k++;
867                                if ((1L << (j - k * SETBITS) & grouping[i][k]) != 0)
868                                        size++;
869                        }
870                        if (size != 1 && !(noroot && size == spp - 1)) {
871                                k = 0;
872                                for (j = 1; j <= spp; j++) {
873                                        if (j == (k + 1) * SETBITS)
874                                                k++;
875                                        if ((1L << (j - k * SETBITS) & grouping[i][k]) != 0)
876                                                putc('*', outfile);
877                                        else
878                                                putc('.', outfile);
879                                        if (j % 10 == 0)
880                                                putc(' ', outfile);
881                                }
882                                for (j = 1; j <= 23 - spp; j++)
883                                        putc(' ', outfile);
884                                fprintf(outfile, "    %5.2f\n", *timesseen[i]);
885                        }
886                }
887        }
888}                               /* printset */
889
890
891void 
892bigsubset(st, n)
893        long           *st;
894        short           n;
895{
896        /* find a maximal subset of st among the groupings */
897        short           i, j, k;
898        long           *su;
899        boolean         max, same;
900
901        su = (long *) Malloc(setsz * sizeof(long));
902        for (i = 0; i < setsz; i++)
903                su[i] = 0;
904        for (i = 0; i < n; i++) {
905                max = true;
906                for (j = 0; j < setsz; j++)
907                        if ((grouping[i][j] & ~st[j]) != 0)
908                                max = false;
909                if (max) {
910                        same = true;
911                        for (j = 0; j < setsz; j++)
912                                if (grouping[i][j] != st[j])
913                                        same = false;
914                        max = !same;
915                }
916                if (max) {
917                        for (j = 0; j < setsz; j++)
918                                if ((su[j] & ~grouping[i][j]) != 0)
919                                        max = false;
920                        if (max) {
921                                same = true;
922                                for (j = 0; j < setsz; j++)
923                                        if (su[j] != grouping[i][j])
924                                                same = false;
925                                max = !same;
926                        }
927                        if (max)
928                                memcpy(su, grouping[i], setsz * sizeof(long));
929                }
930        }
931        memcpy(st, su, setsz * sizeof(long));
932        free(su);
933}                               /* bigsubset */
934
935
936void 
937recontraverse(p, st, n, nextnode)
938        node          **p;
939        long           *st;
940        short           n;
941        short          *nextnode;
942{
943        /* traverse to add next node of consensus tree */
944        short           i, j, k, l;
945        boolean         found, same, zero, zero2;
946        long           *tempset, *st2;
947        node           *q;
948
949        k = l = 0;
950        for (i = 1; i <= spp; i++) {
951                if (i == (l + 1) * SETBITS)
952                        l++;
953                if (((1L << (i - l * SETBITS)) & st[l]) != 0) {
954                        k++;
955                        j = i;
956                }
957        }
958        if (k == 1) {
959                *p = treenode[j - 1];
960                return;
961        }
962        gnu(p);
963        (*nextnode)++;
964        (*p)->index = 0.0;
965        for (i = 0; i < n; i++) {
966                same = true;
967                for (j = 0; j < setsz; j++)
968                        if (grouping[i][j] != st[j])
969                                same = false;
970                if (same)
971                        (*p)->index = *timesseen[i];
972        }
973        tempset = (long *) Malloc(setsz * sizeof(long));
974        memcpy(tempset, st, setsz * sizeof(long));
975        q = *p;
976        st2 = (long *) Malloc(setsz * sizeof(long));
977        memcpy(st2, st, setsz * sizeof(long));
978        zero = true;
979        for (j = 0; j < setsz; j++)
980                if (tempset[j] != 0)
981                        zero = false;
982        if (!zero)
983                bigsubset(tempset, n);
984        zero = zero2 = false;
985        while (!zero && !zero2) {
986                zero = zero2 = true;
987                for (j = 0; j < setsz; j++) {
988                        if (st2[j] != 0)
989                                zero = false;
990                        if (tempset[j] != 0)
991                                zero2 = false;
992                }
993                if (!zero && !zero2) {
994                        gnu(&q->next);
995                        q = q->next;
996                        recontraverse(&q->back, tempset, n, nextnode);
997                        q->back->back = q;
998                        for (j = 0; j < setsz; j++)
999                                st2[j] &= ~tempset[j];
1000                        memcpy(tempset, st2, setsz * sizeof(long));
1001                        found = false;
1002                        i = 1;
1003                        while (!found && i <= n) {
1004                                if (grouping[i - 1]) {
1005                                        same = true;
1006                                        for (j = 0; j < setsz; j++)
1007                                                if (grouping[i - 1][j] != tempset[j])
1008                                                        same = false;
1009                                }
1010                                if (grouping[i - 1] && same)
1011                                        found = true;
1012                                else
1013                                        i++;
1014                        }
1015                        zero = true;
1016                        for (j = 0; j < setsz; j++)
1017                                if (tempset[j] != 0)
1018                                        zero = false;
1019                        if (!zero && !found)
1020                                bigsubset(tempset, n);
1021                }
1022        }
1023        q->next = *p;
1024        free(tempset);
1025        free(st2);
1026}                               /* recontraverse */
1027
1028
1029void 
1030reconstruct(n)
1031        short           n;
1032{
1033        /* reconstruct tree from the subsets */
1034        short           nextnode, i;
1035        long           *s;
1036
1037        nextnode = spp + 1;
1038        s = (long *) Malloc(setsz * sizeof(long));
1039        memcpy(s, fullset, setsz * sizeof(long));
1040        recontraverse(&root, s, n, &nextnode);
1041        free(s);
1042}                               /* reconstruct */
1043
1044
1045void 
1046coordinates(p, tipy)
1047        node           *p;
1048        short          *tipy;
1049{
1050        /* establishes coordinates of nodes */
1051        node           *q, *first, *last;
1052        short           maxx;
1053
1054        if (p->tip) {
1055                p->xcoord = 0;
1056                p->ycoord = *tipy;
1057                p->ymin = *tipy;
1058                p->ymax = *tipy;
1059                (*tipy) += down;
1060                return;
1061        }
1062        q = p->next;
1063        maxx = 0;
1064        while (q != p) {
1065                coordinates(q->back, tipy);
1066                if (!q->back->tip) {
1067                        if (q->back->xcoord > maxx)
1068                                maxx = q->back->xcoord;
1069                }
1070                q = q->next;
1071        }
1072        first = p->next->back;
1073        q = p;
1074        while (q->next != p)
1075                q = q->next;
1076        last = q->back;
1077        p->xcoord = maxx + over;
1078        p->ycoord = (first->ycoord + last->ycoord) / 2;
1079        p->ymin = first->ymin;
1080        p->ymax = last->ymax;
1081}                               /* coordinates */
1082
1083
1084void 
1085drawline(i)
1086        short           i;
1087{
1088        /* draws one row of the tree diagram by moving up tree */
1089        node           *p, *q;
1090        short           n, j;
1091        boolean         extra, done, trif;
1092        node           *r, *first, *last;
1093        boolean         found;
1094
1095        p = root;
1096        q = root;
1097        fprintf(outfile, "  ");
1098        extra = false;
1099        trif = false;
1100        do {
1101                if (!p->tip) {
1102                        found = false;
1103                        r = p->next;
1104                        while (r != p && !found) {
1105                                if (i >= r->back->ymin && i <= r->back->ymax) {
1106                                        q = r->back;
1107                                        found = true;
1108                                } else
1109                                        r = r->next;
1110                        }
1111                        first = p->next->back;
1112                        r = p;
1113                        while (r->next != p)
1114                                r = r->next;
1115                        last = r->back;
1116                }
1117                done = (p->tip || p == q);
1118                n = p->xcoord - q->xcoord;
1119                if (extra) {
1120                        n--;
1121                        extra = false;
1122                }
1123                if (q->ycoord == i && !done) {
1124                        if (trif)
1125                                putc('-', outfile);
1126                        else
1127                                putc('+', outfile);
1128                        trif = false;
1129                        if (!q->tip) {
1130                                for (j = 1; j <= n - 5; j++)
1131                                        putc('-', outfile);
1132                                if (q->index >= 100.0)
1133                                        fprintf(outfile, "%5.1f", q->index);
1134                                else if (q->index >= 10.0)
1135                                        fprintf(outfile, "-%4.1f", q->index);
1136                                else
1137                                        fprintf(outfile, "--%3.1f", q->index);
1138                                extra = true;
1139                                trif = true;
1140                        } else {
1141                                for (j = 1; j < n; j++)
1142                                        putc('-', outfile);
1143                        }
1144                } else if (!p->tip && last->ycoord > i && first->ycoord < i &&
1145                           (i != p->ycoord || p == root)) {
1146                        putc('!', outfile);
1147                        for (j = 1; j < n; j++)
1148                                putc(' ', outfile);
1149                } else {
1150                        for (j = 1; j <= n; j++)
1151                                putc(' ', outfile);
1152                        if (trif)
1153                                trif = false;
1154                }
1155                if (q != p)
1156                        p = q;
1157        } while (!done);
1158        if (p->ycoord == i && p->tip) {
1159                for (j = 0; j < nmlngth; j++)
1160                        putc(nayme[(int) ((long) p->index) - 1][j], outfile);
1161        }
1162        putc('\n', outfile);
1163}                               /* drawline */
1164
1165
1166void 
1167printree()
1168{
1169        /* prints out diagram of the tree */
1170        short           i;
1171        short           tipy;
1172
1173        if (treeprint) {
1174                fprintf(outfile, "\nCONSENSUS TREE:\n");
1175                fprintf(outfile, "the numbers at the forks indicate the number\n");
1176                fprintf(outfile, "of times the group consisting of the species\n");
1177                fprintf(outfile, "which are to the right of that fork occurred\n");
1178                fprintf(outfile, "among the trees, out of %6.2f trees\n", ntrees);
1179                tipy = 1;
1180                coordinates(root, &tipy);
1181                putc('\n', outfile);
1182                for (i = 1; i <= tipy - down; i++)
1183                        drawline(i);
1184                putc('\n', outfile);
1185        }
1186        if (noroot) {
1187                fprintf(outfile, "\n  remember:");
1188                if (didreroot)
1189                        fprintf(outfile, " (though rerooted by outgroup)");
1190                fprintf(outfile, " this is an unrooted tree!\n");
1191        }
1192        putc('\n', outfile);
1193}                               /* printree */
1194
1195
1196#undef down
1197#undef over
1198
1199void 
1200treeout(p)
1201        node           *p;
1202{
1203        /* write out file with representation of final tree */
1204        short           i, n;
1205        double          x;
1206        Char            c;
1207        node           *q;
1208
1209        if (p->tip) {
1210                for (i = 1; i <= nmlngth; i++) {
1211                        if (nayme[(int) ((long) p->index) - 1][i - 1] != ' ')
1212                                n = i;
1213                }
1214                for (i = 0; i < n; i++) {
1215                        c = nayme[(int) ((long) p->index) - 1][i];
1216                        if (c == ' ')
1217                                c = '_';
1218                        putc(c, treefile);
1219                }
1220                col += n;
1221        } else {
1222                putc('(', treefile);
1223                col++;
1224                q = p->next;
1225                while (q != p) {
1226                        treeout(q->back);
1227                        q = q->next;
1228                        if (q == p)
1229                                break;
1230                        putc(',', treefile);
1231                        col++;
1232                        if (col > 60) {
1233                                putc('\n', treefile);
1234                                col = 0;
1235                        }
1236                }
1237                putc(')', treefile);
1238                col++;
1239        }
1240        if (p->tip)
1241                x = ntrees;
1242        else
1243                x = p->index;
1244        if (p == root) {
1245                fprintf(treefile, ";\n");
1246                return;
1247        }
1248        fprintf(treefile, ":%1.3f", x/(double)ntrees);
1249        col += 4;
1250}                               /* treeout */
1251
1252
1253void 
1254consensus()
1255{
1256        short           i;
1257        short           n, n2;
1258        short           tipy;
1259
1260        n2 = 0;
1261        compress(&n);
1262        sort(n);
1263        eliminate(&n, &n2, fullset);
1264        compress(&n);
1265        reconstruct(n);
1266        tipy = 1;
1267        coordinates(root, &tipy);
1268        if (prntsets) {
1269                fprintf(outfile, "\nSets included in the consensus tree\n");
1270                printset(n);
1271                for (i = 0; i < n2; i++) {
1272                        if (!grouping[i]) {
1273                                grouping[i] = (long *) Malloc(setsz * sizeof(long));
1274                                timesseen[i] = (double *) Malloc(sizeof(double));
1275                        }
1276                        memcpy(grouping[i], group2[i], setsz * sizeof(long));
1277                        *timesseen[i] = *times2[i];
1278                }
1279                n = n2;
1280                fprintf(outfile, "\n\nSets NOT included in consensus tree:");
1281                if (n2 == 0)
1282                        fprintf(outfile, " NONE\n");
1283                else {
1284                        putc('\n', outfile);
1285                        printset(n);
1286                }
1287                putc('\n', outfile);
1288        }
1289        printree();
1290        if (progress)
1291                printf("\nOutput written to output file\n\n");
1292        if (trout) {
1293                treeout(root);
1294                if (progress)
1295                        printf("Tree also written onto file\n\n");
1296        }
1297}                               /* consensus */
1298
1299
1300main(argc, argv)
1301        int             argc;
1302        Char           *argv[];
1303{
1304        char            infilename[100], outfilename[100], trfilename[100];
1305#ifdef MAC
1306        macsetup("Consense", "");
1307#endif
1308        openfile(&infile, INFILE, "r", argv[0], infilename);
1309        openfile(&outfile, OUTFILE, "w", argv[0], outfilename);
1310
1311        ibmpc = ibmpc0;
1312        ansi = ansi0;
1313        vt52 = vt520;
1314        didreroot = false;
1315        first = true;
1316        spp = 0;
1317        garbage = NULL;
1318        col = 0;
1319        getoptions();
1320        if (trout)
1321                openfile(&treefile, TREEFILE, "w", argv[0], trfilename);
1322        if (prntsets)
1323                fprintf(outfile, "Species in order: \n\n");
1324        i = 1;
1325        ntrees = 0.0;
1326        maxgrp = 10000;
1327        lasti = -1;
1328        grouping = (long **) Malloc(maxgrp * sizeof(long *));
1329        order = (long **) Malloc(maxgrp * sizeof(long *));
1330        timesseen = (double **) Malloc(maxgrp * sizeof(double *));
1331        while (!anerror && !eof(infile)) {
1332                treeread();
1333                if (first) {
1334                        treenode = (pointarray) Malloc(spp * sizeof(pointarray *));
1335                        nayme = (naym *) Malloc(spp * sizeof(naym));
1336                        initreenode(root);
1337                        setsz = (short) ceil(((double) spp + 1.0) / (double) SETBITS);
1338                        fullset = (long *) Malloc(setsz * sizeof(long));
1339                        for (j = 0; j < setsz; j++)
1340                                fullset[j] = 0L;
1341                        for (j = 0; j < setsz; j++) {
1342                                if (spp + 1 < (j + 1) * SETBITS)
1343                                        fullset[j] = (1L << ((spp + 1) - j * SETBITS)) - 1;
1344                                else
1345                                        fullset[j] = ~0L;
1346                        }
1347                        fullset[0] -= 1;
1348                }
1349                if (goteof)
1350                        continue;
1351                ntrees += trweight;
1352                if (noroot) {
1353                        reroot(treenode[outgrno - 1]);
1354                        didreroot = outgropt;
1355                }
1356                accumulate(root);
1357                first = false;
1358                gdispose(root);
1359                i++;
1360        }
1361        putc('\n', outfile);
1362        group2 = (long **) Malloc(maxgrp * sizeof(long *));
1363        times2 = (double **) Malloc(maxgrp * sizeof(double *));
1364        consensus();
1365
1366        FClose(treefile);
1367        FClose(infile);
1368        FClose(outfile);
1369#ifdef MAC
1370        fixmacfile(outfilename);
1371        fixmacfile(trfilename);
1372#endif
1373        exit(0);
1374}                               /* main */
1375
1376
1377int 
1378eof(f)
1379        FILE           *f;
1380{
1381        register int    ch;
1382
1383        if (feof(f))
1384                return 1;
1385        if (f == stdin)
1386                return 0;
1387        ch = getc(f);
1388        if (ch == EOF)
1389                return 1;
1390        ungetc(ch, f);
1391        return 0;
1392}
1393
1394
1395int 
1396eoln(f)
1397        FILE           *f;
1398{
1399        register int    ch;
1400
1401        ch = getc(f);
1402        if (ch == EOF)
1403                return 1;
1404        ungetc(ch, f);
1405        return (ch == '\n');
1406}
1407
1408void 
1409memerror()
1410{
1411        printf("Error allocating memory\n");
1412        exit(-1);
1413}
1414
1415MALLOCRETURN   *
1416mymalloc(x)
1417        long            x;
1418{
1419        MALLOCRETURN   *mem;
1420        mem = (MALLOCRETURN *) calloc(1, x);
1421        if (!mem)
1422                memerror();
1423        else
1424                return (MALLOCRETURN *) mem;
1425}
Note: See TracBrowser for help on using the repository browser.