Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

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

Line
1	/*
2	* sltree.c Adachi, J. 1995.02.15
3	* Copyright (C) 1992-1995 J. Adachi & M. Hasegawa, All rights reserved.
4	*/
5
6	#define SLDEBUG 0
7
8	#include "protml.h"
9
10
11	Tree *
12	new_stree(maxspc, maxibrnch, numptrn, seqconint)
13	int maxspc, maxibrnch;
14	imatrix seqconint;
15	{
16	int n, i;
17	Tree *tr;
18	Node dp, up;
19
20	tr = (Tree *) malloc(sizeof(Tree));
21	if (tr == NULL) maerror("tr in new_tree().");
22	tr->ebrnchp = (Node *) malloc((unsigned)maxspc sizeof(Node *));
23	if (tr->ebrnchp == NULL) maerror("ebrnchp in new_tree().");
24	tr->ibrnchp = (Node *) malloc((unsigned)maxibrnch sizeof(Node *));
25	if (tr->ibrnchp == NULL) maerror("ibrnchp in new_tree().");
26	tr->bturn = new_ivector(maxspc);
27	for (n = 0; n < maxspc; n++) {
28	tr->bturn[n] = n;
29	dp = (Node *) malloc(sizeof(Node));
30	if (dp == NULL) maerror("dp in new_tree().");
31	up = (Node *) malloc(sizeof(Node));
32	if (up == NULL) maerror("up in new_tree().");
33	dp->isop = NULL;
34	up->isop = NULL;
35	dp->kinp = up;
36	up->kinp = dp;
37	dp->descen = TRUE;
38	up->descen = FALSE;
39	dp->num = n;
40	up->num = n;
41	dp->length = 0.0;
42	up->length = 0.0;
43	dp->lklhdl = 0.0;
44	up->lklhdl = 0.0;
45	dp->paths = new_ivector(maxspc);
46	up->paths = dp->paths;
47	for (i = 0; i < maxspc; i++) dp->paths[i] = 0;
48	dp->paths[n] = 1;
49	dp->eprob = seqconint[n];
50	up->eprob = NULL;
51	dp->iprob = NULL;
52	up->iprob = new_dmatrix(numptrn, Tpmradix);
53	tr->ebrnchp[n] = dp;
54	}
55	for (n = 0; n < maxspc - 1; n++) {
56	tr->ebrnchp[n]->kinp->isop = tr->ebrnchp[n + 1]->kinp;
57	}
58	tr->ebrnchp[maxspc - 1]->kinp->isop = tr->ebrnchp[0]->kinp;
59	for (n = 0; n < maxibrnch; n++) {
60	dp = (Node *) malloc(sizeof(Node));
61	if (dp == NULL) maerror("dp in new_tree().");
62	up = (Node *) malloc(sizeof(Node));
63	if (up == NULL) maerror("up in new_tree().");
64	dp->isop = NULL;
65	up->isop = NULL;
66	dp->kinp = up;
67	up->kinp = dp;
68	dp->descen = TRUE;
69	up->descen = FALSE;
70	dp->num = n + maxspc;
71	up->num = n + maxspc;
72	dp->length = 0.0;
73	up->length = 0.0;
74	dp->lklhdl = 0.0;
75	up->lklhdl = 0.0;
76	dp->paths = new_ivector(maxspc);
77	up->paths = dp->paths;
78	for (i = 0; i < maxspc; i++) dp->paths[i] = 0;
79	dp->eprob = NULL;
80	up->eprob = NULL;
81	dp->iprob = new_dmatrix(numptrn, Tpmradix);
82	up->iprob = new_dmatrix(numptrn, Tpmradix);
83	tr->ibrnchp[n] = dp;
84	}
85	tr->rootp = tr->ebrnchp[maxspc - 1]->kinp;
86
87	return tr;
88	} /_ new_stree /
89
90
91	Infosltree *
92	newinfosltrees(num, maxbrnch)
93	int num;
94	int maxbrnch;
95	{
96	int i;
97	Infosltree *info;
98	double *v;
99
100	info = (Infosltree )malloc((unsigned)num sizeof(Infosltree));
101	if (info == NULL) maerror("1 in newinfosltrees().");
102	v = (double )malloc((unsigned)(num maxbrnch) * sizeof(double));
103	if (v == NULL) maerror("2 in newinfosltrees().");
104	for (i = 0; i < num; i++, v+=maxbrnch)
105	info[i].lengths = v;
106
107	return info;
108	} /_ newinfosltrees /
109
110
111	void
112	insertbranch(ibp, np)
113	Node ibp, np;
114	{
115	np->isop = ibp->isop->isop;
116	ibp->isop->isop = np->kinp;
117	np->kinp->isop = ibp->isop;
118	ibp->isop = np;
119	} /* insertbranch */
120
121
122	void
123	deletebranch(ibp, np)
124	Node ibp, np;
125	{
126	np->kinp->isop->isop = np->isop;
127	ibp->isop = np->kinp->isop;
128	np->kinp->isop = NULL; /* redundancy */
129	np->isop = NULL; /* redundancy */
130	} /* deletebranch */
131
132
133	void
134	movebranch(jbp, ip)
135	Node jbp, ip;
136	{
137	Node *jp;
138
139	jp = jbp->isop;
140	jbp->isop = jp->isop;
141	jp->isop = ip->isop;
142	ip->isop = jp;
143	} /* movebranch */
144
145	void
146	removebranch(jbp, ip)
147	Node jbp, ip;
148	{
149	Node *jp;
150
151	jp = ip->isop;
152	ip->isop = jp->isop;
153	jp->isop = jbp->isop;
154	jbp->isop = jp;
155	} /* removebranch */
156
157
158	void
159	subpathing(np)
160	Node *np;
161	{
162	int i;
163	Node *cp;
164	boolean npi, cpi;
165
166	for (i = 0, npi = np->paths; i < Maxspc; i++) *npi++ = FALSE;
167	for (cp = np->isop; cp != np; cp = cp->isop) {
168	npi = np->paths;
169	cpi = cp->paths;
170	for (i = 0; i < Maxspc; i++, npi++, cpi++) {
171	if (cpi == TRUE) npi = *cpi;
172
173	}
174	}
175	/* for (i=0; i<Maxspc; i++) printf("%1d",np->paths[i]);
176	putchar('\n'); */
177	} /* subpathing */
178
179
180	void
181	copylength(tr, lengs)
182	Tree *tr;
183	dvector lengs;
184	{
185	int i, j;
186	Node **bp;
187
188	bp = tr->ebrnchp;
189	for (i = 0; i < Numspc; i++) {
190	bp[i]->length = lengs[i];
191	bp[i]->kinp->length = lengs[i];
192	}
193	bp = tr->ibrnchp;
194	for (j = 0, i = Numspc; j < Numibrnch; j++, i++) {
195	bp[j]->length = lengs[i];
196	bp[j]->kinp->length = lengs[i];
197	}
198	}
199
200
201	Node*
202	sdml(tr, op)
203	Tree *tr;
204	Node *op;
205	{
206	Node cp, kp, *rp;
207	int i, l, nconv, nconv2;
208	double eps, lendiff;
209
210	/* prtopology(tr); */
211	kp = op->kinp;
212	for (cp = op->isop; cp->isop != op; cp = cp->isop)
213	;
214	rp = cp;
215	do {
216	cp = cp->isop->kinp;
217	if (cp->isop == NULL) { /* external node */
218	/* printf("rmle %3d%3d\n", cp->num+1,cp->descen); */
219	cp = cp->kinp; /* not descen */
220	partelkl(cp);
221	} else { /* internal node */
222	if (cp->kinp->descen != 2) {
223	cp->descen = 2;
224	} else {
225	/* printf("rmli %3d%3d\n", cp->num+1,cp->descen); */
226	prodpart(cp->kinp->isop);
227	partilkl(cp);
228	cp->descen ? (cp->kinp->descen = FALSE)
229	: (cp->kinp->descen = TRUE);
230	}
231	}
232	} while (cp != rp);
233
234
235	prodpart(op->isop);
236	mlibranch(op, 0.1, 5);
237	#if SLDEBUG
238	printf("\n%2s", "");
239	for (i = 0; i < Numbrnch; i++) printf("%5d",i+1); putchar('\n');
240	#endif
241	kp->isop->kinp->descen = kp->isop->isop->kinp->descen = 2;
242
243	for (l = 0, nconv = 0; l < MAXIT; l++) {
244	if (l == 0) eps = 1.0;
245	else if (l == 1) eps = 0.5;
246	else eps = 0.1;
247	#if SLDEBUG
248	printf("%2d", l+1);
249	for (i = 0; i < Numspc; i++)
250	printf("%5.0f",tr->ebrnchp[i]->length*100);
251	for (i = 0; i < Numibrnch; i++)
252	printf("%5.0f",tr->ibrnchp[i]->length*100);
253	putchar('\n');
254	#endif
255	cp = rp;
256	do {
257	cp = cp->isop->kinp;
258	if (!(l == 0 && cp->kinp == op)) prodpart(cp->kinp->isop);
259	if (cp->isop == NULL) { /* external node */
260	cp = cp->kinp; /* not descen */
261	lendiff = cp->length;
262	mlebranch(cp, eps, 5);
263	lendiff = fabs(lendiff - cp->length);
264	lendiff < 0.1 ? (nconv++) : (nconv = 0);
265	/* printf("e%3d%9.3f%9.3f\n", cp->num+1,cp->length,lendiff); */
266	} else { /* internal node */
267	if (cp->descen == 1) {
268	partilkl(cp);
269	} else if (cp->descen == 2 \|\| !cp->descen) {
270	if (cp->descen == 2 ) cp = cp->kinp;
271	lendiff = cp->length;
272	mlibranch(cp, eps, 5);
273	lendiff = fabs(lendiff - cp->length);
274	lendiff < 0.1 ? (nconv++) : (nconv = 0);
275	/* printf("i%3d%9.3f%9.3f\n",
276	cp->num+1,cp->length,lendiff); */
277	}
278	}
279	} while (cp != rp);
280	if (nconv >= 3) break;
281	}
282	kp->isop->descen ?
283	(kp->isop->kinp->descen = 0) : (kp->isop->kinp->descen = 1);
284	kp->isop->isop->descen ?
285	(kp->isop->isop->kinp->descen = 0) : (kp->isop->isop->kinp->descen = 1);
286
287	nconv = nconv2 = 0;
288	Converg = FALSE;
289	for (l = 0, Numit = 1; l < MAXIT; l++, Numit++) {
290	if (l == 0) eps = 1.0;
291	else if (l == 1) eps = 0.5;
292	else if (l == 2) eps = 0.1;
293	else eps = REPSILON;
294	#if SLDEBUG
295	printf("%2d", l+1);
296	for (i = 0; i < Numspc; i++)
297	printf("%5.0f",tr->ebrnchp[i]->length*100);
298	for (i = 0; i < Numibrnch; i++)
299	printf("%5.0f",tr->ibrnchp[i]->length*100);
300	putchar('\n');
301	#endif
302	cp = rp;
303	do {
304	cp = cp->isop->kinp;
305	prodpart(cp->kinp->isop);
306	if (cp->isop == NULL) { /* external node */
307	/* if (Debug) printf("mle %3d%3d\n",cp->num+1,cp->descen); */
308	cp = cp->kinp; /* not descen */
309	lendiff = cp->length;
310	mlebranch(cp, eps, 5);
311	lendiff = fabs(lendiff - cp->length);
312	lendiff < REPSILON ? (nconv++) : (nconv = 0);
313	lendiff < 0.5 ? (nconv2++) : (nconv2 = 0);
314	} else { /* internal node */
315	/* if (Debug) printf("mli %3d%3d\n",cp->num+1,cp->descen); */
316	if (cp->descen) {
317	partilkl(cp);
318	} else {
319	lendiff = cp->length;
320	mlibranch(cp, eps, 5);
321	lendiff = fabs(lendiff - cp->length);
322	lendiff < REPSILON ? (nconv++) : (nconv = 0);
323	lendiff < 0.5 ? (nconv2++) : (nconv2 = 0);
324	}
325	}
326	if (nconv >= Numbrnch) goto convergence;
327	} while (cp != rp);
328	}
329	if (nconv2 >= Numbrnch) Converg = 2;
330	evallkl(cp);
331	return rp;
332
333	convergence:
334	Converg = TRUE;
335	evallkl(cp);
336	return cp;
337	} /_ sdml /
338
339
340	void
341	decomposition(tr, n, infosltrees)
342	Tree *tr;
343	int n;
344	Infosltree *infosltrees;
345	{
346	Node rp, np, ibp, jbp, ip, jp, ibp2, jbp2, maxibp, maxjbp, *op;
347	double maxlklhd, maxaprlkl, res, rate, minres;
348	Infosltree head, tail, xp, yp, *zp;
349	int i, npair, ntree;
350	dvector lengs;
351
352	/* putchar('\n'); */
353
354	rp = tr->rootp;
355	np = tr->ibrnchp[n]->kinp;
356	npair = 0;
357	minres = 1.0e+37;
358	head = &infosltrees[MAXSLBUF];
359	tail = head;
360	for(ibp = rp; ibp->isop != rp; ibp = ibp->isop) {
361	ibp2 = ibp;
362	ip = ibp->isop;
363	/* printf("ip:%3d\n", ip->num+1); */
364	insertbranch(ibp, np);
365	for(jbp = np; jbp != rp; jbp = jbp->isop) {
366	jbp2 = jbp;
367	jp = jbp->isop;
368	if (jbp->isop == rp) tr->rootp = jbp;
369	movebranch(jbp, ip);
370
371	subpathing(np->kinp);
372	/* pathing(tr); */
373	slslength(tr, Distanmat, Maxspc);
374	/* initpartlkl(tr);
375	aproxlkl(tr);
376	praproxlkl(tr);
377	putctopology(tr); */
378
379	res = tr->rssleast;
380	#if 0
381	printf("%4d%4d%4d%4d %.3f\n",
382	n, npair, ip->num+1, jp->num+1, res);
383	#endif
384	if (npair < MAXSLBUF \|\| res < tail->residual) {
385	if (res < minres) minres = res;
386	if (npair < MAXSLBUF) {
387	yp = &infosltrees[npair];
388	} else if (res < tail->residual) {
389	yp = tail;
390	tail = yp->up;
391	}
392	yp->residual = res;
393	yp->ibp = ibp2;
394	yp->jbp = jbp2;
395	lengs = yp->lengths;
396	for (i = 0; i < Numbrnch; i++) lengs[i] = Brnlength[i];
397	for (xp = tail; res < xp->residual; zp = xp, xp = xp->up)
398	;
399	yp->up = xp;
400	if (tail == xp)
401	tail = yp;
402	else
403	zp->up = yp;
404	}
405
406	removebranch(jbp, ip);
407	if (jbp->isop == rp) tr->rootp = rp;
408	if (jbp == rp) rp = jbp->isop;
409	npair++;
410	}
411	deletebranch(ibp, np);
412	}
413
414	if (npair < MAXSLBUF)
415	ntree = npair;
416	else
417	ntree = MAXSLBUF;
418	maxaprlkl = -1.0e+37;
419	for (xp = tail; xp != head; xp = xp->up) {
420	ntree--;
421	rate = xp->residual / minres; /* !? */
422	if (rate < SRSRATE \|\| ntree < MAXSLBUF / 2) {
423	tr->rssleast = rate;
424	ibp = xp->ibp;
425	ip = ibp->isop;
426	insertbranch(ibp, np);
427	jbp = xp->jbp;
428	jp = jbp->isop;
429	if (jbp->isop == rp) tr->rootp = jbp;
430	movebranch(jbp, ip);
431	subpathing(np->kinp);
432	/* pathing(tr); */
433	copylength(tr, xp->lengths);
434	/* lslength(tr, Distanvec, Maxspc); */
435	initpartlkl(tr);
436	aproxlkl(tr);
437	xp->lklaprox = tr->lklhd;
438	#if 0
439	printf("%4d%4d%4d ", ntree, ip->num+1, jp->num+1);
440	praproxlkl(tr);
441	putctopology(tr);
442	#endif
443	removebranch(jbp, ip);
444	if (jbp->isop == rp) tr->rootp = rp;
445	if (jbp == rp) rp = jbp->isop;
446	deletebranch(ibp, np);
447
448	if (tr->lklhd > maxaprlkl) maxaprlkl = tr->lklhd;
449	}
450	}
451
452	/* putchar('\n'); */
453
454	maxlklhd = -1.0e+37;
455	for (xp = tail; xp != head; xp = xp->up) {
456	rate = maxaprlkl - xp->lklaprox; /* !? */
457	if (rate < 20) {
458	tr->rssleast = rate;
459	ibp = xp->ibp;
460	ip = ibp->isop;
461	insertbranch(ibp, np);
462	jbp = xp->jbp;
463	jp = jbp->isop;
464	if (jbp->isop == rp) tr->rootp = jbp;
465	movebranch(jbp, ip);
466	subpathing(np->kinp);
467	/* pathing(tr); */
468	copylength(tr, xp->lengths);
469	/* lslength(tr, Distanvec, Maxspc); */
470	initpartlkl(tr);
471	op = (Node *)mlikelihood(tr);
472	mlvalue(tr, Infotrees);
473	xp->lklaprox = tr->lklhd;
474	#if 0
475	printf("%4d%4d\t", ip->num+1, jp->num+1);
476	praproxlkl(tr);
477	putctopology(tr);
478	#endif
479	/* prtopology(tr);
480	resulttree(tr); */
481	removebranch(jbp, ip);
482	if (jbp->isop == rp) tr->rootp = rp;
483	if (jbp == rp) rp = jbp->isop;
484	deletebranch(ibp, np);
485
486	if (tr->lklhd > maxlklhd) {
487	maxibp = xp->ibp;
488	maxjbp = xp->jbp;
489	maxlklhd = tr->lklhd;
490	}
491	}
492	}
493
494
495	ip = maxibp->isop;
496	insertbranch(maxibp, np);
497	jp = maxjbp->isop;
498	/* printf("%4d%4d%4d\t", n, ip->num + 1, jp->num + 1); */
499	if (maxjbp->isop == rp) tr->rootp = maxjbp;
500	movebranch(maxjbp, ip);
501	subpathing(np->kinp);
502	/* lslength(tr, Distanvec, Maxspc);
503	initpartlkl(tr);
504	op = (Node *)mlikelihood(tr);
505	mlvalue(tr, Infotrees);
506	praproxlkl(tr); */
507	#if 1
508	putctopology(tr);
509	#endif
510	} /* decomposition */
511
512
513	void
514	stardecomp(tr, maxibrnch)
515	Tree *tr;
516	int maxibrnch;
517	{
518	int i;
519	Infosltree *infosltrees;
520
521	infosltrees = (Infosltree *)newinfosltrees(MAXSLBUF + 1, Maxbrnch);
522	infosltrees[MAXSLBUF].lklaprox = -1.0e+36;
523	infosltrees[MAXSLBUF].residual = 0.0;
524	for (i = 0; i < maxibrnch; i++) {
525	Numibrnch++;
526	Numbrnch++;
527	decomposition(tr, i, infosltrees);
528	}
529	} /* stardecomp */
530
531
532	dcube
533	new_dcubesym(nrow, ncol)
534	int nrow;
535	int ncol;
536	/* memory allocate a double cube */
537	{
538	int i, j;
539	dcube c;
540
541	c = (dcube) malloc((unsigned)nrow * sizeof(dmatrix));
542	if (c == NULL) maerror("1 in dcubesym().");
543	c = (dmatrix) malloc((unsigned)(nrow nrow) * sizeof(dvector));
544	if (*c == NULL) maerror("2 in dcubesym().");
545	*c = (dvector) malloc((unsigned)(nrow(nrow+1)/2 * ncol) * sizeof(double));
546	if (**c == NULL) maerror("3 in dcubesym().");
547	for (j = 1; j < nrow; j++) c[0][j] = c[0][j-1] + ncol;
548	for (i = 1; i < nrow; i++) {
549	c[i] = c[i-1] + nrow;
550	for (j = 0; j < i; j++) c[i][j] = c[j][i];
551	c[i][i] = c[i-1][nrow-1] + ncol;
552	for (j = i + 1; j < nrow; j++) c[i][j] = c[i][j-1] + ncol;
553	}
554	return c;
555	}
556
557
558	void
559	free_dcubesym(c)
560	dcube c;
561	{
562	free((char ) *c);
563	free((char ) c);
564	free((char *) c);
565	}
566
567
568	void
569	ystardecomp(tr)
570	Tree *tr;
571	{
572	Node rp, np, op, ibp, jbp, ip, jp, ibp2, jbp2, maxibp, *maxjbp;
573	double lkl, maxlkl, minlkl, suml1, suml2, ldiff, sdlkl, z, minz, maxz, eps;
574	double lk, llk;
575	int i, j, k, l, ii, jj, npair, notu, notu2, nn1, maxi;
576	Node **nodevec;
577	ivector pairnum;
578	dvector pairlkl, pairrel, alklptrn, mlklptrn;
579	dmatrix lklmat, relmat, iprb, jprb;
580	dcube lklcube;
581
582	eps = 1.0; /* !? */
583	rp = tr->rootp;
584	for(ibp = rp, notu = 1; ibp->isop != rp; ibp = ibp->isop, notu++) ;
585	nodevec = (Node *) malloc((unsigned)notu sizeof(Node *));
586	if (nodevec == NULL) maerror("nodevec in xstardecomp().");
587	pairnum = new_ivector(notu);
588	pairlkl = new_dvector(notu);
589	pairrel = new_dvector(notu);
590	lklmat = new_dmatrix(notu, notu);
591	relmat = new_dmatrix(notu, notu);
592	lklcube = new_dcubesym(notu, Numptrn);
593	nn1 = (double)(Numsite / (Numsite-1));
594
595	while (notu > 3) {
596
597	/* printf("notu: %d\n", notu); */
598	/*
599	for (i = 0, k = 1; i < notu; i++) {
600	for (j = i; j < notu; j++) lklcube[i][j][Numptrn-1] = k++;
601	}
602	putchar('\n');
603	for (i = 0; i < notu; i++) {
604	for (j = 0; j < notu; j++) printf("%4.0f", lklcube[i][j][Numptrn-1]);
605	putchar('\n');
606	}
607	*/
608	Numibrnch++; Numbrnch++;
609	for (i = 0; i < notu; i++) pairlkl[i] = -1.0e+37;
610	rp = tr->rootp;
611	for(ip = rp->isop, i = 0; ip != rp; ip = ip->isop, i++) nodevec[i] = ip;
612	/* printf("i: %d\n", i); */
613	nodevec[i] = rp;
614	np = tr->ibrnchp[Numibrnch - 1]->kinp;
615	npair = 0;
616	maxlkl = -1.0e+37;
617	minlkl = 1.0e+37;
618
619	initpartlkl(tr);
620	if (notu < Numspc) {
621	Alklptrn = lklcube[0][1];
622	op = (Node *)mlikelihood(tr);
623	mlvalue(tr, Infotrees);
624	initpartlkl(tr);
625	}
626	regupartlkl(tr); /* !? */
627
628	rp = tr->rootp;
629	for(ibp = rp, ii = 0; ibp->isop != rp; ibp = ibp->isop, ii++) {
630	ibp2 = ibp;
631	ip = ibp->isop;
632	for(jbp = ip, jj = ii + 1; jbp != rp; jbp = jbp->isop, jj++) {
633	jbp2 = jbp;
634	jp = jbp->isop;
635	Alklptrn = lklcube[ii][jj];
636	iprb = ip->iprob;
637	jprb = jp->iprob;
638	for (k = 0, lkl = 0.0; k < Numptrn; k++) {
639	for (l = 0, lk = 0.0; l < Tpmradix; l++) {
640	/*
641	if (notu < Numspc) printf("%3d%3d %3d %3d %20.7e%20.7e\n",
642	ii,jj,k,l,iprb[k][l],jprb[k][l]);
643	*/
644	lk += Freqtpm[l] * iprb[k][l] * jprb[k][l];
645	}
646	/*
647	if (notu < Numspc)
648	printf("%3d%3d %3d %15.10f\n", ii,jj,k,lk);
649	*/
650	llk = log(lk);
651	Alklptrn[k] = llk;
652	lkl += llk * Weight[k];
653	}
654	/* lkl /= Numsite; */
655	lklmat[ii][jj] = lkl;
656	lklmat[jj][ii] = lkl;
657	if (lkl > maxlkl) maxlkl = lkl;
658	if (lkl < minlkl) minlkl = lkl;
659	if (lkl > pairlkl[ii]) { pairnum[ii] = jj; pairlkl[ii] = lkl; }
660	if (lkl > pairlkl[jj]) { pairnum[jj] = ii; pairlkl[jj] = lkl; }
661	/*
662	printf("%-4d%4d%3d %.3f\n", npair,ip->num+1,jp->num+1,lkl);
663	*/
664	npair++;
665	}
666	}
667
668	for (i = 0; i < notu; i++) {
669	mlklptrn = lklcube[i][pairnum[i]];
670	minz = 1.0e+37;
671	for (j = 0; j < notu; j++) {
672	if (i != j) {
673	if (j == pairnum[i]) {
674	relmat[i][j] = 0.0;
675	} else {
676	alklptrn = lklcube[i][j];
677	for (suml1 = suml2 = 0.0, k = 0; k < Numptrn; k++) {
678	ldiff = alklptrn[k] - mlklptrn[k];
679	suml1 += ldiff * Weight[k];
680	suml2 += ldiff * ldiff * Weight[k];
681	}
682	suml1 /= Numsite;
683	sdlkl = sqrt( nn1 * (suml2 - suml1suml1Numsite) );
684	z = (lklmat[i][pairnum[i]] - lklmat[i][j]) / sdlkl;
685	relmat[i][j] = z;
686	if (z < minz) minz = z;
687	}
688	}
689	}
690	pairrel[i] = minz;
691	}
692	printf("\n%2s", "");
693	for (j = 0; j < notu; j++) printf("%4d", nodevec[j]->num+1); putchar('\n');
694	/* printf("\n%2s", "");
695	for (j = 0; j < notu; j++) printf("%4d", j + 1); putchar('\n'); */
696	for (i = 0; i < notu; i++) {
697	printf("%-2d", nodevec[i]->num+1);
698	for (j = 0; j < notu; j++) {
699	i != j ? printf("%4.0f", lklmat[i][j]-minlkl) : printf("%4s", "");
700	} putchar('\n');
701	}
702	printf("\n%2s", "");
703	for (j = 0; j < notu; j++) printf("%4d", nodevec[j]->num+1); putchar('\n');
704	for (j = 0; j < notu; j++) {
705	printf("%-2d", nodevec[j]->num+1);
706	for (i = 0; i < notu; i++) {
707	if (j == i) {
708	printf("%4s", "");
709	} else if (j == pairnum[i]) {
710	if (i == pairnum[j])
711	printf("%4s", "ML");
712	else
713	printf("%4s", "ml");
714	} else {
715	printf("%4.1f", relmat[i][j]);
716	}
717	} putchar('\n');
718	}
719	printf("%2s", "");
720	for (i = 0; i < notu; i++) printf("%4.1f", pairrel[i]); putchar('\n');
721	Numibrnch--; Numbrnch--;
722
723	for (i = 0, maxz = 0.0; i < notu; i++) {
724	j = pairnum[i];
725	if ( i == pairnum[j] && i < j ) {
726	z = (pairrel[i] > pairrel[j] ? pairrel[j] : pairrel[i]);
727	printf("%3d", nodevec[i]->num+1);
728	printf("%3d", nodevec[j]->num+1);
729	printf("%6.2f%6.2f\n", pairrel[i], pairrel[j]);
730	if (z > maxz) {
731	maxz = z;
732	maxi = i;
733	}
734	}
735	}
736	if ( notu == 4 && (maxi == 1 \|\| maxi == 2) ) maxi = 0;
737	for (i = 0, notu2 = notu; i < notu2; i++) {
738	j = pairnum[i];
739	/* if ( i == pairnum[j] && i < j &&
740	((pairrel[i] > eps && pairrel[j] > eps) \|\|
741	(i == maxi && maxz < eps)) ) { */
742	if ( i == pairnum[j] && i<j && (pairrel[i]>eps && pairrel[j]>eps) ) {
743	Numibrnch++; Numbrnch++;
744	np = tr->ibrnchp[Numibrnch - 1]->kinp;
745	ip = nodevec[i];
746	rp = tr->rootp;
747	for(ibp = rp; ibp->isop != ip; ibp = ibp->isop) ;
748	insertbranch(ibp, np);
749	jp = nodevec[j];
750	for(jbp = np; jbp->isop != jp; jbp = jbp->isop) ;
751	/* printf("%3d%3d%4d%3d\t", i, j, ip->num+1, jp->num+1); */
752	if (jbp->isop == rp) tr->rootp = jbp;
753	movebranch(jbp, ip);
754	subpathing(np->kinp);
755	/* lslength(tr, Distanvec, Maxspc);
756	initpartlkl(tr);
757	op = (Node *)mlikelihood(tr);
758	mlvalue(tr, Infotrees);
759	praproxlkl(tr); */
760	putctopology(tr);
761	notu--;
762	if (notu == 3) break;
763	}
764	/* } */
765	}
766	if (notu == notu2) break;
767	if (notu != 3) { putchar('\n'); prtopology(tr); }
768	fflush(stdout);
769
770	} /* while */
771
772	free_dcubesym(lklcube);
773	free_dmatrix(relmat);
774	free_dmatrix(lklmat);
775	free_dvector(pairrel);
776	free_dvector(pairlkl);
777	free_ivector(pairnum);
778	free((Node *) nodevec);
779	} /* ystardecomp */
780
781
782	void
783	xstardecomp(tr)
784	Tree *tr;
785	{
786	Node rp, np, op, ibp, jbp, ip, jp, ibp2, jbp2, maxibp, *maxjbp;
787	double lkl, maxlkl, minlkl, suml1, suml2, ldiff, sdlkl, z, minz, maxz, eps;
788	int nstep, i, j, k, k0, ii, jj, npair, notu, notu2, nn1, maxi, pairord0;
789	Node **nodevec;
790	ivector pairnum, pairord;
791	dvector pairlkl, pairrel, alklptrn, mlklptrn, elenvec, ilenvec;
792	dmatrix lklmat, relmat;
793	dcube lklcube;
794
795	elenvec = new_dvector(Maxspc);
796	ilenvec = new_dvector(Maxibrnch);
797	eps = 0.5; /* !? */
798	rp = tr->rootp;
799	for(ibp = rp, notu = 1; ibp->isop != rp; ibp = ibp->isop, notu++) ;
800	nodevec = (Node *) malloc((unsigned)notu sizeof(Node *));
801	if (nodevec == NULL) maerror("nodevec in xstardecomp().");
802	pairnum = new_ivector(notu);
803	pairord = new_ivector(notu);
804	pairlkl = new_dvector(notu);
805	pairrel = new_dvector(notu);
806	lklmat = new_dmatrix(notu, notu);
807	relmat = new_dmatrix(notu, notu);
808	lklcube = new_dcubesym(notu, Numptrn);
809	nn1 = (double)(Numsite / (Numsite-1));
810	nstep = 0;
811
812	while (notu > 3) {
813
814	if (Verbs_optn) fprintf(stderr, "%d OTUs:", notu);
815	fmlength(tr, Distanmat, Maxspc);
816	initpartlkl(tr);
817	rp = (Node *)mlikelihood(tr);
818	mlvalue(tr, Infotrees);
819	if (Debug_optn) putctopology(tr);
820	printf("#%d\n", nstep);
821	prtopology(Ctree);
822	/* resulttree(Ctree); */
823
824	for (i = 0; i < Numspc; i++) elenvec[i] = tr->ebrnchp[i]->length;
825	for (i = 0; i < Numibrnch; i++) ilenvec[i] = tr->ibrnchp[i]->length;
826	Numibrnch++; Numbrnch++;
827	for (i = 0; i < notu; i++) pairlkl[i] = -1.0e+37;
828	rp = tr->rootp;
829	for(ip = rp->isop, i = 0; ip != rp; ip = ip->isop, i++) nodevec[i] = ip;
830	nodevec[i] = rp;
831	np = tr->ibrnchp[Numibrnch - 1]->kinp;
832	npair = 0;
833	maxlkl = -1.0e+37;
834	minlkl = 1.0e+37;
835	rp = tr->rootp;
836	for(ibp = rp, ii = 0; ibp->isop != rp; ibp = ibp->isop, ii++) {
837	ibp2 = ibp;
838	ip = ibp->isop;
839	insertbranch(ibp, np);
840	for(jbp = np, jj = ii + 1; jbp != rp; jbp = jbp->isop, jj++) {
841	if (Verbs_optn) fprintf(stderr, " %d", npair+1);
842	jbp2 = jbp;
843	jp = jbp->isop;
844	if (jbp->isop == rp) tr->rootp = jbp;
845	movebranch(jbp, ip);
846
847	subpathing(np->kinp);
848	for (i = 0; i < Numspc; i++) {
849	tr->ebrnchp[i]->length = elenvec[i];
850	tr->ebrnchp[i]->kinp->length = elenvec[i];
851	}
852	for (i = 0; i < Numibrnch; i++) {
853	tr->ibrnchp[i]->length = ilenvec[i];
854	tr->ibrnchp[i]->kinp->length = ilenvec[i];
855	}
856	ip->length = ip->kinp->length = ip->length * 0.9;
857	jp->length = jp->kinp->length = jp->length * 0.9;
858	np->length = np->kinp->length = 1.0;
859	Alklptrn = lklcube[ii][jj];
860	op = (Node *)sdml(tr, np);
861	mlvalue(tr, Infotrees);
862	/* putctopology(tr); */
863
864	lkl = tr->lklhd;
865	lklmat[ii][jj] = lkl;
866	lklmat[jj][ii] = lkl;
867	if (lkl > maxlkl) maxlkl = lkl;
868	if (lkl < minlkl) minlkl = lkl;
869	if (lkl > pairlkl[ii]) { pairnum[ii] = jj; pairlkl[ii] = lkl; }
870	if (lkl > pairlkl[jj]) { pairnum[jj] = ii; pairlkl[jj] = lkl; }
871	removebranch(jbp, ip);
872	if (jbp->isop == rp) tr->rootp = rp;
873	if (jbp == rp) rp = jbp->isop;
874	npair++;
875	}
876	deletebranch(ibp, np);
877	}
878
879	for (i = 0; i < notu; i++) {
880	mlklptrn = lklcube[i][pairnum[i]];
881	minz = 1.0e+37;
882	for (j = 0; j < notu; j++) {
883	if (i != j) {
884	if (j == pairnum[i]) {
885	relmat[i][j] = 0.0;
886	} else {
887	alklptrn = lklcube[i][j];
888	for (suml1 = suml2 = 0.0, k = 0; k < Numptrn; k++) {
889	ldiff = alklptrn[k] - mlklptrn[k];
890	suml1 += ldiff * Weight[k];
891	suml2 += ldiff * ldiff * Weight[k];
892	}
893	suml1 /= Numsite;
894	sdlkl = sqrt( nn1 * (suml2 - suml1suml1Numsite) );
895	z = (lklmat[i][pairnum[i]] - lklmat[i][j]) / sdlkl;
896	relmat[i][j] = z;
897	if (z < minz) minz = z;
898	}
899	}
900	}
901	pairrel[i] = minz;
902	}
903
904	if (Info_optn) {
905	printf("\n%2s", "");
906	for (j = 0; j < notu; j++) printf("%4d", nodevec[j]->num+1); putchar('\n');
907	for (i = 0; i < notu; i++) {
908	printf("%-2d", nodevec[i]->num+1);
909	for (j = 0; j < notu; j++) {
910	i != j ? printf("%4.0f", lklmat[i][j]-minlkl) : printf("%4s", "");
911	} putchar('\n');
912	}
913	} /* Info_optn */
914
915	if (Info_optn) {
916	printf("\n%2s", "");
917	for (j = 0; j < notu; j++) printf("%4d", nodevec[j]->num+1); putchar('\n');
918	for (j = 0; j < notu; j++) {
919	printf("%-2d", nodevec[j]->num+1);
920	for (i = 0; i < notu; i++) {
921	if (j == i) {
922	printf("%4s", "");
923	} else if (j == pairnum[i]) {
924	if (i == pairnum[j])
925	printf("%4s", "ML");
926	else
927	printf("%4s", "ml");
928	} else {
929	printf("%4.0f", relmat[i][j] * 10.0);
930	}
931	} putchar('\n');
932	}
933	printf("%2s", "");
934	for (i = 0; i < notu; i++) printf("%4.0f",pairrel[i]*10.0); putchar('\n');
935	} /* Info_optn */
936
937	Numibrnch--; Numbrnch--;
938
939	for (i = 0, maxz = 0.0, pairord0 = notu; i < notu; i++) {
940	j = pairnum[i];
941	if ( i == pairnum[j] && i < j ) {
942	z = (pairrel[i] > pairrel[j] ? pairrel[j] : pairrel[i]);
943	if (Info_optn) {
944	printf("%3d", nodevec[i]->num+1);
945	printf("%3d", nodevec[j]->num+1);
946	printf("%6.2f%6.2f\n", pairrel[i], pairrel[j]);
947	}
948	pairrel[i] = pairrel[j] = z;
949	if (z > maxz) {
950	maxz = z;
951	maxi = i;
952	}
953	if (pairord0 == notu) {
954	pairord0 = i;
955	pairord[i] = notu;
956	} else {
957	for (k = k0 = pairord0; k < notu; k0 = k, k = pairord[k]) {
958	if (z > pairrel[k]) {
959	if (k == pairord0) {
960	pairord[i] = k;
961	pairord0 = i;
962	} else {
963	pairord[i] = k;
964	pairord[k0] = i;
965	}
966	break;
967	}
968	}
969	if (k == notu) {
970	pairord[i] = k;
971	pairord[k0] = i;
972	}
973	}
974	}
975	}
976
977	/*
978	for (k = pairord0; k < notu; k = pairord[k]) {
979	printf("%3d%3d%9.3f\n", k+1, pairnum[k]+1, pairrel[k]);
980	}
981	printf("%3s%3s%9s%3d\n", "", "", "", pairord0+1);
982	for (i = 0; i < notu; i++) {
983	if (i < pairnum[i])
984	printf("%3d%3d%9.3f%3d\n", i+1,pairnum[i]+1,pairrel[i],pairord[i]+1);
985	}
986	*/
987	if ( notu == 4 && (maxi == 1 \|\| maxi == 2) ) maxi = 0;
988	for (i = pairord0, notu2 = notu; i < notu2; i = pairord[i]) {
989	j = pairnum[i];
990	if (pairrel[i] > eps \|\| i == pairord0) {
991	np = tr->ibrnchp[Numibrnch]->kinp;
992	Numibrnch++; Numbrnch++;
993	np->length = np->kinp->length = 10.0;
994	ip = nodevec[i];
995	rp = tr->rootp;
996	for(ibp = rp; ibp->isop != ip; ibp = ibp->isop) ;
997	insertbranch(ibp, np);
998	jp = nodevec[j];
999	for(jbp = np; jbp->isop != jp; jbp = jbp->isop) ;
1000	/* printf("%3d%3d%4d%3d\t", i, j, ip->num+1, jp->num+1); */
1001	if (jbp->isop == rp) tr->rootp = jbp;
1002	movebranch(jbp, ip);
1003	subpathing(np->kinp);
1004	/* putctopology(tr); */
1005	notu--;
1006	if (notu == 3) break;
1007	}
1008	}
1009	if (notu == notu2) break;
1010	fflush(stdout);
1011	if (Verbs_optn) fprintf(stderr, "\n", notu);
1012	nstep++;
1013
1014	} /* while */
1015
1016	reroot(tr, tr->ebrnchp[Maxspc-1]->kinp);
1017	free_dcubesym(lklcube);
1018	free_dmatrix(relmat);
1019	free_dmatrix(lklmat);
1020	free_dvector(pairrel);
1021	free_dvector(pairlkl);
1022	free_ivector(pairord);
1023	free_ivector(pairnum);
1024	free((Node *) nodevec);
1025	free_dvector(elenvec);
1026	free_dvector(ilenvec);
1027	} /* xstardecomp */

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

Download in other formats: