Context Navigation

source: branches/items/GDE/TREEPUZZLE/src/puzzle1.c

Visit:

Last change on this file was 4932, checked in by baderk, 17 years ago
Fixed various typos…
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 140.4 KB

Line
1	/*
2	* puzzle1.c
3	*
4	*
5	* Part of TREE-PUZZLE 5.0 (June 2000)
6	*
7	* (c) 1999-2000 by Heiko A. Schmidt, Korbinian Strimmer,
8	* M. Vingron, and Arndt von Haeseler
9	* (c) 1995-1999 by Korbinian Strimmer and Arndt von Haeseler
10	*
11	* All parts of the source except where indicated are distributed under
12	* the GNU public licence. See http://www.opensource.org for details.
13	*/
14
15
16	#define EXTERN
17
18	#include "puzzle.h"
19	#include "gamma.h"
20
21	void num2quart(uli qnum, int a, int b, int c, int d)
22	{
23	double temp;
24	uli aa, bb, cc, dd;
25	uli lowval=0, highval=0;
26
27	aa=0; bb=1; cc=2; dd=3;
28
29	temp = (double)(24 * qnum);
30	temp = sqrt(temp);
31	temp = sqrt(temp);
32	/* temp = pow(temp, (double)(1/4)); */
33	dd = (uli) floor(temp) + 1;
34	if (dd < 3) dd = 3;
35	lowval = (uli) dd(dd-1)(dd-2)*(dd-3)/24;
36	highval = (uli) (dd+1)dd(dd-1)*(dd-2)/24;
37	if (lowval >= qnum)
38	while ((lowval > qnum)) {
39	dd -= 1; lowval = (uli) dd(dd-1)(dd-2)*(dd-3)/24;
40	}
41	else {
42	while (highval <= qnum) {
43	dd += 1; highval = (uli) (dd+1)dd(dd-1)*(dd-2)/24;
44	}
45	lowval = (uli) dd(dd-1)(dd-2)*(dd-3)/24;
46	}
47	qnum -= lowval;
48	if (qnum > 0) {
49	temp = (double)(6 * qnum);
50	temp = pow(temp, (double)(1/3));
51	cc = (uli) floor(temp);
52	if (cc < 2) cc= 2;
53	lowval = (uli) cc(cc-1)(cc-2)/6;
54	highval = (uli) (cc+1)cc(cc-1)/6;
55	if (lowval >= qnum)
56	while ((lowval > qnum)) {
57	cc -= 1; lowval = (uli) cc(cc-1)(cc-2)/6;
58	}
59	else {
60	while (highval <= qnum) {
61	cc += 1; highval = (uli) (cc+1)cc(cc-1)/6;
62	}
63	lowval = (uli) cc(cc-1)(cc-2)/6;
64	}
65	qnum -= lowval;
66	if (qnum > 0) {
67	temp = (double)(2 * qnum);
68	temp = sqrt(temp);
69	bb = (uli) floor(temp);
70	if (bb < 1) bb= 1;
71	lowval = (uli) bb*(bb-1)/2;
72	highval = (uli) (bb+1)*bb/2;
73	if (lowval >= qnum)
74	while ((lowval > qnum)) {
75	bb -= 1; lowval = (uli) bb*(bb-1)/2;
76	}
77	else {
78	while (highval <= qnum) {
79	bb += 1; highval = (uli) (bb+1)*bb/2;
80	}
81	lowval = (uli) bb*(bb-1)/2;
82	}
83	qnum -= lowval;
84	if (qnum > 0) {
85	aa = (uli) qnum;
86	if (aa < 0) aa= 0;
87	}
88	}
89	}
90	*d = (int)dd;
91	*c = (int)cc;
92	*b = (int)bb;
93	*a = (int)aa;
94	} /* num2quart */
95
96	/******************/
97
98	uli numquarts(int maxspc)
99	{
100	uli tmp;
101	int a, b, c, d;
102
103	if (maxspc < 4)
104	return (uli)0;
105	else {
106	maxspc--;
107	a = maxspc-3;
108	b = maxspc-2;
109	c = maxspc-1;
110	d = maxspc;
111
112	tmp = (uli) 1 + a +
113	(uli) b * (b-1) / 2 +
114	(uli) c * (c-1) * (c-2) / 6 +
115	(uli) d * (d-1) * (d-2) * (d-3) / 24;
116	return (tmp);
117	}
118	} /* numquarts */
119
120	/******************/
121
122	uli quart2num (int a, int b, int c, int d)
123	{
124	uli tmp;
125	if ((a>b) \|\| (b>c) \|\| (c>d)) {
126	fprintf(stderr, "Error PP5 not (%d <= %d <= %d <= %d) !!!\n", a, b, c,
127	d);
128	exit (1);
129	}
130	tmp = (uli) a +
131	(uli) b * (b-1) / 2 +
132	(uli) c * (c-1) * (c-2) / 6 +
133	(uli) d * (d-1) * (d-2) * (d-3) / 24;
134	return (tmp);
135	} /* quart2num */
136
137	/******************/
138
139
140
141	/* flag=0 old allquart binary */
142	/* flag=1 allquart binary */
143	/* flag=2 allquart ACSII */
144	/* flag=3 quartlh binary */
145	/* flag=4 quartlh ASCII */
146
147	void writetpqfheader(int nspec,
148	FILE *ofp,
149	int flag)
150	{ int currspec;
151
152	if (flag == 0) {
153	unsigned long nquart;
154	unsigned long blocklen;
155
156	nquart = numquarts(nspec);
157	/* compute number of bytes */
158	if (nquart % 2 == 0) { /* even number */
159	blocklen = (nquart)/2;
160	} else { /* odd number */
161	blocklen = (nquart + 1)/2;
162	}
163	/* FPRINTF(STDOUTFILE "Writing quartet file: %s\n", filename); */
164	fprintf(ofp, "TREE-PUZZLE\n%s\n\n", VERSION);
165	fprintf(ofp, "species: %d\n", nspec);
166	fprintf(ofp, "quartets: %lu\n", nquart);
167	fprintf(ofp, "bytes: %lu\n\n", blocklen);
168
169
170	/* fwrite(&(quartetinfo[0]), sizeof(char), blocklen, ofp); */
171	}
172
173	if (flag == 1) fprintf(ofp, "##TPQF-BB (TREE-PUZZLE %s)\n%d\n", VERSION, nspec);
174	if (flag == 2) fprintf(ofp, "##TPQF-BA (TREE-PUZZLE %s)\n%d\n", VERSION, nspec);
175	if (flag == 3) fprintf(ofp, "##TPQF-LB (TREE-PUZZLE %s)\n%d\n", VERSION, nspec);
176	if (flag == 4) fprintf(ofp, "##TPQF-LA (TREE-PUZZLE %s)\n%d\n", VERSION, nspec);
177
178	for (currspec=0; currspec<nspec; currspec++) {
179	fputid(ofp, currspec);
180	fprintf(ofp, "\n");
181	} /* for each species */
182	fprintf(ofp, "\n");
183
184	} /* writetpqfheader */
185
186
187
188	void writeallquarts(int nspec,
189	char *filename,
190	unsigned char *quartetinfo)
191	{ unsigned long nquart;
192	unsigned long blocklen;
193	FILE *ofp;
194
195	nquart = numquarts(nspec);
196
197	/* compute number of bytes */
198	if (nquart % 2 == 0) { /* even number */
199	blocklen = (nquart)/2;
200	} else { /* odd number */
201	blocklen = (nquart + 1)/2;
202	}
203
204	FPRINTF(STDOUTFILE "Writing quartet file: %s\n", filename);
205
206	openfiletowrite(&ofp, filename, "all quartets");
207
208	writetpqfheader(nspec, ofp, 0);
209
210	fwrite(&(quartetinfo[0]), sizeof(char), blocklen, ofp);
211	fclose(ofp);
212
213	} /* writeallquart */
214
215
216
217	void readallquarts(int nspec,
218	char *filename,
219	unsigned char *quartetinfo)
220	{ unsigned long nquart;
221	unsigned long blocklen;
222	int currspec;
223	unsigned long dummynquart;
224	unsigned long dummyblocklen;
225	int dummynspec;
226	char dummyversion[128];
227	char dummyname[128];
228	FILE *ifp;
229
230	nquart = numquarts(nspec);
231
232	/* compute number of bytes */
233	if (nquart % 2 == 0) { /* even number */
234	blocklen = (nquart)/2;
235	} else { /* odd number */
236	blocklen = (nquart + 1)/2;
237	}
238
239	/* &(quartetinfo[0] */
240
241	openfiletoread(&ifp, filename, "all quartets");
242
243	fscanf(ifp, "TREE-PUZZLE\n");
244	fscanf(ifp, "%s\n\n", dummyversion);
245
246	fscanf(ifp, "species: %d\n", &dummynspec);
247	fscanf(ifp, "quartets: %lu\n", &dummynquart);
248	fscanf(ifp, "bytes: %lu\n\n", &dummyblocklen);
249
250	if ((nspec != dummynspec) \|\|
251	(nquart != dummynquart) \|\|
252	(blocklen != dummyblocklen)) {
253	FPRINTF(STDOUTFILE "ERROR: Parameters in quartet file %s do not match!!!\n",
254	filename);
255	# if PARALLEL
256	PP_SendDone();
257	MPI_Finalize();
258	# endif /* PARALLEL */
259	exit(1);
260	}
261
262	FPRINTF(STDOUTFILE "Reading quartet file: %s\n", filename);
263	FPRINTF(STDOUTFILE " generated by: TREE-PUZZLE %s\n", dummyversion);
264	FPRINTF(STDOUTFILE " current: TREE-PUZZLE %s\n", VERSION);
265
266	FPRINTF(STDOUTFILE " %d species, %lu quartets, %lu bytes\n",
267	dummynspec, dummynquart, dummyblocklen);
268
269	for (currspec=0; currspec<nspec; currspec++) {
270
271	fscanf(ifp, "%s\n", dummyname);
272	FPRINTF(STDOUTFILE " %3d. %s (", currspec+1, dummyname);
273	fputid(STDOUT, currspec);
274	FPRINTF(STDOUTFILE ")\n");
275
276	} /* for each species */
277	FPRINTF(STDOUTFILE "\n");
278
279	fread(&(quartetinfo[0]), sizeof(char), blocklen, ifp);
280	fclose(ifp);
281
282	} /* writeallquart */
283
284
285
286
287	/******************************************************************************/
288	/* options - file I/O - output */
289	/******************************************************************************/
290
291	/* compute TN parameters according to F84 Ts/Tv ratio */
292	void makeF84model()
293	{
294	double rho, piA, piC, piG, piT, piR, piY, ts, yr;
295
296	piA = Freqtpm[0];
297	piC = Freqtpm[1];
298	piG = Freqtpm[2];
299	piT = Freqtpm[3];
300	piR = piA + piG;
301	piY = piC + piT;
302	if (piCpiTpiR + piApiGpiY == 0.0) {
303	FPRINTF(STDOUTFILE "\n\n\nF84 model not possible ");
304	FPRINTF(STDOUTFILE "(bad combination of base frequencies)\n");
305	tstvf84 = 0.0;
306	return;
307	}
308	rho = (piRpiY(piRpiYtstvf84 - (piApiG + piCpiT)))/
309	(piCpiTpiR + piApiGpiY);
310
311	if(piR == 0.0 \|\| piY == 0.0 \|\| (piR + rho) == 0.0) {
312	FPRINTF(STDOUTFILE "\n\n\nF84 model not possible ");
313	FPRINTF(STDOUTFILE "(bad combination of base frequencies)\n");
314	tstvf84 = 0.0;
315	return;
316	}
317	ts = 0.5 + 0.25rho(1.0/piR + 1.0/piY);
318	yr = (piY + rho)/piY * piR/(piR + rho);
319	if (ts < MINTS \|\| ts > MAXTS) {
320	FPRINTF(STDOUTFILE "\n\n\nF84 model not possible ");
321	FPRINTF(STDOUTFILE "(bad Ts/Tv parameter)\n");
322	tstvf84 = 0.0;
323	return;
324	}
325	if (yr < MINYR \|\| yr > MAXYR) {
326	FPRINTF(STDOUTFILE "\n\n\nF84 model not possible ");
327	FPRINTF(STDOUTFILE "(bad Y/R transition parameter)\n");
328	tstvf84 = 0.0;
329	return;
330	}
331	TSparam = ts;
332	YRparam = yr;
333	optim_optn = FALSE;
334	}
335
336	/* compute number of quartets used in LM analysis */
337	void compnumqts()
338	{
339	if (lmqts == 0) {
340	if (numclust == 4)
341	Numquartets = (uli) clustAclustBclustC*clustD;
342	if (numclust == 3)
343	Numquartets = (uli) clustAclustBclustC*(clustC-1)/2;
344	if (numclust == 2)
345	Numquartets = (uli) clustA(clustA-1)/2 clustB*(clustB-1)/2;
346	if (numclust == 1)
347	Numquartets = (uli) Maxspc(Maxspc-1)(Maxspc-2)*(Maxspc-3)/24;
348	} else {
349	Numquartets = lmqts;
350	}
351	}
352
353	/* set options interactively */
354	void setoptions()
355	{
356	int i, valid;
357	double sumfreq;
358	char ch;
359
360	/* defaults */
361	rhetmode = UNIFORMRATE; /* assume rate homogeneity */
362	numcats = 1;
363	Geta = 0.05;
364	grate_optim = FALSE;
365	fracinv = 0.0;
366	fracinv_optim = FALSE;
367
368	compclock = FALSE; /* compute clocklike branch lengths */
369	locroot = -1; /* search for optimal place of root */
370	qcalg_optn = FALSE; /* don't use sampling of quartets */
371	approxp_optn = TRUE; /* approximate parameter estimates */
372	listqptrees = PSTOUT_NONE; /* list puzzling step trees */
373
374	/* approximate QP quartets? */
375	if (Maxspc <= 6) approxqp = FALSE;
376	else approxqp = TRUE;
377
378	codon_optn = 0; /* use all positions in a codon */
379
380	/* number of puzzling steps */
381	if (Maxspc <= 25) Numtrial = 1000;
382	else if (Maxspc <= 50) Numtrial = 10000;
383	else if (Maxspc <= 75) Numtrial = 25000;
384	else Numtrial = 50000;
385
386	utree_optn = TRUE; /* use first user tree for estimation */
387	outgroup = 0; /* use first taxon as outgroup */
388	sym_optn = FALSE; /* symmetrize doublet frequencies */
389	tstvf84 = 0.0; /* disable F84 model */
390	show_optn = FALSE; /* show unresolved quartets */
391	typ_optn = TREERECON_OPTN; /* tree reconstruction */
392	numclust = 1; /* one clusters in LM analysis */
393	lmqts = 0; /* all quartets in LM analysis */
394	compnumqts();
395	if (Numquartets > 10000) {
396	lmqts = 10000; /* 10000 quartets in LM analysis */
397	compnumqts();
398	}
399
400	do {
401	FPRINTF(STDOUTFILE "\n\n\nGENERAL OPTIONS\n");
402	FPRINTF(STDOUTFILE " b Type of analysis? ");
403	if (typ_optn == TREERECON_OPTN) FPRINTF(STDOUTFILE "Tree reconstruction\n");
404	if (typ_optn == LIKMAPING_OPTN) FPRINTF(STDOUTFILE "Likelihood mapping\n");
405	if (typ_optn == TREERECON_OPTN) {
406	FPRINTF(STDOUTFILE " k Tree search procedure? ");
407	if (puzzlemode == QUARTPUZ) FPRINTF(STDOUTFILE "Quartet puzzling\n");
408	if (puzzlemode == USERTREE) FPRINTF(STDOUTFILE "User defined trees\n");
409	if (puzzlemode == PAIRDIST) FPRINTF(STDOUTFILE "Pairwise distances only (no tree)\n");
410	if (puzzlemode == QUARTPUZ) {
411	FPRINTF(STDOUTFILE " v Approximate quartet likelihood? %s\n",
412	(approxqp ? "Yes" : "No"));
413	FPRINTF(STDOUTFILE " u List unresolved quartets? %s\n",
414	(show_optn ? "Yes" : "No"));
415	FPRINTF(STDOUTFILE " n Number of puzzling steps? %lu\n",
416	Numtrial);
417	FPRINTF(STDOUTFILE " j List puzzling step trees? ");
418	switch (listqptrees) {
419	case PSTOUT_NONE: FPRINTF(STDOUTFILE "No\n"); break;
420	case PSTOUT_ORDER: FPRINTF(STDOUTFILE "Unique topologies\n"); break;
421	case PSTOUT_LISTORDER: FPRINTF(STDOUTFILE "Unique topologies & Chronological list\n"); break;
422	case PSTOUT_LIST: FPRINTF(STDOUTFILE "Chronological list only\n"); break;
423	}
424
425	FPRINTF(STDOUTFILE " o Display as outgroup? ");
426	fputid(STDOUT, outgroup);
427	FPRINTF(STDOUTFILE "\n");
428	}
429	if (puzzlemode == QUARTPUZ \|\| puzzlemode == USERTREE) {
430	FPRINTF(STDOUTFILE " z Compute clocklike branch lengths? ");
431	if (compclock) FPRINTF(STDOUTFILE "Yes\n");
432	else FPRINTF(STDOUTFILE "No\n");
433	}
434	if (compclock)
435	if (puzzlemode == QUARTPUZ \|\| puzzlemode == USERTREE) {
436	FPRINTF(STDOUTFILE " l Location of root? ");
437	if (locroot < 0) FPRINTF(STDOUTFILE "Best place (automatic search)\n");
438	else if (locroot < Maxspc) {
439	FPRINTF(STDOUTFILE "Branch %d (", locroot + 1);
440	fputid(STDOUT, locroot);
441	FPRINTF(STDOUTFILE ")\n");
442	} else FPRINTF(STDOUTFILE "Branch %d (internal branch)\n", locroot + 1);
443	}
444	}
445	if (typ_optn == LIKMAPING_OPTN) {
446	FPRINTF(STDOUTFILE " g Group sequences in clusters? ");
447	if (numclust == 1) FPRINTF(STDOUTFILE "No\n");
448	else FPRINTF(STDOUTFILE "Yes (%d clusters as specified)\n", numclust);
449	FPRINTF(STDOUTFILE " n Number of quartets? ");
450	if (lmqts == 0) FPRINTF(STDOUTFILE "%lu (all possible)\n", Numquartets);
451	else FPRINTF(STDOUTFILE "%lu (random choice)\n", lmqts);
452	}
453	FPRINTF(STDOUTFILE " e Parameter estimates? ");
454	if (approxp_optn) FPRINTF(STDOUTFILE "Approximate (faster)\n");
455	else FPRINTF(STDOUTFILE "Exact (slow)\n");
456	if (!(puzzlemode == USERTREE && typ_optn == TREERECON_OPTN)) {
457	FPRINTF(STDOUTFILE " x Parameter estimation uses? ");
458	if (qcalg_optn) FPRINTF(STDOUTFILE "Quartet sampling + NJ tree\n");
459	else FPRINTF(STDOUTFILE "Neighbor-joining tree\n");
460
461	} else {
462	FPRINTF(STDOUTFILE " x Parameter estimation uses? ");
463	if (utree_optn)
464	FPRINTF(STDOUTFILE "1st input tree\n");
465	else if (qcalg_optn) FPRINTF(STDOUTFILE "Quartet sampling + NJ tree\n");
466	else FPRINTF(STDOUTFILE "Neighbor-joining tree\n");
467	}
468	FPRINTF(STDOUTFILE "SUBSTITUTION PROCESS\n");
469	FPRINTF(STDOUTFILE " d Type of sequence input data? ");
470	if (auto_datatype == AUTO_GUESS) FPRINTF(STDOUTFILE "Auto: ");
471	if (data_optn == NUCLEOTIDE) FPRINTF(STDOUTFILE "Nucleotides\n");
472	if (data_optn == AMINOACID) FPRINTF(STDOUTFILE "Amino acids\n");
473	if (data_optn == BINARY) FPRINTF(STDOUTFILE "Binary states\n");
474	if (data_optn == NUCLEOTIDE && (Maxseqc % 3) == 0 && !SH_optn) {
475	FPRINTF(STDOUTFILE " h Codon positions selected? ");
476	if (codon_optn == 0) FPRINTF(STDOUTFILE "Use all positions\n");
477	if (codon_optn == 1) FPRINTF(STDOUTFILE "Use only 1st positions\n");
478	if (codon_optn == 2) FPRINTF(STDOUTFILE "Use only 2nd positions\n");
479	if (codon_optn == 3) FPRINTF(STDOUTFILE "Use only 3rd positions\n");
480	if (codon_optn == 4) FPRINTF(STDOUTFILE "Use 1st and 2nd positions\n");
481	}
482	FPRINTF(STDOUTFILE " m Model of substitution? ");
483	if (data_optn == NUCLEOTIDE) { /* nucleotides */
484	if (nuc_optn) {
485	if(HKY_optn)
486	FPRINTF(STDOUTFILE "HKY (Hasegawa et al. 1985)\n");
487	else {
488	FPRINTF(STDOUTFILE "TN (Tamura-Nei 1993)\n");
489	FPRINTF(STDOUTFILE " p Constrain TN model to F84 model? ");
490	if (tstvf84 == 0.0)
491	FPRINTF(STDOUTFILE "No\n");
492	else FPRINTF(STDOUTFILE "Yes (Ts/Tv ratio = %.2f)\n", tstvf84);
493	}
494	FPRINTF(STDOUTFILE " t Transition/transversion parameter? ");
495	if (optim_optn)
496	FPRINTF(STDOUTFILE "Estimate from data set\n");
497	else
498	FPRINTF(STDOUTFILE "%.2f\n", TSparam);
499	if (TN_optn) {
500	FPRINTF(STDOUTFILE " r Y/R transition parameter? ");
501	if (optim_optn)
502	FPRINTF(STDOUTFILE "Estimate from data set\n");
503	else
504	FPRINTF(STDOUTFILE "%.2f\n", YRparam);
505	}
506	}
507	if (SH_optn) {
508	FPRINTF(STDOUTFILE "SH (Schoeniger-von Haeseler 1994)\n");
509	FPRINTF(STDOUTFILE " t Transition/transversion parameter? ");
510	if (optim_optn)
511	FPRINTF(STDOUTFILE "Estimate from data set\n");
512	else
513	FPRINTF(STDOUTFILE "%.2f\n", TSparam);
514	}
515	}
516	if (data_optn == NUCLEOTIDE && SH_optn) {
517	FPRINTF(STDOUTFILE " h Doublets defined by? ");
518	if (SHcodon)
519	FPRINTF(STDOUTFILE "1st and 2nd codon positions\n");
520	else
521	FPRINTF(STDOUTFILE "1st+2nd, 3rd+4th, etc. site\n");
522	}
523	if (data_optn == AMINOACID) { /* amino acids */
524	switch (auto_aamodel) {
525	case AUTO_GUESS:
526	FPRINTF(STDOUTFILE "Auto: ");
527	break;
528	case AUTO_DEFAULT:
529	FPRINTF(STDOUTFILE "Def.: ");
530	break;
531	}
532	if (Dayhf_optn) FPRINTF(STDOUTFILE "Dayhoff (Dayhoff et al. 1978)\n");
533	if (Jtt_optn) FPRINTF(STDOUTFILE "JTT (Jones et al. 1992)\n");
534	if (mtrev_optn) FPRINTF(STDOUTFILE "mtREV24 (Adachi-Hasegawa 1996)\n");
535	if (cprev_optn) FPRINTF(STDOUTFILE "cpREV45 (Adachi et al. 2000)\n");
536	if (blosum62_optn) FPRINTF(STDOUTFILE "BLOSUM62 (Henikoff-Henikoff 92)\n");
537	if (vtmv_optn) FPRINTF(STDOUTFILE "VT (Mueller-Vingron 2000)\n");
538	if (wag_optn) FPRINTF(STDOUTFILE "WAG (Whelan-Goldman 2000)\n");
539	}
540	if (data_optn == BINARY) { /* binary states */
541	FPRINTF(STDOUTFILE "Two-state model (Felsenstein 1981)\n");
542	}
543	if (data_optn == AMINOACID)
544	FPRINTF(STDOUTFILE " f Amino acid frequencies? ");
545	else if (data_optn == NUCLEOTIDE && SH_optn)
546	FPRINTF(STDOUTFILE " f Doublet frequencies? ");
547	else if (data_optn == NUCLEOTIDE && nuc_optn)
548	FPRINTF(STDOUTFILE " f Nucleotide frequencies? ");
549	else if (data_optn == BINARY)
550	FPRINTF(STDOUTFILE " f Binary state frequencies? ");
551	FPRINTF(STDOUTFILE "%s\n", (Frequ_optn ? "Estimate from data set" :
552	"Use specified values"));
553	if (data_optn == NUCLEOTIDE && SH_optn)
554	FPRINTF(STDOUTFILE " s Symmetrize doublet frequencies? %s\n",
555	(sym_optn ? "Yes" : "No"));
556
557	FPRINTF(STDOUTFILE "RATE HETEROGENEITY\n");
558	FPRINTF(STDOUTFILE " w Model of rate heterogeneity? ");
559	if (rhetmode == UNIFORMRATE) FPRINTF(STDOUTFILE "Uniform rate\n");
560	if (rhetmode == GAMMARATE ) FPRINTF(STDOUTFILE "Gamma distributed rates\n");
561	if (rhetmode == TWORATE ) FPRINTF(STDOUTFILE "Two rates (1 invariable + 1 variable)\n");
562	if (rhetmode == MIXEDRATE ) FPRINTF(STDOUTFILE "Mixed (1 invariable + %d Gamma rates)\n", numcats);
563
564	if (rhetmode == TWORATE \|\| rhetmode == MIXEDRATE) {
565	FPRINTF(STDOUTFILE " i Fraction of invariable sites? ");
566	if (fracinv_optim) FPRINTF(STDOUTFILE "Estimate from data set");
567	else FPRINTF(STDOUTFILE "%.2f", fracinv);
568	if (fracinv == 0.0 && !fracinv_optim) FPRINTF(STDOUTFILE " (all sites variable)");
569	FPRINTF(STDOUTFILE "\n");
570	}
571	if (rhetmode == GAMMARATE \|\| rhetmode == MIXEDRATE) {
572	FPRINTF(STDOUTFILE " a Gamma distribution parameter alpha? ");
573	if (grate_optim)
574	FPRINTF(STDOUTFILE "Estimate from data set\n");
575	else if (Geta > 0.5)
576	FPRINTF(STDOUTFILE "%.2f (strong rate heterogeneity)\n", (1.0-Geta)/Geta);
577	else FPRINTF(STDOUTFILE "%.2f (weak rate heterogeneity)\n", (1.0-Geta)/Geta);
578	FPRINTF(STDOUTFILE " c Number of Gamma rate categories? %d\n", numcats);
579	}
580
581	FPRINTF(STDOUTFILE "\nQuit [q], confirm [y], or change [menu] settings: ");
582
583	/* read one char */
584	ch = getchar();
585	if (ch != '\n') {
586	do ;
587	while (getchar() != '\n');
588	}
589	ch = (char) tolower((int) ch);
590
591	/* letters in use: a b c d e f g h i j k l m n o p q r s t u v w y x z */
592	/* letters not in use: */
593
594	switch (ch) {
595
596	case '\n': break;
597
598	case 'z': if (typ_optn == TREERECON_OPTN && (puzzlemode == QUARTPUZ \|\| puzzlemode == USERTREE)) {
599	compclock = compclock + 1;
600	if (compclock == 2) compclock = 0;
601	} else {
602	FPRINTF(STDOUTFILE "\n\n\nThis is not a possible option!\n");
603	}
604	break;
605
606	case 'l': if (compclock && typ_optn == TREERECON_OPTN && (puzzlemode == QUARTPUZ \|\| puzzlemode == USERTREE)) {
607	FPRINTF(STDOUTFILE "\n\n\nEnter an invalid branch number to search ");
608	FPRINTF(STDOUTFILE "for the best location!\n");
609	FPRINTF(STDOUTFILE "\nPlace root at branch (1-%d): ",
610	2*Maxspc-3);
611	scanf("%d", &locroot);
612	do ;
613	while (getchar() != '\n');
614	if (locroot < 1 \|\| locroot > 2*Maxspc-3) locroot = 0;
615	locroot = locroot - 1;
616	} else {
617	FPRINTF(STDOUTFILE "\n\n\nThis is not a possible option!\n");
618	}
619	break;
620
621	case 'e': if ((rhetmode == TWORATE \|\| rhetmode == MIXEDRATE) && fracinv_optim) {
622	FPRINTF(STDOUTFILE "\n\n\nInvariable sites estimation needs to be exact!\n");
623	} else {
624	approxp_optn = approxp_optn + 1;
625	if (approxp_optn == 2) approxp_optn = 0;
626	}
627	break;
628
629	case 'w': rhetmode = rhetmode + 1;
630	if (rhetmode == 4) rhetmode = UNIFORMRATE;
631	if (rhetmode == UNIFORMRATE) { /* uniform rate */
632	numcats = 1;
633	Geta = 0.05;
634	grate_optim = FALSE;
635	fracinv = 0.0;
636	fracinv_optim = FALSE;
637	}
638	if (rhetmode == GAMMARATE ) { /* Gamma distributed rates */
639	numcats = 8;
640	Geta = 0.05;
641	grate_optim = TRUE;
642	fracinv = 0.0;
643	fracinv_optim = FALSE;
644	}
645	if (rhetmode == TWORATE ) { /* two rates (1 invariable + 1 variable) */
646	approxp_optn = FALSE;
647	numcats = 1;
648	Geta = 0.05;
649	grate_optim = FALSE;
650	fracinv = 0.0;
651	fracinv_optim = TRUE;
652	}
653	if (rhetmode == MIXEDRATE ) { /* mixed (1 invariable + Gamma rates) */
654	approxp_optn = FALSE;
655	numcats = 8;
656	Geta = 0.05;
657	grate_optim = TRUE;
658	fracinv = 0.0;
659	fracinv_optim = TRUE;
660	}
661	break;
662
663	case 'i': if (rhetmode == TWORATE \|\| rhetmode == MIXEDRATE) {
664	FPRINTF(STDOUTFILE "\n\n\nEnter an invalid value for ");
665	FPRINTF(STDOUTFILE "estimation from data set!\n");
666	FPRINTF(STDOUTFILE "\nFraction of invariable sites among all sites (%.2f-%.2f): ",
667	MINFI, MAXFI);
668	scanf("%lf", &fracinv);
669	do ;
670	while (getchar() != '\n');
671	if (fracinv < MINFI \|\| fracinv > MAXFI) {
672	fracinv_optim = TRUE;
673	fracinv = 0.0;
674	} else {
675	fracinv_optim = FALSE;
676	}
677	} else {
678	FPRINTF(STDOUTFILE "\n\n\nThis is not a possible option!\n");
679	}
680	break;
681
682	case 'a': if (rhetmode == GAMMARATE \|\| rhetmode == MIXEDRATE) {
683	FPRINTF(STDOUTFILE "\n\n\nEnter an invalid value for estimation from data set!\n");
684	FPRINTF(STDOUTFILE "\nGamma distribution parameter alpha (%.2f-%.2f): ",
685	(1.0-MAXGE)/MAXGE, (1.0-MINGE)/MINGE);
686	scanf("%lf", &Geta);
687	do ;
688	while (getchar() != '\n');
689	if (Geta < (1.0-MAXGE)/MAXGE \|\| Geta > (1.0-MINGE)/MINGE) {
690	grate_optim = TRUE;
691	Geta = 0.05;
692	} else {
693	grate_optim = FALSE;
694	Geta = 1.0/(1.0 + Geta);
695	}
696	} else
697	FPRINTF(STDOUTFILE "\n\n\nThis is not a possible option!\n");
698	break;
699
700	case 'c': if (rhetmode == GAMMARATE \|\| rhetmode == MIXEDRATE) {
701	FPRINTF(STDOUTFILE "\n\n\nNumber of Gamma rate categories (%d-%d): ",
702	MINCAT, MAXCAT);
703	scanf("%d", &numcats);
704	do ;
705	while (getchar() != '\n');
706	if (numcats < MINCAT \|\| numcats > MAXCAT) {
707	FPRINTF(STDOUTFILE "\n\n\nThis number of categories is not available!\n");
708	numcats = 4;
709	}
710	} else {
711	FPRINTF(STDOUTFILE "\n\n\nThis is not a possible option!\n");
712	}
713	break;
714
715	case 'h': if (data_optn == NUCLEOTIDE && (Maxseqc % 3) == 0 && !SH_optn) {
716	codon_optn = codon_optn + 1;
717	if (codon_optn == 5) codon_optn = 0;
718	translatedataset();
719	/* reestimate nucleotide frequencies only
720	if user did not specify other values */
721	if (Frequ_optn) estimatebasefreqs();
722
723	} else if (data_optn == NUCLEOTIDE && SH_optn) {
724	if (Maxseqc % 2 != 0 && Maxseqc % 3 == 0) {
725	SHcodon = TRUE;
726	FPRINTF(STDOUTFILE "\n\n\nThis is the only possible option for the data set!\n");
727	}
728	if (Maxseqc % 3 != 0 && Maxseqc % 2 == 0) {
729	SHcodon = FALSE;
730	FPRINTF(STDOUTFILE "\n\n\nThis is the only possible option for the data set!\n");
731	}
732	if (Maxseqc % 2 == 0 && Maxseqc % 3 == 0) {
733	if (SHcodon)
734	SHcodon = FALSE;
735	else
736	SHcodon = TRUE;
737	translatedataset();
738	/* reestimate nucleotide frequencies only
739	if user did not specify other values */
740	if (Frequ_optn) estimatebasefreqs();
741	}
742	} else {
743	FPRINTF(STDOUTFILE "\n\n\nThis is not a possible option!\n");
744	}
745	break;
746
747	case 'x': if (typ_optn == TREERECON_OPTN && puzzlemode == USERTREE) {
748	if (utree_optn) {
749	utree_optn = FALSE;
750	qcalg_optn = FALSE;
751	} else {
752	qcalg_optn = qcalg_optn + 1;
753	if (qcalg_optn == 2) {
754	qcalg_optn = 0;
755	utree_optn = TRUE;
756	}
757	}
758	} else {
759	qcalg_optn = qcalg_optn + 1;
760	if (qcalg_optn == 2) qcalg_optn = 0;
761	}
762	break;
763
764	case 'k': if (typ_optn == TREERECON_OPTN) {
765	puzzlemode = (puzzlemode + 1) % 3;
766	/* puzzlemode = puzzlemode + 1;
767	if (puzzlemode == 3) puzzlemode = 0;
768	xxx */
769	} else {
770	FPRINTF(STDOUTFILE "\n\n\nThis is not a possible option!\n");
771	}
772	break;
773
774	case 'b': typ_optn = (typ_optn + 1) % 2;
775	/* typ_optn = typ_optn + 1;
776	if (typ_optn == 2) typ_optn = TREERECON_OPTN;
777	xxx */
778	break;
779
780	case 'g': if (typ_optn == LIKMAPING_OPTN) {
781	clustA = clustB = clustC = clustD = 0;
782	if (numclust != 1) {
783	numclust = 1;
784	} else {
785	FPRINTF(STDOUTFILE "\n\n\nNumber of clusters (2-4): ");
786	scanf("%d", &numclust);
787	do ;
788	while (getchar() != '\n');
789	if (numclust < 2 \|\| numclust > 4) {
790	numclust = 1;
791	FPRINTF(STDOUTFILE "\n\n\nOnly 2, 3, or 4 ");
792	FPRINTF(STDOUTFILE "clusters possible\n");
793	} else {
794	FPRINTF(STDOUTFILE "\nDistribute all sequences over the ");
795	if (numclust == 2) {
796	FPRINTF(STDOUTFILE "two clusters a and b (At least two\n");
797	FPRINTF(STDOUTFILE "sequences per cluster are necessary), ");
798	}
799	if (numclust == 3) {
800	FPRINTF(STDOUTFILE "three clusters a, b, and c\n");
801	FPRINTF(STDOUTFILE "(At least one sequence in cluster a and b, and at least two\n");
802	FPRINTF(STDOUTFILE "sequences in c are necessary), ");
803	}
804	if (numclust == 4) {
805	FPRINTF(STDOUTFILE "four clusters a, b, c, and d\n");
806	FPRINTF(STDOUTFILE "(At least one sequence per cluster is necessary),\n");
807	}
808	FPRINTF(STDOUTFILE "type x to exclude a sequence:\n\n");
809
810	for (i = 0; i < Maxspc; i++) {
811	valid = FALSE;
812	do {
813	fputid10(STDOUT, i);
814	FPRINTF(STDOUTFILE ": ");
815	/* read one char */
816	ch = getchar();
817	if (ch != '\n') {
818	do ;
819	while (getchar() != '\n');
820	}
821	ch = (char) tolower((int) ch);
822	if (ch == 'a' \|\| ch == 'b' \|\| ch == 'x')
823	valid = TRUE;
824	if (numclust == 3 \|\| numclust == 4)
825	if (ch == 'c') valid = TRUE;
826	if (numclust == 4)
827	if (ch == 'd') valid = TRUE;
828	} while (!valid);
829	if (ch == 'a') {
830	clusterA[clustA] = i;
831	clustA++;
832	}
833	if (ch == 'b') {
834	clusterB[clustB] = i;
835	clustB++;
836	}
837	if (ch == 'c') {
838	clusterC[clustC] = i;
839	clustC++;
840	}
841	if (ch == 'd') {
842	clusterD[clustD] = i;
843	clustD++;
844	}
845	}
846	/* check clusters */
847	valid = TRUE;
848	if (numclust == 4) {
849	if (clustA == 0) {
850	valid = FALSE;
851	numclust = 1;
852	FPRINTF(STDOUTFILE "\n\n\nNo sequence in cluster a\n");
853	}
854	if (clustB == 0) {
855	valid = FALSE;
856	numclust = 1;
857	FPRINTF(STDOUTFILE "\n\n\nNo sequence in cluster b\n");
858	}
859	if (clustC == 0) {
860	valid = FALSE;
861	numclust = 1;
862	FPRINTF(STDOUTFILE "\n\n\nNo sequence in cluster c\n");
863	}
864	if (clustD == 0) {
865	valid = FALSE;
866	numclust = 1;
867	FPRINTF(STDOUTFILE "\n\n\nNo sequence in cluster d\n");
868	}
869	}
870	if (numclust == 3) {
871	if (clustA == 0) {
872	valid = FALSE;
873	numclust = 1;
874	FPRINTF(STDOUTFILE "\n\n\nNo sequence in cluster a\n");
875	}
876	if (clustB == 0) {
877	valid = FALSE;
878	numclust = 1;
879	FPRINTF(STDOUTFILE "\n\n\nNo sequence in cluster b\n");
880	}
881	if (clustC < 2) {
882	valid = FALSE;
883	numclust = 1;
884	if (clustC == 0)
885	FPRINTF(STDOUTFILE "\n\n\nNo sequence in cluster c\n");
886	else
887	FPRINTF(STDOUTFILE "\n\n\nOnly one sequence in cluster c\n");
888	}
889	}
890	if (numclust == 2) {
891	if (clustA < 2) {
892	valid = FALSE;
893	numclust = 1;
894	if (clustA == 0)
895	FPRINTF(STDOUTFILE "\n\n\nNo sequence in cluster a\n");
896	else
897	FPRINTF(STDOUTFILE "\n\n\nOnly one sequence in cluster a\n");
898	}
899	if (clustB < 2) {
900	valid = FALSE;
901	numclust = 1;
902	if (clustB == 0)
903	FPRINTF(STDOUTFILE "\n\n\nNo sequence in cluster b\n");
904	else
905	FPRINTF(STDOUTFILE "\n\n\nOnly one sequence in cluster b\n");
906	}
907	}
908	if (valid) {
909	FPRINTF(STDOUTFILE "\nNumber of sequences in each cluster:\n\n");
910	FPRINTF(STDOUTFILE "Cluster a: %d\n", clustA);
911	FPRINTF(STDOUTFILE "Cluster b: %d\n", clustB);
912	if (numclust > 2)
913	FPRINTF(STDOUTFILE "Cluster c: %d\n", clustC);
914	if (numclust == 4)
915	FPRINTF(STDOUTFILE "Cluster d: %d\n", clustD);
916	FPRINTF(STDOUTFILE "\nExcluded sequences: ");
917	if (numclust == 2) FPRINTF(STDOUTFILE "%d\n",
918	Maxspc-clustA-clustB);
919	if (numclust == 3) FPRINTF(STDOUTFILE "%d\n",
920	Maxspc-clustA-clustB-clustC);
921	if (numclust == 4) FPRINTF(STDOUTFILE "%d\n",
922	Maxspc-clustA-clustB-clustC-clustD);
923
924	}
925	}
926	}
927	/* number of resulting quartets */
928	compnumqts();
929
930	} else {
931	FPRINTF(STDOUTFILE "\n\n\nThis is not a possible option!\n");
932	}
933	break;
934
935	case 'd': if (auto_datatype == AUTO_GUESS) {
936	auto_datatype = AUTO_OFF;
937	guessdata_optn = data_optn;
938	data_optn = 0;
939	} else {
940	data_optn = data_optn + 1;
941	if (data_optn == 3) {
942	auto_datatype = AUTO_GUESS;
943	data_optn = guessdata_optn;
944	}
945	}
946	/* translate characters into format used by ML engine */
947	translatedataset();
948	estimatebasefreqs();
949	break;
950
951	case 'u': if (puzzlemode == QUARTPUZ && typ_optn == TREERECON_OPTN)
952	show_optn = 1 - show_optn;
953	else
954	FPRINTF(STDOUTFILE "\n\n\nThis is not a possible option!\n");
955	break;
956
957	case 'j': if (puzzlemode == QUARTPUZ && typ_optn == TREERECON_OPTN)
958	listqptrees = (listqptrees + 1) % 4;
959	else
960	FPRINTF(STDOUTFILE "\n\n\nThis is not a possible option!\n");
961	break;
962
963	case 'v': if (puzzlemode == QUARTPUZ && typ_optn == TREERECON_OPTN)
964	approxqp = 1 - approxqp;
965	else
966	FPRINTF(STDOUTFILE "\n\n\nThis is not a possible option!\n");
967	break;
968
969	case 'f': if (Frequ_optn) {
970	tstvf84 = 0.0;
971	Frequ_optn = FALSE;
972	sumfreq = 0.0;
973	if (data_optn == AMINOACID)
974	FPRINTF(STDOUTFILE "\n\n\nAmino acid");
975	else if (data_optn == NUCLEOTIDE && SH_optn)
976	FPRINTF(STDOUTFILE "\n\n\nDoublet");
977	else if (data_optn == NUCLEOTIDE && nuc_optn)
978	FPRINTF(STDOUTFILE "\n\n\nNucleotide");
979	else if (data_optn == BINARY)
980	FPRINTF(STDOUTFILE "\n\n\nBinary state");
981	FPRINTF(STDOUTFILE " frequencies (in %%):\n\n");
982	for (i = 0; i < gettpmradix() - 1; i++) {
983	FPRINTF(STDOUTFILE "pi(%s) = ", int2code(i));
984	scanf("%lf", &(Freqtpm[i]));
985	do ;
986	while (getchar() != '\n');
987	Freqtpm[i] = Freqtpm[i]/100.0;
988	if (Freqtpm[i] < 0.0) {
989	FPRINTF(STDOUTFILE "\n\n\nNegative frequency not possible\n");
990	estimatebasefreqs();
991	break;
992	}
993	sumfreq = sumfreq + Freqtpm[i];
994	if (sumfreq > 1.0) {
995	FPRINTF(STDOUTFILE "\n\n\nThe sum of ");
996	FPRINTF(STDOUTFILE "all frequencies exceeds");
997	FPRINTF(STDOUTFILE " 100%%\n");
998	estimatebasefreqs();
999	break;
1000	}
1001	if (i == gettpmradix() - 2)
1002	Freqtpm[i+1] = 1.0 - sumfreq;
1003	}
1004	} else estimatebasefreqs();
1005	break;
1006
1007	case 's': if (data_optn == NUCLEOTIDE && SH_optn) {
1008	sym_optn = 1 - sym_optn;
1009	} else {
1010	FPRINTF(STDOUTFILE "\n\n\nThis is not a possible option!\n");
1011	}
1012	break;
1013
1014	case 'n': if (puzzlemode == QUARTPUZ && typ_optn == TREERECON_OPTN)
1015	{
1016	FPRINTF(STDOUTFILE "\n\n\nNumber of puzzling steps: ");
1017	scanf("%lu", &Numtrial);
1018	do ;
1019	while (getchar() != '\n');
1020	if (Numtrial < 1) {
1021	FPRINTF(STDOUTFILE "\n\n\nThe number of puzzling");
1022	FPRINTF(STDOUTFILE " steps can't be smaller than one\n");
1023	Numtrial = 1000;
1024	}
1025	}
1026	else if (typ_optn == LIKMAPING_OPTN)
1027	{
1028	FPRINTF(STDOUTFILE "\n\nEnter zero to use all possible");
1029	FPRINTF(STDOUTFILE " quartets in the analysis!\n");
1030	FPRINTF(STDOUTFILE "\nNumber of random quartets: ");
1031	scanf("%lu", &lmqts);
1032	do ;
1033	while (getchar() != '\n');
1034
1035	/* compute number of quartets used */
1036	compnumqts();
1037	}
1038	else
1039	{
1040	FPRINTF(STDOUTFILE "\n\n\nThis is not a possible option!\n");
1041	}
1042	break;
1043
1044	case 'o': if (puzzlemode == QUARTPUZ && typ_optn == TREERECON_OPTN) {
1045	FPRINTF(STDOUTFILE "\n\n\nSequence to be displayed as outgroup (1-%d): ",
1046	Maxspc);
1047	scanf("%d", &outgroup);
1048	do ;
1049	while (getchar() != '\n');
1050	if (outgroup < 1 \|\| outgroup > Maxspc) {
1051	FPRINTF(STDOUTFILE "\n\n\nSequences are numbered ");
1052	FPRINTF(STDOUTFILE "from 1 to %d\n",
1053	Maxspc);
1054	outgroup = 1;
1055	}
1056	outgroup = outgroup - 1;
1057	} else {
1058	FPRINTF(STDOUTFILE "\n\n\nThis is not a possible option!\n");
1059	}
1060	break;
1061
1062	case 'm': if (data_optn == NUCLEOTIDE) { /* nucleotide data */
1063	if(HKY_optn && nuc_optn) {
1064	/* HKY -> TN */
1065	tstvf84 = 0.0;
1066	TSparam = 2.0;
1067	YRparam = 0.9;
1068	HKY_optn = FALSE;
1069	TN_optn = TRUE;
1070	optim_optn = TRUE;
1071	nuc_optn = TRUE;
1072	SH_optn = FALSE;
1073	break;
1074	}
1075	if(TN_optn && nuc_optn) {
1076	if (Maxseqc % 2 == 0 \|\| Maxseqc % 3 == 0) {
1077	/* number of chars needs to be a multiple 2 or 3 */
1078	/* TN -> SH */
1079	if (Maxseqc % 2 != 0 && Maxseqc % 3 == 0)
1080	SHcodon = TRUE;
1081	else
1082	SHcodon = FALSE;
1083	tstvf84 = 0.0;
1084	TSparam = 2.0;
1085	YRparam = 1.0;
1086	HKY_optn = TRUE;
1087	TN_optn = FALSE;
1088	optim_optn = TRUE;
1089	nuc_optn = FALSE;
1090	SH_optn = TRUE;
1091	/* translate characters into format */
1092	/* used by ML engine */
1093	translatedataset();
1094	estimatebasefreqs();
1095	} else {
1096	FPRINTF(STDOUTFILE "\n\n\nSH model not ");
1097	FPRINTF(STDOUTFILE "available for the data set!\n");
1098	/* TN -> HKY */
1099	tstvf84 = 0.0;
1100	TSparam = 2.0;
1101	YRparam = 1.0;
1102	HKY_optn = TRUE;
1103	TN_optn = FALSE;
1104	optim_optn = TRUE;
1105	nuc_optn = TRUE;
1106	SH_optn = FALSE;
1107	}
1108	break;
1109	}
1110	if(SH_optn) {
1111	/* SH -> HKY */
1112	tstvf84 = 0.0;
1113	TSparam = 2.0;
1114	YRparam = 1.0;
1115	HKY_optn = TRUE;
1116	TN_optn = FALSE;
1117	optim_optn = TRUE;
1118	nuc_optn = TRUE;
1119	SH_optn = FALSE;
1120	/* translate characters into format */
1121	/* used by ML engine */
1122	translatedataset();
1123	estimatebasefreqs();
1124	break;
1125	}
1126	break;
1127	}
1128	if (data_optn == AMINOACID) { /* amino acid data */
1129	if (auto_aamodel) {
1130	/* AUTO -> Dayhoff */
1131	Dayhf_optn = TRUE;
1132	Jtt_optn = FALSE;
1133	mtrev_optn = FALSE;
1134	cprev_optn = FALSE;
1135	blosum62_optn = FALSE;
1136	vtmv_optn = FALSE;
1137	wag_optn = FALSE;
1138	auto_aamodel = AUTO_OFF;
1139	break;
1140	}
1141	if (Dayhf_optn) {
1142	/* Dayhoff -> JTT */
1143	Dayhf_optn = FALSE;
1144	Jtt_optn = TRUE;
1145	mtrev_optn = FALSE;
1146	cprev_optn = FALSE;
1147	blosum62_optn = FALSE;
1148	vtmv_optn = FALSE;
1149	wag_optn = FALSE;
1150	auto_aamodel = AUTO_OFF;
1151	break;
1152	}
1153	if (Jtt_optn) {
1154	/* JTT -> mtREV */
1155	Dayhf_optn = FALSE;
1156	Jtt_optn = FALSE;
1157	mtrev_optn = TRUE;
1158	cprev_optn = FALSE;
1159	blosum62_optn = FALSE;
1160	vtmv_optn = FALSE;
1161	wag_optn = FALSE;
1162	auto_aamodel = AUTO_OFF;
1163	break;
1164	}
1165	#ifdef CPREV
1166	if (mtrev_optn) {
1167	/* mtREV -> cpREV */
1168	Dayhf_optn = FALSE;
1169	Jtt_optn = FALSE;
1170	mtrev_optn = FALSE;
1171	cprev_optn = TRUE;
1172	blosum62_optn = FALSE;
1173	vtmv_optn = FALSE;
1174	wag_optn = FALSE;
1175	auto_aamodel = AUTO_OFF;
1176	break;
1177	}
1178	#else /* ! CPREV */
1179	if (mtrev_optn) {
1180	/* mtREV -> BLOSUM 62 */
1181	Dayhf_optn = FALSE;
1182	Jtt_optn = FALSE;
1183	mtrev_optn = FALSE;
1184	cprev_optn = FALSE;
1185	blosum62_optn = TRUE;
1186	vtmv_optn = FALSE;
1187	wag_optn = FALSE;
1188	auto_aamodel = AUTO_OFF;
1189	break;
1190	}
1191	#endif /* ! CPREV */
1192
1193	#ifdef CPREV
1194	if (cprev_optn) {
1195	/* cpREV -> BLOSUM 62 */
1196	Dayhf_optn = FALSE;
1197	Jtt_optn = FALSE;
1198	mtrev_optn = FALSE;
1199	cprev_optn = FALSE;
1200	blosum62_optn = TRUE;
1201	vtmv_optn = FALSE;
1202	wag_optn = FALSE;
1203	auto_aamodel = AUTO_OFF;
1204	break;
1205	}
1206	#endif
1207	if (blosum62_optn) {
1208	/* BLOSUM 62 -> VT model */
1209	Dayhf_optn = FALSE;
1210	Jtt_optn = FALSE;
1211	mtrev_optn = FALSE;
1212	cprev_optn = FALSE;
1213	blosum62_optn = FALSE;
1214	vtmv_optn = TRUE;
1215	wag_optn = FALSE;
1216	auto_aamodel = AUTO_OFF;
1217	break;
1218	}
1219	if (vtmv_optn) {
1220	/* VT model -> WAG model */
1221	Dayhf_optn = FALSE;
1222	Jtt_optn = FALSE;
1223	mtrev_optn = FALSE;
1224	cprev_optn = FALSE;
1225	blosum62_optn = FALSE;
1226	vtmv_optn = FALSE;
1227	wag_optn = TRUE;
1228	auto_aamodel = AUTO_OFF;
1229	break;
1230	}
1231	if (wag_optn) {
1232	/* WAG model -> AUTO */
1233	Dayhf_optn = guessDayhf_optn;
1234	Jtt_optn = guessJtt_optn;
1235	mtrev_optn = guessmtrev_optn;
1236	cprev_optn = guesscprev_optn;
1237	blosum62_optn = guessblosum62_optn;
1238	vtmv_optn = guessvtmv_optn;
1239	wag_optn = guesswag_optn;
1240	auto_aamodel = guessauto_aamodel;
1241	break;
1242	}
1243	break;
1244	}
1245	if (data_optn == BINARY) {
1246	FPRINTF(STDOUTFILE "\n\n\nNo other model available!\n");
1247	}
1248	break;
1249
1250	case 't': if (data_optn != NUCLEOTIDE) {
1251	FPRINTF(STDOUTFILE "\n\n\nThis is not a possible option!\n");
1252	} else {
1253	tstvf84 = 0.0;
1254	FPRINTF(STDOUTFILE "\n\n\nEnter an invalid value for ");
1255	FPRINTF(STDOUTFILE "estimation from data set!\n");
1256	FPRINTF(STDOUTFILE "\nTransition/transversion parameter (%.2f-%.2f): ",
1257	MINTS, MAXTS);
1258	scanf("%lf", &TSparam);
1259	do ;
1260	while (getchar() != '\n');
1261	if (TSparam < MINTS \|\| TSparam > MAXTS) {
1262	optim_optn = TRUE;
1263	TSparam = 2.0;
1264	} else {
1265	optim_optn = FALSE;
1266	}
1267	}
1268	break;
1269
1270	case 'q': FPRINTF(STDOUTFILE "\n\n\n");
1271	# if PARALLEL
1272	PP_SendDone();
1273	MPI_Finalize();
1274	# endif /* PARALLEL */
1275	exit(0);
1276
1277	break;
1278
1279	case 'r': if (!(TN_optn && nuc_optn)){
1280	FPRINTF(STDOUTFILE "\n\n\nThis is not a possible option!\n");
1281	} else {
1282	tstvf84 = 0.0;
1283	FPRINTF(STDOUTFILE "\n\n\nEnter an invalid value ");
1284	FPRINTF(STDOUTFILE "for estimation from data set!\n");
1285	FPRINTF(STDOUTFILE "\nY/R transition parameter (%.2f-%.2f): ", MINYR, MAXYR);
1286	scanf("%lf", &YRparam);
1287	do ;
1288	while (getchar() != '\n');
1289	if (YRparam < MINYR \|\| YRparam > MAXYR) {
1290	optim_optn = TRUE;
1291	YRparam = 0.9;
1292	} else if (YRparam == 1.0) {
1293	TN_optn = FALSE;
1294	HKY_optn = TRUE;
1295	if (optim_optn) TSparam = 2.0;
1296	} else {
1297	optim_optn = FALSE;
1298	}
1299	}
1300	break;
1301
1302	case 'p': if (!(TN_optn && nuc_optn)){
1303	FPRINTF(STDOUTFILE "\n\n\nThis is not a possible option!\n");
1304	} else {
1305	FPRINTF(STDOUTFILE "\n\n\nThe F84 model (Felsenstein 1984) is a restricted");
1306	FPRINTF(STDOUTFILE " TN model, and the one\nF84 parameter uniquely");
1307	FPRINTF(STDOUTFILE " determines the two corresponding TN parameters!\n\n");
1308	FPRINTF(STDOUTFILE "F84 expected transition/transversion ratio: ");
1309	scanf("%lf", &tstvf84);
1310	do ;
1311	while (getchar() != '\n');
1312	if (tstvf84 <= 0.0) tstvf84 = 0.0;
1313	else makeF84model();
1314	}
1315	break;
1316
1317	case 'y': break;
1318
1319	default: FPRINTF(STDOUTFILE "\n\n\nThis is not a possible option!\n");
1320	break;
1321	}
1322	} while (ch != 'y');
1323
1324	FPRINTF(STDOUTFILE "\n\n\n");
1325	}
1326
1327	/* open file for reading */
1328	void openfiletoread(FILE **fp, char name[], char descr[])
1329	{
1330	int count = 0;
1331	cvector str;
1332
1333	if ((*fp = fopen(name, "r")) == NULL) {
1334	FPRINTF(STDOUTFILE "\n\n\nPlease enter a file name for the %s: ", descr);
1335	str = mygets();
1336	while ((*fp = fopen(str, "r")) == NULL)
1337	{
1338	count++;
1339	if (count > 10)
1340	{
1341	FPRINTF(STDOUTFILE "\n\n\nToo many trials - quitting ...\n");
1342	exit(1);
1343	}
1344	FPRINTF(STDOUTFILE "File '%s' not found, ", str);
1345	FPRINTF(STDOUTFILE "please enter alternative name: ");
1346	free_cvector(str);
1347	str = mygets();
1348	}
1349	free_cvector(str);
1350	FPRINTF(STDOUTFILE "\n");
1351	}
1352	} /* openfiletoread */
1353
1354
1355	/* open file for writing */
1356	void openfiletowrite(FILE **fp, char name[], char descr[])
1357	{
1358	int count = 0;
1359	cvector str;
1360
1361	if ((*fp = fopen(name, "w")) == NULL) {
1362	FPRINTF(STDOUTFILE "\n\n\nPlease enter a file name for the %s: ", descr);
1363	str = mygets();
1364	while ((*fp = fopen(str, "w")) == NULL)
1365	{
1366	count++;
1367	if (count > 10)
1368	{
1369	FPRINTF(STDOUTFILE "\n\n\nToo many trials - quitting ...\n");
1370	exit(1);
1371	}
1372	FPRINTF(STDOUTFILE "File '%s' not created, ", str);
1373	FPRINTF(STDOUTFILE "please enter other name: ");
1374	free_cvector(str);
1375	str = mygets();
1376	}
1377	free_cvector(str);
1378	FPRINTF(STDOUTFILE "\n");
1379	}
1380	} /* openfiletowrite */
1381
1382
1383	/* open file for appending */
1384	void openfiletoappend(FILE **fp, char name[], char descr[])
1385	{
1386	int count = 0;
1387	cvector str;
1388
1389	if ((*fp = fopen(name, "a")) == NULL) {
1390	FPRINTF(STDOUTFILE "\n\n\nPlease enter a file name for the %s: ", descr);
1391	str = mygets();
1392	while ((*fp = fopen(str, "a")) == NULL)
1393	{
1394	count++;
1395	if (count > 10)
1396	{
1397	FPRINTF(STDOUTFILE "\n\n\nToo many trials - quitting ...\n");
1398	exit(1);
1399	}
1400	FPRINTF(STDOUTFILE "File '%s' not created, ", str);
1401	FPRINTF(STDOUTFILE "please enter other name: ");
1402	free_cvector(str);
1403	str = mygets();
1404	}
1405	free_cvector(str);
1406	FPRINTF(STDOUTFILE "\n");
1407	}
1408	} /* openfiletowrite */
1409
1410
1411	/* close file */
1412	void closefile(FILE *fp)
1413	{
1414	fclose(fp);
1415	} /* closefile */
1416
1417	/* symmetrize doublet frequencies */
1418	void symdoublets()
1419	{
1420	int i, imean;
1421	double mean;
1422
1423	if (data_optn == NUCLEOTIDE && SH_optn && sym_optn) {
1424	/* ML frequencies */
1425	mean = (Freqtpm[1] + Freqtpm[4])/2.0; /* AC CA */
1426	Freqtpm[1] = mean;
1427	Freqtpm[4] = mean;
1428	mean = (Freqtpm[2] + Freqtpm[8])/2.0; /* AG GA */
1429	Freqtpm[2] = mean;
1430	Freqtpm[8] = mean;
1431	mean = (Freqtpm[3] + Freqtpm[12])/2.0; /* AT TA */
1432	Freqtpm[3] = mean;
1433	Freqtpm[12] = mean;
1434	mean = (Freqtpm[6] + Freqtpm[9])/2.0; /* CG GC */
1435	Freqtpm[6] = mean;
1436	Freqtpm[9] = mean;
1437	mean = (Freqtpm[7] + Freqtpm[13])/2.0; /* CT TC */
1438	Freqtpm[7] = mean;
1439	Freqtpm[13] = mean;
1440	mean = (Freqtpm[11] + Freqtpm[14])/2.0; /* GT TG */
1441	Freqtpm[11] = mean;
1442	Freqtpm[14] = mean;
1443
1444	/* base composition of each taxon */
1445	for (i = 0; i < Maxspc; i++) {
1446	imean = (Basecomp[i][1] + Basecomp[i][4])/2; /* AC CA */
1447	Basecomp[i][1] = imean;
1448	Basecomp[i][4] = imean;
1449	imean = (Basecomp[i][2] + Basecomp[i][8])/2; /* AG GA */
1450	Basecomp[i][2] = imean;
1451	Basecomp[i][8] = imean;
1452	imean = (Basecomp[i][3] + Basecomp[i][12])/2; /* AT TA */
1453	Basecomp[i][3] = imean;
1454	Basecomp[i][12] = imean;
1455	imean = (Basecomp[i][6] + Basecomp[i][9])/2; /* CG GC */
1456	Basecomp[i][6] = imean;
1457	Basecomp[i][9] = imean;
1458	imean = (Basecomp[i][7] + Basecomp[i][13])/2; /* CT TC */
1459	Basecomp[i][7] = imean;
1460	Basecomp[i][13] = imean;
1461	imean = (Basecomp[i][11] + Basecomp[i][14])/2; /* GT TG */
1462	Basecomp[i][11] = imean;
1463	Basecomp[i][14] = imean;
1464	}
1465	}
1466	}
1467
1468	/* show Ts/Tv ratio and Ts Y/R ratio */
1469	void computeexpectations()
1470	{
1471	double AlphaYBeta, AlphaRBeta, piR, piY, num, denom, pyr, pur;
1472
1473	if (nuc_optn == TRUE) { /* 4x4 nucs */
1474	piR = Freqtpm[0] + Freqtpm[2];
1475	piY = Freqtpm[1] + Freqtpm[3];
1476	AlphaRBeta = 4.0*TSparam / (1 + YRparam);
1477	AlphaYBeta = AlphaRBeta * YRparam;
1478	tstvratio = (AlphaRBetaFreqtpm[0]Freqtpm[2] +
1479	AlphaYBetaFreqtpm[1]Freqtpm[3])/(piR * piY);
1480	yrtsratio = (AlphaYBetaFreqtpm[1]Freqtpm[3]) /
1481	(AlphaRBetaFreqtpm[0]Freqtpm[2]);
1482	} else { /* 16x16 nucs */
1483	pyr = Freqtpm[1]Freqtpm[3] + Freqtpm[5]Freqtpm[7] +
1484	Freqtpm[9]Freqtpm[11] + Freqtpm[4]Freqtpm[12] +
1485	Freqtpm[5]Freqtpm[13] + Freqtpm[6]Freqtpm[14] +
1486	Freqtpm[7]Freqtpm[15] + Freqtpm[13]Freqtpm[15];
1487	pur = Freqtpm[0]Freqtpm[2] + Freqtpm[4]Freqtpm[6] +
1488	Freqtpm[0]Freqtpm[8] + Freqtpm[1]Freqtpm[9] +
1489	Freqtpm[2]Freqtpm[10] + Freqtpm[8]Freqtpm[10] +
1490	Freqtpm[3]Freqtpm[11] + Freqtpm[12]Freqtpm[14];
1491	num = pyr + pur;
1492	denom = Freqtpm[0]Freqtpm[1] + Freqtpm[1]Freqtpm[2] +
1493	Freqtpm[0]Freqtpm[3] + Freqtpm[2]Freqtpm[3] +
1494	Freqtpm[0]Freqtpm[4] + Freqtpm[1]Freqtpm[5] +
1495	Freqtpm[4]Freqtpm[5] + Freqtpm[2]Freqtpm[6] +
1496	Freqtpm[5]Freqtpm[6] + Freqtpm[3]Freqtpm[7] +
1497	Freqtpm[4]Freqtpm[7] + Freqtpm[6]Freqtpm[7] +
1498	Freqtpm[4]Freqtpm[8] + Freqtpm[5]Freqtpm[9] +
1499	Freqtpm[8]Freqtpm[9] + Freqtpm[6]Freqtpm[10] +
1500	Freqtpm[9]Freqtpm[10] + Freqtpm[7]Freqtpm[11] +
1501	Freqtpm[8]Freqtpm[11] + Freqtpm[10]Freqtpm[11] +
1502	Freqtpm[0]Freqtpm[12] + Freqtpm[8]Freqtpm[12] +
1503	Freqtpm[1]Freqtpm[13] + Freqtpm[9]Freqtpm[13] +
1504	Freqtpm[12]Freqtpm[13] + Freqtpm[2]Freqtpm[14] +
1505	Freqtpm[10]Freqtpm[14] + Freqtpm[13]Freqtpm[14] +
1506	Freqtpm[3]Freqtpm[15] + Freqtpm[11]Freqtpm[15] +
1507	Freqtpm[12]Freqtpm[15] + Freqtpm[14]Freqtpm[15];
1508	tstvratio = 2.0TSparam num/denom;
1509	yrtsratio = pyr/pur;
1510	}
1511	}
1512
1513	/* write ML distance matrix to file */
1514	void putdistance(FILE *fp)
1515	{
1516	int i, j;
1517
1518	fprintf(fp, " %d\n", Maxspc);
1519	for (i = 0; i < Maxspc; i++) {
1520	fputid10(fp, i);
1521	for (j = 0; j < Maxspc; j++) {
1522	fprintf(fp, " %.5f", Distanmat[i][j]/100.0);
1523	/* seven in one row */
1524	if ((j + 1) % 7 == 0 && j+1 != Maxspc)
1525	fprintf(fp, "\n ");
1526	}
1527	fprintf(fp, "\n");
1528	}
1529	}
1530
1531
1532	/* find identical sequences */
1533	void findidenticals(FILE *fp)
1534	{
1535	int i, j, noids;
1536	cvector useqs;
1537
1538	useqs = new_cvector(Maxspc);
1539
1540	for (i = 0; i < Maxspc; i++)
1541	useqs[i] = 0;
1542
1543	noids = TRUE;
1544	for (i = 0; i < Maxspc && noids; i++)
1545	for (j = i + 1; j < Maxspc && noids; j++)
1546	if (Distanmat[i][j] == 0.0) noids = FALSE;
1547
1548	if (noids)
1549	fprintf(fp, " All sequences are unique.\n");
1550	else {
1551	for (i = 0; i < Maxspc; i++) {
1552	noids = TRUE;
1553	for (j = i + 1; j < Maxspc && noids; j++)
1554	if (Distanmat[i][j] == 0.0) noids = FALSE;
1555
1556	if (!noids && useqs[i] == 0) {
1557	fputid(fp, i);
1558	useqs[i] = 1;
1559	for (j = i + 1; j < Maxspc; j++)
1560	if (Distanmat[i][j] == 0.0) {
1561	fprintf(fp, ", ");
1562	fputid(fp, j);
1563	useqs[j] = 1;
1564	}
1565	fprintf(fp, ".\n");
1566	}
1567	}
1568	}
1569	free_cvector(useqs);
1570	}
1571
1572	/* compute average distance */
1573	double averagedist()
1574	{
1575	int i, j;
1576	double sum;
1577
1578	sum = 0.0;
1579	for (i = 0; i < Maxspc; i++)
1580	for (j = i + 1; j < Maxspc; j++)
1581	sum = sum + Distanmat[i][j];
1582
1583	sum = sum / (double) Maxspc / ((double) Maxspc - 1.0) * 2.0;
1584
1585	return sum;
1586	}
1587
1588	/* first lines of EPSF likelihood mapping file */
1589	void initps(FILE *ofp)
1590	{
1591	fprintf(ofp, "%%!PS-Adobe-3.0 EPSF-3.0\n");
1592	fprintf(ofp, "%%%%BoundingBox: 60 210 550 650\n");
1593	fprintf(ofp, "%%%%Pages: 1\n");
1594	# ifndef ALPHA
1595	fprintf(ofp, "%%%%Creator: %s (version %s)\n", PACKAGE, VERSION);
1596	# else
1597	fprintf(ofp, "%%%%Creator: %s (version %s%s)\n", PACKAGE, VERSION, ALPHA);
1598	# endif
1599	fprintf(ofp, "%%%%Title: Likelihood Mapping Analysis\n");
1600	fprintf(ofp, "%%%%CreationDate: %s", asctime(localtime(&Starttime)) );
1601	fprintf(ofp, "%%%%DocumentFonts: Helvetica\n");
1602	fprintf(ofp, "%%%%DocumentNeededFonts: Helvetica\n");
1603	fprintf(ofp, "%%%%EndComments\n");
1604	fprintf(ofp, "%% use inch as unit\n");
1605	fprintf(ofp, "/inch {72 mul} def\n");
1606	fprintf(ofp, "%% triangle side length (3 inch)\n");
1607	fprintf(ofp, "/tl {3 inch mul} def\n");
1608	fprintf(ofp, "%% plot one dot (x-y coordinates on stack)\n");
1609	fprintf(ofp, "/dot {\n");
1610	fprintf(ofp, "newpath\n");
1611	fprintf(ofp, "0.002 tl 0 360 arc %% radius is 0.002 of the triangle length\n");
1612	fprintf(ofp, "closepath\n");
1613	fprintf(ofp, "fill\n");
1614	fprintf(ofp, "} def\n");
1615	fprintf(ofp, "%% preamble\n");
1616	fprintf(ofp, "/Helvetica findfont\n");
1617	fprintf(ofp, "12 scalefont\n");
1618	fprintf(ofp, "setfont\n");
1619	fprintf(ofp, "%% 0/0 for triangle of triangles\n");
1620	fprintf(ofp, "0.9 inch 3 inch translate\n");
1621	fprintf(ofp, "%% first triangle (the one with dots)\n");
1622	fprintf(ofp, "0.6 tl 1.2 tl 0.8660254038 mul translate\n");
1623	fprintf(ofp, "newpath\n");
1624	fprintf(ofp, " 0.0 tl 0.0 tl moveto\n");
1625	fprintf(ofp, " 1.0 tl 0.0 tl lineto\n");
1626	fprintf(ofp, " 0.5 tl 0.8660254038 tl lineto\n");
1627	fprintf(ofp, "closepath\n");
1628	fprintf(ofp, "stroke\n");
1629	}
1630
1631	/* plot one point of likelihood mapping analysis */
1632	void plotlmpoint(FILE *ofp, double w1, double w2)
1633	{
1634	fprintf(ofp,"%.10f tl %.10f tl dot\n",
1635	0.5w1 + w2, w10.8660254038);
1636	}
1637
1638	/* last lines of EPSF likelihood mapping file */
1639	void finishps(FILE *ofp)
1640	{
1641	fprintf(ofp, "stroke\n");
1642	fprintf(ofp, "%% second triangle (the one with 3 basins)\n");
1643	fprintf(ofp, "/secondtriangle {\n");
1644	fprintf(ofp, "newpath\n");
1645	fprintf(ofp, " 0.0 tl 0.0 tl moveto\n");
1646	fprintf(ofp, " 1.0 tl 0.0 tl lineto\n");
1647	fprintf(ofp, " 0.5 tl 0.8660254038 tl lineto\n");
1648	fprintf(ofp, "closepath\n");
1649	fprintf(ofp, "stroke\n");
1650	fprintf(ofp, "newpath\n");
1651	fprintf(ofp, " 0.50 tl 0.2886751346 tl moveto\n");
1652	fprintf(ofp, " 0.50 tl 0.0000000000 tl lineto\n");
1653	fprintf(ofp, "stroke\n");
1654	fprintf(ofp, "newpath\n");
1655	fprintf(ofp, " 0.50 tl 0.2886751346 tl moveto\n");
1656	fprintf(ofp, " 0.25 tl 0.4330127019 tl lineto\n");
1657	fprintf(ofp, "stroke\n");
1658	fprintf(ofp, "newpath\n");
1659	fprintf(ofp, " 0.50 tl 0.2886751346 tl moveto\n");
1660	fprintf(ofp, " 0.75 tl 0.4330127019 tl lineto\n");
1661	fprintf(ofp, "stroke\n");
1662	fprintf(ofp, "0.44 tl 0.5 tl moveto %% up\n");
1663	fprintf(ofp, "(%.1f%%) show\n", (double) ar1*100.0/Numquartets);
1664	fprintf(ofp, "0.25 tl 0.15 tl moveto %% down left\n");
1665	fprintf(ofp, "(%.1f%%) show\n", (double) ar3*100.0/Numquartets);
1666	fprintf(ofp, "0.63 tl 0.15 tl moveto %% down right\n");
1667	fprintf(ofp, "(%.1f%%) show\n", (double) ar2*100.0/Numquartets);
1668	fprintf(ofp, "} def\n");
1669	fprintf(ofp, "%% third triangle (the one with 7 basins)\n");
1670	fprintf(ofp, "/thirdtriangle {\n");
1671	fprintf(ofp, "newpath\n");
1672	fprintf(ofp, " 0.0 tl 0.0 tl moveto\n");
1673	fprintf(ofp, " 1.0 tl 0.0 tl lineto\n");
1674	fprintf(ofp, " 0.5 tl 0.8660254038 tl lineto\n");
1675	fprintf(ofp, "closepath\n");
1676	fprintf(ofp, "stroke\n");
1677	fprintf(ofp, "newpath\n");
1678	fprintf(ofp, " 0.25 tl 0.1443375673 tl moveto\n");
1679	fprintf(ofp, " 0.75 tl 0.1443375673 tl lineto\n");
1680	fprintf(ofp, " 0.50 tl 0.5773502692 tl lineto\n");
1681	fprintf(ofp, "closepath\n");
1682	fprintf(ofp, "stroke\n");
1683	fprintf(ofp, "newpath\n");
1684	fprintf(ofp, " 0.125 tl 0.2165063509 tl moveto\n");
1685	fprintf(ofp, " 0.250 tl 0.1443375673 tl lineto\n");
1686	fprintf(ofp, "stroke\n");
1687	fprintf(ofp, "newpath\n");
1688	fprintf(ofp, " 0.375 tl 0.6495190528 tl moveto\n");
1689	fprintf(ofp, " 0.500 tl 0.5773502692 tl lineto\n");
1690	fprintf(ofp, "stroke\n");
1691	fprintf(ofp, "newpath\n");
1692	fprintf(ofp, " 0.625 tl 0.6495190528 tl moveto\n");
1693	fprintf(ofp, " 0.500 tl 0.5773502692 tl lineto\n");
1694	fprintf(ofp, "stroke\n");
1695	fprintf(ofp, "newpath\n");
1696	fprintf(ofp, " 0.875 tl 0.2165063509 tl moveto\n");
1697	fprintf(ofp, " 0.750 tl 0.1443375673 tl lineto\n");
1698	fprintf(ofp, "stroke\n");
1699	fprintf(ofp, "newpath\n");
1700	fprintf(ofp, " 0.750 tl 0.00 tl moveto\n");
1701	fprintf(ofp, " 0.750 tl 0.1443375673 tl lineto\n");
1702	fprintf(ofp, "stroke\n");
1703	fprintf(ofp, "newpath\n");
1704	fprintf(ofp, " 0.250 tl 0.00 tl moveto\n");
1705	fprintf(ofp, " 0.250 tl 0.1443375673 tl lineto\n");
1706	fprintf(ofp, "stroke\n");
1707	fprintf(ofp, "0.42 tl 0.66 tl moveto %% up\n");
1708	fprintf(ofp, "(%.1f%%) show\n", (double) reg1*100.0/Numquartets);
1709	fprintf(ofp, "0.07 tl 0.05 tl moveto %% down left\n");
1710	fprintf(ofp, "(%.1f%%) show\n", (double) reg3*100.0/Numquartets);
1711	fprintf(ofp, "0.77 tl 0.05 tl moveto %% down right\n");
1712	fprintf(ofp, "(%.1f%%) show\n", (double) reg2*100.0/Numquartets);
1713	fprintf(ofp, "0.43 tl 0.05 tl moveto %% down side\n");
1714	fprintf(ofp, "(%.1f%%) show\n", (double) reg5*100.0/Numquartets);
1715	fprintf(ofp, "0.43 tl 0.28 tl moveto %% center\n");
1716	fprintf(ofp, "(%.1f%%) show\n", (double) reg7*100.0/Numquartets);
1717	fprintf(ofp, "gsave\n");
1718	fprintf(ofp, "-60 rotate\n");
1719	fprintf(ofp, "-0.07 tl 0.77 tl moveto %% right side\n");
1720	fprintf(ofp, "(%.1f%%) show\n", (double) reg4*100.0/Numquartets);
1721	fprintf(ofp, "grestore\n");
1722	fprintf(ofp, "gsave\n");
1723	fprintf(ofp, "60 rotate\n");
1724	fprintf(ofp, "0.4 tl -0.09 tl moveto %% left side\n");
1725	fprintf(ofp, "(%.1f%%) show\n", (double) reg6*100.0/Numquartets);
1726	fprintf(ofp, "grestore\n");
1727	fprintf(ofp, "} def\n");
1728	fprintf(ofp, "%% print the other two triangles\n");
1729	fprintf(ofp, "-0.6 tl -1.2 tl 0.8660254038 mul translate\n");
1730	fprintf(ofp, "secondtriangle\n");
1731	fprintf(ofp, "1.2 tl 0 translate\n");
1732	fprintf(ofp, "thirdtriangle\n");
1733	if (numclust == 4) { /* four cluster analysis */
1734	fprintf(ofp, "%% label corners\n");
1735	fprintf(ofp, "0.375 tl 0.9 tl moveto\n");
1736	fprintf(ofp, "((a,b)-(c,d)) show %% CHANGE HERE IF NECESSARY\n");
1737	fprintf(ofp, "-0.16 tl -0.08 tl moveto\n");
1738	fprintf(ofp, "((a,d)-(b,c)) show %% CHANGE HERE IF NECESSARY\n");
1739	fprintf(ofp, "0.92 tl -0.08 tl moveto\n");
1740	fprintf(ofp, "((a,c)-(b,d)) show %% CHANGE HERE IF NECESSARY\n");
1741	}
1742	if (numclust == 3) { /* three cluster analysis */
1743	fprintf(ofp, "%% label corners\n");
1744	fprintf(ofp, "0.375 tl 0.9 tl moveto\n");
1745	fprintf(ofp, "((a,b)-(c,c)) show %% CHANGE HERE IF NECESSARY\n");
1746	fprintf(ofp, "-0.16 tl -0.08 tl moveto\n");
1747	fprintf(ofp, "((a,c)-(b,c)) show %% CHANGE HERE IF NECESSARY\n");
1748	fprintf(ofp, "0.92 tl -0.08 tl moveto\n");
1749	fprintf(ofp, "((a,c)-(b,c)) show %% CHANGE HERE IF NECESSARY\n");
1750	}
1751	if (numclust == 2) { /* two cluster analysis */
1752	fprintf(ofp, "%% label corners\n");
1753	fprintf(ofp, "0.375 tl 0.9 tl moveto\n");
1754	fprintf(ofp, "((a,a)-(b,b)) show %% CHANGE HERE IF NECESSARY\n");
1755	fprintf(ofp, "-0.16 tl -0.08 tl moveto\n");
1756	fprintf(ofp, "((a,b)-(a,b)) show %% CHANGE HERE IF NECESSARY\n");
1757	fprintf(ofp, "0.92 tl -0.08 tl moveto\n");
1758	fprintf(ofp, "((a,b)-(a,b)) show %% CHANGE HERE IF NECESSARY\n");
1759	}
1760	fprintf(ofp, "showpage\n");
1761	fprintf(ofp, "%%%%EOF\n");
1762	}
1763
1764	/* computes LM point from the three log-likelihood values,
1765	plots the point, and does some statistics */
1766	void makelmpoint(FILE *fp, double b1, double b2, double b3)
1767	{
1768	double w1, w2, w3, temp;
1769	unsigned char qpbranching;
1770	double temp1, temp2, temp3, onethird;
1771	unsigned char discreteweight[3], treebits[3];
1772
1773	onethird = 1.0/3.0;
1774	treebits[0] = (unsigned char) 1;
1775	treebits[1] = (unsigned char) 2;
1776	treebits[2] = (unsigned char) 4;
1777
1778	/* sort in descending order */
1779	qweight[0] = b1;
1780	qweight[1] = b2;
1781	qweight[2] = b3;
1782	sort3doubles(qweight, qworder);
1783
1784	/* compute Bayesian weights */
1785	qweight[qworder[1]] = exp(qweight[qworder[1]]-qweight[qworder[0]]);
1786	qweight[qworder[2]] = exp(qweight[qworder[2]]-qweight[qworder[0]]);
1787	qweight[qworder[0]] = 1.0;
1788	temp = qweight[0] + qweight[1] + qweight[2];
1789	qweight[0] = qweight[0]/temp;
1790	qweight[1] = qweight[1]/temp;
1791	qweight[2] = qweight[2]/temp;
1792
1793	/* plot one point in likelihood mapping triangle */
1794	w1 = qweight[0];
1795	w2 = qweight[1];
1796	w3 = qweight[2];
1797	plotlmpoint(fp, w1, w2);
1798
1799	/* check areas 1,2,3 */
1800	if (treebits[qworder[0]] == 1) ar1++;
1801	else if (treebits[qworder[0]] == 2) ar2++;
1802	else ar3++;
1803
1804	/* check out regions 1,2,3,4,5,6,7 */
1805
1806	/* 100 distribution */
1807	temp1 = 1.0 - qweight[qworder[0]];
1808	sqdiff[0] = temp1*temp1 +
1809	qweight[qworder[1]]*qweight[qworder[1]] +
1810	qweight[qworder[2]]*qweight[qworder[2]];
1811	discreteweight[0] = treebits[qworder[0]];
1812
1813	/* 110 distribution */
1814	temp1 = 0.5 - qweight[qworder[0]];
1815	temp2 = 0.5 - qweight[qworder[1]];
1816	sqdiff[1] = temp1temp1 + temp2temp2 +
1817	qweight[qworder[2]]*qweight[qworder[2]];
1818	discreteweight[1] = treebits[qworder[0]] + treebits[qworder[1]];
1819
1820	/* 111 distribution */
1821	temp1 = onethird - qweight[qworder[0]];
1822	temp2 = onethird - qweight[qworder[1]];
1823	temp3 = onethird - qweight[qworder[2]];
1824	sqdiff[2] = temp1 * temp1 + temp2 * temp2 + temp3 * temp3;
1825	discreteweight[2] = (unsigned char) 7;
1826
1827	/* sort in descending order */
1828	sort3doubles(sqdiff, sqorder);
1829
1830	qpbranching = (unsigned char) discreteweight[sqorder[2]];
1831
1832	if (qpbranching == 1) {
1833	reg1++;
1834	if (w2 < w3) reg1l++;
1835	else reg1r++;
1836	}
1837	if (qpbranching == 2) {
1838	reg2++;
1839	if (w1 < w3) reg2d++;
1840	else reg2u++;
1841	}
1842	if (qpbranching == 4) {
1843	reg3++;
1844	if (w1 < w2) reg3d++;
1845	else reg3u++;
1846	}
1847	if (qpbranching == 3) {
1848	reg4++;
1849	if (w1 < w2) reg4d++;
1850	else reg4u++;
1851	}
1852	if (qpbranching == 6) {
1853	reg5++;
1854	if (w2 < w3) reg5l++;
1855	else reg5r++;
1856	}
1857	if (qpbranching == 5) {
1858	reg6++;
1859	if (w1 < w3) reg6d++;
1860	else reg6u++;
1861	}
1862	if (qpbranching == 7) reg7++;
1863	}
1864
1865	/* print tree statistics */
1866	void printtreestats(FILE *ofp)
1867	{
1868	int i, j, besttree;
1869	double bestlkl, difflkl, difflklps, temp, sum;
1870
1871	/* find best tree */
1872	besttree = 0;
1873	bestlkl = ulkl[0];
1874	for (i = 1; i < numutrees; i++)
1875	if (ulkl[i] > bestlkl) {
1876	besttree = i;
1877	bestlkl = ulkl[i];
1878	}
1879
1880	fprintf(ofp, "\n\nCOMPARISON OF USER TREES (NO CLOCK)\n\n");
1881	fprintf(ofp, "Tree log L difference S.E. Significantly worse\n");
1882	fprintf(ofp, "--------------------------------------------------------\n");
1883	for (i = 0; i < numutrees; i++) {
1884	difflkl = ulkl[besttree]-ulkl[i];
1885	fprintf(ofp, "%2d %10.2f %8.2f ", i+1, ulkl[i], difflkl);
1886	if (i == besttree) {
1887	fprintf(ofp, " <----------------- best tree");
1888	} else {
1889	/* compute variance of Log L differences over sites */
1890	difflklps = difflkl/(double)Maxsite;
1891	sum = 0.0;
1892	for (j = 0; j < Numptrn; j++) {
1893	temp = allsites[besttree][j] - allsites[i][j] - difflklps;
1894	sum += temptempWeight[j];
1895	}
1896	sum = sqrt(fabs(sum/(Maxsite-1.0)*Maxsite));
1897	fprintf(ofp, "%11.2f ", sum);
1898	if (difflkl > 1.96*sum)
1899	fprintf(ofp, "yes");
1900	else
1901	fprintf(ofp, "no");
1902	}
1903	fprintf(ofp, "\n");
1904	}
1905	fprintf(ofp, "\nThis test (5%% significance) follows Kishino and Hasegawa (1989).\n");
1906
1907	if (compclock) {
1908
1909	/* find best tree */
1910	besttree = 0;
1911	bestlkl = ulklc[0];
1912	for (i = 1; i < numutrees; i++)
1913	if (ulklc[i] > bestlkl) {
1914	besttree = i;
1915	bestlkl = ulklc[i];
1916	}
1917
1918	fprintf(ofp, "\n\nCOMPARISON OF USER TREES (WITH CLOCK)\n\n");
1919	fprintf(ofp, "Tree log L difference S.E. Significantly worse\n");
1920	fprintf(ofp, "--------------------------------------------------------\n");
1921	for (i = 0; i < numutrees; i++) {
1922	difflkl = ulklc[besttree]-ulklc[i];
1923	fprintf(ofp, "%2d %10.2f %8.2f ", i+1, ulklc[i], difflkl);
1924	if (i == besttree) {
1925	fprintf(ofp, " <----------------- best tree");
1926	} else {
1927	/* compute variance of Log L differences over sites */
1928	difflklps = difflkl/(double)Maxsite;
1929	sum = 0.0;
1930	for (j = 0; j < Numptrn; j++) {
1931	temp = allsitesc[besttree][j] - allsitesc[i][j] - difflklps;
1932	sum += temptempWeight[j];
1933	}
1934	sum = sqrt(fabs(sum/(Maxsite-1.0)*Maxsite));
1935	fprintf(ofp, "%11.2f ", sum);
1936	if (difflkl > 1.96*sum)
1937	fprintf(ofp, "yes");
1938	else
1939	fprintf(ofp, "no");
1940	}
1941	fprintf(ofp, "\n");
1942	}
1943	fprintf(ofp, "\nThis test (5%% significance) follows Kishino and Hasegawa (1989).\n");
1944	}
1945	}
1946
1947	/* time stamp */
1948	void timestamp(FILE* ofp)
1949	{
1950	double timespan;
1951	double cpuspan;
1952	timespan = difftime(Stoptime, Starttime);
1953	cpuspan = ((double) (Stopcpu - Startcpu) / CLOCKS_PER_SEC);
1954	fprintf(ofp, "\n\nTIME STAMP\n\n");
1955	fprintf(ofp, "Date and time: %s", asctime(localtime(&Starttime)) );
1956	fprintf(ofp, "Runtime (excl. input) : %.0f seconds (= %.1f minutes = %.1f hours)\n",
1957	timespan, timespan/60., timespan/3600.);
1958	fprintf(ofp, "Runtime (incl. input) : %.0f seconds (= %.1f minutes = %.1f hours)\n",
1959	fulltime, fulltime/60., fulltime/3600.);
1960	#ifdef TIMEDEBUG
1961	fprintf(ofp, "CPU time (incl. input): %.0f seconds (= %.1f minutes = %.1f hours)\n\n",
1962	fullcpu, fullcpu/60., fullcpu/3600.);
1963	#endif /* TIMEDEBUG */
1964
1965	}
1966
1967	/* extern int bestrfound; */
1968
1969	/* write output file */
1970	void writeoutputfile(FILE *ofp, int part)
1971	{
1972	int i, fail, df;
1973	uli li;
1974	double pval, delta;
1975
1976	if ((part == WRITEPARAMS) \|\| (part == WRITEALL)) {
1977	# ifndef ALPHA
1978	fprintf(ofp, "TREE-PUZZLE %s\n\n", VERSION);
1979	# else
1980	fprintf(ofp, "TREE-PUZZLE %s%s\n\n", VERSION, ALPHA);
1981	# endif
1982
1983	fprintf(ofp, "Input file name: %s\n",INFILE);
1984	if (puzzlemode == USERTREE) fprintf(ofp, "User tree file name: %s\n",INTREE);
1985
1986
1987	fprintf(ofp, "Type of analysis: ");
1988	if (typ_optn == TREERECON_OPTN) fprintf(ofp, "tree reconstruction\n");
1989	if (typ_optn == LIKMAPING_OPTN) fprintf(ofp, "likelihood mapping\n");
1990	fprintf(ofp, "Parameter estimation: ");
1991	if (approxp_optn) fprintf(ofp, "approximate (faster)\n");
1992	else fprintf(ofp, "accurate (slow)\n");
1993	if (!(puzzlemode == USERTREE && typ_optn == TREERECON_OPTN)) {
1994	fprintf(ofp, "Parameter estimation uses: ");
1995	if (qcalg_optn)
1996	fprintf(ofp, "quartet sampling (for substitution process) + NJ tree (for rate variation)\n");
1997	else
1998	fprintf(ofp, "neighbor-joining tree (for substitution process and rate variation)\n");
1999	} else {
2000	fprintf(ofp, "Parameter estimation uses: ");
2001	if (utree_optn)
2002	fprintf(ofp, "1st user tree (for substitution process and rate variation)\n");
2003	else if (qcalg_optn)
2004	fprintf(ofp, "quartet sampling (for substitution process) + NJ tree (for rate variation)\n");
2005	else
2006	fprintf(ofp, "neighbor-joining tree (for substitution process and rate variation)\n");
2007	}
2008	fprintf(ofp, "\nStandard errors (S.E.) are obtained by the curvature method.\n");
2009	fprintf(ofp, "The upper and lower bounds of an approximate 95%% confidence interval\n");
2010	fprintf(ofp, "for parameter or branch length x are x-1.96S.E. and x+1.96S.E.\n");
2011	fprintf(ofp, "\n\n");
2012	fprintf(ofp, "SEQUENCE ALIGNMENT\n\n");
2013	fprintf(ofp, "Input data: %d sequences with %d ", Maxspc, Maxsite);
2014	if (data_optn == AMINOACID)
2015	fprintf(ofp, "amino acid");
2016	else if (data_optn == NUCLEOTIDE && SH_optn)
2017	fprintf(ofp, "doublet (%d nucleotide)", Maxsite*2);
2018	else if (data_optn == NUCLEOTIDE && nuc_optn)
2019	fprintf(ofp, "nucleotide");
2020	else if (data_optn == BINARY)
2021	fprintf(ofp, "binary state");
2022	fprintf(ofp, " sites");
2023	if (data_optn == NUCLEOTIDE && (Maxseqc % 3) == 0 && !SH_optn) {
2024	if (codon_optn == 1) fprintf(ofp, " (1st codon positions)");
2025	if (codon_optn == 2) fprintf(ofp, " (2nd codon positions)");
2026	if (codon_optn == 3) fprintf(ofp, " (3rd codon positions)");
2027	if (codon_optn == 4) fprintf(ofp, " (1st and 2nd codon positions)");
2028	}
2029	if (data_optn == NUCLEOTIDE && SH_optn) {
2030	if (SHcodon)
2031	fprintf(ofp, " (1st and 2nd codon positions)");
2032	else
2033	fprintf(ofp, " (1st+2nd, 3rd+4th, etc. site)");
2034	}
2035	fprintf(ofp, "\n");
2036	fprintf(ofp, "Number of constant sites: %d (= %.1f%% of all sites)\n",
2037	Numconst, 100.0*fracconst);
2038	fprintf(ofp, "Number of site patterns: %d\n",
2039	Numptrn);
2040	fprintf(ofp, "Number of constant site patterns: %d (= %.1f%% of all site patterns)\n\n\n",
2041	Numconstpat, 100.0*fracconstpat);
2042	fprintf(ofp, "SUBSTITUTION PROCESS\n\n");
2043	fprintf(ofp, "Model of substitution: ");
2044	if (data_optn == NUCLEOTIDE) { /* nucleotides */
2045	if (nuc_optn) {
2046	if(HKY_optn) fprintf(ofp, "HKY (Hasegawa et al. 1985)\n");
2047	else fprintf(ofp, "TN (Tamura-Nei 1993)\n");
2048	fprintf(ofp, "Transition/transversion parameter");
2049	if (optim_optn)
2050	fprintf(ofp, " (estimated from data set)");
2051	fprintf(ofp, ": %.2f", TSparam);
2052	if (optim_optn)
2053	fprintf(ofp, " (S.E. %.2f)", tserr);
2054	fprintf(ofp, "\n");
2055
2056	if (optim_optn && TSparam > MAXTS - 1.0)
2057	fprintf(ofp, "WARNING --- parameter estimate close to internal upper bound!\n");
2058	if (optim_optn && TSparam < MINTS + 0.1)
2059	fprintf(ofp, "WARNING --- parameter estimate close to internal lower bound!\n");
2060
2061	if (TN_optn) {
2062	fprintf(ofp, "Y/R transition parameter");
2063	if (optim_optn)
2064	fprintf(ofp, " (estimated from data set)");
2065	fprintf(ofp, ": %.2f", YRparam);
2066	if (optim_optn)
2067	fprintf(ofp, " (S.E. %.2f)", yrerr);
2068	fprintf(ofp, "\n");
2069
2070	if (optim_optn && YRparam > MAXYR - 0.5)
2071	fprintf(ofp, "WARNING --- parameter estimate close to internal upper bound!\n");
2072	if (optim_optn && YRparam < MINYR + 0.1)
2073	fprintf(ofp, "WARNING --- parameter estimate close to internal lower bound!\n");
2074
2075	}
2076	}
2077	if (SH_optn) {
2078	fprintf(ofp, "SH (Schoeniger-von Haeseler 1994)\n");
2079	fprintf(ofp, "Transition/transversion parameter");
2080	if (optim_optn) fprintf(ofp, " (estimated from data set)");
2081	fprintf(ofp, ": %.2f\n", TSparam);
2082	if (optim_optn)
2083	fprintf(ofp, " (S.E. %.2f)", tserr);
2084	fprintf(ofp, "\n");
2085
2086	if (optim_optn && TSparam > MAXTS - 1.0)
2087	fprintf(ofp, "WARNING --- parameter estimate close to internal upper bound!\n");
2088	if (optim_optn && TSparam < MINTS + 0.1)
2089	fprintf(ofp, "WARNING --- parameter estimate close to internal lower bound!\n");
2090
2091	}
2092	}
2093	if (data_optn == AMINOACID) { /* amino acids */
2094	if (Dayhf_optn) fprintf(ofp, "Dayhoff (Dayhoff et al. 1978)\n");
2095	if (Jtt_optn) fprintf(ofp, "JTT (Jones et al. 1992)\n");
2096	if (mtrev_optn) fprintf(ofp, "mtREV24 (Adachi-Hasegawa 1996)\n");
2097	if (cprev_optn) fprintf(ofp, "cpREV45 (Adachi et al. 2000)\n");
2098	if (blosum62_optn) fprintf(ofp, "BLOSUM 62 (Henikoff-Henikoff 1992)\n");
2099	if (vtmv_optn) fprintf(ofp, "VT (Mueller-Vingron 2000)\n");
2100	if (wag_optn) fprintf(ofp, "WAG (Whelan-Goldman 2000)\n");
2101	}
2102	if (data_optn == BINARY) { /* binary states */
2103	fprintf(ofp, "Two-state model (Felsenstein 1981)\n");
2104	}
2105	if (data_optn == AMINOACID)
2106	fprintf(ofp, "Amino acid ");
2107	else if (data_optn == NUCLEOTIDE && SH_optn)
2108	fprintf(ofp, "Doublet ");
2109	else if (data_optn == NUCLEOTIDE && nuc_optn)
2110	fprintf(ofp, "Nucleotide ");
2111	else if (data_optn == BINARY)
2112	fprintf(ofp, "Binary state ");
2113	fprintf(ofp, "frequencies (");
2114	if (Frequ_optn) fprintf(ofp, "estimated from data set");
2115	else fprintf(ofp, "user specified");
2116	if (data_optn == NUCLEOTIDE && SH_optn && sym_optn)
2117	fprintf(ofp, " and symmetrized");
2118	fprintf(ofp, "):\n\n");
2119	for (i = 0; i < gettpmradix(); i++)
2120	fprintf(ofp, " pi(%s) = %5.1f%%\n",
2121	int2code(i), Freqtpm[i]*100);
2122	if (data_optn == NUCLEOTIDE) {
2123	fprintf(ofp, "\nExpected transition/transversion ratio: %.2f",
2124	tstvratio);
2125	if (tstvf84 == 0.0) fprintf(ofp, "\n");
2126	else fprintf(ofp, " (= F84 parameter)\n");
2127	fprintf(ofp, "Expected pyrimidine transition/purine transition");
2128	fprintf(ofp, " ratio: %.2f\n", yrtsratio);
2129	if (tstvf84 != 0.0 && TN_optn)
2130	fprintf(ofp,
2131	"This TN model is equivalent to a F84 model (Felsenstein 1984).\n");
2132	}
2133	fprintf(ofp, "\n\nRATE HETEROGENEITY\n\n");
2134	fprintf(ofp, "Model of rate heterogeneity: ");
2135	if (rhetmode == UNIFORMRATE) fprintf(ofp, "uniform rate\n");
2136	if (rhetmode == GAMMARATE ) fprintf(ofp, "Gamma distributed rates\n");
2137	if (rhetmode == TWORATE ) fprintf(ofp, "two rates (1 invariable + 1 variable)\n");
2138	if (rhetmode == MIXEDRATE ) fprintf(ofp, "mixed (1 invariable + %d Gamma rates)\n", numcats);
2139	if (rhetmode == TWORATE \|\| rhetmode == MIXEDRATE) {
2140	fprintf(ofp, "Fraction of invariable sites");
2141	if (fracinv_optim) fprintf(ofp, " (estimated from data set)");
2142	fprintf(ofp, ": %.2f", fracinv);
2143	if (fracinv_optim) fprintf(ofp, " (S.E. %.2f)", fierr);
2144	fprintf(ofp, "\n");
2145
2146	if (fracinv_optim && fracinv > MAXFI - 0.05)
2147	fprintf(ofp, "WARNING --- parameter estimate close to internal upper bound!\n");
2148
2149	fprintf(ofp, "Number of invariable sites: %.0f\n", floor(fracinv*Maxsite));
2150	}
2151	if (rhetmode == GAMMARATE \|\| rhetmode == MIXEDRATE) {
2152	fprintf(ofp, "Gamma distribution parameter alpha");
2153	if (grate_optim) fprintf(ofp, " (estimated from data set)");
2154	fprintf(ofp, ": %.2f", (1.0-Geta)/Geta);
2155	if (grate_optim) fprintf(ofp, " (S.E. %.2f)",
2156	geerr/(GetaGeta)); / first order approximation */
2157	fprintf(ofp, "\n");
2158
2159	if (grate_optim && Geta > MAXGE - 0.02)
2160	fprintf(ofp, "WARNING --- parameter estimate close to internal upper bound!\n");
2161	if (grate_optim && Geta < MINGE + 0.01)
2162	fprintf(ofp, "WARNING --- parameter estimate close to internal lower bound!\n");
2163
2164	fprintf(ofp, "Number of Gamma rate categories: %d\n", numcats);
2165	}
2166	if (rhetmode == MIXEDRATE) {
2167	fprintf(ofp, "Total rate heterogeneity (invariable sites + Gamma model): ");
2168	fprintf(ofp, "%.2f", fracinv + Geta - fracinv*Geta);
2169	if (grate_optim && fracinv_optim)
2170	fprintf(ofp, " (S.E. %.2f)", geerr + fierr); /* first order approximation */
2171	else if (grate_optim && !fracinv_optim)
2172	fprintf(ofp, " (S.E. %.2f)", geerr);
2173	else if (!grate_optim && fracinv_optim)
2174	fprintf(ofp, " (S.E. %.2f)", fierr);
2175	fprintf(ofp, "\n");
2176	}
2177	if (rhetmode != UNIFORMRATE) {
2178	fprintf(ofp, "\nRates and their respective probabilities used in the likelihood function:\n");
2179	fprintf(ofp, "\n Category Relative rate Probability\n");
2180	if (rhetmode == TWORATE \|\| rhetmode == MIXEDRATE)
2181	fprintf(ofp, " 0 0.0000 %.4f\n", fracinv);
2182	for (i = 0; i < numcats; i++)
2183	fprintf(ofp, " %d %.4f %.4f\n",
2184	i+1, Rates[i], (1.0-fracinv)/(double) numcats);
2185	}
2186	if (rhetmode == GAMMARATE \|\| rhetmode == MIXEDRATE) {
2187	fprintf(ofp, "\nCategories 1-%d approximate a continuous ", numcats);
2188	fprintf(ofp, "Gamma-distribution with expectation 1\n");
2189	fprintf(ofp, "and variance ");
2190	if (Geta == 1.0) fprintf(ofp, "infinity");
2191	else fprintf(ofp, "%.2f", Geta/(1.0-Geta));
2192	fprintf(ofp, ".\n");
2193	}
2194
2195	if (typ_optn == TREERECON_OPTN && (puzzlemode == QUARTPUZ \|\| puzzlemode == USERTREE))
2196	if (rhetmode != UNIFORMRATE) {
2197	fprintf(ofp, "\nCombination of categories that contributes");
2198	fprintf(ofp, " the most to the likelihood\n");
2199	fprintf(ofp, "(computation done without clock assumption assuming ");
2200	if (puzzlemode == QUARTPUZ) fprintf(ofp, "quartet-puzzling tree");
2201	if (puzzlemode == USERTREE) {
2202	if (utree_optn) fprintf(ofp, "1st user tree");
2203	else fprintf(ofp, "NJ tree");
2204	}
2205	fprintf(ofp, "):\n\n");
2206	if (bestratefound==0) findbestratecombination();
2207	printbestratecombination(ofp);
2208	}
2209
2210	fprintf(ofp, "\n\nSEQUENCE COMPOSITION (SEQUENCES IN INPUT ORDER)\n\n");
2211	fail = FALSE;
2212	fprintf(ofp, " 5%% chi-square test p-value\n");
2213	for (i = 0; i < Maxspc; i++) {
2214	fprintf(ofp, " ");
2215	fputid10(ofp, i);
2216	pval = homogentest(i);
2217	if ( pval < 0.05 ) fprintf(ofp, " failed ");
2218	else fprintf(ofp, " passed ");
2219	if (chi2fail) fail = TRUE;
2220	fprintf(ofp, " %6.2f%% ", pval*100.0);
2221	fprintf(ofp, "\n");
2222	}
2223	fprintf(ofp, "\n");
2224	fprintf(ofp, "The chi-square tests compares the ");
2225	if (data_optn == AMINOACID)
2226	fprintf(ofp, "amino acid");
2227	else if (data_optn == NUCLEOTIDE && SH_optn)
2228	fprintf(ofp, "doublet");
2229	else if (data_optn == NUCLEOTIDE && nuc_optn)
2230	fprintf(ofp, "nucleotide");
2231	else if (data_optn == BINARY)
2232	fprintf(ofp, "binary state");
2233	fprintf(ofp," composition of each sequence\n");
2234	fprintf(ofp, "to the frequency distribution assumed in the maximum likelihood model.\n");
2235	if (fail) {
2236	fprintf(ofp, "\nWARNING: Result of chi-square test may not be valid");
2237	fprintf(ofp, " because of small\nmaximum likelihood frequencies and");
2238	fprintf(ofp, " short sequence length!\n");
2239	}
2240	fprintf(ofp, "\n\nIDENTICAL SEQUENCES\n\n");
2241	fprintf(ofp, "The sequences in each of the following groups are all identical. To speed\n");
2242	fprintf(ofp, "up computation please remove all but one of each group from the data set.\n\n");
2243	findidenticals(ofp);
2244	fprintf(ofp, "\n\nMAXIMUM LIKELIHOOD DISTANCES\n\n");
2245	fprintf(ofp, "Maximum likelihood distances are computed using the ");
2246	fprintf(ofp, "selected model of\nsubstitution and rate heterogeneity.\n\n");
2247	putdistance(ofp);
2248	fprintf(ofp, "\nAverage distance (over all possible pairs of sequences): %.5f\n",
2249	averagedist() / 100.0);
2250
2251
2252	} /* if WRITEPARAMS) \|\| WRITEALL */
2253
2254	if ((part == WRITEREST) \|\| (part == WRITEALL)) {
2255
2256	if (puzzlemode == QUARTPUZ &&typ_optn == TREERECON_OPTN) {
2257	fprintf(ofp, "\n\nBAD QUARTET STATISTICS (SEQUENCES IN INPUT ORDER)\n\n");
2258	for (i = 0; i < Maxspc; i++) {
2259	fprintf(ofp, " ");
2260	fputid10(ofp, i);
2261	if (badqs > 0)
2262	fprintf(ofp, " [%lu] %6.2f%%\n", badtaxon[i], (double) (100 * badtaxon[i]) / (double) badqs);
2263	else
2264	fprintf(ofp, " [%lu]\n", badtaxon[i]);
2265	}
2266	fprintf(ofp, "\nThe number in square brackets indicates how often each sequence is\n");
2267	fprintf(ofp, "involved in one of the %lu completely unresolved quartets of the\n", badqs);
2268	fprintf(ofp, "quartet puzzling tree search.\n");
2269	if (badqs > 0)
2270	fprintf(ofp, "Additionally the according percentages are given.\n");
2271	}
2272
2273	if (typ_optn == TREERECON_OPTN) {
2274
2275	fprintf(ofp, "\n\nTREE SEARCH\n\n");
2276	if (puzzlemode == QUARTPUZ) {
2277	fprintf(ofp, "Quartet puzzling is used to choose from the possible tree topologies\n");
2278	fprintf(ofp, "and to simultaneously infer support values for internal branches.\n\n");
2279	fprintf(ofp, "Number of puzzling steps: %lu\n", Numtrial);
2280	fprintf(ofp, "Analysed quartets: %lu\n", Numquartets);
2281	fprintf(ofp, "Unresolved quartets: %lu (= %.1f%%)\n",
2282	badqs, (double) badqs / (double) Numquartets * 100);
2283	fprintf(ofp, "\nQuartet trees are based on %s maximum likelihood values\n",
2284	(approxqp ? "approximate" : "exact"));
2285	fprintf(ofp, "using the selected model of substitution and rate heterogeneity.\n\n\n");
2286	}
2287	if (puzzlemode == USERTREE) {
2288	fprintf(ofp, "%d tree topologies were specified by the user.\n", numutrees);
2289	}
2290	if (puzzlemode == PAIRDIST) {
2291	fprintf(ofp, "No tree search performed (maximum likelihood distances only).\n");
2292	}
2293
2294	if (puzzlemode == QUARTPUZ) {
2295	fprintf(ofp, "QUARTET PUZZLING TREE\n\n");
2296	fprintf(ofp, "Support for the internal branches of the unrooted quartet puzzling\n");
2297	fprintf(ofp, "tree topology is shown in percent.\n");
2298	if (consincluded == Maxspc - 3)
2299	fprintf(ofp,"\nThis quartet puzzling tree is completely resolved.\n");
2300	else
2301	fprintf(ofp,"\nThis quartet puzzling tree is not completely resolved!\n");
2302	fprintf(ofp, "\n\n");
2303	plotconsensustree(ofp);
2304	fprintf(ofp, "\n\nQuartet puzzling tree (in CLUSTAL W notation):\n\n");
2305	writeconsensustree(ofp);
2306	fprintf(ofp, "\n\nBIPARTITIONS\n\n");
2307	fprintf(ofp, "The following bipartitions occured at least once");
2308	fprintf(ofp, " in all intermediate\ntrees that have been generated ");
2309	fprintf(ofp, "in the %lu puzzling steps:\n\n", Numtrial);
2310	fprintf(ofp, "Bipartitions included in the quartet puzzling tree:\n");
2311	fprintf(ofp,
2312	"(bipartition with sequences in input order : number of times seen)\n\n");
2313	for (li = 0; li < consincluded; li++) {
2314	fprintf(ofp, " ");
2315	printsplit(ofp, splitfreqs[2*li+1]);
2316	fprintf(ofp, " : %lu\n", splitfreqs[2*li]);
2317	}
2318	if (consincluded == 0) fprintf(ofp, " None (no bipartition included)\n");
2319	fprintf(ofp, "\nBipartitions not included in the quartet puzzling tree:\n");
2320	fprintf(ofp,
2321	"(bipartition with sequences in input order : number of times seen)\n\n");
2322
2323	if (consincluded == numbiparts) {
2324	fprintf(ofp, " None (all bipartitions are included)\n");
2325	} else {
2326	/* print first 20 bipartions not included */
2327	for (li = consincluded; (li < numbiparts) && (li < consincluded + 20UL); li++) {
2328	fprintf(ofp, " ");
2329	printsplit(ofp, splitfreqs[2*li+1]);
2330	fprintf(ofp, " : %lu\n", splitfreqs[2*li]);
2331	}
2332	if ((li == consincluded + 20UL) && (li != numbiparts))
2333	fprintf(ofp, "\n(%lu other less frequent bipartitions not shown)\n",
2334	numbiparts - consincluded - 20UL);
2335	}
2336	fprintfsortedpstrees(ofp, psteptreelist, psteptreenum, psteptreesum, 0, 5.0);
2337	}
2338
2339	if (puzzlemode == QUARTPUZ) {
2340	fprintf(ofp, "\n\nMAXIMUM LIKELIHOOD BRANCH LENGTHS ON QUARTET");
2341	fprintf(ofp, " PUZZLING TREE (NO CLOCK)\n\nBranch lengths are computed using");
2342	fprintf(ofp, " the selected model of\nsubstitution and rate heterogeneity.\n\n\n");
2343	clockmode = 0; /* nonclocklike branch lengths */
2344	prtopology(ofp);
2345	fprintf(ofp, "\n");
2346	resulttree(ofp);
2347	fprintf(ofp, "\n\nQuartet puzzling tree with maximum likelihood branch lengths");
2348	fprintf(ofp, "\n(in CLUSTAL W notation):\n\n");
2349	fputphylogeny(ofp);
2350	if (compclock) {
2351	fprintf(ofp, "\n\nMAXIMUM LIKELIHOOD BRANCH LENGTHS OF QUARTET");
2352	fprintf(ofp, " PUZZLING TREE (WITH CLOCK)\n\nBranch lengths are computed using");
2353	fprintf(ofp, " the selected model of\nsubstitution and rate heterogeneity.\n");
2354	fprintf(ofp, "\nRoot located at branch: %d ", locroot+1);
2355	if (rootsearch == 0) fprintf(ofp, "(user specified)\n\n\n");
2356	if (rootsearch == 1) {
2357	fprintf(ofp, "(automatic search)");
2358	if (numbestroot > 1) fprintf(ofp, "- WARNING: %d best locations found! -", numbestroot);
2359	fprintf(ofp, "\n\n");
2360	fprintf(ofp, "If the automatic search misplaces the root please rerun the analysis\n");
2361	fprintf(ofp, "(rename \"outtree\" to \"intree\") and select location of root manually!");
2362	fprintf(ofp, "\n\n\n");
2363	}
2364	if (rootsearch == 2) fprintf(ofp, "(displayed outgroup)\n\n\n");
2365	clockmode = 1; /* clocklike branch lengths */
2366	prtopology(ofp);
2367	fprintf(ofp, "\n");
2368	fprintf(ofp, "\nTree drawn as unrooted tree for better ");
2369	fprintf(ofp, "comparison with non-clock tree!\n");
2370	resulttree(ofp);
2371	fprintf(ofp, "\n");
2372	resultheights(ofp);
2373	fprintf(ofp, "\n\nRooted quartet puzzling tree with clocklike");
2374	fprintf(ofp, " maximum likelihood branch lengths\n");
2375	fprintf(ofp, "(in CLUSTAL W notation):\n\n");
2376	fputrooted(ofp, locroot);
2377	}
2378
2379	if (compclock) {
2380	fprintf(ofp, "\n\nMOLECULAR CLOCK LIKELIHOOD RATIO TEST\n\n");
2381	fprintf(ofp, "log L without clock: %.2f (independent branch parameters: %d)\n",
2382	Ctree->lklhd, Numspc + Numibrnch);
2383	fprintf(ofp, "log L with clock: %.2f (independent branch parameters: %d)\n\n",
2384	Ctree->lklhdc, Numhts + 1);
2385	delta = 2.0*((Ctree->lklhd) - (Ctree->lklhdc));
2386	fprintf(ofp, "Likelihood ratio test statistic delta: %.2f\n", delta);
2387	df = Numspc + Numibrnch - Numhts - 1;
2388	fprintf(ofp, "Degress of freedom of chi-square distribution: %d\n", df);
2389
2390	pval = IncompleteGammaQ(df0.5, delta0.5);
2391
2392	fprintf(ofp, "Critical significance level: %.2f%%\n\n", pval*100.0);
2393	if (pval >= 0.05) {
2394	fprintf(ofp, "The simpler (clocklike) tree can not be rejected on a significance\n");
2395	fprintf(ofp, "level of 5%%. The log-likelihood of the more complex (no clock) tree\n");
2396	fprintf(ofp, "is not significantly increased.\n");
2397	} else {
2398	fprintf(ofp, "The simpler (clocklike) tree is rejected on a significance level\n");
2399	fprintf(ofp, "of 5%%. The log-likelihood of the more complex (no clock) tree is\n");
2400	fprintf(ofp, "significantly increased.\n");
2401	}
2402	fprintf(ofp, "\nPlease take care that the correct root is used!\n");
2403	}
2404
2405	}
2406	}
2407
2408	if (typ_optn == LIKMAPING_OPTN) {
2409
2410	fprintf(ofp, "\n\nLIKELIHOOD MAPPING ANALYSIS\n\n");
2411	fprintf(ofp, "Number of quartets: %lu", Numquartets);
2412	if (lmqts == 0) fprintf(ofp, " (all possible)\n");
2413	else fprintf(ofp, " (random choice)\n");
2414	fprintf(ofp, "\nQuartet trees are based on approximate maximum likelihood values\n");
2415	fprintf(ofp, "using the selected model of substitution and rate heterogeneity.\n\n\n");
2416	if (numclust == 1) {
2417	fprintf(ofp, "Sequences are not grouped in clusters.\n");
2418	} else {
2419	fprintf(ofp, "Sequences are grouped in %d clusters.\n", numclust);
2420	fprintf(ofp, "\nCluster a: %d sequences\n\n", clustA);
2421	for (i = 0; i < clustA; i++) {
2422	fprintf(ofp, " ");
2423	fputid(ofp, clusterA[i]);
2424	fprintf(ofp, "\n");
2425	}
2426	fprintf(ofp, "\nCluster b: %d sequences\n\n", clustB);
2427	for (i = 0; i < clustB; i++) {
2428	fprintf(ofp, " ");
2429	fputid(ofp, clusterB[i]);
2430	fprintf(ofp, "\n");
2431	}
2432	if (numclust > 2) {
2433	fprintf(ofp, "\nCluster c: %d sequences\n\n", clustC);
2434	for (i = 0; i < clustC; i++) {
2435	fprintf(ofp, " ");
2436	fputid(ofp, clusterC[i]);
2437	fprintf(ofp, "\n");
2438	}
2439	}
2440	if (numclust == 4) {
2441	fprintf(ofp, "\nCluster d: %d sequences\n\n", clustD);
2442	for (i = 0; i < clustD; i++) {
2443	fprintf(ofp, " ");
2444	fputid(ofp, clusterD[i]);
2445	fprintf(ofp, "\n");
2446	}
2447	}
2448	fprintf(ofp, "\nQuartets of sequences used in the likelihood");
2449	fprintf(ofp, " mapping analysis are generated\n");
2450	if (numclust == 2)
2451	fprintf(ofp, "by drawing two sequences from cluster a and two from cluster b.");
2452	if (numclust == 3)
2453	fprintf(ofp, "by drawing one sequence from clusters a and b and two from cluster c.");
2454	if (numclust == 4)
2455	fprintf(ofp, "by drawing one sequence from each of the clusters a, b, c, and d.");
2456	}
2457
2458	fprintf(ofp, "\n\nLIKELIHOOD MAPPING STATISTICS\n\n");
2459	fprintf(ofp, "Occupancies of the three areas 1, 2, 3:\n\n");
2460	if (numclust == 4)
2461	fprintf(ofp, " (a,b)-(c,d)\n");
2462	if (numclust == 3)
2463	fprintf(ofp, " (a,b)-(c,c)\n");
2464	if (numclust == 2)
2465	fprintf(ofp, " (a,a)-(b,b)\n");
2466	fprintf(ofp, " /\\\n");
2467	fprintf(ofp, " / \\\n");
2468	fprintf(ofp, " / \\\n");
2469	fprintf(ofp, " / 1 \\\n");
2470	fprintf(ofp, " / \\ / \\\n");
2471	fprintf(ofp, " / \\ / \\\n");
2472	fprintf(ofp, " / \\/ \\\n");
2473	fprintf(ofp, " / 3 : 2 \\\n");
2474	fprintf(ofp, " / : \\\n");
2475	fprintf(ofp, " /__________________\\\n");
2476	if (numclust == 4)
2477	fprintf(ofp, " (a,d)-(b,c) (a,c)-(b,d)\n");
2478	if (numclust == 3)
2479	fprintf(ofp, " (a,c)-(b,c) (a,c)-(b,c)\n");
2480	if (numclust == 2)
2481	fprintf(ofp, " (a,b)-(a,b) (a,b)-(a,b)\n");
2482	fprintf(ofp, "\n");
2483	fprintf(ofp, "Number of quartets in region 1: %lu (= %.1f%%)\n",
2484	ar1, (double) ar1*100.0/Numquartets);
2485	fprintf(ofp, "Number of quartets in region 2: %lu (= %.1f%%)\n",
2486	ar2, (double) ar2*100.0/Numquartets);
2487	fprintf(ofp, "Number of quartets in region 3: %lu (= %.1f%%)\n\n",
2488	ar3, (double) ar3*100.0/Numquartets);
2489	fprintf(ofp, "Occupancies of the seven areas 1, 2, 3, 4, 5, 6, 7:\n\n");
2490	if (numclust == 4)
2491	fprintf(ofp, " (a,b)-(c,d)\n");
2492	if (numclust == 3)
2493	fprintf(ofp, " (a,b)-(c,c)\n");
2494	if (numclust == 2)
2495	fprintf(ofp, " (a,a)-(b,b)\n");
2496	fprintf(ofp, " /\\\n");
2497	fprintf(ofp, " / \\\n");
2498	fprintf(ofp, " / 1 \\\n");
2499	fprintf(ofp, " / \\ / \\\n");
2500	fprintf(ofp, " / /\\ \\\n");
2501	fprintf(ofp, " / 6 / \\ 4 \\\n");
2502	fprintf(ofp, " / / 7 \\ \\\n");
2503	fprintf(ofp, " / \\ /______\\ / \\\n");
2504	fprintf(ofp, " / 3 : 5 : 2 \\\n");
2505	fprintf(ofp, " /__________________\\\n");
2506	if (numclust == 4)
2507	fprintf(ofp, " (a,d)-(b,c) (a,c)-(b,d)\n");
2508	if (numclust == 3)
2509	fprintf(ofp, " (a,c)-(b,c) (a,c)-(b,c)\n");
2510	if (numclust == 2)
2511	fprintf(ofp, " (a,b)-(a,b) (a,b)-(a,b)\n");
2512	fprintf(ofp, "\n");
2513	fprintf(ofp, "Number of quartets in region 1: %lu (= %.1f%%) left: %lu right: %lu\n",
2514	reg1, (double) reg1*100.0/Numquartets, reg1l, reg1r);
2515	fprintf(ofp, "Number of quartets in region 2: %lu (= %.1f%%) bottom: %lu top: %lu\n",
2516	reg2, (double) reg2*100.0/Numquartets, reg2d, reg2u);
2517	fprintf(ofp, "Number of quartets in region 3: %lu (= %.1f%%) bottom: %lu top: %lu\n",
2518	reg3, (double) reg3*100.0/Numquartets, reg3d, reg3u);
2519	fprintf(ofp, "Number of quartets in region 4: %lu (= %.1f%%) bottom: %lu top: %lu\n",
2520	reg4, (double) reg4*100.0/Numquartets, reg4d, reg4u);
2521	fprintf(ofp, "Number of quartets in region 5: %lu (= %.1f%%) left: %lu right: %lu\n",
2522	reg5, (double) reg5*100.0/Numquartets, reg5l, reg5r);
2523	fprintf(ofp, "Number of quartets in region 6: %lu (= %.1f%%) bottom: %lu top: %lu\n",
2524	reg6, (double) reg6*100.0/Numquartets, reg6d, reg6u);
2525	fprintf(ofp, "Number of quartets in region 7: %lu (= %.1f%%)\n",
2526	reg7, (double) reg7*100.0/Numquartets);
2527	}
2528
2529	} /* if WRITEREST) \|\| WRITEALL */
2530	}
2531
2532
2533	#if PARALLEL
2534	void writetimesstat(FILE *ofp)
2535	{
2536	int n;
2537	double cpusum = 0.0;
2538	double wallmax = 0.0;
2539	cputimes[0] = ((double)(cputimestop - cputimestart) / CLOCKS_PER_SEC);
2540	walltimes[0] = difftime(walltimestop, walltimestart);
2541	fullcpu = tarr.fullcpu;
2542	fulltime = tarr.fulltime;
2543	fullcputimes[0] = tarr.fullcpu;
2544	fullwalltimes[0] = tarr.fulltime;
2545	altcputimes[0] = tarr.cpu;
2546	altwalltimes[0] = tarr.time;
2547	fprintf(ofp, "\n\n\nPARALLEL LOAD STATISTICS\n\n");
2548
2549	fprintf(ofp, "The analysis was performed with %d parallel processes (1 master and \n", PP_NumProcs);
2550	fprintf(ofp, "%d worker processes).\n\n", PP_NumProcs-1);
2551	fprintf(ofp, "The following table the distribution of computation to the processes.\n");
2552	fprintf(ofp, "The first column gives the process number, where 0 is the master process.\n");
2553	fprintf(ofp, "The second and third column show the number of quartets computed (3 topologies \n");
2554	fprintf(ofp, "each) and the the number of scheduling blocks the came in. The last two columns \n");
2555	fprintf(ofp, "state the number of puzzling steps done by a process and number of scheduling \n");
2556	fprintf(ofp, "blocks.\n\n");
2557	fprintf(ofp, "process #quartets #chunks #puzzlings #chunks \n");
2558	fprintf(ofp, "-----------------------------------------------\n");
2559	for (n=0; n<PP_NumProcs; n++) {
2560	fprintf(ofp, "%6d %9d %7d %10d %7d \n", n,
2561	quartsent[n], quartsentn[n],
2562	splitsent[n], splitsentn[n]);
2563	} /* for */
2564	fprintf(ofp, "-----------------------------------------------\n");
2565	fprintf(ofp, " Sums: %9d %7d %10d %7d \n",
2566	PP_quartrecved, PP_quartrecvedn,
2567	PP_splitrecved, PP_splitrecvedn);
2568
2569	#ifdef TIMEDEBUG
2570	fprintf(ofp, "\n\nBelow the distribution of computing times (CPU and wallclock) per host is shown.\n");
2571	fprintf(ofp, "The times are shown in seconds, minutes, and hours. At the bottom of the table the\n");
2572	fprintf(ofp, "sum of CPU times and the maximum wallclock time is shown.\n\n");
2573	fprintf(ofp, "process CPU-time[s] [min] [hours] \| wallclock[s] [min] [hours] \n");
2574	fprintf(ofp, "----------------------------------------------------------------------------\n");
2575	for (n=0; n<PP_NumProcs; n++) {
2576
2577	# ifdef TIMEDEBUG
2578	fprintf(ofp, "%6d %11.1f %9.1f %9.1f \| %11.1f %9.1f %9.1f\n", n,
2579	cputimes[n], cputimes[n] /60, cputimes[n] /3600,
2580	walltimes[n], walltimes[n]/60, walltimes[n]/3600);
2581	# endif /* TIMEDEBUG */
2582
2583	fprintf(ofp, "%6d %11.1f %9.1f %9.1f \| %11.1f %9.1f %9.1f\n", n,
2584	fullcputimes[n], fullcputimes[n] /60, fullcputimes[n] /3600,
2585	fullwalltimes[n], fullwalltimes[n]/60, fullwalltimes[n]/3600);
2586
2587	# ifdef TIMEDEBUG
2588	fprintf(ofp, "%6d %11.1f %9.1f %9.1f \| %11.1f %9.1f %9.1f alt\n", n,
2589	altcputimes[n], altcputimes[n] /60, altcputimes[n] /3600,
2590	altwalltimes[n], altwalltimes[n]/60, altwalltimes[n]/3600);
2591	# endif /* TIMEDEBUG */
2592
2593	if (fullwalltimes[n] > wallmax) wallmax=fullwalltimes[n];
2594	cpusum += fullcputimes[n];
2595	} /* for */
2596	fprintf(ofp, "----------------------------------------------------------------------------\n");
2597	fprintf(ofp, "Sum/Max: %11.1f %9.1f %9.1f \| %11.1f %9.1f %9.1f \n",
2598	cpusum, cpusum/60, cpusum/3600, wallmax, wallmax/60, wallmax/3600);
2599	#else /* TIMEDEBUG */
2600	fprintf(ofp, "\n\nBelow the distribution of computing times (wallclock) per host is shown.\n");
2601	fprintf(ofp, "The times are shown in seconds, minutes, and hours. At the bottom of the table the\n");
2602	fprintf(ofp, "the maximum wallclock times is shown.\n\n");
2603	fprintf(ofp, "process wallclock[s] [min] [hours] \n");
2604	fprintf(ofp, "----------------------------------------------------------------------------\n");
2605	for (n=0; n<PP_NumProcs; n++) {
2606
2607	# ifdef TIMEDEBUG
2608	fprintf(ofp, "%6d %11.1f %9.1f %9.1f\n", n,
2609	walltimes[n], walltimes[n]/60, walltimes[n]/3600);
2610	# endif /* TIMEDEBUG */
2611
2612	fprintf(ofp, "%6d %11.1f %9.1f %9.1f\n", n,
2613	fullwalltimes[n], fullwalltimes[n]/60, fullwalltimes[n]/3600);
2614
2615	# ifdef TIMEDEBUG
2616	fprintf(ofp, "%6d %11.1f %9.1f %9.1f alt\n", n,
2617	altwalltimes[n], altwalltimes[n]/60, altwalltimes[n]/3600);
2618	# endif /* TIMEDEBUG */
2619
2620	if (fullwalltimes[n] > wallmax) wallmax=fullwalltimes[n];
2621	cpusum += fullcputimes[n];
2622	} /* for */
2623	fprintf(ofp, "----------------------------------------------------------------------------\n");
2624	fprintf(ofp, "Sum/Max: %11.1f %9.1f %9.1f \n",
2625	wallmax, wallmax/60, wallmax/3600);
2626	#endif /* TIMEDEBUG */
2627
2628	fullcpu = cpusum;
2629	fulltime = wallmax;
2630
2631	} /* writetimesstat */
2632	#endif
2633
2634
2635	/* write current user tree to file */
2636	void writecutree(FILE *ofp, int num)
2637	{
2638	int df;
2639	double pval, delta;
2640
2641
2642	if (typ_optn == TREERECON_OPTN) {
2643
2644	if (puzzlemode == USERTREE) {
2645	fprintf(ofp, "\n\nMAXIMUM LIKELIHOOD BRANCH LENGTHS OF USER");
2646	fprintf(ofp, " DEFINED TREE # %d (NO CLOCK)\n\nBranch lengths are computed using", num);
2647	fprintf(ofp, " the selected model of\nsubstitution and rate heterogeneity.\n\n\n");
2648	clockmode = 0; /* nonclocklike branch lengths */
2649	prtopology(ofp);
2650	fprintf(ofp, "\n");
2651	resulttree(ofp);
2652	fprintf(ofp, "\n\nUnrooted user defined tree with maximum likelihood branch lengths");
2653	fprintf(ofp, "\n(in CLUSTAL W notation):\n\n");
2654	fputphylogeny(ofp);
2655	if (compclock) {
2656	fprintf(ofp, "\n\nMAXIMUM LIKELIHOOD BRANCH LENGTHS OF USER");
2657	fprintf(ofp, " DEFINED TREE # %d (WITH CLOCK)\n\nBranch lengths are computed using", num);
2658	fprintf(ofp, " the selected model of\nsubstitution and rate heterogeneity.\n");
2659	fprintf(ofp, "\nRoot located at branch: %d ", locroot+1);
2660	if (rootsearch == 0) fprintf(ofp, "(user specified)\n\n\n");
2661	if (rootsearch == 1) {
2662	fprintf(ofp, "(automatic search)");
2663	if (numbestroot > 1) fprintf(ofp, "- WARNING: %d best locations found! -", numbestroot);
2664	fprintf(ofp, "\n\n");
2665	fprintf(ofp, "If the automatic search misplaces the root please rerun the analysis\n");
2666	fprintf(ofp, "and select location of root manually!");
2667	fprintf(ofp, "\n\n\n");
2668
2669	}
2670	if (rootsearch == 2) fprintf(ofp, "(displayed outgroup)\n\n\n");
2671	clockmode = 1; /* clocklike branch lengths */
2672	prtopology(ofp);
2673	fprintf(ofp, "\n");
2674	resulttree(ofp);
2675	fprintf(ofp, "\n");
2676	resultheights(ofp);
2677	fprintf(ofp, "\n\nRooted user defined tree with clocklike ");
2678	fprintf(ofp, "maximum likelihood branch lengths\n");
2679	fprintf(ofp, "(in CLUSTAL W notation):\n\n");
2680	fputrooted(ofp, locroot);
2681	}
2682
2683	if (compclock) {
2684	fprintf(ofp, "\n\nMOLECULAR CLOCK LIKELIHOOD RATIO TEST FOR USER TREE # %d\n\n", num);
2685	fprintf(ofp, "log L without clock: %.2f (independent branch parameters: %d)\n",
2686	Ctree->lklhd, Numspc + Numibrnch);
2687	fprintf(ofp, "log L with clock: %.2f (independent branch parameters: %d)\n\n",
2688	Ctree->lklhdc, Numhts + 1);
2689	delta = 2.0*((Ctree->lklhd) - (Ctree->lklhdc));
2690	fprintf(ofp, "Likelihood ratio test statistic delta: %.2f\n", delta);
2691	df = Numspc + Numibrnch - Numhts - 1;
2692	fprintf(ofp, "Degrees of freedom of chi-square distribution: %d\n", df);
2693
2694	pval = IncompleteGammaQ (df0.5, delta0.5);
2695
2696	fprintf(ofp, "Critical significance level: %.2f%%\n\n", pval*100.0);
2697	if (pval >= 0.05) {
2698	fprintf(ofp, "The simpler (clocklike) tree can not be rejected on a significance\n");
2699	fprintf(ofp, "level of 5%%. The log-likelihood of the more complex (no clock) tree\n");
2700	fprintf(ofp, "is not significantly increased.\n");
2701	} else {
2702	fprintf(ofp, "The simpler (clocklike) tree is rejected on a significance level\n");
2703	fprintf(ofp, "of 5%%. The log-likelihood of the more complex (no clock) tree is\n");
2704	fprintf(ofp, "significantly increased.\n");
2705	}
2706	fprintf(ofp, "\nPlease take care that the correct root is used!\n");
2707	}
2708	}
2709	}
2710	}
2711
2712
2713	/******************************************************************************/
2714	/* timer routines */
2715	/******************************************************************************/
2716
2717	/* start timer */
2718	void starttimer()
2719	{
2720	time(&time0);
2721	time1 = time0;
2722	}
2723
2724	/* check remaining time and print message if necessary */
2725	void checktimer(uli numqts)
2726	{
2727	double tc2, mintogo, minutes, hours;
2728
2729	time(&time2);
2730	if ( (time2 - time1) > 900) { /* generate message every 15 minutes */
2731	/* every 900 seconds */
2732	/* percentage of completed quartets */
2733	if (mflag == 0) {
2734	mflag = 1;
2735	FPRINTF(STDOUTFILE "\n");
2736	}
2737	tc2 = 100.*numqts/Numquartets;
2738	mintogo = (100.0-tc2) *
2739	(double) (time2-time0)/60.0/tc2;
2740	hours = floor(mintogo/60.0);
2741	minutes = mintogo - 60.0*hours;
2742	FPRINTF(STDOUTFILE "%.2f%%", tc2);
2743	FPRINTF(STDOUTFILE " completed (remaining");
2744	FPRINTF(STDOUTFILE " time: %.0f", hours);
2745	FPRINTF(STDOUTFILE " hours %.0f", minutes);
2746	FPRINTF(STDOUTFILE " minutes)\n");
2747	fflush(STDOUT);
2748	time1 = time2;
2749	}
2750
2751	}
2752
2753	/* check remaining time and print message if necessary */
2754	void checktimer2(uli numqts, uli all, int flag)
2755	{
2756	double tc2, mintogo, minutes, hours;
2757
2758	static time_t tt1;
2759	static time_t tt2;
2760
2761	if (flag == 1) {
2762	time(&tt1);
2763	time(&tt2);
2764	} else {
2765	time(&tt2);
2766	if ( (tt2 - tt1) > 900) { /* generate message every 15 minutes */
2767	/* every 900 seconds */
2768	/* percentage of completed quartets */
2769	if (mflag == 0) {
2770	mflag = 1;
2771	FPRINTF(STDOUTFILE "\n");
2772	}
2773	tc2 = 100.*numqts/Numquartets;
2774	mintogo = (100.0-tc2) *
2775	(double) (tt2-time0)/60.0/tc2;
2776	hours = floor(mintogo/60.0);
2777	minutes = mintogo - 60.0*hours;
2778	FPRINTF(STDOUTFILE "%.2f%%", tc2);
2779	FPRINTF(STDOUTFILE " completed (remaining");
2780	FPRINTF(STDOUTFILE " time: %.0f", hours);
2781	FPRINTF(STDOUTFILE " hours %.0f", minutes);
2782	FPRINTF(STDOUTFILE " minutes)\n");
2783	fflush(STDOUT);
2784	tt1 = tt2;
2785	}
2786	}
2787	}
2788
2789	void resetqblocktime(timearray_t *ta)
2790	{
2791	ta->quartcpu += ta->quartblockcpu;
2792	ta->quartblockcpu = 0.0;
2793	ta->quarttime += ta->quartblocktime;
2794	ta->quartblocktime = 0.0;
2795	} /* resetqblocktime */
2796
2797
2798	void resetpblocktime(timearray_t *ta)
2799	{
2800	ta->puzzcpu += ta->puzzblockcpu;
2801	ta->puzzblockcpu = 0.0;
2802	ta->puzztime += ta->puzzblocktime;
2803	ta->puzzblocktime = 0.0;
2804	} /* resetpblocktime */
2805
2806
2807	#ifdef TIMEDEBUG
2808	void printtimearr(timearray_t *ta)
2809	{
2810	# if ! PARALLEL
2811	int PP_Myid;
2812	PP_Myid = -1;
2813	# endif
2814	printf("(%2d) MMCPU: %11ld / %11ld \n", PP_Myid, ta->maxcpu, ta->mincpu);
2815	printf("(%2d) CTick: %11.6f [tks] / %11.6f [s] \n", PP_Myid, ta->mincputick, ta->mincputicktime);
2816
2817	printf("(%2d) MMTIM: %11ld / %11ld \n", PP_Myid, ta->maxtime, ta->mintime);
2818
2819	printf("(%2d) Mxblk: %11.6e / %11.6e \n", PP_Myid, ta->maxcpublock, ta->maxtimeblock);
2820	printf("(%2d) Mnblk: %11.6e / %11.6e \n", PP_Myid, ta->mincpublock, ta->mintimeblock);
2821
2822	printf("(%2d) Gnrl: %11.6e / %11.6e \n", PP_Myid, ta->generalcpu, ta->generaltime);
2823	printf("(%2d) Optn: %11.6e / %11.6e \n", PP_Myid, ta->optionscpu, ta->optionstime);
2824	printf("(%2d) Estm: %11.6e / %11.6e \n", PP_Myid, ta->paramestcpu, ta->paramesttime);
2825	printf("(%2d) Qurt: %11.6e / %11.6e \n", PP_Myid, ta->quartcpu, ta->quarttime);
2826	printf("(%2d) QBlk: %11.6e / %11.6e \n", PP_Myid, ta->quartblockcpu, ta->quartblocktime);
2827	printf("(%2d) QMax: %11.6e / %11.6e \n", PP_Myid, ta->quartmaxcpu, ta->quartmaxtime);
2828	printf("(%2d) QMin: %11.6e / %11.6e \n", PP_Myid, ta->quartmincpu, ta->quartmintime);
2829
2830	printf("(%2d) Puzz: %11.6e / %11.6e \n", PP_Myid, ta->puzzcpu, ta->puzztime);
2831	printf("(%2d) PBlk: %11.6e / %11.6e \n", PP_Myid, ta->puzzblockcpu, ta->puzzblocktime);
2832	printf("(%2d) PMax: %11.6e / %11.6e \n", PP_Myid, ta->puzzmaxcpu, ta->puzzmaxtime);
2833	printf("(%2d) PMin: %11.6e / %11.6e \n", PP_Myid, ta->puzzmincpu, ta->puzzmintime);
2834
2835	printf("(%2d) Tree: %11.6e / %11.6e \n", PP_Myid, ta->treecpu, ta->treetime);
2836	printf("(%2d) TBlk: %11.6e / %11.6e \n", PP_Myid, ta->treeblockcpu, ta->treeblocktime);
2837	printf("(%2d) TMax: %11.6e / %11.6e \n", PP_Myid, ta->treemaxcpu, ta->treemaxtime);
2838	printf("(%2d) TMin: %11.6e / %11.6e \n", PP_Myid, ta->treemincpu, ta->treemintime);
2839
2840	printf("(%2d) C/T : %11.6e / %11.6e \n", PP_Myid,
2841	(ta->generalcpu + ta->optionscpu + ta->paramestcpu + ta->quartblockcpu + ta->puzzblockcpu + ta->treeblockcpu),
2842	(ta->generaltime + ta->optionstime + ta->paramesttime + ta->quartblocktime + ta->puzzblocktime + ta->treeblocktime));
2843	printf("(%2d) CPU: %11.6e / Time: %11.6e \n", PP_Myid, ta->cpu, ta->time);
2844	printf("(%2d) aCPU: %11.6e / aTime: %11.6e \n", PP_Myid, ta->fullcpu, ta->fulltime);
2845
2846	} /* printtimearr */
2847	#endif /* TIMEDEBUG */
2848
2849	char *jtype [7];
2850
2851	void inittimearr(timearray_t *ta)
2852	{
2853	clock_t c0, c1, c2;
2854
2855	jtype[OVERALL] = "OVERALL";
2856	jtype[GENERAL] = "GENERAL";
2857	jtype[OPTIONS] = "OPTIONS";
2858	jtype[PARAMEST] = "PARAMeter ESTimation";
2859	jtype[QUARTETS] = "QUARTETS";
2860	jtype[PUZZLING] = "PUZZLING steps";
2861	jtype[TREEEVAL] = "TREE EVALuation";
2862	ta->currentjob = GENERAL;
2863
2864	c1 = clock();
2865	c2 = clock();
2866	while (c1 == c2)
2867	c2 = clock();
2868	ta->mincputick = (double)(c2 - c1);
2869	ta->mincputicktime = ((double)(c2 - c1))/CLOCKS_PER_SEC;
2870
2871	ta->tempcpu = clock();
2872	ta->tempcpustart = ta->tempcpu;
2873	ta->tempfullcpu = ta->tempcpu;
2874	time(&(ta->temptime));
2875	ta->temptimestart = ta->temptime;
2876	ta->tempfulltime = ta->temptime;
2877
2878	c0=0; c1=0; c2=(clock_t)((2 * c1) + 1);;
2879	while (c1 < c2) {
2880	c0 = c1;
2881	c1 = c2;
2882	c2 = (clock_t)((2 * c1) + 1);
2883	}
2884	if (c1 == c2) ta->maxcpu=c0;
2885	if (c1 > c2) ta->maxcpu=c1;
2886
2887	c0=0; c1=0; c2=(clock_t)((2 * c1) - 1);
2888	while (c1 > c2) {
2889	c0 = c1;
2890	c1 = c2;
2891	c2 = (clock_t)((2 * c1) - 1);
2892	}
2893	if (c1 == c2) ta->mincpu=c0;
2894	if (c1 < c2) ta->mincpu=c1;
2895
2896
2897
2898	ta->maxtime = 0;
2899	ta->mintime = 0;
2900
2901	ta->maxcpublock = 0;
2902	ta->mincpublock = DBL_MAX;
2903	ta->maxtimeblock = 0;
2904	ta->mintimeblock = DBL_MAX;
2905
2906	ta->cpu = 0.0;
2907	ta->time = 0.0;
2908
2909	ta->fullcpu = 0.0;
2910	ta->fulltime = 0.0;
2911
2912	ta->generalcpu = 0.0;
2913	ta->optionscpu = 0.0;
2914	ta->paramestcpu = 0.0;
2915	ta->quartcpu = 0.0;
2916	ta->quartblockcpu = 0.0;
2917	ta->quartmaxcpu = 0.0;
2918	ta->quartmincpu = ((double) ta->maxcpu)/CLOCKS_PER_SEC;
2919	ta->puzzcpu = 0.0;
2920	ta->puzzblockcpu = 0.0;
2921	ta->puzzmaxcpu = 0.0;
2922	ta->puzzmincpu = ((double) ta->maxcpu)/CLOCKS_PER_SEC;
2923	ta->treecpu = 0.0;
2924	ta->treeblockcpu = 0.0;
2925	ta->treemaxcpu = 0.0;
2926	ta->treemincpu = ((double) ta->maxcpu)/CLOCKS_PER_SEC;
2927
2928	ta->generaltime = 0.0;
2929	ta->optionstime = 0.0;
2930	ta->paramesttime = 0.0;
2931	ta->quarttime = 0.0;
2932	ta->quartblocktime = 0.0;
2933	ta->quartmaxtime = 0.0;
2934	ta->quartmintime = DBL_MAX;
2935	ta->puzztime = 0.0;
2936	ta->puzzblocktime = 0.0;
2937	ta->puzzmaxtime = 0.0;
2938	ta->puzzmintime = DBL_MAX;
2939	ta->treetime = 0.0;
2940	ta->treeblocktime = 0.0;
2941	ta->treemaxtime = 0.0;
2942	ta->treemintime = DBL_MAX;
2943	} /* inittimearr */
2944
2945
2946	/***************/
2947
2948	void addup(int jobtype, clock_t c1, clock_t c2, time_t t1, time_t t2, timearray_t *ta)
2949	{
2950	double c,
2951	t;
2952
2953	if (t2 != t1) t = difftime(t2, t1);
2954	else t = 0.0;
2955
2956	if (c2 < c1)
2957	c = ((double)(c2 - ta->mincpu))/CLOCKS_PER_SEC +
2958	((double)(ta->maxcpu - c1))/CLOCKS_PER_SEC;
2959	else
2960	c = ((double)(c2 - c1))/CLOCKS_PER_SEC;
2961
2962	if (jobtype != OVERALL) {
2963
2964	if (ta->mincpublock > c) ta->mincpublock = c;
2965	if (ta->maxcpublock < c) ta->maxcpublock = c;
2966	if (ta->mintimeblock > t) ta->mintimeblock = t;
2967	if (ta->maxtimeblock < t) ta->maxtimeblock = t;
2968
2969	switch (jobtype) {
2970	case GENERAL: ta->generalcpu += c;
2971	ta->generaltime += t;
2972	break;
2973	case OPTIONS: ta->optionscpu += c;
2974	ta->optionstime += t;
2975	break;
2976	case PARAMEST: ta->paramestcpu += c;
2977	ta->paramesttime += t;
2978	break;
2979	case QUARTETS: ta->quartblockcpu += c;
2980	ta->quartblocktime += t;
2981	if (ta->quartmincpu > c) ta->quartmincpu = c;
2982	if (ta->quartmaxcpu < c) ta->quartmaxcpu = c;
2983	if (ta->quartmintime > t) ta->quartmintime = t;
2984	if (ta->quartmaxtime < t) ta->quartmaxtime = t;
2985	break;
2986	case PUZZLING: ta->puzzblockcpu += c;
2987	ta->puzzblocktime += t;
2988	if (ta->puzzmincpu > c) ta->puzzmincpu = c;
2989	if (ta->puzzmaxcpu < c) ta->puzzmaxcpu = c;
2990	if (ta->puzzmintime > t) ta->puzzmintime = t;
2991	if (ta->puzzmaxtime < t) ta->puzzmaxtime = t;
2992	break;
2993	case TREEEVAL: ta->treeblockcpu += c;
2994	ta->treeblocktime += t;
2995	if (ta->treemincpu > c) ta->treemincpu = c;
2996	if (ta->treemaxcpu < c) ta->treemaxcpu = c;
2997	if (ta->treemintime > t) ta->treemintime = t;
2998	if (ta->treemaxtime < t) ta->treemaxtime = t;
2999	break;
3000	}
3001	ta->cpu += c;
3002	ta->time += t;
3003
3004	} else {
3005	ta->fullcpu += c;
3006	ta->fulltime += t;
3007	}
3008
3009	# ifdef TIMEDEBUG
3010	{
3011	# if ! PARALLEL
3012	int PP_Myid = -1;
3013	# endif /* !PARALLEL */
3014	printf("(%2d) CPU: +%10.6f / Time: +%10.6f (%s)\n", PP_Myid, c, t, jtype[jobtype]);
3015	printf("(%2d) CPU: %11.6f / Time: %11.6f (%s)\n", PP_Myid, ta->cpu, ta->time, jtype[jobtype]);
3016	printf("(%2d) CPU: %11.6f / Time: %11.6f (%s)\n", PP_Myid, ta->fullcpu, ta->fulltime, jtype[jobtype]);
3017	}
3018	# endif /* TIMEDEBUG */
3019	} /* addup */
3020
3021
3022	/***************/
3023
3024
3025	void addtimes(int jobtype, timearray_t *ta)
3026	{
3027	clock_t tempc;
3028	time_t tempt;
3029
3030	time(&tempt);
3031	tempc = clock();
3032
3033	if ((tempc < ta->tempfullcpu) \|\| (jobtype == OVERALL)) { /* CPU counter overflow for overall time */
3034	addup(OVERALL, ta->tempfullcpu, tempc, ta->tempfulltime, tempt, ta);
3035	ta->tempfullcpu = tempc;
3036	ta->tempfulltime = tempt;
3037	if (jobtype == OVERALL) {
3038	addup(ta->currentjob, ta->tempcpustart, tempc, ta->temptimestart, tempt, ta);
3039	ta->tempcpustart = ta->tempcpu;
3040	ta->tempcpu = tempc;
3041	ta->temptimestart = ta->temptime;
3042	ta->temptime = tempt;
3043	}
3044	}
3045
3046	if((jobtype != ta->currentjob) && (jobtype != OVERALL)) { /* change of job type */
3047	addup(ta->currentjob, ta->tempcpustart, ta->tempcpu, ta->temptimestart, ta->temptime, ta);
3048	ta->tempcpustart = ta->tempcpu;
3049	ta->tempcpu = tempc;
3050	ta->temptimestart = ta->temptime;
3051	ta->temptime = tempt;
3052	ta->currentjob = jobtype;
3053	}
3054
3055	if (tempc < ta->tempcpustart) { /* CPU counter overflow */
3056	addup(jobtype, ta->tempcpustart, tempc, ta->temptimestart, tempt, ta);
3057	ta->tempcpustart = ta->tempcpu;
3058	ta->tempcpu = tempc;
3059	ta->temptimestart = ta->temptime;
3060	ta->temptime = tempt;
3061	}
3062
3063	} /* addtimes */
3064
3065
3066
3067	/******************************************************************************/
3068
3069	/* estimate parameters of substitution process and rate heterogeneity - no tree
3070	n-taxon tree is not needed because of quartet method or NJ tree topology */
3071	void estimateparametersnotree()
3072	{
3073	int it, nump, change;
3074	double TSold, YRold, FIold, GEold;
3075
3076	it = 0;
3077	nump = 0;
3078
3079	/* count number of parameters */
3080	if (data_optn == NUCLEOTIDE && optim_optn) nump++;
3081	if (fracinv_optim \|\| grate_optim) nump++;
3082
3083	do { /* repeat until nothing changes any more */
3084	it++;
3085	change = FALSE;
3086
3087	/* optimize substitution parameters */
3088	if (data_optn == NUCLEOTIDE && optim_optn) {
3089
3090	TSold = TSparam;
3091	YRold = YRparam;
3092
3093
3094	/*
3095	* optimize
3096	*/
3097
3098	FPRINTF(STDOUTFILE "Optimizing missing substitution process parameters\n");
3099	fflush(STDOUT);
3100
3101	if (qcalg_optn) { /* quartet sampling */
3102	optimseqevolparamsq();
3103	} else { /* NJ tree */
3104	tmpfp = tmpfile();
3105	njtree(tmpfp);
3106	rewind(tmpfp);
3107	readusertree(tmpfp);
3108	closefile(tmpfp);
3109	optimseqevolparamst();
3110	}
3111
3112	computedistan(); /* update ML distances */
3113
3114	/* same tolerance as 1D minimization */
3115	if ((fabs(TSparam - TSold) > 3.3*PEPS1) \|\|
3116	(fabs(YRparam - YRold) > 3.3*PEPS1)
3117	) change = TRUE;
3118
3119	}
3120
3121	/* optimize rate heterogeneity variables */
3122	if (fracinv_optim \|\| grate_optim) {
3123
3124	FIold = fracinv;
3125	GEold = Geta;
3126
3127
3128	/*
3129	* optimize
3130	*/
3131
3132	FPRINTF(STDOUTFILE "Optimizing missing rate heterogeneity parameters\n");
3133	fflush(STDOUT);
3134	/* compute NJ tree */
3135	tmpfp = tmpfile();
3136	njtree(tmpfp);
3137	/* use NJ tree topology to estimate parameters */
3138	rewind(tmpfp);
3139	readusertree(tmpfp);
3140	closefile(tmpfp);
3141
3142	optimrateparams();
3143	computedistan(); /* update ML distances */
3144
3145
3146	/* same tolerance as 1D minimization */
3147	if ((fabs(fracinv - FIold) > 3.3*PEPS2) \|\|
3148	(fabs(Geta - GEold) > 3.3*PEPS2)
3149	) change = TRUE;
3150
3151	}
3152
3153	if (nump == 1) return;
3154
3155	} while (it != MAXITS && change);
3156
3157	return;
3158	}
3159
3160
3161	/* estimate parameters of substitution process and rate heterogeneity - tree
3162	same as above but here the n-taxon tree is already in memory */
3163	void estimateparameterstree()
3164	{
3165	int it, nump, change;
3166	double TSold, YRold, FIold, GEold;
3167
3168	it = 0;
3169	nump = 0;
3170
3171	/* count number of parameters */
3172	if (data_optn == NUCLEOTIDE && optim_optn) nump++;
3173	if (fracinv_optim \|\| grate_optim) nump++;
3174
3175	do { /* repeat until nothing changes any more */
3176	it++;
3177	change = FALSE;
3178
3179	/* optimize substitution process parameters */
3180	if (data_optn == NUCLEOTIDE && optim_optn) {
3181
3182	TSold = TSparam;
3183	YRold = YRparam;
3184
3185
3186	/*
3187	* optimize
3188	*/
3189
3190	FPRINTF(STDOUTFILE "Optimizing missing substitution process parameters\n");
3191	fflush(STDOUT);
3192	optimseqevolparamst();
3193	computedistan(); /* update ML distances */
3194
3195
3196	/* same tolerance as 1D minimization */
3197	if ((fabs(TSparam - TSold) > 3.3*PEPS1) \|\|
3198	(fabs(YRparam - YRold) > 3.3*PEPS1)
3199	) change = TRUE;
3200
3201	}
3202
3203	/* optimize rate heterogeneity variables */
3204	if (fracinv_optim \|\| grate_optim) {
3205
3206	FIold = fracinv;
3207	GEold = Geta;
3208
3209
3210	/*
3211	* optimize
3212	*/
3213
3214	FPRINTF(STDOUTFILE "Optimizing missing rate heterogeneity parameters\n");
3215	fflush(STDOUT);
3216	optimrateparams();
3217	computedistan(); /* update ML distances */
3218
3219
3220	/* same tolerance as 1D minimization */
3221	if ((fabs(fracinv - FIold) > 3.3*PEPS2) \|\|
3222	(fabs(Geta - GEold) > 3.3*PEPS2)
3223	) change = TRUE;
3224
3225	}
3226
3227	if (nump == 1) return;
3228
3229	} while (it != MAXITS && change);
3230
3231	return;
3232	}
3233
3234
3235	/******************************************************************************/
3236	/* exported from main */
3237	/******************************************************************************/
3238
3239	void compute_quartlklhds(int a, int b, int c, int d, double d1, double d2, double *d3, int approx)
3240	{
3241	if (approx == APPROX) {
3242
3243	d1 = quartet_alklhd(a,b, c,d); / (a,b)-(c,d) */
3244	d2 = quartet_alklhd(a,c, b,d); / (a,c)-(b,d) */
3245	d3 = quartet_alklhd(a,d, b,c); / (a,d)-(b,c) */
3246
3247	} else /* approx == EXACT */ {
3248
3249	d1 = quartet_lklhd(a,b, c,d); / (a,b)-(c,d) */
3250	d2 = quartet_lklhd(a,c, b,d); / (a,c)-(b,d) */
3251	d3 = quartet_lklhd(a,d, b,c); / (a,d)-(b,c) */
3252
3253	}
3254	}
3255
3256	/***************************************************************/
3257
3258	void recon_tree()
3259	{
3260	int i;
3261	# if ! PARALLEL
3262	int a, b, c;
3263	uli nq;
3264	double tc2, mintogo, minutes, hours;
3265	# endif
3266
3267	/* allocate memory for taxon list of bad quartets */
3268	badtaxon = new_ulivector(Maxspc);
3269	for (i = 0; i < Maxspc; i++) badtaxon[i] = 0;
3270
3271	/* allocate variable used for randomizing input order */
3272	trueID = new_ivector(Maxspc);
3273
3274	/* allocate memory for quartets */
3275	quartetinfo = mallocquartets(Maxspc);
3276
3277	/* prepare for consensus tree analysis */
3278	initconsensus();
3279
3280	if (!(readquart_optn) \|\| (readquart_optn && savequart_optn)) {
3281	/* compute quartets */
3282	FPRINTF(STDOUTFILE "Computing quartet maximum likelihood trees\n");
3283	fflush(STDOUT);
3284	computeallquartets();
3285	}
3286
3287	if (savequart_optn)
3288	writeallquarts(Maxspc, ALLQUART, quartetinfo);
3289	if (readquart_optn) {
3290	int xx1, xx2, xx3, xx4, count;
3291	readallquarts (Maxspc, ALLQUART, quartetinfo);
3292	if (show_optn) { /* list all unresolved quartets */
3293	openfiletowrite(&unresfp, UNRESOLVED, "unresolved quartet trees");
3294	fprintf(unresfp, "List of all completely unresolved quartets:\n\n");
3295	}
3296
3297	/* initialize bad quartet memory */
3298	for (count = 0; count < Maxspc; count++) badtaxon[count] = 0;
3299	badqs = 0;
3300
3301	for (xx4 = 3; xx4 < Maxspc; xx4++)
3302	for (xx3 = 2; xx3 < xx4; xx3++)
3303	for (xx2 = 1; xx2 < xx3; xx2++)
3304	for (xx1 = 0; xx1 < xx2; xx1++) {
3305	if (readquartet(xx1, xx2, xx3, xx4) == 7) {
3306	badqs++;
3307	badtaxon[xx1]++;
3308	badtaxon[xx2]++;
3309	badtaxon[xx3]++;
3310	badtaxon[xx4]++;
3311	if (show_optn) {
3312	fputid10(unresfp, xx1);
3313	fprintf(unresfp, " ");
3314	fputid10(unresfp, xx2);
3315	fprintf(unresfp, " ");
3316	fputid10(unresfp, xx3);
3317	fprintf(unresfp, " ");
3318	fputid (unresfp, xx4);
3319	fprintf(unresfp, "\n");
3320	}
3321	}
3322	} /* end for xx4; for xx3; for xx2; for xx1 */
3323	if (show_optn) /* list all unresolved quartets */
3324	fclose(unresfp);
3325	} /* readquart_optn */
3326
3327	# if PARALLEL
3328	PP_SendAllQuarts(numquarts(Maxspc), quartetinfo);
3329	# endif /* PARALLEL */
3330
3331	FPRINTF(STDOUTFILE "Computing quartet puzzling tree\n");
3332	fflush(STDOUT);
3333
3334	/* start timer - percentage of completed trees */
3335	time(&time0);
3336	time1 = time0;
3337	mflag = 0;
3338
3339	/* open file for chronological list of puzzling step trees */
3340	if((listqptrees == PSTOUT_LIST) \|\| (listqptrees == PSTOUT_LISTORDER))
3341	openfiletowrite(&qptlist, OUTPTLIST, "puzzling step trees (chonological)");
3342
3343	# if PARALLEL
3344	{
3345	PP_SendDoPermutBlock(Numtrial);
3346	}
3347	# else
3348	addtimes(GENERAL, &tarr);
3349	for (Currtrial = 0; Currtrial < Numtrial; Currtrial++) {
3350
3351	/* randomize input order */
3352	chooser(Maxspc, Maxspc, trueID);
3353
3354	/* initialize tree */
3355	inittree();
3356
3357	/* adding all other leafs */
3358	for (i = 3; i < Maxspc; i++) {
3359
3360	/* clear all edgeinfos */
3361	resetedgeinfo();
3362
3363	/* clear counter of quartets */
3364	nq = 0;
3365
3366	/*
3367	* core of quartet puzzling algorithm
3368	*/
3369
3370	for (a = 0; a < nextleaf - 2; a++)
3371	for (b = a + 1; b < nextleaf - 1; b++)
3372	for (c = b + 1; c < nextleaf; c++) {
3373
3374	/* check which two _leaves_ out of a, b, c
3375	are closer related to each other than
3376	to leaf i according to a least squares
3377	fit of the continuous Baysian weights to the
3378	seven trivial "attractive regions". We assign
3379	a score of 1 to all edges between these two leaves
3380	chooseA and chooseB */
3381
3382	checkquartet(a, b, c, i);
3383	incrementedgeinfo(chooseA, chooseB);
3384
3385	nq++;
3386
3387	/* generate message every 15 minutes */
3388
3389	/* check timer */
3390	time(&time2);
3391	if ( (time2 - time1) > 900) {
3392	/* every 900 seconds */
3393	/* percentage of completed trees */
3394	if (mflag == 0) {
3395	FPRINTF(STDOUTFILE "\n");
3396	mflag = 1;
3397	}
3398	tc2 = 100.0*Currtrial/Numtrial +
3399	100.0*nq/Numquartets/Numtrial;
3400	mintogo = (100.0-tc2) *
3401	(double) (time2-time0)/60.0/tc2;
3402	hours = floor(mintogo/60.0);
3403	minutes = mintogo - 60.0*hours;
3404	FPRINTF(STDOUTFILE "%2.2f%%", tc2);
3405	FPRINTF(STDOUTFILE " completed (remaining");
3406	FPRINTF(STDOUTFILE " time: %.0f", hours);
3407	FPRINTF(STDOUTFILE " hours %.0f", minutes);
3408	FPRINTF(STDOUTFILE " minutes)\n");
3409	fflush(STDOUT);
3410	time1 = time2;
3411	}
3412	}
3413
3414	/* find out which edge has the lowest edgeinfo */
3415	minimumedgeinfo();
3416
3417	/* add the next leaf on minedge */
3418	addnextleaf(minedge);
3419	}
3420
3421	/* compute bipartitions of current tree */
3422	computebiparts();
3423	makenewsplitentries();
3424
3425	{
3426	int *ctree, startnode;
3427	char *trstr;
3428	treelistitemtype *treeitem;
3429	ctree = initctree();
3430	copytree(ctree);
3431	startnode = sortctree(ctree);
3432	trstr=sprintfctree(ctree, psteptreestrlen);
3433
3434
3435	treeitem = addtree2list(&trstr, 1, &psteptreelist, &psteptreenum, &psteptreesum);
3436
3437	if((listqptrees == PSTOUT_LIST)
3438	\|\| (listqptrees == PSTOUT_LISTORDER)) {
3439	/* print: order no/# topol per this id/tree id/sum of unique topologies/sum of trees so far */
3440	fprintf(qptlist, "%ld.\t1\t%d\t%d\t%d\t%d\n",
3441	Currtrial + 1, (treeitem).count, (treeitem).id, psteptreenum, psteptreesum);
3442	}
3443
3444	# ifdef VERBOSE1
3445	printf("%s\n", trstr);
3446	printfsortedpstrees(psteptreelist);
3447	# endif
3448	freectree(&ctree);
3449	}
3450
3451
3452
3453	/* free tree before building the next tree */
3454	freetree();
3455
3456	addtimes(PUZZLING, &tarr);
3457	}
3458	# endif /* PARALLEL */
3459
3460	/* close file for list of puzzling step trees */
3461	if((listqptrees == PSTOUT_LIST) \|\| (listqptrees == PSTOUT_LISTORDER))
3462	closefile(qptlist);
3463
3464	if (mflag == 1) FPRINTF(STDOUTFILE "\n");
3465
3466	/* garbage collection */
3467	free(splitcomp);
3468	free_ivector(trueID);
3469
3470	# if ! PARALLEL
3471	free_cmatrix(biparts);
3472	# endif /* PARALLEL */
3473
3474	freequartets();
3475
3476	/* compute majority rule consensus tree */
3477	makeconsensus();
3478
3479	/* write consensus tree to tmp file */
3480	tmpfp = tmpfile();
3481	writeconsensustree(tmpfp);
3482	} /* recon_tree */
3483
3484	/***************************************************************/
3485
3486	void map_lklhd()
3487	{
3488	int i, a, a1, a2, b, b1, b2, c, c1, c2, d;
3489	uli nq;
3490	double logs[3], d1, d2, d3, temp;
3491	ivector qts, mlorder, gettwo;
3492	/* reset variables */
3493	ar1 = ar2 = ar3 = 0;
3494	reg1 = reg2 = reg3 = reg4 = reg5 = reg6 = reg7 = 0;
3495	reg1l = reg1r = reg2u = reg2d = reg3u = reg3d = reg4u =
3496	reg4d = reg5l = reg5r = reg6u = reg6d = 0;
3497
3498	/* place for random quartet */
3499	qts = new_ivector(4);
3500
3501	/* initialize output file */
3502	openfiletowrite(&trifp, TRIANGLE, "Postscript output");
3503	initps(trifp);
3504	FPRINTF(STDOUTFILE "Performing likelihood mapping analysis\n");
3505	fflush(STDOUT);
3506
3507	/* start timer */
3508	starttimer();
3509	nq = 0;
3510	mflag = 0;
3511
3512	addtimes(GENERAL, &tarr);
3513	if (lmqts == 0) { /* all possible quartets */
3514
3515	if (numclust == 4) { /* four-cluster analysis */
3516
3517	for (a = 0; a < clustA; a++)
3518	for (b = 0; b < clustB; b++)
3519	for (c = 0; c < clustC; c++)
3520	for (d = 0; d < clustD; d++) {
3521
3522	nq++;
3523
3524	/* check timer */
3525	checktimer(nq);
3526
3527	/* maximum likelihood values */
3528	/* approximate ML is sufficient */
3529	compute_quartlklhds(clusterA[a],clusterB[b],clusterC[c],clusterD[d],&d1,&d2,&d3, APPROX);
3530
3531	/* draw point for LM analysis */
3532	makelmpoint(trifp, d1, d2, d3);
3533	addtimes(QUARTETS, &tarr);
3534
3535	}
3536	}
3537
3538	if (numclust == 3) { /* three-cluster analysis */
3539
3540	gettwo = new_ivector(2);
3541
3542	for (a = 0; a < clustA; a++)
3543	for (b = 0; b < clustB; b++)
3544	for (c1 = 0; c1 < clustC-1; c1++)
3545	for (c2 = c1+1; c2 < clustC; c2++) {
3546
3547	nq++;
3548
3549	/* check timer */
3550	checktimer(nq);
3551
3552	/* maximum likelihood values */
3553	/* approximate ML is sufficient */
3554	compute_quartlklhds(clusterA[a],clusterB[b],clusterC[c1],clusterC[c2],&d1,&d2,&d3, APPROX);
3555
3556	/* randomize order of d2 and d3 */
3557	if (randominteger(2) == 1) {
3558	temp = d3;
3559	d3 = d2;
3560	d2 = temp;
3561	}
3562
3563	/* draw point for LM analysis */
3564	makelmpoint(trifp, d1, d2, d3);
3565	addtimes(QUARTETS, &tarr);
3566
3567	}
3568	free_ivector(gettwo);
3569	}
3570
3571	if (numclust == 2) { /* two-cluster analysis */
3572
3573	gettwo = new_ivector(2);
3574
3575	for (a1 = 0; a1 < clustA-1; a1++)
3576	for (a2 = a1+1; a2 < clustA; a2++)
3577	for (b1 = 0; b1 < clustB-1; b1++)
3578	for (b2 = b1+1; b2 < clustB; b2++) {
3579
3580	nq++;
3581
3582	/* check timer */
3583	checktimer(nq);
3584
3585	/* maximum likelihood values */
3586	/* approximate ML is sufficient */
3587	compute_quartlklhds(clusterA[a1],clusterA[a2],clusterB[b1],clusterB[b2],&d1,&d2,&d3, APPROX);
3588
3589	/* randomize order of d2 and d3 */
3590	if (randominteger(2) == 1) {
3591	temp = d3;
3592	d3 = d2;
3593	d2 = temp;
3594	}
3595
3596	/* draw point for LM analysis */
3597	makelmpoint(trifp, d1, d2, d3);
3598	addtimes(QUARTETS, &tarr);
3599
3600	}
3601
3602	free_ivector(gettwo);
3603	}
3604
3605	if (numclust == 1) { /* normal likelihood mapping (one cluster) */
3606
3607	mlorder = new_ivector(3);
3608
3609	#if 0
3610	for (i = 3; i < Maxspc; i++)
3611	for (a = 0; a < i - 2; a++)
3612	for (b = a + 1; b < i - 1; b++)
3613	for (c = b + 1; c < i; c++)
3614	for (d = 3; d < Maxspc; d++)
3615	for (c = 2; c < d; c++)
3616	for (b = 1; b < c; b++)
3617	for (a = 0; a < b; a++)
3618	#endif
3619
3620	for (i = 3; i < Maxspc; i++)
3621	for (c = 2; c < i; c++)
3622	for (b = 1; b < c; b++)
3623	for (a = 0; a < b; a++) {
3624
3625	nq++;
3626
3627	/* check timer */
3628	checktimer(nq);
3629
3630	/* maximum likelihood values */
3631	/* approximate ML is sufficient */
3632	compute_quartlklhds(a,b,c,i,&logs[0],&logs[1],&logs[2], APPROX);
3633
3634	/* randomize order */
3635	chooser(3,3,mlorder);
3636	d1 = logs[mlorder[0]];
3637	d2 = logs[mlorder[1]];
3638	d3 = logs[mlorder[2]];
3639
3640	/* draw point for LM analysis */
3641	makelmpoint(trifp, d1, d2, d3);
3642	addtimes(QUARTETS, &tarr);
3643
3644	}
3645	free_ivector(mlorder);
3646	}
3647
3648	} else { /* randomly selected quartets */
3649
3650	if (numclust == 4) { /* four-cluster analysis */
3651
3652	for (lmqts = 0; lmqts < Numquartets; lmqts++) {
3653
3654	nq++;
3655
3656	/* check timer */
3657	checktimer(nq);
3658
3659	/* choose random quartet */
3660	qts[0] = clusterA[ randominteger(clustA) ];
3661	qts[1] = clusterB[ randominteger(clustB) ];
3662	qts[2] = clusterC[ randominteger(clustC) ];
3663	qts[3] = clusterD[ randominteger(clustD) ];
3664
3665	/* maximum likelihood values */
3666	/* approximate ML is sufficient */
3667	compute_quartlklhds(qts[0],qts[1],qts[2],qts[3],&d1,&d2,&d3, APPROX);
3668
3669	/* draw point for LM analysis */
3670	makelmpoint(trifp, d1, d2, d3);
3671	addtimes(QUARTETS, &tarr);
3672
3673	}
3674	}
3675
3676	if (numclust == 3) { /* three-cluster analysis */
3677
3678	gettwo = new_ivector(2);
3679
3680	for (lmqts = 0; lmqts < Numquartets; lmqts++) {
3681
3682	nq++;
3683
3684	/* check timer */
3685	checktimer(nq);
3686
3687	/* choose random quartet */
3688	qts[0] = clusterA[ randominteger(clustA) ];
3689	qts[1] = clusterB[ randominteger(clustB) ];
3690	chooser(clustC, 2, gettwo);
3691	qts[2] = clusterC[gettwo[0]];
3692	qts[3] = clusterC[gettwo[1]];
3693
3694	/* maximum likelihood values */
3695	/* approximate ML is sufficient */
3696	compute_quartlklhds(qts[0],qts[1],qts[2],qts[3],&d1,&d2,&d3, APPROX);
3697
3698	/* order of d2 and d3 is already randomized! */
3699
3700	/* draw point for LM analysis */
3701	makelmpoint(trifp, d1, d2, d3);
3702	addtimes(QUARTETS, &tarr);
3703
3704	}
3705
3706	free_ivector(gettwo);
3707	}
3708
3709	if (numclust == 2) { /* two-cluster analysis */
3710
3711	gettwo = new_ivector(2);
3712
3713	for (lmqts = 0; lmqts < Numquartets; lmqts++) {
3714
3715	nq++;
3716
3717	/* check timer */
3718	checktimer(nq);
3719
3720	/* choose random quartet */
3721	chooser(clustA, 2, gettwo);
3722	qts[0] = clusterA[gettwo[0]];
3723	qts[1] = clusterA[gettwo[1]];
3724	chooser(clustB, 2, gettwo);
3725	qts[2] = clusterB[gettwo[0]];
3726	qts[3] = clusterB[gettwo[1]];
3727
3728	/* maximum likelihood values */
3729	/* approximate ML is sufficient */
3730	compute_quartlklhds(qts[0],qts[1],qts[2],qts[3],&d1,&d2,&d3, APPROX);
3731
3732	/* order of d2 and d3 is already randomized! */
3733
3734	/* draw point for LM analysis */
3735	makelmpoint(trifp, d1, d2, d3);
3736	addtimes(QUARTETS, &tarr);
3737
3738	}
3739	free_ivector(gettwo);
3740	}
3741
3742	if (numclust == 1) { /* normal likelihood mapping (one cluster) */
3743
3744	for (lmqts = 0; lmqts < Numquartets; lmqts++) {
3745
3746	nq++;
3747
3748	/* check timer */
3749	checktimer(nq);
3750
3751	/* choose random quartet */
3752	chooser(Maxspc, 4, qts);
3753
3754	/* maximum likelihood values */
3755	/* approximate ML is sufficient */
3756	compute_quartlklhds(qts[0],qts[1],qts[2],qts[3],&d1,&d2,&d3, APPROX);
3757
3758	/* order of d1, d2, and d3 is already randomized! */
3759
3760	/* draw point for LM analysis */
3761	makelmpoint(trifp, d1, d2, d3);
3762	addtimes(QUARTETS, &tarr);
3763
3764	}
3765	}
3766	}
3767
3768	finishps(trifp);
3769	closefile(trifp);
3770	free_ivector(qts);
3771
3772	} /* map_lklhd */
3773
3774	/***************************************************************/
3775
3776	void setdefaults() {
3777
3778	strcpy(INFILE, INFILEDEFAULT);
3779	strcpy(OUTFILE, OUTFILEDEFAULT);
3780	strcpy(TREEFILE, TREEFILEDEFAULT);
3781	strcpy(INTREE, INTREEDEFAULT);
3782	strcpy(DISTANCES, DISTANCESDEFAULT);
3783	strcpy(TRIANGLE, TRIANGLEDEFAULT);
3784	strcpy(UNRESOLVED, UNRESOLVEDDEFAULT);
3785	strcpy(ALLQUART, ALLQUARTDEFAULT);
3786	strcpy(ALLQUARTLH, ALLQUARTLHDEFAULT);
3787	strcpy(OUTPTLIST, OUTPTLISTDEFAULT);
3788	strcpy(OUTPTORDER, OUTPTORDERDEFAULT);
3789
3790	usebestq_optn = FALSE;
3791	savequartlh_optn = FALSE;
3792	savequart_optn = FALSE;
3793	readquart_optn = FALSE;
3794
3795	randseed = -1; /* to set random random seed */
3796
3797	} /* setdefaults */
3798
3799	/***************************************************************/
3800
3801	void printversion()
3802	{
3803	# if ! PARALLEL
3804	fprintf(stderr, "puzzle (%s) %s\n", PACKAGE, VERSION);
3805	#else
3806	fprintf(stderr, "ppuzzle (%s) %s\n", PACKAGE, VERSION);
3807	# endif
3808	exit (0);
3809	}
3810	/***************************************************************/
3811
3812	void printusage(char *fname)
3813	{
3814	fprintf(stderr, "\n\nUsage: %s [-h] [ Infilename [ UserTreeFilename ] ]\n\n", fname);
3815	# if PARALLEL
3816	PP_SendDone();
3817	MPI_Finalize();
3818	# endif
3819	exit (1);
3820	}
3821
3822	/***************************************************************/
3823
3824	#ifdef HHH
3825	void printusagehhh(char *fname)
3826	{
3827	fprintf(stderr, "\n\nUsage: %s [options] [ Infilename [ UserTreeFilename ] ]\n\n", fname);
3828	fprintf(stderr, " -h - print usage\n");
3829	fprintf(stderr, " -wqf - write quartet file to Infilename.allquart\n");
3830	fprintf(stderr, " -rqf - read quartet file from Infilename.allquart\n");
3831	fprintf(stderr, " -wqlb - write quart lhs to Infilename.allquartlh (binary)\n");
3832	fprintf(stderr, " -wqla - write quart lhs to Infilename.allquartlh (ASCII)\n");
3833	fprintf(stderr, " -bestq - use best quart, no basian weights\n");
3834	fprintf(stderr, " -randseed<#> - use <#> as random number seed, for debug purposes only\n");
3835	# if PARALLEL
3836	PP_SendDone();
3837	MPI_Finalize();
3838	# endif
3839	exit (2);
3840	}
3841	#endif /* HHH */
3842
3843	/***************************************************************/
3844
3845
3846	void scancmdline(int argc, char *argv[])
3847	{
3848	static short infileset = 0;
3849	static short intreefileset = 0;
3850	short flagused;
3851	int n;
3852	int count, dummyint;
3853
3854	for (n = 1; n < *argc; n++) {
3855	# ifdef VERBOSE1
3856	printf("argv[%d] = %s\n", n, (*argv)[n]);
3857	# endif
3858
3859	flagused = FALSE;
3860
3861	# ifdef HHH
3862	dummyint = 0;
3863	count = sscanf((*argv)[n], "-wqlb%n", &dummyint);
3864	if (dummyint == 5) {
3865	savequartlh_optn = TRUE;
3866	saveqlhbin_optn = TRUE;
3867	flagused = TRUE;
3868	}
3869
3870	dummyint = 0;
3871	count = sscanf((*argv)[n], "-wqla%n", &dummyint);
3872	if (dummyint == 5) {
3873	savequartlh_optn = TRUE;
3874	saveqlhbin_optn = FALSE;
3875	flagused = TRUE;
3876	}
3877
3878	dummyint = 0;
3879	count = sscanf((*argv)[n], "-wqf%n", &dummyint);
3880	if (dummyint == 4) {
3881	savequart_optn = TRUE;
3882	flagused = TRUE;
3883	}
3884
3885	dummyint = 0;
3886	count = sscanf((*argv)[n],"-rqf%n", &dummyint);
3887	if (dummyint == 4) {
3888	readquart_optn = TRUE;
3889	flagused = TRUE;
3890	}
3891
3892	dummyint = 0;
3893	count = sscanf((*argv)[n],"-bestq%n", &dummyint);
3894	if (dummyint == 6) {
3895	usebestq_optn = TRUE;
3896	flagused = TRUE;
3897	}
3898
3899	dummyint = 0;
3900	count = sscanf((*argv)[n],"-hhh%n", &dummyint);
3901	if (dummyint==4) {
3902	printusagehhh((*argv)[0]);
3903	flagused = TRUE;
3904	}
3905	# endif /* HHH */
3906
3907	dummyint = 0;
3908	count = sscanf((*argv)[n],"-V%n", &dummyint);
3909	if (dummyint==2) {
3910	printversion((*argv)[0]);
3911	flagused = TRUE;
3912	}
3913
3914	dummyint = 0;
3915	count = sscanf((*argv)[n],"-version%n", &dummyint);
3916	if (dummyint==8) {
3917	printversion((*argv)[0]);
3918	flagused = TRUE;
3919	}
3920
3921	dummyint = 0;
3922	count = sscanf((*argv)[n],"--version%n", &dummyint);
3923	if (dummyint>=4) {
3924	printversion((*argv)[0]);
3925	flagused = TRUE;
3926	}
3927
3928	dummyint = 0;
3929	count = sscanf((*argv)[n],"-h%n", &dummyint);
3930	if (dummyint==2) {
3931	printusage((*argv)[0]);
3932	flagused = TRUE;
3933	}
3934
3935	count = sscanf((*argv)[n],"-randseed%d", &dummyint);
3936	if (count == 1) {
3937	randseed = dummyint;
3938	flagused = TRUE;
3939	}
3940
3941	#if 0
3942	count = sscanf((*argv)[n],"-h%n", &dummyint);
3943	if ((count == 1) && (dummyint>=2)) printusage((*argv)[0]);
3944
3945	count = sscanf((*argv)[n],"-writequarts%n", &dummyint);
3946	if (count == 1) writequartstofile = 1;;
3947
3948	count = sscanf((*argv)[n],"-ws%d", &dummyint);
3949	if (count == 1) windowsize = dummyint;
3950	#endif
3951
3952	if ((*argv)[n][0] != '-') {
3953	if (infileset == 0) {
3954	strcpy(INFILE, (*argv)[n]);
3955	infileset++;
3956	sprintf(OUTFILE ,"%s.%s", INFILE, OUTFILEEXT);
3957	sprintf(TREEFILE ,"%s.%s", INFILE, TREEFILEEXT);
3958	sprintf(DISTANCES ,"%s.%s", INFILE, DISTANCESEXT);
3959	sprintf(TRIANGLE ,"%s.%s", INFILE, TRIANGLEEXT);
3960	sprintf(UNRESOLVED ,"%s.%s", INFILE, UNRESOLVEDEXT);
3961	sprintf(ALLQUART ,"%s.%s", INFILE, ALLQUARTEXT);
3962	sprintf(ALLQUARTLH ,"%s.%s", INFILE, ALLQUARTLHEXT);
3963	sprintf(OUTPTLIST ,"%s.%s", INFILE, OUTPTLISTEXT);
3964	sprintf(OUTPTORDER ,"%s.%s", INFILE, OUTPTORDEREXT);
3965	FPRINTF(STDOUTFILE "Input file: %s\n", INFILE);
3966	flagused = TRUE;
3967	} else {
3968	if (intreefileset == 0) {
3969	strcpy(INTREE, (*argv)[n]);
3970	intreefileset++;
3971	sprintf(OUTFILE ,"%s.%s", INTREE, OUTFILEEXT);
3972	sprintf(TREEFILE ,"%s.%s", INTREE, TREEFILEEXT);
3973	sprintf(DISTANCES ,"%s.%s", INTREE, DISTANCESEXT);
3974	FPRINTF(STDOUTFILE "Usertree file: %s\n", INTREE);
3975	flagused = TRUE;
3976	}
3977	}
3978	}
3979	if (flagused == FALSE) {
3980	fprintf(stderr, "WARNING: commandline parameter %d not recognized (\"%s\")\n", n, (*argv)[n]);
3981	}
3982	flagused = FALSE;
3983	}
3984
3985	} /* scancmdline */
3986
3987
3988	/***************************************************************/
3989
3990	void inputandinit(int argc, char *argv[]) {
3991
3992	int ci;
3993
3994	/* vectors used in QP and LM analysis */
3995	qweight = new_dvector(3);
3996	sqdiff = new_dvector(3);
3997	qworder = new_ivector(3);
3998	sqorder = new_ivector(3);
3999
4000	/* Initialization and parsing of Commandline */
4001	setdefaults();
4002	scancmdline(argc, argv);
4003
4004	/* initialize random numbers generator */
4005	if (randseed >= 0)
4006	fprintf(stderr, "WARNING: random seed set to %d for debugging!\n", randseed);
4007	randseed = initrandom(randseed);
4008
4009	psteptreelist = NULL;
4010	psteptreesum = 0;
4011	bestratefound = 0;
4012
4013	# ifndef ALPHA
4014	FPRINTF(STDOUTFILE "\n\n\nWELCOME TO TREE-PUZZLE %s!\n\n\n", VERSION);
4015	# else
4016	FPRINTF(STDOUTFILE "\n\n\nWELCOME TO TREE-PUZZLE %s%s!\n\n\n", VERSION, ALPHA);
4017	# endif
4018
4019
4020	/* get sequences */
4021	openfiletoread(&seqfp, INFILE, "sequence data");
4022	getsizesites(seqfp);
4023	FPRINTF(STDOUTFILE "\nInput data set contains %d sequences of length %d\n", Maxspc, Maxseqc);
4024	getdataset(seqfp);
4025	closefile(seqfp);
4026	data_optn = guessdatatype();
4027
4028	/* translate characters into format used by ML engine */
4029	nuc_optn = TRUE;
4030	SH_optn = FALSE;
4031	Seqchar = NULL;
4032	translatedataset();
4033
4034	/* estimate base frequencies from data set */
4035	Freqtpm = NULL;
4036	Basecomp = NULL;
4037	estimatebasefreqs();
4038
4039	/* guess model of substitution */
4040	guessmodel();
4041
4042	/* initialize guess variables */
4043	auto_datatype = AUTO_GUESS;
4044	if (data_optn == AMINOACID) auto_aamodel = AUTO_GUESS;
4045	else auto_aamodel = AUTO_DEFAULT;
4046	/* save guessed amino acid options */
4047	guessDayhf_optn = Dayhf_optn;
4048	guessJtt_optn = Jtt_optn;
4049	guessmtrev_optn = mtrev_optn;
4050	guesscprev_optn = cprev_optn;
4051	guessblosum62_optn = blosum62_optn;
4052	guessvtmv_optn = vtmv_optn;
4053	guesswag_optn = wag_optn;
4054	guessauto_aamodel = auto_aamodel;
4055
4056
4057	/* check for user specified tree */
4058	if ((utfp = fopen(INTREE, "r")) != NULL) {
4059	fclose(utfp);
4060	puzzlemode = USERTREE;
4061	} else {
4062	puzzlemode = QUARTPUZ;
4063	}
4064
4065	/* reserve memory for cluster LM analysis */
4066	clusterA = new_ivector(Maxspc);
4067	clusterB = new_ivector(Maxspc);
4068	clusterC = new_ivector(Maxspc);
4069	clusterD = new_ivector(Maxspc);
4070
4071	/* set options interactively */
4072	setoptions();
4073
4074	/* open usertree file right after start */
4075	if (typ_optn == TREERECON_OPTN && puzzlemode == USERTREE) {
4076	openfiletoread(&utfp, INTREE, "user trees");
4077	}
4078
4079	/* start main timer */
4080	time(&Starttime);
4081	Startcpu=clock();
4082	addtimes(OPTIONS, &tarr);
4083
4084	/* symmetrize doublet frequencies if specified */
4085	symdoublets();
4086
4087	/* initialise ML */
4088	mlstart();
4089
4090	/* determine how many usertrees */
4091	if (typ_optn == TREERECON_OPTN && puzzlemode == USERTREE) {
4092	numutrees = 0;
4093	do {
4094	ci = fgetc(utfp);
4095	if ((char) ci == ';') numutrees++;
4096	} while (ci != EOF);
4097	rewind(utfp);
4098	if (numutrees < 1) {
4099	FPRINTF(STDOUTFILE "Unable to proceed (no tree in input tree file)\n\n\n");
4100	exit(1);
4101	}
4102	}
4103
4104	/* check fraction of invariable sites */
4105	if ((rhetmode == TWORATE \|\| rhetmode == MIXEDRATE) && !fracinv_optim)
4106	/* fraction of invariable site was specified manually */
4107	if (fracinv > MAXFI)
4108	fracinv = MAXFI;
4109
4110	addtimes(GENERAL, &tarr);
4111	/* estimate parameters */
4112	if (!(typ_optn == TREERECON_OPTN && puzzlemode == USERTREE)) {
4113	/* no tree present */
4114	estimateparametersnotree();
4115	} else {
4116	if (utree_optn) {
4117	/* use 1st user tree */
4118	readusertree(utfp);
4119	rewind(utfp);
4120	estimateparameterstree();
4121	} else {
4122	/* don't use first user tree */
4123	estimateparametersnotree();
4124	}
4125	}
4126	addtimes(PARAMEST, &tarr);
4127
4128	/* compute expected Ts/Tv ratio */
4129	if (data_optn == NUCLEOTIDE) computeexpectations();
4130
4131	} /* inputandinit */
4132
4133
4134
4135	/***************************************************************/
4136
4137	void evaluatetree(FILE intreefp, FILE outtreefp, int pmode, int utreenum, int maxutree, int *oldlocroot)
4138	{
4139
4140	switch (pmode) {
4141	case QUARTPUZ: /* read QP tree */
4142	readusertree(intreefp);
4143	FPRINTF(STDOUTFILE "Computing maximum likelihood branch lengths (without clock)\n");
4144	fflush(STDOUT);
4145	usertree_lklhd();
4146	findbestratecombination();
4147	break;
4148	case USERTREE: /* read user tree */
4149	readusertree(intreefp);
4150	FPRINTF(STDOUTFILE "Computing maximum likelihood branch lengths (without clock) for tree # %d\n", utreenum+1);
4151	fflush(STDOUT);
4152	usertree_lklhd();
4153	if (maxutree > 1) {
4154	ulkl[utreenum] = Ctree->lklhd;
4155	allsitelkl(Ctree->condlkl, allsites[utreenum]);
4156	}
4157	if (utreenum==0) findbestratecombination();
4158	break;
4159	}
4160
4161
4162	if (compclock) { /* clocklike branch length */
4163	switch (pmode) {
4164	case QUARTPUZ:
4165	FPRINTF(STDOUTFILE "Computing maximum likelihood branch lengths (with clock)\n");
4166	fflush(STDOUT);
4167	break;
4168	case USERTREE:
4169	FPRINTF(STDOUTFILE "Computing maximum likelihood branch lengths (with clock) for tree # %d\n", utreenum+1);
4170	fflush(STDOUT);
4171	break;
4172	}
4173
4174	/* find best place for root */
4175	rootsearch = 0;
4176
4177	if (utreenum==0) locroot = *oldlocroot;
4178	else *oldlocroot = locroot;
4179
4180	if (locroot < 0) {
4181	locroot = findrootedge();
4182	rootsearch = 1;
4183	}
4184	/* if user-specified edge for root does not exist use displayed outgroup */
4185	if (!checkedge(locroot)) {
4186	locroot = outgroup;
4187	rootsearch = 2;
4188	}
4189	/* compute likelihood */
4190	clock_lklhd(locroot);
4191	if (maxutree > 1) {
4192	ulklc[utreenum] = Ctree->lklhdc;
4193	allsitelkl(Ctree->condlkl, allsitesc[utreenum]);
4194	}
4195
4196	}
4197
4198	if (clockmode == 0)
4199	fprintf(outtreefp, "[ lh=%.6f ]", Ctree->lklhd);
4200	else
4201	fprintf(outtreefp, "[ lh=%.6f ]", Ctree->lklhdc);
4202
4203	/* write ML branch length tree to outree file */
4204	clockmode = 0; /* nonclocklike branch lengths */
4205	fputphylogeny(outtreefp);
4206
4207	/* clocklike branch lengths */
4208	if (compclock) {
4209	clockmode = 1;
4210	fputrooted(outtreefp, locroot);
4211	}
4212	} /* evaluatetree */
4213
4214	/***************************************************************/
4215
4216	void memcleanup() {
4217	if (puzzlemode == QUARTPUZ && typ_optn == TREERECON_OPTN) {
4218	free(splitfreqs);
4219	free(splitpatterns);
4220	free(splitsizes);
4221	free_ivector(consconfid);
4222	free_ivector(conssizes);
4223	free_cmatrix(consbiparts);
4224	free_ulivector(badtaxon);
4225	}
4226	free_cmatrix(Identif);
4227	free_dvector(Freqtpm);
4228	free_imatrix(Basecomp);
4229	free_ivector(clusterA);
4230	free_ivector(clusterB);
4231	free_ivector(clusterC);
4232	free_ivector(clusterD);
4233	free_dvector(qweight);
4234	free_dvector(sqdiff);
4235	free_ivector(qworder);
4236	free_ivector(sqorder);
4237	freetreelist(&psteptreelist, &psteptreenum, &psteptreesum);
4238	} /* memcleanup */
4239
4240	/***************************************************************/
4241
4242
4243	/******************************************************************************/
4244	/* main part */
4245	/******************************************************************************/
4246
4247	int main(int argc, char *argv[])
4248	{
4249	int i, oldlocroot=0;
4250
4251	/* start main timer */
4252	time(&walltimestart);
4253	cputimestart = clock();
4254	inittimearr(&tarr);
4255
4256	# if PARALLEL
4257	PP_Init(&argc, &argv);
4258	if (PP_IamSlave) {
4259	slave_main(argc, argv);
4260	} else {
4261	# endif /* PARALLEL */
4262
4263	inputandinit(&argc, &argv);
4264
4265	FPRINTF(STDOUTFILE "Writing parameters to file %s\n", OUTFILE);
4266	openfiletowrite(&ofp, OUTFILE, "general output");
4267	/* openfiletowrite(&ofp, "xxxx", "general output"); */
4268	writeoutputfile(ofp,WRITEPARAMS);
4269	fclose(ofp);
4270
4271
4272	/* write distance matrix */
4273	FPRINTF(STDOUTFILE "Writing pairwise distances to file %s\n", DISTANCES);
4274	openfiletowrite(&dfp, DISTANCES, "pairwise distances");
4275	putdistance(dfp);
4276	closefile(dfp);
4277
4278	# if PARALLEL
4279	PP_SendSizes(Maxspc, Maxsite, numcats, Numptrn, tpmradix, outgroup, fracconst, randseed);
4280	PP_SendData(Seqpat, /* cmatrix */
4281	Alias, Weight, constpat, /* ivector */
4282	Rates, Eval, Freqtpm, /* dvector */
4283	Evec, Ievc, iexp, Distanmat, /* dmatrix */
4284	ltprobr); /* dcube */
4285	# endif /* PARALLEL */
4286	psteptreestrlen = (Maxspc * (int)(1 + log10(Maxspc))) +
4287	(Maxspc * 3);
4288
4289	switch (typ_optn) {
4290	case TREERECON_OPTN: /* tree reconstruction */
4291
4292	if (puzzlemode == QUARTPUZ) { /* quartet puzzling */
4293	recon_tree();
4294	} /* quartet puzzling */
4295	break;
4296
4297	case LIKMAPING_OPTN: /* likelihood mapping */
4298
4299	map_lklhd();
4300	break;
4301	} /* switch typ_optn */
4302
4303
4304	free_cmatrix(Seqchar);
4305	free_cmatrix(seqchars);
4306
4307	/* reserve memory for tree statistics */
4308	if (typ_optn == TREERECON_OPTN && puzzlemode == USERTREE && numutrees > 1) {
4309	ulkl = new_dvector(numutrees);
4310	allsites = new_dmatrix(numutrees,Numptrn);
4311	if (compclock) {
4312	ulklc = new_dvector(numutrees);
4313	allsitesc = new_dmatrix(numutrees,Numptrn);
4314	}
4315	}
4316
4317	/* write puzzling step tree list */
4318	if ((listqptrees == PSTOUT_ORDER) \|\| (listqptrees == PSTOUT_LISTORDER)) {
4319	openfiletowrite(&qptorder, OUTPTORDER, "puzzling step trees (unique)");
4320
4321	fprintfsortedpstrees(qptorder, psteptreelist, psteptreenum, psteptreesum, 1, 0.0);
4322	closefile(qptorder);
4323	}
4324
4325	/* compute ML branch lengths for QP tree and for 1st user tree */
4326	switch(typ_optn) {
4327	case TREERECON_OPTN:
4328
4329	/* open outtree file */
4330	openfiletowrite(&tfp, TREEFILE, "output tree(s)");
4331
4332	addtimes(GENERAL, &tarr);
4333
4334	switch (puzzlemode) {
4335	case QUARTPUZ: /* read QP tree */
4336	rewind(tmpfp);
4337	openfiletowrite(&tfp, TREEFILE, "output tree(s)");
4338	evaluatetree(tmpfp, tfp, puzzlemode, 0, 1, &oldlocroot);
4339	addtimes(TREEEVAL, &tarr);
4340	closefile(tmpfp);
4341	closefile(tfp);
4342
4343	openfiletoappend(&ofp, OUTFILE, "general output");
4344	writeoutputfile(ofp,WRITEREST);
4345	break;
4346	case USERTREE: /* read user tree */
4347	openfiletoappend(&ofp, OUTFILE, "general output");
4348
4349	openfiletowrite(&tfp, TREEFILE, "output tree(s)");
4350	for (i = 0; i < numutrees; i++) {
4351	evaluatetree(utfp, tfp, puzzlemode, i, numutrees, &oldlocroot);
4352	if (i==0) writeoutputfile(ofp,WRITEREST);
4353	writecutree(ofp, i+1);
4354	addtimes(TREEEVAL, &tarr);
4355	}
4356	closefile(tfp);
4357	closefile(utfp);
4358	break;
4359	default:
4360	openfiletoappend(&ofp, OUTFILE, "general output");
4361	writeoutputfile(ofp,WRITEREST);
4362	break;
4363	} /* switch puzzlemode */
4364	break;
4365	default:
4366	openfiletoappend(&ofp, OUTFILE, "general output");
4367	writeoutputfile(ofp,WRITEREST);
4368	break;
4369	} /* switch typ_optn */
4370
4371	/* print tree statistics */
4372	if (typ_optn == TREERECON_OPTN && puzzlemode == USERTREE && numutrees > 1)
4373	printtreestats(ofp);
4374
4375	/* free memory for tree statistics */
4376	if (typ_optn == TREERECON_OPTN && puzzlemode == USERTREE && numutrees > 1) {
4377	free_dvector(ulkl);
4378	free_dmatrix(allsites);
4379	if (compclock) {
4380	free_dvector(ulklc);
4381	free_dmatrix(allsitesc);
4382	}
4383	}
4384
4385	# if PARALLEL
4386	PP_SendDone();
4387	# endif /* PARALLEL */
4388
4389	/* write CPU/Wallclock times and parallel statistics */
4390	time(&walltimestop);
4391	cputimestop = clock();
4392	addtimes(OVERALL, &tarr);
4393	# ifdef TIMEDEBUG
4394	printtimearr(&tarr);
4395	# endif /* TIMEDEBUG */
4396	fullcpu = tarr.fullcpu;
4397	fulltime = tarr.fulltime;
4398
4399	# if PARALLEL
4400	writetimesstat(ofp);
4401	# endif /* PARALLEL */
4402
4403	/* stop timer */
4404	time(&Stoptime);
4405	Stopcpu=clock();
4406	timestamp(ofp);
4407	closefile(ofp);
4408
4409
4410	/* printbestratecombination(stderr); */
4411	mlfinish();
4412
4413	FPRINTF(STDOUTFILE "\nAll results written to disk:\n");
4414	FPRINTF(STDOUTFILE " Puzzle report file: %s\n", OUTFILE);
4415	FPRINTF(STDOUTFILE " Likelihood distances: %s\n", DISTANCES);
4416
4417	if (typ_optn == TREERECON_OPTN && puzzlemode != PAIRDIST)
4418	FPRINTF(STDOUTFILE " Phylip tree file: %s\n", TREEFILE);
4419	if (typ_optn == TREERECON_OPTN && puzzlemode == QUARTPUZ) {
4420	if ((listqptrees == PSTOUT_ORDER) \|\|(listqptrees == PSTOUT_LISTORDER))
4421	FPRINTF(STDOUTFILE " Unique puzzling step trees: %s\n", OUTPTORDER);
4422	if ((listqptrees == PSTOUT_LIST) \|\|(listqptrees == PSTOUT_LISTORDER))
4423	FPRINTF(STDOUTFILE " Puzzling step tree list: %s\n", OUTPTLIST);
4424	}
4425	if (show_optn && typ_optn == TREERECON_OPTN && puzzlemode == QUARTPUZ)
4426	FPRINTF(STDOUTFILE " Unresolved quartets: %s\n", UNRESOLVED);
4427	if (typ_optn == LIKMAPING_OPTN)
4428	FPRINTF(STDOUTFILE " Likelihood mapping diagram: %s\n", TRIANGLE);
4429	FPRINTF(STDOUTFILE "\n");
4430
4431	/* runtime message */
4432	FPRINTF(STDOUTFILE
4433	"The computation took %.0f seconds (= %.1f minutes = %.1f hours)\n",
4434	difftime(Stoptime, Starttime), difftime(Stoptime, Starttime)/60.,
4435	difftime(Stoptime, Starttime)/3600.);
4436	FPRINTF(STDOUTFILE
4437	" including input %.0f seconds (= %.1f minutes = %.1f hours)\n",
4438	fulltime, fulltime/60., fulltime/3600.);
4439	#ifdef TIMEDEBUG
4440	FPRINTF(STDOUTFILE
4441	"and %.0f seconds CPU time (= %.1f minutes = %.1f hours)\n\n",
4442	fullcpu, fullcpu/60., fullcpu/3600.);
4443	#endif /* TIMEDEBUG */
4444
4445	/* free memory */
4446	memcleanup();
4447
4448	# if PARALLEL
4449	} /* !IamSlave */
4450	PP_Finalize();
4451	# endif /* PARALLEL */
4452
4453	return 0;
4454	}
4455
4456
4457	/* compare function for uli - sort largest numbers first */
4458	int ulicmp(const void ap, const void bp)
4459	{
4460	uli a, b;
4461
4462	a = ((uli ) ap);
4463	b = ((uli ) bp);
4464
4465	if (a > b) return -1;
4466	else if (a < b) return 1;
4467	else return 0;
4468	}
4469
4470	/* compare function for int - sort smallest numbers first */
4471	int intcmp(const void ap, const void bp)
4472	{
4473	int a, b;
4474
4475	a = ((int ) ap);
4476	b = ((int ) bp);
4477
4478	if (a < b) return -1;
4479	else if (a > b) return 1;
4480	else return 0;
4481	}

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

Download in other formats: