Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

optimizeModel.c

Visit:

Last change on this file was 10409, checked in by aboeckma, 12 years ago
Updated raxml to current version
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 93.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
28	* of taxa and mixed models".
29	* Bioinformatics 2006; doi: 10.1093/bioinformatics/btl446
30	*/
31
32	#ifndef WIN32
33	#include <unistd.h>
34	#endif
35
36	#include <math.h>
37	#include <time.h>
38	#include <stdlib.h>
39	#include <stdio.h>
40	#include <ctype.h>
41	#include <string.h>
42	#include "axml.h"
43
44
45	static const double MNBRAK_GOLD = 1.618034;
46	static const double MNBRAK_TINY = 1.e-20;
47	static const double MNBRAK_GLIMIT = 100.0;
48	static const double BRENT_ZEPS = 1.e-5;
49	static const double BRENT_CGOLD = 0.3819660;
50
51	extern int optimizeRatesInvocations;
52	extern int optimizeRateCategoryInvocations;
53	extern int optimizeAlphaInvocations;
54	extern int optimizeInvarInvocations;
55	extern double masterTime;
56	extern char ratesFileName[1024];
57	extern char workdir[1024];
58	extern char run_id[128];
59	extern char lengthFileName[1024];
60	extern char lengthFileNameModel[1024];
61	extern char *protModels[NUM_PROT_MODELS];
62
63
64	#ifdef _USE_PTHREADS
65	extern volatile int NumberOfThreads;
66	extern volatile double *reductionBuffer;
67	#endif
68
69	/* TODO remove at some point */
70	#define _DEBUG_MODEL_OPTIMIZATION
71
72	static void brentGeneric(double ax, double bx, double cx, double fb, double tol, double xmin, double result, int numberOfModels,
73	int whichFunction, int rateNumber, tree tr, linkageList ll, double lim_inf, double lim_sup);
74
75	static int brakGeneric(double param, double ax, double bx, double cx, double fa, double fb,
76	double *fc, double lim_inf, double lim_sup,
77	int numberOfModels, int rateNumber, int whichFunction, tree tr, linkageList ll, double modelEpsilon);
78
79	/*******************FUNCTIONS FOOR EXACT MODEL OPTIMIZATION UNDER GTRGAMMA *************************************/
80
81
82	static void setRateModel(tree *tr, int model, double rate, int position)
83	{
84	int
85	states = tr->partitionData[model].states,
86	numRates = (states * states - states) / 2;
87
88	if(tr->partitionData[model].dataType == DNA_DATA)
89	assert(position >= 0 && position < (numRates - 1));
90	else
91	assert(position >= 0 && position < numRates);
92
93	assert(tr->partitionData[model].dataType != BINARY_DATA);
94
95	if(!(tr->partitionData[model].dataType == SECONDARY_DATA \|\|
96	tr->partitionData[model].dataType == SECONDARY_DATA_6 \|\|
97	tr->partitionData[model].dataType == SECONDARY_DATA_7))
98	assert(rate >= RATE_MIN && rate <= RATE_MAX);
99
100	if(tr->partitionData[model].nonGTR)
101	{
102	int
103	i,
104	k = tr->partitionData[model].symmetryVector[position];
105
106	assert(tr->partitionData[model].dataType == SECONDARY_DATA \|\|
107	tr->partitionData[model].dataType == SECONDARY_DATA_6 \|\|
108	tr->partitionData[model].dataType == SECONDARY_DATA_7);
109
110	if(k == -1)
111	tr->partitionData[model].substRates[position] = 0.0;
112	else
113	{
114	if(k == tr->partitionData[model].symmetryVector[numRates - 1])
115	{
116	for(i = 0; i < numRates - 1; i++)
117	if(tr->partitionData[model].symmetryVector[i] == k)
118	tr->partitionData[model].substRates[position] = 1.0;
119	}
120	else
121	{
122	for(i = 0; i < numRates - 1; i++)
123	{
124	if(tr->partitionData[model].symmetryVector[i] == k)
125	tr->partitionData[model].substRates[i] = rate;
126	}
127	}
128	}
129	}
130	else
131	tr->partitionData[model].substRates[position] = rate;
132	}
133
134
135
136
137
138	static linkageList* initLinkageList(int linkList, tree tr)
139	{
140	int
141	k,
142	partitions,
143	numberOfModels = 0,
144	i,
145	pos;
146	linkageList* ll = (linkageList*)rax_malloc(sizeof(linkageList));
147
148	for(i = 0; i < tr->NumberOfModels; i++)
149	{
150	if(linkList[i] > numberOfModels)
151	numberOfModels = linkList[i];
152	}
153
154	numberOfModels++;
155
156	ll->entries = numberOfModels;
157	ll->ld = (linkageData)rax_malloc(sizeof(linkageData) numberOfModels);
158
159
160	for(i = 0; i < numberOfModels; i++)
161	{
162	ll->ld[i].valid = TRUE;
163	partitions = 0;
164
165	for(k = 0; k < tr->NumberOfModels; k++)
166	if(linkList[k] == i)
167	partitions++;
168
169	ll->ld[i].partitions = partitions;
170	ll->ld[i].partitionList = (int)rax_malloc(sizeof(int) partitions);
171
172	for(k = 0, pos = 0; k < tr->NumberOfModels; k++)
173	if(linkList[k] == i)
174	ll->ld[i].partitionList[pos++] = k;
175	}
176
177	return ll;
178	}
179
180
181	static linkageList* initLinkageListGTR(tree *tr)
182	{
183	int
184	i,
185	links = (int)rax_malloc(sizeof(int) * tr->NumberOfModels),
186	firstAA = tr->NumberOfModels + 2,
187	countGTR = 0,
188	countUnlinkedGTR = 0,
189	countOtherModel = 0;
190
191	linkageList*
192	ll;
193
194	for(i = 0; i < tr->NumberOfModels; i++)
195	{
196	if(tr->partitionData[i].dataType == AA_DATA)
197	{
198	if(tr->partitionData[i].protModels == GTR)
199	{
200	if(i < firstAA)
201	firstAA = i;
202	countGTR++;
203	}
204	else
205	{
206	if(tr->partitionData[i].protModels == GTR_UNLINKED)
207	countUnlinkedGTR++;
208	else
209	countOtherModel++;
210	}
211	}
212	}
213
214	/*
215	TODO need to think what we actually want here !
216	Shall we mix everything: linked, unlinked WAG etc?
217	*/
218
219	assert((countGTR > 0 && countOtherModel == 0) \|\| (countGTR == 0 && countOtherModel > 0) \|\| (countGTR == 0 && countOtherModel == 0));
220
221	if(countGTR == 0)
222	{
223	for(i = 0; i < tr->NumberOfModels; i++)
224	links[i] = i;
225	}
226	else
227	{
228	for(i = 0; i < tr->NumberOfModels; i++)
229	{
230	switch(tr->partitionData[i].dataType)
231	{
232	case DNA_DATA:
233	case BINARY_DATA:
234	case GENERIC_32:
235	case GENERIC_64:
236	case SECONDARY_DATA:
237	case SECONDARY_DATA_6:
238	case SECONDARY_DATA_7:
239	links[i] = i;
240	break;
241	case AA_DATA:
242	links[i] = firstAA;
243	break;
244	default:
245	assert(0);
246	}
247	}
248	}
249
250
251	ll = initLinkageList(links, tr);
252
253	rax_free(links);
254
255	return ll;
256	}
257
258
259
260	static void freeLinkageList( linkageList* ll)
261	{
262	int i;
263
264	for(i = 0; i < ll->entries; i++)
265	rax_free(ll->ld[i].partitionList);
266
267	rax_free(ll->ld);
268	rax_free(ll);
269	}
270
271	#define ALPHA_F 0
272	#define INVAR_F 1
273	#define RATE_F 2
274	#define SCALER_F 3
275	#define LXRATE_F 4
276	#define LXWEIGHT_F 5
277
278	static void updateWeights(tree *tr, int model, int rate, double value)
279	{
280	int
281	j;
282
283	double
284	w = 0.0;
285
286	tr->partitionData[model].weightExponents[rate] = value;
287
288	for(j = 0; j < 4; j++)
289	w += exp(tr->partitionData[model].weightExponents[j]);
290
291	for(j = 0; j < 4; j++)
292	tr->partitionData[model].weights[j] = exp(tr->partitionData[model].weightExponents[j]) / w;
293	}
294
295	static void optLG4X_Weights(tree tr, linkageList ll, int numberOfModels, int weightNumber, double modelEpsilon)
296	{
297	int
298	pos,
299	i;
300
301	double
302	lim_inf = -1000000.0,
303	lim_sup = 200.0,
304	startLH = (double )rax_malloc(sizeof(double) * numberOfModels),
305	endLH = (double )rax_malloc(sizeof(double) * numberOfModels),
306	startWeights = (double )rax_malloc(sizeof(double) * numberOfModels * 4),
307	startExponents = (double )rax_malloc(sizeof(double) * numberOfModels * 4),
308	endWeights = (double )rax_malloc(sizeof(double) * numberOfModels),
309	endExponents = (double )rax_malloc(sizeof(double) * numberOfModels),
310	_a = (double )rax_malloc(sizeof(double) * numberOfModels),
311	_b = (double )rax_malloc(sizeof(double) * numberOfModels),
312	_c = (double )rax_malloc(sizeof(double) * numberOfModels),
313	_fa = (double )rax_malloc(sizeof(double) * numberOfModels),
314	_fb = (double )rax_malloc(sizeof(double) * numberOfModels),
315	_fc = (double )rax_malloc(sizeof(double) * numberOfModels),
316	_param = (double )rax_malloc(sizeof(double) * numberOfModels),
317	result = (double )rax_malloc(sizeof(double) * numberOfModels),
318	_x = (double )rax_malloc(sizeof(double) * numberOfModels);
319
320	evaluateGeneric(tr, tr->start);
321
322	#ifdef _DEBUG_MODEL_OPTIMIZATION
323	double
324	initialLH = tr->likelihood;
325	#endif
326
327	assert(weightNumber >= 0 && weightNumber < 4);
328
329	/*
330	at this point here every worker has the traversal data it needs for the
331	search, so we won't re-distribute it he he :-)
332	*/
333
334	for(i = 0, pos = 0; i < ll->entries; i++)
335	{
336	if(ll->ld[i].valid)
337	{
338	int
339	index = ll->ld[i].partitionList[0];
340
341
342	assert(ll->ld[i].partitions == 1);
343
344	memcpy(&startWeights[pos * 4], tr->partitionData[index].weights, 4 * sizeof(double));
345	memcpy(&startExponents[pos * 4], tr->partitionData[index].weightExponents, 4 * sizeof(double));
346
347	_a[pos] = startExponents[pos * 4 + weightNumber] + 0.1;
348	_b[pos] = startExponents[pos * 4 + weightNumber] - 0.1;
349
350	if(_a[pos] < lim_inf)
351	_a[pos] = lim_inf;
352
353	if(_a[pos] > lim_sup)
354	_a[pos] = lim_sup;
355
356	if(_b[pos] < lim_inf)
357	_b[pos] = lim_inf;
358
359	if(_b[pos] > lim_sup)
360	_b[pos] = lim_sup;
361
362	startLH[pos] = tr->perPartitionLH[index];
363	endLH[pos] = unlikely;
364
365	pos++;
366	}
367	}
368	#ifdef _USE_PTHREADS
369	masterBarrier(THREAD_COPY_LG4X_RATES, tr);
370	#endif
371	assert(pos == numberOfModels);
372
373	brakGeneric(_param, _a, _b, _c, _fa, _fb, _fc, lim_inf, lim_sup, numberOfModels, weightNumber, LXWEIGHT_F, tr, ll, modelEpsilon);
374	brentGeneric(_a, _b, _c, _fb, modelEpsilon, _x, result, numberOfModels, LXWEIGHT_F, weightNumber, tr, ll, lim_inf, lim_sup);
375
376	/* now calculate the likelihood with the optimized rate */
377
378	for(i = 0, pos = 0; i < ll->entries; i++)
379	{
380	if(ll->ld[i].valid)
381	{
382	int
383	index = ll->ld[i].partitionList[0];
384
385	assert(ll->ld[i].partitions == 1);
386
387	tr->executeModel[index] = TRUE;
388
389	updateWeights(tr, index, weightNumber, _x[pos]);
390
391	pos++;
392	}
393	}
394
395	#ifdef _USE_PTHREADS
396	masterBarrier(THREAD_COPY_LG4X_RATES, tr);
397	#endif
398
399	evaluateGeneric(tr, tr->start);
400
401	for(i = 0, pos = 0; i < ll->entries; i++)
402	{
403	if(ll->ld[i].valid)
404	{
405	int
406	index = ll->ld[i].partitionList[0];
407
408	assert(ll->ld[i].partitions == 1);
409
410	if(tr->executeModel[index])
411	{
412	endLH[pos] = tr->perPartitionLH[index];
413
414	if(startLH[pos] > endLH[pos])
415	{
416	memcpy(tr->partitionData[index].weights, &startWeights[pos * 4], sizeof(double) * 4);
417	memcpy(tr->partitionData[index].weightExponents, &startExponents[pos * 4], sizeof(double) * 4);
418	}
419	}
420
421	pos++;
422	}
423	}
424
425	assert(pos == numberOfModels);
426
427	#ifdef _USE_PTHREADS
428	masterBarrier(THREAD_COPY_LG4X_RATES, tr);
429	#endif
430
431
432	#ifdef _DEBUG_MODEL_OPTIMIZATION
433	evaluateGeneric(tr, tr->start);
434
435	if(tr->likelihood < initialLH)
436	printf("%f %f\n", tr->likelihood, initialLH);
437	assert(tr->likelihood >= initialLH);
438	#endif
439
440	rax_free(endExponents);
441	rax_free(startExponents);
442	rax_free(startLH);
443	rax_free(endLH);
444	rax_free(startWeights);
445	rax_free(endWeights);
446	rax_free(result);
447	rax_free(_a);
448	rax_free(_b);
449	rax_free(_c);
450	rax_free(_fa);
451	rax_free(_fb);
452	rax_free(_fc);
453	rax_free(_param);
454	rax_free(_x);
455	}
456
457	static void optimizeWeights(tree tr, double modelEpsilon, linkageList ll, int numberOfModels)
458	{
459	int
460	i;
461
462	double
463	initialLH = 0.0,
464	finalLH = 0.0;
465
466	evaluateGeneric(tr, tr->start);
467
468	initialLH = tr->likelihood;
469	//printf("W: %f %f [%f] ->", tr->perPartitionLH[0], tr->perPartitionLH[1], initialLH);
470
471	for(i = 0; i < 4; i++)
472	optLG4X_Weights(tr, ll, numberOfModels, i, modelEpsilon);
473
474	#ifdef _USE_PTHREADS
475	masterBarrier(THREAD_COPY_LG4X_RATES, tr);
476	#endif
477
478	evaluateGenericInitrav(tr, tr->start);
479
480	finalLH = tr->likelihood;
481
482	if(finalLH < initialLH)
483	printf("Final: %f initial: %f\n", finalLH, initialLH);
484	assert(finalLH >= initialLH);
485
486	//printf("%f %f [%f]\n", tr->perPartitionLH[0], tr->perPartitionLH[1], finalLH);
487	}
488
489	static void evaluateChange(tree tr, int rateNumber, double value, double result, boolean converged, int whichFunction, int numberOfModels, linkageList *ll, double modelEpsilon)
490	{
491	int i, k, pos;
492
493	switch(whichFunction)
494	{
495	case ALPHA_F:
496	for(i = 0, pos = 0; i < ll->entries; i++)
497	{
498	int
499	index = ll->ld[i].partitionList[0];
500
501	assert(ll->ld[i].partitions == 1);
502
503	if(ll->ld[i].valid)
504	{
505	if(converged[pos])
506	tr->executeModel[index] = FALSE;
507	else
508	{
509	tr->executeModel[index] = TRUE;
510	tr->partitionData[index].alpha = value[pos];
511	makeGammaCats(tr->partitionData[index].alpha, tr->partitionData[index].gammaRates, 4, tr->useGammaMedian);
512	}
513
514	pos++;
515	}
516	else
517	tr->executeModel[index] = FALSE;
518	}
519
520	assert(pos == numberOfModels);
521
522	#ifdef _USE_PTHREADS
523	{
524	volatile double result;
525
526	masterBarrier(THREAD_OPT_ALPHA, tr);
527	if(tr->NumberOfModels == 1)
528	{
529	for(i = 0, result = 0.0; i < NumberOfThreads; i++)
530	result += reductionBuffer[i];
531	tr->perPartitionLH[0] = result;
532	}
533	else
534	{
535	int j;
536	volatile double partitionResult;
537
538	result = 0.0;
539
540	for(j = 0; j < tr->NumberOfModels; j++)
541	{
542	for(i = 0, partitionResult = 0.0; i < NumberOfThreads; i++)
543	partitionResult += reductionBuffer[i * tr->NumberOfModels + j];
544	result += partitionResult;
545	tr->perPartitionLH[j] = partitionResult;
546	}
547	}
548	}
549	#else
550	evaluateGenericInitrav(tr, tr->start);
551	#endif
552
553	for(i = 0, pos = 0; i < ll->entries; i++)
554	{
555	int
556	index = ll->ld[i].partitionList[0];
557
558	assert(ll->ld[i].partitions == 1);
559
560	if(ll->ld[i].valid)
561	{
562	result[pos] = 0.0;
563
564	assert(tr->perPartitionLH[index] <= 0.0);
565
566	result[pos] -= tr->perPartitionLH[index];
567
568	pos++;
569	}
570
571	tr->executeModel[index] = TRUE;
572	}
573
574	assert(pos == numberOfModels);
575
576	break;
577	case INVAR_F:
578
579	for(i = 0; i < ll->entries; i++)
580	{
581	if(converged[i])
582	{
583	for(k = 0; k < ll->ld[i].partitions; k++)
584	tr->executeModel[ll->ld[i].partitionList[k]] = FALSE;
585	}
586	else
587	{
588	for(k = 0; k < ll->ld[i].partitions; k++)
589	{
590	int index = ll->ld[i].partitionList[k];
591	tr->executeModel[index] = TRUE;
592	tr->partitionData[index].propInvariant = value[i];
593	}
594	}
595	}
596
597	#ifdef _USE_PTHREADS
598	{
599	volatile double result;
600
601	masterBarrier(THREAD_OPT_INVAR, tr);
602	if(tr->NumberOfModels == 1)
603	{
604	for(i = 0, result = 0.0; i < NumberOfThreads; i++)
605	result += reductionBuffer[i];
606	tr->perPartitionLH[0] = result;
607	}
608	else
609	{
610	int j;
611	volatile double partitionResult;
612
613	result = 0.0;
614
615	for(j = 0; j < tr->NumberOfModels; j++)
616	{
617	for(i = 0, partitionResult = 0.0; i < NumberOfThreads; i++)
618	partitionResult += reductionBuffer[i * tr->NumberOfModels + j];
619	result += partitionResult;
620	tr->perPartitionLH[j] = partitionResult;
621	}
622	}
623	}
624	#else
625	evaluateGeneric(tr, tr->start);
626	#endif
627
628	for(i = 0; i < ll->entries; i++)
629	{
630	result[i] = 0.0;
631
632	for(k = 0; k < ll->ld[i].partitions; k++)
633	{
634	int index = ll->ld[i].partitionList[k];
635
636	assert(tr->perPartitionLH[index] <= 0.0);
637
638	result[i] -= tr->perPartitionLH[index];
639	tr->executeModel[index] = TRUE;
640	}
641	}
642
643	break;
644	case RATE_F:
645	for(i = 0, pos = 0; i < ll->entries; i++)
646	{
647	if(ll->ld[i].valid)
648	{
649	if(converged[pos])
650	{
651	for(k = 0; k < ll->ld[i].partitions; k++)
652	tr->executeModel[ll->ld[i].partitionList[k]] = FALSE;
653	}
654	else
655	{
656	for(k = 0; k < ll->ld[i].partitions; k++)
657	{
658	int index = ll->ld[i].partitionList[k];
659	setRateModel(tr, index, value[pos], rateNumber);
660	initReversibleGTR(tr, index);
661	}
662	}
663	pos++;
664	}
665	else
666	{
667	for(k = 0; k < ll->ld[i].partitions; k++)
668	tr->executeModel[ll->ld[i].partitionList[k]] = FALSE;
669	}
670
671	}
672
673	assert(pos == numberOfModels);
674
675	if(tr->useBrLenScaler)
676	determineFullTraversal(tr->start, tr);
677
678	#ifdef _USE_PTHREADS
679	{
680	volatile double result;
681
682	masterBarrier(THREAD_OPT_RATE, tr);
683	if(tr->NumberOfModels == 1)
684	{
685	for(i = 0, result = 0.0; i < NumberOfThreads; i++)
686	result += reductionBuffer[i];
687	tr->perPartitionLH[0] = result;
688	}
689	else
690	{
691	int j;
692	volatile double partitionResult;
693
694	result = 0.0;
695
696	for(j = 0; j < tr->NumberOfModels; j++)
697	{
698	for(i = 0, partitionResult = 0.0; i < NumberOfThreads; i++)
699	partitionResult += reductionBuffer[i * tr->NumberOfModels + j];
700	result += partitionResult;
701	tr->perPartitionLH[j] = partitionResult;
702	}
703	}
704	}
705	#else
706	evaluateGenericInitrav(tr, tr->start);
707	#endif
708
709
710	for(i = 0, pos = 0; i < ll->entries; i++)
711	{
712	if(ll->ld[i].valid)
713	{
714	result[pos] = 0.0;
715	for(k = 0; k < ll->ld[i].partitions; k++)
716	{
717	int index = ll->ld[i].partitionList[k];
718
719	assert(tr->perPartitionLH[index] <= 0.0);
720
721	result[pos] -= tr->perPartitionLH[index];
722
723	}
724	pos++;
725	}
726	for(k = 0; k < ll->ld[i].partitions; k++)
727	{
728	int index = ll->ld[i].partitionList[k];
729	tr->executeModel[index] = TRUE;
730	}
731	}
732
733	assert(pos == numberOfModels);
734	break;
735	case SCALER_F:
736	assert(ll->entries == tr->NumberOfModels);
737	assert(ll->entries == tr->numBranches);
738	for(i = 0; i < ll->entries; i++)
739	{
740	if(converged[i])
741	{
742	for(k = 0; k < ll->ld[i].partitions; k++)
743	tr->executeModel[ll->ld[i].partitionList[k]] = FALSE;
744	}
745	else
746	{
747	for(k = 0; k < ll->ld[i].partitions; k++)
748	{
749	int
750	index = ll->ld[i].partitionList[k];
751
752	tr->executeModel[index] = TRUE;
753	tr->partitionData[index].brLenScaler = value[i];
754	scaleBranches(tr, FALSE);
755	}
756	}
757	}
758
759	#ifdef _USE_PTHREADS
760	{
761	volatile
762	double result;
763
764	/* need to call this here because we need to copy the changed branch lengths
765	(due to changed) scalers into the traversal descriptor */
766
767	determineFullTraversal(tr->start, tr);
768
769	masterBarrier(THREAD_OPT_SCALER, tr);
770	if(tr->NumberOfModels == 1)
771	{
772	for(i = 0, result = 0.0; i < NumberOfThreads; i++)
773	result += reductionBuffer[i];
774	tr->perPartitionLH[0] = result;
775	}
776	else
777	{
778	int j;
779	volatile double partitionResult;
780
781	result = 0.0;
782
783	for(j = 0; j < tr->NumberOfModels; j++)
784	{
785	for(i = 0, partitionResult = 0.0; i < NumberOfThreads; i++)
786	partitionResult += reductionBuffer[i * tr->NumberOfModels + j];
787	result += partitionResult;
788	tr->perPartitionLH[j] = partitionResult;
789	}
790	}
791	}
792	#else
793	evaluateGenericInitrav(tr, tr->start);
794	#endif
795
796	for(i = 0; i < ll->entries; i++)
797	{
798	result[i] = 0.0;
799
800	for(k = 0; k < ll->ld[i].partitions; k++)
801	{
802	int index = ll->ld[i].partitionList[k];
803
804	assert(tr->perPartitionLH[index] <= 0.0);
805
806	result[i] -= tr->perPartitionLH[index];
807	tr->executeModel[index] = TRUE;
808	}
809	}
810	break;
811	case LXRATE_F:
812	{
813	boolean
814	atLeastOnePartition = FALSE;
815
816	assert(rateNumber != -1);
817
818	for(i = 0, pos = 0; i < ll->entries; i++)
819	{
820	int
821	index = ll->ld[i].partitionList[0];
822
823	assert(ll->ld[i].partitions == 1);
824
825	if(ll->ld[i].valid)
826	{
827	if(converged[pos])
828	tr->executeModel[index] = FALSE;
829	else
830	{
831	atLeastOnePartition = TRUE;
832	tr->executeModel[index] = TRUE;
833	tr->partitionData[index].gammaRates[rateNumber] = value[pos];
834	}
835
836	pos++;
837	}
838	else
839	tr->executeModel[index] = FALSE;
840	}
841
842	assert(pos == numberOfModels);
843
844	#ifdef _USE_PTHREADS
845	{
846	volatile double result;
847
848	masterBarrier(THREAD_OPT_LG4X_RATES, tr);
849
850	if(tr->NumberOfModels == 1)
851	{
852	for(i = 0, result = 0.0; i < NumberOfThreads; i++)
853	result += reductionBuffer[i];
854	tr->perPartitionLH[0] = result;
855	}
856	else
857	{
858	int j;
859	volatile double partitionResult;
860
861	result = 0.0;
862
863	for(j = 0; j < tr->NumberOfModels; j++)
864	{
865	for(i = 0, partitionResult = 0.0; i < NumberOfThreads; i++)
866	partitionResult += reductionBuffer[i * tr->NumberOfModels + j];
867	result += partitionResult;
868	tr->perPartitionLH[j] = partitionResult;
869	}
870	}
871	}
872
873	#else
874	evaluateGenericInitrav(tr, tr->start);
875	#endif
876
877	if(atLeastOnePartition)
878	{
879	boolean
880	buffer = (boolean)rax_malloc(tr->NumberOfModels * sizeof(boolean));
881
882	memcpy(buffer, tr->executeModel, sizeof(boolean) * tr->NumberOfModels);
883
884	for(i = 0; i < tr->NumberOfModels; i++)
885	tr->executeModel[i] = FALSE;
886
887	for(i = 0, pos = 0; i < ll->entries; i++)
888	{
889	int
890	index = ll->ld[i].partitionList[0];
891
892	if(ll->ld[i].valid)
893	tr->executeModel[index] = TRUE;
894	}
895
896	optimizeWeights(tr, modelEpsilon, ll, numberOfModels);
897
898	memcpy(tr->executeModel, buffer, sizeof(boolean) * tr->NumberOfModels);
899
900	rax_free(buffer);
901	}
902
903
904	for(i = 0, pos = 0; i < ll->entries; i++)
905	{
906	int
907	index = ll->ld[i].partitionList[0];
908
909	assert(ll->ld[i].partitions == 1);
910
911	if(ll->ld[i].valid)
912	{
913	result[pos] = 0.0;
914
915	assert(tr->perPartitionLH[index] <= 0.0);
916
917	result[pos] -= tr->perPartitionLH[index];
918
919	pos++;
920	}
921
922	tr->executeModel[index] = TRUE;
923	}
924
925	assert(pos == numberOfModels);
926	}
927	break;
928	case LXWEIGHT_F:
929	assert(rateNumber != -1);
930
931	//printf("%d %d\n", tr->executeModel[0], tr->executeModel[1]);
932
933	for(i = 0, pos = 0; i < ll->entries; i++)
934	{
935	int
936	index = ll->ld[i].partitionList[0];
937
938	assert(ll->ld[i].partitions == 1);
939
940	if(ll->ld[i].valid)
941	{
942	if(converged[pos])
943	tr->executeModel[index] = FALSE;
944	else
945	{
946	tr->executeModel[index] = TRUE;
947	updateWeights(tr, index, rateNumber, value[pos]);
948	}
949
950	pos++;
951	}
952	else
953	tr->executeModel[index] = FALSE;
954	}
955
956	assert(pos == numberOfModels);
957
958	#ifdef _USE_PTHREADS
959	{
960	volatile double result;
961
962	masterBarrier(THREAD_OPT_LG4X_RATES, tr);
963
964	if(tr->NumberOfModels == 1)
965	{
966	for(i = 0, result = 0.0; i < NumberOfThreads; i++)
967	result += reductionBuffer[i];
968	tr->perPartitionLH[0] = result;
969	}
970	else
971	{
972	int j;
973	volatile double partitionResult;
974
975	result = 0.0;
976
977	for(j = 0; j < tr->NumberOfModels; j++)
978	{
979	for(i = 0, partitionResult = 0.0; i < NumberOfThreads; i++)
980	partitionResult += reductionBuffer[i * tr->NumberOfModels + j];
981	result += partitionResult;
982	tr->perPartitionLH[j] = partitionResult;
983	}
984	}
985	}
986	#else
987	evaluateGeneric(tr, tr->start);
988	#endif
989
990	//printf("weights Like: %f\n", tr->perPartitionLH[0]);
991
992	for(i = 0, pos = 0; i < ll->entries; i++)
993	{
994	int
995	index = ll->ld[i].partitionList[0];
996
997	assert(ll->ld[i].partitions == 1);
998
999	if(ll->ld[i].valid)
1000	{
1001	result[pos] = 0.0;
1002
1003	assert(tr->perPartitionLH[index] <= 0.0);
1004
1005	result[pos] -= tr->perPartitionLH[index];
1006
1007	pos++;
1008	}
1009
1010	tr->executeModel[index] = TRUE;
1011	}
1012
1013	assert(pos == numberOfModels);
1014
1015	break;
1016	default:
1017	assert(0);
1018	}
1019
1020	}
1021
1022
1023
1024	static void brentGeneric(double ax, double bx, double cx, double fb, double tol, double xmin, double result, int numberOfModels,
1025	int whichFunction, int rateNumber, tree tr, linkageList ll, double lim_inf, double lim_sup)
1026	{
1027	int iter, i;
1028	double
1029	a = (double )rax_malloc(sizeof(double) * numberOfModels),
1030	b = (double )rax_malloc(sizeof(double) * numberOfModels),
1031	d = (double )rax_malloc(sizeof(double) * numberOfModels),
1032	etemp = (double )rax_malloc(sizeof(double) * numberOfModels),
1033	fu = (double )rax_malloc(sizeof(double) * numberOfModels),
1034	fv = (double )rax_malloc(sizeof(double) * numberOfModels),
1035	fw = (double )rax_malloc(sizeof(double) * numberOfModels),
1036	fx = (double )rax_malloc(sizeof(double) * numberOfModels),
1037	p = (double )rax_malloc(sizeof(double) * numberOfModels),
1038	q = (double )rax_malloc(sizeof(double) * numberOfModels),
1039	r = (double )rax_malloc(sizeof(double) * numberOfModels),
1040	tol1 = (double )rax_malloc(sizeof(double) * numberOfModels),
1041	tol2 = (double )rax_malloc(sizeof(double) * numberOfModels),
1042	u = (double )rax_malloc(sizeof(double) * numberOfModels),
1043	v = (double )rax_malloc(sizeof(double) * numberOfModels),
1044	w = (double )rax_malloc(sizeof(double) * numberOfModels),
1045	x = (double )rax_malloc(sizeof(double) * numberOfModels),
1046	xm = (double )rax_malloc(sizeof(double) * numberOfModels),
1047	e = (double )rax_malloc(sizeof(double) * numberOfModels);
1048	boolean converged = (boolean )rax_malloc(sizeof(boolean) * numberOfModels);
1049	boolean allConverged;
1050
1051	for(i = 0; i < numberOfModels; i++)
1052	converged[i] = FALSE;
1053
1054	for(i = 0; i < numberOfModels; i++)
1055	{
1056	e[i] = 0.0;
1057	d[i] = 0.0;
1058	}
1059
1060	for(i = 0; i < numberOfModels; i++)
1061	{
1062	a[i]=((ax[i] < cx[i]) ? ax[i] : cx[i]);
1063	b[i]=((ax[i] > cx[i]) ? ax[i] : cx[i]);
1064	x[i] = w[i] = v[i] = bx[i];
1065	fw[i] = fv[i] = fx[i] = fb[i];
1066	}
1067
1068	for(i = 0; i < numberOfModels; i++)
1069	{
1070	assert(a[i] >= lim_inf && a[i] <= lim_sup);
1071	assert(b[i] >= lim_inf && b[i] <= lim_sup);
1072	assert(x[i] >= lim_inf && x[i] <= lim_sup);
1073	assert(v[i] >= lim_inf && v[i] <= lim_sup);
1074	assert(w[i] >= lim_inf && w[i] <= lim_sup);
1075	}
1076
1077
1078
1079	for(iter = 1; iter <= ITMAX; iter++)
1080	{
1081	allConverged = TRUE;
1082
1083	for(i = 0; i < numberOfModels && allConverged; i++)
1084	allConverged = allConverged && converged[i];
1085
1086	if(allConverged)
1087	{
1088	rax_free(converged);
1089	rax_free(a);
1090	rax_free(b);
1091	rax_free(d);
1092	rax_free(etemp);
1093	rax_free(fu);
1094	rax_free(fv);
1095	rax_free(fw);
1096	rax_free(fx);
1097	rax_free(p);
1098	rax_free(q);
1099	rax_free(r);
1100	rax_free(tol1);
1101	rax_free(tol2);
1102	rax_free(u);
1103	rax_free(v);
1104	rax_free(w);
1105	rax_free(x);
1106	rax_free(xm);
1107	rax_free(e);
1108	return;
1109	}
1110
1111	for(i = 0; i < numberOfModels; i++)
1112	{
1113	if(!converged[i])
1114	{
1115	assert(a[i] >= lim_inf && a[i] <= lim_sup);
1116	assert(b[i] >= lim_inf && b[i] <= lim_sup);
1117	assert(x[i] >= lim_inf && x[i] <= lim_sup);
1118	assert(v[i] >= lim_inf && v[i] <= lim_sup);
1119	assert(w[i] >= lim_inf && w[i] <= lim_sup);
1120
1121	xm[i] = 0.5 * (a[i] + b[i]);
1122	tol2[i] = 2.0 * (tol1[i] = tol * fabs(x[i]) + BRENT_ZEPS);
1123
1124	if(fabs(x[i] - xm[i]) <= (tol2[i] - 0.5 * (b[i] - a[i])))
1125	{
1126	result[i] = -fx[i];
1127	xmin[i] = x[i];
1128	converged[i] = TRUE;
1129	}
1130	else
1131	{
1132	if(fabs(e[i]) > tol1[i])
1133	{
1134	r[i] = (x[i] - w[i]) * (fx[i] - fv[i]);
1135	q[i] = (x[i] - v[i]) * (fx[i] - fw[i]);
1136	p[i] = (x[i] - v[i]) * q[i] - (x[i] - w[i]) * r[i];
1137	q[i] = 2.0 * (q[i] - r[i]);
1138	if(q[i] > 0.0)
1139	p[i] = -p[i];
1140	q[i] = fabs(q[i]);
1141	etemp[i] = e[i];
1142	e[i] = d[i];
1143	if((fabs(p[i]) >= fabs(0.5 * q[i] * etemp[i])) \|\| (p[i] <= q[i] * (a[i]-x[i])) \|\| (p[i] >= q[i] * (b[i] - x[i])))
1144	d[i] = BRENT_CGOLD * (e[i] = (x[i] >= xm[i] ? a[i] - x[i] : b[i] - x[i]));
1145	else
1146	{
1147	d[i] = p[i] / q[i];
1148	u[i] = x[i] + d[i];
1149	if( u[i] - a[i] < tol2[i] \|\| b[i] - u[i] < tol2[i])
1150	d[i] = SIGN(tol1[i], xm[i] - x[i]);
1151	}
1152	}
1153	else
1154	{
1155	d[i] = BRENT_CGOLD * (e[i] = (x[i] >= xm[i] ? a[i] - x[i]: b[i] - x[i]));
1156	}
1157	u[i] = ((fabs(d[i]) >= tol1[i]) ? (x[i] + d[i]): (x[i] +SIGN(tol1[i], d[i])));
1158	}
1159
1160	if(!converged[i])
1161	assert(u[i] >= lim_inf && u[i] <= lim_sup);
1162	}
1163	}
1164
1165	evaluateChange(tr, rateNumber, u, fu, converged, whichFunction, numberOfModels, ll, tol);
1166
1167	for(i = 0; i < numberOfModels; i++)
1168	{
1169	if(!converged[i])
1170	{
1171	if(fu[i] <= fx[i])
1172	{
1173	if(u[i] >= x[i])
1174	a[i] = x[i];
1175	else
1176	b[i] = x[i];
1177
1178	SHFT(v[i],w[i],x[i],u[i]);
1179	SHFT(fv[i],fw[i],fx[i],fu[i]);
1180	}
1181	else
1182	{
1183	if(u[i] < x[i])
1184	a[i] = u[i];
1185	else
1186	b[i] = u[i];
1187
1188	if(fu[i] <= fw[i] \|\| w[i] == x[i])
1189	{
1190	v[i] = w[i];
1191	w[i] = u[i];
1192	fv[i] = fw[i];
1193	fw[i] = fu[i];
1194	}
1195	else
1196	{
1197	if(fu[i] <= fv[i] \|\| v[i] == x[i] \|\| v[i] == w[i])
1198	{
1199	v[i] = u[i];
1200	fv[i] = fu[i];
1201	}
1202	}
1203	}
1204
1205	assert(a[i] >= lim_inf && a[i] <= lim_sup);
1206	assert(b[i] >= lim_inf && b[i] <= lim_sup);
1207	assert(x[i] >= lim_inf && x[i] <= lim_sup);
1208	assert(v[i] >= lim_inf && v[i] <= lim_sup);
1209	assert(w[i] >= lim_inf && w[i] <= lim_sup);
1210	assert(u[i] >= lim_inf && u[i] <= lim_sup);
1211	}
1212	}
1213	}
1214
1215	rax_free(converged);
1216	rax_free(a);
1217	rax_free(b);
1218	rax_free(d);
1219	rax_free(etemp);
1220	rax_free(fu);
1221	rax_free(fv);
1222	rax_free(fw);
1223	rax_free(fx);
1224	rax_free(p);
1225	rax_free(q);
1226	rax_free(r);
1227	rax_free(tol1);
1228	rax_free(tol2);
1229	rax_free(u);
1230	rax_free(v);
1231	rax_free(w);
1232	rax_free(x);
1233	rax_free(xm);
1234	rax_free(e);
1235
1236	printf("\n. Too many iterations in BRENT !");
1237	assert(0);
1238	}
1239
1240
1241
1242	static int brakGeneric(double param, double ax, double bx, double cx, double fa, double fb,
1243	double *fc, double lim_inf, double lim_sup,
1244	int numberOfModels, int rateNumber, int whichFunction, tree tr, linkageList ll, double modelEpsilon)
1245	{
1246	double
1247	ulim = (double )rax_malloc(sizeof(double) * numberOfModels),
1248	u = (double )rax_malloc(sizeof(double) * numberOfModels),
1249	r = (double )rax_malloc(sizeof(double) * numberOfModels),
1250	q = (double )rax_malloc(sizeof(double) * numberOfModels),
1251	fu = (double )rax_malloc(sizeof(double) * numberOfModels),
1252	dum = (double )rax_malloc(sizeof(double) * numberOfModels),
1253	temp = (double )rax_malloc(sizeof(double) * numberOfModels);
1254
1255	int
1256	i,
1257	state = (int )rax_malloc(sizeof(int) * numberOfModels),
1258	endState = (int )rax_malloc(sizeof(int) * numberOfModels);
1259
1260	boolean converged = (boolean )rax_malloc(sizeof(boolean) * numberOfModels);
1261	boolean allConverged;
1262
1263	for(i = 0; i < numberOfModels; i++)
1264	converged[i] = FALSE;
1265
1266	for(i = 0; i < numberOfModels; i++)
1267	{
1268	state[i] = 0;
1269	endState[i] = 0;
1270
1271	u[i] = 0.0;
1272
1273	param[i] = ax[i];
1274
1275	if(param[i] > lim_sup)
1276	param[i] = ax[i] = lim_sup;
1277
1278	if(param[i] < lim_inf)
1279	param[i] = ax[i] = lim_inf;
1280
1281	assert(param[i] >= lim_inf && param[i] <= lim_sup);
1282	}
1283
1284
1285	evaluateChange(tr, rateNumber, param, fa, converged, whichFunction, numberOfModels, ll, modelEpsilon);
1286
1287
1288	for(i = 0; i < numberOfModels; i++)
1289	{
1290	param[i] = bx[i];
1291	if(param[i] > lim_sup)
1292	param[i] = bx[i] = lim_sup;
1293	if(param[i] < lim_inf)
1294	param[i] = bx[i] = lim_inf;
1295
1296	assert(param[i] >= lim_inf && param[i] <= lim_sup);
1297	}
1298
1299	evaluateChange(tr, rateNumber, param, fb, converged, whichFunction, numberOfModels, ll, modelEpsilon);
1300
1301	for(i = 0; i < numberOfModels; i++)
1302	{
1303	if (fb[i] > fa[i])
1304	{
1305	SHFT(dum[i],ax[i],bx[i],dum[i]);
1306	SHFT(dum[i],fa[i],fb[i],dum[i]);
1307	}
1308
1309	cx[i] = bx[i] + MNBRAK_GOLD * (bx[i] - ax[i]);
1310
1311	param[i] = cx[i];
1312
1313	if(param[i] > lim_sup)
1314	param[i] = cx[i] = lim_sup;
1315	if(param[i] < lim_inf)
1316	param[i] = cx[i] = lim_inf;
1317
1318	assert(param[i] >= lim_inf && param[i] <= lim_sup);
1319	}
1320
1321
1322	evaluateChange(tr, rateNumber, param, fc, converged, whichFunction, numberOfModels, ll, modelEpsilon);
1323
1324	while(1)
1325	{
1326	allConverged = TRUE;
1327
1328	for(i = 0; i < numberOfModels && allConverged; i++)
1329	allConverged = allConverged && converged[i];
1330
1331	if(allConverged)
1332	{
1333	for(i = 0; i < numberOfModels; i++)
1334	{
1335	if(ax[i] > lim_sup)
1336	ax[i] = lim_sup;
1337	if(ax[i] < lim_inf)
1338	ax[i] = lim_inf;
1339
1340	if(bx[i] > lim_sup)
1341	bx[i] = lim_sup;
1342	if(bx[i] < lim_inf)
1343	bx[i] = lim_inf;
1344
1345	if(cx[i] > lim_sup)
1346	cx[i] = lim_sup;
1347	if(cx[i] < lim_inf)
1348	cx[i] = lim_inf;
1349	}
1350
1351	rax_free(converged);
1352	rax_free(ulim);
1353	rax_free(u);
1354	rax_free(r);
1355	rax_free(q);
1356	rax_free(fu);
1357	rax_free(dum);
1358	rax_free(temp);
1359	rax_free(state);
1360	rax_free(endState);
1361	return 0;
1362
1363	}
1364
1365	for(i = 0; i < numberOfModels; i++)
1366	{
1367	if(!converged[i])
1368	{
1369	switch(state[i])
1370	{
1371	case 0:
1372	endState[i] = 0;
1373	if(!(fb[i] > fc[i]))
1374	converged[i] = TRUE;
1375	else
1376	{
1377
1378	if(ax[i] > lim_sup)
1379	ax[i] = lim_sup;
1380	if(ax[i] < lim_inf)
1381	ax[i] = lim_inf;
1382	if(bx[i] > lim_sup)
1383	bx[i] = lim_sup;
1384	if(bx[i] < lim_inf)
1385	bx[i] = lim_inf;
1386	if(cx[i] > lim_sup)
1387	cx[i] = lim_sup;
1388	if(cx[i] < lim_inf)
1389	cx[i] = lim_inf;
1390
1391	r[i]=(bx[i]-ax[i])*(fb[i]-fc[i]);
1392	q[i]=(bx[i]-cx[i])*(fb[i]-fa[i]);
1393	u[i]=(bx[i])-((bx[i]-cx[i])q[i]-(bx[i]-ax[i])r[i])/
1394	(2.0*SIGN(MAX(fabs(q[i]-r[i]),MNBRAK_TINY),q[i]-r[i]));
1395
1396	ulim[i]=(bx[i])+MNBRAK_GLIMIT*(cx[i]-bx[i]);
1397
1398	if(u[i] > lim_sup)
1399	u[i] = lim_sup;
1400	if(u[i] < lim_inf)
1401	u[i] = lim_inf;
1402	if(ulim[i] > lim_sup)
1403	ulim[i] = lim_sup;
1404	if(ulim[i] < lim_inf)
1405	ulim[i] = lim_inf;
1406
1407	if ((bx[i]-u[i])*(u[i]-cx[i]) > 0.0)
1408	{
1409	param[i] = u[i];
1410	if(param[i] > lim_sup)
1411	param[i] = u[i] = lim_sup;
1412	if(param[i] < lim_inf)
1413	param[i] = u[i] = lim_inf;
1414	endState[i] = 1;
1415	}
1416	else
1417	{
1418	if ((cx[i]-u[i])*(u[i]-ulim[i]) > 0.0)
1419	{
1420	param[i] = u[i];
1421	if(param[i] > lim_sup)
1422	param[i] = u[i] = lim_sup;
1423	if(param[i] < lim_inf)
1424	param[i] = u[i] = lim_inf;
1425	endState[i] = 2;
1426	}
1427	else
1428	{
1429	if ((u[i]-ulim[i])*(ulim[i]-cx[i]) >= 0.0)
1430	{
1431	u[i] = ulim[i];
1432	param[i] = u[i];
1433	if(param[i] > lim_sup)
1434	param[i] = u[i] = ulim[i] = lim_sup;
1435	if(param[i] < lim_inf)
1436	param[i] = u[i] = ulim[i] = lim_inf;
1437	endState[i] = 0;
1438	}
1439	else
1440	{
1441	u[i]=(cx[i])+MNBRAK_GOLD*(cx[i]-bx[i]);
1442	param[i] = u[i];
1443	endState[i] = 0;
1444	if(param[i] > lim_sup)
1445	param[i] = u[i] = lim_sup;
1446	if(param[i] < lim_inf)
1447	param[i] = u[i] = lim_inf;
1448	}
1449	}
1450	}
1451	}
1452	break;
1453	case 1:
1454	endState[i] = 0;
1455	break;
1456	case 2:
1457	endState[i] = 3;
1458	break;
1459	default:
1460	assert(0);
1461	}
1462	assert(param[i] >= lim_inf && param[i] <= lim_sup);
1463	}
1464	}
1465
1466	evaluateChange(tr, rateNumber, param, temp, converged, whichFunction, numberOfModels, ll, modelEpsilon);
1467
1468	for(i = 0; i < numberOfModels; i++)
1469	{
1470	if(!converged[i])
1471	{
1472	switch(endState[i])
1473	{
1474	case 0:
1475	fu[i] = temp[i];
1476	SHFT(ax[i],bx[i],cx[i],u[i]);
1477	SHFT(fa[i],fb[i],fc[i],fu[i]);
1478	state[i] = 0;
1479	break;
1480	case 1:
1481	fu[i] = temp[i];
1482	if (fu[i] < fc[i])
1483	{
1484	ax[i]=(bx[i]);
1485	bx[i]=u[i];
1486	fa[i]=(fb[i]);
1487	fb[i]=fu[i];
1488	converged[i] = TRUE;
1489	}
1490	else
1491	{
1492	if (fu[i] > fb[i])
1493	{
1494	assert(u[i] >= lim_inf && u[i] <= lim_sup);
1495	cx[i]=u[i];
1496	fc[i]=fu[i];
1497	converged[i] = TRUE;
1498	}
1499	else
1500	{
1501	u[i]=(cx[i])+MNBRAK_GOLD*(cx[i]-bx[i]);
1502	param[i] = u[i];
1503	if(param[i] > lim_sup) {param[i] = u[i] = lim_sup;}
1504	if(param[i] < lim_inf) {param[i] = u[i] = lim_inf;}
1505	state[i] = 1;
1506	}
1507	}
1508	break;
1509	case 2:
1510	fu[i] = temp[i];
1511	if (fu[i] < fc[i])
1512	{
1513	SHFT(bx[i],cx[i],u[i], cx[i]+MNBRAK_GOLD*(cx[i]-bx[i]));
1514	state[i] = 2;
1515	}
1516	else
1517	{
1518	state[i] = 0;
1519	SHFT(ax[i],bx[i],cx[i],u[i]);
1520	SHFT(fa[i],fb[i],fc[i],fu[i]);
1521	}
1522	break;
1523	case 3:
1524	SHFT(fb[i],fc[i],fu[i], temp[i]);
1525	SHFT(ax[i],bx[i],cx[i],u[i]);
1526	SHFT(fa[i],fb[i],fc[i],fu[i]);
1527	state[i] = 0;
1528	break;
1529	default:
1530	assert(0);
1531	}
1532	}
1533	}
1534	}
1535
1536
1537	assert(0);
1538	rax_free(converged);
1539	rax_free(ulim);
1540	rax_free(u);
1541	rax_free(r);
1542	rax_free(q);
1543	rax_free(fu);
1544	rax_free(dum);
1545	rax_free(temp);
1546	rax_free(state);
1547	rax_free(endState);
1548
1549
1550
1551	return(0);
1552	}
1553
1554
1555
1556
1557
1558
1559
1560
1561
1562	static void optInvar(tree tr, double modelEpsilon, linkageList ll)
1563	{
1564	int
1565	i,
1566	k,
1567	numberOfModels = ll->entries;
1568	double lim_inf = INVAR_MIN;
1569	double lim_sup = INVAR_MAX;
1570	double
1571	startLH = (double )rax_malloc(sizeof(double) * numberOfModels),
1572	startInvar = (double )rax_malloc(sizeof(double) * numberOfModels),
1573	endInvar = (double )rax_malloc(sizeof(double) * numberOfModels),
1574	_a = (double )rax_malloc(sizeof(double) * numberOfModels),
1575	_b = (double )rax_malloc(sizeof(double) * numberOfModels),
1576	_c = (double )rax_malloc(sizeof(double) * numberOfModels),
1577	_fa = (double )rax_malloc(sizeof(double) * numberOfModels),
1578	_fb = (double )rax_malloc(sizeof(double) * numberOfModels),
1579	_fc = (double )rax_malloc(sizeof(double) * numberOfModels),
1580	_param = (double )rax_malloc(sizeof(double) * numberOfModels),
1581	result = (double )rax_malloc(sizeof(double) * numberOfModels),
1582	_x = (double )rax_malloc(sizeof(double) * numberOfModels);
1583
1584	evaluateGenericInitrav(tr, tr->start);
1585
1586	#ifdef _USE_PTHREADS
1587	evaluateGeneric(tr, tr->start);
1588	/* to avoid transferring traversal info further on */
1589	#endif
1590
1591	#ifdef _DEBUG_MODEL_OPTIMIZATION
1592	double
1593	initialLH = tr->likelihood;
1594	#endif
1595
1596	for(i = 0; i < numberOfModels; i++)
1597	{
1598	assert(ll->ld[i].valid);
1599
1600	startInvar[i] = tr->partitionData[ll->ld[i].partitionList[0]].propInvariant;
1601	_a[i] = startInvar[i] + 0.1;
1602	_b[i] = startInvar[i] - 0.1;
1603	if(_b[i] < lim_inf)
1604	_b[i] = lim_inf;
1605
1606	startLH[i] = 0.0;
1607
1608	for(k = 0; k < ll->ld[i].partitions; k++)
1609	{
1610	startLH[i] += tr->perPartitionLH[ll->ld[i].partitionList[k]];
1611	/* TODO need to fix the initialization for this assertion not to fail */
1612	/* assert(tr->partitionData[ll->ld[i].partitionList[0]].propInvariant == tr->partitionData[ll->ld[i].partitionList[k]].propInvariant);*/
1613	}
1614	}
1615
1616	brakGeneric(_param, _a, _b, _c, _fa, _fb, _fc, lim_inf, lim_sup, numberOfModels, -1, INVAR_F, tr, ll, modelEpsilon);
1617	brentGeneric(_a, _b, _c, _fb, modelEpsilon, _x, result, numberOfModels, INVAR_F, -1, tr, ll, lim_inf, lim_sup);
1618
1619	for(i = 0; i < numberOfModels; i++)
1620	endInvar[i] = result[i];
1621
1622	for(i = 0; i < numberOfModels; i++)
1623	{
1624	if(startLH[i] > endInvar[i])
1625	{
1626	for(k = 0; k < ll->ld[i].partitions; k++)
1627	tr->partitionData[ll->ld[i].partitionList[k]].propInvariant = startInvar[i];
1628	}
1629	else
1630	{
1631	for(k = 0; k < ll->ld[i].partitions; k++)
1632	tr->partitionData[ll->ld[i].partitionList[k]].propInvariant = _x[i];
1633	}
1634	}
1635
1636	#ifdef _USE_PTHREADS
1637	masterBarrier(THREAD_COPY_INVAR, tr);
1638	#endif
1639
1640	#ifdef _DEBUG_MODEL_OPTIMIZATION
1641	evaluateGenericInitrav(tr, tr->start);
1642
1643	if(tr->likelihood < initialLH)
1644	printf("%f %f\n", tr->likelihood, initialLH);
1645	assert(tr->likelihood >= initialLH);
1646	#endif
1647
1648	rax_free(startLH);
1649	rax_free(startInvar);
1650	rax_free(endInvar);
1651	rax_free(result);
1652	rax_free(_a);
1653	rax_free(_b);
1654	rax_free(_c);
1655	rax_free(_fa);
1656	rax_free(_fb);
1657	rax_free(_fc);
1658	rax_free(_param);
1659	rax_free(_x);
1660
1661	}
1662
1663
1664
1665
1666
1667
1668	/**********************************************************************************************************/
1669	/* ALPHA PARAM ********************************************************************************************/
1670
1671
1672
1673
1674
1675
1676
1677	static void optAlpha(tree tr, double modelEpsilon, linkageList ll, int numberOfModels)
1678	{
1679	int
1680	pos,
1681	i;
1682
1683	double
1684	lim_inf = ALPHA_MIN,
1685	lim_sup = ALPHA_MAX,
1686	endLH = (double )rax_malloc(sizeof(double) * numberOfModels),
1687	startLH = (double )rax_malloc(sizeof(double) * numberOfModels),
1688	startAlpha = (double )rax_malloc(sizeof(double) * numberOfModels),
1689	endAlpha = (double )rax_malloc(sizeof(double) * numberOfModels),
1690	_a = (double )rax_malloc(sizeof(double) * numberOfModels),
1691	_b = (double )rax_malloc(sizeof(double) * numberOfModels),
1692	_c = (double )rax_malloc(sizeof(double) * numberOfModels),
1693	_fa = (double )rax_malloc(sizeof(double) * numberOfModels),
1694	_fb = (double )rax_malloc(sizeof(double) * numberOfModels),
1695	_fc = (double )rax_malloc(sizeof(double) * numberOfModels),
1696	_param = (double )rax_malloc(sizeof(double) * numberOfModels),
1697	result = (double )rax_malloc(sizeof(double) * numberOfModels),
1698	_x = (double )rax_malloc(sizeof(double) * numberOfModels);
1699
1700	evaluateGenericInitrav(tr, tr->start);
1701
1702	#ifdef _DEBUG_MODEL_OPTIMIZATION
1703	double
1704	initialLH = tr->likelihood;
1705	#endif
1706
1707	/*
1708	at this point here every worker has the traversal data it needs for the
1709	search, so we won't re-distribute it he he :-)
1710	*/
1711	for(i = 0, pos = 0; i < ll->entries; i++)
1712	{
1713	if(ll->ld[i].valid)
1714	{
1715	int
1716	index = ll->ld[i].partitionList[0];
1717
1718	assert(ll->ld[i].partitions == 1);
1719
1720	startAlpha[pos] = tr->partitionData[index].alpha;
1721
1722	_a[pos] = startAlpha[pos] + 0.1;
1723	_b[pos] = startAlpha[pos] - 0.1;
1724
1725	if(_a[pos] < lim_inf)
1726	_a[pos] = lim_inf;
1727
1728	if(_a[pos] > lim_sup)
1729	_a[pos] = lim_sup;
1730
1731	if(_b[pos] < lim_inf)
1732	_b[pos] = lim_inf;
1733
1734	if(_b[pos] > lim_sup)
1735	_b[pos] = lim_sup;
1736
1737	startLH[pos] = tr->perPartitionLH[index];
1738	endLH[pos] = unlikely;
1739
1740	pos++;
1741	}
1742	}
1743
1744	brakGeneric(_param, _a, _b, _c, _fa, _fb, _fc, lim_inf, lim_sup, numberOfModels, -1, ALPHA_F, tr, ll, modelEpsilon);
1745	brentGeneric(_a, _b, _c, _fb, modelEpsilon, _x, result, numberOfModels, ALPHA_F, -1, tr, ll, lim_inf, lim_sup);
1746
1747	for(i = 0; i < numberOfModels; i++)
1748	endLH[i] = result[i];
1749
1750	for(i = 0, pos = 0; i < ll->entries; i++)
1751	{
1752	if(ll->ld[i].valid)
1753	{
1754	int
1755	index = ll->ld[i].partitionList[0];
1756
1757	assert(ll->ld[i].partitions == 1);
1758
1759	if(startLH[pos] > endLH[pos])
1760	{
1761	tr->partitionData[index].alpha = startAlpha[pos];
1762	makeGammaCats(tr->partitionData[index].alpha, tr->partitionData[index].gammaRates, 4, tr->useGammaMedian);
1763	}
1764	else
1765	{
1766	tr->partitionData[index].alpha = _x[pos];
1767	makeGammaCats(tr->partitionData[index].alpha, tr->partitionData[index].gammaRates, 4, tr->useGammaMedian);
1768	}
1769
1770	pos++;
1771	}
1772	}
1773
1774	assert(pos == numberOfModels);
1775
1776	#ifdef _USE_PTHREADS
1777	masterBarrier(THREAD_COPY_ALPHA, tr);
1778	#endif
1779
1780
1781	#ifdef _DEBUG_MODEL_OPTIMIZATION
1782	evaluateGenericInitrav(tr, tr->start);
1783
1784	if(tr->likelihood < initialLH)
1785	printf("%f %f\n", tr->likelihood, initialLH);
1786	assert(tr->likelihood >= initialLH);
1787	#endif
1788
1789	rax_free(startLH);
1790	rax_free(endLH);
1791	rax_free(startAlpha);
1792	rax_free(endAlpha);
1793	rax_free(result);
1794	rax_free(_a);
1795	rax_free(_b);
1796	rax_free(_c);
1797	rax_free(_fa);
1798	rax_free(_fb);
1799	rax_free(_fc);
1800	rax_free(_param);
1801	rax_free(_x);
1802	}
1803
1804
1805	/*******************************************************************************************************************/
1806	/*****************RATES ****************************************************************************************/
1807
1808
1809
1810
1811
1812	static void optRate(tree tr, double modelEpsilon, linkageList ll, int numberOfModels, int states, int rateNumber, int numberOfRates)
1813	{
1814	int
1815	l,
1816	k,
1817	j,
1818	pos;
1819
1820	double
1821	lim_inf = RATE_MIN,
1822	lim_sup = RATE_MAX,
1823	startRates = (double )rax_malloc(sizeof(double) * numberOfRates * numberOfModels),
1824	startLH = (double )rax_malloc(sizeof(double) * numberOfModels),
1825	endLH = (double )rax_malloc(sizeof(double) * numberOfModels),
1826	_a = (double )rax_malloc(sizeof(double) * numberOfModels),
1827	_b = (double )rax_malloc(sizeof(double) * numberOfModels),
1828	_c = (double )rax_malloc(sizeof(double) * numberOfModels),
1829	_fa = (double )rax_malloc(sizeof(double) * numberOfModels),
1830	_fb = (double )rax_malloc(sizeof(double) * numberOfModels),
1831	_fc = (double )rax_malloc(sizeof(double) * numberOfModels),
1832	_param = (double )rax_malloc(sizeof(double) * numberOfModels),
1833	result = (double )rax_malloc(sizeof(double) * numberOfModels),
1834	_x = (double )rax_malloc(sizeof(double) * numberOfModels);
1835
1836	assert(states != -1);
1837
1838	evaluateGenericInitrav(tr, tr->start);
1839
1840	#ifdef _DEBUG_MODEL_OPTIMIZATION
1841	double
1842	initialLH = tr->likelihood;
1843	#endif
1844
1845	/*
1846	at this point here every worker has the traversal data it needs for the
1847	search
1848	*/
1849
1850	for(l = 0, pos = 0; l < ll->entries; l++)
1851	{
1852	if(ll->ld[l].valid)
1853	{
1854	endLH[pos] = unlikely;
1855	startLH[pos] = 0.0;
1856
1857	for(j = 0; j < ll->ld[l].partitions; j++)
1858	{
1859	int
1860	index = ll->ld[l].partitionList[j];
1861
1862	startLH[pos] += tr->perPartitionLH[index];
1863
1864	for(k = 0; k < numberOfRates; k++)
1865	startRates[pos * numberOfRates + k] = tr->partitionData[index].substRates[k];
1866	}
1867	pos++;
1868	}
1869	}
1870
1871	assert(pos == numberOfModels);
1872
1873	for(k = 0, pos = 0; k < ll->entries; k++)
1874	{
1875	if(ll->ld[k].valid)
1876	{
1877
1878	/* TODO partition List[0]? */
1879	int
1880	index = ll->ld[k].partitionList[0];
1881
1882
1883	_a[pos] = tr->partitionData[index].substRates[rateNumber] + 0.1;
1884	_b[pos] = tr->partitionData[index].substRates[rateNumber] - 0.1;
1885
1886	if(_a[pos] < lim_inf)
1887	_a[pos] = lim_inf;
1888
1889	if(_a[pos] > lim_sup)
1890	_a[pos] = lim_sup;
1891
1892	if(_b[pos] < lim_inf)
1893	_b[pos] = lim_inf;
1894
1895	if(_b[pos] > lim_sup)
1896	_b[pos] = lim_sup;
1897	pos++;
1898	}
1899	}
1900
1901	assert(pos == numberOfModels);
1902
1903	brakGeneric(_param, _a, _b, _c, _fa, _fb, _fc, lim_inf, lim_sup, numberOfModels, rateNumber, RATE_F, tr, ll, modelEpsilon);
1904
1905	for(k = 0; k < numberOfModels; k++)
1906	{
1907	assert(_a[k] >= lim_inf && _a[k] <= lim_sup);
1908	assert(_b[k] >= lim_inf && _b[k] <= lim_sup);
1909	assert(_c[k] >= lim_inf && _c[k] <= lim_sup);
1910	}
1911
1912	brentGeneric(_a, _b, _c, _fb, modelEpsilon, _x, result, numberOfModels, RATE_F, rateNumber, tr, ll, lim_inf, lim_sup);
1913
1914	for(k = 0; k < numberOfModels; k++)
1915	endLH[k] = result[k];
1916
1917	for(k = 0, pos = 0; k < ll->entries; k++)
1918	{
1919	if(ll->ld[k].valid)
1920	{
1921	if(startLH[pos] > endLH[pos])
1922	{
1923	for(j = 0; j < ll->ld[k].partitions; j++)
1924	{
1925	int
1926	index = ll->ld[k].partitionList[j];
1927
1928	setRateModel(tr, index, startRates[pos * numberOfRates + rateNumber], rateNumber);
1929	//tr->partitionData[index].substRates[rateNumber] = startRates[pos * numberOfRates + rateNumber];
1930	initReversibleGTR(tr, index);
1931	}
1932	}
1933	else
1934	{
1935	for(j = 0; j < ll->ld[k].partitions; j++)
1936	{
1937	int
1938	index = ll->ld[k].partitionList[j];
1939
1940	setRateModel(tr, index, _x[pos], rateNumber);
1941	//tr->partitionData[index].substRates[rateNumber] = _x[pos];
1942	initReversibleGTR(tr, index);
1943	}
1944	}
1945	pos++;
1946	}
1947	}
1948
1949	#ifdef _USE_PTHREADS
1950	masterBarrier(THREAD_COPY_RATES, tr);
1951	#endif
1952	assert(pos == numberOfModels);
1953
1954	rax_free(startLH);
1955	rax_free(endLH);
1956	rax_free(result);
1957	rax_free(_a);
1958	rax_free(_b);
1959	rax_free(_c);
1960	rax_free(_fa);
1961	rax_free(_fb);
1962	rax_free(_fc);
1963	rax_free(_param);
1964	rax_free(_x);
1965	rax_free(startRates);
1966
1967	#ifdef _DEBUG_MODEL_OPTIMIZATION
1968	evaluateGenericInitrav(tr, tr->start);
1969
1970	if(tr->likelihood < initialLH)
1971	printf("%f %f\n", tr->likelihood, initialLH);
1972	assert(tr->likelihood >= initialLH);
1973	#endif
1974
1975	}
1976
1977	static void optRates(tree tr, double modelEpsilon, linkageList ll, int numberOfModels, int states)
1978	{
1979	int
1980	rateNumber,
1981	numberOfRates = ((states * states - states) / 2) - 1;
1982
1983	for(rateNumber = 0; rateNumber < numberOfRates; rateNumber++)
1984	optRate(tr, modelEpsilon, ll, numberOfModels, states, rateNumber, numberOfRates);
1985	}
1986
1987	static boolean AAisGTR(tree *tr)
1988	{
1989	int
1990	i,
1991	count = 0;
1992
1993	for(i = 0; i < tr->NumberOfModels; i++)
1994	{
1995	if(tr->partitionData[i].dataType == AA_DATA)
1996	{
1997	count++;
1998	if(tr->partitionData[i].protModels != GTR)
1999	return FALSE;
2000	}
2001	}
2002
2003	if(count == 0)
2004	return FALSE;
2005
2006	return TRUE;
2007	}
2008
2009	static boolean AAisUnlinkedGTR(tree *tr)
2010	{
2011	int
2012	i,
2013	count = 0;
2014
2015	for(i = 0; i < tr->NumberOfModels; i++)
2016	{
2017	if(tr->partitionData[i].dataType == AA_DATA)
2018	{
2019	count++;
2020	if(tr->partitionData[i].protModels != GTR_UNLINKED)
2021	return FALSE;
2022	}
2023	}
2024
2025	if(count == 0)
2026	return FALSE;
2027
2028	return TRUE;
2029	}
2030
2031	static void optRatesGeneric(tree tr, double modelEpsilon, linkageList ll)
2032	{
2033	int
2034	i,
2035	dnaPartitions = 0,
2036	aaPartitionsLinked = 0,
2037	aaPartitionsUnlinked = 0,
2038	secondaryPartitions = 0,
2039	secondaryModel = -1,
2040	states = -1;
2041
2042	/* assumes homogeneous super-partitions, that either contain DNA or AA partitions !*/
2043	/* does not check whether AA are all linked */
2044
2045	/* first do DNA */
2046
2047	for(i = 0; i < ll->entries; i++)
2048	{
2049	switch(tr->partitionData[ll->ld[i].partitionList[0]].dataType)
2050	{
2051	case DNA_DATA:
2052	states = tr->partitionData[ll->ld[i].partitionList[0]].states;
2053	ll->ld[i].valid = TRUE;
2054	dnaPartitions++;
2055	break;
2056	case BINARY_DATA:
2057	case AA_DATA:
2058	case SECONDARY_DATA:
2059	case SECONDARY_DATA_6:
2060	case SECONDARY_DATA_7:
2061	case GENERIC_32:
2062	case GENERIC_64:
2063	ll->ld[i].valid = FALSE;
2064	break;
2065	default:
2066	assert(0);
2067	}
2068	}
2069
2070	if(dnaPartitions > 0)
2071	optRates(tr, modelEpsilon, ll, dnaPartitions, states);
2072
2073
2074	/* then SECONDARY */
2075
2076	for(i = 0; i < ll->entries; i++)
2077	{
2078	switch(tr->partitionData[ll->ld[i].partitionList[0]].dataType)
2079	{
2080	/* we only have one type of secondary data models in one analysis */
2081	case SECONDARY_DATA_6:
2082	states = tr->partitionData[ll->ld[i].partitionList[0]].states;
2083	secondaryModel = SECONDARY_DATA_6;
2084	ll->ld[i].valid = TRUE;
2085	secondaryPartitions++;
2086	break;
2087	case SECONDARY_DATA_7:
2088	states = tr->partitionData[ll->ld[i].partitionList[0]].states;
2089	secondaryModel = SECONDARY_DATA_7;
2090	ll->ld[i].valid = TRUE;
2091	secondaryPartitions++;
2092	break;
2093	case SECONDARY_DATA:
2094	states = tr->partitionData[ll->ld[i].partitionList[0]].states;
2095	secondaryModel = SECONDARY_DATA;
2096	ll->ld[i].valid = TRUE;
2097	secondaryPartitions++;
2098	break;
2099	case BINARY_DATA:
2100	case AA_DATA:
2101	case DNA_DATA:
2102	case GENERIC_32:
2103	case GENERIC_64:
2104	ll->ld[i].valid = FALSE;
2105	break;
2106	default:
2107	assert(0);
2108	}
2109	}
2110
2111
2112
2113	if(secondaryPartitions > 0)
2114	{
2115	assert(secondaryPartitions == 1);
2116
2117	switch(secondaryModel)
2118	{
2119	case SECONDARY_DATA:
2120	optRates(tr, modelEpsilon, ll, secondaryPartitions, states);
2121	break;
2122	case SECONDARY_DATA_6:
2123	optRates(tr, modelEpsilon, ll, secondaryPartitions, states);
2124	break;
2125	case SECONDARY_DATA_7:
2126	optRates(tr, modelEpsilon, ll, secondaryPartitions, states);
2127	break;
2128	default:
2129	assert(0);
2130	}
2131	}
2132
2133	/* then AA */
2134
2135	if(AAisGTR(tr))
2136	{
2137	for(i = 0; i < ll->entries; i++)
2138	{
2139	switch(tr->partitionData[ll->ld[i].partitionList[0]].dataType)
2140	{
2141	case AA_DATA:
2142	states = tr->partitionData[ll->ld[i].partitionList[0]].states;
2143	ll->ld[i].valid = TRUE;
2144	aaPartitionsLinked++;
2145	break;
2146	case DNA_DATA:
2147	case BINARY_DATA:
2148	case SECONDARY_DATA:
2149	case SECONDARY_DATA_6:
2150	case SECONDARY_DATA_7:
2151	ll->ld[i].valid = FALSE;
2152	break;
2153	default:
2154	assert(0);
2155	}
2156	}
2157
2158	assert(aaPartitionsLinked == 1);
2159
2160	optRates(tr, modelEpsilon, ll, aaPartitionsLinked, states);
2161	}
2162
2163	if(AAisUnlinkedGTR(tr))
2164	{
2165	aaPartitionsUnlinked = 0;
2166
2167	for(i = 0; i < ll->entries; i++)
2168	{
2169	switch(tr->partitionData[ll->ld[i].partitionList[0]].dataType)
2170	{
2171	case AA_DATA:
2172	states = tr->partitionData[ll->ld[i].partitionList[0]].states;
2173	ll->ld[i].valid = TRUE;
2174	aaPartitionsUnlinked++;
2175	break;
2176	case DNA_DATA:
2177	case BINARY_DATA:
2178	case SECONDARY_DATA:
2179	case SECONDARY_DATA_6:
2180	case SECONDARY_DATA_7:
2181	ll->ld[i].valid = FALSE;
2182	break;
2183	default:
2184	assert(0);
2185	}
2186	}
2187
2188	assert(aaPartitionsUnlinked >= 1);
2189
2190	optRates(tr, modelEpsilon, ll, aaPartitionsUnlinked, states);
2191	}
2192	/* then multi-state */
2193
2194	/*
2195	now here we have to be careful, because every multi-state partition can actually
2196	have a distinct number of states, so we will go to every multi-state partition separately,
2197	parallel efficiency for this will suck, but what the hell .....
2198	*/
2199
2200	if(tr->multiStateModel == GTR_MULTI_STATE)
2201	{
2202	for(i = 0; i < ll->entries; i++)
2203	{
2204	switch(tr->partitionData[ll->ld[i].partitionList[0]].dataType)
2205	{
2206	case GENERIC_32:
2207	{
2208	int k;
2209
2210	states = tr->partitionData[ll->ld[i].partitionList[0]].states;
2211
2212	ll->ld[i].valid = TRUE;
2213
2214	for(k = 0; k < ll->entries; k++)
2215	if(k != i)
2216	ll->ld[k].valid = FALSE;
2217
2218	optRates(tr, modelEpsilon, ll, 1, states);
2219	}
2220	break;
2221	case AA_DATA:
2222	case DNA_DATA:
2223	case BINARY_DATA:
2224	case SECONDARY_DATA:
2225	case SECONDARY_DATA_6:
2226	case SECONDARY_DATA_7:
2227	case GENERIC_64:
2228	break;
2229	default:
2230	assert(0);
2231	}
2232	}
2233	}
2234
2235	for(i = 0; i < ll->entries; i++)
2236	ll->ld[i].valid = TRUE;
2237	}
2238
2239	static void optLG4X_Rate(tree tr, double modelEpsilon, linkageList ll, int numberOfModels, int rateNumber)
2240	{
2241	int
2242	pos,
2243	i;
2244
2245	double
2246	lim_inf = LG4X_RATE_MIN,
2247	lim_sup = LG4X_RATE_MAX,
2248	startLH = (double )rax_malloc(sizeof(double) * numberOfModels),
2249	endLH = (double )rax_malloc(sizeof(double) * numberOfModels),
2250	startAlpha = (double )rax_malloc(sizeof(double) * numberOfModels * 4),
2251	startWeights = (double )rax_malloc(sizeof(double) * numberOfModels * 4),
2252	endAlpha = (double )rax_malloc(sizeof(double) * numberOfModels),
2253	_a = (double )rax_malloc(sizeof(double) * numberOfModels),
2254	_b = (double )rax_malloc(sizeof(double) * numberOfModels),
2255	_c = (double )rax_malloc(sizeof(double) * numberOfModels),
2256	_fa = (double )rax_malloc(sizeof(double) * numberOfModels),
2257	_fb = (double )rax_malloc(sizeof(double) * numberOfModels),
2258	_fc = (double )rax_malloc(sizeof(double) * numberOfModels),
2259	_param = (double )rax_malloc(sizeof(double) * numberOfModels),
2260	result = (double )rax_malloc(sizeof(double) * numberOfModels),
2261	_x = (double )rax_malloc(sizeof(double) * numberOfModels);
2262
2263	evaluateGenericInitrav(tr, tr->start);
2264	// printf("Enter lg4x rates: %f\n", tr->likelihood);
2265
2266	#ifdef _DEBUG_MODEL_OPTIMIZATION
2267	double
2268	initialLH = tr->likelihood;
2269	#endif
2270
2271	assert(rateNumber >= 0 && rateNumber < 4);
2272
2273	/*
2274	at this point here every worker has the traversal data it needs for the
2275	search, so we won't re-distribute it he he :-)
2276	*/
2277
2278	for(i = 0, pos = 0; i < ll->entries; i++)
2279	{
2280	if(ll->ld[i].valid)
2281	{
2282	int
2283	index = ll->ld[i].partitionList[0];
2284
2285	assert(ll->ld[i].partitions == 1);
2286
2287	memcpy(&startAlpha[pos * 4], tr->partitionData[index].gammaRates, 4 * sizeof(double));
2288	memcpy(&startWeights[pos * 4], tr->partitionData[index].weights, 4 * sizeof(double));
2289
2290	_a[pos] = startAlpha[pos * 4 + rateNumber] + 0.1;
2291	_b[pos] = startAlpha[pos * 4 + rateNumber] - 0.1;
2292
2293	if(_a[pos] < lim_inf)
2294	_a[pos] = lim_inf;
2295
2296	if(_a[pos] > lim_sup)
2297	_a[pos] = lim_sup;
2298
2299	if(_b[pos] < lim_inf)
2300	_b[pos] = lim_inf;
2301
2302	if(_b[pos] > lim_sup)
2303	_b[pos] = lim_sup;
2304
2305	startLH[pos] = tr->perPartitionLH[index];
2306	endLH[pos] = unlikely;
2307
2308	pos++;
2309	}
2310	}
2311
2312	assert(pos == numberOfModels);
2313
2314	brakGeneric(_param, _a, _b, _c, _fa, _fb, _fc, lim_inf, lim_sup, numberOfModels, rateNumber, LXRATE_F, tr, ll, modelEpsilon);
2315	brentGeneric(_a, _b, _c, _fb, modelEpsilon, _x, result, numberOfModels, LXRATE_F, rateNumber, tr, ll, lim_inf, lim_sup);
2316
2317	/* now calculate the likelihood with the optimized rate */
2318
2319	for(i = 0, pos = 0; i < ll->entries; i++)
2320	{
2321	if(ll->ld[i].valid)
2322	{
2323	int
2324	index = ll->ld[i].partitionList[0];
2325
2326	assert(ll->ld[i].partitions == 1);
2327
2328	tr->executeModel[index] = TRUE;
2329
2330	tr->partitionData[index].gammaRates[rateNumber] = _x[pos];
2331
2332	pos++;
2333	}
2334	}
2335
2336	#ifdef _USE_PTHREADS
2337	masterBarrier(THREAD_COPY_LG4X_RATES, tr);
2338	#endif
2339
2340	evaluateGenericInitrav(tr, tr->start);
2341
2342	for(i = 0, pos = 0; i < ll->entries; i++)
2343	{
2344	if(ll->ld[i].valid)
2345	{
2346	int
2347	index = ll->ld[i].partitionList[0];
2348
2349	assert(ll->ld[i].partitions == 1);
2350
2351	endLH[pos] = tr->perPartitionLH[index];
2352
2353	if(startLH[pos] > endLH[pos])
2354	{
2355	memcpy(tr->partitionData[index].weights, &startWeights[pos * 4], sizeof(double) * 4);
2356	memcpy(tr->partitionData[index].gammaRates, &startAlpha[pos * 4], sizeof(double) * 4);
2357	}
2358
2359	pos++;
2360	}
2361	}
2362
2363	assert(pos == numberOfModels);
2364
2365	#ifdef _USE_PTHREADS
2366	masterBarrier(THREAD_COPY_LG4X_RATES, tr);
2367	#endif
2368
2369
2370	#ifdef _DEBUG_MODEL_OPTIMIZATION
2371	evaluateGenericInitrav(tr, tr->start);
2372
2373	if(tr->likelihood < initialLH)
2374	printf("%f %f\n", tr->likelihood, initialLH);
2375	assert(tr->likelihood >= initialLH);
2376	#endif
2377
2378	rax_free(startLH);
2379	rax_free(endLH);
2380	rax_free(startWeights);
2381	rax_free(startAlpha);
2382	rax_free(endAlpha);
2383	rax_free(result);
2384	rax_free(_a);
2385	rax_free(_b);
2386	rax_free(_c);
2387	rax_free(_fa);
2388	rax_free(_fb);
2389	rax_free(_fc);
2390	rax_free(_param);
2391	rax_free(_x);
2392	}
2393
2394	/*** WARNING: the re-scaling of the branch lengths will only work if
2395	it is done after the rate optimization that modifies fracchange
2396	***/
2397
2398
2399	static void optLG4X(tree tr, double modelEpsilon, linkageList ll, int numberOfModels)
2400	{
2401	int
2402	i;
2403
2404	double
2405	lg4xScaler,
2406	lg4xScalers = (double )rax_calloc(tr->NumberOfModels, sizeof(double)),
2407	modelWeights = (double )rax_calloc(tr->NumberOfModels, sizeof(double)),
2408	wgtsum = 0.0;
2409
2410	for(i = 0; i < 4; i++)
2411	optLG4X_Rate(tr, modelEpsilon, ll, numberOfModels, i);
2412
2413	for(i = 0; i < tr->NumberOfModels; i++)
2414	lg4xScalers[i] = 1.0;
2415
2416	for(i = 0; i < ll->entries; i++)
2417	{
2418	if(ll->ld[i].valid)
2419	{
2420	int
2421	j,
2422	index = ll->ld[i].partitionList[0];
2423
2424	double
2425	averageRate = 0.0;
2426
2427	assert(ll->ld[i].partitions == 1);
2428
2429	for(j = 0; j < 4; j++)
2430	averageRate += tr->partitionData[index].gammaRates[j];
2431
2432	averageRate /= 4.0;
2433
2434	lg4xScalers[index] = averageRate;
2435	}
2436	}
2437
2438	if(tr->NumberOfModels > 1)
2439	{
2440	for(i = 0; i < tr->NumberOfModels; i++)
2441	tr->fracchanges[i] = tr->rawFracchanges[i] * (1.0 / lg4xScalers[i]);
2442	}
2443
2444	for(i = 0; i < tr->cdta->endsite; i++)
2445	{
2446	modelWeights[tr->model[i]] += (double)tr->cdta->aliaswgt[i];
2447	wgtsum += (double)tr->cdta->aliaswgt[i];
2448	}
2449
2450	lg4xScaler = 0.0;
2451
2452	for(i = 0; i < tr->NumberOfModels; i++)
2453	{
2454	double
2455	fraction = modelWeights[i] / wgtsum;
2456
2457	lg4xScaler += (fraction * lg4xScalers[i]);
2458	}
2459
2460	tr->fracchange = tr->rawFracchange * (1.0 / lg4xScaler);
2461
2462	rax_free(lg4xScalers);
2463	rax_free(modelWeights);
2464	}
2465
2466	static void optAlphasGeneric(tree tr, double modelEpsilon, linkageList ll)
2467	{
2468	int
2469	i,
2470	non_LG4X_Partitions = 0,
2471	LG4X_Partitions = 0;
2472
2473	/* assumes homogeneous super-partitions, that either contain DNA or AA partitions !*/
2474	/* does not check whether AA are all linked */
2475
2476	/* first do non-LG4X partitions */
2477
2478	for(i = 0; i < ll->entries; i++)
2479	{
2480	switch(tr->partitionData[ll->ld[i].partitionList[0]].dataType)
2481	{
2482	case DNA_DATA:
2483	case BINARY_DATA:
2484	case SECONDARY_DATA:
2485	case SECONDARY_DATA_6:
2486	case SECONDARY_DATA_7:
2487	case GENERIC_32:
2488	case GENERIC_64:
2489	ll->ld[i].valid = TRUE;
2490	non_LG4X_Partitions++;
2491	break;
2492	case AA_DATA:
2493	if(tr->partitionData[ll->ld[i].partitionList[0]].protModels == LG4X)
2494	{
2495	LG4X_Partitions++;
2496	ll->ld[i].valid = FALSE;
2497	}
2498	else
2499	{
2500	ll->ld[i].valid = TRUE;
2501	non_LG4X_Partitions++;
2502	}
2503	break;
2504	default:
2505	assert(0);
2506	}
2507	}
2508
2509
2510
2511	if(non_LG4X_Partitions > 0)
2512	optAlpha(tr, modelEpsilon, ll, non_LG4X_Partitions);
2513
2514
2515
2516	/* then LG4x partitions */
2517
2518	for(i = 0; i < ll->entries; i++)
2519	{
2520	switch(tr->partitionData[ll->ld[i].partitionList[0]].dataType)
2521	{
2522	case DNA_DATA:
2523	case BINARY_DATA:
2524	case SECONDARY_DATA:
2525	case SECONDARY_DATA_6:
2526	case SECONDARY_DATA_7:
2527	case GENERIC_32:
2528	case GENERIC_64:
2529	ll->ld[i].valid = FALSE;
2530	break;
2531	case AA_DATA:
2532	if(tr->partitionData[ll->ld[i].partitionList[0]].protModels == LG4X)
2533	ll->ld[i].valid = TRUE;
2534	else
2535	ll->ld[i].valid = FALSE;
2536	break;
2537	default:
2538	assert(0);
2539	}
2540	}
2541
2542
2543
2544	if(LG4X_Partitions > 0)
2545	optLG4X(tr, modelEpsilon, ll, LG4X_Partitions);
2546
2547
2548
2549
2550	for(i = 0; i < ll->entries; i++)
2551	ll->ld[i].valid = TRUE;
2552	}
2553
2554
2555	/*******************FUNCTIONS FOR APPROXIMATE MODEL OPTIMIZATION *************************************/
2556
2557
2558
2559
2560
2561
2562	static int catCompare(const void p1, const void p2)
2563	{
2564	rateCategorize rc1 = (rateCategorize )p1;
2565	rateCategorize rc2 = (rateCategorize )p2;
2566
2567	double i = rc1->accumulatedSiteLikelihood;
2568	double j = rc2->accumulatedSiteLikelihood;
2569
2570	if (i > j)
2571	return (1);
2572	if (i < j)
2573	return (-1);
2574	return (0);
2575	}
2576
2577
2578
2579	static void categorizePartition(tree tr, rateCategorize rc, int model, int lower, int upper)
2580	{
2581	int
2582	zeroCounter,
2583	i,
2584	k;
2585
2586	double
2587	diff,
2588	min;
2589
2590	for (i = lower, zeroCounter = 0; i < upper; i++, zeroCounter++)
2591	{
2592	double
2593	temp = tr->cdta->patrat[i];
2594
2595	int
2596	found = 0;
2597
2598	for(k = 0; k < tr->partitionData[model].numberOfCategories; k++)
2599	{
2600	if(temp == rc[k].rate \|\| (fabs(temp - rc[k].rate) < 0.001))
2601	{
2602	found = 1;
2603	tr->cdta->rateCategory[i] = k;
2604	break;
2605	}
2606	}
2607
2608	if(!found)
2609	{
2610	min = fabs(temp - rc[0].rate);
2611	tr->cdta->rateCategory[i] = 0;
2612
2613	for(k = 1; k < tr->partitionData[model].numberOfCategories; k++)
2614	{
2615	diff = fabs(temp - rc[k].rate);
2616
2617	if(diff < min)
2618	{
2619	min = diff;
2620	tr->cdta->rateCategory[i] = k;
2621	}
2622	}
2623	}
2624	}
2625
2626	for(k = 0; k < tr->partitionData[model].numberOfCategories; k++)
2627	tr->partitionData[model].unscaled_perSiteRates[k] = rc[k].rate;
2628	}
2629
2630
2631	#ifdef _USE_PTHREADS
2632
2633	void optRateCatPthreads(tree tr, double lower_spacing, double upper_spacing, double lhs, int n, int tid)
2634	{
2635	int
2636	model;
2637
2638	size_t
2639	localIndex,
2640	i;
2641
2642	for(model = 0; model < tr->NumberOfModels; model++)
2643	{
2644	for(i = tr->partitionData[model].lower, localIndex = 0; i < tr->partitionData[model].upper; i++)
2645	{
2646	if(i % ((size_t)n) == ((size_t)tid))
2647	{
2648
2649	double initialRate, initialLikelihood,
2650	leftLH, rightLH, leftRate, rightRate, v;
2651	const double epsilon = 0.00001;
2652	int k;
2653
2654	tr->cdta->patrat[i] = tr->cdta->patratStored[i];
2655	initialRate = tr->cdta->patrat[i];
2656
2657	initialLikelihood = evaluatePartialGeneric(tr, localIndex, initialRate, model); /* i is real i ??? */
2658
2659
2660	leftLH = rightLH = initialLikelihood;
2661	leftRate = rightRate = initialRate;
2662
2663	k = 1;
2664
2665	while((initialRate - k * lower_spacing > 0.0001) &&
2666	((v = evaluatePartialGeneric(tr, localIndex, initialRate - k * lower_spacing, model))
2667	> leftLH) &&
2668	(fabs(leftLH - v) > epsilon))
2669	{
2670	#ifndef WIN32
2671	if(isnan(v))
2672	assert(0);
2673	#endif
2674
2675	leftLH = v;
2676	leftRate = initialRate - k * lower_spacing;
2677	k++;
2678	}
2679
2680	k = 1;
2681
2682	while(((v = evaluatePartialGeneric(tr, localIndex, initialRate + k * upper_spacing, model)) > rightLH) &&
2683	(fabs(rightLH - v) > epsilon))
2684	{
2685	#ifndef WIN32
2686	if(isnan(v))
2687	assert(0);
2688	#endif
2689	rightLH = v;
2690	rightRate = initialRate + k * upper_spacing;
2691	k++;
2692	}
2693
2694	if(rightLH > initialLikelihood \|\| leftLH > initialLikelihood)
2695	{
2696	if(rightLH > leftLH)
2697	{
2698	tr->cdta->patrat[i] = rightRate;
2699	lhs[i] = rightLH;
2700	}
2701	else
2702	{
2703	tr->cdta->patrat[i] = leftRate;
2704	lhs[i] = leftLH;
2705	}
2706	}
2707	else
2708	lhs[i] = initialLikelihood;
2709
2710	tr->cdta->patratStored[i] = tr->cdta->patrat[i];
2711	localIndex++;
2712	}
2713	}
2714	assert(localIndex == tr->partitionData[model].width);
2715	}
2716	}
2717
2718
2719
2720	#else
2721
2722
2723	static void optRateCatModel(tree tr, int model, double lower_spacing, double upper_spacing, double lhs)
2724	{
2725	int lower = tr->partitionData[model].lower;
2726	int upper = tr->partitionData[model].upper;
2727	int i;
2728	for(i = lower; i < upper; i++)
2729	{
2730	double initialRate, initialLikelihood,
2731	leftLH, rightLH, leftRate, rightRate, v;
2732	const double epsilon = 0.00001;
2733	int k;
2734
2735	tr->cdta->patrat[i] = tr->cdta->patratStored[i];
2736	initialRate = tr->cdta->patrat[i];
2737
2738	initialLikelihood = evaluatePartialGeneric(tr, i, initialRate, model);
2739
2740
2741	leftLH = rightLH = initialLikelihood;
2742	leftRate = rightRate = initialRate;
2743
2744	k = 1;
2745
2746	while((initialRate - k * lower_spacing > 0.0001) &&
2747	((v = evaluatePartialGeneric(tr, i, initialRate - k * lower_spacing, model))
2748	> leftLH) &&
2749	(fabs(leftLH - v) > epsilon))
2750	{
2751	#ifndef WIN32
2752	if(isnan(v))
2753	assert(0);
2754	#endif
2755
2756	leftLH = v;
2757	leftRate = initialRate - k * lower_spacing;
2758	k++;
2759	}
2760
2761	k = 1;
2762
2763	while(((v = evaluatePartialGeneric(tr, i, initialRate + k * upper_spacing, model)) > rightLH) &&
2764	(fabs(rightLH - v) > epsilon))
2765	{
2766	#ifndef WIN32
2767	if(isnan(v))
2768	assert(0);
2769	#endif
2770	rightLH = v;
2771	rightRate = initialRate + k * upper_spacing;
2772	k++;
2773	}
2774
2775	if(rightLH > initialLikelihood \|\| leftLH > initialLikelihood)
2776	{
2777	if(rightLH > leftLH)
2778	{
2779	tr->cdta->patrat[i] = rightRate;
2780	lhs[i] = rightLH;
2781	}
2782	else
2783	{
2784	tr->cdta->patrat[i] = leftRate;
2785	lhs[i] = leftLH;
2786	}
2787	}
2788	else
2789	lhs[i] = initialLikelihood;
2790
2791	tr->cdta->patratStored[i] = tr->cdta->patrat[i];
2792	}
2793
2794	}
2795
2796
2797	#endif
2798
2799
2800
2801	/*
2802	set scaleRates to FALSE everywhere such that
2803	per-site rates are not scaled to obtain an overall mean rate
2804	of 1.0
2805	*/
2806
2807	void updatePerSiteRates(tree *tr, boolean scaleRates)
2808	{
2809	int
2810	i,
2811	model;
2812
2813	if(tr->multiBranch)
2814	{
2815	for(model = 0; model < tr->NumberOfModels; model++)
2816	{
2817	int
2818	lower = tr->partitionData[model].lower,
2819	upper = tr->partitionData[model].upper;
2820
2821	if(scaleRates)
2822	{
2823	double
2824	scaler = 0.0,
2825	accRat = 0.0;
2826
2827	int
2828	accWgt = 0;
2829
2830	for(i = lower; i < upper; i++)
2831	{
2832	int
2833	w = tr->cdta->aliaswgt[i];
2834
2835	double
2836	rate = tr->partitionData[model].unscaled_perSiteRates[tr->cdta->rateCategory[i]];
2837
2838	assert(0 <= tr->cdta->rateCategory[i] && tr->cdta->rateCategory[i] < tr->maxCategories);
2839
2840	accWgt += w;
2841
2842	accRat += (w * rate);
2843	}
2844
2845	accRat /= ((double)accWgt);
2846
2847	scaler = 1.0 / ((double)accRat);
2848
2849	for(i = 0; i < tr->partitionData[model].numberOfCategories; i++)
2850	tr->partitionData[model].perSiteRates[i] = scaler * tr->partitionData[model].unscaled_perSiteRates[i];
2851
2852	accRat = 0.0;
2853
2854	for(i = lower; i < upper; i++)
2855	{
2856	int
2857	w = tr->cdta->aliaswgt[i];
2858
2859	double
2860	rate = tr->partitionData[model].perSiteRates[tr->cdta->rateCategory[i]];
2861
2862	assert(0 <= tr->cdta->rateCategory[i] && tr->cdta->rateCategory[i] < tr->maxCategories);
2863
2864	accRat += (w * rate);
2865	}
2866
2867	accRat /= ((double)accWgt);
2868
2869	if(!(ABS(1.0 - accRat) < 1.0E-5))
2870	printBothOpen("An assertion will fail: accumulated rate categories: %1.40f\n", accRat);
2871
2872	assert(ABS(1.0 - accRat) < 1.0E-5);
2873	}
2874
2875
2876
2877
2878	}
2879	}
2880	else
2881	{
2882	int
2883	accWgt = 0;
2884
2885	double
2886	scaler = 0.0,
2887	accRat = 0.0;
2888
2889	if(scaleRates)
2890	{
2891	for(model = 0, accRat = 0.0, accWgt = 0; model < tr->NumberOfModels; model++)
2892	{
2893	int
2894	localCount = 0,
2895	lower = tr->partitionData[model].lower,
2896	upper = tr->partitionData[model].upper;
2897
2898	for(i = lower, localCount = 0; i < upper; i++, localCount++)
2899	{
2900	int
2901	w = tr->cdta->aliaswgt[i];
2902
2903	double
2904	rate = tr->partitionData[model].unscaled_perSiteRates[tr->cdta->rateCategory[i]];
2905
2906	assert(0 <= tr->cdta->rateCategory[i] && tr->cdta->rateCategory[i] < tr->maxCategories);
2907
2908	accWgt += w;
2909
2910	accRat += (w * rate);
2911	}
2912	}
2913
2914
2915
2916	accRat /= ((double)accWgt);
2917
2918	scaler = 1.0 / ((double)accRat);
2919
2920	for(model = 0; model < tr->NumberOfModels; model++)
2921	{
2922	for(i = 0; i < tr->partitionData[model].numberOfCategories; i++)
2923	tr->partitionData[model].perSiteRates[i] = scaler * tr->partitionData[model].unscaled_perSiteRates[i];
2924	}
2925
2926	for(model = 0, accRat = 0.0; model < tr->NumberOfModels; model++)
2927	{
2928	int
2929	localCount = 0,
2930	lower = tr->partitionData[model].lower,
2931	upper = tr->partitionData[model].upper;
2932
2933	for(i = lower, localCount = 0; i < upper; i++, localCount++)
2934	{
2935	int
2936	w = tr->cdta->aliaswgt[i];
2937
2938	double
2939	rate = tr->partitionData[model].perSiteRates[tr->cdta->rateCategory[i]];
2940
2941	assert(0 <= tr->cdta->rateCategory[i] && tr->cdta->rateCategory[i] < tr->maxCategories);
2942
2943	accRat += (w * rate);
2944	}
2945	}
2946
2947	accRat /= ((double)accWgt);
2948
2949	if(!(ABS(1.0 - accRat) < 1.0E-5))
2950	printBothOpen("An assertion will fail: accumulated rate categories: %1.40f\n", accRat);
2951
2952	assert(ABS(1.0 - accRat) < 1.0E-5);
2953	}
2954
2955
2956
2957	}
2958
2959
2960	#ifdef _USE_PTHREADS
2961	masterBarrier(THREAD_COPY_RATE_CATS, tr);
2962	#endif
2963	}
2964
2965	static void optimizeRateCategories(tree *tr, int _maxCategories)
2966	{
2967	assert(_maxCategories > 0);
2968
2969	if(_maxCategories > 1)
2970	{
2971	double
2972	temp,
2973	lower_spacing,
2974	upper_spacing,
2975	initialLH = tr->likelihood,
2976	ratStored = (double )rax_malloc(sizeof(double) * tr->cdta->endsite),
2977	lhs = (double )rax_malloc(sizeof(double) * tr->cdta->endsite),
2978	oldCategorizedRates = (double )rax_malloc(sizeof(double ) tr->NumberOfModels),
2979	oldUnscaledCategorizedRates = (double )rax_malloc(sizeof(double ) tr->NumberOfModels);
2980
2981	int
2982	i,
2983	k,
2984	maxCategories = _maxCategories,
2985	oldCategory = (int )rax_malloc(sizeof(int) * tr->cdta->endsite),
2986	model,
2987	oldNumbers = (int )rax_malloc(sizeof(int) * tr->NumberOfModels);
2988
2989	assert(isTip(tr->start->number, tr->rdta->numsp));
2990
2991	determineFullTraversal(tr->start, tr);
2992
2993	if(optimizeRateCategoryInvocations == 1)
2994	{
2995	lower_spacing = 0.5 / ((double)optimizeRateCategoryInvocations);
2996	upper_spacing = 1.0 / ((double)optimizeRateCategoryInvocations);
2997	}
2998	else
2999	{
3000	lower_spacing = 0.05 / ((double)optimizeRateCategoryInvocations);
3001	upper_spacing = 0.1 / ((double)optimizeRateCategoryInvocations);
3002	}
3003
3004	if(lower_spacing < 0.001)
3005	lower_spacing = 0.001;
3006
3007	if(upper_spacing < 0.001)
3008	upper_spacing = 0.001;
3009
3010	optimizeRateCategoryInvocations++;
3011
3012	memcpy(oldCategory, tr->cdta->rateCategory, sizeof(int) * tr->cdta->endsite);
3013	memcpy(ratStored, tr->cdta->patratStored, sizeof(double) * tr->cdta->endsite);
3014
3015	for(model = 0; model < tr->NumberOfModels; model++)
3016	{
3017	oldNumbers[model] = tr->partitionData[model].numberOfCategories;
3018
3019	oldCategorizedRates[model] = (double )rax_malloc(sizeof(double) tr->maxCategories);
3020	oldUnscaledCategorizedRates[model] = (double )rax_malloc(sizeof(double) tr->maxCategories);
3021
3022	memcpy(oldCategorizedRates[model], tr->partitionData[model].perSiteRates, tr->maxCategories * sizeof(double));
3023	memcpy(oldUnscaledCategorizedRates[model], tr->partitionData[model].unscaled_perSiteRates, tr->maxCategories * sizeof(double));
3024	}
3025
3026	#ifdef _USE_PTHREADS
3027	tr->lhs = lhs;
3028	tr->lower_spacing = lower_spacing;
3029	tr->upper_spacing = upper_spacing;
3030	masterBarrier(THREAD_RATE_CATS, tr);
3031	#else
3032	for(model = 0; model < tr->NumberOfModels; model++)
3033	optRateCatModel(tr, model, lower_spacing, upper_spacing, lhs);
3034	#endif
3035
3036	for(model = 0; model < tr->NumberOfModels; model++)
3037	{
3038	int
3039	where = 1,
3040	found = 0,
3041	width = tr->partitionData[model].upper - tr->partitionData[model].lower,
3042	upper = tr->partitionData[model].upper,
3043	lower = tr->partitionData[model].lower;
3044
3045	rateCategorize
3046	rc = (rateCategorize )rax_malloc(sizeof(rateCategorize) * width);
3047
3048	for (i = 0; i < width; i++)
3049	{
3050	rc[i].accumulatedSiteLikelihood = 0.0;
3051	rc[i].rate = 0.0;
3052	}
3053
3054	rc[0].accumulatedSiteLikelihood = lhs[lower];
3055	rc[0].rate = tr->cdta->patrat[lower];
3056
3057	tr->cdta->rateCategory[lower] = 0;
3058
3059	for (i = lower + 1; i < upper; i++)
3060	{
3061	temp = tr->cdta->patrat[i];
3062	found = 0;
3063
3064	for(k = 0; k < where; k++)
3065	{
3066	if(temp == rc[k].rate \|\| (fabs(temp - rc[k].rate) < 0.001))
3067	{
3068	found = 1;
3069	rc[k].accumulatedSiteLikelihood += lhs[i];
3070	break;
3071	}
3072	}
3073
3074	if(!found)
3075	{
3076	rc[where].rate = temp;
3077	rc[where].accumulatedSiteLikelihood += lhs[i];
3078	where++;
3079	}
3080	}
3081
3082	qsort(rc, where, sizeof(rateCategorize), catCompare);
3083
3084	if(where < maxCategories)
3085	{
3086	tr->partitionData[model].numberOfCategories = where;
3087	categorizePartition(tr, rc, model, lower, upper);
3088	}
3089	else
3090	{
3091	tr->partitionData[model].numberOfCategories = maxCategories;
3092	categorizePartition(tr, rc, model, lower, upper);
3093	}
3094
3095	rax_free(rc);
3096	}
3097
3098	updatePerSiteRates(tr, TRUE);
3099
3100	evaluateGenericInitrav(tr, tr->start);
3101
3102	if(tr->likelihood < initialLH)
3103	{
3104	for(model = 0; model < tr->NumberOfModels; model++)
3105	{
3106	tr->partitionData[model].numberOfCategories = oldNumbers[model];
3107	memcpy(tr->partitionData[model].perSiteRates, oldCategorizedRates[model], tr->maxCategories * sizeof(double));
3108	memcpy(tr->partitionData[model].unscaled_perSiteRates, oldUnscaledCategorizedRates[model], tr->maxCategories * sizeof(double));
3109	}
3110
3111	memcpy(tr->cdta->patratStored, ratStored, sizeof(double) * tr->cdta->endsite);
3112	memcpy(tr->cdta->rateCategory, oldCategory, sizeof(int) * tr->cdta->endsite);
3113
3114	updatePerSiteRates(tr, FALSE);
3115
3116	evaluateGenericInitrav(tr, tr->start);
3117
3118	/* printf("REVERT: %1.40f %1.40f\n", initialLH, tr->likelihood); */
3119
3120	assert(initialLH == tr->likelihood);
3121	}
3122
3123	for(model = 0; model < tr->NumberOfModels; model++)
3124	{
3125	rax_free(oldCategorizedRates[model]);
3126	rax_free(oldUnscaledCategorizedRates[model]);
3127	}
3128
3129	rax_free(oldCategorizedRates);
3130	rax_free(oldUnscaledCategorizedRates);
3131	rax_free(oldCategory);
3132	rax_free(ratStored);
3133	rax_free(lhs);
3134	rax_free(oldNumbers);
3135	}
3136	}
3137
3138
3139
3140
3141
3142
3143
3144	/*****************************************************************************************************/
3145
3146	void resetBranches(tree *tr)
3147	{
3148	nodeptr p, q;
3149	int nodes, i;
3150
3151	nodes = tr->mxtips + 3 * (tr->mxtips - 2);
3152	p = tr->nodep[1];
3153	while (nodes-- > 0)
3154	{
3155	for(i = 0; i < tr->numBranches; i++)
3156	p->z[i] = defaultz;
3157
3158	q = p->next;
3159	while(q != p)
3160	{
3161	for(i = 0; i < tr->numBranches; i++)
3162	q->z[i] = defaultz;
3163	q = q->next;
3164	}
3165	p++;
3166	}
3167	}
3168
3169	static double fixZ(double z)
3170	{
3171	if(z > zmax)
3172	return zmax;
3173
3174	if(z < zmin)
3175	return zmin;
3176
3177	return z;
3178	}
3179
3180	static double getFracChange(tree *tr, int model)
3181	{
3182	if(tr->NumberOfModels == 1)
3183	return tr->fracchange;
3184	else
3185	return tr->fracchanges[model];
3186	}
3187
3188	void scaleBranches(tree *tr, boolean fromFile)
3189	{
3190	nodeptr
3191	p;
3192
3193	int
3194	model,
3195	i,
3196	nodes,
3197	count = 0;
3198
3199	double
3200	z;
3201
3202	if(!tr->storedBrLens)
3203	tr->storedBrLens = (double )rax_malloc(sizeof(double) (2 * tr->mxtips - 3) * 2);
3204
3205	assert(tr->numBranches == tr->NumberOfModels);
3206
3207	nodes = tr->mxtips + tr->mxtips - 2;
3208	p = tr->nodep[1];
3209
3210	for(i = 1; i <= nodes; i++)
3211	{
3212	p = tr->nodep[i];
3213
3214	if(fromFile)
3215	{
3216	tr->storedBrLens[count] = p->z[0];
3217
3218	for(model = 0; model < tr->NumberOfModels; model++)
3219	{
3220	z = exp(-p->z[model] / getFracChange(tr, model));
3221
3222	z = fixZ(z);
3223
3224	p->z[model] = z;
3225	}
3226	}
3227	else
3228	{
3229	for(model = 0; model < tr->NumberOfModels; model++)
3230	{
3231	z = tr->partitionData[model].brLenScaler * tr->storedBrLens[count];
3232
3233	z = exp(-z / getFracChange(tr, model));
3234
3235	z = fixZ(z);
3236
3237	p->z[model] = z;
3238	}
3239	}
3240	count++;
3241
3242
3243	if(i > tr->mxtips)
3244	{
3245	if(fromFile)
3246	{
3247	tr->storedBrLens[count] = p->next->z[0];
3248
3249	for(model = 0; model < tr->NumberOfModels; model++)
3250	{
3251	z = exp(-p->next->z[model] / getFracChange(tr, model));
3252
3253	z = fixZ(z);
3254
3255	p->next->z[model] = z;
3256
3257	}
3258	}
3259	else
3260	{
3261	for(model = 0; model < tr->NumberOfModels; model++)
3262	{
3263	z = tr->partitionData[model].brLenScaler * tr->storedBrLens[count];
3264	z = exp(-z / getFracChange(tr, model));
3265
3266	z = fixZ(z);
3267
3268	p->next->z[model] = z;
3269	}
3270	}
3271	count++;
3272
3273	if(fromFile)
3274	{
3275	tr->storedBrLens[count] = p->next->next->z[0];
3276
3277	for(model = 0; model < tr->NumberOfModels; model++)
3278	{
3279	z = exp(-p->next->next->z[model] / getFracChange(tr, model));
3280
3281	z = fixZ(z);
3282
3283	p->next->next->z[model] = z;
3284	}
3285	}
3286	else
3287	{
3288	for(model = 0; model < tr->NumberOfModels; model++)
3289	{
3290	z = tr->partitionData[model].brLenScaler * tr->storedBrLens[count];
3291
3292	z = exp(-z / getFracChange(tr, model));
3293
3294	z = fixZ(z);
3295
3296	p->next->next->z[model] = z;
3297	}
3298	}
3299	count++;
3300	}
3301	}
3302
3303	assert(count == (2 * tr->mxtips - 3) * 2);
3304	}
3305
3306
3307	static void printAAmatrix(tree *tr, double epsilon)
3308	{
3309
3310
3311	if(AAisGTR(tr) \|\| AAisUnlinkedGTR(tr))
3312	{
3313	int
3314	model;
3315
3316	for(model = 0; model < tr->NumberOfModels; model++)
3317	{
3318	if(tr->partitionData[model].dataType == AA_DATA)
3319	{
3320	char
3321	gtrFileName[1024];
3322	/epsilonStr[1024];/
3323
3324	FILE
3325	*gtrFile;
3326
3327	double
3328	*rates = tr->partitionData[model].substRates,
3329	*f = tr->partitionData[model].frequencies,
3330	q[20][20];
3331
3332	int
3333	r = 0,
3334	i,
3335	j;
3336
3337	assert(tr->partitionData[model].protModels == GTR \|\| tr->partitionData[model].protModels == GTR_UNLINKED);
3338
3339	/sprintf(epsilonStr, "%f", epsilon);/
3340
3341	strcpy(gtrFileName, workdir);
3342	strcat(gtrFileName, "RAxML_proteinGTRmodel.");
3343	strcat(gtrFileName, run_id);
3344
3345	/*
3346	strcat(gtrFileName, "_");
3347	strcat(gtrFileName, epsilonStr);
3348	*/
3349
3350	strcat(gtrFileName, "_Partition_");
3351	strcat(gtrFileName, tr->partitionData[model].partitionName);
3352
3353	gtrFile = myfopen(gtrFileName, "wb");
3354
3355	for(i = 0; i < 20; i++)
3356	for(j = 0; j < 20; j++)
3357	q[i][j] = 0.0;
3358
3359	for(i = 0; i < 19; i++)
3360	for(j = i + 1; j < 20; j++)
3361	q[i][j] = rates[r++];
3362
3363	for(i = 0; i < 20; i++)
3364	for(j = 0; j <= i; j++)
3365	{
3366	if(i == j)
3367	q[i][j] = 0.0;
3368	else
3369	q[i][j] = q[j][i];
3370	}
3371
3372	for(i = 0; i < 20; i++)
3373	{
3374	for(j = 0; j < 20; j++)
3375	fprintf(gtrFile, "%1.80f ", q[i][j]);
3376
3377	fprintf(gtrFile, "\n");
3378	}
3379	for(i = 0; i < 20; i++)
3380	fprintf(gtrFile, "%1.80f ", f[i]);
3381	fprintf(gtrFile, "\n");
3382
3383	fclose(gtrFile);
3384
3385	if(tr->partitionData[model].protModels == GTR)
3386	printBothOpen("\nPrinted linked AA GTR matrix that achieved an overall improvement of %f log likelihood units for partition %s to file %s\n\n", epsilon, tr->partitionData[model].partitionName, gtrFileName);
3387	else
3388	printBothOpen("\nPrinted unlinked AA GTR matrix that achieved an overall improvement of %f log likelihood units for partition %s to file %s\n\n", epsilon, tr->partitionData[model].partitionName, gtrFileName);
3389	}
3390
3391	}
3392	}
3393	}
3394
3395
3396
3397	static void optScaler(tree tr, double modelEpsilon, linkageList ll)
3398	{
3399	int
3400	i,
3401	k,
3402	numberOfModels = ll->entries;
3403
3404	double
3405	lim_inf = 0.01,
3406	lim_sup = 100.0,
3407	endLH = (double )rax_malloc(sizeof(double) * numberOfModels),
3408	startLH = (double )rax_malloc(sizeof(double) * numberOfModels),
3409	startAlpha = (double )rax_malloc(sizeof(double) * numberOfModels),
3410	endAlpha = (double )rax_malloc(sizeof(double) * numberOfModels),
3411	_a = (double )rax_malloc(sizeof(double) * numberOfModels),
3412	_b = (double )rax_malloc(sizeof(double) * numberOfModels),
3413	_c = (double )rax_malloc(sizeof(double) * numberOfModels),
3414	_fa = (double )rax_malloc(sizeof(double) * numberOfModels),
3415	_fb = (double )rax_malloc(sizeof(double) * numberOfModels),
3416	_fc = (double )rax_malloc(sizeof(double) * numberOfModels),
3417	_param = (double )rax_malloc(sizeof(double) * numberOfModels),
3418	result = (double )rax_malloc(sizeof(double) * numberOfModels),
3419	_x = (double )rax_malloc(sizeof(double) * numberOfModels);
3420
3421
3422	assert(numberOfModels == tr->numBranches);
3423
3424
3425
3426	evaluateGenericInitrav(tr, tr->start);
3427
3428	#ifdef _DEBUG_MODEL_OPTIMIZATION
3429	double
3430	initialLH = tr->likelihood;
3431	#endif
3432
3433	for(i = 0; i < numberOfModels; i++)
3434	{
3435	assert(ll->ld[i].valid);
3436
3437	startAlpha[i] = tr->partitionData[ll->ld[i].partitionList[0]].brLenScaler;
3438	_a[i] = startAlpha[i] + 0.1;
3439	_b[i] = startAlpha[i] - 0.1;
3440	if(_b[i] < lim_inf)
3441	_b[i] = lim_inf;
3442
3443	startLH[i] = 0.0;
3444	endLH[i] = unlikely;
3445
3446	assert(ll->ld[i].partitions == 1);
3447
3448	for(k = 0; k < ll->ld[i].partitions; k++)
3449	startLH[i] += tr->perPartitionLH[ll->ld[i].partitionList[k]];
3450	}
3451
3452	brakGeneric(_param, _a, _b, _c, _fa, _fb, _fc, lim_inf, lim_sup, numberOfModels, -1, SCALER_F, tr, ll, modelEpsilon);
3453	brentGeneric(_a, _b, _c, _fb, modelEpsilon, _x, result, numberOfModels, SCALER_F, -1, tr, ll, lim_inf, lim_sup);
3454
3455	for(i = 0; i < numberOfModels; i++)
3456	endLH[i] = result[i];
3457
3458
3459	for(i = 0; i < numberOfModels; i++)
3460	{
3461	if(startLH[i] > endLH[i])
3462	{
3463	for(k = 0; k < ll->ld[i].partitions; k++)
3464	{
3465	tr->partitionData[ll->ld[i].partitionList[k]].brLenScaler = startAlpha[i];
3466	scaleBranches(tr, FALSE);
3467	}
3468	}
3469	else
3470	{
3471	for(k = 0; k < ll->ld[i].partitions; k++)
3472	{
3473	tr->partitionData[ll->ld[i].partitionList[k]].brLenScaler = _x[i];
3474	scaleBranches(tr, FALSE);
3475	}
3476	}
3477	}
3478
3479
3480	evaluateGenericInitrav(tr, tr->start);
3481
3482	#ifdef _DEBUG_MODEL_OPTIMIZATION
3483	if(tr->likelihood < initialLH)
3484	printf("%f %f\n", tr->likelihood, initialLH);
3485	assert(tr->likelihood >= initialLH);
3486	#endif
3487
3488	rax_free(startLH);
3489	rax_free(endLH);
3490	rax_free(startAlpha);
3491	rax_free(endAlpha);
3492	rax_free(result);
3493	rax_free(_a);
3494	rax_free(_b);
3495	rax_free(_c);
3496	rax_free(_fa);
3497	rax_free(_fb);
3498	rax_free(_fc);
3499	rax_free(_param);
3500	rax_free(_x);
3501
3502	}
3503
3504	static void autoProtein(tree *tr)
3505	{
3506	int
3507	countAutos = 0,
3508	model;
3509
3510	for(model = 0; model < tr->NumberOfModels; model++)
3511	if(tr->partitionData[model].protModels == AUTO)
3512	countAutos++;
3513
3514	if(countAutos > 0)
3515	{
3516	int
3517	i,
3518	numProteinModels = AUTO,
3519	bestIndex = (int)rax_malloc(sizeof(int) * tr->NumberOfModels),
3520	oldIndex = (int)rax_malloc(sizeof(int) * tr->NumberOfModels);
3521
3522	double
3523	startLH,
3524	bestScores = (double)rax_malloc(sizeof(double) * tr->NumberOfModels);
3525
3526	topolRELL_LIST
3527	rl = (topolRELL_LIST )rax_malloc(sizeof(topolRELL_LIST));
3528
3529	initTL(rl, tr, 1);
3530	saveTL(rl, tr, 0);
3531
3532	evaluateGenericInitrav(tr, tr->start);
3533
3534	startLH = tr->likelihood;
3535
3536	for(model = 0; model < tr->NumberOfModels; model++)
3537	{
3538	oldIndex[model] = tr->partitionData[model].autoProtModels;
3539	bestIndex[model] = -1;
3540	bestScores[model] = unlikely;
3541	}
3542
3543	for(i = 0; i < numProteinModels; i++)
3544	{
3545	for(model = 0; model < tr->NumberOfModels; model++)
3546	{
3547	if(tr->partitionData[model].protModels == AUTO)
3548	{
3549	tr->partitionData[model].autoProtModels = i;
3550	initReversibleGTR(tr, model);
3551	}
3552	}
3553
3554	#ifdef _USE_PTHREADS
3555	masterBarrier(THREAD_COPY_RATES, tr);
3556	#endif
3557
3558	resetBranches(tr);
3559	evaluateGenericInitrav(tr, tr->start);
3560	treeEvaluate(tr, 0.5);
3561
3562	for(model = 0; model < tr->NumberOfModels; model++)
3563	{
3564	if(tr->partitionData[model].protModels == AUTO)
3565	{
3566	if(tr->perPartitionLH[model] > bestScores[model])
3567	{
3568	bestScores[model] = tr->perPartitionLH[model];
3569	bestIndex[model] = i;
3570	}
3571	}
3572	}
3573	}
3574
3575	printBothOpen("Automatic protein model assignment algorithm:\n\n");
3576
3577	for(model = 0; model < tr->NumberOfModels; model++)
3578	{
3579	if(tr->partitionData[model].protModels == AUTO)
3580	{
3581	tr->partitionData[model].autoProtModels = bestIndex[model];
3582	initReversibleGTR(tr, model);
3583	printBothOpen("\tPartition: %d best-scoring AA model: %s likelihood %f\n", model, protModels[tr->partitionData[model].autoProtModels], bestScores[model]);
3584	}
3585	}
3586
3587	printBothOpen("\n\n");
3588
3589	#ifdef _USE_PTHREADS
3590	masterBarrier(THREAD_COPY_RATES, tr);
3591	#endif
3592
3593	resetBranches(tr);
3594	evaluateGenericInitrav(tr, tr->start);
3595	treeEvaluate(tr, 2.0);
3596
3597	if(tr->likelihood < startLH)
3598	{
3599	for(model = 0; model < tr->NumberOfModels; model++)
3600	{
3601	if(tr->partitionData[model].protModels == AUTO)
3602	{
3603	tr->partitionData[model].autoProtModels = oldIndex[model];
3604	initReversibleGTR(tr, model);
3605	}
3606	}
3607
3608	#ifdef _USE_PTHREADS
3609	masterBarrier(THREAD_COPY_RATES, tr);
3610	#endif
3611	restoreTL(rl, tr, 0);
3612	evaluateGenericInitrav(tr, tr->start);
3613	}
3614
3615	assert(tr->likelihood >= startLH);
3616
3617	freeTL(rl);
3618	rax_free(rl);
3619
3620	rax_free(oldIndex);
3621	rax_free(bestIndex);
3622	rax_free(bestScores);
3623	}
3624	}
3625
3626
3627	//#define _DEBUG_MOD_OPT
3628
3629	void modOpt(tree tr, analdef adef, boolean resetModel, double likelihoodEpsilon)
3630	{
3631	int i, model, catOpt = 0;
3632	double
3633	inputLikelihood,
3634	currentLikelihood,
3635	modelEpsilon = 0.0001;
3636	linkageList
3637	*alphaList,
3638	*invarList,
3639	*rateList,
3640	*scalerList;
3641	/*
3642	int linkedAlpha[4] = {0, 0, 0, 0};
3643	int linkedInvar[4] = {0, 0, 0, 0};
3644	int linkedRates[4] = {0, 0, 0, 0};
3645	*/
3646	int
3647	unlinked = (int )rax_malloc(sizeof(int) * tr->NumberOfModels),
3648	linked = (int )rax_malloc(sizeof(int) * tr->NumberOfModels);
3649
3650	assert(!adef->useBinaryModelFile);
3651
3652	modelEpsilon = 0.0001;
3653
3654
3655	for(i = 0; i < tr->NumberOfModels; i++)
3656	{
3657	unlinked[i] = i;
3658	linked[i] = 0;
3659	}
3660
3661	alphaList = initLinkageList(unlinked, tr);
3662	invarList = initLinkageList(unlinked, tr);
3663	rateList = initLinkageListGTR(tr);
3664	scalerList = initLinkageList(unlinked, tr);
3665
3666	if(!(adef->mode == CLASSIFY_ML))
3667	tr->start = tr->nodep[1];
3668
3669	if(resetModel)
3670	{
3671
3672
3673	initRateMatrix(tr);
3674
3675	for(model = 0; model < tr->NumberOfModels; model++)
3676	{
3677	if(adef->useInvariant)
3678	{
3679	int lower, upper;
3680	int count = 0;
3681	int total = 0;
3682
3683	lower = tr->partitionData[model].lower;
3684	upper = tr->partitionData[model].upper;
3685
3686	for(i = lower; i < upper; i++)
3687	{
3688	if(tr->invariant[i] < 4)
3689	count += tr->cdta->aliaswgt[i];
3690	total += tr->cdta->aliaswgt[i];
3691	}
3692
3693	tr->partitionData[model].propInvariant = ((double)count)/((double) total);
3694	}
3695
3696	tr->partitionData[model].alpha = 1.0;
3697
3698	initReversibleGTR(tr, model);
3699
3700	makeGammaCats(tr->partitionData[model].alpha, tr->partitionData[model].gammaRates, 4, tr->useGammaMedian);
3701	}
3702	#ifdef _USE_PTHREADS
3703	masterBarrier(THREAD_RESET_MODEL ,tr);
3704	#endif
3705
3706	resetBranches(tr);
3707
3708	evaluateGenericInitrav(tr, tr->start);
3709
3710
3711
3712	treeEvaluate(tr, 0.25);
3713	}
3714
3715	inputLikelihood = tr->likelihood;
3716
3717	evaluateGenericInitrav(tr, tr->start);
3718
3719	assert(inputLikelihood == tr->likelihood);
3720
3721	/* no need for individual models here, just an init on params equal for all partitions*/
3722
3723	do
3724	{
3725	currentLikelihood = tr->likelihood;
3726
3727	#ifdef _DEBUG_MOD_OPT
3728	printf("start: %1.40f\n", currentLikelihood);
3729	#endif
3730
3731	optRatesGeneric(tr, modelEpsilon, rateList);
3732
3733	evaluateGenericInitrav(tr, tr->start);
3734
3735	#ifdef _DEBUG_MOD_OPT
3736	printf("after rates %1.40f\n", tr->likelihood);
3737	#endif
3738
3739	autoProtein(tr);
3740
3741	if(adef->mode != OPTIMIZE_BR_LEN_SCALER)
3742	treeEvaluate(tr, 0.0625);
3743	else
3744	optScaler(tr, modelEpsilon, scalerList);
3745
3746	#ifdef _DEBUG_MOD_OPT
3747	evaluateGenericInitrav(tr, tr->start);
3748	printf("after br-len 1 %1.40f\n", tr->likelihood);
3749	#endif
3750
3751	switch(tr->rateHetModel)
3752	{
3753	case GAMMA_I:
3754	optAlphasGeneric(tr, modelEpsilon, alphaList);
3755
3756	#ifdef _DEBUG_MOD_OPT
3757	evaluateGenericInitrav(tr, tr->start);
3758	printf("after alphas %1.40f\n", tr->likelihood);
3759	#endif
3760
3761	optInvar(tr, modelEpsilon, invarList);
3762
3763	#ifdef _DEBUG_MOD_OPT
3764	evaluateGenericInitrav(tr, tr->start);
3765	printf("after invar %1.40f\n", tr->likelihood);
3766	#endif
3767
3768	if(adef->mode != OPTIMIZE_BR_LEN_SCALER)
3769	treeEvaluate(tr, 0.1);
3770	else
3771	optScaler(tr, modelEpsilon, scalerList);
3772
3773	#ifdef _DEBUG_MOD_OPT
3774	evaluateGenericInitrav(tr, tr->start);
3775	printf("after br-len 2 %1.40f\n", tr->likelihood);
3776	#endif
3777
3778	break;
3779	case GAMMA:
3780	optAlphasGeneric(tr, modelEpsilon, alphaList);
3781
3782	#ifdef _DEBUG_MOD_OPT
3783	evaluateGenericInitrav(tr, tr->start);
3784	printf("after alphas %1.40f\n", tr->likelihood);
3785	#endif
3786
3787
3788	onlyInitrav(tr, tr->start);
3789
3790	evaluateGeneric(tr, tr->start);
3791
3792	if(adef->mode != OPTIMIZE_BR_LEN_SCALER)
3793	treeEvaluate(tr, 0.1);
3794	else
3795	optScaler(tr, modelEpsilon, scalerList);
3796
3797	#ifdef _DEBUG_MOD_OPT
3798	evaluateGenericInitrav(tr, tr->start);
3799	printf("after br-len 3 %1.40f\n", tr->likelihood);
3800	#endif
3801
3802
3803	break;
3804	case CAT:
3805	if(!tr->noRateHet)
3806	{
3807	if(catOpt < 3)
3808	{
3809	evaluateGenericInitrav(tr, tr->start);
3810	optimizeRateCategories(tr, adef->categories);
3811	catOpt++;
3812	}
3813	}
3814
3815	#ifdef _DEBUG_MOD_OPT
3816	evaluateGenericInitrav(tr, tr->start);
3817	printf("after cat-opt %f\n", tr->likelihood);
3818	#endif
3819
3820	break;
3821	default:
3822	assert(0);
3823	}
3824
3825	if(tr->likelihood < currentLikelihood)
3826	{
3827	if(fabs(tr->likelihood - currentLikelihood) > MIN(0.0000001, likelihoodEpsilon))
3828	{
3829	printf("%1.40f %1.40f\n", tr->likelihood, currentLikelihood);
3830	assert(0);
3831	}
3832	}
3833
3834	printAAmatrix(tr, fabs(currentLikelihood - tr->likelihood));
3835	}
3836	while(fabs(currentLikelihood - tr->likelihood) > likelihoodEpsilon);
3837
3838	rax_free(unlinked);
3839	rax_free(linked);
3840	freeLinkageList(alphaList);
3841	freeLinkageList(rateList);
3842	freeLinkageList(invarList);
3843	freeLinkageList(scalerList);
3844	}
3845
3846
3847
3848
3849	/*******************FUNCTIONS FOOR EXACT MODEL OPTIMIZATION UNDER GTRGAMMA *************************************/
3850
3851
3852
3853	static double branchLength(int model, double z, tree tr)
3854	{
3855	double x;
3856
3857	x = z[model];
3858	assert(x > 0);
3859	if (x < zmin)
3860	x = zmin;
3861
3862
3863	assert(x <= zmax);
3864
3865	if(!tr->multiBranch)
3866	x = -log(x) * tr->fracchange;
3867	else
3868	x = -log(x) * tr->fracchanges[model];
3869
3870	return x;
3871
3872	}
3873
3874
3875	static double treeLengthRec(nodeptr p, tree *tr, int model)
3876	{
3877	double
3878	x = branchLength(model, p->z, tr);
3879
3880	if(isTip(p->number, tr->rdta->numsp))
3881	return x;
3882	else
3883	{
3884	double acc = 0;
3885	nodeptr q;
3886
3887	q = p->next;
3888
3889	while(q != p)
3890	{
3891	acc += treeLengthRec(q->back, tr, model);
3892	q = q->next;
3893	}
3894
3895	return acc + x;
3896	}
3897	}
3898
3899	double treeLength(tree *tr, int model)
3900	{
3901	/* printf("SCALER: %f\n", tr->partitionData[model].brLenScaler); */
3902
3903	return treeLengthRec(tr->start->back, tr, model);
3904	}
3905
3906
3907
3908
3909

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

Context Navigation

source: branches/nameserver/GDE/RAxML/optimizeModel.c

Download in other formats: