Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

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

Line
1	/*
2	* mlklhd.c Adachi, J. 1996.03.02
3	* Copyright (C) 1992-1996 J. Adachi & M. Hasegawa. All rights reserved.
4	*/
5
6	#include "protml.h"
7
8	#define RELIDEBUG 0
9	#define QRELIDEBUG 0
10	#define MLEDEBUG 0
11	#define MLEDEBUG2 0
12	#define PRBCHECK 0
13
14
15	double
16	probnormal(z) /* lower probability of normal distribution */
17	double z;
18	{
19	int i;
20	double z2, prev, p, t;
21
22	z2 = z * z;
23	t = p = z * exp(-0.5 * z2) / sqrt(2 * 3.14159265358979323846264);
24	for (i = 3; i < 200; i += 2) {
25	prev = p; t *= z2 / i; p += t;
26	if (p == prev) return 0.5 + p;
27	}
28	return (z > 0);
29	} /* probnormal */
30
31
32	double
33	uprobnormal(z) /* upper probability of normal distribution */
34	double z;
35	{
36	return 1 - probnormal(z);
37	} /* uprobnormal */
38
39
40	#ifdef DEBUG
41	static void
42	prprob(xprob)
43	dmatrix xprob;
44	{
45	int i, j;
46
47	for (i = 0; i < Numptrn; i++) {
48	for (j = 0; j < Tpmradix; j++) printf(" %3.0f", xprob[i][j]*100.0);
49	putchar('\n');
50	}
51	} /* prprob */
52	#endif
53
54
55	void
56	copypart1(op, cp)
57	Node op, cp;
58	{
59	/* op (o) opb <--- cpb (c) cp */
60	int i;
61	dvector opb, cpb;
62
63	opb = *(op->iprob);
64	cpb = *(cp->iprob);
65	for (i = 0; i < Numptrn; i++) {
66	opb++ = cpb++;
67	opb++ = cpb++;
68	opb++ = cpb++;
69	opb++ = cpb++;
70	#ifndef NUC
71	opb++ = cpb++;
72	opb++ = cpb++;
73	opb++ = cpb++;
74	opb++ = cpb++;
75	opb++ = cpb++;
76	opb++ = cpb++;
77	opb++ = cpb++;
78	opb++ = cpb++;
79	opb++ = cpb++;
80	opb++ = cpb++;
81	opb++ = cpb++;
82	opb++ = cpb++;
83	opb++ = cpb++;
84	opb++ = cpb++;
85	opb++ = cpb++;
86	opb++ = cpb++;
87	#endif /* NUC */
88	}
89	} /_ copypart1 /
90
91
92	void
93	prodpart1(op, cp)
94	Node op, cp;
95	{
96	/* op (o) opb <--- cpb (c) cp */
97	int i;
98	dvector opb, cpb;
99
100	opb = *(op->iprob);
101	cpb = *(cp->iprob);
102	for (i = 0; i < Numptrn; i++) {
103	opb++ = *cpb++;
104	opb++ = *cpb++;
105	opb++ = *cpb++;
106	opb++ = *cpb++;
107	#ifndef NUC
108	opb++ = *cpb++;
109	opb++ = *cpb++;
110	opb++ = *cpb++;
111	opb++ = *cpb++;
112	opb++ = *cpb++;
113	opb++ = *cpb++;
114	opb++ = *cpb++;
115	opb++ = *cpb++;
116	opb++ = *cpb++;
117	opb++ = *cpb++;
118	opb++ = *cpb++;
119	opb++ = *cpb++;
120	opb++ = *cpb++;
121	opb++ = *cpb++;
122	opb++ = *cpb++;
123	opb++ = *cpb++;
124	#endif /* NUC */
125	}
126	/* if (Debug) prprob(op->iprob); */
127	} /_ prodpart1 /
128
129
130	void
131	prodpart(op)
132	Node *op;
133	{
134	/* (c) cp
135	* opb <---- cpb
136	* op (o) \|
137	* \| (c) cp'
138	* -- cpb'
139	*
140	* (c) last
141	*/
142	Node *cp;
143	int i;
144	dvector opb, cpb;
145
146	cp = op;
147	while (cp->isop->isop != op) {
148	cp = cp->isop;
149	opb = *(op->iprob);
150	cpb = *(cp->iprob);
151	for (i = 0; i < Numptrn; i++) {
152	opb++ = *cpb++;
153	opb++ = *cpb++;
154	opb++ = *cpb++;
155	opb++ = *cpb++;
156	#ifndef NUC
157	opb++ = *cpb++;
158	opb++ = *cpb++;
159	opb++ = *cpb++;
160	opb++ = *cpb++;
161	opb++ = *cpb++;
162	opb++ = *cpb++;
163	opb++ = *cpb++;
164	opb++ = *cpb++;
165	opb++ = *cpb++;
166	opb++ = *cpb++;
167	opb++ = *cpb++;
168	opb++ = *cpb++;
169	opb++ = *cpb++;
170	opb++ = *cpb++;
171	opb++ = *cpb++;
172	opb++ = *cpb++;
173	#endif /* NUC */
174	}
175	}
176	/* if (Debug) prprob(op->iprob); */
177	} /_ prodpart /
178
179
180	void
181	partilkl(op)
182	Node *op;
183	{
184	/* (i) iprob
185	* op / /
186	* --(a)----------(d) /
187	* oprob <------------
188	*/
189	int i, k;
190	double sum;
191	dmattpmty tprob;
192	dmatrix oprob, cprob;
193	dvector opb, cpb, cpb2, tpb;
194
195	tprobmtrx(op->length, tprob);
196	oprob = op->iprob;
197	cprob = op->kinp->isop->iprob;
198	for (k = 0; k < Numptrn; k++) {
199	opb = oprob[k];
200	cpb = cprob[k];
201	tpb = *tprob;
202	#ifndef NUC
203	for (i = 0; i < Tpmradix; i++) {
204	cpb2 = cpb;
205	sum = tpb++ *cpb2++;
206	sum += tpb++ *cpb2++;
207	sum += tpb++ *cpb2++;
208	sum += tpb++ *cpb2++;
209	sum += tpb++ *cpb2++;
210	sum += tpb++ *cpb2++;
211	sum += tpb++ *cpb2++;
212	sum += tpb++ *cpb2++;
213	sum += tpb++ *cpb2++;
214	sum += tpb++ *cpb2++;
215	sum += tpb++ *cpb2++;
216	sum += tpb++ *cpb2++;
217	sum += tpb++ *cpb2++;
218	sum += tpb++ *cpb2++;
219	sum += tpb++ *cpb2++;
220	sum += tpb++ *cpb2++;
221	sum += tpb++ *cpb2++;
222	sum += tpb++ *cpb2++;
223	sum += tpb++ *cpb2++;
224	sum += tpb++ *cpb2++;
225	opb[i] = sum;
226	}
227	#else /* NUC */
228	opb[0] = tpb[ 0] * cpb[0] + tpb[ 1] * cpb[1]
229	+ tpb[ 2] * cpb[2] + tpb[ 3] * cpb[3];
230	opb[1] = tpb[ 4] * cpb[0] + tpb[ 5] * cpb[1]
231	+ tpb[ 6] * cpb[2] + tpb[ 7] * cpb[3];
232	opb[2] = tpb[ 8] * cpb[0] + tpb[ 9] * cpb[1]
233	+ tpb[10] * cpb[2] + tpb[11] * cpb[3];
234	opb[3] = tpb[12] * cpb[0] + tpb[13] * cpb[1]
235	+ tpb[14] * cpb[2] + tpb[15] * cpb[3];
236	#endif /* NUC */
237	}
238	/* if (Debug) prprob(oprob); */
239	} /_ partilkl /
240
241
242	void
243	partelkl(op)
244	Node *op;
245	{
246	/* op
247	* --(a)----------(d)
248	* oprob <---- dseqi
249	*/
250	int i, k;
251	dmattpmty tprob;
252	dmatrix oprob;
253	ivector dseqi;
254	dvector opb, tbp;
255
256	tprobmtrxt(op->length, tprob);
257	oprob = op->iprob;
258	dseqi = op->kinp->eprob;
259	for (k = 0; k < Numptrn; k++) {
260	opb = oprob[k];
261	if ((i = dseqi[k]) >= 0) {
262	tbp = tprob[i];
263	opb[ 0] = tbp[ 0];
264	opb[ 1] = tbp[ 1];
265	opb[ 2] = tbp[ 2];
266	opb[ 3] = tbp[ 3];
267	#ifndef NUC
268	opb[ 4] = tbp[ 4];
269	opb[ 5] = tbp[ 5];
270	opb[ 6] = tbp[ 6];
271	opb[ 7] = tbp[ 7];
272	opb[ 8] = tbp[ 8];
273	opb[ 9] = tbp[ 9];
274	opb[10] = tbp[10];
275	opb[11] = tbp[11];
276	opb[12] = tbp[12];
277	opb[13] = tbp[13];
278	opb[14] = tbp[14];
279	opb[15] = tbp[15];
280	opb[16] = tbp[16];
281	opb[17] = tbp[17];
282	opb[18] = tbp[18];
283	opb[19] = tbp[19];
284	#endif /* NUC */
285	} else {
286	opb[ 0] = 1.0;
287	opb[ 1] = 1.0;
288	opb[ 2] = 1.0;
289	opb[ 3] = 1.0;
290	#ifndef NUC
291	opb[ 4] = 1.0;
292	opb[ 5] = 1.0;
293	opb[ 6] = 1.0;
294	opb[ 7] = 1.0;
295	opb[ 8] = 1.0;
296	opb[ 9] = 1.0;
297	opb[10] = 1.0;
298	opb[11] = 1.0;
299	opb[12] = 1.0;
300	opb[13] = 1.0;
301	opb[14] = 1.0;
302	opb[15] = 1.0;
303	opb[16] = 1.0;
304	opb[17] = 1.0;
305	opb[18] = 1.0;
306	opb[19] = 1.0;
307	#endif /* NUC */
308	}
309	}
310	} /_ partelkl /
311
312
313	void
314	partelkl2(op)
315	Node *op;
316	{
317	/* op
318	* --(a)----------(d)
319	* oprob <---- dseqi
320	*/
321	int j, k;
322	dmattpmty tprob;
323	dmatrix oprob;
324	ivector dseqi;
325	dvector opb;
326
327	tprobmtrx(op->length, tprob);
328	oprob = op->iprob;
329	dseqi = op->kinp->eprob;
330	for (k = 0; k < Numptrn; k++) {
331	opb = oprob[k];
332	if ((j = dseqi[k]) >= 0) {
333	opb[ 0] = tprob[ 0][j];
334	opb[ 1] = tprob[ 1][j];
335	opb[ 2] = tprob[ 2][j];
336	opb[ 3] = tprob[ 3][j];
337	#ifndef NUC
338	opb[ 4] = tprob[ 4][j];
339	opb[ 5] = tprob[ 5][j];
340	opb[ 6] = tprob[ 6][j];
341	opb[ 7] = tprob[ 7][j];
342	opb[ 8] = tprob[ 8][j];
343	opb[ 9] = tprob[ 9][j];
344	opb[10] = tprob[10][j];
345	opb[11] = tprob[11][j];
346	opb[12] = tprob[12][j];
347	opb[13] = tprob[13][j];
348	opb[14] = tprob[14][j];
349	opb[15] = tprob[15][j];
350	opb[16] = tprob[16][j];
351	opb[17] = tprob[17][j];
352	opb[18] = tprob[18][j];
353	opb[19] = tprob[19][j];
354	#endif /* NUC */
355	} else {
356	opb[ 0] = 1.0;
357	opb[ 1] = 1.0;
358	opb[ 2] = 1.0;
359	opb[ 3] = 1.0;
360	#ifndef NUC
361	opb[ 4] = 1.0;
362	opb[ 5] = 1.0;
363	opb[ 6] = 1.0;
364	opb[ 7] = 1.0;
365	opb[ 8] = 1.0;
366	opb[ 9] = 1.0;
367	opb[10] = 1.0;
368	opb[11] = 1.0;
369	opb[12] = 1.0;
370	opb[13] = 1.0;
371	opb[14] = 1.0;
372	opb[15] = 1.0;
373	opb[16] = 1.0;
374	opb[17] = 1.0;
375	opb[18] = 1.0;
376	opb[19] = 1.0;
377	#endif /* NUC */
378	}
379	}
380	} /_ partelkl /
381
382
383	void
384	initpartlkl(tr)
385	Tree *tr;
386	{
387	Node cp, rp;
388
389	cp = rp = tr->rootp;
390	do {
391	cp = cp->isop->kinp;
392	if (cp->isop == NULL) { /* external node */
393	/* if (Debug) printf("mle %3d %2d\n", cp->num+1, cp->descen); */
394	cp = cp->kinp; /* not descen */
395	partelkl(cp);
396	} else { /* internal node */
397	/* if (Debug) printf("mli %3d %2d\n", cp->num+1, cp->descen); */
398	if (!cp->descen) {
399	prodpart(cp->kinp->isop);
400	partilkl(cp);
401	}
402	}
403	} while (cp != rp);
404	} /_ initpartlkl /
405
406
407	void
408	regupartlkl(tr)
409	Tree *tr;
410	{
411	/* work on star decomposition */
412	Node cp, rp;
413	int i, j, prmax;
414	double pr;
415	dvector prbcv, prbkv;
416	dmatrix prbc, prbk;
417	dmattpmty tprob;
418
419	cp = rp = tr->rootp;
420	do {
421	cp = cp->isop;
422	if (cp->kinp->isop != NULL) { /* internal node */
423	prbc = cp->iprob;
424	prbk = cp->kinp->isop->iprob;
425	printf("length = %12.8f, %12.8f\n", cp->length, cp->kinp->length);
426	tprobmtrx(cp->length, tprob);
427	for (i = 0; i < Numptrn; i++) {
428	prbcv = prbc[i];
429	prbkv = prbk[i];
430	for (j = 1, pr = prbkv[0], prmax = 0; j < Tpmradix; j++) {
431	if (prbkv[j] > pr) {
432	pr = prbkv[j];
433	prmax = j;
434	}
435	}
436	for (j = 0; j < Tpmradix; j++) {
437	prbcv[j] = tprob[j][prmax];
438	}
439	}
440	}
441	} while (cp != rp);
442	} /_ regupartlkl /
443
444
445	void
446	mlibranch(op, eps, nloop)
447	Node op; / op->kinp == descen */
448	double eps;
449	int nloop;
450	{
451	/* (i) cprob
452	* op dist /
453	* (a)----------(b)
454	* / oprob'
455	* (i)
456	* oprob
457	*/
458	boolean conv;
459	int i, j, k, it;
460	double sumlk, sumd1, sumd2, lkld1, lkld2, prod1, prod2, vari;
461	double dist, distold, distdiff, distpre, coef, slk, sd1, sd2;
462	/* double d, d1, d2, diff, dis; */
463	dmattpmty tprob, tdif1, tdif2;
464	dmatrix oprob, cprob;
465	dvector tpb, td1, td2, opb, cpb;
466	#ifdef NUC
467	double cpb0, cpb1, cpb2, cpb3;
468	#endif /* NUC */
469
470	oprob = op->isop->iprob;
471	cprob = op->kinp->isop->iprob;
472	distpre = dist = op->length;
473	/*
474	if (1) {
475	for (k = 0, diff = 0.0; k < Numptrn; k++) {
476	opb = oprob[k];
477	cpb = cprob[k];
478	for (i = 0, d1 = d2 = 0.0; i < Tpmradix; i++) {
479	for (j = 0; j < Tpmradix; j++) {
480	d1 += opb[i] * cpb[j];
481	if (i != j) d2 += Freqtpm[i] * opb[i] * cpb[j];
482	}
483	}
484	d = d2 / d1;
485	diff += d * Weight[k];
486	}
487	if (diff < 1) diff = 1.0;
488	if (diff >= Maxsite) diff = Maxsite - 1.0;
489	dis = -log(1.0 - diff / Maxsite)*100.0;
490	}
491	*/
492	for (it = 0, conv = FALSE, coef = INITCOEFMLE; it < nloop; it++) {
493	tdiffmtrx(dist, tprob, tdif1, tdif2);
494	lkld1 = lkld2 = 0.0;
495	for (k = 0; k < Numptrn; k++) {
496	sumlk = sumd1 = sumd2 = 0.0;
497	opb = oprob[k];
498	cpb = cprob[k];
499	#ifdef NUC
500	cpb0 = cpb[0]; cpb1 = cpb[1]; cpb2 = cpb[2]; cpb3 = cpb[3];
501	#endif /* NUC */
502	for (i = 0; i < Tpmradix; i++) {
503	tpb = tprob[i];
504	td1 = tdif1[i];
505	td2 = tdif2[i];
506	prod1 = Freqtpm[i] * opb[i];
507	#ifndef NUC
508	slk = sd1 = sd2 = 0.0;
509	for (j = 0; j < Tpmradix; j++) {
510	prod2 = cpb[j];
511	slk += prod2 * tpb[j];
512	sd1 += prod2 * td1[j];
513	sd2 += prod2 * td2[j];
514	}
515	sumlk += prod1 * slk;
516	sumd1 += prod1 * sd1;
517	sumd2 += prod1 * sd2;
518	#else /* NUC */
519	slk = cpb0tpb[0] + cpb1tpb[1] + cpb2tpb[2] + cpb3tpb[3];
520	sd1 = cpb0td1[0] + cpb1td1[1] + cpb2td1[2] + cpb3td1[3];
521	sd2 = cpb0td2[0] + cpb1td2[1] + cpb2td2[2] + cpb3td2[3];
522	sumlk += prod1 * slk;
523	sumd1 += prod1 * sd1;
524	sumd2 += prod1 * sd2;
525	#endif /* NUC */
526	}
527	sumd1 /= sumlk;
528	lkld1 += sumd1 * Weight[k];
529	lkld2 += (sumd2 / sumlk - sumd1 * sumd1) * Weight[k];
530	}
531
532	distold = dist;
533	distdiff = - (lkld1 / lkld2);
534	if (lkld1 > 0) { /* NR internal */
535	dist += distdiff;
536	coef = INITCOEFMLE;
537	conv = TRUE;
538	} else {
539	if (lkld2 < 0) {
540	if (dist + distdiff > dist * coef) {
541	dist += distdiff;
542	coef = INITCOEFMLE;
543	conv = TRUE;
544	} else {
545	if (dist < 1.0) {
546	dist = Llimit;
547	coef = INITCOEFMLE;
548	} else {
549	dist *= coef;
550	coef *= INITCOEFMLE;
551	}
552	conv = FALSE;
553	}
554	} else {
555	if (dist < 1.0) {
556	dist = Llimit;
557	coef = INITCOEFMLE;
558	} else {
559	dist *= coef;
560	coef *= INITCOEFMLE;
561	}
562	conv = FALSE;
563	}
564	}
565	if (fabs(lkld1) > 5.0) conv = FALSE;
566	if (dist == distold) conv = TRUE;
567	if (dist < Llimit) { dist = Llimit; conv = TRUE; }
568	if (dist > Ulimit) { dist = Ulimit; conv = TRUE; }
569	#if MLEDEBUG2
570	printf("mli%3d%3d %7.3f %7.3f %9.4f %9.4f %12.5f %9.3f\n",
571	op->num+1,it+1,dist,distold, dist-distold,distdiff,lkld1,lkld2);
572	#endif
573	if (conv && fabs(distold - dist) < eps) break;
574	}
575	op->kinp->length = dist;
576	op->length = dist;
577	vari = 1.0 / fabs(lkld2);
578	op->lklhdl = vari;
579	if (Debug)
580	printf("mli%4d%3d%8.3f%8.3f%12.7f%7.2f%14.6f%14.3f\n",
581	op->num+1,it+1,dist,distpre,distpre-dist, sqrt(vari),lkld1,lkld2);
582
583	if (distold != dist) tprobmtrx(dist, tprob);
584	oprob = op->iprob;
585	#ifdef NUC
586	tpb = tprob[0];
587	#endif /* NUC */
588	for (k = 0; k < Numptrn; k++) {
589	opb = oprob[k];
590	cpb = cprob[k];
591	#ifndef NUC
592	for (i = 0; i < Tpmradix; i++) {
593	tpb = tprob[i];
594	for (j = 0, sumlk = 0.0; j < Tpmradix; j++)
595	sumlk += tpb[j] * cpb[j];
596	opb[i] = sumlk;
597	}
598	#else /* NUC */
599	cpb0 = cpb[0]; cpb1 = cpb[1]; cpb2 = cpb[2]; cpb3 = cpb[3];
600	opb[0] = tpb[ 0]cpb0 + tpb[ 1]cpb1 + tpb[ 2]cpb2 + tpb[ 3]cpb3;
601	opb[1] = tpb[ 4]cpb0 + tpb[ 5]cpb1 + tpb[ 6]cpb2 + tpb[ 7]cpb3;
602	opb[2] = tpb[ 8]cpb0 + tpb[ 9]cpb1 + tpb[10]cpb2 + tpb[11]cpb3;
603	opb[3] = tpb[12]cpb0 + tpb[13]cpb1 + tpb[14]cpb2 + tpb[15]cpb3;
604	#endif /* NUC */
605	}
606	} /_ mlibranch /
607
608
609	void
610	mlebranch(op, eps, nloop)
611	Node op; / op->kinp == descen */
612	double eps;
613	int nloop;
614	{
615	/* op dist
616	* (a)----------(b)
617	* / oprob' dseqi
618	* (i)
619	* oprob
620	*/
621	boolean conv;
622	int i, j, k, it;
623	double sumlk, sumd1, sumd2, lkld1, lkld2, prod, vari;
624	double dist, distold, distdiff, distpre, coef;
625	dmattpmty tprob, tdif1, tdif2;
626	dmatrix oprob;
627	ivector dseqi;
628	dvector opb, tpb, td1, td2;
629	#ifdef NUC
630	double pn0, pn1, pn2, pn3;
631	#endif /* NUC */
632
633	oprob = op->isop->iprob;
634	dseqi = op->kinp->eprob;
635	distpre = dist = op->length;
636
637	for (it = 0, conv = FALSE, coef = INITCOEFMLE; it < nloop; it++) {
638	tdiffmtrx(dist, tprob, tdif1, tdif2);
639	lkld1 = lkld2 = 0.0;
640	for (k = 0; k < Numptrn; k++) {
641	if ((i = dseqi[k]) >= 0) {
642	opb = oprob[k];
643	tpb = tprob[i]; td1 = tdif1[i]; td2 = tdif2[i];
644	#ifndef NUC
645	sumlk = sumd1 = sumd2 = 0.0;
646	for (j = 0; j < Tpmradix; j++) {
647	prod = opb[j];
648	sumlk += prod * tpb[j];
649	sumd1 += prod * td1[j];
650	sumd2 += prod * td2[j];
651	}
652	#else /* NUC */
653	pn0 = opb[0];
654	pn1 = opb[1];
655	pn2 = opb[2];
656	pn3 = opb[3];
657	sumlk = pn0 * tpb[0] + pn1 * tpb[1]
658	+ pn2 * tpb[2] + pn3 * tpb[3];
659	sumd1 = pn0 * td1[0] + pn1 * td1[1]
660	+ pn2 * td1[2] + pn3 * td1[3];
661	sumd2 = pn0 * td2[0] + pn1 * td2[1]
662	+ pn2 * td2[2] + pn3 * td2[3];
663	#endif /* NUC */
664	sumd1 /= sumlk;
665	lkld1 += sumd1 * Weight[k];
666	lkld2 += (sumd2 / sumlk - sumd1 * sumd1) * Weight[k];
667	}
668	}
669
670	distold = dist;
671	distdiff = - (lkld1 / lkld2);
672	if (lkld1 > 0) { /* NR external */
673	dist += distdiff;
674	coef = INITCOEFMLE;
675	conv = TRUE;
676	} else {
677	if (lkld2 < 0) {
678	if (dist + distdiff > dist * coef) {
679	dist += distdiff;
680	coef = INITCOEFMLE;
681	conv = TRUE;
682	} else {
683	if (dist < 0.2) {
684	dist = LOWERLIMIT;
685	coef = INITCOEFMLE;
686	} else {
687	dist *= coef;
688	coef *= INITCOEFMLE;
689	}
690	conv = FALSE;
691	}
692	} else {
693	if (dist < 0.2) {
694	dist = LOWERLIMIT;
695	coef = INITCOEFMLE;
696	} else {
697	dist *= coef;
698	coef *= INITCOEFMLE;
699	}
700	conv = FALSE;
701	}
702	}
703	if (fabs(lkld1) > 5.0) conv = FALSE;
704	if (dist == distold) conv = TRUE;
705	if (dist < LOWERLIMIT) { dist = LOWERLIMIT; conv = TRUE; }
706	if (dist > Ulimit) { dist = Ulimit; conv = TRUE; }
707	#if MLEDEBUG2
708	printf("mle%3d%3d %7.3f %7.3f %9.4f %9.4f %12.5f %9.3f\n",
709	op->num+1,it+1,dist,distold, dist-distold,distdiff,lkld1,lkld2);
710	#endif
711	if (conv && fabs(distold - dist) < eps) break;
712	}
713	op->kinp->length = dist;
714	op->length = dist;
715	vari = 1.0 / fabs(lkld2);
716	op->lklhdl = vari;
717	if (Debug)
718	printf("mle%4d%3d%8.3f%8.3f%12.7f%7.2f%14.6f%14.3f\n",
719	op->num+1,it+1,dist,distpre,distpre-dist,sqrt(vari),lkld1,lkld2);
720
721	tprobmtrxt(dist, tprob); /* if (distold != dist) */
722	oprob = op->iprob;
723	for (k = 0; k < Numptrn; k++) {
724	opb = oprob[k];
725	#ifndef NUC
726	if ((i = dseqi[k]) >= 0) {
727	tpb = tprob[i];
728	for (j = 0; j < Tpmradix; j++)
729	opb[j] = tpb[j];
730	} else {
731	for (j = 0; j < Tpmradix; j++)
732	opb[j] = 1.0; /* !? */
733	}
734	#else /* NUC */
735	if ((i = dseqi[k]) >= 0) {
736	tpb = tprob[i];
737	opb[0] = tpb[0];
738	opb[1] = tpb[1];
739	opb[2] = tpb[2];
740	opb[3] = tpb[3];
741	} else {
742	opb[0] = 1.0;
743	opb[1] = 1.0;
744	opb[2] = 1.0;
745	opb[3] = 1.0;
746	}
747	#endif /* NUC */
748	}
749	} /_ mlebranch /
750
751
752	void
753	mlebranch2(op, eps, nloop)
754	Node op; / op->kinp == descen */
755	double eps;
756	int nloop;
757	{
758	/* op dist
759	* (a)----------(b)
760	* / oprob' dseqi
761	* (i)
762	* oprob
763	*/
764	boolean conv;
765	int i, j, k, it;
766	double sumlk, sumd1, sumd2, lkld1, lkld2, prod, vari;
767	double dist, distold, distdiff, distpre, coef;
768	dmattpmty tprob, tdif1, tdif2;
769	dmatrix oprob;
770	ivector dseqi;
771	dvector opb;
772	#ifdef NUC
773	dvector tpb, td1, td2;
774	double pn0, pn1, pn2, pn3;
775	#endif /* NUC */
776
777	oprob = op->isop->iprob;
778	dseqi = op->kinp->eprob;
779	distpre = dist = op->length;
780
781	for (it = 0, conv = FALSE, coef = INITCOEFMLE; it < nloop; it++) {
782	tdiffmtrx(dist, tprob, tdif1, tdif2);
783	#ifdef NUC
784	tpb = tprob; td1 = tdif1; td2 = *tdif2;
785	#endif /* NUC */
786	lkld1 = lkld2 = 0.0;
787	for (k = 0; k < Numptrn; k++) {
788	if ((j = dseqi[k]) >= 0) {
789	opb = oprob[k];
790	#ifndef NUC
791	sumlk = sumd1 = sumd2 = 0.0;
792	for (i = 0; i < Tpmradix; i++) {
793	prod = Freqtpm[i] * opb[i];
794	sumlk += prod * tprob[i][j];
795	sumd1 += prod * tdif1[i][j];
796	sumd2 += prod * tdif2[i][j];
797	}
798	#else /* NUC */
799	pn0 = Freqtpm[0] * opb[0];
800	pn1 = Freqtpm[1] * opb[1];
801	pn2 = Freqtpm[2] * opb[2];
802	pn3 = Freqtpm[3] * opb[3];
803	sumlk = pn0 * tpb[j ] + pn1 * tpb[j+ 4]
804	+ pn2 * tpb[j+8] + pn3 * tpb[j+12];
805	sumd1 = pn0 * td1[j ] + pn1 * td1[j+ 4]
806	+ pn2 * td1[j+8] + pn3 * td1[j+12];
807	sumd2 = pn0 * td2[j ] + pn1 * td2[j+ 4]
808	+ pn2 * td2[j+8] + pn3 * td2[j+12];
809	#endif /* NUC */
810	sumd1 /= sumlk;
811	lkld1 += sumd1 * Weight[k];
812	lkld2 += (sumd2 / sumlk - sumd1 * sumd1) * Weight[k];
813	}
814	}
815
816	distold = dist;
817	distdiff = - (lkld1 / lkld2);
818	if (lkld1 > 0) { /* NR external */
819	dist += distdiff;
820	coef = INITCOEFMLE;
821	conv = TRUE;
822	} else {
823	if (lkld2 < 0) {
824	if (dist + distdiff > dist * coef) {
825	dist += distdiff;
826	coef = INITCOEFMLE;
827	conv = TRUE;
828	} else {
829	if (dist < 0.2) {
830	dist = LOWERLIMIT;
831	coef = INITCOEFMLE;
832	} else {
833	dist *= coef;
834	coef *= INITCOEFMLE;
835	}
836	conv = FALSE;
837	}
838	} else {
839	if (dist < 0.2) {
840	dist = LOWERLIMIT;
841	coef = INITCOEFMLE;
842	} else {
843	dist *= coef;
844	coef *= INITCOEFMLE;
845	}
846	conv = FALSE;
847	}
848	}
849	if (fabs(lkld1) > 5.0) conv = FALSE;
850	if (dist == distold) conv = TRUE;
851	if (dist < LOWERLIMIT) { dist = LOWERLIMIT; conv = TRUE; }
852	if (dist > Ulimit) { dist = Ulimit; conv = TRUE; }
853	#if MLEDEBUG2
854	printf("mle%3d%3d %7.3f %7.3f %9.4f %9.4f %12.5f %9.3f\n",
855	op->num+1,it+1,dist,distold, dist-distold,distdiff,lkld1,lkld2);
856	#endif
857	if (conv && fabs(distold - dist) < eps) break;
858	}
859	op->kinp->length = dist;
860	op->length = dist;
861	vari = 1.0 / fabs(lkld2);
862	op->lklhdl = vari;
863	if (Debug)
864	printf("mle%4d%3d%8.3f%8.3f%12.7f%7.2f%14.6f%14.3f\n",
865	op->num+1,it+1,dist,distpre,distpre-dist,sqrt(vari),lkld1,lkld2);
866
867	if (distold != dist) tprobmtrx(dist, tprob);
868	oprob = op->iprob;
869	#ifdef NUC
870	tpb = *tprob;
871	#endif /* NUC */
872	for (k = 0; k < Numptrn; k++) {
873	opb = oprob[k];
874	#ifndef NUC
875	if ((j = dseqi[k]) >= 0) {
876	for (i = 0; i < Tpmradix; i++)
877	opb[i] = tprob[i][j];
878	} else {
879	for (i = 0; i < Tpmradix; i++)
880	opb[i] = 1.0; /* !? */
881	}
882	#else /* NUC */
883	if ((j = dseqi[k]) >= 0) {
884	opb[0] = tpb[j];
885	opb[1] = tpb[j+4];
886	opb[2] = tpb[j+8];
887	opb[3] = tpb[j+12];
888	} else {
889	opb[0] = 1.0;
890	opb[1] = 1.0;
891	opb[2] = 1.0;
892	opb[3] = 1.0;
893	}
894	#endif /* NUC */
895	}
896	} /_ mlebranch /
897
898
899	void
900	evallkl(op)
901	Node *op;
902	{
903	/* op
904	* (a)----------( )
905	* / oprob
906	* (i)
907	* iprob
908	*/
909	int i, k;
910	double sumlk, lklhd;
911	dvector opb, ipb;
912	dmatrix oprob, iprob;
913
914	oprob = op->iprob;
915	iprob = op->isop->iprob;
916
917	for (k = 0, lklhd = 0.0; k < Numptrn; k++) {
918	opb = oprob[k];
919	ipb = iprob[k];
920	for (i = 0, sumlk = 0.0; i < Tpmradix; i++) {
921	sumlk += Freqtpm[i] * opb[i] * ipb[i];
922	}
923	sumlk = log(sumlk);
924	Alklptrn[k] = sumlk;
925	lklhd += sumlk * Weight[k];
926	/* printf("%3d %10.5f\n", k, sumlk); */
927	}
928	/* printf("branch:%3d ln L: %12.5f\n", op->kinp->num+1,lklhd); */
929	Ctree->lklhd = lklhd;
930
931	} /_ evallkl /
932
933
934	Node *
935	mlikelihood(tr)
936	Tree *tr;
937	{
938	Node cp, rp;
939	int l, nconv, nconv2;
940	double eps, lendiff;
941	#if MLEDEBUG
942	int i;
943	#endif
944
945	nconv = nconv2 = 0;
946	Converg = FALSE;
947	for (l = 0; l < MAXIT; l++) {
948	Numit = l + 1;
949	if (l == 0) eps = 1.0;
950	else if (l == 1) eps = 0.2;
951	else if (l == 2) eps = 0.1;
952	else eps = Epsilon;
953	if (Debug) printf("ml:%3d\n", l);
954	#if MLEDEBUG2
955	printf("ml:%3d\n", l+1);
956	#endif
957	#if MLEDEBUG
958	if (l == 0) putchar('\n');
959	printf("%2d", l+1);
960	for (i = 0; i < Numspc; i++)
961	printf("%5.0f",tr->ebrnchp[i]->length*100);
962	for (i = 0; i < Numibrnch; i++)
963	printf("%5.0f",tr->ibrnchp[i]->length*100);
964	printf(" %d\n", nconv);
965	#endif
966
967	cp = rp = tr->rootp;
968	do {
969	cp = cp->isop->kinp;
970	prodpart(cp->kinp->isop);
971	if (cp->isop == NULL) { /* external node */
972	/* if (Debug) printf("mle %3d%3d\n",cp->num+1,cp->descen); */
973	cp = cp->kinp; /* not descen */
974	lendiff = cp->length;
975	mlebranch(cp, eps, 5);
976	/* evallkl(cp); */
977	lendiff = fabs(lendiff - cp->length);
978	lendiff < Epsilon ? (nconv++) : (nconv = 0);
979	lendiff < 0.1 ? (nconv2++) : (nconv2 = 0);
980	} else { /* internal node */
981	/* if (Debug) printf("mli %3d%3d\n",cp->num+1,cp->descen); */
982	if (cp->descen) {
983	partilkl(cp);
984	/* mlibranch(cp, eps, 5); */
985	} else {
986	lendiff = cp->length;
987	mlibranch(cp, eps, 5);
988	/* evallkl(cp); */
989	lendiff = fabs(lendiff - cp->length);
990	lendiff < Epsilon ? (nconv++) : (nconv = 0);
991	lendiff < 0.1 ? (nconv2++) : (nconv2 = 0);
992	}
993	}
994	/* if (nconv >= Numbrnch) goto convergence; */
995	} while (cp != rp);
996	if (nconv >= Numbrnch) goto convergence;
997
998	}
999	if (nconv2 >= Numbrnch) Converg = 2;
1000	evallkl(cp);
1001	return rp;
1002
1003	convergence:
1004	Converg = TRUE;
1005	evallkl(cp);
1006	return cp;
1007	} /_ mlikelihood /
1008
1009
1010	void
1011	ribranch(op)
1012	Node op; / op->kinp == descen */
1013	{
1014	/* (i) cprob
1015	* op dist /
1016	* (a)----------(b)
1017	* / oprob'
1018	* (i)
1019	* oprob
1020	*/
1021	int i, j, k;
1022	double dist, sumlk, sumlk0, lsumlk, lsumlk0, slk, lklhd, lklhd0;
1023	double nn1, suml1, suml2, ldiff, sdlkl, dlkl, rel;
1024	dmattpmty tprob;
1025	dmatrix oprob, cprob;
1026	dvector tpb, opb, cpb;
1027	#ifdef NUC
1028	double cpb0, cpb1, cpb2, cpb3;
1029	#endif /* NUC */
1030
1031	oprob = op->isop->iprob;
1032	cprob = op->kinp->isop->iprob;
1033	dist = op->length;
1034
1035	tprobmtrx(dist, tprob);
1036	lklhd = lklhd0 = suml1 = suml2 = 0.0;
1037	for (k = 0; k < Numptrn; k++) {
1038	sumlk = sumlk0 = 0.0;
1039	opb = oprob[k];
1040	cpb = cprob[k];
1041	#ifdef NUC
1042	cpb0 = cpb[0]; cpb1 = cpb[1]; cpb2 = cpb[2]; cpb3 = cpb[3];
1043	#endif /* NUC */
1044	for (i = 0; i < Tpmradix; i++) {
1045	tpb = tprob[i];
1046	#ifndef NUC
1047	for (j = 0, slk = 0.0; j < Tpmradix; j++) {
1048	slk += cpb[j] * tpb[j];
1049	}
1050	#else /* NUC */
1051	slk = cpb0tpb[0] + cpb1tpb[1] + cpb2tpb[2] + cpb3tpb[3];
1052	#endif /* NUC */
1053	sumlk += Freqtpm[i] * opb[i] * slk;
1054	sumlk0 += Freqtpm[i] * opb[i] * cpb[i];
1055	}
1056	lsumlk = log(sumlk);
1057	lsumlk0 = log(sumlk0);
1058	lklhd += lsumlk * Weight[k];
1059	lklhd0 += lsumlk0 * Weight[k];
1060	ldiff = lsumlk - lsumlk0;
1061	suml1 += ldiff * Weight[k];
1062	suml2 += ldiff * ldiff * Weight[k];
1063	#if PRBCHECK
1064	if (sumlk < 0.0 \|\| sumlk0 < 0.0) {
1065	for (i = 0; i < Tpmradix; i++)
1066	printf("%3d %20.10e %20.10e\n", i, opb[i], cpb[i]);
1067	printf("%3d %15.3e %15.3e %9.3f %9.3f %9.3f\n",
1068	k, sumlk, sumlk0, lsumlk, lsumlk0, ldiff);
1069	}
1070	#endif
1071	}
1072	dlkl = lklhd - lklhd0;
1073	suml1 /= Numsite;
1074	nn1 = (double)Numsite / (double)(Numsite-1);
1075	sdlkl = sqrt( nn1 * (suml2 - suml1suml1Numsite) );
1076	Relitrif[op->num - Maxspc] = dlkl / sdlkl;
1077	rel = probnormal(dlkl / sdlkl);
1078	/*
1079	printf("%3d %12.3f %12.3f %9.3f %7.3f %7.3f\n",
1080	op->num+1, lklhd, lklhd0, dlkl, sdlkl, rel);
1081	*/
1082
1083
1084	oprob = op->iprob;
1085	#ifdef NUC
1086	tpb = tprob[0];
1087	#endif /* NUC */
1088	for (k = 0; k < Numptrn; k++) {
1089	opb = oprob[k];
1090	cpb = cprob[k];
1091	#ifndef NUC
1092	for (i = 0; i < Tpmradix; i++) {
1093	tpb = tprob[i];
1094	for (j = 0, sumlk = 0.0; j < Tpmradix; j++)
1095	sumlk += tpb[j] * cpb[j];
1096	opb[i] = sumlk;
1097	}
1098	#else /* NUC */
1099	cpb0 = cpb[0]; cpb1 = cpb[1]; cpb2 = cpb[2]; cpb3 = cpb[3];
1100	opb[0] = tpb[ 0]cpb0 + tpb[ 1]cpb1 + tpb[ 2]cpb2 + tpb[ 3]cpb3;
1101	opb[1] = tpb[ 4]cpb0 + tpb[ 5]cpb1 + tpb[ 6]cpb2 + tpb[ 7]cpb3;
1102	opb[2] = tpb[ 8]cpb0 + tpb[ 9]cpb1 + tpb[10]cpb2 + tpb[11]cpb3;
1103	opb[3] = tpb[12]cpb0 + tpb[13]cpb1 + tpb[14]cpb2 + tpb[15]cpb3;
1104	#endif /* NUC */
1105	}
1106	} /_ ribranch /
1107
1108
1109	Node *
1110	relibranch(op)
1111	Node *op;
1112	{
1113	Node *cp;
1114
1115	cp = op;
1116	do {
1117	cp = cp->isop->kinp;
1118	prodpart(cp->kinp->isop);
1119	if (cp->isop == NULL) { /* external node */
1120	cp = cp->kinp; /* not descen */
1121	partelkl(cp);
1122	} else { /* internal node */
1123	/* if (Debug) printf("mli %3d%3d\n",cp->num+1,cp->descen); */
1124	if (cp->descen) {
1125	partilkl(cp);
1126	} else {
1127	ribranch(cp);
1128	}
1129	}
1130	} while (cp != op);
1131	return op;
1132	} /_ relibranch /
1133
1134
1135	void
1136	mlvalue(tr, infotrs)
1137	Tree *tr;
1138	Infotree *infotrs;
1139	{
1140	int n, k, npara;
1141	double lklhd, tbl, aic, lkl, varilkl, ldiff;
1142
1143	#ifdef DIST
1144	int i, j, i1, i2;
1145	double dist;
1146	dmatrix amt;
1147	imatrix pths;
1148	amt = new_dmatrix(Numpair, Numbrnch);
1149	pths = tr->paths;
1150	for (i = 0, i1 = 0; i1 < (Numspc - 1); i1++) {
1151	for (i2 = i1 + 1; i2 < Numspc; i2++, i++) {
1152	for (j = 0; j < Numbrnch; j++) {
1153	pths[j][i1] != pths[j][i2] ? (amt[i][j]=1.0) : (amt[i][j]=0.0);
1154	}
1155	}
1156	}
1157	for (i = 0, i1 = 0; i1 < (Numspc - 1); i1++) {
1158	for (i2 = i1 + 1; i2 < Numspc; i2++, i++) {
1159	for (j = 0, dist = 0.0; j < Numbrnch; j++) {
1160	if (amt[i][j]) dist += tr->brnchp[j]->length;
1161	} printf("%5.1f",dist);
1162	}
1163	} putchar('\n');
1164	free_dmatrix(amt);
1165	#endif /* DIST */
1166
1167	tbl = 0.0;
1168	for (n = 0; n < Numspc; n++) tbl += tr->ebrnchp[n]->length;
1169	for (n = 0; n < Numibrnch; n++) tbl += tr->ibrnchp[n]->length;
1170
1171	npara = Numspc + Numibrnch;
1172	if (Frequ_optn) npara += Tpmradix - 1;
1173	#ifdef NUC
1174	if (AlphaBeta != 1.0) npara++;
1175	if (AlphaYR != 1.0) npara++;
1176	if (Beta12 != 1.0) npara++;
1177	#endif /* NUC */
1178
1179	lklhd = tr->lklhd;
1180	aic = npara * 2 - 2.0 * lklhd;
1181	lkl = lklhd / Numsite;
1182	for (k = 0, varilkl = 0.0; k < Numptrn; k++) {
1183	ldiff = Alklptrn[k] - lkl;
1184	varilkl += ldiff * ldiff * Weight[k];
1185	}
1186	tr->varilkl = varilkl;
1187	tr->npara = npara;
1188	tr->aic = aic;
1189	tr->tblength = tbl;
1190	infotrs[Cnotree].npara = npara;
1191	infotrs[Cnotree].lklhd = lklhd;
1192	infotrs[Cnotree].lklaprox = 0.0; /* !? */
1193	infotrs[Cnotree].aic = aic;
1194	infotrs[Cnotree].tblength = tbl;
1195	if (Cnotree == 0) {
1196	Maxlkltree = 0;
1197	Minaictree = 0;
1198	Mintbltree = 0;
1199	Maxlkl = lklhd;
1200	Minaic = aic;
1201	Mintbl = tbl;
1202	} else {
1203	if (lklhd > Maxlkl) {
1204	Maxlkltree = Cnotree;
1205	Maxlkl = lklhd;
1206	}
1207	if (aic < Minaic) {
1208	Minaictree = Cnotree;
1209	Minaic = aic;
1210	}
1211	if (tbl < Mintbl) {
1212	Mintbltree = Cnotree;
1213	Mintbl = tbl;
1214	}
1215	}
1216	return;
1217	} /_ mlvalue /
1218
1219
1220	void
1221	reroot(tr, rp)
1222	Tree *tr;
1223	Node *rp;
1224	{
1225	Node cp, op, xp, bp, *ap;
1226	boolean exch_flag;
1227
1228	tr->rootp = rp;
1229	cp = rp;
1230	do {
1231	cp = cp->isop->kinp;
1232	if (cp->isop == NULL) { /* external node */
1233	cp->descen = TRUE;
1234	cp = cp->kinp;
1235	cp->descen = FALSE;
1236	} else { /* internal node */
1237	if (cp->descen != -1) {
1238	cp->descen = TRUE;
1239	cp->kinp->descen = -1;
1240	tr->ibrnchp[cp->num - Maxspc] = cp;
1241	} else {
1242	/* printf("inode %3d\n", cp->kinp->num+1); */
1243	cp->descen = FALSE;
1244	cp = cp->kinp; /* reversed */
1245	op = cp;
1246	do {
1247	exch_flag = FALSE;
1248	/* printf("op:%3d\n", op->num+1); */
1249	for (bp = cp; bp->isop->isop != op; bp = bp->isop) {
1250	xp = bp->isop;
1251	ap = xp->isop;
1252	/* printf("bp:%3d xp:%3d ap:%3d\n",
1253	bp->num+1,xp->num+1,ap->num+1); */
1254	if (ap->num < xp->num) {
1255	xp->isop = ap->isop;
1256	ap->isop = xp;
1257	bp->isop = ap;
1258	exch_flag = TRUE;
1259	}
1260	}
1261	op = bp->isop;
1262	} while (exch_flag);
1263	cp->kinp->num = cp->isop->num;
1264	for (xp = cp->isop; xp != cp; xp = xp->isop) {
1265	xp->num = xp->kinp->num;
1266	}
1267	cp = cp->kinp; /* reversed */
1268	}
1269	}
1270	} while (cp != rp);
1271
1272	/* printf("root: %d\n", cp->kinp->num+1); */
1273	op = cp;
1274	do {
1275	exch_flag = FALSE;
1276	/* printf("op:%3d\n", op->num+1); */
1277	for (bp = cp; bp->isop->isop != op; bp = bp->isop) {
1278	xp = bp->isop;
1279	ap = xp->isop;
1280	/* printf("bp:%3d xp:%3d ap:%3d\n",
1281	bp->num+1,xp->num+1,ap->num+1); */
1282	if (ap->num < xp->num) {
1283	xp->isop = ap->isop;
1284	ap->isop = xp;
1285	bp->isop = ap;
1286	exch_flag = TRUE;
1287	}
1288	}
1289	op = bp->isop;
1290	} while (exch_flag);
1291
1292	cp->num = cp->kinp->num;
1293	for (xp = cp->isop; xp != cp; xp = xp->isop) {
1294	xp->num = xp->kinp->num;
1295	}
1296	} /* reroot */
1297
1298
1299	void
1300	sorttree(tr, rp)
1301	Tree *tr;
1302	Node *rp;
1303	{
1304	Node cp, op, xp, yp;
1305
1306	tr->rootp = rp;
1307	cp = rp;
1308	do {
1309	cp = cp->isop->kinp;
1310	if (cp->isop == NULL) { /* external node */
1311	cp->descen = TRUE;
1312	cp = cp->kinp;
1313	cp->descen = FALSE;
1314	cp->num = 1;
1315	} else { /* internal node */
1316	if (cp->descen != -1) {
1317	cp->descen = TRUE;
1318	cp->kinp->descen = -1;
1319	tr->ibrnchp[cp->num - Maxspc] = cp;
1320	} else {
1321	cp->descen = FALSE;
1322	}
1323	if (!cp->descen) {
1324	op = cp->kinp;
1325	xp = op->isop;
1326	yp = xp->isop;
1327	if (xp->num < yp->num) {
1328	op->isop = yp;
1329	yp->isop = xp;
1330	xp->isop = op;
1331	}
1332	cp->num = xp->num + yp->num;
1333	xp->num = xp->kinp->num;
1334	yp->num = yp->kinp->num;
1335	}
1336	}
1337	} while (cp != rp);
1338	op = cp;
1339	xp = op->isop;
1340	yp = xp->isop;
1341	if (xp->num < yp->num) {
1342	op->isop = yp;
1343	yp->isop = xp;
1344	xp->isop = op;
1345	}
1346	xp->num = xp->kinp->num;
1347	yp->num = yp->kinp->num;
1348	op->num = op->kinp->num;
1349	} /* sorttree */
1350
1351
1352	void
1353	chroot(tr, s1, s2)
1354	Tree *tr;
1355	int s1, s2;
1356	{
1357	Node rp, cp, op, xp, *yp;
1358
1359	if (Outgr_optn == 2) {
1360	} else {
1361	}
1362
1363	rp = tr->ebrnchp[s1]->kinp;
1364	cp = rp;
1365	do {
1366	cp = cp->isop->kinp;
1367	if (cp->isop == NULL) { /* external node */
1368	cp->descen = TRUE;
1369	cp = cp->kinp;
1370	cp->descen = FALSE;
1371	} else { /* internal node */
1372	if (cp->descen != -1) {
1373	cp->descen = TRUE;
1374	cp->kinp->descen = -1;
1375	tr->ibrnchp[cp->num - Maxspc] = cp;
1376	} else {
1377	cp->descen = FALSE;
1378	}
1379	if (!cp->descen) {
1380	op = cp->kinp;
1381	xp = op->isop;
1382	yp = xp->isop;
1383	if (xp->num > yp->num) {
1384	op->isop = yp;
1385	yp->isop = xp;
1386	xp->isop = op;
1387	}
1388	cp->num = op->isop->num;
1389	xp->num = xp->kinp->num;
1390	yp->num = yp->kinp->num;
1391	}
1392	}
1393	} while (cp != rp);
1394	op = cp;
1395	xp = op->isop;
1396	yp = xp->isop;
1397	if (xp->num > yp->num) {
1398	op->isop = yp;
1399	yp->isop = xp;
1400	xp->isop = op;
1401	}
1402	xp->num = xp->kinp->num;
1403	yp->num = yp->kinp->num;
1404	op->num = op->kinp->num;
1405	} /* chroot */
1406
1407
1408	void
1409	noexch(rp, exchstate)
1410	Node *rp;
1411	ivector exchstate;
1412	{
1413	Node *cp;
1414
1415	cp = rp->kinp;
1416	do {
1417	cp = cp->isop->kinp;
1418	if (cp->isop == NULL) { /* external node */
1419	cp = cp->kinp;
1420	} else { /* internal node */
1421	if (exchstate[cp->num - Maxspc] == 1) {
1422	exchstate[cp->num - Maxspc] = 2;
1423	} else if (exchstate[cp->num - Maxspc] == 2) {
1424	exchstate[cp->num - Maxspc] = 0;
1425	} else {
1426	cp = cp->kinp;
1427	}
1428	}
1429	} while (cp != rp);
1430	exchstate[cp->num - Maxspc] = 1;
1431	#if 0
1432	for (j = exchorder[i]; j < Numibrnch; j = exchorder[j]) {
1433	if (exchstate[j] && exchldiff[i] > exchldiff[j]) {
1434	cp = ibrn[j];
1435	if (rp == cp->isop->kinp \|\|
1436	rp == cp->isop->isop->kinp \|\|
1437	rp == cp->kinp->isop \|\|
1438	rp == cp->kinp->isop->kinp \|\|
1439	rp == cp->kinp->isop->isop \|\|
1440	rp == cp->kinp->isop->isop->kinp) {
1441	exchstate[j] = 0;
1442	}
1443	}
1444	}
1445	#endif
1446	} /* noexch */
1447
1448
1449	void
1450	reliml(tr, op, lklorg, mlklptrn, rel)
1451	Tree *tr;
1452	Node *op;
1453	double lklorg;
1454	LPVECTOR mlklptrn;
1455	double *rel;
1456	{
1457	Node cp, kp;
1458	int i, l, nconv, nconv2;
1459	double eps, lendiff, suml1, suml2, ldiff, sdlkl, lkldiff;
1460
1461	/* prtopology(tr); */
1462	kp = op->kinp;
1463	cp = op;
1464	do {
1465	cp = cp->isop->kinp;
1466	if (cp->isop == NULL) { /* external node */
1467	/* printf("rmle %3d%3d\n", cp->num+1,cp->descen); */
1468	cp = cp->kinp; /* not descen */
1469	partelkl(cp);
1470	} else { /* internal node */
1471	if (cp->kinp->descen != 2) {
1472	cp->descen = 2;
1473	} else {
1474	/* printf("rmli %3d%3d\n", cp->num+1,cp->descen); */
1475	prodpart(cp->kinp->isop);
1476	partilkl(cp);
1477	cp->descen ? (cp->kinp->descen = FALSE)
1478	: (cp->kinp->descen = TRUE);
1479	}
1480	}
1481	} while (cp != op);
1482
1483	prodpart(op->isop);
1484	mlibranch(op, 0.1, 5);
1485	#if RELIDEBUG
1486	printf("\n%3s", "");
1487	for (i = 0; i < Numbrnch; i++) printf("%5d",i+1); putchar('\n');
1488	#endif
1489	op->isop->kinp->descen = op->isop->isop->kinp->descen = 2;
1490	kp->isop->kinp->descen = kp->isop->isop->kinp->descen = 2;
1491
1492	for (l = 0, nconv = 0; l < MAXIT; l++) {
1493	if (l == 0) eps = 0.5;
1494	else eps = 0.1;
1495	#if RELIDEBUG
1496	printf("%3d", l+1);
1497	for (i = 0; i < Numspc; i++)
1498	printf("%5.0f",tr->ebrnchp[i]->length*100);
1499	for (i = 0; i < Numibrnch; i++)
1500	printf("%5.0f",tr->ibrnchp[i]->length*100);
1501	putchar('\n');
1502	#endif
1503	cp = op;
1504	do {
1505	cp = cp->isop->kinp;
1506	prodpart(cp->kinp->isop);
1507	if (cp->isop == NULL) { /* external node */
1508	cp = cp->kinp; /* not descen */
1509	lendiff = cp->length;
1510	mlebranch(cp, eps, 5);
1511	lendiff = fabs(lendiff - cp->length);
1512	lendiff < 0.1 ? (nconv++) : (nconv = 0);
1513	/* printf("e%3d%9.3f%9.3f\n", cp->num+1,cp->length,lendiff); */
1514	} else { /* internal node */
1515	if (!cp->descen && cp->kinp->descen)
1516	partilkl(cp);
1517	if (cp->descen == 2 \|\| (cp->descen && !cp->kinp->descen)) {
1518	if (cp->descen ==2 ) cp = cp->kinp;
1519	lendiff = cp->length;
1520	mlibranch(cp, eps, 5);
1521	lendiff = fabs(lendiff - cp->length);
1522	lendiff < 0.1 ? (nconv++) : (nconv = 0);
1523	/* printf("i%3d%9.3f%9.3f\n",
1524	cp->num+1,cp->length,lendiff); */
1525	}
1526	}
1527	} while (cp != op);
1528	if (nconv >= 5) break;
1529	}
1530	op->isop->descen ?
1531	(op->isop->kinp->descen = 0) : (op->isop->kinp->descen = 1);
1532	op->isop->isop->descen ?
1533	(op->isop->isop->kinp->descen = 0) : (op->isop->isop->kinp->descen = 1);
1534	kp->isop->descen ?
1535	(kp->isop->kinp->descen = 0) : (kp->isop->kinp->descen = 1);
1536	kp->isop->isop->descen ?
1537	(kp->isop->isop->kinp->descen = 0) : (kp->isop->isop->kinp->descen = 1);
1538
1539	nconv = nconv2 = 0;
1540	Converg = FALSE;
1541	for (l = 0, Numit = 1; l < MAXIT; l++, Numit++) {
1542	if (l == 0) eps = 0.5;
1543	else if (l == 1) eps = 0.1;
1544	else eps = Epsilon;
1545	#if RELIDEBUG
1546	printf("%3d", l+1);
1547	for (i = 0; i < Numspc; i++)
1548	printf("%5.0f",tr->ebrnchp[i]->length*100);
1549	for (i = 0; i < Numibrnch; i++)
1550	printf("%5.0f",tr->ibrnchp[i]->length*100);
1551	putchar('\n');
1552	#endif
1553	cp = op;
1554	do {
1555	cp = cp->isop->kinp;
1556	prodpart(cp->kinp->isop);
1557	if (cp->isop == NULL) { /* external node */
1558	/* if (Debug) printf("mle %3d%3d\n",cp->num+1,cp->descen); */
1559	cp = cp->kinp; /* not descen */
1560	lendiff = cp->length;
1561	mlebranch(cp, eps, 5);
1562	lendiff = fabs(lendiff - cp->length);
1563	lendiff < Epsilon ? (nconv++) : (nconv = 0);
1564	lendiff < 0.5 ? (nconv2++) : (nconv2 = 0);
1565	} else { /* internal node */
1566	/* if (Debug) printf("mli %3d%3d\n",cp->num+1,cp->descen); */
1567	if (cp->descen) {
1568	partilkl(cp);
1569	} else {
1570	lendiff = cp->length;
1571	mlibranch(cp, eps, 5);
1572	lendiff = fabs(lendiff - cp->length);
1573	lendiff < Epsilon ? (nconv++) : (nconv = 0);
1574	lendiff < 0.5 ? (nconv2++) : (nconv2 = 0);
1575	}
1576	}
1577	if (nconv >= Numbrnch) {
1578	Converg = TRUE;
1579	break;
1580	}
1581	} while (cp != op);
1582	if (Converg) break;
1583	}
1584	if (!Converg && nconv2 >= Numbrnch) Converg = 2;
1585	evallkl(cp);
1586
1587	for (i = 0, suml1 = suml2 = 0.0; i < Numptrn; i++) {
1588	ldiff = Alklptrn[i] - mlklptrn[i];
1589	suml1 += ldiff * Weight[i];
1590	suml2 += ldiff * ldiff * Weight[i];
1591	}
1592	suml1 /= Numsite;
1593	sdlkl = sqrt( (double)Numsite/(Numsite-1) * (suml2-suml1suml1Numsite) );
1594	lkldiff = lklorg - tr->lklhd;
1595	*rel = probnormal(lkldiff / sdlkl);
1596	#if 0
1597	printf("%3d", op->num+1);
1598	printf(" %7.3f %7.3f %7.3f %7.3f", lkldiff,sdlkl,lkldiff/sdlkl,*rel);
1599	printf(" %3s%3d\n",
1600	(Converg == TRUE ? "con" : (Converg == 2 ? "jbc" : "Noc")), Numit);
1601	#endif
1602	return;
1603	} /_ reliml /
1604
1605
1606	void
1607	localbp(reliprob, mlklptrn, rlklptrn, whichml, nb, ns)
1608	dmatrix reliprob;
1609	LPVECTOR mlklptrn;
1610	LPCUBE rlklptrn;
1611	ivector whichml;
1612	int nb, ns;
1613	{
1614	int i, j, k, ib, imax, nsite, nptrn;
1615	double coefrand;
1616	ivector addweight;
1617	imatrix bpn;
1618	dmatrix bpl;
1619
1620	if (Verbs_optn) fprintf(stderr, "bootstraping\n");
1621	addweight = new_ivector(ns);
1622	bpn = new_imatrix(nb, 3);
1623	bpl = new_dmatrix(nb, 3);
1624	coefrand = (double)Numsite / ((double)RANDOM_MAX + 1.0);
1625	for ( j = 0, k = 0; k < Numptrn; k++ ) {
1626	for ( i = 0, imax = Weight[k]; i < imax; i++ ) addweight[j++] = k;
1627	}
1628
1629	for (ib = 0; ib < nb; ib++)
1630	bpn[ib][0] = bpn[ib][1] = bpn[ib][2] = 0;
1631	for (i = 0; i < NUMBOOTSR; i++) {
1632	for (ib = 0; ib < nb; ib++)
1633	bpl[ib][0] = bpl[ib][1] = bpl[ib][2] = 0.0;
1634	for (k = 0; k < ns; k++) {
1635	nsite = (int)( coefrand * (double)rand() ); /* RANDOM */
1636	nptrn = addweight[nsite];
1637	for (ib = 0; ib < nb; ib++) {
1638	bpl[ib][0] += mlklptrn[nptrn];
1639	bpl[ib][1] += rlklptrn[ib][0][nptrn];
1640	bpl[ib][2] += rlklptrn[ib][1][nptrn];
1641	}
1642	}
1643	for (ib = 0; ib < nb; ib++) {
1644	bpl[ib][1] > bpl[ib][2]
1645	? ( bpl[ib][0] > bpl[ib][1] ? bpn[ib][0]++ : bpn[ib][1]++ )
1646	: ( bpl[ib][0] > bpl[ib][2] ? bpn[ib][0]++ : bpn[ib][2]++ );
1647	}
1648	}
1649	for (ib = 0; ib < nb; ib++) {
1650	reliprob[ib][0] = (double)bpn[ib][0] / (double)NUMBOOTSR;
1651	reliprob[ib][1] = (double)bpn[ib][whichml[ib]] / (double)NUMBOOTSR;
1652	}
1653
1654	free_imatrix(bpn);
1655	free_dmatrix(bpl);
1656	free_ivector(addweight);
1657	} /* localbp */
1658
1659
1660	void
1661	reliabranch(tr)
1662	Tree *tr;
1663	{
1664	boolean localconv, vibrate_flag;
1665	int ib, i, j, n, ll, forder, exn, exn0, exb, exb0;
1666	double lklorg, lklold, lklhd1, lklhd2, uldiff, muldiff, rel, minrel, iblen;
1667	Node rp, dp, ap, cp, cp0, cp1, kp, kp0, kp1, ncp;
1668	Node exchbrnch1, exchbrnch2, ebrn, ibrn;
1669	dvector elenvec, ilenvec, exchldiff, exchrelia;
1670	ivector exchstate, exchorder, whichml;
1671	LPVECTOR mlklptrn;
1672	LPCUBE rlklptrn;
1673
1674	/* root kp1 cp [1] [2]
1675	* R D ---(d) ap I dp (a)--- A X Y
1676	* (a)-----(d)
1677	* Z C ---(a) kp0 cp0 (a)--- B Y X
1678	* kp cp1
1679	*
1680	* [1]: A(X) <--> C(Z) , [2]: B(X) <--> C(Z)
1681	*/
1682	elenvec = new_dvector(Maxspc);
1683	ilenvec = new_dvector(Numibrnch);
1684	exchbrnch1 = (Node *)malloc((unsigned)Numibrnch sizeof(Node *));
1685	if (exchbrnch1 == NULL) maerror("exchbrnch1 in reliabranch().");
1686	exchbrnch2 = (Node *)malloc((unsigned)Numibrnch sizeof(Node *));
1687	if (exchbrnch2 == NULL) maerror("exchbrnch2 in reliabranch().");
1688	exchstate = new_ivector(Numibrnch);
1689	exchorder = new_ivector(Numibrnch);
1690	whichml = new_ivector(Numibrnch);
1691	exchrelia = new_dvector(Numibrnch);
1692	exchldiff = new_dvector(Numibrnch+1);
1693	exchldiff[Numibrnch] = - DBL_MAX;
1694	rlklptrn = NEW_LPCUBE(Numibrnch, 2, Numptrn);
1695	ebrn = tr->ebrnchp;
1696	ibrn = tr->ibrnchp;
1697	mlklptrn = Alklptrn;
1698	lklorg = tr->lklhd;
1699	exn = 0;
1700	exb = Numibrnch;
1701
1702	for (ll = 0, localconv = vibrate_flag = FALSE; !localconv; ll++) {
1703	if (ll) printf("%%%d\n", ll);
1704	if (Verbs_optn) fprintf(stderr, "%2d:", ll+1);
1705	for (i = 0; i < Numspc; i++) elenvec[i] = ebrn[i]->length;
1706	for (i = 0; i < Numibrnch; i++) ilenvec[i] = ibrn[i]->length;
1707
1708	for (ib = 0, forder = Numibrnch, muldiff = 0.0; ib < Numibrnch; ib++) {
1709	if (Verbs_optn) fprintf(stderr, " %d", ib+1);
1710	dp = ibrn[ib];
1711	kp0 = ap = dp->kinp;
1712	if (dp->isop->isop->isop != dp \|\| ap->isop->isop->isop != ap) {
1713	Relistat[ib] = -1;
1714	exchstate[ib] = 0;
1715	continue;
1716	}
1717	for (kp = ap->isop; kp->descen \|\| kp == tr->rootp;
1718	kp0 = kp, kp = kp->isop)
1719	;
1720	kp1 = kp->isop;
1721
1722	for (cp=dp->isop, cp0=dp, n=0; cp != dp; cp0=cp, cp=cp->isop, n++) {
1723	cp1 = cp->isop;
1724	kp0->isop = cp; cp->isop = kp1;
1725	cp0->isop = kp; kp->isop = cp1;
1726	for (i = 0; i < Numspc; i++)
1727	ebrn[i]->length = ebrn[i]->kinp->length = elenvec[i];
1728	for (i = 0; i < Numibrnch; i++)
1729	ibrn[i]->length = ibrn[i]->kinp->length = ilenvec[i];
1730	iblen = ap->length * 0.5;
1731	if (iblen < Llimit) iblen = Llimit;
1732	dp->length = ap->length = iblen;
1733	n == 0 ? (Alklptrn=rlklptrn[ib][0]):(Alklptrn=rlklptrn[ib][1]);
1734	reliml(tr, ap, lklorg, mlklptrn, &rel);
1735	n == 0 ? (lklhd1 = tr->lklhd) : (lklhd2 = tr->lklhd);
1736	if (n == 0) {
1737	Relinum[ib][0] = dp->isop->num;
1738	Relinum[ib][1] = dp->isop->isop->num;
1739	ncp = cp;
1740	minrel = rel;
1741	} else if (lklhd2 > lklhd1) {
1742	Relinum[ib][0] = dp->isop->num;
1743	Relinum[ib][1] = dp->isop->isop->num;
1744	ncp = cp;
1745	minrel = rel;
1746	}
1747	kp0->isop = kp; kp->isop = kp1;
1748	cp0->isop = cp; cp->isop = cp1;
1749	} /* for cp */
1750
1751	if (lklhd1 > lklhd2) {
1752	lklorg > lklhd1 ? (Relistat[ib] = 0) : (Relistat[ib] = 1);
1753	uldiff = lklhd1 - lklorg;
1754	whichml[ib] = 1;
1755	} else {
1756	lklorg > lklhd2 ? (Relistat[ib] = 0) : (Relistat[ib] = 2);
1757	uldiff = lklhd2 - lklorg;
1758	whichml[ib] = 2;
1759	}
1760	if (uldiff > 0.0 && uldiff < 0.00001 &&
1761	minrel < 0.5 && minrel > 0.499 &&
1762	fabs(lklhd1 - lklhd2) < 0.00001) { /* attention! */
1763	/* printf("%3d %20.15f %20.15f %20.15f\n",
1764	ib+Numspc+1, uldiff, minrel, fabs(lklhd1-lklhd2)); */
1765	uldiff = 0.0;
1766	minrel = 0.5;
1767	Relistat[ib] = 0;
1768	}
1769	Reliprob[ib][0] = minrel;
1770	Reliprob[ib][1] = 1.0 - minrel;
1771	uldiff > 0 ? (exchstate[ib] = 1) : (exchstate[ib] = 0);
1772	exchorder[ib] = -Numspc-1; /* -1 */
1773	exchldiff[ib] = uldiff;
1774	exchrelia[ib] = 1.0 - minrel;
1775	exchbrnch1[ib] = kp;
1776	exchbrnch2[ib] = ncp;
1777	if (uldiff > muldiff) {
1778	exchorder[ib] = forder;
1779	forder = ib;
1780	muldiff = uldiff;
1781	} else if (uldiff > 0.0) {
1782	for (i = forder; ; i = j) {
1783	j = exchorder[i];
1784	/* printf("exchorder[%3d]:%3d\n",i+1,j+1); */
1785	if (uldiff > exchldiff[j]) {
1786	exchorder[ib] = j;
1787	exchorder[i] = ib;
1788	break;
1789	}
1790	}
1791	}
1792	} /* for ib */
1793	if (Verbs_optn) fprintf(stderr, "\n");
1794	#if 0
1795	printf(" ib RS N1 N2 LBP1 LBP2 ES %3d ldiff relia\n",
1796	forder+Numspc+1);
1797	for (ib = 0; ib < Numibrnch; ib++) {
1798	printf("%3d%3d%3d%3d%6.3f%6.3f%3d%4d%15.10f%14.10f\n",
1799	ib+Numspc+1, Relistat[ib], Relinum[ib][0]+1, Relinum[ib][1]+1,
1800	Reliprob[ib][0], Reliprob[ib][1], exchstate[ib],
1801	exchorder[ib]+Numspc+1, exchldiff[ib], exchrelia[ib]);
1802	} putchar('\n');
1803	#endif
1804	for (i = forder; i < Numibrnch; i = exchorder[i]) {
1805	if (exchstate[i]) noexch(ibrn[i], exchstate);
1806	}
1807
1808	if (Xreli_optn) goto ONCE;
1809
1810	if (muldiff <= 0.0) {
1811	localconv = TRUE;
1812	} else {
1813	prtopology(tr); putchar('\n');
1814	exn0 = exn;
1815	exb0 = exb;
1816	exn = 0;
1817	for (i = forder; i < Numibrnch; i = exchorder[i]) {
1818	if (exchstate[i]) {
1819	exn++;
1820	exb = i;
1821	kp = exchbrnch1[i];
1822	cp = exchbrnch2[i];
1823	/*
1824	printf("%%%3d %3d<->%-3d ln L:%12.3f +%7.3f\n",i+Numspc+1,
1825	kp->kinp->num+1,cp->kinp->num+1,lklorg,exchldiff[i]);
1826	*/
1827	printf("%%%3d %3d<->%-3d ln L:%12.3f +%17.10f\n",i+Numspc+1,
1828	kp->kinp->num+1,cp->kinp->num+1,lklorg,exchldiff[i]);
1829	dp = ibrn[i];
1830	ap = dp->kinp;
1831	kp->isop->isop->isop = cp;
1832	cp->isop->isop->isop = kp;
1833	ncp = cp->isop;
1834	cp->isop = kp->isop;
1835	kp->isop = ncp;
1836	if (vibrate_flag) break;
1837	}
1838	}
1839	fflush(stdout);
1840	Alklptrn = mlklptrn;
1841	pathing(tr);
1842	slslength(tr, Distanmat, Maxspc);
1843	initpartlkl(tr);
1844	rp = (Node *)mlikelihood(tr);
1845	lklold = lklorg;
1846	lklorg = tr->lklhd;
1847	/* printf("%%%3d %3d %3d %3d %20.10f\n",
1848	exn, exn0, exb, exb0, lklorg - lklold); */
1849	if (exn == 1 && exn0 == 1 && exb == exb0) { /* attention! */
1850	if (fabs(lklorg - lklold) < 0.00001) localconv = TRUE;
1851	}
1852	if (lklorg < lklold) vibrate_flag = TRUE;
1853	}
1854	} /* for ll (!localconv) */
1855
1856	ONCE: /* if (Xreli_optn) */
1857
1858	Alklptrn = mlklptrn;
1859	for (i = 0; i < Numspc; i++)
1860	ebrn[i]->length = ebrn[i]->kinp->length = elenvec[i];
1861	for (i = 0; i < Numibrnch; i++)
1862	ibrn[i]->length = ibrn[i]->kinp->length = ilenvec[i];
1863	initpartlkl(tr);
1864	rp = (Node *)mlikelihood(tr);
1865	mlvalue(tr, Infotrees);
1866	prtopology(tr); putchar('\n');
1867	if (Verbs_optn) fputs("rerooting\n", stderr);
1868	reroot(tr, tr->rootp);
1869	putctopology(tr); putchar('\n');
1870	localbp(Reliprob, mlklptrn, rlklptrn, whichml, Numibrnch, Numsite);
1871	FREE_LPCUBE(rlklptrn);
1872	free_dvector(exchrelia);
1873	free_dvector(exchldiff);
1874	free_ivector(exchstate);
1875	free_ivector(exchorder);
1876	free_ivector(whichml);
1877	free(exchbrnch1);
1878	free(exchbrnch2);
1879	free_dvector(elenvec);
1880	free_dvector(ilenvec);
1881	} /_ reliabranch /
1882
1883
1884	void
1885	annealing(tr)
1886	Tree *tr;
1887	{
1888	int i, k;
1889	double lklorg, lklnew, lkldiff, ldiff, suml1, suml2, sdlkl, nn1, z, rel;
1890	double maxprob;
1891	Node rp, dp, ap, cp, cp0, cp1, kp, kp0, *kp1;
1892	LPVECTOR mlklptrn, lklptrn2, lkltemp;
1893
1894	mlklptrn = NEW_LPVECTOR(Numptrn);
1895	lklptrn2 = NEW_LPVECTOR(Numptrn);
1896	for (k = 0; k < Numptrn; k++) mlklptrn[k] = Alklptrn[k];
1897	lkltemp = Alklptrn; Alklptrn = mlklptrn; mlklptrn = lkltemp;
1898	lklorg = tr->lklhd;
1899	nn1 = (double)Numsite / (Numsite-1);
1900	if (Write_optn) putchar('\n');
1901	for (i = 0; i < Numibrnch; i++) {
1902	maxprob = 1.1;
1903	dp = tr->ibrnchp[i];
1904	ap = dp->kinp;
1905	for (kp = ap->isop, kp0 = ap;
1906	kp->descen \|\| kp == tr->rootp;
1907	kp0 = kp, kp = kp->isop)
1908	;
1909	kp1 = kp->isop;
1910	/* printf(" kp:%3d%3d%3d", kp0->num+1, kp->num+1, kp1->num+1); */
1911	for (cp = dp->isop, cp0 = dp; cp != dp; cp0 = cp, cp = cp->isop) {
1912	cp1 = cp->isop;
1913	kp0->isop = cp; cp->isop = kp1;
1914	cp0->isop = kp; kp->isop = cp1;
1915
1916	/* putctopology(tr); */
1917	/* prtopology(tr); */
1918	#if 1
1919	pathing(tr);
1920	slslength(tr, Distanmat, Maxspc);
1921	#endif
1922	initpartlkl(tr);
1923	rp = (Node *)mlikelihood(tr);
1924	lklnew = tr->lklhd;
1925	lkldiff = lklorg - lklnew;
1926	for (k = 0, suml1 = suml2 = 0.0; k < Numptrn; k++) {
1927	ldiff = Alklptrn[k] - mlklptrn[k];
1928	suml1 += ldiff * Weight[k];
1929	suml2 += ldiff * ldiff * Weight[k];
1930	}
1931	suml1 /= Numsite;
1932	sdlkl = sqrt( nn1 * (suml2 - suml1suml1Numsite) );
1933	z = lkldiff / sdlkl;
1934	rel = probnormal(z);
1935	/* printf("%3d%3d%3d%3d",
1936	i+1+Numspc, cp0->num+1, cp->num+1, cp1->num+1); */
1937	if (Write_optn) {
1938	printf("%3d", i+1+Numspc);
1939	printf("%3d%3d ", kp->num+1, cp->num+1);
1940	printf("%3d%3d", dp->isop->num+1, dp->isop->isop->num+1);
1941	printf(" %7.3f %7.3f %7.3f %.3f\n", i,lkldiff,sdlkl,z,rel);
1942	}
1943	if (rel < maxprob) {
1944	maxprob = rel;
1945	Reliprob[i][0] = rel;
1946	Relinum[i][0] = dp->isop->num;
1947	Relinum[i][1] = dp->isop->isop->num;
1948	}
1949
1950	kp0->isop = kp; kp->isop = kp1;
1951	cp0->isop = cp; cp->isop = cp1;
1952	}
1953	}
1954	lkltemp = Alklptrn; Alklptrn = mlklptrn; mlklptrn = lkltemp;
1955	#if 1
1956	pathing(tr);
1957	slslength(tr, Distanmat, Maxspc);
1958	#endif
1959	initpartlkl(tr);
1960	rp = (Node *)mlikelihood(tr);
1961	mlvalue(tr, Infotrees);
1962
1963	FREE_LPVECTOR(mlklptrn);
1964	FREE_LPVECTOR(lklptrn2);
1965
1966	} /_ annealing /
1967
1968
1969	void
1970	qlrsearch(tr)
1971	Tree *tr;
1972	{
1973	boolean localconv;
1974	int ib, i, j, n, ll, forder;
1975	double lklorg, lklold, lklhd1, lklhd2, uldiff, muldiff, rel, minrel, iblen;
1976	Node rp, dp, ap, cp, cp0, cp1, kp, kp0, kp1, ncp;
1977	Node exchbrnch1, exchbrnch2, ebrn, ibrn;
1978	dvector elenvec, ilenvec, exchldiff, exchrelia;
1979	ivector exchstate, exchorder, whichml;
1980	LPVECTOR mlklptrn;
1981	LPCUBE rlklptrn;
1982
1983	boolean exch_flag;
1984	int mm, rnum;
1985	double rmin;
1986
1987	#define RMAX 100.0
1988
1989	/* root kp1 cp [1] [2]
1990	* R D ---(d) ap I dp (a)--- A X Y
1991	* (a)-----(d)
1992	* Z C ---(a) kp0 cp0 (a)--- B Y X
1993	* kp cp1
1994	*
1995	* [1]: A(X) <--> C(Z) , [2]: B(X) <--> C(Z)
1996	*/
1997	elenvec = new_dvector(Maxspc);
1998	ilenvec = new_dvector(Numibrnch);
1999	exchbrnch1 = (Node *)malloc((unsigned)Numibrnch sizeof(Node *));
2000	if (exchbrnch1 == NULL) maerror("exchbrnch1 in reliabranch().");
2001	exchbrnch2 = (Node *)malloc((unsigned)Numibrnch sizeof(Node *));
2002	if (exchbrnch2 == NULL) maerror("exchbrnch2 in reliabranch().");
2003	exchstate = new_ivector(Numibrnch);
2004	exchorder = new_ivector(Numibrnch);
2005	whichml = new_ivector(Numibrnch);
2006	exchrelia = new_dvector(Numibrnch);
2007	exchldiff = new_dvector(Numibrnch+1);
2008	exchldiff[Numibrnch] = - DBL_MAX;
2009	rlklptrn = NEW_LPCUBE(Numibrnch, 2, Numptrn);
2010	ebrn = tr->ebrnchp;
2011	ibrn = tr->ibrnchp;
2012	mlklptrn = Alklptrn;
2013
2014	for (ib = 0; ib < Numibrnch; ib++) {
2015	Relistat[ib] = -1;
2016	Relinum[ib][0] = Relinum[ib][1] = 0;
2017	Reliprob[ib][0] = Reliprob[ib][1] = 0.0;
2018	}
2019
2020	for (ll = 0, localconv = FALSE; !localconv; ll++) {
2021	if (Verbs_optn) fprintf(stderr, "%2d:", ll+1);
2022	Alklptrn = mlklptrn;
2023	initpartlkl(tr);
2024	rp = (Node *)mlikelihood(tr);
2025	rp = (Node *)relibranch(rp);
2026	mlvalue(tr, Infotrees);
2027	reroot(tr, tr->rootp);
2028	mlklptrn = Alklptrn;
2029	lklold = lklorg;
2030	lklorg = tr->lklhd;
2031	putchar('\n');
2032	prtopology(tr);
2033	#if QRELIDEBUG
2034	for (ib = 0; ib < Numibrnch; ib++) {
2035	printf("%3d%9.3f%9.3f\n",
2036	ib+Maxspc+1, Relitrif[ib], probnormal(Relitrif[ib]));
2037	}
2038	#endif
2039	fflush(stdout);
2040
2041	for (i = 0; i < Numspc; i++) elenvec[i] = ebrn[i]->length;
2042	for (i = 0; i < Numibrnch; i++) ilenvec[i] = ibrn[i]->length;
2043
2044	for (exch_flag = FALSE, mm = 0; !exch_flag; mm++) {
2045	if (Verbs_optn) fprintf(stderr, " %d", mm+1);
2046	for (ib = 0, rnum = -1, rmin = RMAX; ib < Numibrnch; ib++) {
2047	if (Relitrif[ib] < rmin) {
2048	rnum = ib;
2049	rmin = Relitrif[ib];
2050	}
2051	}
2052	if (rnum == -1) {
2053	/* printf("%3d %10.3f min\n", rnum+Maxspc+1, rmin); */
2054	localconv = TRUE;
2055	break;
2056	}
2057
2058	ib = rnum;
2059	dp = ibrn[ib];
2060	kp0 = ap = dp->kinp;
2061	for (kp = ap->isop; kp->descen \|\| kp == tr->rootp;
2062	kp0 = kp, kp = kp->isop)
2063	;
2064	kp1 = kp->isop;
2065
2066	for (cp=dp->isop, cp0=dp, n=0; cp != dp; cp0=cp, cp=cp->isop, n++) {
2067	cp1 = cp->isop;
2068	kp0->isop = cp; cp->isop = kp1;
2069	cp0->isop = kp; kp->isop = cp1;
2070	for (i = 0; i < Numspc; i++)
2071	ebrn[i]->length = ebrn[i]->kinp->length = elenvec[i];
2072	for (i = 0; i < Numibrnch; i++)
2073	ibrn[i]->length = ibrn[i]->kinp->length = ilenvec[i];
2074	iblen = ap->length * 0.5;
2075	if (iblen < Llimit) iblen = Llimit;
2076	dp->length = ap->length = iblen;
2077	n == 0 ? (Alklptrn=rlklptrn[ib][0]):(Alklptrn=rlklptrn[ib][1]);
2078	reliml(tr, ap, lklorg, mlklptrn, &rel);
2079	n == 0 ? (lklhd1 = tr->lklhd) : (lklhd2 = tr->lklhd);
2080	if (n == 0) {
2081	if (lklhd1 > lklorg) {
2082	for (i = 0; i < Numspc; i++) elenvec[i]=ebrn[i]->length;
2083	for (i = 0; i < Numibrnch; i++) ilenvec[i]=ibrn[i]->length;
2084	}
2085	Relinum[ib][0] = dp->isop->num;
2086	Relinum[ib][1] = dp->isop->isop->num;
2087	ncp = cp;
2088	minrel = rel;
2089	} else if (lklhd2 > lklhd1) {
2090	if (lklhd2 > lklorg) {
2091	for (i = 0; i < Numspc; i++) elenvec[i]=ebrn[i]->length;
2092	for (i = 0; i < Numibrnch; i++) ilenvec[i]=ibrn[i]->length;
2093	}
2094	Relinum[ib][0] = dp->isop->num;
2095	Relinum[ib][1] = dp->isop->isop->num;
2096	ncp = cp;
2097	minrel = rel;
2098	}
2099	kp0->isop = kp; kp->isop = kp1;
2100	cp0->isop = cp; cp->isop = cp1;
2101	} /* for cp */
2102
2103	if (lklhd1 > lklhd2) {
2104	lklorg > lklhd1 ? (Relistat[ib] = 0) : (Relistat[ib] = 1);
2105	uldiff = lklhd1 - lklorg;
2106	whichml[ib] = 1;
2107	} else {
2108	lklorg > lklhd2 ? (Relistat[ib] = 0) : (Relistat[ib] = 2);
2109	uldiff = lklhd2 - lklorg;
2110	whichml[ib] = 2;
2111	}
2112	Reliprob[ib][0] = minrel;
2113	Reliprob[ib][1] = 1.0 - minrel;
2114	uldiff > 0 ? (exchstate[ib] = 1) : (exchstate[ib] = 0);
2115	exchbrnch1[ib] = kp;
2116	exchbrnch2[ib] = ncp;
2117	#if QRELIDEBUG
2118	printf("%3d %3d %8.3f", mm+1, ib+Maxspc+1, rmin);
2119	printf(" %3d%3d%3d%6.3f%6.3f\n", Relistat[ib], Relinum[ib][0]+1,
2120	Relinum[ib][1]+1, Reliprob[ib][0], Reliprob[ib][1]);
2121	#endif
2122	if (exchstate[ib]) {
2123	kp = exchbrnch1[ib];
2124	cp = exchbrnch2[ib];
2125	printf("\n%%%-3d %3d %3d<->%-3d ln L:%12.3f +%7.3f\n", ll+1,
2126	ib+Numspc+1, kp->kinp->num+1,cp->kinp->num+1,lklorg,uldiff);
2127	dp = ibrn[ib];
2128	ap = dp->kinp;
2129	kp->isop->isop->isop = cp;
2130	cp->isop->isop->isop = kp;
2131	ncp = cp->isop;
2132	cp->isop = kp->isop;
2133	kp->isop = ncp;
2134	Relistat[ib] = 0;
2135	Relinum[ib][0] = Relinum[ib][1] = 0;
2136	Reliprob[ib][0] = Reliprob[ib][1];
2137	Reliprob[ib][1] = 0.0;
2138	for (i = 0; i < Numspc; i++)
2139	ebrn[i]->length = ebrn[i]->kinp->length = elenvec[i];
2140	for (i = 0; i < Numibrnch; i++)
2141	ibrn[i]->length = ibrn[i]->kinp->length = ilenvec[i];
2142	exch_flag = TRUE;
2143	} else {
2144	Relitrif[ib] = RMAX;
2145	}
2146
2147	} /* for ib */
2148	if (Verbs_optn) fprintf(stderr, "\n");
2149
2150	if (Xreli_optn) goto ONCE;
2151
2152	} /* for ll (!localconv) */
2153
2154	ONCE: /* if (Xreli_optn) */
2155
2156	for (i = 0; i < Numspc; i++)
2157	ebrn[i]->length = ebrn[i]->kinp->length = elenvec[i];
2158	for (i = 0; i < Numibrnch; i++)
2159	ibrn[i]->length = ibrn[i]->kinp->length = ilenvec[i];
2160	Alklptrn = mlklptrn;
2161	initpartlkl(tr);
2162	rp = (Node *)mlikelihood(tr);
2163	mlvalue(tr, Infotrees);
2164	reroot(tr, tr->rootp);
2165	putchar('\n');
2166	putctopology(tr);
2167	localbp(Reliprob, mlklptrn, rlklptrn, whichml, Numibrnch, Numsite);
2168	FREE_LPCUBE(rlklptrn);
2169	free_dvector(exchrelia);
2170	free_dvector(exchldiff);
2171	free_ivector(exchstate);
2172	free_ivector(exchorder);
2173	free_ivector(whichml);
2174	free(exchbrnch1);
2175	free(exchbrnch2);
2176	free_dvector(elenvec);
2177	free_dvector(ilenvec);
2178	} /_ qlrsearch /

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

Download in other formats: