Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

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

Line
1	/*
2	* protml.c Adachi, J. 1995.12.01
3	* Copyright (C) 1992-1995 J. Adachi & M. Hasegawa, All rights reserved.
4	*/
5
6	#define NJMLD 0
7	#define BRANCHULIMIT 0
8	#define PTNLKL 0
9	#define TPGRAPH 0
10	#define DISLKL 0
11	#define TREECHECK 0
12	#define FMLEN 0
13	#define COMPCRITERION 0
14	#define TPMZEROCHECK 0
15
16	#define MAIN_MODULE 1
17	#include "protml.h"
18
19	void
20	copyright()
21	{
22	#ifndef NUC
23	fprintf(stderr, "ProtML %s(%s) ", VERSION, DATE );
24	fprintf(stderr, "Maximum Likelihood Inference of Protein Phylogeny\n");
25	#else /* NUC */
26	fprintf(stderr, "NucML %s(%s) ", VERSION, DATE );
27	fprintf(stderr, "Maximum Likelihood Inference of Nucleic Acid Phylogeny\n");
28	#endif /* NUC */
29	fprintf(stderr, "Copyright (C) 1992-1996 J. Adachi & M. Hasegawa. ");
30	fprintf(stderr, "All rights reserved.\n");
31	/* fprintf(stderr, " ProtML comes with NO WARRANTY\n"); */
32	}
33
34
35	void
36	usage()
37	{
38	copyright();
39	fprintf(stderr, "Usage: %s [switches] sequence_file [topology_file]\n",
40	Prog_name);
41	fprintf(stderr, " Help: %s -h\n", Prog_name);
42	}
43
44
45	void
46	helpinfo()
47	{
48	copyright();
49	fprintf(stderr,
50	"Usage: %s [switches] sequence_file [topology_file]\n",Prog_name);
51	fprintf(stderr,
52	"sequence_file = MOLPHY_format \| Sequential(-S) \| Interleaved(-I)\n");
53	fprintf(stderr,
54	"topology_file = users_trees(-u) \| constrained_tree(-e)\n");
55	fprintf(stderr, "Model:\n");
56	#ifndef NUC
57	fprintf(stderr, "-j JTT (default) -jf JTT-F Jones, Taylor & Thornton(1992)\n");
58	fprintf(stderr, "-d Dayhoff -df Dayhoff-F Dayhoff et al.(1978)\n");
59	fprintf(stderr, "-m %-7s -mf %-7s-F Adachi & Hasegawa(1995)\n", MTREVNAME, MTREVNAME);
60	#else /* NUC */
61	fprintf(stderr,
62	"-t n1 n1: Alpha/Beta ratio (default:%.1f) Hasegawa, Kishino & Yano(1985)\n", ALPHABETA);
63	fprintf(stderr,
64	"-t n1,n2 n2: AlphaY/AlphaR ratio (default:%.1f) Tamura & Nei(1993)\n", ALPHAYR);
65	#endif /* NUC */
66	#ifndef NUC
67	fprintf(stderr, "-p Poisson -pf Proportional (-f: with data Frequencies)\n");
68	fprintf(stderr, "-r users RSR -rf users RSR-F (Relative Substitution Rate)\n");
69	#else /* NUC */
70	fprintf(stderr, "-p Proportional -pf Poisson (-f withOUT data Frequencies)\n");
71	fprintf(stderr, "-r users RSR-F -rf users RSR (Relative Substitution Rate)\n");
72	#endif /* NUC */
73	fprintf(stderr, "Search strategy or Mode:\n");
74	fprintf(stderr, "-u Users trees (need users_trees file)\n");
75	fprintf(stderr, "-R local Rearrangement search -RX LBP only\n");
76	fprintf(stderr, "-e Exhaustive search (with/without constrained_tree file)\n");
77	fprintf(stderr, "-s Star decomposition search (may not be the ML tree)\n");
78	fprintf(stderr, "-q Quick add OTUs search (may not be the ML tree)\n");
79	fprintf(stderr, "-D maximum likelihood Distance matrix --> NJDIST\n");
80	fprintf(stderr, "Others:\n");
81	fprintf(stderr, "-n num retained top ranking trees (default -e:%d,-q:%d)\n", MAXALTREES, NUMQLTREES);
82	fprintf(stderr, "-P num Per cent (default -e:%.0f, if possible trees<=945 : %.0f)\n", TBLRATE100, TBLRATE7100);
83	fprintf(stderr, "-b no RELL-BP -M Minimum evolution (with -e)\n");
84	fprintf(stderr, "-S Sequential format -I Interleaved format\n");
85	fprintf(stderr, "-v verbose to stderr -i, -w output some information\n");
86	}
87
88	static void prologue();
89	void pml();
90
91
92	main(argc, argv)
93	int argc;
94	char **argv;
95	{
96	FILE seqfp, tplfp;
97	int i, ch, num;
98	double x;
99	char *cpp, cp;
100	char buf[64];
101	boolean num_flag;
102	extern int Optindex; /* index of next argument */
103	extern char Optargp; / pointer to argument of current option */
104
105	Ct0 = time(NULL);
106	if((Prog_name = strrchr(argv[0], DIR_CHAR)) != NULL )
107	Prog_name++;
108	else
109	Prog_name = argv[0];
110
111	Rrsr_optn = FALSE;
112	Boots_optn = TRUE; /* default with User_optn*/
113	Exhau_optn = FALSE;
114	Const_optn = FALSE;
115	User_optn = FALSE;
116	Aneal_optn = FALSE;
117	Njoin_optn = FALSE;
118	num_flag = FALSE;
119	Percnt_optn = FALSE;
120	Maxaltree = MAXALTREES;
121	Numqltree = NUMQLTREES;
122	#ifndef NUC
123	Jtt_optn = TRUE;
124	Frequ_optn = FALSE;
125	#else /* NUC */
126	Jtt_optn = FALSE;
127	Frequ_optn = TRUE;
128	Tstv_optn = FALSE;
129	AlphaBeta = ALPHABETA;
130	AlphaYR = ALPHAYR;
131	Beta12 = BETA12;
132	Topting = FALSE;
133	#endif /* NUC */
134	Dayhf_optn = FALSE;
135	Mtrev_optn = FALSE;
136	Linesites = LINESITES;
137	while((ch = mygetopt(argc, argv, SWITCHES)) != -1 ) {
138	switch(ch) {
139	#ifndef NUC
140	case 'j':
141	Jtt_optn = TRUE;
142	Dayhf_optn = FALSE;
143	Mtrev_optn = FALSE;
144	Poisn_optn = FALSE;
145	Rrsr_optn = FALSE;
146	break;
147	case 'd':
148	Dayhf_optn = TRUE;
149	Jtt_optn = FALSE;
150	Mtrev_optn = FALSE;
151	Poisn_optn = FALSE;
152	Rrsr_optn = FALSE;
153	break;
154	case 'm':
155	Mtrev_optn = TRUE;
156	Dayhf_optn = FALSE;
157	Jtt_optn = FALSE;
158	Poisn_optn = FALSE;
159	Rrsr_optn = FALSE;
160	break;
161	case 'f':
162	Frequ_optn = TRUE;
163	break;
164	case 'p':
165	Poisn_optn = TRUE;
166	Jtt_optn = FALSE;
167	Dayhf_optn = FALSE;
168	Mtrev_optn = FALSE;
169	Rrsr_optn = FALSE;
170	break;
171	case 'r':
172	Rrsr_optn = TRUE;
173	Jtt_optn = FALSE;
174	Dayhf_optn = FALSE;
175	Mtrev_optn = FALSE;
176	Poisn_optn = FALSE;
177	break;
178	#else /* NUC */
179	case 't':
180	cpp = &Optargp;
181	Tstv_optn = TRUE;
182	if (x = strtod(Optargp, cpp)) {
183	AlphaBeta = x;
184	if (**cpp == ',') {
185	Optargp = ++(*cpp);
186	if (x = strtod(Optargp, cpp)) {
187	AlphaYR = x;
188	if (**cpp == ',') {
189	Optargp = ++(*cpp);
190	if (x = strtod(Optargp, cpp)) {
191	Beta12 = x;
192	}
193	}
194	}
195	}
196	} else {
197	Toptim_optn = TRUE;
198	if (cp = strchr(Optargp, ',')) {
199	Optargp = ++cp;
200	if (x = strtod(Optargp, cpp)) {
201	AlphaYR = x;
202	} else {
203	Toptim_optn += 2;
204	AlphaYR = 1.0;
205	}
206	if (strchr(cp, ',')) Toptim_optn += 4;
207	Beta12 = 1.0;
208	}
209	}
210	break;
211	case 'f':
212	Frequ_optn = FALSE;
213	break;
214	case 'p':
215	Poisn_optn = TRUE;
216	break;
217	case 'r':
218	Rrsr_optn = TRUE;
219	Poisn_optn = FALSE;
220	break;
221	#endif /* NUC */
222	case 'F':
223	Logdet_optn = TRUE;
224	break;
225
226
227	case 'a': Aneal_optn = TRUE;
228	Njoin_optn = FALSE;
229	Exhau_optn = FALSE;
230	User_optn = FALSE;
231	Stard_optn = FALSE;
232	Quick_optn = FALSE;
233	Boots_optn = FALSE; break;
234	case 'e': Exhau_optn = TRUE;
235	Njoin_optn = FALSE;
236	Aneal_optn = FALSE;
237	User_optn = FALSE;
238	Stard_optn = FALSE;
239	Quick_optn = FALSE;
240	Boots_optn = FALSE; break;
241	case 'u': User_optn = TRUE;
242	Njoin_optn = FALSE;
243	Aneal_optn = FALSE;
244	Exhau_optn = FALSE;
245	Stard_optn = FALSE;
246	Quick_optn = FALSE; break;
247	case 'q': Quick_optn = TRUE;
248	Njoin_optn = FALSE;
249	Aneal_optn = FALSE;
250	User_optn = FALSE;
251	Stard_optn = FALSE;
252	Exhau_optn = FALSE; break;
253	case 'Q': Quick1_optn = TRUE;
254	Njoin_optn = FALSE;
255	Aneal_optn = FALSE;
256	User_optn = FALSE;
257	Stard_optn = FALSE;
258	Exhau_optn = FALSE; break;
259	case 's': Stard_optn = TRUE;
260	Njoin_optn = FALSE;
261	Aneal_optn = FALSE;
262	Quick_optn = FALSE;
263	Quick1_optn = FALSE;
264	User_optn = FALSE;
265	Exhau_optn = FALSE; break;
266	case 'N': Njoin_optn = TRUE;
267	Stard_optn = FALSE;
268	Aneal_optn = FALSE;
269	Quick_optn = FALSE;
270	Quick1_optn = FALSE;
271	User_optn = FALSE;
272	Exhau_optn = FALSE; break;
273
274	case 'n': cpp = &Optargp;
275	num_flag = TRUE;
276	if (i = strtol(Optargp, cpp, 10)) num = i; break;
277	case 'o':
278	Outgr_optn = TRUE;
279	cpp = &Optargp;
280	if (i = strtol(Optargp, cpp, 10)) {
281	Outgroup1 = i - 1;
282	printf("Outgroup1 %3d\n", Outgroup1+1);
283	if (**cpp == ',') {
284	Optargp = ++(*cpp);
285	if (i = strtol(Optargp, cpp, 10)) {
286	Outgr_optn = 2;
287	Outgroup2 = i - 1;
288	printf("Outgroup2 %3d\n", Outgroup2+1);
289	} else {
290	fprintf(stderr,"%s: bad out group option \"%s\"\n",
291	Prog_name, Optargp);
292	exit(1);
293	}
294	}
295	} else {
296	fprintf(stderr,"%s: bad out group option \"%s\"\n",
297	Prog_name, Optargp);
298	exit(1);
299	}
300	break;
301	case 'P':
302	cpp = &Optargp;
303	Percnt_optn = TRUE;
304	if (x = strtod(Optargp, cpp)) {
305	Percent = x;
306	/* printf("%8.3f\n", x); */
307	if (x < 0.0) {
308	fprintf(stderr,"%s: bad par cent(-P) open \"%.1f\"\n",
309	Prog_name, x);
310	exit(1);
311	}
312	} else {
313	fprintf(stderr,"%s: bad par cent(-P) open \"%s\"\n",
314	Prog_name, Optargp);
315	exit(1);
316	}
317	break;
318	case 'b': Boots_optn = FALSE; break;
319	case 'R': Relia_optn = TRUE; break;
320	case 'T': Triad_optn = TRUE; break;
321	case 'D': Distn_optn = TRUE; break;
322	case 'A': Aprox_optn = TRUE; break;
323	case 'l':
324	Lklhd_optn = TRUE;
325	Llsfile = new_cvector(strlen(Optargp) + strlen(LLSFILEEXT) + 1);
326	*Llsfile = '\0';
327	if (strcat(Llsfile, Optargp) && strcat(Llsfile, LLSFILEEXT)) {
328	/* printf("Llsfile is \"%s\"\n", Llsfile); */
329	break;
330	}
331	case 'L': Logfl_optn = TRUE; break;
332	case 'i': Info_optn = TRUE; break;
333	case 'M': Mevol_optn = TRUE; break;
334	case 'v': Verbs_optn = TRUE; break;
335	case 'V': Varia_optn = TRUE; break;
336	case 'w': Write_optn = TRUE; break;
337	case 'X': Xreli_optn = TRUE; break;
338	case 'S': Seque_optn = TRUE;
339	Inlvd_optn = FALSE; break;
340	case 'I': Inlvd_optn = TRUE;
341	Seque_optn = FALSE; break;
342	#if 1
343	case 'z': Debug_optn = TRUE; break;
344	case 'Z': Debug = TRUE; break;
345	#endif
346	case 'C': Ctacit_optn = TRUE; break;
347	case 'h':
348	case 'H': helpinfo(); exit(1);
349	default : usage(); exit(1);
350	}
351	}
352	if (Outgr_optn == 1) {
353	printf("Outgroup1 %3d\n", Outgroup1+1);
354	} else if (Outgr_optn == 2) {
355	printf("Outgroup1 %3d Outgroup2 %3d\n", Outgroup1+1, Outgroup2+1);
356	} else {
357	}
358	if (Exhau_optn) {
359	if (Optindex + 2 == argc) Const_optn = TRUE;
360	}
361	#if TPGRAPH
362	#else /* TPGRAPH */
363	if (!Aneal_optn && !Exhau_optn && !User_optn && !Quick_optn &&
364	!Quick1_optn && !Stard_optn && !Distn_optn && !Njoin_optn ) {
365	if (Optindex == argc) {
366	helpinfo(); exit(1); /* default !? */
367	} else if (Optindex + 1 == argc) {
368	helpinfo(); exit(1); /* default !? */
369	} else if (Optindex + 2 == argc) {
370	User_optn = TRUE; /* default !? */
371	}
372	}
373	#endif /* TPGRAPH */
374	if (num_flag) {
375	if (Quick_optn \|\| Quick1_optn) Numqltree = num;
376	if (Exhau_optn) Maxaltree = num;
377	}
378	*Modelname = '\0';
379	#ifndef NUC
380	if (Poisn_optn) {
381	if (!Frequ_optn) strcat(Modelname, "Poisson");
382	else strcat(Modelname, "Proportional");
383	} else {
384	if (Rrsr_optn) {
385	strcat(Modelname, "Read RSR");
386	} else if (Jtt_optn) {
387	strcat(Modelname, "JTT");
388	} else if (Dayhf_optn) {
389	strcat(Modelname, "Dayhoff");
390	} else if (Mtrev_optn) {
391	strcat(Modelname, MTREVNAME);
392	} else {
393	strcat(Modelname, "unknown");
394	}
395	if (Frequ_optn) strcat(Modelname, "-F");
396	}
397	#else /* NUC */
398	if (Poisn_optn) {
399	if (!Frequ_optn) strcat(Modelname, "Poisson");
400	else strcat(Modelname, "Proportional");
401	AlphaBeta = 1.0;
402	AlphaYR = 1.0;
403	} else if (Rrsr_optn) {
404	strcat(Modelname, "Read RSR");
405	if (Frequ_optn) strcat(Modelname, "-F");
406	} else {
407	if (Toptim_optn) {
408	strcat(Modelname, "A/B:opt");
409	if (Toptim_optn == 3 \|\| Toptim_optn == 7)
410	strcat(Modelname, " Ay/Ar:opt");
411	if (Toptim_optn >= 4)
412	strcat(Modelname, " B1/B2:opt");
413	} else {
414	sprintf(buf, "A/B:%.2f\0", AlphaBeta);
415	strcat(Modelname, buf);
416	if (AlphaYR != 1.0) {
417	sprintf(buf, " Ay/Ar:%.2f", AlphaYR);
418	strcat(Modelname, buf);
419	}
420	if (Beta12 != 1.0) {
421	sprintf(buf, " B1/B2:%.2f", Beta12);
422	strcat(Modelname, buf);
423	}
424	}
425	if (Frequ_optn) strcat(Modelname, " F");
426	}
427	#endif /* NUC */
428	/* printf("\"%s\"\n", Modelname); */
429
430	#ifdef DEBUG
431	if (Debug) {
432	printf("argc = %d\n",argc);
433	for(i = 0; i < argc; i++) printf("argv[%d] = %s\n",i,argv[i]);
434	putchar('\n');
435	printf("\nOptindex = %d\n",Optindex);
436	printf("Optargp = %s\n",Optargp);
437	}
438	#endif
439
440	Aneal_mode = 0;
441	if (Optindex == argc) {
442	seqfp = stdin;
443	tplfp = stdin;
444	Aneal_mode = 0;
445	} else if (Optindex + 1 == argc) {
446	if ((seqfp = fopen(argv[Optindex++], "r")) == NULL) {
447	fprintf(stderr,"%s: can't open %s\n",Prog_name,argv[--Optindex]);
448	exit(1);
449	} else {
450	tplfp = seqfp;
451	Aneal_mode = 0;
452	}
453	} else if (Optindex + 2 == argc) {
454	if ((seqfp = fopen(argv[Optindex++], "r")) == NULL) {
455	fprintf(stderr,"%s: can't open %s\n",Prog_name,argv[--Optindex]);
456	exit(1);
457	} else if ((tplfp = fopen(argv[Optindex], "r")) == NULL) {
458	if (*argv[Optindex] == '-') {
459	tplfp = stdin;
460	Aneal_mode = 1;
461	} else {
462	fprintf(stderr,"%s: can't open %s\n",Prog_name,argv[Optindex]);
463	exit(1);
464	}
465	}
466	} else {
467	fprintf(stderr, "%s: Inconsistent number of operand in command line!\n",
468	Prog_name);
469	usage();
470	exit(1);
471	}
472
473	if (Rrsr_optn) {
474	if ((Rtfifp = fopen(RSRINFILE, "r")) == NULL) {
475	fprintf(stderr,
476	"%s: can't open RSR file: %s\n",Prog_name,RSRINFILE);
477	exit(1);
478	}
479	}
480	if (Logfl_optn) {
481	if ((Logfp = fopen(LOGFILE, "w")) == NULL) {
482	fprintf(stderr,
483	"%s: can't open log file: %s\n",Prog_name,LOGFILE);
484	exit(1);
485	}
486	}
487	if (Lklhd_optn) {
488	if ((Lklfp = fopen(Llsfile, "w")) == NULL) {
489	fprintf(stderr,
490	"%s: can't open log likelihood file: %s\n",Prog_name,Llsfile);
491	exit(1);
492	}
493	}
494
495	Epsfile = new_cvector(strlen(EPSFILE) + 1);
496	*Epsfile = '\0';
497	strcat(Epsfile, EPSFILE);
498	if ((Epsfp = fopen(Epsfile, "w")) == NULL) {
499	fprintf(stderr,"%s: can't open eps tree file: %s\n",Prog_name,Epsfile);
500	exit(1);
501	}
502	if (Relia_optn \|\| Njoin_optn) {
503	if ((Tplfp = fopen(TPLFILE, "w")) == NULL) {
504	fprintf(stderr,"%s: can't open tree topology file: %s\n",
505	Prog_name,TPLFILE);
506	exit(1);
507	}
508	}
509	if (User_optn \|\| Stard_optn \|\| Njoin_optn) {
510	if ((Trefp = fopen(TREFILE, "w")) == NULL) {
511	fprintf(stderr,"%s: can't open tree file: %s\n",Prog_name,TREFILE);
512	exit(1);
513	}
514	}
515
516	if (Verbs_optn) { fputc('\n', stderr); copyright(); }
517	Numexe = 1;
518	for (Cnoexe = 0; Cnoexe < Numexe; Cnoexe++) {
519	prologue(seqfp, stdout);
520	if (!Distn_optn) pml(tplfp, stdout);
521	}
522	#if __STDC__ && DIFFTIME
523	if (Verbs_optn) {
524	Ct1 = time(NULL);
525	x = difftime(Ct1, Ct0);
526	fprintf(stderr, "elapse time %02d:%02d:%02d\n",
527	(int)(x/3600), (int)(x/60), (int)x % 60);
528	}
529	#endif
530	if (tplfp != stdin && tplfp != seqfp) fclose(tplfp);
531	if (seqfp != stdin) fclose(seqfp);
532	fclose(Epsfp);
533	if (Relia_optn \|\| Njoin_optn) fclose(Tplfp);
534	if (User_optn \|\| Stard_optn \|\| Njoin_optn) fclose(Trefp);
535	if (Lklhd_optn) fclose(Lklfp);
536	if (Logfl_optn) fclose(Logfp);
537	if (Rrsr_optn) fclose(Rtfifp);
538
539	return 0;
540	}
541
542	void
543	header(ofp, maxspc, numsite, commentp)
544	FILE *ofp;
545	int *maxspc;
546	int *numsite;
547	char **commentp;
548	{
549	char date[32];
550
551	strftime(date, 32, "%x", localtime(&Ct0));
552	fprintf(ofp, "%s %s (%s) %s", Prog_name, VERSION, date, Modelname);
553	fprintf(ofp, " %d OTUs %d sites. %s\n", maxspc, numsite, *commentp);
554	} /_ header /
555
556
557	void
558	headerd(ofp, maxspc, numsite, commentp)
559	FILE *ofp;
560	int *maxspc;
561	int *numsite;
562	char **commentp;
563	{
564
565	fprintf(ofp, "%d %d sites %s %s\n",
566	maxspc, numsite, Modelname, *commentp);
567	} /_ headerd /
568
569
570	#ifdef NUC
571	#include "optimtpm.c"
572	#include "abratio.c"
573	#endif /* NUC */
574
575	#if COMPCRITERION
576	#include "compcrit.c"
577	#endif /* COMPCRITERION */
578
579
580	static void
581	prologue(ifp, ofp)
582	FILE *ifp;
583	FILE *ofp;
584	{
585	int i, j, k;
586	ivector alias;
587	double maxdis;
588	#if TPMZEROCHECK
589	int mini, minj;
590	double dis, minprob;
591	dmattpmty tpm;
592	#endif
593	#if 0
594	dmattpmty tpm;
595	dvectpmty x;
596	double y;
597	#endif
598
599	getsize(ifp, &Maxspc, &Maxsite, &Comment);
600	Numsites = new_ivector(Maxspc);
601	if (Distn_optn) {
602	if (Maxspc < 2) exit(1);
603	} else {
604	if (Maxspc < 3) exit(1);
605	}
606	if ((User_optn \|\| Aneal_optn \|\| Stard_optn \|\| Njoin_optn) && !Ctacit_optn)
607	header(ofp, &Maxspc, &Maxsite, &Comment);
608	if (Distn_optn && Maxspc > 2) headerd(ofp, &Maxspc, &Maxsite, &Comment);
609	Identif = (char *)malloc((unsigned)Maxspc sizeof(char *));
610	if (Identif == NULL) maerror("in prologue, Identif");
611	Sciname = (char *)malloc((unsigned)Maxspc sizeof(char *));
612	if (Sciname == NULL) maerror("in prologue, Sciname");
613	Engname = (char *)malloc((unsigned)Maxspc sizeof(char *));
614	if (Engname == NULL) maerror("in prologue, Engname");
615	Seqchar = new_cmatrix(Maxspc, Maxsite);
616	if (Seque_optn)
617	getseqs(ifp, Identif, Seqchar, Maxspc, Maxsite);
618	else if (Inlvd_optn)
619	getseqi(ifp, Identif, Seqchar, Maxspc, Maxsite);
620	else
621	getseq(ifp, Identif, Sciname, Engname, Seqchar, Maxspc, Maxsite);
622	if (Debug_optn) prsequence(ofp, Identif, Seqchar, Maxspc, Maxsite);
623	getfreqepm(Seqchar, Freqemp, Maxspc, Maxsite);
624	alias = new_ivector(Maxsite);
625	radixsort(Seqchar, alias, Maxspc, Maxsite, &Numptrn);
626	Seqconint = new_imatrix(Maxspc, Numptrn);
627	Weight = new_ivector(Numptrn);
628	condenceseq(Seqchar, alias, Seqconint, Weight, Maxspc, Maxsite, Numptrn);
629	convseq(Seqconint, Maxspc, Numptrn);
630	getnumsites(Seqconint, Numsites, Weight, Maxspc, Numptrn);
631	if (Debug_optn)
632	prcondenceseq(Identif, Seqconint, Weight, Maxspc, Maxsite, Numptrn);
633	free_ivector(alias);
634	checkseq(Seqconint, Maxspc, Numptrn);
635	if (Frequ_optn) convfreq(Freqemp);
636	if (Toptim_optn) Topting = TRUE;
637	#ifdef NUCxxxxxxx
638	if (!Toptim_optn && (AlphaBeta == ALPHABETA)) abratio();
639	#endif /* NUC */
640	tranprobmat();
641	#if TPMZEROCHECK
642	for (k = 1, dis = 1.0; k < 100; k++) {
643	dis = dis / (double)k;
644	tprobmtrx(dis, tpm);
645	for (i = 1, minprob = 2.0; i < Tpmradix; i++) {
646	for (j = 0; j < i; j++) {
647	if (tpm[i][j] < minprob) {
648	minprob = tpm[i][j];
649	mini = i;
650	minj = j;
651	}
652	if (tpm[i][j] < 0.0)
653	printf("%3d %20.15f %2d %2d %25.20f\n",
654	k, dis, i, j, tpm[i][j]);
655	}
656	}
657	printf("%3d %20.15f %2d %2d %25.20f\n", k, dis, mini, minj, minprob);
658	}
659	exit(1);
660	#endif
661	if (Toptim_optn) Topting = FALSE;
662	if (Write_optn) prfreq();
663
664	Distanmat = new_dmatrix(Maxspc, Maxspc);
665	if (Logdet_optn)
666	lddistance(Distanmat, Seqchar, Maxspc, Maxsite);
667	else
668	distance(Distanmat, Seqchar, Maxspc, Maxsite);
669	#if 1
670	if (Triad_optn) tridistance(Distanmat, Seqconint, Weight, Maxspc, Numptrn);
671	#else
672	if (Triad_optn) tridistance2(Distanmat, Seqchar, Maxspc, Maxsite);
673	#endif
674	if (Distn_optn) putdistance(Identif, Sciname, Engname, Distanmat, Maxspc);
675	Maxbrnch = 2 * Maxspc - 3;
676	Maxibrnch = Maxspc - 3;
677	Maxpair = (Maxspc * (Maxspc - 1)) / 2;
678	#if BRANCHULIMIT
679	Llimit = 100.0 / Maxsite;
680	for (i = 1, maxdis = Llimit; i < Maxspc; i++) {
681	for (j = 0; j < i; j++) {
682	if (Distanmat[i][j] > maxdis) maxdis = Distanmat[i][j];
683	}
684	}
685	Ulimit = maxdis;
686	Mlimit = (Ulimit - Llimit) * 0.5;
687	#else /* BRANCHULIMIT */
688	Llimit = LOWERLIMIT;
689	Ulimit = UPPERLIMIT;
690	Mlimit = MIDLIMIT;
691	#endif /* BRANCHULIMIT */
692	#if 0
693	tprobmtrx(1.0, tpm); /* 1 - 10000 */
694	for (i = 0, y = 0.0; i < Tpmradix; i++) {
695	for (j = 0, x[i] = 0.0; j < Tpmradix; j++) {
696	x[i] += tpm[i][j];
697	if (i != j) y += Freqtpm[i] * tpm[i][j];
698	}
699	}
700	putchar('\n');
701	if (Tpmradix == NUMAMI) {
702	for (i = 0; i < Tpmradix; i++) {
703	for (j = 0; j < 10; j++) {
704	printf("%7d", (int)(tpm[i][j]*1.0e6));
705	} putchar('\n');
706	}
707	for (j = 0; j < 10; j++) printf("%7.4f", x[j]); putchar('\n');
708	putchar('\n');
709	for (i = 0; i < Tpmradix; i++) {
710	for (j = 10; j < Tpmradix; j++) {
711	printf("%7d", (int)(tpm[i][j]*1.0e6));
712	} putchar('\n');
713	}
714	for (j = 10; j < Tpmradix; j++) printf("%7.4f", x[j]); putchar('\n');
715	} else {
716	for (i = 0; i < Tpmradix; i++) {
717	for (j = 0; j < Tpmradix; j++) {
718	printf("%7d", (int)(tpm[i][j]*1.0e6));
719	} putchar('\n');
720	}
721	for (j = 0; j < Tpmradix; j++) printf("%7.4f", x[j]); putchar('\n');
722	}
723	printf("%20.10f\n\n", y);
724	exit(1);
725	#endif
726	}
727
728
729	#ifdef BDATE
730	#include "bdate.c"
731	#endif /* BDATE */
732
733	#if PTNLKL
734	#include "ptnlkl.c"
735	#endif /* PTNLKL */
736
737	#if TPGRAPH
738	#include "tpgraph.c"
739	#endif /* TPGRAPH */
740
741
742	#ifndef TPM
743	void
744	pml(ifp, ofp)
745	FILE *ifp;
746	FILE *ofp;
747	{
748	int buftree, cnospc, i;
749	Node rp, poolnode, *addposition;
750	Infoaddtree topaddtree, curp, *insp;
751	char *cltree;
752	double lkla, nt, mt;
753	#if DISLKL
754	int k, l;
755	double dlkl, dislkl;
756	dmatrix dislklhd, iprb, jprb;
757	#endif /* DISLKL */
758
759	buftree = getbuftree(Maxspc, Identif);
760	if (Debug) fprintf(ofp, "buftree: %d\n", buftree);
761	Distanvec = new_dvector(Maxpair);
762	changedistan(Distanmat, Distanvec, Maxspc);
763	Brnlength = new_dvector(Maxbrnch);
764	if (!Cnoexe) getproportion(&Proportion, Distanvec, Maxpair);
765	Epsilon = EPSILON;
766	Numspc = Maxspc;
767	Numbrnch = Maxbrnch;
768	Numpair = Maxpair;
769	Numsite = Maxsite;
770	Numibrnch = Numspc - 3;
771	Converg = TRUE;
772	Numit = 0;
773
774	#if TPGRAPH
775	tpgraph();
776	exit(1);
777	#endif /* TPGRAPH */
778
779	if (User_optn) { /* Users trees MODE */
780
781	if (!Cnoexe)
782	Ctree = (Tree *) new_tree(Maxspc, Maxibrnch, Numptrn, Seqconint);
783	getnumtree(ifp, &Numtree);
784	Infotrees = (Infotree *) newinfotrees(Numtree, buftree);
785	Lklptrn = NEW_LPMATRIX(Numtree, Numptrn);
786	Strtree = new_cvector(buftree);
787	if (Relia_optn) {
788	Relistat = new_ivector(Numibrnch);
789	Reliprob = new_dmatrix(Numibrnch, 2);
790	Relinum = new_imatrix(Numibrnch, 2);
791	fprintf(Tplfp, "%d\n", Numtree);
792	}
793	if (Numtree != 1) fprintf(Trefp, "%d\n", Numtree);
794	for (Cnotree = 0; Cnotree < Numtree; Cnotree++) {
795	if (!Ctacit_optn && !Aprox_optn) printf("#%d\n", Cnotree + 1);
796	getusertree(ifp, Strtree, buftree);
797	if (Debug) puts(Strtree);
798	constructtree(Ctree, Strtree);
799	#if TREECHECK
800	prtopology(Ctree);
801	exit(1);
802	#endif
803	if (Debug) prcurtree(Ctree); /* */
804	Alklptrn = Lklptrn[Cnotree];
805	#ifdef NUC
806	if (Toptim_optn) {
807	Epsilon = REPSILON;
808	optimtpm();
809	tranprobmat();
810	Epsilon = EPSILON;
811	}
812	#endif /* NUC */
813
814	#if FMLEN
815	fmlength(Ctree, Distanmat, Maxspc);
816	#else
817	slslength(Ctree, Distanmat, Maxspc);
818	/* lslength(Ctree, Distanvec, Maxspc); */
819	#endif
820	if (Debug) prcurtree(Ctree);
821	if (Debug_optn) prcurtree(Ctree);
822	if (Aprox_optn) {
823	initpartlkl(Ctree);
824	aproxlkl(Ctree);
825	praproxlkl(Ctree);
826	putctopology(Ctree);
827	} else {
828	initpartlkl(Ctree);
829	if (Ctacit_optn) aproxlkl(Ctree);
830	rp = (Node *)mlikelihood(Ctree);
831	/* rp = (Node )relibranch(rp); /
832	mlvalue(Ctree, Infotrees);
833	if (Debug_optn) putctopology(Ctree);
834	if (Relia_optn) {
835	Epsilon = REPSILON;
836	reliabranch(Ctree);
837	fputctopology(Tplfp, Ctree);
838	Epsilon = EPSILON;
839	}
840	if (!Ctacit_optn) prtopology(Ctree);
841	strctree(Ctree, Infotrees[Cnotree].ltplgy);
842	if (Debug) puts(Infotrees[Cnotree].ltplgy);
843	if (!Ctacit_optn) resulttree(Ctree);
844	if (Cnotree == 0) pstree(Epsfp, Ctree);
845	fputcphylogeny(Trefp, Ctree);
846	#if VARITPM
847	varitpm(Ctree);
848	#endif /* VARITPM */
849	#if COMPCRITERION
850	if (Ctacit_optn) compcrit(Ctree);
851	#endif /* COMPCRITERION */
852	#ifdef BDATE
853	branchdate(Ctree);
854	#endif /* BDATE */
855	#if PTNLKL
856	xpatternlklhd(Ctree);
857	/* patternlklhd(Ctree); */
858	#endif /* PTNLKL */
859	}
860	if (Verbs_optn) {
861	fprintf(stderr," %d", Cnotree + 1);
862	#if __STDC__ && DIFFTIME
863	Ct1 = time(NULL);
864	fprintf(stderr, "(%.0fs)", difftime(Ct1, Ct0));
865	#endif
866	}
867	}
868	if (Ctacit_optn && Numtree == 1) printf("%.10f\n", Ctree->lklhd);
869	if (Relia_optn) {
870	free_ivector(Relistat);
871	free_dmatrix(Reliprob);
872	free_imatrix(Relinum);
873	}
874	if (Lklhd_optn) outlklhd(Lklptrn);
875	if (!Aprox_optn && Numtree > 1 && (!COMPCRITERION \|\| !Ctacit_optn)) {
876	if (Boots_optn && Verbs_optn) fputs("\nbootstraping\n", stderr);
877	if (Boots_optn) bootstrap(Infotrees, Lklptrn);
878	tabletree(Infotrees, Lklptrn);
879	if (Info_optn) tableinfo(Infotrees);
880	}
881	if (Verbs_optn) fputs("\n", stderr);
882	FREE_LPMATRIX(Lklptrn);
883
884	} else if (Aneal_optn) { /* Annealing MODE */
885
886	Relia_optn = TRUE;
887	Ctree = (Tree *) new_tree(Maxspc, Maxibrnch, Numptrn, Seqconint);
888	if (!Aneal_mode) getnumtree(ifp, &Numtree);
889	if (!Numtree) Numtree = 1; /* !? */
890	Infotrees = (Infotree *) newinfotrees(Numtree, buftree);
891	Lklptrn = NEW_LPMATRIX(Numtree, Numptrn);
892	Strtree = new_cvector(buftree);
893	if (Relia_optn) {
894	Reliprob = new_dmatrix(Numibrnch, 2);
895	Relinum = new_imatrix(Numibrnch, 2);
896	}
897	Cnotree = 0;
898
899	if (!Ctacit_optn && !Aprox_optn) printf("\n");
900	getusertree(ifp, Strtree, buftree);
901	if (Debug) puts(Strtree);
902	constructtree(Ctree, Strtree);
903	if (Debug) prcurtree(Ctree); /* */
904	Alklptrn = Lklptrn[Cnotree];
905	#ifdef NUC
906	if (Toptim_optn) {
907	optimtpm();
908	tranprobmat();
909	}
910	#endif /* NUC */
911	#if FMLEN
912	fmlength(Ctree, Distanmat, Maxspc);
913	#else
914	lslength(Ctree, Distanvec, Maxspc);
915	#endif
916	if (Debug) prcurtree(Ctree);
917	if (Debug_optn) prcurtree(Ctree);
918	initpartlkl(Ctree);
919	rp = (Node *)mlikelihood(Ctree);
920	mlvalue(Ctree, Infotrees);
921	if (Debug_optn) putctopology(Ctree);
922
923	if (Relia_optn) annealing(Ctree);
924
925	if (!Ctacit_optn) prtopology(Ctree);
926	strctree(Ctree, Infotrees[Cnotree].ltplgy);
927	if (Debug) puts(Infotrees[Cnotree].ltplgy);
928	if (!Ctacit_optn) resulttree(Ctree);
929
930	if (Ctacit_optn && Numtree == 1) printf("%.10f\n", Ctree->lklhd);
931	if (Verbs_optn) fputs("\n", stderr);
932	if (Lklhd_optn) outlklhd(Lklptrn);
933	if (Relia_optn) {
934	free_dmatrix(Reliprob);
935	free_imatrix(Relinum);
936	}
937	FREE_LPMATRIX(Lklptrn);
938
939	} else if (Exhau_optn) { /* Exhaustive search */
940
941	Numtree = 1;
942	Cnotree = 0;
943	Numspc = Maxspc;
944	Numbrnch = Maxspc;
945	Maxibrnch = Maxspc - 3;
946	/* Infotrees = (Infotree ) newinfotrees(1, buftree); /
947	/* Alklptrn = NEW_LPVECTOR(Numptrn); */
948	if (Maxspc < 4) {
949	fprintf(stderr,"%s: number of OTUs is \"%d\".\n",Prog_name,Maxspc);
950	exit(1);
951	}
952	if (!Const_optn) {
953	for (i = Maxspc * 2 - 5, nt = 1.0; i > 0; i--) nt *= i;
954	for (i = Maxspc - 3, mt = 1.0; i > 0; i--) mt *= i;
955	for (i = Maxspc - 3; i > 0; i--) mt *= 2;
956	nt /= mt;
957	Numtplgy = nt;
958	} else {
959	Numtplgy = 1.0;
960	}
961
962	Ctree = (Tree *) new_njtree(Maxspc, Maxibrnch, Numptrn, Seqconint);
963	enjtree(Ctree, Distanmat, Maxspc, TRUE);
964	pathing(Ctree);
965	Numibrnch = Maxibrnch;
966	slslength(Ctree, Distanmat, Maxspc);
967	Mintbldm = Ctree->tbldis;
968	/* printf("TBL: %7.3f ", Mintbldm); putctopology(Ctree); */
969	free_njtree(Ctree, Maxspc, Maxibrnch);
970
971	Ctree = (Tree *) new_njtree(Maxspc, Maxibrnch, Numptrn, Seqconint);
972	enjtree(Ctree, Distanmat, Maxspc, FALSE);
973	pathing(Ctree);
974	Numibrnch = Maxibrnch;
975	slslength(Ctree, Distanmat, Maxspc);
976	Maxtbldm = Ctree->tbldis;
977	/* printf("TBL: %7.3f ", Maxtbldm); putctopology(Ctree); */
978	free_njtree(Ctree, Maxspc, Maxibrnch);
979
980	Poolorder = new_ivector(Maxspc);
981	Ctree = (Tree *) new_atree(Maxspc, Maxibrnch, Numptrn, Seqconint);
982	poolnode = (Node *)malloc((unsigned)Maxspc sizeof(Node *));
983	if (poolnode == NULL) maerror("poolnode.");
984	addposition = (Node *)malloc((unsigned)Maxspc sizeof(Node *));
985	if (addposition == NULL) maerror("addposition.");
986	if (!Const_optn) { /* without constrained tree */
987	atreeinit(Ctree, poolnode, addposition, Poolorder);
988	} else { /* with constrained tree */
989	Strtree = new_cvector(buftree);
990	getusertree(ifp, Strtree, buftree);
991	if (Debug) puts(Strtree);
992	streeinit(Ctree, Strtree, poolnode, addposition, Poolorder);
993	}
994	Infoaltrees = (Infoaltree *) newinfoaltrees(Maxaltree, buftree);
995
996	if (Percnt_optn) {
997	if (Percent > 100.0)
998	Tblrate = 2.0;
999	else
1000	Tblrate = Percent / 100.0;
1001	} else {
1002	if (Numtplgy <= 945)
1003	Tblrate = TBLRATE7;
1004	else
1005	Tblrate = TBLRATE;
1006	}
1007	Basetbldm = (Maxtbldm - Mintbldm) * Tblrate + Mintbldm;
1008	/* printf("TBL: Min%7.3f max%7.3f bese%7.3f\n",
1009	Mintbldm, Maxtbldm, Basetbldm); */
1010	if (Verbs_optn) {
1011	fprintf(stderr," %.0f possible trees, TBL limit %.1f%%\n",
1012	Numtplgy, Tblrate*100);
1013	if (Maxspc > 9) Numverbs = 10000; else Numverbs = 1000;
1014	}
1015	if (Numtplgy > 655e6) {
1016	fprintf(stderr,"%s: too many possible trees \"%.0f\".\n",
1017	Prog_name, Numtplgy);
1018	exit(1);
1019	}
1020
1021	Ahead.lklaprox = 0.0;
1022	Atail.up = &Ahead;
1023	Tblbin = new_ivector(NUMTBLBIN);
1024	for (i = 0; i < NUMTBLBIN; i++) Tblbin[i] = 0;
1025	Tblunder = Tblover = 0;
1026	Tblcoef = NUMTBLBIN / (Maxtbldm - Mintbldm + 0.00001);
1027	if (addposition[3] == NULL)
1028	autoconstruction(Ctree, 3, poolnode, addposition, Poolorder);
1029	else
1030	wedge(Ctree, 3, poolnode, addposition, Poolorder, addposition[3]);
1031	/* if (Cnotree < Maxaltree) Numaltree = Cnotree; */
1032	if (Verbs_optn) fputc('\n', stderr);
1033	if (!Aprox_optn) tablealtree(Numaltree);
1034	free_ivector(Poolorder);
1035	free_ivector(Tblbin);
1036
1037	} else if (Stard_optn) { /* Star Decomposition MODE */
1038
1039	Numspc = Maxspc;
1040	Numibrnch = 0;
1041	Numbrnch = Maxspc;
1042	Cnotree = 0;
1043	Numtree = MAXSLBUF;
1044	Infotrees = (Infotree *) newinfotrees(Numtree, buftree);
1045	Lklptrn = NEW_LPMATRIX(Numtree, Numptrn);
1046	Ctree = (Tree *) new_stree(Maxspc, Maxibrnch, Numptrn, Seqconint);
1047	if (Debug) prcurtree(Ctree); /* */
1048	Alklptrn = Lklptrn[Cnotree];
1049	#if 0
1050	fmlength(Ctree, Distanmat, Maxspc);
1051	if (Aprox_optn) {
1052	initpartlkl(Ctree);
1053	aproxlkl(Ctree);
1054	praproxlkl(Ctree);
1055	putctopology(Ctree);
1056	} else {
1057	initpartlkl(Ctree);
1058	rp = (Node *)mlikelihood(Ctree);
1059	mlvalue(Ctree, Infotrees);
1060	if (Debug_optn) putctopology(Ctree);
1061	prtopology(Ctree);
1062	resulttree(Ctree);
1063	}
1064	#endif
1065
1066	#if DISLKL
1067	initpartlkl(Ctree);
1068	dislklhd = new_dmatrix(Numspc, Numspc);
1069	for (i = 0; i < Numspc - 1; i++) {
1070	iprb = Ctree->ebrnchp[i]->kinp->iprob;
1071	for (j = i + 1; j < Numspc; j++) {
1072	jprb = Ctree->ebrnchp[j]->kinp->iprob;
1073	for (k = 0, dislkl = 0.0; k < Numptrn; k++) {
1074	for (l = 0, dlkl = 0.0; l < Tpmradix; l++) {
1075	/* dlkl += fabs(iprb[k][l] - jprb[k][l]); */
1076	dlkl += Freqtpm[l] * iprb[k][l] * jprb[k][l];
1077	}
1078	/* dislkl += dlkl * Weight[k]; */
1079	dislkl += log(dlkl) * Weight[k];
1080	}
1081	dislkl = -dislkl;
1082	dislkl /= Numsite;
1083	dislklhd[i][j] = dislkl;
1084	dislklhd[j][i] = dislkl;
1085	}
1086	dislklhd[i][i] = 0.0;
1087	}
1088	dislklhd[Numspc-1][Numspc-1] = 0.0;
1089	/*
1090	for (i = 0; i < Numspc; i++) {
1091	for (j = 0; j < Numspc; j++) {
1092	if (i != j) printf("%5.0f", dislklhd[i][j] * 1000.0);
1093	else printf("%5.3s", Identif[i]);
1094	} putchar('\n');
1095	} putchar('\n');
1096	*/
1097	printf("\n%d %d\n", Numspc, Numsite);
1098	for (i = 0; i < Numspc; i++) {
1099	printf("%s\n", Identif[i]);
1100	for (j = 0; j < Numspc; j++) {
1101	printf(" %15.12f", dislklhd[i][j]);
1102	if ((j + 1) % 5 == 0) putchar('\n');
1103	} if (j % 5 != 0) putchar('\n');
1104	} putchar('\n');
1105
1106	free_dmatrix(dislklhd);
1107	exit(1);
1108	#endif /* DISLKL */
1109
1110	#if 1
1111	xstardecomp(Ctree);
1112	#else
1113	stardecomp(Ctree, Maxibrnch);
1114	#endif
1115	pathing(Ctree);
1116	fmlength(Ctree, Distanmat, Maxspc);
1117	initpartlkl(Ctree);
1118	Alklptrn = Lklptrn[Cnotree];
1119	rp = (Node *)mlikelihood(Ctree);
1120	mlvalue(Ctree, Infotrees);
1121	if (Debug_optn) putctopology(Ctree);
1122	putchar('\n');
1123	prtopology(Ctree);
1124	resulttree(Ctree);
1125	pstree(Epsfp, Ctree);
1126	FREE_LPMATRIX(Lklptrn);
1127
1128	} else if (Njoin_optn) { /* NJ MODE */
1129
1130	Numspc = Maxspc;
1131	Numbrnch = Maxspc;
1132	Maxibrnch = Maxspc - 3;
1133	Numibrnch = 0;
1134	Cnotree = 0;
1135	Numtree = 1;
1136	Infotrees = (Infotree *) newinfotrees(Numtree, buftree);
1137	Lklptrn = NEW_LPMATRIX(Numtree, Numptrn);
1138	Alklptrn = Lklptrn[0];
1139	puts("NJ");
1140
1141	Ctree = (Tree *) new_njtree(Maxspc, Maxibrnch, Numptrn, Seqconint);
1142	#if NJMLD
1143	njmtree(Ctree, Distanmat, Numspc, TRUE);
1144	#else
1145	enjtree(Ctree, Distanmat, Numspc, TRUE);
1146	#endif
1147	Alklptrn = Lklptrn[Cnotree];
1148	initpartlkl(Ctree);
1149	rp = (Node *)mlikelihood(Ctree);
1150	mlvalue(Ctree, Infotrees);
1151	if (Relia_optn) {
1152	Relitrif = new_dvector(Maxibrnch);
1153	Relistat = new_ivector(Numibrnch);
1154	Reliprob = new_dmatrix(Numibrnch, 2);
1155	Relinum = new_imatrix(Numibrnch, 2);
1156	qlrsearch(Ctree);
1157	putchar('\n');
1158	}
1159	prtopology(Ctree);
1160	resulttree(Ctree);
1161	/*
1162	initpartlkl(Ctree);
1163	rp = (Node *)mlikelihood(Ctree);
1164	mlvalue(Ctree, Infotrees);
1165	if (Debug_optn) putctopology(Ctree);
1166	putchar('\n');
1167	prtopology(Ctree);
1168	resulttree(Ctree);
1169	putchar('\n');
1170	putctopology(Ctree);
1171	*/
1172	pstree(Epsfp, Ctree);
1173	fprintf(Tplfp, "%d\n", Numtree);
1174	fputctopology(Tplfp, Ctree);
1175	fprintf(Trefp, "%d\n", Numtree);
1176	fputcphylogeny(Trefp, Ctree);
1177	FREE_LPMATRIX(Lklptrn);
1178	if (Relia_optn) {
1179	free_dvector(Relitrif);
1180	free_ivector(Relistat);
1181	free_dmatrix(Reliprob);
1182	free_imatrix(Relinum);
1183	}
1184	/*
1185	sorttree(Ctree, Ctree->rootp);
1186	prtopology(Ctree);
1187	putsortseq(Ctree);
1188	*/
1189
1190	} else { /* quick MODE */
1191
1192	Numtree = 2;
1193	Lklptrn = NEW_LPMATRIX(Numqltree, Numptrn);
1194	Ctree = (Tree *) new_tree(Maxspc, Maxibrnch, Numptrn, Seqconint);
1195	Strtree = new_cvector(buftree);
1196	Infoqltrees = (Infoqltree *) newinfoqltrees(MAXQLBUF, Maxbrnch);
1197	Qhead = (Infoqltree *) malloc(sizeof(Infoqltree));
1198	if (Qhead == NULL) maerror("Qhead in pml().");
1199	Qhead->lklaprox = 0.0;
1200	Qhead->residual = 0.0;
1201	for (Cnotree = 0; Cnotree < Numqltree; Cnotree++) {
1202	Alklptrn = Lklptrn[Cnotree];
1203	if (Cnotree == 0)
1204	initturn(Ctree);
1205	else
1206	randturn(Ctree);
1207	qtreeinit(Ctree);
1208	for (cnospc = 3; cnospc < Maxspc; cnospc++) {
1209	Numspc = cnospc + 1;
1210	Numibrnch = Numspc - 3;
1211	Numbrnch = Numspc + Numibrnch;
1212	Numpair = (Numspc * (Numspc - 1)) / 2;
1213	convertdistan(Ctree, cnospc+1, Distanmat, Distanvec);
1214	Qtail = Qhead;
1215	roundtree(Ctree, cnospc, Infoqltrees, Qhead, Qtail);
1216	}
1217	#if 0
1218	pathing(Ctree);
1219	fmlength(Ctree, Distanmat, Maxspc);
1220	initpartlkl(Ctree);
1221	rp = (Node *)mlikelihood(Ctree);
1222	/* mlvalue(Ctree, Infotrees); */
1223	prtopology(Ctree);
1224	resulttree(Ctree);
1225	#endif
1226	/* rerootq(Ctree, Numspc); */
1227
1228	for (i = 0; Ctree->bturn[i] != Numspc-1 && i < Numspc; i++)
1229	;
1230	reroot(Ctree, Ctree->ebrnchp[i]->kinp);
1231
1232	if (Aprox_optn) initpartlkl(Ctree);
1233	if (Aprox_optn) praproxlkl(Ctree);
1234	if (Aprox_optn) putctopology(Ctree);
1235	if (Cnotree == 0) {
1236	topaddtree = (Infoaddtree *)newinfoaddtree(buftree);
1237	strctree(Ctree, topaddtree->ltplgy);
1238	topaddtree->lklaprox = Ctree->lklhd;
1239	topaddtree->frequency = 1;
1240	cltree = (char )malloc((unsigned)buftree sizeof(char));
1241	if (cltree == NULL) maerror("cltree in protml.c");
1242	if (Debug) puts(topaddtree->ltplgy);
1243	Numaddtree = 1;
1244	} else {
1245	strctree(Ctree, cltree);
1246	lkla = Ctree->lklhd;
1247	insp = topaddtree;
1248	for (curp = topaddtree; curp != NULL; curp = curp->dp) {
1249	if (strcmp(cltree, curp->ltplgy) == 0) { /* same */
1250	curp->frequency++;
1251	break;
1252	}
1253	if (lkla < curp->lklaprox) insp = curp;
1254	}
1255	if (curp == NULL) {
1256	curp = (Infoaddtree *)newinfoaddtree(buftree);
1257	curp->dp = insp->dp;
1258	insp->dp = curp;
1259	(void)strcpy(curp->ltplgy, cltree);
1260	curp->lklaprox = lkla;
1261	curp->frequency = 1;
1262	Numaddtree++;
1263	}
1264	}
1265	}
1266	if (!Aprox_optn) tableaddtree(topaddtree, Numaddtree);
1267	FREE_LPMATRIX(Lklptrn);
1268	}
1269
1270	free_cvector(Comment);
1271	free_cmatrix(Identif);
1272	free_cmatrix(Sciname);
1273	free_cmatrix(Engname);
1274	free_imatrix(Seqconint);
1275	free_ivector(Weight);
1276	free_cmatrix(Seqchar);
1277	free_dvector(Distanvec);
1278	free_dmatrix(Distanmat);
1279	free_dvector(Brnlength);
1280	}
1281	#endif /* TPM */

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

Download in other formats: