Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

source: tags/initial/GDE/MOLPHY/altree.c

Visit:

Last change on this file was 2, checked in by oldcode, 24 years ago
Initial revision
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 20.4 KB

Line
1	/*
2	* altree.c Adachi, J. 1995.09.22
3	* Copyright (C) 1992-1995 J. Adachi & M. Hasegawa. All rights reserved.
4	*/
5
6	#include "protml.h"
7
8	Tree *
9	new_atree(maxspc, maxibrnch, numptrn, seqconint)
10	int maxspc, maxibrnch, numptrn;
11	imatrix seqconint;
12	{
13	int n, i;
14	Tree *tr;
15	Node dp, up;
16
17	tr = (Tree *) malloc(sizeof(Tree));
18	if (tr == NULL) maerror("tr in new_atree().");
19	tr->ebrnchp = (Node *) malloc((unsigned)maxspc sizeof(Node *));
20	if (tr->ebrnchp == NULL) maerror("ebrnchp in new_atree().");
21	tr->ibrnchp = (Node *) malloc((unsigned)maxibrnch sizeof(Node *));
22	if (tr->ibrnchp == NULL) maerror("ibrnchp in new_atree().");
23	for (n = 0; n < maxspc; n++) {
24	dp = (Node *) malloc(sizeof(Node));
25	if (dp == NULL) maerror("dp in new_atree().");
26	up = (Node *) malloc(sizeof(Node));
27	if (up == NULL) maerror("up in new_atree().");
28	dp->isop = NULL;
29	up->isop = NULL;
30	dp->kinp = up;
31	up->kinp = dp;
32	dp->descen = TRUE;
33	up->descen = FALSE;
34	dp->num = n;
35	up->num = n;
36	dp->length = 0.0;
37	up->length = 0.0;
38	dp->lklhdl = 0.0;
39	up->lklhdl = 0.0;
40	dp->paths = new_ivector(maxspc);
41	up->paths = dp->paths;
42	for (i = 0; i < maxspc; i++) dp->paths[i] = 0;
43	dp->paths[n] = 1;
44	dp->eprob = seqconint[n];
45	dp->iprob = NULL;
46	up->eprob = NULL;
47	up->iprob = new_dmatrix(numptrn, Tpmradix);
48	tr->ebrnchp[n] = dp;
49	}
50	tr->rootp = NULL;
51	return tr;
52	} /_ new_atree /
53
54
55	Node *
56	new_dnode()
57	{
58	Node *dp;
59	int i;
60
61	dp = (Node *) malloc(sizeof(Node));
62	if (dp == NULL) maerror("dp in new_dnode().");
63	dp->isop = NULL;
64	dp->kinp = NULL;
65	dp->descen = TRUE;
66	dp->num = Maxbrnch;
67	dp->length = 0.0;
68	dp->lklhdl = 0.0;
69	dp->paths = new_ivector(Maxspc);
70	for (i = 0; i < Maxspc; i++) dp->paths[i] = 0;
71	dp->eprob = NULL;
72	dp->iprob = new_dmatrix(Numptrn, Tpmradix);
73	return dp;
74	} /_ new_dnode /
75
76
77	Node *
78	new_anode()
79	{
80	Node *ap;
81
82	ap = (Node *) malloc(sizeof(Node));
83	if (ap == NULL) maerror("ap in new_anode().");
84	ap->isop = NULL;
85	ap->kinp = NULL;
86	ap->descen = FALSE;
87	ap->num = Maxbrnch;
88	ap->length = 0.0;
89	ap->lklhdl = 0.0;
90	ap->paths = NULL;
91	ap->eprob = NULL;
92	ap->iprob = new_dmatrix(Numptrn, Tpmradix);
93	return ap;
94	} /_ new_anode /
95
96
97	Node ***
98	new_nodematrix(nrow, ncol)
99	int nrow;
100	int ncol;
101	/* memory allocate a node matrix */
102	{
103	int i;
104	Node ***m;
105
106	m = (Node **) malloc((unsigned)nrow sizeof(Node **));
107	if (m == NULL) maerror("1 in nodematrix().");
108	m = (Node ) malloc((unsigned)(nrow ncol) * sizeof(Node *));
109	if (*m == NULL) maerror("2 in nodematrix().");
110	for (i = 1; i < nrow; i++) m[i] = m[i-1] + ncol;
111	return m;
112	}
113
114
115	void
116	free_nodematrix(m)
117	Node ***m;
118	{
119	free((char ) m);
120	free((char *) m);
121	}
122
123
124	#if 1
125	void
126	aproxlkl(tr)
127	Tree *tr;
128	{
129	int i, k;
130	Node *cp;
131	double sumlk, lklhd;
132	dmatrix prob1, prob2, prob3;
133
134	cp = tr->rootp;
135	lklhd = 0.0;
136	if (cp->isop->isop->isop == cp) { /* binary tree */
137	prob1 = cp->iprob;
138	prob2 = cp->isop->iprob;
139	prob3 = cp->isop->isop->iprob;
140	for (k = 0; k < Numptrn; k++) {
141	sumlk = 0.0;
142	for (i = 0; i < Tpmradix; i++)
143	sumlk += Freqtpm[i] * prob1[k][i] * prob2[k][i] * prob3[k][i];
144	/* printf("%3d%23.20f%10.3f%4d\n", k,sumlk,log(sumlk),Weight[k]); */
145	lklhd += log(sumlk) * Weight[k];
146	}
147	} else { /* multiway tree */
148	prodpart(cp->isop);
149	prob1 = cp->iprob;
150	prob2 = cp->isop->iprob;
151	for (k = 0; k < Numptrn; k++) {
152	sumlk = 0.0;
153	for (i = 0; i < Tpmradix; i++)
154	sumlk += Freqtpm[i] * prob1[k][i] * prob2[k][i];
155	/* printf("%3d%23.20f%10.3f%4d\n", k,sumlk,log(sumlk),Weight[k]); */
156	lklhd += log(sumlk) * Weight[k];
157	}
158	}
159	tr->aproxl = lklhd;
160	tr->lklhd = lklhd; /* ?!?!?! */
161	return;
162	} /_ aproxlkl /
163
164	#else
165	void
166	aproxlkl(tr)
167	Tree *tr;
168	{
169	int i, k;
170	Node cp, rp;
171	double sumlk, lklhd;
172	dmatrix prob1, prob2, prob3;
173
174	cp = rp = tr->rootp;
175	do {
176	cp = cp->isop->kinp;
177	if (cp->isop == NULL) { /* external node */
178	cp = cp->kinp; /* not descen */
179	partelkl(cp);
180	} else { /* internal node */
181	if (!cp->descen) {
182	prodpart(cp->kinp->isop);
183	partilkl(cp);
184	}
185	}
186	} while (cp != rp);
187
188	lklhd = 0.0;
189	if (cp->isop->isop->isop == cp) { /* binary tree */
190	prob1 = cp->iprob;
191	prob2 = cp->isop->iprob;
192	prob3 = cp->isop->isop->iprob;
193	for (k = 0; k < Numptrn; k++) {
194	sumlk = 0.0;
195	for (i = 0; i < Tpmradix; i++)
196	sumlk += Freqtpm[i] * prob1[k][i] * prob2[k][i] * prob3[k][i];
197	/* printf("%3d%23.20f%10.3f%4d\n", k,sumlk,log(sumlk),Weight[k]); */
198	lklhd += log(sumlk) * Weight[k];
199	}
200	} else { /* multiway tree */
201	prodpart(cp->isop);
202	prob1 = cp->iprob;
203	prob2 = cp->isop->iprob;
204	for (k = 0; k < Numptrn; k++) {
205	sumlk = 0.0;
206	for (i = 0; i < Tpmradix; i++)
207	sumlk += Freqtpm[i] * prob1[k][i] * prob2[k][i];
208	/* printf("%3d%23.20f%10.3f%4d\n", k,sumlk,log(sumlk),Weight[k]); */
209	lklhd += log(sumlk) * Weight[k];
210	}
211	}
212	tr->lklhd = lklhd;
213	return;
214	} /_ aproxlkl /
215	#endif
216
217
218	void
219	praproxlkl(tr)
220	Tree *tr;
221	{
222	printf("%.1f\t%.1f\t%d\t",tr->lklhd,tr->tbldis,Cnotree+1); /* offset */
223	} /_ praproxlkl /
224
225
226	void
227	aproxtree(tr, ntr)
228	Tree *tr;
229	int ntr;
230	{
231	int index;
232	double lkl, tbl;
233	Infoaltree cp, bp, *op;
234	static num = 0;
235
236	if (Verbs_optn) {
237	if (ntr % Numverbs == 0) {
238	fprintf(stderr, " %d", ntr+1);
239	#if __STDC__ && DIFFTIME
240	Ct1 = time(NULL);
241	fprintf(stderr, "(%.0fs)", difftime(Ct1, Ct0));
242	#endif
243	}
244	}
245	pathing(tr);
246	slslength(tr, Distanmat, Numspc);
247	tbl = tr->tbldis;
248	index = (int)((tbl - Mintbldm) * Tblcoef);
249	if (index < 0) {
250	Tblunder++;
251	} else if (index < NUMTBLBIN) {
252	Tblbin[index]++;
253	} else {
254	Tblover++;
255	}
256
257	if (!Mevol_optn) {
258	if (tbl > Basetbldm) return;
259	initpartlkl(tr);
260	aproxlkl(tr);
261	}
262	if (Aprox_optn) {
263	praproxlkl(tr);
264	putctopology(tr);
265	}
266	if (Mevol_optn) { /* Minimum Evolution */
267	lkl = -tbl;
268	} else { /* approximate ln likelihood */
269	lkl = tr->aproxl;
270	}
271	if (num < Maxaltree) {
272	Numaltree = num + 1;
273	Infoaltrees[num].lklaprox = lkl;
274	Infoaltrees[num].tbl = tbl;
275	strctree(tr, Infoaltrees[num].ltplgy);
276	op = &Infoaltrees[num];
277	if (Info_optn)
278	printf("%.1f\t%d\t%d\t%s\n", op->lklaprox, ntr,num, op->ltplgy);
279	for (bp = &Atail, cp = bp->up; lkl > cp->lklaprox; bp = cp, cp = cp->up)
280	;
281	op->up = cp;
282	bp->up = op;
283	num++;
284	} else if (lkl > Atail.up->lklaprox) {
285	op = Atail.up;
286	Atail.up = op->up;
287	op->lklaprox = lkl;
288	op->tbl = tbl;
289	strctree(tr, op->ltplgy);
290	if (Info_optn)
291	printf("%.1f\t%d\t%d\t%s\n", op->lklaprox, ntr,num, op->ltplgy);
292	for (bp = &Atail, cp = bp->up; lkl > cp->lklaprox; bp = cp, cp = cp->up)
293	;
294	op->up = cp;
295	bp->up = op;
296	}
297	/* if (ntr % 10000 == 0) putctopology(tr); */
298	} /_ aproxtree /
299
300
301	void autoconstruction();
302
303	void
304	wedge(tr, onode, poolnode, addposition, poolorder, op)
305	Tree *tr;
306	int onode;
307	Node poolnode, addposition;
308	ivector poolorder;
309	Node *op;
310	{
311	Node cp, dp;
312	int i, istart;
313
314	cp = poolnode[onode];
315	if (addposition[onode] != NULL) { /* internal node */
316	if (cp->isop != NULL) { /* binary tree */
317	dp = op->kinp;
318	op->kinp = cp->isop;
319	cp->isop->kinp = op;
320	dp->kinp = cp->isop->isop;
321	cp->isop->isop->kinp = dp;
322	op->paths = cp->isop->paths;
323	cp->isop->isop->paths = dp->paths;
324	for (istart = onode + 1; istart < Numspc; istart++) {
325	if (poolorder[onode] != poolorder[istart])
326	break;
327	}
328	for (i = istart; i < Numspc; i++) {
329	if (op == addposition[i])
330	addposition[i] = dp->kinp;
331	}
332	} else { /* multiway tree */
333	cp->isop = op->isop;
334	op->isop = cp;
335	dp = op->kinp;
336	}
337	} else { /* root */
338	dp = cp->isop->kinp;
339	if (op == tr->rootp)
340	tr->rootp = dp;
341	dp->isop = op->isop;
342	op->isop->isop->isop = dp;
343	op->isop = cp;
344	cp->isop->isop = op;
345	}
346	if (onode < Numspc - 1) {
347	#if 0
348	/* evaluate each subtrees */
349	if (onode < 6) putctopology(tr);
350	#endif
351	/* printf("%3d%3d%3d%3d ",
352	onode,dp->num+1,cp->kinp->num+1,cp->isop->num+1);
353	putctopology(tr); */
354	if (addposition[onode + 1] == NULL)
355	autoconstruction(tr, onode + 1, poolnode, addposition, poolorder);
356	else
357	wedge(tr, onode+1, poolnode, addposition, poolorder, addposition[onode+1]);
358	} else if (onode = Numspc - 1) {
359	/* Numibrnch = Maxibrnch;
360	Numbrnch = Maxbrnch; */
361	aproxtree(tr, Cnotree);
362	Cnotree++;
363	}
364	if (addposition[onode] != NULL) { /* internal node */
365	if (op->isop != cp) { /* binary tree */
366	for (i = istart; i < Numspc; i++) {
367	if (dp->kinp == addposition[i])
368	addposition[i] = op;
369	}
370	cp->isop->kinp = NULL;
371	cp->isop->isop->kinp = NULL;
372	op->kinp = dp;
373	dp->kinp = op;
374	op->paths = dp->paths;
375	} else { /* multiway tree */
376	op->isop = cp->isop;
377	cp->isop = NULL;
378	}
379	} else { /* root */
380	if (dp == tr->rootp)
381	tr->rootp = op;
382	dp->isop->isop->isop = op;
383	op->isop = dp->isop;
384	cp->isop->isop = cp;
385	dp->isop = NULL;
386	}
387	if (op->kinp->isop == NULL) {
388	return;
389	} else {
390	cp = op->kinp->isop;
391	while (cp != op->kinp) {
392	wedge(tr, onode, poolnode, addposition, poolorder, cp);
393	cp = cp->isop;
394	}
395	}
396	} /_ wedge /
397
398
399	void
400	autoconstruction(tr, onode, poolnode, addposition, poolorder)
401	Tree *tr;
402	int onode;
403	Node poolnode, addposition;
404	ivector poolorder;
405	{
406	Node *cp;
407
408	cp = tr->rootp;
409	do {
410	cp = cp->isop;
411	wedge(tr, onode, poolnode, addposition, poolorder, cp);
412	} while (cp != tr->rootp);
413	} /_ autoconstruction /
414
415
416	Node *
417	inbranode(tr, cpp, nenode, numorder, poolnode2, st)
418	Tree *tr;
419	char **cpp;
420	int *nenode;
421	int numorder;
422	Node ***poolnode2;
423	cvector st;
424	{
425	Node xp, yp, *np;
426	int i, k, n, dvg;
427	char *idp, ident[MAXWORD];
428
429	numorder++;
430	if (*cpp == ' ') (cpp)++;
431	if (*cpp == '{') { / internal node of binary tree */
432	if (Debug) printf("inbranode'{': %c\n", **cpp);
433	(*cpp)++;
434	xp = inbranode(tr, cpp, nenode, numorder, poolnode2, st); /first node/
435	dvg = 0;
436	for (n = 0; poolnode2[numorder][n] != NULL; n++)
437	;
438	k = 1;
439	while (*cpp != '}') { / second node and others */
440	np = inbranode(tr, cpp, nenode, numorder, poolnode2, st);
441	np->isop = new_dnode();
442	np->isop->isop = new_anode();
443	np->isop->isop->isop = np;
444	poolnode2[numorder][n + dvg] = np;
445	np->isop->kinp = xp; /* add position */
446	dvg++;
447	Numtplgy = k; / printf("Numtplgy: %.1f\n", Numtplgy); */
448	k += 2;
449	}
450	if (Debug) printf("inbranode'}': %c\n", **cpp);
451	poolnode2[numorder][n + dvg] = NULL;
452	if (dvg < 1) {
453	fprintf(stderr, "ERROR, constrained tree:\n%s\n", st);
454	fprintf(stderr, "redundancy or unnecessary \"{}\" !\n");
455	exit(1);
456	}
457	(*cpp)++;
458	(*nenode)++;
459	return xp;
460	} else if (*cpp == '(') { / internal node of multiway tree */
461	if (Debug) printf("inbranode'(': %c\n", **cpp);
462	(*cpp)++;
463	xp = inbranode(tr, cpp, nenode, numorder, poolnode2, st); /first node/
464	dvg = 0;
465	for (n = 0; poolnode2[numorder][n] != NULL; n++)
466	;
467	if (*cpp != ')') { / second node */
468	np = inbranode(tr, cpp, nenode, numorder, poolnode2, st);
469	np->isop = new_dnode();
470	np->isop->isop = new_anode();
471	np->isop->isop->isop = np;
472	poolnode2[numorder][n + dvg] = np;
473	np->isop->kinp = xp; /* add position */
474	dvg++;
475	}
476	while (*cpp != ')') { / third node and others */
477	yp = np;
478	np = inbranode(tr, cpp, nenode, numorder, poolnode2, st);
479	np->isop = NULL;
480	poolnode2[numorder][n + dvg] = np;
481	np->kinp->isop = yp; /* add position (kinp) */
482	dvg++;
483	}
484	if (Debug) printf("inbranode')': %c\n", **cpp);
485	poolnode2[numorder][n + dvg] = NULL;
486	if (dvg < 1) {
487	fprintf(stderr, "ERROR, constrained tree:\n%s\n", st);
488	fprintf(stderr, "redundancy or unnecessary \"()\" !\n");
489	exit(1);
490	}
491	(*cpp)++;
492	(*nenode)++;
493	return xp;
494
495	} else if (isalnum(*cpp)) { / external node */
496	if (Debug) printf("external: %c\n", **cpp);
497	for (idp = ident; cpp!=' ' && cpp!='(' && **cpp!=')'
498	&& cpp!='{' && cpp!='}'; (*cpp)++) {
499	idp++ = *cpp;
500	if (Debug) putchar(**cpp);
501	}
502	*idp = '\0';
503	if (Debug) putchar('\n');
504	for (i = 0; i < Numspc; i++) {
505	/* puts(Identif[i]); */
506	if (!strcmp(ident, Identif[i])) {
507	return tr->ebrnchp[i]->kinp;
508	}
509	}
510	fprintf(stderr, "ERROR, constrained tree:\n%s\n", st);
511	fprintf(stderr, "ERROR, abnormal identifier(OTU name): %s !\n", ident);
512	exit(1);
513	} else {
514	fprintf(stderr, "ERROR, constrained tree:\n%s\n", st);
515	fprintf(stderr, "ERROR, abnormal character: %s\n", **cpp);
516	exit(1);
517	}
518	return NULL;
519	} /_ inbranode /
520
521
522	void
523	streeinit(tr, strtree, poolnode, addposition, poolorder)
524	Tree *tr;
525	cvector strtree;
526	Node poolnode, addposition;
527	ivector poolorder;
528	{
529	char *chp;
530	int i, j, k, nenode, ninode, dvg, numorder;
531	Node sp0, sp1, sp2, np, ***poolnode2;
532
533	poolnode2 = new_nodematrix(Numspc, Numspc);
534	nenode = 0;
535	numorder = 0;
536	for (i = 0; i < Numspc; i++) {
537	poolnode2[i][0] = NULL;
538	poolnode[i] = NULL;
539	addposition[i] = NULL;
540	poolorder[i] = 0;
541	}
542
543	chp = strtree;
544	if (chp == '{') { / binary */
545	chp++;
546	sp0 = inbranode(tr, &chp, &nenode, numorder, poolnode2, strtree);
547	sp1 = inbranode(tr, &chp, &nenode, numorder, poolnode2, strtree);
548	dvg = 0;
549	k = 3;
550	while (chp != '}') { / fourth node and others */
551	np = inbranode(tr, &chp, &nenode, numorder, poolnode2, strtree);
552	if (*chp == '}') break;
553	np->isop = new_dnode();
554	np->isop->kinp = new_anode();
555	np->isop->kinp->paths = np->isop->paths;
556	np->isop->isop = np;
557	np->isop->kinp->kinp = np->isop;
558	poolnode2[0][dvg] = np;
559	dvg++;
560	Numtplgy *= k;
561	k += 2;
562	}
563	poolnode2[0][dvg] = NULL;
564	sp2 = np;
565	sp0->isop = sp1;
566	sp1->isop = sp2;
567	sp2->isop = sp0;
568	tr->rootp = sp2;
569	} else if (chp == '(') { / multiway */
570	chp++;
571	sp0 = inbranode(tr, &chp, &nenode, numorder, poolnode2, strtree);
572	sp1 = inbranode(tr, &chp, &nenode, numorder, poolnode2, strtree);
573	sp2 = inbranode(tr, &chp, &nenode, numorder, poolnode2, strtree);
574	sp0->isop = sp1;
575	sp1->isop = sp2;
576	sp2->isop = sp0;
577	tr->rootp = sp2;
578	dvg = 0;
579	while (chp != ')') { / fourth node and others */
580	np = inbranode(tr, &chp, &nenode, numorder, poolnode2, strtree);
581	np->isop = NULL;
582	poolnode2[0][dvg] = np;
583	np->kinp->isop = sp2; /* add position (kinp) */
584	dvg++;
585	}
586	poolnode2[0][dvg] = NULL;
587	} else {
588	fprintf(stderr, "ERROR, constrained tree:\n%s\n", strtree);
589	fprintf(stderr, "ERROR, abnormal character: %s\n", *chp);
590	exit(1);
591	}
592
593	if (Debug_optn) {
594	putchar('\n');
595	printf(" %4d%4d%4d\n", sp0->num+1, sp1->num+1, sp2->num+1);
596	for (i = 0; i < Numspc; i++) {
597	printf("%4d", i);
598	if (poolnode2[i][0] != NULL) {
599	for (j = 0; j < Numspc; j++) {
600	if (poolnode2[i][j] != NULL)
601	printf("%4d", poolnode2[i][j]->num+1);
602	else
603	break;
604	} putchar('\n');
605	} else {
606	printf(" nul\n");
607	}
608	}
609	}
610
611	poolnode[0] = sp0;
612	poolnode[1] = sp1;
613	poolnode[2] = sp2;
614	addposition[0] = NULL;
615	addposition[1] = NULL;
616	addposition[2] = NULL;
617	for (k = 3, ninode = 0, i = 0; i < Numspc; i++) {
618	for (j = 0; j < Numspc && poolnode2[i][j] != NULL; j++) {
619	poolorder[k] = i;
620	poolnode[k] = poolnode2[i][j];
621	if (poolnode[k]->isop != NULL) { /* binary tree */
622	poolnode[k]->isop->num = ninode + Maxspc;
623	if (i == 0) { /* root */
624	addposition[k] = NULL;
625	} else { /* not root */
626	addposition[k] = poolnode[k]->isop->kinp;
627	poolnode[k]->isop->kinp = NULL;
628	}
629	tr->ibrnchp[ninode] = poolnode[k]->isop;
630	ninode++;
631	} else { /* multiway tree */
632	addposition[k] = poolnode[k]->kinp->isop;
633	poolnode[k]->kinp->isop = NULL;
634	}
635	k++;
636	}
637	}
638	Numibrnch = ninode;
639	Numbrnch = Numspc + ninode;
640
641	if (Debug_optn) {
642	putchar('\n');
643	for (i = 0; i < k; i++) {
644	printf("%4d ", i);
645	fputid(stdout, Identif[poolnode[i]->kinp->num], 10);
646	printf("%4d", poolnode[i]->kinp->num+1);
647	printf("%4d", poolnode[i]->num+1);
648	if (addposition[i] != NULL) {
649	printf("%4d", addposition[i]->kinp->num+1);
650	printf("%4d", addposition[i]->num+1);
651	} else {
652	fputs(" nul", stdout);
653	fputs(" nul", stdout);
654	}
655	printf(" %4d ", poolorder[i]);
656	if (poolnode[i]->isop != NULL) {
657	np = poolnode[i];
658	while((np = np->isop) != poolnode[i]) {
659	printf("%4d", np->num+1);
660	}
661	} putchar('\n');
662	}
663	printf("Numbrnch =%4d\n", Numbrnch);
664	printf("Numibrnch =%4d\n", Numibrnch);
665	}
666	free_nodematrix(poolnode2);
667	} /_ streeinit /
668
669
670	void
671	atreeinit(tr, poolnode, addposition, poolorder)
672	Tree *tr;
673	Node poolnode, addposition;
674	ivector poolorder;
675	{
676	char *chp;
677	char linebuf[BUFLINE];
678	int i, j, k, nenode, ninode, dvg, numorder, mini;
679	double dis, mindis;
680	Node sp0, sp1, sp2, np;
681	Node ***poolnode2;
682
683	mindis = fabs(Distanmat[0][1] - Distanmat[Numspc-1][1]);
684	mini = 1;
685	for ( i = 2; i < Numspc - 1; i++) {
686	dis = fabs(Distanmat[0][i] - Distanmat[Numspc-1][i]);
687	if (dis < mindis) {
688	mindis = dis;
689	mini = i;
690	}
691	}
692
693	poolnode2 = new_nodematrix(Numspc, Numspc);
694	nenode = 0;
695	numorder = 0;
696	for (i = 0; i < Numspc; i++) {
697	poolnode2[i][0] = NULL;
698	poolnode[i] = NULL;
699	addposition[i] = NULL;
700	poolorder[i] = 0;
701	}
702
703	strcpy(linebuf, Identif[0]);
704	strcat(linebuf, " ");
705	chp = linebuf;
706	sp0 = inbranode(tr, &chp, &nenode, numorder, poolnode2, linebuf);
707	strcpy(linebuf, Identif[mini]);
708	strcat(linebuf, " ");
709	chp = linebuf;
710	sp1 = inbranode(tr, &chp, &nenode, numorder, poolnode2, linebuf);
711	strcpy(linebuf, Identif[Numspc-1]);
712	strcat(linebuf, " ");
713	chp = linebuf;
714	sp2 = inbranode(tr, &chp, &nenode, numorder, poolnode2, linebuf);
715	sp0->isop = sp1;
716	sp1->isop = sp2;
717	sp2->isop = sp0;
718	tr->rootp = sp2;
719	dvg = 0;
720	for ( i = 1; i < Numspc - 1; i++) {
721	if (i != mini) {
722	strcpy(linebuf, Identif[i]);
723	strcat(linebuf, " ");
724	chp = linebuf;
725	np = inbranode(tr, &chp, &nenode, numorder, poolnode2, linebuf);
726	np->isop = new_dnode();
727	np->isop->kinp = new_anode();
728	np->isop->kinp->paths = np->isop->paths;
729	np->isop->isop = np;
730	np->isop->kinp->kinp = np->isop;
731	poolnode2[0][dvg] = np;
732	dvg++;
733	}
734	}
735	poolnode2[0][dvg] = NULL;
736
737	if (Debug_optn) {
738	putchar('\n');
739	printf("%4d%4d%4d\n", sp0->num+1, sp1->num+1, sp2->num+1);
740	for (j = 0; j < Numspc && poolnode2[0][j] != NULL; j++)
741	printf("%4d", poolnode2[0][j]->num+1);
742	putchar('\n');
743	}
744
745	poolnode[0] = sp0;
746	poolnode[1] = sp1;
747	poolnode[2] = sp2;
748	addposition[0] = NULL;
749	addposition[1] = NULL;
750	addposition[2] = NULL;
751	k = 3;
752	ninode = 0;
753	for (j = 0; j < Numspc && poolnode2[0][j] != NULL; j++) {
754	poolorder[k] = 0;
755	poolnode[k] = poolnode2[0][j];
756	poolnode[k]->isop->num = ninode + Maxspc;
757	addposition[k] = NULL;
758	tr->ibrnchp[ninode] = poolnode[k]->isop;
759	ninode++;
760	k++;
761	}
762	Numbrnch = ninode;
763
764	if (Debug_optn) {
765	putchar('\n');
766	for (i = 0; i < k; i++) {
767	printf("%4d ", i);
768	fputid(stdout, Identif[poolnode[i]->kinp->num], 10);
769	printf("%4d", poolnode[i]->kinp->num+1);
770	printf("%4d", poolnode[i]->num+1);
771	if (addposition[i] != NULL) {
772	printf("%4d", addposition[i]->kinp->num+1);
773	printf("%4d", addposition[i]->num+1);
774	} else {
775	fputs(" nul", stdout);
776	fputs(" nul", stdout);
777	}
778	printf(" %4d ", poolorder[i]);
779	if (poolnode[i]->isop != NULL) {
780	np = poolnode[i];
781	while((np = np->isop) != poolnode[i]) {
782	printf("%4d", np->num+1);
783	}
784	} putchar('\n');
785	}
786	printf("Numbrnch =%4d\n", Numbrnch);
787	}
788	free_nodematrix(poolnode2);
789	} /_ atreeinit /
790
791
792	void
793	tablealtree(nt)
794	int nt;
795	{
796	Infoaltree cp, *info;
797	int i, j, k, max;
798
799	info = (Infoaltree *)malloc((unsigned)nt sizeof(Infoaltree *));
800	if (info == NULL) maerror("in tablealtree().");
801	for (cp = Atail.up, i = nt; cp != &Ahead; cp = cp->up) {
802	info[--i] = cp;
803	/* printf("%.1f %s\n", cp->lklaprox, cp->ltplgy); */
804	}
805	printf("%d / %d %s %s \"%s\"", nt, Cnotree, Prog_name, VERSION, Modelname);
806	printf(" %d OTUs %d sites. %s\n", Maxspc, Numsite, Comment);
807
808	if (Mevol_optn) { /* total branch length */
809	printf("# top ranking %1d trees by TBL criterion (Minimum Evolution)\n",
810	Maxaltree);
811	} else { /* approximate ln likelihood */
812	printf("# <= %1d trees (top ranking for approx. ln L)", Maxaltree);
813	printf(" in the top %.1f%% range of TBL\n", Tblrate * 100.0);
814	}
815
816	for (i = 1, max = Tblbin[0]; i < NUMTBLBIN; i++) {
817	if (Tblbin[i] > max) max = Tblbin[i];
818	}
819	printf("# %5s%8s %9s\n", "range", "TBL ", "trees");
820	printf("# < %8.2f %9d\n", Mintbldm, Tblunder);
821	for (i = 0; i < NUMTBLBIN; i++) {
822	printf("# %3.0f%% %8.2f %9d ",
823	(i+1)*100.0/NUMTBLBIN, ((i+1.0)/Tblcoef)+Mintbldm, Tblbin[i]);
824	k = (int)(Tblbin[i] * 50.0 / max);
825	for (j = 0; j < k; j++) putchar('*');
826	putchar('\n');
827	}
828	printf("# %3s %8s %9d\n", "", "over", Tblover);
829
830	putchar('#');
831	if (!Mevol_optn) { /* approximate ln likelihood */
832	printf(" approx. ln L %.1f", info[0]->lklaprox);
833	printf(" ... %.1f", info[nt-1]->lklaprox);
834	printf(" diff %.1f,", info[0]->lklaprox - info[nt-1]->lklaprox);
835	}
836	printf(" TBL %.1f", info[0]->tbl);
837	printf(" ... %.1f", info[nt-1]->tbl);
838	printf(" diff %.1f\n", info[nt-1]->tbl - info[0]->tbl);
839
840	for (i = 0; i < nt; i++) {
841	if (Info_optn ) printf("%.1f ", info[i]->lklaprox);
842	fputs(info[i]->ltplgy, stdout);
843	printf("; %.1f\n", info[0]->lklaprox - info[i]->lklaprox);
844	}
845	} /_ tablealtree /

Note: See TracBrowser for help on using the repository browser.

Download in other formats: