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

void RateAgeACLNMixingMove::rejectCompoundMove( void ) {
    	
    // undo the proposal
	TimeTree& tau = tree->getValue();
	std::vector<double>& nrates = rates->getValue();
	double &rootR = rootRate->getValue();
	
	size_t nn = tau.getNumberOfNodes();
	double c = storedC;
	
	for(size_t i=0; i<nn; i++){
		TopologyNode* node = &tau.getNode(i);
		if(node->isTip() == false){
			double curAge = node->getAge();
			double undoAge = curAge / c;
			tau.setAge( node->getIndex(), undoAge );
		}
	}
	
	size_t nr = nrates.size();
	rootR = rootR * c;
	for(size_t i=0; i<nr; i++){
        double curRt = nrates[i];
        double undoRt = curRt * c;
        nrates[i] = undoRt;
	}
	
#ifdef ASSERTIONS_TREE
    if ( fabs(storedRootAge - tau.getRoot().getAge()) > 1E-8 ) {
        throw RbException("Error while rejecting RateAgeACLNMixingMove proposal: Node ages were not correctly restored!");
    }
#endif
}

void RealNodeContainer::recursiveClampAt(const TopologyNode& from, const ContinuousCharacterData* data, size_t l) {
 
    if (from.isTip())   {
        
        // get taxon index
        size_t index = from.getIndex();
        std::string taxon = tree->getTipNames()[index];
        size_t dataindex = data->getIndexOfTaxon(taxon);
        
        if (data->getCharacter(dataindex,l) != -1000) {
           (*this)[index] = data->getCharacter(dataindex,l);
            clampVector[index] = true;
            //std::cerr << "taxon : " << index << '\t' << taxon << " trait value : " << (*this)[index] << '\n';
        }
        else    {
            std::cerr << "taxon : " << taxon << " is missing for trait " << l+1 << '\n';
        }
    }

    // propagate forward
    size_t numChildren = from.getNumberOfChildren();
    for (size_t i = 0; i < numChildren; ++i) {
        recursiveClampAt(from.getChild(i),data,l);
    }    
}

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;
}

/**
 * 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;
    
}

/**
 * Recursive call to attach ordered interior node times to the time tree psi. Call it initially with the
 * root of the tree.
 */
void HeterogeneousRateBirthDeath::attachTimes(Tree* psi, std::vector<TopologyNode *> &nodes, size_t index, const std::vector<double> &interiorNodeTimes, double originTime )
{
    
    if (index < num_taxa-1)
    {
        // Get the rng
        RandomNumberGenerator* rng = GLOBAL_RNG;
        
        // Randomly draw one node from the list of nodes
        size_t node_index = static_cast<size_t>( floor(rng->uniform01()*nodes.size()) );
        
        // Get the node from the list
        TopologyNode* parent = nodes.at(node_index);
        psi->getNode( parent->getIndex() ).setAge( originTime - interiorNodeTimes[index] );
        
        // Remove the randomly drawn node from the list
        nodes.erase(nodes.begin()+long(node_index));
        
        // Add the left child if an interior node
        TopologyNode* leftChild = &parent->getChild(0);
        if ( !leftChild->isTip() )
        {
            nodes.push_back(leftChild);
        }
        
        // Add the right child if an interior node
        TopologyNode* rightChild = &parent->getChild(1);
        if ( !rightChild->isTip() )
        {
            nodes.push_back(rightChild);
        }
        
        // Recursive call to this function
        attachTimes(psi, nodes, index+1, interiorNodeTimes, originTime);
    }
    
}

void RealNodeContainer::recursiveGetStatsOverTips(const TopologyNode& from, double& e1, double& e2, int& n) const {

    if(from.isTip())   {
        double tmp = (*this)[from.getIndex()];

        n++;
        e1 += tmp;
        e2 += tmp * tmp;
    }
    // propagate forward
    size_t numChildren = from.getNumberOfChildren();
    for (size_t i = 0; i < numChildren; ++i) {
        recursiveGetStatsOverTips(from.getChild(i),e1,e2,n);
    }
    
}

void RealNodeContainer::recursiveGetTipValues(const TopologyNode& from, ContinuousCharacterData& nameToVal) const {
    
    if(from.isTip())   {
        double tmp = (*this)[from.getIndex()];
        std::string name =  tree->getTipNames()[from.getIndex()];
        
        ContinuousTaxonData dataVec = ContinuousTaxonData(name);
        double contObs = tmp;
        dataVec.addCharacter( contObs );
        nameToVal.addTaxonData( dataVec );
        return;
    }
    // propagate forward
    size_t numChildren = from.getNumberOfChildren();
    for (size_t i = 0; i < numChildren; ++i) {
        recursiveGetTipValues(from.getChild(i), nameToVal );
    }
    
}

std::string RealNodeContainer::recursiveGetNewick(const TopologyNode& from) const {

    std::ostringstream s;
    
    if (from.isTip())   {
        s << getTimeTree()->getTipNames()[from.getIndex()] << "_";
//        std::cerr << from.getIndex() << '\t' << getTimeTree()->getTipNames()[from.getIndex()] << "_";
//        std::cerr << (*this)[from.getIndex()] << '\n';
//        exit(1);
    }
    else    {
        s << "(";
        // propagate forward
        size_t numChildren = from.getNumberOfChildren();
        for (size_t i = 0; i < numChildren; ++i) {
            s << recursiveGetNewick(from.getChild(i));
            if (i < numChildren-1)  {
                s << ",";
            }
        }
        s << ")";
    }
    s << (*this)[from.getIndex()];
/*    if (from.isTip() && (! isClamped(from.getIndex()))) {
        std::cerr << "leaf is not clamped\n";
        // get taxon index
        size_t index = from.getIndex();
        std::cerr << "index : " << index << '\n';
        std::string taxon = tree->getTipNames()[index];
        std::cerr << "taxon : " << index << '\t' << taxon << '\n';
        std::cerr << " trait value : " << (*this)[index] << '\n';        
        exit(1);
    }*/
//    if (!from.isRoot()) {
        s << ":";
        s << getTimeTree()->getBranchLength(from.getIndex());
//    }
    
    return s.str();
}

/** 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;

}

/** Perform the move */
double RateAgeACLNMixingMove::performCompoundMove( void ) {
    
    // Get random number generator    
    RandomNumberGenerator* rng     = GLOBAL_RNG;
    
    TimeTree& tau = tree->getValue();
	std::vector<double>& nrates = rates->getValue();
	double &rootR = rootRate->getValue();
	
	
	size_t nn = tau.getNumberOfNodes();
	double u = rng->uniform01();
	double c = exp( epsilon * (u - 0.5) );
	
	for(size_t i=0; i<nn; i++){
		TopologyNode* node = &tau.getNode(i);
		if(node->isTip() == false){
			double curAge = node->getAge();
			double newAge = curAge * c;
			tau.setAge( node->getIndex(), newAge );
			if(node->isRoot()){
				storedRootAge = curAge;
			}
		}
	}
	
	size_t nr = nrates.size();
	rootR = rootR / c;
	for(size_t i=0; i<nrates.size(); i++){
        double curRt = nrates[i];
        double newRt = curRt / c;
        nrates[i] = newRt;
	}
	
	storedC = c;
	double pr = (((double)nn) - (double)nr) * log(c);
	return pr;
}

void PhyloBrownianProcessREML::recursiveComputeLnProbability( const TopologyNode &node, size_t nodeIndex )
{

    // check for recomputation
    if ( node.isTip() == false && dirtyNodes[nodeIndex] )
    {
        // mark as computed
        dirtyNodes[nodeIndex] = false;

        std::vector<double> &p_node  = this->partialLikelihoods[this->activeLikelihood[nodeIndex]][nodeIndex];
        std::vector<double> &mu_node  = this->contrasts[this->activeLikelihood[nodeIndex]][nodeIndex];

        
        // get the number of children
        size_t num_children = node.getNumberOfChildren();
        
        for (size_t j = 1; j < num_children; ++j)
        {
        
            size_t leftIndex = nodeIndex;
            const TopologyNode *left = &node;
            if ( j == 1 )
            {
                left = &node.getChild(0);
                leftIndex = left->getIndex();
                recursiveComputeLnProbability( *left, leftIndex );
            }
            
            const TopologyNode &right = node.getChild(j);
            size_t rightIndex = right.getIndex();
            recursiveComputeLnProbability( right, rightIndex );

            const std::vector<double> &p_left  = this->partialLikelihoods[this->activeLikelihood[leftIndex]][leftIndex];
            const std::vector<double> &p_right = this->partialLikelihoods[this->activeLikelihood[rightIndex]][rightIndex];

            // get the per node and site contrasts
            const std::vector<double> &mu_left  = this->contrasts[this->activeLikelihood[leftIndex]][leftIndex];
            const std::vector<double> &mu_right = this->contrasts[this->activeLikelihood[rightIndex]][rightIndex];

            // get the propagated uncertainties
            double delta_left  = this->contrastUncertainty[this->activeLikelihood[leftIndex]][leftIndex];
            double delta_right = this->contrastUncertainty[this->activeLikelihood[rightIndex]][rightIndex];

            // get the scaled branch lengths
            double v_left  = 0;
            if ( j == 1 )
            {
                v_left = this->computeBranchTime(leftIndex, left->getBranchLength());
            }
            double v_right = this->computeBranchTime(rightIndex, right.getBranchLength());

            // add the propagated uncertainty to the branch lengths
            double t_left  = v_left  + delta_left;
            double t_right = v_right + delta_right;

            // set delta_node = (t_l*t_r)/(t_l+t_r);
            this->contrastUncertainty[this->activeLikelihood[nodeIndex]][nodeIndex] = (t_left*t_right) / (t_left+t_right);

            double stdev = sqrt(t_left+t_right);
            for (int i=0; i<this->numSites; i++)
            {

                mu_node[i] = (mu_left[i]*t_right + mu_right[i]*t_left) / (t_left+t_right);

                // get the site specific rate of evolution
                double standDev = this->computeSiteRate(i) * stdev;

                // compute the contrasts for this site and node
                double contrast = mu_left[i] - mu_right[i];

                // compute the probability for the contrasts at this node
                double lnl_node = RbStatistics::Normal::lnPdf(0, standDev, contrast);

                // sum up the probabilities of the contrasts
                p_node[i] = lnl_node + p_left[i] + p_right[i];

            } // end for-loop over all sites

        } // end for-loop over all children
        
    } // end if we need to compute something for this node.

}

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 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();
//.........这里部分代码省略.........