Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

source: tags/initial/GDE/MOLPHY/mltree.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: 39.6 KB

Line
1	/*
2	* mltree.c Adachi, J. 1996.03.01
3	* Copyright (C) 1992-1996 J. Adachi & M. Hasegawa, All rights reserved.
4	*/
5
6	#include "protml.h"
7
8	#define SLSDEBUG 0
9	#define MLSITE 0
10	#define USERTREEDEBUG 0
11
12
13	Node **
14	new_npvector(n)
15	int n;
16	/* memory allocate a node pointer vector */
17	{
18	Node **v;
19
20	v = (Node *) malloc((unsigned)n sizeof(Node *));
21	if (v == NULL) maerror("in new_npvector().");
22	return v;
23	}
24
25	void
26	free_npvector(v)
27	Node **v;
28	{
29	free((char *) v);
30	}
31
32
33	Tree *
34	new_tree(maxspc, maxibrnch, numptrn, seqconint)
35	int maxspc, maxibrnch;
36	imatrix seqconint;
37	{
38	int n, i;
39	Tree *tr;
40	Node dp, up;
41
42	tr = (Tree *) malloc(sizeof(Tree));
43	if (tr == NULL) maerror("tr in new_tree().");
44	tr->ebrnchp = (Node *) malloc((unsigned)maxspc sizeof(Node *));
45	if (tr->ebrnchp == NULL) maerror("ebrnchp in new_tree().");
46	tr->ibrnchp = (Node *) malloc((unsigned)maxibrnch sizeof(Node *));
47	if (tr->ibrnchp == NULL) maerror("ibrnchp in new_tree().");
48	tr->bturn = new_ivector(maxspc);
49	for (n = 0; n < maxspc; n++) {
50	tr->bturn[n] = n;
51	dp = (Node *) malloc(sizeof(Node));
52	if (dp == NULL) maerror("dp in new_tree().");
53	up = (Node *) malloc(sizeof(Node));
54	if (up == NULL) maerror("up in new_tree().");
55	dp->isop = NULL;
56	up->isop = NULL;
57	dp->kinp = up;
58	up->kinp = dp;
59	dp->descen = TRUE;
60	up->descen = FALSE;
61	dp->num = n;
62	up->num = n;
63	dp->length = 0.0;
64	up->length = 0.0;
65	dp->lklhdl = 0.0;
66	up->lklhdl = 0.0;
67	dp->paths = new_ivector(maxspc);
68	up->paths = dp->paths;
69	for (i = 0; i < maxspc; i++) dp->paths[i] = 0;
70	dp->paths[n] = 1;
71	dp->eprob = seqconint[n];
72	up->eprob = NULL;
73	dp->iprob = NULL;
74	up->iprob = new_dmatrix(numptrn, Tpmradix);
75	tr->ebrnchp[n] = dp;
76	}
77	for (n = 0; n < maxibrnch; n++) {
78	dp = (Node *) malloc(sizeof(Node));
79	if (dp == NULL) maerror("dp in new_tree().");
80	up = (Node *) malloc(sizeof(Node));
81	if (up == NULL) maerror("up in new_tree().");
82	dp->isop = NULL;
83	up->isop = NULL;
84	dp->kinp = up;
85	up->kinp = dp;
86	dp->descen = TRUE;
87	up->descen = FALSE;
88	dp->num = n + maxspc;
89	up->num = n + maxspc;
90	dp->length = 0.0;
91	up->length = 0.0;
92	dp->lklhdl = 0.0;
93	up->lklhdl = 0.0;
94	dp->paths = new_ivector(maxspc);
95	up->paths = dp->paths;
96	for (i = 0; i < maxspc; i++) dp->paths[i] = 0;
97	dp->eprob = NULL;
98	up->eprob = NULL;
99	dp->iprob = new_dmatrix(numptrn, Tpmradix);
100	up->iprob = new_dmatrix(numptrn, Tpmradix);
101	tr->ibrnchp[n] = dp;
102	}
103	tr->rootp = NULL;
104
105	return tr;
106	} /_ new_tree /
107
108
109	int
110	getbuftree(numspc, identif)
111	int numspc;
112	cmatrix identif;
113	{
114	int i, buf;
115
116	buf = numspc * 3; /* - 3 */
117	for (i = 0; i < numspc; i++)
118	buf += strlen(identif[i]);
119	return buf;
120	} /_ getbuftree /
121
122
123	void
124	changedistan(distanmat, distanvec, numspc)
125	dmatrix distanmat;
126	dvector distanvec;
127	int numspc;
128	{
129	int i, j, k;
130
131	for (k = 0, i = 1; i < numspc; i++) {
132	for (j = 0; j < i; j++, k++)
133	distanvec[k] = distanmat[i][j];
134	}
135	}
136
137
138	void
139	getproportion(proportion, distanvec, maxpair)
140	double *proportion;
141	dvector distanvec;
142	int maxpair;
143	{
144	int i;
145	double maxdis;
146
147	maxdis = 0.0;
148	for (i = 0; i < maxpair; i++) {
149	if (distanvec[i] > maxdis) {
150	maxdis = distanvec[i];
151	}
152	}
153	proportion = (double)MAXCOLUMN / (maxdis 3.0);
154	if (proportion > 1.0) proportion = 1.0;
155	if (Debug) printf("Proportion: %.5f maxdis: %.5f\n", *proportion,maxdis);
156	} /* getproportion */
157
158
159	Infotree *
160	newinfotrees(numtree, buftree)
161	int numtree;
162	int buftree;
163	{
164	int i;
165	Infotree *info;
166	cvector m;
167
168	info = (Infotree ) malloc((unsigned)numtree sizeof(Infotree));
169	if (info == NULL) maerror("in newinfotrees().");
170	m = (cvector) malloc((unsigned)(numtree * buftree) * sizeof(char));
171	if (m == NULL) maerror("2 in newinfotrees().");
172	for (i = 0; i < numtree; i++, m+=buftree)
173	info[i].ltplgy = m;
174
175	return info;
176	} /_ newinfotrees /
177
178
179	Infoaltree *
180	newinfoaltrees(numaltree, buftree)
181	int numaltree;
182	int buftree;
183	{
184	int i;
185	Infoaltree *info;
186	cvector m;
187
188	info = (Infoaltree ) malloc((unsigned)numaltree sizeof(Infoaltree));
189	if (info == NULL) maerror("in newinfoaltrees().");
190	m = (cvector) malloc((unsigned)(numaltree * buftree) * sizeof(char));
191	if (m == NULL) maerror("2 in newinfoaltrees().");
192	for (i = 0; i < numaltree; i++, m+=buftree)
193	info[i].ltplgy = m;
194
195	return info;
196	} /_ newinfoaltrees /
197
198
199	void
200	getnumtree(ifp, numtree)
201	FILE *ifp;
202	int *numtree;
203	{
204	char *cp;
205	char line[BUFLINE];
206
207	while (fgets(line, BUFLINE, ifp) != NULL) {
208	/* fputs(line, ofp); */
209	if (sscanf(line, "%d", numtree) == 1) {
210	if (Verbs_optn && User_optn && !Aprox_optn)
211	fprintf(stderr,"%d trees\n", numtree); / stdout */
212	return;
213	} else {
214	if (*line == '#')
215	continue;
216	for (cp = line; isspace(cp) && cp != '\0'; cp++)
217	if (*cp != '\0') {
218	fputs(line, stderr);
219	fprintf(stderr,
220	"\nCan't read number of user trees, bad format.\n");
221	exit(1);
222	} else {
223	continue;
224	}
225	}
226	}
227	fprintf(stderr, "\nCan't read number of user trees.\n");
228	exit(1);
229	} /_ getnumtree /
230
231
232	void
233	getusertree(ifp, strtree, buftree)
234	FILE *ifp;
235	cvector strtree;
236	int buftree;
237	{
238	char line[BUFLINE];
239	char cp, np, xp, ap, *bp;
240	int par, bra, nsp;
241	boolean word;
242
243	word = FALSE;
244	par = bra = nsp = 0;
245	np = strtree;
246	strtree[buftree - 1] = '\0';
247	while (fgets(line, BUFLINE, ifp) != NULL) {
248	if (!isalnum(*line)) word = FALSE;
249	/* printf("%s\n", line); */
250	if (line == '#') continue; / comment line */
251	for (cp = line; *cp != '\0'; ) {
252	switch(*cp) {
253
254	case ';':
255	if (np == strtree) goto next;
256	if (par == 0 && bra == 0) { /* tree */
257	*np = '\0';
258	if (nsp != Numspc) {
259	fprintf(stderr,
260	"ERROR: different number of OTU, %d/%d. \n",
261	nsp, Numspc);
262	fprintf(stderr, "%s\n", strtree);
263	exit(1);
264	}
265	/* printf("%s\n", strtree); */
266	return;
267	} else { /* no tree */
268	if (par != 0) {
269	fprintf(stderr,
270	"ERROR, bad match parenthesis \"()\" !\n");
271	if (par > 0)
272	fprintf(stderr,"')' is %d less than '('.\n",
273	par);
274	if (par < 0)
275	fprintf(stderr,"')' is %d more than '('.\n",
276	abs(par));
277	}
278	if (bra != 0) {
279	fprintf(stderr,
280	"ERROR, bad match brace \"{}\" !\n");
281	if (bra > 0)
282	fprintf(stderr,"'}' is %d less than '{'.\n",
283	par);
284	if (par < 0)
285	fprintf(stderr,"'}' is %d more than '{'.\n",
286	abs(par));
287	}
288	fprintf(stderr, "%s\n", strtree);
289	exit(1);
290	}
291	break;
292	case '(': np++ = cp++; par++; break;
293	case ')': np++ = cp++; par--; break;
294	case '{': np++ = cp++; bra++; break;
295	case '}': np++ = cp++; bra--; break;
296	default:
297	if (isspace(cp) \|\| cp == ',') {
298	cp++;
299	} else { /* identifier */
300	if (isalnum(cp) \|\| cp == '_' \|\| *cp == '-') {
301	if (!word) {
302	np--;
303	if (isalnum(np) \|\| np == '_' \|\| *np == '-') {
304	np++;
305	*np = ' ';
306	}
307	np++;
308	}
309	ap = cp; bp = np;
310	while (isalnum(cp) \|\| cp == '_' \|\| *cp == '-') {
311	np++ = cp++;
312	}
313	if (!word) nsp++;
314	*cp == '\0' ? (word = TRUE) : (word = FALSE);
315	if (Debug) printf("%3d %1d %-10.10s %-20.20s\n",
316	nsp,word,bp,ap);
317	} else {
318	fprintf(stderr, "bad name: %s\n", cp);
319	exit(1);
320	}
321	}
322
323	} /* switch */
324	} /* while */
325	if (par == 0 && bra == 0 && np != strtree) {
326	*np = '\0';
327	/* printf("%s\n", strtree); */
328	return;
329	}
330	if (np > (strtree + buftree - 2)) {
331	fprintf(stderr, "Too long, users tree.\n");
332	fprintf(stderr, "You may forget tree delimiter \";\".\n");
333	fprintf(stderr, "\"%s\"\n", strtree);
334	exit(1);
335	}
336	#if USERTREEDEBUG
337	printf("%s\n", strtree);
338	#endif
339	next: ;
340	}
341	fprintf(stderr, "Can't read users trees of input file.\n");
342	exit(1);
343	} /_ getusertree /
344
345
346	Node *
347	internalnode(tr, cpp, ninode, st)
348	Tree *tr;
349	char **cpp;
350	int *ninode;
351	char *st;
352	{
353	Node xp, np, *rp;
354	int i, j, dvg;
355	char ident[MAXWORD];
356	char *idp;
357
358	/* fprintf(stderr, "in1:%s\n", cpp); /
359
360	if (*cpp == ' ') (cpp)++;
361	if (**cpp == '(') {
362	if (Debug) printf("internal1: %c\n", **cpp);
363	(*cpp)++;
364	xp = internalnode(tr, cpp, ninode, st);
365	xp->isop = xp;
366	dvg = 1;
367	while (**cpp != ')') {
368	dvg++;
369	np = internalnode(tr, cpp, ninode, st);
370	np->isop = xp->isop;
371	xp->isop = np;
372	xp = np;
373	}
374	if (Debug) printf("internal2: %c\n", **cpp);
375	(*cpp)++;
376	if (dvg < 2) {
377	fputs("ERROR, redundancy or unnecessary \"()\" !\n", stderr);
378	fprintf(stderr, "%s\n", st);
379	exit(1);
380	}
381	/* printf("ninode %d\n", ninode + 1); /
382	rp = tr->ibrnchp[*ninode];
383	rp->isop = xp->isop;
384	xp->isop = rp;
385
386	for (j = 0; j < Numspc; j++) rp->paths[j] = 0;
387	for (xp = rp->isop; xp != rp; xp = xp->isop) {
388	for (j = 0; j < Numspc; j++) {
389	if (xp->paths[j] == 1) rp->paths[j] = 1;
390	}
391	}
392	/*
393	if (Debug) {
394	for (j = 0; j < Numspc; j++) printf("%2d",rp->paths[j]);
395	putchar('\n');
396	}
397	*/
398	(*ninode)++;
399	if (*ninode > Maxibrnch) {
400	fputs("ERROR, too many number of internal branch!\n", stderr);
401	fputs("This tree may be rooted tree.\n", stderr);
402	fprintf(stderr, "%s\n", st);
403	exit(1);
404	}
405
406	/* fprintf(stderr, "in2:%s\n", cpp); /
407
408	return rp->kinp;
409	} else if (isalnum(**cpp)) {
410	if (Debug) printf("external: %c\n", **cpp);
411	for (idp = ident; cpp!=' ' && cpp!='(' && *cpp!=')'; (cpp)++) {
412	idp++ = *cpp;
413	if (Debug) putchar(**cpp);
414	}
415	*idp = '\0';
416	if (Debug) putchar('\n');
417
418	for (i = 0; i < Numspc; i++) {
419	/* puts(Identif[i]); */
420	if (!strcmp(ident, Identif[i])) {
421	return tr->ebrnchp[i]->kinp;
422	}
423	}
424	fprintf(stderr, "ERROR, abnormal identifier(OTU name): %s !\n", ident);
425	exit(1);
426	} else {
427	fprintf(stderr, "ERROR, abnormal tree topology character: %s\n", *cpp);
428	exit(1);
429	}
430	return NULL;
431	} /_ internalnode /
432
433
434	void
435	constructtree(tr, strtree)
436	Tree *tr;
437	cvector strtree;
438	{
439	char *cp;
440	int ninode, dvg;
441	Node xp, np;
442
443	#if USERTREEDEBUG
444	printf("%s\n", strtree);
445	#endif
446	ninode = 0;
447	cp = strtree;
448	if (*cp == '(') {
449	if (Debug) printf("roottre0: %c\n", *cp);
450	cp++;
451	xp = internalnode(tr, &cp, &ninode, strtree);
452	xp->isop = xp;
453	dvg = 1;
454	while (*cp != ')') {
455	dvg++;
456	np = internalnode(tr, &cp, &ninode, strtree);
457	np->isop = xp->isop;
458	xp->isop = np;
459	xp = np;
460	}
461	if (Debug) printf("roottre0: %c\n", *cp);
462	if (dvg < 3) {
463	fputs("ERROR: this tree is not unroot tree!\n", stderr);
464	fprintf(stderr, "%s\n", strtree);
465	exit(1);
466	}
467	tr->rootp = xp;
468	Numibrnch = ninode;
469	Numbrnch = Numspc + ninode;
470	} else {
471	fprintf(stderr, "ERROR users tree:\n");
472	fprintf(stderr, "%s\n", strtree);
473	exit(1);
474	}
475	} /_ constructtree /
476
477
478	void
479	prcurtree(tr)
480	Tree *tr;
481	{
482	Node dp, up, *cp;
483	int i, j;
484	double sum;
485
486	/* printf("\nStructure of Tree\n"); */
487	printf("\n%4s%5s%5s%5s%7s%7s%37s\n",
488	"num","kinp","isop","isop","descen","length", "namesp");
489	for (i = 0; i < Numspc; i++) {
490	dp = tr->ebrnchp[i];
491	up = dp->kinp;
492	printf("%4d", dp->num+1);
493	printf("%5d", up->num+1);
494	if (dp->isop == NULL) printf("%5s", "null");
495	else printf("%5d", dp->isop->num+1);
496	if (up->isop == NULL) printf("%5s", "null");
497	else printf("%5d", up->isop->num+1);
498	printf("%3d", dp->descen);
499	printf("%3d", up->descen);
500	printf("%8.3f", dp->length);
501	printf("%8.3f", up->length);
502	printf("%5d%5d", dp->eprob[0],dp->eprob[Numptrn-1]);
503	for (sum = 0.0, j = 0; j < Tpmradix; j++) sum += up->iprob[0][j];
504	printf("%5.0f", sum * 1000.0);
505	for (sum = 0.0, j = 0; j < Tpmradix; j++) sum += up->iprob[Numptrn-1][j];
506	printf("%5.0f", sum * 1000.0);
507	printf(" %s", Identif[dp->num]);
508	putchar('\n');
509	}
510	for (i = 0; i < Numibrnch; i++) {
511	dp = tr->ibrnchp[i];
512	up = dp->kinp;
513	printf("%4d", dp->num+1);
514	printf("%5d", up->num+1);
515	if (dp->isop == NULL) printf("%5s", "null");
516	else printf("%5d", dp->isop->num+1);
517	if (up->isop == NULL) printf("%5s", "null");
518	else printf("%5d", up->isop->num+1);
519	printf("%3d", dp->descen);
520	printf("%3d", up->descen);
521	printf("%8.3f", dp->length);
522	printf("%8.3f", up->length);
523	for (sum = 0.0, j = 0; j < Tpmradix; j++) sum += dp->iprob[0][j];
524	printf("%5.0f", sum * 1000.0);
525	for (sum = 0.0, j = 0; j < Tpmradix; j++) sum += dp->iprob[Numptrn-1][j];
526	printf("%5.0f", sum * 1000.0);
527	for (sum = 0.0, j = 0; j < Tpmradix; j++) sum += up->iprob[0][j];
528	printf("%5.0f", sum * 1000.0);
529	for (sum = 0.0, j = 0; j < Tpmradix; j++) sum += up->iprob[Numptrn-1][j];
530	printf("%5.0f", sum * 1000.0);
531	for (cp = dp->isop; cp != dp; cp = cp->isop) {
532	printf("%5d", cp->num+1);
533	}
534	putchar('\n');
535	}
536	dp = tr->rootp->isop;
537	printf("%4d", dp->num+1);
538	for (cp = dp->isop; cp != dp; cp = cp->isop) {
539	printf("%5d", cp->num+1);
540	}
541	putchar('\n');
542	} /_ prcurtree /
543
544
545	void
546	pathing(tr)
547	Tree *tr;
548	{
549	int i, j;
550	ivector pthsd, pthsa;
551	Node rp, cp, xp, dp;
552
553	cp = rp = tr->rootp;
554	do {
555	cp = cp->isop->kinp;
556	if (cp->isop == NULL) { /* external node */
557	/*
558	pthsd = cp->paths;
559	*/
560	cp = cp->kinp;
561	/*
562	pthsa = cp->paths;
563	for (i = 0; i < Numspc; i++) {
564	pthsd[i] == 1 ? (pthsa[i] = 0) : (pthsa[i] = 1);
565	}
566	*/
567	} else { /* internal node */
568	if (!cp->descen) { /* ascent */
569	dp = cp->kinp;
570	pthsd = dp->paths;
571	for (i = 0; i < Numspc; i++) pthsd[i] = 0;
572	for (xp = dp->isop; xp != dp; xp = xp->isop) {
573	pthsa = xp->kinp->paths;
574	for (i = 0; i < Numspc; i++) {
575	if (pthsa[i]) pthsd[i] = 1;
576	}
577	}
578	/*
579	pthsa = cp->paths;
580	for (i = 0; i < Numspc; i++) {
581	pthsd[i] == 1 ? (pthsa[i] = 0) : (pthsa[i] = 1);
582	}
583	*/
584	}
585	}
586	} while (cp != rp);
587	/* */
588	if (Debug_optn) {
589	puts("\npath of internal nodes");
590	for (i = 0; i < Maxibrnch; i++) {
591	pthsd = tr->ibrnchp[i]->paths;
592	for (j = 0; j < Numspc; j++) printf("%2d", pthsd[j]); putchar('\n');
593	}
594	}
595	/* */
596	} /_ pathing /
597
598
599	#if 0
600	static void
601	leastsquares(am, ym)
602	dmatrix am;
603	dvector ym;
604	{
605	int i, j, k;
606	double pivot, element;
607	dmatrix im;
608
609	im = new_dmatrix(Numbrnch, Numbrnch);
610	for (i = 0; i < Numbrnch; i++) {
611	for (j = 0; j < Numbrnch; j++)
612	im[i][j] = 0.0;
613	im[i][i] = 1.0;
614	}
615	for (k = 0; k < Numbrnch; k++) {
616	pivot = am[k][k];
617	ym[k] /= pivot;
618	for (j = 0; j < Numbrnch; j++) {
619	am[k][j] /= pivot;
620	im[k][j] /= pivot;
621	}
622	for (i = 0; i < Numbrnch; i++) {
623	if (k != i) {
624	element = am[i][k];
625	ym[i] -= element * ym[k];
626	for (j = 0; j < Numbrnch; j++) {
627	am[i][j] -= element * am[k][j];
628	im[i][j] -= element * im[k][j];
629	}
630	}
631	}
632	}
633	/*
634	putchar('\n');
635	for (i = 0; i < Numbrnch; i++) {
636	for (j = 0; j < Numbrnch; j++)
637	printf("%10.6f", im[i][j]);
638	putchar('\n');
639	}
640	*/
641	free_dmatrix(im);
642	} /_ leastsquares /
643	#endif
644
645
646	static void
647	luequation(amat, yvec, size)
648	dmatrix amat;
649	dvector yvec;
650	int size;
651	{
652	/* SOLVE THE LINEAR SET OF EQUATIONS ON LU DECOMPOSITION */
653	double eps = 1.0e-20; /* ! */
654	int i, j, k, l, maxi, idx;
655	double sum, tmp, maxb, aw;
656	dvector wk;
657	ivector index;
658
659	wk = new_dvector(size);
660	index = new_ivector(size);
661	aw = 1.0;
662	for (i = 0; i < size; i++) {
663	maxb = 0.0;
664	for (j = 0; j < size; j++) {
665	if (fabs(amat[i][j]) > maxb)
666	maxb = fabs(amat[i][j]);
667	}
668	if (maxb == 0.0) {
669	fprintf(stderr, "luequation: singular matrix\n");
670	exit(1);
671	}
672	wk[i] = 1.0 / maxb;
673	}
674	for (j = 0; j < size; j++) {
675	for (i = 0; i < j; i++) {
676	sum = amat[i][j];
677	for (k = 0; k < i; k++)
678	sum -= amat[i][k] * amat[k][j];
679	amat[i][j] = sum;
680	}
681	maxb = 0.0;
682	for (i = j; i < size; i++) {
683	sum = amat[i][j];
684	for (k = 0; k < j; k++)
685	sum -= amat[i][k] * amat[k][j];
686	amat[i][j] = sum;
687	tmp = wk[i] * fabs(sum);
688	if (tmp >= maxb) {
689	maxb = tmp;
690	maxi = i;
691	}
692	}
693	if (j != maxi) {
694	for (k = 0; k < size; k++) {
695	tmp = amat[maxi][k];
696	amat[maxi][k] = amat[j][k];
697	amat[j][k] = tmp;
698	}
699	aw = -aw;
700	wk[maxi] = wk[j];
701	}
702	index[j] = maxi;
703	if (amat[j][j] == 0.0)
704	amat[j][j] = eps;
705	if (j != size - 1) {
706	tmp = 1.0 / amat[j][j];
707	for (i = j + 1; i < size; i++)
708	amat[i][j] *= tmp;
709	}
710	}
711
712	l = -1;
713	for (i = 0; i < size; i++) {
714	idx = index[i];
715	sum = yvec[idx];
716	yvec[idx] = yvec[i];
717	if (l != -1) {
718	for (j = l; j < i; j++)
719	sum -= amat[i][j] * yvec[j];
720	} else if (sum != 0.0)
721	l = i;
722	yvec[i] = sum;
723	}
724	for (i = size - 1; i >= 0; i--) {
725	sum = yvec[i];
726	for (j = i + 1; j < size; j++)
727	sum -= amat[i][j] * yvec[j];
728	yvec[i] = sum / amat[i][i];
729	}
730	free_ivector(index);
731	free_dvector(wk);
732	} /_ luequation /
733
734
735	void
736	lslength(tr, distanvec, numspc)
737	Tree *tr;
738	dvector distanvec;
739	int numspc;
740	{
741	int i, i1, j, j1, j2, k, numibrnch, numbrnch, numpair;
742	double sum, leng, alllen, rss;
743	ivector pths;
744	dmatrix atmt, atamt;
745	Node ebp, ibp;
746
747	if (Debug) printf("numspc = %d\n", numspc);
748	numibrnch = Numibrnch;
749	numbrnch = numspc + numibrnch;
750	numpair = (numspc * (numspc - 1)) / 2;
751	atmt = new_dmatrix(numbrnch, numpair);
752	atamt = new_dmatrix(numbrnch, numbrnch);
753	ebp = tr->ebrnchp;
754	ibp = tr->ibrnchp;
755	if (Debug) {
756	putchar('\n');
757	for (i = 0; i < numspc; i++) {
758	for (j = 0; j < numspc; j++) printf("%2d",ebp[i]->paths[j]);
759	putchar('\n');
760	}
761	for (i = 0; i < numibrnch; i++) {
762	for (j = 0; j < numspc; j++) printf("%2d",ibp[i]->paths[j]);
763	putchar('\n');
764	}
765	}
766
767	for (i = 0; i < numspc; i++) {
768	for (j1 = 1, j = 0; j1 < numspc; j1++) {
769	if (j1 == i) {
770	for (j2 = 0; j2 < j1; j2++, j++) {
771	atmt[i][j] = 1.0;
772	}
773	} else {
774	for (j2 = 0; j2 < j1; j2++, j++) {
775	if (j2 == i)
776	atmt[i][j] = 1.0;
777	else
778	atmt[i][j] = 0.0;
779	}
780	}
781	}
782	}
783	for (i1 = 0, i = numspc; i1 < numibrnch; i1++, i++) {
784	pths = ibp[i1]->paths;
785	for (j1 = 1, j = 0; j1 < numspc; j1++) {
786	for (j2 = 0; j2 < j1; j2++, j++) {
787	if (pths[j1] != pths[j2])
788	atmt[i][j] = 1.0;
789	else
790	atmt[i][j] = 0.0;
791	}
792	}
793	}
794
795	if (Debug) {
796	putchar('\n');
797	for (i = 0; i < numpair; i++) {
798	for (j = 0; j < numbrnch; j++)
799	printf("%2.0f",atmt[j][i]); /* !? */
800	printf("%6.1f\n",distanvec[i]);
801	}
802	}
803	#ifdef DIST
804	putchar('\n');
805	for (i = 0; i < numpair; i++)
806	printf("%5.1f",distanvec[i]);
807	putchar('\n');
808	#endif
809	for (i = 0; i < numbrnch; i++) {
810	for (j = 0; j <= i; j++) {
811	for (k = 0, sum = 0.0; k < numpair; k++)
812	sum += atmt[i][k] * atmt[j][k];
813	atamt[i][j] = sum;
814	atamt[j][i] = sum;
815	}
816	}
817	for (i = 0; i < numbrnch; i++) {
818	for (k = 0, sum = 0.0; k < numpair; k++)
819	sum += atmt[i][k] * distanvec[k];
820	Brnlength[i] = sum;
821	}
822	#if 0
823	putchar('\n');
824	for (i = 0; i < numbrnch; i++) {
825	for (j = 0; j < numbrnch; j++)
826	printf("%5.1f", atamt[i][j]);
827	printf("%7.3f",Brnlength[i]);
828	putchar('\n');
829	}
830	#endif
831	luequation(atamt, Brnlength, numbrnch);
832	#if 0
833	putchar('\n');
834	for (i = 0; i < numbrnch; i++) {
835	for (j = 0; j < numbrnch; j++)
836	printf("%5.1f", atamt[i][j]);
837	printf("%7.3f",Brnlength[i]);
838	putchar('\n');
839	}
840	#endif
841	if (Debug) {
842	for (i = 0; i < numbrnch; i++) {
843	printf("%7.3f",Brnlength[i]);
844	} putchar('\n');
845	}
846	for (i = 0, rss = 0.0; i < numpair; i++) {
847	sum = distanvec[i];
848	for (j = 0; j < numbrnch; j++) {
849	if (atmt[j][i] == 1.0 && Brnlength[j] > 0.0) /* caution */
850	sum -= Brnlength[j];
851	}
852	rss += sum * sum;
853	}
854	tr->rssleast = sqrt(rss); /* / numpair !? */
855	alllen = 0.0;
856	for (i = 0; i < numspc; i++) {
857	leng = Brnlength[i];
858	/* if (leng > 0.0) leng = 0.0 caution */
859	alllen += leng;
860	if (leng < Llimit) leng = Llimit;
861	if (leng > Ulimit) leng = Ulimit;
862	ebp[i]->length = leng;
863	ebp[i]->kinp->length = leng;
864	Brnlength[i] = leng;
865	}
866	for (i = 0, j = numspc; i < numibrnch; i++, j++) {
867	leng = Brnlength[j];
868	/* if (leng > 0.0) leng = 0.0 caution */
869	alllen += leng;
870	if (leng < Llimit) leng = Llimit;
871	if (leng > Ulimit) leng = Ulimit;
872	ibp[i]->length = leng;
873	ibp[i]->kinp->length = leng;
874	Brnlength[j] = leng;
875	}
876	tr->tbldis = alllen;
877	free_dmatrix(atmt);
878	free_dmatrix(atamt);
879	} /_ lslength /
880
881
882	void
883	slslength(tr, dmat, ns) /* simplify least square method (not equal LS) */
884	Tree *tr;
885	dmatrix dmat;
886	int ns;
887	{
888	int i, j, k, m, ne, ni, l, r, ll, rr;
889	int numibrnch, numbrnch;
890	double sumc, suml, sumr, leng, alllen, rss, coef;
891	ivector pths, lpaths, rpaths, npaths;
892	Node ebp, ibp, dp, ap, *xp;
893
894	lpaths = new_ivector(ns);
895	rpaths = new_ivector(ns);
896	npaths = new_ivector(ns);
897	numibrnch = Numibrnch;
898	numbrnch = ns + numibrnch;
899	ebp = tr->ebrnchp;
900	ibp = tr->ibrnchp;
901	alllen = 0.0;
902
903	#if SLSDEBUG
904	for (i = 0; i < ns; i++) {
905	for (j = 0; j < ns; j++) printf("%6.0f", dmat[i][j]*1000);
906	putchar('\n');
907	} putchar('\n');
908	for (i = 0; i < ns; i++) {
909	for (j = 0, sumc = 0.0; j < ns; j++) sumc += dmat[i][j];
910	printf("%6.0f", sumc*1000);
911	} putchar('\n'); putchar('\n');
912	#endif /* SLSDEBUG */
913
914	for ( ne = 0; ne < numbrnch; ne++) {
915
916	for (k = 0; k < ns; k++) {
917	npaths[k] = rpaths[k] = 0;
918	}
919	if (ne < ns) { /* internal */
920	l = 1;
921	ll = 1;
922	dp = ebp[ne];
923	ap = ebp[ne]->kinp;
924	pths = dp->paths;
925	for (k = 0; k < ns; k++) lpaths[k] = pths[k];
926	npaths[ne] = 1;
927	} else { /* external */
928	ni = ne - ns;
929	dp = ibp[ni];
930	ap = ibp[ni]->kinp;
931	pths = dp->paths;
932	for (k = 0, l = 0; k < ns; k++) {
933	lpaths[k] = 0;
934	if (pths[k]) l++;
935	}
936	for (xp = dp->isop, j = 1; xp != dp; xp = xp->isop, j++) {
937	pths = xp->paths;
938	for (k = 0, m = 0; k < ns; k++) if (pths[k]) m++;
939	if (xp->descen) {
940	m = ns - m;
941	for (k = 0; k < ns; k++) {
942	if (!pths[k]) {
943	lpaths[k] = j;
944	npaths[k] = m;
945	}
946	}
947	} else {
948	for (k = 0; k < ns; k++) {
949	if (pths[k]) {
950	lpaths[k] = j;
951	npaths[k] = m;
952	}
953	}
954	}
955	}
956	ll = j - 1;
957	}
958	r = ns - l;
959	for (xp = ap->isop, j = 1; xp != ap; xp = xp->isop, j++) {
960	pths = xp->paths;
961	for (k = 0, m = 0; k < ns; k++) if (pths[k]) m++;
962	if (xp->descen) {
963	m = ns - m;
964	for (k = 0; k < ns; k++) {
965	if (!pths[k]) {
966	rpaths[k] = j;
967	npaths[k] = m;
968	}
969	}
970	} else {
971	for (k = 0; k < ns; k++) {
972	if (pths[k]) {
973	rpaths[k] = j;
974	npaths[k] = m;
975	}
976	}
977	}
978	}
979	rr = j - 1;
980
981	#if SLSDEBUG
982	for (k = 0; k < ns; k++) printf("%2d", lpaths[k]); putchar('\n');
983	for (k = 0; k < ns; k++) printf("%2d", rpaths[k]); putchar('\n');
984	for (k = 0; k < ns; k++) printf("%2d", npaths[k]); putchar('\n');
985	#endif /* SLSDEBUG */
986
987	sumc = suml = sumr = 0.0;
988	for (i = 0; i < ns - 1; i++) {
989	for (j = i + 1; j < ns; j++) {
990	coef = (double)npaths[i] * (double)npaths[j];
991	if (lpaths[i]) {
992	if (rpaths[j]) {
993	sumc += dmat[i][j] / coef;
994	#if SLSDEBUG
995	printf("c1:%3d%3d%8.0f%7.0f\n",
996	i+1,j+1,sumc1000,dmat[i][j]/coef1000);
997	#endif /* SLSDEBUG */
998	} else if (lpaths[i] != lpaths[j]) {
999	suml += dmat[i][j] / coef;
1000	#if SLSDEBUG
1001	printf("l :%3d%3d%8.0f%7.0f\n",
1002	i+1,j+1,suml1000,dmat[i][j]/coef1000);
1003	#endif /* SLSDEBUG */
1004	}
1005	} else {
1006	if (!rpaths[j]) {
1007	sumc += dmat[i][j] / coef;
1008	#if SLSDEBUG
1009	printf("c2:%3d%3d%8.0f%7.0f\n",
1010	i+1,j+1,sumc1000,dmat[i][j]/coef1000);
1011	#endif /* SLSDEBUG */
1012	} else if (rpaths[i] != rpaths[j]) {
1013	sumr += dmat[i][j] / coef;
1014	#if SLSDEBUG
1015	printf("r :%3d%3d%8.0f%7.0f\n",
1016	i+1,j+1,sumr1000,dmat[i][j]/coef1000);
1017	#endif /* SLSDEBUG */
1018	}
1019	}
1020	}
1021	}
1022
1023	#if SLSDEBUG
1024	printf("%3s:%3d",(ne < ns ? "ext" : "int"), ne+1);
1025	printf(" %3d%3d%3d%3d", l, r, ll, rr);
1026	printf(" %5.2f%5.2f", (ne < ns ? 0 : (double)rr / (double)(ll-1)),
1027	(double)ll / (double)(rr-1));
1028	printf("%7.0f%7.0f%7.0f", sumc1000, suml1000, sumr*1000);
1029	#endif /* SLSDEBUG */
1030	if (ne >= ns) suml *= (double)rr / (double)(ll-1);
1031	sumr *= (double)ll / (double)(rr-1);
1032	leng = ( sumc - sumr - suml ) / (double)(ll * rr);
1033	#if SLSDEBUG
1034	printf(" leng:%9.5f\n", leng);
1035	#endif /* SLSDEBUG */
1036
1037	alllen += leng;
1038
1039	if (leng < 0.0) leng = 0.0; /* caution!? */
1040	if (ne < ns) /* external */
1041	pths = ebp[ne]->paths;
1042	else /* internal */
1043	pths = ibp[ni]->paths;
1044	for (i = 1; i < ns ; i++) {
1045	for (j = 0; j < i; j++) {
1046	if (pths[i] != pths[j]) dmat[i][j] -= leng;
1047	}
1048	}
1049
1050	if (leng > Ulimit) leng = Ulimit;
1051	if (ne < ns) { /* external */
1052	if (leng < LOWERLIMIT) leng = LOWERLIMIT;
1053	ebp[ne]->length = ebp[ne]->kinp->length = leng;
1054	} else { /* internal */
1055	if (leng < Llimit) leng = Llimit;
1056	ibp[ni]->length = ibp[ni]->kinp->length = leng;
1057	}
1058	Brnlength[ne] = leng; /* !? */
1059
1060	} /* for ne */
1061
1062	for (i = 1, rss = 0.0; i < ns ; i++) {
1063	for (j = 0; j < i; j++) {
1064	rss += dmat[i][j] * dmat[i][j];
1065	#if SLSDEBUG
1066	printf("%3d%3d%9.3f\n", i+1, j+1, dmat[i][j]);
1067	#endif /* SLSDEBUG */
1068	}
1069	}
1070	tr->rssleast = sqrt(rss); /* / numpair !? */
1071	#if SLSDEBUG
1072	printf("rss: %9.3f%9.3f\n", rss, tr->rssleast);
1073	#endif /* SLSDEBUG */
1074	for (i = 1; i < ns ; i++) {
1075	for (j = 0; j < i; j++) dmat[i][j] = dmat[j][i];
1076	}
1077
1078	tr->tbldis = alllen;
1079	free_ivector(lpaths);
1080	free_ivector(rpaths);
1081	free_ivector(npaths);
1082	} /_ slslength /
1083
1084
1085	#define DEBUGFM 0
1086	#define LSLEN 0
1087	void
1088	fmlength(tr, dmat, ns)
1089	Tree *tr;
1090	dmatrix dmat;
1091	int ns;
1092	{
1093	int i, j, k, no, nx, ny, m;
1094	double sy, x, xav, alllen, rss;
1095	ivector otux, path;
1096	dvector xv, yv;
1097	dmatrix dism;
1098	Node *otup, cp, rp, dp, *xp;
1099	#if LSLEN
1100	ivector otun;
1101	#endif /* LSLEN */
1102
1103	dism = new_dmatrix(ns, ns);
1104	for (i = 0; i < ns; i++) {
1105	for (j = 0; j < ns; j++) dism[i][j] = dmat[i][j];
1106	}
1107	otux = new_ivector(ns);
1108	xv = new_dvector(ns);
1109	yv = new_dvector(ns);
1110	otup = (Node **)new_npvector(ns);
1111	#if LSLEN
1112	otun = new_ivector(ns);
1113	#endif /* LSLEN */
1114	for (i = 0; i < ns; i++) {
1115	otux[i] = 0;
1116	otup[i] = tr->ebrnchp[i]->kinp;
1117	#if LSLEN
1118	otun[i] = 1;
1119	#endif /* LSLEN */
1120	}
1121	alllen = 0.0;
1122
1123	no = ns;
1124	cp = rp = tr->rootp;
1125	do {
1126	cp = cp->isop->kinp;
1127	if (cp->isop == NULL) { /* external node */
1128	cp = cp->kinp; /* not descen */
1129	} else { /* internal node */
1130	if (!cp->descen) {
1131	dp = cp->kinp;
1132	m = ns;
1133	for (xp = dp->isop, nx = 0; xp != dp; xp = xp->isop, nx++) {
1134	i = xp->num;
1135	otux[i] = i;
1136	if (i < m) m = i;
1137	}
1138	ny = no - nx;
1139	for (xp = dp->isop, sy = 0; xp != dp; xp = xp->isop) {
1140	i = xp->num;
1141	for (j = 0, xv[i] = 0, yv[i] = 0; j < ns; j++) {
1142	if (otup[j] != NULL) {
1143	if (otux[j])
1144	xv[i] += dism[i][j];
1145	else
1146	yv[i] += dism[i][j];
1147	}
1148	}
1149	yv[i] /= ny;
1150	sy += yv[i];
1151	}
1152	#if DEBUGFM
1153	printf("\nint. baranch %3d", cp->num+1);
1154	printf(" m:%3d nx:%3d ny:%3d no:%3d\n", m+1, nx, ny, no);
1155	printf("%3s", "");
1156	for (j = 0; j < ns; j++) printf("%4d", j+1); putchar('\n');
1157	for (i = 0; i < ns; i++) {
1158	printf("%3d", i+1);
1159	for (j = 0; j < ns; j++) printf("%4.0f", dism[i][j] * 10);
1160	putchar('\n');
1161	}
1162	printf("\n%3s", "otp");
1163	for (j = 0; j < ns; j++) printf("%4d", otup[j] != NULL ? 1 : 0);
1164	putchar('\n');
1165	printf("%3s", "otx");
1166	for (j = 0; j < ns; j++) printf("%4d", otux[j]); putchar('\n');
1167	printf("%3s", "xv");
1168	for (j = 0; j < ns; j++) printf("%4.0f", xv[j]*10);
1169	putchar('\n');
1170	printf("%3s", "yv");
1171	for (j = 0; j < ns; j++) printf("%4.0f", yv[j]*10);
1172	putchar('\n');
1173	#endif /* DEBUGFM */
1174	for (xp = dp->isop, xav = 0; xp != dp; xp = xp->isop) {
1175	i = xp->num;
1176	x = (xv[i]+yv[i])/nx - (sy-yv[i])/(nx*(nx-1));
1177	alllen += x; /* caution!? */
1178	if (x < Llimit) x = Llimit;
1179	xp->length = xp->kinp->length = x;
1180	xav += x;
1181	Brnlength[xp->kinp->num] = x;
1182	#if DEBUGFM
1183	printf("b-len %3d%9.4f%9.4f%9.4f\n", xp->kinp->num+1
1184	x, (xv[i]+yv[i])/nx, (sy-yv[i])/(nx*(nx-1)));
1185	#endif /* DEBUGFM */
1186	if (i != m) {
1187	for (j = 0; j < ns; j++) {
1188	if (otup[j] != NULL && j != m) {
1189	dism[m][j] += dism[i][j];
1190	dism[j][m] += dism[j][i];
1191	}
1192	#if DEBUGFM
1193	dism[i][j] = dism[j][i] = 0.0;
1194	#endif /* DEBUGFM */
1195	}
1196	otup[i] = NULL;
1197	}
1198	otux[i] = 0;
1199	if (xp->kinp->isop != NULL) xp->num = xp->kinp->num;
1200	}
1201
1202	for (j = 0; j < ns; j++) {
1203	if (otup[j] != NULL && j != m) {
1204	dism[m][j] = (dism[m][j] - xav)/nx;
1205	dism[j][m] = (dism[j][m] - xav)/nx;
1206	}
1207	xv[j] = yv[j] = 0.0;
1208	}
1209	cp->num = m;
1210	no -= (nx - 1);
1211
1212	}
1213	}
1214	} while (cp != rp);
1215
1216	xp = rp; nx = 0;
1217	do {
1218	xp = xp->isop;
1219	i = xp->num;
1220	otux[i] = 1;
1221	nx++;
1222	} while (xp != rp);
1223	ny = no - nx;
1224	xp = rp; sy = 0.0;
1225	do {
1226	xp = xp->isop;
1227	i = xp->num;
1228	for (j = 0, xv[i] = 0; j < ns; j++) {
1229	if (otup[j] != NULL && otux[j]) xv[i] += dism[i][j];
1230	}
1231	sy += xv[i];
1232	} while (xp != rp);
1233	sy /= 2;
1234	#if DEBUGFM
1235	printf("root nx:%3d ny:%3d no:%3d\n", nx, ny, no);
1236	printf("%3s", "");
1237	for (j = 0; j < ns; j++) printf("%4d", j+1); putchar('\n');
1238	for (i = 0; i < ns; i++) {
1239	printf("%3d", i+1);
1240	for (j = 0; j < ns; j++) printf("%4.0f", dism[i][j] * 10);
1241	putchar('\n');
1242	}
1243	printf("\n%3s", "otp");
1244	for (j = 0; j < ns; j++) printf("%4d", otup[j] != NULL ? 1 : 0);
1245	putchar('\n');
1246	printf("%3s", "otx");
1247	for (j = 0; j < ns; j++) printf("%4d", otux[j]); putchar('\n');
1248	printf("%3s", "xv");
1249	for (j = 0; j < ns; j++) printf("%4.0f", xv[j] * 10); putchar('\n');
1250	printf("sy:%9.5f\n", sy);
1251	#endif /* DEBUGFM */
1252	xp = rp;
1253	do {
1254	xp = xp->isop;
1255	i = xp->num;
1256	x = ( xv[i] - (sy-xv[i])/(nx-2) ) / (nx-1);
1257	if (xp->kinp->isop == NULL) {
1258	if (x < LOWERLIMIT) x = LOWERLIMIT;
1259	} else {
1260	if (x < Llimit) x = Llimit;
1261	}
1262	xp->length = xp->kinp->length = x;
1263	alllen += x;
1264	Brnlength[xp->kinp->num] = x;
1265	#if DEBUGFM
1266	printf("b-len %3d%9.4f%9.4f%9.4f%9.4f\n",
1267	xp->kinp->unm+1, x, xv[i], sy-xv[i], (sy-xv[i])/(nx-2) );
1268	#endif /* DEBUGFM */
1269	if (xp->kinp->isop != NULL) xp->num = xp->kinp->num;
1270	} while (xp != rp);
1271
1272	#if DEBUGFM
1273	for (j = 0; j < Numbrnch; j++) printf("%3d %9.4f\n", j+1, Brnlength[j]);
1274	#endif /* DEBUGFM */
1275
1276	for (i = 0; i < ns-1; i++) {
1277	for (j = i+1; j < ns; j++) dism[i][j] = dmat[i][j];
1278	}
1279	for (i = 0; i < ns; i++) {
1280	x = tr->ebrnchp[i]->length;
1281	path = tr->ebrnchp[i]->paths;
1282	for (j = 0; j < ns; j++) {
1283	if (path[j]) {
1284	for (k = 0; k < ns; k++) {
1285	if (path[k] == 0) {
1286	if (j < k) {
1287	dism[j][k] -= x;
1288	#if 0
1289	printf("%3d%3d%3d%9.4f%9.4f%9.4f\n",
1290	i+1,j+1,k+1,x,dism[j][k],dmat[j][k]);
1291	#endif
1292	} else {
1293	dism[k][j] -= x;
1294	#if 0
1295	printf("%3d%3d%3d%9.4f%9.4f%9.4f\n",
1296	i+1,j+1,k+1,x,dism[k][j],dmat[k][j]);
1297	#endif
1298	}
1299	}
1300	}
1301	}
1302	}
1303	}
1304	for (i = 0; i < Numibrnch; i++) {
1305	x = tr->ibrnchp[i]->length;
1306	if (x == Llimit) x = 0.0; /* caution */
1307	path = tr->ibrnchp[i]->paths;
1308	for (j = 0; j < ns; j++) {
1309	if (path[j]) {
1310	for (k = 0; k < ns; k++) {
1311	if (path[k] == 0) {
1312	if (j < k) {
1313	dism[j][k] -= x;
1314	#if 0
1315	printf("%3d%3d%3d%9.4f%9.4f%9.4f\n",
1316	i+ns+1,j+1,k+1,x,dism[j][k],dmat[j][k]);
1317	#endif
1318	} else {
1319	dism[k][j] -= x;
1320	#if 0
1321	printf("%3d%3d%3d%9.4f%9.4f%9.4f\n",
1322	i+ns+1,j+1,k+1,x,dism[k][j],dmat[k][j]);
1323	#endif
1324	}
1325	}
1326	}
1327	}
1328	}
1329	}
1330	for (i = 0, rss = 0.0; i < ns-1; i++) {
1331	for (j = i+1; j < ns; j++) {
1332	rss += dism[i][j] * dism[i][j];
1333	#if 0
1334	printf("%3d%3d%9.4f%9.4f%9.4f\n",
1335	i+1, j+1, dism[i][j], dmat[i][j], rss);
1336	#endif
1337	}
1338	}
1339
1340	tr->rssleast = rss;
1341	tr->tbldis = alllen;
1342	free_dmatrix(dism);
1343	free_ivector(otux);
1344	free_dvector(xv);
1345	free_dvector(yv);
1346	free_npvector(otup);
1347	#if LSLEN
1348	free_ivector(otun);
1349	#endif /* LSLEN */
1350	} /_ fmlength /
1351
1352
1353	extern double probnormal();
1354
1355	void
1356	resulttree(tr)
1357	Tree *tr;
1358	{
1359	boolean llimit, rtrif;
1360	int n, ne;
1361	double reli;
1362	Node ep, ip;
1363
1364	if (Relitrif != NULL) rtrif = TRUE; else rtrif = FALSE;
1365	printf("\nNo.%-3d", Cnotree + 1); /* offset */
1366	printf("%9s%7s%5s", "ext.","branch","S.E.");
1367	printf("%6s%7s%5s", "int.","branch","S.E.");
1368	if (Relia_optn) printf(" %5s%7s%6s", "LBP ", "2nd ", "pair");
1369	/* if (Info_optn) printf("%14s", "Likelihood"); */
1370	putchar('\n');
1371	#if 0
1372	for (n = 0; n < 28; n++) putchar('-');
1373	for (n = 0; n < 5; n++) putchar(' ');
1374	for (n = 0; n < 16; n++) putchar('-'); putchar('\n');
1375	#endif
1376	llimit = FALSE;
1377	for (n = 0; n < Numspc; n++) {
1378	ep = tr->ebrnchp[n];
1379	ne = ep->num;
1380	fputid(stdout, Identif[ne], 10);
1381	fputs(" ", stdout);
1382	printf("%3d", ne + 1); /* offset */
1383	if (ep->length == LOWERLIMIT)
1384	/* printf("%7s %-5s", "lower", "limit"); */
1385	printf("%7.2f%6.4s", ep->length, "----");
1386	else if (ep->length == Ulimit)
1387	printf("%7s %-5s", "UPPER", "LIMIT");
1388	else
1389	printf("%7.2f%6.2f", ep->length, sqrt(ep->kinp->lklhdl));
1390	if (n < Numibrnch) {
1391	ip = tr->ibrnchp[n];
1392	printf("%5d", n + 1 + Numspc); /* offset */
1393	if (ip->length == Llimit) {
1394	printf("%7s %-5s", "lower", "limit");
1395	llimit = TRUE;
1396	} else if (ip->length == Ulimit)
1397	printf("%7s %-5s", "UPPER", "LIMIT");
1398	else
1399	printf("%7.2f%6.2f", ip->length, sqrt(ip->kinp->lklhdl));
1400	if (Relia_optn && Relistat[n] >= 0) {
1401	if (Reliprob[n][0] != 1.0)
1402	printf(" %.3f", Reliprob[n][0]);
1403	else
1404	printf(" %.1f ", Reliprob[n][0]);
1405	(Relistat[n] > 0) ? putchar('*') : putchar(' ');
1406	if (Reliprob[n][0] != 1.0)
1407	printf(" %.3f", Reliprob[n][1]);
1408	else
1409	printf(" %.1f ", Reliprob[n][1]);
1410	printf(" %3d&%-3d", Relinum[n][0]+1, Relinum[n][1]+1);
1411	}
1412	#if 0
1413	reli = probnormal( ip->length / sqrt(ip->kinp->lklhdl) );
1414	if (ip->length == Llimit)
1415	printf("%6s", "----");
1416	else
1417	printf("%6.2f", reli);
1418	#endif
1419	/* if (Info_optn) printf("%14.3f", ip->lklhdl); */
1420	putchar('\n');
1421	} else {
1422	if (n == Numspc - 3) {
1423	printf("%7s%11.2f", "TBL :",tr->tblength * 1.0);
1424	printf("%7s %d", "iter:", Numit);
1425	if (!Converg)
1426	printf("%s", " non convergence!");
1427	else if (Converg == 2)
1428	printf("%s", " just before convergence");
1429	putchar('\n');
1430	} else if (n == Numspc - 2) {
1431	printf("%7s%11.2f +- %.2f\n", "ln L:",
1432	tr->lklhd, sqrt(tr->varilkl));
1433	} else if (n == Numspc - 1) {
1434	printf("%7s%11.2f", "AIC :",tr->aic);
1435	if (llimit) printf("%14s %.3f", "lower limit:",Llimit);
1436	putchar('\n');
1437	} else {
1438	putchar('\n');
1439	}
1440	}
1441	}
1442	} /_ resulttree /
1443
1444
1445	void
1446	bootstrap(infotrs, lklptrn)
1447	Infotree *infotrs;
1448	LPMATRIX lklptrn;
1449	{
1450	int i, j, k, n, imax, maxtree, nsite, nptrn, same, allsame;
1451	double lklmax, coefrand;
1452	ivector addweight;
1453	dvector boots;
1454
1455	addweight = new_ivector(Numsite);
1456	for ( j = 0, k = 0; k < Numptrn; k++ ) {
1457	for ( i = 0, imax = Weight[k]; i < imax; i++ )
1458	addweight[j++] = k;
1459	}
1460	if (j != Numsite) {
1461	fputs("ERROR in bootstrap(), sum Weight != Numsite\n", stderr);
1462	exit(1);
1463	}
1464
1465	coefrand = (double)Numsite / ((double)RANDOM_MAX + 1.0);
1466	boots = new_dvector(Numtree);
1467	allsame = 0;
1468	for (n = 0; n < Numtree; n++) infotrs[n].bsprob = 0;
1469
1470	#if MLSITE
1471	for (n = 0; n < Numtree; n++) infotrs[n].mlsite = 0;
1472	for (k = 0; k < Numptrn; k++) {
1473	for (n = 1, maxtree = 0, lklmax = lklptrn[0][k]; n < Numtree; n++) {
1474	if (lklptrn[n][k] > lklmax) {
1475	lklmax = lklptrn[n][k];
1476	maxtree = n;
1477	}
1478	}
1479	infotrs[maxtree].mlsite += Weight[k];
1480	}
1481	#endif
1482
1483	for (i = 0; i < NUMBOOTS; i++) {
1484	if (Verbs_optn && i % 100 == 99) fprintf(stderr, " %d", i+1);
1485	for (n = 0; n < Numtree; n++)
1486	boots[n] = 0.0;
1487	for (k = 0; k < Numsite; k++) {
1488	nsite = (int)( coefrand * (double)rand() ); /* RANDOM */
1489	nptrn = addweight[nsite];
1490	for (n = 0; n < Numtree; n++)
1491	boots[n] += lklptrn[n][nptrn];
1492	}
1493	maxtree = 0;
1494	lklmax = boots[0];
1495	same = 0;
1496	for (n = 1; n < Numtree; n++) {
1497	if (boots[n] >= lklmax) {
1498	if (boots[n] > lklmax) {
1499	maxtree = n;
1500	lklmax = boots[n];
1501	same = 0;
1502	} else {
1503	same++;
1504	}
1505	}
1506	}
1507	allsame += same;
1508	infotrs[maxtree].bsprob++;
1509	if (Debug_optn && i < 10)
1510	printf("%3d%4d%10.2f%6d%6d\n", i, maxtree, lklmax, nptrn, nsite);
1511	}
1512	if (allsame > 0)
1513	printf("\nsame bootstrap likelihood occured %d times\n", allsame);
1514	free_dvector(boots);
1515	free_ivector(addweight);
1516	} /_ bootstrap /
1517
1518
1519	extern double probnormal();
1520	extern double uprobnormal();
1521
1522	void
1523	tabletree(infotrs, lklptrn)
1524	Infotree *infotrs;
1525	LPMATRIX lklptrn;
1526	{
1527	int i, k, maxi, n;
1528	double ldiff, suml1, suml2, sdlkl, nn1, z;
1529	Infotree *info;
1530	LPVECTOR mlklptrn, alklptrn;
1531
1532	maxi = 52;
1533	printf("\n%s %s %s", Prog_name, VERSION, Modelname);
1534	printf(" %d trees %d OTUs %d sites. %s\n\n",
1535	Numtree, Numspc, Numsite, Comment);
1536	printf("%4s%9s%11s%6s%6s%6s%10s",
1537	"Tree","ln L","Diff ln L","S.E.","#Para","AIC","Diff AIC");
1538	/* if (Mevol_optn) */
1539	printf("%7s", " TBL ");
1540	if (Boots_optn) printf("%8s", "RELL-BP");
1541	putchar('\n');
1542	for (i = 0; i < maxi; i++) putchar('-');
1543	/* if (Mevol_optn) */
1544	{ for (i = 0; i < 7; i++) putchar('-'); }
1545	if (Boots_optn) { for (i = 0; i < 8; i++) putchar('-'); }
1546	putchar('\n');
1547	mlklptrn = lklptrn[Maxlkltree];
1548	nn1 = (double)Numsite / (double)(Numsite-1);
1549	for (n = 0; n < Numtree; n++) {
1550	info = &infotrs[n];
1551	alklptrn = lklptrn[n];
1552	printf("%-4d%11.1f%8.1f", n+1, info->lklhd, info->lklhd - Maxlkl); /* offset */
1553	if(n == Maxlkltree) {
1554	fputs(" <-best",stdout);
1555	} else {
1556	for (suml1 = suml2 = 0.0, k = 0; k < Numptrn; k++) {
1557	ldiff = alklptrn[k] - mlklptrn[k];
1558	suml1 += ldiff * Weight[k];
1559	suml2 += ldiff * ldiff * Weight[k];
1560	}
1561	suml1 /= Numsite;
1562	sdlkl = sqrt( nn1 * (suml2 - suml1suml1Numsite) );
1563	printf("%7.1f", sdlkl);
1564	#if 0
1565	for (suml2 = 0.0, k = 0; k < Numptrn; k++) {
1566	ldiff = alklptrn[k] - mlklptrn[k] - suml1;
1567	suml2 += ldiff * ldiff * Weight[k];
1568	}
1569	sdlkl = sqrt( nn1 * suml2 );
1570	printf("%7.1f", sdlkl);
1571	#endif
1572	}
1573	printf("%5d%10.1f%7.1f", info->npara, info->aic, info->aic - Minaic);
1574	/* if (Mevol_optn) { */
1575	if (n == Mintbltree) {
1576	printf("%6s ", "ME");
1577	} else {
1578	printf("%7.1f", info->tblength - Mintbl);
1579	}
1580	/* } */
1581	if (Boots_optn) printf("%8.4f", (double)info->bsprob / NUMBOOTS);
1582	#if MLSITE
1583	if (Boots_optn) printf("%4d", info->mlsite);
1584	#endif
1585	#if 0
1586	if(n == Maxlkltree) {
1587	fputs(" Base",stdout);
1588	} else {
1589	z = (info->lklhd - Maxlkl) / sdlkl;
1590	printf(" %.3f", probnormal(z)); /* NORMAL */
1591	}
1592	#endif
1593	putchar('\n');
1594	}
1595	} /_ tabletree /
1596
1597
1598	void
1599	tableinfo(infotrs)
1600	Infotree *infotrs;
1601	{
1602	int n;
1603	Infotree *info;
1604
1605	putchar('\n');
1606	for (n = 0; n < Numtree; n++) {
1607	info = &infotrs[n];
1608	if (Boots_optn) printf("%.4f\t", (double) info->bsprob / NUMBOOTS);
1609	printf("%.1f\t%d\t", info->lklhd - Maxlkl, n+1); /* offset */
1610	puts(info->ltplgy);
1611	}
1612	} /_ tableinfo /
1613
1614	#if 0
1615	void
1616	rerootq(tr, numspc)
1617	Tree *tr;
1618	int numspc;
1619	{
1620	Node cp, rp, op, xp, *yp;
1621	int i;
1622
1623	for (i = 0; i < numspc; i++) {
1624	if (tr->bturn[i] == numspc - 1)
1625	break;
1626	}
1627	rp = tr->ebrnchp[i]->kinp;
1628	tr->rootp = rp;
1629	cp = rp;
1630	do {
1631	cp = cp->isop->kinp;
1632	if (cp->isop == NULL) { /* external node */
1633	cp->descen = TRUE;
1634	cp = cp->kinp;
1635	cp->descen = FALSE;
1636	} else { /* internal node */
1637	if (cp->descen != -1) {
1638	cp->descen = TRUE;
1639	cp->kinp->descen = -1;
1640	tr->ibrnchp[cp->num - Maxspc] = cp;
1641	} else {
1642	cp->descen = FALSE;
1643	}
1644	if (!cp->descen) {
1645	op = cp->kinp;
1646	xp = op->isop;
1647	yp = xp->isop;
1648	if (xp->num > yp->num) {
1649	op->isop = yp;
1650	yp->isop = xp;
1651	xp->isop = op;
1652	}
1653	cp->num = op->isop->num;
1654	xp->num = xp->kinp->num;
1655	yp->num = yp->kinp->num;
1656	}
1657	}
1658	} while (cp != rp);
1659	op = cp;
1660	xp = op->isop;
1661	yp = xp->isop;
1662	if (xp->num > yp->num) {
1663	op->isop = yp;
1664	yp->isop = xp;
1665	xp->isop = op;
1666	}
1667	xp->num = xp->kinp->num;
1668	yp->num = yp->kinp->num;
1669	op->num = op->kinp->num;
1670	} /* rerootq */
1671	#endif
1672
1673
1674	void
1675	outlklhd(lklptrn)
1676	LPMATRIX lklptrn;
1677	{
1678	int i, j, k, l, maxk;
1679
1680	#if 1
1681	fprintf(Lklfp, "%d %d ", Numtree, Numsite);
1682	header(Lklfp, &Maxspc, &Maxsite, &Comment);
1683	for (i = 0; i < Numtree; i++) {
1684	fprintf(Lklfp, "# %d %.1f %s;\n",
1685	i+1, Infotrees[i].lklhd, Infotrees[i].ltplgy);
1686	for (j = 0, l = 0; j < Numptrn; j++) {
1687	for (k = 0, maxk = Weight[j]; k < maxk; k++, l++) {
1688	fprintf(Lklfp, "%16.8e", lklptrn[i][j]);
1689	if ((l+1) % 5 == 0) fputc('\n', Lklfp);
1690	}
1691	}
1692	if (l % 5 != 0) fputc('\n', Lklfp);
1693	}
1694	#else
1695	for (j = 0, l = 0; j < Numptrn; j++) {
1696	for (k = 0, maxk = Weight[j]; k < maxk; k++, l++) {
1697	for (i = 1; i < Numtree; i++) {
1698	fprintf(Lklfp, "%15.8f", lklptrn[0][j] - lklptrn[i][j]);
1699	}
1700	fputc('\n', Lklfp);
1701	}
1702	}
1703	#endif
1704	} /* outlklhd */
1705
1706
1707	void
1708	putsortseq(tr)
1709	Tree *tr;
1710	{
1711	Node rp, cp;
1712	int i, j, k;
1713
1714	printf("%d %d %s\n", Maxspc, Maxsite, Comment);
1715	cp = rp = tr->rootp;
1716	do {
1717	cp = cp->isop->kinp;
1718	if (cp->isop == NULL) { /* external node */
1719	i = cp->num;
1720	printf("%s", Identif[i]);
1721	if (Sciname[i]) printf(" %s", Sciname[i]);
1722	putchar('\n');
1723	for (k = 0; k < Maxsite; k += 60) {
1724	for (j = k; j < k + 60 && j < Maxsite; j++)
1725	putchar(int2acid(Seqchar[i][j]));
1726	putchar('\n');
1727	}
1728	cp = cp->kinp;
1729	}
1730	} while (cp != rp);
1731	} /* putsortseq */

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

Download in other formats: