src/best.c

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