6 * This file originally contributed by:
9 * Department of Biomathematics
10 * University of California, Los Angeles
11 * Los Angeles, CA 90095
14 * With important contributions by
16 * Maxim Teslenko (maxim.teslenko@nrm.se)
18 * and by many users (run 'acknowledgments' to see more info)
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.
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).
38 const char* const svnRevisionMbbeagleC = "$Rev: 1008 $"; /* Revision keyword which is expended/updated by svn on each commit/update */
40 /* Functions and variables defined in mcmc.c that are not exported in mcmc.h */
41 void LaunchLogLikeForDivision(int chain, int d, MrBFlt* lnL);
43 void FlipCondLikeSpace (ModelInfo *m, int chain, int nodeIndex);
44 void FlipNodeScalerSpace (ModelInfo *m, int chain, int nodeIndex);
45 void FlipSiteScalerSpace (ModelInfo *m, int chain);
46 void FlipTiProbsSpace (ModelInfo *m, int chain, int nodeIndex);
47 void ResetSiteScalers (ModelInfo *m, int chain);
48 void CopySiteScalers (ModelInfo *m, int chain);
50 int TreeCondLikes_Beagle_No_Rescale (Tree *t, int division, int chain);
51 int TreeCondLikes_Beagle_Rescale_All (Tree *t, int division, int chain);
54 extern int numLocalChains;
57 #if defined (BEAGLE_ENABLED)
58 /*------------------------------------------------------------------------
60 | InitBeagleInstance: create and initialize a beagle instance
62 -------------------------------------------------------------------------*/
63 int InitBeagleInstance (ModelInfo *m, int division)
65 int i, j, k, c, s, *inStates, numPartAmbigTips;
68 BeagleInstanceDetails details;
69 long preferedFlags, requiredFlags;
72 if (m->useBeagle == NO)
75 /* at least one eigen buffer needed */
76 if (m->nCijkParts == 0)
79 /* allocate memory used by beagle */
80 m->logLikelihoods = (MrBFlt *) SafeCalloc ((numLocalChains)*m->numChars, sizeof(MrBFlt));
81 m->inRates = (MrBFlt *) SafeCalloc (m->numGammaCats, sizeof(MrBFlt));
82 m->branchLengths = (MrBFlt *) SafeCalloc (2*numLocalTaxa, sizeof(MrBFlt));
83 m->tiProbIndices = (int *) SafeCalloc (2*numLocalTaxa, sizeof(int));
84 m->inWeights = (MrBFlt *) SafeCalloc (m->numGammaCats*m->nCijkParts, sizeof(MrBFlt));
85 m->bufferIndices = (int *) SafeCalloc (m->nCijkParts, sizeof(int));
86 m->eigenIndices = (int *) SafeCalloc (m->nCijkParts, sizeof(int));
87 m->childBufferIndices = (int *) SafeCalloc (m->nCijkParts, sizeof(int));
88 m->childTiProbIndices = (int *) SafeCalloc (m->nCijkParts, sizeof(int));
89 m->cumulativeScaleIndices = (int *) SafeCalloc (m->nCijkParts, sizeof(int));
92 if (m->numStates != m->numModelStates)
93 numPartAmbigTips = numLocalTaxa;
96 for (i=0; i<numLocalTaxa; i++)
98 if (m->isPartAmbig[i] == YES)
103 if (beagleResourceNumber >= 0 && beagleResourceNumber != 99)
105 resource = beagleResourceNumber;
106 beagleResourceCount = 1;
108 else if (beagleResourceCount != 0)
110 resource = beagleResource[beagleInstanceCount % beagleResourceCount];
112 preferedFlags = beagleFlags;
116 if (beagleScalingScheme == MB_BEAGLE_SCALE_ALWAYS)
117 requiredFlags |= BEAGLE_FLAG_SCALERS_LOG; //BEAGLE_FLAG_SCALERS_RAW;
119 /* TODO: allocate fewer buffers when nCijkParts > 1 */
120 /* create beagle instance */
121 m->beagleInstance = beagleCreateInstance(numLocalTaxa,
122 m->numCondLikes * m->nCijkParts,
123 numLocalTaxa - numPartAmbigTips,
126 (numLocalChains + 1) * m->nCijkParts,
127 m->numTiProbs*m->nCijkParts,
129 m->numScalers * m->nCijkParts,
130 (beagleResourceCount == 0 ? NULL : &resource),
131 (beagleResourceCount == 0 ? 0 : 1),
136 if (m->beagleInstance < 0)
138 MrBayesPrint ("%s Failed to start BEAGLE instance\n", spacer);
143 MrBayesPrint ("\n%s Using BEAGLE resource %i for division %d:", spacer, details.resourceNumber, division+1);
144 # if defined (THREADS_ENABLED)
145 MrBayesPrint (" (%s)\n", (tryToUseThreads ? "threaded" : "non-threaded"));
147 MrBayesPrint (" (non-threaded)\n");
149 MrBayesPrint ("%s Rsrc Name : %s\n", spacer, details.resourceName);
150 MrBayesPrint ("%s Impl Name : %s\n", spacer, details.implName);
151 MrBayesPrint ("%s Flags:", spacer);
152 BeaglePrintFlags(details.flags);
154 beagleInstanceCount++;
157 /* initialize tip data */
158 inStates = (int *) SafeMalloc (m->numChars * sizeof(int));
161 inPartials = (double *) SafeMalloc (m->numChars * m->numModelStates * sizeof(double));
164 for (i=0; i<numLocalTaxa; i++)
166 if (m->isPartAmbig[i] == NO)
168 charBits = m->parsSets[i];
169 for (c=0; c<m->numChars; c++)
171 for (s=j=0; s<m->numModelStates; s++)
173 if (IsBitSet(s, charBits))
179 if (j == m->numModelStates)
183 charBits += m->nParsIntsPerSite;
185 beagleSetTipStates(m->beagleInstance, i, inStates);
187 else /* if (m->isPartAmbig == YES) */
190 charBits = m->parsSets[i];
191 for (c=0; c<m->numChars; c++)
193 for (s=0; s<m->numModelStates; s++)
195 if (IsBitSet(s%m->numStates, charBits))
196 inPartials[k++] = 1.0;
198 inPartials[k++] = 0.0;
200 charBits += m->nParsIntsPerSite;
202 beagleSetTipPartials(m->beagleInstance, i, inPartials);
212 /*-----------------------------------------------------------------
214 | LaunchBEAGLELogLikeForDivision: calculate the log likelihood
215 | of the new state of the chain for a single division
217 -----------------------------------------------------------------*/
218 void LaunchBEAGLELogLikeForDivision(int chain, int d, ModelInfo* m, Tree* tree, MrBFlt* lnL)
220 int i, rescaleFreqNew;
224 if (beagleScalingScheme == MB_BEAGLE_SCALE_ALWAYS)
227 # if defined (DEBUG_MB_BEAGLE_FLOW)
228 MrBayesPrint ("ALWAYS RESCALING\n");
230 /* Flip and copy or reset site scalers */
231 FlipSiteScalerSpace(m, chain);
232 if (m->upDateAll == YES) {
233 for (i=0; i<m->nCijkParts; i++) {
234 beagleResetScaleFactors(m->beagleInstance, m->siteScalerIndex[chain] + i);
238 CopySiteScalers(m, chain);
240 TreeTiProbs_Beagle(tree, d, chain);
241 TreeCondLikes_Beagle(tree, d, chain);
242 TreeLikelihood_Beagle(tree, d, chain, lnL, (chainId[chain] % chainParams.numChains));
245 { /* MB_BEAGLE_SCALE_DYNAMIC */
247 /* This flag is only valid within this block */
248 m->rescaleBeagleAll = NO;
249 TreeTiProbs_Beagle(tree, d, chain);
250 if (m->succesCount[chain] > 1000)
252 m->succesCount[chain] = 10;
253 m->rescaleFreq[chain]++; /* increase rescaleFreq independent of whether we accept or reject new state */
254 m->rescaleFreqOld = rescaleFreqNew = m->rescaleFreq[chain];
255 for (i=0; i<tree->nIntNodes; i++)
257 p = tree->intDownPass[i];
258 if (p->upDateCl == YES) {
259 /* flip to the new workspace since TreeCondLikes_Beagle_Rescale_All() does not do it for
260 (p->upDateCl == YES) since it assumes that TreeCondLikes_Beagle_No_Rescale() did it */
261 FlipCondLikeSpace (m, chain, p->index);
267 if (beagleScalingFrequency != 0 &&
268 m->beagleComputeCount[chain] % (beagleScalingFrequency) == 0)
270 m->rescaleFreqOld = rescaleFreqNew = m->rescaleFreq[chain];
271 for (i=0; i<tree->nIntNodes; i++)
273 p = tree->intDownPass[i];
274 if (p->upDateCl == YES) {
275 /* flip to the new workspace since TreeCondLikes_Beagle_Rescale_All() does not do it for
276 (p->upDateCl == YES) since it assumes that TreeCondLikes_Beagle_No_Rescale() did it */
277 FlipCondLikeSpace (m, chain, p->index);
283 TreeCondLikes_Beagle_No_Rescale(tree, d, chain);
285 /* Check if likelihood is valid */
286 if (TreeLikelihood_Beagle(tree, d, chain, lnL, (chainId[chain] % chainParams.numChains)) == BEAGLE_ERROR_FLOATING_POINT)
288 m->rescaleFreqOld = rescaleFreqNew = m->rescaleFreq[chain];
289 if (rescaleFreqNew > 1 && m->succesCount[chain] < 40)
291 if (m->succesCount[chain] < 10)
293 if (m->succesCount[chain] < 4)
295 rescaleFreqNew -= rescaleFreqNew >> 3; /* <== we cut up to 12,5% of rescaleFreq */
296 if (m->succesCount[chain] < 2)
298 rescaleFreqNew -= rescaleFreqNew >> 3;
299 /* 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*/
300 /* 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 */
301 /* 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*/
302 m->rescaleFreqOld -= m->rescaleFreqOld >> 3;
304 m->rescaleFreqOld -= m->rescaleFreqOld >> 3;
306 m->rescaleFreqOld = (m->rescaleFreqOld ? m->rescaleFreqOld:1);
307 m->recalculateScalers = YES;
312 rescaleFreqNew = (rescaleFreqNew ? rescaleFreqNew : 1);
314 m->succesCount[chain] = 0;
316 # if defined (DEBUG_MB_BEAGLE_FLOW)
317 MrBayesPrint ("NUMERICAL RESCALING\n");
320 m->rescaleBeagleAll = YES;
321 FlipSiteScalerSpace(m, chain);
322 isScalerNode = m->isScalerNode[chain];
324 ResetScalersPartition (isScalerNode, tree, rescaleFreqNew);
325 for (i=0; i<m->nCijkParts; i++)
327 beagleResetScaleFactors(m->beagleInstance, m->siteScalerIndex[chain] + i);
330 TreeCondLikes_Beagle_Rescale_All (tree, d, chain);
331 if (TreeLikelihood_Beagle(tree, d, chain, lnL, (chainId[chain] % chainParams.numChains)) == BEAGLE_ERROR_FLOATING_POINT)
333 if (rescaleFreqNew > 1)
335 /* Swap back scalers which were swapped in TreeCondLikes_Beagle_Rescale_All() */
336 for (i=0; i<tree->nIntNodes; i++)
338 p = tree->intDownPass[i];
339 if (isScalerNode[p->index] == YES)
340 FlipNodeScalerSpace (m, chain, p->index);
342 rescaleFreqNew -= rescaleFreqNew >> 3; /* <== we cut up to 12,5% of rescaleFreq */
343 rescaleFreqNew--; /* we cut extra 1 of rescaleFreq */
347 m->rescaleFreq[chain] = rescaleFreqNew;
351 /* Count number of evaluations */
352 m->beagleComputeCount[chain]++;
356 void recalculateScalers(int chain)
358 int i, d, rescaleFreqNew;
363 for (d=0; d<numCurrentDivisions; d++)
365 m = &modelSettings[d];
366 if (m->recalculateScalers == YES)
368 m->recalculateScalers = NO;
369 tree = GetTree(m->brlens, chain, state[chain]);
371 rescaleFreqNew = m->rescaleFreq[chain];
372 isScalerNode = m->isScalerNode[chain];
374 ResetScalersPartition (isScalerNode, tree, rescaleFreqNew);
375 for (i=0; i<m->nCijkParts; i++) {
376 beagleResetScaleFactors(m->beagleInstance, m->siteScalerIndex[chain] + i);
378 /* here it does not matter if we flip CL space or not */
379 TreeCondLikes_Beagle_Rescale_All (tree, d, chain);
385 void BeagleAddGPUDevicesToList(int **newResourceList, int *beagleResourceCount)
387 BeagleResourceList* beagleResources;
390 beagleResources = beagleGetResourceList();
391 if (*newResourceList == NULL) {
392 *newResourceList = (int*) SafeCalloc(sizeof(int), beagleResources->length);
395 for (i = 0; i < beagleResources->length; i++) {
396 if (beagleResources->list[i].supportFlags & BEAGLE_FLAG_PROCESSOR_GPU) {
397 (*newResourceList)[gpuCount] = i;
401 *beagleResourceCount = gpuCount;
405 void BeagleRemoveGPUDevicesFromList(int **beagleResource, int *beagleResourceCount)
407 *beagleResourceCount = 0;
413 | BeaglePrintResources: outputs the available BEAGLE resources
416 void BeaglePrintResources()
419 BeagleResourceList* beagleResources;
421 beagleResources = beagleGetResourceList();
422 MrBayesPrint ("%s Available resources reported by beagle library:\n", spacer);
423 for (i=0; i<beagleResources->length; i++)
425 MrBayesPrint ("\tResource %i:\n", i);
426 MrBayesPrint ("\tName: %s\n", beagleResources->list[i].name);
429 MrBayesPrint ("\tDesc: %s\n", beagleResources->list[i].description);
431 MrBayesPrint ("\tFlags:");
432 BeaglePrintFlags(beagleResources->list[i].supportFlags);
433 MrBayesPrint ("\n\n");
435 MrBayesPrint ("%s BEAGLE vesion: %s\n", spacer, beagleGetVersion());
439 int BeagleCheckFlagCompatability(long inFlags)
441 if (inFlags & BEAGLE_FLAG_PROCESSOR_GPU) {
442 if (inFlags & BEAGLE_FLAG_VECTOR_SSE) {
443 MrBayesPrint ("%s Simultaneous use of GPU and SSE not available.\n", spacer);
446 if (inFlags & BEAGLE_FLAG_THREADING_OPENMP) {
447 MrBayesPrint ("%s Simultaneous use of GPU and OpenMP not available.\n", spacer);
456 /*-------------------
458 | BeaglePrintFlags: outputs beagle instance details
460 ______________________*/
461 void BeaglePrintFlags(long inFlags)
464 char *names[] = { "PROCESSOR_CPU",
470 "COMPUTATION_ASYNCH",
485 long flags[] = { BEAGLE_FLAG_PROCESSOR_CPU,
486 BEAGLE_FLAG_PROCESSOR_GPU,
487 BEAGLE_FLAG_PROCESSOR_FPGA,
488 BEAGLE_FLAG_PROCESSOR_CELL,
489 BEAGLE_FLAG_PRECISION_DOUBLE,
490 BEAGLE_FLAG_PRECISION_SINGLE,
491 BEAGLE_FLAG_COMPUTATION_ASYNCH,
492 BEAGLE_FLAG_COMPUTATION_SYNCH,
493 BEAGLE_FLAG_EIGEN_REAL,
494 BEAGLE_FLAG_EIGEN_COMPLEX,
495 BEAGLE_FLAG_SCALING_MANUAL,
496 BEAGLE_FLAG_SCALING_AUTO,
497 BEAGLE_FLAG_SCALING_ALWAYS,
498 BEAGLE_FLAG_SCALING_DYNAMIC,
499 BEAGLE_FLAG_SCALERS_RAW,
500 BEAGLE_FLAG_SCALERS_LOG,
501 BEAGLE_FLAG_VECTOR_NONE,
502 BEAGLE_FLAG_VECTOR_SSE,
503 BEAGLE_FLAG_THREADING_NONE,
504 BEAGLE_FLAG_THREADING_OPENMP
507 for (i=k=0; i<20; i++)
509 if (inFlags & flags[i])
511 if (k%4 == 0 && k > 0)
512 MrBayesPrint ("\n%s ", spacer);
513 MrBayesPrint (" %s", names[i]);
520 #if defined(THREADS_ENABLED)
521 void *LaunchThreadLogLikeForDivision(void *arguments)
525 LaunchStruct* launchStruct;
527 launchStruct = (LaunchStruct*) arguments;
528 chain = launchStruct->chain;
529 d = launchStruct->division;
530 lnL = launchStruct->lnL;
531 LaunchLogLikeForDivision(chain, d, lnL);
536 MrBFlt LaunchLogLikeForAllDivisionsInParallel(int chain)
541 LaunchStruct* launchValues;
548 /* TODO Initialize only once */
549 threads = (pthread_t*) SafeMalloc(sizeof(pthread_t) * numCurrentDivisions);
550 launchValues = (LaunchStruct*) SafeMalloc(sizeof(LaunchStruct) * numCurrentDivisions);
551 wait = (int*) SafeMalloc(sizeof(int) * numCurrentDivisions);
553 /* Cycle through divisions and recalculate tis and cond likes as necessary. */
554 /* Code below does not try to avoid recalculating ti probs for divisions */
555 /* that could share ti probs with other divisions. */
556 for (d=0; d<numCurrentDivisions; d++)
559 # if defined (BEST_MPI_ENABLED)
560 if (isDivisionActive[d] == NO)
563 m = &modelSettings[d];
565 if (m->upDateCl == YES)
567 launchValues[d].chain = chain;
568 launchValues[d].division = d;
569 launchValues[d].lnL = &(m->lnLike[2*chain + state[chain]]);
571 threadError = pthread_create(&threads[d], NULL,
572 LaunchThreadLogLikeForDivision,
573 (void*) &launchValues[d]);
574 assert (0 == threadError);
583 for (d = 0; d < numCurrentDivisions; d++)
588 threadError = pthread_join(threads[d], NULL);
589 assert (0 == threadError);
591 m = &modelSettings[d];
592 chainLnLike += m->lnLike[2*chain + state[chain]];
595 /* TODO Free these once */
605 int ScheduleLogLikeForAllDivisions()
608 int divisionsToLaunch = 0;
611 if (numCurrentDivisions < 2) {
615 for (d=0; d<numCurrentDivisions; d++) {
616 m = &modelSettings[d];
617 if (m->upDateCl == YES) {
622 return (divisionsToLaunch > 1);
626 /*----------------------------------------------------------------
628 | TreeCondLikes_Beagle: This routine updates all conditional
629 | (partial) likelihoods of a beagle instance while doing no rescaling.
630 | That potentialy can make final liklihood bad then calculation with rescaling needs to be done.
632 -----------------------------------------------------------------*/
633 int TreeCondLikes_Beagle_No_Rescale (Tree *t, int division, int chain)
635 int i, j, cumulativeScaleIndex;
636 BeagleOperation operations;
639 unsigned chil1Step, chil2Step;
642 m = &modelSettings[division];
643 isScalerNode = m->isScalerNode[chain];
645 for (i=0; i<t->nIntNodes; i++)
647 p = t->intDownPass[i];
649 /* check if conditional likelihoods need updating */
650 if (p->upDateCl == YES)
652 /* flip to the new workspace */
653 FlipCondLikeSpace (m, chain, p->index);
655 /* update conditional likelihoods */
656 operations.destinationPartials = m->condLikeIndex[chain][p->index ];
657 operations.child1Partials = m->condLikeIndex[chain][p->left->index ];
658 operations.child1TransitionMatrix = m->tiProbsIndex [chain][p->left->index ];
659 operations.child2Partials = m->condLikeIndex[chain][p->right->index];
660 operations.child2TransitionMatrix = m->tiProbsIndex [chain][p->right->index];
662 /* All partials for tips are the same across omega categories, thus we are doing the following two if statments.*/
663 if (p->left->left== NULL)
668 if (p->right->left== NULL)
673 operations.destinationScaleWrite = BEAGLE_OP_NONE;
674 cumulativeScaleIndex = BEAGLE_OP_NONE;
675 if (isScalerNode[p->index] == YES)
677 operations.destinationScaleRead = m->nodeScalerIndex[chain][p->index];
681 operations.destinationScaleRead = BEAGLE_OP_NONE;
684 for (j=0; j<m->nCijkParts; j++)
686 beagleUpdatePartials(m->beagleInstance,
689 cumulativeScaleIndex);
691 operations.destinationPartials++;
692 operations.child1Partials+=chil1Step;
693 operations.child1TransitionMatrix++;
694 operations.child2Partials+=chil2Step;
695 operations.child2TransitionMatrix++;
697 if (isScalerNode[p->index] == YES)
698 operations.destinationScaleRead++;
707 /*----------------------------------------------------------------
709 | TreeCondLikes_Beagle: This routine updates all conditional
710 | (partial) likelihoods of a beagle instance while rescaling at every node.
711 | Note: all nodes get recalculated, not only tached by move.
713 -----------------------------------------------------------------*/
714 int TreeCondLikes_Beagle_Rescale_All (Tree *t, int division, int chain)
716 int i, j, cumulativeScaleIndex;
717 BeagleOperation operations;
720 unsigned chil1Step, chil2Step;
723 m = &modelSettings[division];
724 isScalerNode = m->isScalerNode[chain];
726 for (i=0; i<t->nIntNodes; i++)
728 p = t->intDownPass[i];
730 if (p->upDateCl == NO) {
732 /* flip to the new workspace */
733 FlipCondLikeSpace (m, chain, p->index);
736 /* update conditional likelihoods */
737 operations.destinationPartials = m->condLikeIndex[chain][p->index ];
738 operations.child1Partials = m->condLikeIndex[chain][p->left->index ];
739 operations.child1TransitionMatrix = m->tiProbsIndex [chain][p->left->index ];
740 operations.child2Partials = m->condLikeIndex[chain][p->right->index];
741 operations.child2TransitionMatrix = m->tiProbsIndex [chain][p->right->index];
743 /* All partials for tips are the same across omega catigoris, thus we are doing the following two if statments.*/
744 if (p->left->left== NULL)
749 if (p->right->left== NULL)
754 operations.destinationScaleRead = BEAGLE_OP_NONE;
755 if (isScalerNode[p->index] == YES)
757 FlipNodeScalerSpace (m, chain, p->index);
758 operations.destinationScaleWrite = m->nodeScalerIndex[chain][p->index];
759 cumulativeScaleIndex = m->siteScalerIndex[chain];
763 operations.destinationScaleWrite = BEAGLE_OP_NONE;
764 cumulativeScaleIndex = BEAGLE_OP_NONE;
769 for (j=0; j<m->nCijkParts; j++)
771 beagleUpdatePartials(m->beagleInstance,
774 cumulativeScaleIndex);
776 operations.destinationPartials++;
777 operations.child1Partials+=chil1Step;
778 operations.child1TransitionMatrix++;
779 operations.child2Partials+=chil2Step;
780 operations.child2TransitionMatrix++;
782 if (isScalerNode[p->index] == YES) {
783 operations.destinationScaleWrite++;
784 cumulativeScaleIndex++;
794 /*----------------------------------------------------------------
796 | TreeCondLikes_Beagle: This routine updates all conditional
797 | (partial) likelihoods of a beagle instance.
799 -----------------------------------------------------------------*/
800 int TreeCondLikes_Beagle (Tree *t, int division, int chain)
802 int i, j, destinationScaleRead, cumulativeScaleIndex;
803 BeagleOperation operations;
806 unsigned chil1Step, chil2Step;
808 m = &modelSettings[division];
810 for (i=0; i<t->nIntNodes; i++)
812 p = t->intDownPass[i];
814 /* check if conditional likelihoods need updating */
815 if (p->upDateCl == YES)
817 /* remove old scalers */
818 if (m->scalersSet[chain][p->index] == YES && m->upDateAll == NO)
820 destinationScaleRead = m->nodeScalerIndex[chain][p->index];
821 cumulativeScaleIndex = m->siteScalerIndex[chain];
822 for (j=0; j<m->nCijkParts; j++)
824 beagleRemoveScaleFactors(m->beagleInstance,
825 &destinationScaleRead,
827 cumulativeScaleIndex);
828 destinationScaleRead++;
829 cumulativeScaleIndex++;
833 /* flip to the new workspace */
834 FlipCondLikeSpace (m, chain, p->index);
835 FlipNodeScalerSpace (m, chain, p->index);
836 m->scalersSet[chain][p->index] = NO;
838 /* update conditional likelihoods */
839 operations.destinationPartials = m->condLikeIndex[chain][p->index ];
840 operations.child1Partials = m->condLikeIndex[chain][p->left->index ];
841 operations.child1TransitionMatrix = m->tiProbsIndex [chain][p->left->index ];
842 operations.child2Partials = m->condLikeIndex[chain][p->right->index];
843 operations.child2TransitionMatrix = m->tiProbsIndex [chain][p->right->index];
845 /* All partials for tips are the same across omega catigoris, thus we are doing the following two if statments.*/
846 if (p->left->left== NULL && p->left->right== NULL)
851 if (p->right->left== NULL && p->right->right== NULL)
856 if (p->scalerNode == YES)
858 m->scalersSet[chain][p->index] = YES;
859 operations.destinationScaleWrite = m->nodeScalerIndex[chain][p->index];
860 cumulativeScaleIndex = m->siteScalerIndex[chain];
864 operations.destinationScaleWrite = BEAGLE_OP_NONE;
865 cumulativeScaleIndex = BEAGLE_OP_NONE;
867 operations.destinationScaleRead = BEAGLE_OP_NONE;
869 for (j=0; j<m->nCijkParts; j++)
871 beagleUpdatePartials(m->beagleInstance,
874 cumulativeScaleIndex);
876 operations.destinationPartials++;
877 operations.child1Partials+=chil1Step;
878 operations.child1TransitionMatrix++;
879 operations.child2Partials+=chil2Step;
880 operations.child2TransitionMatrix++;
881 if (p->scalerNode == YES)
883 operations.destinationScaleWrite++;
884 cumulativeScaleIndex++;
895 /**---------------------------------------------------------------------------
897 | TreeLikelihood_Beagle: Accumulate the log likelihoods calculated by Beagle
900 ---------------------------------------- -------------------------------------*/
901 int TreeLikelihood_Beagle (Tree *t, int division, int chain, MrBFlt *lnL, int whichSitePats)
903 int i, j, c = 0, nStates, hasPInvar, beagleReturn;
904 MrBFlt *swr, s01, s10, probOn, probOff, covBF[40], pInvar=0.0, *bs, freq, likeI, lnLikeI, diff, *omegaCatFreq;
905 CLFlt *clInvar=NULL, *nSitesOfPat;
906 double *nSitesOfPat_Beagle;
911 # if defined (MB_PRINT_DYNAMIC_RESCALE_FAIL_STAT)
912 static unsigned countBeagleDynamicFail=0;
913 static unsigned countALL=0;
919 /* find model settings and nStates, pInvar, invar cond likes */
920 m = &modelSettings[division];
922 nStates = m->numModelStates;
923 if (m->pInvar == NULL)
930 pInvar = *(GetParamVals (m->pInvar, chain, state[chain]));
931 clInvar = m->invCondLikes;
934 /* find base frequencies */
935 bs = GetParamSubVals (m->stateFreq, chain, state[chain]);
937 /* if covarion model, adjust base frequencies */
938 if (m->switchRates != NULL)
940 /* find the stationary frequencies */
941 swr = GetParamVals(m->switchRates, chain, state[chain]);
944 probOn = s01 / (s01 + s10);
945 probOff = 1.0 - probOn;
947 /* now adjust the base frequencies; on-state stored first in cond likes */
948 for (j=0; j<nStates/2; j++)
950 covBF[j] = bs[j] * probOn;
951 covBF[j+nStates/2] = bs[j] * probOff;
954 /* finally set bs pointer to adjusted values */
958 /* TODO Really only need to check if state frequencies have changed */
959 if (m->upDateCijk == YES)
961 /* set base frequencies in BEAGLE instance */
962 for (i=0; i<m->nCijkParts; i++)
963 beagleSetStateFrequencies(m->beagleInstance,
964 m->cijkIndex[chain] + i,
968 /* find category frequencies */
970 freq = 1.0 / m->numGammaCats;
972 freq = (1.0 - pInvar) / m->numGammaCats;
974 /* TODO: cat weights only need to be set when they change */
975 /* set category frequencies in beagle instance */
976 if (m->numOmegaCats > 1)
978 omegaCatFreq = GetParamSubVals(m->omega, chain, state[chain]);
979 for (i=0; i<m->nCijkParts; i++)
981 for (j=0; j<m->numGammaCats; j++)
982 m->inWeights[j] = freq * omegaCatFreq[i];
983 beagleSetCategoryWeights(m->beagleInstance,
984 m->cijkIndex[chain] + i,
988 else if (hasPInvar == YES)
990 for (i=0; i<m->numGammaCats; i++)
991 m->inWeights[i] = freq;
992 beagleSetCategoryWeights(m->beagleInstance,
997 /* find nSitesOfPat */
998 nSitesOfPat = numSitesOfPat + (whichSitePats*numCompressedChars) + m->compCharStart;
1000 /* TODO: pattern weights only need to be set when they change */
1001 /* set pattern weights in beagle instance if using dynamic reweighting */
1002 if (chainParams.weightScheme[0] + chainParams.weightScheme[1] > ETA)
1004 nSitesOfPat_Beagle = (double *) SafeMalloc (m->numChars * sizeof(double));
1005 for (c=0; c<m->numChars; c++)
1006 nSitesOfPat_Beagle[c] = numSitesOfPat[m->compCharStart + c];
1007 beagleSetPatternWeights(m->beagleInstance,
1008 nSitesOfPat_Beagle);
1009 SafeFree ((void **)(&nSitesOfPat_Beagle));
1012 /* find root log likelihoods and scalers */
1013 for (i=0; i<m->nCijkParts; i++)
1015 m->bufferIndices[i] = m->condLikeIndex[chain][p->index] + i;
1016 m->eigenIndices[i] = m->cijkIndex[chain] + i;
1017 m->cumulativeScaleIndices[i] = m->siteScalerIndex[chain] + i;
1018 if (t->isRooted == NO)
1020 m->childBufferIndices[i] = m->condLikeIndex [chain][p->anc->index];
1021 m->childTiProbIndices[i] = m->tiProbsIndex [chain][p->index] + i;
1028 /* get root log likelihoods */
1029 if (t->isRooted == YES)
1031 beagleReturn = beagleCalculateRootLogLikelihoods(m->beagleInstance,
1035 m->cumulativeScaleIndices,
1042 beagleReturn = beagleCalculateEdgeLogLikelihoods(m->beagleInstance,
1044 m->childBufferIndices,
1045 m->childTiProbIndices,
1050 m->cumulativeScaleIndices,
1056 # if defined (MB_PRINT_DYNAMIC_RESCALE_FAIL_STAT)
1059 if (beagleReturn == BEAGLE_ERROR_FLOATING_POINT)
1061 # if defined (MB_PRINT_DYNAMIC_RESCALE_FAIL_STAT)
1062 countBeagleDynamicFail++;
1063 MrBayesPrint ("DEBUG INFO (not an error) countBeagleDynamicFail:%d countALL:%d\n", countBeagleDynamicFail, countALL);
1065 return beagleReturn;
1067 assert (beagleReturn == BEAGLE_SUCCESS);
1068 m->succesCount[chain]++;
1070 /* accumulate logs across sites */
1071 if (hasPInvar == NO)
1073 if (m->dataType == RESTRICTION)
1075 beagleGetSiteLogLikelihoods(m->beagleInstance, m->logLikelihoods);
1078 for (c=0; c<m->numDummyChars; c++)
1080 pUnobserved += exp((double)m->logLikelihoods[c]);
1082 /* correct for absent characters */
1083 (*lnL) -= log (1-pUnobserved) * (m->numUncompressedChars);
1084 for (; c<m->numChars; c++)
1086 (*lnL) += m->logLikelihoods[c] * nSitesOfPat[c];
1089 /* already done, just check for numerical errors */
1090 assert ((*lnL) == (*lnL));
1094 /* has invariable category */
1095 beagleGetSiteLogLikelihoods(m->beagleInstance,
1098 for (c=0; c<m->numChars; c++)
1101 for (j=0; j<nStates; j++)
1102 likeI += (*(clInvar++)) * bs[j];
1105 lnLikeI = log(likeI * pInvar);
1106 diff = lnLikeI - m->logLikelihoods[c];
1111 (*lnL) += m->logLikelihoods[c] * nSitesOfPat[c];
1112 else if (diff > 200.0)
1113 (*lnL) += lnLikeI * nSitesOfPat[c];
1116 (*lnL) += (m->logLikelihoods[c] + log(1.0 + exp(diff))) * nSitesOfPat[c];
1119 /* check for numerical errors */
1120 assert ((*lnL) == (*lnL));
1128 /*----------------------------------------------------------------
1130 | TreeTiProbs_Beagle: This routine updates all transition
1131 | probability matrices of a beagle instance.
1133 -----------------------------------------------------------------*/
1134 int TreeTiProbs_Beagle (Tree *t, int division, int chain)
1137 MrBFlt correctionFactor, theRate, baseRate, *catRate, length;
1141 /* get model settings */
1142 m = &modelSettings[division];
1144 /* find the correction factor to make branch lengths
1145 in terms of expected number of substitutions per character */
1146 correctionFactor = 1.0;
1147 if (m->dataType == DNA || m->dataType == RNA)
1149 if (m->nucModelId == NUCMODEL_DOUBLET)
1150 correctionFactor = 2.0;
1151 else if (m->nucModelId == NUCMODEL_CODON)
1152 correctionFactor = 3.0;
1155 /* get rate multipliers (for gamma & partition specific rates) */
1157 baseRate = GetRate (division, chain);
1159 /* compensate for invariable sites if appropriate */
1160 if (m->pInvar != NULL)
1161 baseRate /= (1.0 - (*GetParamVals(m->pInvar, chain, state[chain])));
1163 /* get category rates for gamma */
1164 if (m->shape == NULL)
1167 catRate = GetParamSubVals (m->shape, chain, state[chain]);
1169 /* get effective category rates */
1170 for (k=0; k<m->numGammaCats; k++)
1171 m->inRates[k] = baseRate * catRate[k] * correctionFactor;
1173 /* TODO: only need to set category rates when they change */
1174 /* set category rates */
1175 beagleSetCategoryRates(m->beagleInstance, m->inRates);
1177 /* get ti prob indices and branch lengths to update */
1178 for (i=count=0; i<t->nNodes; i++)
1180 p = t->allDownPass[i];
1182 /* check if transition probs need updating */
1183 if (p->upDateTi == YES)
1185 /* flip transition probability */
1186 FlipTiProbsSpace (m, chain, p->index);
1189 if (m->cppEvents != NULL)
1191 length = GetParamSubVals (m->cppEvents, chain, state[chain])[p->index];
1193 else if (m->tk02BranchRates != NULL)
1195 length = GetParamSubVals (m->tk02BranchRates, chain, state[chain])[p->index];
1197 else if (m->igrBranchRates != NULL)
1199 length = GetParamSubVals (m->igrBranchRates, chain, state[chain])[p->index];
1201 else if (m->mixedBrchRates != NULL)
1203 length = GetParamSubVals (m->mixedBrchRates, chain, state[chain])[p->index];
1208 /* numerical errors might ensue if we allow very large or very small branch lengths, which might
1209 occur in relaxed clock models; an elegant solution would be to substitute the stationary
1210 probs and initial probs but for now we truncate lengths at small or large values */
1211 if (length > BRLENS_MAX)
1212 length = BRLENS_MAX;
1213 else if (length < BRLENS_MIN)
1214 length = BRLENS_MIN;
1216 m->branchLengths[count] = length;
1219 m->tiProbIndices[count] = m->tiProbsIndex[chain][p->index];
1224 /* TODO: only need to update branches that have changed */
1225 /* calculate transition probabilities */
1227 for (i=0; i<m->nCijkParts; i++)
1229 beagleUpdateTransitionMatrices(m->beagleInstance,
1230 m->cijkIndex[chain] + i,
1236 for (j=0; j<count; j++)
1237 m->tiProbIndices[j]++;
1241 /* return success */
1245 #endif // BEAGLE_ENABLED
1248 void BeagleNotLinked()
1250 MrBayesPrint ("%s BEAGLE library is not linked to this executable.\n", spacer);
1254 void BeagleThreadsNotLinked()
1256 MrBayesPrint ("%s Pthreads library is not linked to this executable.\n", spacer);