Context Navigation

mbbeagle.c

Visit:

Last change on this file was 10416, checked in by aboeckma, 11 years ago
Added mr bayes (no menu yet)
File size: 41.1 KB

Line
1	/*
2	* MrBayes 3
3	*
4	* (c) 2002-2010
5	*
6	* This file originally contributed by:
7	*
8	* Marc A. Suchard
9	* Department of Biomathematics
10	* University of California, Los Angeles
11	* Los Angeles, CA 90095
12	* msuchard@ucla.edu
13	*
14	* With important contributions by
15	*
16	* Maxim Teslenko (maxim.teslenko@nrm.se)
17	*
18	* and by many users (run 'acknowledgements' to see more info)
19	*
20	* This program is free software; you can redistribute it and/or
21	* modify it under the terms of the GNU General Public License
22	* as published by the Free Software Foundation; either version 2
23	* of the License, or (at your option) any later version.
24	*
25	* This program is distributed in the hope that it will be useful,
26	* but WITHOUT ANY WARRANTY; without even the implied warranty of
27	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
28	* GNU General Public License for more details (www.gnu.org).
29	*
30	*/
31
32	#include <stdio.h>
33	#include <stdlib.h>
34	#include <time.h>
35	#include <math.h>
36	#include <string.h>
37	#include <ctype.h>
38	#include <assert.h>
39	#include "mb.h"
40	#include "mbbeagle.h"
41	#include "globals.h"
42	#include "utils.h"
43	#include "mcmc.h"
44	#include "model.h"
45
46	const char* const svnRevisionMbbeagleC="$Rev: 469 $"; /* Revision keyword which is expended/updated by svn on each commit/update*/
47
48	/* Functions and variables defined in mcmc.c that are not exported in mcmc.h */
49	void LaunchLogLikeForDivision(int chain, int d, MrBFlt* lnL);
50
51	void FlipCondLikeSpace (ModelInfo *m, int chain, int nodeIndex);
52	void FlipNodeScalerSpace (ModelInfo *m, int chain, int nodeIndex);
53	void FlipSiteScalerSpace (ModelInfo *m, int chain);
54	void ResetSiteScalers (ModelInfo *m, int chain);
55	void CopySiteScalers (ModelInfo *m, int chain);
56
57	int TreeCondLikes_Beagle (Tree *t, int division, int chain);
58	int TreeLikelihood_Beagle (Tree t, int division, int chain, MrBFlt lnL, int whichSitePats);
59	int TreeTiProbs_Beagle (Tree *t, int division, int chain);
60
61	int TreeCondLikes_Beagle_No_Rescale (Tree *t, int division, int chain);
62	int TreeCondLikes_Beagle_Rescale_All (Tree *t, int division, int chain);
63
64	MrBFlt GetRate (int division, int chain);
65	void FlipTiProbsSpace (ModelInfo *m, int chain, int nodeIndex);
66
67	extern int *chainId;
68	extern int numLocalChains;
69
70	//#define DEBUG_MB_BEAGLE_FLOW
71
72	#if defined (BEAGLE_ENABLED)
73	/*------------------------------------------------------------------------
74	\|
75	\| InitBeagleInstance: create and initialize a beagle instance
76	\|
77	-------------------------------------------------------------------------*/
78	int InitBeagleInstance (ModelInfo *m, int division)
79	{
80	int i, j, k, c, s, *inStates, numPartAmbigTips;
81	double *inPartials;
82	SafeLong *charBits;
83	BeagleInstanceDetails details;
84	long preferedFlags, requiredFlags;
85	int resource;
86
87	if (m->useBeagle == NO)
88	return ERROR;
89
90	/* at least one eigen buffer needed */
91	if (m->nCijkParts == 0)
92	m->nCijkParts = 1;
93
94	/* allocate memory used by beagle */
95	m->logLikelihoods = (MrBFlt ) SafeCalloc ((numLocalChains)m->numChars, sizeof(MrBFlt));
96	m->inRates = (MrBFlt *) SafeCalloc (m->numGammaCats, sizeof(MrBFlt));
97	m->branchLengths = (MrBFlt ) SafeCalloc (2numLocalTaxa, sizeof(MrBFlt));
98	m->tiProbIndices = (int ) SafeCalloc (2numLocalTaxa, sizeof(int));
99	m->inWeights = (MrBFlt ) SafeCalloc (m->numGammaCatsm->nCijkParts, sizeof(MrBFlt));
100	m->bufferIndices = (int *) SafeCalloc (m->nCijkParts, sizeof(int));
101	m->eigenIndices = (int *) SafeCalloc (m->nCijkParts, sizeof(int));
102	m->childBufferIndices = (int *) SafeCalloc (m->nCijkParts, sizeof(int));
103	m->childTiProbIndices = (int *) SafeCalloc (m->nCijkParts, sizeof(int));
104	m->cumulativeScaleIndices = (int *) SafeCalloc (m->nCijkParts, sizeof(int));
105
106	numPartAmbigTips = 0;
107	if (m->numStates != m->numModelStates)
108	numPartAmbigTips = numLocalTaxa;
109	else
110	{
111	for (i=0; i<numLocalTaxa; i++)
112	{
113	if (m->isPartAmbig[i] == YES)
114	numPartAmbigTips++;
115	}
116	}
117
118	if (beagleResourceCount == 0)
119	{
120	preferedFlags = beagleFlags;
121	}
122	else
123	{
124	resource = beagleResource[beagleInstanceCount % beagleResourceCount];
125	preferedFlags = beagleFlags;
126	}
127
128	requiredFlags = 0L;
129
130	if (beagleScalingScheme == MB_BEAGLE_SCALE_ALWAYS)
131	requiredFlags \|= BEAGLE_FLAG_SCALERS_LOG; //BEAGLE_FLAG_SCALERS_RAW;
132
133
134	/* TODO: allocate fewer buffers when nCijkParts > 1 */
135	/* create beagle instance */
136	m->beagleInstance = beagleCreateInstance(numLocalTaxa,
137	m->numCondLikes * m->nCijkParts,
138	numLocalTaxa - numPartAmbigTips,
139	m->numModelStates,
140	m->numChars,
141	(numLocalChains + 1) * m->nCijkParts,
142	m->numTiProbs*m->nCijkParts,
143	m->numGammaCats,
144	m->numScalers * m->nCijkParts,
145	(beagleResourceCount == 0 ? NULL : &resource),
146	(beagleResourceCount == 0 ? 0 : 1),
147	preferedFlags,
148	requiredFlags,
149	&details);
150
151	if (m->beagleInstance < 0)
152	{
153	MrBayesPrint ("%s Failed to start BEAGLE instance\n", spacer);
154	return (ERROR);
155	}
156	else
157	{
158	MrBayesPrint( "\n%s Using BEAGLE resource %i for division %d:", spacer, details.resourceNumber, division+1);
159	#if defined (THREADS_ENABLED)
160	MrBayesPrint( " (%s)\n", (tryToUseThreads ? "threaded" : "non-threaded"));
161	#else
162	MrBayesPrint( " (non-threaded)\n");
163	#endif
164	MrBayesPrint( "%s Rsrc Name : %s\n", spacer, details.resourceName);
165	MrBayesPrint( "%s Impl Name : %s\n", spacer, details.implName);
166	MrBayesPrint( "%s Flags:", spacer);
167	BeaglePrintFlags(details.flags);
168	MrBayesPrint( "\n");
169	beagleInstanceCount++;
170	}
171
172	/* initialize tip data */
173	inStates = (int ) SafeMalloc (m->numChars sizeof(int));
174	if (!inStates)
175	return ERROR;
176	inPartials = (double ) SafeMalloc (m->numChars m->numModelStates * sizeof(double));
177	if (!inPartials)
178	return ERROR;
179	for (i=0; i<numLocalTaxa; i++)
180	{
181	if (m->isPartAmbig[i] == NO)
182	{
183	charBits = m->parsSets[i];
184	for (c=0; c<m->numChars; c++)
185	{
186	for (s=j=0; s<m->numModelStates; s++)
187	{
188	if (IsBitSet(s, charBits))
189	{
190	inStates[c] = s;
191	j++;
192	}
193	}
194	if (j == m->numModelStates)
195	inStates[c] = j;
196	else
197	assert(j==1);
198	charBits += m->nParsIntsPerSite;
199	}
200	beagleSetTipStates(m->beagleInstance, i, inStates);
201	}
202	else /* if (m->isPartAmbig == YES) */
203	{
204	k = 0;
205	charBits = m->parsSets[i];
206	for (c=0; c<m->numChars; c++)
207	{
208	for (s=0; s<m->numModelStates; s++)
209	{
210	if (IsBitSet(s%m->numStates, charBits))
211	inPartials[k++] = 1.0;
212	else
213	inPartials[k++] = 0.0;
214	}
215	charBits += m->nParsIntsPerSite;
216	}
217	beagleSetTipPartials(m->beagleInstance, i, inPartials);
218	}
219	}
220	free (inStates);
221	free (inPartials);
222
223	return NO_ERROR;
224	}
225
226
227
228
229	/*-----------------------------------------------------------------
230	\|
231	\| LaunchBEAGLELogLikeForDivision: calculate the log likelihood
232	\| of the new state of the chain for a single division
233	\|
234	-----------------------------------------------------------------*/
235	void LaunchBEAGLELogLikeForDivision(int chain, int d, ModelInfo* m, Tree* tree, MrBFlt* lnL) {
236
237	int i, rescaleFreqNew;
238	int *isScalerNode;
239	TreeNode *p;
240
241	if (beagleScalingScheme == MB_BEAGLE_SCALE_ALWAYS)
242	{
243
244	#if defined (DEBUG_MB_BEAGLE_FLOW)
245	printf("ALWAYS RESCALING\n");
246	#endif
247	/* Flip and copy or reset site scalers */
248	FlipSiteScalerSpace(m, chain);
249	if (m->upDateAll == YES) {
250	for (i=0; i<m->nCijkParts; i++) {
251	beagleResetScaleFactors(m->beagleInstance, m->siteScalerIndex[chain] + i);
252	}
253	}
254	else
255	CopySiteScalers(m, chain);
256
257	TreeTiProbs_Beagle(tree, d, chain);
258	TreeCondLikes_Beagle(tree, d, chain);
259	TreeLikelihood_Beagle(tree, d, chain, lnL, (chainId[chain] % chainParams.numChains));
260	}
261	else
262	{ /* MB_BEAGLE_SCALE_DYNAMIC */
263
264	/* This flag is only valid within this block */
265	m->rescaleBeagleAll = NO;
266	TreeTiProbs_Beagle(tree, d, chain);
267	if( m->succesCount[chain] > 1000 )
268	{
269	m->succesCount[chain] = 10;
270	m->rescaleFreq[chain]++; /* increase rescaleFreq independent of whether we accept or reject new state*/
271	m->rescaleFreqOld = rescaleFreqNew = m->rescaleFreq[chain];
272	for (i=0; i<tree->nIntNodes; i++)
273	{
274	p = tree->intDownPass[i];
275	if ( p->upDateCl == YES ) {
276	/* flip to the new workspace since TreeCondLikes_Beagle_Rescale_All() does not do it for
277	(p->upDateCl == YES) since it assumes that TreeCondLikes_Beagle_No_Rescale() did it */
278	FlipCondLikeSpace (m, chain, p->index);
279	}
280	}
281	goto rescale_all;
282	}
283
284	if( beagleScalingFrequency != 0 &&
285	m->beagleComputeCount[chain] % (beagleScalingFrequency) == 0 )
286	{
287	m->rescaleFreqOld = rescaleFreqNew = m->rescaleFreq[chain];
288	for (i=0; i<tree->nIntNodes; i++)
289	{
290	p = tree->intDownPass[i];
291	if ( p->upDateCl == YES ) {
292	/* flip to the new workspace since TreeCondLikes_Beagle_Rescale_All() does not do it for (p->upDateCl == YES) since it assumes that TreeCondLikes_Beagle_No_Rescale() did it*/
293	FlipCondLikeSpace (m, chain, p->index);
294	}
295	}
296	goto rescale_all;
297	}
298
299	TreeCondLikes_Beagle_No_Rescale(tree, d, chain);
300
301	/* Check if likelihood is valid */
302	if( TreeLikelihood_Beagle(tree, d, chain, lnL, (chainId[chain] % chainParams.numChains)) == BEAGLE_ERROR_FLOATING_POINT )
303	{
304	m->rescaleFreqOld = rescaleFreqNew = m->rescaleFreq[chain];
305	if(rescaleFreqNew > 1 && m->succesCount[chain] < 40)
306	{
307	if( m->succesCount[chain] < 10 )
308	{
309	if( m->succesCount[chain] < 4 )
310	{
311	rescaleFreqNew-= rescaleFreqNew >> 3; /* <== we cut up to 12,5% of rescaleFreq */
312	if( m->succesCount[chain] < 2 )
313	{
314	rescaleFreqNew-= rescaleFreqNew >> 3;
315	/* to avoid situation when we may stack at high rescaleFreq when new states do not get accepted because of low liklihood but there proposed frequency is high we reduce rescaleFreq even if we reject the last move*/
316	/* basically the higher probability of proposing of low liklihood state which needs smaller rescaleFreq would lead to higher probability of hitting this code which should reduce rescaleFreqOld thus reduce further probability of hitting this code */
317	/* at some point this negative feedback mechanism should get in balance with the mechanism of periodically increasing rescaleFreq when long sequence of successes is achieved*/
318	m->rescaleFreqOld-= m->rescaleFreqOld >> 3;
319	}
320	m->rescaleFreqOld-= m->rescaleFreqOld >> 3;
321	m->rescaleFreqOld--;
322	m->rescaleFreqOld = ( m->rescaleFreqOld ? m->rescaleFreqOld:1);
323	m->recalculateScalers = YES;
324	recalcScalers = YES;
325	}
326	}
327	rescaleFreqNew--;
328	rescaleFreqNew = ( rescaleFreqNew ? rescaleFreqNew:1);
329	}
330	m->succesCount[chain] = 0;
331	rescale_all:
332	#if defined (DEBUG_MB_BEAGLE_FLOW)
333	printf("NUMERICAL RESCALING\n");
334	#endif
335
336	m->rescaleBeagleAll = YES;
337	FlipSiteScalerSpace(m, chain);
338	isScalerNode = m->isScalerNode[chain];
339	while_loop:
340	ResetScalersPartition ( isScalerNode, tree, rescaleFreqNew );
341	for (i=0; i<m->nCijkParts; i++)
342	{
343	beagleResetScaleFactors(m->beagleInstance, m->siteScalerIndex[chain] + i);
344	}
345
346	TreeCondLikes_Beagle_Rescale_All (tree, d, chain);
347	if( TreeLikelihood_Beagle(tree, d, chain, lnL, (chainId[chain] % chainParams.numChains)) == BEAGLE_ERROR_FLOATING_POINT )
348	{
349	if( rescaleFreqNew > 1 )
350	{
351	/Swap back scalers which were swapped in TreeCondLikes_Beagle_Rescale_All() /
352	for (i=0; i<tree->nIntNodes; i++)
353	{
354	p = tree->intDownPass[i];
355	if ( isScalerNode[p->index] == YES)
356	FlipNodeScalerSpace (m, chain, p->index);
357	}
358	rescaleFreqNew -= rescaleFreqNew >> 3; /<== we cut up to 12,5% of rescaleFreq /
359	rescaleFreqNew--; /* we cut extra 1 of rescaleFreq */
360	goto while_loop;
361	}
362	}
363	m->rescaleFreq[chain] = rescaleFreqNew;
364	}
365	}
366
367	/* Count number of evaluations */
368	m->beagleComputeCount[chain]++;
369	}
370
371
372	void recalculateScalers(int chain)
373	{
374	int i, d, rescaleFreqNew;
375	int *isScalerNode;
376	ModelInfo* m;
377	Tree *tree;
378
379	for (d=0; d<numCurrentDivisions; d++)
380	{
381	m = &modelSettings[d];
382	if( m->recalculateScalers == YES)
383	{
384	m->recalculateScalers = NO;
385	tree = GetTree(m->brlens, chain, state[chain]);
386
387	rescaleFreqNew = m->rescaleFreq[chain];
388	isScalerNode = m->isScalerNode[chain];
389
390	ResetScalersPartition ( isScalerNode, tree, rescaleFreqNew );
391	for (i=0; i<m->nCijkParts; i++) {
392	beagleResetScaleFactors(m->beagleInstance, m->siteScalerIndex[chain] + i);
393	}
394	/* here it does not matter if we flip CL space or not */
395	TreeCondLikes_Beagle_Rescale_All (tree, d, chain);
396	}
397	}
398	}
399
400	void BeagleAddGPUDevicesToList(int *newResourceList, int beagleResourceCount) {
401	BeagleResourceList* beagleResources;
402	int i, gpuCount;
403
404	beagleResources = beagleGetResourceList();
405	if (*newResourceList == NULL) {
406	newResourceList = (int) SafeCalloc(sizeof(int), beagleResources->length);
407	}
408	gpuCount = 0;
409	for (i = 0; i < beagleResources->length; i++) {
410	if (beagleResources->list[i].supportFlags & BEAGLE_FLAG_PROCESSOR_GPU) {
411	(*newResourceList)[gpuCount] = i;
412	gpuCount++;
413	}
414	}
415	*beagleResourceCount = gpuCount;
416	}
417
418	void BeagleRemoveGPUDevicesFromList(int *beagleResource, int beagleResourceCount) {
419	*beagleResourceCount = 0;
420	}
421
422	/*-----
423	\|
424	\| BeaglePrintResources: outputs the available BEAGLE resources
425	\|
426	----------*/
427
428	void BeaglePrintResources()
429	{
430	int i;
431	BeagleResourceList* beagleResources;
432
433	beagleResources = beagleGetResourceList();
434	MrBayesPrint ("%s Available resources reported by beagle library:\n", spacer);
435	for (i=0; i<beagleResources->length; i++)
436	{
437	MrBayesPrint ("\tResource %i:\n", i);
438	MrBayesPrint ("\tName: %s\n", beagleResources->list[i].name);
439	if (i > 0)
440	{
441	MrBayesPrint ("\tDesc: %s\n", beagleResources->list[i].description);
442	}
443	MrBayesPrint("\tFlags:");
444	BeaglePrintFlags(beagleResources->list[i].supportFlags);
445	MrBayesPrint("\n\n");
446	}
447	}
448
449	int BeagleCheckFlagCompatability(long inFlags) {
450	if (inFlags & BEAGLE_FLAG_PROCESSOR_GPU) {
451	if (inFlags & BEAGLE_FLAG_VECTOR_SSE) {
452	MrBayesPrint("%s Simultaneous use of GPU and SSE not available.\n", spacer);
453	return NO;
454	}
455	if (inFlags & BEAGLE_FLAG_THREADING_OPENMP) {
456	MrBayesPrint("%s Simultaneous use of GPU and OpenMP not available.\n", spacer);
457	return NO;
458	}
459	}
460
461	return YES;
462	}
463
464
465	/*-------------------
466	\|
467	\| BeaglePrintFlags: outputs beagle instance details
468	\|
469	______________________*/
470	void BeaglePrintFlags(long inFlags)
471	{
472	int i, k;
473	char *names[] = { "PROCESSOR_CPU",
474	"PROCESSOR_GPU",
475	"PROCESSOR_FPGA",
476	"PROCESSOR_CELL",
477	"PRECISION_DOUBLE",
478	"PRECISION_SINGLE",
479	"COMPUTATION_ASYNCH",
480	"COMPUTATION_SYNCH",
481	"EIGEN_REAL",
482	"EIGEN_COMPLEX",
483	"SCALING_MANUAL",
484	"SCALING_AUTO",
485	"SCALING_ALWAYS",
486	"SCALING_DYNAMIC",
487	"SCALERS_RAW",
488	"SCALERS_LOG",
489	"VECTOR_NONE",
490	"VECTOR_SSE",
491	"THREADING_NONE",
492	"THREADING_OPENMP"
493	};
494	long flags[] = { BEAGLE_FLAG_PROCESSOR_CPU,
495	BEAGLE_FLAG_PROCESSOR_GPU,
496	BEAGLE_FLAG_PROCESSOR_FPGA,
497	BEAGLE_FLAG_PROCESSOR_CELL,
498	BEAGLE_FLAG_PRECISION_DOUBLE,
499	BEAGLE_FLAG_PRECISION_SINGLE,
500	BEAGLE_FLAG_COMPUTATION_ASYNCH,
501	BEAGLE_FLAG_COMPUTATION_SYNCH,
502	BEAGLE_FLAG_EIGEN_REAL,
503	BEAGLE_FLAG_EIGEN_COMPLEX,
504	BEAGLE_FLAG_SCALING_MANUAL,
505	BEAGLE_FLAG_SCALING_AUTO,
506	BEAGLE_FLAG_SCALING_ALWAYS,
507	BEAGLE_FLAG_SCALING_DYNAMIC,
508	BEAGLE_FLAG_SCALERS_RAW,
509	BEAGLE_FLAG_SCALERS_LOG,
510	BEAGLE_FLAG_VECTOR_NONE,
511	BEAGLE_FLAG_VECTOR_SSE,
512	BEAGLE_FLAG_THREADING_NONE,
513	BEAGLE_FLAG_THREADING_OPENMP
514	};
515
516	for (i=k=0; i<20; i++)
517	{
518	if (inFlags & flags[i])
519	{
520	if (k%4 == 0 && k > 0)
521	MrBayesPrint("\n%s ", spacer);
522	MrBayesPrint (" %s", names[i]);
523	k++;
524	}
525	}
526	}
527
528	int ScheduleLogLikeForAllDivisions() {
529	int d;
530	int divisionsToLaunch = 0;
531	ModelInfo *m;
532
533	if (numCurrentDivisions < 2) {
534	return 0;
535	}
536
537	for (d=0; d<numCurrentDivisions; d++) {
538	m = &modelSettings[d];
539	if (m->upDateCl == YES) {
540	divisionsToLaunch++;
541	}
542	}
543	return (divisionsToLaunch > 1);
544	}
545
546	#if defined(THREADS_ENABLED)
547	void LaunchThreadLogLikeForDivision(void arguments) {
548	int d, chain;
549	MrBFlt *lnL;
550	LaunchStruct* launchStruct;
551
552	launchStruct = (LaunchStruct*) arguments;
553	chain = launchStruct->chain;
554	d = launchStruct->division;
555	lnL = launchStruct->lnL;
556	LaunchLogLikeForDivision(chain, d, lnL);
557	return 0;
558	}
559
560	MrBFlt LaunchLogLikeForAllDivisionsInParallel(int chain) {
561	int d;
562	int threadError;
563	pthread_t* threads;
564	LaunchStruct* launchValues;
565	int* wait;
566	ModelInfo* m;
567	MrBFlt chainLnLike;
568
569	chainLnLike = 0.0;
570
571	/* TODO Initialize only once */
572	threads = (pthread_t) SafeMalloc(sizeof(pthread_t) numCurrentDivisions);
573	launchValues = (LaunchStruct) SafeMalloc(sizeof(LaunchStruct) numCurrentDivisions);
574	wait = (int) SafeMalloc(sizeof(int) numCurrentDivisions);
575
576	/* Cycle through divisions and recalculate tis and cond likes as necessary. */
577	/* Code below does not try to avoid recalculating ti probs for divisions */
578	/* that could share ti probs with other divisions. */
579	for (d=0; d<numCurrentDivisions; d++)
580	{
581
582	#if defined (BEST_MPI_ENABLED)
583	if (isDivisionActive[d] == NO)
584	continue;
585	#endif
586	m = &modelSettings[d];
587
588	if (m->upDateCl == YES)
589	{
590	launchValues[d].chain = chain;
591	launchValues[d].division = d;
592	launchValues[d].lnL = &(m->lnLike[2*chain + state[chain]]);
593	/* Fork */
594	threadError = pthread_create(&threads[d], NULL,
595	LaunchThreadLogLikeForDivision,
596	(void*) &launchValues[d]);
597	assert(0 == threadError);
598	wait[d] = 1;
599	}
600	else
601	{
602	wait[d] = 0;
603	}
604	}
605
606	for (d = 0; d < numCurrentDivisions; d++)
607	{
608	/* Join */
609	if (wait[d])
610	{
611	threadError = pthread_join(threads[d], NULL);
612	assert(0 == threadError);
613	}
614	m = &modelSettings[d];
615	chainLnLike += m->lnLike[2*chain + state[chain]];
616	}
617
618	/* TODO Free these once */
619	free(threads);
620	free(launchValues);
621	free(wait);
622
623	return chainLnLike;
624	}
625	#endif
626
627	/*----------------------------------------------------------------
628	\|
629	\| TreeCondLikes_Beagle: This routine updates all conditional
630	\| (partial) likelihoods of a beagle instance while doing no rescaling.
631	\| That potentialy can make final liklihood bad then calculation with rescaling needs to be done.
632	\|
633	-----------------------------------------------------------------*/
634	int TreeCondLikes_Beagle_No_Rescale (Tree *t, int division, int chain)
635	{
636	int i, j, cumulativeScaleIndex;
637	BeagleOperation operations;
638	TreeNode *p;
639	ModelInfo *m;
640	unsigned chil1Step, chil2Step;
641	int *isScalerNode;
642
643	m = &modelSettings[division];
644	isScalerNode = m->isScalerNode[chain];
645
646	for (i=0; i<t->nIntNodes; i++)
647	{
648	p = t->intDownPass[i];
649
650	/* check if conditional likelihoods need updating */
651	if (p->upDateCl == YES)
652	{
653	/* flip to the new workspace */
654	FlipCondLikeSpace (m, chain, p->index);
655
656	/* update conditional likelihoods */
657	operations.destinationPartials = m->condLikeIndex[chain][p->index ];
658	operations.child1Partials = m->condLikeIndex[chain][p->left->index ];
659	operations.child1TransitionMatrix = m->tiProbsIndex [chain][p->left->index ];
660	operations.child2Partials = m->condLikeIndex[chain][p->right->index];
661	operations.child2TransitionMatrix = m->tiProbsIndex [chain][p->right->index];
662
663	/* All partials for tips are the same across omega catigoris, thus we are doing the following two if statments.*/
664	if(p->left->left== NULL )
665	chil1Step=0;
666	else
667	chil1Step=1;
668
669	if(p->right->left== NULL )
670	chil2Step=0;
671	else
672	chil2Step=1;
673
674	operations.destinationScaleWrite = BEAGLE_OP_NONE;
675	cumulativeScaleIndex = BEAGLE_OP_NONE;
676	if ( isScalerNode[p->index] == YES )
677	{
678	operations.destinationScaleRead = m->nodeScalerIndex[chain][p->index];
679	}
680	else
681	{
682	operations.destinationScaleRead = BEAGLE_OP_NONE;
683	}
684
685	for (j=0; j<m->nCijkParts; j++)
686	{
687	beagleUpdatePartials(m->beagleInstance,
688	&operations,
689	1,
690	cumulativeScaleIndex);
691
692	operations.destinationPartials++;
693	operations.child1Partials+=chil1Step;
694	operations.child1TransitionMatrix++;
695	operations.child2Partials+=chil2Step;
696	operations.child2TransitionMatrix++;
697
698	if ( isScalerNode[p->index] == YES )
699	operations.destinationScaleRead++;
700	}
701	}
702	}
703
704	return NO_ERROR;
705	}
706
707
708
709
710	/*----------------------------------------------------------------
711	\|
712	\| TreeCondLikes_Beagle: This routine updates all conditional
713	\| (partial) likelihoods of a beagle instance while rescaling at every node.
714	\| Note: all nodes get recalculated, not only tached by move.
715	\|
716	-----------------------------------------------------------------*/
717	int TreeCondLikes_Beagle_Rescale_All (Tree *t, int division, int chain)
718	{
719	int i, j, cumulativeScaleIndex;
720	BeagleOperation operations;
721	TreeNode *p;
722	ModelInfo *m;
723	unsigned chil1Step, chil2Step;
724	int *isScalerNode;
725
726	m = &modelSettings[division];
727	isScalerNode = m->isScalerNode[chain];
728
729	for (i=0; i<t->nIntNodes; i++)
730	{
731	p = t->intDownPass[i];
732
733	if (p->upDateCl == NO ) {
734	//p->upDateCl = YES;
735	/* flip to the new workspace */
736	FlipCondLikeSpace (m, chain, p->index);
737	}
738
739	/* update conditional likelihoods */
740	operations.destinationPartials = m->condLikeIndex[chain][p->index ];
741	operations.child1Partials = m->condLikeIndex[chain][p->left->index ];
742	operations.child1TransitionMatrix = m->tiProbsIndex [chain][p->left->index ];
743	operations.child2Partials = m->condLikeIndex[chain][p->right->index];
744	operations.child2TransitionMatrix = m->tiProbsIndex [chain][p->right->index];
745
746	/* All partials for tips are the same across omega catigoris, thus we are doing the following two if statments.*/
747	if(p->left->left== NULL)
748	chil1Step=0;
749	else
750	chil1Step=1;
751
752	if(p->right->left== NULL)
753	chil2Step=0;
754	else
755	chil2Step=1;
756
757	operations.destinationScaleRead = BEAGLE_OP_NONE;
758	if ( isScalerNode[p->index] == YES )
759	{
760	FlipNodeScalerSpace (m, chain, p->index);
761	operations.destinationScaleWrite = m->nodeScalerIndex[chain][p->index];
762	cumulativeScaleIndex = m->siteScalerIndex[chain];
763	}
764	else
765	{
766	operations.destinationScaleWrite = BEAGLE_OP_NONE;
767	cumulativeScaleIndex = BEAGLE_OP_NONE;
768	}
769
770
771
772	for (j=0; j<m->nCijkParts; j++)
773	{
774	beagleUpdatePartials(m->beagleInstance,
775	&operations,
776	1,
777	cumulativeScaleIndex);
778
779	operations.destinationPartials++;
780	operations.child1Partials+=chil1Step;
781	operations.child1TransitionMatrix++;
782	operations.child2Partials+=chil2Step;
783	operations.child2TransitionMatrix++;
784
785	if ( isScalerNode[p->index] == YES ){
786	operations.destinationScaleWrite++;
787	cumulativeScaleIndex++;
788	}
789
790	}
791	}
792
793	return NO_ERROR;
794	}
795
796
797	/*----------------------------------------------------------------
798	\|
799	\| TreeCondLikes_Beagle: This routine updates all conditional
800	\| (partial) likelihoods of a beagle instance.
801	\|
802	-----------------------------------------------------------------*/
803	int TreeCondLikes_Beagle (Tree *t, int division, int chain)
804	{
805	int i, j, destinationScaleRead, cumulativeScaleIndex;
806	BeagleOperation operations;
807	TreeNode *p;
808	ModelInfo *m;
809	unsigned chil1Step, chil2Step;
810
811	m = &modelSettings[division];
812
813	for (i=0; i<t->nIntNodes; i++)
814	{
815	p = t->intDownPass[i];
816
817	/* check if conditional likelihoods need updating */
818	if (p->upDateCl == YES)
819	{
820	/* remove old scalers */
821	if (m->scalersSet[chain][p->index] == YES && m->upDateAll == NO)
822	{
823	destinationScaleRead = m->nodeScalerIndex[chain][p->index];
824	cumulativeScaleIndex = m->siteScalerIndex[chain];
825	for (j=0; j<m->nCijkParts; j++)
826	{
827	beagleRemoveScaleFactors(m->beagleInstance,
828	&destinationScaleRead,
829	1,
830	cumulativeScaleIndex);
831	destinationScaleRead++;
832	cumulativeScaleIndex++;
833	}
834	}
835
836	/* flip to the new workspace */
837	FlipCondLikeSpace (m, chain, p->index);
838	FlipNodeScalerSpace (m, chain, p->index);
839	m->scalersSet[chain][p->index] = NO;
840
841	/* update conditional likelihoods */
842	operations.destinationPartials = m->condLikeIndex[chain][p->index ];
843	operations.child1Partials = m->condLikeIndex[chain][p->left->index ];
844	operations.child1TransitionMatrix = m->tiProbsIndex [chain][p->left->index ];
845	operations.child2Partials = m->condLikeIndex[chain][p->right->index];
846	operations.child2TransitionMatrix = m->tiProbsIndex [chain][p->right->index];
847
848	/* All partials for tips are the same across omega catigoris, thus we are doing the following two if statments.*/
849	if(p->left->left== NULL && p->left->right== NULL)
850	chil1Step=0;
851	else
852	chil1Step=1;
853
854	if(p->right->left== NULL && p->right->right== NULL)
855	chil2Step=0;
856	else
857	chil2Step=1;
858
859	if (p->scalerNode == YES)
860	{
861	m->scalersSet[chain][p->index] = YES;
862	operations.destinationScaleWrite = m->nodeScalerIndex[chain][p->index];
863	cumulativeScaleIndex = m->siteScalerIndex[chain];
864	}
865	else
866	{
867	operations.destinationScaleWrite = BEAGLE_OP_NONE;
868	cumulativeScaleIndex = BEAGLE_OP_NONE;
869	}
870	operations.destinationScaleRead = BEAGLE_OP_NONE;
871
872	for (j=0; j<m->nCijkParts; j++)
873	{
874	beagleUpdatePartials(m->beagleInstance,
875	&operations,
876	1,
877	cumulativeScaleIndex);
878
879	operations.destinationPartials++;
880	operations.child1Partials+=chil1Step;
881	operations.child1TransitionMatrix++;
882	operations.child2Partials+=chil2Step;
883	operations.child2TransitionMatrix++;
884	if (p->scalerNode == YES)
885	{
886	operations.destinationScaleWrite++;
887	cumulativeScaleIndex++;
888	}
889	}
890
891	}
892	} /* end of for */
893
894	return NO_ERROR;
895	}
896
897
898
899
900
901	/**---------------------------------------------------------------------------
902	\|
903	\| TreeLikelihood_Beagle: Accumulate the log likelihoods calculated by Beagle
904	\| at the root.
905	\|
906	---------------------------------------- -------------------------------------*/
907	int TreeLikelihood_Beagle (Tree t, int division, int chain, MrBFlt lnL, int whichSitePats)
908
909	{
910	int i, j, c = 0, nStates, hasPInvar, beagleReturn;
911	MrBFlt swr, s01, s10, probOn, probOff, covBF[40], pInvar=0.0, bs, freq, likeI, lnLikeI, diff, *omegaCatFreq;
912	CLFlt clInvar=NULL, nSitesOfPat;
913	double *nSitesOfPat_Beagle;
914	TreeNode *p;
915	ModelInfo *m;
916	double pUnobserved;
917
918	#if defined (MB_PRINT_DYNAMIC_RESCALE_FAIL_STAT)
919	static unsigned countBeagleDynamicFail=0;
920	static unsigned countALL=0;
921	#endif
922
923	/* find root node */
924	p = t->root->left;
925
926	/* find model settings and nStates, pInvar, invar cond likes */
927	m = &modelSettings[division];
928
929	nStates = m->numModelStates;
930	if (m->pInvar == NULL)
931	{
932	hasPInvar = NO;
933	}
934	else
935	{
936	hasPInvar = YES;
937	pInvar = *(GetParamVals (m->pInvar, chain, state[chain]));
938	clInvar = m->invCondLikes;
939	}
940
941	/* find base frequencies */
942	bs = GetParamSubVals (m->stateFreq, chain, state[chain]);
943
944	/* if covarion model, adjust base frequencies */
945	if (m->switchRates != NULL)
946	{
947	/* find the stationary frequencies */
948	swr = GetParamVals(m->switchRates, chain, state[chain]);
949	s01 = swr[0];
950	s10 = swr[1];
951	probOn = s01 / (s01 + s10);
952	probOff = 1.0 - probOn;
953
954	/* now adjust the base frequencies; on-state stored first in cond likes */
955	for (j=0; j<nStates/2; j++)
956	{
957	covBF[j] = bs[j] * probOn;
958	covBF[j+nStates/2] = bs[j] * probOff;
959	}
960
961	/* finally set bs pointer to adjusted values */
962	bs = covBF;
963	}
964
965	if (m->upDateCijk == YES) { /* TODO Really only need to check if state frequencies have changed */
966
967	/* set base frequencies in BEAGLE instance */
968	for (i=0; i<m->nCijkParts; i++)
969	beagleSetStateFrequencies(m->beagleInstance,
970	m->cijkIndex[chain] + i,
971	bs);
972	}
973
974	/* find category frequencies */
975	if (hasPInvar == NO)
976	freq = 1.0 / m->numGammaCats;
977	else
978	freq = (1.0 - pInvar) / m->numGammaCats;
979
980	/* TODO: cat weights only need to be set when they change */
981	/* set category frequencies in beagle instance */
982	if (m->numOmegaCats > 1)
983	{
984	omegaCatFreq = GetParamSubVals(m->omega, chain, state[chain]);
985	for (i=0; i<m->nCijkParts; i++)
986	{
987	for (j=0; j<m->numGammaCats; j++)
988	m->inWeights[j] = freq * omegaCatFreq[i];
989	beagleSetCategoryWeights(m->beagleInstance,
990	m->cijkIndex[chain] + i,
991	m->inWeights);
992	}
993	}
994	else if (hasPInvar == YES)
995	{
996	for (i=0; i<m->numGammaCats; i++)
997	m->inWeights[i] = freq;
998	beagleSetCategoryWeights(m->beagleInstance,
999	m->cijkIndex[chain],
1000	m->inWeights);
1001	}
1002
1003	/* find nSitesOfPat */
1004	nSitesOfPat = numSitesOfPat + (whichSitePats*numCompressedChars) + m->compCharStart;
1005
1006	/* TODO: pattern weights only need to be set when they change */
1007	/* set pattern weights in beagle instance if using dynamic reweighting */
1008	if (chainParams.weightScheme[0] + chainParams.weightScheme[1] > ETA)
1009	{
1010	nSitesOfPat_Beagle = (double ) SafeMalloc (m->numChars sizeof(double));
1011	for (c=0; c<m->numChars; c++)
1012	nSitesOfPat_Beagle[c] = numSitesOfPat[m->compCharStart + c];
1013	beagleSetPatternWeights(m->beagleInstance,
1014	nSitesOfPat_Beagle);
1015	SafeFree ((void **)(&nSitesOfPat_Beagle));
1016	}
1017
1018	/* find root log likelihoods and scalers */
1019	for (i=0; i<m->nCijkParts; i++)
1020	{
1021	m->bufferIndices[i] = m->condLikeIndex[chain][p->index] + i;
1022	m->eigenIndices[i] = m->cijkIndex[chain] + i;
1023	m->cumulativeScaleIndices[i] = m->siteScalerIndex[chain] + i;
1024	if (t->isRooted == NO)
1025	{
1026	m->childBufferIndices[i] = m->condLikeIndex [chain][p->anc->index];
1027	m->childTiProbIndices[i] = m->tiProbsIndex [chain][p->index] + i;
1028	}
1029	}
1030
1031	/* reset lnL */
1032	*lnL = 0.0;
1033
1034	/* get root log likelihoods */
1035	if (t->isRooted == YES)
1036	{
1037	beagleReturn = beagleCalculateRootLogLikelihoods(m->beagleInstance,
1038	m->bufferIndices,
1039	m->eigenIndices,
1040	m->eigenIndices,
1041	m->cumulativeScaleIndices,
1042	m->nCijkParts,
1043	lnL);
1044
1045	}
1046	else
1047	{
1048	beagleReturn = beagleCalculateEdgeLogLikelihoods(m->beagleInstance,
1049	m->bufferIndices,
1050	m->childBufferIndices,
1051	m->childTiProbIndices,
1052	NULL,
1053	NULL,
1054	m->eigenIndices,
1055	m->eigenIndices,
1056	m->cumulativeScaleIndices,
1057	m->nCijkParts,
1058	lnL,
1059	NULL,
1060	NULL);
1061	}
1062	#if defined (MB_PRINT_DYNAMIC_RESCALE_FAIL_STAT)
1063	countALL++;
1064	#endif
1065	if( beagleReturn == BEAGLE_ERROR_FLOATING_POINT )
1066	{
1067	#if defined (MB_PRINT_DYNAMIC_RESCALE_FAIL_STAT)
1068	countBeagleDynamicFail++;
1069	printf("#####DEBUG INFO (it is not an error)############## countBeagleDynamicFail:%d countALL:%d\n",countBeagleDynamicFail,countALL);
1070	#endif
1071	return beagleReturn;
1072	}
1073	assert(beagleReturn == BEAGLE_SUCCESS);
1074	m->succesCount[chain]++;
1075
1076	/* accumulate logs across sites */
1077	if (hasPInvar == NO)
1078	{
1079	if( m->dataType == RESTRICTION )
1080	{
1081	beagleGetSiteLogLikelihoods(m->beagleInstance, m->logLikelihoods);
1082	(*lnL) = 0.0;
1083	pUnobserved = 0.0;
1084	for (c=0; c<m->numDummyChars; c++)
1085	{
1086	pUnobserved += exp((double)m->logLikelihoods[c]);
1087	}
1088	/* correct for absent characters */
1089	(lnL) -= log( 1-pUnobserved ) (m->numUncompressedChars);
1090	for (; c<m->numChars; c++)
1091	{
1092	(lnL) += m->logLikelihoods[c] nSitesOfPat[c];
1093	}
1094	}
1095	/* already done, just check for numerical errors */
1096	assert ((lnL) == (lnL));
1097	}
1098	else
1099	{
1100	/* has invariable category */
1101	beagleGetSiteLogLikelihoods(m->beagleInstance,
1102	m->logLikelihoods);
1103	(*lnL) = 0.0;
1104	for (c=0; c<m->numChars; c++)
1105	{
1106	likeI = 0.0;
1107	for (j=0; j<nStates; j++)
1108	likeI += ((clInvar++)) bs[j];
1109	if (likeI != 0.0)
1110	{
1111	lnLikeI = log(likeI * pInvar);
1112	diff = lnLikeI - m->logLikelihoods[c];
1113	}
1114	else
1115	diff = -1000.0;
1116	if (diff < -200.0)
1117	(lnL) += m->logLikelihoods[c] nSitesOfPat[c];
1118	else if (diff > 200.0)
1119	(lnL) += lnLikeI nSitesOfPat[c];
1120	else
1121	{
1122	(lnL) += (m->logLikelihoods[c] + log(1.0 + exp(diff))) nSitesOfPat[c];
1123	}
1124
1125	/* check for numerical errors */
1126	assert((lnL) == (lnL));
1127	}
1128	}
1129
1130	return (NO_ERROR);
1131	}
1132
1133
1134
1135
1136
1137	/*----------------------------------------------------------------
1138	\|
1139	\| TreeTiProbs_Beagle: This routine updates all transition
1140	\| probability matrices of a beagle instance.
1141	\|
1142	-----------------------------------------------------------------*/
1143	int TreeTiProbs_Beagle (Tree *t, int division, int chain)
1144
1145	{
1146
1147	int i, j, k, count;
1148	MrBFlt correctionFactor, theRate, baseRate, *catRate, length;
1149	TreeNode *p;
1150	ModelInfo *m;
1151
1152	/* get model settings */
1153	m = &modelSettings[division];
1154
1155	/* find the correction factor to make branch lengths
1156	in terms of expected number of substitutions per character */
1157	correctionFactor = 1.0;
1158	if (m->dataType == DNA \|\| m->dataType == RNA)
1159	{
1160	if (m->nucModelId == NUCMODEL_DOUBLET)
1161	correctionFactor = 2.0;
1162	else if (m->nucModelId == NUCMODEL_CODON)
1163	correctionFactor = 3.0;
1164	}
1165
1166	/* get rate multipliers (for gamma & partition specific rates) */
1167	theRate = 1.0;
1168	baseRate = GetRate (division, chain);
1169
1170	/* compensate for invariable sites if appropriate */
1171	if (m->pInvar != NULL)
1172	baseRate /= ( 1.0 - ( *GetParamVals(m->pInvar, chain, state[chain])));
1173
1174	/* get category rates for gamma */
1175	if (m->shape == NULL)
1176	catRate = &theRate;
1177	else
1178	catRate = GetParamSubVals (m->shape, chain, state[chain]);
1179
1180	/* get effective category rates */
1181	for (k=0; k<m->numGammaCats; k++)
1182	m->inRates[k] = baseRate * catRate[k] * correctionFactor;
1183
1184	/* TODO: only need to set category rates when they change */
1185	/* set category rates */
1186	beagleSetCategoryRates(m->beagleInstance, m->inRates);
1187
1188	/* get ti prob indices and branch lengths to update */
1189	for (i=count=0; i<t->nNodes; i++)
1190	{
1191	p = t->allDownPass[i];
1192
1193	/* check if transition probs need updating */
1194	if (p->upDateTi == YES)
1195	{
1196	/* flip transition probability */
1197	FlipTiProbsSpace (m, chain, p->index);
1198
1199	/* find length */
1200	if (m->cppEvents != NULL)
1201	{
1202	length = GetParamSubVals (m->cppEvents, chain, state[chain])[p->index];
1203	}
1204	else if (m->tk02BranchRates != NULL)
1205	{
1206	length = GetParamSubVals (m->tk02BranchRates, chain, state[chain])[p->index];
1207	}
1208	else if (m->igrBranchRates != NULL)
1209	{
1210	length = GetParamSubVals (m->igrBranchRates, chain, state[chain])[p->index];
1211	}
1212	else
1213	length = p->length;
1214
1215	/* numerical errors might ensue if we allow very large or very small branch lengths, which might
1216	occur in relaxed clock models; an elegant solution would be to substitute the stationary
1217	probs and initial probs but for now we truncate lengths at small or large values */
1218	if (length > BRLENS_MAX)
1219	length = BRLENS_MAX;
1220	else if (length < BRLENS_MIN)
1221	length = BRLENS_MIN;
1222
1223	m->branchLengths[count] = length;
1224
1225	/* find index */
1226	m->tiProbIndices[count] = m->tiProbsIndex[chain][p->index];
1227	count++;
1228	}
1229	}
1230
1231	/* TODO: only need to update branches that have changed */
1232	/* calculate transition probabilities */
1233	if (count > 0) {
1234	for (i=0; i<m->nCijkParts; i++)
1235	{
1236	beagleUpdateTransitionMatrices(m->beagleInstance,
1237	m->cijkIndex[chain] + i,
1238	m->tiProbIndices,
1239	NULL,
1240	NULL,
1241	m->branchLengths,
1242	count);
1243	for (j=0; j<count; j++)
1244	m->tiProbIndices[j]++;
1245	}
1246	}
1247
1248	/* return success */
1249	return NO_ERROR;
1250	}
1251
1252	#endif // BEAGLE_ENABLED
1253
1254	void BeagleNotLinked()
1255	{
1256	MrBayesPrint("%s BEAGLE library is not linked to this executable.\n", spacer);
1257	}
1258
1259	void BeagleThreadsNotLinked()
1260	{
1261	MrBayesPrint("%s Pthreads library is not linked to this executable.\n", spacer);
1262	}

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

Context Navigation

source: trunk/GDE/MrBAYES/mrbayes_3.2.1/mbbeagle.c

Download in other formats: