source: tags/arb_5.0/TOOLS/arb_dnarates.c

Last change on this file was 5823, checked in by westram, 15 years ago
  • more uniform function names
    • EXP_create_experiment_rel_experiment_data → EXP_find_or_create_experiment_rel_exp_data
    • EXP_find_experiment_rel_experiment_data → EXP_find_experiment_rel_exp_data
    • GBT_create_SAI → GBT_find_or_create_SAI
    • GBT_create_item → GBT_find_or_create_item_rel_item_data
    • GBT_create_species → GBT_find_or_create_species
    • GBT_create_species_rel_species_data → GBT_find_or_create_species_rel_species_data
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 55.9 KB
Line 
1#define programName     "DNAml_rates"
2#define programVersion  "1.0.0"
3#define programDate     "April 11, 1992"
4
5/*  Maximum likelihood site rate calculation, Gary Olsen, 1991, 1992.
6 *
7 *  Portions based on the program dnaml version 3.3 by Joseph Felsenstein
8 *
9 *  Copyright notice from dnaml:
10 *
11 *  version 3.3. (c) Copyright 1986, 1990 by the University of Washington
12 *  and Joseph Felsenstein.  Written by Joseph Felsenstein.  Permission is
13 *  granted to copy and use this program provided no fee is charged for it
14 *  and provided that this copyright notice is not removed.
15 */
16
17/* Conversion to C by Gary Olsen, 1991 */
18
19#include <stdio.h>
20#include <stdlib.h>
21#include <unistd.h>
22/* #include <malloc.h> */
23#include <string.h>
24#include <math.h>
25#include "DNAml_rates_1_0.h"
26
27#include <aw_awars.hxx>
28#include <arbdb.h>
29#include <arbdbt.h>
30
31#ifndef ARB_ASSERT_H
32#include <arb_assert.h>
33#endif
34#define assert(bed) arb_assert(bed)
35
36
37/*  Global variables */
38
39xarray
40*usedxtip, *freextip;
41
42double
43dLidki[maxpatterns],      /*  change in pattern i likelihood with rate */
44    bestki[maxpatterns],      /*  best rate for pattern i */
45    bestLi[maxpatterns],      /*  best likelihood found for pattern i */
46    patrate[maxpatterns];     /*  current rate for pattern i */
47
48double
49xi, xv, ttratio,          /*  transition/transversion info */
50    freqa, freqc, freqg, freqt,  /*  base frequencies */
51    freqr, freqy, invfreqr, invfreqy,
52    freqar, freqcy, freqgr, freqty,
53    fracchange;    /*  random matching fraquency (in a sense) */
54
55int
56info[maxsites+1],         /*  number of informative nucleotides */
57    patsite[maxsites+1],      /*  site corresponding to pattern */
58    pattern[maxsites+1],      /*  pattern number corresponding to site */
59    patweight[maxsites+1],    /*  weight of given pattern */
60    weight[maxsites+1];       /*  weight of sequence site */
61
62int
63categs,        /*  number of rate categories */
64    endsite,       /*  number of unique sequence patterns */
65    mininfo,       /*  minimum number of informative sequences for rate est */
66    numsp,         /*  number of species */
67    sites,         /*  number of input sequence positions */
68    weightsum;     /*  sum of weights of positions in analysis */
69
70boolean
71anerror,       /*  error flag */
72    freqsfrom,     /*  use empirical base frequencies */
73    interleaved,   /*  input data are in interleaved format */
74    printdata,     /*  echo data to output stream */
75    writefile,     /*  put weight and rate data in file */
76    userweights;   /*  use user-supplied position weight mask */
77
78char
79*y[maxsp+1];                    /*  sequence data array */
80
81
82#if UseStdin       /*  Use standard input */
83#  define  INFILE  stdin
84#else              /*  Use file named "infile" */
85#  define  INFILE  fp
86FILE  *fp;
87#endif
88
89#if DebugData
90FILE *debug;
91#endif
92
93
94/*=======================================================================*/
95/*                              PROGRAM                                  */
96/*=======================================================================*/
97
98#if 0
99void hang(msg) char *msg; {printf("Hanging around: %s\n", msg); while(1);}
100#endif
101
102void getnums ()         /* input number of species, number of sites */
103{ /* getnums */
104    printf("\n%s, version %s, %s\n\n",
105           programName,
106           programVersion,
107           programDate);
108    printf("Portions based on Joseph Felsenstein's Nucleic acid sequence\n");
109    printf("Maximum Likelihood method, version 3.3\n\n");
110
111    if (fscanf(INFILE, "%d %d", &numsp, &sites) != 2) {
112        printf("ERROR: Problem reading number of species and sites\n");
113        anerror = TRUE;
114        return;
115    }
116    printf("%d Species, %d Sites\n\n", numsp, sites);
117
118    if (numsp > maxsp) {
119        printf("ERROR: Too many species; adjust program constants\n");
120        anerror = TRUE;
121    }
122    else if (numsp < 4) {
123        printf("ERROR: Too few species\n");
124        anerror = TRUE;
125    }
126
127    if (sites > maxsites) {
128        printf("ERROR: Too many sites; adjust program constants\n");
129        anerror = TRUE;
130    }
131    else if (sites < 1) {
132        printf("ERROR: Too few sites\n");
133        anerror = TRUE;
134    }
135} /* getnums */
136
137
138boolean digit (ch) int ch; {return (ch >= '0' && ch <= '9'); }
139
140
141boolean white (ch) int ch; { return (ch == ' ' || ch == '\n' || ch == '\t'); }
142
143
144void uppercase (chptr)
145     int *chptr;
146     /* convert character to upper case -- either ASCII or EBCDIC */
147{ /* uppercase */
148    int  ch;
149    ch = *chptr;
150    if ((ch >= 'a' && ch <= 'i') || (ch >= 'j' && ch <= 'r')
151        || (ch >= 's' && ch <= 'z'))
152        *chptr = ch + 'A' - 'a';
153} /* uppercase */
154
155
156int base36 (ch)
157     int ch;
158{ /* base36 */
159    if      (ch >= '0' && ch <= '9') return (ch - '0');
160    else if (ch >= 'A' && ch <= 'I') return (ch - 'A' + 10);
161    else if (ch >= 'J' && ch <= 'R') return (ch - 'J' + 19);
162    else if (ch >= 'S' && ch <= 'Z') return (ch - 'S' + 28);
163    else if (ch >= 'a' && ch <= 'i') return (ch - 'a' + 10);
164    else if (ch >= 'j' && ch <= 'r') return (ch - 'j' + 19);
165    else if (ch >= 's' && ch <= 'z') return (ch - 's' + 28);
166    else return -1;
167} /* base36 */
168
169
170int itobase36 (i)
171     int  i;
172{ /* itobase36 */
173    if      (i <  0) return '?';
174    else if (i < 10) return (i      + '0');
175    else if (i < 19) return (i - 10 + 'A');
176    else if (i < 28) return (i - 19 + 'J');
177    else if (i < 36) return (i - 28 + 'S');
178    else return '?';
179} /* itobase36 */
180
181
182int findch (c)
183     int  c;
184{
185    int ch;
186
187    while ((ch = getc(INFILE)) != EOF && ch != c) ;
188    return  ch;
189}
190
191
192void inputweights ()
193/* input the character weights 0, 1, 2 ... 9, A, B, ... Y, Z */
194{ /* inputweights */
195    int i, ch, wi;
196
197    for (i = 2; i <= nmlngth; i++)  (void) getc(INFILE);
198    weightsum = 0;
199    i = 1;
200    while (i <= sites) {
201        ch = getc(INFILE);
202        wi = base36(ch);
203        if (wi >= 0)
204            weightsum += weight[i++] = wi;
205        else if (! white(ch)) {
206            printf("ERROR: Bad weight character: '%c'", ch);
207            printf("       Weights must be a digit or a letter.\n");
208            anerror = TRUE;
209            return;
210        }
211    }
212
213    if (findch('\n') == EOF) {      /* skip to end of line */
214        printf("ERROR: Missing newline at end of weight data\n");
215        anerror = TRUE;
216    }
217} /* inputweights */
218
219
220void getoptions ()
221{ /* getoptions */
222    int  ch, i, extranum;
223
224    categs      =     0;  /*  Number of rate categories */
225    freqsfrom   = FALSE;  /*  Use empirical base frequencies */
226    interleaved =  TRUE;  /*  By default, data format is interleaved */
227    mininfo     = MIN_INFO; /*  Default minimum number of informative seqs */
228    printdata   = FALSE;  /*  Don't echo data to output stream */
229    ttratio     =   2.0;  /*  Transition/transversion rate ratio */
230    userweights = FALSE;  /*  User-defined position weights */
231    writefile   = FALSE;  /*  Do not write to file */
232    extranum    =     0;
233
234    while ((ch = getc(INFILE)) != '\n' && ch != EOF) {
235        uppercase (& ch);
236        switch (ch) {
237            case '1' : printdata    = ! printdata; break;
238            case 'C' : categs       = -1; extranum++; break;
239            case 'F' : freqsfrom    = TRUE; break;
240            case 'I' : interleaved  = ! interleaved; break;
241            case 'L' : break;
242            case 'M' : mininfo      = 0; extranum++; break;
243            case 'T' : ttratio      = -1.0; extranum++; break;
244            case 'U' : break;
245            case 'W' : userweights  = TRUE; weightsum = 0; extranum++; break;
246            case 'Y' : writefile    = ! writefile; break;
247            case ' ' : break;
248            case '\t': break;
249            default  :
250                printf("ERROR: Bad option character: '%c'\n", ch);
251                anerror = TRUE;
252                return;
253        }
254    }
255
256    if (ch == EOF) {
257        printf("ERROR: End-of-file in options list\n");
258        anerror = TRUE;
259        return;
260    }
261
262    /*  process lines with auxiliary data */
263
264    while (extranum-- && ! anerror) {
265        ch = getc(INFILE);
266        uppercase (& ch);
267        switch (ch) {
268
269            case 'C':
270                if (categs >= 0) {
271                    printf("ERROR: Unexpected Categories data\n");
272                    anerror = TRUE;
273                }
274                else if (fscanf(INFILE,"%d",&categs) != 1 || findch('\n')==EOF) {
275                    printf("ERROR: Problem reading number of rate categories\n");
276                    anerror = TRUE;
277                }
278                else if (categs < 1 || categs > maxcategories) {
279                    printf("ERROR: Bad number of rate categories: %d\n", categs);
280                    anerror = TRUE;
281                }
282                break;
283
284            case 'M':   /*  Minimum infomative sequences  */
285                if (mininfo > 0) {
286                    printf("ERROR: Unexpected Min informative residues data\n");
287                    anerror = TRUE;
288                }
289                else if (fscanf(INFILE,"%d",&mininfo)!=1 || findch('\n')==EOF) {
290                    printf("ERROR: Problem reading min informative residues\n");
291                    anerror = TRUE;
292                }
293                else if (mininfo < 2 || mininfo > numsp) {
294                    printf("ERROR: Bad number for informative residues: %d\n",
295                           mininfo);
296                    anerror = TRUE;
297                }
298                break;
299
300            case 'T':   /*  Transition/transversion ratio  */
301                if (ttratio >= 0.0) {
302                    printf("ERROR: Unexpected Transition/transversion data\n");
303                    anerror = TRUE;
304                }
305                else if (fscanf(INFILE,"%lf",&ttratio)!=1 || findch('\n')==EOF) {
306                    printf("ERROR: Problem reading transition/transversion data\n");
307                    anerror = TRUE;
308                }
309                break;
310
311            case 'W':    /*  Weights  */
312                if (! userweights || weightsum > 0) {
313                    printf("ERROR: Unexpected Weights data\n");
314                    anerror = TRUE;
315                }
316                else {
317                    inputweights();
318                }
319                break;
320
321            default:
322                printf("ERROR: Auxiliary data line starts with '%c'\n", ch);
323                anerror = TRUE;
324                break;
325        }
326    }
327
328    if (anerror) return;
329
330    if (categs < 0) {
331        printf("ERROR: Category data missing from input\n");
332        anerror = TRUE;
333        return;
334    }
335
336    if (mininfo <= 0) {
337        printf("ERROR: Minimum informative residues missing from input\n");
338        anerror = TRUE;
339        return;
340    }
341    else {
342        printf("There must be at least %d informative residues per column\n\n", mininfo);
343    }
344
345    if (ttratio < 0.0) {
346        printf("ERROR: Transition/transversion data missing from input\n");
347        anerror = TRUE;
348        return;
349    }
350
351    if (! userweights) {
352        for (i = 1; i <= sites; i++)  weight[i] = 1;
353        weightsum = sites;
354    }
355    else if (weightsum < 1) {
356        printf("ERROR: Weight data invalid or missing from input\n");
357        anerror = TRUE;
358        return;
359    }
360
361} /* getoptions */
362
363
364void getbasefreqs ()
365{ /* getbasefreqs */
366    double  suma, sumb;
367
368    if (freqsfrom)  printf("Empirical ");
369    printf("Base Frequencies:\n\n");
370
371    if (! freqsfrom) {
372        if (fscanf(INFILE, "%lf%lf%lf%lf", &freqa, &freqc, &freqg, &freqt) != 4
373            || findch('\n') == EOF) {
374            printf("ERROR: Problem reading user base frequencies\n");
375            anerror = TRUE;
376        }
377    }
378
379    printf("   A    %10.5f\n", freqa);
380    printf("   C    %10.5f\n", freqc);
381    printf("   G    %10.5f\n", freqg);
382    printf("  T(U)  %10.5f\n\n", freqt);
383    freqr = freqa + freqg;
384    invfreqr = 1.0/freqr;
385    freqar = freqa * invfreqr;
386    freqgr = freqg * invfreqr;
387    freqy = freqc + freqt;
388    invfreqy = 1.0/freqy;
389    freqcy = freqc * invfreqy;
390    freqty = freqt * invfreqy;
391    printf("Transition/transversion ratio = %10.6f\n\n", ttratio);
392    suma = ttratio*freqr*freqy - (freqa*freqg + freqc*freqt);
393    sumb = freqa*freqgr + freqc*freqty;
394    xi = suma/(suma+sumb);
395    xv = 1.0 - xi;
396    ttratio = xi / xv;
397    if (xi <= 0.0) {
398        printf("WARNING: This transition/transversion ratio is\n");
399        printf("         impossible with these base frequencies!\n");
400        xi = 3/5;
401        xv = 2/5;
402        printf("Transition/transversion parameter reset\n\n");
403    }
404    printf("(Transition/transversion parameter = %10.6f)\n\n", xi/xv);
405    fracchange = xi*(2*freqa*freqgr + 2*freqc*freqty)
406        + xv*(1.0 - freqa*freqa - freqc*freqc
407              - freqg*freqg - freqt*freqt);
408} /* getbasefreqs */
409
410
411void getyspace ()
412{ /* getyspace */
413    long size;
414    int  i;
415    char *y0;
416
417    size = 4 * (sites/4 + 1);
418    if (! (y0 = malloc((unsigned) (sizeof(char) * size * (numsp+1))))) {
419        printf("ERROR: Unable to obtain space for data array\n");
420        anerror = TRUE;
421    }
422    else {
423        for (i = 0; i <= numsp; i++) {
424            y[i] = y0;
425            y0 += size;
426        }
427    }
428} /* getyspace */
429
430
431void setuptree (tr, numSp)
432     tree  *tr;
433     int    numSp;
434{ /* setuptree */
435    int     i, j;
436    nodeptr p = 0, q;
437
438    for (i = 1; i <= numSp; i++) {   /*  Set-up tips */
439        if ((anerror = !(p = (nodeptr) malloc((unsigned) sizeof(node))))) break;
440        p->x      = (xarray *) NULL;
441        p->tip    = (char *) NULL;
442        p->number = i;
443        p->next   = (node *) NULL;
444        p->back   = (node *) NULL;
445        tr->nodep[i] = p;
446    }
447
448    for (i = numSp+1; i <= 2*numSp-1 && ! anerror; i++) { /* Internal nodes */  /* was : 2*numSp-2 (ralf) */
449        q = (node *) NULL;
450        for (j = 1; j <= 3; j++) {
451            if ((anerror = !(p = (nodeptr) malloc((unsigned) sizeof(node))))) break;
452            p->x      = (xarray *) NULL;
453            p->tip    = (char *) NULL;
454            p->number = i;
455            p->next   = q;
456            p->back   = (node *) NULL;
457            q = p;
458        }
459        if (anerror) break;
460        p->next->next->next = p;
461        tr->nodep[i] = p;
462    }
463
464    tr->likelihood = unlikely;
465    tr->start      = tr->nodep[1];
466    tr->mxtips     = numSp;
467    tr->ntips      = 0;
468    tr->nextnode   = 0;
469    tr->opt_level  = 0;
470    tr->smoothed   = FALSE;
471    if (anerror) printf("ERROR: Unable to obtain sufficient memory");
472} /* setuptree */
473
474
475void freeTreeNode(p)   /*  Free a tree node (sector) */
476     nodeptr  p;
477{ /* freeTreeNode */
478    if (p) {
479        if (p->x) {
480            if (p->x->a) free((char *) p->x->a);
481            free ((char *) p->x);
482        }
483        free ((char *) p);
484    }
485} /* freeTreeNode */
486
487
488void freeTree (tr)
489     tree  *tr;
490{ /* freeTree */
491    int  i;
492    nodeptr  p, q;
493
494    for (i = 1; i <= tr->mxtips; i++) freeTreeNode(tr->nodep[i]);
495
496    for (i = tr->mxtips+1; i <= 2*(tr->mxtips)-2; i++) {
497        if ((p = tr->nodep[i])) {
498            if ((q = p->next)) {
499                freeTreeNode(q->next);
500                freeTreeNode(q);
501            }
502            freeTreeNode(p);
503        }
504    }
505} /* freeTree */
506
507
508void getdata (tr)
509     tree  *tr;
510     /* read sequences */
511{ /* getdata */
512    int  i, j, k, l, basesread, basesnew, ch;
513    int  meaning[256];          /*  meaning of input characters */
514    char *nameptr;
515    boolean  allread, firstpass;
516
517    for (i = 0; i <= 255; i++) meaning[i] = 0;
518    meaning['A'] =  1;
519    meaning['B'] = 14;
520    meaning['C'] =  2;
521    meaning['D'] = 13;
522    meaning['G'] =  4;
523    meaning['H'] = 11;
524    meaning['K'] = 12;
525    meaning['M'] =  3;
526    meaning['N'] = 15;
527    meaning['O'] = 15;
528    meaning['R'] =  5;
529    meaning['S'] =  6;
530    meaning['T'] =  8;
531    meaning['U'] =  8;
532    meaning['V'] =  7;
533    meaning['W'] =  9;
534    meaning['X'] = 15;
535    meaning['Y'] = 10;
536    meaning['?'] = 15;
537    meaning['-'] = 15;
538
539    basesread = basesnew = 0;
540
541    allread = FALSE;
542    firstpass = TRUE;
543    ch = ' ';
544
545    while (! allread) {
546        for (i = 1; i <= numsp; i++) {          /*  Read data line */
547            if (firstpass) {                      /*  Read species names */
548                j = 1;
549                while (white(ch = getc(INFILE))) {  /*  Skip blank lines */
550                    if (ch == '\n')  j = 1;  else  j++;
551                }
552
553                if (j > nmlngth) {
554                    printf("ERROR: Blank name for species %d; ", i);
555                    printf("check number of species,\n");
556                    printf("       number of sites, and interleave option.\n");
557                    anerror = TRUE;
558                    return;
559                }
560
561                nameptr = tr->nodep[i]->name;
562                for (k = 1; k < j; k++)  *nameptr++ = ' ';
563
564                while (ch != '\n' && ch != EOF) {
565                    if (ch == '_' || white(ch))  ch = ' ';
566                    *nameptr++ = ch;
567                    if (++j > nmlngth) break;
568                    ch = getc(INFILE);
569                }
570
571                while (j++ <= nmlngth) *nameptr++ = ' ';
572                *nameptr = '\0';                      /*  add null termination */
573
574                if (ch == EOF) {
575                    printf("ERROR: End-of-file in name of species %d\n", i);
576                    anerror = TRUE;
577                    return;
578                }
579            }    /* if (firstpass) */
580
581            j = basesread;
582            while ((j < sites) && ((ch = getc(INFILE)) != EOF)
583                   && ((! interleaved) || (ch != '\n'))) {
584                uppercase (& ch);
585                if (meaning[ch] || ch == '.') {
586                    j++;
587                    if (ch == '.') {
588                        if (i != 1) ch = y[1][j];
589                        else {
590                            printf("ERROR: Dot (.) found at site %d of sequence 1\n", j);
591                            anerror = TRUE;
592                            return;
593                        }
594                    }
595                    y[i][j] = ch;
596                }
597                else if (white(ch) || digit(ch)) ;
598                else {
599                    printf("ERROR: Bad base (%c) at site %d of sequence %d\n",
600                           ch, j, i);
601                    anerror = TRUE;
602                    return;
603                }
604            }
605
606            if (ch == EOF) {
607                printf("ERROR: End-of-file at site %d of sequence %d\n", j, i);
608                anerror = TRUE;
609                return;
610            }
611
612            if (! firstpass && (j == basesread)) i--;        /* no data on line */
613            else if (i == 1) basesnew = j;
614            else if (j != basesnew) {
615                printf("ERROR: Sequences out of alignment\n");
616                printf("       %d (instead of %d) residues read in sequence %d\n",
617                       j - basesread, basesnew - basesread, i);
618                anerror = TRUE;
619                return;
620            }
621
622            while (ch != '\n' && ch != EOF) ch = getc(INFILE); /* flush line */
623        }                                                  /* next sequence */
624        firstpass = FALSE;
625        basesread = basesnew;
626        allread = (basesread >= sites);
627    }
628
629    /*  Print listing of sequence alignment */
630
631    if (printdata) {
632        j = nmlngth - 5 + ((sites + ((sites-1)/10))/2);
633        if (j < nmlngth - 1) j = nmlngth - 1;
634        if (j > 37) j = 37;
635        printf("Name"); for (i=1;i<=j;i++) putchar(' '); printf("Sequences\n");
636        printf("----"); for (i=1;i<=j;i++) putchar(' '); printf("---------\n");
637        putchar('\n');
638
639        for (i = 1; i <= sites; i += 60) {
640            l = i + 59;
641            if (l > sites) l = sites;
642
643            if (userweights) {
644                printf("Weights   ");
645                for (j = 11; j <= nmlngth+3; j++) putchar(' ');
646                for (k = i; k <= l; k++) {
647                    putchar(itobase36(weight[k]));
648                    if (((k % 10) == 0) && ((k % 60) != 0)) putchar(' ');
649                }
650                putchar('\n');
651            }
652
653            for (j = 1; j <= numsp; j++) {
654                printf("%s   ", tr->nodep[j]->name);
655                for (k = i; k <= l; k++) {
656                    ch = y[j][k];
657                    if ((j > 1) && (ch == y[1][k])) ch = '.';
658                    putchar(ch);
659                    if (((k % 10) == 0) && ((k % 60) != 0)) putchar(' ');
660                }
661                putchar('\n');
662            }
663            putchar('\n');
664        }
665    }
666
667    /*  Convert characters to meanings */
668
669    for (i = 1; i <= sites; i++)  info[i] = 0;
670
671    for (j = 1; j <= numsp; j++) {
672        for (i = 1; i <= sites; i++) {
673            if ((y[j][i] = meaning[(int)y[j][i]]) != 15) info[i]++;
674        }
675    }
676
677    for (i = 1; i <= sites; i++)  if (info[i] < MIN_INFO)  weight[i] = 0;
678
679} /* getdata */
680
681
682void sitesort ()
683/* Shell sort keeping sites, weights in same order */
684{ /* sitesort */
685    int  gap, i, j, jj, jg, k;
686    boolean  flip, tied;
687
688    for (gap = sites/2; gap > 0; gap /= 2) {
689        for (i = gap + 1; i <= sites; i++) {
690            j = i - gap;
691
692            do {
693                jj = patsite[j];
694                jg = patsite[j+gap];
695                flip = FALSE;
696                tied = TRUE;
697                for (k = 1; tied && (k <= numsp); k++) {
698                    flip = (y[k][jj] >  y[k][jg]);
699                    tied = (y[k][jj] == y[k][jg]);
700                }
701                if (flip) {
702                    patsite[j]     = jg;
703                    patsite[j+gap] = jj;
704                    j -= gap;
705                }
706            } while (flip && (j > 0));
707
708        }
709    }
710} /* sitesort */
711
712
713void sitecombcrunch ()
714/* combine sites that have identical patterns (and nonzero weight) */
715{ /* sitecombcrunch */
716    int  i, sitei, j, sitej, k;
717    boolean  tied;
718
719    i = 0;
720    patsite[0] = patsite[1];
721    patweight[0] = 0;
722
723    for (j = 1; j <= sites; j++) {
724        tied = TRUE;
725        sitei = patsite[i];
726        sitej = patsite[j];
727
728        for (k = 1; tied && (k <= numsp); k++)
729            tied = (y[k][sitei] == y[k][sitej]);
730
731        if (tied) {
732            patweight[i] += weight[sitej];
733        }
734        else {
735            if (patweight[i] > 0) i++;
736            patweight[i] = weight[sitej];
737            patsite[i] = sitej;
738        }
739
740        pattern[sitej] = i;
741    }
742
743    endsite = i;
744    if (patweight[i] > 0) endsite++;
745} /* sitecombcrunch */
746
747
748void makeweights ()
749/* make up weights vector to avoid duplicate computations */
750{ /* makeweights */
751    int  i;
752
753    for (i = 1; i <= sites; i++)  patsite[i] = i;
754    sitesort();
755    sitecombcrunch();
756    if (endsite > maxpatterns) {
757        printf("ERROR:  Too many patterns in data\n");
758        printf("        Increase maxpatterns to at least %d\n", endsite);
759        anerror = TRUE;
760    }
761    else {
762        printf("Analyzing %d distinct data patterns (columns)\n\n", endsite);
763    }
764} /* makeweights */
765
766
767void makevalues (tr)
768     tree  *tr;
769     /* set up fractional likelihoods at tips */
770{ /* makevalues */
771    nodeptr  p;
772    int  i, j;
773
774    for (i = 1; i <= numsp; i++) {    /* Pack and move tip data */
775        for (j = 0; j < endsite; j++)
776            y[i-1][j] = y[i][patsite[j]];
777
778        p = tr->nodep[i];
779        p->tip = &(y[i-1][0]);
780    }
781
782} /* makevalues */
783
784
785void empiricalfreqs (tr)
786     tree  *tr;
787     /* Get empirical base frequencies from the data */
788{ /* empiricalfreqs */
789    double  sum, suma, sumc, sumg, sumt, wj, fa, fc, fg, ft;
790    int  i, j, k, code;
791    char *yptr;
792
793    freqa = 0.25;
794    freqc = 0.25;
795    freqg = 0.25;
796    freqt = 0.25;
797    for (k = 1; k <= 8; k++) {
798        suma = 0.0;
799        sumc = 0.0;
800        sumg = 0.0;
801        sumt = 0.0;
802        for (i = 1; i <= numsp; i++) {
803            yptr = tr->nodep[i]->tip;
804            for (j = 0; j < endsite; j++) {
805                code = *yptr++;
806                fa = freqa * ( code       & 1);
807                fc = freqc * ((code >> 1) & 1);
808                fg = freqg * ((code >> 2) & 1);
809                ft = freqt * ((code >> 3) & 1);
810                wj = patweight[j] / (fa + fc + fg + ft);
811                suma += wj * fa;
812                sumc += wj * fc;
813                sumg += wj * fg;
814                sumt += wj * ft;
815            }
816        }
817        sum = suma + sumc + sumg + sumt;
818        freqa = suma / sum;
819        freqc = sumc / sum;
820        freqg = sumg / sum;
821        freqt = sumt / sum;
822    }
823} /* empiricalfreqs */
824
825
826void getinput (tr)
827     tree  *tr;
828{ /* getinput */
829    getnums();                        if (anerror) return;
830    getoptions();                     if (anerror) return;
831    if (! freqsfrom) getbasefreqs();  if (anerror) return;
832    getyspace();                      if (anerror) return;
833    setuptree(tr, numsp);             if (anerror) return;
834    getdata (tr);                     if (anerror) return;
835    makeweights();                    if (anerror) return;
836    makevalues (tr);                  if (anerror) return;
837    if (freqsfrom) {
838        empiricalfreqs (tr);            if (anerror) return;
839        getbasefreqs();
840    }
841} /* getinput */
842
843
844xarray *setupxarray ()
845{ /* setupxarray */
846    xarray  *x;
847    xtype  *data;
848
849    x = (xarray *) malloc((unsigned) sizeof(xarray));
850    if (x) {
851        data = (xtype *) malloc((unsigned) (4 * endsite * sizeof(xtype)));
852        if (data) {
853            x->a = data;
854            x->c = data += endsite;
855            x->g = data += endsite;
856            x->t = data +  endsite;
857            x->prev = x->next = x;
858            x->owner = (node *) NULL;
859        }
860        else {
861            free ((char *) x);
862            return (xarray *) NULL;
863        }
864    }
865    return x;
866} /* setupxarray */
867
868
869void linkxarray (req, min, freexptr, usedxptr)
870     int  req, min;
871     xarray **freexptr, **usedxptr;
872     /*  Link a set of xarrays */
873{ /* linkxarray */
874    xarray  *first, *prev, *x;
875    int  i;
876
877    first = prev = (xarray *) NULL;
878    i = 0;
879
880    do {
881        x = setupxarray();
882        if (x) {
883            if (! first) first = x;
884            else {
885                prev->next = x;
886                x->prev = prev;
887            }
888            prev = x;
889            i++;
890        }
891        else {
892            printf("ERROR: Failure to get xarray memory.\n");
893            if (i < min) anerror = TRUE;
894        }
895    } while ((i < req) && x);
896
897    if (first) {
898        first->prev = prev;
899        prev->next = first;
900    }
901
902    *freexptr = first;
903    *usedxptr = (xarray *) NULL;
904} /* linkxarray */
905
906
907void setupnodex (tr)
908     tree  *tr;
909{ /* setupnodex */
910    nodeptr  p;
911    int  i;
912
913    for (i = tr->mxtips + 1; (i <= 2*(tr->mxtips) - 2) && ! anerror; i++) {
914        p = tr->nodep[i];
915        if ((anerror = !(p->x = setupxarray()))) break;
916    }
917} /* setupnodex */
918
919
920xarray *getxtip (p)
921     nodeptr  p;
922{ /* getxtip */
923    xarray  *new;
924    boolean  splice;
925
926    if (! p) return (xarray *) NULL;
927
928    splice = FALSE;
929
930    if (p->x) {
931        new = p->x;
932        if (new == new->prev) ;             /* linked to self; leave it */
933        else if (new == usedxtip) usedxtip = usedxtip->next; /* at head */
934        else if (new == usedxtip->prev) ;   /* already at tail */
935        else {                              /* move to tail of list */
936            new->prev->next = new->next;
937            new->next->prev = new->prev;
938            splice = TRUE;
939        }
940    }
941
942    else if (freextip) {
943        p->x = new = freextip;
944        new->owner = p;
945        if (new->prev != new) {            /* not only member of freelist */
946            new->prev->next = new->next;
947            new->next->prev = new->prev;
948            freextip = new->next;
949        }
950        else
951            freextip = (xarray *) NULL;
952
953        splice = TRUE;
954    }
955
956    else if (usedxtip) {
957        usedxtip->owner->x = (xarray *) NULL;
958        p->x = new = usedxtip;
959        new->owner = p;
960        usedxtip = usedxtip->next;
961    }
962
963    else {
964        printf ("ERROR: Unable to locate memory for a tip.\n");
965        anerror = TRUE;
966        exit(0);
967    }
968
969    if (splice) {
970        if (usedxtip) {                  /* list is not empty */
971            usedxtip->prev->next = new;
972            new->prev = usedxtip->prev;
973            usedxtip->prev = new;
974            new->next = usedxtip;
975        }
976        else
977            usedxtip = new->prev = new->next = new;
978    }
979
980    return  new;
981} /* getxtip */
982
983
984xarray *getxnode (p)
985     nodeptr  p;
986     /* Ensure that internal node p has memory */
987{ /* getxnode */
988    nodeptr  s;
989
990    if (! (p->x)) {  /*  Move likelihood array on this node to sector p */
991        if ((s = p->next)->x || (s = s->next)->x) {
992            p->x = s->x;
993            s->x = (xarray *) NULL;
994        }
995        else {
996            printf("ERROR in getxnode: Unable to locate memory at internal node.");
997            exit(0);
998        }
999    }
1000    return  p->x;
1001} /* getxnode */
1002
1003
1004void newview (p)                      /*  Update likelihoods at node */
1005     nodeptr  p;
1006{ /* newview */
1007    double   z1, lz1, xvlz1, z2, lz2, xvlz2,
1008        zz1, zv1, fx1r, fx1y, fx1n, suma1, sumg1, sumc1, sumt1,
1009        zz2, zv2, fx2r, fx2y, fx2n, ki, tempi, tempj;
1010    nodeptr  q, r;
1011    xtype   *x1a, *x1c, *x1g, *x1t, *x2a, *x2c, *x2g, *x2t,
1012        *x3a, *x3c, *x3g, *x3t;
1013    int  i;
1014
1015    if (p->tip) {             /*  Make sure that data are at tip */
1016        int code;
1017        char *yptr;
1018
1019        if (p->x) return;       /*  They are already there */
1020        (void) getxtip(p);      /*  They are not, so get memory */
1021        x3a = &(p->x->a[0]);    /*  Move tip data to xarray */
1022        x3c = &(p->x->c[0]);
1023        x3g = &(p->x->g[0]);
1024        x3t = &(p->x->t[0]);
1025        yptr = p->tip;
1026        for (i = 0; i < endsite; i++) {
1027            code = *yptr++;
1028            *x3a++ =  code       & 1;
1029            *x3c++ = (code >> 1) & 1;
1030            *x3g++ = (code >> 2) & 1;
1031            *x3t++ = (code >> 3) & 1;
1032        }
1033        return;
1034    }
1035
1036    /*  Internal node needs update */
1037
1038    q = p->next->back;
1039    r = p->next->next->back;
1040
1041    while ((! p->x) || (! q->x) || (! r->x)) {
1042        if (! q->x) newview(q);
1043        if (! r->x) newview(r);
1044        if (! p->x)  (void) getxnode(p);
1045    }
1046
1047    x1a = &(q->x->a[0]);
1048    x1c = &(q->x->c[0]);
1049    x1g = &(q->x->g[0]);
1050    x1t = &(q->x->t[0]);
1051    z1 = q->z;
1052    lz1 = (z1 > zmin) ? log(z1) : log(zmin);
1053    xvlz1 = xv * lz1;
1054
1055    x2a = &(r->x->a[0]);
1056    x2c = &(r->x->c[0]);
1057    x2g = &(r->x->g[0]);
1058    x2t = &(r->x->t[0]);
1059    z2 = r->z;
1060    lz2 = (z2 > zmin) ? log(z2) : log(zmin);
1061    xvlz2 = xv * lz2;
1062
1063    x3a = &(p->x->a[0]);
1064    x3c = &(p->x->c[0]);
1065    x3g = &(p->x->g[0]);
1066    x3t = &(p->x->t[0]);
1067
1068    { double  *rptr;
1069
1070        rptr = &(patrate[0]);
1071        for (i = 0; i < endsite; i++) {
1072            ki = *rptr++;
1073
1074            zz1 = exp(ki *   lz1);
1075            zv1 = exp(ki * xvlz1);
1076            fx1r = freqa * *x1a + freqg * *x1g;
1077            fx1y = freqc * *x1c + freqt * *x1t;
1078            fx1n = fx1r + fx1y;
1079            tempi = fx1r * invfreqr;
1080            tempj = zv1 * (tempi-fx1n) + fx1n;
1081            suma1 = zz1 * (*x1a++ - tempi) + tempj;
1082            sumg1 = zz1 * (*x1g++ - tempi) + tempj;
1083            tempi = fx1y * invfreqy;
1084            tempj = zv1 * (tempi-fx1n) + fx1n;
1085            sumc1 = zz1 * (*x1c++ - tempi) + tempj;
1086            sumt1 = zz1 * (*x1t++ - tempi) + tempj;
1087
1088            zz2 = exp(ki *   lz2);
1089            zv2 = exp(ki * xvlz2);
1090            fx2r = freqa * *x2a + freqg * *x2g;
1091            fx2y = freqc * *x2c + freqt * *x2t;
1092            fx2n = fx2r + fx2y;
1093            tempi = fx2r * invfreqr;
1094            tempj = zv2 * (tempi-fx2n) + fx2n;
1095            *x3a++ = suma1 * (zz2 * (*x2a++ - tempi) + tempj);
1096            *x3g++ = sumg1 * (zz2 * (*x2g++ - tempi) + tempj);
1097            tempi = fx2y * invfreqy;
1098            tempj = zv2 * (tempi-fx2n) + fx2n;
1099            *x3c++ = sumc1 * (zz2 * (*x2c++ - tempi) + tempj);
1100            *x3t++ = sumt1 * (zz2 * (*x2t++ - tempi) + tempj);
1101        }
1102    }
1103} /* newview */
1104
1105
1106void hookup (p, q, z)
1107     nodeptr  p, q;
1108     double   z;
1109{ /* hookup */
1110    p->back = q;
1111    q->back = p;
1112    p->z = q->z = z;
1113} /* hookup */
1114
1115
1116void initrav (p)
1117     nodeptr  p;
1118{ /* initrav */
1119    if (! p->tip) {
1120        initrav(p->next->back);
1121        initrav(p->next->next->back);
1122        newview(p);
1123    }
1124} /* initrav */
1125
1126
1127/*=======================================================================*/
1128/*                         Read a tree from a file                       */
1129/*=======================================================================*/
1130
1131int treeFinishCom ()
1132{ /* treeFinishCom */
1133    int      ch;
1134    boolean  inquote;
1135
1136    inquote = FALSE;
1137    while ((ch = getc(INFILE)) != EOF && (inquote || ch != ']')) {
1138        if (ch == '[' && ! inquote) {             /* comment; find its end */
1139            if ((ch = treeFinishCom()) == EOF)  break;
1140        }
1141        else if (ch == '\'') inquote = ! inquote;  /* start or end of quote */
1142    }
1143
1144    return  ch;
1145} /* treeFinishCom */
1146
1147
1148int treeGetCh ()
1149/* get next nonblank, noncomment character */
1150{ /* treeGetCh */
1151    int  ch;
1152
1153    while ((ch = getc(INFILE)) != EOF) {
1154        if (white(ch)) ;
1155        else if (ch == '[') {                   /* comment; find its end */
1156            if ((ch = treeFinishCom()) == EOF)  break;
1157        }
1158        else  break;
1159    }
1160
1161    return  ch;
1162} /* treeGetCh */
1163
1164
1165void treeFlushLabel()
1166{
1167    int     ch;
1168    boolean done;
1169
1170    if ((ch = treeGetCh()) == EOF)  return;
1171    done = (ch == ':' || ch == ',' || ch == ')'  || ch == '[' || ch == ';');
1172
1173    if (!done) {
1174        boolean quoted = (ch == '\'');
1175        if (quoted) ch = getc(INFILE);
1176
1177        while (! done) {
1178            if (quoted) {
1179                if ((ch = findch('\'')) == EOF)  return;      /* find close quote */
1180                ch = getc(INFILE);                            /* check next char */
1181                if (ch != '\'') done = TRUE;                  /* not doubled quote */
1182            }
1183            else if (ch == ':' || ch == ',' || ch == ')'  || ch == '['
1184                     || ch == ';' || ch == '\n' || ch == EOF) {
1185                done = TRUE;
1186            }
1187            if (! done)  done = ((ch = getc(INFILE)) == EOF);
1188        }
1189    }
1190
1191    if (ch != EOF)  (void) ungetc(ch, INFILE);
1192} /* treeFlushLabel */
1193
1194
1195int  findTipName (tr, ch)
1196     tree  *tr;
1197     int    ch;
1198{ /* findTipName */
1199    nodeptr  q;
1200    char  *nameptr, str[nmlngth+1];
1201    int  i, n;
1202    boolean  found, quoted, done;
1203
1204    if ((quoted = (ch == '\'')))  ch = getc(INFILE);
1205    done = FALSE;
1206    i = 0;
1207
1208    do {
1209        if (quoted) {
1210            if (ch == '\'') {
1211                ch = getc(INFILE);
1212                if (ch != '\'') done = TRUE;
1213            }
1214            else if (ch == EOF)
1215                done = TRUE;
1216            else if (ch == '\n' || ch == '\t')
1217                ch = ' ';
1218        }
1219        else if (ch == ':' || ch == ','  || ch == ')'  || ch == '['
1220                 || ch == '\n' || ch == EOF)
1221            done = TRUE;
1222        else if (ch == '_' || ch == '\t')
1223            ch = ' ';
1224
1225        if (! done) {
1226            if (i < nmlngth)  str[i++] = ch;
1227            ch = getc(INFILE);
1228        }
1229    } while (! done);
1230
1231    if (ch == EOF) {
1232        printf("ERROR: End-of-file in tree species name\n");
1233        return  0;
1234    }
1235
1236    (void) ungetc(ch, INFILE);
1237    while (i < nmlngth)  str[i++] = ' ';     /*  Pad name */
1238
1239    n = 1;
1240    do {
1241        q = tr->nodep[n];
1242        if (! (q->back)) {          /*  Only consider unused tips */
1243            i = 0;
1244            nameptr = q->name;
1245            do {found = str[i] == *nameptr++;}  while (found && (++i < nmlngth));
1246        }
1247        else
1248            found = FALSE;
1249    } while ((! found) && (++n <= tr->mxtips));
1250
1251    if (! found) {
1252        i = nmlngth;
1253        do {str[i] = '\0';} while (i-- && (str[i] <= ' '));
1254        printf("ERROR: Cannot find data for tree species: %s\n", str);
1255    }
1256
1257    return  (found ? n : 0);
1258} /* findTipName */
1259
1260
1261double processLength ()
1262{ /* processLength */
1263    double  branch;
1264    int     ch;
1265    char    string[41];
1266
1267    ch = treeGetCh();                            /*  Skip comments */
1268    if (ch != EOF)  (void) ungetc(ch, INFILE);
1269
1270    if (fscanf(INFILE, "%lf", &branch) != 1) {
1271        printf("ERROR: Problem reading branch length in processLength:\n");
1272        if (fscanf(INFILE, "%40s", string) == 1)  printf("%s\n", string);
1273        anerror = TRUE;
1274        branch = 0.0;
1275    }
1276
1277    return  branch;
1278} /* processLength */
1279
1280
1281void  treeFlushLen ()
1282{ /* treeFlushLen */
1283    int  ch;
1284
1285    if ((ch = treeGetCh()) == ':')
1286        (void) processLength();
1287    else if (ch != EOF)
1288        (void) ungetc(ch, INFILE);
1289
1290} /* treeFlushLen */
1291
1292
1293void  treeNeedCh (c1, where)
1294     int    c1;
1295     char  *where;
1296{ /* treeNeedCh */
1297    int c2, i;
1298
1299    if ((c2 = treeGetCh()) == c1)  return;
1300    /* assert(c2 == c1); */ /* stop here with debugger  */
1301
1302    printf("ERROR: Missing '%c' %s tree; ", c1, where);
1303    if (c2 == EOF)
1304        printf("End-of-File");
1305    else {
1306        putchar('\'');
1307        for (i = 24; i-- && (c2 != EOF); c2 = getc(INFILE))  putchar(c2);
1308        putchar('\'');
1309    }
1310    printf(" found instead\n");
1311    anerror = TRUE;
1312} /* treeNeedCh */
1313
1314
1315void  addElementLen (tr, p)
1316     tree    *tr;
1317     nodeptr  p;
1318{ /* addElementLen */
1319    double   z, branch;
1320    nodeptr  q;
1321    int      n, ch;
1322
1323    if ((ch = treeGetCh()) == '(') {     /*  A new internal node */
1324        n = (tr->nextnode)++;
1325        if (n > 2*(tr->mxtips) - 2) {
1326            if (tr->rooted || n > 2*(tr->mxtips) - 1) {
1327                printf("ERROR: Too many internal nodes.  Is tree rooted?\n");
1328                printf("       Deepest splitting should be a trifurcation.\n");
1329                anerror = TRUE;
1330                return;
1331            }
1332            else {
1333                tr->rooted = TRUE;
1334            }
1335        }
1336        q = tr->nodep[n];
1337        assert(q);
1338        addElementLen(tr, q->next);        if (anerror)  return;
1339        treeNeedCh(',', "in");             if (anerror)  return;
1340        assert(q->next);
1341        addElementLen(tr, q->next->next);  if (anerror)  return;
1342        treeNeedCh(')', "in");             if (anerror)  return;
1343        treeFlushLabel();                  if (anerror)  return;
1344    }
1345
1346    else {                               /*  A new tip */
1347        n = findTipName(tr, ch);
1348        if (n <= 0) {anerror = TRUE; return; }
1349        q = tr->nodep[n];
1350        if (tr->start->number > n)  tr->start = q;
1351        (tr->ntips)++;
1352    }                                  /* End of tip processing */
1353
1354    treeNeedCh(':', "in");             if (anerror)  return;
1355    branch = processLength();          if (anerror)  return;
1356    z = exp(-branch / fracchange);
1357    if (z > zmax)  z = zmax;
1358    hookup(p, q, z);
1359} /* addElementLen */
1360
1361
1362void uprootTree (tr, p)
1363     tree   *tr;
1364     nodeptr p;
1365{ /* uprootTree */
1366    nodeptr  q, r, s;
1367    int  n;
1368
1369    if (p->tip || p->back) {
1370        printf("ERROR: Unable to uproot tree.\n");
1371        printf("       Inappropriate node marked for removal.\n");
1372        anerror = TRUE;
1373        return;
1374    }
1375
1376    n = --(tr->nextnode);               /* last internal node added */
1377    if (n != tr->mxtips + tr->ntips - 1) {
1378        printf("ERROR: Unable to uproot tree.  Inconsistent\n");
1379        printf("       number of tips and nodes for rooted tree.\n");
1380        anerror = TRUE;
1381        return;
1382    }
1383
1384    q = p->next->back;                  /* remove p from tree */
1385    r = p->next->next->back;
1386    hookup(q, r, q->z * r->z);
1387
1388    q = tr->nodep[n];
1389    r = q->next;
1390    s = q->next->next;
1391    if (tr->ntips > 2 && p != q && p != r && p != s) {
1392        hookup(p,             q->back, q->z);   /* move connections to p */
1393        hookup(p->next,       r->back, r->z);
1394        hookup(p->next->next, s->back, s->z);
1395    }
1396
1397    q->back = r->back = s->back = (nodeptr) NULL;
1398    tr->rooted = FALSE;
1399} /* uprootTree */
1400
1401
1402void treeReadLen (tr)
1403     tree  *tr;
1404{ /* treeReadLen */
1405    nodeptr  p;
1406    int  i, ch;
1407
1408    for (i = 1; i <= tr->mxtips; i++) tr->nodep[i]->back = (node *) NULL;
1409    tr->start       = tr->nodep[tr->mxtips];
1410    tr->ntips       = 0;
1411    tr->nextnode    = tr->mxtips + 1;
1412    tr->opt_level   = 0;
1413    tr->smoothed    = FALSE;
1414    tr->rooted      = FALSE;
1415
1416    p = tr->nodep[(tr->nextnode)++];
1417    treeNeedCh('(', "at start of");                  if (anerror)  return;
1418    addElementLen(tr, p);                            if (anerror)  return;
1419    treeNeedCh(',', "in");                           if (anerror)  return;
1420    addElementLen(tr, p->next);                      if (anerror)  return;
1421    if (! tr->rooted) {
1422        if ((ch = treeGetCh()) == ',') {        /*  An unrooted format */
1423            addElementLen(tr, p->next->next);            if (anerror)  return;
1424        }
1425        else {                                  /*  A rooted format */
1426            p->next->next->back = (nodeptr) NULL;
1427            tr->rooted = TRUE;
1428            if (ch != EOF)  (void) ungetc(ch, INFILE);
1429        }
1430    }
1431    treeNeedCh(')', "in");
1432    if (anerror)  {
1433        printf("(This error also happens if the last species in the tree is unmarked)\n");
1434        return;
1435    }
1436
1437
1438    treeFlushLabel();                                if (anerror)  return;
1439    treeFlushLen();                                  if (anerror)  return;
1440    treeNeedCh(';', "at end of");                    if (anerror)  return;
1441
1442    if (tr->rooted)  uprootTree(tr, p->next->next);  if (anerror)  return;
1443    tr->start = p->next->next->back;  /* This is start used by treeString */
1444
1445    initrav(tr->start);
1446    initrav(tr->start->back);
1447} /* treeReadLen */
1448
1449/*=======================================================================*/
1450/*                           End of Tree Reading                         */
1451/*=======================================================================*/
1452
1453
1454double evaluate (tr, p)
1455     tree    *tr;
1456     nodeptr  p;
1457{ /* evaluate */
1458    double   sum, z, lz, xvlz,
1459        ki, zz, zv, fx1a, fx1c, fx1g, fx1t, fx1r, fx1y, fx2r, fx2y,
1460        suma, sumb, sumc, term;
1461    double  *log_f, *rptr;
1462    xtype   *x1a, *x1c, *x1g, *x1t, *x2a, *x2c, *x2g, *x2t;
1463    nodeptr  q;
1464    int  i, *wptr;
1465
1466    q = p->back;
1467    while ((! p->x) || (! q->x)) {
1468        if (! (p->x)) newview(p);
1469        if (! (q->x)) newview(q);
1470    }
1471
1472    x1a = &(p->x->a[0]);
1473    x1c = &(p->x->c[0]);
1474    x1g = &(p->x->g[0]);
1475    x1t = &(p->x->t[0]);
1476
1477    x2a = &(q->x->a[0]);
1478    x2c = &(q->x->c[0]);
1479    x2g = &(q->x->g[0]);
1480    x2t = &(q->x->t[0]);
1481
1482    z = p->z;
1483    if (z < zmin) z = zmin;
1484    lz = log(z);
1485    xvlz = xv * lz;
1486
1487    wptr = &(patweight[0]);
1488    rptr = &(patrate[0]);
1489    log_f = tr->log_f;
1490    sum = 0.0;
1491
1492    for (i = 0; i < endsite; i++) {
1493        fx1a  = freqa * *x1a++;
1494        fx1g  = freqg * *x1g++;
1495        fx1c  = freqc * *x1c++;
1496        fx1t  = freqt * *x1t++;
1497        suma  = fx1a * *x2a + fx1c * *x2c + fx1g * *x2g + fx1t * *x2t;
1498        fx2r  = freqa * *x2a++ + freqg * *x2g++;
1499        fx2y  = freqc * *x2c++ + freqt * *x2t++;
1500        fx1r  = fx1a + fx1g;
1501        fx1y  = fx1c + fx1t;
1502        sumc  = (fx1r + fx1y) * (fx2r + fx2y);
1503        sumb  = fx1r * fx2r * invfreqr + fx1y * fx2y * invfreqy;
1504        suma -= sumb;
1505        sumb -= sumc;
1506
1507        ki = *rptr++;
1508        zz = exp(ki *   lz);
1509        zv = exp(ki * xvlz);
1510
1511        term = log(zz * suma + zv * sumb + sumc);
1512        sum += *wptr++ * term;
1513        *log_f++ = term;
1514    }
1515
1516    tr->likelihood = sum;
1517    return  sum;
1518} /* evaluate */
1519
1520
1521void dli_dki (p)  /*  d(Li)/d(ki) */
1522     nodeptr  p;
1523{ /* dli_dki */
1524    double   z, lz, xvlz;
1525    double   ki, fx1a, fx1c, fx1g, fx1t, fx1r, fx1y, fx2r, fx2y,
1526        suma, sumb, sumc;
1527    double  *rptr;
1528    xtype   *x1a, *x1c, *x1g, *x1t, *x2a, *x2c, *x2g, *x2t;
1529    nodeptr  q;
1530    int  i, *wptr;
1531
1532
1533    q = p->back;
1534    while ((! p->x) || (! q->x)) {
1535        if (! p->x) newview(p);
1536        if (! q->x) newview(q);
1537    }
1538
1539    x1a = &(p->x->a[0]);
1540    x1c = &(p->x->c[0]);
1541    x1g = &(p->x->g[0]);
1542    x1t = &(p->x->t[0]);
1543
1544    x2a = &(q->x->a[0]);
1545    x2c = &(q->x->c[0]);
1546    x2g = &(q->x->g[0]);
1547    x2t = &(q->x->t[0]);
1548
1549    z = p->z;
1550    if (z < zmin) z = zmin;
1551    lz = log(z);
1552    xvlz = xv * lz;
1553
1554    rptr = &(patrate[0]);
1555    wptr = &(patweight[0]);
1556
1557    for (i = 0; i < endsite; i++) {
1558        fx1a  = freqa * *x1a++;
1559        fx1g  = freqg * *x1g++;
1560        fx1c  = freqc * *x1c++;
1561        fx1t  = freqt * *x1t++;
1562        suma  = fx1a * *x2a + fx1c * *x2c + fx1g * *x2g + fx1t * *x2t;
1563        fx2r  = freqa * *x2a++ + freqg * *x2g++;
1564        fx2y  = freqc * *x2c++ + freqt * *x2t++;
1565        fx1r  = fx1a + fx1g;
1566        fx1y  = fx1c + fx1t;
1567        sumc  = (fx1r + fx1y) * (fx2r + fx2y);
1568        sumb  = fx1r * fx2r * invfreqr + fx1y * fx2y * invfreqy;
1569        suma -= sumb;
1570        sumb -= sumc;
1571
1572        ki    = *rptr++;
1573        suma *= exp(ki *   lz);
1574        sumb *= exp(ki * xvlz);
1575        dLidki[i] += *wptr++ * lz * (suma + sumb*xv);
1576    }
1577} /* dli_dki */
1578
1579
1580void spanSubtree (p)
1581     nodeptr  p;
1582{ /* spanSubtree */
1583    dli_dki (p);
1584
1585    if (! p->tip) {
1586        spanSubtree(p->next->back);
1587        spanSubtree(p->next->next->back);
1588    }
1589} /* spanSubtree */
1590
1591
1592void findSiteRates (tr, ki_min, ki_max, d_ki, max_error)
1593     tree    *tr;
1594     double   ki_min, ki_max, d_ki, max_error;
1595{ /* findSiteRates */
1596    double  inv_d_ki, ki;
1597    int     i;
1598
1599    if (ki_min <= 0.0 || ki_max <= ki_min) {
1600        printf("ERROR: Bad rate value limits to findSiteRates\n");
1601        anerror = TRUE;
1602        return;
1603    }
1604    else if (d_ki <= 1.0) {
1605        printf("ERROR: Bad rate step to findSiteRates\n");
1606        anerror = TRUE;
1607        return;
1608    }
1609
1610    for (i = 0; i < endsite; i++) {
1611        bestki[i] = 1.0;                       /*  dummy initial rates */
1612        bestLi[i] = unlikely;
1613    }
1614
1615    for (ki = ki_min; ki <= ki_max; ki *= d_ki) {
1616        for (i = 0; i < endsite; i++)  patrate[i] = ki;
1617        initrav(tr->start);
1618        initrav(tr->start->back);
1619        (void) evaluate(tr, tr->start->back);
1620        for (i = 0; i < endsite; i++) {
1621            if (tr->log_f[i] > bestLi[i]) {
1622                bestki[i] = ki;
1623                bestLi[i] = tr->log_f[i];
1624            }
1625        }
1626    }
1627
1628    for (i = 0; i < endsite; i++)  patrate[i] = bestki[i];
1629    initrav(tr->start);
1630    initrav(tr->start->back);
1631
1632    while (d_ki > 1.0 + max_error) {
1633        for (i = 0; i < endsite; i++)  dLidki[i] = 0.0;
1634        spanSubtree(tr->start->back);
1635        if (! tr->start->tip) {
1636            spanSubtree(tr->start->next->back);
1637            spanSubtree(tr->start->next->next->back);
1638        }
1639
1640        d_ki = sqrt(d_ki);
1641        inv_d_ki = 1.0/d_ki;
1642        for (i = 0; i < endsite; i++) {
1643            if (dLidki[i] > 0.0) {
1644                patrate[i] *= d_ki;
1645                if (patrate[i] > ki_max) patrate[i] = ki_max;
1646            }
1647            else {
1648                patrate[i] *= inv_d_ki;
1649                if (patrate[i] < ki_min) patrate[i] = ki_min;
1650            }
1651        }
1652
1653        initrav(tr->start);
1654        initrav(tr->start->back);
1655    }
1656} /* findSiteRates */
1657
1658
1659double  subtreeLength (p)
1660     nodeptr  p;
1661{ /* subtreeLength */
1662    double sum;
1663
1664    sum = -fracchange * log(p->z);
1665    if (! p->tip) {
1666        sum += subtreeLength(p->next->back);
1667        sum += subtreeLength(p->next->next->back);
1668    }
1669
1670    return sum;
1671} /* subtreeLength */
1672
1673
1674double  treeLength(tr)
1675     tree  *tr;
1676{ /* treeLength */
1677    double sum;
1678
1679    sum = subtreeLength(tr->start->back);
1680    if (! tr->start->tip) {
1681        sum += subtreeLength(tr->start->next->back);
1682        sum += subtreeLength(tr->start->next->next->back);
1683    }
1684
1685    return sum;
1686} /* treeLength */
1687
1688
1689void categorize (Sites, Categs, Weight, Pattern, Patrate,
1690                 categrate, sitecateg)
1691     int Sites;
1692     int Categs;
1693     int Weight[];                                        /* one based */
1694int      Pattern[];                                       /* one based */
1695double   Patrate[];                                       /* zero based */
1696double   categrate[];                                     /* zero based */
1697int      sitecateg[];                                     /* one based */
1698{ /* categorize */
1699    double  ki, min_1, min_2, max_1, max_2, a, b;
1700    int  i, k;
1701
1702    min_1 = 1.0E37;
1703    min_2 = 1.0E37;
1704    max_1 = 0.0;
1705    max_2 = 0.0;
1706    for (i = 1; i <= Sites; i++) {
1707        if (Weight[i] > 0) {
1708            ki = Patrate[Pattern[i]];
1709            if (ki < min_2) {
1710                if (ki < min_1) {
1711                    if (ki < 0.995 * min_1)  min_2 = min_1;
1712                    min_1 = ki;
1713                }
1714                else if (ki > 1.005 * min_1) {
1715                    min_2 = ki;
1716                }
1717            }
1718            else if (ki > max_2) {
1719                if (ki > max_1) {
1720                    if (ki > 1.005 * max_1)  max_2 = max_1;
1721                    max_1 = ki;
1722                }
1723                else if (ki < 0.995 * max_1) {
1724                    max_2 = ki;
1725                }
1726            }
1727        }
1728    }
1729
1730    a = (Categs - 3.0)/log(max_2/min_2);
1731    b = - a * log(min_2) + 2.0;
1732
1733    categrate[0] = min_1;
1734    for (k = 1; k <= Categs-2; k++)  categrate[k] = min_2 * exp((k-1)/a);
1735    categrate[Categs-1] = max_1;
1736
1737    for (i = 1; i <= Sites; i++) {
1738        if (Weight[i] > 0) {
1739            ki = Patrate[Pattern[i]];
1740            if      (ki < 0.99 * min_2) sitecateg[i] = 1;
1741            else if (ki > 1.00 * max_2) sitecateg[i] = Categs;
1742            else sitecateg[i] = nint(a * log(Patrate[Pattern[i]]) + b);
1743        }
1744        else
1745            sitecateg[i] = Categs;
1746    }
1747} /* categorize */
1748
1749
1750char *arb_filter;
1751char *alignment_name;
1752
1753GBDATA *gb_main;
1754/** Get the calling filter, needed to expand weights afterwards */
1755void getArbFilter(){
1756    GB_begin_transaction(gb_main);
1757    arb_filter = GBT_read_string(gb_main,AWAR_GDE_EXPORT_FILTER);
1758    alignment_name = GBT_get_default_alignment(gb_main);
1759    GB_commit_transaction(gb_main);
1760}
1761
1762void writeToArb(){
1763    char    sai_name[1024];
1764    char    type_info[1024];
1765    char    category_string[1024];
1766    char   *p;
1767    long    ali_len;
1768    int     i,k;
1769    int     ali_pos;
1770    float  *rates;              /* rates to export */
1771    char   *cats;               /* categories */
1772    GBDATA *gb_sai;
1773    GBDATA *gb_data;
1774
1775
1776    GB_begin_transaction(gb_main);
1777
1778    /* First create new SAI */
1779    sprintf(sai_name, "POS_VAR_BY_ML_%i",getpid());
1780    printf("Writing '%s'\n", sai_name);
1781    ali_len = GBT_get_alignment_len(gb_main,alignment_name);
1782    cats    = (char *)GB_calloc(ali_len,sizeof(char));
1783    rates   = (float *)GB_calloc(ali_len,sizeof(float));
1784
1785    /* fill in rates and categories */
1786    {
1787        double  categrate[maxcategories]; /* rate of a given category */
1788        int     sitecateg[maxsites+1];    /* category of a given site */
1789
1790        categorize(sites, categs, weight, pattern, patrate,categrate, sitecateg);
1791
1792        i = 1;                  /* thanks to pascal */
1793        for ( ali_pos = 0; ali_pos < ali_len; ali_pos++){
1794            if (arb_filter[ali_pos] == '0'){
1795                cats[ali_pos] = '.';
1796                rates[ali_pos] = KI_MAX;
1797                continue; /* filter says not written */
1798            }
1799            if (weight[i] > 0) {
1800                rates[ali_pos] = patrate[pattern[i]];
1801                cats[ali_pos] = itobase36(categs - sitecateg[i]);
1802            }else{
1803                rates[i] = KI_MAX;      /* infinite rate */
1804                cats[ali_pos] = '.';
1805            }
1806            i++;
1807        }
1808
1809        /* write categories */
1810        p    = category_string;
1811        p[0] = 0; /* if no categs */
1812        for (k = 1; k <= categs; k ++) {
1813            sprintf(p," %G", categrate[categs-k]);
1814            p += strlen(p);
1815        }
1816    }
1817
1818
1819    sprintf(type_info, "PVML: Positional Variability by ML (Olsen)");
1820    gb_sai = GBT_find_or_create_SAI(gb_main,sai_name);
1821    if (!gb_sai){
1822        GB_print_error();
1823    }else{
1824        gb_data = GBT_add_data(gb_sai, alignment_name, "rates", GB_FLOATS);
1825        GB_write_floats(gb_data, rates,ali_len);
1826        gb_data = GBT_add_data(gb_sai, alignment_name, "data", GB_STRING);
1827        GB_write_string(gb_data, cats);
1828        gb_data = GBT_add_data(gb_sai, alignment_name, "_CATEGORIES", GB_STRING);
1829        GB_write_string(gb_data,category_string);
1830        gb_data = GBT_add_data(gb_sai, alignment_name, "_TYPE", GB_STRING);
1831        GB_write_string(gb_data, type_info);
1832    }
1833
1834    GB_commit_transaction(gb_main);
1835}
1836
1837void openArb(){
1838    gb_main = GB_open(":","rw");
1839    if (!gb_main){
1840        GB_warning("Cannot find ARB server");
1841        exit(-1);;
1842    }
1843}
1844
1845void closeArb(){
1846    GB_close(gb_main);
1847}
1848
1849void wrfile (outfile, Sites, Categs, Weight, categrate, sitecateg)
1850     FILE   *outfile;
1851     int     Sites;
1852     int     Categs;
1853     int     Weight[];      /* one based */
1854double  categrate[];   /* zero based */
1855int     sitecateg[];   /* one based */
1856{ /* wrfile */
1857
1858
1859    int  i, k, l;
1860
1861
1862    for (k = 1; k <= Sites; k += 60) {
1863        l = k + 59;
1864        if (l > Sites)  l = Sites;
1865        fprintf(outfile, "%s  ", k == 1 ? "Weights   " : "          ");
1866
1867        for (i = k; i <= l; i++) {
1868            putc(itobase36(Weight[i]), outfile);
1869            if (((i % 10) == 0) && ((i % 60) != 0)) putc(' ', outfile);
1870        }
1871
1872        putc('\n', outfile);
1873    }
1874    for (k = 1; k <= Categs; k += 7) {
1875        l = k + 6;
1876        if (l > Categs)  l = Categs;
1877        if (k == 1)  fprintf(outfile, "C %2d", Categs);
1878        else         fprintf(outfile, "    ");
1879
1880        for (i = k-1; i < l; i++)  fprintf(outfile, " %9.5f", categrate[i]);
1881
1882        putc('\n', outfile);
1883    }
1884
1885    for (k = 1; k <= Sites; k += 60) {
1886        l = k + 59;
1887        if (l > Sites)  l = Sites;
1888        fprintf(outfile, "%s  ", k == 1 ? "Categories" : "          ");
1889
1890        for (i = k; i <= l; i++) {
1891            putc(itobase36(sitecateg[i]), outfile);
1892            if (((i % 10) == 0) && ((i % 60) != 0)) putc(' ', outfile);
1893        }
1894
1895        putc('\n', outfile);
1896    }
1897
1898} /* wrfile */
1899
1900
1901void summarize (treenum)
1902     int    treenum;
1903{ /* summarize */
1904    int  i;
1905
1906    printf("  Site      Rate\n");
1907    printf("  ----   ---------\n");
1908
1909    for (i = 1; i <= sites; i++) {
1910        if (weight[i] > 0) {
1911            printf("%6d %11.4f\n", i, patrate[pattern[i]]);
1912        }
1913        else
1914            printf("%6d   Undefined\n", i);
1915    }
1916
1917    putchar('\n');
1918
1919    if (categs > 1) {
1920        double  categrate[maxcategories]; /* rate of a given category */
1921        int     sitecateg[maxsites+1];    /* category of a given site */
1922
1923        categorize(sites, categs, weight, pattern, patrate,categrate, sitecateg);
1924
1925        wrfile(stdout, sites, categs, weight, categrate, sitecateg);
1926        putchar('\n');
1927
1928        if (writefile) {
1929            char  filename[512];
1930            FILE *outfile;
1931
1932            if (treenum <= 0)
1933                (void) sprintf(filename, "%s.%d", "weight_rate", getpid());
1934            else
1935                (void) sprintf(filename, "%s_%2d.%d", "weight_rate",
1936                               treenum, getpid());
1937
1938            outfile = fopen(filename, "w");
1939            if (outfile) {
1940                wrfile(outfile, sites, categs, weight, categrate, sitecateg);
1941                (void) fclose(outfile);
1942                printf("Weights and categories also writen to %s\n\n", filename);
1943            }
1944        }
1945    }
1946} /* summarize */
1947
1948
1949void makeUserRates (tr)
1950     tree   *tr;
1951{ /* makeUserRates */
1952    double  tree_length;
1953    int     numtrees, which, i;
1954
1955    if (fscanf(INFILE, "%d", &numtrees) != 1 || findch('\n') == EOF) {
1956        printf("ERROR: Problem reading number of user trees\n");
1957        anerror = TRUE;
1958        return;
1959    }
1960
1961    printf("User-defined %s:\n\n", (numtrees == 1) ? "tree" : "trees");
1962
1963    for (which = 1; which <= numtrees; which++) {
1964        for (i = 0; i < endsite; i++)  patrate[i] = 1.0;
1965        treeReadLen(tr);                                 if (anerror) break;
1966        tree_length = treeLength(tr);
1967
1968        printf("%d taxon user-supplied tree read\n\n", tr->ntips);
1969        printf("Total length of tree branches = %8.6f\n\n", tree_length);
1970
1971        findSiteRates(tr, 0.5/tree_length, KI_MAX, RATE_STEP, MAX_ERROR);
1972        if (anerror) break;
1973
1974        summarize(numtrees == 1 ? 0 : which);            if (anerror) break;
1975    }
1976
1977} /* makeUserRates */
1978
1979
1980int main ()
1981{ /* Maximum Likelihood Site Rate */
1982    tree   curtree, *tr;
1983
1984#   if DebugData
1985    debug = fopen("ml_rates.debug","w");
1986#   endif
1987
1988#   if ! UseStdin
1989    fp = fopen("infile", "rt");
1990#   endif
1991
1992#if defined(DEBUG)
1993    {
1994        int i;
1995        for (i = 0; i<(2*maxsp-1); ++i) {
1996            curtree.nodep[i] = 0;
1997        }
1998    }
1999#endif /* DEBUG */
2000
2001    openArb();
2002    getArbFilter();
2003
2004    tr = & curtree;
2005    getinput(tr);                               if (anerror)  return 1;
2006    linkxarray(3, 3, & freextip, & usedxtip);   if (anerror)  return 1;
2007    setupnodex(tr);                             if (anerror)  return 1;
2008    makeUserRates(tr);                          if (anerror)  return 1;
2009
2010    writeToArb();
2011    closeArb();
2012
2013    freeTree(tr);
2014
2015#   if ! UseStdin
2016    (void) fclose(fp);
2017#   endif
2018
2019#   if DebugData
2020    fclose(debug);
2021#   endif
2022
2023    return 0;
2024} /* Maximum Likelihood Site Rate */
Note: See TracBrowser for help on using the repository browser.