1 #include "fluctuate.h"
2
3 #ifndef MODELLIKE_INCLUDE
4 #include "fluc_modellike.h"
5 #endif
6
7 #ifdef DMALLOC_FUNC_CHECK
8 #include "/usr/local/include/dmalloc.h"
9 #endif
10
11 /***************************************************************************
12 * ALPHA *
13 * version 1.40. (c) Copyright 1986, 1991, 1992 by the University of *
14 * Washington and Joseph Felsenstein. Written by Joseph Felsenstein, *
15 * Mary K. Kuhner and Jon A. Yamato, with some additional grunt work by *
16 * Sean T. Lamont. Permission is granted to copy and use this program *
17 * provided no fee is charged for it and provided that this copyright *
18 * notice is not removed. *
19 * *
20 ***************************************************************************/
21
22 /* This program implements a Hastings-Metropolis Monte Carlo
23 Maximum Likelihood sampling method for phylogenetic trees
24 without recombination. */
25
26 /* Time, 'tyme', in this tree is measured from the tips to the root.
27 I.e. the tips are at tyme '0', and the root node has the largest
28 value for 'tyme'. */
29
30 /* Only long chains included in joint estimates */
31
32 /* Modified to accept multiple loci. */
33
34 extern double **growi, **theti, **lntheti, xinterval;
35 extern treerec ***sum;
36
37 FILE *infile, *outfile, *treefile, *bestree, *simlog,
38 *parmfile, *thetafile, *intree;
39 long numseq, numnodes, sites, rootnum, categs, apps, inseed, col,
40 slidenum, numloci, numtrees, locus, totchains, holding;
41 boolean **sametree;
42 tree *curtree;
43 double locus_ttratio, freqa, freqc, freqg, freqt,
44 lambda, probsum, theta0, growth0, branch0;
45 long *category,*weight;
46 double *rate, *ne_ratio, *mu_ratio, *theta_ratio;
47 double *probcat;
48 contribarr *contribution;
49 node **freenodes;
50 node **slidenodes;
51 rlrec *alogf;
52 char gch;
53 dnadata *dna;
54 option_struct *op;
55 longer seed;
56
57 /* the following are used in site aliasing (makesiteptr) */
58 long *siteptr;
59
60 /* the following are for reading in parameters (readparmfile) */
61 char *booltokens[NUMBOOL] = {"interleaved","printdata","progress",
62 "print-trees","freqs-from-data","categories",
63 "watterson", "usertree", "autocorrelation",
64 "interactive"},
65 *numbertokens[NUMNUMBER] = {"ttratio","short-chains",
66 "short-steps","short-inc","long-chains",
67 "long-inc","long-steps","growth-rate"};
68
69 /* "Nearly global" variables for maketree: */
70
71 double sumweightrat, oldlikelihood, watttheta;
72 double *weightrat;
73 valrec *tbl;
74 long slid, slacc, indecks, chaintype;
75
76
openfile(FILE ** fp,char * filename,char * mode,char * application,char * perm)77 void openfile(FILE **fp, char *filename, char *mode, char *application,
78 char *perm)
79 {
80 FILE *of;
81 char file[100];
82
83 strcpy(file,filename);
84 while (1){
85 of = fopen(file,mode);
86 if (of)
87 break;
88 else {
89 switch (*mode){
90 case 'r':
91 fprintf(stdout,"%s: can't read %s\n",application,file);
92 file[0] = '\0';
93 while (file[0] =='\0'){
94 fprintf(stdout,"Please enter a new filename>");
95 fgets(file,100,stdin);
96 }
97 break;
98 case 'w':
99 fprintf(stdout,"%s: can't write %s\n",application,file);
100 file[0] = '\0';
101 while (file[0] =='\0'){
102 fprintf(stdout,"Please enter a new filename>");
103 fgets(file,100,stdin);
104 }
105 break;
106 default:
107 fprintf(ERRFILE,"%s: unknown file mode for %s\n",
108 application, file);
109 exit(-1);
110 break;
111 }
112 }
113 }
114 *fp=of;
115 if (perm != NULL)
116 strcpy(perm,file);
117 } /* openfile */
118
randum()119 double randum()
120 /* Mary's version--faster but needs 32 bits. Loops have been unrolled
121 for speed. */
122 {
123 long newseed0, newseed1, newseed2;
124
125 newseed0 = 1549*seed[0];
126 newseed1 = newseed0/2048;
127 newseed0 &= 2047;
128 newseed2 = newseed1/2048;
129 newseed1 &= 2047;
130 newseed1 += 1549*seed[1]
131 + 812*seed[0];
132 newseed2 += newseed1/2048;
133 newseed1 &= 2047;
134 newseed2 += 1549*seed[2]
135 + 812*seed[1];
136
137 seed[0] = newseed0;
138 seed[1] = newseed1;
139 seed[2] = newseed2 & 1023;
140 return (((seed[0]/2048.0 + seed[1])/2048.0 + seed[2])/1024.0);
141 } /* randum */
142
readseedfile(long * inseedptr)143 boolean readseedfile(long *inseedptr)
144 /* read in the seed from the seed file */
145 {
146 FILE *seedfile;
147 boolean success;
148
149 seedfile = fopen("seedfile","r");
150 if (seedfile) {
151 success = fscanf(seedfile, "%ld%*[^\n]", inseedptr);
152 fclose(seedfile);
153 return(success);
154 }
155 return(false); /* didn't find seedfile */
156 } /* readseedfile */
157
printtymelist()158 void printtymelist()
159 /* prints the entire tymelist: a debug function */
160 {
161 long i;
162 tlist *t;
163 long limit;
164
165 t = curtree->tymelist;
166 fprintf(ERRFILE,"TYMELIST BEGINS\n");
167 while (t != NULL) {
168 fprintf(ERRFILE,"%3ld age%8.6f branches %3ld--",
169 t->eventnode->number, t->age, t->numbranch);
170 limit = t->numbranch;
171 for (i = 0; i < limit; i++)
172 fprintf(ERRFILE," %3ld", t->branchlist[i]->number);
173 fprintf(ERRFILE,"\n");
174 t = t->succ;
175 }
176 fprintf(ERRFILE,"TYMELIST ENDS\n");
177 } /* printtymelist */
178
179 /* The next two functions are utility function to rescale time
180 when growth is non-zero */
181
to_real(double tyme)182 double to_real (double tyme)
183 /* convert magic time to real time */
184 {
185 double result;
186
187 /* first convert to Joe-Time */
188 tyme *= -1.0;
189 /* then convert to real time */
190 result = -log(1.0-growth0*tyme)/growth0;
191 /* now convert back to Mary-Time */
192 result *= -1.0;
193
194 /* if we're, theoretically, at infinity,
195 don't calculate f(infinity), just return it! */
196 if(tyme >= curtree->root->tyme) return(curtree->root->tyme);
197
198 return(result);
199 } /* to_real */
200
to_magic(double tyme)201 double to_magic (double tyme)
202 /* convert real time to magic tyme */
203 {
204 double result;
205
206 /* if we're, theoretically, at infinity,
207 don't calculate f(infinity), just return it! */
208 if(tyme >= curtree->root->tyme) return(curtree->root->tyme);
209
210 /* first convert to Joe-Time */
211 tyme *= -1.0;
212 /* then convert to magic time */
213 result = (1.0-exp(-growth0*tyme))/growth0;
214 /* now convert back to Mary-Time */
215 result *= -1.0;
216
217 return(result);
218 } /* to_magic */
219
220 /* end time-rescalers */
221
222 /*****************************************************************
223 * lengthof() returns the length of the branch "rootward" of the *
224 * passed node */
lengthof(node * p)225 double lengthof(node *p)
226 {
227 return fabs(p->tyme - p->back->tyme);
228 } /* lengthof */
229
230 /*********************************************************************
231 * howlong() returns the amount of tyme passing between 2 entries in *
232 * the tymelist. Magical tyme is returned if growth != 0. */
howlong(tlist * t)233 double howlong(tlist *t)
234 {
235 double tlength, tyme1, tyme2;
236
237 /* don't try this if growth rate is 0 */
238 if (!growth0 || !op->growthused) return(t->age - t->eventnode->tyme);
239
240 /* don't bother if the answer will be zero anyway */
241 if (t->age==0.0) return(0.0);
242 /* if we're on the theoretically infinitely long branch,
243 don't calculate f(infinity), just return it! */
244 if (t->age >= curtree->root->tyme) return(rootlength);
245
246 /* otherwise transform to 'magical: fixed theta' time */
247
248 tyme1 = to_magic(t->age);
249 tyme2 = to_magic(t->eventnode->tyme);
250 tlength = tyme1 - tyme2;
251
252 return(tlength);
253 } /* howlong */
254
findtop(node * p)255 node *findtop(node *p)
256 /* findtop returns the first 'top' it finds, for a given node the _same_
257 'top' will always be found */
258 {
259 while (!p->top)
260 p = p->next;
261 return(p);
262 } /* findtop */
263
VarMalloc(node * p,boolean allokate)264 void VarMalloc(node *p, boolean allokate)
265 /* callocs or frees the 'x' [basewise Lnlikelihood] field of a node */
266 {
267 long i;
268
269 if (allokate) {
270 if (p->x == NULL) {
271 p->x = (phenotype)calloc(sites,sizeof(ratelike));
272 p->x[0] = (ratelike)calloc(sites*categs,sizeof(sitelike));
273 for(i = 1; i < sites; i++)
274 p->x[i] = p->x[0] + i * categs;
275 }
276 }
277 else {
278 if (p->x != NULL) {
279 free(p->x[0]);
280 free(p->x);
281 }
282 p->x = NULL;
283 }
284 } /* VarMalloc */
285
286 /* "newnode" & "freenode" are paired memory managers for tree nodes.
287 They maintain a linked list of unused nodes which should be faster
288 to use than "calloc" & "free" (at least for recombination) */
newnode(node ** p)289 void newnode(node **p)
290 {
291 long i;
292
293 i = 0;
294 while (i < numnodes + 3) {
295 if (freenodes[i] != NULL) {
296 *p = freenodes[i];
297 freenodes[i] = NULL;
298 return;
299 }
300 i++;
301 }
302 fprintf(ERRFILE,"newnode failed!\n");
303 exit(-1);
304 } /* newnode */
305
freenode(node * p)306 void freenode(node *p)
307 {
308 long i;
309
310 i = 0;
311 while (i < numnodes + 3) {
312 if (freenodes[i] == NULL) {
313 freenodes[i] = p;
314 return;
315 }
316 i++;
317 }
318 fprintf(ERRFILE,"freenode failed!\n");
319 exit(-1);
320 } /* freenode */
321 /* END of treenode allocation functions */
322
newtymenode(tlist ** t)323 void newtymenode(tlist **t)
324 {
325 *t = (tlist *)calloc(1,sizeof(tlist));
326 (*t)->branchlist = (node **)calloc(numseq,sizeof(node *));
327 } /* newtymenode */
328
freetymenode(tlist * t)329 void freetymenode(tlist *t)
330 {
331 free(t->branchlist);
332 free(t);
333 } /* freetymenode*/
334
freetymelist(tlist * t)335 void freetymelist(tlist *t)
336 {
337 if (t->succ != NULL) freetymelist(t->succ);
338 freetymenode(t);
339 } /* freetymelist */
340
hookup(node * p,node * q)341 void hookup(node *p, node *q)
342 {
343 p->back = q;
344 q->back = p;
345 } /* hookup */
346
atr(node * p)347 void atr(node *p)
348 /* "atr" prints out a text representation of a tree. Pass
349 curtree->root->back for normal results */
350 {
351 if (p->back == curtree->root) {
352 fprintf(ERRFILE,"next node is root\n");
353 fprintf(ERRFILE,"Node %4ld length %12.6f tyme %15.10f",
354 p->back->number, lengthof(p), p->back->tyme);
355 fprintf(ERRFILE," --> %4ld\n",p->number);
356 }
357 fprintf(ERRFILE,"Node %4ld length %12.6f tyme %23.18f -->",
358 p->number, lengthof(p), p->tyme);
359 if (p->top && p->back->top) fprintf(ERRFILE,"TWO TOPS HERE!!!!");
360 if (!p->tip) {
361 if (!p->next->top) fprintf(ERRFILE,"%4ld",p->next->back->number);
362 if (!p->next->next->top) fprintf(ERRFILE,"%4ld",
363 p->next->next->back->number);
364 fprintf(ERRFILE,"\n");
365 if (!p->next->top) atr(p->next->back);
366 if (!p->next->next->top)
367 atr(p->next->next->back);
368 }
369 else fprintf(ERRFILE,"\n");
370 } /* atr */
371
372 /* The next set of functions [zerocollis/onecollis/twocollis] compute
373 the chance that there are 0/1/2 coalescences [respectively]
374 in an interval of length "length", with "numl" active lineages, and
375 "numother" inactive lineages */
zerocollis(long numl,long numother,double length)376 double zerocollis(long numl, long numother, double length)
377 {
378 double result;
379
380 result = exp(-(numl * (numl - 1) + numl * numother * 2) * (length / theta0));
381
382 if (result == 0.0 && length != 0.0) result = EPSILON;
383
384 return(result);
385 } /* zerocollis */
386
387
onecollis(long numl,long numother,double length)388 double onecollis(long numl, long numother, double length)
389 {
390 double expon1, expon2, result;
391
392 expon1 = -((numl - 1) * numother * 2 + (numl - 1) * (numl - 2)) *
393 (length / theta0);
394 expon2 = -(numl * numother * 2 + numl * (numl - 1)) * (length / theta0);
395
396 result = (numl * (numl - 1.0) / (numother * 2 + (numl - 1) * 2) *
397 (exp(expon1) - exp(expon2)));
398
399 if (result == 0.0 && length != 0.0) result = EPSILON;
400
401 return(result);
402 } /* onecollis */
403
404
twocollis(long numother,double length)405 double twocollis(long numother, double length)
406 /* For this case "numl" is assumed to be equal to 3 */
407 {
408 double expon1, expon2, expon3, result;
409
410 expon1 = numother * -2 * (length / theta0);
411 expon2 = -(numother * 4 + 2.0) * (length / theta0);
412 expon3 = -(numother * 6 + 6.0) * (length / theta0);
413
414 result = (6.0 / (numother + 1) *
415 (1.0 / (numother * 4 + 6) * (exp(expon1) - exp(expon3)) -
416 1.0 / (numother * 2 + 4) * (exp(expon2) - exp(expon3))));
417
418 if (result == 0.0 && length != 0.0) result = EPSILON;
419
420 return(result);
421 } /* twocollis */
422 /* End of coalescence functions */
423
gettymenode(long target)424 tlist *gettymenode(long target)
425 /* Return a pointer to the tymelist entry whose 'eventnode' has
426 the number of 'target'. */
427 {
428 tlist *t;
429
430 if (target == curtree->root->number) {
431 return(NULL);
432 }
433 t = curtree->tymelist;
434 if (target <= numseq) return(t); /* it's a tip! */
435 while (t != NULL) {
436 if (t->eventnode->number == target) return(t);
437 t = t->succ;
438 }
439 fprintf(ERRFILE,"In gettymenode, failed to find node%12ld\n", target);
440 fprintf(ERRFILE,"CATASTROPHIC ERROR\n");
441 exit(-1);
442 } /* gettymenode */
443
gettyme(node * p,node * daughter1,node * daughter2,node * ans)444 tlist *gettyme(node *p, node *daughter1, node *daughter2,
445 node *ans)
446 /* Return a pointer to the tymelist entry which encompasses the time
447 into which you wish to place node "p".
448 tipwards/upper bound: "daughter1" and "daughter2"
449 rootward/lower bound: "ans" */
450 {
451 boolean found;
452 tlist *t, *b1, *b2, *before, *after;
453
454 /* first establish a tipward bound on the search */
455 before = curtree->tymelist;
456 found = false;
457 b1 = gettymenode(daughter1->number);
458 b2 = gettymenode(daughter2->number);
459 while (true) {
460 if ((before == b1) || (before == b2)) {
461 if ((found) || (b1 == b2)) break;
462 found = true;
463 }
464 before = before->succ;
465 }
466 /* now establish a rootward bound on the search */
467 after = gettymenode(ans->number);
468 /* begin the search at the tipward bound */
469 t = before;
470 found = false;
471 while (t != after && !found) {
472 if (t->age >= p->tyme)
473 found = true;
474 else
475 t = t->succ;
476 }
477 if (!found)
478 /* prime^.tyme is tied with after, so goes directly in front of it */
479 t = t->prev;
480 return(t);
481 } /* gettyme */
482
inserttymelist(node * prime)483 void inserttymelist(node *prime)
484 /* inserts 2 entries into the tymelist:
485 "prime" and "primans" [prime->back]. */
486 {
487 tlist *t, *temp;
488 /* d[2] is primans' daughter, d[0] and [1] parent's daughters. */
489 node *parent, *q, *primans, *d[3];
490 long i, j;
491 d[0] = d[1] = d[2] = parent = primans = NULL; /* just to be careful */
492
493 newtymenode(&t);
494 /* find daughters and parents */
495 /* this complicated mess is needed because prime must be correctly
496 bounded, not by primans (which is not yet in the tymelist), but by
497 the parent of primans */
498 q = prime;
499 j = 0;
500 for (i = 1; i <= 3; i++) {
501 if (q->top) {
502 primans = q->back;
503 if (primans->next->top) {
504 parent = primans->next->back;
505 d[2] = primans->next->next->back;
506 } else {
507 parent = primans->next->next->back;
508 d[2] = primans->next->back;
509 }
510 } else {
511 d[j] = q->back;
512 j++;
513 }
514 q = q->next;
515 }
516 /* insert prime */
517 t->eventnode = prime;
518 temp = gettyme(prime, d[0], d[1], parent);
519 t->succ = temp->succ;
520 t->prev = temp;
521 if (temp->succ != NULL)
522 temp->succ->prev = t;
523 temp->succ = t;
524 if (t->succ != NULL)
525 t->age = t->succ->eventnode->tyme;
526 else
527 t->age = t->prev->age;
528 t->prev->age = t->eventnode->tyme;
529 /* insert primans */
530 newtymenode(&t);
531 t->eventnode = primans;
532 temp = gettyme(primans, prime, d[2], parent);
533 t->succ = temp->succ;
534 t->prev = temp;
535 if (temp->succ != NULL)
536 temp->succ->prev = t;
537 temp->succ = t;
538 if (t->succ != NULL)
539 t->age = t->succ->eventnode->tyme;
540 else
541 t->age = t->prev->age;
542 t->prev->age = t->eventnode->tyme;
543 } /* inserttymelist */
544
subtymelist(node * ndonor,node * nrecip)545 void subtymelist(node *ndonor, node *nrecip)
546 /* takes out 2 entries from the tymelist:
547 "ndonor" and "nrecip" [which must be tipward/above ndonor] */
548 {
549 long i, j, limit;
550 tlist *d, *r, *t;
551 node *badbranch, *p;
552 boolean found;
553
554 badbranch = NULL; /* just to be careful */
555
556 i = 0;
557
558 r = gettymenode(nrecip->number);
559 d = gettymenode(ndonor->number);
560 p = nrecip;
561 for (j = 1; j <= 3; j++) {
562 p = p->next;
563 if (p->back->number == ndonor->number)
564 badbranch = p;
565 }
566
567 t = r;
568
569 while (t != d) {
570 limit = t->numbranch;
571 for (i = 1; i <= limit; i++) {
572 if (t->branchlist[i - 1] == badbranch) {
573 j = i;
574 t->numbranch--;
575 }
576 }
577 for (i = j; i <= t->numbranch; i++)
578 t->branchlist[i - 1] = t->branchlist[i];
579 t = t->succ;
580 }
581 p = ndonor;
582 p = findtop(p);
583 badbranch = p;
584 found = true;
585 while (t != NULL && found) {
586 found = false;
587 for (i = 1; i <= t->numbranch; i++) {
588 if (t->branchlist[i - 1] == badbranch) {
589 j = i;
590 t->numbranch--;
591 found = true;
592 }
593 }
594 if (found) {
595 for (i = j; i <= t->numbranch; i++)
596 t->branchlist[i - 1] = t->branchlist[i];
597 }
598 t = t->succ;
599 }
600 r->prev->succ = r->succ;
601 r->succ->prev = r->prev;
602 r->prev->age = r->age;
603 freetymenode(r);
604 d->prev->succ = d->succ;
605 if (d->succ != NULL)
606 d->succ->prev = d->prev;
607 d->prev->age = d->age;
608 freetymenode(d);
609 } /* subtymelist */
610
ltov(node * p)611 void ltov(node *p)
612 /* ltov recalculates the proper "v" value of a branch, from
613 the tymes at either end of the branch */
614 {
615 p->v = 1.0 - exp(-(lengthof(p) / dna->fracchange));
616 p->back->v = p->v;
617 } /* ltov */
618
getnums()619 void getnums()
620 {
621 /* input number of sequences, number of sites */
622 fprintf(outfile, "\n");
623 fscanf(infile, "%ld%ld", &numseq, &sites);
624 fprintf(outfile, "%4ld Sequences, %4ld Sites\n", numseq, sites);
625 numnodes = numseq * 2 - 1; /*number of nodes in tree, excluding root*/
626 rootnum = numnodes + 3;
627 setupdata(&dna, sites, numseq);
628 freenodes = (node **)calloc(2,sizeof(node *));
629 /* number of internal nodes in tree is numseq-1 */
630 slidenodes = (node **)calloc(numseq-1,sizeof(node *));
631 } /* getnums */
632
633 /* boolcheck(), booleancheck(), numbercheck(), and readparmfile()
634 are used in reading the parameter file "parmfile" */
boolcheck(char ch)635 long boolcheck(char ch)
636 {
637 ch = toupper(ch);
638 if (ch == 'F') return 0;
639 if (ch == 'T') return 1;
640 return -1;
641 } /* boolcheck */
642
booleancheck(char * var,char * value)643 boolean booleancheck(char *var, char *value)
644 {
645 long i, j, check;
646
647 check = boolcheck(value[0]);
648 if(check == -1) return false;
649
650 for(i = 0; i < NUMBOOL; i++) {
651 if(!strcmp(var,booltokens[i])) {
652 if(i == 0) op->interleaved = (boolean)(check);
653 if(i == 1) op->printdata = (boolean)(check);
654 if(i == 2) op->progress = (boolean)(check);
655 if(i == 3) op->treeprint = (boolean)(check);
656 if(i == 4) {
657 op->freqsfrom = (boolean)(check);
658 if(!op->freqsfrom) {
659 strtok(value,":");
660 freqa = (double)atof(strtok(NULL,";"));
661 freqc = (double)atof(strtok(NULL,";"));
662 freqg = (double)atof(strtok(NULL,";"));
663 freqt = (double)atof(strtok(NULL,";"));
664 }
665 }
666 if(i == 5) {
667 op->ctgry = (boolean)(check);
668 if(op->ctgry) {
669 strtok(value,":");
670 categs = (long)atof(strtok(NULL,";"));
671 rate = (double *)realloc(rate,categs*sizeof(double));
672 probcat = (double *)realloc(probcat,categs*sizeof(double));
673 for(j = 0; j < categs; j++) {
674 rate[j] = (double)atof(strtok(NULL,";"));
675 probcat[j] = (double)atof(strtok(NULL,";"));
676 }
677 }
678 }
679 if(i == 6) {
680 op->watt = (boolean)(check);
681 if (!op->watt) {
682 strtok(value,":");
683 theta0 = (double)atof(strtok(NULL,";"));
684 }
685 }
686 if(i == 7) op->usertree = (boolean)(check);
687 if(i == 8) {
688 op->autocorr = (boolean)(check);
689 if (op->autocorr) {
690 strtok(value,":");
691 lambda = 1.0 / (double)atof(strtok(NULL,";"));
692 }
693 }
694 if(i == 9) op->interact = (boolean)(check);
695 return true;
696 }
697 }
698 return false;
699 } /* booleancheck */
700
numbercheck(char * var,char * value)701 boolean numbercheck(char *var, char *value)
702 {
703 long i;
704
705 for(i = 0; i < NUMNUMBER; i++) {
706 if(!strcmp(var,numbertokens[i])) {
707 if(i == 0) locus_ttratio = atof(value);
708 if(i == 1) op->numchains[0] = atol(value);
709 if(i == 2) op->steps[0] = atol(value);
710 if(i == 3) op->increm[0] = atol(value);
711 if(i == 4) op->numchains[1] = atol(value);
712 if(i == 5) op->increm[1] = atol(value);
713 if(i == 6) op->steps[1] = atol(value);
714 if(i == 7) {growth0 = atof(value); op->growthused = true;}
715 return true;
716 }
717 }
718 return false;
719 } /* numbercheck */
720
readparmfile()721 void readparmfile()
722 {
723 char fileline[LINESIZE],parmvar[LINESIZE],varvalue[LINESIZE];
724
725 parmfile = fopen("parmfile","r");
726
727 if(parmfile) {
728 while(fgets(fileline, LINESIZE, parmfile) != NULL) {
729 if(fileline[0] == '#') continue;
730 if(!strncmp(fileline,"end",3)) break;
731 strcpy(parmvar,strtok(fileline,"="));
732 strcpy(varvalue,strtok(NULL,"\n"));
733 /* now to process... */
734 if(!booleancheck(parmvar,varvalue))
735 if(!numbercheck(parmvar,varvalue)) {
736 fprintf(ERRFILE,
737 "Inappropiate entry in parmfile: %s\n", fileline);
738 exit(-1);
739 }
740 }
741 } else
742 if (!menu) {
743 fprintf(simlog,"Parameter file (parmfile) missing\n");
744 exit(-1);
745 }
746 } /* readparmfile */
747 /* END parameter file read */
748
getoptions()749 void getoptions()
750 /* interactively set options using a very basic menu */
751 {
752 boolean done, done1, done2;
753 char ch;
754 long i, j;
755 char input[LINESIZE];
756
757 rate = (double *)calloc(1,sizeof(double));
758 probcat = (double *)calloc(1,sizeof(double));
759
760 /* first some multiple rate-categories code stuff */
761 op->ctgry = false;
762 rate[0] = 1.0;
763 probcat[0] = 1.0;
764 categs = 1;
765 lambda = 1.0;
766 op->autocorr = false; /* false if categs == 1 */
767 /* end categories code stuff */
768
769 fscanf(infile,"%ld",&numloci);
770 ne_ratio = (double *)calloc(numloci,sizeof(double));
771 for(i = 0; i < numloci; i++) ne_ratio[i] = 1.0;
772 op->same_ne = true;
773 mu_ratio = (double *)calloc(numloci,sizeof(double));
774 for(i = 0; i < numloci; i++) mu_ratio[i] = 1.0;
775 op->same_mu = true;
776 theta_ratio = (double *)calloc(numloci,sizeof(double));
777
778 /* default initializations */
779 holding = 0;
780 op->interleaved = false;
781 op->printdata = false;
782 op->progress = true;
783 op->treeprint = false;
784 locus_ttratio = 2.0;
785 op->freqsfrom = true;
786 op->watt = false;
787 op->usertree = false;
788 op->interact = false;
789 theta0 = 1.0;
790 op->growthused = false;
791 growth0 = 0.0;
792 op->numchains[0] = 10;
793 op->increm[0] = 20;
794 op->steps[0] = 200;
795 op->numchains[1] = 2;
796 op->increm[1] = 20;
797 op->steps[1] = 20000;
798 /* end defaults */
799
800 readparmfile();
801 fprintf(outfile, "\nFluctuating population size HMMC ML method");
802 fprintf(outfile, " method, version 1.4\n\n");
803 if (!readseedfile(&inseed)) {
804 if (!menu) {
805 fprintf(ERRFILE,"Failed to find seedfile, aborting.\n");
806 exit(-1);
807 } else {
808 printf("Random number seed (must be odd)?\n");
809 scanf("%ld%*[^\n]", &inseed);
810 getchar();
811 }
812 }
813 /* now insure that the randum number seed is of the form 4n+1 */
814 /* inseed = inseed * 4 + 1; */
815 for (i = 0; i <= 2; i++)
816 seed[i] = 0;
817 i = 0;
818 for (i = 0; i <= 2; i++) {
819 seed[i] = inseed & 2047;
820 inseed /= 2048;
821 if (inseed == 0) break;
822 i++;
823 }
824 if (!menu) {
825 for (i=0; i<numloci; i++)
826 theta_ratio[i] = ne_ratio[i] * mu_ratio[i];
827 return;
828 }
829 putchar('\n');
830 do {
831 printf("\n%s", op->ansi ? "\033[2J\033[H" : "\n");
832 printf("Hastings-Metropolis Markov Chain Monte Carlo");
833 printf(" method, version 1.4\n\n");
834 printf("INPUT/OUTPUT FORMATS\n");
835 printf(" I Input sequences interleaved? %s\n",
836 op->interleaved ? "Yes" : "No, sequential");
837 printf(" E Echo the data at start of run? %s\n",
838 op->printdata ? "Yes" : "No");
839 printf(" P Print indications of progress of run? %s\n",
840 op->progress ? "Yes" : "No");
841 printf(" G Print out genealogies? %s\n",
842 op->treeprint ? "Yes" : "No");
843 printf(" Q Allow interactive design of output? %s\n",
844 op->interact ? "Yes" : "No");
845 printf("MODEL PARAMETERS\n");
846 printf(" T Transition/transversion ratio:");
847 printf(" %8.4f\n",locus_ttratio);
848 printf(" F Use empirical base frequencies? %s\n",
849 op->freqsfrom ? "Yes" : "No");
850 printf(" C One category of substitution rates?");
851 if (!op->ctgry || categs == 1)
852 printf(" Yes\n");
853 else {
854 printf(" %ld categories\n", categs);
855 printf(" R Rates at adjacent sites correlated?");
856 if (!op->autocorr)
857 printf(" No, they are independent\n");
858 else
859 printf(" Yes, mean block length =%6.1f\n", 1.0 / lambda);
860 }
861 printf(" W Use Watterson estimate of theta?");
862 if (op->watt)
863 printf(" Yes\n");
864 else
865 printf(" No, initial theta = %6.4f\n", theta0);
866 if (op->growthused) {
867 printf(" H Population can change in size? Yes\n");
868 printf(" V Rate of change parameter for growth: %e\n",
869 growth0);
870 }
871 else
872 printf(" H Population can change in size? No\n");
873 printf(" A Allow parameters to change?");
874 if (holding == 0) printf(" Yes\n");
875 else if (holding == 1) printf(" No, theta fixed\n");
876 else if (holding == 2) printf(" No, growth-rate fixed\n");
877 else printf(" Unknown option!!!\n");
878 printf(" U Use user tree in file \"intree\" ? %s\n",
879 op->usertree ? "Yes" : "No, construct a random tree");
880 if (numloci > 1) {
881 printf("MULTIPLE LOCI\n");
882 printf(" Z Population size equal among loci?");
883 if (op->same_ne)
884 printf(" Yes\n");
885 else
886 printf(" No\n");
887 printf(" M Mutation rate equal among loci?");
888 if (op->same_mu)
889 printf(" Yes\n");
890 else
891 printf(" No\n");
892 }
893 printf("SEARCH STRATEGY\n");
894 printf(" S Number of short chains to run? %6ld\n", op->numchains[0]);
895 if (op->numchains[0] > 0) {
896 printf(" 1 Short sampling increment? %6ld\n",
897 op->increm[0]);
898 printf(" 2 Number of steps along short chains? %6ld\n",
899 op->steps[0]);
900 }
901 printf(" L Number of long chains to run? %6ld\n", op->numchains[1]);
902 if (op->numchains[1] > 0) {
903 printf(" 3 Long sampling increment? %6ld\n",
904 op->increm[1]);
905 printf(" 4 Number of steps along long chains? %6ld\n",
906 op->steps[1]);
907 }
908 printf("\n");
909 printf("Are these settings correct? (type Y or the letter for one to change)\n");
910 fgets(input,LINESIZE,stdin);
911 ch = input[0];
912 ch = toupper(ch);
913 done = (ch == 'Y');
914 if (!done) {
915 ch = toupper(ch);
916 if (strchr("AZFGHTIE1234CWVUSLPRQM",ch) != NULL){
917 switch (ch) {
918
919 case 'A':
920 holding += 1;
921 if (holding > 2) holding = 0;
922 break;
923
924 case 'Z':
925 op->same_ne = !op->same_ne;
926 if (!op->same_ne) {
927 printf("Enter relative population size for each locus");
928 printf(" in input order?\n");
929 i = 0;
930 do {
931 scanf("%lf", &ne_ratio[i]);
932 if (ne_ratio[i] <= 0.0) {
933 printf(" Ratios must be positive, please reenter\n");
934 } else i++;
935 } while (i < numloci);
936 } else for(i = 0; i < numloci; i++) ne_ratio[i] = 1.0;
937 break;
938
939 case 'M':
940 op->same_mu = !op->same_mu;
941 if (!op->same_mu) {
942 printf("Enter relative mutation rate for each locus");
943 printf(" in input order?\n");
944 i = 0;
945 do {
946 scanf("%lf", &mu_ratio[i]);
947 if (mu_ratio[i] <= 0.0) {
948 printf(" Ratios must be positive, please reenter\n");
949 } else i++;
950 } while (i < numloci);
951 } else for(i = 0; i < numloci; i++) mu_ratio[i] = 1.0;
952 break;
953
954 case 'S':
955 do {
956 printf("How many Short Chains?\n");
957 fgets(input,LINESIZE,stdin);
958 op->numchains[0] = atoi(input);
959 if (op->numchains[0] < 0)
960 printf("Must be non-negative\n");
961 } while (op->numchains[0] < 0);
962 break;
963
964 case 'H':
965 op->growthused = !op->growthused;
966 break;
967
968 case 'V':
969 if (!op->growthused) break;
970 printf("What parameter value for growth?\n");
971 fgets(input,LINESIZE,stdin);
972 growth0 = atof(input);
973 break;
974
975 case 'L':
976 do {
977 printf("How many Long Chains?\n");
978 fgets(input,LINESIZE,stdin);
979 op->numchains[1] = atoi(input);
980 if (op->numchains[1] < 0)
981 printf("Must be non-negative\n");
982 } while (op->numchains[1] < 0);
983 break;
984
985 case 'C':
986 op->ctgry = !op->ctgry;
987 if (!op->ctgry)
988 op->autocorr = false;
989 if (op->ctgry) {
990 do {
991 printf("Number of categories ?");
992 fgets(input,LINESIZE,stdin);
993 categs = atoi(input);
994 } while (categs < 1);
995 free(rate);
996 free(probcat);
997 printf("Rate for each category? (use a space to");
998 printf(" separate)\n");
999 rate = (double *)calloc(categs,sizeof(double));
1000 probcat = (double *)calloc(categs,sizeof(double));
1001 for (j = 0; j < categs; j++)
1002 scanf("%lf*[^\n]", &rate[j]);
1003
1004 getchar();
1005 do {
1006 printf("Probability for each category?");
1007 printf(" (use a space to separate)\n");
1008 for (j = 0; j < categs; j++)
1009 scanf("%lf", &probcat[j]);
1010 scanf("%*[^\n]");
1011 getchar();
1012 done2 = true;
1013 probsum = 0.0;
1014 for (j = 0; j < categs; j++)
1015 probsum += probcat[j];
1016 if (fabs(1.0 - probsum) > 0.001) {
1017 done2 = false;
1018 printf("Probabilities must add up to");
1019 printf(" 1.0, plus or minus 0.001.\n");
1020 }
1021 } while (!done2);
1022 }
1023 break;
1024
1025 case 'R':
1026 op->autocorr = !op->autocorr;
1027 if (op->autocorr) {
1028 do {
1029 printf("Mean block length of sites having the same ");
1030 printf("rate (greater than 1)?\n");
1031 scanf("%lf%*[^\n]", &lambda);
1032 getchar();
1033 } while (lambda <= 1.0);
1034 lambda = 1.0 / lambda;
1035 }
1036 break;
1037
1038 case 'F':
1039 op->freqsfrom = !op->freqsfrom;
1040 if (!op->freqsfrom) {
1041 printf("Base frequencies for A, C, G, T/U (use blanks to separate)?\n");
1042 scanf("%lf%lf%lf%lf", &freqa, &freqc, &freqg, &freqt);
1043 scanf("%*[^\n]");
1044 }
1045 break;
1046
1047 case 'T':
1048 do {
1049 printf("Transition/transversion ratio?\n");
1050 fgets(input,LINESIZE,stdin);
1051 locus_ttratio = atof(input);
1052 } while (locus_ttratio < 0.0);
1053 break;
1054
1055 case 'I':
1056 op->interleaved = !op->interleaved;
1057 break;
1058
1059 case 'W':
1060 op->watt = !op->watt;
1061 if (!op->watt) {
1062 do {
1063 printf("Initial theta estimate?\n");
1064 fgets(input,LINESIZE,stdin);
1065 theta0 = atof(input);
1066 } while (theta0 <= 0.0);
1067 }
1068 break;
1069
1070 case 'U':
1071 op->usertree = !op->usertree;
1072 break;
1073
1074 case 'E':
1075 op->printdata = !op->printdata;
1076 break;
1077
1078 case 'P':
1079 op->progress = !op->progress;
1080 break;
1081
1082 case 'G':
1083 op->treeprint = !op->treeprint;
1084 break;
1085
1086 case 'Q':
1087 op->interact = !op->interact;
1088 break;
1089
1090 case '1':
1091 done1 = false;
1092 while (!done1) {
1093 printf("How often to sample trees?\n");
1094 fgets(input,LINESIZE,stdin);
1095 op->increm[0] = atoi(input);
1096 if (op->increm[0] > 0)
1097 done1 = true;
1098 else
1099 printf("Must be positive\n");
1100 }
1101 break;
1102
1103 case '2':
1104 done1 = false;
1105 while (!done1) {
1106 printf("How many short steps?\n");
1107 fgets(input,LINESIZE,stdin);
1108 op->steps[0] = atoi(input);
1109 if (op->steps[0] > 0)
1110 done1 = true;
1111 else
1112 printf("Must be a positive integer\n");
1113 }
1114 break;
1115
1116 case '3':
1117 done1 = false;
1118 while (!done1) {
1119 printf("How often to sample trees?\n");
1120 fgets(input,LINESIZE,stdin);
1121 op->increm[1] = atoi(input);
1122 if (op->increm[1] > 0)
1123 done1 = true;
1124 else
1125 printf("Must be positive\n");
1126 }
1127 break;
1128
1129 case '4':
1130 done1 = false;
1131 while (!done1) {
1132 printf("How many long steps?\n");
1133 fgets(input,LINESIZE,stdin);
1134 op->steps[1] = atoi(input);
1135 if (op->steps[1] > 0)
1136 done1 = true;
1137 else
1138 printf("Must be a positive integer\n");
1139 }
1140 break;
1141
1142 default:
1143 fprintf(stderr,"Impossible option %c detected!\n",ch);
1144 break;
1145
1146 }
1147 } else
1148 printf("Not a possible option!\n");
1149 }
1150 } while (!done);
1151 for (i=0; i<numloci; i++) {
1152 theta_ratio[i] = ne_ratio[i] * mu_ratio[i];
1153 }
1154 } /* getoptions */
1155
firstinit()1156 void firstinit()
1157 /* initialization for things that are recorded over multiple loci */
1158 {
1159 long i;
1160
1161 totchains = op->numchains[0] + op->numchains[1];
1162
1163 numtrees = MAX(op->steps[0]/op->increm[0],op->steps[1]/op->increm[1]);
1164
1165 sametree = (boolean **)calloc(numloci,sizeof(boolean *));
1166 sametree[0] = (boolean *)calloc(numloci*numtrees,sizeof(boolean));
1167 for(i = 1; i < numloci; i++)
1168 sametree[i] = sametree[0] + i*numtrees;
1169
1170 model_alloc();
1171
1172 } /* firstinit */
1173
locusinit()1174 void locusinit()
1175 /* initialization of things that are specific to one locus */
1176 {
1177 long i, j;
1178
1179 getnums();
1180
1181 if ((op->increm[0] < 0) || (op->increm[1] < 0)) {
1182 fprintf(ERRFILE,"Error in input sampling increment");
1183 fprintf(ERRFILE," increment set to 10\n");
1184 if (op->increm[0] < 0)
1185 op->increm[0] = 10;
1186 if (op->increm[1] < 0)
1187 op->increm[1] = 10;
1188 }
1189
1190 for (i = 0; i < 1+op->numchains[1]; i++)
1191 for (j = 0; j < numtrees; j++) {
1192 sum[locus][i][j].kk =
1193 (long *)calloc(numseq,sizeof(long));
1194 sum[locus][i][j].kend =
1195 (double *)calloc(numseq,sizeof(double));
1196 }
1197
1198 } /* locusinit */
1199
inputoptions()1200 void inputoptions()
1201 {
1202 char ch;
1203 long i, extranum;
1204
1205 category = (long *)calloc(sites,sizeof(long));
1206 weight = (long *)calloc(sites,sizeof(long));
1207
1208 for (i = 0; i < sites; i++)
1209 category[i] = 1,
1210 weight[i] = 1;
1211 extranum = 0;
1212 while (!(eoln(infile))) {
1213 ch = getc(infile);
1214 if (ch == '\n')
1215 ch = ' ';
1216 ch = isupper(ch) ? ch : toupper(ch);
1217 if (ch == 'C')
1218 extranum++;
1219 else if (ch != ' ') {
1220 printf("BAD OPTION CHARACTER: %c\n", ch);
1221 exit(-1);
1222 }
1223 }
1224 fscanf(infile, "%*[^\n]");
1225 getc(infile);
1226 for (i = 1; i <= extranum; i++) {
1227 ch = getc(infile);
1228 if (ch == '\n')
1229 ch = ' ';
1230 ch = isupper(ch) ? ch : toupper(ch);
1231 if (ch != 'W'){
1232 printf("ERROR: INCORRECT AUXILIARY OPTIONS LINE WHICH STARTS WITH %c\n",
1233 ch);
1234 exit(-1);
1235 }
1236 }
1237 if (categs <= 1)
1238 return;
1239 fprintf(outfile, "\nSite category Rate of change Probability\n");
1240 for (i = 1; i <= categs; i++)
1241 fprintf(outfile, "%12ld%13.3f%13.3f\n", i, rate[i - 1], probcat[i - 1]);
1242 putc('\n', outfile);
1243 } /* inputoptions */
1244
setuptree()1245 void setuptree()
1246 {
1247 long i, j;
1248 node *p, *q;
1249
1250 curtree = (tree *)calloc(1,sizeof(tree));
1251 curtree->nodep = (node **)calloc(numnodes + 3,sizeof(node *));
1252 alogf = (rlrec *)calloc(1,sizeof(rlrec));
1253 alogf->val = (double *)calloc(sites,sizeof(double));
1254 newtymenode(&curtree->tymelist);
1255 for (i = 0; i < 2; i++)
1256 freenodes[i] = NULL;
1257
1258 /* make the tips first */
1259 for (i = 0; i < numseq; i++) {
1260 curtree->nodep[i] = (node *)calloc(1,sizeof(node));
1261 curtree->nodep[i]->tip = true;
1262 curtree->nodep[i]->number = i + 1;
1263 VarMalloc(curtree->nodep[i],true);
1264 }
1265
1266 /* now make the interior nodes and the freenodes, but not the root */
1267 for (i = numseq; i < numnodes+2; i++) {
1268 q = NULL;
1269 for (j = 0; j < 3; j++) {
1270 p = (node *)calloc(1,sizeof(node));
1271 if (p == NULL) {
1272 fprintf(ERRFILE,"tree setup fails, allocate more space\n");
1273 exit(-1);
1274 }
1275 p->number = i + 1;
1276 p->tip = false;
1277 /* initialize the following pointers to NULL
1278 space will be allocated as appropiate in procedure
1279 orient() */
1280 p->x = NULL;
1281 /* end NULL assignments */
1282 p->next = q;
1283 q = p;
1284 }
1285 p->next->next->next = p; /* close up the chain into a loop */
1286 curtree->nodep[i] = p;
1287 if (i >= numnodes) {
1288 freenodes[i - numnodes] = p;
1289 /* do memory allocation for initial freenodes now, orient
1290 only covers nodes in initial tree */
1291 p->top = true;
1292 VarMalloc(p,true);
1293 }
1294 }
1295
1296 /* now make the root */
1297 curtree->nodep[rootnum - 1] = (node *)calloc(1,sizeof(node));
1298 curtree->nodep[rootnum - 1]->tip = true;
1299 curtree->nodep[rootnum - 1]->number = rootnum;
1300 VarMalloc(curtree->nodep[rootnum - 1],true);
1301 strncpy(curtree->nodep[rootnum-1]->nayme,"ROOT",4);
1302 /* guarantee that the root node contributes nothing to the likelihood
1303 of a tree (since its supposed to be at the end of a theoretically
1304 infinite root branch) */
1305 for (i = 0; i < sites; i++) {
1306 for (j = 0; j < categs; j++) {
1307 curtree->nodep[rootnum - 1]->x[i][j][baseA] = 1.0;
1308 curtree->nodep[rootnum - 1]->x[i][j][baseC] = 1.0;
1309 curtree->nodep[rootnum - 1]->x[i][j][baseG] = 1.0;
1310 curtree->nodep[rootnum - 1]->x[i][j][baseT] = 1.0;
1311 }
1312 }
1313 curtree->likelihood = NEGMAX;
1314 } /* setuptree */
1315
freetree()1316 void freetree()
1317 /* we do not free the following arrays:
1318 sum, theti, lntheti, fixed, numout, ne_ratio, mu_ratio, theta_ratio */
1319 {
1320 long i;
1321 node *p;
1322
1323 free(alogf->val);
1324 free(alogf);
1325 free(category);
1326 free(weight);
1327 freetymelist(curtree->tymelist);
1328 /* free the tips */
1329 for(i = 0; i < numseq; i++) {
1330 VarMalloc(curtree->nodep[i],false);
1331 free(curtree->nodep[i]);
1332 }
1333 /* free internal nodes including slidenodes */
1334 for(i = numseq; i < numnodes + 2; i++) {
1335 p = curtree->nodep[i];
1336 VarMalloc(p,false);
1337 VarMalloc(p->next,false);
1338 VarMalloc(p->next->next,false);
1339 free(p->next->next);
1340 free(p->next);
1341 free(p);
1342 }
1343 free(slidenodes);
1344 free(freenodes);
1345 /* free the root node */
1346 VarMalloc(curtree->nodep[rootnum-1],false);
1347 free(curtree->nodep[rootnum-1]);
1348 /* free the tree */
1349 free(curtree->nodep);
1350 /* free the aliases */
1351 free(siteptr);
1352 /* free the working arrays */
1353 free(weightrat);
1354 free(tbl);
1355 free(contribution[0]);
1356 free(contribution);
1357 } /* freetree */
1358
sitecompare(long site1,long site2)1359 boolean sitecompare(long site1, long site2)
1360 {
1361 long i;
1362
1363 for(i = 0; i < numseq; i++)
1364 if(dna->seqs[i][site1] != dna->seqs[i][site2])
1365 return false;
1366
1367 return true;
1368 } /* sitecompare */
1369
makesiteptr()1370 void makesiteptr()
1371 /* create the siteptr array: -1 means do a likelihood calculation
1372 for this site; a number >= 0 means use the site with that number
1373 for x (likelihood) values */
1374 {
1375 long whichsite, i;
1376 boolean found;
1377
1378 siteptr = (long *)calloc(sites,sizeof(long));
1379
1380 siteptr[0] = -1;
1381
1382 for(whichsite = 1; whichsite < sites; whichsite++) {
1383 found = false;
1384 for(i = whichsite - 1; i >= 0; i--) {
1385 if(sitecompare(i,whichsite)) {
1386 siteptr[whichsite] = i;
1387 found = true;
1388 break;
1389 }
1390 }
1391 if (found) continue;
1392 siteptr[whichsite] = -1;
1393 }
1394 } /* makesiteptr */
1395
getinput()1396 void getinput()
1397 {
1398 /* reads the input data */
1399 inputoptions();
1400 dna->freqa = freqa;
1401 dna->freqc = freqc;
1402 dna->freqg = freqg;
1403 dna->freqt = freqt;
1404 setuptree();
1405 getdata(curtree, dna, op, infile, outfile);
1406 makesiteptr();
1407 makevalues(dna, categs, curtree);
1408 if (op->freqsfrom) {
1409 empiricalfreqs(curtree, dna, weight);
1410 }
1411 getbasefreqs(dna, op, locus_ttratio, outfile);
1412 } /* getinput */
1413
watterson()1414 double watterson()
1415 {
1416 /* estimate theta using method of Watterson */
1417 long i, j, kn;
1418 boolean varies;
1419 double watter;
1420
1421 kn = 0;
1422 for (i = 0; i < sites; i++) {
1423 varies = false;
1424 for (j = 1; j < numseq; j++) {
1425 if (dna->seqs[j][i] != dna->seqs[0][i])
1426 varies = true;
1427 }
1428 if (varies)
1429 kn++;
1430 }
1431 watter = 0.0;
1432 if (kn > 0) {
1433 for (i = 1; i < numseq; i++)
1434 watter += 1.0 / i;
1435 watter = kn / (sites * watter);
1436 return watter;
1437 }
1438 fprintf(outfile, "Warning: There are no variable sites");
1439 fprintf(outfile, " in this data set.\n\n");
1440 if (menu) printf("Warning: There are no variable sites in this data set.\n");
1441 else {
1442 fprintf(simlog, "Warning: There are no variable sites");
1443 fprintf(simlog, " in this data set.\n\n");
1444 }
1445 exit(-1);
1446 } /* watterson */
1447
orient(node * p)1448 void orient(node *p)
1449 {
1450 tlist *t, *u;
1451
1452 t = curtree->tymelist;
1453
1454 if (p->tip) {
1455 p->top = true;
1456 p->tyme = 0.0;
1457 t->eventnode = p;
1458 t->branchlist[p->number - 1] = p;
1459 return;
1460 }
1461
1462 p->top = true;
1463 curtree->nodep[p->number-1] = p; /* insure that curtree->nodep points
1464 to nodes with info */
1465
1466 /* since p is a top nodelet, it needs to actually store
1467 likelihood information, x is a NULL pointer
1468 in all other non-tip nodelets */
1469 VarMalloc(p,true);
1470
1471 p->next->top = false;
1472 p->next->next->top = false;
1473
1474 orient(p->next->back);
1475 orient(p->next->next->back);
1476 p->tyme = p->next->length + p->next->back->tyme;
1477 p->next->tyme = p->tyme;
1478 p->next->next->tyme = p->tyme;
1479 if (p->number == curtree->root->back->number) {
1480 p->back->top = false;
1481 p->back->tyme = rootlength;
1482 }
1483 newtymenode(&u);
1484 u->eventnode = p;
1485 while (t != NULL) {
1486 if (u->eventnode->tyme < t->eventnode->tyme) {
1487 u->prev = t->prev;
1488 t->prev = u;
1489 u->succ = t;
1490 u->prev->succ = u;
1491 break;
1492 }
1493 if (t->succ != NULL)
1494 t = t->succ;
1495 else {
1496 t->succ = u;
1497 u->prev = t;
1498 u->succ = NULL;
1499 break;
1500 }
1501 }
1502 } /* orient */
1503
finishsetup(node * p)1504 void finishsetup(node *p)
1505 {
1506 if (p->tip) {
1507 ltov(p);
1508 return;
1509 }
1510 ltov(p);
1511 finishsetup(p->next->back);
1512 finishsetup(p->next->next->back);
1513 return;
1514 } /* finishsetup */
1515
initbranchlist()1516 void initbranchlist()
1517 {
1518 tlist *t;
1519 node *p, *q;
1520 long i, j, k, n;
1521
1522 t = curtree->tymelist;
1523 n = numseq;
1524 t->numbranch = n;
1525 t->age = t->succ->eventnode->tyme;
1526 t = t->succ;
1527 for (i = 0; i < (numnodes - numseq); i++) {
1528 /* for each interior node, do...assumes at least 3 tips */
1529 n--;
1530 t->numbranch = n;
1531 if (n == 1)
1532 t->age = t->eventnode->tyme + rootlength;
1533 else
1534 t->age = t->succ->eventnode->tyme;
1535 p = t->eventnode->next->back;
1536 q = t->eventnode->next->next->back;
1537 k = 0;
1538 for (j = 0; j < t->prev->numbranch ; j++) {
1539 /* for the number of branches above the coalescent node, do...*/
1540 if (t->prev->branchlist[j] != p && t->prev->branchlist[j] != q) {
1541 t->branchlist[k] = t->prev->branchlist[j];
1542 k++;
1543 }
1544 }
1545 t->branchlist[t->numbranch - 1] = t->eventnode;
1546 t = t->succ;
1547 }
1548 /* initialize the slidelist, assume that curtree->nodep[numseq] through
1549 curtree->nodep[numnodes-1] point to all the interior nodes of the
1550 initial tree (one node will be curtree->root->back and so ineligible
1551 to be slid) */
1552 i = 0;
1553 slidenum = 0;
1554 for (j = numseq; j < numnodes; j++) {
1555 if (!(curtree->nodep[j]->back->number == rootnum)) {
1556 slidenodes[i] = curtree->nodep[j];
1557 i++;
1558 slidenum++;
1559 }
1560 }
1561 } /* initbranchlist */
1562
inittable()1563 void inittable()
1564 {
1565 long i;
1566 tbl = (valrec *)calloc(categs,sizeof(valrec));
1567 /* Define a lookup table. Precompute values and store them in a table */
1568 for (i = 0; i < categs; i++) {
1569 tbl[i].rat_xi = rate[i] * dna->xi;
1570 tbl[i].rat_xv = rate[i] * dna->xv;
1571 }
1572 } /* inittable */
1573
initweightrat()1574 void initweightrat()
1575 {
1576 long i;
1577 weightrat = (double *)calloc(sites,sizeof(double));
1578 sumweightrat = 0.0;
1579 for (i = 0; i < sites; i++) {
1580 weightrat[i] = weight[i] * rate[category[i] - 1];
1581 sumweightrat += weightrat[i];
1582 }
1583 } /* initweightrat */
1584
treeout(node * p,long s,FILE ** usefile)1585 void treeout(node *p, long s, FILE **usefile)
1586 {
1587 /* write out file with representation of final tree */
1588 long i, n, w;
1589 char c;
1590 double x;
1591
1592 if (p->tip) {
1593 n = 0;
1594 for (i = 1; i <= NMLNGTH; i++) {
1595 if (p->nayme[i - 1] != ' ')
1596 n = i;
1597 }
1598 for (i = 0; i < n; i++) {
1599 c = p->nayme[i];
1600 if (c == ' ')
1601 c = '_';
1602 putc(c, *usefile);
1603 }
1604 col += n;
1605 } else {
1606 putc('(', *usefile);
1607 col++;
1608 treeout(p->next->back, s, usefile);
1609 putc(',', *usefile);
1610 col++;
1611 if (col > 45) {
1612 putc('\n', *usefile);
1613 col = 0;
1614 }
1615 treeout(p->next->next->back, s, usefile);
1616 putc(')', *usefile);
1617 col++;
1618 }
1619 if (p->v >= 1.0)
1620 x = -1.0;
1621 else
1622 x = lengthof(p);
1623 if (x > 0.0)
1624 w = (long)(0.4343 * log(x));
1625 else if (x == 0.0)
1626 w = 0;
1627 else
1628 w = (long)(0.4343 * log(-x)) + 1;
1629 if (w < 0)
1630 w = 0;
1631 if (p == curtree->root->back)
1632 putc(';', *usefile);
1633 else {
1634 fprintf(*usefile, ":%*.10f", (int)(w + 7), x);
1635 col += w + 8;
1636 }
1637 } /* treeout */
1638
evaluate(tree * tr,boolean first)1639 void evaluate(tree *tr, boolean first)
1640 {
1641 double temp, sum2, sumc, sumterm, lterm, termcheck;
1642 contribarr like, nulike, term, clai;
1643 long i, j, k;
1644 node *p;
1645 sitelike x1;
1646
1647 like = (double *)calloc(categs,sizeof(double));
1648 nulike = (double *)calloc(categs,sizeof(double));
1649 term = (double *)calloc(categs,sizeof(double));
1650 clai = (double *)calloc(categs,sizeof(double));
1651
1652 temp = 0.0;
1653 p = tr->root->back;
1654
1655 for (i = 0; i < sites; i++) {
1656 termcheck = 0.0;
1657 for (j = 0; j < categs; j++) {
1658 memcpy((void *)x1, (void *)p->x[i][j], sizeof(sitelike));
1659 term[j] = dna->freqa * x1[baseA] + dna->freqc * x1[baseC] +
1660 dna->freqg * x1[baseG] + dna->freqt * x1[baseT];
1661 termcheck += term[j];
1662 }
1663 if (!termcheck) {
1664 fprintf(ERRFILE,"Encountered tree incompatible with data\n");
1665 if(first) {
1666 fprintf(ERRFILE,"starting tree needs to be legal\n");
1667 exit(-1);
1668 }
1669 curtree->likelihood = NEGMAX;
1670 return;
1671 }
1672 sumterm = 0.0;
1673 for (j = 0; j < categs; j++)
1674 sumterm += probcat[j] * term[j];
1675 lterm = log(sumterm);
1676 for (j = 0; j < categs; j++)
1677 clai[j] = term[j] / sumterm;
1678 memcpy((void *)contribution[i], (void *)clai, categs*sizeof(double));
1679 if (!op->autocorr)
1680 alogf->val[i] = lterm;
1681 temp += weight[i] * lterm;
1682 }
1683 for (j = 0; j < categs; j++)
1684 like[j] = 1.0;
1685 for (i = 0; i < sites; i++) {
1686 sumc = 0.0;
1687 for (k = 1; k <= categs; k++)
1688 sumc += probcat[k - 1] * like[k - 1];
1689 sumc *= lambda;
1690 memcpy((void *)clai, (void *)contribution[i], categs*sizeof(double));
1691 for (j = 0; j < categs; j++)
1692 nulike[j] = ((1.0 - lambda) * like[j] + sumc) * clai[j];
1693 memcpy((void *)like, (void *)nulike, categs*sizeof(double));
1694 }
1695 sum2 = 0.0;
1696 for (i = 0; i < categs; i++)
1697 sum2 += probcat[i] * like[i];
1698 temp += log(sum2);
1699 curtree->likelihood = temp;
1700 free(like);
1701 free(nulike);
1702 free(term);
1703 free(clai);
1704 } /* evaluate */
1705
nuview(node * p)1706 void nuview(node *p)
1707 {
1708 long i, j;
1709 double w1, w2, lw1, lw2, yy1, yy2, ww1zz1, vv1zz1, ww2zz2, vv2zz2,
1710 vv1zz1_sumr1, vv2zz2_sumr2, vv1zz1_sumy1, vv2zz2_sumy2, sum1, sum2,
1711 sumr1, sumr2, sumy1, sumy2;
1712 node *q, *r;
1713 sitelike xx1, xx2, xx3;
1714
1715 q = p->next->back;
1716 r = p->next->next->back;
1717
1718 w1 = 1.0 - q->v;
1719 w2 = 1.0 - r->v;
1720 if (w1 > 0.0) {
1721 lw1 = log(w1);
1722 for (i = 0; i < categs; i++) {
1723 tbl[i].ww1 = exp(tbl[i].rat_xi * lw1);
1724 tbl[i].zz1 = exp(tbl[i].rat_xv * lw1);
1725 tbl[i].ww1zz1 = tbl[i].ww1 * tbl[i].zz1;
1726 tbl[i].vv1zz1 = (1.0 - tbl[i].ww1) * tbl[i].zz1;
1727 }
1728 }
1729 if (w2 > 0.0) {
1730 lw2 = log(w2);
1731 for (i = 0; i < categs; i++) {
1732 tbl[i].ww2 = exp(tbl[i].rat_xi * lw2);
1733 tbl[i].zz2 = exp(tbl[i].rat_xv * lw2);
1734 tbl[i].ww2zz2 = tbl[i].ww2 * tbl[i].zz2;
1735 tbl[i].vv2zz2 = (1.0 - tbl[i].ww2) * tbl[i].zz2;
1736 }
1737 }
1738 for (i = 0; i < sites; i++) {
1739 for (j = 0; j < categs; j++) {
1740 if(siteptr[i] == -1) { /* if we need to calculate this site */
1741 if (w1 <= 0.0) {
1742 ww1zz1 = 0.0;
1743 vv1zz1 = 0.0;
1744 yy1 = 1.0;
1745 } else {
1746 ww1zz1 = tbl[j].ww1zz1;
1747 vv1zz1 = tbl[j].vv1zz1;
1748 yy1 = 1.0 - tbl[j].zz1;
1749 }
1750 if (w2 <= 0.0) {
1751 ww2zz2 = 0.0;
1752 vv2zz2 = 0.0;
1753 yy2 = 1.0;
1754 } else {
1755 ww2zz2 = tbl[j].ww2zz2;
1756 vv2zz2 = tbl[j].vv2zz2;
1757 yy2 = 1.0 - tbl[j].zz2;
1758 }
1759 memcpy((void *)xx1, (void *)q->x[i][j], sizeof(sitelike));
1760 memcpy((void *)xx2, (void *)r->x[i][j], sizeof(sitelike));
1761 sum1 = yy1 * (dna->freqa * xx1[baseA] + dna->freqc * xx1[baseC] +
1762 dna->freqg * xx1[baseG] + dna->freqt * xx1[baseT]);
1763 sum2 = yy2 * (dna->freqa * xx2[baseA] + dna->freqc * xx2[baseC] +
1764 dna->freqg * xx2[baseG] + dna->freqt * xx2[baseT]);
1765 sumr1 = dna->freqar * xx1[baseA] + dna->freqgr * xx1[baseG];
1766 sumr2 = dna->freqar * xx2[baseA] + dna->freqgr * xx2[baseG];
1767 sumy1 = dna->freqcy * xx1[baseC] + dna->freqty * xx1[baseT];
1768 sumy2 = dna->freqcy * xx2[baseC] + dna->freqty * xx2[baseT];
1769 vv1zz1_sumr1 = vv1zz1 * sumr1;
1770 vv2zz2_sumr2 = vv2zz2 * sumr2;
1771 vv1zz1_sumy1 = vv1zz1 * sumy1;
1772 vv2zz2_sumy2 = vv2zz2 * sumy2;
1773 xx3[baseA] = (sum1 + ww1zz1 * xx1[baseA] + vv1zz1_sumr1) *
1774 (sum2 + ww2zz2 * xx2[baseA] + vv2zz2_sumr2);
1775 xx3[baseC] = (sum1 + ww1zz1 * xx1[baseC] + vv1zz1_sumy1) *
1776 (sum2 + ww2zz2 * xx2[baseC] + vv2zz2_sumy2);
1777 xx3[baseG] = (sum1 + ww1zz1 * xx1[baseG] + vv1zz1_sumr1) *
1778 (sum2 + ww2zz2 * xx2[baseG] + vv2zz2_sumr2);
1779 xx3[baseT] = (sum1 + ww1zz1 * xx1[baseT] + vv1zz1_sumy1) *
1780 (sum2 + ww2zz2 * xx2[baseT] + vv2zz2_sumy2);
1781 memcpy((void *)p->x[i][j], (void *)xx3, sizeof(sitelike));
1782 }
1783 else {
1784 /* this site is just like site #(siteptr[i]), use its values */
1785 memcpy((void *)p->x[i][j], (void *)p->x[siteptr[i]][j], sizeof(sitelike));
1786 }
1787 }
1788 }
1789 } /* nuview */
1790
update(node * p)1791 void update(node *p)
1792 {
1793 if (!p->tip)
1794 nuview(p);
1795 } /* update */
1796
smooth(node * p)1797 void smooth(node *p)
1798 {
1799 if (!p->tip) {
1800 if (!p->next->top)
1801 smooth(p->next->back);
1802 if (!p->next->next->top)
1803 smooth(p->next->next->back);
1804 }
1805 update(p);
1806 } /* smooth */
1807
localsmooth(node * p)1808 void localsmooth(node *p)
1809 {
1810 if (p->number != curtree->root->number) {
1811 p = findtop(p);
1812 nuview(p);
1813 }
1814 if (p->number != curtree->root->number)
1815 localsmooth(p->back);
1816 } /* localsmooth */
1817
testratio()1818 boolean testratio()
1819 {
1820 /* decide to accept or not */
1821 double test, x;
1822
1823 if(curtree->likelihood == NEGMAX)
1824 return false;
1825 test = curtree->likelihood - oldlikelihood;
1826 if (test >= 1.0)
1827 return true;
1828 else {
1829 x = log(randum());
1830 if (test >= x)
1831 return true;
1832 else
1833 return false;
1834 }
1835 } /* testratio */
1836
seekch(char c)1837 void seekch(char c) /* use only in reading file intree! */
1838 {
1839 if (gch == c)
1840 return;
1841 do {
1842 if (eoln(intree)) {
1843 fscanf(intree, "%*[^\n]");
1844 getc(intree);
1845 }
1846 gch = getc(intree);
1847 if (gch == '\n')
1848 gch = ' ';
1849 } while (gch != c);
1850 } /* seekch */
1851
getch(char * c)1852 void getch(char *c) /* use only in reading file intree! */
1853 {
1854 /* get next nonblank character */
1855 do {
1856 if (eoln(intree)) {
1857 fscanf(intree, "%*[^\n]");
1858 getc(intree);
1859 }
1860 *c = getc(intree);
1861 if (*c == '\n')
1862 *c = ' ';
1863 } while (*c == ' ');
1864 } /* getch */
1865
processlength(node * p)1866 void processlength(node *p)
1867 {
1868 long digit;
1869 double valyew, divisor;
1870 boolean pointread;
1871
1872 pointread = false;
1873 valyew = 0.0;
1874 divisor = 1.0;
1875 getch(&gch);
1876 digit = gch - '0';
1877 while (((unsigned long)digit <= 9) || gch == '.'){
1878 if (gch == '.')
1879 pointread = true;
1880 else {
1881 valyew = valyew * 10.0 + digit;
1882 if (pointread)
1883 divisor *= 10.0;
1884 }
1885 getch(&gch);
1886 digit = gch - '0';
1887 }
1888 p->length = valyew / divisor;
1889 p->back->length = p->length;
1890 } /* processlength */
1891
addelement(node * p,long * nextnode)1892 void addelement(node *p, long *nextnode)
1893 {
1894 node *q;
1895 long i, n;
1896 boolean found;
1897 char str[NMLNGTH];
1898
1899 getch(&gch);
1900 if (gch == '(') {
1901 (*nextnode)++;
1902 q = curtree->nodep[(*nextnode) - 1];
1903 hookup(p, q);
1904 addelement(q->next,nextnode);
1905 seekch(',');
1906 addelement(q->next->next, nextnode);
1907 seekch(')');
1908 getch(&gch);
1909 } else {
1910 for (i = 0; i < NMLNGTH; i++)
1911 str[i] = ' ';
1912 n = 1;
1913 do {
1914 if (gch == '_')
1915 gch = ' ';
1916 str[n - 1] = gch;
1917 if (eoln(intree)) {
1918 fscanf(intree, "%*[^\n]");
1919 getc(intree);
1920 }
1921 gch = getc(intree);
1922 if (gch == '\n')
1923 gch = ' ';
1924 n++;
1925 } while (gch != ':' && gch != ',' && gch != ')' && n <= NMLNGTH);
1926 n = 1;
1927 do {
1928 found = true;
1929 for (i = 0; i < NMLNGTH; i++)
1930 found = (found && str[i] == curtree->nodep[n - 1]->nayme[i]);
1931 if (!found)
1932 n++;
1933 } while (!(n > numseq || found));
1934 if (n > numseq) {
1935 printf("Cannot find sequence: ");
1936 for (i = 0; i < NMLNGTH; i++)
1937 putchar(str[i]);
1938 putchar('\n');
1939 }
1940 hookup(curtree->nodep[n - 1], p);
1941 }
1942 if (gch == ':')
1943 processlength(p);
1944 } /* addelement */
1945
treeread()1946 void treeread()
1947 {
1948 long nextnode;
1949 node *p;
1950
1951 curtree->root = curtree->nodep[rootnum - 1];
1952 getch(&gch);
1953 if (gch == '(') {
1954 nextnode = numseq + 1;
1955 p = curtree->nodep[nextnode - 1];
1956 addelement(p, &nextnode);
1957 seekch(',');
1958 addelement(p->next, &nextnode);
1959 hookup(p->next->next, curtree->nodep[rootnum - 1]);
1960 p->next->next->length = rootlength;
1961 curtree->nodep[rootnum - 1]->length = p->next->next->length;
1962 ltov(curtree->nodep[rootnum - 1]);
1963 }
1964 fscanf(intree, "%*[^\n]");
1965 getc(intree);
1966 } /* treeread */
1967
treevaluate()1968 void treevaluate()
1969 {
1970
1971 smooth(curtree->root->back);
1972 smooth(curtree->root);
1973 evaluate(curtree,true);
1974 } /* treevaluate */
1975
localevaluate(node * p,node * pansdaught)1976 void localevaluate(node *p, node *pansdaught)
1977 /* routine assumes that p points to the only 'top' nodelet
1978 in node 'p' */
1979 {
1980
1981 /* first update all daughters and p itself */
1982 if (!p->next->back->tip)
1983 nuview(p->next->back);
1984 if (!p->next->next->back->tip)
1985 nuview(p->next->next->back);
1986 nuview(p);
1987 if (!pansdaught->tip)
1988 nuview(pansdaught);
1989 /* now update the rest of the tree */
1990 localsmooth(p->back);
1991 evaluate(curtree,false);
1992 } /* localevaluate */
1993
copynode(node * source,node * target)1994 void copynode(node *source, node *target)
1995 /* copies source node to target node */
1996 {
1997 long i, j;
1998
1999 for (i = 1; i <= 3; i++) {
2000 /* NEVER! target->next := source->next; */
2001 target->back = source->back;
2002 target->tip = source->tip;
2003 /* but NOT target->number := source->number; */
2004 if (source->x != NULL) {
2005 VarMalloc(target,true);
2006 for (j = 0; j < sites; j++) {
2007 memcpy((void *)target->x[j], (void *)source->x[j], categs*sizeof(sitelike));
2008 }
2009 }
2010 else
2011 VarMalloc(target,false);
2012 memcpy((void *)target->nayme, (void *)source->nayme, sizeof(naym));
2013 target->top = source->top;
2014 target->v = source->v;
2015 target->tyme = source->tyme;
2016 target->length = source->length;
2017 source = source->next;
2018 target = target->next;
2019 }
2020 } /* copynode */
2021
2022 /* joinnode and constructree are used for constructing a rather bad
2023 starting tree if the user doesn't provide one */
joinnode(double length,node * p,node * q)2024 void joinnode(double length, node *p, node *q)
2025 {
2026 hookup(p,q);
2027 p->length = length;
2028 q->length = length;
2029 ltov(p);
2030 } /* joinnode */
2031
constructtree(long numtips,double branchlength)2032 void constructtree(long numtips, double branchlength)
2033 {
2034 long i, j, nextnode;
2035 double height;
2036 node *p, *q;
2037
2038 curtree->root = curtree->nodep[rootnum - 1];
2039 nextnode = numseq;
2040 p = curtree->root;
2041 q = curtree->nodep[nextnode];
2042
2043 p->back = q;
2044 q->back = p;
2045 p->length = rootlength;
2046 q->length = rootlength;
2047 ltov(p);
2048
2049 height = (numtips - 1) * branchlength;
2050 p->tyme = rootlength + height;
2051 for (i = 0; i < numtips - 1; i++) {
2052 p = curtree->nodep[i];
2053 q = curtree->nodep[nextnode]->next;
2054 joinnode(height,p,q);
2055 q = q->next;
2056 if (i != numtips-2) {
2057 nextnode++;
2058 p = curtree->nodep[nextnode];
2059 joinnode(branchlength,p,q);
2060 height -= branchlength;
2061 }
2062 else {
2063 p = curtree->nodep[numtips - 1];
2064 joinnode(height,p,q);
2065 }
2066 for (j = 0; j < 3; j++)
2067 q->tyme = height;
2068 }
2069 } /* constructtree */
2070 /* End bad starting tree construction */
2071
updateslide(node * p,boolean wanted)2072 void updateslide(node *p, boolean wanted)
2073 /* pass me FALSE only if sure that the node is invalid */
2074 {
2075 boolean valid, found;
2076 node *q;
2077 long j, k;
2078
2079 k = 0; /* just to be careful */
2080
2081 valid = true;
2082 q = p;
2083 if (!wanted)
2084 valid = false;
2085 else {
2086 if (q->tip)
2087 valid = false;
2088 else {
2089 q =findtop(q);
2090 if (q->back->tip)
2091 valid = false;
2092 }
2093 }
2094 found = false;
2095 j = 1;
2096 while (!found && j <= slidenum) {
2097 if (slidenodes[j - 1]->number == p->number) {
2098 found = true;
2099 k = j;
2100 }
2101 j++;
2102 }
2103 if (valid && !found) {
2104 slidenum++;
2105 slidenodes[slidenum - 1] = p;
2106 }
2107 if (valid || !found)
2108 return;
2109 while (k < slidenum) {
2110 slidenodes[k - 1] = slidenodes[k];
2111 k++;
2112 }
2113 slidenum--;
2114 } /* updateslide */
2115
rebuildbranch()2116 void rebuildbranch()
2117 {
2118 tlist *t;
2119 node *p;
2120 long i, k;
2121 boolean done;
2122
2123 t = curtree->tymelist->succ;
2124 done = false;
2125 do {
2126 if (t->succ == NULL) done = true;
2127 p = t->eventnode;
2128 k = 1;
2129 p = findtop(p);
2130 for (i = 0; i < t->prev->numbranch; i++) {
2131 if (t->prev->branchlist[i] != p->next->back &&
2132 t->prev->branchlist[i] != p->next->next->back) {
2133 t->branchlist[k - 1] = t->prev->branchlist[i];
2134 k++;
2135 }
2136 }
2137 t->numbranch = t->prev->numbranch - 1;
2138 t->branchlist[t->numbranch - 1] = p;
2139 t = t->succ;
2140 } while (!done);
2141 } /* rebuildbranch */
2142
2143 /****************************************************************
2144 * setlength() returns true if it succeeds in setting a length; *
2145 * false otherwise */
setlength(long numl,long numother,double tstart,double tlength,node * p)2146 boolean setlength(long numl, long numother, double tstart,
2147 double tlength, node *p)
2148 {
2149 double x, e1, r;
2150
2151 r = randum();
2152 e1 = (numl - 1.0) * 2 + numother * 2;
2153 x = -(theta0 / e1) *
2154 log(1 - r * (1 - exp(-(e1 * tlength / theta0))));
2155 if ((unsigned)x > tlength) {
2156 fprintf(ERRFILE,"disaster in setlength\n");
2157 return(false);
2158 }
2159 if (!growth0 || !op->growthused) p->tyme = tstart + x;
2160 else {
2161 tstart = to_magic(tstart);
2162 x = tstart + x;
2163 /* if the time is now impossible to big... */
2164 if (x*growth0 < -1.0) return(false);
2165 else p->tyme = to_real(x);
2166 }
2167
2168 return(true);
2169 } /* setlength */
2170
2171 #define TRY_SOLVE 20 /* number of times to iteratively solve for
2172 length of new branch */
2173 /*********************************************************************
2174 * setlength2() returns true if it succeeds in setting both lengths; *
2175 * false otherwise */
setlength2(long numother,double tstart,double tlength,node * p,node * q)2176 boolean setlength2(long numother, double tstart, double tlength,
2177 node *p, node *q)
2178 {
2179 long i;
2180 double x, xmin, xmax, r, xnew, e1, e2, norm;
2181
2182 r = randum();
2183 e1 = exp(numother * -2 * tlength / theta0);
2184 e2 = exp(-((numother * 2 + 2) * tlength / theta0));
2185 norm = -3 * e1 / ((numother + 1) * twocollis(numother, tlength));
2186 xmin = 0.0;
2187 xmax = tlength;
2188 for (i = 0; i < TRY_SOLVE; i++) {
2189 x = (xmax + xmin) / 2.0;
2190 xnew = norm *
2191 (1.0 / (numother * 2 + 3) *
2192 (exp(-((numother * 4 + 6) * x / theta0)) - 1) -
2193 e2 / (numother + 2) * (exp(-((numother * 2 + 4) * x / theta0)) - 1));
2194 if (xnew > r)
2195 xmax = x;
2196 else
2197 xmin = x;
2198 }
2199 if ((unsigned)x > tlength) {
2200 fprintf(ERRFILE,"disaster in setlength2\n");
2201 return(false);
2202 }
2203 if (!growth0 || !op->growthused) p->tyme = tstart + x;
2204 else {
2205 tstart = to_magic(tstart) + x;
2206 /* if the time is now impossible to big... */
2207 if (growth0*tstart < -1.0) return(false);
2208 else
2209 p->tyme = to_real(tstart);
2210 }
2211
2212
2213 return(setlength(2L, numother, p->tyme, tlength - x, q));
2214 } /* setlength2 */
2215
updatebranch(node * oldans,node * oldp,node * prime)2216 void updatebranch(node *oldans, node *oldp, node *prime)
2217 {
2218 subtymelist(oldans, oldp);
2219 inserttymelist(prime);
2220 rebuildbranch();
2221 } /* updatebranch */
2222
counttymelist(tlist * first,tlist * last)2223 long counttymelist(tlist *first, tlist *last)
2224 {
2225 long count;
2226 tlist *t;
2227
2228 count = 0;
2229 for (t = first; t != last; t = t->succ)
2230 count++;
2231 return(count+1);
2232 } /* counttymelist */
2233
accept_tree(node * prime,node * primans,node * ans,node * pansdaught,node * newbr[],tlist * tplus)2234 boolean accept_tree(node *prime, node *primans, node *ans,
2235 node *pansdaught, node *newbr[], tlist *tplus)
2236 {
2237 node *p;
2238 long i, leftout;
2239
2240 /* set up the phylogeny */
2241 if (prime->tyme < tplus->eventnode->tyme) {
2242 /* case 1 */
2243 hookup(newbr[0], prime->next);
2244 hookup(newbr[1], prime->next->next);
2245 hookup(newbr[2], pansdaught);
2246 } else {
2247 /* case 2 */
2248 leftout = (long)(randum() * 3.0) + 1;
2249 p = prime->next;
2250 for (i = 1; i <= 3; i++) {
2251 if (i != leftout) {
2252 hookup(newbr[i - 1], p);
2253 p = p->next;
2254 } else hookup(newbr[i - 1], pansdaught);
2255 }
2256 }
2257 if (primans->next == pansdaught)
2258 hookup(primans->next->next, ans);
2259 else hookup(primans->next, ans);
2260
2261 prime->next->tyme = prime->tyme;
2262 prime->next->next->tyme = prime->tyme;
2263 primans->next->tyme = primans->tyme;
2264 primans->next->next->tyme = primans->tyme;
2265
2266 for (i = 0; i <= 2; i++) ltov(newbr[i]);
2267 ltov(prime);
2268 ltov(ans);
2269 /* get acceptance ratio */
2270 if (newbr[1]->tyme == ans->tyme) return(true);
2271
2272 localevaluate(prime,pansdaught->back);
2273 return (testratio());
2274
2275 } /* accept_tree */
2276
slide()2277 boolean slide()
2278 {
2279 /* Local variables for slide: */
2280
2281 node *prime, *oldp, *oldans, *primans, *pansdaught, *ans, *p,
2282 *oldbr[3], *newbr[3];
2283 long i, j, k, cline, numother, numintervals;
2284 double chance, tlength, normalizer, **coll2, **coll3, c[3];
2285 tlist *t, *tstart, *tend, *tplus;
2286 boolean accept_slide, done, skipped, succeeded;
2287
2288 newbr[0] = newbr[1] = newbr[2] = NULL; /* just to be careful */
2289
2290 do { /* There exists a chance that randum returns a 1 */
2291 i = (long)(randum() * slidenum) + 1;
2292 } while (i == slidenum + 1);
2293 oldp = slidenodes[i - 1];
2294 oldp = findtop(oldp);
2295 oldans = oldp->back;
2296 /* copy old nodes to new */
2297 newnode(&prime);
2298 copynode(oldp, prime);
2299 newnode(&primans);
2300 copynode(oldans, primans);
2301 hookup(prime, primans);
2302 /* name and connect nodes */
2303 oldbr[0] = prime->next->back;
2304 oldbr[1] = prime->next->next->back;
2305 if (!primans->next->top) {
2306 pansdaught = primans->next;
2307 oldbr[2] = primans->next->back;
2308 ans = primans->next->next->back;
2309 } else {
2310 pansdaught = primans->next->next;
2311 oldbr[2] = primans->next->next->back;
2312 ans = primans->next->back;
2313 }
2314 /* tymelist sort the three branches' tips in arbitrary order */
2315 j = 1;
2316 for (i = 0; i <= 2; i++) {
2317 if (oldbr[i]->tip) {
2318 newbr[j - 1] = oldbr[i];
2319 j++;
2320 }
2321 }
2322 /* remainder of tree */
2323 if (j <= 3) {
2324 t = curtree->tymelist->succ;
2325 done = false;
2326 while (!done) {
2327 for (i = 0; i <= 2; i++) {
2328 if (oldbr[i]->number == t->eventnode->number) {
2329 newbr[j - 1] = oldbr[i];
2330 j++;
2331 if (j > 3) done = true;
2332 }
2333 }
2334 t = t->succ;
2335 if (t == NULL) {
2336 printf("ERROR IN TYMESORT\n");
2337 exit(-1);
2338 }
2339 }
2340 }
2341 tplus = gettymenode(newbr[2]->number);
2342 skipped = false; /* needed for frees in zero length case */
2343 succeeded = true; /* records success of setlength operations */
2344 /* zero length branches are a special case */
2345 if (newbr[1]->tyme == ans->tyme) {
2346 prime->tyme = newbr[1]->tyme;
2347 primans->tyme = newbr[1]->tyme;
2348 coll2 = NULL; /* just to be careful */
2349 coll3 = NULL;
2350 skipped = true;
2351 } else {
2352 /* initialize probability arrays for state (ie. 1,2 or 3 branches
2353 present) */
2354 tstart = gettymenode(newbr[1]->number);
2355 tend = gettymenode(ans->number);
2356 numintervals = counttymelist(tstart,tend);
2357 coll2 = (double **)calloc(2,sizeof(double *));
2358 coll2[0] = (double *)calloc(numintervals,sizeof(double));
2359 coll2[1] = (double *)calloc(numintervals,sizeof(double));
2360 coll3 = (double **)calloc(3,sizeof(double *));
2361 coll3[0] = (double *)calloc(numintervals,sizeof(double));
2362 coll3[1] = (double *)calloc(numintervals,sizeof(double));
2363 coll3[2] = (double *)calloc(numintervals,sizeof(double));
2364 t = tstart;
2365 i = 0;
2366 cline = 1;
2367 numother = tstart->numbranch - 2;
2368 /* initialize 2-array */
2369 coll2[0][i] = 0.0;
2370 coll2[1][i] = 1.0;
2371 /* fill up 2-array */
2372 while (t != tplus) {
2373 i++;
2374 tlength = howlong(t);
2375 if ((numother * 6 + 6) * tlength / theta0 > 20.0)
2376 tlength = 20 * theta0 / (numother * 6 + 6);
2377 coll2[0][i] = coll2[0][i - 1] * zerocollis(1L, numother, tlength)
2378 + coll2[1][i - 1] * onecollis(2L, numother, tlength);
2379 coll2[1][i] = coll2[1][i - 1] * zerocollis(2L, numother, tlength);
2380 normalizer = coll2[0][i] + coll2[1][i];
2381 coll2[0][i] /= normalizer;
2382 coll2[1][i] /= normalizer;
2383 if (normalizer == 0.0) {
2384 fprintf(ERRFILE,"Encountered machine precision limits!\n");
2385 exit(-1);
2386 }
2387 t = t->succ;
2388 if (t->eventnode->number != oldp->number &&
2389 t->eventnode->number != oldans->number)
2390 numother--;
2391 }
2392 if (newbr[2]->tyme >= ans->tyme) {
2393 /* case 2 zero length */
2394 primans->tyme = tplus->eventnode->tyme;
2395 } else {
2396 /* initialize 3-array */
2397 j = 0;
2398 numother--;
2399 coll3[0][j] = 0.0;
2400 coll3[1][j] = coll2[0][i];
2401 coll3[2][j] = coll2[1][i];
2402 /* fill up 3-array */
2403 done = false;
2404 while (!done) {
2405 j++;
2406 tlength = howlong(t);
2407 if ((numother * 6 + 6) * tlength / theta0 > 20.0)
2408 tlength = 20 * theta0 / (numother * 6 + 6);
2409 coll3[0][j] = coll3[0][j - 1] * zerocollis(1L, numother, tlength) +
2410 coll3[1][j - 1] * onecollis(2L, numother, tlength) +
2411 coll3[2][j - 1] * twocollis(numother, tlength);
2412 coll3[1][j] = coll3[1][j - 1] * zerocollis(2L, numother, tlength) +
2413 coll3[2][j - 1] * onecollis(3L, numother, tlength);
2414 coll3[2][j] = coll3[2][j - 1] * zerocollis(3L, numother, tlength);
2415 normalizer = coll3[0][j] + coll3[1][j] + coll3[2][j];
2416 coll3[0][j] /= normalizer;
2417 coll3[1][j] /= normalizer;
2418 coll3[2][j] /= normalizer;
2419 if (normalizer == 0.0) {
2420 fprintf(ERRFILE,"Encountered machine precision limits!\n");
2421 exit(-1);
2422 }
2423 if (t->succ == tend) {
2424 done = true;
2425 break;
2426 }
2427 t = t->succ;
2428 if (t->eventnode->number != oldp->number &&
2429 t->eventnode->number != oldans->number)
2430 numother--;
2431 }
2432 /* now find out when prime and primans collide */
2433 k = j;
2434 while (cline != 3 && k != 0 && t != NULL) {
2435 tlength = howlong(t);
2436 if ((numother * 6 + 6) * tlength / theta0 > 20.0)
2437 tlength = 20 * theta0 / (numother * 6 + 6);
2438 chance = randum();
2439 if (cline == 1) {
2440 c[0] = coll3[0][k-1] * zerocollis(1L, numother, tlength);
2441 c[1] = coll3[1][k-1] * onecollis(2L, numother, tlength);
2442 c[2] = coll3[2][k-1] * twocollis(numother, tlength);
2443 normalizer = c[0] + c[1] + c[2];
2444 c[0] /= normalizer;
2445 c[1] /= normalizer;
2446 if (chance > c[0]) {
2447 chance -= c[0];
2448 if (chance > c[1]) { /* two collisions */
2449 cline += 2;
2450 if (!setlength2(numother, t->eventnode->tyme, tlength, prime,
2451 primans)) succeeded = false;
2452 } else { /* one collision */
2453 cline++;
2454 if (!setlength(2L, numother, t->eventnode->tyme,
2455 tlength, primans)) succeeded = false;
2456 }
2457 }
2458 } else { /* cline must equal 2 */
2459 c[0] = coll3[1][k-1] * zerocollis(2L, numother, tlength);
2460 c[1] = coll3[2][k-1] * onecollis(3L, numother, tlength);
2461 normalizer = c[0] + c[1];
2462 c[0] /= normalizer;
2463 if (chance > c[0]) {
2464 cline++;
2465 if (!setlength(3L, numother, t->eventnode->tyme, tlength, prime))
2466 succeeded = false;
2467 }
2468 }
2469 if (t == NULL) {
2470 fprintf(ERRFILE,"ERROR in slide, ran off tymelist!\n");
2471 exit(-1);
2472 }
2473 if (t->eventnode->number != oldp->number &&
2474 t->eventnode->number != oldans->number)
2475 numother++;
2476 k--;
2477 t = t->prev;
2478 }
2479 cline--; /* A lineage dies! */
2480 numother++;
2481 if (cline == 0) {
2482 printf("ERROR in slide, no lineages left!");
2483 printf(" cline = 0, too few uncollisions\n");
2484 printf(" in loop %12ld\n", indecks);
2485 exit(-1);
2486 }
2487 }
2488 k = i;
2489 t = tplus->prev;
2490 while (cline == 1 && k != 0 && t != NULL) {
2491 tlength = howlong(t);
2492 if ((numother * 6 + 6) * tlength / theta0 > 20.0)
2493 tlength = 20 * theta0 / (numother * 6 + 6);
2494 chance = randum();
2495 c[0] = coll2[0][k-1] * zerocollis(1L, numother, tlength);
2496 c[1] = coll2[1][k-1] * onecollis(2L, numother, tlength);
2497 normalizer = c[0] + c[1];
2498 c[0] /= normalizer;
2499 if (chance > c[0]) {
2500 cline++;
2501 if (!setlength(2L, numother, t->eventnode->tyme, tlength,prime))
2502 succeeded = false;
2503 } else {
2504 k--;
2505 t = t->prev;
2506 }
2507 if (t->eventnode->number != oldp->number &&
2508 t->eventnode->number != oldans->number)
2509 numother++;
2510 }
2511 }
2512 if (succeeded) /* all the setlength operations worked */
2513 accept_slide = accept_tree(prime, primans, ans, pansdaught, newbr,
2514 tplus);
2515 else /* a setlength failed, garbage present! */
2516 accept_slide = false;
2517
2518 slid++;
2519 if (!skipped) {
2520 free(coll2[0]);
2521 free(coll2[1]);
2522 free(coll2);
2523 free(coll3[0]);
2524 free(coll3[1]);
2525 free(coll3[2]);
2526 free(coll3);
2527 }
2528 if (accept_slide) {
2529 slacc++;
2530 updatebranch(oldans,oldp,prime);
2531 updateslide(oldp, false);
2532 freenode(oldp);
2533 updateslide(oldans, false);
2534 freenode(oldans);
2535 updateslide(primans, true);
2536 updateslide(prime, true);
2537 return (accept_slide);
2538 }
2539 hookup(oldbr[0], oldp->next);
2540 hookup(oldbr[1], oldp->next->next);
2541 if (!oldans->next->top) {
2542 hookup(oldbr[2], oldans->next);
2543 hookup(oldans->next->next, ans);
2544 } else {
2545 hookup(oldbr[2], oldans->next->next);
2546 hookup(oldans->next, ans);
2547 }
2548 p = oldp;
2549 for (i = 1; i <= 2; i++) {
2550 p = p->next;
2551 p->back->v = p->v;
2552 }
2553 p = oldans;
2554 for (i = 1; i <= 2; i++) {
2555 p = p->next;
2556 p->back->v = p->v;
2557 }
2558 localevaluate(oldp,oldbr[2]);
2559 freenode(prime);
2560 freenode(primans);
2561
2562 curtree->likelihood = oldlikelihood;
2563
2564 return (accept_slide);
2565 } /* slide */
2566
maketree()2567 void maketree()
2568 {
2569 long incrprog, i, metout, progout;
2570 double bestlike;
2571 boolean chainend, runend, changetree;
2572 char *chainlit[2];
2573
2574 chainlit[0] = "Short";
2575 chainlit[1] = "Long";
2576
2577 contribution = (contribarr *)calloc(sites,sizeof(contribarr));
2578 contribution[0] = (double *)calloc(sites*categs,sizeof(double));
2579 for (i=1;i<sites;i++)
2580 contribution[i] = contribution[0] + i*categs;
2581
2582 inittable();
2583 initweightrat();
2584 getc(infile);
2585 fprintf(outfile, "Watterson estimate of theta is %12.8f\n", watttheta);
2586 if (op->usertree)
2587 treeread();
2588 else {
2589 branch0 = watttheta/numseq;
2590 constructtree(numseq, branch0);
2591 }
2592 orient(curtree->root->back);
2593 finishsetup(curtree->root->back);
2594 initbranchlist();
2595 treevaluate();
2596 bestlike = NEGMAX;
2597 growi[locus][0] = growth0;
2598 theti[locus][0] = theta0;
2599 lntheti[locus][0] = log(theta0);
2600 runend = false;
2601 /* We're going to start sampling thetas with tree 10, and resample after
2602 10 more trees have been outputted. */
2603 /**********************************/
2604 /* Begin Hastings-Metropolis loop */
2605 /**********************************/
2606 for (apps = 0; apps < totchains; apps++) {
2607 if (apps >= op->numchains[0]) chaintype = 1;
2608 else chaintype = 0;
2609 if (op->progress) {
2610 printf("%s chain %ld ",chainlit[chaintype],
2611 ((chaintype == 0) ? (apps + 1) : (apps + 1 - op->numchains[0])));
2612 fflush(stdout);
2613 }
2614 metout = op->increm[chaintype] - 1;
2615 incrprog = (long)(op->steps[chaintype] / 10.0);
2616 progout = incrprog - 1;
2617 op->numout[chaintype] = 0;
2618 xinterval = 0.0; /* initialize largest interval */
2619 slacc = 0;
2620 slid = 0;
2621 chainend = false;
2622 changetree = true; /* just to be careful */
2623
2624 for (indecks=0; indecks < op->steps[chaintype]; indecks++) {
2625 oldlikelihood = curtree->likelihood;
2626 col = 0; /* column number, used in treeout */
2627 if (slide()) changetree = true;
2628 if (indecks == op->steps[chaintype] - 1) { /* end of chain? */
2629 chainend = true;
2630 if (apps == totchains - 1) /* end of run? */
2631 runend = true;
2632 }
2633 if (curtree->likelihood > bestlike) {
2634 if (onebestree) {
2635 FClose(bestree);
2636 bestree = fopen("bestree","w+");
2637 fprintf(bestree, "Chain #%2ld (%s) Step:%8ld\n",apps+1,
2638 chainlit[chaintype], indecks+1);
2639 treeout(curtree->root->back, 1L, &bestree);
2640 fprintf(bestree, " [%12.10f]\n", curtree->likelihood);
2641 bestlike = curtree->likelihood;
2642 }
2643 else {
2644 fprintf(bestree, "Chain #%2ld (%s) Step:%8ld\n",apps+1,
2645 chainlit[chaintype], indecks+1);
2646 treeout(curtree->root->back, 1L, &bestree);
2647 fprintf(bestree, " [%12.10f]\n", curtree->likelihood);
2648 bestlike = curtree->likelihood;
2649 }
2650 }
2651 if (indecks == metout) {
2652 if (op->numout[chaintype] == 0) sametree[locus][0] = false;
2653 else sametree[locus][op->numout[chaintype]] = !changetree;
2654 changetree = false;
2655 op->numout[chaintype]++;
2656 scoretree(apps);
2657 metout += op->increm[chaintype];
2658 if (op->treeprint && apps == totchains-1) {
2659 fprintf(treefile,"\nlocus = %ld, chain = %ld, tree = %ld,",
2660 locus,apps,indecks);
2661 fprintf(treefile,"likelihood = %e\n",curtree->likelihood);
2662 treeout(curtree->root->back, 1L, &treefile);
2663 }
2664 }
2665 if (op->progress && indecks == progout) {
2666 printf(".");
2667 fflush(stdout);
2668 progout += incrprog;
2669 }
2670 }
2671 if (op->progress) printf("\nAccepted %ld/%ld rearrangements\n",slacc,slid);
2672 if (!op->growthused) {
2673 theta0 = coal_singlechain(apps, chainend, runend);
2674 } else {
2675 fluc_estimate(apps,runend);
2676 theta0 = theti[locus][apps+1];
2677 growth0 = growi[locus][apps+1];
2678 }
2679 }
2680 if(slacc == 0) {
2681 fprintf(outfile,"WARNING--no proposed trees ever accepted\n");
2682 fprintf(ERRFILE,"WARNING--no proposed trees ever accepted\n");
2683 }
2684 } /* maketree */
2685
finalfree()2686 void finalfree()
2687 /* free everything at end of program; a debugging convenience to check
2688 for memory leaks */
2689 {
2690 free(rate);
2691 free(probcat);
2692 free(ne_ratio);
2693 free(mu_ratio);
2694 free(theta_ratio);
2695 free(sametree[0]);
2696 free(sametree);
2697 free(curtree);
2698 freedata(dna);
2699 } /* finalfree */
2700
main(int argc,char * argv[])2701 int main(int argc, char *argv[])
2702 { /* Fluctuate */
2703 char infilename[100],outfilename[100],trfilename[100],logfilename[100],
2704 thfilename[100],intrfilename[100],bestrfilename[100];
2705
2706 double clearseed;
2707 long i;
2708
2709 #ifdef MAC
2710 /*macsetup("Fluctuate","Fluctuate");*/
2711 argv[0] = "Fluctuate";
2712 #endif
2713
2714 /* Open various filenames. */
2715
2716 openfile(&infile,INFILE,"r",argv[0],infilename);
2717 if (!menu)
2718 {
2719 openfile(&simlog,"simlog","w",argv[0],logfilename);
2720 }
2721 if (thetaout) openfile(&thetafile,"thetafile","w",argv[0],thfilename);
2722 openfile(&outfile,OUTFILE,"w",argv[0],outfilename);
2723
2724 op = (option_struct *)calloc(1,sizeof(option_struct));
2725 op->ibmpc = IBMCRT;
2726 op->ansi = ANSICRT;
2727 getoptions();
2728 /* start up randum numbers */
2729 for (i = 1; i <= 1000; i++)
2730 clearseed = randum();
2731 if (op->usertree)
2732 openfile(&intree,INTREE,"r",argv[0],intrfilename);
2733 if (op->treeprint)
2734 openfile(&treefile,TREEFILE,"w",argv[0],trfilename);
2735 openfile(&bestree,"bestree","w",argv[0],bestrfilename);
2736 firstinit();
2737 for (locus = 0; locus < numloci; locus++) {
2738 if (op->progress) printf("Locus %ld\n",locus+1);
2739 fprintf(outfile,"\n---------------------------------------\n");
2740 fprintf(outfile, "Locus %ld\n",locus+1);
2741 fprintf(outfile,"---------------------------------------\n");
2742 fprintf(outfile,"Assumed relative Ne = %f\n",ne_ratio[locus]);
2743 fprintf(outfile,"Assumed relative mu = %f\n",mu_ratio[locus]);
2744 if (holding)
2745 if (holding == 1)
2746 fprintf(outfile,"Theta was held constant.\n");
2747 else
2748 fprintf(outfile,"Growth-rate was held constant.\n");
2749 locusinit();
2750 getinput();
2751 watttheta = watterson();
2752 if (op->watt)
2753 theta0 = watttheta;
2754 maketree();
2755 freetree();
2756 }
2757 if (!op->growthused) coal_curveplot();
2758 if (numloci > 1) {
2759 fluc_locus_estimate();
2760 }
2761 modelfree();
2762 finalfree();
2763 FClose(infile);
2764 FClose(outfile);
2765 FClose(treefile);
2766 FClose(bestree);
2767 FClose(simlog);
2768 FClose(parmfile);
2769 FClose(thetafile);
2770 #ifdef MAC
2771 fixmacfile(outfilename);
2772 fixmacfile(trfilename);
2773 fixmacfile(bestrfilename);
2774 fixmacfile(intrfilename);
2775 fixmacfile(thfilename);
2776 fixmacfile(logfilename);
2777 #endif
2778 printf("PROGRAM DONE\n");
2779 exit(0);
2780 } /* Fluctuate */
2781
eof(FILE * f)2782 int eof(FILE *f)
2783 {
2784 register int ch;
2785
2786 if (feof(f))
2787 return 1;
2788 if (f == stdin)
2789 return 0;
2790 ch = getc(f);
2791 if (ch == EOF)
2792 return 1;
2793 ungetc(ch, f);
2794 return 0;
2795 } /* eof */
2796
eoln(FILE * f)2797 int eoln(FILE *f)
2798 {
2799 register int ch;
2800
2801 ch = getc(f);
2802 if (ch == EOF)
2803 return 1;
2804 ungetc(ch, f);
2805 return (ch == '\n');
2806 } /* eoln */
2807
memerror()2808 void memerror()
2809 {
2810 printf("Error allocating memory\n");
2811 exit(-1);
2812 } /* memerror */
2813