25 #define programName "DNAml_rates"
26 #define programVersion "1.0.0"
27 #define programDate "April 11, 1992"
36 #define assert(bed) arb_assert(bed)
89 printf(
"\n%s, version %s, %s\n\n",
93 printf(
"Portions based on Joseph Felsenstein's Nucleic acid sequence\n");
94 printf(
"Maximum Likelihood method, version 3.3\n\n");
96 if (fscanf(INFILE,
"%d %d", &
numsp, &
sites) != 2) {
97 printf(
"ERROR: Problem reading number of species and sites\n");
101 printf(
"%d Species, %d Sites\n\n",
numsp,
sites);
104 printf(
"ERROR: Too many species; adjust program constants\n");
107 else if (
numsp < 4) {
108 printf(
"ERROR: Too few species (need at least 4)\n");
113 printf(
"ERROR: Too many sites; adjust program constants\n");
116 else if (
sites < 1) {
117 printf(
"ERROR: Too few sites\n");
122 inline int min(
int i,
int j) {
return i<j ? i : j; }
123 inline bool digit(
int ch) {
return ch >=
'0' && ch <=
'9'; }
124 inline bool white(
int ch) {
return ch ==
' ' || ch ==
'\n' || ch ==
'\t'; }
129 if ((ch >=
'a' && ch <=
'i') ||
130 (ch >=
'j' && ch <=
'r') ||
131 (ch >=
's' && ch <=
'z'))
133 *chptr = ch +
'A' -
'a';
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;
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';
162 while ((ch = getc(INFILE)) != EOF && ch != c) ;
170 for (
int i = 2; i <=
nmlngth; i++) (
void) getc(INFILE);
175 int ch = getc(INFILE);
180 else if (!
white(ch)) {
181 printf(
"ERROR: Bad weight character: '%c'", ch);
182 printf(
" Weights must be a digit or a letter.\n");
188 if (
findch(
'\n', INFILE) == EOF) {
189 printf(
"ERROR: Missing newline at end of weight data\n");
207 while ((ch = getc(INFILE)) !=
'\n' && ch != EOF) {
211 case 'C':
categs = -1; extranum++;
break;
215 case 'M':
mininfo = 0; extranum++;
break;
216 case 'T':
ttratio = -1.0; extranum++;
break;
223 printf(
"ERROR: Bad option character: '%c'\n", ch);
230 printf(
"ERROR: End-of-file in options list\n");
237 while (extranum-- && !
anerror) {
244 printf(
"ERROR: Unexpected Categories data\n");
247 else if (fscanf(INFILE,
"%d", &
categs) != 1 ||
findch(
'\n', INFILE)==EOF) {
248 printf(
"ERROR: Problem reading number of rate categories\n");
252 printf(
"ERROR: Bad number of rate categories: %d\n",
categs);
259 printf(
"ERROR: Unexpected Min informative residues data\n");
262 else if (fscanf(INFILE,
"%d", &
mininfo)!=1 ||
findch(
'\n', INFILE)==EOF) {
263 printf(
"ERROR: Problem reading min informative residues\n");
266 else if (mininfo < 2 || mininfo >
numsp) {
267 printf(
"ERROR: Bad number for informative residues: %d\n",
275 printf(
"ERROR: Unexpected Transition/transversion data\n");
278 else if (fscanf(INFILE,
"%lf", &
ttratio)!=1 ||
findch(
'\n', INFILE)==EOF) {
279 printf(
"ERROR: Problem reading transition/transversion data\n");
286 printf(
"ERROR: Unexpected Weights data\n");
295 printf(
"ERROR: Auxiliary data line starts with '%c'\n", ch);
304 printf(
"ERROR: Category data missing from input\n");
310 printf(
"ERROR: Minimum informative residues missing from input\n");
315 printf(
"There must be at least %d informative residues per column\n\n",
mininfo);
319 printf(
"ERROR: Transition/transversion data missing from input\n");
329 printf(
"ERROR: Weight data invalid or missing from input\n");
339 printf(
"Base Frequencies:\n\n");
343 ||
findch(
'\n', INFILE) == EOF) {
344 printf(
"ERROR: Problem reading user base frequencies\n");
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);
361 printf(
"Transition/transversion ratio = %10.6f\n\n",
ttratio);
364 xi = suma/(suma+sumb);
368 printf(
"WARNING: This transition/transversion ratio is\n");
369 printf(
" impossible with these base frequencies!\n");
372 printf(
"Transition/transversion parameter reset\n\n");
374 printf(
"(Transition/transversion parameter = %10.6f)\n\n",
xi/
xv);
377 - freqg*freqg - freqt*freqt);
382 long size = 4 * (
sites/4 + 1);
383 char *y0 = (
char*)malloc((
unsigned) (
sizeof(
char) * size * (
numsp+1)));
386 printf(
"ERROR: Unable to obtain space for data array\n");
390 for (
int i = 0; i <=
numsp; i++) {
401 for (
int i = 1; i <= numSp; i++) {
403 if (!p) {
anerror =
true;
break; }
414 for (
int i = numSp+1; i <= nodes && !
anerror; i++) {
416 for (
int j = 1; j <= 3; j++) {
418 if (!p) {
anerror =
true;
break; }
441 if (
anerror) printf(
"ERROR: Unable to obtain sufficient memory");
462 for (
int i = leafs+1; i <= nodes; i++) {
481 for (
int i = 0; i <= 255; i++) meaning[i] = 0;
508 bool allread =
false;
509 bool firstpass =
true;
514 for (
int i = 1; i <=
numsp; i++) {
517 while (
white(ch = getc(INFILE))) {
518 if (ch ==
'\n') j = 1;
else j++;
522 printf(
"ERROR: Blank name for species %d; ", i);
523 printf(
"check number of species,\n");
524 printf(
" number of sites, and interleave option.\n");
530 for (
int k = 1; k < j; k++) *nameptr++ =
' ';
532 while (ch !=
'\n' && ch != EOF) {
533 if (ch ==
'_' ||
white(ch)) ch =
' ';
539 while (j++ <=
nmlngth) *nameptr++ =
' ';
543 printf(
"ERROR: End-of-file in name of species %d\n", i);
550 while ((j <
sites) && ((ch = getc(INFILE)) != EOF)
553 if (meaning[ch] || ch ==
'.') {
556 if (i != 1) ch =
y[1][j];
558 printf(
"ERROR: Dot (.) found at site %d of sequence 1\n", j);
567 printf(
"ERROR: Bad base (%c) at site %d of sequence %d\n",
575 printf(
"ERROR: End-of-file at site %d of sequence %d\n", j, i);
580 if (! firstpass && (j == basesread)) i--;
581 else if (i == 1) basesnew = j;
582 else if (j != basesnew) {
583 printf(
"ERROR: Sequences out of alignment\n");
584 printf(
" %d (instead of %d) residues read in sequence %d\n",
585 j - basesread, basesnew - basesread, i);
590 while (ch !=
'\n' && ch != EOF) ch = getc(INFILE);
593 basesread = basesnew;
594 allread = (basesread >=
sites);
603 printf(
"Name");
for (
int i=1; i<=j; i++) putchar(
' '); printf(
"Sequences\n");
604 printf(
"----");
for (
int i=1; i<=j; i++) putchar(
' '); printf(
"---------\n");
607 for (
int i = 1; i <=
sites; i += 60) {
613 for (j = 11; j <=
nmlngth+3; j++) putchar(
' ');
614 for (
int k = i; k <= l; k++) {
616 if (((k % 10) == 0) && ((k % 60) != 0)) putchar(
' ');
621 for (j = 1; j <=
numsp; j++) {
623 for (
int k = i; k <= l; k++) {
625 if ((j > 1) && (ch ==
y[1][k])) ch =
'.';
627 if (((k % 10) == 0) && ((k % 60) != 0)) putchar(
' ');
637 for (
int i = 1; i <=
sites; i++)
info[i] = 0;
639 for (
int j = 1; j <=
numsp; j++) {
640 for (
int i = 1; i <=
sites; i++) {
641 if ((
y[j][i] = meaning[(
int)
y[j][i]]) != 15)
info[i]++;
652 for (
int gap =
sites/2; gap > 0; gap /= 2) {
653 for (
int i = gap + 1; i <=
sites; i++) {
664 for (
int k = 1; tied && (k <=
numsp); k++) {
665 flip = (
y[k][jj] >
y[k][jg]);
666 tied = (
y[k][jj] ==
y[k][jg]);
674 while (flip && (j > 0));
686 for (
int j = 1; j <=
sites; j++) {
691 for (
int k = 1; tied && (k <=
numsp); k++) {
692 tied = (
y[k][sitei] ==
y[k][sitej]);
718 printf(
"ERROR: Too many patterns in data\n");
719 printf(
" Increase maxpatterns to at least %d\n",
endsite);
723 printf(
"Analyzing %d distinct data patterns (columns)\n\n",
endsite);
730 for (
int i = 1; i <=
numsp; i++) {
731 for (
int j = 0; j <
endsite; j++) {
735 p->
tip = &(
y[i-1][0]);
748 for (
int k = 1; k <= 8; k++) {
755 for (
int i = 1; i <=
numsp; i++) {
757 for (
int j = 0; j <
endsite; j++) {
760 double fa =
freqa * (code & 1);
761 double fc =
freqc * ((code >> 1) & 1);
762 double fg =
freqg * ((code >> 2) & 1);
763 double ft =
freqt * ((code >> 3) & 1);
765 double wj =
patweight[j] / (fa + fc + fg + ft);
774 double sum = suma + sumc + sumg + sumt;
835 if (! first) first = x;
844 printf(
"ERROR: Failure to get xarray memory.\n");
847 }
while ((i < req) && x);
864 if (!p->
x) {
anerror =
true;
break; }
875 if (new_xarray == new_xarray->
prev) ;
889 new_xarray->
owner = p;
890 if (new_xarray->
prev != new_xarray) {
906 new_xarray->
owner = p;
911 printf (
"ERROR: Unable to locate memory for a tip.\n");
936 if ((s = p->
next)->x || (s = s->
next)->x) {
941 printf(
"ERROR in getxnode: Unable to locate memory at internal node.");
962 for (
int i = 0; i <
endsite; i++) {
965 *x3c++ = (code >> 1) & 1;
966 *x3g++ = (code >> 2) & 1;
967 *x3t++ = (code >> 3) & 1;
977 while ((! p->
x) || (! q->
x) || (! r->
x)) {
989 double lz1 = (z1 >
zmin) ? log(z1) : log(
zmin);
990 double xvlz1 =
xv * lz1;
998 double lz2 = (z2 >
zmin) ? log(z2) : log(
zmin);
999 double xvlz2 =
xv * lz2;
1008 for (
int i = 0; i <
endsite; i++) {
1009 double ki = *rptr++;
1011 double zz1 = exp(ki * lz1);
1012 double zv1 = exp(ki * xvlz1);
1016 double fx1n = fx1r + fx1y;
1019 double tempj = zv1 * (tempi-fx1n) + fx1n;
1021 double suma1 = zz1 * (*x1a++ - tempi) + tempj;
1022 double sumg1 = zz1 * (*x1g++ - tempi) + tempj;
1025 tempj = zv1 * (tempi-fx1n) + fx1n;
1027 double sumc1 = zz1 * (*x1c++ - tempi) + tempj;
1028 double sumt1 = zz1 * (*x1t++ - tempi) + tempj;
1030 double zz2 = exp(ki * lz2);
1031 double zv2 = exp(ki * xvlz2);
1035 double fx2n = fx2r + fx2y;
1038 tempj = zv2 * (tempi-fx2n) + fx2n;
1040 *x3a++ = suma1 * (zz2 * (*x2a++ - tempi) + tempj);
1041 *x3g++ = sumg1 * (zz2 * (*x2g++ - tempi) + tempj);
1044 tempj = zv2 * (tempi-fx2n) + fx2n;
1046 *x3c++ = sumc1 * (zz2 * (*x2c++ - tempi) + tempj);
1047 *x3t++ = sumt1 * (zz2 * (*x2t++ - tempi) + tempj);
1076 while ((ch = getc(INFILE)) != EOF && (inquote || ch !=
']')) {
1077 if (ch ==
'[' && ! inquote) {
1080 else if (ch ==
'\'') inquote = !
inquote;
1091 while ((ch = getc(INFILE)) != EOF) {
1093 else if (ch ==
'[') {
1104 if (ch == EOF)
return;
1106 bool done = (ch ==
':' || ch ==
',' || ch ==
')' || ch ==
'[' || ch ==
';');
1108 bool quoted = (ch ==
'\'');
1109 if (quoted) ch = getc(INFILE);
1113 if ((ch =
findch(
'\'', INFILE)) == EOF)
return;
1115 if (ch !=
'\'') done =
true;
1117 else if (ch ==
':' || ch ==
',' || ch ==
')' || ch ==
'['
1118 || ch ==
';' || ch ==
'\n' || ch == EOF) {
1121 if (! done) done = ((ch = getc(INFILE)) == EOF);
1125 if (ch != EOF) (void) ungetc(ch, INFILE);
1130 bool quoted = (ch ==
'\'');
1131 if (quoted) ch = getc(INFILE);
1141 if (ch !=
'\'') done =
true;
1143 else if (ch == EOF) {
1146 else if (ch ==
'\n' || ch ==
'\t') {
1150 else if (ch ==
':' || ch ==
',' || ch ==
')' || ch ==
'[' || ch ==
'\n' || ch == EOF) {
1153 else if (ch ==
'_' || ch ==
'\t') {
1158 if (i <
nmlngth) str[i++] = ch;
1164 printf(
"ERROR: End-of-file in tree species name\n");
1168 (void) ungetc(ch, INFILE);
1169 while (i <
nmlngth) str[i++] =
' ';
1177 char *nameptr = q->
name;
1178 do { found = str[i] == *nameptr++; }
while (found && (++i <
nmlngth));
1183 }
while ((! found) && (++n <= tr->
mxtips));
1187 do { str[i] =
'\0'; }
while (i-- && (str[i] <=
' '));
1188 printf(
"ERROR: Cannot find data for tree species: %s\n", str);
1191 return found ? n : 0;
1197 if (ch != EOF) (void) ungetc(ch, INFILE);
1200 if (fscanf(INFILE,
"%lf", &branch) != 1) {
1201 printf(
"ERROR: Problem reading branch length in processLength:\n");
1204 if (fscanf(INFILE,
"%40s", str) == 1) printf(
"%s\n", str);
1218 else if (ch != EOF) (void) ungetc(ch, INFILE);
1222 static void treeNeedCh(
int c1,
const char *where, FILE *INFILE) {
1225 printf(
"ERROR: Missing '%c' %s tree; ", c1, where);
1227 printf(
"End-of-File");
1231 for (
int i = 24; i-- && (c2 != EOF); c2 = getc(INFILE)) putchar(c2);
1234 printf(
" found instead\n");
1244 if (n > 2*(tr->
mxtips) - 2) {
1246 printf(
"ERROR: Too many internal nodes. Is tree rooted?\n");
1247 printf(
" Deepest splitting should be a trifurcation.\n");
1267 if (n <= 0) {
anerror =
true;
return; }
1283 printf(
"ERROR: Unable to uproot tree.\n");
1284 printf(
" Inappropriate node marked for removal.\n");
1290 printf(
"ERROR: Unable to uproot tree. Inconsistent\n");
1291 printf(
" number of tips and nodes for rooted tree.\n");
1304 if (tr->
ntips > 2 && p != q && p != r && p != s) {
1340 p->next->next->back =
NULp;
1342 if (ch != EOF) (void) ungetc(ch, INFILE);
1347 printf(
"(This error also happens if the last species in the tree is unmarked)\n");
1372 while ((! p->
x) || (! q->
x)) {
1391 double xvlz =
xv * lz;
1395 double *log_f = tr->
log_f;
1399 for (
int i = 0; i <
endsite; i++) {
1400 double fx1a =
freqa * *x1a++;
1401 double fx1g =
freqg * *x1g++;
1402 double fx1c =
freqc * *x1c++;
1403 double fx1t =
freqt * *x1t++;
1405 double suma = fx1a * *x2a + fx1c * *x2c + fx1g * *x2g + fx1t * *x2t;
1407 double fx2r =
freqa * *x2a++ +
freqg * *x2g++;
1408 double fx2y =
freqc * *x2c++ +
freqt * *x2t++;
1409 double fx1r = fx1a + fx1g;
1410 double fx1y = fx1c + fx1t;
1412 double sumc = (fx1r + fx1y) * (fx2r + fx2y);
1418 double ki = *rptr++;
1419 double zz = exp(ki * lz);
1420 double zv = exp(ki * xvlz);
1422 double term = log(zz * suma + zv * sumb + sumc);
1423 sum += *wptr++ * term;
1436 while ((! p->
x) || (! q->
x)) {
1455 double xvlz =
xv * lz;
1460 for (
int i = 0; i <
endsite; i++) {
1461 double fx1a =
freqa * *x1a++;
1462 double fx1g =
freqg * *x1g++;
1463 double fx1c =
freqc * *x1c++;
1464 double fx1t =
freqt * *x1t++;
1466 double suma = fx1a * *x2a + fx1c * *x2c + fx1g * *x2g + fx1t * *x2t;
1468 double fx2r =
freqa * *x2a++ +
freqg * *x2g++;
1469 double fx2y =
freqc * *x2c++ +
freqt * *x2t++;
1470 double fx1r = fx1a + fx1g;
1471 double fx1y = fx1c + fx1t;
1473 double sumc = (fx1r + fx1y) * (fx2r + fx2y);
1479 double ki = *rptr++;
1481 suma *= exp(ki * lz);
1482 sumb *= exp(ki * xvlz);
1484 dLidki[i] += *wptr++ * lz * (suma + sumb*
xv);
1499 if (ki_min <= 0.0 || ki_max <= ki_min) {
1500 printf(
"ERROR: Bad rate value limits to findSiteRates\n");
1504 else if (d_ki <= 1.0) {
1505 printf(
"ERROR: Bad rate step to findSiteRates\n");
1510 for (
int i = 0; i <
endsite; i++) {
1515 for (
double ki = ki_min; ki <= ki_max; ki *= d_ki) {
1520 for (
int i = 0; i <
endsite; i++) {
1532 while (d_ki > 1.0 + max_error) {
1541 double inv_d_ki = 1.0/d_ki;
1542 for (
int i = 0; i <
endsite; i++) {
1589 double min_1 = 1.0E37;
1590 double min_2 = 1.0E37;
1594 for (
int i = 1; i <= Sites; i++) {
1595 if (Weight[i] > 0) {
1596 double ki = Patrate[Pattern[i]];
1599 if (ki < 0.995 * min_1) min_2 = min_1;
1602 else if (ki > 1.005 * min_1) {
1606 else if (ki > max_2) {
1608 if (ki > 1.005 * max_1) max_2 = max_1;
1611 else if (ki < 0.995 * max_1) {
1618 double a = (Categs - 3.0)/log(max_2/min_2);
1619 double b = - a * log(min_2) + 2.0;
1621 categrate[0] = min_1;
1622 for (
int k = 1; k <= Categs-2; k++) categrate[k] = min_2 * exp((k-1)/a);
1623 if (Categs>0) categrate[Categs-1] = max_1;
1625 for (
int i = 1; i <= Sites; i++) {
1626 if (Weight[i] > 0) {
1627 double ki = Patrate[Pattern[i]];
1628 if (ki < 0.99 * min_2) sitecateg[i] = 1;
1629 else if (ki > 1.00 * max_2) sitecateg[i] = Categs;
1630 else sitecateg[i] =
nint(a * log(Patrate[Pattern[i]]) + b);
1633 sitecateg[i] = Categs;
1653 if (!gb_sai_count) {
1671 printf(
"Writing '%s'\n", sai_name);
1687 char *cats = ARB_calloc<char>(ali_len+1);
1688 float *rates = ARB_calloc<float>(ali_len);
1689 char category_string[1024];
1693 long filter_len = strlen(arb_filter);
1694 if (filter_len != ali_len) {
1695 error =
GBS_global_string(
"Filter length (%li) does not match alignment length (%li)",
1696 filter_len, ali_len);
1708 for (
int ali_pos = 0; ali_pos < ali_len; ali_pos++) {
1709 if (arb_filter[ali_pos] ==
'0') {
1710 cats[ali_pos] =
'.';
1720 cats[ali_pos] =
'.';
1725 int unfiltered_sites = i-1;
1726 if (unfiltered_sites !=
sites) {
1727 error =
GBS_global_string(
"Filter positions (%i) do not match input sequence positions (%i)",
1728 unfiltered_sites,
sites);
1733 char *p = category_string;
1736 for (
int k = 1; k <=
categs; k ++) {
1737 sprintf(p,
" %G", categrate[
categs-k]);
1760 GB_write_string(gb_data,
"PVML: Positional Variability by ML (Olsen)");
1770 fprintf(stderr,
"Error in arb_dnarates: %s\n", error);
1777 gb_main =
GB_open(dbname,
"rw");
1779 if (strcmp(dbname,
":") == 0) {
1792 GB_warningf(
"Error saving '%s': %s", saveAs, error);
1808 for (
int k = 1; k <= Sites; k += 60) {
1809 int j =
min(k + 59, Sites);
1811 fprintf(outfile,
"%s ", k == 1 ?
"Weights " :
" ");
1813 for (
int i = k; i <= j; i++) {
1815 if (((i % 10) == 0) && ((i % 60) != 0)) putc(
' ', outfile);
1818 putc(
'\n', outfile);
1820 for (
int k = 1; k <= Categs; k += 7) {
1821 int j =
min(k + 6, Categs);
1823 if (k == 1) fprintf(outfile,
"C %2d", Categs);
1824 else fprintf(outfile,
" ");
1826 for (
int i = k-1; i < j; i++) fprintf(outfile,
" %9.5f", categrate[i]);
1828 putc(
'\n', outfile);
1831 for (
int k = 1; k <= Sites; k += 60) {
1832 int j =
min(k + 59, Sites);
1834 fprintf(outfile,
"%s ", k == 1 ?
"Categories" :
" ");
1836 for (
int i = k; i <= j; i++) {
1838 if (((i % 10) == 0) && ((i % 60) != 0)) putc(
' ', outfile);
1841 putc(
'\n', outfile);
1848 printf(
" Site Rate\n");
1849 printf(
" ---- ---------\n");
1851 for (
int i = 1; i <=
sites; i++) {
1856 printf(
"%6d Undefined\n", i);
1874 (void) sprintf(filename,
"%s.%d",
"weight_rate", getpid());
1877 (void) sprintf(filename,
"%s_%2d.%d",
"weight_rate", treenum, getpid());
1880 FILE *outfile = fopen(filename,
"w");
1883 (void) fclose(outfile);
1884 printf(
"Weights and categories also written to %s\n\n", filename);
1893 if (fscanf(INFILE,
"%d", &numtrees) != 1 ||
findch(
'\n', INFILE) == EOF) {
1894 printf(
"ERROR: Problem reading number of user trees\n");
1899 printf(
"User-defined %s:\n\n", (numtrees == 1) ?
"tree" :
"trees");
1901 for (
int which = 1; which <= numtrees; which++) {
1909 printf(
"%d taxon user-supplied tree read\n\n", tr->
ntips);
1910 printf(
"Total length of tree branches = %8.6f\n\n", tree_length);
1920 inline bool is_char(
const char *name,
char c) {
return name[0] == c && !name[1]; }
1925 const char *dbname =
":";
1926 const char *dbsavename =
NULp;
1929 const char *inputname =
NULp;
1930 FILE *infile =
NULp;
1933 case 3: error =
"'dbname' may only be used together with 'dbsavename'";
break;
1937 dbsavename = argv[3];
1940 inputname = argv[1];
1941 infile =
wantSTDIN(inputname) ? stdin : fopen(inputname,
"rt");
1942 if (!infile) error =
GB_IO_error(
"reading", inputname);
1946 case 1: error =
"missing arguments";
break;
1948 default : error =
"too many arguments";
break;
1952 fprintf(stderr,
"arb_dnarates: Error: %s\n", error);
1957 fputs(
"Usage: arb_dnarates input [dbname dbsavename]\n"
1958 " Expects phylip sequence data as 'input'.\n"
1959 " Reads from STDIN if 'input' is '-'.\n"
1960 " (e.g. cat data.phyl | arb_dnarates - ...\n"
1961 " or arb_dnarates data.phyl ...)\n"
1962 " Expects a 'dbname' or a running ARB DB server.\n"
1964 " - Saves calculated SAI to server (POS_VAR_BY_ML_...)\n"
1965 " Using 'dbname' does only make sense for unittests!\n"
1982 tree *tr = &curtree;
1995 if (
wantSTDIN(inputname)) fclose(infile);
bool wantSTDIN(const char *iname)
GB_ERROR GB_begin_transaction(GBDATA *gbd)
GBDATA * GB_open(const char *path, const char *opent)
static void empiricalfreqs(tree *tr)
static double evaluate(tree *tr, nodeptr p)
GB_ERROR GB_commit_transaction(GBDATA *gbd)
void GB_warning(const char *message)
static void getbasefreqs(FILE *INFILE)
long GB_read_int(GBDATA *gbd)
static void findSiteRates(tree *tr, double ki_min, double ki_max, double d_ki, double max_error)
GB_ERROR GB_save(GBDATA *gb, const char *path, const char *savetype)
static int patsite[maxsites+1]
GB_ERROR GB_write_string(GBDATA *gbd, const char *s)
static double bestki[maxpatterns]
static double dLidki[maxpatterns]
GB_ERROR GB_end_transaction(GBDATA *gbd, GB_ERROR error)
static void initrav(nodeptr p)
GB_ERROR GB_IO_error(const char *action, const char *filename)
static void treeReadLen(tree *tr, FILE *INFILE)
static void newview(nodeptr p)
const char * GBS_global_string(const char *templat,...)
long GBT_get_alignment_len(GBDATA *gb_main, const char *aliname)
static char * alignment_name
static void hookup(nodeptr p, nodeptr q, double z)
static void treeNeedCh(int c1, const char *where, FILE *INFILE)
static int treeFinishCom(FILE *INFILE)
static xarray * getxtip(nodeptr p)
static int get_next_SAI_count()
#define AWAR_GDE_EXPORT_FILTER
static double processLength(FILE *INFILE)
static void uppercase(int *chptr)
static void openArb(const char *dbname)
static void saveArb(const char *saveAs)
static void getArbFilter()
int ARB_main(int argc, char *argv[])
GBDATA * GBT_find_SAI(GBDATA *gb_main, const char *name)
NOT4PERL GBDATA * GBT_add_data(GBDATA *species, const char *ali_name, const char *key, GB_TYPES type) __ATTR__DEPRECATED_TODO("better use GBT_create_sequence_data()")
CONSTEXPR_INLINE int leafs_2_nodes(int leafs, TreeModel model)
static void getdata(tree *tr, FILE *INFILE)
GB_ERROR GB_await_error()
static void makevalues(tree *tr)
char * GBT_read_string(GBDATA *gb_container, const char *fieldpath)
static int weight[maxsites+1]
static void setuptree(tree *tr, const int numSp)
void GB_warningf(const char *templat,...)
static xarray * getxnode(nodeptr p)
static xarray * setupxarray()
bool is_char(const char *name, char c)
static void treeFlushLen(FILE *INFILE)
static double patrate[maxpatterns]
static int treeGetCh(FILE *INFILE)
static void error(const char *msg)
static void uprootTree(tree *tr, nodeptr p)
static void summarize(int treenum)
static void wrfile(FILE *outfile, int Sites, int Categs, int Weight[], double categrate[], int sitecateg[])
static void inputweights(FILE *INFILE)
static void freeTreeNode(nodeptr p)
static void makeUserRates(tree *tr, FILE *INFILE)
static int base36(int ch)
static void linkxarray(int req, int min, xarray **freexptr, xarray **usedxptr)
static void dli_dki(nodeptr p)
static GBDATA * create_next_SAI()
GB_ERROR GB_write_int(GBDATA *gbd, long i)
static void categorize(int Sites, int Categs, int Weight[], int Pattern[], double Patrate[], double categrate[], int sitecateg[])
static double bestLi[maxpatterns]
static double subtreeLength(nodeptr p)
static void freeTree(tree *tr)
fputs(TRACE_PREFIX, stderr)
static void makeweights()
static int findch(int c, FILE *INFILE)
static void getoptions(FILE *INFILE)
static int pattern[maxsites+1]
static int itobase36(int i)
static double treeLength(tree *tr)
static void spanSubtree(nodeptr p)
GB_ERROR GB_write_floats(GBDATA *gbd, const float *f, long size)
static void sitecombcrunch()
char * GBT_get_default_alignment(GBDATA *gb_main)
static void treeFlushLabel(FILE *INFILE)
static void getnums(FILE *INFILE)
double log_f[maxpatterns]
GBDATA * GBT_find_or_create_SAI(GBDATA *gb_main, const char *name)
static void getinput(tree *tr, FILE *INFILE)
GBDATA * GB_search(GBDATA *gbd, const char *fieldpath, GB_TYPES create)
static int patweight[maxsites+1]
static void setupnodex(tree *tr)
static int info[maxsites+1]
static int findTipName(tree *tr, int ch, FILE *INFILE)
static void addElementLen(tree *tr, nodeptr p, FILE *INFILE)
void GB_close(GBDATA *gbd)
GB_write_int const char s