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 = &params[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