Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

source: tags/initial/GDE/PHYLIP/protpars.c

Visit:

Last change on this file was 2, checked in by oldcode, 24 years ago
Initial revision
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 58.7 KB

Line
1	#include "phylip.h"
2
3	/*
4	* version 3.52c. (c) Copyright 1993 by Joseph Felsenstein. Written by Joseph
5	* Felsenstein, Akiko Fuseki, Sean Lamont, and Andrew Keeffe. Permission is
6	* granted to copy and use this program provided no fee is charged for it and
7	* provided that this copyright notice is not removed.
8	*/
9
10	#define nmlngth 10 /* number of characters in species name */
11	#define maxtrees 50 /* maximum number of tied trees stored */
12	#define maxuser 8 /* maximum number of user-defined trees */
13	#define down 2
14
15	#define ibmpc0 false
16	#define ansi0 true
17	#define vt520 false
18
19
20	typedef long *steptr;
21	typedef enum {
22	ala, arg, asn, asp, cys, gln, glu, gly, his, ileu, leu, lys, met, phe, pro,
23	ser1, ser2, thr, trp, tyr, val, del, stop, asx, glx, ser, unk, quest
24	} aas;
25	typedef long sitearray[3];
26	typedef sitearray *seqptr;
27
28	/* nodes will form a binary tree */
29
30	typedef struct node { /* describes a tip species or an ancestor */
31	struct node next, back; /* pointers to nodes */
32	long index; /* number of the node */
33	boolean tip, bottom; /* present species are tips of tree */
34	aas seq; / the sequence */
35	seqptr siteset;/* temporary storage for aa's */
36	steptr numsteps; /* bookkeeps steps */
37	long xcoord, ycoord, ymin; /* used by printree */
38
39	long ymax;
40	} node;
41
42	typedef node **pointptr;
43	typedef long longer[6];
44
45	typedef struct gseq {
46	seqptr seq;
47	struct gseq *next;
48	} gseq;
49
50
51	Static node *root;
52	Static FILE infile, outfile, *treefile;
53	Static long spp, nonodes, chars, inseed, outgrno, col, datasets, ith, i,
54	j, l, jumb, njumble;
55	/*
56	* spp = number of species nonodes = number of nodes in tree chars = number
57	* of sites in actual sequences outgrno indicates outgroup
58	*/
59	Static boolean jumble, usertree, weights, thresh, trout, outgropt, printdata,
60	progress, treeprint, stepbox, ancseq, mulsets, interleaved,
61	ibmpc, vt52, ansi, firstset;
62	Static long fullset, fulldel;
63	Static steptr weight;
64	Static pointptr treenode; /* pointers to all nodes in tree */
65	Static Char *nayme; / names of species */
66	Static double threshold;
67	Static steptr threshwt;
68	Static longer seed;
69	Static long *enterorder;
70	Static sitearray translate[(long) quest - (long) ala + 1];
71	Static long **fsteps;
72	Static long **bestrees;
73	Static gseq *garbage;
74	Static node temp, temp1;
75	Char ch;
76	aas tmpa;
77
78	/* Local variables for maketree, propagated globally for c version: */
79	long nextree, which, minwhich;
80	double like, bestyet, bestlike, minsteps, bstlike2;
81	boolean lastrearr, recompute;
82	node *there;
83	double nsteps[maxuser];
84	long *place;
85	boolean *names;
86
87
88
89	openfile(fp, filename, mode, application, perm)
90	FILE **fp;
91	char *filename;
92	char *mode;
93	char *application;
94	char *perm;
95	{
96	FILE *of;
97	char file[100];
98	strcpy(file, filename);
99	while (1) {
100	of = fopen(file, mode);
101	if (of)
102	break;
103	else {
104	switch (*mode) {
105	case 'r':
106	printf("%s: can't read %s\n", application, file);
107	file[0] = '\0';
108	while (file[0] == '\0') {
109	printf("Please enter a new filename>");
110	gets(file);
111	}
112	break;
113	case 'w':
114	case 'a':
115	printf("%s: can't write %s\n", application, file);
116	file[0] = '\0';
117	while (file[0] == '\0') {
118	printf("Please enter a new filename>");
119	gets(file);
120	}
121	break;
122	}
123	}
124	}
125	*fp = of;
126	if (perm != NULL)
127	strcpy(perm, file);
128	}
129
130
131	void
132	gnu(p)
133	gseq **p;
134	{
135	/*
136	* this and the following are do-it-yourself garbage collectors. Make
137	* a new node or pull one off the garbage list
138	*/
139	if (garbage != NULL) {
140	*p = garbage;
141	garbage = garbage->next;
142	} else {
143	p = (gseq ) Malloc(sizeof(gseq));
144	(p)->seq = (seqptr) Malloc(chars sizeof(sitearray));
145	}
146	(*p)->next = NULL;
147	} /* gnu */
148
149
150	void
151	chuck(p)
152	gseq *p;
153	{
154	/* collect garbage on p -- put it on front of garbage list */
155	p->next = garbage;
156	garbage = p;
157	} /* chuck */
158
159
160	void
161	setup()
162	{
163	/* set up set table to get aasets from aas */
164	aas a, b;
165	long s;
166
167	for (a = ala; (long) a <= (long) stop; a = (aas) ((long) a + 1))
168	translate[(long) a - (long) ala][0] = 1L << ((long) a);
169	translate[0]
170	[1] = (1L << ((long) ala)) \| (1L << ((long) asp)) \| (1L << ((long) glu)) \|
171	(1L << ((long) gly)) \| (1L << ((long) ser1)) \| (1L << ((long) pro)) \|
172	(1L << ((long) thr)) \| (1L << ((long) val));
173	translate[(long) arg - (long) ala]
174	[1] = (1L << ((long) arg)) \| (1L << ((long) cys)) \| (1L << ((long) gln)) \|
175	(1L << ((long) gly)) \| (1L << ((long) his)) \| (1L << ((long) ileu)) \|
176	(1L << ((long) leu)) \| (1L << ((long) lys)) \| (1L << ((long) met)) \|
177	(1L << ((long) pro)) \| (1L << ((long) ser2)) \| (1L << ((long) thr)) \|
178	(1L << ((long) trp)) \| (1L << ((long) stop));
179	translate[(long) asn - (long) ala]
180	[1] = (1L << ((long) asn)) \| (1L << ((long) asp)) \| (1L << ((long) his)) \|
181	(1L << ((long) ileu)) \| (1L << ((long) lys)) \| (1L << ((long) ser2)) \|
182	(1L << ((long) thr)) \| (1L << ((long) tyr));
183	translate[(long) asp - (long) ala]
184	[1] = (1L << ((long) ala)) \| (1L << ((long) asp)) \| (1L << ((long) asn)) \|
185	(1L << ((long) glu)) \| (1L << ((long) gly)) \| (1L << ((long) his)) \|
186	(1L << ((long) tyr)) \| (1L << ((long) val));
187	translate[(long) cys - (long) ala]
188	[1] = (1L << ((long) arg)) \| (1L << ((long) cys)) \| (1L << ((long) gly)) \|
189	(1L << ((long) phe)) \| (1L << ((long) ser1)) \| (1L << ((long) ser2)) \|
190	(1L << ((long) trp)) \| (1L << ((long) tyr)) \| (1L << ((long) stop));
191	translate[(long) gln - (long) ala]
192	[1] = (1L << ((long) arg)) \| (1L << ((long) gln)) \| (1L << ((long) glu)) \|
193	(1L << ((long) his)) \| (1L << ((long) leu)) \| (1L << ((long) lys)) \|
194	(1L << ((long) pro)) \| (1L << ((long) stop));
195	translate[(long) glu - (long) ala]
196	[1] = (1L << ((long) ala)) \| (1L << ((long) asp)) \| (1L << ((long) gln)) \|
197	(1L << ((long) glu)) \| (1L << ((long) gly)) \| (1L << ((long) lys)) \|
198	(1L << ((long) val)) \| (1L << ((long) stop));
199	translate[(long) gly - (long) ala]
200	[1] = (1L << ((long) ala)) \| (1L << ((long) arg)) \| (1L << ((long) asp)) \|
201	(1L << ((long) cys)) \| (1L << ((long) glu)) \| (1L << ((long) gly)) \|
202	(1L << ((long) ser2)) \| (1L << ((long) trp)) \| (1L << ((long) val)) \|
203	(1L << ((long) stop));
204	translate[(long) his - (long) ala]
205	[1] = (1L << ((long) arg)) \| (1L << ((long) asn)) \| (1L << ((long) asp)) \|
206	(1L << ((long) gln)) \| (1L << ((long) his)) \| (1L << ((long) leu)) \|
207	(1L << ((long) pro)) \| (1L << ((long) tyr));
208	translate[(long) ileu - (long) ala]
209	[1] = (1L << ((long) arg)) \| (1L << ((long) asn)) \| (1L << ((long) ileu)) \|
210	(1L << ((long) leu)) \| (1L << ((long) lys)) \| (1L << ((long) met)) \|
211	(1L << ((long) phe)) \| (1L << ((long) ser2)) \| (1L << ((long) thr)) \|
212	(1L << ((long) val));
213	translate[(long) leu - (long) ala]
214	[1] = (1L << ((long) arg)) \| (1L << ((long) gln)) \| (1L << ((long) his)) \|
215	(1L << ((long) ileu)) \| (1L << ((long) leu)) \| (1L << ((long) met)) \|
216	(1L << ((long) phe)) \| (1L << ((long) pro)) \| (1L << ((long) ser1)) \|
217	(1L << ((long) trp)) \| (1L << ((long) val)) \| (1L << ((long) stop));
218	translate[(long) lys - (long) ala]
219	[1] = (1L << ((long) arg)) \| (1L << ((long) asn)) \| (1L << ((long) gln)) \|
220	(1L << ((long) glu)) \| (1L << ((long) ileu)) \| (1L << ((long) lys)) \|
221	(1L << ((long) met)) \| (1L << ((long) thr)) \| (1L << ((long) stop));
222	translate[(long) met - (long) ala]
223	[1] = (1L << ((long) arg)) \| (1L << ((long) ileu)) \| (1L << ((long) leu)) \|
224	(1L << ((long) lys)) \| (1L << ((long) met)) \| (1L << ((long) val)) \|
225	(1L << ((long) thr));
226	translate[(long) phe - (long) ala]
227	[1] = (1L << ((long) cys)) \| (1L << ((long) ileu)) \| (1L << ((long) leu)) \|
228	(1L << ((long) phe)) \| (1L << ((long) ser1)) \| (1L << ((long) tyr)) \|
229	(1L << ((long) val));
230	translate[(long) pro - (long) ala]
231	[1] = (1L << ((long) ala)) \| (1L << ((long) arg)) \| (1L << ((long) gln)) \|
232	(1L << ((long) his)) \| (1L << ((long) leu)) \| (1L << ((long) pro)) \|
233	(1L << ((long) ser1)) \| (1L << ((long) thr));
234	translate[(long) ser1 - (long) ala]
235	[1] = (1L << ((long) ala)) \| (1L << ((long) cys)) \| (1L << ((long) leu)) \|
236	(1L << ((long) phe)) \| (1L << ((long) pro)) \| (1L << ((long) ser1)) \|
237	(1L << ((long) thr)) \| (1L << ((long) trp)) \| (1L << ((long) tyr)) \|
238	(1L << ((long) stop));
239	translate[(long) ser2 - (long) ala]
240	[1] = (1L << ((long) arg)) \| (1L << ((long) asn)) \| (1L << ((long) cys)) \|
241	(1L << ((long) gly)) \| (1L << ((long) ileu)) \| (1L << ((long) ser2)) \|
242	(1L << ((long) thr));
243	translate[(long) thr - (long) ala]
244	[1] = (1L << ((long) ala)) \| (1L << ((long) arg)) \| (1L << ((long) asn)) \|
245	(1L << ((long) ileu)) \| (1L << ((long) lys)) \| (1L << ((long) met)) \|
246	(1L << ((long) pro)) \| (1L << ((long) ser1)) \| (1L << ((long) ser2)) \|
247	(1L << ((long) thr));
248	translate[(long) trp - (long) ala]
249	[1] = (1L << ((long) arg)) \| (1L << ((long) cys)) \| (1L << ((long) gly)) \|
250	(1L << ((long) leu)) \| (1L << ((long) ser1)) \| (1L << ((long) stop)) \|
251	(1L << ((long) trp));
252	translate[(long) tyr - (long) ala]
253	[1] = (1L << ((long) asn)) \| (1L << ((long) asp)) \| (1L << ((long) cys)) \|
254	(1L << ((long) his)) \| (1L << ((long) phe)) \| (1L << ((long) ser1)) \|
255	(1L << ((long) stop)) \| (1L << ((long) tyr));
256	translate[(long) val - (long) ala]
257	[1] = (1L << ((long) ala)) \| (1L << ((long) asp)) \| (1L << ((long) glu)) \|
258	(1L << ((long) gly)) \| (1L << ((long) ileu)) \| (1L << ((long) leu)) \|
259	(1L << ((long) met)) \| (1L << ((long) phe)) \| (1L << ((long) val));
260	translate[(long) stop - (long) ala]
261	[1] = (1L << ((long) arg)) \| (1L << ((long) cys)) \| (1L << ((long) gln)) \|
262	(1L << ((long) glu)) \| (1L << ((long) gly)) \| (1L << ((long) leu)) \|
263	(1L << ((long) lys)) \| (1L << ((long) ser1)) \| (1L << ((long) trp)) \|
264	(1L << ((long) tyr)) \| (1L << ((long) stop));
265	translate[(long) del - (long) ala][1] = 1L << ((long) del);
266	fulldel = (1L << ((long) stop + 1)) - (1L << ((long) ala));
267	fullset = fulldel & (~(1L << ((long) del)));
268	translate[(long) asx - (long) ala]
269	[0] = (1L << ((long) asn)) \| (1L << ((long) asp));
270	translate[(long) glx - (long) ala]
271	[0] = (1L << ((long) gln)) \| (1L << ((long) glu));
272	translate[(long) ser - (long) ala]
273	[0] = (1L << ((long) ser1)) \| (1L << ((long) ser2));
274	translate[(long) unk - (long) ala][0] = fullset;
275	translate[(long) quest - (long) ala][0] = fulldel;
276	translate[(long) asx - (long) ala]
277	[1] = translate[(long) asn - (long) ala]
278	[1] \| translate[(long) asp - (long) ala][1];
279	translate[(long) glx - (long) ala]
280	[1] = translate[(long) gln - (long) ala]
281	[1] \| translate[(long) glu - (long) ala][1];
282	translate[(long) ser - (long) ala]
283	[1] = translate[(long) ser1 - (long) ala]
284	[1] \| translate[(long) ser2 - (long) ala][1];
285	translate[(long) unk - (long) ala][1] = fullset;
286	translate[(long) quest - (long) ala][1] = fulldel;
287	for (a = ala; (long) a <= (long) quest; a = (aas) ((long) a + 1)) {
288	s = 0;
289	for (b = ala; (long) b <= (long) stop; b = (aas) ((long) b + 1)) {
290	if (((1L << ((long) b)) & translate[(long) a - (long) ala][1]) != 0)
291	s \|= translate[(long) b - (long) ala][1];
292	}
293	translate[(long) a - (long) ala][2] = s;
294	}
295	} /* setup */
296
297
298	double
299	randum(seed)
300	long *seed;
301	{
302	/* random number generator -- slow but machine independent */
303	long i, j, k, sum;
304	longer mult, newseed;
305	double x;
306
307	mult[0] = 13;
308	mult[1] = 24;
309	mult[2] = 22;
310	mult[3] = 6;
311	for (i = 0; i <= 5; i++)
312	newseed[i] = 0;
313	for (i = 0; i <= 5; i++) {
314	sum = newseed[i];
315	k = i;
316	if (i > 3)
317	k = 3;
318	for (j = 0; j <= k; j++)
319	sum += mult[j] * seed[i - j];
320	newseed[i] = sum;
321	for (j = i; j <= 4; j++) {
322	newseed[j + 1] += newseed[j] / 64;
323	newseed[j] &= 63;
324	}
325	}
326	memcpy(seed, newseed, sizeof(longer));
327	seed[5] &= 3;
328	x = 0.0;
329	for (i = 0; i <= 5; i++)
330	x = x / 64.0 + seed[i];
331	x /= 4.0;
332	return x;
333	} /* randum */
334
335
336
337	void
338	uppercase(ch)
339	Char *ch;
340	{
341	/* convert ch to upper case -- either ASCII or EBCDIC */
342	ch = (islower(ch) ? toupper(ch) : (ch));
343	} /* uppercase */
344
345
346	void
347	inputnumbers()
348	{
349	/* input the numbers of species and of characters */
350	fscanf(infile, "%ld%ld", &spp, &chars);
351	if (printdata)
352	fprintf(outfile, "%2ld species, %3ld sites\n", spp, chars);
353	if (printdata)
354	putc('\n', outfile);
355	nonodes = spp * 2 - 1;
356	} /* inputnumbers */
357
358	void
359	getoptions()
360	{
361	/* interactively set options */
362	long i, inseed0;
363	Char ch;
364	boolean done, done1;
365
366	fprintf(outfile, "\nProtein parsimony algorithm, version %s\n\n", VERSION);
367	putchar('\n');
368	jumble = false;
369	njumble = 1;
370	outgrno = 1;
371	outgropt = false;
372	thresh = false;
373	trout = true;
374	usertree = false;
375	weights = false;
376	printdata = false;
377	progress = true;
378	treeprint = true;
379	stepbox = false;
380	ancseq = false;
381	interleaved = true;
382	for (;;) {
383	printf(ansi ? "\033[2J\033[H" : vt52 ? "\033E\033H" : "\n");
384	printf("\nProtein parsimony algorithm, version %s\n\n", VERSION);
385	printf("Setting for this run:\n");
386	printf(" U Search for best tree? %s\n",
387	(usertree ? "No, use user trees in input file" : "Yes"));
388	if (!usertree) {
389	printf(" J Randomize input order of sequences?");
390	if (jumble)
391	printf(" Yes (seed =%8ld,%3ld times)\n", inseed0, njumble);
392	else
393	printf(" No. Use input order\n");
394	}
395	printf(" O Outgroup root?");
396	if (outgropt)
397	printf(" Yes, at sequence number%3ld\n", outgrno);
398	else
399	printf(" No, use as outgroup species%3ld\n", outgrno);
400	printf(" T Use Threshold parsimony?");
401	if (thresh)
402	printf(" Yes, count steps up to%4.1f per site\n", threshold);
403	else
404	printf(" No, use ordinary parsimony\n");
405	printf(" M Analyze multiple data sets?");
406	if (mulsets)
407	printf(" Yes, %2ld sets\n", datasets);
408	else
409	printf(" No\n");
410	printf(" I Input sequences interleaved? %s\n",
411	(interleaved ? "Yes" : "No"));
412	printf(" 0 Terminal type (IBM PC, VT52, ANSI)? %s\n",
413	(ibmpc ? "IBM PC" : ansi ? "ANSI" : vt52 ? "VT52" : "(none)"));
414	printf(" 1 Print out the data at start of run %s\n",
415	(printdata ? "Yes" : "No"));
416	printf(" 2 Print indications of progress of run %s\n",
417	(progress ? "Yes" : "No"));
418	printf(" 3 Print out tree %s\n",
419	(treeprint ? "Yes" : "No"));
420	printf(" 4 Print out steps in each site %s\n",
421	(stepbox ? "Yes" : "No"));
422	printf(" 5 Print sequences at all nodes of tree %s\n",
423	(ancseq ? "Yes" : "No"));
424	printf(" 6 Write out trees onto tree file? %s\n",
425	(trout ? "Yes" : "No"));
426	printf(
427	"\nAre these settings correct? (type Y or the letter for one to change)\n");
428	scanf("%c%*[^\n]", &ch);
429	getchar();
430	uppercase(&ch);
431	if (ch == 'Y')
432	break;
433	if (strchr("JOTUMI1234560", ch)) {
434	switch (ch) {
435
436	case 'J':
437	jumble = !jumble;
438	if (jumble) {
439	printf("Random number seed (must be odd)?\n");
440	scanf("%ld%*[^\n]", &inseed);
441	getchar();
442	inseed0 = inseed;
443	for (i = 0; i <= 5; i++)
444	seed[i] = 0;
445	i = 0;
446	do {
447	seed[i] = inseed & 63;
448	inseed /= 64;
449	i++;
450	} while (inseed != 0);
451	printf("Number of times to jumble?\n");
452	scanf("%ld%*[^\n]", &njumble);
453	getchar();
454	} else
455	njumble = 1;
456	break;
457
458	case 'O':
459	outgropt = !outgropt;
460	if (outgropt) {
461	done1 = true;
462	do {
463	printf("Type number of the outgroup:\n");
464	scanf("%ld%*[^\n]", &outgrno);
465	getchar();
466	done1 = (outgrno >= 1 && outgrno <= spp);
467	if (!done1) {
468	printf("BAD OUTGROUP NUMBER: %4ld\n", outgrno);
469	printf(" Must be in range 1 -%2ld\n", spp);
470	}
471	} while (done1 != true);
472	} else
473	outgrno = 1;
474	break;
475
476	case 'T':
477	thresh = !thresh;
478	if (thresh) {
479	done1 = false;
480	do {
481	printf("What will be the threshold value?\n");
482	scanf("%lf%*[^\n]", &threshold);
483	getchar();
484	done1 = (threshold >= 1.0);
485	if (!done1)
486	printf("BAD THRESHOLD VALUE: it must be greater than 1\n");
487	else
488	threshold = (long) (threshold * 10.0 + 0.5) / 10.0;
489	} while (done1 != true);
490	}
491	break;
492
493	case 'M':
494	mulsets = !mulsets;
495	if (mulsets) {
496	done1 = false;
497	do {
498	printf("How many data sets?\n");
499	scanf("%ld%*[^\n]", &datasets);
500	getchar();
501	done1 = (datasets >= 1);
502	if (!done1)
503	printf("BAD DATA SETS NUMBER: it must be greater than 1\n");
504	} while (done1 != true);
505	}
506	break;
507
508	case 'I':
509	interleaved = !interleaved;
510	break;
511
512	case 'U':
513	usertree = !usertree;
514	break;
515
516	case '0':
517	if (ibmpc) {
518	ibmpc = false;
519	vt52 = true;
520	} else {
521	if (vt52) {
522	vt52 = false;
523	ansi = true;
524	} else if (ansi)
525	ansi = false;
526	else
527	ibmpc = true;
528	}
529	break;
530
531	case '1':
532	printdata = !printdata;
533	break;
534
535	case '2':
536	progress = !progress;
537	break;
538
539	case '3':
540	treeprint = !treeprint;
541	break;
542
543	case '4':
544	stepbox = !stepbox;
545	break;
546
547	case '5':
548	ancseq = !ancseq;
549	break;
550
551	case '6':
552	trout = !trout;
553	break;
554	}
555	} else
556	printf("Not a possible option!\n");
557	}
558	} /* getoptions */
559
560
561	void
562	doinit()
563	{
564	/* initializes variables */
565	long i;
566	node p, q;
567
568	inputnumbers();
569	getoptions();
570	treenode = (pointptr) Malloc(nonodes * sizeof(node *));
571	for (i = 0; i < (spp); i++) {
572	treenode[i] = (node *) Malloc(sizeof(node));
573	treenode[i]->numsteps = (steptr) Malloc(chars * sizeof(long));
574	treenode[i]->siteset = (seqptr) Malloc(chars * sizeof(sitearray));
575	treenode[i]->seq = (aas ) Malloc(chars sizeof(aas));
576	}
577	for (i = spp; i < (nonodes); i++) {
578	q = NULL;
579	for (j = 1; j <= 3; j++) {
580	p = (node *) Malloc(sizeof(node));
581	p->numsteps = (steptr) Malloc(chars * sizeof(long));
582	p->siteset = (seqptr) Malloc(chars * sizeof(sitearray));
583	p->seq = (aas ) Malloc(chars sizeof(aas));
584	p->next = q;
585	q = p;
586	}
587	p->next->next->next = p;
588	treenode[i] = p;
589	}
590	} /* doinit */
591
592	void
593	inputweights()
594	{
595	/* input the character weights, 0-9 and A-Z for weights 10 - 35 */
596	Char ch;
597	long i;
598
599	for (i = 1; i < nmlngth; i++)
600	ch = getc(infile);
601	for (i = 0; i < (chars); i++) {
602	do {
603	if (eoln(infile)) {
604	fscanf(infile, "%*[^\n]");
605	getc(infile);
606	}
607	ch = getc(infile);
608	} while (ch == ' ');
609	weight[i] = 1;
610	if (isdigit(ch))
611	weight[i] = ch - '0';
612	else if (isalpha(ch)) {
613	uppercase(&ch);
614	if (ch >= 'A' && ch <= 'I')
615	weight[i] = ch - 55;
616	else if (ch >= 'J' && ch <= 'R')
617	weight[i] = ch - 55;
618	else
619	weight[i] = ch - 55;
620	} else {
621	printf("BAD WEIGHT CHARACTER: %c\n", ch);
622	exit(-1);
623	}
624	}
625	fscanf(infile, "%*[^\n]");
626	getc(infile);
627	weights = true;
628	} /* inputweights */
629
630	void
631	printweights()
632	{
633	/* print out the weights of sites */
634	long i, j, k;
635
636	fprintf(outfile, " Sites are weighted as follows:\n");
637	fprintf(outfile, " ");
638	for (i = 0; i <= 9; i++)
639	fprintf(outfile, "%3ld", i);
640	fprintf(outfile, "\n *---------------------------------\n");
641	for (j = 0; j <= (chars / 10); j++) {
642	fprintf(outfile, "%5ld! ", j * 10);
643	for (i = 0; i <= 9; i++) {
644	k = j * 10 + i;
645	if (k > 0 && k <= chars)
646	fprintf(outfile, "%3ld", weight[k - 1]);
647	else
648	fprintf(outfile, " ");
649	}
650	putc('\n', outfile);
651	}
652	putc('\n', outfile);
653	} /* printweights */
654
655	void
656	inputoptions()
657	{
658	/* input the information on the options */
659	Char ch;
660	long extranum, i, cursp, curchs, FORLIM;
661
662	if (!firstset) {
663	if (eoln(infile)) {
664	fscanf(infile, "%*[^\n]");
665	getc(infile);
666	}
667	fscanf(infile, "%ld%ld", &cursp, &curchs);
668	if (cursp != spp) {
669	printf("\nERROR: INCONSISTENT NUMBER OF SPECIES IN DATA SET %4ld\n", ith);
670	exit(-1);
671	}
672	chars = curchs;
673	}
674	extranum = 0;
675	while (!(eoln(infile))) {
676	ch = getc(infile);
677	uppercase(&ch);
678	if (ch == 'W')
679	extranum++;
680	else if (ch != ' ') {
681	printf("BAD OPTION CHARACTER: %c\n", ch);
682	exit(-1);
683	}
684	}
685	fscanf(infile, "%*[^\n]");
686	getc(infile);
687	for (i = 0; i < (chars); i++)
688	weight[i] = 1;
689	for (i = 1; i <= extranum; i++) {
690	ch = getc(infile);
691	uppercase(&ch);
692	if (ch != 'W') {
693	printf("ERROR: INCORRECT AUXILIARY OPTIONS LINE");
694	printf(" WHICH STARTS WITH %c\n", ch);
695	exit(-1);
696	}
697	if (ch == 'W')
698	inputweights();
699	}
700	if (weights)
701	printweights();
702	if (!thresh)
703	threshold = spp * 3.0;
704	for (i = 0; i < (chars); i++) {
705	weight[i] *= 10;
706	threshwt[i] = (long) (threshold * weight[i] + 0.5);
707	}
708	} /* inputoptions */
709
710	void
711	inputdata()
712	{
713	/* input the names and sequences for each species */
714	long i, j, k, l, aasread, aasnew;
715	Char charstate;
716	boolean allread, done;
717	aas aa; /* temporary amino acid for input */
718	long FORLIM, FORLIM1;
719
720	if (progress)
721	putchar('\n');
722	j = nmlngth + (chars + (chars - 1) / 10) / 2 - 5;
723	if (j < nmlngth - 1)
724	j = nmlngth - 1;
725	if (j > 37)
726	j = 37;
727	if (printdata) {
728	fprintf(outfile, "Name");
729	for (i = 1; i <= j; i++)
730	putc(' ', outfile);
731	fprintf(outfile, "Sequences\n");
732	fprintf(outfile, "----");
733	for (i = 1; i <= j; i++)
734	putc(' ', outfile);
735	fprintf(outfile, "---------\n\n");
736	}
737	aasread = 0;
738	allread = false;
739	while (!(allread)) {
740	allread = true;
741	if (eoln(infile)) {
742	fscanf(infile, "%*[^\n]");
743	getc(infile);
744	}
745	i = 1;
746	while (i <= spp) {
747	if ((interleaved && aasread == 0) \|\| !interleaved) {
748	for (j = 0; j < nmlngth; j++) {
749	if (eof(infile) \|\| eoln(infile)) {
750	printf("ERROR: END-OF-LINE OR END-OF-FILE");
751	printf(" IN THE MIDDLE OF A SPECIES NAME\n");
752	exit(-1);
753	}
754	nayme[i - 1][j] = getc(infile);
755	}
756	}
757	j = interleaved ? aasread : 0;
758	done = false;
759	while (!done && !eof(infile)) {
760	if (interleaved)
761	done = true;
762	while (j < chars && !(eoln(infile) \|\| eof(infile))) {
763	charstate = getc(infile);
764	if (charstate == ' ' \|\| (charstate >= '0' && charstate <= '9'))
765	continue;
766	uppercase(&charstate);
767	if ((!isalpha(charstate) && charstate != '.' && charstate != '?' &&
768	charstate != '-' && charstate != '*') \|\| charstate == 'J' \|\|
769	charstate == 'O' \|\| charstate == 'U') {
770	printf("WARNING -- BAD AMINO ACID:%c", charstate);
771	printf(" AT POSITION%5ld OF SPECIES %3ld\n", j, i);
772	exit(-1);
773	}
774	j++;
775	if (charstate == '.') {
776	treenode[i - 1]->seq[j - 1] = treenode[0]->seq[j - 1];
777	memcpy(treenode[i - 1]->siteset[j - 1],
778	treenode[0]->siteset[j - 1], sizeof(sitearray));
779	continue;
780	}
781	aa = (charstate == 'A') ? ala :
782	(charstate == 'B') ? asx :
783	(charstate == 'C') ? cys :
784	(charstate == 'D') ? asp :
785	(charstate == 'E') ? glu :
786	(charstate == 'F') ? phe :
787	(charstate == 'G') ? gly : aa;
788	aa = (charstate == 'H') ? his :
789	(charstate == 'I') ? ileu :
790	(charstate == 'K') ? lys :
791	(charstate == 'L') ? leu :
792	(charstate == 'M') ? met :
793	(charstate == 'N') ? asn :
794	(charstate == 'P') ? pro :
795	(charstate == 'Q') ? gln :
796	(charstate == 'R') ? arg : aa;
797	aa = (charstate == 'S') ? ser :
798	(charstate == 'T') ? thr :
799	(charstate == 'V') ? val :
800	(charstate == 'W') ? trp :
801	(charstate == 'X') ? unk :
802	(charstate == 'Y') ? tyr :
803	(charstate == 'Z') ? glx :
804	(charstate == '*') ? stop :
805	(charstate == '?') ? quest :
806	(charstate == '-') ? del : aa;
807
808	treenode[i - 1]->seq[j - 1] = aa;
809	memcpy(treenode[i - 1]->siteset[j - 1],
810	translate[(long) aa - (long) ala], sizeof(sitearray));
811	}
812	if (interleaved)
813	continue;
814	if (j < chars) {
815	fscanf(infile, "%*[^\n]");
816	getc(infile);
817	} else if (j == chars)
818	done = true;
819	}
820	if (interleaved && i == 1)
821	aasnew = j;
822	fscanf(infile, "%*[^\n]");
823	getc(infile);
824	if ((interleaved && j != aasnew) \|\| ((!interleaved) && j != chars)) {
825	printf("ERROR: SEQUENCES OUT OF ALIGNMENT\n");
826	exit(-1);
827	}
828	i++;
829	}
830	if (interleaved) {
831	aasread = aasnew;
832	allread = (aasread == chars);
833	} else
834	allread = (i > spp);
835	}
836	if (printdata) {
837	for (i = 1; i <= ((chars - 1) / 60 + 1); i++) {
838	for (j = 1; j <= (spp); j++) {
839	for (k = 0; k < nmlngth; k++)
840	putc(nayme[j - 1][k], outfile);
841	fprintf(outfile, " ");
842	l = i * 60;
843	if (l > chars)
844	l = chars;
845	for (k = (i - 1) * 60 + 1; k <= l; k++) {
846	if (j > 1 && treenode[j - 1]->seq[k - 1] == treenode[0]->seq[k - 1])
847	charstate = '.';
848	else {
849	tmpa = treenode[j - 1]->seq[k - 1];
850	charstate = (tmpa == ala) ? 'A' :
851	(tmpa == asx) ? 'B' :
852	(tmpa == cys) ? 'C' :
853	(tmpa == asp) ? 'D' :
854	(tmpa == glu) ? 'E' :
855	(tmpa == phe) ? 'F' :
856	(tmpa == gly) ? 'G' :
857	(tmpa == his) ? 'H' :
858	(tmpa == ileu) ? 'I' :
859	(tmpa == lys) ? 'K' :
860	(tmpa == leu) ? 'L' : charstate;
861	charstate = (tmpa == met) ? 'M' :
862	(tmpa == asn) ? 'N' :
863	(tmpa == pro) ? 'P' :
864	(tmpa == gln) ? 'Q' :
865	(tmpa == arg) ? 'R' :
866	(tmpa == ser) ? 'S' :
867	(tmpa == ser1) ? 'S' :
868	(tmpa == ser2) ? 'S' : charstate;
869	charstate = (tmpa == thr) ? 'T' :
870	(tmpa == val) ? 'V' :
871	(tmpa == trp) ? 'W' :
872	(tmpa == unk) ? 'X' :
873	(tmpa == tyr) ? 'Y' :
874	(tmpa == glx) ? 'Z' :
875	(tmpa == del) ? '-' :
876	(tmpa == stop) ? '*' :
877	(tmpa == quest) ? '?' : charstate;
878	}
879	putc(charstate, outfile);
880	if (k % 10 == 0 && k % 60 != 0)
881	putc(' ', outfile);
882	}
883	putc('\n', outfile);
884	}
885	putc('\n', outfile);
886	}
887	putc('\n', outfile);
888	}
889	putc('\n', outfile);
890	} /* inputdata */
891
892	void
893	makevalues()
894	{
895	/* set up fractional likelihoods at tips */
896	long i, j;
897	node *p;
898
899	for (i = 1; i <= nonodes; i++) {
900	treenode[i - 1]->back = NULL;
901	treenode[i - 1]->tip = (i <= spp);
902	treenode[i - 1]->index = i;
903	for (j = 0; j < (chars); j++)
904	treenode[i - 1]->numsteps[j] = 0;
905	if (i > spp) {
906	p = treenode[i - 1]->next;
907	while (p != treenode[i - 1]) {
908	p->back = NULL;
909	p->tip = false;
910	p->index = i;
911	for (j = 0; j < (chars); j++)
912	p->numsteps[j] = 0;
913	p = p->next;
914	}
915	}
916	}
917	} /* makevalues */
918
919	void
920	doinput()
921	{
922	/* reads the input data */
923	inputoptions();
924	inputdata();
925	makevalues();
926	} /* doinput */
927
928
929
930	void
931	fillin(p, left, rt)
932	node p, left, *rt;
933	{
934	/*
935	* sets up for each node in the tree the aa set for site m at that
936	* point and counts the changes. The program spends much of its time
937	* in this PROCEDURE
938	*/
939	boolean counted;
940	aas aa;
941	long s;
942	sitearray ls, rs, qs;
943	long i, j, k, m, n;
944
945	for (m = 0; m < chars; m++) {
946	k = 0;
947	if (left != NULL)
948	memcpy(ls, left->siteset[m], sizeof(sitearray));
949	if (rt != NULL)
950	memcpy(rs, rt->siteset[m], sizeof(sitearray));
951	if (left == NULL) {
952	n = rt->numsteps[m];
953	memcpy(qs, rs, sizeof(sitearray));
954	} else if (rt == NULL) {
955	n = left->numsteps[m];
956	memcpy(qs, ls, sizeof(sitearray));
957	} else {
958	n = left->numsteps[m] + rt->numsteps[m];
959	counted = false;
960	for (i = 0; i <= 5; i++) {
961	if (k < 3) {
962	switch (i) {
963
964	case 0:
965	s = ls[0] & rs[0];
966	break;
967
968	case 1:
969	s = (ls[0] & rs[1]) \| (ls[1] & rs[0]);
970	break;
971
972	case 2:
973	s = (ls[0] & rs[2]) \| (ls[1] & rs[1]) \| (ls[2] & rs[0]);
974	break;
975
976	case 3:
977	s = ls[0] \| (ls[1] & rs[2]) \| (ls[2] & rs[1]) \| rs[0];
978	break;
979
980	case 4:
981	s = ls[1] \| (ls[2] & rs[2]) \| rs[1];
982	break;
983
984	case 5:
985	s = ls[2] \| rs[2];
986	break;
987	}
988	if (counted \|\| s != 0) {
989	qs[k] = s;
990	k++;
991	counted = true;
992	} else if (!counted)
993	n += weight[m];
994	}
995	}
996	}
997	for (i = 0; i <= 1; i++) {
998	for (aa = ala; (long) aa <= (long) stop; aa = (aas) ((long) aa + 1)) {
999	if (((1L << ((long) aa)) & qs[i]) != 0) {
1000	for (j = i + 1; j <= 2; j++) {
1001	#if 0
1002	if (translate[(long) aa - (long) ala][j - i] & ~qs[j])
1003	printf("%i %i %x %x\n", i, j, qs[j], translate[(long) aa - (long) ala][j - i]);
1004	#endif
1005	qs[j] \|= translate[(long) aa - (long) ala][j - i];
1006	}
1007	}
1008	}
1009	}
1010	p->numsteps[m] = n;
1011	memcpy(p->siteset[m], qs, sizeof(sitearray));
1012	}
1013	} /* fillin */
1014
1015	void
1016	preorder(p)
1017	node *p;
1018	{
1019	/*
1020	* recompute number of steps in preorder taking both ancestoral and
1021	* descendent steps into account
1022	*/
1023	if (p != NULL && !p->tip) {
1024	fillin(p->next, p->next->next->back, p->back);
1025	fillin(p->next->next, p->back, p->next->back);
1026	preorder(p->next->back);
1027	preorder(p->next->next->back);
1028	}
1029	} /* preorder */
1030
1031
1032	void
1033	add(below, newtip, newfork)
1034	node below, newtip, *newfork;
1035	{
1036	/*
1037	* inserts the nodes newfork and its left descendant, newtip, to the
1038	* tree. below becomes newfork's right descendant
1039	*/
1040	if (below != treenode[below->index - 1])
1041	below = treenode[below->index - 1];
1042	if (below->back != NULL)
1043	below->back->back = newfork;
1044	newfork->back = below->back;
1045	below->back = newfork->next->next;
1046	newfork->next->next->back = below;
1047	newfork->next->back = newtip;
1048	newtip->back = newfork->next;
1049	if (root == below)
1050	root = newfork;
1051	root->back = NULL;
1052	if (recompute) {
1053	fillin(newfork, newfork->next->back, newfork->next->next->back);
1054	preorder(newfork);
1055	if (newfork != root)
1056	preorder(newfork->back);
1057	}
1058	} /* add */
1059
1060	void
1061	re_move(item, fork)
1062	node item, fork;
1063	{
1064	/*
1065	* removes nodes item and its ancestor, fork, from the tree. the new
1066	* descendant of fork's ancestor is made to be fork's second
1067	* descendant (other than item). Also returns pointers to the
1068	* deleted nodes, item and fork
1069	*/
1070	node p, q, *other;
1071
1072	if ((*item)->back == NULL) {
1073	*fork = NULL;
1074	return;
1075	}
1076	fork = treenode[(item)->back->index - 1];
1077	if ((item) == (fork)->next->back)
1078	other = (*fork)->next->next->back;
1079	else
1080	other = (*fork)->next->back;
1081	if (root == *fork)
1082	root = other;
1083	p = (*item)->back->next->back;
1084	q = (*item)->back->next->next->back;
1085	if (p != NULL)
1086	p->back = q;
1087	if (q != NULL)
1088	q->back = p;
1089	(*fork)->back = NULL;
1090	p = (*fork)->next;
1091	do {
1092	p->back = NULL;
1093	p = p->next;
1094	} while (p != (*fork));
1095	(*item)->back = NULL;
1096	if (recompute) {
1097	preorder(other);
1098	if (other != root)
1099	preorder(other->back);
1100	}
1101	} /* re_move */
1102
1103	void
1104	evaluate(r)
1105	node *r;
1106	{
1107	/*
1108	* determines the number of steps needed for a tree. this is the
1109	* minimum number of steps needed to evolve sequences on this tree
1110	*/
1111	long i, steps, term;
1112	double sum;
1113
1114	sum = 0.0;
1115	for (i = 0; i < (chars); i++) {
1116	steps = r->numsteps[i];
1117	if (steps <= threshwt[i])
1118	term = steps;
1119	else
1120	term = threshwt[i];
1121	sum += term;
1122	if (usertree && which <= maxuser)
1123	fsteps[which - 1][i] = term;
1124	}
1125	if (usertree && which <= maxuser) {
1126	nsteps[which - 1] = sum;
1127	if (which == 1) {
1128	minwhich = 1;
1129	minsteps = sum;
1130	} else if (sum < minsteps) {
1131	minwhich = which;
1132	minsteps = sum;
1133	}
1134	}
1135	like = -sum;
1136	} /* evaluate */
1137
1138	void
1139	postorder(p)
1140	node *p;
1141	{
1142	/*
1143	* traverses a binary tree, calling PROCEDURE fillin at a node's
1144	* descendants before calling fillin at the node
1145	*/
1146	if (p->tip)
1147	return;
1148	postorder(p->next->back);
1149	postorder(p->next->next->back);
1150	fillin(p, p->next->back, p->next->next->back);
1151	} /* postorder */
1152
1153	void
1154	reroot(outgroup)
1155	node *outgroup;
1156	{
1157	/* reorients tree, putting outgroup in desired position. */
1158	node p, q;
1159
1160	if (outgroup->back->index == root->index)
1161	return;
1162	p = root->next;
1163	q = root->next->next;
1164	p->back->back = q->back;
1165	q->back->back = p->back;
1166	p->back = outgroup;
1167	q->back = outgroup->back;
1168	outgroup->back->back = q;
1169	outgroup->back = p;
1170	} /* reroot */
1171
1172
1173
1174	void
1175	savetraverse(p, pos, found)
1176	node *p;
1177	long *pos;
1178	boolean *found;
1179	{
1180	/* sets BOOLEANs that indicate which way is down */
1181	p->bottom = true;
1182	if (p->tip)
1183	return;
1184	p->next->bottom = false;
1185	savetraverse(p->next->back, pos, found);
1186	p->next->next->bottom = false;
1187	savetraverse(p->next->next->back, pos, found);
1188	} /* savetraverse */
1189
1190	void
1191	savetree(pos, found)
1192	long *pos;
1193	boolean *found;
1194	{
1195	/*
1196	* record in place where each species has to be added to reconstruct
1197	* this tree
1198	*/
1199	long i, j;
1200	node *p;
1201	boolean done;
1202
1203	reroot(treenode[outgrno - 1]);
1204	savetraverse(root, pos, found);
1205	for (i = 0; i < (nonodes); i++)
1206	place[i] = 0;
1207	place[root->index - 1] = 1;
1208	for (i = 1; i <= (spp); i++) {
1209	p = treenode[i - 1];
1210	while (place[p->index - 1] == 0) {
1211	place[p->index - 1] = i;
1212	while (!p->bottom)
1213	p = p->next;
1214	p = p->back;
1215	}
1216	if (i > 1) {
1217	place[i - 1] = place[p->index - 1];
1218	j = place[p->index - 1];
1219	done = false;
1220	while (!done) {
1221	place[p->index - 1] = spp + i - 1;
1222	while (!p->bottom)
1223	p = p->next;
1224	p = p->back;
1225	done = (p == NULL);
1226	if (!done)
1227	done = (place[p->index - 1] != j);
1228	}
1229	}
1230	}
1231	} /* savetree */
1232
1233	void
1234	findtree(pos, found)
1235	long *pos;
1236	boolean *found;
1237	{
1238	/*
1239	* finds tree given by ARRAY place in ARRAY bestrees by binary search
1240	*/
1241	long i, lower, upper;
1242	boolean below, done;
1243
1244	below = false;
1245	lower = 1;
1246	upper = nextree - 1;
1247	*found = false;
1248	while (!(*found) && lower <= upper) {
1249	(*pos) = (lower + upper) / 2;
1250	i = 3;
1251	done = false;
1252	while (!done) {
1253	done = (i > spp);
1254	if (!done)
1255	done = (place[i - 1] != bestrees[(*pos) - 1][i - 1]);
1256	if (!done)
1257	i++;
1258	}
1259	(*found) = (i > spp);
1260	below = (place[i - 1] < bestrees[(*pos) - 1][i - 1]);
1261	if (*found)
1262	break;
1263	if (below)
1264	upper = (*pos) - 1;
1265	else
1266	lower = (*pos) + 1;
1267	}
1268	if (!(*found) && !below)
1269	(*pos)++;
1270	} /* findtree */
1271
1272	void
1273	addtree(pos, found)
1274	long *pos;
1275	boolean *found;
1276	{
1277	/*
1278	* puts tree from ARRAY place in its proper position in ARRAY
1279	* bestrees
1280	*/
1281	long i, FORLIM;
1282
1283	for (i = nextree - 1; i >= (*pos); i--)
1284	memcpy(bestrees[i], bestrees[i - 1], spp * sizeof(long));
1285	for (i = 0; i < (spp); i++)
1286	bestrees[(*pos) - 1][i] = place[i];
1287	nextree++;
1288	} /* addtree */
1289
1290	void
1291	tryadd(p, item, nufork)
1292	node p, item, *nufork;
1293	{
1294	/*
1295	* temporarily adds one fork and one tip to the tree. if the location
1296	* where they are added yields greater "likelihood" than other
1297	* locations tested up to that time, then keeps that location as
1298	* there
1299	*/
1300	/* Local variables for tryadd: */
1301	long pos;
1302	boolean found;
1303	node rute, q;
1304
1305	if (p == root)
1306	fillin(temp, *item, p);
1307	else {
1308	fillin(temp1, *item, p);
1309	fillin(temp, temp1, p->back);
1310	}
1311	evaluate(temp);
1312	if (lastrearr) {
1313	if (like < bestlike) {
1314	if ((item) == (nufork)->next->next->back) {
1315	q = (*nufork)->next;
1316	(nufork)->next = (nufork)->next->next;
1317	(*nufork)->next->next = q;
1318	q->next = (*nufork);
1319	}
1320	} else if (like >= bstlike2) {
1321	recompute = false;
1322	add(p, (item), (nufork));
1323	rute = root->next->back;
1324	savetree(&pos, &found);
1325	reroot(rute);
1326	if (like > bstlike2) {
1327	bestlike = bstlike2 = like;
1328	pos = 1;
1329	nextree = 1;
1330	addtree(&pos, &found);
1331	} else {
1332	pos = 0;
1333	findtree(&pos, &found);
1334	if (!found) {
1335	if (nextree <= maxtrees)
1336	addtree(&pos, &found);
1337	}
1338	}
1339	re_move(item, nufork);
1340	recompute = true;
1341	}
1342	}
1343	if (like > bestyet) {
1344	bestyet = like;
1345	there = p;
1346	}
1347	} /* tryadd */
1348
1349	void
1350	addpreorder(p, item, nufork)
1351	node p, item, *nufork;
1352	{
1353	/*
1354	* traverses a binary tree, calling PROCEDURE tryadd at a node before
1355	* calling tryadd at its descendants
1356	*/
1357
1358	if (p == NULL)
1359	return;
1360	tryadd(p, &item, &nufork);
1361	if (!p->tip) {
1362	addpreorder(p->next->back, item, nufork);
1363	addpreorder(p->next->next->back, item, nufork);
1364	}
1365	} /* addpreorder */
1366
1367	void
1368	tryrearr(p, success)
1369	node *p;
1370	boolean *success;
1371	{
1372	/*
1373	* evaluates one rearrangement of the tree. if the new tree has
1374	* greater "likelihood" than the old one sets success := TRUE and
1375	* keeps the new tree. otherwise, restores the old tree
1376	*/
1377	node frombelow, whereto, forknode, q;
1378	double oldlike;
1379
1380	if (p->back == NULL)
1381	return;
1382	forknode = treenode[p->back->index - 1];
1383	if (forknode->back == NULL)
1384	return;
1385	oldlike = bestyet;
1386	if (p->back->next->next == forknode)
1387	frombelow = forknode->next->next->back;
1388	else
1389	frombelow = forknode->next->back;
1390	whereto = treenode[forknode->back->index - 1];
1391	if (whereto->next->back == forknode)
1392	q = whereto->next->next->back;
1393	else
1394	q = whereto->next->back;
1395	fillin(temp1, frombelow, q);
1396	fillin(temp, temp1, p);
1397	fillin(temp1, temp, whereto->back);
1398	evaluate(temp1);
1399	if (like <= oldlike) {
1400	if (p == forknode->next->next->back) {
1401	q = forknode->next;
1402	forknode->next = forknode->next->next;
1403	forknode->next->next = q;
1404	q->next = forknode;
1405	}
1406	} else {
1407	recompute = false;
1408	re_move(&p, &forknode);
1409	fillin(whereto, whereto->next->back, whereto->next->next->back);
1410	recompute = true;
1411	add(whereto, p, forknode);
1412	*success = true;
1413	bestyet = like;
1414	}
1415	} /* tryrearr */
1416
1417	void
1418	repreorder(p, success)
1419	node *p;
1420	boolean *success;
1421	{
1422	/*
1423	* traverses a binary tree, calling PROCEDURE tryrearr at a node
1424	* before calling tryrearr at its descendants
1425	*/
1426	if (p == NULL)
1427	return;
1428	tryrearr(p, success);
1429	if (!p->tip) {
1430	repreorder(p->next->back, success);
1431	repreorder(p->next->next->back, success);
1432	}
1433	} /* repreorder */
1434
1435	void
1436	rearrange(r)
1437	node **r;
1438	{
1439	/*
1440	* traverses the tree (preorder), finding any local rearrangement
1441	* which decreases the number of steps. if traversal succeeds in
1442	* increasing the tree's "likelihood", PROCEDURE rearrange runs
1443	* traversal again
1444	*/
1445	boolean success = true;
1446	while (success) {
1447	success = false;
1448	repreorder(*r, &success);
1449	}
1450	} /* rearrange */
1451
1452
1453	void
1454	findch(c)
1455	Char c;
1456	{
1457	/* scan forward until find character c */
1458	boolean done;
1459
1460	done = false;
1461	while (!done) {
1462	if (c == ',') {
1463	if (ch == '(' \|\| ch == ')' \|\| ch == ';') {
1464	printf("\nERROR IN USER TREE:");
1465	printf(" UNMATCHED PARENTHESIS OR MISSING COMMA\n");
1466	exit(-1);
1467	} else if (ch == ',')
1468	done = true;
1469	} else if (c == ')') {
1470	if (ch == '(' \|\| ch == ',' \|\| ch == ';') {
1471	printf("\nERROR IN USER TREE:");
1472	printf(" UNMATCHED PARENTHESIS OR NOT BIFURCATED NODE\n");
1473	exit(-1);
1474	} else {
1475	if (ch == ')')
1476	done = true;
1477	}
1478	} else if (c == ';') {
1479	if (ch != ';') {
1480	printf("\nERROR IN USER TREE:");
1481	printf(" UNMATCHED PARENTHESIS OR MISSING SEMICOLON\n");
1482	exit(-1);
1483	} else
1484	done = true;
1485	}
1486	if ((done && ch == ')') \|\| !done) {
1487	if (eoln(infile)) {
1488	fscanf(infile, "%*[^\n]");
1489	getc(infile);
1490	}
1491	ch = getc(infile);
1492	}
1493	}
1494	} /* findch */
1495
1496	void
1497	addelement(p, nextnode, lparens, names)
1498	node **p;
1499	long nextnode, lparens;
1500	boolean *names;
1501	{
1502	/* recursive procedure adds nodes to user-defined tree */
1503	node *q;
1504	long i, n;
1505	boolean found;
1506	Char str[nmlngth];
1507
1508	do {
1509	if (eoln(infile)) {
1510	fscanf(infile, "%*[^\n]");
1511	getc(infile);
1512	}
1513	ch = getc(infile);
1514	} while (ch == ' ');
1515	if (ch == '(') {
1516	if ((*lparens) >= spp - 1) {
1517	printf("\nERROR IN USER TREE: TOO MANY LEFT PARENTHESES\n");
1518	exit(-1);
1519	}
1520	(*nextnode)++;
1521	(*lparens)++;
1522	q = treenode[(*nextnode) - 1];
1523	addelement(&q->next->back, nextnode, lparens, names);
1524	q->next->back->back = q->next;
1525	findch(',');
1526	addelement(&q->next->next->back, nextnode, lparens, names);
1527	q->next->next->back->back = q->next->next;
1528	findch(')');
1529	*p = q;
1530	return;
1531	}
1532	for (i = 0; i < nmlngth; i++)
1533	str[i] = ' ';
1534	n = 1;
1535	do {
1536	if (ch == '_')
1537	ch = ' ';
1538	str[n - 1] = ch;
1539	if (eoln(infile)) {
1540	fscanf(infile, "%*[^\n]");
1541	getc(infile);
1542	}
1543	ch = getc(infile);
1544	n++;
1545	} while (ch != ',' && ch != ')' && ch != ':' && n <= nmlngth);
1546	n = 1;
1547	do {
1548	found = true;
1549	for (i = 0; i < nmlngth; i++)
1550	found = (found && str[i] == nayme[n - 1][i]);
1551	if (found) {
1552	if (names[n - 1] == false) {
1553	*p = treenode[n - 1];
1554	names[n - 1] = true;
1555	} else {
1556	printf("\nERROR IN USER TREE: DUPLICATE NAME FOUND -- ");
1557	for (i = 0; i < nmlngth; i++)
1558	putchar(nayme[n - 1][i]);
1559	putchar('\n');
1560	exit(-1);
1561	}
1562	} else
1563	n++;
1564	} while (!(n > spp \|\| found));
1565	if (n <= spp)
1566	return;
1567	printf("CANNOT FIND SPECIES: ");
1568	for (i = 0; i < nmlngth; i++)
1569	putchar(str[i]);
1570	putchar('\n');
1571	} /* addelement */
1572
1573	void
1574	treeread()
1575	{
1576	/* read in user-defined tree and set it up */
1577	long nextnode, lparens;
1578	long i;
1579
1580	root = treenode[spp];
1581	nextnode = spp;
1582	root->back = NULL;
1583	for (i = 0; i < (spp); i++)
1584	names[i] = false;
1585	lparens = 0;
1586	addelement(&root, &nextnode, &lparens, names);
1587	findch(';');
1588	if (progress)
1589	printf("\n\n");
1590	fscanf(infile, "%*[^\n]");
1591	getc(infile);
1592	} /* treeread */
1593
1594
1595	void
1596	coordinates(p, tipy)
1597	node *p;
1598	long *tipy;
1599	{
1600	/* establishes coordinates of nodes */
1601	node q, first, *last;
1602
1603	if (p->tip) {
1604	p->xcoord = 0;
1605	p->ycoord = (*tipy);
1606	p->ymin = (*tipy);
1607	p->ymax = (*tipy);
1608	(*tipy) += down;
1609	return;
1610	}
1611	q = p->next;
1612	do {
1613	coordinates(q->back, tipy);
1614	q = q->next;
1615	} while (p != q);
1616	first = p->next->back;
1617	q = p->next;
1618	while (q->next != p)
1619	q = q->next;
1620	last = q->back;
1621	p->xcoord = last->ymax - first->ymin;
1622	p->ycoord = (first->ycoord + last->ycoord) / 2;
1623	p->ymin = first->ymin;
1624	p->ymax = last->ymax;
1625	} /* coordinates */
1626
1627	void
1628	drawline(i, scale)
1629	long i;
1630	double scale;
1631	{
1632	/* draws one row of the tree diagram by moving up tree */
1633	node p, q, r, first, *last;
1634	long n, j;
1635	boolean extra, done;
1636
1637	p = root;
1638	q = root;
1639	extra = false;
1640	if (i == p->ycoord && p == root) {
1641	if (p->index - spp >= 10)
1642	fprintf(outfile, "-%2ld", p->index - spp);
1643	else
1644	fprintf(outfile, "--%ld", p->index - spp);
1645	extra = true;
1646	} else
1647	fprintf(outfile, " ");
1648	do {
1649	if (!p->tip) {
1650	r = p->next;
1651	done = false;
1652	do {
1653	if (i >= r->back->ymin && i <= r->back->ymax) {
1654	q = r->back;
1655	done = true;
1656	}
1657	r = r->next;
1658	} while (!(done \|\| r == p));
1659	first = p->next->back;
1660	r = p->next;
1661	while (r->next != p)
1662	r = r->next;
1663	last = r->back;
1664	}
1665	done = (p == q);
1666	n = (long) (scale * (p->xcoord - q->xcoord) + 0.5);
1667	if (n < 3 && !q->tip)
1668	n = 3;
1669	if (extra) {
1670	n--;
1671	extra = false;
1672	}
1673	if (q->ycoord == i && !done) {
1674	putc('+', outfile);
1675	if (!q->tip) {
1676	for (j = 1; j <= n - 2; j++)
1677	putc('-', outfile);
1678	if (q->index - spp >= 10)
1679	fprintf(outfile, "%2ld", q->index - spp);
1680	else
1681	fprintf(outfile, "-%ld", q->index - spp);
1682	extra = true;
1683	} else {
1684	for (j = 1; j < n; j++)
1685	putc('-', outfile);
1686	}
1687	} else if (!p->tip) {
1688	if (last->ycoord > i && first->ycoord < i && i != p->ycoord) {
1689	putc('!', outfile);
1690	for (j = 1; j < n; j++)
1691	putc(' ', outfile);
1692	} else {
1693	for (j = 1; j <= n; j++)
1694	putc(' ', outfile);
1695	}
1696	} else {
1697	for (j = 1; j <= n; j++)
1698	putc(' ', outfile);
1699	}
1700	if (p != q)
1701	p = q;
1702	} while (!done);
1703	if (p->ycoord == i && p->tip) {
1704	for (j = 0; j < nmlngth; j++)
1705	putc(nayme[p->index - 1][j], outfile);
1706	}
1707	putc('\n', outfile);
1708	} /* drawline */
1709
1710	void
1711	printree()
1712	{
1713	/* prints out diagram of the tree */
1714	/* Local variables for printree: */
1715	long tipy;
1716	double scale;
1717	long i;
1718
1719	putc('\n', outfile);
1720	if (!treeprint)
1721	return;
1722	putc('\n', outfile);
1723	tipy = 1;
1724	coordinates(root, &tipy);
1725	scale = 1.5;
1726	putc('\n', outfile);
1727	for (i = 1; i <= (tipy - down); i++)
1728	drawline(i, scale);
1729	fprintf(outfile, "\n remember:");
1730	if (outgropt)
1731	fprintf(outfile, " (although rooted by outgroup)");
1732	fprintf(outfile, " this is an unrooted tree!\n\n");
1733	} /* printree */
1734
1735
1736
1737	void
1738	ancestset(a, b, c, d, k)
1739	long a, b, c, d;
1740	long *k;
1741	{
1742	/* sets up the aa set array. */
1743	aas aa;
1744	long s, sa, sb;
1745	long i, j, m, n;
1746	boolean counted;
1747
1748	counted = false;
1749	*k = 0;
1750	for (i = 0; i <= 5; i++) {
1751	if (*k < 3) {
1752	s = 0;
1753	if (i > 3)
1754	n = i - 3;
1755	else
1756	n = 0;
1757	for (j = n; j <= (i - n); j++) {
1758	if (j < 3)
1759	sa = a[j];
1760	else
1761	sa = fullset;
1762	for (m = n; m <= (i - j - n); m++) {
1763	if (m < 3)
1764	sb = sa & b[m];
1765	else
1766	sb = sa;
1767	if (i - j - m < 3)
1768	sb &= c[i - j - m];
1769	s \|= sb;
1770	}
1771	}
1772	if (counted \|\| s != 0) {
1773	d[*k] = s;
1774	(*k)++;
1775	counted = true;
1776	}
1777	}
1778	}
1779	for (i = 0; i <= 1; i++) {
1780	for (aa = ala; (long) aa <= (long) stop; aa = (aas) ((long) aa + 1)) {
1781	if (((1L << ((long) aa)) & d[i]) != 0) {
1782	for (j = i + 1; j <= 2; j++)
1783	d[j] \|= translate[(long) aa - (long) ala][j - i];
1784	}
1785	}
1786	}
1787	} /* ancestset */
1788
1789
1790	void
1791	hyprint(b1, b2, bottom, r, nonzero, maybe, nothing)
1792	long b1, b2;
1793	boolean bottom, nonzero, *maybe;
1794	node *r;
1795	sitearray nothing;
1796	{
1797	/* print out states in sites b1 through b2 at node */
1798	long i;
1799	boolean dot;
1800	Char ch;
1801	aas aa;
1802
1803	if (*bottom) {
1804	if (!outgropt)
1805	fprintf(outfile, " ");
1806	else
1807	fprintf(outfile, "root ");
1808	} else
1809	fprintf(outfile, "%3ld ", r->back->index - spp);
1810	if (r->tip) {
1811	for (i = 0; i < nmlngth; i++)
1812	putc(nayme[r->index - 1][i], outfile);
1813	} else
1814	fprintf(outfile, "%4ld ", r->index - spp);
1815	if (*bottom)
1816	fprintf(outfile, " ");
1817	else if (*nonzero)
1818	fprintf(outfile, " yes ");
1819	else if (*maybe)
1820	fprintf(outfile, " maybe ");
1821	else
1822	fprintf(outfile, " no ");
1823	for (i = b1 - 1; i < b2; i++) {
1824	aa = r->seq[i];
1825	switch (aa) {
1826
1827	case ala:
1828	ch = 'A';
1829	break;
1830
1831	case asx:
1832	ch = 'B';
1833	break;
1834
1835	case cys:
1836	ch = 'C';
1837	break;
1838
1839	case asp:
1840	ch = 'D';
1841	break;
1842
1843	case glu:
1844	ch = 'E';
1845	break;
1846
1847	case phe:
1848	ch = 'F';
1849	break;
1850
1851	case gly:
1852	ch = 'G';
1853	break;
1854
1855	case his:
1856	ch = 'H';
1857	break;
1858
1859	case ileu:
1860	ch = 'I';
1861	break;
1862
1863	case lys:
1864	ch = 'K';
1865	break;
1866
1867	case leu:
1868	ch = 'L';
1869	break;
1870
1871	case met:
1872	ch = 'M';
1873	break;
1874
1875	case asn:
1876	ch = 'N';
1877	break;
1878
1879	case pro:
1880	ch = 'P';
1881	break;
1882
1883	case gln:
1884	ch = 'Q';
1885	break;
1886
1887	case arg:
1888	ch = 'R';
1889	break;
1890
1891	case ser:
1892	ch = 'S';
1893	break;
1894
1895	case ser1:
1896	ch = 'S';
1897	break;
1898
1899	case ser2:
1900	ch = 'S';
1901	break;
1902
1903	case thr:
1904	ch = 'T';
1905	break;
1906
1907	case trp:
1908	ch = 'W';
1909	break;
1910
1911	case tyr:
1912	ch = 'Y';
1913	break;
1914
1915	case val:
1916	ch = 'V';
1917	break;
1918
1919	case glx:
1920	ch = 'Z';
1921	break;
1922
1923	case del:
1924	ch = '-';
1925	break;
1926
1927	case stop:
1928	ch = '*';
1929	break;
1930
1931	case unk:
1932	ch = 'X';
1933	break;
1934
1935	case quest:
1936	ch = '?';
1937	break;
1938	}
1939	if (!(*bottom))
1940	dot = (r->siteset[i][0] == treenode[r->back->index - 1]->siteset[i][0]
1941	\|\| ((r->siteset[i][0] &
1942	(~((1L << ((long) ser1)) \| (1L << ((long) ser2)) \|
1943	(1L << ((long) ser))))) == 0 &&
1944	(treenode[r->back->index - 1]->siteset[i][0] &
1945	(~((1L << ((long) ser1)) \| (1L << ((long) ser2)) \|
1946	(1L << ((long) ser))))) == 0));
1947	else
1948	dot = false;
1949	if (dot)
1950	putc('.', outfile);
1951	else
1952	putc(ch, outfile);
1953	if ((i + 1) % 10 == 0)
1954	putc(' ', outfile);
1955	}
1956	putc('\n', outfile);
1957	} /* hyprint */
1958
1959	void
1960	hyptrav(r, hypset, b1, b2, k, bottom, nothing)
1961	node *r;
1962	sitearray *hypset;
1963	long b1, b2, *k;
1964	boolean *bottom;
1965	sitearray nothing;
1966	{
1967	boolean maybe, nonzero;
1968	long i;
1969	aas aa;
1970	long anc, hset;
1971	gseq ancset, temparray;
1972
1973	gnu(&ancset);
1974	gnu(&temparray);
1975	maybe = false;
1976	nonzero = false;
1977	for (i = b1 - 1; i < b2; i++) {
1978	if (!r->tip) {
1979	ancestset(hypset[i], r->next->back->siteset[i],
1980	r->next->next->back->siteset[i], temparray->seq[i], k);
1981	memcpy(r->siteset[i], temparray->seq[i], sizeof(sitearray));
1982	}
1983	if (!(*bottom))
1984	anc = treenode[r->back->index - 1]->siteset[i][0];
1985	if (!r->tip) {
1986	hset = r->siteset[i][0];
1987	r->seq[i] = quest;
1988	for (aa = ala; (long) aa <= (long) stop; aa = (aas) ((long) aa + 1)) {
1989	if (hset == 1L << ((long) aa))
1990	r->seq[i] = aa;
1991	}
1992	if (hset == ((1L << ((long) asn)) \| (1L << ((long) asp))))
1993	r->seq[i] = asx;
1994	if (hset == ((1L << ((long) gln)) \| (1L << ((long) gly))))
1995	r->seq[i] = glx;
1996	if (hset == ((1L << ((long) ser1)) \| (1L << ((long) ser2))))
1997	r->seq[i] = ser;
1998	if (hset == fullset)
1999	r->seq[i] = unk;
2000	}
2001	nonzero = (nonzero \|\| (r->siteset[i][0] & anc) == 0);
2002	maybe = (maybe \|\| r->siteset[i][0] != anc);
2003	}
2004	hyprint(b1, b2, bottom, r, &nonzero, &maybe, nothing);
2005	*bottom = false;
2006	if (!r->tip) {
2007	memcpy(temparray->seq, r->next->back->siteset, chars * sizeof(sitearray));
2008	for (i = b1 - 1; i < b2; i++)
2009	ancestset(hypset[i], r->next->next->back->siteset[i], nothing,
2010	ancset->seq[i], k);
2011	hyptrav(r->next->back, ancset->seq, b1, b2, k, bottom, nothing);
2012	for (i = b1 - 1; i < b2; i++)
2013	ancestset(hypset[i], temparray->seq[i], nothing, ancset->seq[i], k);
2014	hyptrav(r->next->next->back, ancset->seq, b1, b2, k, bottom, nothing);
2015	}
2016	chuck(temparray);
2017	chuck(ancset);
2018	} /* hyptrav */
2019
2020	void
2021	hypstates(k)
2022	long *k;
2023	{
2024	/* fill in and describe states at interior nodes */
2025	boolean bottom;
2026	sitearray nothing;
2027	long i, n;
2028	seqptr hypset;
2029
2030	fprintf(outfile, "\nFrom To Any Steps? State at upper node\n");
2031	fprintf(outfile, " ");
2032	fprintf(outfile, "( . means same as in the node below it on tree)\n\n");
2033	memcpy(nothing, translate[(long) quest - (long) ala], sizeof(sitearray));
2034	hypset = (seqptr) Malloc(chars * sizeof(sitearray));
2035	for (i = 0; i < (chars); i++)
2036	memcpy(hypset[i], nothing, sizeof(sitearray));
2037	bottom = true;
2038	for (i = 1; i <= ((chars - 1) / 40 + 1); i++) {
2039	putc('\n', outfile);
2040	n = i * 40;
2041	if (n > chars)
2042	n = chars;
2043	bottom = true;
2044	hyptrav(root, hypset, i * 40 - 39, n, k, &bottom, nothing);
2045	}
2046	free(hypset);
2047	} /* hypstates */
2048
2049	void
2050	treeout(p)
2051	node *p;
2052	{
2053	/* write out file with representation of final tree */
2054	long i, n;
2055	Char c;
2056
2057	if (p->tip) {
2058	n = 0;
2059	for (i = 1; i <= nmlngth; i++) {
2060	if (nayme[p->index - 1][i - 1] != ' ')
2061	n = i;
2062	}
2063	for (i = 0; i < n; i++) {
2064	c = nayme[p->index - 1][i];
2065	if (c == ' ')
2066	c = '_';
2067	putc(c, treefile);
2068	}
2069	col += n;
2070	} else {
2071	putc('(', treefile);
2072	col++;
2073	treeout(p->next->back);
2074	putc(',', treefile);
2075	col++;
2076	if (col > 65) {
2077	putc('\n', treefile);
2078	col = 0;
2079	}
2080	treeout(p->next->next->back);
2081	putc(')', treefile);
2082	col++;
2083	}
2084	if (p != root)
2085	return;
2086	if (nextree > 2)
2087	fprintf(treefile, "[%6.4f];\n", 1.0 / (nextree - 1));
2088	else
2089	fprintf(treefile, ";\n");
2090	} /* treeout */
2091
2092	void
2093	describe()
2094	{
2095	/*
2096	* prints ancestors, steps and table of numbers of steps in each site
2097	*/
2098	long i, j, k;
2099
2100	if (treeprint)
2101	fprintf(outfile, "\nrequires a total of %10.3f\n", like / -10);
2102	if (stepbox) {
2103	putc('\n', outfile);
2104	if (weights)
2105	fprintf(outfile, " weighted");
2106	fprintf(outfile, "steps in each position:\n");
2107	fprintf(outfile, " ");
2108	for (i = 0; i <= 9; i++)
2109	fprintf(outfile, "%4ld", i);
2110	fprintf(outfile, "\n *-----------------------------------------\n");
2111	for (i = 0; i <= (chars / 10); i++) {
2112	fprintf(outfile, "%5ld", i * 10);
2113	putc('!', outfile);
2114	for (j = 0; j <= 9; j++) {
2115	k = i * 10 + j;
2116	if (k == 0 \|\| k > chars)
2117	fprintf(outfile, " ");
2118	else
2119	fprintf(outfile, "%4ld", root->numsteps[k - 1] / 10);
2120	}
2121	putc('\n', outfile);
2122	}
2123	}
2124	if (ancseq) {
2125	hypstates(&k);
2126	putc('\n', outfile);
2127	}
2128	putc('\n', outfile);
2129	if (trout) {
2130	col = 0;
2131	treeout(root);
2132	}
2133	} /* describe */
2134
2135
2136	void
2137	maketree()
2138	{
2139	/*
2140	* constructs a binary tree from the pointers in treenode. adds each
2141	* node at location which yields highest "likelihood" then rearranges
2142	* the tree for greatest "likelihood"
2143	*/
2144	long i, j, k, numtrees, num;
2145	double gotlike, wt, sumw, sum, sum2, sd;
2146	node item, nufork, *dummy;
2147	double TEMP;
2148
2149	if (!usertree) {
2150	for (i = 1; i <= (spp); i++)
2151	enterorder[i - 1] = i;
2152	if (jumble) {
2153	for (i = 0; i < (spp); i++) {
2154	j = (long) (randum(seed) * spp) + 1;
2155	k = enterorder[j - 1];
2156	enterorder[j - 1] = enterorder[i];
2157	enterorder[i] = k;
2158	}
2159	}
2160	root = treenode[enterorder[0] - 1];
2161	recompute = true;
2162	add(treenode[enterorder[0] - 1], treenode[enterorder[1] - 1],
2163	treenode[spp]);
2164	if (progress) {
2165	printf("Adding species:\n");
2166	printf(" ");
2167	for (i = 0; i < nmlngth; i++)
2168	putchar(nayme[enterorder[0] - 1][i]);
2169	printf("\n ");
2170	for (i = 0; i < nmlngth; i++)
2171	putchar(nayme[enterorder[1] - 1][i]);
2172	putchar('\n');
2173	}
2174	lastrearr = false;
2175	for (i = 3; i <= (spp); i++) {
2176	bestyet = -1.0e6;
2177	there = root;
2178	item = treenode[enterorder[i - 1] - 1];
2179	nufork = treenode[spp + i - 2];
2180	addpreorder(root, item, nufork);
2181	add(there, item, nufork);
2182	like = bestyet;
2183	rearrange(&root);
2184	if (progress) {
2185	printf(" ");
2186	for (j = 0; j < nmlngth; j++)
2187	putchar(nayme[enterorder[i - 1] - 1][j]);
2188	putchar('\n');
2189	}
2190	lastrearr = (i == spp);
2191	if (lastrearr) {
2192	if (progress) {
2193	printf("\nDoing global rearrangements\n");
2194	printf(" !");
2195	for (j = 1; j <= nonodes; j++)
2196	putchar('-');
2197	printf("!\n");
2198	}
2199	bestlike = bestyet;
2200	if (jumb == 1) {
2201	bstlike2 = bestlike;
2202	nextree = 1;
2203	}
2204	do {
2205	if (progress)
2206	printf(" ");
2207	gotlike = bestlike;
2208	for (j = 0; j < (nonodes); j++) {
2209	bestyet = -1.0e6;
2210	item = treenode[j];
2211	if (item != root) {
2212	nufork = treenode[treenode[j]->back->index - 1];
2213	re_move(&item, &nufork);
2214	there = root;
2215	addpreorder(root, item, nufork);
2216	add(there, item, nufork);
2217	}
2218	if (progress) {
2219	putchar('.');
2220	fflush(stdout);
2221	}
2222	}
2223	if (progress)
2224	putchar('\n');
2225	} while (bestlike > gotlike);
2226	}
2227	}
2228	if (progress)
2229	putchar('\n');
2230	for (i = spp - 1; i >= 1; i--)
2231	re_move(&treenode[i], &dummy);
2232	if (jumb == njumble) {
2233	if (treeprint) {
2234	putc('\n', outfile);
2235	if (nextree == 2)
2236	fprintf(outfile, "One most parsimonious tree found:\n");
2237	else
2238	fprintf(outfile, "%6ld trees in all found\n", nextree - 1);
2239	}
2240	if (nextree > maxtrees + 1) {
2241	if (treeprint)
2242	fprintf(outfile, "here are the first%4ld of them\n", (long) maxtrees);
2243	nextree = maxtrees + 1;
2244	}
2245	if (treeprint)
2246	putc('\n', outfile);
2247	recompute = false;
2248	for (i = 0; i <= (nextree - 2); i++) {
2249	root = treenode[0];
2250	add(treenode[0], treenode[1], treenode[spp]);
2251	for (j = 3; j <= (spp); j++)
2252	add(treenode[bestrees[i][j - 1] - 1], treenode[j - 1],
2253	treenode[spp + j - 2]);
2254	reroot(treenode[outgrno - 1]);
2255	postorder(root);
2256	evaluate(root);
2257	printree();
2258	describe();
2259	for (j = 1; j < (spp); j++)
2260	re_move(&treenode[j], &dummy);
2261	}
2262	}
2263	} else {
2264	fscanf(infile, "%ld%*[^\n]", &numtrees);
2265	getc(infile);
2266	if (treeprint) {
2267	fprintf(outfile, "User-defined tree");
2268	if (numtrees > 1)
2269	putc('s', outfile);
2270	fprintf(outfile, ":\n\n\n\n");
2271	}
2272	which = 1;
2273	while (which <= numtrees) {
2274	treeread();
2275	if (outgropt)
2276	reroot(treenode[outgrno - 1]);
2277	postorder(root);
2278	evaluate(root);
2279	printree();
2280	describe();
2281	which++;
2282	}
2283	putc('\n', outfile);
2284	if (numtrees > 1 && chars > 1) {
2285	fprintf(outfile, "Tree Steps Diff Steps Its S.D.");
2286	fprintf(outfile, " Significantly worse?\n\n");
2287	if (numtrees > maxuser)
2288	num = maxuser;
2289	else
2290	num = numtrees;
2291	which = 1;
2292	while (which <= num) {
2293	fprintf(outfile, "%3ld%10.1f", which, nsteps[which - 1] / 10);
2294	if (minwhich == which)
2295	fprintf(outfile, " <------ best\n");
2296	else {
2297	sumw = 0.0;
2298	sum = 0.0;
2299	sum2 = 0.0;
2300	for (j = 0; j < (chars); j++) {
2301	if (weight[j] > 0) {
2302	wt = weight[j] / 10.0;
2303	sumw += wt;
2304	sum += (fsteps[which - 1][j] -
2305	fsteps[minwhich - 1][j]) / 10.0;
2306	TEMP = (fsteps[which - 1][j] -
2307	fsteps[minwhich - 1][j]) / 10.0;
2308	sum2 += TEMP * TEMP / wt;
2309	}
2310	}
2311	TEMP = sum / sumw;
2312	sd = sqrt(sumw / (sumw - 1.0) * (sum2 - TEMP * TEMP));
2313	fprintf(outfile, "%10.1f%12.4f",
2314	(nsteps[which - 1] - minsteps) / 10, sd);
2315	if (sum > 1.95996 * sd)
2316	fprintf(outfile, " Yes\n");
2317	else
2318	fprintf(outfile, " No\n");
2319	}
2320	which++;
2321	}
2322	fprintf(outfile, "\n\n");
2323	}
2324	}
2325	if (jumb == njumble && progress) {
2326	printf("Output written to output file\n\n");
2327	if (trout)
2328	printf("Trees also written onto file\n\n");
2329	}
2330	} /* maketree */
2331
2332
2333	main(argc, argv)
2334	int argc;
2335	Char *argv[];
2336	{ /* Protein parsimony by uphill search */
2337	char infilename[100], outfilename[100], trfilename[100];
2338	#ifdef MAC
2339	macsetup("Protpars", "");
2340	#endif
2341	openfile(&infile, INFILE, "r", argv[0], infilename);
2342	openfile(&outfile, OUTFILE, "w", argv[0], outfilename);
2343
2344	ibmpc = ibmpc0;
2345	ansi = ansi0;
2346	vt52 = vt520;
2347	garbage = NULL;
2348	mulsets = false;
2349	datasets = 1;
2350	firstset = true;
2351	setup();
2352	doinit();
2353	if (usertree) {
2354	fsteps = (long *) Malloc(maxuser sizeof(long *));
2355	for (j = 1; j <= maxuser; j++)
2356	fsteps[j - 1] = (long ) Malloc(chars sizeof(long));
2357	}
2358	bestrees = (long *) Malloc(maxtrees sizeof(long *));
2359	for (j = 1; j <= maxtrees; j++)
2360	bestrees[j - 1] = (long ) Malloc(spp sizeof(long));
2361	nayme = (Char *) Malloc(spp sizeof(Char *));
2362	for (j = 1; j <= spp; j++)
2363	nayme[j - 1] = (Char ) Malloc(nmlngth sizeof(Char));
2364	enterorder = (long ) Malloc(spp sizeof(long));
2365	names = (boolean ) Malloc(spp sizeof(boolean));
2366	place = (long ) Malloc(nonodes sizeof(long));
2367	weight = (steptr) Malloc(chars * sizeof(long));
2368	threshwt = (steptr) Malloc(chars * sizeof(long));
2369	temp = (node *) Malloc(sizeof(node));
2370	temp->numsteps = (steptr) Malloc(chars * sizeof(long));
2371	temp->siteset = (seqptr) Malloc(chars * sizeof(sitearray));
2372	temp->seq = (aas ) Malloc(chars sizeof(aas));
2373	temp1 = (node *) Malloc(sizeof(node));
2374	temp1->numsteps = (steptr) Malloc(chars * sizeof(long));
2375	temp1->siteset = (seqptr) Malloc(chars * sizeof(sitearray));
2376	temp1->seq = (aas ) Malloc(chars sizeof(aas));
2377	if (trout)
2378	openfile(&treefile, TREEFILE, "w", argv[0], trfilename);
2379	for (ith = 1; ith <= datasets; ith++) {
2380	doinput();
2381	if (ith == 1)
2382	firstset = false;
2383	if (datasets > 1) {
2384	fprintf(outfile, "Data set # %ld:\n\n", ith);
2385	if (progress)
2386	printf("Data set # %ld:\n\n", ith);
2387	}
2388	for (jumb = 1; jumb <= njumble; jumb++)
2389	maketree();
2390	}
2391	FClose(infile);
2392	FClose(outfile);
2393	FClose(treefile);
2394	#ifdef MAC
2395	fixmacfile(outfilename);
2396	fixmacfile(trfilename);
2397	#endif
2398	exit(0);
2399	} /* Protein parsimony by uphill search */
2400
2401
2402	int
2403	eof(f)
2404	FILE *f;
2405	{
2406	register int ch;
2407
2408	if (feof(f))
2409	return 1;
2410	if (f == stdin)
2411	return 0;
2412	ch = getc(f);
2413	if (ch == EOF)
2414	return 1;
2415	ungetc(ch, f);
2416	return 0;
2417	}
2418
2419
2420	int
2421	eoln(f)
2422	FILE *f;
2423	{
2424	register int ch;
2425
2426	ch = getc(f);
2427	if (ch == EOF)
2428	return 1;
2429	ungetc(ch, f);
2430	return (ch == '\n');
2431	}
2432
2433	void
2434	memerror()
2435	{
2436	printf("Error allocating memory\n");
2437	exit(-1);
2438	}
2439
2440	MALLOCRETURN *
2441	mymalloc(x)
2442	long x;
2443	{
2444	MALLOCRETURN *mem;
2445	mem = (MALLOCRETURN *) malloc(x);
2446	if (!mem)
2447	memerror();
2448	else
2449	return (MALLOCRETURN *) mem;
2450	}

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

Download in other formats: