1 /*
2 * ml2.c
3 *
4 *
5 * Part of TREE-PUZZLE 5.2 (July 2004)
6 *
7 * (c) 2003-2004 by Heiko A. Schmidt, Korbinian Strimmer, and Arndt von Haeseler
8 * (c) 1999-2003 by Heiko A. Schmidt, Korbinian Strimmer,
9 * M. Vingron, and Arndt von Haeseler
10 * (c) 1995-1999 by Korbinian Strimmer and Arndt von Haeseler
11 *
12 * All parts of the source except where indicated are distributed under
13 * the GNU public licence. See http://www.opensource.org for details.
14 *
15 * ($Id$)
16 *
17 */
18
19 #ifdef HAVE_CONFIG_H
20 # include <config.h>
21 #endif
22
23
24 #define EXTERN extern
25
26 /* prototypes */
27 #include <stdio.h>
28 #include <stdlib.h>
29 #include <math.h>
30 #include <ctype.h>
31 #include <string.h>
32 #include "util.h"
33 #include "ml.h"
34
35 #define STDOUT stdout
36 #ifndef PARALLEL /* because printf() runs significantly faster */
37 /* than fprintf(stdout) on an Apple McIntosh */
38 /* (HS) */
39 # define FPRINTF printf
40 # define STDOUTFILE
41 #else
42 # define FPRINTF fprintf
43 # define STDOUTFILE STDOUT,
44 void PP_Finalize();
45 #endif
46
47 #define CHAR_PER_TREENODE 53
48
49 /* prototypes for two functions of puzzle2.c */
50 void fputid10(FILE *, int);
51 int fputid(FILE *, int);
52
53
54 /******************************************************************************/
55 /* user tree input */
56 /******************************************************************************/
57
58 /* read/convert branch length from string stream itfp */
freadedgelength(double * length,char * ch,FILE * itfp)59 double freadedgelength(double *length, char *ch, FILE *itfp)
60 {
61 long digit, ordzero;
62 int dotread;
63 double divisor;
64 int minusread;
65
66 ordzero = '0';
67 dotread = FALSE;
68 minusread = FALSE;
69 *length = 0.0;
70 divisor = 1.0;
71 *ch = fgetc(itfp);
72 digit = (long)(*ch - ordzero);
73 while ( ((digit <= 9) && (digit >= 0)) || *ch == '.' || *ch == '-' ) {
74 if (*ch == '.' )
75 dotread = TRUE;
76 else if (*ch == '-' )
77 minusread = TRUE;
78 else {
79 *length = *length * 10.0 + digit;
80 if (dotread)
81 divisor *= 10.0;
82 }
83 ch++;
84 *ch = fgetc(itfp);
85 digit = (long)(*ch - ordzero);
86 }
87
88 *length *= divisor;
89
90 if (minusread)
91 *length = -(*length);
92
93 return (*length);
94
95 } /* freadedgelength */
96
97
98 /* read/convert branch length from string *ch */
sreadedgelength(double * length,char ** ch)99 double sreadedgelength(double *length, char **ch)
100 {
101 long digit, ordzero;
102 int dotread;
103 double divisor;
104 int minusread;
105
106 ordzero = '0';
107 dotread = FALSE;
108 minusread = FALSE;
109 *length = 0.0;
110 divisor = 1.0;
111 /* **ch = fgetc(itfp); */
112 if (**ch == ':') (*ch)++;
113 digit = (long)((**ch) - ordzero);
114 while ( ((digit <= 9) && (digit >= 0)) || **ch == '.' || **ch == '-' ) {
115 if (**ch == '.' )
116 dotread = TRUE;
117 else if (**ch == '-' )
118 minusread = TRUE;
119 else {
120 *length = *length * 10.0 + digit;
121 if (dotread)
122 divisor *= 10.0;
123 }
124 (*ch)++;
125 /* *ch = fgetc(itfp); */
126 digit = (long)((**ch) - ordzero);
127 }
128
129 *length *= 100;
130 *length /= divisor;
131
132 if (minusread)
133 *length = -(*length);
134
135 return (*length);
136
137 } /* sreadedgelength */
138
139
fremovecomment(FILE * itfp,char * ci)140 void fremovecomment(FILE *itfp, char *ci)
141 {
142 int comment;
143
144 if (*ci == '[') comment = 1;
145 else {
146 comment = 0;
147 return;
148 }
149
150 do {
151 *ci = fgetc(itfp);
152 if (*ci == EOF) {
153 fprintf(STDOUT, "\n\n\nUnable to proceed (missing ';' or ',' in tree)\n\n\n");
154 exit(1);
155 }
156 if ((*ci == ']') && comment) {
157 comment = 0;
158 *ci = fgetc(itfp);
159 }
160 if (*ci == '[') comment = 1;
161 } while (comment);
162
163 } /* fremovecomment */
164
165
166
167
168
169 /* read user tree, drop all blanks, tabs, and newlines.
170 Drop edgelengths (after :) but keep internal
171 labels. Check whether all pairs of brackets match. */
getusertree(FILE * itfp,cvector tr,int maxlen,int usebranch)172 void getusertree(FILE *itfp, cvector tr, int maxlen, int usebranch)
173 {
174 int n;
175 int brac;
176 char ci;
177
178 /* look for opening bracket */
179 n = 0;
180 brac = 0;
181 do {
182 ci = fgetc(itfp);
183 if (ci == EOF) {
184 fprintf(STDOUT, "\n\n\nUnable to proceed (missing start bracket in tree)\n\n\n");
185 # if PARALLEL
186 PP_Finalize();
187 # endif
188 exit(1);
189 }
190 if (ci == '[') fremovecomment(itfp, &ci);
191 } while ((ci != '('));
192 tr[n] = ci;
193 brac++;
194
195 do {
196 /* get next character (skip blanks, newlines, and tabs) */
197 do {
198 ci = fgetc(itfp);
199 if (ci == EOF) {
200 fprintf(STDOUT, "\n\n\nUnable to proceed (no more characters in tree)\n\n\n");
201 # if PARALLEL
202 PP_Finalize();
203 # endif
204 exit(1);
205 }
206 if (ci == '[') fremovecomment(itfp, &ci);
207 } while (isspace((int) ci));
208
209 /* remove edgelengths, if not needed (HAS) */
210 if (ci == ':') { /* skip characters until a ,) appears */
211 if (!usebranch) {
212 do {
213 ci = fgetc(itfp);
214 if (ci == EOF) {
215 fprintf(STDOUT, "\n\n\nUnable to proceed (missing ';' or ',' in tree)\n\n\n");
216 # if PARALLEL
217 PP_Finalize();
218 # endif
219 exit(1);
220 }
221 if (ci == '[') fremovecomment(itfp, &ci);
222 } while ((ci != ',' && ci != ')') );
223
224 } else { /* if usebranch */
225 n++;
226 tr[n] = ci; /* add ':' to string */
227 do {
228 ci = fgetc(itfp);
229 if (ci == EOF) {
230 fprintf(STDOUT, "\n\n\nUnable to proceed (missing ';' or ',' in tree)\n\n\n");
231 exit(1);
232 }
233 if (ci == '[') fremovecomment(itfp, &ci);
234 if(!isspace((int) ci) && (ci != ',') && (ci != ')')) {
235 n++;
236 tr[n] = ci;
237 };
238 } while ((ci != ',' && ci != ')') );
239 }
240 }
241
242 if (ci == '(') {
243 brac++;
244 }
245 if (ci == ')') {
246 brac--;
247 }
248
249 n++;
250 tr[n] = ci;
251
252 } while ((ci != ';') && (n != maxlen-2));
253
254 if (n == maxlen-2) {
255 fprintf(STDOUT, "\n\n\nUnable to proceed (tree description too long)\n\n\n");
256 fprintf(STDOUT, ">>%s<<\n\n\n", tr);
257
258 # if PARALLEL
259 PP_Finalize();
260 # endif
261 exit(1);
262 }
263
264 if (brac != 0) {
265 fprintf(STDOUT, "\n\n\nUnable to proceed (brackets don't match in tree)\n\n\n");
266 # if PARALLEL
267 PP_Finalize();
268 # endif
269 exit(1);
270 }
271
272 n++;
273 tr[n] = '\0';
274 } /* getusertree */
275
276
internalnode(Tree * tr,char ** chpp,int * ninode)277 Node *internalnode(Tree *tr, char **chpp, int *ninode)
278 {
279 Node *xp;
280 Node *np;
281 Node *rp;
282 int i, j;
283 int dvg;
284 int ff;
285 int stop;
286 int numc;
287 char ident[100];
288 char idcomp[11];
289 char *idp;
290 double blen;
291
292 (*chpp)++;
293 if (**chpp == '(') { /* process subgroup */
294
295 xp = internalnode(tr, chpp, ninode);
296 xp->isop = xp;
297 dvg = 1;
298 while (**chpp != ')') {
299 if (**chpp == '\0') {
300 fprintf(STDOUT, "\n\n\nUnable to proceed (unexpected end of tree)\n\n\n");
301 # if PARALLEL
302 PP_Finalize();
303 # endif
304 exit(1);
305 }
306 dvg++;
307 /* insert edges around node */
308 np = internalnode(tr, chpp, ninode);
309 np->isop = xp->isop;
310 xp->isop = np;
311 xp = np;
312 }
313 /* closing bracket reached */
314 /* read edgelength !!!! (HAS) */
315
316 (*chpp)++;
317 if (dvg < 2) {
318 fprintf(STDOUT, "\n\n\nUnable to proceed (only one OTU inside pair of brackets)\n\n\n");
319 # if PARALLEL
320 PP_Finalize();
321 # endif
322 exit(1);
323 }
324
325 if ((*ninode) >= Maxspc-3) { /* all internal nodes already used */
326 fprintf(STDOUT, "\n\n\nUnable to proceed (no unrooted tree)\n\n\n");
327 # if PARALLEL
328 PP_Finalize();
329 # endif
330 exit(1);
331 }
332
333 rp = tr->ibrnchp[*ninode];
334 rp->isop = xp->isop;
335 xp->isop = rp;
336
337 for (j = 0; j < Numspc; j++)
338 rp->paths[j] = 0;
339 xp = rp->isop;
340 while (xp != rp) {
341 for (j = 0; j < Numspc; j++) {
342 if (xp->paths[j] == 1)
343 rp->paths[j] = 1;
344 }
345 xp = xp->isop;
346 }
347 (*ninode)++;
348
349 if ((**chpp) == ',' || (**chpp) == ')') return rp->kinp;
350 if ((**chpp) == '\0') {
351 fprintf(STDOUT, "\n\n\nUnable to proceed (unexpected end of tree)\n\n\n");
352 # if PARALLEL
353 PP_Finalize();
354 # endif
355 exit(1);
356 }
357
358 if ((**chpp) != ':') {
359
360 /* read internal label into rp->label (max. 20 characters) */
361 rp->label = new_cvector(21);
362 (rp->label)[0] = '\0';
363 for (numc = 0; ((numc < 20) && ((**chpp) != ':')); numc++) {
364 if ((**chpp) == ',' || (**chpp) == ')') return rp->kinp;
365 if ((**chpp) == '\0') {
366 fprintf(STDOUT, "\n\n\nUnable to proceed (unexpected end of tree)\n\n\n");
367 exit(1);
368 }
369 (rp->label)[numc] = **chpp;
370 (rp->label)[numc+1] = '\0';
371 (*chpp)++;
372 }
373
374 /* skip the rest of the internal label */
375 while (((**chpp) != ',') && ((**chpp) != ')') && ((**chpp) != ':')) {
376 (*chpp)++;
377 if ((**chpp) == '\0') {
378 fprintf(STDOUT, "\n\n\nUnable to proceed (unexpected end of tree)\n\n\n");
379 exit(1);
380 }
381 }
382 } /* if ! ':' */
383
384 /* read edgelength !!!! (HAS) */
385 if (**chpp == ':') {
386 (void) sreadedgelength(&blen, chpp);
387 rp->lengthext=blen;
388 rp->lengthset=1;
389 rp->kinp->lengthext=blen;
390 rp->kinp->lengthset=1;
391 }
392
393
394 return rp->kinp;
395
396 } else { /* process species names */
397 /* read species name */
398 for (idp = ident;
399 **chpp != ',' && **chpp != ')' && **chpp != ':' && **chpp != '\0';
400 (*chpp)++) {
401 *idp++ = **chpp;
402 }
403 *idp = '\0';
404
405 /* look for internal number */
406 idcomp[10] = '\0';
407
408 for (i = 0; i < Maxspc; i++) {
409 ff = 0;
410 stop = FALSE;
411 do {
412 idcomp[ff] = Identif[i][ff];
413 ff++;
414 if (idcomp[ff-1] == ' ') stop = TRUE;
415 } while (!stop && (ff != 10));
416 if (stop) idcomp[ff-1] = '\0';
417
418 if (!strcmp(ident, idcomp)) {
419 if (usedtaxa[i]) {
420 fprintf(STDOUT, "\n\n\nUnable to proceed (multiple occurence of sequence '");
421 fprintf(STDOUT, "%s' in tree)\n\n\n", ident);
422 # if PARALLEL
423 PP_Finalize();
424 # endif
425 exit(1);
426 }
427 usedtaxa[i] = TRUE;
428
429 /* read edgelength !!!! (HAS) */
430 if (**chpp == ':') {
431 (void) sreadedgelength(&blen, chpp);
432 tr->ebrnchp[i]->lengthext=blen;
433 tr->ebrnchp[i]->lengthset=1;
434 tr->ebrnchp[i]->kinp->lengthext=blen;
435 tr->ebrnchp[i]->kinp->lengthset=1;
436 }
437
438 return tr->ebrnchp[i]->kinp;
439 }
440 }
441
442 fprintf(STDOUT, "\n\n\nUnable to proceed (unknown sequence '%s' in tree)\n\n\n", ident);
443 # if PARALLEL
444 PP_Finalize();
445 # endif
446 exit(1);
447 }
448 return NULL; /* never returned but without some compilers complain */
449 } /* internalnode */
450
451
452 /* make tree structure, the tree description may contain internal
453 labels but no edge lengths */
constructtree(Tree * tr,cvector strtree,int usebranch)454 void constructtree(Tree *tr, cvector strtree, int usebranch)
455 {
456 char *chp; /* char ptr to tree string */
457 int ninode; /* number of internal nodes < Maxspc-3 */
458 int i; /* counter */
459 int dvg; /* degree vertex ? */
460 int numc; /* char counter */
461 Node *xp; /* last read Node/Subtree */
462 Node *np; /* newly read Node/Subtree */
463
464 ninode = 0;
465 chp = strtree;
466 usedtaxa = new_ivector(Maxspc);
467 for (i = 0; i < Maxspc; i++) usedtaxa[i] = FALSE;
468
469 xp = internalnode(tr, &chp, &ninode);
470 xp->isop = xp;
471 dvg = 1;
472 while (*chp != ')') { /* look for closing bracket */
473 if (*chp == '\0') {
474 fprintf(STDOUT, "\n\n\nUnable to proceed (unexpected end of tree)\n\n\n");
475 # if PARALLEL
476 PP_Finalize();
477 # endif
478 exit(1);
479 }
480 dvg++;
481 /* insert edges around node */
482 np = internalnode(tr, &chp, &ninode);
483 np->isop = xp->isop;
484 xp->isop = np;
485 xp = np;
486 }
487
488 for (i = 0; i < Maxspc; i++)
489 if (usedtaxa[i] == FALSE) {
490 fprintf(STDOUT, "\n\n\nUnable to proceed (sequences missing in tree)\n\n\n");
491 # if PARALLEL
492 PP_Finalize();
493 # endif
494 exit(1);
495 }
496
497 /* closing bracket reached */
498 if (dvg < 3) {
499 fprintf(STDOUT, "\n\n\nUnable to proceed (no unrooted tree)\n\n\n");
500 # if PARALLEL
501 PP_Finalize();
502 # endif
503 exit(1);
504 }
505 tr->rootp = xp;
506 Numibrnch = ninode;
507 Numbrnch = Numspc + ninode;
508
509 chp++;
510 if (*chp == ';' || *chp == '\0') {
511 free_ivector(usedtaxa);
512 return;
513 }
514
515 /* copy last internal label (max. 20 characters) */
516 xp->label = new_cvector(21);
517 (xp->label)[0] = *chp;
518 (xp->label)[1] = '\0';
519 for (numc = 1; numc < 20; numc++) {
520 chp++;
521 if (*chp == ';' || *chp == '\0') {
522 free_ivector(usedtaxa);
523 return;
524 } else {
525 (xp->label)[numc] = *chp;
526 (xp->label)[numc+1] = '\0';
527 }
528 }
529 free_ivector(usedtaxa);
530 return;
531 } /* constructtree */
532
533
534 /* remove possible basal bifurcation */
removebasalbif(cvector strtree)535 void removebasalbif(cvector strtree)
536 {
537 int n, c, brak, cutflag, h;
538
539 /* check how many OTUs on basal level */
540 n = 0;
541 c = 0;
542 brak = 0;
543 do {
544 if (strtree[n] == '(') brak++;
545 if (strtree[n] == ')') brak--;
546
547 if (strtree[n] == ',' && brak == 1) c++; /* number of commas in outer bracket */
548
549 n++;
550 } while (strtree[n] != '\0');
551
552 /* if only 1 OTU inside outer bracket stop now */
553 if (c == 0) {
554 fprintf(STDOUT, "\n\n\nUnable to proceed (Only 1 OTU inside outer bracket in tree)\n\n\n");
555 # if PARALLEL
556 PP_Finalize();
557 # endif
558 exit(1);
559 }
560
561 /* if only 2 OTUs inside outer bracket delete second pair of
562 brackets from the right to remove basal bifurcation */
563
564 if (c == 1) {
565
566 n = 0;
567 brak = 0;
568 cutflag = 0; /* not yet cutted */
569 h = 0;
570 do {
571 if (strtree[n] == '(') brak++;
572 if (strtree[n] == ')') brak--;
573
574 if (brak == 2 && cutflag == 0) cutflag = 1; /* cutting */
575 if (brak == 1 && cutflag == 1) {
576 cutflag = 2; /* cutted */
577 /* leave out internal label */
578 do {
579 h++;
580 } while (strtree[n+h] != ')' && strtree[n+h] != ',');
581
582 }
583
584 if (cutflag == 1) strtree[n] = strtree[n+1];
585 if (cutflag == 2) strtree[n-1] = strtree[n+h];
586
587 n++;
588 } while (strtree[n] != '\0');
589 }
590 } /* removebasalbif */
591
592
makeusertree(FILE * itfp,int usebranchlen)593 void makeusertree(FILE *itfp, int usebranchlen)
594 {
595 cvector strtree;
596
597 strtree = new_cvector(CHAR_PER_TREENODE*2*Maxspc); /* for treefile */
598 getusertree(itfp, strtree, CHAR_PER_TREENODE*2*Maxspc, usebranchlen);
599 removebasalbif(strtree);
600 constructtree(Ctree, strtree, usebranchlen);
601 free_cvector(strtree);
602 } /* makeusertree */
603
604
605 /******************************************************************************/
606 /* memory organisation for maximum likelihood tree */
607 /******************************************************************************/
608
609 /* initialise new tree */
new_tree(int maxspc,int numptrn,cmatrix seqconint)610 Tree *new_tree(int maxspc, int numptrn, cmatrix seqconint)
611 {
612 int n, i, maxibrnch;
613 Tree *tr;
614 Node *dp, *up;
615
616 maxibrnch = maxspc - 3;
617 heights = (Node **) calloc((size_t)(maxspc-2), sizeof(Node *));
618 if (heights == NULL) maerror("heights in new_tree");
619 tr = (Tree *) calloc((size_t)1, sizeof(Tree));
620 if (tr == NULL) maerror("tr in new_tree");
621 tr->ebrnchp = (Node **) calloc((size_t)maxspc, sizeof(Node *));
622 if (tr->ebrnchp == NULL) maerror("ebrnchp in new_tree");
623 tr->ibrnchp = (Node **) calloc((size_t)maxibrnch, sizeof(Node *));
624 if (tr->ibrnchp == NULL) maerror("ibrnchp in new_tree");
625 tr->condlkl = new_dmatrix(numcats, numptrn);
626 for (n = 0; n < maxspc; n++) {
627 dp = (Node *) calloc((size_t) 1, sizeof(Node));
628 if (dp == NULL) maerror("dp in new_tree");
629 up = (Node *) calloc((size_t) 1, sizeof(Node));
630 if (up == NULL) maerror("up in new_tree");
631 dp->isop = NULL;
632 up->isop = NULL;
633 dp->kinp = up;
634 up->kinp = dp;
635 dp->descen = TRUE;
636 up->descen = FALSE;
637 dp->number = n;
638 up->number = n;
639 dp->length = 0.0;
640 up->length = 0.0;
641 dp->lengthc = 0.0;
642 up->lengthc = 0.0;
643 dp->lengthext = 0.0;
644 up->lengthext = 0.0;
645 dp->lengthset = FALSE;
646 up->lengthset = FALSE;
647 dp->varlen = 0.0;
648 up->varlen = 0.0;
649 dp->paths = new_ivector(maxspc);
650 up->paths = dp->paths;
651 for (i = 0; i < maxspc; i++) dp->paths[i] = 0;
652 dp->paths[n] = 1;
653 dp->eprob = seqconint[n];
654 up->eprob = NULL;
655 dp->partials = NULL;
656 up->partials = new_dcube(numcats, numptrn, tpmradix);
657 tr->ebrnchp[n] = dp;
658 up->label = NULL;
659 dp->label = NULL;
660 }
661 for (n = 0; n < maxibrnch; n++) {
662 dp = (Node *) calloc((size_t) 1, sizeof(Node));
663 if (dp == NULL) maerror("dp in new_tree");
664 up = (Node *) calloc((size_t) 1, sizeof(Node));
665 if (up == NULL) maerror("up in new_tree");
666 dp->isop = NULL;
667 up->isop = NULL;
668 dp->kinp = up;
669 up->kinp = dp;
670 dp->descen = TRUE;
671 up->descen = FALSE;
672 dp->number = n;
673 up->number = n;
674 dp->length = 0.0;
675 up->length = 0.0;
676 dp->lengthc = 0.0;
677 up->lengthc = 0.0;
678 dp->lengthext = 0.0;
679 up->lengthext = 0.0;
680 dp->lengthset = FALSE;
681 up->lengthset = FALSE;
682 dp->varlen = 0.0;
683 up->varlen = 0.0;
684 dp->paths = new_ivector(maxspc);
685 up->paths = dp->paths;
686 for (i = 0; i < maxspc; i++) dp->paths[i] = 0;
687 dp->eprob = NULL;
688 up->eprob = NULL;
689 dp->partials = new_dcube(numcats, numptrn, tpmradix);
690 up->partials = new_dcube(numcats, numptrn, tpmradix);
691 tr->ibrnchp[n] = dp;
692 up->label = NULL;
693 dp->label = NULL;
694 }
695 tr->rootp = NULL;
696
697 /*
698 * reserve memory for lengths of the tree branches
699 * and for the distance matrix as a vector
700 * (needed for LS estimation of tree branch lengths)
701 */
702
703 Brnlength = new_dvector(2 * maxspc - 3);
704 Distanvec = new_dvector((maxspc * (maxspc - 1)) / 2);
705
706 return tr;
707 } /* new_tree */
708
709
710 /* initialise quartet tree */
new_quartet(int numptrn,cmatrix seqconint)711 Tree *new_quartet(int numptrn, cmatrix seqconint)
712 {
713 int n, i;
714 Tree *tr;
715 Node *dp, *up;
716
717 heights = (Node **) calloc((size_t) 2, sizeof(Node *));
718 if (heights == NULL) maerror("heights in new_quartet");
719 /* reserve memory for tree */
720 tr = (Tree *) calloc((size_t) 1, sizeof(Tree));
721 if (tr == NULL) maerror("tr in new_quartet");
722 tr->ebrnchp = (Node **) calloc((size_t) 4, sizeof(Node *));
723 if (tr->ebrnchp == NULL) maerror("ebrnchp in new_quartet");
724 tr->ibrnchp = (Node **) calloc((size_t) 1, sizeof(Node *));
725 if (tr->ibrnchp == NULL) maerror("ibrnchp in new_quartet");
726 tr->condlkl = new_dmatrix(numcats, numptrn);
727 /* reserve memory for nodes */
728 for (n = 0; n < 4; n++) {
729 dp = (Node *) calloc((size_t) 1, sizeof(Node));
730 if (dp == NULL) maerror("dp in new_quartet");
731 up = (Node *) calloc((size_t) 1, sizeof(Node));
732 if (up == NULL) maerror("dp in new_quartet");
733 dp->isop = NULL;
734 dp->kinp = up;
735 up->kinp = dp;
736 dp->descen = TRUE;
737 up->descen = FALSE;
738 dp->number = n;
739 up->number = n;
740 dp->length = 0.0;
741 up->length = 0.0;
742 dp->lengthc = 0.0;
743 up->lengthc = 0.0;
744 dp->lengthext = 0.0;
745 up->lengthext = 0.0;
746 dp->lengthset = FALSE;
747 up->lengthset = FALSE;
748 dp->varlen = 0.0;
749 up->varlen = 0.0;
750 dp->paths = new_ivector(4);
751 up->paths = dp->paths;
752 for (i = 0; i < 4; i++) dp->paths[i] = 0;
753 dp->paths[n] = 1;
754 dp->eprob = seqconint[n]; /* make quartet (0,1)-(2,3) as default */
755 up->eprob = NULL;
756 dp->partials = NULL;
757 up->partials = new_dcube(numcats, numptrn, tpmradix);
758 tr->ebrnchp[n] = dp;
759 }
760
761 /* reserve memory for internal branch */
762 dp = (Node *) calloc((size_t) 1, sizeof(Node));
763 if (dp == NULL) maerror("dp in new_quartet");
764 up = (Node *) calloc((size_t) 1, sizeof(Node));
765 if (up == NULL) maerror("dp in new_quartet");
766 dp->isop = tr->ebrnchp[3]->kinp; /* connect internal branch */
767 up->isop = tr->ebrnchp[0]->kinp;
768 dp->kinp = up;
769 up->kinp = dp;
770 dp->descen = TRUE;
771 up->descen = FALSE;
772 dp->number = 0;
773 up->number = 0;
774 dp->length = 0.0;
775 up->length = 0.0;
776 dp->lengthc = 0.0;
777 up->lengthc = 0.0;
778 dp->lengthext = 0.0;
779 up->lengthext = 0.0;
780 dp->lengthset = FALSE;
781 up->lengthset = FALSE;
782 dp->varlen = 0.0;
783 up->varlen = 0.0;
784 dp->paths = new_ivector(4);
785 up->paths = dp->paths;
786 up->paths[0] = 0;
787 up->paths[1] = 0;
788 up->paths[2] = 1;
789 up->paths[3] = 1;
790 dp->eprob = NULL;
791 up->eprob = NULL;
792 dp->partials = new_dcube(numcats, numptrn, tpmradix);
793 up->partials = new_dcube(numcats, numptrn, tpmradix);
794 tr->ibrnchp[0] = dp;
795
796 /* place root */
797 tr->rootp = up;
798
799 /* connect external branches */
800 tr->ebrnchp[0]->kinp->isop = tr->ebrnchp[1]->kinp;
801 tr->ebrnchp[1]->kinp->isop = tr->rootp;
802 tr->ebrnchp[3]->kinp->isop = tr->ebrnchp[2]->kinp;
803 tr->ebrnchp[2]->kinp->isop = tr->rootp->kinp;
804
805 /*
806 * reserve memory for lengths of the five branches
807 * of a quartet and for the six possible distances
808 * (needed for LS estimation of branch lengths)
809 */
810 Brnlength = new_dvector(NUMQBRNCH);
811 Distanvec = new_dvector(NUMQSPC*(NUMQSPC-1)/2);
812
813 return tr;
814 } /* new_quartet */
815
816
817 /* free tree memory */
free_tree(Tree * tr,int taxa)818 void free_tree(Tree *tr, int taxa)
819 {
820 int n;
821 Node *dp, *up;
822
823 free(heights);
824 free_dmatrix(tr->condlkl);
825 for (n = 0; n < taxa; n++) {
826 dp = tr->ebrnchp[n];
827 up = dp->kinp;
828 free_ivector(dp->paths);
829 free_dcube(up->partials);
830 free(dp);
831 free(up);
832 }
833 free(tr->ebrnchp);
834 for (n = 0; n < (taxa-3); n++) {
835 dp = tr->ibrnchp[n];
836 up = dp->kinp;
837 free_dcube(dp->partials);
838 free_dcube(up->partials);
839 free_ivector(dp->paths);
840 if (dp->label != NULL) free_cvector(dp->label);
841 if (up->label != NULL) free_cvector(up->label);
842 free(dp);
843 free(up);
844 }
845 free(tr->ibrnchp);
846 free(tr);
847 free_dvector(Brnlength); /* branch lengths (for LS estimation) */
848 free_dvector(Distanvec); /* distances (for LS estimation) */
849 } /* free_tree */
850
851
852 /* make (a,b)-(c,d) quartet
853
854 a ---+ +--- c
855 +-----+
856 b ---+ +--- d
857
858 species numbers range from 0 to Maxspc - 1 */
859
make_quartet(int a,int b,int c,int d)860 void make_quartet(int a, int b, int c, int d)
861 {
862 /* place sequences */
863 Ctree->ebrnchp[0]->eprob = Seqpat[a];
864 Ctree->ebrnchp[1]->eprob = Seqpat[b];
865 Ctree->ebrnchp[2]->eprob = Seqpat[c];
866 Ctree->ebrnchp[3]->eprob = Seqpat[d];
867
868 /* make distance vector */
869 Distanvec[0] = Distanmat[b][a];
870 Distanvec[1] = Distanmat[c][a];
871 Distanvec[2] = Distanmat[c][b];
872 Distanvec[3] = Distanmat[d][a];
873 Distanvec[4] = Distanmat[d][b];
874 Distanvec[5] = Distanmat[d][c];
875 } /* make_quartet */
876
877 /* write distance matrix as vector */
changedistan(dmatrix distanmat,dvector distanvec,int numspc)878 void changedistan(dmatrix distanmat, dvector distanvec, int numspc)
879 {
880 int i, j, k;
881
882 for (k = 0, i = 1; i < numspc; i++) {
883 for (j = 0; j < i; j++, k++)
884 distanvec[k] = distanmat[i][j];
885 }
886 } /* changedistan */
887
888
889 /******************************************************************************/
890 /* computation of maximum likelihood tree */
891 /******************************************************************************/
892
893
894 /* compute the likelihood for (a,b)-(c,d) quartet */
quartet_lklhd(int a,int b,int c,int d)895 double quartet_lklhd(int a, int b, int c, int d)
896 {
897 /* reserve memory for quartet if necessary */
898 if (mlmode != ML_QUART) { /* no quartet tree */
899 if (Ctree != NULL)
900 free_tree(Ctree, Numspc);
901 Ctree = new_quartet(Numptrn, Seqpat);
902 Numbrnch = NUMQBRNCH;
903 Numibrnch = NUMQIBRNCH;
904 Numspc = NUMQSPC;
905 mlmode = ML_QUART;
906 }
907
908 /* make (a,b)-(c,d) quartet */
909 make_quartet(a,b,c,d);
910
911 clockmode = 0; /* nonclocklike branch lengths */
912
913 /* least square estimate for branch length */
914 lslength(Ctree, Distanvec, Numspc, Numibrnch, Brnlength);
915
916 /* compute likelihood */
917 Ctree->lklhd = optlkl(Ctree);
918
919 return Ctree->lklhd;
920 } /* quartet_lklhd */
921
922
923 /* compute the approximate likelihood for (a,b)-(c,d) quartet */
quartet_alklhd(int a,int b,int c,int d)924 double quartet_alklhd(int a, int b, int c, int d)
925 {
926 /* reserve memory for quartet if necessary */
927 if (mlmode != ML_QUART) { /* no quartet tree */
928 if (Ctree != NULL)
929 free_tree(Ctree, Numspc);
930 Ctree = new_quartet(Numptrn, Seqpat);
931 Numbrnch = NUMQBRNCH;
932 Numibrnch = NUMQIBRNCH;
933 Numspc = NUMQSPC;
934 mlmode = ML_QUART;
935 }
936
937 /* make (a,b)-(c,d) quartet */
938 make_quartet(a,b,c,d);
939
940 clockmode = 0; /* nonclocklike branch lengths */
941
942 /* least square estimate for branch length */
943 lslength(Ctree, Distanvec, Numspc, Numibrnch, Brnlength);
944
945 /* compute likelihood */
946 Ctree->lklhd = treelkl(Ctree);
947
948 return Ctree->lklhd;
949 } /* quartet_alklhd */
950
951
952 /* read usertree from file to memory */
readusertree(FILE * ifp,int usebranch)953 void readusertree(FILE *ifp, int usebranch)
954 {
955 /* reserve memory for tree if necessary */
956 if (mlmode != ML_TREE) { /* no tree */
957 if (Ctree != NULL)
958 free_tree(Ctree, Numspc);
959 Ctree = new_tree(Maxspc, Numptrn, Seqpat);
960 Numbrnch = 2*Maxspc-3;
961 Numibrnch = Maxspc-3;
962 Numspc = Maxspc;
963 mlmode = ML_TREE;
964 }
965
966 /* read tree */
967 makeusertree(ifp, usebranch);
968 } /* readusertree */
969
970
971 /* compute the likelihood of a usertree */
usertree_lklhd(int usebranch,int parallel)972 double usertree_lklhd(int usebranch, int parallel)
973 {
974 /* be sure to have a usertree in memory and
975 to have pairwise distances computed */
976
977 clockmode = 0; /* nonclocklike branch lengths */
978
979 if(usebranch) {
980 /* set preset barnch lengths from usertree file */
981 setextlength(Ctree, Numspc, Numibrnch, Brnlength);
982
983 /* compute likelihood */
984 Ctree->lklhd = treelkl(Ctree);
985
986 } else {
987 /* least square estimate for branch length */
988 changedistan(Distanmat, Distanvec, Numspc);
989 /*epe*/
990 # if PARALLEL
991 if (parallel)
992 lslength_par(Ctree, Distanvec, Numspc, Numibrnch, Brnlength);
993 else
994 # endif
995 lslength(Ctree, Distanvec, Numspc, Numibrnch, Brnlength);
996
997
998 /* compute and optimize likelihood */
999 Ctree->lklhd = optlkl(Ctree);
1000 }
1001
1002 return Ctree->lklhd;
1003 } /* usertree_lklhd */
1004
1005
1006 /* compute the approximate likelihood of a usertree */
usertree_alklhd(int usebranch,int parallel)1007 double usertree_alklhd(int usebranch, int parallel)
1008 {
1009 /* be sure to have a usertree in memory and
1010 to have pairwise distances computed */
1011
1012 clockmode = 0; /* nonclocklike branch lengths */
1013
1014 if(usebranch) {
1015 /* set preset barnch lengths from usertree file */
1016 setextlength(Ctree, Numspc, Numibrnch, Brnlength);
1017
1018 } else {
1019 /* least square estimate for branch length */
1020 changedistan(Distanmat, Distanvec, Numspc);
1021 /*epe*/
1022 # if PARALLEL
1023 if (parallel)
1024 lslength_par(Ctree, Distanvec, Numspc, Numibrnch, Brnlength);
1025 else
1026 # endif
1027 lslength(Ctree, Distanvec, Numspc, Numibrnch, Brnlength);
1028
1029 }
1030
1031 /* compute likelihood */
1032 Ctree->lklhd = treelkl(Ctree);
1033
1034 return Ctree->lklhd;
1035 } /* usertree_alklhd */
1036
1037
1038 /* preparation for ML analysis */
mlstart()1039 void mlstart()
1040 {
1041 /* number of states and code length */
1042 tpmradix = gettpmradix();
1043
1044 /* declare variables */
1045 Eval = new_dvector(tpmradix);
1046 Evec = new_dmatrix(tpmradix,tpmradix);
1047 Ievc = new_dmatrix(tpmradix,tpmradix);
1048 iexp = new_dmatrix(tpmradix,tpmradix);
1049 Alias = new_ivector(Maxsite);
1050
1051 /* process sequence information */
1052 evaluateseqs();
1053 bestrate = new_ivector(Numptrn);
1054
1055 /* compute transition probability matrix */
1056 tranprobmat();
1057
1058 /* non-zero rate categories */
1059 Rates = new_dvector(numcats);
1060 updaterates();
1061 ltprobr = new_dcube(numcats, tpmradix,tpmradix);
1062
1063 /* compute distance matrix */
1064 Distanmat = new_dmatrix(Maxspc, Maxspc);
1065 initdistan();
1066
1067 /* initialize tree pointer for quartet tree */
1068 mlmode = ML_QUART;
1069 Ctree = new_quartet(Numptrn, Seqpat);
1070 Numbrnch = NUMQBRNCH;
1071 Numibrnch = NUMQIBRNCH;
1072 Numspc = NUMQSPC;
1073
1074 /* computing ML distances */
1075 /*epe computedistan() moved to inputandinit;*/
1076
1077 } /* mlstart */
1078
1079
1080 /* recompute ml distances for quartet only */
distupdate(int a,int b,int c,int d)1081 void distupdate(int a, int b, int c, int d)
1082 {
1083 /* update distance matrix */
1084 /* consider only entries relevant to quartet */
1085 Distanmat[a][b] = mldistance(a, b, Distanmat[a][b]);
1086 Distanmat[b][a] = Distanmat[a][b];
1087 Distanmat[a][c] = mldistance(a, c, Distanmat[a][c]);
1088 Distanmat[c][a] = Distanmat[a][c];
1089 Distanmat[a][d] = mldistance(a, d, Distanmat[a][d]);
1090 Distanmat[d][a] = Distanmat[a][d];
1091 Distanmat[b][c] = mldistance(b, c, Distanmat[b][c]);
1092 Distanmat[c][b] = Distanmat[b][c];
1093 Distanmat[b][d] = mldistance(b, d, Distanmat[b][d]);
1094 Distanmat[d][b] = Distanmat[b][d];
1095 Distanmat[c][d] = mldistance(c, d, Distanmat[c][d]);
1096 Distanmat[d][c] = Distanmat[c][d];
1097 } /* distupdate */
1098
1099
1100 /* cleanup after ML analysis */
mlfinish()1101 void mlfinish()
1102 {
1103 if (Ctree != NULL)
1104 free_tree(Ctree, Numspc);
1105 free_ivector(bestrate);
1106 free_ivector(Alias);
1107 free_cmatrix(Seqpat);
1108 free_ivector(constpat);
1109 free_ivector(Weight);
1110 free_dmatrix(Distanmat);
1111 free_dvector(Eval);
1112 free_dmatrix(Evec);
1113 free_dmatrix(Ievc);
1114 free_dvector(Rates);
1115 free_dcube(ltprobr);
1116 free_dmatrix(iexp);
1117 } /* mlfinish */
1118
1119
1120 /******************************************************************************/
1121 /* tree output */
1122 /******************************************************************************/
1123
1124
1125 #define MAXOVER 50
1126 #define MAXLENG 30
1127 #define MAXCOLUMN 80
1128
1129
prbranch(Node * up,int depth,int m,int maxm,ivector umbrella,ivector column,FILE * outfp)1130 void prbranch(Node *up, int depth, int m, int maxm,
1131 ivector umbrella, ivector column, FILE *outfp)
1132 {
1133 int i, num, n, maxn, lim;
1134 Node *cp;
1135 char bch;
1136
1137 if ((int)((clockmode ? up->lengthc : up->length) * Proportion) >= MAXOVER) {
1138 column[depth] = MAXLENG;
1139 bch = '+';
1140 } else {
1141 column[depth] = (int)((clockmode ? up->lengthc : up->length) * Proportion) + 3;
1142 bch = '-';
1143 }
1144
1145 if (up->isop == NULL) { /* external branch */
1146 num = up->number + 1; /* offset */
1147 if (m == 1) umbrella[depth - 1] = TRUE;
1148 for (i = 0; i < depth; i++) {
1149 if (umbrella[i])
1150 fprintf(outfp, "%*c", column[i], ':');
1151 else
1152 fprintf(outfp, "%*c", column[i], ' ');
1153 }
1154 if (m == maxm)
1155 umbrella[depth - 1] = FALSE;
1156 for (i = 0, lim = column[depth] - 3; i < lim; i++)
1157 fputc(bch, outfp);
1158 fprintf(outfp, "-%d ", num);
1159
1160 fputid(outfp, up->number);
1161
1162
1163 fputc('\n', outfp);
1164 fputc(' ', outfp);
1165 return;
1166 }
1167
1168 num = up->number + 1 + Numspc; /* offset, internal branch */
1169 for (cp = up->isop, maxn = 0; cp != up; cp = cp->isop, maxn++)
1170 ;
1171 for (cp = up->isop, n = 1; cp != up; cp = cp->isop, n++) {
1172 prbranch(cp->kinp, depth + 1, n, maxn, umbrella, column, outfp);
1173 if (m == 1 && n == maxn / 2) umbrella[depth - 1] = TRUE;
1174 if (n != maxn) {
1175 for (i = 0; i < depth; i++) {
1176 if (umbrella[i])
1177 fprintf(outfp, "%*c", column[i], ':');
1178 else
1179 fprintf(outfp, "%*c", column[i], ' ');
1180 }
1181 if (n == maxn / 2) { /* internal branch */
1182 for (i = 0, lim = column[depth] - 3; i < lim; i++)
1183 fputc(bch, outfp);
1184 if (num < 10)
1185 fprintf(outfp, "--%d", num);
1186 else if (num < 100)
1187 fprintf(outfp, "-%2d", num);
1188 else
1189 fprintf(outfp, "%3d", num);
1190 } else {
1191 if (umbrella[depth])
1192 fprintf(outfp, "%*c", column[depth], ':');
1193 else
1194 fprintf(outfp, "%*c", column[depth], ' ');
1195 }
1196 fputc('\n', outfp);
1197 fputc(' ', outfp);
1198 }
1199 if (m == maxm) umbrella[depth - 1] = FALSE;
1200 }
1201 return;
1202 } /* prbranch */
1203
1204
getproportion(double * proportion,dvector distanvec,int numspc)1205 void getproportion(double *proportion, dvector distanvec, int numspc)
1206 {
1207 int i, maxpair;
1208 double maxdis;
1209
1210 maxpair = (numspc*(numspc-1))/2;
1211
1212 maxdis = 0.0;
1213 for (i = 0; i < maxpair; i++) {
1214 if (distanvec[i] > maxdis) {
1215 maxdis = distanvec[i];
1216 }
1217 }
1218 *proportion = (double) MAXCOLUMN / (maxdis * 3.0);
1219 if (*proportion > 1.0) *proportion = 1.0;
1220 } /* getproportion */
1221
1222
prtopology(FILE * outfp)1223 void prtopology(FILE *outfp)
1224 {
1225 int n, maxn, depth;
1226 ivector umbrella;
1227 ivector column;
1228 Node *cp, *rp;
1229
1230 getproportion(&Proportion, Distanvec, Numspc);
1231
1232 umbrella = new_ivector(Numspc);
1233 column = new_ivector(Numspc);
1234
1235 for (n = 0; n < Numspc; n++) {
1236 umbrella[n] = FALSE;
1237 column[n] = 3;
1238 }
1239 column[0] = 1;
1240
1241 fputc(' ', outfp);
1242
1243 /* original code: rp = Ctree->rootp */
1244 /* but we want to print the first group in the
1245 trichotomy as outgroup at the bottom! */
1246 rp = Ctree->rootp->isop;
1247
1248 for (maxn = 1, cp = rp->isop; cp != rp; cp = cp->isop, maxn++)
1249 ;
1250 depth = 1;
1251 n = 0;
1252
1253 cp = rp;
1254 do {
1255 cp = cp->isop;
1256 n++;
1257 prbranch(cp->kinp, depth, n, maxn, umbrella, column, outfp);
1258 if (cp != rp) fprintf(outfp, "%*c\n ", column[0], ':');
1259 } while (cp != rp);
1260
1261 free_ivector(umbrella);
1262 free_ivector(column);
1263 } /* prtopology */
1264
1265
1266 /* print unrooted tree file with branch lengths */
fputphylogeny(FILE * fp)1267 void fputphylogeny(FILE *fp)
1268 {
1269 Node *cp, *rp;
1270 int n;
1271
1272 cp = rp = Ctree->rootp;
1273 putc('(', fp);
1274 n = 1;
1275 do {
1276 cp = cp->isop->kinp;
1277 if (cp->isop == NULL) { /* external node */
1278 if (n > 60) {
1279 fprintf(fp, "\n");
1280 n = 2;
1281 }
1282 n += fputid(fp, cp->number);
1283 fprintf(fp, ":%.5f", ((clockmode ? cp->lengthc : cp->length))*0.01);
1284 n += 7;
1285 cp = cp->kinp;
1286 } else { /* internal node */
1287 if (cp->descen) {
1288 if (n > 60) {
1289 fprintf(fp, "\n");
1290 n = 1;
1291 }
1292 putc('(', fp);
1293 n++;
1294 } else {
1295 putc(')', fp);
1296 n++;
1297 if (n > 60) {
1298 fprintf(fp, "\n");
1299 n = 1;
1300 }
1301 /* internal label */
1302 if (cp->kinp->label != NULL) {
1303 fprintf(fp, "%s", cp->kinp->label);
1304 n += strlen(cp->kinp->label);
1305 }
1306 fprintf(fp, ":%.5f", ((clockmode ? cp->lengthc : cp->length))*0.01);
1307 n += 7;
1308 }
1309 }
1310 if (!cp->descen && !cp->isop->descen && cp != rp) {
1311 putc(',', fp); /* not last subtree */
1312 n++;
1313 }
1314 } while (cp != rp);
1315 fprintf(fp, ")");
1316 /* internal label */
1317 if (cp->label != NULL)
1318 fprintf(fp, "%s", cp->label);
1319 fprintf(fp, ";\n");
1320 } /* fputphylogeny */
1321
1322
resulttree(FILE * outfp,int usebranch)1323 void resulttree(FILE *outfp, int usebranch)
1324 {
1325 int n, ne, closeflag;
1326 Node *ep, *ip;
1327 double blen;
1328
1329 closeflag = FALSE;
1330
1331 if (clockmode) {
1332 fprintf(outfp, "\n branch length nc/c");
1333 fprintf(outfp, " branch length nc/c (= non-clock/clock)\n");
1334 } else {
1335 fprintf(outfp, "\n branch length S.E.");
1336 fprintf(outfp, " branch length S.E.\n");
1337 }
1338 for (n = 0; n < Numspc; n++) {
1339 ep = Ctree->ebrnchp[n];
1340 ne = ep->number;
1341 fputid10(outfp, ne);
1342 fputs(" ", outfp);
1343 fprintf(outfp, "%3d", ne + 1);
1344 blen = (clockmode ? ep->lengthc : ep->length);
1345 fprintf(outfp, "%9.5f", blen*0.01);
1346 if (blen < 5.0*MINARC || blen > 0.95*MAXARC) closeflag = TRUE;
1347 if (clockmode)
1348 fprintf(outfp, "%9.3f", (ep->length)/(ep->lengthc));
1349 else
1350 fprintf(outfp, "%9.5f", 0.01*sqrt(ep->kinp->varlen));
1351 if (n < Numibrnch) {
1352 ip = Ctree->ibrnchp[n];
1353 fprintf(outfp, "%8d", n + 1 + Numspc);
1354 blen = (clockmode ? ip->lengthc : ip->length);
1355 fprintf(outfp, "%9.5f", blen*0.01);
1356 if (blen < 5.0*MINARC || blen > 0.95*MAXARC) closeflag = TRUE;
1357 if (clockmode)
1358 fprintf(outfp, "%9.3f", (ip->length)/(ip->lengthc));
1359 else
1360 fprintf(outfp, "%9.5f", 0.01*sqrt(ip->kinp->varlen));
1361 fputc('\n', outfp);
1362 } else {
1363 if (n == Numspc - 3) {
1364 fputc('\n', outfp);
1365 } else if (n == Numspc - 2) {
1366 if (clockmode) {
1367 if (usebranch) {
1368 fprintf(outfp, " Branch lengths set by user (S.E.=0.0)!\n");
1369 } else {
1370 if (!Convergc)
1371 fprintf(outfp, " No convergence after %d iterations!\n", Numitc);
1372 else
1373 fprintf(outfp, " %d iterations until convergence\n", Numitc);
1374 }
1375 } else {
1376 if (usebranch) {
1377 fprintf(outfp, " Branch lengths set by user (S.E.=0.0)!\n");
1378 } else {
1379 if (!Converg)
1380 fprintf(outfp, " No convergence after %d iterations!\n", Numit);
1381 else
1382 fprintf(outfp, " %d iterations until convergence\n", Numit);
1383 }
1384 }
1385 } else if (n == Numspc - 1) {
1386 fprintf(outfp, " log L: %.2f\n", (clockmode ? Ctree->lklhdc : Ctree->lklhd));
1387 } else {
1388 fputc('\n', outfp);
1389 }
1390 }
1391 }
1392 if(closeflag)
1393 fprintf(outfp, "\nWARNING --- at least one branch length is close to an internal boundary!\n");
1394 } /* resulttree */
1395
1396
1397 /******************************************************************************/
1398 /* Neighbor-joining tree */
1399 /******************************************************************************/
1400
1401
1402 /* compute NJ tree and write to file */
njtree(FILE * fp)1403 void njtree(FILE *fp)
1404 {
1405 /* reserve memory for tree if necessary */
1406 if (mlmode != NJ_TREE) { /* no tree */
1407 if (Ctree != NULL)
1408 free_tree(Ctree, Numspc);
1409 Ctree = new_tree(Maxspc, Numptrn, Seqpat);
1410 Numbrnch = 2*Maxspc-3;
1411 Numibrnch = Maxspc-3;
1412 Numspc = Maxspc;
1413 mlmode = NJ_TREE;
1414 }
1415
1416 /* construct NJ tree from distance matrix */
1417 njdistantree(Ctree);
1418
1419 fputphylogeny(fp);
1420 } /* njtree */
1421
1422
1423 /* construct NJ tree from distance matrix */
njdistantree(Tree * tr)1424 void njdistantree(Tree *tr)
1425 {
1426 int i, j, otui=0, otuj=0, otuk, nsp2, cinode, step, restsp, k;
1427 double dij, bix, bjx, bkx, sij, smax, dnsp2;
1428 dvector r;
1429 dmatrix distan;
1430 Node **psotu, *cp, *ip, *jp, *kp;
1431
1432 distan = new_dmatrix(Maxspc,Maxspc);
1433 for (i = 0; i < Maxspc; i++)
1434 for (j = 0; j < Maxspc; j++)
1435 distan[i][j] = Distanmat[i][j];
1436
1437 nsp2 = Maxspc - 2;
1438 dnsp2 = 1.0 / nsp2;
1439
1440 r = new_dvector(Maxspc);
1441
1442 psotu = (Node **) calloc((size_t)Maxspc, sizeof(Node *));
1443 if (psotu == NULL) maerror("psotu in njdistantree");
1444
1445 /* external branches are start OTUs */
1446 for (i = 0; i < Maxspc; i++)
1447 psotu[i] = tr->ebrnchp[i]->kinp;
1448
1449 restsp = Maxspc;
1450 cinode = 0; /* counter for internal nodes */
1451
1452 for (step = 0; restsp > 3; step++) { /* NJ clustering steps */
1453
1454 for (i = 0; i < Maxspc; i++) {
1455 if (psotu[i] != NULL) {
1456 for (j = 0, r[i] = 0.0; j < Maxspc; j++)
1457 if (psotu[j] != NULL)
1458 r[i] += distan[i][j];
1459 }
1460 }
1461
1462 smax = -1.0;
1463 for (i = 0; i < Maxspc-1; i++) {
1464 if (psotu[i] != NULL) {
1465
1466 for (j = i+1; j < Maxspc; j++) {
1467 if (psotu[j] != NULL)
1468 {
1469 sij = ( r[i] + r[j] ) * dnsp2 - distan[i][j];
1470
1471 if (sij > smax) {
1472 smax = sij;
1473 otui = i;
1474 otuj = j;
1475 }
1476 }
1477 }
1478 }
1479 }
1480
1481 /* new pair: otui and otuj */
1482
1483 dij = distan[otui][otuj];
1484 bix = (dij + r[otui]/nsp2 - r[otuj]/nsp2) * 0.5;
1485 bjx = dij - bix;
1486
1487 cp = tr->ibrnchp[cinode];
1488
1489 ip = psotu[otui];
1490 jp = psotu[otuj];
1491 cp->isop = ip;
1492 ip->isop = jp;
1493 jp->isop = cp;
1494 ip->length = bix;
1495 jp->length = bjx;
1496 ip->kinp->length = ip->length;
1497 jp->kinp->length = jp->length;
1498
1499 cp = cp->kinp;
1500
1501 for (k = 0; k < Maxspc; k++)
1502 {
1503 if (psotu[k] != NULL && k != otui && k != otuj)
1504 {
1505 dij = (distan[otui][k] + distan[otuj][k] - distan[otui][otuj]) * 0.5;
1506 distan[otui][k] = dij;
1507 distan[k][otui] = dij;
1508 }
1509 }
1510 distan[otui][otui] = 0.0;
1511
1512 psotu[otui] = cp;
1513 psotu[otuj] = NULL;
1514
1515 cinode++;
1516
1517 restsp--;
1518 nsp2--;
1519 dnsp2 = 1.0 / nsp2;
1520 }
1521
1522 otui = otuj = otuk = -1;
1523 for (i = 0; i < Maxspc; i++)
1524 {
1525 if (psotu[i] != NULL) {
1526 if (otui == -1) otui = i;
1527 else if (otuj == -1) otuj = i;
1528 else otuk = i;
1529 }
1530 }
1531 bix = (distan[otui][otuj] + distan[otui][otuk] - distan[otuj][otuk]) * 0.5;
1532 bjx = distan[otui][otuj] - bix;
1533 bkx = distan[otui][otuk] - bix;
1534 ip = psotu[otui];
1535 jp = psotu[otuj];
1536 kp = psotu[otuk];
1537 ip->isop = jp;
1538 jp->isop = kp;
1539 kp->isop = ip;
1540 ip->length = bix;
1541 jp->length = bjx;
1542 kp->length = bkx;
1543 ip->kinp->length = ip->length;
1544 jp->kinp->length = jp->length;
1545 kp->kinp->length = kp->length;
1546
1547 tr->rootp = kp;
1548
1549 free_dvector(r);
1550 free_dmatrix(distan);
1551 free((Node *) psotu);
1552 } /* njdistantree */
1553
1554 /******************************************************************************/
1555 /* find best assignment of rate categories */
1556 /******************************************************************************/
1557
1558 /* find best assignment of rate categories */
1559 #if 0
1560 void findbestratecombination(Tree *ctree, /* Ctree */
1561 cmatrix Seqpat,
1562 dvector Freqtpm,
1563 double fracinv,
1564 ivector constpat,
1565 ivector bestrate,
1566 int bestratefound,
1567 int ncats) /* numcats */
1568 #endif
findbestratecombination()1569 void findbestratecombination()
1570 {
1571 int k, u;
1572 double bestvalue, fv2;
1573 dvector catprob;
1574 dmatrix cdl;
1575
1576 /* cdl = ctree->condlkl; */
1577 cdl = Ctree->condlkl;
1578
1579
1580 /* catprob = new_dvector(ncats+1); */
1581 catprob = new_dvector(numcats+1);
1582 /* fv2 = (1.0-fracinv)/(double) ncats; */
1583 fv2 = (1.0-fracinv)/(double) numcats;
1584
1585 for (k = 0; k < Numptrn; k++) {
1586
1587 /* zero rate */
1588 if (constpat[k] == TRUE)
1589 catprob[0] = fracinv*Freqtpm[(int) Seqpat[0][k]];
1590 else
1591 catprob[0] = 0.0;
1592
1593 /* non-zero-rates */
1594 /* for (u = 1; u < ncats+1; u++) { */
1595 for (u = 1; u < numcats+1; u++) {
1596 catprob[u] = fv2*cdl[u-1][k];
1597 }
1598
1599 /* find best */
1600 bestvalue = catprob[0];
1601 bestrate[k] = 0;
1602 /* for (u = 1; u < ncats+1; u++) */
1603 for (u = 1; u < numcats+1; u++)
1604 if (catprob[u] >= bestvalue) {
1605 bestvalue = catprob[u];
1606 bestrate[k] = u;
1607 }
1608 } /* for all column patterns */
1609
1610 free_dvector(catprob);
1611 bestratefound = 1;
1612 } /* findbestratecombination */
1613
1614
1615 /* print best assignment of rate categories */
printbestratecombination(FILE * fp)1616 void printbestratecombination(FILE *fp)
1617 {
1618 int s, k;
1619
1620 #ifndef USE_WINDOWS
1621 for (s = 0; s < Maxsite; s++) {
1622 #else
1623 for (s = 0; s < alimaxsite; s++) {
1624 #endif
1625 k = Alias[s];
1626 fprintf(fp, "%2d", bestrate[k]);
1627 if ((s+1) % 30 == 0)
1628 fprintf(fp, "\n");
1629 else if ((s+1) % 10 == 0)
1630 fprintf(fp, " ");
1631 }
1632 if (s % 70 != 0)
1633 fprintf(fp, "\n");
1634 } /* printbestratecombination */
1635
1636
1637
1638
1639 #if 0
1640 # define MLUNIFORMRATE 0
1641 # define MLGAMMARATE 1
1642 # define MLTWORATE 2
1643 # define MLMIXEDRATE 3
1644 #endif
1645 void printbestratecombinationtofile(FILE *fp, int rhettype)
1646 {
1647 int s;
1648 char name[11];
1649
1650 switch(rhettype) {
1651 case GAMMARATE:
1652 sprintf(name, "%drates", numcats);
1653 break;;
1654 case MIXEDRATE:
1655 sprintf(name, "1inv+%d", numcats);
1656 break;;
1657 case TWORATE:
1658 sprintf(name, "1inv+1var");
1659 break;;
1660 default: /* UNIFORMRATE */
1661 sprintf(name, "uniform");
1662 }
1663 #ifndef USE_WINDOWS
1664 fprintf(fp, " 1 %d %d \n", Maxsite, numcats);
1665 fprintf(fp, "%-10s", name);
1666 for (s = 0; s < Maxsite; s++) {
1667 #else
1668 fprintf(fp, " 1 %d %d \n", alimaxsite, numcats);
1669 fprintf(fp, "%-10s", name);
1670 for (s = 0; s < alimaxsite; s++) {
1671 #endif
1672 if (rhettype != UNIFORMRATE)
1673 fprintf(fp, " %.5f", Rates[bestrate[Alias[s]] - 1]);
1674 else
1675 fprintf(fp, " 1.0");
1676 }
1677 fprintf(fp, "\n");
1678
1679 } /* printbestratecombination */
1680 #if 0
1681 # undef MLUNIFORMRATE
1682 # undef MLGAMMARATE
1683 # undef MLTWORATE
1684 # undef MLMIXEDRATE
1685 #endif
1686
1687
1688 /******************************************************************************/
1689 /* computation of clocklike branch lengths */
1690 /******************************************************************************/
1691
1692 /* checks wether e is a valid edge specification */
1693 int checkedge(int e)
1694 {
1695 /* there are Numspc external branches:
1696 0 - Numspc-1
1697 there are Numibrnch internal branches:
1698 Numspc - Numspc+Numibrnch-1
1699 */
1700
1701 if (e < 0) return FALSE;
1702 if (e < Numspc+Numibrnch) return TRUE;
1703 else return FALSE;
1704 } /* checkedge */
1705
1706 /* print topology of subtree */
1707 void fputsubstree(FILE *fp, Node *ip)
1708 {
1709 Node *cp;
1710
1711 if (ip->isop == NULL) { /* terminal nodes */
1712 numtc += fputid(fp, ip->number);
1713 } else {
1714 cp = ip;
1715 fprintf(fp, "(");
1716 numtc += 1;
1717 do {
1718 cp = cp->isop->kinp;
1719 if (cp->isop == NULL) { /* external node */
1720 numtc += fputid(fp, cp->number);
1721 fprintf(fp, ":%.5f", (cp->lengthc)*0.01);
1722 numtc += 7;
1723 cp = cp->kinp;
1724 } else { /* internal node */
1725 if (cp->height > 0.0) {
1726 fprintf(fp, "(");
1727 numtc += 1;
1728 } else if (cp->height < 0.0) {
1729 fprintf(fp, ")");
1730 numtc += 1;
1731 if (numtc > 60) {
1732 fprintf(fp, "\n");
1733 numtc = 1;
1734 }
1735 /* internal label */
1736 if (cp->kinp->label != NULL) {
1737 fprintf(fp, "%s", cp->kinp->label);
1738 numtc += strlen(cp->kinp->label);
1739 }
1740 if (numtc > 60) {
1741 fprintf(fp, "\n");
1742 numtc = 1;
1743 }
1744 fprintf(fp, ":%.5f", (cp->lengthc)*0.01);
1745 numtc += 6;
1746 if (numtc > 60) {
1747 fprintf(fp, "\n");
1748 numtc = 1;
1749 }
1750 }
1751 }
1752 if (cp->height <= 0.0 && cp->isop->height <= 0.0 &&
1753 cp->isop != ip) {
1754 putc(',', fp); /* not last subtree */
1755 numtc += 1;
1756 if (numtc > 60) {
1757 fprintf(fp, "\n");
1758 numtc = 1;
1759 }
1760 }
1761 } while (cp->isop != ip);
1762 fprintf(fp, ")");
1763 numtc += 1;
1764 }
1765 if (numtc > 60) {
1766 fprintf(fp, "\n");
1767 numtc = 1;
1768 }
1769
1770 } /* fputsubstree */
1771
1772 /* print rooted tree file */
1773 void fputrooted(FILE *fp, int e)
1774 {
1775 Node *rootbr;
1776
1777 /* to be called only after clocklike branch
1778 lengths have been computed */
1779
1780 /* pointer to root branch */
1781 if (e < Numspc) rootbr = Ctree->ebrnchp[e];
1782 else rootbr = Ctree->ibrnchp[e - Numspc];
1783
1784 fprintf(fp, "(");
1785 numtc = 2;
1786 fputsubstree(fp, rootbr);
1787 /* internal label */
1788 if (rootbr->label != NULL) {
1789 fprintf(fp, "%s", rootbr->label);
1790 numtc += strlen(rootbr->label);
1791 }
1792 if (numtc > 60) {
1793 fprintf(fp, "\n");
1794 numtc = 1;
1795 }
1796 fprintf(fp, ":%.5f,", (hroot - rootbr->height)*0.01);
1797 numtc += 7;
1798 if (numtc > 60) {
1799 fprintf(fp, "\n");
1800 numtc = 1;
1801 }
1802 fputsubstree(fp, rootbr->kinp);
1803 /* internal label */
1804 if (rootbr->kinp->label != NULL) {
1805 fprintf(fp, "%s", rootbr->kinp->label);
1806 numtc += strlen(rootbr->kinp->label);
1807 }
1808 if (numtc > 60) {
1809 fprintf(fp, "\n");
1810 numtc = 1;
1811 }
1812 fprintf(fp, ":%.5f);\n", (hroot - rootbr->kinp->height)*0.01);
1813 } /* fputrooted */
1814
1815 /* finds heights in subtree */
1816 void findheights(Node *ip)
1817 {
1818 Node *cp, *rp;
1819
1820 if (ip->isop != NULL) { /* forget terminal nodes */
1821
1822 cp = ip;
1823
1824 /* initialise node */
1825 cp->height = 1.0; /* up */
1826 rp = cp;
1827 while (rp->isop != cp) {
1828 rp = rp->isop;
1829 rp->height = -1.0; /* down */
1830 }
1831
1832 do {
1833 cp = cp->isop->kinp;
1834 if (cp->isop == NULL) { /* external node */
1835 cp = cp->kinp;
1836 } else { /* internal node */
1837 if (cp->height == 0.0) { /* node not yet visited */
1838 cp->height = 1.0; /* up */
1839 rp = cp;
1840 while (rp->isop != cp) {
1841 rp = rp->isop;
1842 rp->height = -1.0; /* down */
1843 }
1844 } else if (cp->kinp->height == 1.0) {
1845 /* cp->kinp is next height pointer */
1846 heights[Numhts] = cp->kinp;
1847 Numhts++;
1848 }
1849 }
1850 } while (cp->isop != ip);
1851 /* ip is last height pointer */
1852 heights[Numhts] = ip;
1853 Numhts++;
1854 }
1855 } /* findheights */
1856
1857
1858 /* initialise clocklike branch lengths (with root on edge e) */
1859 void initclock(int e)
1860 {
1861 int n, h, count;
1862 Node *cp, *rp;
1863 double sum, minh, aveh, len;
1864
1865 /* be sure to have a Ctree in memory and
1866 to have pairwise distances computed */
1867
1868 clockmode = 1; /* clocklike branch lengths */
1869
1870 /* least square estimate for branch length */
1871 changedistan(Distanmat, Distanvec, Numspc);
1872 lslength(Ctree, Distanvec, Numspc, Numibrnch, Brnlength);
1873
1874 /* pointer to root branch */
1875 if (e < Numspc) rootbr = Ctree->ebrnchp[e];
1876 else rootbr = Ctree->ibrnchp[e - Numspc];
1877
1878 /* clear all heights */
1879 for (n = 0; n < Numspc; n++) {
1880 Ctree->ebrnchp[n]->height = 0.0;
1881 Ctree->ebrnchp[n]->kinp->height = 0.0;
1882 Ctree->ebrnchp[n]->varheight = 0.0;
1883 Ctree->ebrnchp[n]->kinp->varheight = 0.0;
1884 if (n < Numibrnch) {
1885 Ctree->ibrnchp[n]->height = 0.0;
1886 Ctree->ibrnchp[n]->kinp->height = 0.0;
1887 Ctree->ibrnchp[n]->varheight = 0.0;
1888 Ctree->ibrnchp[n]->kinp->varheight = 0.0;
1889 }
1890 }
1891
1892 /* collect pointers to height nodes */
1893 Numhts = 0;
1894 findheights(rootbr); /* one side */
1895 findheights(rootbr->kinp); /* other side */
1896
1897 /* assign preliminary approximate heights and
1898 corresponding branch lengths */
1899 for (h = 0; h < Numhts; h++) {
1900
1901 cp = rp = heights[h];
1902 sum = 0;
1903 count = 0;
1904 minh = 0.0;
1905 while (rp->isop != cp) {
1906 count++;
1907 rp = rp->isop;
1908 sum += rp->lengthc + rp->kinp->height;
1909 if (rp->kinp->height > minh) minh = rp->kinp->height;
1910 }
1911 aveh = sum / (double) count;
1912 if (aveh < minh + MINARC) aveh = minh + MINARC;
1913 cp->height = aveh;
1914 rp = cp;
1915 while (rp->isop != cp) {
1916 rp = rp->isop;
1917 len = aveh - rp->kinp->height;
1918 rp->kinp->lengthc = len;
1919 rp->lengthc = len;
1920 }
1921
1922 }
1923 if (rootbr->height > rootbr->kinp->height) minh = rootbr->height;
1924 else minh = rootbr->kinp->height;
1925 aveh = 0.5*(rootbr->lengthc + rootbr->height + rootbr->kinp->height);
1926 if (aveh < minh + MINARC) aveh = minh + MINARC;
1927 hroot = aveh;
1928 maxhroot = RMHROOT*hroot; /* maximal possible hroot */
1929 len = (hroot - rootbr->height) + (hroot - rootbr->kinp->height);
1930 rootbr->lengthc = len;
1931 rootbr->kinp->lengthc = len;
1932 } /* initclock */
1933
1934 /* approximate likelihood under the constaining assumption of
1935 clocklike branch lengths (with root on edge e) */
1936 double clock_alklhd(int e)
1937 {
1938 initclock(e);
1939 Ctree->lklhdc = treelkl(Ctree);
1940
1941 return Ctree->lklhdc;
1942 } /* clock_alklhd */
1943
1944 /* log-likelihood given height ht at node pointed to by chep */
1945 double heightlkl(double ht)
1946 {
1947 Node *rp;
1948 double len;
1949
1950 /* adjust branch lengths */
1951 chep->height = ht;
1952 /* descendent branches */
1953 rp = chep;
1954 while (rp->isop != chep) {
1955 rp = rp->isop;
1956 len = chep->height - rp->kinp->height;
1957 rp->kinp->lengthc = len;
1958 rp->lengthc = len;
1959 }
1960 /* upward branch */
1961 if (chep == rootbr || chep->kinp == rootbr) {
1962 len = (hroot - chep->height) + (hroot - chep->kinp->height);
1963 chep->lengthc = len;
1964 chep->kinp->lengthc = len;
1965 } else {
1966 rp = chep->kinp;
1967 while (rp->isop->height <= 0.0)
1968 rp = rp->isop;
1969 chep->lengthc = rp->isop->height - chep->height;
1970 chep->kinp->lengthc = rp->isop->height - chep->height;
1971 }
1972
1973 /* compute likelihood */
1974 Ctree->lklhdc = treelkl(Ctree);
1975
1976 return -(Ctree->lklhdc); /* we use a minimizing procedure */
1977 } /* heightlkl */
1978
1979 /* optimize current height */
1980 void optheight(void)
1981 {
1982 double he, fx, f2x, minh, maxh, len;
1983 Node *rp;
1984
1985 /* current height */
1986 he = chep->height;
1987
1988 /* minimum */
1989 minh = 0.0;
1990 rp = chep;
1991 while (rp->isop != chep) {
1992 rp = rp->isop;
1993 if (rp->kinp->height > minh)
1994 minh = rp->kinp->height;
1995 }
1996 minh += MINARC;
1997
1998 /* maximum */
1999 if (chep == rootbr || chep->kinp == rootbr) {
2000 maxh = hroot;
2001 } else {
2002 rp = chep->kinp;
2003 while (rp->isop->height <= 0.0)
2004 rp = rp->isop;
2005 maxh = rp->isop->height;
2006 }
2007 maxh -= MINARC;
2008
2009 /* check borders for height */
2010 if (he < minh) he = minh;
2011 if (he > maxh) he = maxh;
2012
2013 /* optimization */
2014 if (!(he == minh && he == maxh))
2015 he = onedimenmin(minh, he, maxh, heightlkl, EPSILON_HEIGHTS, &fx, &f2x);
2016
2017 /* variance of height */
2018 f2x = fabs(f2x);
2019 if (1.0/(maxhroot*maxhroot) < f2x)
2020 chep->varheight = 1.0/f2x;
2021 else
2022 chep->varheight = maxhroot*maxhroot;
2023
2024 /* adjust branch lengths */
2025 chep->height = he;
2026 /* descendent branches */
2027 rp = chep;
2028 while (rp->isop != chep) {
2029 rp = rp->isop;
2030 len = chep->height - rp->kinp->height;
2031 rp->kinp->lengthc = len;
2032 rp->lengthc = len;
2033 }
2034 /* upward branch */
2035 if (chep == rootbr || chep->kinp == rootbr) {
2036 len = (hroot - chep->height) + (hroot - chep->kinp->height);
2037 chep->lengthc = len;
2038 chep->kinp->lengthc = len;
2039 } else {
2040 rp = chep->kinp;
2041 while (rp->isop->height <= 0.0)
2042 rp = rp->isop;
2043 chep->lengthc = rp->isop->height - chep->height;
2044 chep->kinp->lengthc = rp->isop->height - chep->height;
2045 }
2046 } /* optheight */
2047
2048 /* log-likelihood given height ht at root */
2049 double rheightlkl(double ht)
2050 {
2051 double len;
2052
2053 /* adjust branch lengths */
2054 hroot = ht;
2055 len = (hroot - rootbr->height) + (hroot - rootbr->kinp->height);
2056 rootbr->lengthc = len;
2057 rootbr->kinp->lengthc = len;
2058
2059 /* compute likelihood */
2060 Ctree->lklhdc = treelkl(Ctree);
2061
2062 return -(Ctree->lklhdc); /* we use a minimizing procedure */
2063 } /* rheightlkl */
2064
2065 /* optimize height of root */
2066 void optrheight(void)
2067 {
2068 double he, fx, f2x, minh, len;
2069
2070 /* current height */
2071 he = hroot;
2072
2073 /* minimum */
2074 if (rootbr->height > rootbr->kinp->height)
2075 minh = rootbr->height;
2076 else
2077 minh = rootbr->kinp->height;
2078 minh += MINARC;
2079
2080 /* check borders for height */
2081 if (he < minh) he = minh;
2082 if (he > maxhroot) he = maxhroot;
2083
2084 /* optimization */
2085 he = onedimenmin(minh, he, maxhroot, rheightlkl, EPSILON_HEIGHTS, &fx, &f2x);
2086
2087 /* variance of height of root */
2088 f2x = fabs(f2x);
2089 if (1.0/(maxhroot*maxhroot) < f2x)
2090 varhroot = 1.0/f2x;
2091 else
2092 varhroot = maxhroot*maxhroot;
2093
2094 /* adjust branch lengths */
2095 hroot = he;
2096 len = (hroot - rootbr->height) + (hroot - rootbr->kinp->height);
2097 rootbr->lengthc = len;
2098 rootbr->kinp->lengthc = len;
2099 } /* optrheight */
2100
2101 /* exact likelihood under the constaining assumption of
2102 clocklike branch lengths (with root on edge e) */
2103 double clock_lklhd(int e)
2104 {
2105 int h, nconv;
2106 double old;
2107
2108 Numitc = 0;
2109 Convergc = FALSE;
2110
2111 initclock(e);
2112
2113 do {
2114
2115 Numitc++;
2116 nconv = 0;
2117
2118 /* optimize height of root */
2119 old = hroot;
2120 optrheight();
2121 if (fabs(old - hroot) < EPSILON_HEIGHTS) nconv++;
2122
2123 /* optimize height of nodes */
2124 for (h = Numhts-1; h >= 0; h--) {
2125
2126 /* pointer chep to current height node */
2127 chep = heights[h];
2128
2129 /* store old value */
2130 old = chep->height;
2131
2132 /* find better height */
2133 optheight();
2134
2135 /* converged ? */
2136 if (fabs(old - chep->height) < EPSILON_HEIGHTS) nconv++;
2137 }
2138
2139 if (nconv == Numhts+1) Convergc = TRUE;
2140
2141 } while (Numitc < MAXIT && !Convergc);
2142
2143 /* compute final likelihood */
2144 Ctree->lklhdc = treelkl(Ctree);
2145
2146 return Ctree->lklhdc;
2147 } /* clock_lklhd */
2148
2149 /* find out the edge containing the root */
2150 int findrootedge()
2151 {
2152 int e, ebest;
2153 double logbest, logtest;
2154
2155 /* compute the likelihood for all edges and take the edge with
2156 best likelihood (using approximate ML) */
2157
2158 ebest = 0;
2159 logbest = clock_alklhd(0);
2160 numbestroot = 1;
2161 for (e = 1; e < Numspc+Numibrnch; e++) {
2162 logtest = clock_alklhd(e);
2163 if (logtest > logbest) {
2164 ebest = e;
2165 logbest = logtest;
2166 numbestroot = 1;
2167 } else if (logtest == logbest) {
2168 numbestroot++;
2169 }
2170 }
2171
2172 return ebest;
2173 } /* findrootedge */
2174
2175 /* show heights and corresponding standard errors */
2176 void resultheights(FILE *fp)
2177 {
2178 int h, num;
2179 Node *cp;
2180
2181 fprintf(fp, " height S.E. of node common to branches\n");
2182 for (h = 0; h < Numhts; h++) {
2183 fprintf(fp, "%.5f %.5f ", (heights[h]->height)*0.01,
2184 sqrt(heights[h]->varheight)*0.01);
2185 cp = heights[h];
2186 do {
2187 num = (cp->number) + 1;
2188 if (cp->kinp->isop != NULL) num += Numspc; /* internal branch */
2189 fprintf(fp, "%d ", num);
2190 cp = cp->isop;
2191 } while (cp != heights[h]);
2192 fprintf(fp, "\n");
2193
2194 }
2195 fprintf(fp, "%.5f %.5f of root at branch %d\n",
2196 hroot*0.01, sqrt(varhroot)*0.01, locroot+1);
2197 } /* resultheights */
2198
2199