C++ TopologyNode::getParent方法代码示例

std::set<TopologyNode*> SpeciesNarrowExchangeProposal::getOldestSubtreesNodesInPopulation( Tree &tau, TopologyNode &n )
{
    
    // I need all the oldest nodes/subtrees that have the same tips.
    // Those nodes need to be scaled too.
    
    // get the beginning and ending age of the population
    double max_age = -1.0;
    if ( n.isRoot() == false )
    {
        max_age = n.getParent().getAge();
    }
    
    // get all the taxa from the species tree that are descendants of node i
    std::vector<TopologyNode*> species_taxa;
    TreeUtilities::getTaxaInSubtree( &n, species_taxa );
    
    // get all the individuals
    std::set<TopologyNode*> individualTaxa;
    for (size_t i = 0; i < species_taxa.size(); ++i)
    {
        const std::string &name = species_taxa[i]->getName();
        std::vector<TopologyNode*> ind = tau.getTipNodesWithSpeciesName( name );
        for (size_t j = 0; j < ind.size(); ++j)
        {
            individualTaxa.insert( ind[j] );
        }
    }
    
    // create the set of the nodes within this population
    std::set<TopologyNode*> nodesInPopulationSet;
    
    // now go through all nodes in the gene
    while ( individualTaxa.empty() == false )
    {
        std::set<TopologyNode*>::iterator it = individualTaxa.begin();
        individualTaxa.erase( it );
        
        TopologyNode *geneNode = *it;
        
        // add this node to our list of node we need to scale, if:
        // a) this is the root node
        // b) this is not the root and the age of the parent node is larger than the parent's age of the species node
        if ( geneNode->isRoot() == false && ( max_age == -1.0 || max_age > geneNode->getParent().getAge() ) )
        {
            // push the parent to our current list
            individualTaxa.insert( &geneNode->getParent() );
        }
        else if ( geneNode->isTip() == false )
        {
            // add this node if it is within the age of our population
            nodesInPopulationSet.insert( geneNode );
        }
        
    }
    
    return nodesInPopulationSet;
}

void MultivariateBrownianPhyloProcess::recursiveSimulate(const TopologyNode& from)  {
    
    size_t index = from.getIndex();
    if (from.isRoot())    {
        
        std::vector<double>& val = (*value)[index];
        for (size_t i=0; i<getDim(); i++)   {
            val[i] = 0;
        }
    }
    
    else    {
        
        // x ~ normal(x_up, sigma^2 * branchLength)

        std::vector<double>& val = (*value)[index];
                
        sigma->getValue().drawNormalSampleCovariance((*value)[index]);

        size_t upindex = from.getParent().getIndex();
        std::vector<double>& upval = (*value)[upindex];

        for (size_t i=0; i<getDim(); i++)   {
            val[i] += upval[i];
        }        
    }
    
    // propagate forward
    size_t numChildren = from.getNumberOfChildren();
    for (size_t i = 0; i < numChildren; ++i) {
        recursiveSimulate(from.getChild(i));
    }
    
}

void PhyloBrownianProcessREML::recursivelyFlagNodeDirty( const TopologyNode &n )
{
    
    // we need to flag this node and all ancestral nodes for recomputation
    size_t index = n.getIndex();
    
    // if this node is already dirty, the also all the ancestral nodes must have been flagged as dirty
    if ( !dirtyNodes[index] )
    {
        // the root doesn't have an ancestor
        if ( !n.isRoot() )
        {
            recursivelyFlagNodeDirty( n.getParent() );
        }
        
        // set the flag
        dirtyNodes[index] = true;
        
        // if we previously haven't touched this node, then we need to change the active likelihood pointer
        if ( changedNodes[index] == false )
        {
            activeLikelihood[index] = (activeLikelihood[index] == 0 ? 1 : 0);
            changedNodes[index] = true;
        }
        
    }
    
}

double BrownianPhyloProcess::recursiveLnProb( const TopologyNode& from ) {
    
    double lnProb = 0.0;
    size_t index = from.getIndex();
    double val = (*value)[index];

    if (! from.isRoot()) {
        
        // x ~ normal(x_up, sigma^2 * branchLength)
        
        size_t upindex = from.getParent().getIndex();
        double standDev = sigma->getValue() * sqrt(from.getBranchLength());
        double mean = (*value)[upindex] + drift->getValue() * from.getBranchLength();
        lnProb += RbStatistics::Normal::lnPdf(val, standDev, mean);
    }
    
    // propagate forward
    size_t numChildren = from.getNumberOfChildren();
    
    for (size_t i = 0; i < numChildren; ++i) {
        lnProb += recursiveLnProb(from.getChild(i));
    }
    
    return lnProb;
    
}

void BrownianPhyloProcess::recursiveSimulate(const TopologyNode& from)  {
    
    size_t index = from.getIndex();
    
    if (! from.isRoot())    {
        
        // x ~ normal(x_up, sigma^2 * branchLength)
        
        size_t upindex = from.getParent().getIndex();
        double standDev = sigma->getValue() * sqrt(from.getBranchLength());
        double mean = (*value)[upindex] + drift->getValue() * from.getBranchLength();

        // simulate the new Val
        RandomNumberGenerator* rng = GLOBAL_RNG;
        (*value)[index] = RbStatistics::Normal::rv( mean, standDev, *rng);
     
    }
    
    // propagate forward
    size_t numChildren = from.getNumberOfChildren();
    for (size_t i = 0; i < numChildren; ++i) {
        recursiveSimulate(from.getChild(i));
    }
    
}

void AutocorrelatedBranchMatrixDistribution::recursiveSimulate(const TopologyNode& node, RbVector< RateMatrix > *values, const std::vector< double > &scaledParent) {
    
    // get the index
    size_t nodeIndex = node.getIndex();
    
    // first we simulate our value
    RandomNumberGenerator* rng = GLOBAL_RNG;
    // do we keep our parents values?
    double u = rng->uniform01();
    if ( u < changeProbability->getValue() ) {
        // change
        
        // draw a new value for the base frequencies
        std::vector<double> newParent = RbStatistics::Dirichlet::rv(scaledParent, *rng);
        std::vector<double> newScaledParent = newParent;
        
        // compute the new scaled parent
        std::vector<double>::iterator end = newScaledParent.end();
        for (std::vector<double>::iterator it = newScaledParent.begin(); it != end; ++it) {
            (*it) *= alpha->getValue();
        }
        
        RateMatrix_GTR rm = RateMatrix_GTR( newParent.size() );
        RbPhylogenetics::Gtr::computeRateMatrix( exchangeabilityRates->getValue(), newParent, &rm );
        
        uniqueBaseFrequencies.push_back( newParent );
        uniqueMatrices.push_back( rm );
        matrixIndex[nodeIndex] = uniqueMatrices.size()-1;
        values->insert(nodeIndex, rm);
        
        size_t numChildren = node.getNumberOfChildren();
        if ( numChildren > 0 ) {
            
            for (size_t i = 0; i < numChildren; ++i) {
                const TopologyNode& child = node.getChild(i);
                recursiveSimulate(child,values,newScaledParent);
            }
        }
        
    }
    else {
        // no change
        size_t parentIndex = node.getParent().getIndex();
        values->insert(nodeIndex, uniqueMatrices[ matrixIndex[ parentIndex ] ]);
        
        size_t numChildren = node.getNumberOfChildren();
        if ( numChildren > 0 ) {
            
            for (size_t i = 0; i < numChildren; ++i) {
                const TopologyNode& child = node.getChild(i);
                recursiveSimulate(child,values,scaledParent);
            }
        }
    }
    
}

/**
 * Perform the proposal.
 *
 * A Beta-simplex proposal randomly changes some values of a simplex, although the other values
 * change too because of the renormalization.
 * First, some random indices are drawn. Then, the proposal draws a new somplex
 *   u ~ Beta(val[index] * alpha)
 * where alpha is the tuning parameter.The new value is set to u.
 * The simplex is then renormalized.
 *
 * \return The hastings ratio.
 */
double SubtreeScaleProposal::doProposal( void )
{
    // Get random number generator
    RandomNumberGenerator* rng     = GLOBAL_RNG;
    
    TimeTree& tau = variable->getValue();
    
    // pick a random node which is not the root and neither the direct descendant of the root
    TopologyNode* node;
    do {
        double u = rng->uniform01();
        size_t index = size_t( std::floor(tau.getNumberOfNodes() * u) );
        node = &tau.getNode(index);
    } while ( node->isRoot() || node->isTip() );
    
    TopologyNode& parent = node->getParent();
    
    // we need to work with the times
    double parent_age  = parent.getAge();
    double my_age      = node->getAge();
    
    // now we store all necessary values
    storedNode = node;
    storedAge = my_age;
    
    // lower bound
    double min_age = 0.0;
    TreeUtilities::getOldestTip(&tau, node, min_age);
    
    // draw new ages and compute the hastings ratio at the same time
    double my_new_age = min_age + (parent_age - min_age) * rng->uniform01();
    
    double scalingFactor = my_new_age / my_age;
    
    size_t nNodes = node->getNumberOfNodesInSubtree(false);
    
    // rescale the subtrees
    TreeUtilities::rescaleSubtree(&tau, node, scalingFactor );
    
    if (min_age != 0.0)
    {
        for (size_t i = 0; i < tau.getNumberOfTips(); i++)
        {
            if (tau.getNode(i).getAge() < 0.0) {
                return RbConstants::Double::neginf;
            }
        }
    }
    
    // compute the Hastings ratio
    double lnHastingsratio = (nNodes > 1 ? log( scalingFactor ) * (nNodes-1) : 0.0 );
    
    return lnHastingsratio;
    
}

double MultivariateBrownianPhyloProcess::recursiveLnProb( const TopologyNode& from ) {
    
    double lnProb = 0.0;
    size_t index = from.getIndex();
    std::vector<double> val = (*value)[index];
    
    if (! from.isRoot()) {
        
        if (1)  {
//        if (dirtyNodes[index])  {

            // x ~ normal(x_up, sigma^2 * branchLength)

            size_t upindex = from.getParent().getIndex();
            std::vector<double> upval = (*value)[upindex];

            const MatrixReal& om = sigma->getValue().getInverse();

            double s2 = 0;
            for (size_t i = 0; i < getDim(); i++) {
                double tmp = 0;
                for (size_t j = 0; j < getDim(); j++) {
                    tmp += om[i][j] * (val[j] - upval[j]);
                }
                s2 += (val[i] - upval[i]) * tmp;
            }

            double logprob = 0;
            logprob -= 0.5 * s2 / from.getBranchLength();
            logprob -= 0.5 * (sigma->getValue().getLogDet() + sigma->getValue().getDim() * log(from.getBranchLength()));
            nodeLogProbs[index] = logprob;
            dirtyNodes[index] = false;
        }
        lnProb += nodeLogProbs[index];
    }
    
    // propagate forward
    size_t numChildren = from.getNumberOfChildren();
    
    for (size_t i = 0; i < numChildren; ++i) {
        lnProb += recursiveLnProb(from.getChild(i));
    }
    
    return lnProb;
    
}

/**
 * Compute the diversity of the tree at time t.
 *
 * \param[in]    t      time at which we want to know the diversity.
 *
 * \return The diversity (number of species in the reconstructed tree).
 */
int PiecewiseConstantSerialSampledBirthDeathProcess::survivors(double t) const
{

    const std::vector<TopologyNode*>& nodes = value->getNodes();

    int survivors = 0;
    for (std::vector<TopologyNode*>::const_iterator it = nodes.begin(); it != nodes.end(); ++it)
    {
        TopologyNode* n = *it;
        if ( n->getAge() < t )
        {
            if ( n->isRoot() || n->getParent().getAge() > t )
            {
                survivors++;
            }
        }
    }

    return survivors;
}

/** Perform the move */
double SubtreeScale::performSimpleMove( void ) {
    
    // Get random number generator    
    RandomNumberGenerator* rng     = GLOBAL_RNG;
    
    TimeTree& tau = variable->getValue();
    
    // pick a random node which is not the root and neither the direct descendant of the root
    TopologyNode* node;
    do {
        double u = rng->uniform01();
        size_t index = size_t( std::floor(tau.getNumberOfNodes() * u) );
        node = &tau.getNode(index);
    } while ( node->isRoot() || node->isTip() );
    
    TopologyNode& parent = node->getParent();
    
    // we need to work with the times
    double parent_age  = parent.getAge();
    double my_age      = node->getAge();
    
    // now we store all necessary values
    storedNode = node;
    storedAge = my_age;
        
    // draw new ages and compute the hastings ratio at the same time
    double my_new_age = parent_age * rng->uniform01();
    
    double scalingFactor = my_new_age / my_age;
    
    size_t nNodes = node->getNumberOfNodesInSubtree(false);
    
    // rescale the subtrees
    TreeUtilities::rescaleSubtree(&tau, node, scalingFactor );
    
    // compute the Hastings ratio
    double lnHastingsratio = (nNodes > 1 ? log( scalingFactor ) * (nNodes-1) : 0.0 );
    
    return lnHastingsratio;
}

/**
 * Perform the proposal.
 *
 * A Uniform-simplex proposal randomly changes some values of a simplex, although the other values
 * change too because of the renormalization.
 * First, some random indices are drawn. Then, the proposal draws a new somplex
 *   u ~ Uniform(val[index] * alpha)
 * where alpha is the tuning parameter.The new value is set to u.
 * The simplex is then renormalized.
 *
 * \return The hastings ratio.
 */
double NodeTimeSlideUniformProposal::doProposal( void )
{

    // Get random number generator
    RandomNumberGenerator* rng     = GLOBAL_RNG;

    Tree& tau = variable->getValue();

    // pick a random node which is not the root and neithor the direct descendant of the root
    TopologyNode* node;
    do {
        double u = rng->uniform01();
        size_t index = size_t( std::floor(tau.getNumberOfNodes() * u) );
        node = &tau.getNode(index);
    } while ( node->isRoot() || node->isTip() );

    TopologyNode& parent = node->getParent();

    // we need to work with the times
    double parent_age  = parent.getAge();
    double my_age      = node->getAge();
    double child_Age   = node->getChild( 0 ).getAge();
    if ( child_Age < node->getChild( 1 ).getAge())
    {
        child_Age = node->getChild( 1 ).getAge();
    }

    // now we store all necessary values
    storedNode = node;
    storedAge = my_age;

    // draw new ages and compute the hastings ratio at the same time
    double my_new_age = (parent_age-child_Age) * rng->uniform01() + child_Age;

    // set the age
    tau.getNode(node->getIndex()).setAge( my_new_age );

    return 0.0;

}

std::vector<TopologyNode*> TreeNodeAgeUpdateProposal::getNodesInPopulation( Tree &tau, TopologyNode &n )
{

    // I need all the oldest nodes/subtrees that have the same tips.
    // Those nodes need to be scaled too.

    // get the beginning and ending age of the population
    double max_age = -1.0;
    if ( n.isRoot() == false )
    {
        max_age = n.getParent().getAge();
    }

    // get all the taxa from the species tree that are descendants of node i
    double min_age_left = n.getChild(0).getAge();
    std::vector<TopologyNode*> speciesTaxa_left;
    TreeUtilities::getTaxaInSubtree( &n.getChild(0), speciesTaxa_left );

    // get all the individuals
    std::set<TopologyNode*> individualTaxa_left;
    for (size_t i = 0; i < speciesTaxa_left.size(); ++i)
    {
        const std::string &name = speciesTaxa_left[i]->getName();
        std::vector<TopologyNode*> ind = tau.getTipNodesWithSpeciesName( name );
        for (size_t j = 0; j < ind.size(); ++j)
        {
            individualTaxa_left.insert( ind[j] );
        }
    }

    // create the set of the nodes within this population
    std::set<TopologyNode*> nodesInPopulationSet;

    // now go through all nodes in the gene
    while ( individualTaxa_left.empty() == false )
    {
        // get the first element
        std::set<TopologyNode*>::iterator it = individualTaxa_left.begin();

        // store the pointer
        TopologyNode *geneNode = *it;

        // and now remove the element from the list
        individualTaxa_left.erase( it );

        // add this node to our list of node we need to scale, if:
        // a) this is the root node
        // b) this is not the root and the age of the parent node is larger than the parent's age of the species node
        if ( geneNode->getAge() > min_age_left && geneNode->getAge() < max_age && geneNode->isTip() == false )
        {
            // add this node if it is within the age of our population
            nodesInPopulationSet.insert( geneNode );
        }

        if ( geneNode->isRoot() == false && ( max_age == -1.0 || max_age > geneNode->getParent().getAge() ) )
        {
            // push the parent to our current list
            individualTaxa_left.insert( &geneNode->getParent() );
        }

    }

    // get all the taxa from the species tree that are descendants of node i
    double min_age_right = n.getChild(1).getAge();
    std::vector<TopologyNode*> speciesTaxa_right;
    TreeUtilities::getTaxaInSubtree( &n.getChild(1), speciesTaxa_right );

    // get all the individuals
    std::set<TopologyNode*> individualTaxa_right;
    for (size_t i = 0; i < speciesTaxa_right.size(); ++i)
    {
        const std::string &name = speciesTaxa_right[i]->getName();
        std::vector<TopologyNode*> ind = tau.getTipNodesWithSpeciesName( name );
        for (size_t j = 0; j < ind.size(); ++j)
        {
            individualTaxa_right.insert( ind[j] );
        }
    }

    // now go through all nodes in the gene
    while ( individualTaxa_right.empty() == false )
    {
        // get the first element
        std::set<TopologyNode*>::iterator it = individualTaxa_right.begin();

        // store the pointer
        TopologyNode *geneNode = *it;

        // and now remove the element from the list
        individualTaxa_right.erase( it );

        // add this node to our list of node we need to scale, if:
        // a) this is the root node
        // b) this is not the root and the age of the parent node is larger than the parent's age of the species node
        if ( geneNode->getAge() > min_age_right && geneNode->getAge() < max_age && geneNode->isTip() == false )
        {
            // add this node if it is within the age of our population
            nodesInPopulationSet.insert( geneNode );
        }

//.........这里部分代码省略.........

/** Perform the move */
void RateAgeBetaShift::performMcmcMove( double lHeat, double pHeat )
{
    
    // Get random number generator
    RandomNumberGenerator* rng     = GLOBAL_RNG;
    
    Tree& tau = tree->getValue();
    RbOrderedSet<DagNode*> affected;
    tree->getAffectedNodes( affected );
    
    double oldLnLike = 0.0;
    bool checkLikelihoodShortcuts = rng->uniform01() < 0.001;
    if ( checkLikelihoodShortcuts == true )
    {
        for (RbOrderedSet<DagNode*>::iterator it = affected.begin(); it != affected.end(); ++it)
        {
            (*it)->touch();
            oldLnLike += (*it)->getLnProbability();
        }
    }
    
    // pick a random node which is not the root and neithor the direct descendant of the root
    TopologyNode* node;
    size_t nodeIdx = 0;
    do {
        double u = rng->uniform01();
        nodeIdx = size_t( std::floor(tau.getNumberOfNodes() * u) );
        node = &tau.getNode(nodeIdx);
    } while ( node->isRoot() || node->isTip() ); 
    
    TopologyNode& parent = node->getParent();
    
    // we need to work with the times
    double parent_age  = parent.getAge();
    double my_age      = node->getAge();
    double child_Age   = node->getChild( 0 ).getAge();
    if ( child_Age < node->getChild( 1 ).getAge())
    {
        child_Age = node->getChild( 1 ).getAge();
    }
    
    // now we store all necessary values
    storedNode = node;
    storedAge = my_age;
    
    
    storedRates[nodeIdx] = rates[nodeIdx]->getValue();
    for (size_t i = 0; i < node->getNumberOfChildren(); i++)
    {
        size_t childIdx = node->getChild(i).getIndex();
        storedRates[childIdx] = rates[childIdx]->getValue();
    }
    
    
    // draw new ages and compute the hastings ratio at the same time
    double m = (my_age-child_Age) / (parent_age-child_Age);
    double a = delta * m + 1.0;
    double b = delta * (1.0-m) + 1.0;
    double new_m = RbStatistics::Beta::rv(a, b, *rng);
    double my_new_age = (parent_age-child_Age) * new_m + child_Age;
    
    // compute the Hastings ratio
    double forward = RbStatistics::Beta::lnPdf(a, b, new_m);
    double new_a = delta * new_m + 1.0;
    double new_b = delta * (1.0-new_m) + 1.0;
    double backward = RbStatistics::Beta::lnPdf(new_a, new_b, m);
    
    // set the age
    tau.getNode(nodeIdx).setAge( my_new_age );
    
    // touch the tree so that the likelihoods are getting stored
    tree->touch();
    
    // get the probability ratio of the tree
    double treeProbRatio = tree->getLnProbabilityRatio();
    
    
    // set the rates
    double pa = node->getParent().getAge();
    double my_new_rate =(pa - my_age) * storedRates[nodeIdx] / (pa - my_new_age);
    
    // now we set the new value
    // this will automcatically call a touch
    rates[nodeIdx]->setValue( new double( my_new_rate ) );
    
    // get the probability ratio of the new rate
    double ratesProbRatio = rates[nodeIdx]->getLnProbabilityRatio();
    
    for (size_t i = 0; i < node->getNumberOfChildren(); i++)
    {
        size_t childIdx = node->getChild(i).getIndex();
        double a = node->getChild(i).getAge();
        double child_new_rate = (my_age - a) * storedRates[childIdx] / (my_new_age - a);
        
        // now we set the new value
        // this will automcatically call a touch
        rates[childIdx]->setValue( new double( child_new_rate ) );

        // get the probability ratio of the new rate
//.........这里部分代码省略.........

/**
 * Perform the proposal.
 *
 * \return The hastings ratio.
 */
double SpeciesNarrowExchangeProposal::doProposal( void )
{
    
    // empty the previous vectors
    storedGeneTreeNodes.clear();
    storedOldBrothers.clear();
    
    // Get random number generator
    RandomNumberGenerator* rng     = GLOBAL_RNG;
    
    Tree& tau = speciesTree->getValue();
    
    // pick a random node which is not the root and neithor a direct descendant of the root
    TopologyNode* node;
    do {
        double u = rng->uniform01();
        size_t index = size_t( std::floor(tau.getNumberOfNodes() * u) );
        node = &tau.getNode(index);
    } while ( node->isRoot() || node->getParent().isRoot() );
    
    TopologyNode& parent = node->getParent();
    TopologyNode& grandparent = parent.getParent();
    TopologyNode* uncle = &grandparent.getChild( 0 );
    // check if we got the correct child
    if ( uncle == &parent )
    {
        uncle = &grandparent.getChild( 1 );
    }
    TopologyNode* brother = &parent.getChild( 0 );
    // check if we got the correct child
    if ( brother == node )
    {
        brother = &parent.getChild( 1 );
    }
    
    // we need to work with the times
    double parent_age   = parent.getAge();
    double uncles_age   = uncle->getAge();
    
    if( uncles_age < parent_age )
    {
        failed = false;
        
        double lnHastingsRatio = 0.0;
        
        // now we store all necessary values
        storedChoosenNode   = node;
        storedUncle         = uncle;
        
        // now we need to find for each gene tree the nodes that need to be moved as well
        // only nodes that have a coalescent event within the lifetime of the parents populations
        // from lineages belonging to the chosen node with lineages belonging to the brother population
        // need to be changed
        for ( size_t i=0; i<geneTrees.size(); ++i )
        {
            // get the i-th gene tree
            Tree& geneTree = geneTrees[i]->getValue();
            
            std::vector<TopologyNode*> nodes = getNodesToChange(geneTree, *node, *brother );
            
            // get the set of nodes in my uncles populations
            // these are the nodes that are possible re-attachment points
            std::set<TopologyNode*> new_siblings = getOldestSubtreesNodesInPopulation(geneTree, *uncle);
            std::set<TopologyNode*> old_siblings = getOldestSubtreesNodesInPopulation(geneTree, *brother);
            
            for (size_t j=0; j<nodes.size(); ++j)
            {
                
                TopologyNode *the_gene_node = nodes[i];

                // first we need to compute the backward probability
                std::set<TopologyNode*> old_candidate_siblings = getPossibleSiblings(the_gene_node, old_siblings);
                
                // add the backward probability to the hastings ratio
                lnHastingsRatio += log( old_siblings.size() );
                
                // then we need to compute the forward probability
                std::set<TopologyNode*> new_candidate_siblings = getPossibleSiblings(the_gene_node, new_siblings);
                
                // add the forward  probability to the hastings ratio
                lnHastingsRatio += log( new_candidate_siblings.size() );
                
                // actually pick a new sibling
                size_t new_index = size_t( floor(rng->uniform01() * new_candidate_siblings.size() ) );
                std::set<TopologyNode*>::iterator it = new_candidate_siblings.begin();
                std::advance(it,new_index);
                TopologyNode *new_child = *it;
                
                // store nodes
                storedGeneTreeNodes.push_back( the_gene_node );
                TopologyNode &the_parent = the_gene_node->getParent();
                TopologyNode *old_brother = &the_parent.getChild( 0 );
                if ( old_brother == the_gene_node )
                {
                    old_brother = &the_parent.getChild( 1 );
//.........这里部分代码省略.........

/** Perform the move */
void RateAgeBetaShift::performMove( double heat, bool raiseLikelihoodOnly )
{

    // Get random number generator
    RandomNumberGenerator* rng     = GLOBAL_RNG;

    TimeTree& tau = tree->getValue();

    // pick a random node which is not the root and neithor the direct descendant of the root
    TopologyNode* node;
    size_t nodeIdx = 0;
    do {
        double u = rng->uniform01();
        nodeIdx = size_t( std::floor(tau.getNumberOfNodes() * u) );
        node = &tau.getNode(nodeIdx);
    } while ( node->isRoot() || node->isTip() );

    TopologyNode& parent = node->getParent();

    // we need to work with the times
    double parent_age  = parent.getAge();
    double my_age      = node->getAge();
    double child_Age   = node->getChild( 0 ).getAge();
    if ( child_Age < node->getChild( 1 ).getAge())
    {
        child_Age = node->getChild( 1 ).getAge();
    }

    // now we store all necessary values
    storedNode = node;
    storedAge = my_age;


    storedRates[nodeIdx] = rates[nodeIdx]->getValue();
    for (size_t i = 0; i < node->getNumberOfChildren(); i++)
    {
        size_t childIdx = node->getChild(i).getIndex();
        storedRates[childIdx] = rates[childIdx]->getValue();
    }


    // draw new ages and compute the hastings ratio at the same time
    double m = (my_age-child_Age) / (parent_age-child_Age);
    double a = delta * m + 1.0;
    double b = delta * (1.0-m) + 1.0;
    double new_m = RbStatistics::Beta::rv(a, b, *rng);
    double my_new_age = (parent_age-child_Age) * new_m + child_Age;

    // compute the Hastings ratio
    double forward = RbStatistics::Beta::lnPdf(a, b, new_m);
    double new_a = delta * new_m + 1.0;
    double new_b = delta * (1.0-new_m) + 1.0;
    double backward = RbStatistics::Beta::lnPdf(new_a, new_b, m);

    // set the age
    tau.setAge( node->getIndex(), my_new_age );
    tree->touch();

    double treeProbRatio = tree->getLnProbabilityRatio();


    // set the rates
    rates[nodeIdx]->setValue( new double((node->getParent().getAge() - my_age) * storedRates[nodeIdx] / (node->getParent().getAge() - my_new_age)));
    double ratesProbRatio = rates[nodeIdx]->getLnProbabilityRatio();

    for (size_t i = 0; i < node->getNumberOfChildren(); i++)
    {
        size_t childIdx = node->getChild(i).getIndex();
        rates[childIdx]->setValue( new double((my_age - node->getChild(i).getAge()) * storedRates[childIdx] / (my_new_age - node->getChild(i).getAge())));
        ratesProbRatio += rates[childIdx]->getLnProbabilityRatio();

    }

    std::set<DagNode*> affected;
    tree->getAffectedNodes( affected );
    double lnProbRatio = 0;
    for (std::set<DagNode*>::iterator it = affected.begin(); it != affected.end(); ++it)
    {
        (*it)->touch();
        lnProbRatio += (*it)->getLnProbabilityRatio();
    }

    if ( fabs(lnProbRatio) > 1E-6 ) {
//        throw RbException("Likelihood shortcut computation failed in rate-age-proposal.");
        std::cout << "Likelihood shortcut computation failed in rate-age-proposal." << std::endl;
    }

    double hastingsRatio = backward - forward;
    double lnAcceptanceRatio = treeProbRatio + ratesProbRatio + hastingsRatio;

    if (lnAcceptanceRatio >= 0.0)
    {
        numAccepted++;

        tree->keep();
        rates[nodeIdx]->keep();
        for (size_t i = 0; i < node->getNumberOfChildren(); i++)
        {
            size_t childIdx = node->getChild(i).getIndex();
//.........这里部分代码省略.........