Context Navigation

source: branches/port5/GDE/RAxML/axml.c

Visit:

Last change on this file was 5262, checked in by westram, 17 years ago
update to RAxML 7.0.3
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 193.1 KB

Line
1	/* RAxML-VI-HPC (version 2.2) a program for sequential and parallel estimation of phylogenetic trees
2	* Copyright August 2006 by Alexandros Stamatakis
3	*
4	* Partially derived from
5	* fastDNAml, a program for estimation of phylogenetic trees from sequences by Gary J. Olsen
6	*
7	* and
8	*
9	* Programs of the PHYLIP package by Joe Felsenstein.
10	*
11	* This program is free software; you may redistribute it and/or modify its
12	* under the terms of the GNU General Public License as published by the Free
13	* Software Foundation; either version 2 of the License, or (at your option)
14	* any later version.
15	*
16	* This program is distributed in the hope that it will be useful, but
17	* WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18	* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19	* for more details.
20	*
21	*
22	* For any other enquiries send an Email to Alexandros Stamatakis
23	* Alexandros.Stamatakis@epfl.ch
24	*
25	* When publishing work that is based on the results from RAxML-VI-HPC please cite:
26	*
27	* Alexandros Stamatakis:"RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models".
28	* Bioinformatics 2006; doi: 10.1093/bioinformatics/btl446
29	*/
30
31	#ifdef WIN32
32	#include <direct.h>
33	#endif
34
35	#ifndef WIN32
36	#include <sys/times.h>
37	#include <sys/types.h>
38	#include <sys/time.h>
39	#include <unistd.h>
40	#endif
41
42	#include <math.h>
43	#include <time.h>
44	#include <stdlib.h>
45	#include <stdio.h>
46	#include <ctype.h>
47	#include <string.h>
48
49	#ifdef PARALLEL
50	#include <mpi.h>
51	#endif
52
53
54
55	#ifdef _USE_PTHREADS
56	#include <pthread.h>
57	#endif
58
59
60	#include "axml.h"
61	#include "globalVariables.h"
62
63
64	/*************** UTILITY FUNCTIONS ************************/
65
66
67
68
69
70	double gettime(void)
71	{
72	#ifdef WIN32
73	time_t tp;
74	struct tm localtm;
75	tp = time(NULL);
76	localtm = *localtime(&tp);
77	return 60.0*localtm.tm_min + localtm.tm_sec;
78	#else
79	struct timeval ttime;
80	gettimeofday(&ttime , NULL);
81	return ttime.tv_sec + ttime.tv_usec * 0.000001;
82	#endif
83	}
84
85	int gettimeSrand(void)
86	{
87	#ifdef WIN32
88	time_t tp;
89	struct tm localtm;
90	tp = time(NULL);
91	localtm = *localtime(&tp);
92	return 246060localtm.tm_yday + 6060localtm.tm_hour + 60localtm.tm_min + localtm.tm_sec;
93	#else
94	struct timeval ttime;
95	gettimeofday(&ttime , NULL);
96	return ttime.tv_sec + ttime.tv_usec;
97	#endif
98	}
99
100	double randum (long *seed)
101	/* random number generator, modified to use 12 bit chunks */
102	{ /* randum */
103	long sum, mult0, mult1, seed0, seed1, seed2, newseed0, newseed1, newseed2;
104	double res;
105
106	mult0 = 1549;
107	seed0 = *seed & 4095;
108	sum = mult0 * seed0;
109	newseed0 = sum & 4095;
110	sum >>= 12;
111	seed1 = (*seed >> 12) & 4095;
112	mult1 = 406;
113	sum += mult0 * seed1 + mult1 * seed0;
114	newseed1 = sum & 4095;
115	sum >>= 12;
116	seed2 = (*seed >> 24) & 255;
117	sum += mult0 * seed2 + mult1 * seed1;
118	newseed2 = sum & 255;
119
120	*seed = newseed2 << 24 \| newseed1 << 12 \| newseed0;
121	res = 0.00390625 * (newseed2 + 0.000244140625 * (newseed1 + 0.000244140625 * newseed0));
122
123	return res;
124	} /* randum */
125
126	static int filexists(char *filename)
127	{
128	FILE *fp;
129	int res;
130	fp = fopen(filename,"r");
131
132	if(fp)
133	{
134	res = 1;
135	fclose(fp);
136	}
137	else
138	res = 0;
139
140	return res;
141	}
142
143	/******************* END UTILITY FUNCTIONS ******************/
144
145
146	/****************************some functions for the likelihood computation **************************/
147
148	#ifdef WIN32
149	boolean isTip(int number, int maxTips)
150	#else
151	inline boolean isTip(int number, int maxTips)
152	#endif
153	{
154	assert(number > 0);
155
156	if(number <= maxTips)
157	return TRUE;
158	else
159	return FALSE;
160	}
161
162
163	#ifdef WIN32
164	double getLikelihoodArray(int number, int mxtips, double *xVector)
165	#else
166	inline double getLikelihoodArray(int number, int mxtips, double *xVector)
167	#endif
168	{
169	return (xVector[number - mxtips - 1]);
170	}
171
172	#ifdef WIN32
173	int getScalingArray(int number, int endsite, int mxtips, int scalingArray)
174	#else
175	inline int getScalingArray(int number, int endsite, int mxtips, int scalingArray)
176	#endif
177	{
178	return &(scalingArray[endsite * (number - mxtips - 1)]);
179	}
180
181
182	#ifdef _MULTI_GENE
183	void getxsnode (nodeptr p, int model)
184	{
185	assert(p->xs[model] \|\| p->next->xs[model] \|\| p->next->next->xs[model]);
186	assert(p->xs[model] + p->next->xs[model] + p->next->next->xs[model] == 1);
187
188	assert(p == p->next->next->next);
189
190	p->xs[model] = 1;
191
192	if(p->next->xs[model])
193	{
194	p->next->xs[model] = 0;
195	return;
196	}
197	else
198	{
199	p->next->next->xs[model] = 0;
200	return;
201	}
202
203	assert(0);
204	}
205
206	#endif
207
208	void getxnode (nodeptr p)
209	{
210	nodeptr s;
211
212	if ((s = p->next)->x \|\| (s = s->next)->x)
213	{
214	p->x = s->x;
215	s->x = 0;
216	}
217
218	assert(p->x);
219	}
220
221
222
223
224
225	void hookup (nodeptr p, nodeptr q, double *z, int numBranches)
226	{
227	int i;
228
229	p->back = q;
230	q->back = p;
231
232	for(i = 0; i < numBranches; i++)
233	p->z[i] = q->z[i] = z[i];
234	}
235
236	void hookupDefault (nodeptr p, nodeptr q, int numBranches)
237	{
238	int i;
239
240	p->back = q;
241	q->back = p;
242
243	for(i = 0; i < numBranches; i++)
244	p->z[i] = q->z[i] = defaultz;
245	}
246
247
248	/****************************some functions for the likelihood computation **************************/
249
250
251
252
253	/*********************reading and initializing input ****************/
254
255	static void getnums (rawdata *rdta)
256	{
257	if (fscanf(INFILE, "%d %d", & rdta->numsp, & rdta->sites) != 2)
258	{
259	if(processID == 0)
260	printf("ERROR: Problem reading number of species and sites\n");
261	errorExit(-1);
262	}
263
264	if (rdta->numsp < 4)
265	{
266	if(processID == 0)
267	printf("TOO FEW SPECIES\n");
268	errorExit(-1);
269	}
270
271	if (rdta->sites < 1)
272	{
273	if(processID == 0)
274	printf("TOO FEW SITES\n");
275	errorExit(-1);
276	}
277
278	return;
279	}
280
281
282	static boolean digitchar (int ch)
283	{
284	return (ch >= '0' && ch <= '9');
285	}
286
287
288	boolean whitechar (int ch)
289	{
290	return (ch == ' ' \|\| ch == '\n' \|\| ch == '\t' \|\| ch == '\r'); /* PC-LINEBREAK*/
291	}
292
293
294	static void uppercase (int *chptr)
295	{
296	int ch;
297
298	ch = *chptr;
299	if ((ch >= 'a' && ch <= 'i') \|\| (ch >= 'j' && ch <= 'r')
300	\|\| (ch >= 's' && ch <= 'z'))
301	*chptr = ch + 'A' - 'a';
302	}
303
304
305
306
307	static void getyspace (rawdata *rdta)
308	{
309	long size;
310	int i;
311	char *y0;
312
313	if (! (rdta->y = (char *) malloc((rdta->numsp + 1) sizeof(char *))))
314	{
315	printf("ERROR: Unable to obtain space for data array pointers\n");
316	exit(-1);
317	}
318
319	size = 4 * (rdta->sites / 4 + 1);
320
321	if (! (y0 = (char ) malloc((rdta->numsp + 1) size * sizeof(char))))
322	{
323	printf("ERROR: Unable to obtain space for data array\n");
324	exit(-1);
325	}
326
327	rdta->y0 = y0;
328
329	for (i = 0; i <= rdta->numsp; i++)
330	{
331	rdta->y[i] = y0;
332	y0 += size;
333	}
334
335	return;
336	}
337
338	static boolean setupTree (tree tr, analdef adef)
339	{
340	nodeptr p0, p, q;
341	int i, j, tips, inter;
342
343	tr->expArray = (int*)NULL;
344	tr->likelihoodArray = (double*)NULL;
345	tr->sumBuffer = (double*)NULL;
346
347	tr->bigCutoff = FALSE;
348	tr->mixedData = FALSE;
349
350	tr->ib = (insertionBranch *)NULL;
351	tr->ip = (insertionPoints *)NULL;
352	tr->numberOfTipsForInsertion = 0;
353
354	for(i = 0; i < NUM_BRANCHES; i++)
355	tr->partitionContributions[i] = -1.0;
356
357
358	if(!adef->useMultipleModel)
359	tr->NumberOfModels = 1;
360
361	if(adef->grouping)
362	tr->grouped = TRUE;
363	else
364	tr->grouped = FALSE;
365
366	if(adef->constraint)
367	tr->constrained = TRUE;
368	else
369	tr->constrained = FALSE;
370
371	tr->treeID = 0;
372
373	tips = tr->mxtips;
374	inter = tr->mxtips - 1;
375
376	tr->tipVectorDNA = (double )malloc(tr->NumberOfModels 64 * sizeof(double));
377	tr->tipVectorAA = (double )malloc(tr->NumberOfModels 460 * sizeof(double));
378
379	tr->EV_DNA = (double )malloc(tr->NumberOfModels 16 * sizeof(double));
380	tr->EV_AA = (double )malloc(tr->NumberOfModels 400 * sizeof(double));
381
382	tr->EI_DNA = (double )malloc(tr->NumberOfModels 12 * sizeof(double));
383	tr->EI_AA = (double )malloc(tr->NumberOfModels 380 * sizeof(double));
384
385	tr->EIGN_DNA = (double )malloc(tr->NumberOfModels 3 * sizeof(double));
386	tr->EIGN_AA = (double )malloc(tr->NumberOfModels 19 * sizeof(double));
387
388	tr->frequencies_DNA = (double )malloc(tr->NumberOfModels 4 * sizeof(double));
389	tr->frequencies_AA = (double )malloc(tr->NumberOfModels 20 * sizeof(double));
390
391
392	tr->initialRates_DNA = (double )malloc(tr->NumberOfModels 5 * sizeof(double));
393	tr->initialRates_AA = (double )malloc(tr->NumberOfModels 190 * sizeof(double));
394
395	tr->xVector = (double *)malloc(tr->mxtips sizeof(double *));
396	tr->yVector = (char *) malloc((tr->mxtips + 1) sizeof(char *));
397
398	tr->gammaRates = (double )malloc(tr->NumberOfModels 4 * sizeof(double));
399	tr->alphas = (double )malloc(tr->NumberOfModels sizeof(double));
400	tr->invariants = (double )malloc(tr->NumberOfModels sizeof(double));
401	tr->fracchanges = (double )malloc(tr->NumberOfModels sizeof(double));
402	tr->likelihoods = (double )malloc(adef->multipleRuns sizeof(double));
403	tr->treeStringLength = tr->mxtips * (nmlngth+128) + 256 + tr->mxtips * 2;
404	tr->tree_string = (char )malloc(tr->treeStringLength sizeof(char));
405	/TODO, must that be so long ?/
406
407
408
409
410	/* tr->ti = (traversalInfo )malloc(sizeof(traversalInfo) tr->mxtips);*/
411
412	tr->td[0].count = 0;
413	tr->td[0].ti = (traversalInfo )malloc(sizeof(traversalInfo) tr->mxtips);
414
415	for(i = 0; i < tr->NumberOfModels; i++)
416	tr->fracchanges[i] = -1.0;
417	tr->fracchange = -1.0;
418
419
420	#ifdef _MULTI_GENE
421	tr->doMulti = 0;
422	{
423	int k;
424
425	for(k = 0; k < tr->numBranches; k++)
426	{
427	tr->td[k].count = 0;
428	tr->td[k].ti = (traversalInfo )malloc(sizeof(traversalInfo) tr->mxtips);
429	}
430	}
431	#endif
432
433
434	tr->constraintVector = (int )malloc((2 tr->mxtips) * sizeof(int));
435
436	tr->nameList = (char *)malloc(sizeof(char ) * (tips + 1));
437
438	if (!(p0 = (nodeptr) malloc((tips + 3inter) sizeof(node))))
439	{
440	printf("ERROR: Unable to obtain sufficient tree memory\n");
441	return FALSE;
442	}
443
444	if (!(tr->nodep = (nodeptr ) malloc((2tr->mxtips) * sizeof(nodeptr))))
445	{
446	printf("ERROR: Unable to obtain sufficient tree memory, too\n");
447	return FALSE;
448	}
449
450	tr->nodep[0] = (node ) NULL; / Use as 1-based array */
451
452	for (i = 1; i <= tips; i++)
453	{
454	p = p0++;
455	p->x = 0;
456	p->number = i;
457	p->next = p;
458	p->back = (node *) NULL;
459	#ifdef _MULTI_GENE
460	{
461	int k;
462
463	for(k = 0; k < tr->numBranches; k++)
464	{
465	p->xs[k] = 0;
466	p->backs[k] = (nodeptr)NULL;
467	}
468	}
469	#endif
470	tr->nodep[i] = p;
471	}
472
473	for (i = tips + 1; i <= tips + inter; i++)
474	{
475	q = (node *) NULL;
476	for (j = 1; j <= 3; j++)
477	{
478	p = p0++;
479	p->x = 0;
480	p->number = i;
481	p->next = q;
482
483	p->back = (node *) NULL;
484	#ifdef _MULTI_GENE
485	{
486	int k;
487
488	for(k = 0; k < tr->numBranches; k++)
489	{
490	p->xs[k] = 0;
491	p->backs[k] = (nodeptr)NULL;
492	}
493	}
494	#endif
495	q = p;
496	}
497	p->next->next->next = p;
498	tr->nodep[i] = p;
499	}
500
501	tr->likelihood = unlikely;
502	tr->start = (node *) NULL;
503	tr->ntips = 0;
504	tr->nextnode = 0;
505	tr->prelabeled = TRUE;
506	tr->smoothed = FALSE;
507
508	return TRUE;
509	}
510
511
512
513
514	static boolean getdata (analdef adef, rawdata rdta, tree *tr)
515	{
516	int i, j, basesread, basesnew, ch, my_i, meaning;
517	int meaningAA[256], meaningDNA[256];
518	boolean allread, firstpass;
519	char buffer[300];
520	int len;
521	unsigned long total = 0;
522	unsigned long gaps = 0;
523	int gapValueAA, gapValueDNA;
524
525	for (i = 0; i <= 255; i++)
526	{
527	meaningAA[i] = -1;
528	meaningDNA[i] = -1;
529	}
530
531	/* AA data */
532
533	meaningAA['A'] = 0; /* alanine */
534	meaningAA['R'] = 1; /* arginine */
535	meaningAA['N'] = 2; /* asparagine*/
536	meaningAA['D'] = 3; /* aspartic */
537	meaningAA['C'] = 4; /* cysteine */
538	meaningAA['Q'] = 5; /* glutamine */
539	meaningAA['E'] = 6; /* glutamic */
540	meaningAA['G'] = 7; /* glycine */
541	meaningAA['H'] = 8; /* histidine */
542	meaningAA['I'] = 9; /* isoleucine */
543	meaningAA['L'] = 10; /* leucine */
544	meaningAA['K'] = 11; /* lysine */
545	meaningAA['M'] = 12; /* methionine */
546	meaningAA['F'] = 13; /* phenylalanine */
547	meaningAA['P'] = 14; /* proline */
548	meaningAA['S'] = 15; /* serine */
549	meaningAA['T'] = 16; /* threonine */
550	meaningAA['W'] = 17; /* tryptophan */
551	meaningAA['Y'] = 18; /* tyrosine */
552	meaningAA['V'] = 19; /* valine */
553	meaningAA['B'] = 20;/* asparagine, aspartic 2 and 3*/
554	meaningAA['Z'] = 21;/21 glutamine glutamic 5 and 6/
555	meaningAA['X'] = meaningAA['?'] = meaningAA[''] = meaningAA['-'] = 22; / all = 1.0 */
556	gapValueAA = 22;
557
558	/* DNA data */
559
560	meaningDNA['A'] = 1;
561	meaningDNA['B'] = 14;
562	meaningDNA['C'] = 2;
563	meaningDNA['D'] = 13;
564	meaningDNA['G'] = 4;
565	meaningDNA['H'] = 11;
566	meaningDNA['K'] = 12;
567	meaningDNA['M'] = 3;
568	meaningDNA['N'] = 15;
569	meaningDNA['O'] = 15;
570	meaningDNA['R'] = 5;
571	meaningDNA['S'] = 6;
572	meaningDNA['T'] = 8;
573	meaningDNA['U'] = 8;
574	meaningDNA['V'] = 7;
575	meaningDNA['W'] = 9;
576	meaningDNA['X'] = 15;
577	meaningDNA['Y'] = 10;
578	meaningDNA['-'] = 15;
579	meaningDNA['?'] = 15;
580	gapValueDNA = 15;
581
582	/*******************************************************************/
583
584	basesread = basesnew = 0;
585
586	allread = FALSE;
587	firstpass = TRUE;
588	ch = ' ';
589
590	while (! allread)
591	{
592	for (i = 1; i <= tr->mxtips; i++)
593	{
594	if (firstpass)
595	{
596	ch = getc(INFILE);
597	while(ch == ' ' \|\| ch == '\n' \|\| ch == '\t' \|\| ch == '\r') /* PC-LINEBREAK*/
598	{
599	ch = getc(INFILE);
600	}
601	my_i = 0;
602
603	do
604	{
605	buffer[my_i] = ch;
606	ch = getc(INFILE);
607	my_i++;
608	if(my_i >= nmlngth)
609	{
610	if(processID == 0)
611	{
612	printf("Taxon Name to long at taxon %d, adapt constant nmlngth in\n", i);
613	printf("axml.h, current setting %d\n", nmlngth);
614	}
615	errorExit(-1);
616	}
617	}
618	while(ch != ' ' && ch != '\n' && ch != '\t' && ch != '\r'); /* PC-LINEBREAK*/
619
620	buffer[my_i] = '\0';
621	len = strlen(buffer) + 1;
622	tr->nameList[i] = (char )malloc(sizeof(char) len);
623	strcpy(tr->nameList[i], buffer);
624	}
625
626	j = basesread;
627	while ((j < rdta->sites) && ((ch = getc(INFILE)) != EOF) && (ch != '\n') && (ch != '\r')) /* PC-LINEBREAK*/
628	{
629	uppercase(& ch);
630
631	assert(tr->dataVector[j + 1] != -1);
632
633	switch(tr->dataVector[j + 1])
634	{
635	case DNA_DATA:
636	meaning = meaningDNA[ch];
637	break;
638	case AA_DATA:
639	meaning = meaningAA[ch];
640	break;
641	default:
642	assert(0);
643	}
644
645	if ((meaning != -1) \|\| ch == '.')
646	{
647	j++;
648	if (ch == '.')
649	{
650	if (i != 1)
651	ch = rdta->y[1][j];
652	else
653	{
654	printf("ERROR: Dot (.) found at site %d of sequence 1\n", j + 1);
655	return FALSE;
656	}
657	}
658	rdta->y[i][j] = ch;
659	}
660	else
661	{
662	if(whitechar(ch) \|\| digitchar(ch)) ;
663	else
664	{
665	printf("ERROR: Bad base (%c) at site %d of sequence %d\n",
666	ch, j + 1, i);
667	return FALSE;
668	}
669	}
670	}
671
672	if (ch == EOF)
673	{
674	printf("ERROR: End-of-file at site %d of sequence %d\n", j + 1, i);
675	return FALSE;
676	}
677
678	if (! firstpass && (j == basesread))
679	i--;
680	else
681	{
682	if (i == 1)
683	basesnew = j;
684	else
685	if (j != basesnew)
686	{
687	printf("ERROR: Sequences out of alignment\n");
688	printf("%d (instead of %d) residues read in sequence %d %s\n",
689	j - basesread, basesnew - basesread, i, tr->nameList[i]);
690	return FALSE;
691	}
692	}
693	while (ch != '\n' && ch != EOF && ch != '\r') ch = getc(INFILE); /* flush line // PC-LINEBREAK*/
694	}
695
696	firstpass = FALSE;
697	basesread = basesnew;
698	allread = (basesread >= rdta->sites);
699	}
700
701	for(j = 1; j <= tr->mxtips; j++)
702	for(i = 1; i <= rdta->sites; i++)
703	{
704	assert(tr->dataVector[i] != -1);
705
706	switch(tr->dataVector[i])
707	{
708	case DNA_DATA:
709	meaning = meaningDNA[rdta->y[j][i]];
710	if(meaning == gapValueDNA)
711	gaps++;
712	break;
713	case AA_DATA:
714	meaning = meaningAA[rdta->y[j][i]];
715	if(meaning == gapValueAA)
716	gaps++;
717	break;
718	default:
719	assert(0);
720	}
721
722	total++;
723	rdta->y[j][i] = meaning;
724	}
725
726	adef->gapyness = (double)gaps / (double)total;
727
728	return TRUE;
729	}
730
731
732
733	static void inputweights (rawdata *rdta)
734	{
735	int i, w, fres;
736	FILE *weightFile;
737	int wv = (int )malloc(sizeof(int) * rdta->sites + 1);
738
739	weightFile = fopen(weightFileName, "r");
740	if (!weightFile)
741	{
742	if(processID == 0)
743	printf( "Could not open weight file: %s\n", weightFileName);
744	errorExit(-1);
745	}
746
747	i = 1;
748
749	while((fres = fscanf(weightFile,"%d", &w)) != EOF)
750	{
751	if(!fres)
752	{
753	if(processID == 0)
754	printf("error reading weight file probably encountered a non-integer weight value\n");
755	errorExit(-1);
756	}
757	wv[i] = w;
758	i++;
759	}
760
761	if(i != (rdta->sites + 1))
762	{
763	if(processID == 0)
764	printf("number %d of weights not equal to number %d of alignment columns\n", i, rdta->sites);
765	errorExit(-1);
766	}
767
768	for(i = 1; i <= rdta->sites; i++)
769	rdta->wgt[i] = wv[i];
770
771	fclose(weightFile);
772	free(wv);
773	}
774
775
776
777	static void getinput(analdef adef, rawdata rdta, cruncheddata cdta, tree tr)
778	{
779	int i;
780
781	getnums(rdta);
782
783	tr->mxtips = rdta->numsp;
784	rdta->wgt = (int ) malloc((rdta->sites + 1) sizeof(int));
785	rdta->wgt2 = (int ) malloc((rdta->sites + 1) sizeof(int));
786	cdta->alias = (int ) malloc((rdta->sites + 1) sizeof(int));
787	cdta->aliaswgt = (int ) malloc((rdta->sites + 1) sizeof(int));
788	cdta->rateCategory = (int ) malloc((rdta->sites + 1) sizeof(int));
789	tr->model = (int *) calloc((rdta->sites + 1), sizeof(int));
790	tr->dataVector = (int ) malloc((rdta->sites + 1) sizeof(int));
791	cdta->wr = (double ) malloc((rdta->sites + 1) sizeof(double));
792	cdta->wr2 = (double ) malloc((rdta->sites + 1) sizeof(double));
793	cdta->patrat = (double ) malloc((rdta->sites + 1) sizeof(double));
794	cdta->patratStored = (double ) malloc((rdta->sites + 1) sizeof(double));
795
796
797	if(!adef->useWeightFile)
798	{
799	for (i = 1; i <= rdta->sites; i++)
800	rdta->wgt[i] = 1;
801	}
802	else
803	inputweights(rdta);
804
805	tr->multiBranch = 0;
806	tr->numBranches = 1;
807
808	if(adef->useMultipleModel)
809	{
810	int ref;
811
812	parsePartitions(adef, rdta, tr);
813
814	for(i = 1; i <= rdta->sites; i++)
815	{
816	ref = tr->model[i];
817	tr->dataVector[i] = tr->partitionData[ref].dataType;
818	}
819	}
820	else
821	{
822	int dataType;
823
824	tr->partitionData = (pInfo*)malloc(sizeof(pInfo));
825	tr->partitionData[0].partitionName = (char)malloc(128 sizeof(char));
826	strcpy(tr->partitionData[0].partitionName, "No Name Provided");
827
828	tr->partitionData[0].protModels = adef->proteinMatrix;
829	tr->partitionData[0].protFreqs = adef->protEmpiricalFreqs;
830
831
832	tr->NumberOfModels = 1;
833
834	if(adef->model == M_PROTCAT \|\| adef->model == M_PROTGAMMA)
835	dataType = AA_DATA;
836	else
837	dataType = DNA_DATA;
838
839	/* INIT data-type, model, dataVector for good */
840	/* those values will be constant throughout the */
841	/* inference process */
842
843	tr->partitionData[0].dataType = dataType;
844
845	for(i = 0; i <= rdta->sites; i++)
846	{
847	tr->dataVector[i] = dataType;
848	tr->model[i] = 0;
849	}
850	}
851
852	#ifdef _MULTI_GENE
853	{
854	int i;
855
856	tr->startVector = (nodeptr )malloc(sizeof(nodeptr) tr->NumberOfModels);
857	tr->tipMissing = (char *)malloc(sizeof(char ) * (tr->mxtips + 1));
858
859	for(i = 0; i <= tr->mxtips; i++)
860	tr->tipMissing[i] = (char )malloc(sizeof(char) (tr->NumberOfModels));
861	}
862	#endif
863
864
865	getyspace(rdta);
866	setupTree(tr, adef);
867
868	if(!getdata(adef, rdta, tr))
869	{
870	printf("Problem reading alignment file \n");
871	errorExit(1);
872	}
873
874	return;
875	}
876
877
878
879
880	static void sitesort(rawdata rdta, cruncheddata cdta, tree tr, analdef adef)
881	{
882	int gap, i, j, jj, jg, k, n, nsp;
883	int index, category, *superCategory;
884	boolean flip, tied;
885	char **data;
886
887	if(adef->useMultipleModel)
888	{
889	superCategory = tr->dataVector;
890	category = tr->model;
891	}
892	else
893	{
894	category = (int*)NULL;
895	superCategory = (int*)NULL;
896	}
897
898	index = cdta->alias;
899	data = rdta->y;
900	n = rdta->sites;
901	nsp = rdta->numsp;
902	index[0] = -1;
903
904
905	if((adef->mode != OPTIMIZE_RATES))
906	{
907	for (gap = n / 2; gap > 0; gap /= 2)
908	{
909	for (i = gap + 1; i <= n; i++)
910	{
911	j = i - gap;
912
913	do
914	{
915	jj = index[j];
916	jg = index[j+gap];
917	if(adef->useMultipleModel)
918	{
919	assert(superCategory[jj] != -1 &&
920	superCategory[jg] != -1 &&
921	category[jj] != -1 &&
922	category[jg] != -1);
923
924	if(superCategory[jj] > superCategory[jg])
925	{
926	flip = TRUE;
927	tied = FALSE;
928	}
929	else
930	{
931	flip = (category[jj] > category[jg]);
932	tied = (category[jj] == category[jg]);
933	}
934	}
935	else
936	{
937	flip = 0;
938	tied = 1;
939	}
940
941	for (k = 1; (k <= nsp) && tied; k++)
942	{
943	flip = (data[k][jj] > data[k][jg]);
944	tied = (data[k][jj] == data[k][jg]);
945	}
946
947	if (flip)
948	{
949	index[j] = jg;
950	index[j+gap] = jj;
951	j -= gap;
952	}
953	}
954	while (flip && (j > 0));
955	}
956	}
957	}
958	}
959
960
961	static void sitecombcrunch (rawdata rdta, cruncheddata cdta, tree tr, analdef adef)
962	{
963	int i, sitei, j, sitej, k;
964	boolean tied;
965	int *aliasModel;
966	int *aliasSuperModel;
967
968
969
970	if(adef->useMultipleModel)
971	{
972	aliasSuperModel = (int)malloc(sizeof(int) (rdta->sites + 1));
973	aliasModel = (int)malloc(sizeof(int) (rdta->sites + 1));
974	}
975	else
976	{
977	aliasModel = (int*)NULL;
978	aliasSuperModel = (int*)NULL;
979	}
980
981	i = 0;
982	cdta->alias[0] = cdta->alias[1];
983	cdta->aliaswgt[0] = 0;
984
985	for (j = 1; j <= rdta->sites; j++)
986	{
987	sitei = cdta->alias[i];
988	sitej = cdta->alias[j];
989	if(adef->mode == OPTIMIZE_RATES)
990	tied = 0;
991	else
992	{
993	if(adef->useMultipleModel)
994	{
995	tied = (tr->model[sitei] == tr->model[sitej]);
996	if(tied)
997	assert(tr->dataVector[sitei] == tr->dataVector[sitej]);
998	}
999	else
1000	tied = 1;
1001	}
1002
1003	for (k = 1; tied && (k <= rdta->numsp); k++)
1004	tied = (rdta->y[k][sitei] == rdta->y[k][sitej]);
1005
1006	if (tied)
1007	{
1008	cdta->aliaswgt[i] += rdta->wgt2[sitej];
1009	if(adef->useMultipleModel)
1010	{
1011	aliasModel[i] = tr->model[sitej];
1012	aliasSuperModel[i] = tr->dataVector[sitej];
1013	}
1014	}
1015	else
1016	{
1017	if (cdta->aliaswgt[i] > 0) i++;
1018	cdta->aliaswgt[i] = rdta->wgt2[sitej];
1019	cdta->alias[i] = sitej;
1020	if(adef->useMultipleModel)
1021	{
1022	aliasModel[i] = tr->model[sitej];
1023	aliasSuperModel[i] = tr->dataVector[sitej];
1024	}
1025	}
1026	}
1027
1028	cdta->endsite = i;
1029	if (cdta->aliaswgt[i] > 0) cdta->endsite++;
1030
1031	if(adef->useMultipleModel)
1032	{
1033	for(i = 0; i <= rdta->sites; i++)
1034	{
1035	tr->model[i] = aliasModel[i];
1036	tr->dataVector[i] = aliasSuperModel[i];
1037	}
1038	}
1039
1040	if(adef->useMultipleModel)
1041	{
1042	free(aliasModel);
1043	free(aliasSuperModel);
1044	}
1045
1046	/*for(i = 0; i < cdta->endsite; i++)
1047	{
1048	printf("%d %d %d\n", i, tr->model[i], tr->dataVector[i]);
1049	}*/
1050	}
1051
1052
1053	static boolean makeweights (analdef adef, rawdata rdta, cruncheddata cdta, tree tr)
1054	{
1055	int i;
1056
1057
1058	for (i = 1; i <= rdta->sites; i++)
1059	rdta->wgt2[i] = rdta->wgt[i];
1060
1061	for (i = 1; i <= rdta->sites; i++)
1062	cdta->alias[i] = i;
1063
1064	sitesort(rdta, cdta, tr, adef);
1065	sitecombcrunch(rdta, cdta, tr, adef);
1066
1067	return TRUE;
1068	}
1069
1070
1071
1072
1073	static boolean makevalues(rawdata rdta, cruncheddata cdta, tree tr, analdef adef)
1074	{
1075	int i, j, model, fullSites = 0, modelCounter;
1076
1077	char y = (char )malloc(rdta->numsp * cdta->endsite * sizeof(char));
1078	char yBUF = (char )malloc(rdta->numsp * cdta->endsite * sizeof(char));
1079
1080	for (i = 1; i <= rdta->numsp; i++)
1081	for (j = 0; j < cdta->endsite; j++)
1082	y[((i - 1) * cdta->endsite) + j] = rdta->y[i][cdta->alias[j]];
1083
1084	free(rdta->y0);
1085	free(rdta->y);
1086
1087	rdta->y0 = y;
1088	memcpy(yBUF, y, rdta->numsp * cdta->endsite * sizeof(char));
1089	rdta->yBUF = yBUF;
1090
1091	if(!adef->useMultipleModel)
1092	tr->NumberOfModels = 1;
1093
1094	if(adef->useMultipleModel)
1095	{
1096	int perm = (int)malloc(sizeof(int) * tr->NumberOfModels);
1097	pInfo partitionData = (pInfo)malloc(sizeof(pInfo) * tr->NumberOfModels);
1098
1099	tr->partitionData[0].lower = 0;
1100
1101	model = tr->model[0];
1102	modelCounter = 0;
1103	perm[modelCounter] = model;
1104	i = 1;
1105
1106	while(i < cdta->endsite)
1107	{
1108	if(tr->model[i] != model)
1109	{
1110	tr->partitionData[modelCounter].upper = i;
1111	tr->partitionData[modelCounter + 1].lower = i;
1112
1113	model = tr->model[i];
1114	perm[modelCounter + 1] = model;
1115	modelCounter++;
1116	}
1117	i++;
1118	}
1119
1120
1121	tr->partitionData[tr->NumberOfModels - 1].upper = cdta->endsite;
1122
1123	memcpy(partitionData, tr->partitionData, tr->NumberOfModels * sizeof(pInfo));
1124
1125	for(i = 0; i < tr->NumberOfModels; i++)
1126	{
1127	tr->partitionData[i].dataType = partitionData[perm[i]].dataType;
1128	tr->partitionData[i].protModels = partitionData[perm[i]].protModels;
1129	tr->partitionData[i].protFreqs = partitionData[perm[i]].protFreqs;
1130	}
1131
1132	model = tr->model[0];
1133	modelCounter = 0;
1134	tr->model[0] = modelCounter;
1135	i = 1;
1136
1137	while(i < cdta->endsite)
1138	{
1139	if(tr->model[i] != model)
1140	{
1141	model = tr->model[i];
1142	modelCounter++;
1143	tr->model[i] = modelCounter;
1144	}
1145	else
1146	tr->model[i] = modelCounter;
1147	i++;
1148	}
1149
1150	for(i = 0; i < (cdta->endsite - 1); i++)
1151	{
1152	if(tr->dataVector[i] != tr->dataVector[i + 1])
1153	{
1154	tr->mixedData = TRUE;
1155	break;
1156	}
1157	}
1158
1159	free(perm);
1160	free(partitionData);
1161	}
1162	else
1163	{
1164	tr->partitionData[0].lower = 0;
1165	tr->partitionData[0].upper = cdta->endsite;
1166	}
1167
1168	#ifdef _USE_PTHREADS
1169	/*
1170	TODO-MIX Have to seriously think about whether
1171	to implement this or not
1172	*/
1173
1174	if(tr->mixedData)
1175	{
1176	printf("No Pthreads implementation for mixed data available right now \n");
1177	assert(0);
1178	}
1179	#endif
1180
1181
1182	tr->rdta = rdta;
1183	tr->cdta = cdta;
1184
1185	tr->invariant = (int )malloc(cdta->endsite sizeof(int));
1186	tr->originalDataVector = (int )malloc(cdta->endsite sizeof(int));
1187	tr->originalModel = (int )malloc(cdta->endsite sizeof(int));
1188	tr->originalWeights = (int )malloc(cdta->endsite sizeof(int));
1189
1190	memcpy(tr->originalModel, tr->model, cdta->endsite * sizeof(int));
1191	memcpy(tr->originalDataVector, tr->dataVector, cdta->endsite * sizeof(int));
1192	memcpy(tr->originalWeights, tr->cdta->aliaswgt, cdta->endsite * sizeof(int));
1193
1194	tr->originalCrunchedLength = tr->cdta->endsite;
1195	for(i = 0; i < tr->cdta->endsite; i++)
1196	fullSites += tr->cdta->aliaswgt[i];
1197
1198	tr->fullSites = fullSites;
1199
1200	for(i = 0; i < rdta->numsp; i++)
1201	tr->yVector[i + 1] = &(rdta->y0[tr->originalCrunchedLength * i]);
1202
1203	return TRUE;
1204	}
1205
1206
1207
1208
1209
1210
1211
1212	static int sequenceSimilarity(char tipJ, char tipK, int n)
1213	{
1214	int i;
1215
1216	for(i = 0; i < n; i++)
1217	if(tipJ++ != tipK++)
1218	return 0;
1219
1220	return 1;
1221	}
1222
1223	static void checkSequences(tree tr, rawdata rdta, analdef *adef)
1224	{
1225	int n = tr->mxtips + 1;
1226	int i, j;
1227	int omissionList = (int )malloc(n * sizeof(int));
1228	int undeterminedList = (int )malloc((rdta->sites + 1)* sizeof(int));
1229	int modelList = (int )malloc((rdta->sites + 1)* sizeof(int));
1230	int count = 0;
1231	int countNameDuplicates = 0;
1232	int countUndeterminedColumns = 0;
1233	int countOnlyGaps = 0;
1234	int modelCounter = 1;
1235	char undetermined_AA = 22;
1236	char undetermined_DNA = 15;
1237	char tipI, tipJ;
1238	FILE *f;
1239
1240
1241	if(processID == 0)
1242	f = fopen(infoFileName, "a");
1243	else
1244	f = (FILE *)NULL;
1245
1246	for(i = 1; i < n; i++)
1247	omissionList[i] = 0;
1248
1249	for(i = 0; i < rdta->sites + 1; i++)
1250	undeterminedList[i] = 0;
1251
1252	for(i = 1; i < n; i++)
1253	{
1254	for(j = i + 1; j < n; j++)
1255	if(strcmp(tr->nameList[i], tr->nameList[j]) == 0)
1256	{
1257	countNameDuplicates++;
1258	if(processID == 0)
1259	{
1260	printf("Sequence names of taxon %d and %d are identical, they are both called %s\n", i, j, tr->nameList[i]);
1261	fprintf(f, "Sequence names of taxon %d and %d are identical, they are both called %s\n", i, j, tr->nameList[i]);
1262	}
1263	}
1264	}
1265
1266	if(countNameDuplicates > 0)
1267	{
1268	if(processID == 0)
1269	{
1270	printf("ERROR: Found %d taxa that had equal names in the alignment, exiting...\n", countNameDuplicates);
1271	fprintf(f, "ERROR: Found %d taxa that had equal names in the alignment, exiting...\n", countNameDuplicates);
1272	fclose(f);
1273	}
1274	errorExit(-1);
1275	}
1276
1277	for(i = 1; i < n; i++)
1278	{
1279	j = 1;
1280
1281	while(j <= rdta->sites)
1282	{
1283	if(tr->dataVector[j] == DNA_DATA && rdta->y[i][j] != undetermined_DNA)
1284	break;
1285	if(tr->dataVector[j] == AA_DATA && rdta->y[i][j] != undetermined_AA)
1286	break;
1287	j++;
1288	}
1289
1290	if(j == (rdta->sites + 1))
1291	{
1292	if(processID == 0)
1293	{
1294	printf("ERROR: Sequence %s consists entirely of undetermined values which will be treated as missing data\n",
1295	tr->nameList[i]);
1296	fprintf(f, "ERROR: Sequence %s consists entirely of undetermined values which will be treated as missing data\n",
1297	tr->nameList[i]);
1298	}
1299	countOnlyGaps++;
1300	}
1301
1302	}
1303
1304	if(countOnlyGaps > 0)
1305	{
1306	if(processID == 0)
1307	{
1308	printf("ERROR: Found %d sequences that consist entirely of undetermined values, exiting...\n", countOnlyGaps);
1309	fprintf(f, "ERROR: Found %d sequences that consist entirely of undetermined values, exiting...\n", countOnlyGaps);
1310	fclose(f);
1311	}
1312	errorExit(-1);
1313	}
1314
1315	for(i = 0; i <= rdta->sites; i++)
1316	modelList[i] = -1;
1317
1318	for(i = 1; i <= rdta->sites; i++)
1319	{
1320	j = 1;
1321
1322	while(j < n)
1323	{
1324	if(tr->dataVector[i] == DNA_DATA && rdta->y[j][i] != undetermined_DNA)
1325	break;
1326	if(tr->dataVector[i] == AA_DATA && rdta->y[j][i] != undetermined_AA)
1327	break;
1328	j++;
1329	}
1330
1331	if(j == n)
1332	{
1333	undeterminedList[i] = 1;
1334	if(processID == 0)
1335	{
1336	printf("IMPORTANT WARNING: Alignment column %d contains only undetermined values which will be treated as missing data\n",
1337	i);
1338	fprintf(f, "IMPORTANT WARNING: Alignment column %d contains only undetermined values which will be treated as missing data\n", i);
1339	}
1340	countUndeterminedColumns++;
1341	}
1342	else
1343	{
1344	if(adef->useMultipleModel)
1345	{
1346	modelList[modelCounter] = tr->model[i];
1347	modelCounter++;
1348	}
1349	}
1350	}
1351
1352
1353	for(i = 1; i < n; i++)
1354	{
1355	if(omissionList[i] == 0)
1356	{
1357	tipI = &(rdta->y[i][1]);
1358
1359	for(j = i + 1; j < n; j++)
1360	{
1361	if(omissionList[j] == 0)
1362	{
1363	tipJ = &(rdta->y[j][1]);
1364	if(sequenceSimilarity(tipI, tipJ, rdta->sites))
1365	{
1366	if(processID == 0)
1367	{
1368	printf("\n\nIMPORTANT WARNING: Sequences %s and %s are exactly identical\n", tr->nameList[i], tr->nameList[j]);
1369	fprintf(f, "\n\nIMPORTANT WARNING: Sequences %s and %s are exactly identical\n", tr->nameList[i], tr->nameList[j]);
1370	}
1371	omissionList[j] = 1;
1372	count++;
1373	}
1374	}
1375	}
1376	}
1377	}
1378
1379	if(count > 0 \|\| countUndeterminedColumns > 0)
1380	{
1381	char noDupFile[2048];
1382	char noDupModels[2048];
1383
1384	if(count > 0)
1385	{
1386	if(processID == 0)
1387	{
1388	printf("\n");
1389
1390	printf("IMPORTANT WARNING\n");
1391
1392	printf("Found %d %s that %s exactly identical to other sequences in the alignment.\n", count, (count == 1)?"sequence":"sequences", (count == 1)?"is":"are");
1393	printf("Normally they should be excluded from the analysis.\n\n");
1394
1395	fprintf(f, "\n");
1396
1397	fprintf(f, "IMPORTANT WARNING\n");
1398
1399	fprintf(f, "Found %d %s that %s exactly identical to other sequences in the alignment.\n", count, (count == 1)?"sequence":"sequences", (count == 1)?"is":"are");
1400	fprintf(f, "Normally they should be excluded from the analysis.\n\n");
1401	}
1402	}
1403
1404	if(countUndeterminedColumns > 0)
1405	{
1406	if(processID == 0)
1407	{
1408	printf("\n");
1409
1410	printf("IMPORTANT WARNING\n");
1411
1412	printf("Found %d %s that %s only undetermined values which will be treated as missing data.\n",
1413	countUndeterminedColumns, (countUndeterminedColumns == 1)?"column":"columns", (countUndeterminedColumns == 1)?"contains":"contain");
1414	printf("Normally these columns should be excluded from the analysis.\n\n");
1415
1416	fprintf(f, "\n");
1417
1418	fprintf(f, "IMPORTANT WARNING\n");
1419
1420	fprintf(f, "Found %d %s that %s only undetermined values which will be treated as missing data.\n",
1421	countUndeterminedColumns, (countUndeterminedColumns == 1)?"column":"columns", (countUndeterminedColumns == 1)?"contains":"contain");
1422	fprintf(f, "Normally these columns should be excluded from the analysis.\n\n");
1423	}
1424	}
1425
1426	strcpy(noDupFile, seq_file);
1427	strcat(noDupFile, ".reduced");
1428
1429	strcpy(noDupModels, modelFileName);
1430	strcat(noDupModels, ".reduced");
1431
1432	if(processID == 0)
1433	{
1434
1435	if(adef->useMultipleModel && !filexists(noDupModels) && countUndeterminedColumns)
1436	{
1437	FILE *newFile = fopen(noDupModels, "w");
1438
1439	printf("\nJust in case you might need it, a mixed model file with \n");
1440	printf("model assignments for undetermined columns removed is printed to file %s\n",noDupModels);
1441
1442	fprintf(f, "\nJust in case you might need it, a mixed model file with \n");
1443	fprintf(f, "model assignments for undetermined columns removed is printed to file %s\n",noDupModels);
1444
1445
1446	for(i = 0; i < tr->NumberOfModels; i++)
1447	{
1448	boolean modelStillExists = FALSE;
1449
1450	for(j = 1; (j <= rdta->sites) && (!modelStillExists); j++)
1451	{
1452	if(modelList[j] == i)
1453	modelStillExists = TRUE;
1454	}
1455
1456	if(modelStillExists)
1457	{
1458	char *protModels[10] = {"DAYHOFF", "DCMUT", "JTT", "MTREV", "WAG", "RTREV", "CPREV", "VT", "BLOSUM62", "MTMAM"};
1459	int k = 1;
1460	int lower, upper;
1461	int parts = 0;
1462
1463
1464	switch(tr->partitionData[i].dataType)
1465	{
1466	case AA_DATA:
1467	{
1468	char AAmodel[1024];
1469
1470	strcpy(AAmodel, protModels[tr->partitionData[i].protModels]);
1471	if(tr->partitionData[i].protFreqs)
1472	strcat(AAmodel, "F");
1473
1474	fprintf(newFile, "%s, ", AAmodel);
1475	}
1476	break;
1477	case DNA_DATA:
1478	fprintf(newFile, "DNA, ");
1479	break;
1480	default:
1481	assert(0);
1482	}
1483
1484	fprintf(newFile, "%s = ", tr->partitionData[i].partitionName);
1485
1486	while(k <= rdta->sites)
1487	{
1488	if(modelList[k] == i)
1489	{
1490	lower = k;
1491	while((modelList[k + 1] == i) && (k <= rdta->sites))
1492	k++;
1493	upper = k;
1494
1495	if(lower == upper)
1496	{
1497	if(parts == 0)
1498	fprintf(newFile, "%d", lower);
1499	else
1500	fprintf(newFile, ",%d", lower);
1501	}
1502	else
1503	{
1504	if(parts == 0)
1505	fprintf(newFile, "%d-%d", lower, upper);
1506	else
1507	fprintf(newFile, ",%d-%d", lower, upper);
1508	}
1509	parts++;
1510	}
1511	k++;
1512	}
1513	fprintf(newFile, "\n");
1514	}
1515	}
1516	fclose(newFile);
1517	}
1518	else
1519	{
1520	if(adef->useMultipleModel)
1521	{
1522	printf("\n A mixed model file with model assignments for undetermined\n");
1523	printf("columns removed has already been printed to file %s\n",noDupModels);
1524
1525	fprintf(f, "\n A mixed model file with model assignments for undetermined\n");
1526	fprintf(f, "columns removed has already been printed to file %s\n",noDupModels);
1527	}
1528	}
1529
1530
1531	if(!filexists(noDupFile))
1532	{
1533	FILE *newFile;
1534
1535	printf("Just in case you might need it, an alignment file with \n");
1536	if(count && !countUndeterminedColumns)
1537	printf("sequence duplicates removed is printed to file %s\n", noDupFile);
1538	if(!count && countUndeterminedColumns)
1539	printf("undetermined columns removed is printed to file %s\n", noDupFile);
1540	if(count && countUndeterminedColumns)
1541	printf("sequence duplicates and undetermined columns removed is printed to file %s\n", noDupFile);
1542
1543	fprintf(f, "Just in case you might need it, an alignment file with \n");
1544	if(count && !countUndeterminedColumns)
1545	fprintf(f, "sequence duplicates removed is printed to file %s\n", noDupFile);
1546	if(!count && countUndeterminedColumns)
1547	fprintf(f, "undetermined columns removed is printed to file %s\n", noDupFile);
1548	if(count && countUndeterminedColumns)
1549	fprintf(f, "sequence duplicates and undetermined columns removed is printed to file %s\n", noDupFile);
1550
1551	newFile = fopen(noDupFile, "w");
1552
1553	fprintf(newFile, "%d %d\n", tr->mxtips - count, rdta->sites - countUndeterminedColumns);
1554
1555	for(i = 1; i < n; i++)
1556	{
1557	if(!omissionList[i])
1558	{
1559	fprintf(newFile, "%s ", tr->nameList[i]);
1560	tipI = &(rdta->y[i][1]);
1561
1562	for(j = 0; j < rdta->sites; j++)
1563	{
1564	if(undeterminedList[j + 1] == 0)
1565	{
1566	switch(tr->dataVector[j + 1])
1567	{
1568	case AA_DATA:
1569	fprintf(newFile, "%c", inverseMeaningPROT[tipI[j]]);
1570	break;
1571	case DNA_DATA:
1572	fprintf(newFile, "%c", inverseMeaningDNA[tipI[j]]);
1573	break;
1574	default:
1575	assert(0);
1576	}
1577	}
1578	}
1579
1580	fprintf(newFile, "\n");
1581	}
1582	}
1583
1584	fclose(newFile);
1585	}
1586	else
1587	{
1588	if(count && !countUndeterminedColumns)
1589	printf("An alignment file with sequence duplicates removed has already\n");
1590	if(!count && countUndeterminedColumns)
1591	printf("An alignment file with undetermined columns removed has already\n");
1592	if(count && countUndeterminedColumns)
1593	printf("An alignment file with undetermined columns and sequence duplicates removed has already\n");
1594
1595	printf("been printed to file %s\n", noDupFile);
1596
1597	if(count && !countUndeterminedColumns)
1598	fprintf(f, "An alignment file with sequence duplicates removed has already\n");
1599	if(!count && countUndeterminedColumns)
1600	fprintf(f, "An alignment file with undetermined columns removed has already\n");
1601	if(count && countUndeterminedColumns)
1602	fprintf(f, "An alignment file with undetermined columns and sequence duplicates removed has already\n");
1603
1604	fprintf(f, "been printed to file %s\n", noDupFile);
1605	}
1606	}
1607	}
1608
1609
1610	free(undeterminedList);
1611	free(omissionList);
1612	free(modelList);
1613	if(processID == 0)
1614	fclose(f);
1615	}
1616
1617
1618
1619	static float dist(int i, int j, const int sites, const float nDouble, char **y)
1620	{
1621	int k, count;
1622	char *tipI = &(y[i + 1][1]);
1623	char *tipJ = &(y[j + 1][1]);
1624
1625	for(k = 0, count = 0; k < sites; k ++)
1626	if(tipI[k] == tipJ[k])
1627	count++;
1628
1629	return (((float)count) * nDouble);
1630	}
1631
1632	static void distArray(int i, const int sites, const float nDouble, int n, float ref, int omitted, char **y)
1633	{
1634	int k, l;
1635	char *tipI = &(y[i + 1][1]);
1636
1637	for(l = 0; l < n; l++)
1638	{
1639	if((!omitted[l]) && (l != i))
1640	{
1641	char *tipJ = &(y[l + 1][1]);
1642	int count = 0;
1643	for(k = 0, count = 0; k < sites; k ++)
1644	if(tipI[k] == tipJ[k])
1645	count++;
1646	ref[l] = (((float)count) * nDouble);
1647	}
1648	}
1649	}
1650
1651
1652
1653	static int qtCompare(const void p1, const void p2)
1654	{
1655	qtData rc1 = (qtData )p1;
1656	qtData rc2 = (qtData )p2;
1657
1658	float i = rc1->val;
1659	float j = rc2->val;
1660
1661	if (i > j)
1662	return (-1);
1663	if (i < j)
1664	return (1);
1665	return (0);
1666	}
1667
1668
1669	static qtList * clusterQT_LARGE(int n, float thres, int ccc, rawdata rdta)
1670	{
1671	int clusterCount;
1672	int i, j;
1673	int omitted, current, *best;
1674	qtList clusters = (qtList)NULL;
1675	const float nDouble = 1.0 / (float)(rdta->sites);
1676	double t = gettime();
1677	const int sites = rdta->sites;
1678	char **y = rdta->y;
1679	float *ref;
1680	qtData *candidates;
1681
1682	candidates = (qtData )malloc(sizeof(qtData) n);
1683	clusters = (qtList )malloc(sizeof(qtList) n);
1684	omitted = (int*)calloc(n, sizeof(int));
1685	current = (int)malloc(sizeof(int) n);
1686	best = (int)malloc(sizeof(int) n);
1687	ref = (float)malloc(sizeof(float) n);
1688	clusterCount = 0;
1689
1690
1691	for(i = 0; i < n; i++)
1692	{
1693	if(!omitted[i])
1694	{
1695	int entCount = 0;
1696	int countCandidates = 0;
1697
1698	current[entCount++] = i;
1699	omitted[i] = 1;
1700
1701	distArray(i, sites, nDouble, n, ref, omitted, y);
1702
1703	for(j = 0; j < n; j++)
1704	{
1705	if(!omitted[j] && i != j)
1706	{
1707	float temp;
1708
1709	if((temp = ref[j]) >= thres)
1710	{
1711	candidates[countCandidates].val = temp;
1712	candidates[countCandidates].number = j;
1713	countCandidates++;
1714	}
1715	}
1716	}
1717
1718	if(countCandidates > 0)
1719	{
1720	qsort(candidates, countCandidates, sizeof(qtData), qtCompare);
1721
1722	for(j = 0; j < countCandidates; j++)
1723	{
1724	int k;
1725
1726	for(k = 0; k < entCount; k++)
1727	if(dist(current[k], candidates[j].number, sites, nDouble, y) < thres)
1728	break;
1729
1730	if(k == entCount)
1731	{
1732	current[entCount++] = candidates[j].number;
1733	omitted[candidates[j].number] = 1;
1734	}
1735	}
1736	}
1737
1738	clusters[clusterCount].entries = (int )malloc(sizeof(int) entCount);
1739	memcpy(clusters[clusterCount].entries, current, entCount * sizeof(int));
1740	clusters[clusterCount++].count = entCount;
1741	}
1742	}
1743
1744	printf("Time %f\n", gettime() - t);
1745	printf("FOUND %d Clusters\n", clusterCount);
1746
1747
1748	if(1)
1749	{
1750	int ver = 0;
1751	int check = 0;
1752	int total = 0;
1753	for(i = 0; i < n; i++)
1754	ver += i;
1755
1756	for(i = 0; i < clusterCount; i++)
1757	{
1758	{
1759	int k;
1760	for(j = 0; j < clusters[i].count; j++)
1761	for(k = 0; k < clusters[i].count; k++)
1762	assert(dist(clusters[i].entries[j], clusters[i].entries[k],sites, nDouble, y) >= thres);
1763	}
1764
1765	for(j = 0; j < clusters[i].count; j++)
1766	{
1767	check += clusters[i].entries[j];
1768	total++;
1769	}
1770	}
1771	assert(ver == check);
1772	printf("Total: %d\n", total);
1773	}
1774
1775
1776
1777	for(i = 0; i < clusterCount; i++)
1778	{
1779	float max = 0.0;
1780	int length = clusters[i].count;
1781	int pos = -1;
1782	int *c = clusters[i].entries;
1783	int buf;
1784
1785	if(length > 2)
1786	{
1787	for(j = 0; j < length; j++)
1788	{
1789	int k;
1790	float avg = 0.0;
1791
1792	for(k = 0; k < length; k++)
1793	{
1794	if(j != k)
1795	avg += dist(c[j], c[k], sites, nDouble, y);
1796	}
1797
1798	if(avg > max)
1799	{
1800	max = avg;
1801	pos = j;
1802	}
1803	}
1804	if(pos > 0)
1805	{
1806	buf = c[0];
1807	c[0] = c[pos];
1808	c[pos] = buf;
1809	}
1810	}
1811	for(j = 0; j < length; j++)
1812	c[j] = c[j] + 1;
1813	}
1814
1815	free(candidates);
1816	free(omitted);
1817	free(current);
1818	free(best);
1819	free(ref);
1820	*ccc = clusterCount;
1821	return clusters;
1822	}
1823
1824
1825
1826	static qtList * clusterQT(float *d, int n, float thres, int ccc)
1827	{
1828	int clusterCount;
1829	int i, j;
1830	int omitted, current, *best;
1831	int total = 0;
1832	qtList clusters = (qtList)NULL;
1833	double t = gettime();
1834
1835	clusters = (qtList )malloc(sizeof(qtList) n);
1836	omitted = (int*)calloc(n, sizeof(int));
1837	current = (int)malloc(sizeof(int) n);
1838	best = (int)malloc(sizeof(int) n);
1839
1840	clusterCount = 0;
1841
1842	while(1)
1843	{
1844	int max = -1;
1845	int maxPos = -1;
1846
1847	for(i = 0; i < n; i++)
1848	{
1849	if(!omitted[i])
1850	{
1851	int entCount = 0;
1852	int inSet = (int )calloc(n, sizeof(int));
1853	boolean aboveThres = TRUE;
1854
1855	current[entCount++] = i;
1856	inSet[i] = 1;
1857
1858	while(aboveThres)
1859	{
1860	float dm = -1.0;
1861	int dmPos = -1;
1862
1863	for(j = 0; j < n; j++)
1864	if(i != j && (!omitted[j]) && (!inSet[j]) && d[i][j] > dm)
1865	{
1866	dm = d[i][j];
1867	dmPos = j;
1868	}
1869
1870	if(dmPos == -1)
1871	aboveThres = FALSE;
1872	else
1873	{
1874	for(j = 0; j < entCount && aboveThres; j++)
1875	if(d[current[j]][dmPos] < thres)
1876	aboveThres = FALSE;
1877
1878	if(aboveThres)
1879	{
1880	current[entCount++] = dmPos;
1881	inSet[dmPos] = 1;
1882	}
1883	}
1884	}
1885
1886	if(entCount > max)
1887	{
1888	max = entCount;
1889	maxPos = i;
1890	memcpy(best, current, entCount * sizeof(int));
1891	}
1892	free(inSet);
1893	}
1894	}
1895
1896	if(maxPos == -1)
1897	break;
1898
1899	clusters[clusterCount].entries = (int )malloc(sizeof(int) max);
1900	memcpy(clusters[clusterCount].entries, best, max * sizeof(int));
1901
1902	for(i = 0; i < max; i++)
1903	omitted[best[i]] = 1;
1904
1905	clusters[clusterCount++].count = max;
1906	}
1907
1908	printf("Time %f\n", gettime() - t);
1909	printf("FOUND %d Clusters\n", clusterCount);
1910
1911	if(1)
1912	{
1913	int ver = 0;
1914	int check = 0;
1915	for(i = 0; i < n; i++)
1916	ver += i;
1917
1918	for(i = 0; i < clusterCount; i++)
1919	{
1920	/printf("Cluster %d:", i);/
1921
1922	{
1923	int k;
1924	for(j = 0; j < clusters[i].count; j++)
1925	for(k = 0; k < clusters[i].count; k++)
1926	assert(d[clusters[i].entries[j]][clusters[i].entries[k]] >= thres);
1927	}
1928
1929	for(j = 0; j < clusters[i].count; j++)
1930	{
1931	check += clusters[i].entries[j];
1932	/printf("%d ", clusters[i].entries[j]);/
1933	total++;
1934	}
1935	/printf("\n");/
1936	}
1937	assert(ver == check);
1938
1939	/printf("TOTAL: %d\n", total);/
1940	}
1941
1942	for(i = 0; i < clusterCount; i++)
1943	{
1944	float max = 0.0;
1945	int length = clusters[i].count;
1946	int pos = -1;
1947	int *c = clusters[i].entries;
1948	int buf;
1949
1950	if(length > 2)
1951	{
1952	for(j = 0; j < length; j++)
1953	{
1954	int k;
1955	float avg = 0.0;
1956	for(k = 0; k < length; k++)
1957	{
1958	if(j != k)
1959	avg += d[c[j]][c[k]];
1960	}
1961
1962	if(avg > max)
1963	{
1964	max = avg;
1965	pos = j;
1966	}
1967	}
1968	/printf("Cluster %d length %d avg %f\n", i, length, max);/
1969
1970	if(pos > 0)
1971	{
1972	/printf("Cluster %d siwtching %d <-> %d\n", i, 0, pos);/
1973	buf = c[0];
1974	c[0] = c[pos];
1975	c[pos] = buf;
1976	}
1977	}
1978	for(j = 0; j < length; j++)
1979	c[j] = c[j] + 1;
1980	}
1981
1982	free(omitted);
1983	free(current);
1984	free(best);
1985	*ccc = clusterCount;
1986	return clusters;
1987	}
1988
1989	static void reduceBySequenceSimilarity(tree tr, rawdata rdta, analdef *adef)
1990	{
1991	int n = tr->mxtips + 1;
1992	int i, j;
1993	int omissionList = (int )malloc(n * sizeof(int));
1994	int undeterminedList = (int )malloc((rdta->sites + 1)* sizeof(int));
1995	int modelList = (int )malloc((rdta->sites + 1)* sizeof(int));
1996	int countNameDuplicates = 0;
1997	int countUndeterminedColumns = 0;
1998	int countOnlyGaps = 0;
1999	int modelCounter = 1;
2000	char buf[16], outName[1024];
2001	char undetermined_AA = 22;
2002	char undetermined_DNA = 15;
2003	char *tipI;
2004	qtList clusters = (qtList)NULL;
2005	FILE f, assoc;
2006	int numberOfClusters = 0;
2007	int nonTrivial = 0;
2008
2009	strcpy(outName, workdir);
2010	strcat(outName, "RAxML_reducedList.");
2011	strcat(outName, run_id);
2012
2013	if(processID == 0)
2014	f = fopen(infoFileName, "a");
2015	else
2016	f = (FILE *)NULL;
2017
2018
2019
2020	for(i = 1; i < n; i++)
2021	omissionList[i] = 0;
2022
2023	for(i = 0; i < rdta->sites + 1; i++)
2024	undeterminedList[i] = 0;
2025
2026	for(i = 1; i < n; i++)
2027	{
2028	for(j = i + 1; j < n; j++)
2029	if(strcmp(tr->nameList[i], tr->nameList[j]) == 0)
2030	{
2031	countNameDuplicates++;
2032	if(processID == 0)
2033	{
2034	printf("Sequence names of taxon %d and %d are identical, they are both called %s\n", i, j, tr->nameList[i]);
2035	fprintf(f, "Sequence names of taxon %d and %d are identical, they are both called %s\n", i, j, tr->nameList[i]);
2036	}
2037	}
2038	}
2039
2040	if(countNameDuplicates > 0)
2041	{
2042	if(processID == 0)
2043	{
2044	printf("ERROR: Found %d taxa that had equal names in the alignment, exiting...\n", countNameDuplicates);
2045	fprintf(f, "ERROR: Found %d taxa that had equal names in the alignment, exiting...\n", countNameDuplicates);
2046	fclose(f);
2047	}
2048	errorExit(-1);
2049	}
2050
2051	for(i = 1; i < n; i++)
2052	{
2053	j = 1;
2054
2055	while(j <= rdta->sites)
2056	{
2057	if(tr->dataVector[j] == DNA_DATA && rdta->y[i][j] != undetermined_DNA)
2058	break;
2059	if(tr->dataVector[j] == AA_DATA && rdta->y[i][j] != undetermined_AA)
2060	break;
2061	j++;
2062	}
2063
2064	if(j == (rdta->sites + 1))
2065	{
2066	if(processID == 0)
2067	{
2068	printf("ERROR: Sequence %s consists entirely of undetermined values which will be treated as missing data\n", tr->nameList[i]);
2069	fprintf(f, "ERROR: Sequence %s consists entirely of undetermined values which will be treated as missing data\n", tr->nameList[i]);
2070	}
2071	countOnlyGaps++;
2072	}
2073
2074	}
2075
2076	if(countOnlyGaps > 0)
2077	{
2078	if(processID == 0)
2079	{
2080	printf("ERROR: Found %d sequences that consist entirely of undetermined values, exiting...\n", countOnlyGaps);
2081	fprintf(f, "ERROR: Found %d sequences that consist entirely of undetermined values, exiting...\n", countOnlyGaps);
2082	fclose(f);
2083	}
2084	errorExit(-1);
2085	}
2086
2087	for(i = 0; i <= rdta->sites; i++)
2088	modelList[i] = -1;
2089
2090	for(i = 1; i <= rdta->sites; i++)
2091	{
2092	j = 1;
2093
2094	while(j < n)
2095	{
2096	if(tr->dataVector[i] == DNA_DATA && rdta->y[j][i] != undetermined_DNA)
2097	break;
2098	if(tr->dataVector[i] == AA_DATA && rdta->y[j][i] != undetermined_AA)
2099	break;
2100	j++;
2101	}
2102
2103	if(j == n)
2104	{
2105	undeterminedList[i] = 1;
2106	if(processID == 0)
2107	{
2108	printf("IMPORTANT WARNING: Alignment column %d contains only undetermined values which will be treated as missing data\n", i);
2109	fprintf(f, "IMPORTANT WARNING: Alignment column %d contains only undetermined values which will be treated as missing data\n", i);
2110	}
2111	countUndeterminedColumns++;
2112	}
2113	else
2114	{
2115	if(adef->useMultipleModel)
2116	{
2117	modelList[modelCounter] = tr->model[i];
2118	modelCounter++;
2119	}
2120	}
2121	}
2122
2123	switch(adef->similarityFilterMode)
2124	{
2125	case SMALL_DATA:
2126	{
2127	float **d;
2128	int n = tr->mxtips;
2129	int i, j;
2130	double t = gettime();
2131	float nDouble = 1.0 / (float)(rdta->sites);
2132	int sites = rdta->sites;
2133	char tipI, tipJ;
2134
2135
2136	d = (float *)malloc(sizeof(float ) * n);
2137	for(i = 0; i < n; i++)
2138	d[i] = (float )malloc(sizeof(float) n);
2139
2140	for(i = 0; i < n; i++)
2141	{
2142	d[i][i] = 1.0;
2143	tipI = &(rdta->y[i + 1][1]);
2144	for(j = i + 1; j < n; j++)
2145	{
2146	int k;
2147	int count = 0;
2148	tipJ = &(rdta->y[j + 1][1]);
2149	for(k = 0; k < sites; k++)
2150	if(tipJ[k] == tipI[k])
2151	count++;
2152
2153	d[i][j] = ((float)count * nDouble);
2154	d[j][i] = d[i][j];
2155	}
2156	}
2157
2158	printf("DistMat %f\n", gettime() - t);
2159
2160	t = gettime();
2161	clusters = clusterQT(d, n, (float)(adef->sequenceSimilarity), &numberOfClusters);
2162	printf("QT %f %d\n", gettime() - t, numberOfClusters);
2163	}
2164	break;
2165	case LARGE_DATA:
2166	{
2167	double t;
2168
2169	t = gettime();
2170	clusters = clusterQT_LARGE(tr->mxtips, (float)(adef->sequenceSimilarity), &numberOfClusters, rdta);
2171	printf("QT %f %d\n", gettime() - t, numberOfClusters);
2172	}
2173	break;
2174	default:
2175	assert(0);
2176	}
2177
2178	assoc = fopen(outName, "w");
2179
2180	for(i = 0; i < numberOfClusters; i++)
2181	{
2182	int length = clusters[i].count;
2183	int *c = clusters[i].entries;
2184	int j;
2185
2186	if(length > 1)
2187	{
2188	fprintf(assoc, "%s:%s", tr->nameList[c[0]], tr->nameList[c[1]]);
2189	for(j = 2; j < length; j++)
2190	fprintf(assoc, ",%s", tr->nameList[c[j]]);
2191	fprintf(assoc, "\n");
2192
2193	nonTrivial++;
2194	}
2195	}
2196
2197	fclose(assoc);
2198
2199
2200	if(nonTrivial > 0 \|\| countUndeterminedColumns > 0)
2201	{
2202	char noDupFile[2048];
2203	char noDupModels[2048];
2204
2205	if(nonTrivial > 0)
2206	{
2207	if(processID == 0)
2208	{
2209	printf("\n");
2210
2211	printf("Found %d non-trival clusters, reduction to %d sequences\n", nonTrivial, numberOfClusters);
2212
2213	fprintf(f, "\n");
2214
2215	fprintf(f, "Found %d non-trival clusters, reduction to %d sequences\n", nonTrivial, numberOfClusters);
2216	}
2217	}
2218
2219	if(countUndeterminedColumns > 0)
2220	{
2221	if(processID == 0)
2222	{
2223	printf("\n");
2224
2225	printf("IMPORTANT WARNING\n");
2226
2227	printf("Found %d %s that %s only undetermined values which will be treated as missing data.\n",
2228	countUndeterminedColumns, (countUndeterminedColumns == 1)?"column":"columns", (countUndeterminedColumns == 1)?"contains":"contain");
2229	printf("Normally these columns should be excluded from the analysis.\n\n");
2230
2231	fprintf(f, "\n");
2232
2233	fprintf(f, "IMPORTANT WARNING\n");
2234
2235	fprintf(f, "Found %d %s that %s only undetermined values which will be treated as missing data.\n",
2236	countUndeterminedColumns, (countUndeterminedColumns == 1)?"column":"columns", (countUndeterminedColumns == 1)?"contains":"contain");
2237	fprintf(f, "Normally these columns should be excluded from the analysis.\n\n");
2238	}
2239	}
2240
2241	sprintf(buf, "%f", adef->sequenceSimilarity);
2242
2243	strcpy(noDupFile, seq_file);
2244	strcat(noDupFile, ".reducedBy.");
2245	strcat(noDupFile, buf);
2246
2247
2248	strcpy(noDupModels, modelFileName);
2249	strcat(noDupModels, ".reducedBy.");
2250	strcat(noDupModels, buf);
2251
2252
2253	if(processID == 0)
2254	{
2255	if(adef->useMultipleModel && !filexists(noDupModels) && countUndeterminedColumns)
2256	{
2257	FILE *newFile = fopen(noDupModels, "w");
2258
2259	printf("\nJust in case you might need it, a mixed model file with \n");
2260	printf("model assignments for undetermined columns removed is printed to file %s\n",noDupModels);
2261
2262	fprintf(f, "\nJust in case you might need it, a mixed model file with \n");
2263	fprintf(f, "model assignments for undetermined columns removed is printed to file %s\n",noDupModels);
2264
2265
2266	for(i = 0; i < tr->NumberOfModels; i++)
2267	{
2268	boolean modelStillExists = FALSE;
2269
2270	for(j = 1; (j <= rdta->sites) && (!modelStillExists); j++)
2271	{
2272	if(modelList[j] == i)
2273	modelStillExists = TRUE;
2274	}
2275
2276	if(modelStillExists)
2277	{
2278	char *protModels[10] = {"DAYHOFF", "DCMUT", "JTT", "MTREV", "WAG", "RTREV", "CPREV", "VT", "BLOSUM62", "MTMAM"};
2279	int k = 1;
2280	int lower, upper;
2281	int parts = 0;
2282
2283	switch(tr->partitionData[i].dataType)
2284	{
2285	case AA_DATA:
2286	{
2287	char AAmodel[1024];
2288
2289	strcpy(AAmodel, protModels[tr->partitionData[i].protModels]);
2290	if(tr->partitionData[i].protFreqs)
2291	strcat(AAmodel, "F");
2292
2293	fprintf(newFile, "%s, ", AAmodel);
2294	}
2295	break;
2296	case DNA_DATA:
2297	fprintf(newFile, "DNA, ");
2298	break;
2299	default:
2300	assert(0);
2301	}
2302
2303	fprintf(newFile, "%s = ", tr->partitionData[i].partitionName);
2304
2305
2306	while(k <= rdta->sites)
2307	{
2308	if(modelList[k] == i)
2309	{
2310	lower = k;
2311	while((modelList[k + 1] == i) && (k <= rdta->sites))
2312	k++;
2313	upper = k;
2314
2315	if(lower == upper)
2316	{
2317	if(parts == 0)
2318	fprintf(newFile, "%d", lower);
2319	else
2320	fprintf(newFile, ",%d", lower);
2321	}
2322	else
2323	{
2324	if(parts == 0)
2325	fprintf(newFile, "%d-%d", lower, upper);
2326	else
2327	fprintf(newFile, ",%d-%d", lower, upper);
2328	}
2329	parts++;
2330	}
2331	k++;
2332	}
2333	fprintf(newFile, "\n");
2334	}
2335	}
2336	fclose(newFile);
2337	}
2338	else
2339	{
2340	if(adef->useMultipleModel)
2341	{
2342	printf("\n A mixed model file with model assignments for undetermined\n");
2343	printf("columns removed has already been printed to file %s\n",noDupModels);
2344
2345	fprintf(f, "\n A mixed model file with model assignments for undetermined\n");
2346	fprintf(f, "columns removed has already been printed to file %s\n",noDupModels);
2347	}
2348	}
2349
2350
2351	if(!filexists(noDupFile))
2352	{
2353	FILE *newFile;
2354
2355	printf("Just in case you might need it, an alignment file with \n");
2356	if(nonTrivial && !countUndeterminedColumns)
2357	printf("similar sequences removed is printed to file %s\n", noDupFile);
2358	if(!nonTrivial && countUndeterminedColumns)
2359	printf("undetermined columns removed is printed to file %s\n", noDupFile);
2360	if(nonTrivial && countUndeterminedColumns)
2361	printf("similar sequences and undetermined columns removed is printed to file %s\n", noDupFile);
2362
2363	fprintf(f, "Just in case you might need it, an alignment file with \n");
2364	if(nonTrivial && !countUndeterminedColumns)
2365	fprintf(f, "similar sequences removed is printed to file %s\n", noDupFile);
2366	if(!nonTrivial && countUndeterminedColumns)
2367	fprintf(f, "undetermined columns removed is printed to file %s\n", noDupFile);
2368	if(nonTrivial && countUndeterminedColumns)
2369	fprintf(f, "similar sequences and undetermined columns removed is printed to file %s\n", noDupFile);
2370
2371	newFile = fopen(noDupFile, "w");
2372
2373	fprintf(newFile, "%d %d\n", numberOfClusters, rdta->sites - countUndeterminedColumns);
2374
2375	for(i = 0; i < numberOfClusters; i++)
2376	{
2377
2378	fprintf(newFile, "%s ", tr->nameList[clusters[i].entries[0]]);
2379	tipI = &(rdta->y[clusters[i].entries[0]][1]);
2380
2381	for(j = 0; j < rdta->sites; j++)
2382	{
2383	if(undeterminedList[j + 1] == 0)
2384	{
2385	switch(tr->dataVector[j + 1])
2386	{
2387	case AA_DATA:
2388	fprintf(newFile, "%c", inverseMeaningPROT[tipI[j]]);
2389	break;
2390	case DNA_DATA:
2391	fprintf(newFile, "%c", inverseMeaningDNA[tipI[j]]);
2392	break;
2393	default:
2394	assert(0);
2395	}
2396	}
2397	}
2398
2399	fprintf(newFile, "\n");
2400	}
2401	fclose(newFile);
2402	}
2403	else
2404	{
2405	if(nonTrivial && !countUndeterminedColumns)
2406	printf("An alignment file with similar sequences removed has already\n");
2407	if(!nonTrivial && countUndeterminedColumns)
2408	printf("An alignment file with undetermined columns removed has already\n");
2409	if(nonTrivial && countUndeterminedColumns)
2410	printf("An alignment file with undetermined columns and similar sequences removed has already\n");
2411
2412	printf("been printed to file %s\n", noDupFile);
2413
2414	if(nonTrivial && !countUndeterminedColumns)
2415	fprintf(f, "An alignment file with similar sequences removed has already\n");
2416	if(!nonTrivial && countUndeterminedColumns)
2417	fprintf(f, "An alignment file with undetermined columns removed has already\n");
2418	if(nonTrivial && countUndeterminedColumns)
2419	fprintf(f, "An alignment file with undetermined columns and similar sequences removed has already\n");
2420
2421	fprintf(f, "been printed to file %s\n", noDupFile);
2422	}
2423	}
2424	}
2425
2426
2427	free(undeterminedList);
2428	free(omissionList);
2429	free(modelList);
2430	if(processID == 0)
2431	fclose(f);
2432	}
2433
2434
2435
2436
2437	static void generateBS(tree tr, analdef adef)
2438	{
2439	int i, j, k, w;
2440	int count;
2441	char outName[1024], buf[16];
2442	FILE *of;
2443
2444	assert(adef->boot != 0);
2445
2446	for(i = 0; i < adef->multipleRuns; i++)
2447	{
2448	makeboot(adef, tr);
2449
2450	count = 0;
2451	for(j = 0; j < tr->cdta->endsite; j++)
2452	count += tr->cdta->aliaswgt[j];
2453
2454	assert(count == tr->rdta->sites);
2455
2456	strcpy(outName, workdir);
2457	strcat(outName, seq_file);
2458	strcat(outName, ".BS");
2459	sprintf(buf, "%d", i);
2460	strcat(outName, buf);
2461	printf("Printing replicate %d to %s\n", i, outName);
2462
2463	of = fopen(outName, "w");
2464
2465	fprintf(of, "%d %d\n", tr->mxtips, count);
2466
2467	for(j = 1; j <= tr->mxtips; j++)
2468	{
2469	char *tip = tr->yVector[tr->nodep[j]->number];
2470	fprintf(of, "%s ", tr->nameList[j]);
2471
2472	for(k = 0; k < tr->cdta->endsite; k++)
2473	{
2474	switch(tr->dataVector[k])
2475	{
2476	case DNA_DATA:
2477	for(w = 0; w < tr->cdta->aliaswgt[k]; w++)
2478	fprintf(of, "%c", inverseMeaningDNA[tip[k]]);
2479	break;
2480	case AA_DATA:
2481	for(w = 0; w < tr->cdta->aliaswgt[k]; w++)
2482	fprintf(of, "%c", inverseMeaningPROT[tip[k]]);
2483	break;
2484	default:
2485	assert(0);
2486	}
2487	}
2488
2489	fprintf(of, "\n");
2490	}
2491	fclose(of);
2492	}
2493	}
2494
2495
2496
2497
2498
2499	static void splitMultiGene(tree tr, rawdata rdta)
2500	{
2501	int i, l;
2502	int n = rdta->sites + 1;
2503	int modelFilter = (int )malloc(sizeof(int) * n);
2504	int length, k;
2505	char *tip;
2506	FILE *outf;
2507	char outFileName[2048];
2508	char buf[16];
2509
2510	for(i = 0; i < tr->NumberOfModels; i++)
2511	{
2512	strcpy(outFileName, seq_file);
2513	sprintf(buf, "%d", i);
2514	strcat(outFileName, ".GENE.");
2515	strcat(outFileName, buf);
2516	outf = fopen(outFileName, "w");
2517	length = 0;
2518	for(k = 1; k < n; k++)
2519	{
2520	if(tr->model[k] == i)
2521	{
2522	modelFilter[k] = 1;
2523	length++;
2524	}
2525	else
2526	modelFilter[k] = -1;
2527	}
2528
2529	fprintf(outf, "%d %d\n", rdta->numsp, length);
2530
2531	for(l = 1; l <= rdta->numsp; l++)
2532	{
2533	fprintf(outf, "%s ", tr->nameList[l]);
2534
2535	tip = &(rdta->y[l][0]);
2536
2537	for(k = 1; k < n; k++)
2538	{
2539	if(modelFilter[k] == 1)
2540	{
2541	switch(tr->dataVector[k])
2542	{
2543	case AA_DATA:
2544	fprintf(outf, "%c", inverseMeaningPROT[tip[k]]);
2545	break;
2546	case DNA_DATA:
2547	fprintf(outf, "%c", inverseMeaningDNA[tip[k]]);
2548	break;
2549	default:
2550	assert(0);
2551	}
2552	}
2553	}
2554	fprintf(outf, "\n");
2555
2556	}
2557
2558	fclose(outf);
2559
2560	printf("Wrote individual gene/partition alignment to file %s\n", outFileName);
2561	}
2562
2563	free(modelFilter);
2564	printf("Wrote all %d individual gene/partition alignments\n", tr->NumberOfModels);
2565	printf("Exiting normally\n");
2566	}
2567
2568
2569
2570	void calculateModelOffsets(tree *tr)
2571	{
2572	int
2573	i,
2574	patterns,
2575	currentOffset = 0,
2576	dnaSpan,
2577	aaSpan;
2578
2579	switch(tr->rateHetModel)
2580	{
2581	case CAT:
2582	dnaSpan = 4;
2583	aaSpan = 20;
2584	break;
2585	case GAMMA:
2586	case GAMMA_I:
2587	dnaSpan = 16;
2588	aaSpan = 80;
2589	break;
2590	default:
2591	assert(0);
2592	}
2593
2594	tr->partitionData[0].modelOffset = currentOffset;
2595
2596	for(i = 1; i < tr->NumberOfModels; i++)
2597	{
2598	patterns = tr->partitionData[i - 1].upper - tr->partitionData[i - 1].lower;
2599	switch(tr->partitionData[i - 1].dataType)
2600	{
2601	case AA_DATA:
2602	currentOffset += aaSpan * patterns;
2603	break;
2604	case DNA_DATA:
2605	currentOffset += dnaSpan * patterns;
2606	break;
2607	default:
2608	assert(0);
2609	}
2610
2611	tr->partitionData[i].modelOffset = currentOffset;
2612	}
2613	}
2614
2615
2616	void allocNodex (tree tr, analdef adef)
2617	{
2618	nodeptr p;
2619	int i;
2620
2621	assert(tr->expArray == (int*)NULL);
2622	assert(tr->likelihoodArray == (double*)NULL);
2623	assert(tr->sumBuffer == (double *)NULL);
2624
2625	#ifdef _LOCAL_DATA
2626	tr->currentModel = adef->model;
2627	masterBarrier(THREAD_ALLOC_LIKELIHOOD, tr);
2628	#else
2629	{
2630	int span;
2631
2632	tr->expArray = (int )malloc(tr->cdta->endsite tr->mxtips * sizeof(int));
2633
2634	if(tr->mixedData)
2635	{
2636	tr->numberOfProteinPositions = 0;
2637	tr->numberOfNucleotidePositions = 0;
2638	for(i = 0; i < tr->cdta->endsite; i++)
2639	{
2640	switch(tr->dataVector[i])
2641	{
2642	case AA_DATA:
2643	tr->numberOfProteinPositions++;
2644	break;
2645	case DNA_DATA:
2646	tr->numberOfNucleotidePositions++;
2647	break;
2648	default:
2649	assert(0);
2650	}
2651	}
2652
2653	switch(adef->model)
2654	{
2655	case M_PROTCAT:
2656	case M_GTRCAT:
2657	span = tr->numberOfNucleotidePositions * 4 + tr->numberOfProteinPositions * 20;
2658	tr->likelihoodArray = (double )malloc(tr->mxtips span * sizeof(double));
2659	tr->sumBuffer = (double )malloc(span sizeof(double));
2660	break;
2661	case M_PROTGAMMA:
2662	case M_GTRGAMMA:
2663	span = tr->numberOfNucleotidePositions * 16 + tr->numberOfProteinPositions * 80;
2664	tr->likelihoodArray = (double )malloc(tr->mxtips span * sizeof(double));
2665	tr->sumBuffer = (double )malloc(span sizeof(double));
2666	break;
2667	default:
2668	assert(0);
2669	}
2670
2671	calculateModelOffsets(tr);
2672	/printf("DNA %d AA %d\n", tr->numberOfNucleotidePositions, tr->numberOfProteinPositions); /
2673	}
2674	else
2675	{
2676	switch(adef->model)
2677	{
2678	case M_PROTCAT:
2679	span = 20 * tr->cdta->endsite;
2680	tr->likelihoodArray = (double )malloc(span tr->mxtips * sizeof(double));
2681	tr->sumBuffer = (double )malloc(span sizeof(double));
2682	break;
2683	case M_PROTGAMMA:
2684	span = 80 * tr->cdta->endsite;
2685	tr->likelihoodArray = (double )malloc(span tr->mxtips * sizeof(double));
2686	tr->sumBuffer = (double )malloc(span sizeof(double));
2687	break;
2688	case M_GTRGAMMA:
2689	span = 16 * tr->cdta->endsite;
2690	tr->likelihoodArray = (double )malloc(span tr->mxtips * sizeof(double));
2691	tr->sumBuffer = (double )malloc(span sizeof(double));
2692	break;
2693	case M_GTRCAT:
2694	span = 4 * tr->cdta->endsite;
2695	tr->likelihoodArray = (double )malloc(span tr->mxtips * sizeof(double));
2696	tr->sumBuffer = (double )malloc(span sizeof(double));
2697	break;
2698	default:
2699	assert(0);
2700	}
2701	}
2702
2703	for(i = 0; i < tr->mxtips; i++)
2704	tr->xVector[i] = &(tr->likelihoodArray[i * span]);
2705	}
2706	#endif
2707
2708	for (i = tr->mxtips + 1; (i <= 2*(tr->mxtips) - 2); i++)
2709	{
2710	p = tr->nodep[i];
2711	p->x = 1;
2712	}
2713	}
2714
2715
2716	void freeNodex(tree *tr)
2717	{
2718	nodeptr p;
2719	int i;
2720
2721	#ifdef _LOCAL_DATA
2722	masterBarrier(THREAD_FREE_LIKELIHOOD, tr);
2723	#else
2724	free(tr->expArray);
2725	free(tr->likelihoodArray);
2726	free(tr->sumBuffer);
2727	tr->expArray = (int*)NULL;
2728	tr->likelihoodArray = (double*)NULL;
2729	tr->sumBuffer = (double*)NULL;
2730	#endif
2731
2732	for (i = tr->mxtips + 1; (i <= 2*(tr->mxtips) - 2); i++)
2733	{
2734	p = tr->nodep[i];
2735	while(!p->x)
2736	p = p->next;
2737	p->x = 0;
2738	p->next->x = 0;
2739	p->next->next->x = 0;
2740	}
2741	}
2742
2743	static void initAdef(analdef *adef)
2744	{
2745
2746	adef->bootstrapBranchLengths = FALSE;
2747	adef->model = M_GTRCAT;
2748	adef->max_rearrange = 21;
2749	adef->stepwidth = 5;
2750	adef->initial = adef->bestTrav = 10;
2751	adef->initialSet = FALSE;
2752	adef->restart = FALSE;
2753	adef->mode = BIG_RAPID_MODE;
2754	adef->categories = 25;
2755	adef->boot = 0;
2756	adef->rapidBoot = 0;
2757	adef->useWeightFile = FALSE;
2758	adef->checkpoints = 0;
2759	adef->startingTreeOnly = 0;
2760	adef->useMixedModel = 0;
2761	adef->multipleRuns = 1;
2762	adef->useMultipleModel = FALSE;
2763	adef->likelihoodEpsilon = 0.1;
2764	adef->constraint = FALSE;
2765	adef->grouping = FALSE;
2766	adef->randomStartingTree = FALSE;
2767	adef->categorizeGamma = FALSE;
2768	adef->parsimonySeed = 0;
2769	adef->proteinMatrix = JTT;
2770	adef->protEmpiricalFreqs = 0;
2771	adef->outgroup = FALSE;
2772	adef->useInvariant = FALSE;
2773	adef->sequenceSimilarity = 1.0;
2774	adef->permuteTreeoptimize = FALSE;
2775	adef->useInvariant = FALSE;
2776	adef->allInOne = FALSE;
2777	adef->multiBoot = 0;
2778	adef->likelihoodTest = FALSE;
2779	adef->reallyThoroughBoot = FALSE;
2780	adef->perGeneBranchLengths = FALSE;
2781	adef->treeLength = FALSE;
2782	adef->computePerSiteLLs = FALSE;
2783	adef->generateBS = FALSE;
2784	adef->bootStopOnly = 0;
2785	adef->bootStopping = FALSE;
2786	adef->gapyness = 0.0;
2787	adef->similarityFilterMode = 0;
2788	adef->bootstopCutoff = 0.0;
2789	adef->useExcludeFile = FALSE;
2790	adef->userProteinModel = FALSE;
2791	adef->externalAAMatrix = (double*)NULL;
2792	adef->rapidML_Addition = FALSE;
2793	adef->computeELW = FALSE;
2794	#ifdef _VINCENT
2795	adef->optimizeBSmodel = TRUE;
2796	#endif
2797	}
2798
2799
2800
2801
2802	static int modelExists(char model, analdef adef)
2803	{
2804	int i;
2805	char *protModels[10] = {"DAYHOFF", "DCMUT", "JTT", "MTREV", "WAG", "RTREV", "CPREV", "VT", "BLOSUM62", "MTMAM"};
2806	char thisModel[1024];
2807
2808
2809	/********* DNA ********************/
2810
2811	if(strcmp(model, "GTRGAMMAI\0") == 0)
2812	{
2813	adef->model = M_GTRGAMMA;
2814	adef->useInvariant = TRUE;
2815	return 1;
2816	}
2817
2818	if(strcmp(model, "GTRGAMMA\0") == 0)
2819	{
2820	adef->model = M_GTRGAMMA;
2821	return 1;
2822	}
2823
2824	if(strcmp(model, "GTRCAT\0") == 0)
2825	{
2826	adef->model = M_GTRCAT;
2827	return 1;
2828	}
2829
2830	if(strcmp(model, "GTRMIX\0") == 0)
2831	{
2832	adef->model = M_GTRCAT;
2833	adef->useMixedModel = 1;
2834	return 1;
2835	}
2836
2837	if(strcmp(model, "GTRMIXI\0") == 0)
2838	{
2839	adef->model = M_GTRCAT;
2840	adef->useMixedModel = 1;
2841	adef->useInvariant = TRUE;
2842	return 1;
2843	}
2844
2845
2846	if(strcmp(model, "GTRCAT_GAMMA\0") == 0)
2847	{
2848	adef->model = M_GTRCAT;
2849	adef->useMixedModel = 1;
2850	adef->categorizeGamma = TRUE;
2851	return 1;
2852	}
2853
2854	if(strcmp(model, "GTRCAT_GAMMAI\0") == 0)
2855	{
2856	adef->model = M_GTRCAT;
2857	adef->useMixedModel = 1;
2858	adef->categorizeGamma = TRUE;
2859	adef->useInvariant = TRUE;
2860	return 1;
2861	}
2862
2863	/************* AA GTR ******************/
2864
2865	if(strcmp(model, "PROTCATGTR\0") == 0)
2866	{
2867	adef->model = M_PROTCAT;
2868	adef->proteinMatrix = GTR;
2869	return 1;
2870	}
2871	if(strcmp(model, "PROTMIXGTR\0") == 0)
2872	{
2873	adef->model = M_PROTCAT;
2874	adef->proteinMatrix = GTR;
2875	adef->useMixedModel = 1;
2876	return 1;
2877	}
2878	if(strcmp(model, "PROTGAMMAGTR\0") == 0)
2879	{
2880	adef->model = M_PROTGAMMA;
2881	adef->proteinMatrix = GTR;
2882	return 1;
2883	}
2884
2885	if(strcmp(model, "PROTCAT_GAMMAGTR\0") == 0)
2886	{
2887	adef->model = M_PROTCAT;
2888	adef->proteinMatrix = GTR;
2889	adef->useMixedModel = 1;
2890	adef->categorizeGamma = TRUE;
2891	return 1;
2892	}
2893
2894	if(strcmp(model, "PROTCAT_GAMMAIGTR\0") == 0)
2895	{
2896	adef->model = M_PROTCAT;
2897	adef->proteinMatrix = GTR;
2898	adef->useMixedModel = 1;
2899	adef->categorizeGamma = TRUE;
2900	adef->useInvariant = TRUE;
2901	return 1;
2902	}
2903
2904	/**************** AA **********************/
2905
2906	for(i = 0; i < 10; i++)
2907	{
2908	/* check CAT */
2909
2910	strcpy(thisModel, "PROTCAT");
2911	strcat(thisModel, protModels[i]);
2912
2913	if(strcmp(model, thisModel) == 0)
2914	{
2915	adef->model = M_PROTCAT;
2916	adef->proteinMatrix = i;
2917	return 1;
2918	}
2919
2920	/* check CATF */
2921
2922	strcpy(thisModel, "PROTCAT");
2923	strcat(thisModel, protModels[i]);
2924	strcat(thisModel, "F");
2925
2926	if(strcmp(model, thisModel) == 0)
2927	{
2928	adef->model = M_PROTCAT;
2929	adef->proteinMatrix = i;
2930	adef->protEmpiricalFreqs = 1;
2931	return 1;
2932	}
2933
2934	/**************check MIX **********************/
2935
2936	strcpy(thisModel, "PROTMIX");
2937	strcat(thisModel, protModels[i]);
2938
2939	if(strcmp(model, thisModel) == 0)
2940	{
2941	adef->model = M_PROTCAT;
2942	adef->proteinMatrix = i;
2943	adef->useMixedModel = 1;
2944	return 1;
2945	}
2946
2947	/check MIXI /
2948
2949	strcpy(thisModel, "PROTMIXI");
2950	strcat(thisModel, protModels[i]);
2951
2952	if(strcmp(model, thisModel) == 0)
2953	{
2954	adef->model = M_PROTCAT;
2955	adef->proteinMatrix = i;
2956	adef->useMixedModel = 1;
2957	adef->useInvariant = TRUE;
2958	return 1;
2959	}
2960
2961
2962	/* check MIXmodelF */
2963
2964	strcpy(thisModel, "PROTMIX");
2965	strcat(thisModel, protModels[i]);
2966	strcat(thisModel, "F");
2967
2968	if(strcmp(model, thisModel) == 0)
2969	{
2970	adef->model = M_PROTCAT;
2971	adef->proteinMatrix = i;
2972	adef->useMixedModel = 1;
2973	adef->protEmpiricalFreqs = 1;
2974	return 1;
2975	}
2976
2977	/* check MIXImodelF */
2978
2979	strcpy(thisModel, "PROTMIXI");
2980	strcat(thisModel, protModels[i]);
2981	strcat(thisModel, "F");
2982
2983	if(strcmp(model, thisModel) == 0)
2984	{
2985	adef->model = M_PROTCAT;
2986	adef->proteinMatrix = i;
2987	adef->useMixedModel = 1;
2988	adef->protEmpiricalFreqs = 1;
2989	adef->useInvariant = TRUE;
2990	return 1;
2991	}
2992
2993	/**************check GAMMA **********************/
2994
2995	strcpy(thisModel, "PROTGAMMA");
2996	strcat(thisModel, protModels[i]);
2997
2998	if(strcmp(model, thisModel) == 0)
2999	{
3000	adef->model = M_PROTGAMMA;
3001	adef->proteinMatrix = i;
3002	return 1;
3003	}
3004
3005	/check GAMMAI/
3006
3007	strcpy(thisModel, "PROTGAMMAI");
3008	strcat(thisModel, protModels[i]);
3009
3010	if(strcmp(model, thisModel) == 0)
3011	{
3012	adef->model = M_PROTGAMMA;
3013	adef->proteinMatrix = i;
3014	adef->useInvariant = TRUE;
3015	return 1;
3016	}
3017
3018
3019	/* check GAMMAmodelF */
3020
3021	strcpy(thisModel, "PROTGAMMA");
3022	strcat(thisModel, protModels[i]);
3023	strcat(thisModel, "F");
3024
3025	if(strcmp(model, thisModel) == 0)
3026	{
3027	adef->model = M_PROTGAMMA;
3028	adef->proteinMatrix = i;
3029	adef->protEmpiricalFreqs = 1;
3030	return 1;
3031	}
3032
3033	/* check GAMMAImodelF */
3034
3035	strcpy(thisModel, "PROTGAMMAI");
3036	strcat(thisModel, protModels[i]);
3037	strcat(thisModel, "F");
3038
3039	if(strcmp(model, thisModel) == 0)
3040	{
3041	adef->model = M_PROTGAMMA;
3042	adef->proteinMatrix = i;
3043	adef->protEmpiricalFreqs = 1;
3044	adef->useInvariant = TRUE;
3045	return 1;
3046	}
3047
3048	/**************check CAT_GAMMA **********************/
3049
3050
3051	strcpy(thisModel, "PROTCAT_GAMMA");
3052	strcat(thisModel, protModels[i]);
3053
3054	if(strcmp(model, thisModel) == 0)
3055	{
3056	adef->model = M_PROTCAT;
3057	adef->useMixedModel = 1;
3058	adef->categorizeGamma = TRUE;
3059	adef->proteinMatrix = i;
3060	return 1;
3061	}
3062
3063	/* check CAT_GAMMAI */
3064
3065	strcpy(thisModel, "PROTCAT_GAMMAI");
3066	strcat(thisModel, protModels[i]);
3067
3068	if(strcmp(model, thisModel) == 0)
3069	{
3070	adef->model = M_PROTCAT;
3071	adef->useMixedModel = 1;
3072	adef->categorizeGamma = TRUE;
3073	adef->proteinMatrix = i;
3074	adef->useInvariant = TRUE;
3075	return 1;
3076	}
3077
3078	/* check CAT_GAMMAmodelF */
3079
3080	strcpy(thisModel, "PROTCAT_GAMMA");
3081	strcat(thisModel, protModels[i]);
3082	strcat(thisModel, "F");
3083
3084	if(strcmp(model, thisModel) == 0)
3085	{
3086	adef->model = M_PROTCAT;
3087	adef->useMixedModel = 1;
3088	adef->categorizeGamma = TRUE;
3089	adef->proteinMatrix = i;
3090	adef->protEmpiricalFreqs = 1;
3091	return 1;
3092	}
3093
3094	/check CAT_GAMMAImodelF /
3095
3096	strcpy(thisModel, "PROTCAT_GAMMAI");
3097	strcat(thisModel, protModels[i]);
3098	strcat(thisModel, "F");
3099
3100	if(strcmp(model, thisModel) == 0)
3101	{
3102	adef->model = M_PROTCAT;
3103	adef->useMixedModel = 1;
3104	adef->categorizeGamma = TRUE;
3105	adef->proteinMatrix = i;
3106	adef->protEmpiricalFreqs = 1;
3107	adef->useInvariant = TRUE;
3108	return 1;
3109	}
3110
3111
3112
3113	}
3114
3115	/*********************************************************************************/
3116
3117
3118
3119	return 0;
3120	}
3121
3122
3123
3124	static int mygetopt(int argc, char *argv, char opts, int optind, char *optarg)
3125	{
3126	static int sp = 1;
3127	register int c;
3128	register char *cp;
3129
3130	if(sp == 1)
3131	{
3132	if(optind >= argc \|\| argv[optind][0] != '-' \|\| argv[*optind][1] == '\0')
3133	return -1;
3134	}
3135	else
3136	{
3137	if(strcmp(argv[*optind], "--") == 0)
3138	{
3139	optind = optind + 1;
3140	return -1;
3141	}
3142	}
3143
3144	c = argv[*optind][sp];
3145	if(c == ':' \|\| (cp=strchr(opts, c)) == 0)
3146	{
3147	printf(": illegal option -- %c \n", c);
3148	if(argv[*optind][++sp] == '\0')
3149	{
3150	optind = optind + 1;
3151	sp = 1;
3152	}
3153	return('?');
3154	}
3155	if(*++cp == ':')
3156	{
3157	if(argv[*optind][sp+1] != '\0')
3158	{
3159	optarg = &argv[optind][sp+1];
3160	optind = optind + 1;
3161	}
3162	else
3163	{
3164	optind = optind + 1;
3165	if(*optind >= argc)
3166	{
3167	printf(": option requires an argument -- %c\n", c);
3168	sp = 1;
3169	return('?');
3170	}
3171	else
3172	{
3173	optarg = argv[optind];
3174	optind = optind + 1;
3175	}
3176	}
3177	sp = 1;
3178	}
3179	else
3180	{
3181	if(argv[*optind][++sp] == '\0')
3182	{
3183	sp = 1;
3184	optind = optind + 1;
3185	}
3186	*optarg = 0;
3187	}
3188	return(c);
3189	}
3190
3191	static void checkOutgroups(tree tr, analdef adef)
3192	{
3193	if(adef->outgroup)
3194	{
3195	boolean found;
3196	int i, j;
3197
3198	if(tr->numberOfOutgroups != 1 && adef->mode == MEHRING_ALGO)
3199	{
3200	printf("Error, you must specify exactly one sequence via \"-o\" to \n");
3201	printf("to run the sequence position determination algorithm\n");
3202	exit(-1);
3203	}
3204
3205	for(j = 0; j < tr->numberOfOutgroups; j++)
3206	{
3207	found = FALSE;
3208	for(i = 1; (i <= tr->mxtips) && !found; i++)
3209	{
3210	if(strcmp(tr->nameList[i], tr->outgroups[j]) == 0)
3211	{
3212	tr->outgroupNums[j] = i;
3213	found = TRUE;
3214	}
3215	}
3216	if(!found)
3217	{
3218	printf("Error, the outgroup name \"%s\" you specified can not be found in the alignment, exiting ....\n", tr->outgroups[j]);
3219	errorExit(-1);
3220	}
3221	}
3222	}
3223
3224	}
3225
3226	static void parseOutgroups(char outgr[2048], tree *tr)
3227	{
3228	int count = 1, i, k;
3229	char name[nmlngth];
3230
3231	i = 0;
3232	while(outgr[i] != '\0')
3233	{
3234	if(outgr[i] == ',')
3235	count++;
3236	i++;
3237	}
3238
3239	tr->numberOfOutgroups = count;
3240
3241	tr->outgroups = (char *)malloc(sizeof(char ) * count);
3242
3243	for(i = 0; i < tr->numberOfOutgroups; i++)
3244	tr->outgroups[i] = (char )malloc(sizeof(char) nmlngth);
3245
3246	tr->outgroupNums = (int )malloc(sizeof(int) count);
3247
3248	i = 0;
3249	k = 0;
3250	count = 0;
3251	while(outgr[i] != '\0')
3252	{
3253	if(outgr[i] == ',')
3254	{
3255	name[k] = '\0';
3256	strcpy(tr->outgroups[count], name);
3257	count++;
3258	k = 0;
3259	}
3260	else
3261	{
3262	name[k] = outgr[i];
3263	k++;
3264	}
3265	i++;
3266	}
3267
3268	name[k] = '\0';
3269	strcpy(tr->outgroups[count], name);
3270
3271	/*for(i = 0; i < tr->numberOfOutgroups; i++)
3272	printf("%d %s \n", i, tr->outgroups[i]);*/
3273
3274
3275	/printf("%s \n", name);/
3276	}
3277
3278
3279	/********************************* OUTGROUP STUFF END *******************************************************/
3280	static void printVersionInfo(void)
3281	{
3282	printf("\nThis is %s version %s released by Alexandros Stamatakis in %s\n\n", programName, programVersion, programDate);
3283	}
3284
3285	static void printMinusFUsage(void)
3286	{
3287	printf("\n");
3288	printf(" \"-f a\": rapid Bootstrap analysis and search for best-scoring ML tree in one program run\n");
3289	printf(" \"-f b\": draw bipartition information on a tree provided with \"-t\" based on multiple trees\n");
3290	printf(" (e.g. form a bootstrap) in a file specifed by \"-z\"\n");
3291	printf(" \"-f c\": check if the alignment can be properly read by RAxML\n");
3292	printf(" \"-f d\": new rapid hill-climbing \n");
3293	printf(" \"-f e\": optimize model+branch lengths for given input tree under GAMMA/GAMMAI only\n");
3294
3295	/*printf(" \"-f f\": optimize individual per-site evolutionary rates
3296	on a fixed input tree and compute sliding window tree lengths\n");*/
3297
3298	printf(" \"-f g\": compute per site log Likelihoods for one ore more trees passed via\n");
3299	printf(" \"-z\" and write them to a file that can be read by CONSEL\n");
3300	printf(" \"-f h\": compute log likelihood test (SH-test) between best tree passed via \"-t\"\n");
3301	printf(" and a bunch of other trees passed via \"-z\" \n");
3302	printf(" \"-f i\": perform a really thorough bootstrap, refinement of final BS tree under GAMMA and a\n");
3303	printf(" more exhaustive algorithm\n");
3304	printf(" \"-f j\": generate a bunch of bootstrapped alignment files from an original alignemnt file\n");
3305
3306	/* printf(" \"-f k\": \n"); */
3307
3308	printf(" \"-f m\": Compare bipartitions between two bunches of trees passed via \"-t\" and \"-z\" \n");
3309	printf(" respectively. This will return the Pearson correlation between all bipartitions found\n");
3310	printf(" in the two tree files. A file called RAxML_bipartitionFrequencies.outpuFileName\n");
3311	printf(" will be printed that contains the pair-wise bipartition frequencies of the two sets\n");
3312	printf(" \"-f n\": Compute the log likelihood score of all trees contained in a tree file provided by\n");
3313	printf(" \"-z\" under GAMMA or GAMMA+P-Invar\n");
3314	printf(" \"-f o\": old and slower rapid hill-climbing \n");
3315	printf(" \"-f p\": perform pure stepwise MP addition of new sequences to an incomplete starting tree\n");
3316
3317	/* printf(" \"-f r\": optimize individual per-site evolutionary rates on a fixed input tree\n");*/
3318
3319	printf(" \"-f s\": split up a multi-gene partitioned alignment into the respective subalignments \n");
3320	printf(" \"-f t\": do randomized tree searches on one fixed starting tree\n");
3321	printf(" \"-f w\": compute ELW test on a bunch of trees passed via \"-z\" \n");
3322	printf("\n");
3323	printf(" DEFAULT: new rapid hill climbing\n");
3324	printf("\n");
3325	}
3326
3327
3328	static void printREADME(void)
3329	{
3330	printVersionInfo();
3331	printf("\n");
3332	printf("Please also consult the RAxML-manual\n");
3333	printf("To report bugs send an email to Alexandros.Stamatakis@epfl.ch\n\n\n");
3334
3335	printf("raxmlHPC[-MPI\|-PTHREADS] -s sequenceFileName -n outputFileName -m substitutionModel\n");
3336	printf(" [-a weightFileName] [-b bootstrapRandomNumberSeed] [-c numberOfCategories]\n");
3337	printf(" [-d] [-e likelihoodEpsilon] [-E excludeFileName] [-f a\|b\|c\|d\|e\|g\|h\|i\|j\|m\|n\|o\|p\|s\|t\|w]\n");
3338	printf(" [-g groupingFileName] [-h] [-i initialRearrangementSetting] [-j] [-k] \n");
3339	printf(" [-l sequenceSimilarityThreshold] [-L sequenceSimilarityThreshold] [-M]\n");
3340	printf(" [-o outGroupName1[,outGroupName2[,...]]] [-p parsimonyRandomSeed] [-P proteinModel]\n");
3341	printf(" [-q multipleModelFileName] [-r binaryConstraintTree] [-t userStartingTree]\n");
3342	printf(" [-T numberOfThreads] [-u multiBootstrapSearches] [-v][-w workingDirectory]\n");
3343	printf(" [-x rapidBootstrapRandomNumberSeed][-y][-z multipleTreesFile] [-#\|-N numberOfRuns]\n");
3344	printf("\n");
3345	printf(" -a Specify a column weight file name to assign individual weights to each column of \n");
3346	printf(" the alignment. Those weights must be integers separated by any type and number \n");
3347	printf(" of whitespaces whithin a separate file, see file \"example_weights\" for an example.\n");
3348	printf("\n");
3349	printf(" -b Specify an integer number (random seed) and turn on bootstrapping\n");
3350	printf("\n");
3351	printf(" DEFAULT: OFF\n");
3352	printf("\n");
3353	printf(" -c Specify number of distinct rate catgories for RAxML when modelOfEvolution\n");
3354	printf(" is set to GTRCAT or GTRMIX\n");
3355	printf(" Individual per-site rates are categorized into numberOfCategories rate \n");
3356	printf(" categories to accelerate computations. \n");
3357	printf("\n");
3358	printf(" DEFAULT: 25\n");
3359	printf("\n");
3360	printf(" -d start ML optimization from random starting tree \n");
3361	printf("\n");
3362	printf(" DEFAULT: OFF\n");
3363	printf("\n");
3364	printf(" -e set model optimization precision in log likelihood units for final\n");
3365	printf(" optimization of tree topology under MIX/MIXI or GAMMA/GAMMAI\n");
3366	printf("\n");
3367	printf(" DEFAULT: 0.1 for models not using proportion of invariant sites estimate\n");
3368	printf(" 0.001 for models using proportion of invariant sites estimate\n");
3369	printf("\n");
3370	printf(" -E specify an exclude file name, that contains a specification of alignment positions you wish to exclude.\n");
3371	printf(" Format is similar to Nexus, the file shall contain entries like \"100-200 300-400\", to exclude a\n");
3372	printf(" single column write, e.g., \"100-100\", if you use a mixed model, an appropriatly adapted model file\n");
3373	printf(" will be written.\n");
3374	printf("\n");
3375	printf(" -f select algorithm:\n");
3376
3377	printMinusFUsage();
3378
3379	printf("\n");
3380	printf(" -g specify the file name of a multifurcating constraint tree\n");
3381	printf(" this tree does not need to be comprehensive, i.e. must not contain all taxa\n");
3382	printf("\n");
3383	printf(" -h Display this help message.\n");
3384	printf("\n");
3385	printf(" -i Initial rearrangement setting for the subsequent application of topological \n");
3386	printf(" changes phase\n");
3387	printf("\n");
3388	printf(" DEFAULT: determined by program\n");
3389	printf("\n");
3390	printf(" -j Specifies if checkpoints will be written by the program. If checkpoints \n");
3391	printf(" (intermediate tree topologies) shall be written by the program specify \"-j\"\n");
3392	printf("\n");
3393	printf(" DEFAULT: OFF\n");
3394	printf("\n");
3395	printf(" -k Specifies that bootstrapped trees should be printed with branch lengths.\n");
3396	printf(" The bootstraps will run a bit longer, because model parameters will be optimized\n");
3397	printf(" at the end of each run. Use with CATMIX/PROTMIX or GAMMA/GAMMAI.\n");
3398	printf("\n");
3399	printf(" DEFAULT: OFF\n");
3400	printf("\n");
3401	printf(" -l Specify a threshold for sequence similarity clustering. RAxML will then print out an alignment\n");
3402	printf(" to a file called sequenceFileName.reducedBy.threshold that only contains sequences <= the\n");
3403	printf(" specified thresold that must be between 0.0 and 1.0. RAxML uses the QT-clustering algorithm \n");
3404	printf(" to perform this task. In addition, a file called RAxML_reducedList.outputFileName will be written\n");
3405	printf(" that contains clustering information.\n");
3406	printf("\n");
3407	printf(" DEFAULT: OFF\n");
3408	printf("\n");
3409	printf(" -L Same functionality as \"-l\" above, but uses a less exhasutive and thus faster clustering algorithm\n");
3410	printf(" This is intended for very large datasets with more than 20,000-30,000 sequences\n");
3411	printf("\n");
3412	printf(" DEFAULT: OFF\n");
3413	printf("\n");
3414	printf(" -m Model of Nucleotide or Amino Acid Substitution: \n");
3415	printf("\n");
3416	printf(" NUCLEOTIDES:\n\n");
3417	printf(" \"-m GTRCAT\" : GTR + Optimization of substitution rates + Optimization of site-specific\n");
3418	printf(" evolutionary rates which are categorized into numberOfCategories distinct \n");
3419	printf(" rate categories for greater computational efficiency\n");
3420	printf(" if you do a multiple analysis with \"-#\" or \"-N\" but without bootstrapping the program\n");
3421	printf(" will use GTRMIX instead\n");
3422	printf(" \"-m GTRGAMMA\" : GTR + Optimization of substitution rates + GAMMA model of rate \n");
3423	printf(" heterogeneity (alpha parameter will be estimated)\n");
3424	printf(" \"-m GTRMIX\" : Inference of the tree under GTRCAT\n");
3425	printf(" and thereafter evaluation of the final tree topology under GTRGAMMA\n");
3426	printf(" \"-m GTRCAT_GAMMA\" : Inference of the tree with site-specific evolutionary rates.\n");
3427	printf(" However, here rates are categorized using the 4 discrete GAMMA rates.\n");
3428	printf(" Evaluation of the final tree topology under GTRGAMMA\n");
3429	printf(" \"-m GTRGAMMAI\" : Same as GTRGAMMA, but with estimate of proportion of invariable sites \n");
3430	printf(" \"-m GTRMIXI\" : Same as GTRMIX, but with estimate of proportion of invariable sites \n");
3431	printf(" \"-m GTRCAT_GAMMAI\" : Same as GTRCAT_GAMMA, but with estimate of proportion of invariable sites \n");
3432	printf("\n");
3433	printf(" AMINO ACIDS:\n\n");
3434	printf(" \"-m PROTCATmatrixName[F]\" : specified AA matrix + Optimization of substitution rates + Optimization of site-specific\n");
3435	printf(" evolutionary rates which are categorized into numberOfCategories distinct \n");
3436	printf(" rate categories for greater computational efficiency\n");
3437	printf(" if you do a multiple analysis with \"-#\" or \"-N\" but without bootstrapping the program\n");
3438	printf(" will use PROTMIX... instead\n");
3439	printf(" \"-m PROTGAMMAmatrixName[F]\" : specified AA matrix + Optimization of substitution rates + GAMMA model of rate \n");
3440	printf(" heterogeneity (alpha parameter will be estimated)\n");
3441	printf(" \"-m PROTMIXmatrixName[F]\" : Inference of the tree under specified AA matrix + CAT\n");
3442	printf(" and thereafter evaluation of the final tree topology under specified AA matrix + GAMMA\n");
3443	printf(" \"-m PROTCAT_GAMMAmatrixName[F]\" : Inference of the tree under specified AA matrix and site-specific evolutionary rates.\n");
3444	printf(" However, here rates are categorized using the 4 discrete GAMMA rates.\n");
3445	printf(" Evaluation of the final tree topology under specified AA matrix + GAMMA\n");
3446	printf(" \"-m PROTGAMMAImatrixName[F]\" : Same as PROTGAMMAmatrixName[F], but with estimate of proportion of invariable sites \n");
3447	printf(" \"-m PROTMIXImatrixName[F]\" : Same as PROTMIXmatrixName[F], but with estimate of proportion of invariable sites \n");
3448	printf(" \"-m PROTCAT_GAMMAImatrixName[F]\" : Same as PROTCAT_GAMMAmatrixName[F], but with estimate of proportion of invariable sites \n");
3449	printf("\n");
3450	printf(" Available AA substitution models: DAYHOFF, DCMUT, JTT, MTREV, WAG, RTREV, CPREV, VT, BLOSUM62, MTMAM, GTR\n");
3451	printf(" With the optional \"F\" appendix you can specify if you want to use empirical base frequencies\n");
3452	printf(" Please not that for mixed models you can in addition specify the per-gene AA model in\n");
3453	printf(" the mixed model file (see manual for details)\n");
3454	printf("\n");
3455	printf(" -M Switch on estimation of individual per-partition branch lengths. Only has effect when used in combination with \"-q\"\n");
3456	printf(" Branch lengths for individual partitions will be printed to separate files\n");
3457	printf(" A weighted average of the branch lengths is computed by using the respective partition lengths\n");
3458	printf("\n"),
3459	printf(" DEFAULT: OFF\n");
3460	printf("\n");
3461	printf(" -n Specifies the name of the output file.\n");
3462	printf("\n");
3463	printf(" -o Specify the name of a single outgrpoup or a comma-separated list of outgroups, eg \"-o Rat\" \n");
3464	printf(" or \"-o Rat,Mouse\", in case that multiple outgroups are not monophyletic the first name \n");
3465	printf(" in the list will be selected as outgroup, don't leave spaces between taxon names!\n");
3466	printf("\n");
3467	printf(" -q Specify the file name which contains the assignment of models to alignment\n");
3468	printf(" partitions for multiple models of substitution. For the syntax of this file\n");
3469	printf(" please consult the manual.\n");
3470	printf("\n");
3471	printf(" -p Specify a random number seed for the parsimony inferences. This allows you to reproduce your results\n");
3472	printf(" and will help me debug the program. This option HAS NO EFFECT in the parallel MPI version\n");
3473	printf("\n");
3474	printf(" -P Specify the file name of a user-defined AA (Protein) substitution model. This file must contain\n");
3475	printf(" 420 entries, the first 400 being the AA substitution rates (this must be a symmetric matrix) and the\n");
3476	printf(" last 20 are the empirical base frequencies\n");
3477	printf("\n");
3478	printf(" -r Specify the file name of a binary constraint tree.\n");
3479	printf(" this tree does not need to be comprehensive, i.e. must not contain all taxa\n");
3480	printf("\n");
3481	printf(" -s Specify the name of the alignment data file in PHYLIP format\n");
3482	printf("\n");
3483	printf(" -t Specify a user starting tree file name in Newick format\n");
3484	printf("\n");
3485	printf(" -T PTHREADS VERSION ONLY! Specify the number of threads you want to run.\n");
3486	printf(" Make sure to set \"-T\" to at most the number of CPUs you have on your machine,\n");
3487	printf(" otherwise, there will be a huge performance decrease!\n");
3488	printf("\n");
3489	printf(" -u Specify the number of multiple BS searches per replicate\n");
3490	printf(" to obtain better ML trees for each replicate\n");
3491	printf("\n");
3492	printf(" DEFAULT: One ML search per BS replicate\n");
3493	printf("\n");
3494	printf(" -v Display version information\n");
3495	printf("\n");
3496	printf(" -w Name of the working directory where RAxML will write its output files\n");
3497	printf("\n");
3498	printf(" DEFAULT: current directory\n");
3499	printf("\n");
3500	printf(" -x Specify an integer number (random seed) and turn on rapid bootstrapping\n");
3501	printf("\n");
3502	printf(" -y If you want to only compute a parsimony starting tree with RAxML specify \"-y\",\n");
3503	printf(" the program will exit after computation of the starting tree\n");
3504	printf("\n");
3505	printf(" DEFAULT: OFF\n");
3506	printf("\n");
3507	printf(" -z Specify the file name of a file containing multiple trees e.g. from a bootstrap\n");
3508	printf(" that shall be used to draw bipartition values onto a tree provided with \"-t\",\n");
3509	printf(" It can also be used to compute per site log likelihoods in combination with \"-f g\"\n");
3510	printf(" and to read a bunch of trees for a couple of other options (\"-f h\", \"-f m\", \"-f n\").\n");
3511	printf("\n");
3512	printf(" -#\|-N Specify the number of alternative runs on distinct starting trees\n");
3513	printf(" In combination with the \"-b\" option, this will invoke a multiple boostrap analysis\n");
3514	printf(" Note that \"-N\" has been added as an alternative since \"-#\" sometimes caused problems\n");
3515	printf(" with certain MPI job submission systems, since \"-#\" is often used to start comments\n");
3516	printf("\n");
3517	printf(" DEFAULT: 1 single analysis\n");
3518	printf("\n\n\n\n");
3519
3520	}
3521
3522
3523	static void get_args(int argc, char argv[], analdef adef, tree *tr)
3524	{
3525	int optind = 1;
3526	int c;
3527	boolean bad_opt=FALSE;
3528	char aut[256];
3529	char buf[2048];
3530	char *optarg;
3531	char model[2048] = "";
3532	char outgroups[2048] = "";
3533	char modelChar;
3534	double likelihoodEpsilon, sequenceSimilarity;
3535	int nameSet = 0,
3536	alignmentSet = 0,
3537	multipleRuns = 0,
3538	constraintSet = 0,
3539	treeSet = 0,
3540	groupSet = 0,
3541	modelSet = 0,
3542	treesSet = 0,
3543	multipleBoots = 0;
3544	long parsimonySeed = 0;
3545	run_id[0] = 0;
3546	workdir[0] = 0;
3547	seq_file[0] = 0;
3548	tree_file[0] = 0;
3549	model[0] = 0;
3550	weightFileName[0] = 0;
3551	modelFileName[0] = 0;
3552
3553	/********* tr inits ************/
3554
3555	#ifdef _USE_PTHREADS
3556	NumberOfThreads = 0;
3557	#endif
3558
3559	tr->doCutoff = TRUE;
3560
3561	/******* tr inits end***********/
3562
3563	#ifdef _VINCENT
3564	while(!bad_opt &&
3565	((c = mygetopt(argc,argv,"T:E:N:u:l:x:X:z:g:r:e:a:b:c:f:i:m:t:w:s:n:o:L:B:q:#:p:vdyjhkM", &optind, &optarg))!=-1))
3566	#else
3567	while(!bad_opt &&
3568	((c = mygetopt(argc,argv,"T:E:N:u:l:x:z:g:r:e:a:b:c:f:i:m:t:w:s:n:o:L:B:P:q:#:p:vdyjhkM", &optind, &optarg))!=-1))
3569	#endif
3570	{
3571	switch(c)
3572	{
3573	case 'T':
3574	#ifdef _USE_PTHREADS
3575	sscanf(optarg,"%d", &NumberOfThreads);
3576	#else
3577	if(processID == 0)
3578	{
3579	printf("Option -T does not have any effect with the sequential or parallel MPI version.\n");
3580	printf("It is used to specify the number of threads for the Pthreads-based parallelization\n");
3581	}
3582	#endif
3583	break;
3584	case 'P':
3585	strcpy(proteinModelFileName, optarg);
3586	adef->userProteinModel = TRUE;
3587	parseProteinModel(adef);
3588	break;
3589	case 'E':
3590	strcpy(excludeFileName, optarg);
3591	adef->useExcludeFile = TRUE;
3592	break;
3593	case 'M':
3594	adef->perGeneBranchLengths = TRUE;
3595	break;
3596	case 'u':
3597	sscanf(optarg,"%d", &multipleBoots);
3598	adef->multiBoot = multipleBoots;
3599	break;
3600	case 'o':
3601	strcpy(outgroups, optarg);
3602	parseOutgroups(outgroups, tr);
3603	adef->outgroup = TRUE;
3604	break;
3605	case 'k':
3606	adef->bootstrapBranchLengths = TRUE;
3607	break;
3608	case 'z':
3609	strcpy(bootStrapFile, optarg);
3610	treesSet = 1;
3611	break;
3612	case 'd':
3613	adef->randomStartingTree = TRUE;
3614	break;
3615	case 'g':
3616	strcpy(tree_file, optarg);
3617	adef->grouping = TRUE;
3618	adef->restart = TRUE;
3619	groupSet = 1;
3620	break;
3621	case 'r':
3622	strcpy(tree_file, optarg);
3623	adef->restart = TRUE;
3624	adef->constraint = TRUE;
3625	constraintSet = 1;
3626	break;
3627	case 'e':
3628	sscanf(optarg,"%lf", &likelihoodEpsilon);
3629	adef->likelihoodEpsilon = likelihoodEpsilon;
3630	break;
3631	case 'q':
3632	strcpy(modelFileName,optarg);
3633	adef->useMultipleModel = TRUE;
3634	break;
3635	case 'p':
3636	sscanf(optarg,"%ld", &parsimonySeed);
3637	adef->parsimonySeed = parsimonySeed;
3638	break;
3639	case 'N':
3640	case '#':
3641	/* TODO include auto in readme */
3642	if(sscanf(optarg,"%d", &multipleRuns) > 0)
3643	{
3644	adef->multipleRuns = multipleRuns;
3645	}
3646	else
3647	{
3648	if((sscanf(optarg,"%s", aut) > 0) &&
3649	(
3650	(strcmp(aut, "auto") == 0) \|\|
3651	(strcmp(aut, "Auto") == 0) \|\|
3652	(strcmp(aut, "AUTO") == 0) \|\|
3653	(strcmp(aut, "automatic") == 0) \|\|
3654	(strcmp(aut, "Automatic") == 0) \|\|
3655	(strcmp(aut, "AUTOMATIC") == 0)
3656	)
3657	)
3658	{
3659	adef->bootStopping = TRUE;
3660	adef->multipleRuns = 1000;
3661	}
3662	else
3663	{
3664	if(processID == 0)
3665	{
3666	printf("Use -# or -N option either with an integer, e.g., -# 100 or with -# auto\n");
3667	printf("or -N 100 or -N auto respectively, note that auto will not work for the\n");
3668	printf("MPI-based parallel version\n");
3669	}
3670	errorExit(0);
3671	}
3672	}
3673	break;
3674	case 'v':
3675	printVersionInfo();
3676	errorExit(0);
3677	case 'y':
3678	adef->startingTreeOnly = 1;
3679	break;
3680	case 'h':
3681	printREADME();
3682	errorExit(0);
3683	case 'j':
3684	adef->checkpoints = 1;
3685	break;
3686	case 'a':
3687	strcpy(weightFileName,optarg);
3688	adef->useWeightFile = TRUE;
3689	break;
3690	case 'b':
3691	sscanf(optarg,"%ld", &adef->boot);
3692	break;
3693	case 'x':
3694	sscanf(optarg,"%ld", &adef->rapidBoot);
3695	#ifdef _VINCENT
3696	adef->optimizeBSmodel = FALSE;
3697	#endif
3698	break;
3699	#ifdef _VINCENT
3700	case 'X':
3701	sscanf(optarg,"%ld", &adef->rapidBoot);
3702	adef->optimizeBSmodel = TRUE;
3703	break;
3704	#endif
3705	case 'c':
3706	sscanf(optarg, "%d", &adef->categories);
3707	break;
3708	case 'l':
3709	sscanf(optarg,"%lf", &sequenceSimilarity);
3710	adef->sequenceSimilarity = sequenceSimilarity;
3711	adef->mode = SEQUENCE_SIMILARITY_FILTER;
3712	adef->similarityFilterMode = SMALL_DATA;
3713	break;
3714	case 'L':
3715	sscanf(optarg,"%lf", &sequenceSimilarity);
3716	adef->sequenceSimilarity = sequenceSimilarity;
3717	adef->mode = SEQUENCE_SIMILARITY_FILTER;
3718	adef->similarityFilterMode = LARGE_DATA;
3719	break;
3720	case 'B':
3721	/* TODO include in readme */
3722	sscanf(optarg,"%lf", &(adef->bootstopCutoff));
3723	if(adef->bootstopCutoff <= 0.0)
3724	{
3725	printf("ERROR BootstopCutoff was set to %f, but must be greater than 0.0\n",
3726	adef->bootstopCutoff);
3727	errorExit(-1);
3728	}
3729	if(adef->bootstopCutoff == 0.5)
3730	{
3731	printf("\n\nWARNING: BootstopCutoff was set to %f, this is equivalent to default\n", adef->bootstopCutoff);
3732	printf("Bootstopping without the \"-B\" option. Are you sure that this is \n");
3733	printf("what you want to do?\n\n");
3734	}
3735	if(adef->bootstopCutoff > 0.5)
3736	{
3737	printf("ERROR BootstopCutoff was set to %f, but must be smaller or equal to 0.5\n",
3738	adef->bootstopCutoff);
3739	errorExit(-1);
3740	}
3741	break;
3742	case 'f':
3743	sscanf(optarg, "%c", &modelChar);
3744	switch(modelChar)
3745	{
3746	case 'a':
3747	adef->allInOne = TRUE;
3748	adef->mode = BIG_RAPID_MODE;
3749	tr->doCutoff = TRUE;
3750	break;
3751	case 'b':
3752	adef->mode = CALC_BIPARTITIONS;
3753	break;
3754	case 'c':
3755	adef->mode = CHECK_ALIGNMENT;
3756	break;
3757	case 'd':
3758	adef->mode = BIG_RAPID_MODE;
3759	tr->doCutoff = TRUE;
3760	break;
3761	case 'e':
3762	adef->mode = TREE_EVALUATION;
3763	break;
3764	case 'f':
3765	/* TODO include in readme */
3766	adef->mode = OPTIMIZE_RATES;
3767	adef->treeLength = TRUE;
3768	break;
3769	case 'g':
3770	adef->mode = OPTIMIZE_RATES;
3771	adef->computePerSiteLLs = TRUE;
3772	break;
3773	case 'h':
3774	adef->mode = TREE_EVALUATION;
3775	adef->likelihoodTest = TRUE;
3776	break;
3777	case 'i':
3778	adef->reallyThoroughBoot = TRUE;
3779	break;
3780	case 'j':
3781	adef->generateBS = TRUE;
3782	break;
3783	case 'k':
3784	/* TODO include in readme */
3785	adef->bootStopOnly = 1;
3786	break;
3787	case 'l':
3788	/* TODO include in readme */
3789	adef->bootStopOnly = 2;
3790	break;
3791	case 'm':
3792	adef->bootStopOnly = 3;
3793	break;
3794	case 'n':
3795	adef->bootStopOnly = 4;
3796	break;
3797	case 'o':
3798	adef->mode = BIG_RAPID_MODE;
3799	tr->doCutoff = FALSE;
3800	break;
3801	case 'p':
3802	adef->mode = PARSIMONY_ADDITION;
3803	break;
3804	case 'q':
3805	/* TODO include in README */
3806	adef->mode = MEHRING_ALGO;
3807	break;
3808	case 'r':
3809	adef->mode = OPTIMIZE_RATES;
3810	break;
3811	case 's':
3812	adef->mode = SPLIT_MULTI_GENE;
3813	break;
3814	case 't':
3815	adef->mode = BIG_RAPID_MODE;
3816	tr->doCutoff = TRUE;
3817	adef->permuteTreeoptimize = TRUE;
3818	break;
3819	case 'u':
3820	/* TODO readme */
3821	adef->mode = ARNDT_MODE;
3822	break;
3823	case 'v':
3824	/* TODO README */
3825	adef->rapidML_Addition = TRUE;
3826	break;
3827	case 'w':
3828	adef->computeELW = TRUE;
3829	break;
3830	default:
3831	{
3832	if(processID == 0)
3833	{
3834	printf("Error select one of the following algorithms via -f :\n");
3835	printMinusFUsage();
3836	}
3837	errorExit(-1);
3838	}
3839	}
3840	break;
3841	case 'i':
3842	sscanf(optarg, "%d", &adef->initial);
3843	adef->initialSet = TRUE;
3844	break;
3845	case 'n':
3846	strcpy(run_id,optarg);
3847	nameSet = 1;
3848	break;
3849	case 'w':
3850	strcpy(workdir,optarg);
3851	break;
3852	case 't':
3853	strcpy(tree_file, optarg);
3854	adef->restart = TRUE;
3855	treeSet = 1;
3856	break;
3857	case 's':
3858	strcpy(seq_file, optarg);
3859	alignmentSet = 1;
3860	break;
3861	case 'm':
3862	strcpy(model,optarg);
3863	if(modelExists(model, adef) == 0)
3864	{
3865	if(processID == 0)
3866	{
3867	printf("Model %s does not exist\n\n", model);
3868	printf("For DNA data use: GTRCAT or GTRGAMMA or\n");
3869	printf(" GTRMIX or GTRMIXI or\n");
3870	printf(" GTRGAMMAI or GTRCAT_GAMMAI or\n");
3871	printf(" GTRCAT_GAMMA\n\n");
3872	printf("For AA data use: PROTCATmatrixName[F] or PROTGAMMAmatrixName[F] or\n");
3873	printf(" PROTMIXmatrixName[F] or PROTMIXImatrixName[F] or\n");
3874	printf(" PROTGAMMAImatrixName[F] or PROTCAT_GAMMAImatrixName[F] or\n");
3875	printf(" PROTCAT_GAMMAImatrixName[F]\n\n");
3876	printf("The AA substitution matrix can be one of the following: \n");
3877	printf("DAYHOFF, DCMUT, JTT, MTREV, WAG, RTREV, CPREV, VT, BLOSUM62, MTMAM, GTR\n\n");
3878	printf("With the optional \"F\" appendix you can specify if you want to use empirical base frequencies\n");
3879	printf("Please note that for mixed models you can in addition specify the per-gene model in\n");
3880	printf("the mixed model file (see manual for details)\n");
3881	}
3882	errorExit(-1);
3883	}
3884	else
3885	modelSet = 1;
3886	break;
3887	default:
3888	errorExit(-1);
3889	}
3890	}
3891
3892	#ifdef _USE_PTHREADS
3893	if(NumberOfThreads < 2)
3894	{
3895	printf("\nThe number of threads is currently set to %d\n", NumberOfThreads);
3896	printf("Specify the number of threads to run via -T numberOfThreads\n");
3897	printf("NumberOfThreads must be set to an integer value greater than 1\n\n");
3898	errorExit(-1);
3899	}
3900	#endif
3901
3902
3903	if(adef->computeELW)
3904	{
3905	if(processID == 0)
3906	{
3907	if(adef->boot == 0)
3908	{
3909	printf("Error, you must specify a bootstrap seed via \"-b\" to compute ELW statistics\n");
3910	errorExit(-1);
3911	}
3912
3913	if(adef->multipleRuns < 2)
3914	{
3915	printf("Error, you must specify the number of BS replicates via \"-#\" or \"-N\" to compute ELW statistics\n");
3916	printf("it should be larger than 1, recommended setting is 100\n");
3917	errorExit(-1);
3918	}
3919
3920	if(!treesSet)
3921	{
3922	printf("Error, you must specify an input file containing several candidate trees\n");
3923	printf("via \"-z\" to compute ELW statistics.\n");
3924	errorExit(-1);
3925	}
3926
3927	if(!(adef->model == M_PROTGAMMA \|\| adef->model == M_GTRGAMMA))
3928	{
3929	printf("Error ELW test can only be conducted undetr GAMMA or GAMMA+P-Invar models\n");
3930	errorExit(-1);
3931	}
3932	}
3933	}
3934
3935	if(adef->mode == MEHRING_ALGO && !(adef->restart && adef->outgroup))
3936	{
3937	if(processID == 0)
3938	{
3939	printf("\nTo use the sequence position determination algorithm you have to specify a starting tree with \"-t\" \n");
3940	printf("and a taxon to be re-inserted with \"-o\" \n");
3941	errorExit(-1);
3942	}
3943	}
3944
3945
3946
3947	if(((!adef->boot) && (!adef->rapidBoot)) && adef->bootStopping)
3948	{
3949	if(processID == 0)
3950	{
3951	printf("Can't use automatic bootstopping without actually doing a Bootstrap\n");
3952	printf("Specify either -x randomNumberSeed (rapid) or -b randomNumberSeed (standard)\n");
3953	errorExit(-1);
3954	}
3955	}
3956
3957	if(adef->boot && adef->rapidBoot)
3958	{
3959	if(processID == 0)
3960	{
3961	printf("Can't use standard and rapid BOOTSTRAP simultaneously\n");
3962	errorExit(-1);
3963	}
3964	}
3965
3966	if(adef->rapidBoot && !(adef->mode == MEHRING_ALGO))
3967	{
3968	if(processID == 0 && (adef->restart \|\| treesSet))
3969	{
3970	printf("Error, starting tree(s) will be ignored by rapid Bootstrapping\n");
3971	errorExit(-1);
3972	}
3973
3974	if(processID == 0 && (groupSet \|\| constraintSet))
3975	{
3976	printf("Error, constraint tree will be ignored by rapid Bootstrapping\n");
3977	errorExit(-1);
3978	}
3979	}
3980
3981	if(adef->allInOne && (adef->rapidBoot == 0))
3982	{
3983	if(processID == 0)
3984	{
3985	printf("Error, to carry out an ML search after a rapid BS inference you must specify a random number seed with -x\n");
3986	errorExit(-1);
3987	}
3988	}
3989
3990	if(adef->mode == SEQUENCE_SIMILARITY_FILTER)
3991	{
3992	if(processID == 0)
3993	{
3994	if(adef->sequenceSimilarity <= 0.0 \|\| adef->sequenceSimilarity >= 1.0)
3995	{
3996	printf("\n ERROR: sequence similarity must be > 0.0 and < 1.0, exiting ...\n");
3997	errorExit(-1);
3998	}
3999	}
4000	}
4001
4002	if(adef->mode == OPTIMIZE_RATES)
4003	{
4004	if(adef->treeLength && !(adef->model == M_GTRGAMMA \|\| adef->model == M_PROTGAMMA))
4005	{
4006	if(processID == 0)
4007	printf("\n ERROR: Tree-Length-based sliding window approach only allowed under GAMMA model of rate heterogeneity!\n");
4008	errorExit(-1);
4009	}
4010
4011	if(adef->computePerSiteLLs)
4012	{
4013	if(!(adef->model == M_GTRGAMMA \|\| adef->model == M_PROTGAMMA))
4014	{
4015	if(processID == 0)
4016	printf("\n ERROR: Computation of per-site log LHs is only allowed under GAMMA model of rate heterogeneity!\n");
4017	errorExit(-1);
4018	}
4019	if(!treesSet)
4020	{
4021	if(processID == 0)
4022	printf("\n ERROR: For Computation of per-site log LHs you need to specify several input trees with \"-z\"\n");
4023	errorExit(-1);
4024	}
4025	}
4026
4027	if(!adef->restart)
4028	{
4029	if(!adef->computePerSiteLLs && processID == 0)
4030	{
4031	if(!adef->treeLength)
4032	printf("\n You need to specify an input tree with \"-t\" to optimize rates using \"-f r\"\n");
4033	else
4034	printf("\n You need to specify an input tree with \"-t\" to optimize rates and compute the sliding window tree length using \"-f f\"\n"); errorExit(-1);
4035	}
4036	}
4037	}
4038
4039	if(adef->bootstrapBranchLengths && (adef->model == M_GTRCAT \|\| adef->model == M_PROTCAT) && (!adef->useMixedModel))
4040	{
4041	if(processID == 0)
4042	{
4043	printf("\nWARNING: you want to print out the branch lengths of your bootstrapped trees\n");
4044	printf("WARNING: However you have currently chosen one of the CAT models where the branch lengths\n");
4045	printf("WARNING: are essentially meaningless, you should better use CATMIX/PROTMIX instead\n");
4046	}
4047	}
4048
4049	if(adef->mode == SPLIT_MULTI_GENE && (!adef->useMultipleModel))
4050	{
4051	if(processID == 0)
4052	{
4053	printf("\n Error, you are trying to split a multi-gene alignment into individual genes with the \"-f s\" option\n");
4054	printf("Without specifying a multiple model file with \"-q modelFileName\" \n");
4055	}
4056	errorExit(-1);
4057	}
4058
4059	if(adef->mode == CALC_BIPARTITIONS && !treesSet)
4060	{
4061	if(processID == 0)
4062	printf("\n Error, in bipartition computation mode you must specify a file containing multiple trees with the \"-z\" option\n");
4063	errorExit(-1);
4064	}
4065
4066	if(adef->mode == CALC_BIPARTITIONS && !adef->restart)
4067	{
4068	if(processID == 0)
4069	printf("\n Error, in bipartition computation mode you must specify a tree on which bipartition information will be drawn with the \"-t\" option\n");
4070	errorExit(-1);
4071	}
4072
4073	if(!modelSet)
4074	{
4075	if(processID == 0)
4076	printf("\n Error, you must specify a model of substitution with the \"-m\" option\n");
4077	errorExit(-1);
4078	}
4079
4080
4081
4082	if(adef->useMultipleModel && (adef->model == M_PROTGAMMA \|\| adef->model == M_PROTCAT) && (adef->proteinMatrix == GTR))
4083	{
4084	if(processID == 0)
4085	printf("\n Error GTR model of AA substiution in combination with mixed models is currently not implemented\n");
4086	errorExit(-1);
4087	}
4088
4089
4090
4091	if(!adef->restart && adef->mode == PARSIMONY_ADDITION)
4092	{
4093	if(processID == 0)
4094	{
4095	printf("\n You need to specify an incomplete binary input tree with \"-t\" to execute \n");
4096	printf(" RAxML MP stepwise addition with \"-f p\"\n");
4097	}
4098	errorExit(-1);
4099	}
4100
4101	if(adef->restart && adef->randomStartingTree)
4102	{
4103	if(processID == 0)
4104	{
4105	if(adef->constraint)
4106	{
4107	printf("\n Error you specified a binary constraint tree with -r AND the computation\n");
4108	printf("of a random starting tree with -d for the same run\n");
4109	}
4110	else
4111	{
4112	if(adef->grouping)
4113	{
4114	printf("\n Error you specified a multifurcating constraint tree with -g AND the computation\n");
4115	printf("of a random starting tree with -d for the same run\n");
4116	}
4117	else
4118	{
4119	printf("\n Error you specified a starting tree with -t AND the computation\n");
4120	printf("of a random starting tree with -d for the same run\n");
4121	}
4122	}
4123	}
4124	errorExit(-1);
4125	}
4126
4127	if(treeSet && constraintSet)
4128	{
4129	if(processID == 0)
4130	printf("\n Error you specified a binary constraint tree AND a starting tree for the same run\n");
4131	errorExit(-1);
4132	}
4133
4134
4135	if(treeSet && groupSet)
4136	{
4137	if(processID == 0)
4138	printf("\n Error you specified a multifurcating constraint tree AND a starting tree for the same run\n");
4139	errorExit(-1);
4140	}
4141
4142
4143	if(groupSet && constraintSet)
4144	{
4145	if(processID == 0)
4146	printf("\n Error you specified a bifurcating constraint tree AND a multifurcating constraint tree for the same run\n");
4147	errorExit(-1);
4148	}
4149
4150	if(adef->restart && adef->startingTreeOnly)
4151	{
4152	if(processID == 0)
4153	{
4154	printf("\n Error conflicting options: you want to compute only a parsimony starting tree with -y\n");
4155	printf(" while you actually specified a starting tree with -t %s\n", tree_file);
4156	}
4157	errorExit(-1);
4158	}
4159
4160	if(adef->mode == TREE_EVALUATION && (!adef->restart))
4161	{
4162	if(processID == 0)
4163	printf("\n Error: please specify a treefile for the tree you want to evaluate with -t\n");
4164	errorExit(-1);
4165	}
4166
4167	#ifdef PARALLEL
4168
4169	if(adef->mode == SPLIT_MULTI_GENE)
4170	{
4171	if(processID == 0)
4172	printf("Multi gene alignment splitting (-f s) not implemented for the MPI-Version\n");
4173	errorExit(-1);
4174	}
4175
4176	if(adef->mode == TREE_EVALUATION)
4177	{
4178	if(processID == 0)
4179	printf("Tree Evaluation mode (-f e) noot implemented for the MPI-Version\n");
4180	errorExit(-1);
4181	}
4182
4183	if(adef->mode == CALC_BIPARTITIONS)
4184	{
4185	if(processID == 0)
4186	printf("Computation of bipartitions (-f b) not implemented for the MPI-Version\n");
4187	errorExit(-1);
4188	}
4189
4190	if(adef->multipleRuns == 1)
4191	{
4192	if(processID == 0)
4193	{
4194	printf("Error: you are running the parallel MPI program but only want to compute one tree\n");
4195	printf("For the MPI version you must specify a number of trees greater than 1 with the -# or -N option\n");
4196	}
4197	errorExit(-1);
4198	}
4199	#endif
4200
4201
4202
4203	if(adef->mode == TREE_EVALUATION && (adef->model == M_GTRCAT \|\| adef->model == M_PROTCAT))
4204	{
4205	if(processID == 0)
4206	{
4207	printf("\n Error: No tree evaluation with GTRCAT/PROTCAT possible\n");
4208	printf("the GTRCAT likelihood values are instable at present and should not\n");
4209	printf("be used to compare trees based on ML values\n");
4210	}
4211	errorExit(-1);
4212	}
4213
4214	if(!nameSet)
4215	{
4216	if(processID == 0)
4217	printf("\n Error: please specify a name for this run with -n\n");
4218	errorExit(-1);
4219	}
4220
4221	if(! alignmentSet)
4222	{
4223	if(processID == 0)
4224	printf("\n Error: please specify an alignment for this run with -s\n");
4225	errorExit(-1);
4226	}
4227
4228
4229	#ifdef WIN32
4230	if(workdir[0]==0 \|\| workdir[0] != '\\')
4231	{
4232	getcwd(buf,sizeof(buf));
4233	if( buf[strlen(buf)-1] != '\\') strcat(buf,"\\");
4234	strcat(buf,workdir);
4235	if( buf[strlen(buf)-1] != '\\') strcat(buf,"\\");
4236	strcpy(workdir,buf);
4237	}
4238	#else
4239	if(workdir[0]==0 \|\| workdir[0] != '/')
4240	{
4241	getcwd(buf,sizeof(buf));
4242	if( buf[strlen(buf)-1] != '/') strcat(buf,"/");
4243	strcat(buf,workdir);
4244	if( buf[strlen(buf)-1] != '/') strcat(buf,"/");
4245	strcpy(workdir,buf);
4246	}
4247	#endif
4248
4249	return;
4250	}
4251
4252
4253
4254
4255	void errorExit(int e)
4256	{
4257	#ifdef PARALLEL
4258	MPI_Status msgStatus;
4259	int i, dummy;
4260
4261	if(processID == 0)
4262	{
4263	for(i = 1; i < numOfWorkers; i++)
4264	MPI_Send(&dummy, 1, MPI_INT, i, FINALIZE, MPI_COMM_WORLD);
4265
4266	MPI_Finalize();
4267	exit(e);
4268	}
4269	else
4270	{
4271	MPI_Recv(&dummy, 1, MPI_INT, 0, FINALIZE, MPI_COMM_WORLD, &msgStatus);
4272	MPI_Finalize();
4273	exit(e);
4274	}
4275	#else
4276	exit(e);
4277	#endif
4278	}
4279
4280
4281
4282	static void makeFileNames(void)
4283	{
4284	int infoFileExists = 0;
4285	#ifdef PARALLEL
4286	MPI_Status msgStatus;
4287	#endif
4288
4289	strcpy(permFileName, workdir);
4290	strcpy(resultFileName, workdir);
4291	strcpy(logFileName, workdir);
4292	strcpy(checkpointFileName, workdir);
4293	strcpy(infoFileName, workdir);
4294	strcpy(randomFileName, workdir);
4295	strcpy(bootstrapFileName, workdir);
4296	strcpy(bipartitionsFileName, workdir);
4297	strcpy(ratesFileName, workdir);
4298	strcpy(lengthFileName, workdir);
4299	strcpy(lengthFileNameModel, workdir);
4300	strcpy( perSiteLLsFileName, workdir);
4301
4302	strcat(permFileName, "RAxML_parsimonyTree.");
4303	strcat(resultFileName, "RAxML_result.");
4304	strcat(logFileName, "RAxML_log.");
4305	strcat(checkpointFileName, "RAxML_checkpoint.");
4306	strcat(infoFileName, "RAxML_info.");
4307	strcat(randomFileName, "RAxML_randomTree.");
4308	strcat(bootstrapFileName, "RAxML_bootstrap.");
4309	strcat(bipartitionsFileName, "RAxML_bipartitions.");
4310	strcat(ratesFileName, "RAxML_perSiteRates.");
4311	strcat(lengthFileName, "RAxML_treeLength.");
4312	strcat(lengthFileNameModel, "RAxML_treeLengthModel.");
4313	strcat( perSiteLLsFileName, "RAxML_perSiteLLs.");
4314
4315	strcat(permFileName, run_id);
4316	strcat(resultFileName, run_id);
4317	strcat(logFileName, run_id);
4318	strcat(checkpointFileName, run_id);
4319	strcat(infoFileName, run_id);
4320	strcat(randomFileName, run_id);
4321	strcat(bootstrapFileName, run_id);
4322	strcat(bipartitionsFileName, run_id);
4323	strcat(ratesFileName, run_id);
4324	strcat(lengthFileName, run_id);
4325	strcat(lengthFileNameModel, run_id);
4326	strcat(perSiteLLsFileName, run_id);
4327
4328	if(processID == 0)
4329	{
4330	infoFileExists = filexists(infoFileName);
4331
4332	#ifdef PARALLEL
4333	{
4334	int i;
4335
4336	for(i = 1; i < numOfWorkers; i++)
4337	MPI_Send(&infoFileExists, 1, MPI_INT, i, FINALIZE, MPI_COMM_WORLD);
4338	}
4339	#endif
4340
4341	if(infoFileExists)
4342	{
4343	printf("RAxML output files with the run ID <%s> already exist \n", run_id);
4344	printf("in directory %s ...... exiting\n", workdir);
4345	#ifdef PARALLEL
4346	MPI_Finalize();
4347	exit(-1);
4348	#else
4349	exit(-1);
4350	#endif
4351	}
4352	}
4353	#ifdef PARALLEL
4354	else
4355	{
4356	MPI_Recv(&infoFileExists, 1, MPI_INT, 0, FINALIZE, MPI_COMM_WORLD, &msgStatus);
4357	if(infoFileExists)
4358	{
4359	MPI_Finalize();
4360	exit(-1);
4361	}
4362	}
4363	#endif
4364	}
4365
4366
4367	static void readData(analdef adef, rawdata rdta, cruncheddata cdta, tree tr)
4368	{
4369	INFILE = fopen(seq_file, "r");
4370
4371	if (!INFILE)
4372	{
4373	if(processID == 0)
4374	printf( "Could not open sequence file: %s\n", seq_file);
4375	errorExit(-1);
4376	}
4377	getinput(adef, rdta, cdta, tr);
4378
4379	fclose(INFILE);
4380	}
4381
4382
4383
4384	/*********************reading and initializing input ****************/
4385
4386
4387	/******************PRINTING various INFO ************************************/
4388
4389
4390	static void printModelAndProgramInfo(tree tr, analdef adef, int argc, char *argv[])
4391	{
4392	if(processID == 0)
4393	{
4394	int i, model;
4395	FILE *infoFile = fopen(infoFileName, "a");
4396	char modelType[128];
4397
4398	if(adef->useInvariant)
4399	strcpy(modelType, "GAMMA+P-Invar");
4400	else
4401	strcpy(modelType, "GAMMA");
4402
4403	printf("\n\nYou are using %s version %s released by Alexandros Stamatakis in %s\n", programName, programVersion, programDate);
4404	fprintf(infoFile, "\n\nYou are using %s version %s released by Alexandros Stamatakis in %s\n", programName, programVersion, programDate);
4405
4406	if(adef->mode == OPTIMIZE_RATES)
4407	{
4408	printf("\nAlignment has %d columns\n\n", tr->cdta->endsite);
4409	fprintf(infoFile, "\nAlignment has %d columns\n\n", tr->cdta->endsite);
4410	}
4411	else
4412	{
4413	printf("\nAlignment has %d distinct alignment patterns\n\n", tr->cdta->endsite);
4414	fprintf(infoFile, "\nAlignment has %d distinct alignment patterns\n\n", tr->cdta->endsite);
4415	}
4416
4417	if(adef->useInvariant)
4418	{
4419	printf("Found %d invariant alignment patterns that correspond to %d columns \n", tr->numberOfInvariableColumns, tr->weightOfInvariableColumns);
4420	fprintf(infoFile, "Found %d invariant alignment patterns that correspond to %d columns \n", tr->numberOfInvariableColumns, tr->weightOfInvariableColumns);
4421	}
4422
4423	printf("Proportion of gaps and completely undetermined characters in this alignment: %f\n", adef->gapyness);
4424	fprintf(infoFile, "Proportion of gaps and completely undetermined characters in this alignment: %f\n",
4425	adef->gapyness);
4426
4427	switch(adef->mode)
4428	{
4429	case ARNDT_MODE:
4430	printf("Arndt-Mode\n");
4431	fprintf(infoFile, "Arndt-Mode\n");
4432	break;
4433	case TREE_EVALUATION :
4434	printf("\nRAxML Model Optimization up to an accuracy of %f log likelihood units\n\n",
4435	adef->likelihoodEpsilon);
4436	fprintf(infoFile, "\nRAxML Model Optimization up to an accuracy of %f log likelihood units\n\n",
4437	adef->likelihoodEpsilon);
4438	break;
4439	case BIG_RAPID_MODE:
4440	if(adef->rapidBoot)
4441	{
4442	if(adef->allInOne)
4443	{
4444	printf("\nRAxML rapid bootstrapping and subsequent ML search\n\n");
4445	fprintf(infoFile, "\nRAxML rapid bootstrapping and subsequent ML search\n\n");
4446	}
4447	else
4448	{
4449	printf("\nRAxML rapid bootstrapping algorithm\n\n");
4450	fprintf(infoFile, "\nRAxML rapid bootstrapping algorithm\n\n");
4451	}
4452	}
4453	else
4454	{
4455	printf("\nRAxML rapid hill-climbing mode\n\n");
4456	fprintf(infoFile, "\nRAxML rapid hill-climbing mode\n\n");
4457	}
4458	break;
4459	case CALC_BIPARTITIONS:
4460	printf("\nRAxML Bipartition Computation: Drawing support values from trees in file %s onto tree in file %s\n\n",
4461	bootStrapFile, tree_file);
4462	fprintf(infoFile, "\nRAxML Bipartition Computation: Drawing support values from trees in file %s onto tree in file %s\n\n",
4463	bootStrapFile, tree_file);
4464	fclose(infoFile);
4465	return;
4466	case OPTIMIZE_RATES:
4467	if(!(adef->treeLength \|\| adef->computePerSiteLLs))
4468	{
4469	printf("\nRAxML optimization of per-site evolutionary rates\n\n");
4470	fprintf(infoFile,"\nRAxML optimization of per-site evolutionary rates\n\n");
4471	}
4472	if(adef->treeLength)
4473	{
4474	printf("\nRAxML optimization of per-site evolutionary rates and tree-length sliding window\n\n");
4475	fprintf(infoFile,"\nRAxML optimization of per-site evolutionary rates and tree-length sliding window\n\n");
4476	}
4477	if(adef->computePerSiteLLs)
4478	{
4479	printf("\nRAxML computation of per-site log likelihoods\n\n");
4480	fprintf(infoFile,"\nRAxML computation of per-site log likelihoods\n\n");
4481	}
4482	fclose(infoFile);
4483	return;
4484	case PARSIMONY_ADDITION:
4485	printf("\nRAxML stepwise MP addition to incomplete starting tree\n\n");
4486	fprintf(infoFile,"\nRAxML stepwise MP addition to incomplete starting tree\n\n");
4487	fclose(infoFile);
4488	return;
4489	case MEHRING_ALGO:
4490	printf("\nRAxML single-sequence position determination algorithm\n\n");
4491	fprintf(infoFile,"\nRAxML single-sequence position determination algorithm\n\n");
4492	break;
4493	default:
4494	printf("Oups, forgot to implement mode description %d exiting\n", adef->mode);
4495	exit(-1);
4496	}
4497
4498	if(tr->NumberOfModels > 1)
4499	{
4500	if(adef->perGeneBranchLengths)
4501	{
4502	printf("Partitioned Data Mode: Using %d distinct models/partitions with individual per partition branch length optimization\n", tr->NumberOfModels);
4503	fprintf(infoFile, "Partitioned Data Mode: Using %d distinct models/partitions with individual per partition branch length optimization\n", tr->NumberOfModels);
4504	printf( "\n\n");
4505	fprintf(infoFile, "\n\n");
4506	}
4507	else
4508	{
4509	printf("Partitioned Data Mode: Using %d distinct models/partitions with joint branch length optimization\n",
4510	tr->NumberOfModels);
4511	fprintf(infoFile, "Partitioned Data Mode: Using %d distinct models/partitions with joint branch length optimization\n",
4512	tr->NumberOfModels);
4513	printf( "\n\n");
4514	fprintf(infoFile, "\n\n");
4515	}
4516	}
4517
4518	if(adef->rapidBoot)
4519	{
4520	if(adef->allInOne)
4521	{
4522	printf("\nExecuting %d rapid bootstrap inferences and thereafter a thorough ML search \n\n", adef->multipleRuns);
4523	fprintf(infoFile, "\nExecuting %d rapid bootstrap inferences and thereafter a thorough ML search \n\n", adef->multipleRuns);
4524	}
4525	else
4526	{
4527	printf("\nExecuting %d rapid bootstrap inferences\n\n", adef->multipleRuns);
4528	fprintf(infoFile, "\nExecuting %d rapid bootstrap inferences\n\n", adef->multipleRuns);
4529	}
4530	}
4531	else
4532	{
4533	if(adef->boot)
4534	{
4535	if(adef->multipleRuns > 1)
4536	{
4537	printf("Executing %d non-parametric bootstrap inferences\n\n", adef->multipleRuns);
4538	fprintf(infoFile, "Executing %d non-parametric bootstrap inferences\n\n", adef->multipleRuns);
4539	}
4540	else
4541	{
4542	printf("Executing %d non-parametric bootstrap inference\n\n", adef->multipleRuns);
4543	fprintf(infoFile, "Executing %d non-parametric bootstrap inference\n\n", adef->multipleRuns);
4544	}
4545	}
4546	else
4547	{
4548	char treeType[1024];
4549
4550	if(adef->restart)
4551	strcpy(treeType, "user-specifed");
4552	else
4553	{
4554	if(adef->randomStartingTree)
4555	strcpy(treeType, "distinct complete random");
4556	else
4557	strcpy(treeType, "distinct randomized MP");
4558	}
4559
4560
4561	if(adef->multipleRuns > 1)
4562	{
4563	printf("Executing %d inferences on the original alignment using %d %s trees\n\n",
4564	adef->multipleRuns, adef->multipleRuns, treeType);
4565	fprintf(infoFile, "Executing %d inferences on the original alignment using %d %s trees\n\n",
4566	adef->multipleRuns, adef->multipleRuns, treeType);
4567	}
4568	else
4569	{
4570	printf("Executing %d inference on the original alignment using a %s tree\n\n",
4571	adef->multipleRuns, treeType);
4572	fprintf(infoFile, "Executing %d inference on the original alignment using a %s tree\n\n",
4573	adef->multipleRuns, treeType);
4574	}
4575	}
4576	}
4577
4578	if(tr->rateHetModel == GAMMA \|\| tr->rateHetModel == GAMMA_I)
4579	{
4580	printf( "All free model parameters will be estimated by RAxML\n");
4581	printf( "%s model of rate heteorgeneity, ML estimate of alpha-parameter\n", modelType);
4582	printf( "%s Model parameters will be estimated up to an accuracy of %2.10f Log Likelihood units\n\n",
4583	modelType, adef->likelihoodEpsilon);
4584	fprintf(infoFile, "All free model parameters will be estimated by RAxML\n");
4585	fprintf(infoFile, "%s model of rate heterogeneity, ML estimate of alpha-parameter\n", modelType);
4586	fprintf(infoFile, "%s Model parameters will be estimated up to an accuracy of %2.10f Log Likelihood units\n\n",
4587	modelType, adef->likelihoodEpsilon);
4588	}
4589	else
4590	{
4591	if(adef->useMixedModel)
4592	{
4593	printf( "All free model parameters will be estimated by RAxML\n");
4594	printf( "ML estimate of %d per site rate categories\n", adef->categories);
4595	printf( "Likelihood of final tree will be evaluated and optimized under %s\n", modelType);
4596	printf( "Final model parameters will be estimated up to an accuracy of %2.10f Log Likelihood units\n\n",
4597	adef->likelihoodEpsilon);
4598
4599	fprintf(infoFile, "All free model parameters will be estimated by RAxML\n");
4600	fprintf(infoFile, "ML estimate of %d per site rate categories\n", adef->categories);
4601	fprintf(infoFile, "Likelihood of final tree will be evaluated and optimized under %s\n", modelType);
4602	fprintf(infoFile, "Final model parameters will be estimated up to an accuracy of %2.10f Log Likelihood units\n\n",
4603	adef->likelihoodEpsilon);
4604	}
4605	else
4606	{
4607	printf( "Approximation of rate heterogeneity only!\n");
4608	printf( "All free model parameters will be estimated by RAxML\n");
4609	printf( "ML estimate of %d per site rate categories\n", adef->categories);
4610	printf( "WARNING: CAT-based likelihood values should NEVER be used to COMPARE trees!\n\n");
4611
4612	fprintf(infoFile, "Approximation of rate heterogeneity only!\n");
4613	fprintf(infoFile, "All free model parameters will be estimated by RAxML\n");
4614	fprintf(infoFile, "ML estimate of %d per site rate categories\n", adef->categories);
4615	fprintf(infoFile, "WARNING: CAT-based likelihood values should NEVER be used to COMPARE trees!\n\n");
4616	}
4617	}
4618
4619
4620	for(model = 0; model < tr->NumberOfModels; model++)
4621	{
4622	printf("Partition: %d\n", model);
4623	printf("Name: %s\n", tr->partitionData[model].partitionName);
4624
4625	fprintf(infoFile, "Partition: %d\n", model);
4626	fprintf(infoFile, "Name: %s\n", tr->partitionData[model].partitionName);
4627
4628	switch(tr->partitionData[model].dataType)
4629	{
4630	case DNA_DATA:
4631	printf("DataType: DNA\n");
4632	printf("Substitution Matrix: GTR\n");
4633
4634	if(adef->boot == 0)
4635	{
4636	printf("Empirical Base Frequencies:\n");
4637	printf("pi(A): %f pi(C): %f pi(G): %f pi(T): %f",
4638	tr->frequencies_DNA[model * 4 + 0], tr->frequencies_DNA[model * 4 + 1],
4639	tr->frequencies_DNA[model * 4 + 2], tr->frequencies_DNA[model * 4 + 3]);
4640	}
4641	else
4642	{
4643	printf("Empirical Base Frequencies will not be printed for Bootstrapping\n");
4644	}
4645
4646	fprintf(infoFile, "DataType: DNA\n");
4647	fprintf(infoFile, "Substitution Matrix: GTR\n");
4648
4649	if(adef->boot == 0)
4650	{
4651	fprintf(infoFile, "Empirical Base Frequencies:\n");
4652	fprintf(infoFile, "pi(A): %f pi(C): %f pi(G): %f pi(T): %f",
4653	tr->frequencies_DNA[model * 4 + 0], tr->frequencies_DNA[model * 4 + 1],
4654	tr->frequencies_DNA[model * 4 + 2], tr->frequencies_DNA[model * 4 + 3]);
4655	}
4656	else
4657	{
4658	fprintf(infoFile, "Empirical Base Frequencies will not be printed for Bootstrapping\n");
4659	}
4660	break;
4661	case AA_DATA:
4662	{
4663	char *protStrings[10] = {"DAYHOFF", "DCMUT", "JTT", "MTREV", "WAG", "RTREV", "CPREV", "VT", "BLOSUM62", "MTMAM"};
4664	char basesPROT[20] = {'A', 'R', 'N' , 'D', 'C', 'Q','E','G','H','I','L','K','M','F','P','S','T','W','Y','V'};
4665
4666	assert(tr->partitionData[model].protModels >= 0 && tr->partitionData[model].protModels < 10);
4667
4668	printf("DataType: AA\n");
4669	printf("Substitution Matrix: %s\n", protStrings[tr->partitionData[model].protModels]);
4670	printf("%s Base Frequencies:\n", (tr->partitionData[model].protFreqs == 1)?"Empirical":"Fixed");
4671
4672	if(adef->boot == 0)
4673	{
4674	int k;
4675
4676	for(k = 0; k < 20; k++)
4677	{
4678	if(k % 4 == 0 && k > 0)
4679	printf("\n");
4680
4681	printf("pi(%c): %f ", basesPROT[k], tr->frequencies_AA[model * 20 + k]);
4682	}
4683	}
4684	else
4685	{
4686	printf("Base Frequencies will not be printed for Bootstrapping\n");
4687	}
4688
4689	fprintf(infoFile, "DataType: AA\n");
4690	fprintf(infoFile, "Substitution Matrix: %s\n", protStrings[tr->partitionData[model].protModels]);
4691	fprintf(infoFile, "%s Base Frequencies:\n", (tr->partitionData[model].protFreqs == 1)?"Empirical":"Fixed");
4692
4693	if(adef->boot == 0)
4694	{
4695	int k;
4696
4697	for(k = 0; k < 20; k++)
4698	{
4699	if(k % 4 == 0)
4700	fprintf(infoFile, "\n");
4701
4702	fprintf(infoFile, "pi(%c): %f ", basesPROT[k], tr->frequencies_AA[model * 20 + k]);
4703	}
4704	}
4705	else
4706	{
4707	fprintf(infoFile, "Base Frequencies will not be printed for Bootstrapping\n");
4708	}
4709
4710
4711	}
4712	break;
4713	default:
4714	assert(0);
4715	}
4716
4717	printf("\n\n\n");
4718	fprintf(infoFile,"\n\n\n");
4719	}
4720
4721	printf("\n");
4722	fprintf(infoFile, "\n");
4723
4724	fprintf(infoFile,"RAxML was called as follows:\n\n");
4725	for(i = 0; i < argc; i++)
4726	fprintf(infoFile,"%s ", argv[i]);
4727	fprintf(infoFile,"\n\n\n");
4728
4729	fclose(infoFile);
4730	}
4731	}
4732
4733	void printResult(tree tr, analdef adef, boolean finalPrint)
4734	{
4735	FILE *logFile;
4736	char temporaryFileName[1024] = "", treeID[64] = "";
4737
4738	strcpy(temporaryFileName, resultFileName);
4739
4740	switch(adef->mode)
4741	{
4742	case TREE_EVALUATION:
4743
4744
4745	Tree2String(tr->tree_string, tr, tr->start->back, TRUE, TRUE, FALSE, FALSE, finalPrint, adef, SUMMARIZE_LH);
4746
4747	logFile = fopen(temporaryFileName, "w");
4748	fprintf(logFile, "%s", tr->tree_string);
4749	fclose(logFile);
4750
4751	if(adef->perGeneBranchLengths)
4752	printTreePerGene(tr, adef, temporaryFileName, "w");
4753
4754
4755	break;
4756	case BIG_RAPID_MODE:
4757	if(!adef->boot)
4758	{
4759	if(adef->multipleRuns > 1)
4760	{
4761	sprintf(treeID, "%d", tr->treeID);
4762	strcat(temporaryFileName, ".RUN.");
4763	strcat(temporaryFileName, treeID);
4764	}
4765
4766	if((adef->model == M_GTRCAT \|\| adef->model == M_PROTCAT) && (adef->useMixedModel == 0))
4767	{
4768	Tree2String(tr->tree_string, tr, tr->start->back, FALSE, TRUE, FALSE, FALSE, finalPrint, adef, NO_BRANCHES);
4769	logFile = fopen(temporaryFileName, "w");
4770	fprintf(logFile, "%s", tr->tree_string);
4771	fclose(logFile);
4772	}
4773	else
4774	{
4775	if(finalPrint)
4776	{
4777	Tree2String(tr->tree_string, tr, tr->start->back, TRUE, TRUE, FALSE, FALSE, finalPrint, adef,
4778	SUMMARIZE_LH);
4779
4780	logFile = fopen(temporaryFileName, "w");
4781	fprintf(logFile, "%s", tr->tree_string);
4782	fclose(logFile);
4783
4784	if(adef->perGeneBranchLengths)
4785	printTreePerGene(tr, adef, temporaryFileName, "w");
4786	}
4787	else
4788	{
4789	Tree2String(tr->tree_string, tr, tr->start->back, FALSE, TRUE, FALSE, FALSE, finalPrint, adef,
4790	NO_BRANCHES);
4791	logFile = fopen(temporaryFileName, "w");
4792	fprintf(logFile, "%s", tr->tree_string);
4793	fclose(logFile);
4794	}
4795	}
4796	}
4797	break;
4798	default:
4799	printf("FATAL ERROR call to printResult from undefined STATE %d\n", adef->mode);
4800	exit(-1);
4801	break;
4802	}
4803	}
4804
4805	void printBootstrapResult(tree tr, analdef adef, boolean finalPrint)
4806	{
4807	if(processID == 0)
4808	{
4809	FILE *logFile;
4810
4811	if(adef->mode == BIG_RAPID_MODE && (adef->boot \|\| adef->rapidBoot))
4812	{
4813	#ifndef PARALLEL
4814	if(adef->bootstrapBranchLengths)
4815	{
4816	Tree2String(tr->tree_string, tr, tr->start->back, TRUE, TRUE, FALSE, FALSE, finalPrint, adef, SUMMARIZE_LH);
4817	logFile = fopen(bootstrapFileName, "a");
4818	fprintf(logFile, "%s", tr->tree_string);
4819	fclose(logFile);
4820	if(adef->perGeneBranchLengths)
4821	printTreePerGene(tr, adef, bootstrapFileName, "a");
4822	}
4823	else
4824	{
4825	Tree2String(tr->tree_string, tr, tr->start->back, FALSE, TRUE, FALSE, FALSE, finalPrint, adef, NO_BRANCHES);
4826	logFile = fopen(bootstrapFileName, "a");
4827	fprintf(logFile, "%s", tr->tree_string);
4828	fclose(logFile);
4829	}
4830	#else
4831	logFile = fopen(bootstrapFileName, "a");
4832	fprintf(logFile, "%s", tr->tree_string);
4833	fclose(logFile);
4834	#endif
4835	}
4836	else
4837	{
4838	printf("FATAL ERROR in printBootstrapResult\n");
4839	exit(-1);
4840	}
4841	}
4842	}
4843
4844
4845
4846	void printBipartitionResult(tree tr, analdef adef, boolean finalPrint)
4847	{
4848	if(processID == 0 \|\| adef->allInOne)
4849	{
4850	FILE *logFile;
4851
4852	Tree2String(tr->tree_string, tr, tr->start->back, FALSE, TRUE, FALSE, TRUE, finalPrint, adef, NO_BRANCHES);
4853	logFile = fopen(bipartitionsFileName, "a");
4854	fprintf(logFile, "%s", tr->tree_string);
4855	fclose(logFile);
4856	}
4857	}
4858
4859
4860
4861	void printLog(tree tr, analdef adef, boolean finalPrint)
4862	{
4863	FILE *logFile;
4864	char temporaryFileName[1024] = "", checkPoints[1024] = "", treeID[64] = "";
4865	double lh, t;
4866
4867	lh = tr->likelihood;
4868	t = gettime() - masterTime;
4869
4870	strcpy(temporaryFileName, logFileName);
4871	strcpy(checkPoints, checkpointFileName);
4872
4873	switch(adef->mode)
4874	{
4875	case TREE_EVALUATION:
4876	logFile = fopen(temporaryFileName, "a");
4877
4878	printf("%f %f\n", t, lh);
4879	fprintf(logFile, "%f %f\n", t, lh);
4880
4881	fclose(logFile);
4882	break;
4883	case BIG_RAPID_MODE:
4884	if(adef->boot \|\| adef->rapidBoot)
4885	{
4886	/* testing only printf("%f %f\n", t, lh);*/
4887	/* NOTHING PRINTED so far */
4888	}
4889	else
4890	{
4891	if(adef->multipleRuns > 1)
4892	{
4893	sprintf(treeID, "%d", tr->treeID);
4894	strcat(temporaryFileName, ".RUN.");
4895	strcat(temporaryFileName, treeID);
4896
4897	strcat(checkPoints, ".RUN.");
4898	strcat(checkPoints, treeID);
4899	}
4900
4901
4902	if(!adef->checkpoints)
4903	{
4904	logFile = fopen(temporaryFileName, "a");
4905	#ifndef PARALLEL
4906	printf("%f %f\n", t, lh);
4907	#endif
4908	fprintf(logFile, "%f %f\n", t, lh);
4909
4910	fclose(logFile);
4911	}
4912	else
4913	{
4914	logFile = fopen(temporaryFileName, "a");
4915	#ifndef PARALLEL
4916	printf("%f %f %d\n", t, lh, tr->checkPointCounter);
4917	#endif
4918	fprintf(logFile, "%f %f %d\n", t, lh, tr->checkPointCounter);
4919
4920	fclose(logFile);
4921
4922	strcat(checkPoints, ".");
4923
4924	sprintf(treeID, "%d", tr->checkPointCounter);
4925	strcat(checkPoints, treeID);
4926
4927	Tree2String(tr->tree_string, tr, tr->start->back, FALSE, TRUE, FALSE, FALSE, finalPrint, adef, NO_BRANCHES);
4928
4929	logFile = fopen(checkPoints, "a");
4930	fprintf(logFile, "%s", tr->tree_string);
4931	fclose(logFile);
4932
4933	tr->checkPointCounter++;
4934	}
4935	}
4936	break;
4937	default:
4938	printf("FATAL ERROR call to printLog from undefined STATE %d\n", adef->mode);
4939	exit(-1);
4940	break;
4941	}
4942	}
4943
4944
4945
4946	void printStartingTree(tree tr, analdef adef, boolean finalPrint)
4947	{
4948	if(adef->boot)
4949	{
4950	/* not printing starting trees for bootstrap */
4951	}
4952	else
4953	{
4954	FILE *treeFile;
4955	char temporaryFileName[1024] = "", treeID[64] = "";
4956
4957	Tree2String(tr->tree_string, tr, tr->start->back, FALSE, TRUE, FALSE, FALSE, finalPrint, adef, NO_BRANCHES);
4958
4959	if(adef->randomStartingTree)
4960	strcpy(temporaryFileName, randomFileName);
4961	else
4962	strcpy(temporaryFileName, permFileName);
4963
4964	if(adef->multipleRuns > 1)
4965	{
4966	sprintf(treeID, "%d", tr->treeID);
4967	strcat(temporaryFileName, ".RUN.");
4968	strcat(temporaryFileName, treeID);
4969	}
4970
4971	treeFile = fopen(temporaryFileName, "a");
4972	fprintf(treeFile, "%s", tr->tree_string);
4973	fclose(treeFile);
4974	}
4975	}
4976
4977	void writeInfoFile(analdef adef, tree tr, double t)
4978	{
4979	if(processID == 0)
4980	{
4981	FILE *infoFile = fopen(infoFileName, "a");
4982
4983	switch(adef->mode)
4984	{
4985	case TREE_EVALUATION:
4986	break;
4987	case BIG_RAPID_MODE:
4988	if(adef->boot \|\| adef->rapidBoot)
4989	{
4990	if(!adef->initialSet)
4991	{
4992	fprintf(infoFile, "Bootstrap[%d]: Time %f bootstrap likelihood %f, best rearrangement setting %d\n", tr->treeID, t, tr->likelihood, adef->bestTrav);
4993	printf("Bootstrap[%d]: Time %f bootstrap likelihood %f, best rearrangement setting %d\n", tr->treeID, t, tr->likelihood, adef->bestTrav);
4994	}
4995	else
4996	{
4997	fprintf(infoFile, "Bootstrap[%d]: Time %f bootstrap likelihood %f\n", tr->treeID, t, tr->likelihood);
4998	printf("Bootstrap[%d]: Time %f bootstrap likelihood %f\n", tr->treeID, t, tr->likelihood);
4999	}
5000	}
5001	else
5002	{
5003	if((adef->model == M_GTRCAT \|\| adef->model == M_PROTCAT) && !adef->useMixedModel)
5004	{
5005	if(!adef->initialSet)
5006	fprintf(infoFile, "Inference[%d]: Time %f CAT-likelihood %f, best rearrangement setting %d\n", tr->treeID, t, tr->likelihood, adef->bestTrav);
5007	else
5008	fprintf(infoFile, "Inference[%d]: Time %f CAT-likelihood %f\n", tr->treeID, t, tr->likelihood);
5009	}
5010	else
5011	{
5012	int model;
5013	char modelType[128];
5014
5015	if(adef->useInvariant)
5016	strcpy(modelType, "GAMMA+P-Invar");
5017	else
5018	strcpy(modelType, "GAMMA");
5019
5020	if(!adef->initialSet)
5021	fprintf(infoFile, "Inference[%d]: Time %f %s-likelihood %f, best rearrangement setting %d, ",
5022	tr->treeID, t, modelType, tr->likelihood, adef->bestTrav);
5023	else
5024	fprintf(infoFile, "Inference[%d]: Time %f %s-likelihood %f, ",
5025	tr->treeID, t, modelType, tr->likelihood);
5026
5027	for(model = 0; model < tr->NumberOfModels; model++)
5028	{
5029	fprintf(infoFile, "alpha[%d]: %f ", model, tr->alphas[model]);
5030	if(adef->useInvariant)
5031	fprintf(infoFile, "invar[%d]: %f ", model, tr->invariants[model]);
5032	#ifndef PARALLEL
5033	if(adef->model == M_GTRCAT \|\| adef->model == M_GTRGAMMA)
5034	{
5035	int k;
5036
5037	fprintf(infoFile, "rates[%d] ac ag at cg ct gt: ",model);
5038	for(k = 0; k < DNA_RATES; k++)
5039	fprintf(infoFile, "%f ", tr->initialRates_DNA[model * DNA_RATES + k]);
5040	}
5041	fprintf(infoFile, "1.0 ");
5042	#endif
5043	}
5044
5045	fprintf(infoFile, "\n");
5046	}
5047	}
5048	break;
5049	default:
5050	assert(0);
5051	}
5052
5053	fclose(infoFile);
5054	}
5055	}
5056
5057	static void finalizeInfoFile(tree tr, analdef adef)
5058	{
5059	if(processID == 0)
5060	{
5061	FILE *infoFile = fopen(infoFileName, "a");
5062	double t;
5063	int model;
5064
5065	t = gettime() - masterTime;
5066
5067	switch(adef->mode)
5068	{
5069	case TREE_EVALUATION :
5070	printf("\n\nOverall Time for Tree Evaluation %f\n", t);
5071	printf("Final GAMMA likelihood: %f\n", tr->likelihood);
5072
5073	fprintf(infoFile, "\n\nOverall Time for Tree Evaluation %f\n", t);
5074	fprintf(infoFile, "Final GAMMA likelihood: %f\n", tr->likelihood);
5075
5076	{
5077	int
5078	params,
5079	paramsBrLen;
5080
5081	if(tr->NumberOfModels == 1)
5082	{
5083	if(adef->useInvariant)
5084	{
5085	params = 1 /* INVAR / + 5 / RATES / + 3 / freqs / + 1 / alpha */;
5086	paramsBrLen = 1 /* INVAR / + 5 / RATES / + 3 / freqs / + 1 / alpha */ +
5087	(2 * tr->mxtips - 3);
5088	}
5089	else
5090	{
5091	params = 5 /* RATES / + 3 / freqs / + 1 / alpha */;
5092	paramsBrLen = 5 /* RATES / + 3 / freqs / + 1 / alpha */ +
5093	(2 * tr->mxtips - 3);
5094	}
5095	}
5096	else
5097	{
5098	if(tr->multiBranch)
5099	{
5100	if(adef->useInvariant)
5101	{
5102	params = tr->NumberOfModels * (1 /* INVAR / + 5 / RATES / + 3 / freqs / + 1 / alpha */);
5103	paramsBrLen = tr->NumberOfModels * (1 /* INVAR / + 5 / RATES / + 3 / freqs / + 1 / alpha */ +
5104	(2 * tr->mxtips - 3));
5105	}
5106	else
5107	{
5108	params = tr->NumberOfModels * (5 /* RATES / + 3 / freqs / + 1 / alpha */);
5109	paramsBrLen = tr->NumberOfModels * (5 /* RATES / + 3 / freqs / + 1 / alpha */ +
5110	(2 * tr->mxtips - 3));
5111	}
5112	}
5113	else
5114	{
5115	if(adef->useInvariant)
5116	{
5117	params = tr->NumberOfModels * (1 /* INVAR / + 5 / RATES / + 3 / freqs / + 1 / alpha */);
5118	paramsBrLen = tr->NumberOfModels * (1 /* INVAR / + 5 / RATES / + 3 / freqs / + 1 / alpha */)
5119	+ (2 * tr->mxtips - 3);
5120	}
5121	else
5122	{
5123	params = tr->NumberOfModels * (5 /* RATES / + 3 / freqs / + 1 / alpha */);
5124	paramsBrLen = tr->NumberOfModels * (5 /* RATES / + 3 / freqs / + 1 / alpha */)
5125	+ (2 * tr->mxtips - 3);
5126	}
5127
5128	}
5129	}
5130
5131	if(!tr->mixedData && tr->partitionData[0].dataType == DNA_DATA)
5132	{
5133	printf("Number of free parameters for AIC-TEST(BR-LEN): %d\n", paramsBrLen);
5134	printf("Number of free parameters for AIC-TEST(NO-BR-LEN): %d\n", params);
5135	fprintf(infoFile, "Number of free parameters for AIC-TEST(BR-LEN): %d\n", paramsBrLen);
5136	fprintf(infoFile, "Number of free parameters for AIC-TEST(NO-BR-LEN): %d\n", params);
5137	}
5138
5139	}
5140
5141	printf("\n\n");
5142	fprintf(infoFile, "\n\n");
5143
5144	for(model = 0; model < tr->NumberOfModels; model++)
5145	{
5146	double tl;
5147	char typeOfData[1024];
5148
5149	switch(tr->partitionData[model].dataType)
5150	{
5151	case AA_DATA:
5152	strcpy(typeOfData,"AA");
5153	break;
5154	case DNA_DATA:
5155	strcpy(typeOfData,"DNA");
5156	break;
5157	default:
5158	assert(0);
5159	}
5160
5161	fprintf(infoFile, "Model Parameters of Partition %d, Name: %s, Type of Data: %s\n",
5162	model, tr->partitionData[model].partitionName, typeOfData);
5163	fprintf(infoFile, "alpha: %f\n", tr->alphas[model]);
5164
5165	printf("Model Parameters of Partition %d, Name: %s, Type of Data: %s\n",
5166	model, tr->partitionData[model].partitionName, typeOfData);
5167	printf("alpha: %f\n", tr->alphas[model]);
5168
5169	if(adef->useInvariant)
5170	{
5171	fprintf(infoFile, "invar: %f\n", tr->invariants[model]);
5172	printf("invar: %f\n", tr->invariants[model]);
5173	}
5174
5175	if(adef->perGeneBranchLengths)
5176	tl = treeLength(tr, model);
5177	else
5178	tl = treeLength(tr, 0);
5179
5180	fprintf(infoFile, "Tree-Length: %f\n", tl);
5181	printf("Tree-Length: %f\n", tl);
5182
5183
5184
5185	switch(tr->partitionData[model].dataType)
5186	{
5187	case AA_DATA:
5188	break;
5189	case DNA_DATA:
5190	{
5191	int k;
5192	char *names[6] = {"a<->c", "a<->g", "a<->t", "c<->g", "c<->t", "g<->t"};
5193	for(k = 0; k < DNA_RATES; k++)
5194	{
5195	fprintf(infoFile, "rate %s: %f\n", names[k], tr->initialRates_DNA[model * DNA_RATES + k]);
5196	printf("rate %s: %f\n", names[k], tr->initialRates_DNA[model * DNA_RATES + k]);
5197	}
5198
5199	fprintf(infoFile, "rate %s: %f\n", names[5], 1.0);
5200	printf("rate %s: %f\n", names[5], 1.0);
5201	}
5202	break;
5203	default:
5204	assert(0);
5205	}
5206
5207	fprintf(infoFile, "\n");
5208	printf("\n");
5209	}
5210
5211	printf("Final tree written to: %s\n", resultFileName);
5212	printf("Execution Log File written to: %s\n", logFileName);
5213
5214	fprintf(infoFile, "Final tree written to: %s\n", resultFileName);
5215	fprintf(infoFile, "Execution Log File written to: %s\n", logFileName);
5216
5217	break;
5218	case BIG_RAPID_MODE:
5219	if(adef->boot)
5220	{
5221	printf("\n\nOverall Time for %d Bootstraps %f\n", adef->multipleRuns, t);
5222	printf("\n\nAverage Time per Bootstrap %f\n", (double)(t/((double)adef->multipleRuns)));
5223	printf("All %d bootstrapped trees written to: %s\n", adef->multipleRuns, bootstrapFileName);
5224
5225	fprintf(infoFile, "\n\nOverall Time for %d Bootstraps %f\n", adef->multipleRuns, t);
5226	fprintf(infoFile, "Average Time per Bootstrap %f\n", (double)(t/((double)adef->multipleRuns)));
5227	fprintf(infoFile, "\n\nAll %d bootstrapped trees written to: %s\n", adef->multipleRuns, bootstrapFileName);
5228	}
5229	else
5230	{
5231	if(adef->multipleRuns > 1)
5232	{
5233	double avgLH = 0;
5234	double bestLH = unlikely;
5235	int i, bestI = 0;
5236
5237	for(i = 0; i < adef->multipleRuns; i++)
5238	{
5239	avgLH += tr->likelihoods[i];
5240	if(tr->likelihoods[i] > bestLH)
5241	{
5242	bestLH = tr->likelihoods[i];
5243	bestI = i;
5244	}
5245	}
5246	avgLH /= ((double)adef->multipleRuns);
5247
5248	printf("\n\nOverall Time for %d Inferences %f\n", adef->multipleRuns, t);
5249	printf("Average Time per Inference %f\n", (double)(t/((double)adef->multipleRuns)));
5250	printf("Average Likelihood : %f\n", avgLH);
5251	printf("\n");
5252	printf("Best Likelihood in run number %d: likelihood %f\n\n", bestI, bestLH);
5253
5254	if(adef->checkpoints)
5255	printf("Checkpoints written to: %s.RUN.%d.* to %d.*\n", checkpointFileName, 0, adef->multipleRuns - 1);
5256	if(!adef->restart)
5257	{
5258	if(adef->randomStartingTree)
5259	printf("Random starting trees written to: %s.RUN.%d to %d\n", randomFileName, 0, adef->multipleRuns - 1);
5260	else
5261	printf("Parsimony starting trees written to: %s.RUN.%d to %d\n", permFileName, 0, adef->multipleRuns - 1);
5262	}
5263	printf("Final trees written to: %s.RUN.%d to %d\n", resultFileName, 0, adef->multipleRuns - 1);
5264	printf("Execution Log Files written to: %s.RUN.%d to %d\n", logFileName, 0, adef->multipleRuns - 1);
5265	printf("Execution information file written to: %s\n", infoFileName);
5266
5267
5268	fprintf(infoFile, "\n\nOverall Time for %d Inferences %f\n", adef->multipleRuns, t);
5269	fprintf(infoFile, "Average Time per Inference %f\n", (double)(t/((double)adef->multipleRuns)));
5270	fprintf(infoFile, "Average Likelihood : %f\n", avgLH);
5271	fprintf(infoFile, "\n");
5272	fprintf(infoFile, "Best Likelihood in run number %d: likelihood %f\n\n", bestI, bestLH);
5273	if(adef->checkpoints)
5274	fprintf(infoFile, "Checkpoints written to: %s.RUN.%d.* to %d.*\n", checkpointFileName, 0, adef->multipleRuns - 1);
5275	if(!adef->restart)
5276	{
5277	if(adef->randomStartingTree)
5278	fprintf(infoFile, "Random starting trees written to: %s.RUN.%d to %d\n", randomFileName, 0, adef->multipleRuns - 1);
5279	else
5280	fprintf(infoFile, "Parsimony starting trees written to: %s.RUN.%d to %d\n", permFileName, 0, adef->multipleRuns - 1);
5281	}
5282	fprintf(infoFile, "Final trees written to: %s.RUN.%d to %d\n", resultFileName, 0, adef->multipleRuns - 1);
5283	fprintf(infoFile, "Execution Log Files written to: %s.RUN.%d to %d\n", logFileName, 0, adef->multipleRuns - 1);
5284	fprintf(infoFile, "Execution information file written to: %s\n", infoFileName);
5285
5286	}
5287	else
5288	{
5289	printf("\n\nOverall Time for 1 Inference %f\n", t);
5290	printf("Likelihood : %f\n", tr->likelihood);
5291	printf("\n\n");
5292
5293	if(adef->checkpoints)
5294	printf("Checkpoints written to: %s.*\n", checkpointFileName);
5295	if(!adef->restart)
5296	{
5297	if(adef->randomStartingTree)
5298	printf("Random starting tree written to: %s\n", randomFileName);
5299	else
5300	printf("Parsimony starting tree written to: %s\n", permFileName);
5301	}
5302	printf("Final tree written to: %s\n", resultFileName);
5303	printf("Execution Log File written to: %s\n", logFileName);
5304	printf("Execution information file written to: %s\n",infoFileName);
5305
5306
5307
5308	fprintf(infoFile, "\n\nOverall Time for 1 Inference %f\n", t);
5309	fprintf(infoFile, "Likelihood : %f\n", tr->likelihood);
5310	fprintf(infoFile, "\n\n");
5311
5312	if(adef->checkpoints)
5313	fprintf(infoFile, "Checkpoints written to: %s.*\n", checkpointFileName);
5314	if(!adef->restart)
5315	{
5316	if(adef->randomStartingTree)
5317	fprintf(infoFile, "Random starting tree written to: %s\n", randomFileName);
5318	else
5319	fprintf(infoFile, "Parsimony starting tree written to: %s\n", permFileName);
5320	}
5321	fprintf(infoFile, "Final tree written to: %s\n", resultFileName);
5322	fprintf(infoFile, "Execution Log File written to: %s\n", logFileName);
5323	fprintf(infoFile, "Execution information file written to: %s\n",infoFileName);
5324
5325	}
5326	}
5327
5328	break;
5329	case CALC_BIPARTITIONS:
5330	printf("\n\nTime for Computation of Bipartitions %f\n", t);
5331	printf("Tree with bipartitions written to file: %s\n", bipartitionsFileName);
5332	printf("Execution information file written to : %s\n",infoFileName);
5333
5334
5335	fprintf(infoFile, "\n\nTime for Computation of Bipartitions %f\n", t);
5336	fprintf(infoFile, "Tree with bipartitions written to file: %s\n", bipartitionsFileName);
5337
5338
5339	break;
5340	case OPTIMIZE_RATES:
5341	if(! (adef->computePerSiteLLs \|\| adef->treeLength))
5342	{
5343	printf("\n\nTime for Optimization of per-site rates %f\n", t);
5344	printf("Optimized rates written to file: %s\n", ratesFileName);
5345	printf("Execution information file written to : %s\n",infoFileName);
5346
5347
5348	fprintf(infoFile, "\n\nTime for Optimization of per-site rates %f\n", t);
5349	fprintf(infoFile, "Optimized rates written to file: %s\n", ratesFileName);
5350	}
5351
5352	if(adef->treeLength)
5353	{
5354	printf("\n\nTime for Optimization of per-site rates and sliding window tree length %f\n", t);
5355	printf("Optimized rates written to file: %s, \n sliding window data written to %s and %s\n", ratesFileName, lengthFileName, lengthFileNameModel);
5356	printf("Execution information file written to : %s\n",infoFileName);
5357
5358
5359	fprintf(infoFile, "\n\nTime for Optimization of per-site rates and sliding window tree length %f\n", t);
5360	fprintf(infoFile, "Optimized rates written to file: %s, \n sliding window data written to %s and %s\n", ratesFileName, lengthFileName, lengthFileNameModel);
5361	}
5362
5363	if(adef->computePerSiteLLs)
5364	{
5365	printf("\n\nTime for Optimization of per-site log likelihoods %f\n", t);
5366	printf("Per-site Log Likelihoods written to File %s in Tree-Puzzle format\n", perSiteLLsFileName);
5367	printf("Execution information file written to : %s\n",infoFileName);
5368
5369
5370	fprintf(infoFile, "\n\nTime for Optimization of per-site log likelihoods %f\n", t);
5371	fprintf(infoFile, "Per-site Log Likelihoods written to File %s in Tree-Puzzle format\n", perSiteLLsFileName);
5372	}
5373
5374	break;
5375	case PARSIMONY_ADDITION:
5376	printf("\n\nTime for MP stepwise addition %f\n", t);
5377	printf("Execution information file written to : %s\n",infoFileName);
5378	printf("Complete parsimony tree written to: %s\n", permFileName);
5379
5380
5381
5382	fprintf(infoFile, "\n\nTime for MP stepwise addition %f\n", t);
5383	fprintf(infoFile, "Complete parsimony tree written to: %s\n", permFileName);
5384
5385
5386	break;
5387	default:
5388	assert(0);
5389	}
5390	fclose(infoFile);
5391	}
5392
5393	}
5394
5395
5396
5397	/******************PRINTING various INFO ************************************/
5398
5399	/************************************************************************************/
5400
5401	static void computeLHTest(tree tr, analdef adef, char *bootStrapFileName)
5402	{
5403	int numberOfTrees = 0, i;
5404	char ch;
5405	double bestLH, currentLH, weightSum = 0.0;
5406	double bestVector, otherVector;
5407
5408
5409	bestVector = (double)malloc(sizeof(double) tr->cdta->endsite);
5410	otherVector = (double)malloc(sizeof(double) tr->cdta->endsite);
5411
5412	for(i = 0; i < tr->cdta->endsite; i++)
5413	weightSum += (double)(tr->cdta->aliaswgt[i]);
5414
5415	modOpt(tr, adef);
5416	printf("Model optimization, best Tree: %f\n", tr->likelihood);
5417	bestLH = tr->likelihood;
5418
5419
5420	evaluateGenericInitrav(tr, tr->start);
5421
5422	evaluateGenericVector(tr, tr->start, bestVector);
5423
5424	INFILE = fopen(bootStrapFileName, "r");
5425	while((ch = getc(INFILE)) != EOF)
5426	{
5427	if(ch == ';')
5428	numberOfTrees++;
5429	}
5430	rewind(INFILE);
5431
5432	printf("Found %d trees in File %s\n", numberOfTrees, bootStrapFileName);
5433
5434	for(i = 0; i < numberOfTrees; i++)
5435	{
5436	treeReadLen(INFILE, tr, adef);
5437	treeEvaluate(tr, 2);
5438	tr->start = tr->nodep[1];
5439
5440	evaluateGenericInitrav(tr, tr->start);
5441
5442	currentLH = tr->likelihood;
5443	if(currentLH > bestLH)
5444	{
5445	printf("Better tree found %d at %f\n", i, currentLH);
5446	/exit(1);/
5447	}
5448	/printf("Tree %d %f\n",i, tr->likelihood);/
5449
5450	evaluateGenericVector(tr, tr->start, otherVector);
5451
5452	{
5453	int j;
5454	double temp, wtemp, sum, sum2, sd;
5455
5456	sum = 0.0;
5457	sum2 = 0.0;
5458
5459	for (j = 0; j < tr->cdta->endsite; j++)
5460	{
5461	temp = bestVector[j] - otherVector[j];
5462	wtemp = tr->cdta->aliaswgt[j] * temp;
5463	sum += wtemp;
5464	sum2 += wtemp * temp;
5465	}
5466
5467	sd = sqrt( weightSum * (sum2 - sum*sum / weightSum)
5468	/ (weightSum - 1) );
5469
5470	printf("Tree: %d Likelihood: %f D(LH): %f SD: %f Significantly Worse: %s\n", i, currentLH, currentLH - bestLH, sd, (sum > 1.95996 * sd) ? "Yes" : " No");
5471	}
5472	}
5473
5474	fclose(INFILE);
5475
5476	free(bestVector);
5477	free(otherVector);
5478	exit(0);
5479	}
5480
5481	static void computePerSiteLLs(tree tr, analdef adef, char *bootStrapFileName)
5482	{
5483	int numberOfTrees = 0, i, j;
5484	char ch;
5485	double *otherVector;
5486	FILE *tlf;
5487
5488	tlf = fopen( perSiteLLsFileName, "w");
5489
5490	otherVector = (double)malloc(sizeof(double) tr->cdta->endsite);
5491
5492	allocNodex(tr, adef);
5493
5494	INFILE = fopen(bootStrapFileName, "r");
5495	while((ch = getc(INFILE)) != EOF)
5496	{
5497	if(ch == ';')
5498	numberOfTrees++;
5499	}
5500	rewind(INFILE);
5501
5502	printf("Found %d trees in File %s\n", numberOfTrees, bootStrapFileName);
5503
5504	fprintf(tlf, " %d %d\n", numberOfTrees, tr->cdta->endsite);
5505
5506	for(i = 0; i < numberOfTrees; i++)
5507	{
5508	treeReadLen(INFILE, tr, adef);
5509	if(i == 0)
5510	modOpt(tr, adef);
5511	else
5512	treeEvaluate(tr, 2);
5513
5514	printf("Tree %d: %f\n", i, tr->likelihood);
5515
5516	tr->start = tr->nodep[1];
5517
5518	evaluateGenericInitrav(tr, tr->start);
5519
5520	evaluateGenericVector(tr, tr->start, otherVector);
5521
5522	fprintf(tlf, "tr%d\t", i + 1);
5523	for(j = 0; j < tr->cdta->endsite; j++)
5524	{
5525	fprintf(tlf, "%f ", otherVector[j]);
5526	}
5527	fprintf(tlf, "\n");
5528	}
5529
5530	fclose(INFILE);
5531	fclose(tlf);
5532	free(otherVector);
5533	}
5534
5535	#ifdef _USE_PTHREADS
5536
5537	#ifndef _MAC
5538	#include <sched.h>
5539
5540	static void pinThread2Cpu(int myTid)
5541	{
5542	char *coreSteppingStr;
5543	int myCore, len;
5544	cpu_set_t cpuMask;
5545
5546	coreSteppingStr = getenv("SCHEDULE");
5547	len = coreSteppingStr ? strlen(coreSteppingStr) : 0;
5548
5549	myCore = myTid;
5550
5551	if (myTid < len)
5552	{
5553	if ((coreSteppingStr[myTid] >= '0') && (coreSteppingStr[myTid] <= '9'))
5554	myCore = coreSteppingStr[myTid] - '0';
5555
5556	if ((coreSteppingStr[myTid] >= 'a') && (coreSteppingStr[myTid] <= 'f'))
5557	myCore = coreSteppingStr[myTid] - 'a' + 10;
5558
5559	if ((coreSteppingStr[myTid] >= 'A') && (coreSteppingStr[myTid] <= 'F'))
5560	myCore = coreSteppingStr[myTid] - 'A' + 10;
5561	}
5562
5563	CPU_ZERO(&cpuMask);
5564	CPU_SET(myCore, &cpuMask);
5565
5566	if (sched_setaffinity(0, sizeof(cpuMask), &cpuMask))
5567	{
5568	printf("Error while scheduling Thread #%d to CPU %d\n",
5569	myTid, myCore);
5570	exit(1);
5571	}
5572
5573	/* printf("Scheduled Thread #%d to logical CPU %d\n", myTid, myCore); */
5574	return;
5575	}
5576
5577	#endif
5578
5579	typedef struct {
5580	tree *tr;
5581	int threadNumber;
5582	} threadData;
5583
5584
5585	static void calcBounds(int tid, const int n, int start, int end, int l, int u)
5586	{
5587	int span = end - start;
5588
5589	/* LTD */
5590	/* assert(span % n == 0); */
5591
5592	if(span % n == 0)
5593	span = span / n;
5594	else
5595	span = 1 + (span / n);
5596
5597	l = start + tid span;
5598	if(tid == n - 1)
5599	*u = end;
5600	else
5601	u = l + span;
5602	}
5603
5604	#ifdef _LOCAL_DATA
5605
5606	static void strided_Bounds(int tid, int endsite, int n, int startIndex, int endIndex)
5607	{
5608	int endsiteL = endsite - tid;
5609
5610	if(endsiteL % n == 0)
5611	endsiteL = endsiteL / n;
5612	else
5613	endsiteL = 1 + (endsiteL / n);
5614
5615	*startIndex = 0;
5616	*endIndex = endsiteL;
5617	}
5618
5619	static void collectDouble(double dest, double source, const int totallength, const int stride, const int offset)
5620	{
5621	int
5622	i = 0,
5623	k = offset;
5624
5625	for(; k < totallength; i++, k += stride)
5626	dest[k] = source[i];
5627	}
5628
5629
5630
5631
5632	static void strideTips(char dest, char source, const int totallength, const int stride,
5633	const int offset, const int mxtips, int strideLength)
5634	{
5635	int i, j, k;
5636
5637	assert(offset < stride);
5638
5639	for(i = 0; i < mxtips; i++)
5640	{
5641	char *d = &dest[i + 1][0];
5642	char *s = &source[i + 1][0];
5643
5644	for(k = 0, j = offset; j < totallength; j += stride, k++)
5645	{
5646	assert(k < strideLength);
5647	d[k] = s[j];
5648	}
5649	}
5650
5651	}
5652
5653	static void strideInt(int dest, int source, const int totallength, const int stride, const int offset)
5654	{
5655	int i, k,
5656	*d = &dest[0];
5657
5658	for(i = offset, k = 0; i < totallength; i += stride, k++)
5659	d[k] = source[i];
5660	}
5661
5662	static void strideDouble(double dest, double source, const int totallength, const int stride, const int offset)
5663	{
5664	int i = offset;
5665	double *d = dest;
5666
5667	for(; i < totallength; i += stride, d++)
5668	*d = source[i];
5669	}
5670
5671
5672	static void stridePartitionData(tree *localTree, int tid, int n, int length)
5673	{
5674	int
5675	i,
5676	endsite,
5677	dummy;
5678
5679	strided_Bounds(tid, length, n, &dummy, &endsite);
5680
5681	/* printf("%d %d\n", endsite, tid); */
5682
5683	for(i = 0; i < localTree->NumberOfModels; i++)
5684	{
5685	localTree->strided_partitionData[i].dataType = localTree->partitionData[i].dataType;
5686	localTree->strided_partitionData[i].protModels = localTree->strided_partitionData[i].protModels;
5687	localTree->strided_partitionData[i].protFreqs = localTree->strided_partitionData[i].protFreqs;
5688	}
5689
5690	if(localTree->NumberOfModels > 1)
5691	{
5692	int i, model;
5693
5694	localTree->strided_partitionData[0].lower = 0;
5695
5696	model = localTree->strided_model[0];
5697	i = 1;
5698
5699	while(i < endsite)
5700	{
5701	if(localTree->strided_model[i] != model)
5702	{
5703	localTree->strided_partitionData[model].upper = i;
5704	localTree->strided_partitionData[model + 1].lower = i;
5705	model = localTree->strided_model[i];
5706	}
5707	i++;
5708	}
5709
5710
5711	localTree->strided_partitionData[localTree->NumberOfModels - 1].upper = endsite;
5712
5713	/*
5714	for(i = 0; i < localTree->NumberOfModels; i++)
5715	printf("%d %d %d\n", tid, localTree->strided_partitionData[i].lower, localTree->strided_partitionData[i].upper);
5716	*/
5717	}
5718	else
5719	{
5720	localTree->strided_partitionData[0].lower = 0;
5721	localTree->strided_partitionData[0].upper = endsite;
5722	}
5723
5724	}
5725
5726
5727	inline static void sendTraversalInfo(tree localTree, tree tr)
5728	{
5729	/* the one below is a hack we are re-assigning the local pointer to the global one
5730	the memcpy version below is just for testing and preparing the
5731	fine-grained MPI BlueGene version */
5732
5733	if(1)
5734	{
5735	localTree->td[0] = tr->td[0];
5736	}
5737	else
5738	{
5739	localTree->td[0].count = tr->td[0].count;
5740	memcpy(localTree->td[0].ti, tr->td[0].ti, localTree->td[0].count * sizeof(traversalInfo));
5741	}
5742	}
5743
5744	#endif
5745
5746
5747
5748
5749	static void execFunction(tree tr, tree localTree, const int startIndex, const int endIndex,
5750	const int parsimonyStartIndex, const int parsimonyEndIndex, int tid, int n)
5751	{
5752	double result, dlnLdlz, d2lnLdlz2;
5753	int parsimonyResult;
5754
5755	/* new */
5756	int currentJob;
5757
5758	currentJob = threadJob >> 16;
5759
5760	/* new */
5761	switch(currentJob)
5762	/switch(threadJob) /
5763	{
5764	case THREAD_NEWVIEW:
5765	#ifdef _LOCAL_DATA
5766	/* send */
5767	sendTraversalInfo(localTree, tr);
5768
5769	newviewIterative(localTree, startIndex, endIndex);
5770	#else
5771	newviewIterative(tr, startIndex, endIndex);
5772	#endif
5773	break;
5774
5775	/*****************************************************/
5776
5777	case THREAD_EVALUATE:
5778	#ifdef _LOCAL_DATA
5779	/* send */
5780	sendTraversalInfo(localTree, tr);
5781	result = evaluateIterative(localTree, startIndex, endIndex);
5782	#else
5783	result = evaluateIterative(tr, startIndex, endIndex);
5784	#endif
5785
5786	/* receive */
5787	reductionBuffer[tid] = result;
5788	break;
5789
5790	/*****************************************************/
5791
5792	case THREAD_SUM_MAKENEWZ:
5793	#ifdef _LOCAL_DATA
5794	/* send */
5795	sendTraversalInfo(localTree, tr);
5796	makenewzIterative(localTree, startIndex, endIndex);
5797	#else
5798	makenewzIterative(tr, startIndex, endIndex);
5799	#endif
5800
5801	break;
5802
5803	/*****************************************************/
5804
5805	case THREAD_MAKENEWZ:
5806	#ifdef _LOCAL_DATA
5807
5808	/* send */
5809
5810	localTree->modelNumber = tr->modelNumber;
5811	localTree->coreLZ = tr->coreLZ;
5812
5813	if(localTree->multiBranch)
5814	execCore(localTree, &dlnLdlz, &d2lnLdlz2, localTree->strided_partitionData[localTree->modelNumber].lower,
5815	localTree->strided_partitionData[localTree->modelNumber].upper, localTree->modelNumber);
5816	else
5817	execCore(localTree, &dlnLdlz, &d2lnLdlz2, startIndex, endIndex, localTree->modelNumber);
5818	#else
5819	if(tr->multiBranch)
5820	{
5821	int u, l;
5822	int start = tr->partitionData[tr->modelNumber].lower;
5823	int end = tr->partitionData[tr->modelNumber].upper;
5824
5825	calcBounds(tid, n, start, end, &l, &u);
5826
5827
5828	execCore(tr, &dlnLdlz, &d2lnLdlz2, l, u, tr->modelNumber);
5829	}
5830	else
5831	{
5832	execCore(tr, &dlnLdlz, &d2lnLdlz2, startIndex, endIndex, tr->modelNumber);
5833	}
5834	#endif
5835
5836
5837	/* receive */
5838	reductionBuffer[tid] = dlnLdlz;
5839	reductionBufferTwo[tid] = d2lnLdlz2;
5840	break;
5841
5842	/*****************************************************/
5843
5844	case THREAD_SUM_MAKENEWZ_PARTITION:
5845	{
5846	int u, l, start, end;
5847
5848	#ifdef _LOCAL_DATA
5849	/* TODO */
5850	assert(0);
5851	/* only required for a rarely used, undocumented function */
5852	#endif
5853
5854	start = tr->partitionData[tr->modelNumber].lower;
5855	end = tr->partitionData[tr->modelNumber].upper;
5856	calcBounds(tid, n, start, end, &l, &u);
5857
5858	makenewzIterativePartition(tr, l, u, tr->modelNumber);
5859	}
5860	break;
5861
5862	/*****************************************************/
5863
5864	case THREAD_MAKENEWZ_PARTITION:
5865	{
5866	int u, l, start, end;
5867
5868	#ifdef _LOCAL_DATA
5869	/* TODO */
5870	assert(0);
5871	/* only required for a rarely used, undocumented function */
5872	#endif
5873
5874
5875	start = tr->partitionData[tr->modelNumber].lower;
5876	end = tr->partitionData[tr->modelNumber].upper;
5877	calcBounds(tid, n, start, end, &l, &u);
5878
5879	execCorePartition(tr, &dlnLdlz, &d2lnLdlz2, l, u, tr->modelNumber);
5880	reductionBuffer[tid] = dlnLdlz;
5881	reductionBufferTwo[tid] = d2lnLdlz2;
5882	}
5883	break;
5884
5885	/*****************************************************/
5886
5887	case THREAD_NEWVIEW_PARTITION:
5888	#ifdef _LOCAL_DATA
5889	/* send */
5890	localTree->modelNumber = tr->modelNumber;
5891	sendTraversalInfo(localTree, tr);
5892
5893	newviewIterativePartition(localTree, localTree->strided_partitionData[localTree->modelNumber].lower,
5894	localTree->strided_partitionData[localTree->modelNumber].upper, localTree->modelNumber);
5895	#else
5896	{
5897	int u, l;
5898	int start = tr->partitionData[tr->modelNumber].lower;
5899	int end = tr->partitionData[tr->modelNumber].upper;
5900	calcBounds(tid, n, start, end, &l, &u);
5901
5902	newviewIterativePartition(tr, l, u, tr->modelNumber);
5903	}
5904	#endif
5905	break;
5906
5907	/*****************************************************/
5908
5909	case THREAD_EVALUATE_PARTITION:
5910	#ifdef _LOCAL_DATA
5911	/* send */
5912
5913	localTree->modelNumber = tr->modelNumber;
5914	sendTraversalInfo(localTree, tr);
5915
5916
5917
5918
5919	result = evaluateIterativePartition(localTree, localTree->strided_partitionData[localTree->modelNumber].lower,
5920	localTree->strided_partitionData[localTree->modelNumber].upper,
5921	localTree->modelNumber);
5922	#else
5923	{
5924	int u, l;
5925	int start = tr->partitionData[tr->modelNumber].lower;
5926	int end = tr->partitionData[tr->modelNumber].upper;
5927	calcBounds(tid, n, start, end, &l, &u);
5928	result = evaluateIterativePartition(tr, l, u, tr->modelNumber);
5929	}
5930	#endif
5931
5932	/* receive */
5933	reductionBuffer[tid] = result;
5934	break;
5935
5936	/*****************************************************/
5937
5938	case THREAD_RATE_CATS:
5939	#ifdef _LOCAL_DATA
5940
5941	/* send */
5942
5943	sendTraversalInfo(localTree, tr);
5944	localTree->lower_spacing = tr->lower_spacing;
5945	localTree->upper_spacing = tr->upper_spacing;
5946
5947	optRateCat_LOCAL(localTree, startIndex, endIndex,
5948	localTree->lower_spacing, localTree->upper_spacing, localTree->strided_lhs);
5949
5950	/* receive */
5951
5952	collectDouble(tr->cdta->patrat, localTree->strided_patrat, tr->cdta->endsite, n, tid);
5953	collectDouble(tr->cdta->patratStored, localTree->strided_patratStored, tr->cdta->endsite, n, tid);
5954	collectDouble(tr->lhs, localTree->strided_lhs, tr->cdta->endsite, n, tid);
5955	#else
5956	{
5957	int i;
5958	for(i = 0; i < tr->cdta->endsite; i++)
5959	if(i % n == tid)
5960	optRateCat(tr, i, tr->lower_spacing, tr->upper_spacing, tr->lhs);
5961	}
5962	#endif
5963	break;
5964
5965	/*****************************************************/
5966
5967	case THREAD_NEWVIEW_PARSIMONY:
5968
5969	#ifdef _LOCAL_DATA
5970	/* send */
5971	sendTraversalInfo(localTree, tr);
5972	newviewParsimonyIterative(localTree, parsimonyStartIndex, parsimonyEndIndex);
5973	#else
5974	newviewParsimonyIterative(tr, parsimonyStartIndex, parsimonyEndIndex);
5975	#endif
5976
5977	break;
5978
5979	/*****************************************************/
5980
5981	case THREAD_EVALUATE_PARSIMONY:
5982
5983	#ifdef _LOCAL_DATA
5984	/* send */
5985	sendTraversalInfo(localTree, tr);
5986	parsimonyResult = evaluateParsimonyIterative(localTree, parsimonyStartIndex, parsimonyEndIndex);
5987	#else
5988	parsimonyResult = evaluateParsimonyIterative(tr, parsimonyStartIndex, parsimonyEndIndex);
5989	#endif
5990
5991	/* receive */
5992	reductionBufferParsimony[tid] = parsimonyResult;
5993
5994	break;
5995
5996	/*****************************************************/
5997
5998	case THREAD_EVALUATE_VECTOR:
5999	#ifdef _LOCAL_DATA
6000	sendTraversalInfo(localTree, tr);
6001	evaluateGenericVectorIterative(localTree, startIndex, endIndex);
6002
6003	collectDouble(tr->siteLL_Vector, localTree->strided_siteLL_Vector, tr->cdta->endsite, n, tid);
6004	/* TODO */
6005	/* assert(0);*/
6006	/* rarely used function */
6007	#else
6008
6009	evaluateGenericVectorIterative(tr, startIndex, endIndex);
6010	#endif
6011	break;
6012
6013	/*****************************************************/
6014
6015	case THREAD_CATEGORIZE:
6016	#ifdef _LOCAL_DATA
6017	{
6018	int i;
6019
6020	/* send */
6021	sendTraversalInfo(localTree, tr);
6022
6023	for(i = 0; i < localTree->NumberOfModels; i++)
6024	{
6025	localTree->strided_patrat[i * 4] = localTree->gammaRates[i * 4];
6026	localTree->strided_patrat[i * 4 + 1] = localTree->gammaRates[i * 4 + 1];
6027	localTree->strided_patrat[i * 4 + 2] = localTree->gammaRates[i * 4 + 2];
6028	localTree->strided_patrat[i * 4 + 3] = localTree->gammaRates[i * 4 + 3];
6029	assert(i * 4 + 3 < localTree->originalCrunchedLength);
6030	}
6031
6032	localTree->NumberOfCategories = 4 * localTree->NumberOfModels;
6033	categorizeIterative(localTree, startIndex, endIndex);
6034
6035	for(i = startIndex; i < endIndex; i++)
6036	{
6037	double temp, wtemp;
6038	temp = localTree->gammaRates[localTree->strided_rateCategory[i]];
6039	localTree->strided_wr[i] = wtemp = temp * localTree->strided_aliaswgt[i];
6040	localTree->strided_wr2[i] = temp * wtemp;
6041	}
6042	}
6043	#else
6044	categorizeIterative(tr, startIndex, endIndex);
6045	#endif
6046	break;
6047
6048	/*****************************************************/
6049
6050
6051	#ifdef _LOCAL_DATA
6052	case THREAD_PREPARE_PARSIMONY:
6053	/printf("THREAD_PREPARE_PARSIMONY\n"); /
6054	if(tid > 0)
6055	{
6056	localTree->parsimonyLength = tr->parsimonyLength;
6057	memcpy(localTree->partitionData , tr->partitionData, sizeof(pInfo) * localTree->NumberOfModels);
6058	}
6059
6060	strideInt(localTree->strided_model, tr->model,
6061	localTree->originalCrunchedLength, n, tid);
6062	strideInt(localTree->strided_dataVector, tr->dataVector,
6063	localTree->originalCrunchedLength, n, tid);
6064	stridePartitionData(localTree, tid, n, localTree->parsimonyLength);
6065
6066	strideTips(localTree->strided_yVector, tr->yVector, localTree->originalCrunchedLength, n, tid,
6067	localTree->mxtips, localTree->strideLength);
6068	strideInt(localTree->strided_aliaswgt, tr->cdta->aliaswgt,
6069	localTree->originalCrunchedLength, n, tid);
6070
6071	localTree->mySpan = 1 + (localTree->parsimonyLength / n);
6072	localTree->parsimonyData = (parsimonyVector )malloc(sizeof(parsimonyVector)
6073	localTree->mxtips * localTree->mySpan);
6074
6075	break;
6076
6077	/*****************************************************/
6078
6079	case THREAD_FINISH_PARSIMONY:
6080	free(localTree->parsimonyData);
6081
6082	if(tid > 0)
6083	{
6084	localTree->cdta->endsite = tr->cdta->endsite;
6085	memcpy(localTree->partitionData , tr->partitionData, sizeof(pInfo) * localTree->NumberOfModels);
6086	}
6087	strideInt(localTree->strided_model, tr->model,
6088	localTree->originalCrunchedLength, n, tid);
6089	strideInt(localTree->strided_dataVector, tr->dataVector,
6090	localTree->originalCrunchedLength, n, tid);
6091	stridePartitionData(localTree, tid, n, localTree->cdta->endsite);
6092
6093	strideTips(localTree->strided_yVector, tr->yVector, localTree->originalCrunchedLength, n, tid,
6094	localTree->mxtips, localTree->strideLength);
6095	strideInt(localTree->strided_aliaswgt, tr->cdta->aliaswgt,
6096	localTree->originalCrunchedLength, n, tid);
6097
6098	break;
6099
6100	/*****************************************************/
6101
6102	case THREAD_ALLOC_LIKELIHOOD:
6103	/printf("THREAD_ALLOC_LIKELIHOOD\n");/
6104	{
6105	int span, i;
6106
6107	localTree->likelihoodFunction = tr->likelihoodFunction;
6108	localTree->currentModel = tr->currentModel;
6109	localTree->cdta->endsite = tr->cdta->endsite;
6110	localTree->mySpan = 1 + (localTree->cdta->endsite / n);
6111
6112	localTree->expArray = (int )malloc(localTree->mySpan localTree->mxtips * sizeof(int));
6113
6114	if(localTree->mixedData)
6115	{
6116	assert(0);
6117	}
6118	else
6119	{
6120	switch(localTree->currentModel)
6121	{
6122	case M_PROTCAT:
6123	span = 20 * localTree->mySpan;
6124	localTree->sumBuffer = (double )malloc(20 localTree->strideLength * sizeof(double));
6125	break;
6126	case M_PROTGAMMA:
6127	span = 80 * localTree->mySpan;
6128	localTree->sumBuffer = (double )malloc(80 localTree->strideLength * sizeof(double));
6129	break;
6130	case M_GTRGAMMA:
6131	span = 16 * localTree->mySpan;
6132	localTree->sumBuffer = (double )malloc(16 localTree->strideLength * sizeof(double));
6133	break;
6134	case M_GTRCAT:
6135	span = 4 * localTree->mySpan;
6136	localTree->sumBuffer = (double )malloc(4 localTree->strideLength * sizeof(double));
6137	break;
6138	default:
6139	assert(0);
6140	}
6141	localTree->likelihoodArray = (double )malloc(span localTree->mxtips * sizeof(double));
6142	}
6143
6144
6145	for(i = 0; i < localTree->mxtips; i++)
6146	localTree->xVector[i] = &(localTree->likelihoodArray[i * span]);
6147	}
6148	break;
6149
6150	/*****************************************************/
6151
6152	case THREAD_FREE_LIKELIHOOD:
6153	free(localTree->expArray);
6154	free(localTree->likelihoodArray);
6155	free(localTree->sumBuffer);
6156	localTree->expArray = (int*)NULL;
6157	localTree->likelihoodArray = (double*)NULL;
6158	localTree->sumBuffer = (double*)NULL;
6159	break;
6160
6161	/*****************************************************/
6162
6163	case THREAD_COPY_REVERSIBLE:
6164	if(tid > 0)
6165	{
6166	memcpy(localTree->tipVectorDNA, tr->tipVectorDNA, localTree->NumberOfModels * 64 * sizeof(double));
6167	memcpy(localTree->tipVectorAA, tr->tipVectorAA, localTree->NumberOfModels * 460 * sizeof(double));
6168
6169	memcpy(localTree->EV_DNA, tr->EV_DNA, localTree->NumberOfModels * 16 * sizeof(double));
6170	memcpy(localTree->EV_AA, tr->EV_AA, localTree->NumberOfModels * 400 * sizeof(double));
6171
6172	memcpy(localTree->EI_DNA, tr->EI_DNA, localTree->NumberOfModels * 12 * sizeof(double));
6173	memcpy(localTree->EI_AA, tr->EI_AA, localTree->NumberOfModels * 380 * sizeof(double));
6174
6175	memcpy(localTree->EIGN_DNA, tr->EIGN_DNA, localTree->NumberOfModels * 3 * sizeof(double));
6176	memcpy(localTree->EIGN_AA, tr->EIGN_AA, localTree->NumberOfModels * 19 * sizeof(double));
6177	}
6178	break;
6179
6180	/*****************************************************/
6181
6182	case THREAD_COPY_RATE_CATS:
6183	if(tid > 0)
6184	localTree->NumberOfCategories = tr->NumberOfCategories;
6185
6186	strideInt(localTree->strided_rateCategory, tr->cdta->rateCategory,
6187	localTree->originalCrunchedLength, n, tid);
6188
6189	memcpy(localTree->strided_patrat, tr->cdta->patrat, localTree->originalCrunchedLength * sizeof(double));
6190
6191	strideDouble(localTree->strided_patratStored, tr->cdta->patratStored,
6192	localTree->originalCrunchedLength, n, tid);
6193	strideDouble(localTree->strided_wr, tr->cdta->wr,
6194	localTree->originalCrunchedLength, n, tid);
6195	strideDouble(localTree->strided_wr2, tr->cdta->wr2,
6196	localTree->originalCrunchedLength, n, tid);
6197
6198	break;
6199
6200	/*****************************************************/
6201
6202	case THREAD_COPY_GAMMA_RATES:
6203	if(tid > 0)
6204	memcpy(localTree->gammaRates, tr->gammaRates, localTree->NumberOfModels * 4 * sizeof(double));
6205	break;
6206
6207	/*****************************************************/
6208
6209	case THREAD_COPY_INVARIANTS:
6210	if(tid > 0)
6211	memcpy(localTree->invariants, tr->invariants, localTree->NumberOfModels * sizeof(double));
6212	break;
6213
6214	/*****************************************************/
6215
6216	case THREAD_COPY_INIT_MODEL:
6217	/* printf("THREAD_COPY_INIT_MODEL\n"); */
6218	if(tid > 0)
6219	{
6220	localTree->NumberOfCategories = tr->NumberOfCategories;
6221	localTree->likelihoodFunction = tr->likelihoodFunction;
6222	localTree->cdta->endsite = tr->cdta->endsite;
6223
6224	memcpy(localTree->tipVectorDNA, tr->tipVectorDNA, localTree->NumberOfModels * 64 * sizeof(double));
6225	memcpy(localTree->tipVectorAA, tr->tipVectorAA, localTree->NumberOfModels * 460 * sizeof(double));
6226
6227	memcpy(localTree->EV_DNA, tr->EV_DNA, localTree->NumberOfModels * 16 * sizeof(double));
6228	memcpy(localTree->EV_AA, tr->EV_AA, localTree->NumberOfModels * 400 * sizeof(double));
6229
6230	memcpy(localTree->EI_DNA, tr->EI_DNA, localTree->NumberOfModels * 12 * sizeof(double));
6231	memcpy(localTree->EI_AA, tr->EI_AA, localTree->NumberOfModels * 380 * sizeof(double));
6232
6233	memcpy(localTree->EIGN_DNA, tr->EIGN_DNA, localTree->NumberOfModels * 3 * sizeof(double));
6234	memcpy(localTree->EIGN_AA, tr->EIGN_AA, localTree->NumberOfModels * 19 * sizeof(double));
6235
6236	memcpy(localTree->frequencies_DNA, tr->frequencies_DNA, localTree->NumberOfModels * 4 * sizeof(double));
6237	memcpy(localTree->frequencies_AA, tr->frequencies_AA, localTree->NumberOfModels * 20 * sizeof(double));
6238
6239	memcpy(localTree->invariants, tr->invariants, localTree->NumberOfModels * sizeof(double));
6240
6241	memcpy(localTree->gammaRates, tr->gammaRates, localTree->NumberOfModels * 4 * sizeof(double));
6242
6243	memcpy(localTree->partitionData , tr->partitionData, sizeof(pInfo) * localTree->NumberOfModels);
6244	}
6245
6246	strideInt(localTree->strided_model, tr->model,
6247	localTree->originalCrunchedLength, n, tid);
6248	strideInt(localTree->strided_dataVector, tr->dataVector,
6249	localTree->originalCrunchedLength, n, tid);
6250	stridePartitionData(localTree, tid, n, localTree->cdta->endsite);
6251
6252	strideInt(localTree->strided_rateCategory, tr->cdta->rateCategory,
6253	localTree->originalCrunchedLength, n, tid);
6254
6255	strideInt(localTree->strided_aliaswgt, tr->cdta->aliaswgt,
6256	localTree->originalCrunchedLength, n, tid);
6257
6258
6259	strideInt(localTree->strided_invariant, tr->invariant,
6260	localTree->originalCrunchedLength, n, tid);
6261
6262	memcpy(localTree->strided_patrat, tr->cdta->patrat, localTree->originalCrunchedLength * sizeof(double));
6263
6264	strideDouble(localTree->strided_patratStored, tr->cdta->patratStored,
6265	localTree->originalCrunchedLength, n, tid);
6266	strideDouble(localTree->strided_wr, tr->cdta->wr,
6267	localTree->originalCrunchedLength, n, tid);
6268	strideDouble(localTree->strided_wr2, tr->cdta->wr2,
6269	localTree->originalCrunchedLength, n, tid);
6270
6271	strideTips(localTree->strided_yVector, tr->yVector, localTree->originalCrunchedLength, n, tid,
6272	localTree->mxtips, localTree->strideLength);
6273	break;
6274
6275	/*****************************************************/
6276
6277	case THREAD_NEXT_REPLICATE:
6278	/* printf("THREAD_NEXT_REPLICATE\n"); */
6279	if(tid > 0)
6280	{
6281	localTree->cdta->endsite = tr->cdta->endsite;
6282	memcpy(localTree->partitionData , tr->partitionData, sizeof(pInfo) * localTree->NumberOfModels);
6283	}
6284
6285	strideInt(localTree->strided_model, tr->model,
6286	localTree->originalCrunchedLength, n, tid);
6287	strideInt(localTree->strided_dataVector, tr->dataVector,
6288	localTree->originalCrunchedLength, n, tid);
6289	stridePartitionData(localTree, tid, n, localTree->cdta->endsite);
6290
6291	strideInt(localTree->strided_aliaswgt, tr->cdta->aliaswgt,
6292	localTree->originalCrunchedLength, n, tid);
6293	strideInt(localTree->strided_rateCategory, tr->cdta->rateCategory,
6294	localTree->originalCrunchedLength, n, tid);
6295
6296	strideDouble(localTree->strided_wr, tr->cdta->wr,
6297	localTree->originalCrunchedLength, n, tid);
6298	strideDouble(localTree->strided_wr2, tr->cdta->wr2,
6299	localTree->originalCrunchedLength, n, tid);
6300
6301	memcpy(localTree->strided_patrat, tr->cdta->patrat, localTree->originalCrunchedLength * sizeof(double));
6302
6303	strideTips(localTree->strided_yVector, tr->yVector, localTree->originalCrunchedLength, n, tid,
6304	localTree->mxtips, localTree->strideLength);
6305
6306	break;
6307	#endif
6308	default:
6309	assert(0);
6310	}
6311	}
6312
6313
6314
6315	void masterBarrier(int jobType, tree *tr)
6316	{
6317	const int n = NumberOfThreads;
6318	int startIndex, endIndex, i, sum,
6319	parsimonyStartIndex, parsimonyEndIndex;
6320
6321	/* new */
6322	jobCycle = !jobCycle;
6323	threadJob = (jobType << 16) + jobCycle;
6324
6325	/*
6326	old
6327	threadJob = jobType;
6328	jobCycle = !jobCycle;
6329	*/
6330
6331	#ifdef _LOCAL_DATA
6332	strided_Bounds(0, tr->cdta->endsite, n, &startIndex, &endIndex);
6333	strided_Bounds(0, tr->parsimonyLength, n, &parsimonyStartIndex, &parsimonyEndIndex);
6334	#else
6335	calcBounds(0, n, 0, tr->parsimonyLength, &parsimonyStartIndex, &parsimonyEndIndex);
6336	calcBounds(0, n, 0, tr->cdta->endsite, &startIndex, &endIndex);
6337	#endif
6338
6339	execFunction(tr, tr, startIndex, endIndex, parsimonyStartIndex, parsimonyEndIndex, 0, n);
6340
6341	do
6342	{
6343	for(i = 1, sum = 1; i < n; i++)
6344	sum += barrierBuffer[i];
6345	}
6346	while(sum < n);
6347
6348	for(i = 1; i < n; i++)
6349	barrierBuffer[i] = 0;
6350	/threadJob = -1; /
6351	}
6352
6353
6354	#ifdef _LOCAL_DATA
6355
6356	static void allocStrides(tree *tr)
6357	{
6358	int i;
6359
6360	if(tr->numBranches < NUM_BRANCHES)
6361	{
6362	printf("PERFORMANCE WARNING: for optimal efficiency on this dataset\n");
6363	printf("set NUM_BRANCHES to %d in file axml.h an re-compile\n", tr->numBranches);
6364	}
6365
6366	tr->strideLength = 1 + (tr->originalCrunchedLength / NumberOfThreads);
6367
6368	tr->strided_y0 = (char )malloc(tr->strideLength tr->mxtips * sizeof(char));
6369	tr->strided_yVector = (char *)malloc((tr->mxtips + 1) sizeof(char *));
6370
6371	for(i = 0; i < tr->mxtips; i++)
6372	tr->strided_yVector[i + 1] = &(tr->strided_y0[tr->strideLength * i]);
6373
6374	tr->strided_aliaswgt = (int )malloc(sizeof(int) tr->strideLength);
6375	tr->strided_invariant = (int )malloc(sizeof(int) tr->strideLength);
6376	tr->strided_model = (int )malloc(sizeof(int) tr->strideLength);
6377	tr->strided_rateCategory = (int )malloc(sizeof(int) tr->strideLength);
6378	tr->strided_dataVector = (int )malloc(sizeof(int) tr->strideLength);
6379
6380	tr->strided_wr = (double )malloc(sizeof(double) tr->strideLength);
6381	tr->strided_wr2 = (double )malloc(sizeof(double) tr->strideLength);
6382	tr->strided_siteLL_Vector = (double )malloc(sizeof(double) tr->strideLength);
6383
6384	/* this is a bit ugly here */
6385
6386	tr->strided_patrat = (double )malloc(sizeof(double) tr->originalCrunchedLength);
6387
6388
6389
6390	tr->strided_patratStored = (double )malloc(sizeof(double) tr->strideLength);
6391	tr->strided_lhs = (double )malloc(sizeof(double) tr->strideLength);
6392
6393	tr->strided_partitionData = (pInfo)malloc(sizeof(pInfo) tr->NumberOfModels);
6394	}
6395
6396	#endif
6397
6398	static void likelihoodThread(void tData)
6399	{
6400	threadData td = (threadData)tData;
6401	tree *tr = td->tr;
6402	tree localTree = (tree )malloc(sizeof(tree));
6403
6404	int
6405	parsimonyStartIndex,
6406	parsimonyEndIndex,
6407	startIndex,
6408	endIndex,
6409	myCycle = 0;
6410
6411	const int n = NumberOfThreads;
6412	const int tid = td->threadNumber;
6413
6414	#ifdef _LOCAL_DATA
6415	cruncheddata cdta = (cruncheddata )malloc(sizeof(cruncheddata));
6416	#endif
6417
6418	#ifndef _MAC
6419	pinThread2Cpu(tid);
6420	#endif
6421
6422	#ifdef _LOCAL_DATA
6423	localTree->expArray = (int*)NULL;
6424	localTree->likelihoodArray = (double*)NULL;
6425	localTree->sumBuffer = (double*)NULL;
6426	localTree->cdta = cdta;
6427	localTree->mixedData = tr->mixedData;
6428	localTree->NumberOfModels = tr->NumberOfModels;
6429	localTree->mxtips = tr->mxtips;
6430	localTree->originalCrunchedLength = tr->originalCrunchedLength;
6431	localTree->multiBranch = tr->multiBranch;
6432	localTree->numBranches = tr->numBranches;
6433
6434	localTree->tipVectorDNA = (double )malloc(localTree->NumberOfModels 64 * sizeof(double));
6435	localTree->tipVectorAA = (double )malloc(localTree->NumberOfModels 460 * sizeof(double));
6436
6437	localTree->EV_DNA = (double )malloc(localTree->NumberOfModels 16 * sizeof(double));
6438	localTree->EV_AA = (double )malloc(localTree->NumberOfModels 400 * sizeof(double));
6439
6440	localTree->EI_DNA = (double )malloc(localTree->NumberOfModels 12 * sizeof(double));
6441	localTree->EI_AA = (double )malloc(localTree->NumberOfModels 380 * sizeof(double));
6442
6443	localTree->EIGN_DNA = (double )malloc(localTree->NumberOfModels 3 * sizeof(double));
6444	localTree->EIGN_AA = (double )malloc(localTree->NumberOfModels 19 * sizeof(double));
6445
6446	localTree->frequencies_DNA = (double )malloc(localTree->NumberOfModels 4 * sizeof(double));
6447	localTree->frequencies_AA = (double )malloc(localTree->NumberOfModels 20 * sizeof(double));
6448
6449	localTree->initialRates_DNA = (double )malloc(localTree->NumberOfModels 5 * sizeof(double));
6450	localTree->initialRates_AA = (double )malloc(localTree->NumberOfModels 190 * sizeof(double));
6451
6452	localTree->xVector = (double *)malloc(localTree->mxtips sizeof(double *));
6453
6454	localTree->gammaRates = (double )malloc(localTree->NumberOfModels 4 * sizeof(double));
6455	localTree->invariants = (double )malloc(localTree->NumberOfModels sizeof(double));
6456	localTree->model = (int ) malloc(localTree->originalCrunchedLength sizeof(int));
6457
6458	localTree->partitionData = (pInfo)malloc(sizeof(pInfo) localTree->NumberOfModels);
6459
6460	localTree->td[0].count = 0;
6461	localTree->td[0].ti = (traversalInfo )malloc(sizeof(traversalInfo) localTree->mxtips);
6462
6463	localTree->NumberOfCategories = tr->NumberOfCategories;
6464
6465	allocStrides(localTree);
6466
6467
6468	#endif
6469
6470	printf("\nThis is RAxML Worker Pthread Number: %d\n", tid);
6471
6472	while(1)
6473	{
6474	/*
6475	old
6476	while (myCycle == jobCycle);
6477	myCycle = !myCycle;
6478	*/
6479
6480	/* new */
6481	while (myCycle == threadJob);
6482	myCycle = threadJob;
6483
6484	#ifdef _LOCAL_DATA
6485	strided_Bounds(tid, localTree->cdta->endsite, n, &startIndex, &endIndex);
6486	strided_Bounds(tid, localTree->parsimonyLength, n, &parsimonyStartIndex, &parsimonyEndIndex);
6487	#else
6488	calcBounds(tid, n, 0, tr->cdta->endsite, &startIndex, &endIndex);
6489	calcBounds(tid, n, 0, tr->parsimonyLength, &parsimonyStartIndex, &parsimonyEndIndex);
6490	#endif
6491
6492	execFunction(tr, localTree, startIndex, endIndex, parsimonyStartIndex, parsimonyEndIndex, tid, n);
6493
6494	barrierBuffer[tid] = 1;
6495	}
6496
6497	return (void*)NULL;
6498	}
6499
6500	static void startPthreads(tree *tr)
6501	{
6502	pthread_t *threads;
6503	pthread_attr_t attr;
6504	int rc, t;
6505	threadData *tData;
6506
6507
6508	/* TODO pthread_attr_getstackaddr and pthread_attr_setstackaddr */
6509
6510	jobCycle = 0;
6511	/* old */
6512	/* threadJob = -1; */
6513	/* new */
6514	threadJob = 0;
6515
6516	printf("\nThis is the RAxML Master Pthread\n");
6517
6518	pthread_attr_init(&attr);
6519	pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_DETACHED);
6520
6521	threads = (pthread_t )malloc(NumberOfThreads sizeof(pthread_t));
6522	tData = (threadData )malloc(NumberOfThreads sizeof(threadData));
6523	reductionBuffer = (double )malloc(sizeof(double) NumberOfThreads);
6524	reductionBufferTwo = (double )malloc(sizeof(double) NumberOfThreads);
6525	reductionBufferParsimony = (int )malloc(sizeof(int) NumberOfThreads);
6526	barrierBuffer = (int )malloc(sizeof(int) NumberOfThreads);
6527
6528	#ifdef _LOCAL_DATA
6529	allocStrides(tr);
6530	#endif
6531
6532	for(t = 0; t < NumberOfThreads; t++)
6533	barrierBuffer[t] = 0;
6534
6535	#ifndef _MAC
6536	pinThread2Cpu(0);
6537	#endif
6538
6539	for(t = 1; t < NumberOfThreads; t++)
6540	{
6541	tData[t].tr = tr;
6542	tData[t].threadNumber = t;
6543	rc = pthread_create(&threads[t], &attr, likelihoodThread, (void *)(&tData[t]));
6544	if(rc)
6545	{
6546	printf("ERROR; return code from pthread_create() is %d\n", rc);
6547	exit(-1);
6548	}
6549	}
6550	}
6551
6552
6553
6554	#endif
6555
6556
6557	/*************************************************************************************************************************************************************/
6558
6559	typedef struct {
6560	double lh;
6561	int tree;
6562	double weight;
6563	} elw;
6564
6565	static int elwCompare(const void p1, const void p2)
6566	{
6567	elw rc1 = (elw )p1;
6568	elw rc2 = (elw )p2;
6569
6570	double i = rc1->weight;
6571	double j = rc2->weight;
6572
6573	if (i > j)
6574	return (-1);
6575	if (i < j)
6576	return (1);
6577	return (0);
6578	}
6579
6580
6581
6582
6583
6584
6585	static void computeELW(tree tr, analdef adef, char *bootStrapFileName)
6586	{
6587	int
6588	numberOfTrees = 0,
6589	i, k;
6590	char ch;
6591
6592	/*
6593	double
6594	bestLH = unlikely,
6595	elwSum = 0.0;
6596	*/
6597
6598	FILE *infoFile;
6599	int originalRateCategories = (int)malloc(tr->cdta->endsite * sizeof(int));
6600	int originalInvariant = (int)malloc(tr->cdta->endsite * sizeof(int));
6601	long startSeed;
6602	double **lhs;
6603	double **lhweights;
6604	elw *bootweights;
6605
6606	infoFile = fopen(infoFileName, "a");
6607
6608	initModel(tr, tr->rdta, tr->cdta, adef);
6609	allocNodex(tr, adef);
6610
6611	INFILE = fopen(bootStrapFileName, "r");
6612	while((ch = getc(INFILE)) != EOF)
6613	{
6614	if(ch == ';')
6615	numberOfTrees++;
6616	}
6617	rewind(INFILE);
6618
6619	if(numberOfTrees < 2)
6620	{
6621	printf("Error, there is only one tree in file %s which you want to use to conduct an ELW test\n", bootStrapFileName);
6622
6623	exit(-1);
6624	}
6625
6626	printf("\n\nFound %d trees in File %s\n\n", numberOfTrees, bootStrapFileName);
6627	fprintf(infoFile, "\n\nFound %d trees in File %s\n\n", numberOfTrees, bootStrapFileName);
6628
6629	bootweights = (elw )malloc(sizeof(elw) numberOfTrees);
6630
6631	lhs = (double *)malloc(sizeof(double ) * numberOfTrees);
6632
6633	for(k = 0; k < numberOfTrees; k++)
6634	lhs[k] = (double )malloc(sizeof(double) adef->multipleRuns);
6635
6636	lhweights = (double *)malloc(sizeof(double ) * numberOfTrees);
6637
6638	for(k = 0; k < numberOfTrees; k++)
6639	lhweights[k] = (double )malloc(sizeof(double) adef->multipleRuns);
6640
6641
6642	treeReadLen(INFILE, tr, adef);
6643	modOpt(tr, adef);
6644	rewind(INFILE);
6645
6646	/*
6647	This is for testing only !
6648	for(i = 0; i < numberOfTrees; i++)
6649	{
6650	treeReadLen(INFILE, tr, adef);
6651	treeEvaluate(tr, 2.0);
6652	bootweights[i].lh = tr->likelihood;
6653	}
6654	rewind(INFILE);
6655	*/
6656
6657	printf("Model optimization, first Tree: %f\n", tr->likelihood);
6658	fprintf(infoFile, "Model optimization, first Tree: %f\n", tr->likelihood);
6659
6660
6661	memcpy(originalRateCategories, tr->cdta->rateCategory, sizeof(int) * tr->cdta->endsite);
6662	memcpy(originalInvariant, tr->invariant, sizeof(int) * tr->cdta->endsite);
6663
6664	assert(adef->boot > 0);
6665	/* this is ugly, should be passed as param to computenextreplicate() */
6666	startSeed = adef->boot;
6667
6668
6669	for(i = 0; i < numberOfTrees; i++)
6670	{
6671	treeReadLen(INFILE, tr, adef);
6672	resetBranches(tr);
6673	adef->rapidBoot = startSeed;
6674
6675	for(k = 0; k < adef->multipleRuns; k++)
6676	{
6677	computeNextReplicate(tr, adef, originalRateCategories, originalInvariant);
6678
6679	if(k == 0)
6680	treeEvaluate(tr, 2.0);
6681	else
6682	treeEvaluate(tr, 0.5);
6683	/printf("%d %d %f\n", i, k, tr->likelihood);/
6684	lhs[i][k] = tr->likelihood;
6685	}
6686
6687	reductionCleanup(tr, adef, originalRateCategories, originalInvariant);
6688	}
6689
6690
6691
6692	for(k = 0; k < adef->multipleRuns; k++)
6693	{
6694	double best = unlikely;
6695	double sum = 0.0;
6696
6697	for(i = 0; i < numberOfTrees; i++)
6698	if(lhs[i][k] > best)
6699	best = lhs[i][k];
6700
6701	for(i = 0; i < numberOfTrees; i++)
6702	lhweights[i][k] = exp(lhs[i][k] - best);
6703
6704	for(i = 0; i < numberOfTrees; i++)
6705	sum += lhweights[i][k];
6706
6707	for(i = 0; i < numberOfTrees; i++)
6708	lhweights[i][k] = lhweights[i][k] / sum;
6709
6710	}
6711
6712
6713
6714	for(i = 0; i < numberOfTrees; i++)
6715	{
6716	double sum = 0.0;
6717
6718	for(k = 0; k < adef->multipleRuns; k++)
6719	sum += lhweights[i][k];
6720
6721	bootweights[i].weight = sum / ((double)adef->multipleRuns);
6722	bootweights[i].tree = i;
6723	}
6724
6725	qsort(bootweights, numberOfTrees, sizeof(elw), elwCompare);
6726
6727
6728	{
6729	double sum = 0.0;
6730
6731	/printf("Tree\t Posterior Probability \t Cumulative posterior probability \t Original Likelihood\n");/
6732	printf("Tree\t Posterior Probability \t Cumulative posterior probability\n");
6733	fprintf(infoFile, "Tree\t Posterior Probability \t Cumulative posterior probability\n");
6734	for(i = 0; i < numberOfTrees; i++)
6735	{
6736	sum += bootweights[i].weight;
6737	/printf("%d\t\t %f \t\t %f \t\t\t %f\n", bootweights[i].tree, bootweights[i].weight, sum, bootweights[i].lh);/
6738	printf("%d\t\t %f \t\t %f\n", bootweights[i].tree, bootweights[i].weight, sum);
6739	fprintf(infoFile, "%d\t\t %f \t\t %f\n", bootweights[i].tree, bootweights[i].weight, sum);
6740	}
6741	}
6742
6743	free(originalRateCategories);
6744	free(originalInvariant);
6745
6746	fclose(INFILE);
6747	fclose(infoFile);
6748	exit(0);
6749	}
6750
6751
6752
6753
6754
6755	int main (int argc, char *argv[])
6756	{
6757	rawdata *rdta;
6758	cruncheddata *cdta;
6759	tree *tr;
6760	analdef *adef;
6761
6762	#ifdef PARALLEL
6763	MPI_Init(&argc, &argv);
6764	MPI_Comm_rank(MPI_COMM_WORLD, &processID);
6765	MPI_Comm_size(MPI_COMM_WORLD, &numOfWorkers);
6766	if(processID == 0)
6767	printf("\nThis is the RAxML MPI Master process\n");
6768	else
6769	printf("\nThis is the RAxML MPI Worker Process Number: %d\n", processID);
6770	#else
6771	processID = 0;
6772	#endif
6773
6774	masterTime = gettime();
6775
6776	adef = (analdef *)malloc(sizeof(analdef));
6777	rdta = (rawdata *)malloc(sizeof(rawdata));
6778	cdta = (cruncheddata *)malloc(sizeof(cruncheddata));
6779	tr = (tree *)malloc(sizeof(tree));
6780
6781	initAdef(adef);
6782	get_args(argc,argv, adef, tr);
6783
6784	if(adef->model == M_PROTCAT \|\| adef->model == M_GTRCAT)
6785	tr->rateHetModel = CAT;
6786	else
6787	{
6788	if(adef->useInvariant)
6789	tr->rateHetModel = GAMMA_I;
6790	else
6791	tr->rateHetModel = GAMMA;
6792	}
6793
6794	/*
6795	This is a very ugly numerical bug fix, that intends to avoid the unaesthetic phenomena
6796	that can occur during model param optimization due to the dependency between parameters
6797	alpha and invar which are NOT independent from each other.
6798	When using P-Invar set likelihood epsilon to a lower value!
6799
6800	TODO-MIX this is very ugly !
6801
6802	*/
6803
6804	if(adef->useInvariant && adef->likelihoodEpsilon > 0.001)
6805	adef->likelihoodEpsilon = 0.001;
6806
6807	readData(adef, rdta, cdta, tr);
6808
6809	checkOutgroups(tr, adef);
6810	makeFileNames();
6811
6812	if(adef->useExcludeFile)
6813	{
6814	handleExcludeFile(tr, adef, rdta);
6815	exit(0);
6816	}
6817
6818	if(adef->mode != SEQUENCE_SIMILARITY_FILTER)
6819	{
6820	checkSequences(tr, rdta, adef);
6821	}
6822	else
6823	{
6824	reduceBySequenceSimilarity(tr, rdta, adef);
6825	exit(0);
6826	}
6827
6828	if(adef->mode == SPLIT_MULTI_GENE)
6829	{
6830	splitMultiGene(tr, rdta);
6831	exit(0);
6832	}
6833
6834	if(adef->mode == CHECK_ALIGNMENT)
6835	{
6836	printf("Alignment format can be read by RAxML \n");
6837	exit(0);
6838	}
6839
6840	makeweights(adef, rdta, cdta, tr);
6841	makevalues(rdta, cdta, tr, adef);
6842
6843	if(adef->generateBS)
6844	{
6845	generateBS(tr, adef);
6846	exit(0);
6847	}
6848
6849	#ifdef _USE_PTHREADS
6850	startPthreads(tr);
6851	#endif
6852
6853	if(adef->computeELW)
6854	computeELW(tr, adef, bootStrapFile);
6855
6856	if(adef->boot)
6857	makeboot(adef, tr);
6858
6859	initModel(tr, rdta, cdta, adef);
6860
6861	printModelAndProgramInfo(tr, adef, argc, argv);
6862
6863	if(adef->bootStopOnly > 0)
6864	{
6865	computeBootStopOnly(tr, adef, bootStrapFile);
6866	exit(0);
6867	}
6868
6869	switch(adef->mode)
6870	{
6871	case ARNDT_MODE:
6872	printf("OPT_ARNDT\n");
6873	getStartingTree(tr, adef);
6874	optimizeArndt(tr, adef);
6875	break;
6876	case MEHRING_ALGO:
6877	getStartingTree(tr, adef);
6878	determineSequencePosition(tr, adef);
6879	break;
6880	case PARSIMONY_ADDITION:
6881	getStartingTree(tr, adef);
6882	printStartingTree(tr, adef, TRUE);
6883	break;
6884	case OPTIMIZE_RATES:
6885	if(adef->computePerSiteLLs)
6886	computePerSiteLLs(tr, adef, bootStrapFile);
6887	else
6888	optimizeRatesOnly(tr, adef);
6889	break;
6890	case TREE_EVALUATION:
6891	getStartingTree(tr, adef);
6892	if(adef->likelihoodTest)
6893	computeLHTest(tr, adef, bootStrapFile);
6894	else
6895	{
6896	modOpt(tr, adef);
6897	printLog(tr, adef, TRUE);
6898	printResult(tr, adef, TRUE);
6899	break;
6900	}
6901	case CALC_BIPARTITIONS:
6902	calcBipartitions(tr, adef, tree_file, bootStrapFile);
6903	break;
6904	case BIG_RAPID_MODE:
6905	if(adef->boot)
6906	doBootstrap(tr, adef, rdta, cdta);
6907	else
6908	{
6909	if(adef->rapidBoot)
6910	{
6911	#ifdef _VINCENT
6912	doAllInOneVincent(tr, adef);
6913	#else
6914	doAllInOne(tr, adef);
6915	#endif
6916	}
6917	else
6918	doInference(tr, adef, rdta, cdta);
6919	}
6920	break;
6921	default:
6922	assert(0);
6923	}
6924
6925	finalizeInfoFile(tr, adef);
6926
6927	#ifdef PARALLEL
6928	MPI_Finalize();
6929	#endif
6930
6931	return 0;
6932	}

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

Download in other formats: