source: branches/profile/TOOLS/arb_dnarates.cxx

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