Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

puzzle1.c

Visit:

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

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

Context Navigation

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

Download in other formats: