source: tags/svn.1.5.4/TOOLS/arb_dnarates.cxx

Last change on this file was 8309, checked in by westram, 14 years ago
  • moved much code into static scope

(partly reverted by [8310])

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