1 /*
2 * MrBayes 3
3 *
4 * (c) 2002-2013
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 'acknowledgments' 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 "bayes.h"
33 #include "mbbeagle.h"
34 #include "mcmc.h"
35 #include "model.h"
36 #include "utils.h"
37
38 // #define DEBUG_MB_BEAGLE_MULTIPART
39 // #define DEBUG_MB_BEAGLE_MULTIPART_SITELNL
40
41 const char* const svnRevisionMbbeagleC = "$Rev$"; /* Revision keyword which is expended/updated by svn on each commit/update */
42
43 /* Functions and variables defined in mcmc.c that are not exported in mcmc.h */
44 void LaunchLogLikeForDivision(int chain, int d, MrBFlt* lnL);
45
46 void FlipCondLikeSpace (ModelInfo *m, int chain, int nodeIndex);
47 void FlipNodeScalerSpace (ModelInfo *m, int chain, int nodeIndex);
48 void FlipSiteScalerSpace (ModelInfo *m, int chain);
49 void FlipTiProbsSpace (ModelInfo *m, int chain, int nodeIndex);
50 void ResetSiteScalers (ModelInfo *m, int chain);
51 void CopySiteScalers (ModelInfo *m, int chain);
52
53
54 extern int *chainId;
55 extern int numLocalChains;
56
57
58 #if defined (BEAGLE_ENABLED)
59 /*------------------------------------------------------------------------
60 |
61 | InitBeagleInstance: create and initialize a beagle instance
62 |
63 -------------------------------------------------------------------------*/
InitBeagleInstance(ModelInfo * m,int division)64 int InitBeagleInstance (ModelInfo *m, int division)
65 {
66 int i, j, k, c, s, *inStates, numPartAmbigTips;
67 double *inPartials;
68 BitsLong *charBits;
69
70 if (m->beagleInstance == -99) {
71 MrBayesPrint ("\n%s Error: Consecutive MCMC runs not currently supported with BEAGLE. Please restart MrBayes.\n\n", spacer);
72 exit(1);
73 }
74
75
76 if (m->useBeagle == NO)
77 return ERROR;
78
79 /* at least one eigen buffer needed */
80 if (m->nCijkParts == 0)
81 m->nCijkParts = 1;
82
83 /* allocate memory used by beagle */
84 m->logLikelihoods = (MrBFlt *) SafeCalloc ((numLocalChains)*m->numChars, sizeof(MrBFlt));
85 m->inRates = (MrBFlt *) SafeCalloc (m->numRateCats, sizeof(MrBFlt));
86 m->branchLengths = (MrBFlt *) SafeCalloc (2*numLocalTaxa, sizeof(MrBFlt));
87 m->tiProbIndices = (int *) SafeCalloc (2*numLocalTaxa, sizeof(int));
88 m->inWeights = (MrBFlt *) SafeCalloc (m->numRateCats*m->nCijkParts, sizeof(MrBFlt));
89 m->bufferIndices = (int *) SafeCalloc (m->nCijkParts, sizeof(int));
90 m->eigenIndices = (int *) SafeCalloc (m->nCijkParts, sizeof(int));
91 m->childBufferIndices = (int *) SafeCalloc (m->nCijkParts, sizeof(int));
92 m->childTiProbIndices = (int *) SafeCalloc (m->nCijkParts, sizeof(int));
93 m->cumulativeScaleIndices = (int *) SafeCalloc (m->nCijkParts, sizeof(int));
94 m->operations = (BeagleOperation *) SafeCalloc(m->numCondLikes*m->nCijkParts, sizeof(BeagleOperation));
95 m->scaleFactorsOps = (int *) SafeCalloc (2*numLocalTaxa*m->nCijkParts, sizeof(int));
96
97 numPartAmbigTips = 0;
98 if (m->numStates != m->numModelStates)
99 numPartAmbigTips = numLocalTaxa;
100 else
101 {
102 for (i=0; i<numLocalTaxa; i++)
103 {
104 #if defined(BEAGLE_V3_ENABLED)
105 if (m->isPartAmbig[i] == YES || beagleAllFloatTips)
106 #else
107 if (m->isPartAmbig[i] == YES)
108 #endif
109 {
110 numPartAmbigTips++;
111 }
112 }
113 }
114
115 m->beagleInstance = createBeagleInstance(m, m->nCijkParts, m->numRateCats, m->numModelStates, m->numCondLikes, m->numScalers, m->numChars, m->numTiProbs, numPartAmbigTips, division);
116
117 if (m->beagleInstance < 0)
118 return ERROR;
119
120 /* initialize tip data */
121 inStates = (int *) SafeMalloc (m->numChars * sizeof(int));
122 if (!inStates)
123 return ERROR;
124 inPartials = (double *) SafeMalloc (m->numChars * m->numModelStates * sizeof(double));
125 if (!inPartials)
126 return ERROR;
127 for (i=0; i<numLocalTaxa; i++)
128 {
129 #if defined(BEAGLE_V3_ENABLED)
130 if (m->isPartAmbig[i] == NO && !(beagleAllFloatTips))
131 #else
132 if (m->isPartAmbig[i] == NO)
133 #endif
134 {
135 charBits = m->parsSets[i];
136 for (c=0; c<m->numChars; c++)
137 {
138 for (s=j=0; s<m->numModelStates; s++)
139 {
140 if (IsBitSet(s, charBits))
141 {
142 inStates[c] = s;
143 j++;
144 }
145 }
146 if (j == m->numModelStates)
147 inStates[c] = j;
148 else
149 assert (j==1);
150 charBits += m->nParsIntsPerSite;
151 }
152 beagleSetTipStates(m->beagleInstance, i, inStates);
153 }
154 else /* if (m->isPartAmbig == YES) */
155 {
156 k = 0;
157 charBits = m->parsSets[i];
158 for (c=0; c<m->numChars; c++)
159 {
160 for (s=0; s<m->numModelStates; s++)
161 {
162 if (IsBitSet(s%m->numStates, charBits))
163 inPartials[k++] = 1.0;
164 else
165 inPartials[k++] = 0.0;
166 }
167 charBits += m->nParsIntsPerSite;
168 }
169 beagleSetTipPartials(m->beagleInstance, i, inPartials);
170 }
171 }
172 free (inStates);
173 free (inPartials);
174
175 return NO_ERROR;
176 }
177
178
createBeagleInstance(ModelInfo * m,int nCijkParts,int numRateCats,int numModelStates,int numCondLikes,int numScalers,int numChars,int numTiProbs,int numPartAmbigTips,int division)179 int createBeagleInstance(ModelInfo *m, int nCijkParts, int numRateCats, int numModelStates, int numCondLikes, int numScalers, int numChars, int numTiProbs, int numPartAmbigTips, int division)
180 {
181 int resource, resourceCount, beagleInstance;
182 BeagleInstanceDetails details;
183 long preferredFlags, requiredFlags;
184
185 #if defined (BEAGLE_V3_ENABLED)
186 int i;
187 Tree* t;
188 #endif
189
190 resourceCount = beagleResourceCount;
191
192 preferredFlags = beagleFlags;
193 requiredFlags = 0L;
194
195 if (beagleScalingScheme == MB_BEAGLE_SCALE_ALWAYS)
196 {
197 requiredFlags |= BEAGLE_FLAG_SCALERS_LOG;
198 }
199
200 if (beagleResourceNumber >= 0 && beagleResourceNumber != 99)
201 {
202 resource = beagleResourceNumber;
203 resourceCount = 1;
204 }
205 else if (beagleResourceCount != 0)
206 {
207 # if defined (MPI_ENABLED)
208 resource = beagleResource[(beagleInstanceCount + proc_id) % beagleResourceCount];
209 # else
210 resource = beagleResource[beagleInstanceCount % beagleResourceCount];
211 # endif
212 }
213 #if defined (BEAGLE_V3_ENABLED)
214 else if (beagleResourceNumber == 99)
215 {
216 int numInstancePartitions = 1;
217 if (division == -1)
218 {
219 numInstancePartitions = numCurrentDivisions;
220 }
221
222 MrBayesPrint ("\n%s Running benchmarks to automatically select fastest BEAGLE resource... ", spacer);
223
224 long benchmarkFlags = BEAGLE_BENCHFLAG_SCALING_NONE;
225 if (beagleScalingScheme == MB_BEAGLE_SCALE_ALWAYS)
226 {
227 benchmarkFlags = BEAGLE_BENCHFLAG_SCALING_ALWAYS;
228 }
229 else if (beagleScalingScheme == MB_BEAGLE_SCALE_DYNAMIC)
230 {
231 benchmarkFlags = BEAGLE_BENCHFLAG_SCALING_DYNAMIC;
232 }
233
234 // select fastest resource
235 BeagleBenchmarkedResourceList* rBList;
236 rBList = beagleGetBenchmarkedResourceList(
237 numLocalTaxa,
238 numLocalTaxa - numPartAmbigTips,
239 numModelStates,
240 numChars,
241 numRateCats,
242 NULL,
243 0,
244 preferredFlags,
245 requiredFlags,
246 nCijkParts,
247 numInstancePartitions,
248 0,
249 benchmarkFlags);
250
251 # if defined (DEBUG_MB_BEAGLE_FLOW)
252 fprintf(stdout, "Resource benchmarks:\n");
253 fprintf(stdout, "\ttipCount : %d\n", numLocalTaxa);
254 fprintf(stdout, "\tcompactBufferCount : %d\n", numLocalTaxa - numPartAmbigTips);
255 fprintf(stdout, "\tstateCount : %d\n", numModelStates);
256 fprintf(stdout, "\tpatternCount : %d\n", numChars);
257 fprintf(stdout, "\tcategoryCount : %d\n", numRateCats);
258 fprintf(stdout, "\teigenModelCount : %d\n", nCijkParts);
259 fprintf(stdout, "\tpartitionCount : %d\n", numInstancePartitions);
260 fprintf(stdout, "\tPreferred Flags:");
261 BeaglePrintFlags(preferredFlags);
262 fprintf(stdout, "\n");
263 fprintf(stdout, "\tRequired Flags:");
264 BeaglePrintFlags(requiredFlags);
265 fprintf(stdout, "\n");
266 fprintf(stdout, "\tBenchmark Flags: ");
267 (benchmarkFlags & BEAGLE_BENCHFLAG_SCALING_ALWAYS ? fprintf(stdout, "BEAGLE_BENCHFLAG_SCALING_ALWAYS") : fprintf(stdout, "BEAGLE_BENCHFLAG_SCALING_DYNAMIC"));
268 fprintf(stdout, "\n");
269 fprintf(stdout, "\n");
270
271
272
273 if (rBList != NULL)
274 {
275 for (int i = 0; i < rBList->length; i++)
276 {
277 fprintf(stdout, "\tResource %i:\n\t\tName : %s\n", i, rBList->list[i].name);
278 fprintf(stdout, "\t\tDesc : %s\n", rBList->list[i].description);
279 fprintf(stdout, "\t\tSupport Flags:");
280 BeaglePrintFlags(rBList->list[i].supportFlags);
281 fprintf(stdout, "\n");
282 fprintf(stdout, "\t\tRequired Flags:");
283 BeaglePrintFlags(rBList->list[i].requiredFlags);
284 fprintf(stdout, "\n");
285 fprintf(stdout, "\t\tBenchmark Results:\n");
286 fprintf(stdout, "\t\t\tNmbr : %d\n", rBList->list[i].number);
287 fprintf(stdout, "\t\t\tImpl : %s\n", rBList->list[i].implName);
288 fprintf(stdout, "\t\t\tFlags:");
289 BeaglePrintFlags(rBList->list[i].benchedFlags);
290 fprintf(stdout, "\n");
291 fprintf(stdout, "\t\t\tPerf : %.4f ms (%.2fx CPU)\n", rBList->list[i].benchmarkResult, rBList->list[i].performanceRatio);
292 }
293 }
294 fprintf(stdout, "\n");
295 # endif
296
297 if (rBList != NULL)
298 {
299 double fastestTime = rBList->list[0].benchmarkResult;
300 resource = rBList->list[0].number;
301 resourceCount = 1;
302
303 for (int i = 1; i < rBList->length; i++)
304 {
305 if (rBList->list[i].benchmarkResult < fastestTime)
306 {
307 fastestTime = rBList->list[i].benchmarkResult;
308 resource = rBList->list[i].number;
309 }
310 }
311 }
312 }
313 #endif /* BEAGLE_V3_ENABLED */
314
315 /* TODO: allocate fewer buffers when nCijkParts > 1 */
316 /* create beagle instance */
317 beagleInstance = beagleCreateInstance(numLocalTaxa,
318 numCondLikes * nCijkParts,
319 numLocalTaxa - numPartAmbigTips,
320 numModelStates,
321 numChars,
322 ((numLocalChains + 1) * nCijkParts) * numCurrentDivisions,
323 numTiProbs*nCijkParts,
324 numRateCats,
325 numScalers * nCijkParts,
326 (resourceCount == 0 ? NULL : &resource),
327 (resourceCount == 0 ? 0 : 1),
328 preferredFlags,
329 requiredFlags,
330 &details);
331
332 if (beagleInstance < 0)
333 {
334 MrBayesPrint ("%s Failed to start BEAGLE instance\n", spacer);
335 }
336 else
337 {
338 if (division < 0)
339 {
340 /* do not use multipartition mode (division < 0) for CPU resource */
341 if (details.flags & BEAGLE_FLAG_FRAMEWORK_CPU)
342 {
343 MrBayesPrint ("\n%s Selected resource is the CPU, changing to multi-instance BEAGLE mode\n", spacer);
344 beagleFinalizeInstance(beagleInstance);
345 return -1;
346 }
347 else
348 {
349 MrBayesPrint ("\n%s Using BEAGLE v%s resource %i for %d division%s:\n", spacer, beagleGetVersion(), details.resourceNumber, numCurrentDivisions, (numCurrentDivisions > 1 ? "s" : ""));
350 }
351 }
352 else
353 {
354 MrBayesPrint ("\n%s Using BEAGLE v%s resource %i for division %d:\n", spacer, beagleGetVersion(), details.resourceNumber, division+1);
355 }
356 MrBayesPrint ("%s Rsrc Name : %s\n", spacer, details.resourceName);
357 MrBayesPrint ("%s Impl Name : %s\n", spacer, details.implName);
358 MrBayesPrint ("%s Flags:", spacer);
359 BeaglePrintFlags(details.flags);
360 MrBayesPrint ("%s MODEL STATES: %d", spacer, numModelStates);
361 MrBayesPrint ("\n");
362 beagleInstanceCount++;
363 }
364
365 #if defined (BEAGLE_V3_ENABLED)
366 /* use level-order traversal with CUDA implementation or OpenCL with multi-partition */
367 if(((details.flags & BEAGLE_FLAG_FRAMEWORK_CUDA) && division < 1 ) ||
368 ((details.flags & BEAGLE_FLAG_FRAMEWORK_OPENCL) && division < 0))
369 {
370 for (i=0; i<(numTrees * 2 * numGlobalChains); i++)
371 {
372 t = mcmcTree[i];
373 if (t != NULL)
374 {
375 if ((t->intDownPassLevel = (TreeNode **) SafeCalloc (numTaxa, sizeof (TreeNode *))) == NULL)
376 {
377 free (t->nodes);
378 free (t->allDownPass);
379 free (t);
380 return -1;
381 }
382 t->levelPassEnabled = 1;
383 GetDownPass(t);
384 }
385 }
386 }
387 /* set max number of CPU threads to be used by BEAGLE with CPU implementation */
388 if ((details.flags & BEAGLE_FLAG_FRAMEWORK_CPU) && beagleThreadCount > 1 && beagleThreadCount != 99)
389 {
390 beagleSetCPUThreadCount(beagleInstance, beagleThreadCount);
391 }
392 #endif /* BEAGLE_V3_ENABLED */
393
394 return beagleInstance;
395 }
396
397
398 /*-----------------------------------------------------------------
399 |
400 | LaunchBEAGLELogLikeForDivision: calculate the log likelihood
401 | of the new state of the chain for a single division
402 |
403 -----------------------------------------------------------------*/
LaunchBEAGLELogLikeForDivision(int chain,int d,ModelInfo * m,Tree * tree,MrBFlt * lnL)404 void LaunchBEAGLELogLikeForDivision(int chain, int d, ModelInfo* m, Tree* tree, MrBFlt* lnL)
405 {
406 int i, rescaleFreqNew;
407 int *isScalerNode;
408 TreeNode *p;
409
410 if (beagleScalingScheme == MB_BEAGLE_SCALE_ALWAYS)
411 {
412
413 # if defined (DEBUG_MB_BEAGLE_FLOW)
414 MrBayesPrint ("ALWAYS RESCALING\n");
415 # endif
416 #if defined (DEBUG_MB_BEAGLE_MULTIPART)
417 printf("m->upDateAll %d\n", m->upDateAll);
418 #endif
419 /* Flip and copy or reset site scalers */
420 FlipSiteScalerSpace(m, chain);
421 if (m->upDateAll == YES) {
422 for (i=0; i<m->nCijkParts; i++) {
423 beagleResetScaleFactors(m->beagleInstance, m->siteScalerIndex[chain] + i);
424 }
425 }
426 else
427 CopySiteScalers(m, chain);
428
429 TreeTiProbs_Beagle(tree, d, chain);
430 TreeCondLikes_Beagle_Always_Rescale(tree, d, chain);
431 TreeLikelihood_Beagle(tree, d, chain, lnL, (chainId[chain] % chainParams.numChains));
432 }
433 else
434 { /* MB_BEAGLE_SCALE_DYNAMIC */
435
436 /* This flag is only valid within this block */
437 m->rescaleBeagleAll = NO;
438 TreeTiProbs_Beagle(tree, d, chain);
439 if (m->successCount[chain] > 1000)
440 {
441 m->successCount[chain] = 10;
442 m->rescaleFreq[chain]++; /* increase rescaleFreq independent of whether we accept or reject new state */
443 m->rescaleFreqOld = rescaleFreqNew = m->rescaleFreq[chain];
444 for (i=0; i<tree->nIntNodes; i++)
445 {
446 p = tree->intDownPass[i];
447 if (p->upDateCl == YES) {
448 /* flip to the new workspace since TreeCondLikes_Beagle_Rescale_All() does not do it for
449 (p->upDateCl == YES) since it assumes that TreeCondLikes_Beagle_No_Rescale() did it */
450 FlipCondLikeSpace (m, chain, p->index);
451 }
452 }
453 goto rescale_all;
454 }
455
456 if (beagleScalingFrequency != 0 &&
457 m->beagleComputeCount[chain] % (beagleScalingFrequency) == 0)
458 {
459 m->rescaleFreqOld = rescaleFreqNew = m->rescaleFreq[chain];
460 for (i=0; i<tree->nIntNodes; i++)
461 {
462 p = tree->intDownPass[i];
463 if (p->upDateCl == YES) {
464 /* flip to the new workspace since TreeCondLikes_Beagle_Rescale_All() does not do it for
465 (p->upDateCl == YES) since it assumes that TreeCondLikes_Beagle_No_Rescale() did it */
466 FlipCondLikeSpace (m, chain, p->index);
467 }
468 }
469 goto rescale_all;
470 }
471
472 TreeCondLikes_Beagle_No_Rescale(tree, d, chain);
473
474 /* Check if likelihood is valid */
475 if (TreeLikelihood_Beagle(tree, d, chain, lnL, (chainId[chain] % chainParams.numChains)) == BEAGLE_ERROR_FLOATING_POINT)
476 {
477 m->rescaleFreqOld = rescaleFreqNew = m->rescaleFreq[chain];
478 if (rescaleFreqNew > 1 && m->successCount[chain] < 40)
479 {
480 if (m->successCount[chain] < 10)
481 {
482 if (m->successCount[chain] < 4)
483 {
484 rescaleFreqNew -= rescaleFreqNew >> 3; /* <== we cut up to 12,5% of rescaleFreq */
485 if (m->successCount[chain] < 2)
486 {
487 rescaleFreqNew -= rescaleFreqNew >> 3;
488 /* 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*/
489 /* 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 */
490 /* 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*/
491 m->rescaleFreqOld -= m->rescaleFreqOld >> 3;
492 }
493 m->rescaleFreqOld -= m->rescaleFreqOld >> 3;
494 m->rescaleFreqOld--;
495 m->rescaleFreqOld = (m->rescaleFreqOld ? m->rescaleFreqOld:1);
496 m->recalculateScalers = YES;
497 recalcScalers = YES;
498 }
499 }
500 rescaleFreqNew--;
501 rescaleFreqNew = (rescaleFreqNew ? rescaleFreqNew : 1);
502 }
503 m->successCount[chain] = 0;
504 rescale_all:
505 # if defined (DEBUG_MB_BEAGLE_FLOW)
506 MrBayesPrint ("NUMERICAL RESCALING\n");
507 # endif
508
509 m->rescaleBeagleAll = YES;
510 FlipSiteScalerSpace(m, chain);
511 isScalerNode = m->isScalerNode[chain];
512 while_loop:
513 ResetScalersPartition (isScalerNode, tree, rescaleFreqNew);
514 for (i=0; i<m->nCijkParts; i++)
515 {
516 beagleResetScaleFactors(m->beagleInstance, m->siteScalerIndex[chain] + i);
517 }
518
519 TreeCondLikes_Beagle_Rescale_All (tree, d, chain);
520 if (TreeLikelihood_Beagle(tree, d, chain, lnL, (chainId[chain] % chainParams.numChains)) == BEAGLE_ERROR_FLOATING_POINT)
521 {
522 if (rescaleFreqNew > 1)
523 {
524 /* Swap back scalers which were swapped in TreeCondLikes_Beagle_Rescale_All() */
525 for (i=0; i<tree->nIntNodes; i++)
526 {
527 p = tree->intDownPass[i];
528 if (isScalerNode[p->index] == YES)
529 FlipNodeScalerSpace (m, chain, p->index);
530 }
531 rescaleFreqNew -= rescaleFreqNew >> 3; /* <== we cut up to 12,5% of rescaleFreq */
532 rescaleFreqNew--; /* we cut extra 1 of rescaleFreq */
533 goto while_loop;
534 }
535 }
536 m->rescaleFreq[chain] = rescaleFreqNew;
537 }
538 }
539
540 /* Count number of evaluations */
541 m->beagleComputeCount[chain]++;
542 }
543
544
recalculateScalers(int chain)545 void recalculateScalers(int chain)
546 {
547 int i, d, rescaleFreqNew;
548 int *isScalerNode;
549 ModelInfo* m;
550 Tree *tree;
551 #if defined (BEAGLE_V3_ENABLED)
552 int divisionCount;
553 int *divisions;
554 #endif
555
556 if (modelSettings[0].useBeagleMultiPartitions == NO)
557 {
558 for (d=0; d<numCurrentDivisions; d++)
559 {
560 m = &modelSettings[d];
561 if (m->recalculateScalers == YES)
562 {
563 m->recalculateScalers = NO;
564 tree = GetTree(m->brlens, chain, state[chain]);
565
566 rescaleFreqNew = m->rescaleFreq[chain];
567 isScalerNode = m->isScalerNode[chain];
568
569 ResetScalersPartition (isScalerNode, tree, rescaleFreqNew);
570 for (i=0; i<m->nCijkParts; i++) {
571 beagleResetScaleFactors(m->beagleInstance, m->siteScalerIndex[chain] + i);
572 }
573 /* here it does not matter if we flip CL space or not */
574 TreeCondLikes_Beagle_Rescale_All (tree, d, chain);
575 }
576 }
577 }
578 #if defined (BEAGLE_V3_ENABLED)
579 else
580 {
581 divisions = (int *) SafeCalloc (numCurrentDivisions, sizeof(int));
582 divisionCount = 0;
583
584 for (d=0; d<numCurrentDivisions; d++)
585 {
586 m = &modelSettings[d];
587 if (m->recalculateScalers == YES)
588 {
589 m->recalculateScalers = NO;
590 tree = GetTree(m->brlens, chain, state[chain]);
591
592 rescaleFreqNew = m->rescaleFreq[chain];
593 isScalerNode = m->isScalerNode[chain];
594
595 ResetScalersPartition (isScalerNode, tree, rescaleFreqNew);
596 for (i=0; i<m->nCijkParts; i++) {
597 beagleResetScaleFactorsByPartition(m->beagleInstance, m->siteScalerIndex[chain] + i, m->divisionIndex);
598 }
599
600 divisions[divisionCount++] = d;
601 }
602 }
603
604 /* here it does not matter if we flip CL space or not */
605 TreeCondLikes_BeagleMultiPartition_Rescale_All(divisions, divisionCount, chain);
606
607 free(divisions);
608 }
609 #endif
610 }
611
612
BeagleAddGPUDevicesToList(int ** newResourceList,int * beagleResourceCount)613 void BeagleAddGPUDevicesToList(int **newResourceList, int *beagleResourceCount)
614 {
615 BeagleResourceList* beagleResources;
616 int i, gpuCount;
617
618 beagleResources = beagleGetResourceList();
619 if (*newResourceList == NULL) {
620 *newResourceList = (int*) SafeCalloc(sizeof(int), beagleResources->length);
621 }
622 gpuCount = 0;
623 for (i = 0; i < beagleResources->length; i++) {
624 if (beagleResources->list[i].supportFlags & BEAGLE_FLAG_PROCESSOR_GPU) {
625 (*newResourceList)[gpuCount] = i;
626 gpuCount++;
627 }
628 }
629 *beagleResourceCount = gpuCount;
630 }
631
632
BeagleRemoveGPUDevicesFromList(int ** beagleResource,int * beagleResourceCount)633 void BeagleRemoveGPUDevicesFromList(int **beagleResource, int *beagleResourceCount)
634 {
635 *beagleResourceCount = 0;
636 }
637
638
639 /*-----
640 |
641 | BeaglePrintResources: outputs the available BEAGLE resources
642 |
643 ----------*/
BeaglePrintResources()644 void BeaglePrintResources()
645 {
646 int i;
647 BeagleResourceList* beagleResources;
648
649 beagleResources = beagleGetResourceList();
650 MrBayesPrint ("%s Available resources reported by beagle library:\n", spacer);
651 for (i=0; i<beagleResources->length; i++)
652 {
653 MrBayesPrint ("\tResource %i:\n", i);
654 MrBayesPrint ("\tName: %s\n", beagleResources->list[i].name);
655 if (i > 0)
656 {
657 MrBayesPrint ("\tDesc: %s\n", beagleResources->list[i].description);
658 }
659 MrBayesPrint ("\tFlags:");
660 BeaglePrintFlags(beagleResources->list[i].supportFlags);
661 MrBayesPrint ("\n\n");
662 }
663 MrBayesPrint ("%s BEAGLE version: %s\n", spacer, beagleGetVersion());
664 }
665
666
BeagleCheckFlagCompatability(long inFlags)667 int BeagleCheckFlagCompatability(long inFlags)
668 {
669 if (inFlags & BEAGLE_FLAG_PROCESSOR_GPU) {
670 if (inFlags & BEAGLE_FLAG_VECTOR_SSE) {
671 MrBayesPrint ("%s Simultaneous use of GPU and SSE not available.\n", spacer);
672 return NO;
673 }
674 if (inFlags & BEAGLE_FLAG_THREADING_OPENMP) {
675 MrBayesPrint ("%s Simultaneous use of GPU and OpenMP not available.\n", spacer);
676 return NO;
677 }
678 }
679
680 return YES;
681 }
682
683
684 /*-------------------
685 |
686 | BeaglePrintFlags: outputs beagle instance details
687 |
688 ______________________*/
BeaglePrintFlags(long inFlags)689 void BeaglePrintFlags(long inFlags)
690 {
691 int i, k;
692 char *names[] = { "PROCESSOR_CPU",
693 "PROCESSOR_GPU",
694 "PROCESSOR_FPGA",
695 "PROCESSOR_CELL",
696 "PRECISION_DOUBLE",
697 "PRECISION_SINGLE",
698 "COMPUTATION_ASYNCH",
699 "COMPUTATION_SYNCH",
700 "EIGEN_REAL",
701 "EIGEN_COMPLEX",
702 "SCALING_MANUAL",
703 "SCALING_AUTO",
704 "SCALING_ALWAYS",
705 "SCALING_DYNAMIC",
706 "SCALERS_RAW",
707 "SCALERS_LOG",
708 "VECTOR_NONE",
709 "VECTOR_SSE",
710 "THREADING_NONE",
711 #if defined (BEAGLE_V3_ENABLED)
712 "THREADING_CPP",
713 #endif
714 "THREADING_OPENMP"
715 };
716 long flags[] = { BEAGLE_FLAG_PROCESSOR_CPU,
717 BEAGLE_FLAG_PROCESSOR_GPU,
718 BEAGLE_FLAG_PROCESSOR_FPGA,
719 BEAGLE_FLAG_PROCESSOR_CELL,
720 BEAGLE_FLAG_PRECISION_DOUBLE,
721 BEAGLE_FLAG_PRECISION_SINGLE,
722 BEAGLE_FLAG_COMPUTATION_ASYNCH,
723 BEAGLE_FLAG_COMPUTATION_SYNCH,
724 BEAGLE_FLAG_EIGEN_REAL,
725 BEAGLE_FLAG_EIGEN_COMPLEX,
726 BEAGLE_FLAG_SCALING_MANUAL,
727 BEAGLE_FLAG_SCALING_AUTO,
728 BEAGLE_FLAG_SCALING_ALWAYS,
729 BEAGLE_FLAG_SCALING_DYNAMIC,
730 BEAGLE_FLAG_SCALERS_RAW,
731 BEAGLE_FLAG_SCALERS_LOG,
732 BEAGLE_FLAG_VECTOR_NONE,
733 BEAGLE_FLAG_VECTOR_SSE,
734 BEAGLE_FLAG_THREADING_NONE,
735 #if defined (BEAGLE_V3_ENABLED)
736 BEAGLE_FLAG_THREADING_CPP,
737 #endif
738 BEAGLE_FLAG_THREADING_OPENMP
739 };
740
741 int flagCount = 20;
742
743 #if defined (BEAGLE_V3_ENABLED)
744 flagCount = 21;
745 #endif
746
747 for (i=k=0; i<flagCount; i++)
748 {
749 if (inFlags & flags[i])
750 {
751 if (k%4 == 0 && k > 0)
752 MrBayesPrint ("\n%s ", spacer);
753 MrBayesPrint (" %s", names[i]);
754 k++;
755 }
756 }
757 }
758
ScheduleLogLikeForAllDivisions()759 int ScheduleLogLikeForAllDivisions()
760 {
761 int d;
762 int divisionsToLaunch = 0;
763 ModelInfo *m;
764
765 if (numCurrentDivisions < 2) {
766 return 0;
767 }
768
769 for (d=0; d<numCurrentDivisions; d++) {
770 m = &modelSettings[d];
771 if (m->upDateCl == YES) {
772 divisionsToLaunch++;
773 }
774 }
775
776 return (divisionsToLaunch > 1);
777 }
778
779
780 /*----------------------------------------------------------------
781 |
782 | TreeCondLikes_Beagle: This routine updates all conditional
783 | (partial) likelihoods of a beagle instance while doing no rescaling.
784 | That potentialy can make final liklihood bad then calculation with rescaling needs to be done.
785 |
786 -----------------------------------------------------------------*/
TreeCondLikes_Beagle_No_Rescale(Tree * t,int division,int chain)787 int TreeCondLikes_Beagle_No_Rescale (Tree *t, int division, int chain)
788 {
789 int i, j, cumulativeScaleIndex, op;
790 TreeNode *p;
791 ModelInfo *m;
792 unsigned chil1Step, chil2Step;
793 int *isScalerNode;
794
795 m = &modelSettings[division];
796 isScalerNode = m->isScalerNode[chain];
797
798 op = 0;
799
800 for (i=0; i<t->nIntNodes; i++)
801 {
802
803 #if defined (BEAGLE_V3_ENABLED)
804 if (t->levelPassEnabled)
805 {
806 p = t->intDownPassLevel[i];
807 }
808 else
809 #endif
810 {
811 p = t->intDownPass[i];
812 }
813
814 /* check if conditional likelihoods need updating */
815 if (p->upDateCl == YES)
816 {
817 /* flip to the new workspace */
818 FlipCondLikeSpace (m, chain, p->index);
819
820 /* update conditional likelihoods */
821 m->operations[op].destinationPartials = m->condLikeIndex[chain][p->index ];
822 m->operations[op].child1Partials = m->condLikeIndex[chain][p->left->index ];
823 m->operations[op].child1TransitionMatrix = m->tiProbsIndex [chain][p->left->index ];
824 m->operations[op].child2Partials = m->condLikeIndex[chain][p->right->index];
825 m->operations[op].child2TransitionMatrix = m->tiProbsIndex [chain][p->right->index];
826
827 /* All partials for tips are the same across omega categories, thus we are doing the following two if statments.*/
828 if (p->left->left== NULL)
829 chil1Step=0;
830 else
831 chil1Step=1;
832
833 if (p->right->left== NULL)
834 chil2Step=0;
835 else
836 chil2Step=1;
837
838 m->operations[op].destinationScaleWrite = BEAGLE_OP_NONE;
839 cumulativeScaleIndex = BEAGLE_OP_NONE;
840 if (isScalerNode[p->index] == YES)
841 {
842 m->operations[op].destinationScaleRead = m->nodeScalerIndex[chain][p->index];
843 }
844 else
845 {
846 m->operations[op].destinationScaleRead = BEAGLE_OP_NONE;
847 }
848
849 for (j=1; j<m->nCijkParts; j++)
850 {
851 m->operations[op+1].destinationPartials = m->operations[op].destinationPartials + 1;
852 m->operations[op+1].child1Partials = m->operations[op].child1Partials + chil1Step;
853 m->operations[op+1].child1TransitionMatrix = m->operations[op].child1TransitionMatrix + 1;
854 m->operations[op+1].child2Partials = m->operations[op].child2Partials + chil2Step;
855 m->operations[op+1].child2TransitionMatrix = m->operations[op].child2TransitionMatrix + 1;
856
857 m->operations[op+1].destinationScaleWrite = BEAGLE_OP_NONE;
858 if ( isScalerNode[p->index] == YES )
859 {
860 m->operations[op+1].destinationScaleRead = m->operations[op].destinationScaleRead + 1;
861 }
862 else
863 {
864 m->operations[op+1].destinationScaleRead = BEAGLE_OP_NONE;
865 }
866 op++;
867 }
868 op++;
869 }
870 }
871
872 beagleUpdatePartials(m->beagleInstance,
873 m->operations,
874 op,
875 cumulativeScaleIndex);
876
877 return NO_ERROR;
878 }
879
880
881 /*----------------------------------------------------------------
882 |
883 | TreeCondLikes_Beagle: This routine updates all conditional
884 | (partial) likelihoods of a beagle instance while rescaling at every node.
885 | Note: all nodes get recalculated.
886 |
887 -----------------------------------------------------------------*/
TreeCondLikes_Beagle_Rescale_All(Tree * t,int division,int chain)888 int TreeCondLikes_Beagle_Rescale_All (Tree *t, int division, int chain)
889 {
890 int i, j, cumulativeScaleIndex, op, opJ, opPrev;
891 TreeNode *p;
892 ModelInfo *m;
893 unsigned chil1Step, chil2Step;
894 int *isScalerNode;
895
896 m = &modelSettings[division];
897 isScalerNode = m->isScalerNode[chain];
898
899 op = 0;
900
901 for (i=0; i<t->nIntNodes; i++)
902 {
903
904 #if defined (BEAGLE_V3_ENABLED)
905 if (t->levelPassEnabled)
906 {
907 p = t->intDownPassLevel[i];
908 }
909 else
910 #endif
911 {
912 p = t->intDownPass[i];
913 }
914
915 if (p->upDateCl == NO)
916 {
917 //p->upDateCl = YES;
918 /* flip to the new workspace */
919 FlipCondLikeSpace (m, chain, p->index);
920 }
921
922 /* update conditional likelihoods */
923 m->operations[op].destinationPartials = m->condLikeIndex[chain][p->index ];
924 m->operations[op].child1Partials = m->condLikeIndex[chain][p->left->index ];
925 m->operations[op].child1TransitionMatrix = m->tiProbsIndex [chain][p->left->index ];
926 m->operations[op].child2Partials = m->condLikeIndex[chain][p->right->index];
927 m->operations[op].child2TransitionMatrix = m->tiProbsIndex [chain][p->right->index];
928
929 /* All partials for tips are the same across omega catigoris, thus we are doing the following two if statments.*/
930 if (p->left->left== NULL)
931 chil1Step=0;
932 else
933 chil1Step=1;
934
935 if (p->right->left== NULL)
936 chil2Step=0;
937 else
938 chil2Step=1;
939
940 m->operations[op].destinationScaleRead = BEAGLE_OP_NONE;
941 if (isScalerNode[p->index] == YES)
942 {
943 FlipNodeScalerSpace (m, chain, p->index);
944 m->operations[op].destinationScaleWrite = m->nodeScalerIndex[chain][p->index];
945 cumulativeScaleIndex = m->siteScalerIndex[chain];
946 }
947 else
948 {
949 m->operations[op].destinationScaleWrite = BEAGLE_OP_NONE;
950 cumulativeScaleIndex = BEAGLE_OP_NONE;
951 }
952
953
954
955 for (j=1; j<m->nCijkParts; j++)
956 {
957 opJ = op + j * t->nIntNodes;
958 opPrev = op + (j-1) * t->nIntNodes;
959
960 m->operations[opJ].destinationPartials = m->operations[opPrev].destinationPartials + 1;
961 m->operations[opJ].child1Partials = m->operations[opPrev].child1Partials + chil1Step;
962 m->operations[opJ].child1TransitionMatrix = m->operations[opPrev].child1TransitionMatrix + 1;
963 m->operations[opJ].child2Partials = m->operations[opPrev].child2Partials + chil2Step;
964 m->operations[opJ].child2TransitionMatrix = m->operations[opPrev].child2TransitionMatrix + 1;
965
966 m->operations[opJ].destinationScaleRead = BEAGLE_OP_NONE;
967 if ( isScalerNode[p->index] == YES )
968 {
969 m->operations[opJ].destinationScaleWrite = m->operations[opPrev].destinationScaleWrite + 1;
970 }
971 else
972 {
973 m->operations[opJ].destinationScaleWrite = BEAGLE_OP_NONE;
974 }
975 }
976 op++;
977 }
978
979 cumulativeScaleIndex = m->siteScalerIndex[chain];
980 for (j=0; j<m->nCijkParts; j++)
981 {
982 beagleUpdatePartials(m->beagleInstance,
983 &m->operations[j*t->nIntNodes],
984 t->nIntNodes,
985 cumulativeScaleIndex);
986 cumulativeScaleIndex++;
987 }
988
989 return NO_ERROR;
990 }
991
992
993 /*----------------------------------------------------------------
994 |
995 | TreeCondLikes_Beagle: This routine updates and rescales
996 | all conditional (partial) likelihoods of a beagle instance.
997 |
998 -----------------------------------------------------------------*/
TreeCondLikes_Beagle_Always_Rescale(Tree * t,int division,int chain)999 int TreeCondLikes_Beagle_Always_Rescale (Tree *t, int division, int chain)
1000 {
1001 int i, j, op, opJ, opPrev, scaleOp;
1002 int destinationScaleRead, cumulativeScaleIndex;
1003 TreeNode *p;
1004 ModelInfo *m;
1005 unsigned chil1Step, chil2Step;
1006
1007 m = &modelSettings[division];
1008
1009 op = 0;
1010 scaleOp = 0;
1011
1012 for (i=0; i<t->nIntNodes; i++)
1013 {
1014
1015 #if defined (BEAGLE_V3_ENABLED)
1016 if (t->levelPassEnabled)
1017 {
1018 p = t->intDownPassLevel[i];
1019 }
1020 else
1021 #endif
1022 {
1023 p = t->intDownPass[i];
1024 }
1025
1026 /* check if conditional likelihoods need updating */
1027 if (p->upDateCl == YES)
1028 {
1029 /* remove old scalers */
1030 if (m->upDateAll == NO)
1031 {
1032 destinationScaleRead = m->nodeScalerIndex[chain][p->index];
1033 for (j=0; j<m->nCijkParts; j++)
1034 {
1035 m->scaleFactorsOps[scaleOp+j*t->nIntNodes] = destinationScaleRead;
1036 destinationScaleRead++;
1037 }
1038 scaleOp++;
1039 }
1040
1041 /* flip to the new workspace */
1042 FlipCondLikeSpace (m, chain, p->index);
1043 FlipNodeScalerSpace (m, chain, p->index);
1044
1045 /* update conditional likelihoods */
1046 m->operations[op].destinationPartials = m->condLikeIndex[chain][p->index ];
1047 m->operations[op].child1Partials = m->condLikeIndex[chain][p->left->index ];
1048 m->operations[op].child1TransitionMatrix = m->tiProbsIndex [chain][p->left->index ];
1049 m->operations[op].child2Partials = m->condLikeIndex[chain][p->right->index];
1050 m->operations[op].child2TransitionMatrix = m->tiProbsIndex [chain][p->right->index];
1051
1052 /* All partials for tips are the same across omega catigoris, thus we are doing the following two if statments.*/
1053 if (p->left->left== NULL && p->left->right== NULL)
1054 chil1Step=0;
1055 else
1056 chil1Step=1;
1057
1058 if (p->right->left== NULL && p->right->right== NULL)
1059 chil2Step=0;
1060 else
1061 chil2Step=1;
1062
1063 m->operations[op].destinationScaleWrite = m->nodeScalerIndex[chain][p->index];
1064 m->operations[op].destinationScaleRead = BEAGLE_OP_NONE;
1065
1066 for (j=1; j<m->nCijkParts; j++)
1067 {
1068 opJ = op + j * t->nIntNodes;
1069 opPrev = op + (j-1) * t->nIntNodes;
1070
1071 m->operations[opJ].destinationPartials = m->operations[opPrev].destinationPartials + 1;
1072 m->operations[opJ].child1Partials = m->operations[opPrev].child1Partials + chil1Step;
1073 m->operations[opJ].child1TransitionMatrix = m->operations[opPrev].child1TransitionMatrix + 1;
1074 m->operations[opJ].child2Partials = m->operations[opPrev].child2Partials + chil2Step;
1075 m->operations[opJ].child2TransitionMatrix = m->operations[opPrev].child2TransitionMatrix + 1;
1076
1077 m->operations[opJ].destinationScaleWrite = m->operations[opPrev].destinationScaleWrite + 1;
1078 m->operations[opJ].destinationScaleRead = BEAGLE_OP_NONE;
1079 }
1080 #if defined (DEBUG_MB_BEAGLE_MULTIPART)
1081 printf("%d %d ", op, m->operations[op].destinationPartials);
1082 printf("%d ", m->operations[op].destinationScaleWrite);
1083 printf("%d ", m->operations[op].destinationScaleRead);
1084 printf("%d ", m->operations[op].child1Partials);
1085 printf("%d ", m->operations[op].child1TransitionMatrix);
1086 printf("%d ", m->operations[op].child2Partials);
1087 printf("%d\n", m->operations[op].child2TransitionMatrix);
1088 #endif
1089 op++;
1090 }
1091 } /* end of for */
1092
1093 cumulativeScaleIndex = m->siteScalerIndex[chain];
1094 for (j=0; j<m->nCijkParts; j++)
1095 {
1096 if (m->upDateAll == NO)
1097 {
1098 beagleRemoveScaleFactors(m->beagleInstance,
1099 &m->scaleFactorsOps[j*t->nIntNodes],
1100 scaleOp,
1101 cumulativeScaleIndex);
1102 }
1103
1104 beagleUpdatePartials(m->beagleInstance,
1105 &m->operations[j*t->nIntNodes],
1106 op,
1107 cumulativeScaleIndex);
1108 cumulativeScaleIndex++;
1109 }
1110
1111 return NO_ERROR;
1112 }
1113
1114
1115 /**---------------------------------------------------------------------------
1116 |
1117 | TreeLikelihood_Beagle: Accumulate the log likelihoods calculated by Beagle
1118 | at the root.
1119 |
1120 ---------------------------------------- -------------------------------------*/
TreeLikelihood_Beagle(Tree * t,int division,int chain,MrBFlt * lnL,int whichSitePats)1121 int TreeLikelihood_Beagle (Tree *t, int division, int chain, MrBFlt *lnL, int whichSitePats)
1122 {
1123 int i, j, c = 0, nStates, hasPInvar, beagleReturn;
1124 MrBFlt *swr, s01, s10, probOn, probOff, covBF[40], pInvar=0.0, *bs, freq, likeI, lnLikeI, diff, *omegaCatFreq;
1125 CLFlt *clInvar=NULL, *nSitesOfPat;
1126 double *nSitesOfPat_Beagle;
1127 TreeNode *p;
1128 ModelInfo *m;
1129 double pUnobserved;
1130
1131 # if defined (MB_PRINT_DYNAMIC_RESCALE_FAIL_STAT)
1132 static unsigned countBeagleDynamicFail=0;
1133 static unsigned countALL=0;
1134 # endif
1135
1136 /* find root node */
1137 p = t->root->left;
1138
1139 /* find model settings and nStates, pInvar, invar cond likes */
1140 m = &modelSettings[division];
1141
1142 nStates = m->numModelStates;
1143 if (m->pInvar == NULL)
1144 {
1145 hasPInvar = NO;
1146 }
1147 else
1148 {
1149 hasPInvar = YES;
1150 pInvar = *(GetParamVals (m->pInvar, chain, state[chain]));
1151 clInvar = m->invCondLikes;
1152 }
1153
1154 /* find base frequencies */
1155 bs = GetParamSubVals (m->stateFreq, chain, state[chain]);
1156
1157 /* if covarion model, adjust base frequencies */
1158 if (m->switchRates != NULL)
1159 {
1160 /* find the stationary frequencies */
1161 swr = GetParamVals(m->switchRates, chain, state[chain]);
1162 s01 = swr[0];
1163 s10 = swr[1];
1164 probOn = s01 / (s01 + s10);
1165 probOff = 1.0 - probOn;
1166
1167 /* now adjust the base frequencies; on-state stored first in cond likes */
1168 for (j=0; j<nStates/2; j++)
1169 {
1170 covBF[j] = bs[j] * probOn;
1171 covBF[j+nStates/2] = bs[j] * probOff;
1172 }
1173
1174 /* finally set bs pointer to adjusted values */
1175 bs = covBF;
1176 }
1177
1178 /* TODO Really only need to check if state frequencies have changed */
1179 if (m->upDateCijk == YES)
1180 {
1181 /* set base frequencies in BEAGLE instance */
1182 for (i=0; i<m->nCijkParts; i++)
1183 beagleSetStateFrequencies(m->beagleInstance,
1184 m->cijkIndex[chain] + i,
1185 bs);
1186 }
1187
1188 /* find category frequencies */
1189 if (hasPInvar == NO)
1190 freq = 1.0 / m->numRateCats;
1191 else
1192 freq = (1.0 - pInvar) / m->numRateCats;
1193 #if defined (DEBUG_MB_BEAGLE_MULTIPART)
1194 printf("freq = %f\n", freq);
1195 #endif
1196 /* TODO: cat weights only need to be set when they change */
1197 /* set category frequencies in beagle instance */
1198 if (m->numOmegaCats > 1)
1199 {
1200 omegaCatFreq = GetParamSubVals(m->omega, chain, state[chain]);
1201 for (i=0; i<m->nCijkParts; i++)
1202 {
1203 for (j=0; j<m->numRateCats; j++)
1204 m->inWeights[j] = freq * omegaCatFreq[i];
1205 beagleSetCategoryWeights(m->beagleInstance,
1206 m->cijkIndex[chain] + i,
1207 m->inWeights);
1208 }
1209 }
1210 else if (hasPInvar == YES)
1211 {
1212 for (i=0; i<m->numRateCats; i++)
1213 m->inWeights[i] = freq;
1214 beagleSetCategoryWeights(m->beagleInstance,
1215 m->cijkIndex[chain],
1216 m->inWeights);
1217 }
1218
1219 /* find nSitesOfPat */
1220 nSitesOfPat = numSitesOfPat + (whichSitePats*numCompressedChars) + m->compCharStart;
1221
1222 /* TODO: pattern weights only need to be set when they change */
1223 /* set pattern weights in beagle instance if using dynamic reweighting */
1224 if (chainParams.weightScheme[0] + chainParams.weightScheme[1] > ETA)
1225 {
1226 nSitesOfPat_Beagle = (double *) SafeMalloc (m->numChars * sizeof(double));
1227 for (c=0; c<m->numChars; c++)
1228 nSitesOfPat_Beagle[c] = numSitesOfPat[m->compCharStart + c];
1229 beagleSetPatternWeights(m->beagleInstance,
1230 nSitesOfPat_Beagle);
1231 SAFEFREE (nSitesOfPat_Beagle);
1232 }
1233
1234 /* find root log likelihoods and scalers */
1235 for (i=0; i<m->nCijkParts; i++)
1236 {
1237 m->bufferIndices[i] = m->condLikeIndex[chain][p->index] + i;
1238 m->eigenIndices[i] = m->cijkIndex[chain] + i;
1239 m->cumulativeScaleIndices[i] = m->siteScalerIndex[chain] + i;
1240 if (t->isRooted == NO)
1241 {
1242 m->childBufferIndices[i] = m->condLikeIndex [chain][p->anc->index];
1243 m->childTiProbIndices[i] = m->tiProbsIndex [chain][p->index] + i;
1244 }
1245 }
1246 #if defined (DEBUG_MB_BEAGLE_MULTIPART)
1247 printf("bufferIndices = %d, eigenIndices = %d, cumulativeScaleIndices = %d\n", m->bufferIndices[0], m->eigenIndices[0], m->cumulativeScaleIndices[0]);
1248 #endif
1249 /* reset lnL */
1250 *lnL = 0.0;
1251
1252 /* get root log likelihoods */
1253 if (t->isRooted == YES)
1254 {
1255 beagleReturn = beagleCalculateRootLogLikelihoods(m->beagleInstance,
1256 m->bufferIndices,
1257 m->eigenIndices,
1258 m->eigenIndices,
1259 m->cumulativeScaleIndices,
1260 m->nCijkParts,
1261 lnL);
1262
1263 }
1264 else
1265 {
1266 beagleReturn = beagleCalculateEdgeLogLikelihoods(m->beagleInstance,
1267 m->bufferIndices,
1268 m->childBufferIndices,
1269 m->childTiProbIndices,
1270 NULL,
1271 NULL,
1272 m->eigenIndices,
1273 m->eigenIndices,
1274 m->cumulativeScaleIndices,
1275 m->nCijkParts,
1276 lnL,
1277 NULL,
1278 NULL);
1279 }
1280
1281 # if defined (MB_PRINT_DYNAMIC_RESCALE_FAIL_STAT)
1282 countALL++;
1283 # endif
1284 if (*lnL > DBL_MAX || *lnL < -DBL_MAX) {
1285 beagleReturn = BEAGLE_ERROR_FLOATING_POINT;
1286 }
1287 if (beagleReturn == BEAGLE_ERROR_FLOATING_POINT)
1288 {
1289 # if defined (MB_PRINT_DYNAMIC_RESCALE_FAIL_STAT)
1290 countBeagleDynamicFail++;
1291 MrBayesPrint ("DEBUG INFO (not an error) countBeagleDynamicFail:%d countALL:%d\n", countBeagleDynamicFail, countALL);
1292 # endif
1293 return beagleReturn;
1294 }
1295 assert (beagleReturn == BEAGLE_SUCCESS);
1296 m->successCount[chain]++;
1297
1298 #if defined (DEBUG_MB_BEAGLE_MULTIPART_SITELNL)
1299 beagleGetSiteLogLikelihoods(m->beagleInstance, m->logLikelihoods);
1300 printf("lnL = ");
1301 for (c=0; c<m->numChars; c++)
1302 {
1303 printf("[%d] %f ", c, m->logLikelihoods[c]);
1304 }
1305 printf("\n");
1306 #endif
1307
1308 /* accumulate logs across sites */
1309 if (hasPInvar == NO)
1310 {
1311 if (m->dataType == RESTRICTION)
1312 {
1313 beagleGetSiteLogLikelihoods(m->beagleInstance, m->logLikelihoods);
1314 (*lnL) = 0.0;
1315 pUnobserved = 0.0;
1316 for (c=0; c<m->numDummyChars; c++)
1317 {
1318 pUnobserved += exp((double)m->logLikelihoods[c]);
1319 }
1320 /* correct for absent characters */
1321 (*lnL) -= log (1-pUnobserved) * (m->numUncompressedChars);
1322 for (; c<m->numChars; c++)
1323 {
1324 (*lnL) += m->logLikelihoods[c] * nSitesOfPat[c];
1325 }
1326 }
1327 /* already done, just check for numerical errors */
1328 assert ((*lnL) == (*lnL));
1329 }
1330 else
1331 {
1332 /* has invariable category */
1333 beagleGetSiteLogLikelihoods(m->beagleInstance,
1334 m->logLikelihoods);
1335 (*lnL) = 0.0;
1336 for (c=0; c<m->numChars; c++)
1337 {
1338 likeI = 0.0;
1339 for (j=0; j<nStates; j++)
1340 likeI += (*(clInvar++)) * bs[j];
1341 if (likeI != 0.0)
1342 {
1343 lnLikeI = log(likeI * pInvar);
1344 diff = lnLikeI - m->logLikelihoods[c];
1345 }
1346 else
1347 diff = -1000.0;
1348 if (diff < -200.0)
1349 (*lnL) += m->logLikelihoods[c] * nSitesOfPat[c];
1350 else if (diff > 200.0)
1351 (*lnL) += lnLikeI * nSitesOfPat[c];
1352 else
1353 {
1354 (*lnL) += (m->logLikelihoods[c] + log(1.0 + exp(diff))) * nSitesOfPat[c];
1355 }
1356
1357 /* check for numerical errors */
1358 assert ((*lnL) == (*lnL));
1359 }
1360 }
1361
1362 return (NO_ERROR);
1363 }
1364
1365
1366 /*----------------------------------------------------------------
1367 |
1368 | TreeTiProbs_Beagle: This routine updates all transition
1369 | probability matrices of a beagle instance.
1370 |
1371 -----------------------------------------------------------------*/
TreeTiProbs_Beagle(Tree * t,int division,int chain)1372 int TreeTiProbs_Beagle (Tree *t, int division, int chain)
1373 {
1374 int i, j, k, count;
1375 MrBFlt correctionFactor, theRate, baseRate, *catRate, length;
1376 TreeNode *p;
1377 ModelInfo *m;
1378
1379 /* get model settings */
1380 m = &modelSettings[division];
1381
1382 /* find the correction factor to make branch lengths
1383 in terms of expected number of substitutions per character */
1384 correctionFactor = 1.0;
1385 if (m->dataType == DNA || m->dataType == RNA)
1386 {
1387 if (m->nucModelId == NUCMODEL_DOUBLET)
1388 correctionFactor = 2.0;
1389 else if (m->nucModelId == NUCMODEL_CODON)
1390 correctionFactor = 3.0;
1391 }
1392
1393 /* get rate multipliers (for gamma & partition specific rates) */
1394 theRate = 1.0;
1395 baseRate = GetRate (division, chain);
1396
1397 /* compensate for invariable sites if appropriate */
1398 if (m->pInvar != NULL)
1399 baseRate /= (1.0 - (*GetParamVals(m->pInvar, chain, state[chain])));
1400
1401 /* get category rates for gamma */
1402 if (m->shape == NULL)
1403 catRate = &theRate;
1404 else
1405 catRate = GetParamSubVals (m->shape, chain, state[chain]);
1406
1407 /* get effective category rates */
1408 for (k=0; k<m->numRateCats; k++)
1409 m->inRates[k] = baseRate * catRate[k] * correctionFactor;
1410
1411 /* TODO: only need to set category rates when they change */
1412 /* set category rates */
1413 beagleSetCategoryRates(m->beagleInstance, m->inRates);
1414
1415 /* get ti prob indices and branch lengths to update */
1416 for (i=count=0; i<t->nNodes; i++)
1417 {
1418 p = t->allDownPass[i];
1419
1420 /* check if transition probs need updating */
1421 if (p->upDateTi == YES)
1422 {
1423 /* flip transition probability */
1424 FlipTiProbsSpace (m, chain, p->index);
1425
1426 /* find length */
1427 if (m->cppEvents != NULL)
1428 {
1429 length = GetParamSubVals (m->cppEvents, chain, state[chain])[p->index];
1430 }
1431 else if (m->tk02BranchRates != NULL)
1432 {
1433 length = GetParamSubVals (m->tk02BranchRates, chain, state[chain])[p->index];
1434 }
1435 else if (m->igrBranchRates != NULL)
1436 {
1437 length = GetParamSubVals (m->igrBranchRates, chain, state[chain])[p->index];
1438 }
1439 else if (m->mixedBrchRates != NULL)
1440 {
1441 length = GetParamSubVals (m->mixedBrchRates, chain, state[chain])[p->index];
1442 }
1443 else
1444 length = p->length;
1445
1446 /* numerical errors might ensue if we allow very large or very small branch lengths, which might
1447 occur in relaxed clock models; an elegant solution would be to substitute the stationary
1448 probs and initial probs but for now we truncate lengths at small or large values */
1449 if (length > BRLENS_MAX)
1450 length = BRLENS_MAX;
1451 else if (length < BRLENS_MIN)
1452 length = BRLENS_MIN;
1453
1454 m->branchLengths[count] = length;
1455
1456 /* find index */
1457 m->tiProbIndices[count] = m->tiProbsIndex[chain][p->index];
1458 #if defined (DEBUG_MB_BEAGLE_MULTIPART)
1459 printf("%d %f ", count, m->branchLengths[count]);
1460 printf(" %d", m->tiProbIndices[count]);
1461 printf(" %d", 0);
1462 printf(" %d\n", m->cijkIndex[chain]);
1463 #endif
1464 count++;
1465 }
1466 }
1467
1468 /* TODO: only need to update branches that have changed */
1469 /* calculate transition probabilities */
1470 if (count > 0) {
1471 for (i=0; i<m->nCijkParts; i++)
1472 {
1473 beagleUpdateTransitionMatrices(m->beagleInstance,
1474 m->cijkIndex[chain] + i,
1475 m->tiProbIndices,
1476 NULL,
1477 NULL,
1478 m->branchLengths,
1479 count);
1480 for (j=0; j<count; j++)
1481 m->tiProbIndices[j]++;
1482 }
1483 }
1484
1485 /* return success */
1486 return NO_ERROR;
1487 }
1488
1489 #if defined (BEAGLE_V3_ENABLED)
1490
1491 /*------------------------------------------------------------------------
1492 |
1493 | InitBeagleMultiPartitionInstance: create and initialize a beagle instance for multiple partitions
1494 |
1495 -------------------------------------------------------------------------*/
InitBeagleMultiPartitionInstance()1496 int InitBeagleMultiPartitionInstance ()
1497 {
1498 int i, j, k, c, s, d, *inStates, numPartAmbigTips, sizePD;
1499 int nCijkParts, numRateCats, numModelStates, numCondLikes, numScalers;
1500 int numChars, numTiProbs;
1501 int *anyDivPartAmbigTip;
1502 double *inPartials;
1503 ModelInfo *m, *mPrev;
1504 BitsLong *charBits;
1505 int beagleInstance;
1506
1507 m = &modelSettings[0];
1508
1509 if (m->beagleInstance == -99) {
1510 MrBayesPrint ("\n%s Error: Consecutive MCMC runs not currently supported with BEAGLE. Please restart MrBayes.\n\n", spacer);
1511 exit(1);
1512 }
1513
1514 if (m->useBeagle == NO)
1515 return ERROR;
1516
1517 /* check if divisions can be merged into single instance */
1518 for (d=1; d<numCurrentDivisions; d++)
1519 {
1520 mPrev = m;
1521 m = &modelSettings[d];
1522 if (m->useBeagle == NO ||
1523 m->nCijkParts != mPrev->nCijkParts ||
1524 m->numRateCats != mPrev->numRateCats ||
1525 m->numModelStates != mPrev->numModelStates ||
1526 m->numCondLikes != mPrev->numCondLikes ||
1527 m->numScalers != mPrev->numScalers ||
1528 m->numTiProbs != mPrev->numTiProbs)
1529 {
1530 #ifdef DEBUG_MB_BEAGLE_MULTIPART
1531 printf("m->useBeagle = %d\n m->nCijkParts = %d | mPrev->nCijkParts = %d\n m->numRateCats = %d | mPrev->numRateCats = %d\n m->numModelStates = %d | mPrev->numModelStates = %d\n m->numCondLikes = %d | mPrev->numCondLikes = %d\n m->numScalers = %d | mPrev->numScalers = %d\n m->numTiProbs = %d | mPrev->numTiProbs = %d\n",
1532 m->useBeagle,
1533 m->nCijkParts , mPrev->nCijkParts ,
1534 m->numRateCats , mPrev->numRateCats ,
1535 m->numModelStates, mPrev->numModelStates,
1536 m->numCondLikes , mPrev->numCondLikes ,
1537 m->numScalers , mPrev->numScalers ,
1538 m->numTiProbs , mPrev->numTiProbs);
1539 #endif
1540 return ERROR;
1541 }
1542 }
1543
1544 anyDivPartAmbigTip = (int *) SafeCalloc (numLocalTaxa * numCurrentDivisions, sizeof(int));
1545
1546 m = &modelSettings[0];
1547
1548 nCijkParts = (m->nCijkParts == 0 ? 1 : m->nCijkParts);
1549 numRateCats = m->numRateCats;
1550 numModelStates = m->numModelStates;
1551 numCondLikes = m->numCondLikes;
1552 numScalers = m->numScalers;
1553 numChars = 0;
1554 numTiProbs = 0;
1555
1556 sizePD = nCijkParts*numCurrentDivisions;
1557 m->cijkIndicesAll = (int *) SafeCalloc (sizePD*2*numLocalTaxa, sizeof(int));
1558 m->categoryRateIndicesAll = (int *) SafeCalloc (sizePD*2*numLocalTaxa, sizeof(int));
1559 m->operationsAll = (BeagleOperationByPartition *) SafeCalloc(sizePD*m->numCondLikes, sizeof(BeagleOperationByPartition));
1560 m->tiProbIndices = (int *) SafeCalloc (sizePD*2*numLocalTaxa, sizeof(int));
1561 m->branchLengths = (MrBFlt *) SafeCalloc (sizePD*2*numLocalTaxa, sizeof(MrBFlt));
1562 m->bufferIndices = (int *) SafeCalloc (sizePD, sizeof(int));
1563 m->childBufferIndices = (int *) SafeCalloc (sizePD, sizeof(int));
1564 m->childTiProbIndices = (int *) SafeCalloc (sizePD, sizeof(int));
1565 m->eigenIndices = (int *) SafeCalloc (sizePD, sizeof(int));
1566 m->cumulativeScaleIndices = (int *) SafeCalloc (sizePD, sizeof(int));
1567
1568 for (d=0; d<numCurrentDivisions; d++)
1569 {
1570 m = &modelSettings[d];
1571
1572 /* at least one eigen buffer needed */
1573 if (m->nCijkParts == 0)
1574 m->nCijkParts = 1;
1575
1576 /* allocate memory used by beagle */
1577 m->inRates = (MrBFlt *) SafeCalloc (m->numRateCats, sizeof(MrBFlt));
1578 m->inWeights = (MrBFlt *) SafeCalloc (m->numRateCats*m->nCijkParts, sizeof(MrBFlt));
1579 m->operationsByPartition = (BeagleOperationByPartition *) SafeCalloc(m->numCondLikes*m->nCijkParts, sizeof(BeagleOperationByPartition));
1580 m->scaleFactorsOps = (int *) SafeCalloc (2*numLocalTaxa*m->nCijkParts, sizeof(int));
1581
1582 numChars += m->numChars;
1583 numTiProbs += m->numTiProbs;
1584 }
1585
1586 m = &modelSettings[0];
1587 m->numCharsAll = numChars;
1588 m->logLikelihoodsAll = (MrBFlt *) SafeCalloc ((numLocalChains)*m->numCharsAll, sizeof(MrBFlt));
1589
1590 numPartAmbigTips = 0;
1591 if (m->numStates != m->numModelStates)
1592 numPartAmbigTips = numLocalTaxa;
1593 else
1594 {
1595 for (i=0; i<numLocalTaxa; i++)
1596 {
1597 #if defined(BEAGLE_V3_ENABLED)
1598 if (!beagleAllFloatTips)
1599 #else
1600 if (1)
1601 #endif
1602 {
1603 for (d=0; d<numCurrentDivisions; d++)
1604 {
1605 m = &modelSettings[d];
1606 if (m->isPartAmbig[i] == YES)
1607 {
1608 anyDivPartAmbigTip[i] = YES;
1609 numPartAmbigTips++;
1610 break;
1611 }
1612 }
1613 }
1614 else
1615 {
1616 anyDivPartAmbigTip[i] = YES;
1617 numPartAmbigTips++;
1618 }
1619 }
1620 }
1621
1622 beagleInstance = createBeagleInstance(m, nCijkParts, numRateCats, numModelStates, numCondLikes, numScalers, numChars, numTiProbs, numPartAmbigTips, -1);
1623
1624 if (beagleInstance < 0)
1625 return ERROR;
1626
1627 for (d=0; d<numCurrentDivisions; d++)
1628 {
1629 m = &modelSettings[d];
1630 m->beagleInstance = beagleInstance;
1631 m->useBeagleMultiPartitions = YES;
1632 }
1633
1634 /* initialize tip data */
1635 inStates = (int *) SafeMalloc (numChars * sizeof(int));
1636 if (!inStates)
1637 return ERROR;
1638 inPartials = (double *) SafeMalloc (numChars * numModelStates * sizeof(double));
1639 if (!inPartials)
1640 return ERROR;
1641 for (i=0; i<numLocalTaxa; i++)
1642 {
1643 if (anyDivPartAmbigTip[i] == NO)
1644 {
1645 k=0;
1646 for (d=0; d<numCurrentDivisions; d++)
1647 {
1648 m = &modelSettings[d];
1649 charBits = m->parsSets[i];
1650 for (c=0; c<m->numChars; c++)
1651 {
1652 for (s=j=0; s<m->numModelStates; s++)
1653 {
1654 if (IsBitSet(s, charBits))
1655 {
1656 inStates[k] = s;
1657 j++;
1658 }
1659 }
1660 if (j == m->numModelStates)
1661 inStates[k] = j;
1662 else
1663 assert (j==1);
1664 charBits += m->nParsIntsPerSite;
1665 k++;
1666 }
1667 }
1668 beagleSetTipStates(beagleInstance, i, inStates);
1669 }
1670 else /* if (m->isPartAmbig == YES) */
1671 {
1672 k = 0;
1673 for (d=0; d<numCurrentDivisions; d++)
1674 {
1675 m = &modelSettings[d];
1676 charBits = m->parsSets[i];
1677 for (c=0; c<m->numChars; c++)
1678 {
1679 for (s=0; s<m->numModelStates; s++)
1680 {
1681 if (IsBitSet(s%m->numStates, charBits))
1682 inPartials[k++] = 1.0;
1683 else
1684 inPartials[k++] = 0.0;
1685 }
1686 charBits += m->nParsIntsPerSite;
1687 }
1688 }
1689 beagleSetTipPartials(beagleInstance, i, inPartials);
1690 }
1691 }
1692
1693 free (anyDivPartAmbigTip);
1694 free (inStates);
1695 free (inPartials);
1696
1697 return NO_ERROR;
1698 }
1699
1700
LaunchBEAGLELogLikeMultiPartition(int * divisions,int divisionCount,int chain,MrBFlt * lnL)1701 void LaunchBEAGLELogLikeMultiPartition(int* divisions, int divisionCount, int chain, MrBFlt* lnL)
1702 {
1703 int i, d, dIndex, rescaleDivisionCount, norescaleDivisionCount, beagleReturn;
1704 int *isScalerNode, *rescaleDivisions, *norescaleDivisions;
1705 MrBFlt rescaleLnL;
1706 Tree *tree;
1707 TreeNode *p;
1708 ModelInfo *m;
1709
1710 if (beagleScalingScheme == MB_BEAGLE_SCALE_ALWAYS)
1711 {
1712
1713 # if defined (DEBUG_MB_BEAGLE_FLOW)
1714 MrBayesPrint ("ALWAYS RESCALING\n");
1715 # endif
1716
1717 for (d=0; d<divisionCount; d++)
1718 {
1719 dIndex = divisions[d];
1720 m = &modelSettings[dIndex];
1721
1722 /* Flip and copy or reset site scalers */
1723 FlipSiteScalerSpace(m, chain);
1724
1725 /* rescale if there is at least one partition requiring it */
1726 if (m->upDateAll == YES)
1727 {
1728 for (i=0; i<m->nCijkParts; i++)
1729 {
1730 beagleResetScaleFactorsByPartition(m->beagleInstance, m->siteScalerIndex[chain] + i, m->divisionIndex);
1731 }
1732 }
1733 else
1734 {
1735 CopySiteScalers(m, chain);
1736 }
1737 }
1738
1739 TreeTiProbs_BeagleMultiPartition(divisions, divisionCount, chain);
1740 TreeCondLikes_BeagleMultiPartition_Always_Rescale(divisions, divisionCount, chain);
1741 /* TODO integrate only divisions being updated */
1742 TreeLikelihood_BeagleMultiPartition(divisions, divisionCount, chain, lnL, (chainId[chain] % chainParams.numChains));
1743 }
1744 else
1745 { /* MB_BEAGLE_SCALE_DYNAMIC */
1746
1747 (*lnL) = 0.0;
1748 rescaleLnL = 0.0;
1749
1750 rescaleDivisions = (int *) SafeCalloc (numCurrentDivisions, sizeof(int));
1751 norescaleDivisions = (int *) SafeCalloc (numCurrentDivisions, sizeof(int));
1752 rescaleDivisionCount = 0;
1753 norescaleDivisionCount = 0;
1754
1755 TreeTiProbs_BeagleMultiPartition(divisions, divisionCount, chain);
1756
1757 for (d=0; d<divisionCount; d++)
1758 {
1759 dIndex = divisions[d];
1760 norescaleDivisions[d] = dIndex;
1761 m = &modelSettings[dIndex];
1762 tree = GetTree(m->brlens, chain, state[chain]);
1763
1764 /* This flag is only valid within this block */
1765 m->rescaleBeagleAll = NO;
1766
1767 if (m->successCount[chain] > 1000)
1768 {
1769 m->successCount[chain] = 10;
1770 m->rescaleFreq[chain]++; /* increase rescaleFreq independent of whether we accept or reject new state */
1771 m->rescaleFreqOld = m->rescaleFreqNew = m->rescaleFreq[chain];
1772 for (i=0; i<tree->nIntNodes; i++)
1773 {
1774 p = tree->intDownPass[i];
1775 if (p->upDateCl == YES) {
1776 /* flip to the new workspace since TreeCondLikes_Beagle_Rescale_All() does not do it for
1777 (p->upDateCl == YES) since it assumes that TreeCondLikes_Beagle_No_Rescale() did it */
1778 FlipCondLikeSpace (m, chain, p->index);
1779 }
1780 }
1781 rescaleDivisions[rescaleDivisionCount++] = dIndex;
1782 #if defined (DEBUG_MB_BEAGLE_MULTIPART)
1783 printf ("NUMERICAL RESCALING FOR DIVISION %d due to successes\n", dIndex);
1784 #endif
1785
1786 norescaleDivisions[d] = -1;
1787 }
1788 else if (beagleScalingFrequency != 0 && m->beagleComputeCount[chain] % (beagleScalingFrequency) == 0)
1789 {
1790 m->rescaleFreqOld = m->rescaleFreqNew = m->rescaleFreq[chain];
1791 for (i=0; i<tree->nIntNodes; i++)
1792 {
1793 p = tree->intDownPass[i];
1794 if (p->upDateCl == YES) {
1795 /* flip to the new workspace since TreeCondLikes_Beagle_Rescale_All() does not do it for
1796 (p->upDateCl == YES) since it assumes that TreeCondLikes_Beagle_No_Rescale() did it */
1797 FlipCondLikeSpace (m, chain, p->index);
1798 }
1799 }
1800 rescaleDivisions[rescaleDivisionCount++] = dIndex;
1801 #if defined (DEBUG_MB_BEAGLE_MULTIPART)
1802 printf ("NUMERICAL RESCALING FOR DIVISION %d due to beagleScalingFrequency\n", dIndex);
1803 #endif
1804 norescaleDivisions[d] = -1;
1805 }
1806 }
1807
1808 if (rescaleDivisionCount < divisionCount)
1809 {
1810 i = 0;
1811 for (d=0; d<divisionCount; d++)
1812 {
1813 dIndex = norescaleDivisions[d];
1814 if (dIndex != -1)
1815 {
1816 norescaleDivisions[i++] = dIndex;
1817 }
1818 }
1819 norescaleDivisionCount = i;
1820
1821 TreeCondLikes_BeagleMultiPartition_No_Rescale(norescaleDivisions, norescaleDivisionCount, chain);
1822
1823 beagleReturn = TreeLikelihood_BeagleMultiPartition(norescaleDivisions, norescaleDivisionCount, chain, lnL, (chainId[chain] % chainParams.numChains));
1824
1825 /* Check if likelihood is valid */
1826 if (beagleReturn == BEAGLE_ERROR_FLOATING_POINT)
1827 {
1828 (*lnL)=0;
1829 for (d=0; d<norescaleDivisionCount; d++)
1830 {
1831 dIndex = norescaleDivisions[d];
1832 m = &modelSettings[dIndex];
1833 if (m->lnLike[2*chain + state[chain]] > DBL_MAX ||
1834 m->lnLike[2*chain + state[chain]] < -DBL_MAX ||
1835 m->lnLike[2*chain + state[chain]] != m->lnLike[2*chain + state[chain]])
1836 {
1837
1838 m->rescaleFreqOld = m->rescaleFreqNew = m->rescaleFreq[chain];
1839 if (m->rescaleFreqNew > 1 && m->successCount[chain] < 40)
1840 {
1841 if (m->successCount[chain] < 10)
1842 {
1843 if (m->successCount[chain] < 4)
1844 {
1845 m->rescaleFreqNew -= m->rescaleFreqNew >> 3; /* <== we cut up to 12,5% of rescaleFreq */
1846 if (m->successCount[chain] < 2)
1847 {
1848 m->rescaleFreqNew -= m->rescaleFreqNew >> 3;
1849 /* 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*/
1850 /* 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 */
1851 /* 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*/
1852 m->rescaleFreqOld -= m->rescaleFreqOld >> 3;
1853 }
1854 m->rescaleFreqOld -= m->rescaleFreqOld >> 3;
1855 m->rescaleFreqOld--;
1856 m->rescaleFreqOld = (m->rescaleFreqOld ? m->rescaleFreqOld:1);
1857 m->recalculateScalers = YES;
1858 recalcScalers = YES;
1859 }
1860 }
1861 m->rescaleFreqNew--;
1862 m->rescaleFreqNew = (m->rescaleFreqNew ? m->rescaleFreqNew : 1);
1863 }
1864 m->successCount[chain] = 0;
1865 rescaleDivisions[rescaleDivisionCount++] = dIndex;
1866 norescaleDivisions[d] = -1;
1867 }
1868 else
1869 {
1870 (*lnL)+= m->lnLike[2*chain + state[chain]];
1871 }
1872 }
1873 // /* rescaling all divisions */
1874 // for (d=0; d<divisionCount; d++)
1875 // {
1876 // rescaleDivisions[d] = divisions[d];
1877 // }
1878 // rescaleDivisionCount = divisionCount;
1879 }
1880 }
1881
1882 if (rescaleDivisionCount > 0)
1883 {
1884 for (d=0; d<rescaleDivisionCount; d++)
1885 {
1886 dIndex = rescaleDivisions[d];
1887 m = &modelSettings[dIndex];
1888
1889 # if defined (DEBUG_MB_BEAGLE_FLOW)
1890 MrBayesPrint ("NUMERICAL RESCALING FOR DIVISION %d\n", dIndex);
1891 # endif
1892
1893 m->rescaleBeagleAll = YES;
1894 FlipSiteScalerSpace(m, chain);
1895 }
1896
1897 while_loop:
1898 for (d=0; d<rescaleDivisionCount; d++)
1899 {
1900 dIndex = rescaleDivisions[d];
1901 m = &modelSettings[dIndex];
1902
1903 isScalerNode = m->isScalerNode[chain];
1904 ResetScalersPartition (isScalerNode, tree, m->rescaleFreqNew);
1905 for (i=0; i<m->nCijkParts; i++)
1906 {
1907 beagleResetScaleFactorsByPartition(m->beagleInstance, m->siteScalerIndex[chain] + i, m->divisionIndex);
1908 }
1909 }
1910
1911 TreeCondLikes_BeagleMultiPartition_Rescale_All (rescaleDivisions, rescaleDivisionCount, chain);
1912 if (TreeLikelihood_BeagleMultiPartition(rescaleDivisions, rescaleDivisionCount, chain, &rescaleLnL, (chainId[chain] % chainParams.numChains)) == BEAGLE_ERROR_FLOATING_POINT)
1913 {
1914
1915 for (d=0; d<rescaleDivisionCount; d++)
1916 {
1917 dIndex = rescaleDivisions[d];
1918 m = &modelSettings[dIndex];
1919 tree = GetTree(m->brlens, chain, state[chain]);
1920 isScalerNode = m->isScalerNode[chain];
1921
1922 if (m->rescaleFreqNew > 1)
1923 {
1924 /* Swap back scalers which were swapped in TreeCondLikes_Beagle_Rescale_All() */
1925 for (i=0; i<tree->nIntNodes; i++)
1926 {
1927 p = tree->intDownPass[i];
1928 if (isScalerNode[p->index] == YES)
1929 FlipNodeScalerSpace (m, chain, p->index);
1930 }
1931 m->rescaleFreqNew -= m->rescaleFreqNew >> 3; /* <== we cut up to 12,5% of rescaleFreq */
1932 m->rescaleFreqNew--; /* we cut extra 1 of rescaleFreq */
1933 }
1934 else
1935 {
1936 m->rescaleFreq[chain] = m->rescaleFreqNew;
1937 rescaleDivisions[d] = -1;
1938 }
1939 }
1940
1941 i = 0;
1942 for (d=0; d<rescaleDivisionCount; d++)
1943 {
1944 dIndex = rescaleDivisions[d];
1945 if (dIndex != -1)
1946 {
1947 rescaleDivisions[i++] = dIndex;
1948 }
1949 }
1950 rescaleDivisionCount = i;
1951
1952 if (rescaleDivisionCount > 0)
1953 {
1954 goto while_loop;
1955 }
1956 }
1957
1958 for (d=0; d<rescaleDivisionCount; d++)
1959 {
1960 dIndex = rescaleDivisions[d];
1961 m = &modelSettings[dIndex];
1962 m->rescaleFreq[chain] = m->rescaleFreqNew;
1963 }
1964
1965 (*lnL) += rescaleLnL;
1966 }
1967
1968 for (d=0; d<divisionCount; d++)
1969 {
1970 dIndex = divisions[d];
1971 m = &modelSettings[dIndex];
1972 /* Count number of evaluations */
1973 m->beagleComputeCount[chain]++;
1974 }
1975
1976 free(rescaleDivisions);
1977 free(norescaleDivisions);
1978
1979 }
1980
1981 return;
1982 }
1983
1984
1985 /*----------------------------------------------------------------
1986 |
1987 | TreeTiProbs_BeagleMultiPartition: This routine updates all transition
1988 | probability matrices of a beagle instance across all divisions.
1989 |
1990 -----------------------------------------------------------------*/
TreeTiProbs_BeagleMultiPartition(int * divisions,int divisionCount,int chain)1991 int TreeTiProbs_BeagleMultiPartition (int* divisions, int divisionCount, int chain)
1992 {
1993 int i, j, k, d, dIndex, count, divisionOffset, divisionCijkIndex;
1994 MrBFlt correctionFactor, theRate, baseRate, *catRate, length, *branchLengthsAll;
1995 TreeNode *p;
1996 Tree *t;
1997 ModelInfo *m;
1998 int *cijkIndicesAll, *tiProbIndicesAll, *categoryRateIndicesAll;
1999
2000 cijkIndicesAll = modelSettings[0].cijkIndicesAll;
2001 categoryRateIndicesAll = modelSettings[0].categoryRateIndicesAll;
2002 tiProbIndicesAll = modelSettings[0].tiProbIndices;
2003 branchLengthsAll = modelSettings[0].branchLengths;
2004
2005 count = 0;
2006 divisionOffset = 0;
2007 for (d=0; d<divisionCount; d++)
2008 {
2009 dIndex = divisions[d];
2010
2011 /* get model settings */
2012 m = &modelSettings[dIndex];
2013 t = GetTree(m->brlens, chain, state[chain]);
2014 divisionOffset = m->numTiProbs * m->nCijkParts * m->divisionIndex;
2015 #if defined (DEBUG_MB_BEAGLE_MULTIPART)
2016 printf("m->upDateAll %d\n", m->upDateAll);
2017 #endif
2018 /* find the correction factor to make branch lengths
2019 in terms of expected number of substitutions per character */
2020 correctionFactor = 1.0;
2021 if (m->dataType == DNA || m->dataType == RNA)
2022 {
2023 if (m->nucModelId == NUCMODEL_DOUBLET)
2024 correctionFactor = 2.0;
2025 else if (m->nucModelId == NUCMODEL_CODON)
2026 correctionFactor = 3.0;
2027 }
2028
2029 /* get rate multipliers (for gamma & partition specific rates) */
2030 theRate = 1.0;
2031 baseRate = GetRate (dIndex, chain);
2032
2033 /* compensate for invariable sites if appropriate */
2034 if (m->pInvar != NULL)
2035 baseRate /= (1.0 - (*GetParamVals(m->pInvar, chain, state[chain])));
2036
2037 /* get category rates for gamma */
2038 if (m->shape == NULL)
2039 catRate = &theRate;
2040 else
2041 catRate = GetParamSubVals (m->shape, chain, state[chain]);
2042
2043 /* get effective category rates */
2044 for (k=0; k<m->numRateCats; k++)
2045 m->inRates[k] = baseRate * catRate[k] * correctionFactor;
2046
2047 /* TODO: only need to set category rates when they change */
2048 /* set category rates */
2049 beagleSetCategoryRatesWithIndex(m->beagleInstance, m->divisionIndex, m->inRates);
2050
2051 divisionCijkIndex = m->cijkIndex[chain] + (numLocalChains + 1) * m->nCijkParts * m->divisionIndex;
2052
2053 /* get ti prob indices and branch lengths to update */
2054 for (i=0; i<t->nNodes; i++)
2055 {
2056 p = t->allDownPass[i];
2057
2058 /* check if transition probs need updating */
2059 if (p->upDateTi == YES)
2060 {
2061 /* flip transition probability */
2062 FlipTiProbsSpace (m, chain, p->index);
2063
2064 /* find length */
2065 if (m->cppEvents != NULL)
2066 {
2067 length = GetParamSubVals (m->cppEvents, chain, state[chain])[p->index];
2068 }
2069 else if (m->tk02BranchRates != NULL)
2070 {
2071 length = GetParamSubVals (m->tk02BranchRates, chain, state[chain])[p->index];
2072 }
2073 else if (m->igrBranchRates != NULL)
2074 {
2075 length = GetParamSubVals (m->igrBranchRates, chain, state[chain])[p->index];
2076 }
2077 else if (m->mixedBrchRates != NULL)
2078 {
2079 length = GetParamSubVals (m->mixedBrchRates, chain, state[chain])[p->index];
2080 }
2081 else
2082 length = p->length;
2083
2084 /* numerical errors might ensue if we allow very large or very small branch lengths, which might
2085 occur in relaxed clock models; an elegant solution would be to substitute the stationary
2086 probs and initial probs but for now we truncate lengths at small or large values */
2087 if (length > BRLENS_MAX)
2088 length = BRLENS_MAX;
2089 else if (length < BRLENS_MIN)
2090 length = BRLENS_MIN;
2091
2092 branchLengthsAll[count] = length;
2093 tiProbIndicesAll[count] = m->tiProbsIndex[chain][p->index] + divisionOffset;
2094 categoryRateIndicesAll[count] = m->divisionIndex;
2095 cijkIndicesAll[count] = divisionCijkIndex;
2096 #if defined (DEBUG_MB_BEAGLE_MULTIPART)
2097 {
2098 printf("%2d %f ", count, branchLengthsAll[count]);
2099 printf(" %3d", tiProbIndicesAll[count]);
2100 printf(" %d", categoryRateIndicesAll[count]);
2101 printf(" %f %f %f %f", m->inRates[0], m->inRates[1], m->inRates[2], m->inRates[3]);
2102 printf(" %d\n", cijkIndicesAll[count]);
2103 }
2104 #endif
2105 count++;
2106 }
2107 }
2108 }
2109
2110 /* TODO: only need to update branches that have changed */
2111 /* calculate transition probabilities */
2112 if (count > 0)
2113 {
2114 for (i=1; i<m->nCijkParts; i++)
2115 {
2116 for (j=0; j<count; j++)
2117 {
2118 k = j + i*count;
2119 branchLengthsAll[k] = branchLengthsAll[j];
2120 tiProbIndicesAll[k] = tiProbIndicesAll[j] + i;
2121 categoryRateIndicesAll[k] = categoryRateIndicesAll[j];
2122 cijkIndicesAll[k] = cijkIndicesAll[j] + i;
2123 }
2124 }
2125
2126 beagleUpdateTransitionMatricesWithMultipleModels(m->beagleInstance,
2127 cijkIndicesAll,
2128 categoryRateIndicesAll,
2129 tiProbIndicesAll,
2130 NULL,
2131 NULL,
2132 branchLengthsAll,
2133 count * m->nCijkParts);
2134
2135 }
2136
2137 /* return success */
2138 return NO_ERROR;
2139 }
2140
2141
2142 /*----------------------------------------------------------------
2143 |
2144 | TreeCondLikes_Beagle: This routine updates all conditional
2145 | (partial) likelihoods of a beagle instance across all divisions while doing no rescaling.
2146 | That potentialy can make final liklihood bad then calculation with rescaling needs to be done.
2147 |
2148 -----------------------------------------------------------------*/
TreeCondLikes_BeagleMultiPartition_No_Rescale(int * divisions,int divisionCount,int chain)2149 int TreeCondLikes_BeagleMultiPartition_No_Rescale (int* divisions, int divisionCount, int chain)
2150 {
2151 int i, j, d, dIndex, divisionOffset;
2152 int opJ, opPrev, opCountMax, opCountTotal;
2153 Tree *t;
2154 TreeNode *p;
2155 ModelInfo *m;
2156 unsigned chil1Step, chil2Step;
2157 int *isScalerNode;
2158 BeagleOperationByPartition *operationsAll;
2159
2160 m = &modelSettings[0];
2161 operationsAll = m->operationsAll;
2162 opCountMax = 0;
2163
2164 for (d=0; d<divisionCount; d++)
2165 {
2166 dIndex = divisions[d];
2167 m = &modelSettings[dIndex];
2168 t = GetTree(m->brlens, chain, state[chain]);
2169 divisionOffset = m->numTiProbs * m->nCijkParts * m->divisionIndex;
2170
2171 isScalerNode = m->isScalerNode[chain];
2172
2173 m->opCount = 0;
2174
2175 for (i=0; i<t->nIntNodes; i++)
2176 {
2177
2178 if (t->levelPassEnabled)
2179 {
2180 p = t->intDownPassLevel[i];
2181 }
2182 else
2183 {
2184 p = t->intDownPass[i];
2185 }
2186
2187 /* check if conditional likelihoods need updating */
2188 if (p->upDateCl == YES)
2189 {
2190 /* flip to the new workspace */
2191 FlipCondLikeSpace (m, chain, p->index);
2192
2193 /* update conditional likelihoods */
2194 m->operationsByPartition[m->opCount].destinationPartials = m->condLikeIndex[chain][p->index ];
2195 m->operationsByPartition[m->opCount].child1Partials = m->condLikeIndex[chain][p->left->index ];
2196 m->operationsByPartition[m->opCount].child1TransitionMatrix = m->tiProbsIndex [chain][p->left->index ] + divisionOffset;
2197 m->operationsByPartition[m->opCount].child2Partials = m->condLikeIndex[chain][p->right->index];
2198 m->operationsByPartition[m->opCount].child2TransitionMatrix = m->tiProbsIndex [chain][p->right->index] + divisionOffset;
2199
2200 /* All partials for tips are the same across omega categories, thus we are doing the following two if statments.*/
2201 if (p->left->left== NULL)
2202 chil1Step=0;
2203 else
2204 chil1Step=1;
2205
2206 if (p->right->left== NULL)
2207 chil2Step=0;
2208 else
2209 chil2Step=1;
2210
2211 m->operationsByPartition[m->opCount].destinationScaleWrite = BEAGLE_OP_NONE;
2212 m->operationsByPartition[m->opCount].cumulativeScaleIndex = BEAGLE_OP_NONE;
2213 if (isScalerNode[p->index] == YES)
2214 {
2215 m->operationsByPartition[m->opCount].destinationScaleRead = m->nodeScalerIndex[chain][p->index];
2216 }
2217 else
2218 {
2219 m->operationsByPartition[m->opCount].destinationScaleRead = BEAGLE_OP_NONE;
2220 }
2221
2222 for (j=1; j<m->nCijkParts; j++)
2223 {
2224 opJ = m->opCount + j * t->nIntNodes;
2225 opPrev = m->opCount + (j-1) * t->nIntNodes;
2226
2227 m->operationsByPartition[opJ].destinationPartials = m->operationsByPartition[opPrev].destinationPartials + 1;
2228 m->operationsByPartition[opJ].child1Partials = m->operationsByPartition[opPrev].child1Partials + chil1Step;
2229 m->operationsByPartition[opJ].child1TransitionMatrix = m->operationsByPartition[opPrev].child1TransitionMatrix + 1;
2230 m->operationsByPartition[opJ].child2Partials = m->operationsByPartition[opPrev].child2Partials + chil2Step;
2231 m->operationsByPartition[opJ].child2TransitionMatrix = m->operationsByPartition[opPrev].child2TransitionMatrix + 1;
2232
2233 m->operationsByPartition[opJ].destinationScaleWrite = BEAGLE_OP_NONE;
2234 m->operationsByPartition[opJ].cumulativeScaleIndex = BEAGLE_OP_NONE;
2235 if ( isScalerNode[p->index] == YES )
2236 {
2237 m->operationsByPartition[opJ].destinationScaleRead = m->operationsByPartition[opPrev].destinationScaleRead + 1;
2238 }
2239 else
2240 {
2241 m->operationsByPartition[opJ].destinationScaleRead = BEAGLE_OP_NONE;
2242 }
2243 }
2244 m->opCount++;
2245 }
2246 }
2247 if (m->opCount > opCountMax)
2248 opCountMax = m->opCount;
2249 }
2250
2251 opCountTotal = 0;
2252 for (i=0; i<opCountMax; i++)
2253 {
2254 for (j=0; j<m->nCijkParts; j++)
2255 {
2256 opJ = i + j*t->nIntNodes;
2257 for (d=0; d<divisionCount; d++)
2258 {
2259 dIndex = divisions[d];
2260 m = &modelSettings[dIndex];
2261 if (i<m->opCount)
2262 {
2263 operationsAll[opCountTotal].destinationPartials = m->operationsByPartition[opJ].destinationPartials;
2264 operationsAll[opCountTotal].destinationScaleWrite = m->operationsByPartition[opJ].destinationScaleWrite;
2265 operationsAll[opCountTotal].destinationScaleRead = m->operationsByPartition[opJ].destinationScaleRead;
2266 operationsAll[opCountTotal].child1Partials = m->operationsByPartition[opJ].child1Partials;
2267 operationsAll[opCountTotal].child1TransitionMatrix = m->operationsByPartition[opJ].child1TransitionMatrix;
2268 operationsAll[opCountTotal].child2Partials = m->operationsByPartition[opJ].child2Partials;
2269 operationsAll[opCountTotal].child2TransitionMatrix = m->operationsByPartition[opJ].child2TransitionMatrix;
2270 operationsAll[opCountTotal].partition = m->divisionIndex;
2271 operationsAll[opCountTotal].cumulativeScaleIndex = m->operationsByPartition[opJ].cumulativeScaleIndex;
2272 opCountTotal++;
2273 }
2274 }
2275 }
2276 }
2277
2278 beagleUpdatePartialsByPartition(m->beagleInstance,
2279 operationsAll,
2280 opCountTotal);
2281
2282 return NO_ERROR;
2283 }
2284
2285
2286 /*----------------------------------------------------------------
2287 |
2288 | TreeCondLikes_Beagle: This routine updates all conditional
2289 | (partial) likelihoods of a beagle instance across all divisions while rescaling at every node.
2290 | Note: all nodes get recalculated, not only tached by move.
2291 |
2292 -----------------------------------------------------------------*/
TreeCondLikes_BeagleMultiPartition_Rescale_All(int * divisions,int divisionCount,int chain)2293 int TreeCondLikes_BeagleMultiPartition_Rescale_All (int* divisions, int divisionCount, int chain)
2294 {
2295 int i, j, d, dIndex, divisionOffset;
2296 int opJ, opPrev, opCountMax, opCountTotal;
2297 Tree *t;
2298 TreeNode *p;
2299 ModelInfo *m;
2300 unsigned chil1Step, chil2Step;
2301 int *isScalerNode;
2302 BeagleOperationByPartition *operationsAll;
2303
2304 m = &modelSettings[0];
2305 operationsAll = m->operationsAll;
2306 opCountMax = 0;
2307
2308 t = GetTree(m->brlens, chain, state[chain]);
2309
2310 for (d=0; d<divisionCount; d++)
2311 {
2312 dIndex = divisions[d];
2313 m = &modelSettings[dIndex];
2314 t = GetTree(m->brlens, chain, state[chain]);
2315 divisionOffset = m->numTiProbs * m->nCijkParts * m->divisionIndex;
2316
2317 isScalerNode = m->isScalerNode[chain];
2318
2319 m->opCount = 0;
2320
2321 for (i=0; i<t->nIntNodes; i++)
2322 {
2323
2324 if (t->levelPassEnabled)
2325 {
2326 p = t->intDownPassLevel[i];
2327 }
2328 else
2329 {
2330 p = t->intDownPass[i];
2331 }
2332
2333 if (p->upDateCl == NO)
2334 {
2335 //p->upDateCl = YES;
2336 /* flip to the new workspace */
2337 FlipCondLikeSpace (m, chain, p->index);
2338 }
2339
2340 /* update conditional likelihoods */
2341 m->operationsByPartition[m->opCount].destinationPartials = m->condLikeIndex[chain][p->index ];
2342 m->operationsByPartition[m->opCount].child1Partials = m->condLikeIndex[chain][p->left->index ];
2343 m->operationsByPartition[m->opCount].child1TransitionMatrix = m->tiProbsIndex [chain][p->left->index ] + divisionOffset;
2344 m->operationsByPartition[m->opCount].child2Partials = m->condLikeIndex[chain][p->right->index];
2345 m->operationsByPartition[m->opCount].child2TransitionMatrix = m->tiProbsIndex [chain][p->right->index] + divisionOffset;
2346
2347 /* All partials for tips are the same across omega catigoris, thus we are doing the following two if statments.*/
2348 if (p->left->left== NULL)
2349 chil1Step=0;
2350 else
2351 chil1Step=1;
2352
2353 if (p->right->left== NULL)
2354 chil2Step=0;
2355 else
2356 chil2Step=1;
2357
2358 m->operationsByPartition[m->opCount].destinationScaleRead = BEAGLE_OP_NONE;
2359 if (isScalerNode[p->index] == YES)
2360 {
2361 FlipNodeScalerSpace (m, chain, p->index);
2362 m->operationsByPartition[m->opCount].destinationScaleWrite = m->nodeScalerIndex[chain][p->index];
2363 m->operationsByPartition[m->opCount].cumulativeScaleIndex = m->siteScalerIndex[chain];
2364 }
2365 else
2366 {
2367 m->operationsByPartition[m->opCount].destinationScaleWrite = BEAGLE_OP_NONE;
2368 m->operationsByPartition[m->opCount].cumulativeScaleIndex = BEAGLE_OP_NONE;
2369 }
2370
2371 for (j=1; j<m->nCijkParts; j++)
2372 {
2373 opJ = m->opCount + j * t->nIntNodes;
2374 opPrev = m->opCount + (j-1) * t->nIntNodes;
2375
2376 m->operationsByPartition[opJ].destinationPartials = m->operationsByPartition[opPrev].destinationPartials + 1;
2377 m->operationsByPartition[opJ].child1Partials = m->operationsByPartition[opPrev].child1Partials + chil1Step;
2378 m->operationsByPartition[opJ].child1TransitionMatrix = m->operationsByPartition[opPrev].child1TransitionMatrix + 1;
2379 m->operationsByPartition[opJ].child2Partials = m->operationsByPartition[opPrev].child2Partials + chil2Step;
2380 m->operationsByPartition[opJ].child2TransitionMatrix = m->operationsByPartition[opPrev].child2TransitionMatrix + 1;
2381
2382 m->operationsByPartition[opJ].destinationScaleRead = BEAGLE_OP_NONE;
2383 if ( isScalerNode[p->index] == YES )
2384 {
2385 m->operationsByPartition[opJ].destinationScaleWrite = m->operationsByPartition[opPrev].destinationScaleWrite + 1;
2386 m->operationsByPartition[opJ].cumulativeScaleIndex = m->operationsByPartition[opPrev].cumulativeScaleIndex + 1;
2387 }
2388 else
2389 {
2390 m->operationsByPartition[opJ].destinationScaleWrite = BEAGLE_OP_NONE;
2391 m->operationsByPartition[opJ].cumulativeScaleIndex = BEAGLE_OP_NONE;
2392 }
2393 }
2394 m->opCount++;
2395 }
2396 if (m->opCount > opCountMax)
2397 opCountMax = m->opCount;
2398 }
2399
2400
2401 for (j=0; j<m->nCijkParts; j++)
2402 {
2403 opCountTotal = 0;
2404 for (i=0; i<opCountMax; i++)
2405 {
2406 opJ = i + j*t->nIntNodes;
2407 for (d=0; d<divisionCount; d++)
2408 {
2409 dIndex = divisions[d];
2410 m = &modelSettings[dIndex];
2411 if (i<m->opCount)
2412 {
2413 operationsAll[opCountTotal].destinationPartials = m->operationsByPartition[opJ].destinationPartials;
2414 operationsAll[opCountTotal].destinationScaleWrite = m->operationsByPartition[opJ].destinationScaleWrite;
2415 operationsAll[opCountTotal].destinationScaleRead = m->operationsByPartition[opJ].destinationScaleRead;
2416 operationsAll[opCountTotal].child1Partials = m->operationsByPartition[opJ].child1Partials;
2417 operationsAll[opCountTotal].child1TransitionMatrix = m->operationsByPartition[opJ].child1TransitionMatrix;
2418 operationsAll[opCountTotal].child2Partials = m->operationsByPartition[opJ].child2Partials;
2419 operationsAll[opCountTotal].child2TransitionMatrix = m->operationsByPartition[opJ].child2TransitionMatrix;
2420 operationsAll[opCountTotal].partition = m->divisionIndex;
2421 operationsAll[opCountTotal].cumulativeScaleIndex = m->operationsByPartition[opJ].cumulativeScaleIndex;
2422 opCountTotal++;
2423 }
2424 }
2425 }
2426 beagleUpdatePartialsByPartition(m->beagleInstance,
2427 operationsAll,
2428 opCountTotal);
2429 }
2430
2431 return NO_ERROR;
2432 }
2433
2434
2435 /*----------------------------------------------------------------
2436 |
2437 | TreeCondLikes_BeagleMultiPartition: This routine updates and rescales all conditional
2438 | (partial) likelihoods of a beagle instance across all divisions.
2439 |
2440 -----------------------------------------------------------------*/
TreeCondLikes_BeagleMultiPartition_Always_Rescale(int * divisions,int divisionCount,int chain)2441 int TreeCondLikes_BeagleMultiPartition_Always_Rescale (int* divisions, int divisionCount, int chain)
2442 {
2443 int i, j, d, dIndex, destinationScaleRead, cumulativeScaleIndex;
2444 int opJ, opPrev, opCountMax, opCountTotal, divisionOffset, scaleOp;
2445 Tree *t;
2446 TreeNode *p;
2447 ModelInfo *m;
2448 unsigned chil1Step, chil2Step;
2449 BeagleOperationByPartition *operationsAll;
2450
2451 m = &modelSettings[0];
2452 operationsAll = m->operationsAll;
2453 opCountMax = 0;
2454
2455 for (d=0; d<divisionCount; d++)
2456 {
2457 dIndex = divisions[d];
2458
2459 /* get model settings */
2460 m = &modelSettings[dIndex];
2461 t = GetTree(m->brlens, chain, state[chain]);
2462 divisionOffset = m->numTiProbs * m->nCijkParts * m->divisionIndex;
2463
2464 m->opCount = 0;
2465 scaleOp = 0;
2466
2467 for (i=0; i<t->nIntNodes; i++)
2468 {
2469
2470 if (t->levelPassEnabled)
2471 {
2472 p = t->intDownPassLevel[i];
2473 }
2474 else
2475 {
2476 p = t->intDownPass[i];
2477 }
2478
2479
2480 /* check if conditional likelihoods need updating */
2481 if (p->upDateCl == YES)
2482 {
2483
2484 /* remove old scalers */
2485 if (m->upDateAll == NO)
2486 {
2487 destinationScaleRead = m->nodeScalerIndex[chain][p->index];
2488 for (j=0; j<m->nCijkParts; j++)
2489 {
2490 m->scaleFactorsOps[scaleOp+j*t->nIntNodes] = destinationScaleRead;
2491 destinationScaleRead++;
2492 }
2493 scaleOp++;
2494 }
2495
2496
2497 /* flip to the new workspace */
2498 FlipCondLikeSpace (m, chain, p->index);
2499 FlipNodeScalerSpace (m, chain, p->index);
2500
2501 /* update conditional likelihoods */
2502 m->operationsByPartition[m->opCount].destinationPartials = m->condLikeIndex[chain][p->index ];
2503 m->operationsByPartition[m->opCount].child1Partials = m->condLikeIndex[chain][p->left->index ];
2504 m->operationsByPartition[m->opCount].child1TransitionMatrix = m->tiProbsIndex [chain][p->left->index ] + divisionOffset;
2505 m->operationsByPartition[m->opCount].child2Partials = m->condLikeIndex[chain][p->right->index];
2506 m->operationsByPartition[m->opCount].child2TransitionMatrix = m->tiProbsIndex [chain][p->right->index] + divisionOffset;
2507
2508 /* All partials for tips are the same across omega catigoris, thus we are doing the following two if statments.*/
2509 if (p->left->left== NULL && p->left->right== NULL)
2510 chil1Step=0;
2511 else
2512 chil1Step=1;
2513
2514 if (p->right->left== NULL && p->right->right== NULL)
2515 chil2Step=0;
2516 else
2517 chil2Step=1;
2518
2519 m->operationsByPartition[m->opCount].destinationScaleWrite = m->nodeScalerIndex[chain][p->index];
2520 m->operationsByPartition[m->opCount].cumulativeScaleIndex = m->siteScalerIndex[chain];
2521 m->operationsByPartition[m->opCount].destinationScaleRead = BEAGLE_OP_NONE;
2522
2523 for (j=1; j<m->nCijkParts; j++)
2524 {
2525 opJ = m->opCount + j * t->nIntNodes;
2526 opPrev = m->opCount + (j-1) * t->nIntNodes;
2527
2528 m->operationsByPartition[opJ].destinationPartials = m->operationsByPartition[opPrev].destinationPartials + 1;
2529 m->operationsByPartition[opJ].child1Partials = m->operationsByPartition[opPrev].child1Partials + chil1Step;
2530 m->operationsByPartition[opJ].child1TransitionMatrix = m->operationsByPartition[opPrev].child1TransitionMatrix + 1;
2531 m->operationsByPartition[opJ].child2Partials = m->operationsByPartition[opPrev].child2Partials + chil2Step;
2532 m->operationsByPartition[opJ].child2TransitionMatrix = m->operationsByPartition[opPrev].child2TransitionMatrix + 1;
2533
2534 m->operationsByPartition[opJ].destinationScaleWrite = m->operationsByPartition[opPrev].destinationScaleWrite + 1;
2535 m->operationsByPartition[opJ].cumulativeScaleIndex = m->operationsByPartition[opPrev].cumulativeScaleIndex + 1;
2536 m->operationsByPartition[opJ].destinationScaleRead = BEAGLE_OP_NONE;
2537 }
2538 m->opCount++;
2539 }
2540 } /* end of nIntNodes for */
2541
2542 if (m->opCount > opCountMax)
2543 {
2544 opCountMax = m->opCount;
2545 }
2546
2547 if (m->upDateAll == NO)
2548 {
2549 cumulativeScaleIndex = m->siteScalerIndex[chain];
2550 for (j=0; j<m->nCijkParts; j++)
2551 {
2552 beagleRemoveScaleFactorsByPartition(m->beagleInstance,
2553 &m->scaleFactorsOps[j*t->nIntNodes],
2554 scaleOp,
2555 cumulativeScaleIndex,
2556 m->divisionIndex);
2557 cumulativeScaleIndex++;
2558 }
2559 }
2560 }
2561
2562 for (j=0; j<m->nCijkParts; j++)
2563 {
2564 opCountTotal = 0;
2565 for (i=0; i<opCountMax; i++)
2566 {
2567 opJ = i + j*t->nIntNodes;
2568 for (d=0; d<divisionCount; d++)
2569 {
2570 dIndex = divisions[d];
2571 m = &modelSettings[dIndex];
2572 if (i<m->opCount)
2573 {
2574 operationsAll[opCountTotal].destinationPartials = m->operationsByPartition[opJ].destinationPartials;
2575 operationsAll[opCountTotal].destinationScaleWrite = m->operationsByPartition[opJ].destinationScaleWrite;
2576 operationsAll[opCountTotal].destinationScaleRead = m->operationsByPartition[opJ].destinationScaleRead;
2577 operationsAll[opCountTotal].child1Partials = m->operationsByPartition[opJ].child1Partials;
2578 operationsAll[opCountTotal].child1TransitionMatrix = m->operationsByPartition[opJ].child1TransitionMatrix;
2579 operationsAll[opCountTotal].child2Partials = m->operationsByPartition[opJ].child2Partials;
2580 operationsAll[opCountTotal].child2TransitionMatrix = m->operationsByPartition[opJ].child2TransitionMatrix;
2581 operationsAll[opCountTotal].partition = m->divisionIndex;
2582 operationsAll[opCountTotal].cumulativeScaleIndex = m->operationsByPartition[opJ].cumulativeScaleIndex;
2583
2584 #if defined (DEBUG_MB_BEAGLE_MULTIPART)
2585 {
2586 printf("%2d %3d ", opCountTotal, operationsAll[opCountTotal].destinationPartials);
2587 printf("%3d ", operationsAll[opCountTotal].destinationScaleWrite);
2588 printf("%3d ", operationsAll[opCountTotal].destinationScaleRead);
2589 printf("%3d ", operationsAll[opCountTotal].child1Partials);
2590 printf("%3d ", operationsAll[opCountTotal].child1TransitionMatrix);
2591 printf("%3d ", operationsAll[opCountTotal].child2Partials);
2592 printf("%3d ", operationsAll[opCountTotal].child2TransitionMatrix);
2593 printf("%d ", operationsAll[opCountTotal].partition);
2594 printf("%3d\n", operationsAll[opCountTotal].cumulativeScaleIndex);
2595 }
2596 #endif
2597
2598 opCountTotal++;
2599 }
2600 }
2601 }
2602 beagleUpdatePartialsByPartition(m->beagleInstance,
2603 operationsAll,
2604 opCountTotal);
2605 }
2606
2607 return NO_ERROR;
2608 }
2609
2610 /**---------------------------------------------------------------------------
2611 |
2612 | TreeLikelihood_BeagleMultiPartition: Accumulate the log likelihoods calculated by Beagle
2613 | at the root across all divisions.
2614 |
2615 ---------------------------------------- -------------------------------------*/
TreeLikelihood_BeagleMultiPartition(int * divisions,int divisionCount,int chain,MrBFlt * lnL,int whichSitePats)2616 int TreeLikelihood_BeagleMultiPartition (int* divisions, int divisionCount, int chain, MrBFlt *lnL, int whichSitePats)
2617 {
2618 int i, j, d, c = 0, nStates, beagleReturn, site, dIndex, divisionOffset;
2619 int hasPInvar, hasAnyPInvar, hasAnyDataRestriction;
2620 MrBFlt *swr, s01, s10, probOn, probOff, covBF[40], pInvar=0.0, *bs, freq, likeI, lnLikeI, diff, *omegaCatFreq, *lnLDiv;
2621 CLFlt *clInvar=NULL, *nSitesOfPat;
2622 double *nSitesOfPat_Beagle;
2623 TreeNode *p;
2624 Tree *t, *tZero;
2625 ModelInfo *m;
2626 double pUnobserved;
2627
2628 # if defined (MB_PRINT_DYNAMIC_RESCALE_FAIL_STAT)
2629 static unsigned countBeagleDynamicFail=0;
2630 static unsigned countALL=0;
2631 # endif
2632
2633 hasAnyPInvar = NO;
2634 hasAnyDataRestriction = NO;
2635
2636 m = &modelSettings[0];
2637
2638 for (d=0; d<divisionCount; d++)
2639 {
2640 dIndex = divisions[d];
2641
2642 /* find model settings and nStates, pInvar, invar cond likes */
2643 m = &modelSettings[dIndex];
2644
2645 divisionOffset = (numLocalChains + 1) * m->nCijkParts * m->divisionIndex;
2646
2647 t = GetTree(m->brlens, chain, state[chain]);
2648
2649 /* find root node */
2650 p = t->root->left;
2651
2652 nStates = m->numModelStates;
2653 if (m->pInvar == NULL)
2654 {
2655 hasPInvar = NO;
2656 }
2657 else
2658 {
2659 hasPInvar = YES;
2660 hasAnyPInvar = YES;
2661 pInvar = *(GetParamVals (m->pInvar, chain, state[chain]));
2662
2663 #if defined (DEBUG_MB_BEAGLE_MULTIPART)
2664 printf("pInvar[%d] = %f\n", dIndex, *(GetParamVals (m->pInvar, chain, state[chain])));
2665 #endif
2666 }
2667
2668 if (m->dataType == RESTRICTION)
2669 hasAnyDataRestriction = YES;
2670
2671 /* find base frequencies */
2672 bs = GetParamSubVals (m->stateFreq, chain, state[chain]);
2673
2674 /* if covarion model, adjust base frequencies */
2675 if (m->switchRates != NULL)
2676 {
2677 /* find the stationary frequencies */
2678 swr = GetParamVals(m->switchRates, chain, state[chain]);
2679 s01 = swr[0];
2680 s10 = swr[1];
2681 probOn = s01 / (s01 + s10);
2682 probOff = 1.0 - probOn;
2683
2684 /* now adjust the base frequencies; on-state stored first in cond likes */
2685 for (j=0; j<nStates/2; j++)
2686 {
2687 covBF[j] = bs[j] * probOn;
2688 covBF[j+nStates/2] = bs[j] * probOff;
2689 }
2690
2691 /* finally set bs pointer to adjusted values */
2692 bs = covBF;
2693 }
2694
2695 /* TODO Really only need to check if state frequencies have changed */
2696 if (m->upDateCijk == YES)
2697 {
2698
2699 /* set base frequencies in BEAGLE instance */
2700 for (i=0; i<m->nCijkParts; i++)
2701 beagleSetStateFrequencies(m->beagleInstance,
2702 m->cijkIndex[chain] + i + divisionOffset,
2703 bs);
2704 }
2705
2706 /* find category frequencies */
2707 if (hasPInvar == NO)
2708 freq = 1.0 / m->numRateCats;
2709 else
2710 freq = (1.0 - pInvar) / m->numRateCats;
2711
2712 #if defined (DEBUG_MB_BEAGLE_MULTIPART)
2713 printf("freq = %f\n", freq);
2714 #endif
2715
2716 /* TODO: cat weights only need to be set when they change */
2717 /* set category frequencies in beagle instance */
2718 if (m->numOmegaCats > 1)
2719 {
2720 omegaCatFreq = GetParamSubVals(m->omega, chain, state[chain]);
2721 for (i=0; i<m->nCijkParts; i++)
2722 {
2723 for (j=0; j<m->numRateCats; j++)
2724 m->inWeights[j] = freq * omegaCatFreq[i];
2725 beagleSetCategoryWeights(m->beagleInstance,
2726 m->cijkIndex[chain] + i + divisionOffset,
2727 m->inWeights);
2728 }
2729 }
2730 else if (hasPInvar == YES)
2731 {
2732 for (i=0; i<m->numRateCats; i++)
2733 m->inWeights[i] = freq;
2734 beagleSetCategoryWeights(m->beagleInstance,
2735 m->cijkIndex[chain] + divisionOffset,
2736 m->inWeights);
2737 }
2738 }
2739
2740 /* TODO: pattern weights only need to be set when they change */
2741 /* set pattern weights in beagle instance if using dynamic reweighting */
2742 if (chainParams.weightScheme[0] + chainParams.weightScheme[1] > ETA)
2743 {
2744 nSitesOfPat_Beagle = (double *) SafeMalloc (m->numCharsAll * sizeof(double));
2745 site = 0;
2746 for (d=0; d<numCurrentDivisions; d++)
2747 {
2748 for (c=0; c < modelSettings[d].numChars; c++)
2749 {
2750 nSitesOfPat_Beagle[site++] = numSitesOfPat[modelSettings[d].compCharStart + c];
2751 }
2752 }
2753 beagleSetPatternWeights(m->beagleInstance,
2754 nSitesOfPat_Beagle);
2755 SAFEFREE (nSitesOfPat_Beagle);
2756 }
2757
2758 tZero = GetTree(modelSettings[0].brlens, chain, state[chain]);
2759
2760 /* find root log likelihoods and scalers */
2761 j=0;
2762 for (i=0; i<m->nCijkParts; i++)
2763 {
2764 for (d=0; d<divisionCount; d++)
2765 {
2766 dIndex = divisions[d];
2767 m = &modelSettings[dIndex];
2768 t = GetTree(m->brlens, chain, state[chain]);
2769 p = t->root->left;
2770
2771
2772 divisionOffset = (numLocalChains + 1) * m->nCijkParts * m->divisionIndex;
2773
2774 modelSettings[0].bufferIndices[j] = m->condLikeIndex[chain][p->index] + i;
2775 modelSettings[0].eigenIndices[j] = m->cijkIndex[chain] + i + divisionOffset;
2776 modelSettings[0].cumulativeScaleIndices[j] = m->siteScalerIndex[chain] + i;
2777
2778 #if defined (DEBUG_MB_BEAGLE_MULTIPART)
2779 printf("dIndex = %d, bufferIndices = %d, eigenIndices = %d, cumulativeScaleIndices = %d\n", dIndex, modelSettings[0].bufferIndices[j], modelSettings[0].eigenIndices[j], modelSettings[0].cumulativeScaleIndices[j]);
2780 #endif
2781 if (tZero->isRooted == NO)
2782 {
2783 divisionOffset = m->numTiProbs * m->nCijkParts * m->divisionIndex;
2784
2785 modelSettings[0].childBufferIndices[j] = m->condLikeIndex[chain][p->anc->index];
2786 modelSettings[0].childTiProbIndices[j] = m->tiProbsIndex [chain][p->index] + i + divisionOffset;
2787
2788 #if defined (DEBUG_MB_BEAGLE_MULTIPART)
2789 printf("childBufferIndices = %d, childTiProbIndices = %d\n", modelSettings[0].childBufferIndices[j], modelSettings[0].childTiProbIndices[j]);
2790 #endif
2791 }
2792 j++;
2793 }
2794
2795 }
2796
2797 /* reset lnL */
2798 *lnL = 0.0;
2799
2800 /* get root log likelihoods */
2801 if (tZero->isRooted == YES)
2802 {
2803 beagleReturn = beagleCalculateRootLogLikelihoodsByPartition(modelSettings[0].beagleInstance,
2804 modelSettings[0].bufferIndices,
2805 modelSettings[0].eigenIndices,
2806 modelSettings[0].eigenIndices,
2807 modelSettings[0].cumulativeScaleIndices,
2808 divisions,
2809 divisionCount,
2810 modelSettings[0].nCijkParts,
2811 modelSettings[0].logLikelihoodsAll, /* assuming numChars >= numDivisions */
2812 lnL);
2813
2814 }
2815 else
2816 {
2817 beagleReturn = beagleCalculateEdgeLogLikelihoodsByPartition(
2818 modelSettings[0].beagleInstance,
2819 modelSettings[0].bufferIndices,
2820 modelSettings[0].childBufferIndices,
2821 modelSettings[0].childTiProbIndices,
2822 NULL,
2823 NULL,
2824 modelSettings[0].eigenIndices,
2825 modelSettings[0].eigenIndices,
2826 modelSettings[0].cumulativeScaleIndices,
2827 divisions,
2828 divisionCount,
2829 modelSettings[0].nCijkParts,
2830 modelSettings[0].logLikelihoodsAll,
2831 lnL,
2832 NULL,
2833 NULL,
2834 NULL,
2835 NULL);
2836
2837 }
2838
2839 # if defined (MB_PRINT_DYNAMIC_RESCALE_FAIL_STAT)
2840 countALL++;
2841 # endif
2842
2843 if (*lnL > DBL_MAX || *lnL < -DBL_MAX) {
2844 beagleReturn = BEAGLE_ERROR_FLOATING_POINT;
2845 }
2846
2847 for (d=0; d<divisionCount; d++)
2848 {
2849 dIndex = divisions[d];
2850 m = &modelSettings[dIndex];
2851 m->lnLike[2*chain + state[chain]] = modelSettings[0].logLikelihoodsAll[d];
2852 if (m->lnLike[2*chain + state[chain]] > DBL_MAX ||
2853 m->lnLike[2*chain + state[chain]] < -DBL_MAX ||
2854 m->lnLike[2*chain + state[chain]] != m->lnLike[2*chain + state[chain]])
2855 {
2856 beagleReturn = BEAGLE_ERROR_FLOATING_POINT;
2857 }
2858 else
2859 {
2860 m->successCount[chain]++;
2861 }
2862 }
2863
2864 if (beagleReturn == BEAGLE_ERROR_FLOATING_POINT)
2865 {
2866 # if defined (MB_PRINT_DYNAMIC_RESCALE_FAIL_STAT)
2867 countBeagleDynamicFail++;
2868 MrBayesPrint ("DEBUG INFO (not an error) countBeagleDynamicFail:%d countALL:%d\n", countBeagleDynamicFail, countALL);
2869 # endif
2870 }
2871
2872 #if defined (DEBUG_MB_BEAGLE_MULTIPART_SITELNL)
2873 beagleGetSiteLogLikelihoods(modelSettings[0].beagleInstance, modelSettings[0].logLikelihoodsAll);
2874 printf("lnL = ");
2875 site = 0;
2876 for (d=0; d<numCurrentDivisions; d++)
2877 {
2878 for (c=0; c<modelSettings[d].numChars; c++)
2879 {
2880 printf("[%d] %f ", c, modelSettings[0].logLikelihoodsAll[site]);
2881 site++;
2882 }
2883 }
2884 printf("\n");
2885 #endif
2886
2887 /* accumulate logs across sites */
2888 if (hasAnyPInvar || hasAnyDataRestriction)
2889 {
2890 beagleGetSiteLogLikelihoods(modelSettings[0].beagleInstance, modelSettings[0].logLikelihoodsAll);
2891
2892 (*lnL) = 0.0;
2893 for (d=0; d<divisionCount; d++)
2894 {
2895 dIndex = divisions[d];
2896 m = &modelSettings[dIndex];
2897 lnLDiv = &m->lnLike[2*chain + state[chain]];
2898
2899 if (!(*lnLDiv > DBL_MAX || *lnLDiv < -DBL_MAX || *lnLDiv != *lnLDiv))
2900 {
2901 /* find nSitesOfPat */
2902 nSitesOfPat = numSitesOfPat + (whichSitePats*numCompressedChars) + m->compCharStart;
2903
2904 (*lnLDiv) = 0.0;
2905
2906 site = 0;
2907 for (i=0; i<dIndex; i++) {
2908 site += modelSettings[i].numChars;
2909 }
2910 if (m->pInvar == NULL)
2911 {
2912 if (m->dataType == RESTRICTION)
2913 {
2914 pUnobserved = 0.0;
2915 for (c=0; c<m->numDummyChars; c++)
2916 {
2917 pUnobserved += exp((double)modelSettings[0].logLikelihoodsAll[site]);
2918 site++;
2919 }
2920 /* correct for absent characters */
2921 (*lnLDiv) -= log (1-pUnobserved) * (m->numUncompressedChars);
2922 for (; c<m->numChars; c++)
2923 {
2924 (*lnLDiv) += modelSettings[0].logLikelihoodsAll[site] * nSitesOfPat[c];
2925 site++;
2926 }
2927 }
2928 }
2929 else
2930 {
2931 /* has invariable category */
2932 pInvar = *(GetParamVals (m->pInvar, chain, state[chain]));
2933 clInvar = m->invCondLikes;
2934
2935 /* find base frequencies */
2936 bs = GetParamSubVals (m->stateFreq, chain, state[chain]);
2937
2938 /* if covarion model, adjust base frequencies */
2939 if (m->switchRates != NULL)
2940 {
2941 /* find the stationary frequencies */
2942 swr = GetParamVals(m->switchRates, chain, state[chain]);
2943 s01 = swr[0];
2944 s10 = swr[1];
2945 probOn = s01 / (s01 + s10);
2946 probOff = 1.0 - probOn;
2947
2948 /* now adjust the base frequencies; on-state stored first in cond likes */
2949 for (j=0; j<nStates/2; j++)
2950 {
2951 covBF[j] = bs[j] * probOn;
2952 covBF[j+nStates/2] = bs[j] * probOff;
2953 }
2954
2955 /* finally set bs pointer to adjusted values */
2956 bs = covBF;
2957 }
2958
2959 for (c=0; c<m->numChars; c++)
2960 {
2961 likeI = 0.0;
2962 for (j=0; j<nStates; j++)
2963 likeI += (*(clInvar++)) * bs[j];
2964 if (likeI != 0.0)
2965 {
2966 lnLikeI = log(likeI * pInvar);
2967 diff = lnLikeI - modelSettings[0].logLikelihoodsAll[site];
2968 }
2969 else
2970 diff = -1000.0;
2971 if (diff < -200.0)
2972 (*lnLDiv) += modelSettings[0].logLikelihoodsAll[site] * nSitesOfPat[c];
2973 else if (diff > 200.0)
2974 (*lnLDiv) += lnLikeI * nSitesOfPat[c];
2975 else
2976 {
2977 (*lnLDiv) += (modelSettings[0].logLikelihoodsAll[site] + log(1.0 + exp(diff))) * nSitesOfPat[c];
2978 }
2979 site++;
2980 }
2981 }
2982 (*lnL) += m->lnLike[2*chain + state[chain]];
2983 }
2984 }
2985 /* check for numerical errors */
2986 assert ((*lnL) == (*lnL));
2987 }
2988
2989 return beagleReturn;
2990 }
2991
2992 #endif /* BEAGLE_V3_ENABLED */
2993
2994 #endif /* BEAGLE_ENABLED */
2995
2996
BeagleNotLinked()2997 void BeagleNotLinked()
2998 {
2999 MrBayesPrint ("%s BEAGLE library is not linked to this executable.\n", spacer);
3000 }
3001
3002
BeagleThreadsNotAvailable()3003 void BeagleThreadsNotAvailable()
3004 {
3005 MrBayesPrint ("%s BEAGLE CPU threading requires v3.1 and higher of the library.\n", spacer);
3006 }
3007