1 /*
2 * Best 2.2
3 *
4 * This file contains the functions
5 * for calculating the probability of
6 * gene trees given the species tree
7 * and the prior probability of the
8 * species tree
9 *
10 * Liang Liu
11 * Department of Statistics
12 * The Ohio State University
13 * Columbus, Ohio
14 *
15 * liuliang@stat.ohio-state.edu
16 */
17
18 #include "bayes.h"
19 #include "best.h"
20 #include "command.h"
21 #include "mcmc.h"
22 #include "model.h"
23 #include "proposal.h"
24 #include "utils.h"
25
26 /****************************** Local functions converted by Fredrik from BEST code *****************************/
27 int CompareDepths (const void *x, const void *y);
28 int CompareDoubles (const void *x, const void *y);
29 int CompareNodes (const void *x, const void *y);
30 int CompareNodesByX (const void *x, const void *y);
31 int GetSpeciesTreeFromMinDepths (Tree* speciesTree, double *depthMatrix);
32 int GetDepthMatrix(Tree * speciesTree, double *depthMatrix);
33 int GetMeanDist (Tree *speciesTree, double *depthMatrix, double *mean);
34 int GetMinDepthMatrix (Tree **geneTrees, int numGeneTrees, double *depthMatrix);
35 void LineagesIn (TreeNode* geneTreeNode, TreeNode* speciesTreeNode);
36 double LnPriorProbGeneTree (Tree *geneTree, double mu, Tree *speciesTree, double *popSizePtr);
37 double LnProposalProbSpeciesTree (Tree *speciesTree, double *depthMatrix, double expRate);
38 void MapGeneTreeToSpeciesTree (Tree *geneTree, Tree *speciesTree);
39 int ModifyDepthMatrix (double expRate, double *depthMatrix, RandLong *seed);
40
41 /* Global BEST variables */
42 BitsLong **speciesPairSets;
43 double *depthMatrix;
44
45 /* Allocate variables used by best code during mcmc */
AllocateBestChainVariables(void)46 void AllocateBestChainVariables (void)
47 {
48 int i, j, index, numUpperTriang, nLongsNeeded;
49
50 // Free if by mistake variables are already allocated
51 if (memAllocs[ALLOC_BEST] == YES)
52 FreeBestChainVariables ();
53
54 // Allocate space for upper triangular pair sets
55 numUpperTriang = (numSpecies * (numSpecies-1)) / 2;
56 nLongsNeeded = ((numSpecies - 1) / nBitsInALong) + 1;
57 speciesPairSets = (BitsLong **) SafeCalloc (numUpperTriang, sizeof(BitsLong *));
58 speciesPairSets[0] = (BitsLong *) SafeCalloc (numUpperTriang*nLongsNeeded, sizeof(BitsLong));
59 for (i=1; i<numUpperTriang; i++)
60 speciesPairSets[i] = speciesPairSets[0] + i*nLongsNeeded;
61
62 // Set upper triangular pair partitions once and for all
63 index = 0;
64 for (i=0; i<numSpecies; i++) {
65 for (j=i+1; j<numSpecies; j++) {
66 SetBit(i, speciesPairSets[index]);
67 SetBit(j, speciesPairSets[index]);
68 index++;
69 }
70 }
71
72 /* allocate species for depthMatrix */
73 depthMatrix = SafeCalloc (numUpperTriang, sizeof(double));
74
75 memAllocs[ALLOC_BEST] = YES;
76 }
77
78
79 /** Compare function (Depth struct) for qsort */
CompareDepths(const void * x,const void * y)80 int CompareDepths (const void *x, const void *y) {
81
82 if ((*((Depth *)(x))).depth < (*((Depth *)(y))).depth)
83 return -1;
84 else if ((*((Depth *)(x))).depth > (*((Depth *)(y))).depth)
85 return 1;
86 else
87 return 0;
88 }
89
90
91 /** Compare function (doubles) for qsort */
CompareDoubles(const void * x,const void * y)92 int CompareDoubles (const void *x, const void *y) {
93
94 if (*((double *)(x)) < *((double *)(y)))
95 return -1;
96 else if (*((double *)(x)) > *((double *)(y)))
97 return 1;
98 else
99 return 0;
100 }
101
102
103 /** Compare function (TreeNode struct) for qsort */
CompareNodes(const void * x,const void * y)104 int CompareNodes (const void *x, const void *y) {
105
106 if ((*((TreeNode **)(x)))->nodeDepth < (*((TreeNode**)(y)))->nodeDepth)
107 return -1;
108 else if ((*((TreeNode **)(x)))->nodeDepth > (*((TreeNode**)(y)))->nodeDepth)
109 return 1;
110 else
111 return 0;
112 }
113
114
115 /** Compare function (TreeNode struct; sort by x, then by nodeDepth) for qsort */
CompareNodesByX(const void * x,const void * y)116 int CompareNodesByX (const void *x, const void *y) {
117
118 if ((*((TreeNode **)(x)))->x < (*((TreeNode**)(y)))->x)
119 return -1;
120 else if ((*((TreeNode **)(x)))->x > (*((TreeNode**)(y)))->x)
121 return 1;
122 else {
123 if ((*((TreeNode **)(x)))->nodeDepth < (*((TreeNode**)(y)))->nodeDepth)
124 return -1;
125 else if ((*((TreeNode **)(x)))->nodeDepth > (*((TreeNode**)(y)))->nodeDepth)
126 return 1;
127 else
128 return 0;
129 }
130 }
131
132
133 /**-----------------------------------------------------------------
134 |
135 | FillSpeciesTreeParams: Fill in species trees (start value)
136 |
137 ------------------------------------------------------------------*/
FillSpeciesTreeParams(RandLong * seed,int fromChain,int toChain)138 int FillSpeciesTreeParams (RandLong *seed, int fromChain, int toChain)
139 {
140 int i, k, chn, numGeneTrees, freeBestChainVars;
141 Param *p;
142 Tree *speciesTree, **geneTrees;
143
144 // Allocate space for global best model variables used in this function, in case they are not allocated
145 if (memAllocs[ALLOC_BEST] == NO)
146 {
147 freeBestChainVars = YES;
148 AllocateBestChainVariables();
149 }
150 else
151 freeBestChainVars = NO;
152
153 // Use global variable numTopologies to calculate number of gene trees
154 // There is one topology for the species tree, the other ones are gene trees
155 // The number of current divisions is not safe because one gene tree can have
156 // several partitions, for instance if we assign the different genes on the
157 // mitochondrion different substitution models, or if we assign different rates
158 // to the codon site positions in a sequence
159 numGeneTrees = numTopologies - 1;
160 geneTrees = (Tree **) SafeCalloc (numGeneTrees, sizeof(Tree*));
161
162 // Build species trees for state 0
163 for (chn=fromChain; chn<toChain; chn++)
164 {
165 for (k=0; k<numParams; k++)
166 {
167 p = ¶ms[k];
168 if (p->paramType == P_SPECIESTREE)
169 {
170 // Find species tree and gene trees
171 speciesTree = GetTree(p, chn, 0);
172 for (i=0; i<p->nSubParams; i++)
173 geneTrees[i] = GetTree(p->subParams[i], chn, 0);
174
175 // Get minimum depth matrix for species tree
176 GetMinDepthMatrix (geneTrees, numGeneTrees, depthMatrix);
177
178 // Get a species tree from min depth matrix
179 GetSpeciesTreeFromMinDepths(speciesTree, depthMatrix);
180
181 assert (IsSpeciesTreeConsistent(speciesTree, chn) == YES);
182
183 // Label the tips
184 if (LabelTree (speciesTree, speciesNameSets[speciespartitionNum].names) == ERROR)
185 {
186 FreeBestChainVariables();
187 return (ERROR);
188 }
189 }
190 }
191 }
192
193 // Free gene trees
194 free (geneTrees);
195
196 // Free best model variables if appropriate
197 if (freeBestChainVars == YES)
198 FreeBestChainVariables();
199
200 return (NO_ERROR);
201 }
202
203
204 /**-----------------------------------------------------------------
205 |
206 | FreeBestChainVariables: Free best variables used during an mcmc
207 | run.
208 |
209 ------------------------------------------------------------------*/
FreeBestChainVariables(void)210 void FreeBestChainVariables(void)
211 {
212 if (memAllocs[ALLOC_BEST] == YES) {
213 free (speciesPairSets[0]);
214 free (speciesPairSets);
215 speciesPairSets = NULL;
216 }
217
218 free (depthMatrix);
219 depthMatrix = NULL;
220
221 memAllocs[ALLOC_BEST] = NO;
222 }
223
224
225 /**---------------------------------------------------------------------
226 |
227 | GetDepthMatrix:
228 |
229 | This algorithm calculates the upper triangular depth matrix for the
230 | species tree. Time complexity O(n^2).
231 |
232 | @param speciesTree The species tree (in)
233 | @param depthMatrix The minimum depth matrix, upper triangular array (out)
234 | @returns Returns ERROR or NO_ERROR
235 ----------------------------------------------------------------------*/
GetDepthMatrix(Tree * speciesTree,double * depthMatrix)236 int GetDepthMatrix (Tree *speciesTree, double *depthMatrix) {
237
238 int i, left, right, numUpperTriang, index, nLongsNeeded, freeBitsets;
239 double maxDepth;
240 TreeNode *p;
241
242 // Make sure we have bitfields allocated and set
243 if (speciesTree->bitsets == NULL)
244 {
245 AllocateTreePartitions(speciesTree);
246 freeBitsets = YES;
247 }
248 else
249 {
250 ResetTreePartitions(speciesTree); // just in case
251 freeBitsets = NO;
252 }
253
254 // Calculate number of values in the upper triangular matrix
255 numUpperTriang = numSpecies * (numSpecies - 1) / 2;
256
257 // Number of longs needed in a bitfield representing a species set
258 nLongsNeeded = ((numSpecies -1) / nBitsInALong) + 1;
259
260 // Set all cells to max
261 maxDepth = speciesTree->root->left->nodeDepth; // root depth
262 for (i=0; i<numUpperTriang; i++)
263 depthMatrix[i] = maxDepth;
264
265 // Loop over interior nodes
266 for (i=0; i<speciesTree->nIntNodes; i++)
267 {
268 p = speciesTree->intDownPass[i];
269 for (left = FirstTaxonInPartition(p->left->partition, nLongsNeeded); left < numSpecies; left = NextTaxonInPartition(left, p->left->partition, nLongsNeeded))
270 {
271 for (right = FirstTaxonInPartition(p->right->partition, nLongsNeeded); right < numSpecies; right = NextTaxonInPartition(right, p->right->partition, nLongsNeeded))
272 {
273 index = UpperTriangIndex(left, right, numSpecies);
274 depthMatrix[index] = p->nodeDepth;
275 }
276 }
277 }
278
279 // Free partitions if appropriate
280 if (freeBitsets == YES)
281 FreeTreePartitions(speciesTree);
282
283 return (NO_ERROR);
284 }
285
286
287 /**---------------------------------------------------------------------
288 |
289 | GetMeanDist
290 |
291 | This algorithm calculates the mean distance between a distance matrix
292 | and the minimum depths that define a species tree.
293 |
294 | @param speciesTree The species tree (in)
295 | @param minDepthMatrix The minimum depth matrix, upper triangular array (in)
296 | @param mean The mean distance (out)
297 | @returns Returns ERROR or NO_ERROR
298 ----------------------------------------------------------------------*/
GetMeanDist(Tree * speciesTree,double * minDepthMatrix,double * mean)299 int GetMeanDist (Tree *speciesTree, double *minDepthMatrix, double *mean) {
300
301 int i, left, right, index, nLongsNeeded, freeBitsets;
302 double dist, minDist=0.0, distSum;
303 TreeNode *p;
304
305 // Make sure we have bitfields allocated and set
306 if (speciesTree->bitsets == NULL)
307 {
308 AllocateTreePartitions(speciesTree);
309 freeBitsets = YES;
310 }
311 else
312 {
313 ResetTreePartitions(speciesTree); // just in case
314 freeBitsets = NO;
315 }
316
317 // Number of longs needed in a bitfield representing a species set
318 nLongsNeeded = ((numSpecies -1) / nBitsInALong) + 1;
319
320 // Loop over interior nodes
321 distSum = 0.0;
322 for (i=0; i<speciesTree->nIntNodes; i++)
323 {
324 p = speciesTree->intDownPass[i];
325 p->x = 0;
326 while ((left=FirstTaxonInPartition(p->left->partition, nLongsNeeded)) < numSpecies)
327 {
328 while ((right=FirstTaxonInPartition(p->right->partition, nLongsNeeded)) < numSpecies)
329 {
330 p->x++;
331 index = UpperTriangIndex(left, right, numSpecies);
332 dist = depthMatrix[index] - p->nodeDepth;
333 if (p->x == 1)
334 minDist = dist;
335 else if (dist < minDist)
336 minDist = dist;
337 ClearBit(right, p->right->partition);
338 }
339 ClearBit(left, p->left->partition);
340 }
341 distSum += minDist;
342 }
343
344 (*mean) = distSum / speciesTree->nIntNodes;
345
346 // Reset partitions that were destroyed above or free partitions, as appropriate
347 if (freeBitsets == YES)
348 FreeTreePartitions(speciesTree);
349 else
350 ResetTreePartitions(speciesTree);
351
352 return (NO_ERROR);
353 }
354
355
356 /**---------------------------------------------------------------------
357 |
358 | GetMinDepthMatrix: converted from GetMinDists.
359 |
360 | This algorithm scans the gene trees and calculates the minimum depth
361 | (height) separating species across gene trees. The complexity of the
362 | original algorithm was O(mn^3), where m is the number of gene trees and
363 | n is the number of taxa in each gene tree. I think this algorithm has
364 | complexity that is better on average, but the difference is small.
365 |
366 | I have rewritten the algorithm also to show alternative techniques that
367 | could be used in this and other BEST algorithms.
368 |
369 | @param geneTrees The gene trees (in)
370 | @param depthMatrix The minimum depth matrix, upper triangular array (out)
371 | @returns Returns ERROR or NO_ERROR
372 ----------------------------------------------------------------------*/
GetMinDepthMatrix(Tree ** geneTrees,int numGeneTrees,double * depthMatrix)373 int GetMinDepthMatrix (Tree **geneTrees, int numGeneTrees, double *depthMatrix) {
374
375 int i, j, w, nLongsNeeded, numUpperTriang, index, trace=0;
376 double maxDepth;
377 TreeNode *p;
378 BitsLong **speciesSets;
379
380 // Allocate space for species partitions
381 nLongsNeeded = ((numSpecies -1) / nBitsInALong) + 1; // number of longs needed in a bitfield representing a species set
382 speciesSets = (BitsLong **) SafeCalloc ((2*(size_t)numLocalTaxa-1), sizeof(BitsLong *));
383 speciesSets[0] = (BitsLong *) SafeCalloc ((2*(size_t)numLocalTaxa-1)*nLongsNeeded, sizeof(BitsLong));
384 for (i=1; i<2*numLocalTaxa-1; i++)
385 speciesSets[i] = speciesSets[0] + i*nLongsNeeded;
386
387 // Set tip species partitions once and for all
388 for (i=0; i<numLocalTaxa; i++)
389 SetBit(speciespartitionId[i][speciespartitionNum]-1, speciesSets[i]);
390
391 // Set initial max depth for upper triangular matrix
392 numUpperTriang = (numSpecies * (numSpecies - 1)) / 2;
393 maxDepth = geneTrees[0]->root->left->nodeDepth;
394 for (i=0; i<numUpperTriang; i++)
395 depthMatrix[i] = maxDepth;
396
397 // Now we are ready to cycle over gene trees
398 for (w=0; w<numGeneTrees; w++)
399 {
400 if (trace) {
401 printf ("\nGene %d\n",w);
402 ShowTree(geneTrees[w]);
403 }
404
405 // Set species sets for interior nodes. O(n)
406 for (i=0; i<geneTrees[w]->nIntNodes; i++) {
407 p = geneTrees[w]->intDownPass[i];
408 for (j=0; j<nLongsNeeded; j++)
409 speciesSets[p->index][j] = speciesSets[p->left->index][j] | speciesSets[p->right->index][j];
410 }
411
412 // Now order the interior nodes in terms of node depth. We rely on the fact that the
413 // ordered sequence is a valid downpass sequence. O(log n).
414 qsort((void *)(geneTrees[w]->intDownPass), (size_t)(geneTrees[w]->nIntNodes), sizeof(TreeNode *), CompareNodes);
415
416 // Finally find the minimum for each cell in the upper triangular matrix
417 // This is the time critical step with complexity O(n^3) in the simplest
418 // algorithm version. This algorithm should do a little better in most cases.
419 for (i=0; i<numUpperTriang; i++) {
420
421 // Find shallowest node that has the pair
422 for (j=0; j<geneTrees[w]->nIntNodes; j++) {
423 p = geneTrees[w]->intDownPass[j];
424
425 // Because nodes are ordered in time, if this test is true then we cannot beat the minimum
426 if (p->nodeDepth > depthMatrix[i])
427 break;
428
429 // Check whether the node is a candidate minimum for the species pair
430 // If the test is true, we know from the test above that p->nodeDepth is
431 // either a tie or the new minimum
432 if (IsPartNested(speciesPairSets[i], speciesSets[p->index], nLongsNeeded) == YES) {
433 depthMatrix[i] = p->nodeDepth;
434 break;
435 }
436 }
437 }
438 } // Next gene tree
439
440 if (trace)
441 {
442 index = 0;
443 printf ("Mindepth matrix\n");
444 for (i=0;i<numSpecies;i++) {
445 for (j=0; j<i; j++)
446 printf (" ");
447 for (j=i+1;j<numSpecies;j++) {
448 printf ("%.6f ",depthMatrix[index]);
449 index++;
450 }
451 printf ("\n");
452 }
453 printf ("\n");
454 }
455
456 free (speciesSets[0]);
457 free (speciesSets);
458
459 return (NO_ERROR);
460 }
461
462
463 /**---------------------------------------------------------------------
464 |
465 | GetSpeciesTreeFromMinDepths: converted from GetConstraints, Startsptree,
466 | and MaximumTree.
467 |
468 | This is a clustering algorithm based on minimum depths for species pairs.
469 | It reduces an n choose 2 upper triangular min depth matrix to an array
470 | of n-1 node depths, which fit onto a tree.
471 |
472 | @param speciesTree The species tree to be filled (out)
473 | @param depthMatrix The min depth matrix, upper triangular array (in)
474 | @returns Returns NO_ERROR if success, ERROR if negative brlens occur
475 ----------------------------------------------------------------------*/
GetSpeciesTreeFromMinDepths(Tree * speciesTree,double * depthMatrix)476 int GetSpeciesTreeFromMinDepths (Tree* speciesTree, double *depthMatrix) {
477
478 int i, j, numUpperTriang, nLongsNeeded, index, nextNodeIndex;
479 Depth *minDepth;
480 PolyTree *polyTree;
481 PolyNode *p, *q, *r, *u, *qPrev, *rPrev;
482
483 nLongsNeeded = ((numSpecies - 1) / nBitsInALong) + 1;
484 numUpperTriang = numSpecies*(numSpecies - 1) / 2;
485 minDepth = (Depth *) SafeCalloc (numUpperTriang, sizeof(Depth));
486
487 // Convert depthMatrix to an array of Depth structs
488 index = 0;
489 for (i=0; i<numSpecies; i++) {
490 for (j=i+1; j<numSpecies; j++) {
491 minDepth[index].depth = depthMatrix[index];
492 minDepth[index].pairSet = speciesPairSets[index];
493 index++;
494 }
495 }
496
497 // Sort the array of distance structs (O(log n^2))
498 qsort((void *)(minDepth), (size_t)(numUpperTriang), sizeof(Depth), CompareDepths);
499
500 // The algorithm below reduces the upper triangular matrix (n choose 2) to an n-1
501 // array in O(n^2log(n)) time. We build the tree at the same time, since we can
502 // find included pairs in the tree in log(n) time. We use a polytomous tree for this.
503
504 // Allocate space for polytomous tree and set up partitions
505 polyTree = AllocatePolyTree(numSpecies);
506 AllocatePolyTreePartitions(polyTree);
507
508 // Build initial tree (a bush)
509 polyTree->isRooted = YES;
510 polyTree->isClock = YES;
511 polyTree->root = &polyTree->nodes[2*numSpecies-2];
512 for (i=0; i<numSpecies; i++) {
513 p = &polyTree->nodes[i];
514 p->index = i;
515 p->depth = 0.0;
516 p->left = NULL;
517 if (i<numSpecies-1)
518 p->sib = &polyTree->nodes[i+1];
519 else
520 p->sib = NULL;
521 p->anc = polyTree->root;
522 }
523 p = polyTree->root;
524 p->index = 2*numSpecies - 2;
525 p->left = &polyTree->nodes[0];
526 p->sib = NULL;
527 p->anc = NULL;
528 p->depth = -1.0;
529 polyTree->nNodes = numSpecies + 1;
530 polyTree->nIntNodes = 1;
531 GetPolyDownPass(polyTree);
532 ResetPolyTreePartitions(polyTree); /* set bitsets (partitions) for initial tree */
533
534 // Resolve bush using sorted depth structs
535 nextNodeIndex = numSpecies;
536 for (i=0; i<numUpperTriang; i++) {
537
538 // Find tip corresponding to first taxon in pair
539 p = &polyTree->nodes[FirstTaxonInPartition(minDepth[i].pairSet, nLongsNeeded)];
540
541 // Descend tree until we find a node within which the pair set is nested
542 do {
543 p = p->anc;
544 }
545 while (!IsPartNested(minDepth[i].pairSet, p->partition, nLongsNeeded));
546
547 if (p->left->sib->sib != NULL) {
548 // This node is still a polytomy
549
550 // Find left and right descendants of new node
551 qPrev = NULL;
552 for (q=p->left; IsSectionEmpty(q->partition, minDepth[i].pairSet, nLongsNeeded); q=q->sib)
553 qPrev = q;
554 rPrev = q;
555 for (r=q->sib; IsSectionEmpty(r->partition, minDepth[i].pairSet, nLongsNeeded); r=r->sib)
556 rPrev = r;
557
558 // Introduce the new node
559 u = &polyTree->nodes[nextNodeIndex];
560 u->index = nextNodeIndex;
561 nextNodeIndex++;
562 polyTree->nIntNodes++;
563 polyTree->nNodes++;
564 u->left = q;
565 u->anc = p;
566 if (p->left == q)
567 p->left = u;
568 else
569 qPrev->sib = u;
570 // former upstream sibling to r should point to r->sib
571 if (rPrev == q)
572 u->sib = r->sib;
573 else
574 rPrev->sib = r->sib;
575 if (q->sib == r)
576 u->sib = r->sib;
577 else
578 u->sib = q->sib;
579 u->depth = minDepth[i].depth; // because minDepth structs are sorted, we know this is the min depth
580 assert (u->depth > 0.0);
581
582 // Create new taxon set with bitfield operations
583 for (j=0; j<nLongsNeeded; j++)
584 u->partition[j] = q->partition[j] | r->partition[j];
585
586 // Patch the tree together with the new node added
587 q->sib = r;
588 r->sib = NULL;
589 q->anc = u;
590 r->anc = u;
591 }
592 else if (p == polyTree->root && p->depth < 0.0) {
593 // This is the first time we hit the root of the tree && it is resolved
594 p->depth = minDepth[i].depth;
595 assert (p->depth > 0.0);
596 }
597 // other cases should not be added to tree
598 }
599
600 // Make sure we have a complete species tree
601 assert (polyTree->nIntNodes == numSpecies - 1);
602
603 // Set traversal sequences
604 GetPolyDownPass(polyTree);
605
606 // Set branch lengths from node depths (not done automatically for us)
607 // Make sure all branch lengths are nonnegative (we can have 0.0 brlens, they
608 // should not be problematic in a species tree; they occur when there are
609 // ties in the min depth matrix that have not been modified by the move)
610 for (i=0; i<polyTree->nNodes; i++) {
611 p = polyTree->allDownPass[i];
612 if (p->anc == NULL)
613 p->length = 0.0;
614 else
615 p->length = p->anc->depth - p->depth;
616 if (p->length < 0.0) {
617 FreePolyTree(polyTree);
618 free (minDepth);
619 return (ERROR);
620 }
621 }
622
623 // Copy to species tree from polytomous tree
624 CopyToSpeciesTreeFromPolyTree (speciesTree, polyTree);
625
626 // Free locally allocated variables
627 FreePolyTree(polyTree);
628 free (minDepth);
629
630 return(NO_ERROR);
631 }
632
633
634 /**---------------------------------------------------------------------------------------
635 |
636 | IsSpeciesTreeConsistent: Called when user tries to set a species tree or when
637 | attempting to use a species tree from a check point as starting value.
638 |
639 -----------------------------------------------------------------------------------------*/
IsSpeciesTreeConsistent(Tree * speciesTree,int chain)640 int IsSpeciesTreeConsistent (Tree *speciesTree, int chain)
641 {
642 int i, answer, numGeneTrees, numUpperTriang, freeBestVars;
643 double *speciesTreeDepthMatrix;
644 Tree **geneTrees;
645
646 freeBestVars = NO;
647 if (memAllocs[ALLOC_BEST] == NO)
648 {
649 AllocateBestChainVariables();
650 freeBestVars = YES;
651 }
652
653 numGeneTrees = numTrees - 1;
654 geneTrees = (Tree **) SafeCalloc (numGeneTrees, sizeof(Tree *));
655 for (i=0; i<numTrees-1; i++)
656 geneTrees[i] = GetTreeFromIndex(i, chain, state[chain]);
657
658 numUpperTriang = numSpecies * (numSpecies - 1) / 2;
659 speciesTreeDepthMatrix = (double *) SafeCalloc (numUpperTriang, sizeof(double));
660
661 GetMinDepthMatrix(geneTrees, numGeneTrees, depthMatrix);
662 GetDepthMatrix(speciesTree, speciesTreeDepthMatrix);
663
664 for (i=0; i<numUpperTriang; i++)
665 {
666 if (depthMatrix[i] < speciesTreeDepthMatrix[i])
667 break;
668 }
669
670 if (i == numUpperTriang)
671 answer = YES;
672 else
673 answer = NO;
674
675 if (answer == NO)
676 ShowNodes(speciesTree->root, 0, YES);
677
678 if (freeBestVars == YES)
679 FreeBestChainVariables();
680
681 free (speciesTreeDepthMatrix);
682 free (geneTrees);
683
684 return answer;
685 }
686
687
688 /**---------------------------------------------------------------------------------------
689 |
690 | LineagesIn: Recursive function to get number of gene tree lineages coming into each
691 | branch of the species tree (in ->x of speciestree nodes). We also assign each gene
692 | tree lineage to the corresponding species tree lineage (in ->x of genetree nodes).
693 | Finally, number of coalescent events is recorded (in ->y of speciestree nodes).
694 | Time complexity is O(n).
695 |
696 -----------------------------------------------------------------------------------------*/
LineagesIn(TreeNode * geneTreeNode,TreeNode * speciesTreeNode)697 void LineagesIn (TreeNode *geneTreeNode, TreeNode *speciesTreeNode)
698 {
699 int nLongsNeeded;
700
701 if (geneTreeNode->nodeDepth < speciesTreeNode->nodeDepth) {
702 // climb up species tree
703 if (speciesTreeNode->left == NULL) {
704 assert (geneTreeNode->left == NULL);
705 speciesTreeNode->x++;
706 }
707 else {
708 nLongsNeeded = (numSpecies - 1) / nBitsInALong + 1;
709 speciesTreeNode->x++;
710 if (IsPartNested(geneTreeNode->partition, speciesTreeNode->left->partition, nLongsNeeded) == YES)
711 LineagesIn (geneTreeNode, speciesTreeNode->left);
712 else if (IsPartNested(geneTreeNode->partition, speciesTreeNode->right->partition, nLongsNeeded) == YES)
713 LineagesIn (geneTreeNode, speciesTreeNode->right);
714 }
715 }
716 else {
717 // climb up gene tree
718 if (geneTreeNode->left != NULL)
719 LineagesIn(geneTreeNode->left, speciesTreeNode);
720 if (geneTreeNode->right != NULL)
721 LineagesIn(geneTreeNode->right, speciesTreeNode);
722 if (geneTreeNode->left == NULL) {
723 speciesTreeNode->x++;
724 assert (speciesTreeNode->left == NULL);
725 }
726 else {
727 speciesTreeNode->y++;
728 }
729 geneTreeNode->x = speciesTreeNode->index;
730 }
731 }
732
733
734 /**-----------------------------------------------------------------
735 |
736 | LnSpeciesTreeProb: Wrapper for LnJointGeneTreeSpeciesTreePr to
737 | free calling functions from retrieving gene and species trees.
738 |
739 ------------------------------------------------------------------*/
LnSpeciesTreeProb(int chain)740 double LnSpeciesTreeProb(int chain)
741 {
742 int i, numGeneTrees;
743 double lnProb;
744 Tree **geneTrees, *speciesTree;
745 ModelInfo *m;
746
747 m = &modelSettings[0];
748
749 speciesTree = GetTree(m->speciesTree, chain, state[chain]);
750
751 numGeneTrees = m->speciesTree->nSubParams;
752 geneTrees = (Tree **) SafeCalloc (numGeneTrees, sizeof(Tree *));
753
754 for (i=0; i<m->speciesTree->nSubParams; i++)
755 geneTrees[i] = GetTree(m->speciesTree->subParams[i], chain, state[chain]);
756
757 lnProb = LnJointGeneTreeSpeciesTreePr(geneTrees, numGeneTrees, speciesTree, chain);
758
759 free (geneTrees);
760
761 return lnProb;
762 }
763
764
765 /**-----------------------------------------------------------------
766 |
767 | LnJointGeneTreeSpeciesTreePr: Converted from LnJointGenetreePr,
768 | SPLogLike, SPLogPrior.
769 |
770 | In this function we calculate the entire probability of the species
771 | tree, including its probability given its priors, and the probability
772 | of the gene trees given the species tree.
773 |
774 ------------------------------------------------------------------*/
LnJointGeneTreeSpeciesTreePr(Tree ** geneTrees,int numGeneTrees,Tree * speciesTree,int chain)775 double LnJointGeneTreeSpeciesTreePr(Tree **geneTrees, int numGeneTrees, Tree *speciesTree, int chain)
776 {
777 double lnPrior, lnLike, clockRate, mu, *popSizePtr, sR, eR, sF;
778 int i;
779 ModelInfo *m;
780 ModelParams *mp;
781
782 // Get model info for species tree
783 m = &modelSettings[speciesTree->relParts[0]];
784
785 // Get model params for species tree
786 mp = &modelParams[speciesTree->relParts[0]];
787
788 // Get popSize ptr
789 popSizePtr = GetParamVals(m->popSize, chain, state[chain]);
790
791 // Get clock rate
792 if (speciesTree->isCalibrated == YES)
793 clockRate = *GetParamVals(m->clockRate, chain, state[chain]);
794 else
795 clockRate = 1.0;
796
797 // Calculate probability of gene trees given species tree
798 // ShowNodes(speciesTree->root, 0, YES);
799 lnLike = 0.0;
800 mu = clockRate;
801 for (i=0; i<numGeneTrees; i++) {
802 lnLike += LnPriorProbGeneTree(geneTrees[i], mu, speciesTree, popSizePtr);
803 }
804
805 // Calculate probability of species tree given its priors
806 if (strcmp(mp->speciesTreeBrlensPr, "Birthdeath") == 0) {
807 sR = *GetParamVals(m->speciationRates, chain, state[chain]);
808 eR = *GetParamVals(m->extinctionRates, chain, state[chain]);
809 sF = mp->sampleProb;
810 lnPrior = 0.0;
811 LnBirthDeathPriorPr(speciesTree, clockRate, &lnPrior, sR, eR, mp->sampleStrat, sF);
812 }
813 else
814 lnPrior = 0.0;
815
816 // The population size is taken care of elsewhere
817
818 return lnLike + lnPrior;
819 }
820
821
822 /**-----------------------------------------------------------------
823 |
824 | LnPriorProbGeneTree: Calculate the prior probability of a gene
825 | tree.
826 |
827 ------------------------------------------------------------------*/
LnPriorProbGeneTree(Tree * geneTree,double mu,Tree * speciesTree,double * popSizePtr)828 double LnPriorProbGeneTree (Tree *geneTree, double mu, Tree *speciesTree, double *popSizePtr)
829 {
830 int i, k, index, nEvents, trace=0;
831 double N, lnProb, ploidyFactor, theta, timeInterval;
832 TreeNode *p, *q=NULL, *r;
833 ModelParams *mp;
834
835 // Get model params
836 mp = &modelParams[speciesTree->relParts[0]];
837
838 // Find ploidy setting
839 if (strcmp(mp->ploidy, "Diploid") == 0)
840 ploidyFactor = 4.0;
841 else if (strcmp(mp->ploidy, "Haploid") == 0)
842 ploidyFactor = 2.0;
843 else /* if (strcmp(mp->ploidy, "Zlinked") == 0) */
844 ploidyFactor = 3.0;
845
846 // Initialize species tree with theta in d
847 for (i=0; i<speciesTree->nNodes-1; i++) {
848 p = speciesTree->allDownPass[i];
849 if (strcmp(mp->popVarPr, "Equal") != 0)
850 N = popSizePtr[p->index];
851 else
852 N = popSizePtr[0];
853 p->d = ploidyFactor * N * mu;
854 }
855
856 // Map gene tree to species tree
857 MapGeneTreeToSpeciesTree(geneTree, speciesTree);
858
859 // Sort gene tree interior nodes first by speciestree branch on which they coalesce, then in time order
860 qsort((void *)(geneTree->intDownPass), (size_t)(geneTree->nIntNodes), sizeof(TreeNode *), CompareNodesByX);
861
862 // Debug output of qsort result
863 if (trace) {
864 printf ("index -- x -- nodeDepth for gene tree\n");
865 for (i=0; i<geneTree->nIntNodes; i++)
866 printf ("%d -- %d -- %e\n", geneTree->intDownPass[i]->index, geneTree->intDownPass[i]->x, geneTree->intDownPass[i]->nodeDepth);
867 }
868
869 // Now calculate probability after making sure species tree nodes appear in index order
870 // (the order does not have to be a correct downpass sequence)
871 for (i=0; i<speciesTree->memNodes; i++)
872 {
873 p = &(speciesTree->nodes[i]);
874 speciesTree->allDownPass[p->index] = p;
875 }
876 index = 0;
877 lnProb = 0.0;
878 for (i=0; i<speciesTree->nNodes-1; i++) {
879
880 p = speciesTree->allDownPass[i];
881
882 // Get theta
883 theta = p->d;
884
885 // Get number of events
886 nEvents = p->y;
887
888 // Calculate probability
889 lnProb += nEvents * log (2.0 / theta);
890
891 for (k=p->x; k > p->x - p->y; k--) {
892
893 q = geneTree->intDownPass[index];
894 assert (q->x == p->index);
895
896 if (k == p->x)
897 timeInterval = q->nodeDepth - p->nodeDepth;
898 else {
899 r = geneTree->intDownPass[index-1];
900 timeInterval = q->nodeDepth - r->nodeDepth;
901 }
902
903 lnProb -= (k * (k - 1) * timeInterval) / theta;
904 index++;
905 }
906
907 if (p->x - p->y > 1) {
908
909 if (nEvents == 0)
910 timeInterval = p->anc->nodeDepth - p->nodeDepth;
911 else
912 {
913 /* FIXME: q == NULL if above loop not run (from clang static analyzer) */
914 timeInterval = p->anc->nodeDepth - q->nodeDepth;
915 }
916
917 assert (p->anc->anc != NULL);
918 assert (timeInterval >= 0.0);
919
920 k = p->x - p->y;
921 lnProb -= (k * (k - 1) * timeInterval) / theta;
922 }
923 }
924
925 // Restore downpass sequences (probably not necessary for gene tree, but may be if some
926 // code relies on intDownPass and allDownPass to be in same order)
927 GetDownPass(speciesTree);
928 GetDownPass(geneTree);
929
930 // Free space
931 FreeTreePartitions(speciesTree);
932 FreeTreePartitions(geneTree);
933
934 return lnProb;
935 }
936
937
938 /**---------------------------------------------------------------------
939 |
940 | LnProposalProbSpeciesTree:
941 |
942 | This algorithm calculates the probability of proposing a particular
943 | species tree given a distance matrix modified using a scheme based on
944 | truncated exponential distributions with rate expRate.
945 |
946 | @param speciesTree The species tree (in)
947 | @param depthMatrix The minimum depth matrix, upper triangular array (in)
948 | @param expRate Rate of truncated exponential distribution
949 | @returns Returns probability of proposing the species tree
950 ----------------------------------------------------------------------*/
LnProposalProbSpeciesTree(Tree * speciesTree,double * depthMatrix,double expRate)951 double LnProposalProbSpeciesTree (Tree *speciesTree, double *depthMatrix, double expRate)
952 {
953 int i, left, right, index, nLongsNeeded, freeBitsets;
954 double dist, normConst=1.0, negLambdaX=0.0, eNegLambdaX, density, prob,
955 sumDensRatio, prodProb, lnProb;
956 TreeNode *p;
957
958 // Make sure we have bitfields allocated and set
959 if (speciesTree->bitsets == NULL)
960 freeBitsets = YES;
961 else
962 freeBitsets = NO;
963 AllocateTreePartitions(speciesTree);
964
965 // Number of longs needed in a bitfield representing a species set
966 nLongsNeeded = ((numSpecies -1) / nBitsInALong) + 1;
967
968 // Loop over interior nodes
969 lnProb = 0.0;
970 for (i=0; i<speciesTree->nIntNodes; i++)
971 {
972 p = speciesTree->intDownPass[i];
973 p->x = 0;
974 sumDensRatio = 0.0;
975 prodProb = 1.0;
976 for (left = FirstTaxonInPartition(p->left->partition, nLongsNeeded); left < numSpecies; left = NextTaxonInPartition(left, p->left->partition, nLongsNeeded))
977 {
978 for (right = FirstTaxonInPartition(p->right->partition, nLongsNeeded); right < numSpecies; right = NextTaxonInPartition(right, p->right->partition, nLongsNeeded))
979 {
980 p->x++;
981 index = UpperTriangIndex(left, right, numSpecies); assert (index < numSpecies*(numSpecies - 1) / 2);
982 dist = depthMatrix[index] - p->nodeDepth; // distance between depth matrix entry and actual species-tree node
983 normConst = 1.0 - exp(-expRate * depthMatrix[index]); // normalization constant because of truncation of exp distribution
984 negLambdaX = - expRate * dist;
985 eNegLambdaX = exp(negLambdaX);
986 density = expRate * eNegLambdaX / normConst; // density for x == dist, f(dist)
987 prob = (1.0 - eNegLambdaX) / normConst; // cumulative prob for x <= dist, F(dist)
988 sumDensRatio += density / prob; // warning: if dist==0, prob is ZERO!
989 prodProb *= prob;
990 }
991 }
992 if (p->x == 1)
993 lnProb += log(expRate) + negLambdaX - log(normConst);
994 else
995 lnProb += log(sumDensRatio * prodProb);
996 }
997
998 // to avoid lnProposalProb is NaN at initial steps
999 if (lnProb != lnProb) lnProb = 0.0;
1000
1001 // Free partitions if appropriate
1002 if (freeBitsets == YES)
1003 FreeTreePartitions(speciesTree);
1004
1005 return (lnProb);
1006 }
1007
1008
1009 /**-----------------------------------------------------------------
1010 |
1011 | MapGeneTreeToSpeciesTree: Fold gene tree into species tree. We
1012 | are going to use ->x of gene tree to give index of the
1013 | corresponding node in the species tree. ->x in the species
1014 | tree will give the number of lineages into the corresponding
1015 | branch, and ->y will give the number of coalescent events on
1016 | that branch.
1017 |
1018 ------------------------------------------------------------------*/
MapGeneTreeToSpeciesTree(Tree * geneTree,Tree * speciesTree)1019 void MapGeneTreeToSpeciesTree (Tree *geneTree, Tree *speciesTree)
1020 {
1021 int i, j, nLongsNeeded, trace=0;
1022 TreeNode *p;
1023
1024 // Initialize species partitions for both gene tree and species tree
1025 // This will set the partitions to reflect the partitions in the tree itself,
1026 // which is OK for the species tree, but we want the gene tree partitions to
1027 // reflect the species partitions and not the gene partitions, so we need to
1028 // set them here
1029 AllocateTreePartitions(geneTree);
1030 AllocateTreePartitions(speciesTree);
1031 nLongsNeeded = (numSpecies - 1) / nBitsInALong + 1;
1032 for (i=0; i<geneTree->nNodes-1; i++) {
1033 p = geneTree->allDownPass[i];
1034 ClearBits(p->partition, nLongsNeeded);
1035 if (p->left == NULL)
1036 SetBit(speciespartitionId[p->index][speciespartitionNum]-1, p->partition);
1037 else {
1038 for (j=0; j<nLongsNeeded; j++)
1039 p->partition[j] = p->left->partition[j] | p->right->partition[j];
1040 }
1041 }
1042 // Species tree partitions already set by call to AllocateTreePartitions
1043
1044 // Reset ->x and ->y of species tree (->x of gene tree does not need to be initialized)
1045 for (i=0; i<speciesTree->nNodes; i++)
1046 {
1047 p = speciesTree->allDownPass[i];
1048 p->x = 0;
1049 p->y = 0;
1050 }
1051
1052 // Call recursive function to match gene tree and species tree
1053 LineagesIn(geneTree->root->left, speciesTree->root->left);
1054
1055 if (trace) {
1056 printf ("index -- x -- y for species tree\n");
1057 for (i=0; i<speciesTree->nNodes-1; i++)
1058 printf ("%-2d -- %d -- %d\n", speciesTree->allDownPass[i]->index, speciesTree->allDownPass[i]->x, speciesTree->allDownPass[i]->y);
1059 }
1060
1061 if (trace) {
1062 printf ("index -- x -- nodeDepth for gene tree\n");
1063 for (i=0; i<geneTree->nIntNodes; i++)
1064 printf ("%-2d -- %d -- %e\n", geneTree->intDownPass[i]->index, geneTree->intDownPass[i]->x, geneTree->intDownPass[i]->nodeDepth);
1065 }
1066
1067 // Free space
1068 FreeTreePartitions(speciesTree);
1069 FreeTreePartitions(geneTree);
1070 }
1071
1072
1073 /**---------------------------------------------------------------------
1074 |
1075 | ModifyDepthMatrix:
1076 |
1077 | This algorithm uses a truncated exponential distribution to modify
1078 | a depth matrix.
1079 |
1080 | @param expRate The rate of the exponential distribution (in)
1081 | @param depthMatrix The minimum depth matrix to be modified, upper triangular array (in/out)
1082 | @param seed Pointer to seed for random number generator (in/ut)
1083 | @returns Returns ERROR or NO_ERROR
1084 ----------------------------------------------------------------------*/
ModifyDepthMatrix(double expRate,double * depthMatrix,RandLong * seed)1085 int ModifyDepthMatrix (double expRate, double *depthMatrix, RandLong *seed)
1086 {
1087 int i, numUpperTriang;
1088 double u, interval, delta;
1089
1090 numUpperTriang = numSpecies * (numSpecies - 1) / 2;
1091 for (i=0; i<numUpperTriang; i++)
1092 {
1093 interval = depthMatrix[i];
1094 u = RandomNumber(seed);
1095 delta = log (1.0 - u*(1.0 - exp(-expRate*interval))) / (-expRate);
1096 assert (delta <= interval);
1097 depthMatrix[i] -= delta;
1098 }
1099
1100 return (NO_ERROR);
1101 }
1102
1103
1104 /**-----------------------------------------------------------------
1105 |
1106 | Move_GeneTree1: Propose a new gene tree using ExtSPRClock
1107 |
1108 | @param param The parameter (gene tree) to change
1109 | @param chain The chain number
1110 | @param seed Pointer to the seed of the random number gen.
1111 | @param lnPriorRatio Pointer to the log prior ratio (out)
1112 | @param lnProposalRatio Pointer to the log proposal (Hastings) ratio (out)
1113 | @param mvp Pointer to tuning parameter(s)
1114 ------------------------------------------------------------------*/
Move_GeneTree1(Param * param,int chain,RandLong * seed,MrBFlt * lnPriorRatio,MrBFlt * lnProposalRatio,MrBFlt * mvp)1115 int Move_GeneTree1 (Param *param, int chain, RandLong *seed, MrBFlt *lnPriorRatio, MrBFlt *lnProposalRatio, MrBFlt *mvp)
1116 {
1117 int i, numGeneTrees, numUpperTriang;
1118 double newLnProb, oldLnProb, backwardLnProposalProb, forwardLnProposalProb,
1119 *oldMinDepths, *modMinDepths, forwardLambda, backwardLambda, mean;
1120 Tree *newSpeciesTree, *oldSpeciesTree, **geneTrees;
1121 ModelInfo *m;
1122
1123 // Calculate number of gene trees
1124 numGeneTrees = numTopologies - 1;
1125
1126 // Get model settings
1127 m = &modelSettings[param->relParts[0]];
1128
1129 // Get species tree (this trick is possible because we always copy tree params)
1130 newSpeciesTree = GetTree (m->speciesTree, chain, state[chain]);
1131 oldSpeciesTree = GetTree (m->speciesTree, chain, state[chain] ^ 1);
1132
1133 // Get gene trees
1134 geneTrees = (Tree **) SafeCalloc (2*numGeneTrees, sizeof(Tree *));
1135 for (i=0; i<m->speciesTree->nSubParams; i++) {
1136 geneTrees[i] = GetTree(m->speciesTree->subParams[i], chain, state[chain]);
1137 }
1138
1139 // Allocate space for depth matrix copy
1140 numUpperTriang = numSpecies * (numSpecies - 1) / 2;
1141 oldMinDepths = (double *) SafeCalloc (2*numUpperTriang, sizeof(double));
1142 modMinDepths = oldMinDepths + numUpperTriang;
1143
1144 // Get min depth matrix for old gene trees
1145 GetMinDepthMatrix(geneTrees, numTopologies-1, depthMatrix);
1146
1147 // Save a copy
1148 for (i=0; i<numUpperTriang; i++)
1149 oldMinDepths[i] = depthMatrix[i];
1150
1151 // Get forward lambda
1152 GetMeanDist(oldSpeciesTree, depthMatrix, &mean);
1153 // if (mean < 1E-6) mean = 1E-6;
1154 forwardLambda = 1.0 / mean;
1155
1156 // Calculate joint probability of old gene trees and old species tree
1157 oldLnProb = LnJointGeneTreeSpeciesTreePr(geneTrees, numGeneTrees, oldSpeciesTree, chain);
1158
1159 // Modify the picked gene tree using code from a regular MrBayes move
1160 Move_ExtSPRClock(param, chain, seed, lnPriorRatio, lnProposalRatio, mvp);
1161
1162 // Update the min depth matrix
1163 GetMinDepthMatrix(geneTrees, numTopologies-1, depthMatrix);
1164
1165 // Copy the min depth matrix
1166 for (i=0; i<numUpperTriang; i++)
1167 modMinDepths[i] = depthMatrix[i];
1168
1169 // Modify the min depth matrix
1170 ModifyDepthMatrix (forwardLambda, modMinDepths, seed);
1171
1172 // Get a new species tree
1173 if (GetSpeciesTreeFromMinDepths (newSpeciesTree, modMinDepths) == ERROR) {
1174 abortMove = YES;
1175 free (geneTrees);
1176 free (oldMinDepths);
1177 return (NO_ERROR);
1178 }
1179
1180 // Calculate joint probability of new gene trees and new species tree
1181 newLnProb = LnJointGeneTreeSpeciesTreePr(geneTrees, numGeneTrees, newSpeciesTree, chain);
1182
1183 // Get backward lambda
1184 GetMeanDist(newSpeciesTree, depthMatrix, &mean);
1185 // if (mean < 1E-6) mean = 1E-6;
1186 backwardLambda = 1.0 / mean;
1187
1188 // Get proposal probability of old species tree
1189 backwardLnProposalProb = LnProposalProbSpeciesTree (oldSpeciesTree, oldMinDepths, backwardLambda);
1190
1191 // Get proposal probability of new species tree
1192 forwardLnProposalProb = LnProposalProbSpeciesTree (newSpeciesTree, depthMatrix, forwardLambda);
1193
1194 // Update prior ratio taking species tree into account
1195 (*lnPriorRatio) += (newLnProb - oldLnProb);
1196
1197 // Update proposal ratio based on this move
1198 (*lnProposalRatio) += (backwardLnProposalProb - forwardLnProposalProb);
1199
1200 // Free allocated memory
1201 free (geneTrees);
1202 free (oldMinDepths);
1203
1204 return (NO_ERROR);
1205 }
1206
1207
1208 /**-----------------------------------------------------------------
1209 |
1210 | Move_GeneTree2: Propose a new gene tree using NNIClock
1211 |
1212 | @param param The parameter to change
1213 | @param chain The chain number
1214 | @param seed Pointer to the seed of the random number gen.
1215 | @param lnPriorRatio Pointer to the log prior ratio (out)
1216 | @param lnProposalRatio Pointer to the log proposal (Hastings) ratio (out)
1217 | @param mvp Pointer to tuning parameter(s)
1218 ------------------------------------------------------------------*/
Move_GeneTree2(Param * param,int chain,RandLong * seed,MrBFlt * lnPriorRatio,MrBFlt * lnProposalRatio,MrBFlt * mvp)1219 int Move_GeneTree2 (Param *param, int chain, RandLong *seed, MrBFlt *lnPriorRatio, MrBFlt *lnProposalRatio, MrBFlt *mvp)
1220 {
1221 int i, numGeneTrees, numUpperTriang;
1222 double newLnProb, oldLnProb, backwardLnProposalProb, forwardLnProposalProb,
1223 *oldMinDepths, *modMinDepths, forwardLambda, backwardLambda, mean;
1224 Tree *newSpeciesTree, *oldSpeciesTree, **geneTrees;
1225 ModelInfo *m;
1226
1227 // Calculate number of gene trees
1228 numGeneTrees = numTopologies - 1;
1229
1230 // Get model settings
1231 m = &modelSettings[param->relParts[0]];
1232
1233 // Get species tree (this trick is possible because we always copy tree params)
1234 newSpeciesTree = GetTree (m->speciesTree, chain, state[chain]);
1235 oldSpeciesTree = GetTree (m->speciesTree, chain, state[chain] ^ 1);
1236
1237 // Get gene trees
1238 geneTrees = (Tree **) SafeCalloc (2*numGeneTrees, sizeof(Tree *));
1239 for (i=0; i<m->speciesTree->nSubParams; i++) {
1240 geneTrees[i] = GetTree(m->speciesTree->subParams[i], chain, state[chain]);
1241 }
1242
1243 // Allocate space for depth matrix copy
1244 numUpperTriang = numSpecies * (numSpecies - 1) / 2;
1245 oldMinDepths = (double *) SafeCalloc (2*numUpperTriang, sizeof(double));
1246 modMinDepths = oldMinDepths + numUpperTriang;
1247
1248 // Get min depth matrix for old gene trees
1249 GetMinDepthMatrix(geneTrees, numTopologies-1, depthMatrix);
1250
1251 // Save a copy
1252 for (i=0; i<numUpperTriang; i++)
1253 oldMinDepths[i] = depthMatrix[i];
1254
1255 // Get forward lambda
1256 GetMeanDist(oldSpeciesTree, depthMatrix, &mean);
1257 // if (mean < 1E-6) mean = 1E-6;
1258 forwardLambda = 1.0 / mean;
1259
1260 // Calculate joint probability of old gene trees and old species tree
1261 oldLnProb = LnJointGeneTreeSpeciesTreePr(geneTrees, numGeneTrees, oldSpeciesTree, chain);
1262
1263 // Modify the picked gene tree using code from a regular MrBayes move (no tuning parameter, so passing on mvp is OK)
1264 Move_NNIClock(param, chain, seed, lnPriorRatio, lnProposalRatio, mvp);
1265
1266 // Update the min depth matrix
1267 GetMinDepthMatrix(geneTrees, numTopologies-1, depthMatrix);
1268
1269 // Copy the min depth matrix
1270 for (i=0; i<numUpperTriang; i++)
1271 modMinDepths[i] = depthMatrix[i];
1272
1273 // Modify the min depth matrix
1274 ModifyDepthMatrix (forwardLambda, modMinDepths, seed);
1275
1276 // Get a new species tree
1277 if (GetSpeciesTreeFromMinDepths (newSpeciesTree, modMinDepths) == ERROR) {
1278 abortMove = YES;
1279 free (geneTrees);
1280 free (oldMinDepths);
1281 return (NO_ERROR);
1282 }
1283
1284 // Calculate joint probability of new gene trees and new species tree
1285 newLnProb = LnJointGeneTreeSpeciesTreePr(geneTrees, numGeneTrees, newSpeciesTree, chain);
1286
1287 // Get backward lambda
1288 GetMeanDist(newSpeciesTree, depthMatrix, &mean);
1289 // if (mean < 1E-6) mean = 1E-6;
1290 backwardLambda = 1.0 / mean;
1291
1292 // Get proposal probability of old species tree
1293 backwardLnProposalProb = LnProposalProbSpeciesTree (oldSpeciesTree, oldMinDepths, backwardLambda);
1294
1295 // Get proposal probability of new species tree
1296 forwardLnProposalProb = LnProposalProbSpeciesTree (newSpeciesTree, depthMatrix, forwardLambda);
1297
1298 // Update prior ratio taking species tree into account
1299 (*lnPriorRatio) += (newLnProb - oldLnProb);
1300
1301 // Update proposal ratio based on this move
1302 (*lnProposalRatio) += (backwardLnProposalProb - forwardLnProposalProb);
1303
1304 // Free allocated memory
1305 free (geneTrees);
1306 free (oldMinDepths);
1307
1308 return (NO_ERROR);
1309 }
1310
1311
1312 /**-----------------------------------------------------------------
1313 |
1314 | Move_GeneTree3: Propose a new gene tree using ParsSPRClock
1315 |
1316 | @param param The parameter to change
1317 | @param chain The chain number
1318 | @param seed Pointer to the seed of the random number gen.
1319 | @param lnPriorRatio Pointer to the log prior ratio (out)
1320 | @param lnProposalRatio Pointer to the log proposal (Hastings) ratio (out)
1321 | @param mvp Pointer to tuning parameter(s)
1322 ------------------------------------------------------------------*/
Move_GeneTree3(Param * param,int chain,RandLong * seed,MrBFlt * lnPriorRatio,MrBFlt * lnProposalRatio,MrBFlt * mvp)1323 int Move_GeneTree3 (Param *param, int chain, RandLong *seed, MrBFlt *lnPriorRatio, MrBFlt *lnProposalRatio, MrBFlt *mvp)
1324 {
1325 int i, numGeneTrees, numUpperTriang;
1326 double newLnProb, oldLnProb, backwardLnProposalProb, forwardLnProposalProb,
1327 *oldMinDepths, *modMinDepths, forwardLambda, backwardLambda, mean;
1328 Tree *newSpeciesTree, *oldSpeciesTree, **geneTrees;
1329 ModelInfo *m;
1330
1331 // Calculate number of gene trees
1332 numGeneTrees = numTopologies - 1;
1333
1334 // Get model settings
1335 m = &modelSettings[param->relParts[0]];
1336
1337 // Get species tree (this trick is possible because we always copy tree params)
1338 newSpeciesTree = GetTree (m->speciesTree, chain, state[chain]);
1339 oldSpeciesTree = GetTree (m->speciesTree, chain, state[chain] ^ 1);
1340
1341 // Get gene trees
1342 geneTrees = (Tree **) SafeCalloc (2*numGeneTrees, sizeof(Tree *));
1343 for (i=0; i<m->speciesTree->nSubParams; i++) {
1344 geneTrees[i] = GetTree(m->speciesTree->subParams[i], chain, state[chain]);
1345 }
1346
1347 // Allocate space for depth matrix copy
1348 numUpperTriang = numSpecies * (numSpecies - 1) / 2;
1349 oldMinDepths = (double *) SafeCalloc (2*numUpperTriang, sizeof(double));
1350 modMinDepths = oldMinDepths + numUpperTriang;
1351
1352 // Get min depth matrix for old gene trees
1353 GetMinDepthMatrix(geneTrees, numTopologies-1, depthMatrix);
1354
1355 // Save a copy
1356 for (i=0; i<numUpperTriang; i++)
1357 oldMinDepths[i] = depthMatrix[i];
1358
1359 // Get forward lambda
1360 GetMeanDist(oldSpeciesTree, depthMatrix, &mean);
1361 // if (mean < 1E-6) mean = 1E-6;
1362 forwardLambda = 1.0 / mean;
1363
1364 // Calculate joint probability of old gene trees and old species tree
1365 oldLnProb = LnJointGeneTreeSpeciesTreePr(geneTrees, numGeneTrees, oldSpeciesTree, chain);
1366
1367 // Modify the picked gene tree using code from a regular MrBayes move
1368 Move_ParsSPRClock(param, chain, seed, lnPriorRatio, lnProposalRatio, mvp);
1369
1370 // Update the min depth matrix
1371 GetMinDepthMatrix(geneTrees, numTopologies-1, depthMatrix);
1372
1373 // Copy the min depth matrix
1374 for (i=0; i<numUpperTriang; i++)
1375 modMinDepths[i] = depthMatrix[i];
1376
1377 // Modify the min depth matrix
1378 ModifyDepthMatrix (forwardLambda, modMinDepths, seed);
1379
1380 // Get a new species tree
1381 if (GetSpeciesTreeFromMinDepths (newSpeciesTree, modMinDepths) == ERROR) {
1382 abortMove = YES;
1383 free (geneTrees);
1384 free (oldMinDepths);
1385 return (NO_ERROR);
1386 }
1387
1388 // Calculate joint probability of new gene trees and new species tree
1389 newLnProb = LnJointGeneTreeSpeciesTreePr(geneTrees, numGeneTrees, newSpeciesTree, chain);
1390
1391 // Get backward lambda
1392 GetMeanDist(newSpeciesTree, depthMatrix, &mean);
1393 // if (mean < 1E-6) mean = 1E-6;
1394 backwardLambda = 1.0 / mean;
1395
1396 // Get proposal probability of old species tree
1397 backwardLnProposalProb = LnProposalProbSpeciesTree (oldSpeciesTree, oldMinDepths, backwardLambda);
1398
1399 // Get proposal probability of new species tree
1400 forwardLnProposalProb = LnProposalProbSpeciesTree (newSpeciesTree, depthMatrix, forwardLambda);
1401
1402 // Update prior ratio taking species tree into account
1403 (*lnPriorRatio) += (newLnProb - oldLnProb);
1404
1405 // Update proposal ratio based on this move
1406 (*lnProposalRatio) += (backwardLnProposalProb - forwardLnProposalProb);
1407
1408 // Free allocated memory
1409 free (geneTrees);
1410 free (oldMinDepths);
1411
1412 return (NO_ERROR);
1413 }
1414
1415
1416 /*-----------------------------------------------------------------------------------
1417 |
1418 | Move_NodeSliderGeneTree: Move the position of one (root or nonroot) node in a
1419 | gene tree inside a species tree.
1420 |
1421 -------------------------------------------------------------------------------------*/
Move_NodeSliderGeneTree(Param * param,int chain,RandLong * seed,MrBFlt * lnPriorRatio,MrBFlt * lnProposalRatio,MrBFlt * mvp)1422 int Move_NodeSliderGeneTree (Param *param, int chain, RandLong *seed, MrBFlt *lnPriorRatio, MrBFlt *lnProposalRatio, MrBFlt *mvp)
1423 {
1424 int i, *nEvents;
1425 MrBFlt window, minDepth, maxDepth, oldDepth, newDepth,
1426 oldLeftLength=0.0, oldRightLength=0.0, clockRate,
1427 oldPLength=0.0, lambda=0.0, nu=0.0, igrvar=0.0,
1428 *brlens=NULL, *tk02Rate=NULL, *igrRate=NULL, *popSizePtr;
1429 TreeNode *p, *q;
1430 ModelInfo *m;
1431 Tree *geneTree, *speciesTree;
1432 Param *subParm;
1433
1434 window = mvp[0]; /* window size */
1435
1436 m = &modelSettings[param->relParts[0]];
1437
1438 /* get gene tree and species tree */
1439 geneTree = GetTree (param, chain, state[chain]);
1440 speciesTree = GetTree (m->speciesTree, chain, state[chain]);
1441
1442 /* get population size(s) */
1443 popSizePtr = GetParamVals(m->popSize, chain, state[chain]);
1444
1445 /* get clock rate */
1446 if (m->clockRate == NULL)
1447 clockRate = 1.0;
1448 else
1449 clockRate = *GetParamVals(m->clockRate, chain, state[chain]);
1450
1451 /* pick a node to be changed */
1452 p = geneTree->intDownPass[(int)(RandomNumber(seed)*geneTree->nIntNodes)];
1453
1454 # if defined (DEBUG_CSLIDER)
1455 printf ("Before node slider (gene tree):\n");
1456 printf ("Picked branch with index %d and depth %f\n", p->index, p->nodeDepth);
1457 if (p->anc->anc == NULL)
1458 printf ("Old clock rate: %f\n", clockRate);
1459 ShowNodes (t->root, 0, t->isRooted);
1460 getchar();
1461 # endif
1462
1463 /* get gene tree prior prob before move */
1464 (*lnPriorRatio) -= LnPriorProbGeneTree(geneTree, clockRate, speciesTree, popSizePtr);
1465
1466 /* store values needed later for prior calculation (relaxed clocks) */
1467 oldPLength = p->length;
1468 if (p->left != NULL)
1469 {
1470 oldLeftLength = p->left->length;
1471 oldRightLength = p->right->length;
1472 }
1473 else
1474 oldLeftLength = oldRightLength = 0.0;
1475
1476 /* find species tree branch to which the gene tree node belongs */
1477 MapGeneTreeToSpeciesTree(geneTree, speciesTree);
1478 q = NULL;
1479 for (i=0; i<speciesTree->nNodes-1; i++)
1480 {
1481 q = speciesTree->allDownPass[i];
1482 if (p->x == q->index)
1483 break;
1484 }
1485 assert (q != NULL && p->x == q->index);
1486
1487 /* determine lower and upper bound */
1488 /* FIXME: p->left == NULL? (from clang static analyzer) */
1489 minDepth = p->left->nodeDepth + POS_MIN;
1490 if (p->right->nodeDepth + POS_MIN > minDepth)
1491 minDepth = p->right->nodeDepth + POS_MIN;
1492 if (q->nodeDepth + POS_MIN > minDepth)
1493 minDepth = q->nodeDepth + POS_MIN;
1494 if (p->anc->anc == NULL)
1495 maxDepth = TREEHEIGHT_MAX;
1496 else
1497 maxDepth = p->anc->nodeDepth - POS_MIN;
1498
1499 /* abort if impossible */
1500 if (minDepth >= maxDepth)
1501 {
1502 abortMove = YES;
1503 return (NO_ERROR);
1504 }
1505
1506 if (maxDepth-minDepth < window)
1507 {
1508 window = maxDepth-minDepth;
1509 }
1510
1511 /* pick the new node depth */
1512 oldDepth = p->nodeDepth;
1513 newDepth = oldDepth + (RandomNumber(seed) - 0.5) * window;
1514
1515 /* reflect the new node depth */
1516 while (newDepth < minDepth || newDepth > maxDepth)
1517 {
1518 if (newDepth < minDepth)
1519 newDepth = 2.0 * minDepth - newDepth;
1520 if (newDepth > maxDepth)
1521 newDepth = 2.0 * maxDepth - newDepth;
1522 }
1523 p->nodeDepth = newDepth;
1524
1525 /* determine new branch lengths around p and set update of transition probabilities */
1526 if (p->left != NULL)
1527 {
1528 p->left->length = p->nodeDepth - p->left->nodeDepth;
1529 assert (p->left->length >= POS_MIN);
1530 p->left->upDateTi = YES;
1531 p->right->length = p->nodeDepth - p->right->nodeDepth;
1532 assert (p->right->length >= POS_MIN);
1533 p->right->upDateTi = YES;
1534 }
1535 if (p->anc->anc != NULL)
1536 {
1537 p->length = p->anc->nodeDepth - p->nodeDepth;
1538 assert (p->length >= POS_MIN);
1539 p->upDateTi = YES;
1540 }
1541
1542 /* set flags for update of cond likes from p and down to root */
1543 q = p;
1544 while (q->anc != NULL)
1545 {
1546 q->upDateCl = YES;
1547 q = q->anc;
1548 }
1549
1550 /* calculate proposal ratio */
1551 (*lnProposalRatio) = 0.0;
1552
1553 /* calculate prior ratio */
1554 (*lnPriorRatio) += LnPriorProbGeneTree (geneTree, clockRate, speciesTree, popSizePtr);
1555
1556 /* adjust proposal and prior ratio for relaxed clock models */
1557 for (i=0; i<param->nSubParams; i++)
1558 {
1559 subParm = param->subParams[i];
1560 if (subParm->paramType == P_CPPEVENTS)
1561 {
1562 nEvents = subParm->nEvents[2*chain+state[chain]];
1563 lambda = *GetParamVals (modelSettings[subParm->relParts[0]].cppRate, chain, state[chain]);
1564
1565 /* proposal ratio */
1566 if (p->left != NULL)
1567 {
1568 (*lnProposalRatio) += nEvents[p->left->index ] * log (p->left->length / oldLeftLength);
1569 (*lnProposalRatio) += nEvents[p->right->index] * log (p->right->length / oldRightLength);
1570 }
1571 if (p->anc->anc != NULL)
1572 (*lnProposalRatio) += nEvents[p->index] * log (p->length / oldPLength);
1573
1574 /* prior ratio */
1575 if (p->anc->anc == NULL) // two branches changed in same direction
1576 (*lnPriorRatio) += lambda * (2.0 * (oldDepth - newDepth));
1577 else if (p->left != NULL) // two branches changed in one direction, one branch in the other direction
1578 (*lnPriorRatio) += lambda * (oldDepth - newDepth);
1579 else /* if (p->left == NULL) */ // one branch changed
1580 (*lnPriorRatio) += lambda * (newDepth - oldDepth);
1581
1582 /* update effective evolutionary lengths */
1583 if (UpdateCppEvolLengths (subParm, p, chain) == ERROR)
1584 {
1585 abortMove = YES;
1586 return (NO_ERROR);
1587 }
1588 }
1589 else if ( subParm->paramType == P_TK02BRANCHRATES ||
1590 (subParm->paramType == P_MIXEDBRCHRATES && *GetParamIntVals(subParm, chain, state[chain]) == RCL_TK02))
1591 {
1592 if (subParm->paramType == P_TK02BRANCHRATES)
1593 nu = *GetParamVals (modelSettings[subParm->relParts[0]].tk02var, chain, state[chain]);
1594 else
1595 nu = *GetParamVals (modelSettings[subParm->relParts[0]].mixedvar, chain, state[chain]);
1596 tk02Rate = GetParamVals (subParm, chain, state[chain]);
1597 brlens = GetParamSubVals (subParm, chain, state[chain]);
1598
1599 /* no proposal ratio effect */
1600
1601 /* prior ratio */
1602 if (p->left != NULL)
1603 {
1604 (*lnPriorRatio) -= LnProbTK02LogNormal (tk02Rate[p->index], nu*oldLeftLength, tk02Rate[p->left->index]);
1605 (*lnPriorRatio) -= LnProbTK02LogNormal (tk02Rate[p->index], nu*oldRightLength, tk02Rate[p->right->index]);
1606 (*lnPriorRatio) += LnProbTK02LogNormal (tk02Rate[p->index], nu*p->left->length, tk02Rate[p->left->index]);
1607 (*lnPriorRatio) += LnProbTK02LogNormal (tk02Rate[p->index], nu*p->right->length, tk02Rate[p->right->index]);
1608 }
1609 if (p->anc->anc != NULL)
1610 {
1611 (*lnPriorRatio) -= LnProbTK02LogNormal (tk02Rate[p->anc->index], nu*oldPLength, tk02Rate[p->index]);
1612 (*lnPriorRatio) += LnProbTK02LogNormal (tk02Rate[p->anc->index], nu*p->length, tk02Rate[p->index]);
1613 }
1614
1615 /* update effective evolutionary lengths */
1616 if (p->left != NULL)
1617 {
1618 brlens[p->left->index] = p->left->length * (tk02Rate[p->left->index]+tk02Rate[p->index])/2.0;
1619 brlens[p->right->index] = p->right->length * (tk02Rate[p->right->index]+tk02Rate[p->index])/2.0;
1620 }
1621 if (p->anc->anc != NULL)
1622 {
1623 brlens[p->index] = p->length * (tk02Rate[p->index]+tk02Rate[p->anc->index])/2.0;
1624 }
1625 }
1626 else if ( subParm->paramType == P_IGRBRANCHRATES ||
1627 (subParm->paramType == P_MIXEDBRCHRATES && *GetParamIntVals(subParm, chain, state[chain]) == RCL_IGR))
1628 {
1629 if (subParm->paramType == P_IGRBRANCHRATES)
1630 igrvar = *GetParamVals (modelSettings[subParm->relParts[0]].igrvar, chain, state[chain]);
1631 else
1632 igrvar = *GetParamVals (modelSettings[subParm->relParts[0]].mixedvar, chain, state[chain]);
1633 igrRate = GetParamVals (subParm, chain, state[chain]);
1634 brlens = GetParamSubVals (subParm, chain, state[chain]);
1635
1636 if (p->left != NULL)
1637 {
1638 (*lnPriorRatio) -= LnProbGamma (oldLeftLength/igrvar, oldLeftLength/igrvar, igrRate[p->left->index ]);
1639 (*lnPriorRatio) -= LnProbGamma (oldRightLength/igrvar, oldRightLength/igrvar, igrRate[p->right->index]);
1640 (*lnPriorRatio) += LnProbGamma (p->left->length/igrvar, p->left->length/igrvar, igrRate[p->left->index ]);
1641 (*lnPriorRatio) += LnProbGamma (p->right->length/igrvar, p->right->length/igrvar, igrRate[p->right->index]);
1642 }
1643 if (p->anc->anc != NULL)
1644 {
1645 (*lnPriorRatio) -= LnProbGamma (oldPLength/igrvar, oldPLength/igrvar, igrRate[p->index]);
1646 (*lnPriorRatio) += LnProbGamma (p->length /igrvar, p->length /igrvar, igrRate[p->index]);
1647 }
1648
1649 if (p->left != NULL)
1650 {
1651 brlens[p->left->index ] = igrRate[p->left->index ] * p->left->length;
1652 brlens[p->right->index] = igrRate[p->right->index] * p->right->length;
1653 }
1654 if (p->anc->anc != NULL)
1655 {
1656 brlens[p->index] = igrRate[p->index] * p->length;
1657 }
1658 }
1659 }
1660
1661 # if defined (DEBUG_CSLIDER)
1662 printf ("After node slider (gene tree):\n");
1663 printf ("Old depth: %f -- New depth: %f -- LnPriorRatio %f -- LnProposalRatio %f\n",
1664 oldDepth, newDepth, (*lnPriorRatio), (*lnProposalRatio));
1665 ShowNodes (t->root, 0, t->isRooted);
1666 getchar();
1667 # endif
1668
1669 return (NO_ERROR);
1670
1671 }
1672
1673
1674 /*------------------------------------------------------------------
1675 |
1676 | Move_SpeciesTree: Propose a new species tree
1677 |
1678 ------------------------------------------------------------------*/
Move_SpeciesTree(Param * param,int chain,RandLong * seed,MrBFlt * lnPriorRatio,MrBFlt * lnProposalRatio,MrBFlt * mvp)1679 int Move_SpeciesTree (Param *param, int chain, RandLong *seed, MrBFlt *lnPriorRatio, MrBFlt *lnProposalRatio, MrBFlt *mvp)
1680 {
1681 int i, numGeneTrees, numUpperTriang;
1682 double newLnProb, oldLnProb, backwardLnProposalProb, forwardLnProposalProb, *modMinDepths,
1683 forwardLambda, backwardLambda, lambdaDiv, mean;
1684 Tree *newSpeciesTree, *oldSpeciesTree, **geneTrees;
1685 ModelInfo *m;
1686
1687 /* get tuning parameter (lambda divider) */
1688 lambdaDiv = mvp[0];
1689
1690 /* calculate number of gene trees */
1691 numGeneTrees = param->nSubParams;
1692
1693 /* get model settings */
1694 m = &modelSettings[param->relParts[0]];
1695
1696 /* get new and old species trees */
1697 newSpeciesTree = GetTree (m->speciesTree, chain, state[chain]);
1698 oldSpeciesTree = GetTree (m->speciesTree, chain, state[chain] ^ 1);
1699
1700 /* get gene trees */
1701 geneTrees = (Tree **) SafeCalloc (numGeneTrees, sizeof(Tree*));
1702 for (i=0; i<param->nSubParams; i++)
1703 geneTrees[i] = GetTree(param->subParams[i], chain, state[chain]);
1704
1705 /* get minimum depth matrix */
1706 GetMinDepthMatrix(geneTrees, numGeneTrees, depthMatrix);
1707
1708 /* get forward lambda */
1709 GetMeanDist(oldSpeciesTree, depthMatrix, &mean);
1710 forwardLambda = 1.0 / (mean * lambdaDiv);
1711
1712 /* make a copy for modification */
1713 numUpperTriang = numSpecies * (numSpecies - 1) / 2;
1714 modMinDepths = (double *) SafeCalloc (numUpperTriang, sizeof(double));
1715 for (i=0; i<numUpperTriang; i++)
1716 modMinDepths[i] = depthMatrix[i];
1717
1718 /* modify minimum depth matrix */
1719 ModifyDepthMatrix (forwardLambda, modMinDepths, seed);
1720
1721 /* construct a new species tree from the modified constraints */
1722 if (GetSpeciesTreeFromMinDepths(newSpeciesTree, modMinDepths) == ERROR) {
1723 abortMove = YES;
1724 free (modMinDepths);
1725 free (geneTrees);
1726 return (NO_ERROR);
1727 }
1728
1729 /* get lambda for back move */
1730 GetMeanDist(newSpeciesTree, depthMatrix, &mean);
1731 backwardLambda = 1.0 / (mean * lambdaDiv);
1732
1733 /* calculate proposal ratio */
1734 backwardLnProposalProb = LnProposalProbSpeciesTree (oldSpeciesTree, depthMatrix, backwardLambda);
1735 forwardLnProposalProb = LnProposalProbSpeciesTree (newSpeciesTree, depthMatrix, forwardLambda);
1736 (*lnProposalRatio) = backwardLnProposalProb - forwardLnProposalProb;
1737
1738 # if defined (BEST_MPI_ENABLED)
1739 // Broadcast the proposed species tree to all processors if MPI version
1740 # endif
1741
1742 # if defined (BEST_MPI_ENABLED)
1743 // Let each processor calculate the ln probability ratio of its current gene tree(s)
1744 // given the new and old species tree in the MPI version
1745
1746 // Assemble the ln probability ratios across the processors and to lnPriorRatio
1747 # else
1748 /* calculate the ln probability ratio of the current gene trees
1749 given the new and old species trees */
1750 newLnProb = LnJointGeneTreeSpeciesTreePr(geneTrees, numGeneTrees, newSpeciesTree, chain);
1751 oldLnProb = LnJointGeneTreeSpeciesTreePr(geneTrees, numGeneTrees, oldSpeciesTree, chain);
1752 # endif
1753
1754 /* set (*lnPriorRatio) to ln probability ratio */
1755 (*lnPriorRatio) = (newLnProb - oldLnProb);
1756
1757 /* free allocated space */
1758 free (modMinDepths);
1759 free (geneTrees);
1760
1761 return (NO_ERROR);
1762 }
1763
1764
1765 /** Show upper triangular matrix */
ShowUpperTriangMatrix(double * values,int squareSize)1766 void ShowUpperTriangMatrix (double *values, int squareSize)
1767 {
1768 int i, j, index;
1769
1770 index = 0;
1771 printf ("Upper triang matrix:\n");
1772 for (i=0; i<squareSize; i++) {
1773 for (j=0; j<i; j++)
1774 printf (" ");
1775 for (j=i+1; j<squareSize; j++) {
1776 printf ("%.6f ", values[index]);
1777 index++;
1778 }
1779 printf ("\n");
1780 }
1781 printf ("\n");
1782 }
1783
1784