Context Navigation

source: branches/stable/GDE/TREEPUZZLE/src/puzzle2.c

Visit:

Last change on this file was 5782, checked in by westram, 16 years ago
fixed spelling of 'occurrence'
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 68.4 KB

Line
1	/*
2	* puzzle2.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 extern
17
18	#include "puzzle.h"
19	#include <string.h>
20
21	#if PARALLEL
22	# include "sched.h"
23	#endif /* PARALLEL */
24
25
26	/******************************************************************************/
27	/* sequences */
28	/******************************************************************************/
29
30	/* read ten characters of current line as identifier */
31	void readid(FILE *infp, int t)
32	{
33	int i, j, flag, ci;
34
35	for (i = 0; i < 10; i++) {
36	ci = fgetc(infp);
37	if (ci == EOF \|\| !isprint(ci)) {
38	FPRINTF(STDOUTFILE "\n\n\nUnable to proceed (no name for sequence %d)\n\n\n", t+1);
39	exit(1);
40	}
41	Identif[t][i] = (char) ci;
42	}
43	/* convert leading blanks in taxon name to underscores */
44	flag = FALSE;
45	for (i = 9; i > -1; i--) {
46	if (flag == FALSE) {
47	if (Identif[t][i] != ' ') flag = TRUE;
48	} else {
49	if (Identif[t][i] == ' ') Identif[t][i] = '_';
50	}
51	}
52	/* check whether this name is already used */
53	for (i = 0; i < t; i++) { /* compare with all other taxa */
54	flag = TRUE; /* assume identity */
55	for (j = 0; (j < 10) && (flag == TRUE); j++)
56	if (Identif[t][j] != Identif[i][j])
57	flag = FALSE;
58	if (flag) {
59	FPRINTF(STDOUTFILE "\n\n\nUnable to proceed (multiple occurrences of sequence name '");
60	fputid(STDOUT, t);
61	FPRINTF(STDOUTFILE "')\n\n\n");
62	exit(1);
63	}
64	}
65	}
66
67	/* read next allowed character */
68	char readnextcharacter(FILE *ifp, int notu, int nsite)
69	{
70	char c;
71
72	/* ignore blanks and control characters except newline */
73	do {
74	if (fscanf(ifp, "%c", &c) != 1) {
75	FPRINTF(STDOUTFILE "\n\n\nUnable to proceed (missing character at position %d in sequence '", nsite + 1);
76	fputid(STDOUT, notu);
77	FPRINTF(STDOUTFILE "')\n\n\n");
78	exit(1);
79	}
80	} while (c == ' ' \|\| (iscntrl((int) c) && c != '\n'));
81	return c;
82	}
83
84	/* skip rest of the line */
85	void skiprestofline(FILE* ifp, int notu, int nsite)
86	{
87	int ci;
88
89	/* read chars until the first newline */
90	do{
91	ci = fgetc(ifp);
92	if (ci == EOF) {
93	FPRINTF(STDOUTFILE "Unable to proceed (missing newline at position %d in sequence '", nsite + 1);
94	fputid(STDOUT, notu);
95	FPRINTF(STDOUTFILE "')\n\n\n");
96	exit(1);
97	}
98	} while ((char) ci != '\n');
99	}
100
101	/* skip control characters and blanks */
102	void skipcntrl(FILE *ifp, int notu, int nsite)
103	{
104	int ci;
105
106	/* read over all control characters and blanks */
107	do {
108	ci = fgetc(ifp);
109	if (ci == EOF) {
110	FPRINTF(STDOUTFILE "\n\n\nUnable to proceed (missing character at position %d in sequence '", nsite + 1);
111	fputid(STDOUT, notu);
112	FPRINTF(STDOUTFILE "')\n\n\n");
113	exit(1);
114	}
115	} while (iscntrl(ci) \|\| (char) ci == ' ');
116	/* go one character back */
117	if (ungetc(ci, ifp) == EOF) {
118	FPRINTF(STDOUTFILE "\n\n\nUnable to proceed (positioning error at position %d in sequence '", nsite + 1);
119	fputid(STDOUT, notu);
120	FPRINTF(STDOUTFILE "')\n\n\n");
121	exit(1);
122	}
123	}
124
125	/* read sequences of one data set */
126	void getseqs(FILE *ifp)
127	{
128	int notu, nsite, endofline, linelength, i;
129	char c;
130
131	seqchars = new_cmatrix(Maxspc, Maxseqc);
132	/* read all characters */
133	nsite = 0; /* next site to be read */
134	while (nsite < Maxseqc) {
135	/* read first taxon */
136	notu = 0;
137	/* go to next true line */
138	skiprestofline(ifp, notu, nsite);
139	skipcntrl(ifp, notu, nsite);
140	if (nsite == 0) readid(ifp, notu);
141	endofline = FALSE;
142	linelength = 0;
143	do {
144	c = readnextcharacter(ifp, notu, nsite + linelength);
145	if (c == '\n') endofline = TRUE;
146	else if (c == '.') {
147	FPRINTF(STDOUTFILE "\n\n\nUnable to proceed (invalid character '.' at position ");
148	FPRINTF(STDOUTFILE "%d in first sequence)\n\n\n", nsite + linelength + 1);
149	exit(1);
150	} else if (nsite + linelength < Maxseqc) {
151	/* change to upper case */
152	seqchars[notu][nsite + linelength] = (char) toupper((int) c);
153	linelength++;
154	} else {
155	endofline = TRUE;
156	skiprestofline(ifp, notu, nsite + linelength);
157	}
158	} while (!endofline);
159	if (linelength == 0) {
160	FPRINTF(STDOUTFILE "\n\n\nUnable to proceed (line with length 0 at position %d in sequence '", nsite + 1);
161	fputid(STDOUT, notu);
162	FPRINTF(STDOUTFILE "')\n\n\n");
163	exit(1);
164	}
165	/* read other taxa */
166	for (notu = 1; notu < Maxspc; notu++) {
167	/* go to next true line */
168	if (notu != 1) skiprestofline(ifp, notu, nsite);
169	skipcntrl(ifp, notu, nsite);
170	if (nsite == 0) readid(ifp, notu);
171	for (i = nsite; i < nsite + linelength; i++) {
172	c = readnextcharacter(ifp, notu, i);
173	if (c == '\n') { /* too short */
174	FPRINTF(STDOUTFILE "\n\n\nUnable to proceed (line to short at position %d in sequence '", i + 1);
175	fputid(STDOUT, notu);
176	FPRINTF(STDOUTFILE "')\n\n\n");
177	exit(1);
178	} else if (c == '.') {
179	seqchars[notu][i] = seqchars[0][i];
180	} else {
181	/* change to upper case */
182	seqchars[notu][i] = (char) toupper((int) c);
183	}
184	}
185	}
186	nsite = nsite + linelength;
187	}
188	}
189
190	/* initialize identifer array */
191	void initid(int t)
192	{
193	int i, j;
194
195	Identif = new_cmatrix(t, 10);
196	for (i = 0; i < t; i++)
197	for (j = 0; j < 10; j++)
198	Identif[i][j] = ' ';
199	}
200
201	/* print identifier of specified taxon in full 10 char length */
202	void fputid10(FILE *ofp, int t)
203	{
204	int i;
205
206	for (i = 0; i < 10; i++) fputc(Identif[t][i], ofp);
207	}
208
209	/* print identifier of specified taxon up to first space */
210	int fputid(FILE *ofp, int t)
211	{
212	int i;
213
214	i = 0;
215	while (Identif[t][i] != ' ' && i < 10) {
216	fputc(Identif[t][i], ofp);
217	i++;
218	}
219	return i;
220	}
221
222	/* read first line of sequence data set */
223	void getsizesites(FILE *ifp)
224	{
225	if (fscanf(ifp, "%d", &Maxspc) != 1) {
226	FPRINTF(STDOUTFILE "\n\n\nUnable to proceed (missing number of sequences)\n\n\n");
227	exit(1);
228	}
229	if (fscanf(ifp, "%d", &Maxseqc) != 1) {
230	FPRINTF(STDOUTFILE "\n\n\nUnable to proceed (missing number of sites)\n\n\n");
231	exit(1);
232	}
233
234	if (Maxspc < 4) {
235	FPRINTF(STDOUTFILE "\n\n\nUnable to proceed (less than 4 sequences)\n\n\n");
236	exit(1);
237	}
238	if (Maxspc > 257) {
239	FPRINTF(STDOUTFILE "\n\n\nUnable to proceed (more than 257 sequences)\n\n\n");
240	exit(1);
241	}
242	if (Maxseqc < 1) {
243	FPRINTF(STDOUTFILE "\n\n\nUnable to proceed (no sequence sites)\n\n\n");
244	exit(1);
245	}
246	Maxbrnch = 2*Maxspc - 3;
247	}
248
249	/* read one data set - PHYLIP interleaved */
250	void getdataset(FILE *ifp)
251	{
252	initid(Maxspc);
253	getseqs(ifp);
254	}
255
256	/* guess data type */
257	int guessdatatype()
258	{
259	uli numnucs, numchars, numbins;
260	int notu, nsite;
261	char c;
262
263	/* count A, C, G, T, U, N */
264	numnucs = 0;
265	numchars = 0;
266	numbins = 0;
267	for (notu = 0; notu < Maxspc; notu++)
268	for (nsite = 0; nsite < Maxseqc; nsite++) {
269	c = seqchars[notu][nsite];
270	if (c == 'A' \|\| c == 'C' \|\| c == 'G' \|\|
271	c == 'T' \|\| c == 'U' \|\| c == 'N') numnucs++;
272	if (c != '-' && c != '?') numchars++;
273	if (c == '0' \|\| c == '1') numbins++;
274	}
275	if (numchars == 0) numchars = 1;
276	/* more than 85 % frequency means nucleotide data */
277	if ((double) numnucs / (double) numchars > 0.85) return 0;
278	else if ((double) numbins / (double) numchars > 0.2) return 2;
279	else return 1;
280	}
281
282	/* translate characters into format used by ML engine */
283	void translatedataset()
284	{
285	int notu, sn, co;
286	char c;
287	cvector code;
288
289
290	/* determine Maxsite - number of ML sites per taxon */
291	if (data_optn == 0 && SH_optn) {
292	if (SHcodon)
293	Maxsite = Maxseqc / 3;
294	else
295	Maxsite = Maxseqc / 2; /* assume doublets */
296
297	} else
298	Maxsite = Maxseqc;
299	if (data_optn == 0 && (Maxsite % 3) == 0 && !SH_optn) {
300	if (codon_optn == 1 \|\| codon_optn == 2 \|\| codon_optn == 3)
301	Maxsite = Maxsite / 3; /* only one of the three codon positions */
302	if (codon_optn == 4)
303	Maxsite = 2(Maxsite / 3); / 1st + 2nd codon positions */
304	}
305
306	/* reserve memory */
307	if (Seqchar != NULL) free_cmatrix(Seqchar);
308	Seqchar = new_cmatrix(Maxspc, Maxsite);
309
310	/* code length */
311	if (data_optn == 0 && SH_optn)
312	code = new_cvector(2);
313	else
314	code = new_cvector(1);
315
316	/* decode characters */
317	if (data_optn == 0 && SH_optn) { /* SH doublets */
318
319	for (notu = 0; notu < Maxspc; notu++) {
320	for (sn = 0; sn < Maxsite; sn++) {
321	for (co = 0; co < 2; co++) {
322	if (SHcodon)
323	c = seqchars[notu][sn*3 + co];
324	else
325	c = seqchars[notu][sn*2 + co];
326	code[co] = c;
327	}
328	Seqchar[notu][sn] = code2int(code);
329	}
330	}
331
332	} else if (!(data_optn == 0 && (Maxseqc % 3) == 0)) { /* use all */
333
334	for (notu = 0; notu < Maxspc; notu++) {
335	for (sn = 0; sn < Maxsite; sn++) {
336	code[0] = seqchars[notu][sn];
337	Seqchar[notu][sn] = code2int(code);
338	}
339	}
340
341	} else { /* codons */
342
343	for (notu = 0; notu < Maxspc; notu++) {
344	for (sn = 0; sn < Maxsite; sn++) {
345	if (codon_optn == 1 \|\| codon_optn == 2 \|\| codon_optn == 3)
346	code[0] = seqchars[notu][sn*3+codon_optn-1];
347	else if (codon_optn == 4) {
348	if ((sn % 2) == 0)
349	code[0] = seqchars[notu][(sn/2)*3];
350	else
351	code[0] = seqchars[notu][((sn-1)/2)*3+1];
352	} else
353	code[0] = seqchars[notu][sn];
354	Seqchar[notu][sn] = code2int(code);
355	}
356	}
357
358	}
359	free_cvector(code);
360	}
361
362	/* estimate mean base frequencies from translated data set */
363	void estimatebasefreqs()
364	{
365	int tpmradix, i, j;
366	uli all, *gene;
367
368	tpmradix = gettpmradix();
369
370	if (Freqtpm != NULL) free_dvector(Freqtpm);
371	Freqtpm = new_dvector(tpmradix);
372
373	if (Basecomp != NULL) free_imatrix(Basecomp);
374	Basecomp = new_imatrix(Maxspc, tpmradix);
375
376	gene = (uli ) malloc((unsigned) ((tpmradix + 1) sizeof(uli)));
377	if (gene == NULL) maerror("gene in estimatebasefreqs");
378
379	for (i = 0; i < tpmradix + 1; i++) gene[i] = 0;
380	for (i = 0; i < Maxspc; i++)
381	for (j = 0; j < tpmradix; j++) Basecomp[i][j] = 0;
382	for (i = 0; i < Maxspc; i++)
383	for (j = 0; j < Maxsite; j++) {
384	gene[(int) Seqchar[i][j]]++;
385	if (Seqchar[i][j] != tpmradix) Basecomp[i][(int) Seqchar[i][j]]++;
386	}
387
388	all = Maxspc * Maxsite - gene[tpmradix];
389	if (all != 0) { /* normal case */
390	for (i = 0; i < tpmradix; i++)
391	Freqtpm[i] = (double) gene[i] / (double) all;
392	} else { /* pathological case with no unique character in data set */
393	for (i = 0; i < tpmradix; i++)
394	Freqtpm[i] = 1.0 / (double) tpmradix;
395	}
396
397	free(gene);
398
399	Frequ_optn = TRUE;
400	}
401
402	/* guess model of substitution */
403	void guessmodel()
404	{
405	double c1, c2, c3, c4, c5, c6;
406	dvector f;
407	dmatrix a;
408	int i;
409
410	Dayhf_optn = FALSE;
411	Jtt_optn = TRUE;
412	mtrev_optn = FALSE;
413	cprev_optn = FALSE;
414	blosum62_optn = FALSE;
415	vtmv_optn = FALSE;
416	wag_optn = FALSE;
417	TSparam = 2.0;
418	YRparam = 1.0;
419	optim_optn = TRUE;
420	HKY_optn = TRUE;
421	TN_optn = FALSE;
422
423	if (data_optn == 1) { /* amino acids */
424
425	/* chi2 fit to amino acid frequencies */
426
427	f = new_dvector(20);
428	a = new_dmatrix(20,20);
429	/* chi2 distance Dayhoff */
430	dyhfdata(a, f);
431	c1 = 0;
432	for (i = 0; i < 20; i++)
433	c1 = c1 + (Freqtpm[i]-f[i])*(Freqtpm[i]-f[i]);
434	/* chi2 distance JTT */
435	jttdata(a, f);
436	c2 = 0;
437	for (i = 0; i < 20; i++)
438	c2 = c2 + (Freqtpm[i]-f[i])*(Freqtpm[i]-f[i]);
439	/* chi2 distance mtREV */
440	mtrevdata(a, f);
441	c3 = 0;
442	for (i = 0; i < 20; i++)
443	c3 = c3 + (Freqtpm[i]-f[i])*(Freqtpm[i]-f[i]);
444	/* chi2 distance VT */
445	vtmvdata(a, f);
446	c4 = 0;
447	for (i = 0; i < 20; i++)
448	c4 = c4 + (Freqtpm[i]-f[i])*(Freqtpm[i]-f[i]);
449	/* chi2 distance WAG */
450	wagdata(a, f);
451	c5 = 0;
452	for (i = 0; i < 20; i++)
453	c5 = c5 + (Freqtpm[i]-f[i])*(Freqtpm[i]-f[i]);
454	/* chi2 distance cpREV */
455	cprev45data(a, f);
456	c6 = 0;
457	for (i = 0; i < 20; i++)
458	c6 = c6 + (Freqtpm[i]-f[i])*(Freqtpm[i]-f[i]);
459
460	free_dvector(f);
461	free_dmatrix(a);
462
463	#ifndef CPREV
464	if ((c1 < c2) && (c1 < c3) && (c1 < c4) && (c1 < c5)) {
465	/* c1 -> Dayhoff */
466	Dayhf_optn = TRUE;
467	Jtt_optn = FALSE;
468	mtrev_optn = FALSE;
469	cprev_optn = FALSE;
470	vtmv_optn = FALSE;
471	wag_optn = FALSE;
472	FPRINTF(STDOUTFILE "(consists very likely of amino acids encoded on nuclear DNA)\n");
473	} else {
474	if ((c2 < c3) && (c2 < c4) && (c2 < c5)) {
475	/* c2 -> JTT */
476	Dayhf_optn = FALSE;
477	Jtt_optn = TRUE;
478	mtrev_optn = FALSE;
479	cprev_optn = FALSE;
480	vtmv_optn = FALSE;
481	wag_optn = FALSE;
482	FPRINTF(STDOUTFILE "(consists very likely of amino acids encoded on nuclear DNA)\n");
483	} else {
484	if ((c3 < c4) && (c3 < c5)) {
485	/* c3 -> mtREV */
486	Dayhf_optn = FALSE;
487	Jtt_optn = FALSE;
488	mtrev_optn = TRUE;
489	cprev_optn = FALSE;
490	vtmv_optn = FALSE;
491	wag_optn = FALSE;
492	FPRINTF(STDOUTFILE "(consists very likely of amino acids encoded on mtDNA)\n");
493	} else {
494	if ((c4 < c5)) {
495	/* c4 -> VT */
496	Dayhf_optn = FALSE;
497	Jtt_optn = FALSE;
498	mtrev_optn = FALSE;
499	cprev_optn = FALSE;
500	vtmv_optn = TRUE;
501	wag_optn = FALSE;
502	FPRINTF(STDOUTFILE "(consists very likely of amino acids encoded on nuclear DNA)\n");
503	} else {
504	/* c5 -> WAG */
505	Dayhf_optn = FALSE;
506	Jtt_optn = FALSE;
507	mtrev_optn = FALSE;
508	cprev_optn = FALSE;
509	vtmv_optn = FALSE;
510	wag_optn = TRUE;
511	FPRINTF(STDOUTFILE "(consists very likely of amino acids encoded on nuclear DNA)\n");
512	} /* if c4 else c5 */
513	} /* if c3 else c4 */
514	} /* if c2 */
515	} /* if c1 */
516
517	#else /* CPREV */
518
519	if ((c1 < c2) && (c1 < c3) && (c1 < c4) && (c1 < c5) && (c1 < c6)) {
520	/* c1 -> Dayhoff */
521	Dayhf_optn = TRUE;
522	Jtt_optn = FALSE;
523	mtrev_optn = FALSE;
524	cprev_optn = FALSE;
525	vtmv_optn = FALSE;
526	wag_optn = FALSE;
527	FPRINTF(STDOUTFILE "(consists very likely of amino acids encoded on nuclear DNA)\n");
528	} else {
529	if ((c2 < c3) && (c2 < c4) && (c2 < c5) && (c2 < c6)) {
530	/* c2 -> JTT */
531	Dayhf_optn = FALSE;
532	Jtt_optn = TRUE;
533	mtrev_optn = FALSE;
534	cprev_optn = FALSE;
535	vtmv_optn = FALSE;
536	wag_optn = FALSE;
537	FPRINTF(STDOUTFILE "(consists very likely of amino acids encoded on nuclear DNA)\n");
538	} else {
539	if ((c3 < c4) && (c3 < c5) && (c3 < c6)) {
540	/* c3 -> mtREV */
541	Dayhf_optn = FALSE;
542	Jtt_optn = FALSE;
543	mtrev_optn = TRUE;
544	cprev_optn = FALSE;
545	vtmv_optn = FALSE;
546	wag_optn = FALSE;
547	FPRINTF(STDOUTFILE "(consists very likely of amino acids encoded on mtDNA)\n");
548	} else {
549	if ((c4 < c5) && (c4 < c6)) {
550	/* c4 -> VT */
551	Dayhf_optn = FALSE;
552	Jtt_optn = FALSE;
553	mtrev_optn = FALSE;
554	cprev_optn = FALSE;
555	vtmv_optn = TRUE;
556	wag_optn = FALSE;
557	FPRINTF(STDOUTFILE "(consists very likely of amino acids encoded on nuclear DNA)\n");
558	} else {
559	if (c5 < c6) {
560	/* c5 -> WAG */
561	Dayhf_optn = FALSE;
562	Jtt_optn = FALSE;
563	mtrev_optn = FALSE;
564	cprev_optn = FALSE;
565	vtmv_optn = FALSE;
566	wag_optn = TRUE;
567	FPRINTF(STDOUTFILE "(consists very likely of amino acids encoded on nuclear DNA)\n");
568	} else {
569	/* if (c6) */
570	/* c6 -> cpREV */
571	Dayhf_optn = FALSE;
572	Jtt_optn = FALSE;
573	mtrev_optn = FALSE;
574	cprev_optn = TRUE;
575	vtmv_optn = FALSE;
576	wag_optn = FALSE;
577	FPRINTF(STDOUTFILE "(consists very likely of amino acids encoded on cpDNA)\n");
578	} /* if c5 else c6 */
579	} /* if c4 else c5 */
580	} /* if c3 else c4 */
581	} /* if c2 */
582	} /* if c1 */
583	#endif /* CPREV */
584
585	} else if (data_optn == 0) {
586	FPRINTF(STDOUTFILE "(consists very likely of nucleotides)\n");
587	} else {
588	FPRINTF(STDOUTFILE "(consists very likely of binary state data)\n");
589	}
590	} /* guessmodel */
591
592
593	/******************************************************************************/
594	/* functions for representing and building puzzling step trees */
595	/******************************************************************************/
596
597	/* initialize tree with the following starting configuration
598
599	2
600	0 +------- C(=2)
601	A(=0) -----+
602	+------- B(=1)
603	1
604	*/
605	void inittree()
606	{
607	int i;
608
609	/* allocate the memory for the whole tree */
610
611	/* allocate memory for vector with all the edges of the tree */
612	edge = (ONEEDGE *) calloc(Maxbrnch, sizeof(ONEEDGE) );
613	if (edge == NULL) maerror("edge in inittree");
614
615	/* allocate memory for vector with edge numbers of leaves */
616	edgeofleaf = (int *) calloc(Maxspc, sizeof(int) );
617	if (edgeofleaf == NULL) maerror("edgeofleaf in inittree");
618
619	/* allocate memory for all the edges the edge map */
620	for (i = 0; i < Maxbrnch; i++) {
621	edge[i].edgemap = (int *) calloc(Maxbrnch, sizeof(int) );
622	if (edge[i].edgemap == NULL) maerror("edgemap in inittree");
623	}
624
625	/* number all edges */
626	for (i = 0; i < Maxbrnch; i++) edge[i].numedge = i;
627
628	/* initialize tree */
629
630	nextedge = 3;
631	nextleaf = 3;
632
633	/* edge maps */
634	(edge[0].edgemap)[0] = 0; /* you are on the right edge */
635	(edge[0].edgemap)[1] = 4; /* go down left for leaf 1 */
636	(edge[0].edgemap)[2] = 5; /* go down right for leaf 2 */
637	(edge[1].edgemap)[0] = 1; /* go up for leaf 0 */
638	(edge[1].edgemap)[1] = 0; /* you are on the right edge */
639	(edge[1].edgemap)[2] = 3; /* go up/down right for leaf 2 */
640	(edge[2].edgemap)[0] = 1; /* go up for leaf 0 */
641	(edge[2].edgemap)[1] = 2; /* go up/down left for leaf 1 */
642	(edge[2].edgemap)[2] = 0; /* you are on the right edge */
643
644	/* interconnection */
645	edge[0].up = NULL;
646	edge[0].downleft = &edge[1];
647	edge[0].downright = &edge[2];
648	edge[1].up = &edge[0];
649	edge[1].downleft = NULL;
650	edge[1].downright = NULL;
651	edge[2].up = &edge[0];
652	edge[2].downleft = NULL;
653	edge[2].downright = NULL;
654
655	/* edges of leaves */
656	edgeofleaf[0] = 0;
657	edgeofleaf[1] = 1;
658	edgeofleaf[2] = 2;
659	} /* inittree */
660
661	/* add next leaf on the specified edge */
662	void addnextleaf(int dockedge)
663	{
664	int i;
665
666	if (dockedge >= nextedge) {
667	/* Trying to add leaf nextleaf to nonexisting edge dockedge */
668	FPRINTF(STDOUTFILE "\n\n\nHALT: PLEASE REPORT ERROR F TO DEVELOPERS\n\n\n");
669	exit(1);
670	}
671
672	if (nextleaf >= Maxspc) {
673	/* Trying to add leaf nextleaf to a tree with Maxspc leaves */
674	FPRINTF(STDOUTFILE "\n\n\nHALT: PLEASE REPORT ERROR G TO DEVELOPERS\n\n\n");
675	exit(1);
676	}
677
678	/* necessary change in edgeofleaf if dockedge == edgeofleaf[0] */
679	if (edgeofleaf[0] == dockedge) edgeofleaf[0] = nextedge;
680
681	/* adding nextedge to the tree */
682	edge[nextedge].up = edge[dockedge].up;
683	edge[nextedge].downleft = &edge[dockedge];
684	edge[nextedge].downright = &edge[nextedge+1];
685	edge[dockedge].up = &edge[nextedge];
686
687	if (edge[nextedge].up != NULL) {
688	if ( ((edge[nextedge].up)->downleft) == &edge[dockedge] )
689	(edge[nextedge].up)->downleft = &edge[nextedge];
690	else
691	(edge[nextedge].up)->downright = &edge[nextedge];
692	}
693
694	/* adding nextedge + 1 to the tree */
695	edge[nextedge+1].up = &edge[nextedge];
696	edge[nextedge+1].downleft = NULL;
697	edge[nextedge+1].downright = NULL;
698	edgeofleaf[nextleaf] = nextedge+1;
699
700	/* the two new edges get info about the old edges */
701	/* nextedge */
702	for (i = 0; i < nextedge; i++) {
703	switch ( (edge[dockedge].edgemap)[i] ) {
704
705	/* down right changes to down left */
706	case 5: (edge[nextedge].edgemap)[i] = 4;
707	break;
708
709	/* null changes to down left */
710	case 0: (edge[nextedge].edgemap)[i] = 4;
711	break;
712
713	default: (edge[nextedge].edgemap)[i] =
714	(edge[dockedge].edgemap)[i];
715	break;
716	}
717	}
718
719	/* nextedge + 1 */
720	for (i = 0; i < nextedge; i++) {
721	switch ( (edge[dockedge].edgemap)[i] ) {
722
723	/* up/down left changes to up */
724	case 2: (edge[nextedge+1].edgemap)[i] = 1;
725	break;
726
727	/* up/down right changes to up */
728	case 3: (edge[nextedge+1].edgemap)[i] = 1;
729	break;
730
731	/* down left changes to up/down left */
732	case 4: (edge[nextedge+1].edgemap)[i] = 2;
733	break;
734
735	/* down right changes to up/down left */
736	case 5: (edge[nextedge+1].edgemap)[i] = 2;
737	break;
738
739	/* null changes to up/down left */
740	case 0: (edge[nextedge+1].edgemap)[i] = 2;
741	break;
742
743	/* up stays up */
744	default: (edge[nextedge+1].edgemap)[i] =
745	(edge[dockedge].edgemap)[i];
746	break;
747	}
748	}
749
750	/* dockedge */
751	for (i = 0; i < nextedge; i++) {
752	switch ( (edge[dockedge].edgemap)[i] ) {
753
754	/* up/down right changes to up */
755	case 3: (edge[dockedge].edgemap)[i] = 1;
756	break;
757
758	/* up/down left changes to up */
759	case 2: (edge[dockedge].edgemap)[i] = 1;
760	break;
761
762	default: break;
763	}
764	}
765
766	/* all edgemaps are updated for the two new edges */
767	/* nextedge */
768	(edge[nextedge].edgemap)[nextedge] = 0;
769	(edge[nextedge].edgemap)[nextedge+1] = 5; /* down right */
770
771	/* nextedge + 1 */
772	(edge[nextedge+1].edgemap)[nextedge] = 1; /* up */
773	(edge[nextedge+1].edgemap)[nextedge+1] = 0;
774
775	/* all other edges */
776	for (i = 0; i < nextedge; i++) {
777	(edge[i].edgemap)[nextedge] = (edge[i].edgemap)[dockedge];
778	(edge[i].edgemap)[nextedge+1] = (edge[i].edgemap)[dockedge];
779	}
780
781	/* an extra for dockedge */
782	(edge[dockedge].edgemap)[nextedge] = 1; /* up */
783	(edge[dockedge].edgemap)[nextedge+1] = 3; /* up/down right */
784
785	nextleaf++;
786	nextedge = nextedge + 2;
787	} /* addnextleaf */
788
789
790	/* free memory (to be called after inittree) */
791	void freetree()
792	{
793	int i;
794
795	for (i = 0; i < 2 * Maxspc - 3; i++) free(edge[i].edgemap);
796	free(edge);
797	free(edgeofleaf);
798	} /* freetree */
799
800	/* writes OTU sitting on edge ed */
801	void writeOTU(FILE *outfp, int ed)
802	{
803	int i;
804
805	/* test whether we are on a leaf */
806	if (edge[ed].downright == NULL && edge[ed].downleft == NULL) {
807	for (i = 1; i < nextleaf; i++) {
808	if (edgeofleaf[i] == ed) { /* i is the leaf of ed */
809	column += fputid(outfp, trueID[i]);
810	return;
811	}
812	}
813	}
814
815	/* we are NOT on a leaf */
816	fprintf(outfp, "(");
817	column++;
818	writeOTU(outfp, edge[ed].downleft->numedge);
819	fprintf(outfp, ",");
820	column++;
821	column++;
822	if (column > 55) {
823	column = 2;
824	fprintf(outfp, "\n ");
825	}
826	writeOTU(outfp, edge[ed].downright->numedge);
827	fprintf(outfp, ")");
828	column++;
829	} /* writeOTU */
830
831	/* write tree */
832	void writetree(FILE *outfp)
833	{
834	column = 1;
835	fprintf(outfp, "(");
836	column += fputid(outfp, trueID[0]) + 3;
837	fprintf(outfp, ",");
838	writeOTU(outfp, edge[edgeofleaf[0]].downleft->numedge);
839	column++;
840	column++;
841	fprintf(outfp, ",");
842	writeOTU(outfp, edge[edgeofleaf[0]].downright->numedge);
843	fprintf(outfp, ");\n");
844	} /* writetree */
845
846
847	/* clear all edgeinfos */
848	void resetedgeinfo()
849	{
850	int i;
851
852	for (i = 0; i < nextedge; i++)
853	edge[i].edgeinfo = 0;
854	} /* resetedgeinfo */
855
856	/* increment all edgeinfo between leaf A and B */
857	void incrementedgeinfo(int A, int B)
858	{
859	int curredge, finaledge, nextstep;
860
861	if (A == B) return;
862
863	finaledge = edgeofleaf[B];
864
865	curredge = edgeofleaf[A];
866	edge[curredge].edgeinfo = edge[curredge].edgeinfo + 1;
867
868	while (curredge != finaledge) {
869	nextstep = (edge[curredge].edgemap)[finaledge];
870	switch (nextstep) {
871
872	/* up */
873	case 1: curredge = (edge[curredge].up)->numedge;
874	break;
875
876	/* up/down left */
877	case 2: curredge = ((edge[curredge].up)->downleft)->numedge;
878	break;
879
880	/* up/down right */
881	case 3: curredge = ((edge[curredge].up)->downright)->numedge;
882	break;
883
884	/* down left */
885	case 4: curredge = (edge[curredge].downleft)->numedge;
886	break;
887
888	/* down right */
889	case 5: curredge = (edge[curredge].downright)->numedge;
890	break;
891
892	}
893	edge[curredge].edgeinfo = edge[curredge].edgeinfo + 1;
894	}
895	} /* incrementedgeinfo */
896
897	/* checks which edge has the lowest edgeinfo
898	if there are several edges with the same lowest edgeinfo,
899	one of them will be selected randomly */
900	void minimumedgeinfo()
901	{
902	int i, k, howmany, randomnum;
903
904	howmany = 1;
905	minedge = 0;
906	mininfo = edge[0].edgeinfo;
907	for (i = 1; i < nextedge; i++)
908	if (edge[i].edgeinfo <= mininfo) {
909	if (edge[i].edgeinfo == mininfo) {
910	howmany++;
911	} else {
912	minedge = i;
913	mininfo = edge[i].edgeinfo;
914	howmany = 1;
915	}
916	}
917
918	if (howmany > 1) { /* draw random edge */
919	randomnum = randominteger(howmany) + 1; /* 1 to howmany */
920	i = -1;
921	for (k = 0; k < randomnum; k++) {
922	do {
923	i++;
924	} while (edge[i].edgeinfo != mininfo);
925	minedge = i;
926	}
927	}
928	} /* minimumedgeinfo */
929
930
931
932
933	/*******************************************/
934	/* tree sorting */
935	/*******************************************/
936
937	/* compute address of the 4 int (sort key) in the 4 int node */
938	int ct_sortkeyaddr(int addr)
939	{
940	int a, res;
941	a = addr % 4;
942	res = addr - a + 3;
943	return res;
944	}
945
946
947	/**********/
948
949	/* compute address of the next edge pointer in a 4 int node (0->1->2->0) */
950	int ct_nextedgeaddr(int addr)
951	{
952	int a, res;
953	a = addr % 4;
954	if ( a == 2 ) { res = addr - 2; }
955	else { res = addr + 1; }
956	return res;
957	}
958
959
960	/**********/
961
962	/* compute address of 1st edge of a 4 int node from node number */
963	int ct_1stedge(int node)
964	{
965	int res;
966	res = 4 * node;
967	return res;
968	}
969
970
971	/**********/
972
973	/* compute address of 2nd edge of a 4 int node from node number */
974	int ct_2ndedge(int node)
975	{
976	int res;
977	res = 4 * node +1;
978	return res;
979	}
980
981
982	/**********/
983
984	/* compute address of 3rd edge of a 4 int node from node number */
985	int ct_3rdedge(int node)
986	{
987	int res;
988	res = 4 * node +2;
989	return res;
990	}
991
992
993	/**********/
994
995	/* check whether node 'node' is a leaf (2nd/3rd edge pointer = -1) */
996	int ct_isleaf(int node, int *ctree)
997	{
998	return (ctree[ct_3rdedge(node)] < 0);
999	}
1000
1001
1002	/**********/
1003
1004	/* compute node number of 4 int node from an edge addr. */
1005	int ct_addr2node(int addr)
1006	{
1007	int a, res;
1008	a = addr % 4;
1009	res = (int) ((addr - a) / 4);
1010	return res;
1011	}
1012
1013
1014	/**********/
1015
1016	/* print graph pointers for checking */
1017	void printctree(int *ctree)
1018	{
1019	int n;
1020	for (n=0; n < 2*Maxspc; n++) {
1021	printf("n[%3d] = (%3d.%2d, %3d.%2d, %3d.%2d \| %3d)\n", n,
1022	(int) ctree[ct_1stedge(n)]/4,
1023	(int) ctree[ct_1stedge(n)]%4,
1024	(int) ctree[ct_2ndedge(n)]/4,
1025	(int) ctree[ct_2ndedge(n)]%4,
1026	(int) ctree[ct_3rdedge(n)]/4,
1027	(int) ctree[ct_3rdedge(n)]%4,
1028	ctree[ct_3rdedge(n)+1]);
1029	}
1030	printf("\n");
1031	} /* printctree */
1032
1033
1034	/**********/
1035
1036	/* allocate memory for ctree 3 ints pointer plus 1 check byte */
1037	int *initctree()
1038	{
1039	int *snodes;
1040	int n;
1041
1042	snodes = (int ) malloc(4 2 * Maxspc * sizeof(int));
1043	if (snodes == NULL) maerror("snodes in copytree");
1044
1045	for (n=0; n<(4 * 2 * Maxspc); n++) {
1046	snodes[n]=-1;
1047	}
1048	return snodes;
1049	}
1050
1051
1052	/**********/
1053
1054	/* free memory of a tree for sorting */
1055	void freectree(int **snodes)
1056	{
1057	free(*snodes);
1058	*snodes = NULL;
1059	}
1060
1061
1062	/**********/
1063
1064	/* copy subtree recursively */
1065	void copyOTU(int ctree, / tree array struct */
1066	int ct_nextnode, / next free node */
1067	int ct_curredge, /* currende edge to add subtree */
1068	int ct_nextleaf, / next free leaf (0-maxspc) */
1069	int ed) /* edge in puzzling step tree */
1070	{
1071	int i, nextcurredge;
1072
1073	/* test whether we are on a leaf */
1074	if (edge[ed].downright == NULL && edge[ed].downleft == NULL) {
1075	for (i = 1; i < nextleaf; i++) {
1076	if (edgeofleaf[i] == ed) { /* i is the leaf of ed */
1077	nextcurredge = ct_1stedge(*ct_nextleaf);
1078	ctree[ct_curredge] = nextcurredge;
1079	ctree[nextcurredge] = ct_curredge;
1080	ctree[ct_sortkeyaddr(nextcurredge)] = trueID[i];
1081	(*ct_nextleaf)++;
1082	return;
1083	}
1084	}
1085	}
1086
1087	/* we are NOT on a leaf */
1088	nextcurredge = ct_1stedge(*ct_nextnode);
1089	ctree[ct_curredge] = nextcurredge;
1090	ctree[nextcurredge] = ct_curredge;
1091	(*ct_nextnode)++;
1092	nextcurredge = ct_nextedgeaddr(nextcurredge);
1093	copyOTU(ctree, ct_nextnode, nextcurredge,
1094	ct_nextleaf, edge[ed].downleft->numedge);
1095
1096	nextcurredge = ct_nextedgeaddr(nextcurredge);
1097	copyOTU(ctree, ct_nextnode, nextcurredge,
1098	ct_nextleaf, edge[ed].downright->numedge);
1099	}
1100
1101
1102	/**********/
1103
1104	/* copy treestructure to sorting structure */
1105	void copytree(int *ctree)
1106	{
1107	int ct_curredge;
1108	int ct_nextleaf;
1109	int ct_nextnode;
1110
1111	ct_nextnode = Maxspc;
1112	ct_curredge = ct_1stedge(ct_nextnode);
1113	ct_nextleaf = 1;
1114
1115	ctree[ct_1stedge(0)] = ct_curredge;
1116	ctree[ct_curredge] = ct_1stedge(0);
1117	ctree[ct_sortkeyaddr(0)] = trueID[0];
1118
1119	ct_nextnode++;
1120
1121	ct_curredge = ct_nextedgeaddr(ct_curredge);
1122	copyOTU(ctree, &ct_nextnode, ct_curredge,
1123	&ct_nextleaf, edge[edgeofleaf[0]].downleft->numedge);
1124
1125	ct_curredge = ct_nextedgeaddr(ct_curredge);
1126	copyOTU(ctree, &ct_nextnode, ct_curredge,
1127	&ct_nextleaf, edge[edgeofleaf[0]].downright->numedge);
1128	}
1129
1130
1131	/**********/
1132
1133	/* sort subtree from edge recursively by indices */
1134	int sortOTU(int edge, int *ctree)
1135	{
1136	int key1, key2;
1137	int edge1, edge2;
1138	int tempedge;
1139
1140	if (ctree[ct_2ndedge((int) (edge / 4))] < 0)
1141	return ctree[ct_sortkeyaddr(edge)];
1142
1143	edge1 = ctree[ct_nextedgeaddr(edge)];
1144	edge2 = ctree[ct_nextedgeaddr(ct_nextedgeaddr(edge))];
1145
1146	/* printf ("visiting [%5d] -> [%5d], [%5d]\n", edge, edge1, edge2); */
1147	/* printf ("visiting [%2d.%2d] -> [%2d.%2d], [%2d.%2d]\n",
1148	(int)(edge/4), edge%4, (int)(edge1/4), edge1%4,
1149	(int)(edge2/4), edge2%4); */
1150
1151	key1 = sortOTU(edge1, ctree);
1152	key2 = sortOTU(edge2, ctree);
1153
1154	if (key2 < key1) {
1155	tempedge = ctree[ctree[edge1]];
1156	ctree[ctree[edge1]] = ctree[ctree[edge2]];
1157	ctree[ctree[edge2]] = tempedge;
1158	tempedge = ctree[edge1];
1159	ctree[edge1] = ctree[edge2];
1160	ctree[edge2] = tempedge;
1161	ctree[ct_sortkeyaddr(edge)] = key2;
1162
1163	} else {
1164	ctree[ct_sortkeyaddr(edge)] = key1;
1165	}
1166	return ctree[ct_sortkeyaddr(edge)];
1167	}
1168
1169
1170	/**********/
1171
1172	/* sort ctree recursively by indices */
1173	int sortctree(int *ctree)
1174	{
1175	int n, startnode=-1;
1176	for(n=0; n<Maxspc; n++) {
1177	if (ctree[ct_sortkeyaddr(n*4)] == 0)
1178	startnode = n;
1179	}
1180	sortOTU(ctree[startnode * 4], ctree);
1181	return startnode;
1182	}
1183
1184
1185	/**********/
1186
1187	/* print recursively subtree of edge of sorted tree ctree */
1188	void printfsortOTU(int edge, int *ctree)
1189	{
1190	int edge1, edge2;
1191
1192	if (ctree[ct_2ndedge((int) (edge / 4))] < 0) {
1193	printf("%d", ctree[ct_sortkeyaddr(edge)]);
1194	return;
1195	}
1196
1197	edge1 = ctree[ct_nextedgeaddr(edge)];
1198	edge2 = ctree[ct_nextedgeaddr(ct_nextedgeaddr(edge))];
1199
1200	printf("(");
1201	printfsortOTU(edge1, ctree);
1202	printf(",");
1203	printfsortOTU(edge2, ctree);
1204	printf(")");
1205
1206	}
1207
1208
1209	/**********/
1210
1211	/* print recursively sorted tree ctree */
1212	int printfsortctree(int *ctree)
1213	{
1214	int n, startnode=-1;
1215	for(n=0; n<Maxspc; n++) {
1216	if (ctree[ct_sortkeyaddr(n*4)] == 0)
1217	startnode = n;
1218	}
1219	printf ("(%d,", ctree[ct_sortkeyaddr(startnode*4)]);
1220	printfsortOTU(ctree[startnode * 4], ctree);
1221	printf (");\n");
1222	return startnode;
1223	}
1224
1225
1226	/**********/
1227
1228	/* print recursively subtree of edge of sorted tree ctree to string */
1229	void sprintfOTU(char str, int len, int edge, int *ctree)
1230	{
1231	int edge1, edge2;
1232
1233	if (ctree[ct_2ndedge((int) (edge / 4))] < 0) {
1234	len+=sprintf(&(str[len]), "%d", ctree[ct_sortkeyaddr(edge)]);
1235	return;
1236	}
1237
1238	edge1 = ctree[ct_nextedgeaddr(edge)];
1239	edge2 = ctree[ct_nextedgeaddr(ct_nextedgeaddr(edge))];
1240
1241	sprintf(&(str[*len]), "(");
1242	(*len)++;
1243	sprintfOTU(str, len, edge1, ctree);
1244	sprintf(&(str[*len]), ",");
1245	(*len)++;
1246	sprintfOTU(str, len, edge2, ctree);
1247	sprintf(&(str[*len]), ")");
1248	(*len)++;
1249	}
1250
1251	/**********/
1252
1253	/* print recursively sorted tree ctree to string */
1254	char sprintfctree(int ctree, int strglen)
1255	{
1256	char *treestr,
1257	*tmpptr;
1258	int n,
1259	len=0,
1260	startnode=-1;
1261	treestr = (char ) malloc(strglen sizeof(char));
1262	tmpptr = treestr;
1263	for(n=0; n<Maxspc; n++) {
1264	if (ctree[ct_sortkeyaddr(n*4)] == 0)
1265	startnode = n;
1266	}
1267	len+=sprintf (&(tmpptr[len]), "(%d,", ctree[ct_sortkeyaddr(startnode*4)]);
1268	sprintfOTU(tmpptr, &len, ctree[startnode * 4], ctree);
1269	len+=sprintf (&(tmpptr[len]), ");");
1270	return treestr;
1271	}
1272
1273
1274	/**********/
1275
1276
1277	/***********************************************/
1278	/* establish and handle a list of sorted trees */
1279	/***********************************************/
1280
1281	int itemcount;
1282
1283	/* initialize structure */
1284	treelistitemtype inittreelist(int treenum)
1285	{
1286	*treenum = 0;
1287	return NULL;
1288	}
1289
1290
1291	/**********/
1292
1293	/* malloc new tree list item */
1294	treelistitemtype *gettreelistitem()
1295	{
1296	treelistitemtype *tmpptr;
1297	tmpptr = (treelistitemtype *)malloc(sizeof(treelistitemtype));
1298	if (tmpptr == NULL) maerror("item of intermediate tree structures");
1299	(*tmpptr).pred = NULL;
1300	(*tmpptr).succ = NULL;
1301	(*tmpptr).tree = NULL;
1302	(*tmpptr).count = 0;
1303	(*tmpptr).idx = itemcount++;
1304	return tmpptr;
1305	}
1306
1307	/**********/
1308
1309	/* free whole tree list */
1310	void freetreelist(treelistitemtype **list,
1311	int *numitems,
1312	int *numsum)
1313	{
1314	treelistitemtype *current;
1315	treelistitemtype *next;
1316	current = *list;
1317	while (current != NULL) {
1318	next = (*current).succ;
1319	if ((*current).tree != NULL) {
1320	free ((*current).tree);
1321	(*current).tree = NULL;
1322	}
1323	free(current);
1324	current = next;
1325	}
1326	*list = NULL;
1327	*numitems = 0;
1328	*numsum = 0;
1329	} /* freetreelist */
1330
1331
1332	/**********/
1333
1334	/* add tree to the tree list */
1335	treelistitemtype addtree2list(char tree, / sorted tree string */
1336	int numtrees, /* how many occurred, e.g. in parallel */
1337	treelistitemtype *list, / addr. of tree list */
1338	int *numitems,
1339	int *numsum)
1340	{
1341	treelistitemtype *tmpptr = NULL;
1342	treelistitemtype *newptr = NULL;
1343	int result;
1344	int done = 0;
1345
1346	if ((*list == NULL) \|\| (numitems == 0)) {
1347	newptr = gettreelistitem();
1348	(newptr).tree = tree;
1349	*tree = NULL;
1350	(newptr).id = numitems;
1351	(*newptr).count = numtrees;
1352	*numitems = 1;
1353	*numsum = numtrees;
1354	*list = newptr;
1355	} else {
1356	tmpptr = *list;
1357	while(done == 0) {
1358	result = strcmp( (tmpptr).tree, tree);
1359	if (result==0) {
1360	free(tree); tree = NULL;
1361	(*tmpptr).count += numtrees;
1362	*numsum += numtrees;
1363	done = 1;
1364	newptr = tmpptr;
1365	} else { if (result < 0) {
1366	if ((*tmpptr).succ != NULL)
1367	tmpptr = (*tmpptr).succ;
1368	else {
1369	newptr = gettreelistitem();
1370	(newptr).tree = tree;
1371	*tree = NULL;
1372	(newptr).id = numitems;
1373	(*newptr).count = numtrees;
1374	(*newptr).pred = tmpptr;
1375	(*tmpptr).succ = newptr;
1376	(*numitems)++;
1377	*numsum += numtrees;
1378	done = 1;
1379	}
1380	} else { /* result < 0 */
1381	newptr = gettreelistitem();
1382	(newptr).tree = tree;
1383	*tree = NULL;
1384	(newptr).id = numitems;
1385	(*newptr).count = numtrees;
1386	(*newptr).succ = tmpptr;
1387	(newptr).pred = (tmpptr).pred;
1388	(*tmpptr).pred = newptr;
1389	*numsum += numtrees;
1390
1391	if ((*newptr).pred != NULL) {
1392	((newptr).pred).succ = newptr;
1393	} else {
1394	*list = newptr;
1395	}
1396	(*numitems)++;
1397	done = 1;
1398	} /* end if result < 0 */
1399	} /* end if result != 0 */
1400	} /* while searching in list */
1401	} /* if list empty, else */
1402	return (newptr);
1403	} /* addtree2list */
1404
1405
1406	/**********/
1407
1408	/* resort list of trees by number of occurrences for output */
1409	void sortbynum(treelistitemtype list, treelistitemtype *sortlist)
1410	{
1411	treelistitemtype *tmpptr = NULL;
1412	treelistitemtype *curr = NULL;
1413	treelistitemtype *next = NULL;
1414	int xchange = 1;
1415
1416	if (list == NULL) fprintf(stderr, "Grrrrrrrrr>>>>\n");
1417	tmpptr = list;
1418	*sortlist = list;
1419	while (tmpptr != NULL) {
1420	(tmpptr).sortnext = (tmpptr).succ;
1421	(tmpptr).sortlast = (tmpptr).pred;
1422	tmpptr = (*tmpptr).succ;
1423	}
1424
1425	while (xchange > 0) {
1426	curr = *sortlist;
1427	xchange = 0;
1428	if (curr == NULL) fprintf(stderr, "Grrrrrrrrr>>>>\n");
1429	while((*curr).sortnext != NULL) {
1430	next = (*curr).sortnext;
1431	if ((curr).count >= (next).count)
1432	curr = (*curr).sortnext;
1433	else {
1434	if ((*curr).sortlast != NULL)
1435	(((curr).sortlast)).sortnext = next;
1436	if (*sortlist == curr)
1437	*sortlist = next;
1438	(next).sortlast = (curr).sortlast;
1439
1440	if ((*next).sortnext != NULL)
1441	(((next).sortnext)).sortlast = curr;
1442	(curr).sortnext = (next).sortnext;
1443
1444	(*curr).sortlast = next;
1445	(*next).sortnext = curr;
1446
1447	xchange++;
1448	}
1449	}
1450	}
1451	} /* sortbynum */
1452
1453
1454	/**********/
1455
1456	/* print puzzling step tree structures for checking */
1457	void printfpstrees(treelistitemtype *list)
1458	{
1459	char ch;
1460	treelistitemtype *tmpptr = NULL;
1461	tmpptr = list;
1462	ch = '-';
1463	while (tmpptr != NULL) {
1464	printf ("%c[%2d] %5d %s\n", ch, (tmpptr).idx, (tmpptr).count, (*tmpptr).tree);
1465	tmpptr = (*tmpptr).succ;
1466	ch = ' ';
1467	}
1468	}
1469
1470	/**********/
1471
1472	/* print sorted puzzling step tree structure with names */
1473	void fprintffullpstree(FILE outf, char treestr)
1474	{
1475	int count = 0;
1476	int idnum = 0;
1477	int n;
1478	for(n=0; treestr[n] != '\0'; n++){
1479	while(isdigit((int)treestr[n])){
1480	idnum = (10 * idnum) + ((int)treestr[n]-48);
1481	n++;
1482	count++;
1483	}
1484	if (count > 0){
1485	# ifdef USEQUOTES
1486	fprintf(outf, "'");
1487	# endif
1488	(void)fputid(outf, idnum);
1489	# ifdef USEQUOTES
1490	fprintf(outf, "'");
1491	# endif
1492	count = 0;
1493	idnum = 0;
1494	}
1495	fprintf(outf, "%c", treestr[n]);
1496	}
1497	}
1498
1499
1500	/**********/
1501
1502	/* print sorted puzzling step tree structures with names */
1503	void fprintfsortedpstrees(FILE *output,
1504	treelistitemtype list, / tree list */
1505	int itemnum, /* order number */
1506	int itemsum, /* number of trees */
1507	int comment, /* with statistics, or puzzle report ? */
1508	float cutoff) /* cutoff percentage */
1509	{
1510	treelistitemtype *tmpptr = NULL;
1511	treelistitemtype *slist = NULL;
1512	int num = 1;
1513	float percent;
1514
1515	if (list == NULL) fprintf(stderr, "Grrrrrrrrr>>>>\n");
1516	sortbynum(list, &slist);
1517
1518	tmpptr = slist;
1519	while (tmpptr != NULL) {
1520	percent = (float)(100.0 * (*tmpptr).count / itemsum);
1521	if ((cutoff == 0.0) \|\| (cutoff <= percent)) {
1522	if (comment)
1523	fprintf (output, "[ %d. %d %.2f %d %d %d ]", num++, (tmpptr).count, percent, (tmpptr).id, itemnum, itemsum);
1524	else {
1525	if (num == 1){
1526	fprintf (output, "\n");
1527	fprintf (output, "The following tree(s) occured in more than %.2f%% of the %d puzzling steps.\n", cutoff, itemsum);
1528	fprintf (output, "The trees are orderd descending by the number of occurrences.\n");
1529	fprintf (output, "\n");
1530	fprintf (output, "\n occurrences ID Phylip tree\n");
1531	}
1532	fprintf (output, "%2d. %5d %6.2f%% %5d ", num++, (tmpptr).count, percent, (tmpptr).id);
1533	}
1534	fprintffullpstree(output, (*tmpptr).tree);
1535	fprintf (output, "\n");
1536	}
1537	tmpptr = (*tmpptr).sortnext;
1538	}
1539
1540	if (!comment) {
1541	fprintf (output, "\n");
1542	switch(num) {
1543	case 1: fprintf (output, "There were no tree topologies (out of %d) occuring with a percentage >= %.2f%% of the %d puzzling steps.\n", itemnum, cutoff, itemsum); break;
1544	case 2: fprintf (output, "There was one tree topology (out of %d) occuring with a percentage >= %.2f%%.\n", itemnum, cutoff); break;
1545	default: fprintf (output, "There were %d tree topologies (out of %d) occuring with a percentage >= %.2f%%.\n", num-1, itemnum, cutoff); break;
1546	}
1547	fprintf (output, "\n");
1548	fprintf (output, "\n");
1549	}
1550
1551	} /* fprintfsortedpstrees */
1552
1553	/**********/
1554
1555	/* print sorted tree topologies for checking */
1556	void printfsortedpstrees(treelistitemtype *list)
1557	{
1558	treelistitemtype *tmpptr = NULL;
1559	treelistitemtype *slist = NULL;
1560
1561	sortbynum(list, &slist);
1562
1563	tmpptr = slist;
1564	while (tmpptr != NULL) {
1565	printf ("[%2d] %5d %s\n", (tmpptr).idx, (tmpptr).count, (*tmpptr).tree);
1566	tmpptr = (*tmpptr).sortnext;
1567	}
1568	} /* printfsortedpstrees */
1569
1570
1571	/*******************************************/
1572	/* end of tree sorting */
1573	/*******************************************/
1574
1575
1576
1577	/******************************************************************************/
1578	/* functions for computing the consensus tree */
1579	/******************************************************************************/
1580
1581	/* prepare for consensus tree analysis */
1582	void initconsensus()
1583	{
1584	# if ! PARALLEL
1585	biparts = new_cmatrix(Maxspc-3, Maxspc);
1586	# endif /* PARALLEL */
1587
1588	if (Maxspc % 32 == 0)
1589	splitlength = Maxspc/32;
1590	else splitlength = (Maxspc + 32 - (Maxspc % 32))/32;
1591	numbiparts = 0; /* no pattern stored so far */
1592	maxbiparts = 0; /* no memory reserved so far */
1593	splitfreqs = NULL;
1594	splitpatterns = NULL;
1595	splitsizes = NULL;
1596	splitcomp = (uli ) malloc(splitlength sizeof(uli) );
1597	if (splitcomp == NULL) maerror("splitcomp in initconsensus");
1598	}
1599
1600	/* prototype needed for recursive function */
1601	void makepart(int i, int curribrnch);
1602
1603	/* recursive function to get bipartitions */
1604	void makepart(int i, int curribrnch)
1605	{
1606	int j;
1607
1608	if ( edge[i].downright == NULL \|\|
1609	edge[i].downleft == NULL) { /* if i is leaf */
1610
1611	/* check out what leaf j sits on this edge i */
1612	for (j = 1; j < Maxspc; j++) {
1613	if (edgeofleaf[j] == i) {
1614	biparts[curribrnch][trueID[j]] = '*';
1615	return;
1616	}
1617	}
1618	} else { /* still on inner branch */
1619	makepart(edge[i].downleft->numedge, curribrnch);
1620	makepart(edge[i].downright->numedge, curribrnch);
1621	}
1622	}
1623
1624	/* compute bipartitions of tree of current puzzling step */
1625	void computebiparts()
1626	{
1627	int i, j, curribrnch;
1628
1629	curribrnch = -1;
1630
1631	for (i = 0; i < Maxspc - 3; i++)
1632	for (j = 0; j < Maxspc; j++)
1633	biparts[i][j] = '.';
1634
1635	for (i = 0; i < Maxbrnch; i++) {
1636	if (!( edgeofleaf[0] == i \|\|
1637	edge[i].downright == NULL \|\|
1638	edge[i].downleft == NULL) ) { /* check all inner branches */
1639	curribrnch++;
1640	makepart(i, curribrnch);
1641
1642	/* make sure that the root is always a '' /
1643	if (biparts[curribrnch][outgroup] == '.') {
1644	for (j = 0; j < Maxspc; j++) {
1645	if (biparts[curribrnch][j] == '.')
1646	biparts[curribrnch][j] = '*';
1647	else
1648	biparts[curribrnch][j] = '.';
1649	}
1650	}
1651	}
1652	}
1653	}
1654
1655	/* print out the bipartition n of all different splitpatterns */
1656	void printsplit(FILE *fp, uli n)
1657	{
1658	int i, j, col;
1659	uli z;
1660
1661	col = 0;
1662	for (i = 0; i < splitlength; i++) {
1663	z = splitpatterns[n*splitlength + i];
1664	for (j = 0; j < 32 && col < Maxspc; j++) {
1665	if (col % 10 == 0 && col != 0) fprintf(fp, " ");
1666	if (z & 1) fprintf(fp, ".");
1667	else fprintf(fp, "*");
1668	z = (z >> 1);
1669	col++;
1670	}
1671	}
1672	}
1673
1674	/* make new entries for new different bipartitions and count frequencies */
1675	void makenewsplitentries()
1676	{
1677	int i, j, bpc, identical, idflag, bpsize;
1678	uli nextentry, obpc;
1679
1680	/* where the next entry would be in splitpatterns */
1681	nextentry = numbiparts;
1682
1683	for (bpc = 0; bpc < Maxspc - 3; bpc++) { /* for every new bipartition */
1684	/* convert bipartition into a more compact format */
1685	bpsize = 0;
1686	for (i = 0; i < splitlength; i++) {
1687	splitcomp[i] = 0;
1688	for (j = 0; j < 32; j++) {
1689	splitcomp[i] = splitcomp[i] >> 1;
1690	if (i*32 + j < Maxspc)
1691	if (biparts[bpc][i*32 + j] == '.') {
1692	/* set highest bit */
1693	splitcomp[i] = (splitcomp[i] \| 2147483648UL);
1694	bpsize++; /* count the '.' */
1695	}
1696	}
1697	}
1698	/* compare to the old patterns */
1699	identical = FALSE;
1700	for (obpc = 0; (obpc < numbiparts) && (!identical); obpc++) {
1701	/* compare first partition size */
1702	if (splitsizes[obpc] == bpsize) idflag = TRUE;
1703	else idflag = FALSE;
1704	/* if size is identical compare whole partition */
1705	for (i = 0; (i < splitlength) && idflag; i++)
1706	if (splitcomp[i] != splitpatterns[obpc*splitlength + i])
1707	idflag = FALSE;
1708	if (idflag) identical = TRUE;
1709	}
1710	if (identical) { /* if identical increase frequency */
1711	splitfreqs[2*(obpc-1)]++;
1712	} else { /* create new entry */
1713	if (nextentry == maxbiparts) { /* reserve more memory */
1714	maxbiparts = maxbiparts + 2*Maxspc;
1715	splitfreqs = (uli *) myrealloc(splitfreqs,
1716	2maxbiparts sizeof(uli) );
1717	/* 2x: splitfreqs contains also an index (sorting!) */
1718	if (splitfreqs == NULL) maerror("splitfreqs in makenewsplitentries");
1719	splitpatterns = (uli *) myrealloc(splitpatterns,
1720	splitlengthmaxbiparts sizeof(uli) );
1721	if (splitpatterns == NULL) maerror("splitpatterns in makenewsplitentries");
1722	splitsizes = (int *) myrealloc(splitsizes,
1723	maxbiparts * sizeof(int) );
1724	if (splitsizes == NULL) maerror("splitsizes in makenewsplitentries");
1725	}
1726	splitfreqs[2nextentry] = 1; / frequency */
1727	splitfreqs[2nextentry+1] = nextentry; / index for sorting */
1728	for (i = 0; i < splitlength; i++)
1729	splitpatterns[nextentry*splitlength + i] = splitcomp[i];
1730	splitsizes[nextentry] = bpsize;
1731	nextentry++;
1732	}
1733	}
1734	numbiparts = nextentry;
1735	}
1736
1737	/* general remarks:
1738
1739	- every entry in consbiparts is one node of the consensus tree
1740	- for each node one has to know which taxa and which other nodes
1741	are directly descending from it
1742	- for every taxon/node number there is a flag that shows
1743	whether it descends from the node or not
1744	- '0' means that neither a taxon nor another node with the
1745	corresponding number decends from the node
1746	'1' means that the corresponding taxon descends from the node
1747	'2' means that the corresponding node descends from the node
1748	'3' means that the corresponding taxon and node descends from the node
1749	*/
1750
1751	/* copy bipartition n of all different splitpatterns to consbiparts[k] */
1752	void copysplit(uli n, int k)
1753	{
1754	int i, j, col;
1755	uli z;
1756
1757	col = 0;
1758	for (i = 0; i < splitlength; i++) {
1759	z = splitpatterns[n*splitlength + i];
1760	for (j = 0; j < 32 && col < Maxspc; j++) {
1761	if (z & 1) consbiparts[k][col] = '1';
1762	else consbiparts[k][col] = '0';
1763	z = (z >> 1);
1764	col++;
1765	}
1766	}
1767	}
1768
1769	/* compute majority rule consensus tree */
1770	void makeconsensus()
1771	{
1772	int i, j, k, size, subnode;
1773	char chari, charj;
1774
1775	/* sort bipartition frequencies */
1776	qsort(splitfreqs, numbiparts, 2*sizeof(uli), ulicmp);
1777	/* how many bipartitions are included in the consensus tree */
1778	consincluded = 0;
1779	for (i = 0; i < numbiparts && i == consincluded; i++) {
1780	if (2splitfreqs[2i] > Numtrial) consincluded = i + 1;
1781	}
1782
1783	/* collect all info about majority rule consensus tree */
1784	/* the +1 is due to the edge with the root */
1785	consconfid = new_ivector(consincluded + 1);
1786	conssizes = new_ivector(2*consincluded + 2);
1787	consbiparts = new_cmatrix(consincluded + 1, Maxspc);
1788
1789	for (i = 0; i < consincluded; i++) {
1790	/* copy partition to consbiparts */
1791	copysplit(splitfreqs[2*i+1], i);
1792	/* frequency in percent (rounded to integer) */
1793	consconfid[i] = (int) floor(100.0splitfreqs[2i]/Numtrial + 0.5);
1794	/* size of partition */
1795	conssizes[2i] = splitsizes[splitfreqs[2i+1]];
1796	conssizes[2*i+1] = i;
1797	}
1798	for (i = 0; i < Maxspc; i++) consbiparts[consincluded][i] = '1';
1799	consbiparts[consincluded][outgroup] = '0';
1800	consconfid[consincluded] = 100;
1801	conssizes[2*consincluded] = Maxspc - 1;
1802	conssizes[2*consincluded + 1] = consincluded;
1803
1804	/* sort bipartitions according to cluster size */
1805	qsort(conssizes, consincluded + 1, 2*sizeof(int), intcmp);
1806
1807	/* reconstruct consensus tree */
1808	for (i = 0; i < consincluded; i++) { /* try every node */
1809	size = conssizes[2i]; / size of current node */
1810	for (j = i + 1; j < consincluded + 1; j++) {
1811
1812	/* compare only with nodes with more descendants */
1813	if (size == conssizes[2*j]) continue;
1814
1815	/* check whether node i is a subnode of j */
1816	subnode = FALSE;
1817	for (k = 0; k < Maxspc && !subnode; k++) {
1818	chari = consbiparts[ conssizes[2*i+1] ][k];
1819	if (chari != '0') {
1820	charj = consbiparts[ conssizes[2*j+1] ][k];
1821	if (chari == charj \|\| charj == '3') subnode = TRUE;
1822	}
1823	}
1824
1825	/* if i is a subnode of j change j accordingly */
1826	if (subnode) {
1827	/* remove subnode i from j */
1828	for (k = 0; k < Maxspc; k++) {
1829	chari = consbiparts[ conssizes[2*i+1] ][k];
1830	if (chari != '0') {
1831	charj = consbiparts[ conssizes[2*j+1] ][k];
1832	if (chari == charj)
1833	consbiparts[ conssizes[2*j+1] ][k] = '0';
1834	else if (charj == '3') {
1835	if (chari == '1')
1836	consbiparts[ conssizes[2*j+1] ][k] = '2';
1837	else if (chari == '2')
1838	consbiparts[ conssizes[2*j+1] ][k] = '1';
1839	else {
1840	/* Consensus tree [1] */
1841	FPRINTF(STDOUTFILE "\n\n\nHALT: PLEASE REPORT ERROR H TO DEVELOPERS\n\n\n");
1842	exit(1);
1843	}
1844	} else {
1845	/* Consensus tree [2] */
1846	FPRINTF(STDOUTFILE "\n\n\nHALT: PLEASE REPORT ERROR I TO DEVELOPERS\n\n\n");
1847	exit(1);
1848	}
1849	}
1850	}
1851	/* add link to subnode i in node j */
1852	charj = consbiparts[ conssizes[2j+1] ][ conssizes[2i+1] ];
1853	if (charj == '0')
1854	consbiparts[ conssizes[2j+1] ][ conssizes[2i+1] ] = '2';
1855	else if (charj == '1')
1856	consbiparts[ conssizes[2j+1] ][ conssizes[2i+1] ] = '3';
1857	else {
1858	/* Consensus tree [3] */
1859	FPRINTF(STDOUTFILE "\n\n\nHALT: PLEASE REPORT ERROR J TO DEVELOPERS\n\n\n");
1860	exit(1);
1861	}
1862	}
1863	}
1864	}
1865	}
1866
1867	/* prototype for recursion */
1868	void writenode(FILE *treefile, int node);
1869
1870	/* write node (writeconsensustree) */
1871	void writenode(FILE *treefile, int node)
1872	{
1873	int i, first;
1874
1875	fprintf(treefile, "(");
1876	column++;
1877	/* write descending nodes */
1878	first = TRUE;
1879	for (i = 0; i < Maxspc; i++) {
1880	if (consbiparts[node][i] == '2' \|\|
1881	consbiparts[node][i] == '3') {
1882	if (first) first = FALSE;
1883	else {
1884	fprintf(treefile, ",");
1885	column++;
1886	}
1887	if (column > 60) {
1888	column = 2;
1889	fprintf(treefile, "\n");
1890	}
1891	/* write node i */
1892	writenode(treefile, i);
1893
1894	/* reliability value as internal label */
1895	fprintf(treefile, "%d", consconfid[i]);
1896
1897	column = column + 3;
1898	}
1899	}
1900	/* write descending taxa */
1901	for (i = 0; i < Maxspc; i++) {
1902	if (consbiparts[node][i] == '1' \|\|
1903	consbiparts[node][i] == '3') {
1904	if (first) first = FALSE;
1905	else {
1906	fprintf(treefile, ",");
1907	column++;
1908	}
1909	if (column > 60) {
1910	column = 2;
1911	fprintf(treefile, "\n");
1912	}
1913	column += fputid(treefile, i);
1914	}
1915	}
1916	fprintf(treefile, ")");
1917	column++;
1918	}
1919
1920	/* write consensus tree */
1921	void writeconsensustree(FILE *treefile)
1922	{
1923	int i, first;
1924
1925	column = 1;
1926	fprintf(treefile, "(");
1927	column += fputid(treefile, outgroup) + 2;
1928	fprintf(treefile, ",");
1929	/* write descending nodes */
1930	first = TRUE;
1931	for (i = 0; i < Maxspc; i++) {
1932	if (consbiparts[consincluded][i] == '2' \|\|
1933	consbiparts[consincluded][i] == '3') {
1934	if (first) first = FALSE;
1935	else {
1936	fprintf(treefile, ",");
1937	column++;
1938	}
1939	if (column > 60) {
1940	column = 2;
1941	fprintf(treefile, "\n");
1942	}
1943	/* write node i */
1944	writenode(treefile, i);
1945
1946	/* reliability value as internal label */
1947	fprintf(treefile, "%d", consconfid[i]);
1948
1949	column = column + 3;
1950	}
1951	}
1952	/* write descending taxa */
1953	for (i = 0; i < Maxspc; i++) {
1954	if (consbiparts[consincluded][i] == '1' \|\|
1955	consbiparts[consincluded][i] == '3') {
1956	if (first) first = FALSE;
1957	else {
1958	fprintf(treefile, ",");
1959	column++;
1960	}
1961	if (column > 60) {
1962	column = 2;
1963	fprintf(treefile, "\n");
1964	}
1965	column += fputid(treefile, i);
1966	}
1967	}
1968	fprintf(treefile, ");\n");
1969	}
1970
1971	/* prototype for recursion */
1972	void nodecoordinates(int node);
1973
1974	/* establish node coordinates (plotconsensustree) */
1975	void nodecoordinates(int node)
1976	{
1977	int i, ymin, ymax, xcoordinate;
1978
1979	/* first establish coordinates of descending nodes */
1980	for (i = 0; i < Maxspc; i++) {
1981	if (consbiparts[node][i] == '2' \|\|
1982	consbiparts[node][i] == '3')
1983	nodecoordinates(i);
1984	}
1985
1986	/* then establish coordinates of descending taxa */
1987	for (i = 0; i < Maxspc; i++) {
1988	if (consbiparts[node][i] == '1' \|\|
1989	consbiparts[node][i] == '3') {
1990	/* y-coordinate of taxon i */
1991	ycortax[i] = ytaxcounter;
1992	ytaxcounter = ytaxcounter - 2;
1993	}
1994	}
1995
1996	/* then establish coordinates of this node */
1997	ymin = 2*Maxspc - 2;
1998	ymax = 0;
1999	xcoordinate = 0;
2000	for (i = 0; i < Maxspc; i++) {
2001	if (consbiparts[node][i] == '2' \|\|
2002	consbiparts[node][i] == '3') {
2003	if (ycor[i] > ymax) ymax = ycor[i];
2004	if (ycor[i] < ymin) ymin = ycor[i];
2005	if (xcor[i] > xcoordinate) xcoordinate = xcor[i];
2006	}
2007	}
2008	for (i = 0; i < Maxspc; i++) {
2009	if (consbiparts[node][i] == '1' \|\|
2010	consbiparts[node][i] == '3') {
2011	if (ycortax[i] > ymax) ymax = ycortax[i];
2012	if (ycortax[i] < ymin) ymin = ycortax[i];
2013	}
2014	}
2015	ycormax[node] = ymax;
2016	ycormin[node] = ymin;
2017	ycor[node] = (int) floor(0.5*(ymax + ymin) + 0.5);
2018	if (xcoordinate == 0) xcoordinate = 9;
2019	xcor[node] = xcoordinate + 4;
2020	}
2021
2022	/* prototype for recursion */
2023	void drawnode(int node, int xold);
2024
2025	/* drawnode (plotconsensustree) */
2026	void drawnode(int node, int xold)
2027	{
2028	int i, j;
2029	char buf[4];
2030
2031	/* first draw vertical line */
2032	for (i = ycormin[node] + 1; i < ycormax[node]; i++)
2033	treepict[xcor[node]][i] = ':';
2034
2035	/* then draw descending nodes */
2036	for (i = 0; i < Maxspc; i++) {
2037	if (consbiparts[node][i] == '2' \|\|
2038	consbiparts[node][i] == '3')
2039	drawnode(i, xcor[node]);
2040	}
2041
2042	/* then draw descending taxa */
2043	for (i = 0; i < Maxspc; i++) {
2044	if (consbiparts[node][i] == '1' \|\|
2045	consbiparts[node][i] == '3') {
2046	treepict[xcor[node]][ycortax[i]] = ':';
2047	for (j = xcor[node] + 1; j < xsize-10; j++)
2048	treepict[j][ycortax[i]] = '-';
2049	for (j = 0; j < 10; j++)
2050	treepict[xsize-10+j][ycortax[i]] = Identif[i][j];
2051	}
2052	}
2053
2054	/* then draw internal edge with consensus value */
2055	treepict[xold][ycor[node]] = ':';
2056	treepict[xcor[node]][ycor[node]] = ':';
2057	for (i = xold + 1; i < xcor[node]-3; i++)
2058	treepict[i][ycor[node]] = '-';
2059	sprintf(buf, "%d", consconfid[node]);
2060	if (consconfid[node] == 100) {
2061	treepict[xcor[node]-3][ycor[node]] = buf[0];
2062	treepict[xcor[node]-2][ycor[node]] = buf[1];
2063	treepict[xcor[node]-1][ycor[node]] = buf[2];
2064	} else {
2065	treepict[xcor[node]-3][ycor[node]] = '-';
2066	treepict[xcor[node]-2][ycor[node]] = buf[0];
2067	treepict[xcor[node]-1][ycor[node]] = buf[1];
2068	}
2069	}
2070
2071	/* plot consensus tree */
2072	void plotconsensustree(FILE *plotfp)
2073	{
2074	int i, j, yroot, startree;
2075
2076	/* star tree or no star tree */
2077	if (consincluded == 0) {
2078	startree = TRUE;
2079	consincluded = 1; /* avoids problems with malloc */
2080	} else
2081	startree = FALSE;
2082
2083	/* memory for x-y-coordinates of each bipartition */
2084	xcor = new_ivector(consincluded);
2085	ycor = new_ivector(consincluded);
2086	ycormax = new_ivector(consincluded);
2087	ycormin = new_ivector(consincluded);
2088	if (startree) consincluded = 0; /* avoids problems with malloc */
2089
2090	/* y-coordinates of each taxon */
2091	ycortax = new_ivector(Maxspc);
2092	ycortax[outgroup] = 0;
2093
2094	/* establish coordinates */
2095	ytaxcounter = 2*Maxspc - 2;
2096
2097	/* first establish coordinates of descending nodes */
2098	for (i = 0; i < Maxspc; i++) {
2099	if (consbiparts[consincluded][i] == '2' \|\|
2100	consbiparts[consincluded][i] == '3')
2101	nodecoordinates(i);
2102	}
2103
2104	/* then establish coordinates of descending taxa */
2105	for (i = 0; i < Maxspc; i++) {
2106	if (consbiparts[consincluded][i] == '1' \|\|
2107	consbiparts[consincluded][i] == '3') {
2108	/* y-coordinate of taxon i */
2109	ycortax[i] = ytaxcounter;
2110	ytaxcounter = ytaxcounter - 2;
2111	}
2112	}
2113
2114	/* then establish length of root edge and size of whole tree */
2115	yroot = 0;
2116	xsize = 0;
2117	for (i = 0; i < Maxspc; i++) {
2118	if (consbiparts[consincluded][i] == '2' \|\|
2119	consbiparts[consincluded][i] == '3') {
2120	if (ycor[i] > yroot) yroot = ycor[i];
2121	if (xcor[i] > xsize) xsize = xcor[i];
2122	}
2123	}
2124	for (i = 0; i < Maxspc; i++) {
2125	if (consbiparts[consincluded][i] == '1' \|\|
2126	consbiparts[consincluded][i] == '3') {
2127	if (ycortax[i] > yroot) yroot = ycortax[i];
2128	}
2129	}
2130	if (xsize == 0) xsize = 9;
2131	/* size in x direction inclusive one blank on the left */
2132	xsize = xsize + 6;
2133
2134	/* change all x-labels so that (0,0) is down-left */
2135	for (i = 0; i < consincluded; i++)
2136	xcor[i] = xsize-1-xcor[i];
2137
2138	/* draw tree */
2139	treepict = new_cmatrix(xsize, 2*Maxspc-1);
2140	for (i = 0; i < xsize; i++)
2141	for (j = 0; j < 2*Maxspc-1; j++)
2142	treepict[i][j] = ' ';
2143
2144	/* draw root */
2145	for (i = 1; i < yroot; i++)
2146	treepict[1][i] = ':';
2147	treepict[1][0] = ':';
2148	for (i = 2; i < xsize - 10; i++)
2149	treepict[i][0] = '-';
2150	for (i = 0; i < 10; i++)
2151	treepict[xsize-10+i][0] = Identif[outgroup][i];
2152
2153	/* then draw descending nodes */
2154	for (i = 0; i < Maxspc; i++) {
2155	if (consbiparts[consincluded][i] == '2' \|\|
2156	consbiparts[consincluded][i] == '3')
2157	drawnode(i, 1);
2158	}
2159
2160	/* then draw descending taxa */
2161	for (i = 0; i < Maxspc; i++) {
2162	if (consbiparts[consincluded][i] == '1' \|\|
2163	consbiparts[consincluded][i] == '3') {
2164	treepict[1][ycortax[i]] = ':';
2165	for (j = 2; j < xsize-10; j++)
2166	treepict[j][ycortax[i]] = '-';
2167	for (j = 0; j < 10; j++)
2168	treepict[xsize-10+j][ycortax[i]] = Identif[i][j];
2169	}
2170	}
2171
2172	/* plot tree */
2173	for (i = 2*Maxspc-2; i > -1; i--) {
2174	for (j = 0; j < xsize; j++)
2175	fputc(treepict[j][i], plotfp);
2176	fputc('\n', plotfp);
2177	}
2178
2179	free_ivector(xcor);
2180	free_ivector(ycor);
2181	free_ivector(ycormax);
2182	free_ivector(ycormin);
2183	free_ivector(ycortax);
2184	free_cmatrix(treepict);
2185	}
2186
2187
2188
2189	/******************************************************************************/
2190	/* storing and evaluating quartet branching information */
2191	/******************************************************************************/
2192
2193	/* general remarks:
2194
2195	for a quartet with the taxa a, b, c, d there are
2196	three possible binary trees:
2197
2198	1) (a,b)-(c,d)
2199	2) (a,c)-(b,d)
2200	3) (a,d)-(b,c)
2201
2202	For every quartet information about its branching structure is
2203	stored. With the functions readquartet and writequartet
2204	this information can be accessed. For every quartet (a,b,c,d)
2205	with a < b < c < d (taxa) the branching information is encoded
2206	using 4 bits:
2207
2208	value 8 4 2 1
2209	+-------------+-------------+-------------+-------------+
2210	\| not used \| tree 3 \| tree 2 \| tree 1 \|
2211	+-------------+-------------+-------------+-------------+
2212
2213	If the branching structure of the taxa corresponds to one of the
2214	three trees the corresponding bit is set. If the branching structure
2215	is unclear because two of the three trees have the same maximum
2216	likelihood value the corresponding two bits are set. If the branching
2217	structure is completely unknown all the bits are set (the highest
2218	bit is always cleared because it is not used).
2219
2220	*/
2221
2222	/* allocate memory for quartets */
2223	unsigned char *mallocquartets(int taxa)
2224	{
2225	uli nc, numch;
2226	unsigned char *qinfo;
2227
2228	/* compute number of quartets */
2229	Numquartets = (uli) taxa(taxa-1)(taxa-2)*(taxa-3)/24;
2230	if (Numquartets % 2 == 0) { /* even number */
2231	numch = Numquartets/2;
2232	} else { /* odd number */
2233	numch = (Numquartets + 1)/2;
2234	}
2235	/* allocate memory */
2236	qinfo = (unsigned char ) malloc(numch sizeof(unsigned char) );
2237	if (qinfo == NULL) maerror("quartetinfo in mallocquartets");
2238	for (nc = 0; nc < numch; nc++) qinfo[nc] = 0;
2239	return(qinfo);
2240	}
2241
2242	/* free quartet memory */
2243	void freequartets()
2244	{
2245	free(quartetinfo);
2246	}
2247
2248	/* read quartet info - a < b < c < d */
2249	unsigned char readquartet(int a, int b, int c, int d)
2250	{
2251	uli qnum;
2252
2253	qnum = (uli) a
2254	+ (uli) b*(b-1)/2
2255	+ (uli) c(c-1)(c-2)/6
2256	+ (uli) d(d-1)(d-2)*(d-3)/24;
2257	if (qnum % 2 == 0) { /* even number */
2258	/* bits 0 to 3 */
2259	return (quartetinfo[qnum/2] & (unsigned char) 15);
2260	} else { /* odd number */
2261	/* bits 4 to 7 */
2262	return ((quartetinfo[(qnum-1)/2] & (unsigned char) 240)>>4);
2263	}
2264	}
2265
2266	/* write quartet info - a < b < c < d, 0 <= info <= 15 */
2267	void writequartet(int a, int b, int c, int d, unsigned char info)
2268	{
2269	uli qnum;
2270
2271	qnum = (uli) a
2272	+ (uli) b*(b-1)/2
2273	+ (uli) c(c-1)(c-2)/6
2274	+ (uli) d(d-1)(d-2)*(d-3)/24;
2275	if (qnum % 2 == 0) { /* even number */
2276	/* bits 0 to 3 */
2277	quartetinfo[qnum/2] =
2278	((quartetinfo[qnum/2] & (unsigned char) 240) \|
2279	(info & (unsigned char) 15));
2280	} else { /* odd number */
2281	/* bits 4 to 7 */
2282	quartetinfo[(qnum-1)/2] =
2283	((quartetinfo[(qnum-1)/2] & (unsigned char) 15) \|
2284	((info & (unsigned char) 15)<<4));
2285	}
2286	}
2287
2288	/* prototypes */
2289	void openfiletowrite(FILE **, char[], char[]);
2290	void closefile(FILE *);
2291
2292	/* sorts three doubles in descending order */
2293	void sort3doubles(dvector num, ivector order)
2294	{
2295	if (num[0] > num[1]) {
2296	if(num[2] > num[0]) {
2297	order[0] = 2;
2298	order[1] = 0;
2299	order[2] = 1;
2300	} else if (num[2] < num[1]) {
2301	order[0] = 0;
2302	order[1] = 1;
2303	order[2] = 2;
2304	} else {
2305	order[0] = 0;
2306	order[1] = 2;
2307	order[2] = 1;
2308	}
2309	} else {
2310	if(num[2] > num[1]) {
2311	order[0] = 2;
2312	order[1] = 1;
2313	order[2] = 0;
2314	} else if (num[2] < num[0]) {
2315	order[0] = 1;
2316	order[1] = 0;
2317	order[2] = 2;
2318	} else {
2319	order[0] = 1;
2320	order[1] = 2;
2321	order[2] = 0;
2322	}
2323	}
2324	}
2325
2326	/* checks out all possible quartets */
2327	void computeallquartets()
2328	{
2329	double onethird;
2330	uli nq;
2331	unsigned char treebits[3];
2332	FILE *lhfp;
2333	# if ! PARALLEL
2334	int a, b, c, i;
2335	double qc2, mintogo, minutes, hours, temp;
2336	double temp1, temp2, temp3;
2337	unsigned char discreteweight[3];
2338	# endif
2339
2340	onethird = 1.0/3.0;
2341	treebits[0] = (unsigned char) 1;
2342	treebits[1] = (unsigned char) 2;
2343	treebits[2] = (unsigned char) 4;
2344
2345	if (show_optn) { /* list all unresolved quartets */
2346	openfiletowrite(&unresfp, UNRESOLVED, "unresolved quartet trees");
2347	fprintf(unresfp, "List of all completely unresolved quartets:\n\n");
2348	}
2349
2350	nq = 0;
2351	badqs = 0;
2352
2353	/* start timer - percentage of completed quartets */
2354	time(&time0);
2355	time1 = time0;
2356	mflag = 0;
2357
2358	# if PARALLEL
2359	{
2360	schedtype sched;
2361	int flag;
2362	MPI_Status stat;
2363	int dest = 1;
2364	uli qaddr =0;
2365	uli qamount=0;
2366	int qblocksent = 0;
2367	int apr;
2368	uli sq, noq;
2369	initsched(&sched, numquarts(Maxspc), PP_NumProcs-1, 4);
2370	qamount=sgss(&sched);
2371	while (qamount > 0) {
2372	if (PP_emptyslave()) {
2373	PP_RecvQuartBlock(0, &sq, &noq, quartetinfo, &apr);
2374	qblocksent -= noq;
2375	}
2376	dest = PP_getslave();
2377	PP_SendDoQuartBlock(dest, qaddr, qamount, (approxqp ? APPROX : EXACT));
2378	qblocksent += qamount;
2379	qaddr += qamount;
2380	qamount=sgss(&sched);
2381
2382	MPI_Iprobe(MPI_ANY_SOURCE, PP_QUARTBLOCKSPECS, PP_Comm, &flag, &stat);
2383	while (flag) {
2384	PP_RecvQuartBlock(0, &sq, &noq, quartetinfo, &apr);
2385	qblocksent -= noq;
2386	MPI_Iprobe(MPI_ANY_SOURCE, PP_QUARTBLOCKSPECS, PP_Comm, &flag, &stat);
2387	}
2388	}
2389	while (qblocksent > 0) {
2390	PP_RecvQuartBlock(0, &sq, &noq, quartetinfo, &apr);
2391	qblocksent -= noq;
2392	}
2393	}
2394	# else /* PARALLEL */
2395
2396	addtimes(GENERAL, &tarr);
2397	if (savequartlh_optn) {
2398	openfiletowrite(&lhfp, ALLQUARTLH, "all quartet likelihoods");
2399	if (saveqlhbin_optn) writetpqfheader(Maxspc, lhfp, 3);
2400	else writetpqfheader(Maxspc, lhfp, 4);
2401	}
2402
2403	for (i = 3; i < Maxspc; i++)
2404	for (c = 2; c < i; c++)
2405	for (b = 1; b < c; b++)
2406	for (a = 0; a < b; a++) {
2407	nq++;
2408
2409	/* generate message every 15 minutes */
2410	/* check timer */
2411	time(&time2);
2412	if ( (time2 - time1) > 900) {
2413	/* every 900 seconds */
2414	/* percentage of completed quartets */
2415	if (mflag == 0) {
2416	FPRINTF(STDOUTFILE "\n");
2417	mflag = 1;
2418	}
2419	qc2 = 100.*nq/Numquartets;
2420	mintogo = (100.0-qc2) *
2421	(double) (time2-time0)/60.0/qc2;
2422	hours = floor(mintogo/60.0);
2423	minutes = mintogo - 60.0*hours;
2424	FPRINTF(STDOUTFILE "%.2f%%", qc2);
2425	FPRINTF(STDOUTFILE " completed (remaining");
2426	FPRINTF(STDOUTFILE " time: %.0f", hours);
2427	FPRINTF(STDOUTFILE " hours %.0f", minutes);
2428	FPRINTF(STDOUTFILE " minutes)\n");
2429	fflush(STDOUT);
2430	time1 = time2;
2431	}
2432
2433	/* maximum likelihood values */
2434
2435	/* exact or approximate maximum likelihood values */
2436	compute_quartlklhds(a,b,c,i,&qweight[0],&qweight[1],&qweight[2], (approxqp ? APPROX : EXACT));
2437
2438	if (savequartlh_optn) {
2439	if (saveqlhbin_optn)
2440	fwrite(qweight, sizeof(double), 3, lhfp);
2441	else
2442	fprintf(lhfp, "(%d,%d,%d,%d)\t%f\t%f\t%f\n", a, b, c, i,
2443	qweight[0], qweight[1], qweight[2]);
2444	}
2445
2446	/* sort in descending order */
2447	sort3doubles(qweight, qworder);
2448
2449	if (usebestq_optn) {
2450	sqorder[2] = 2;
2451	discreteweight[sqorder[2]] = treebits[qworder[0]];
2452	if (qweight[qworder[0]] == qweight[qworder[1]]) {
2453	discreteweight[sqorder[2]] = discreteweight[sqorder[2]] \|\| treebits[qworder[1]];
2454	if (qweight[qworder[1]] == qweight[qworder[2]]) {
2455	discreteweight[sqorder[2]] = discreteweight[sqorder[2]] \|\| treebits[qworder[2]];
2456	discreteweight[sqorder[2]] = 7;
2457	}
2458	}
2459	} else {
2460
2461	/* compute Bayesian weights */
2462	qweight[qworder[1]] = exp(qweight[qworder[1]]-qweight[qworder[0]]);
2463	qweight[qworder[2]] = exp(qweight[qworder[2]]-qweight[qworder[0]]);
2464	qweight[qworder[0]] = 1.0;
2465	temp = qweight[0] + qweight[1] + qweight[2];
2466	qweight[0] = qweight[0]/temp;
2467	qweight[1] = qweight[1]/temp;
2468	qweight[2] = qweight[2]/temp;
2469
2470	/* square deviations */
2471	temp1 = 1.0 - qweight[qworder[0]];
2472	sqdiff[0] = temp1 * temp1 +
2473	qweight[qworder[1]] * qweight[qworder[1]] +
2474	qweight[qworder[2]] * qweight[qworder[2]];
2475	discreteweight[0] = treebits[qworder[0]];
2476
2477	temp1 = 0.5 - qweight[qworder[0]];
2478	temp2 = 0.5 - qweight[qworder[1]];
2479	sqdiff[1] = temp1 * temp1 + temp2 * temp2 +
2480	qweight[qworder[2]] * qweight[qworder[2]];
2481	discreteweight[1] = treebits[qworder[0]] + treebits[qworder[1]];
2482
2483	temp1 = onethird - qweight[qworder[0]];
2484	temp2 = onethird - qweight[qworder[1]];
2485	temp3 = onethird - qweight[qworder[2]];
2486	sqdiff[2] = temp1 * temp1 + temp2 * temp2 + temp3 * temp3;
2487	discreteweight[2] = (unsigned char) 7;
2488
2489	/* sort in descending order */
2490	sort3doubles(sqdiff, sqorder);
2491	}
2492
2493	/* determine best discrete weight */
2494	writequartet(a, b, c, i, discreteweight[sqorder[2]]);
2495
2496	/* counting completely unresolved quartets */
2497	if (discreteweight[sqorder[2]] == 7) {
2498	badqs++;
2499	badtaxon[a]++;
2500	badtaxon[b]++;
2501	badtaxon[c]++;
2502	badtaxon[i]++;
2503	if (show_optn) {
2504	fputid10(unresfp, a);
2505	fprintf(unresfp, " ");
2506	fputid10(unresfp, b);
2507	fprintf(unresfp, " ");
2508	fputid10(unresfp, c);
2509	fprintf(unresfp, " ");
2510	fputid(unresfp, i);
2511	fprintf(unresfp, "\n");
2512	}
2513	}
2514	addtimes(QUARTETS, &tarr);
2515	}
2516	if (savequartlh_optn) {
2517	closefile(lhfp);
2518	}
2519	if (show_optn)
2520	closefile(unresfp);
2521	if (mflag == 1)
2522	FPRINTF(STDOUTFILE "\n");
2523	# endif /* PARALLEL */
2524
2525	}
2526
2527	/* check the branching structure between the leaves (not the taxa!)
2528	A, B, C, and I (A, B, C, I don't need to be ordered). As a result,
2529	the two leaves that are closer related to each other than to leaf I
2530	are found in chooseA and chooseB. If the branching structure is
2531	not uniquely defined, ChooseA and ChooseB are chosen randomly
2532	from the possible taxa */
2533	void checkquartet(int A, int B, int C, int I)
2534	{
2535	int i, j, a, b, taxon[5], leaf[5], ipos;
2536	unsigned char qresult;
2537	int notunique = FALSE;
2538
2539	/* The relationship between leaves and taxa is defined by trueID */
2540	taxon[1] = trueID[A]; /* taxon number */
2541	leaf[1] = A; /* leaf number */
2542	taxon[2] = trueID[B];
2543	leaf[2] = B;
2544	taxon[3] = trueID[C];
2545	leaf[3] = C;
2546	taxon[4] = trueID[I];
2547	leaf[4] = I;
2548
2549	/* sort for taxa */
2550	/* Source: Numerical Recipes (PIKSR2.C) */
2551	for (j = 2; j <= 4; j++) {
2552	a = taxon[j];
2553	b = leaf[j];
2554	i = j-1;
2555	while (i > 0 && taxon[i] > a) {
2556	taxon[i+1] = taxon[i];
2557	leaf[i+1] = leaf[i];
2558	i--;
2559	}
2560	taxon[i+1] = a;
2561	leaf[i+1] = b;
2562	}
2563
2564	/* where is leaf I ? */
2565	ipos = 1;
2566	while (leaf[ipos] != I) ipos++;
2567
2568	/* look at sequence quartet */
2569	qresult = readquartet(taxon[1], taxon[2], taxon[3], taxon[4]);
2570
2571	/* chooseA and chooseB */
2572	do {
2573	switch (qresult) {
2574
2575	/* one single branching structure */
2576
2577	/* 001 */
2578	case 1: if (ipos == 1 \|\| ipos == 2) {
2579	chooseA = leaf[3];
2580	chooseB = leaf[4];
2581	} else {
2582	chooseA = leaf[1];
2583	chooseB = leaf[2];
2584	}
2585	notunique = FALSE;
2586	break;
2587
2588	/* 010 */
2589	case 2: if (ipos == 1 \|\| ipos == 3) {
2590	chooseA = leaf[2];
2591	chooseB = leaf[4];
2592	} else {
2593	chooseA = leaf[1];
2594	chooseB = leaf[3];
2595	}
2596	notunique = FALSE;
2597	break;
2598
2599	/* 100 */
2600	case 4: if (ipos == 1 \|\| ipos == 4) {
2601	chooseA = leaf[2];
2602	chooseB = leaf[3];
2603	} else {
2604	chooseA = leaf[1];
2605	chooseB = leaf[4];
2606	}
2607	notunique = FALSE;
2608	break;
2609
2610	/* two possible branching structures */
2611
2612	/* 011 */
2613	case 3: if (randominteger(2)) qresult = 1;
2614	else qresult = 2;
2615	notunique = TRUE;
2616	break;
2617
2618	/* 101 */
2619	case 5: if (randominteger(2)) qresult = 1;
2620	else qresult = 4;
2621	notunique = TRUE;
2622	break;
2623
2624	/* 110 */
2625	case 6: if (randominteger(2)) qresult = 2;
2626	else qresult = 4;
2627	notunique = TRUE;
2628	break;
2629
2630	/* three possible branching structures */
2631
2632	/* 111 */
2633	case 7: qresult = (1 << randominteger(3)); /* 1, 2, or 4 */
2634	notunique = TRUE;
2635	break;
2636
2637	default: /* Program error [checkquartet] */
2638	#if PARALLEL
2639	FPRINTF(STDOUTFILE "\n\n\n(%2d)HALT: PLEASE REPORT ERROR K-PARALLEL TO DEVELOPERS (%d,%d,%d,%d) = %ld\n\n\n",
2640	PP_Myid, taxon[1], taxon[2], taxon[3], taxon[4],
2641	quart2num(taxon[1], taxon[2], taxon[3], taxon[4]));
2642	#else
2643	FPRINTF(STDOUTFILE "\n\n\nHALT: PLEASE REPORT ERROR K TO DEVELOPERS\n\n\n");
2644	#endif
2645
2646	}
2647	} while (notunique);
2648
2649	return;
2650	}
2651

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

Download in other formats: