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

void GibbsPruneAndRegraft::findNewBrothers(std::vector<TopologyNode *> &b, TopologyNode &p, TopologyNode *n) {
    
    // security check that I'm not a tip
    if (!n->isTip() && &p != n) 
    {
        // check the first child
        TopologyNode* child = &n->getChild( 0 );
        if ( child->getAge() < p.getAge() ) 
        {
            b.push_back( child );
        } 
        else 
        {
            findNewBrothers(b, p, child);
        }
        
        // check the second child
        child = &n->getChild( 1 );
        if ( child->getAge() < p.getAge() ) 
        {
            b.push_back( child );
        } 
        else 
        {
            findNewBrothers(b, p, child);
        }
    }
}

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
}

/**
 * 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 RootTimeSlideUniformProposal::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 = &tau.getRoot();
    
    // we need to work with the times
    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
    storedAge = my_age;
    
    // draw new ages and compute the hastings ratio at the same time
    double my_new_age = (origin->getValue() - child_Age) * rng->uniform01() + child_Age;
    
    // set the age
    node->setAge( my_new_age );
    
    return 0.0;
    
}

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

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

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

std::set<TopologyNode*> SpeciesNarrowExchangeProposal::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 the node
    double min_age = n.getAge();
    std::vector<TopologyNode*> species_taxa;
    TreeUtilities::getTaxaInSubtree( &n, species_taxa );
    
    // get all the individuals
    std::set<TopologyNode*> individual_taxa;
    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)
        {
            individual_taxa.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 ( individual_taxa.empty() == false )
    {
        std::set<TopologyNode*>::iterator it = individual_taxa.begin();
        individual_taxa.erase( it );
        
        TopologyNode *gene_node = *it;
        
        if ( gene_node->getAge() < min_age )
        {
            // the age of the node is younger than the populations start age
            // -> add the node to the current working set
            individual_taxa.insert( &gene_node->getParent() );
            
        }
        else if ( gene_node->getAge() < max_age || max_age == -1.0 )
        {
            // the age of the node is within the population
            // -> add the node to the current return set
            nodesInPopulationSet.insert( gene_node );
            
            // if this is not the root then we need to add the parent node to the working set
            if ( gene_node->isRoot() == false )
            {
                individual_taxa.insert( &gene_node->getParent() );
                
            }
            
        }
        
    }
    
    return nodesInPopulationSet;
}

void AbstractRootedTreeDistribution::simulateTree( void )
{

    // the time tree object (topology & times)
    Tree *psi = new Tree();

    // internally we treat unrooted topologies the same as rooted
    psi->setRooted( true );

    // create the tip nodes
    std::vector<TopologyNode*> nodes;
    for (size_t i=0; i<num_taxa; ++i)
    {

        // create the i-th taxon
        TopologyNode* node = new TopologyNode( taxa[i], i );

        // set the age of this tip node
        node->setAge( taxa[i].getAge() );

        if (node->getAge() > 0)
        {
            node->setFossil(true);
        }

        // add the new node to the list
        nodes.push_back( node );

    }


    double ra = root_age->getValue();
    double present = ra;

    // we need a sorted vector of constraints, starting with the smallest
    std::vector<Clade> sorted_clades;

    std::vector<Clade> constraints;

    for (size_t i = 0; i < constraints.size(); ++i)
    {
        if (constraints[i].getAge() > ra)
        {
            throw RbException("Cannot simulate tree: clade constraints are older than the root age.");
        }

        // set the ages of each of the taxa in the constraint
        for (size_t j = 0; j < constraints[i].size(); ++j)
        {
            for (size_t k = 0; k < num_taxa; k++)
            {
                if ( taxa[k].getName() == constraints[i].getTaxonName(j) )
                {
                    constraints[i].setTaxonAge(j, taxa[k].getAge());
                    break;
                }
            }
        }

        if ( constraints[i].size() > 1 && constraints[i].size() < num_taxa )
        {
            sorted_clades.push_back( constraints[i] );
        }
    }

    // create a clade that contains all species
    Clade all_species = Clade(taxa);
    all_species.setAge( ra );
    sorted_clades.push_back(all_species);

    // next sort the clades
    std::sort(sorted_clades.begin(),sorted_clades.end());

    // remove duplicates
    std::vector<Clade> tmp;
    tmp.push_back( sorted_clades[0] );
    for (size_t i = 1; i < sorted_clades.size(); ++i)
    {
        Clade &a = tmp[tmp.size()-1];
        Clade &b = sorted_clades[i];

        if ( a.size() != b.size() )
        {
            tmp.push_back( sorted_clades[i] );
        }
        else
        {
            bool equal = true;
            for (size_t j = 0; j < a.size(); ++j)
            {
                if ( a.getTaxon(j) != b.getTaxon(j) )
                {
                    equal = false;
                    break;
                }
            }
            if ( equal == false )
            {
                tmp.push_back( sorted_clades[i] );
            }
//.........这里部分代码省略.........

void StateDependentSpeciationExtinctionProcess::recursivelyDrawJointConditionalAncestralStates(const TopologyNode &node, std::vector<size_t>& startStates, std::vector<size_t>& endStates)
{
    
    size_t node_index = node.getIndex();
    
    if ( node.isTip() == true )
    {
        const AbstractHomologousDiscreteCharacterData& data = static_cast<TreeDiscreteCharacterData*>(this->value)->getCharacterData();
        const AbstractDiscreteTaxonData& taxon_data = data.getTaxonData( node.getName() );
        
        const DiscreteCharacterState &char_state = taxon_data.getCharacter(0);
        
        // we need to treat ambiguous state differently
        if ( char_state.isAmbiguous() == false )
        {
            endStates[node_index] = char_state.getStateIndex();
        }
        else
        {
            // initialize the conditional likelihoods for this branch
            state_type branch_conditional_probs = std::vector<double>(2 * num_states, 0);
            size_t start_state = startStates[node_index];
            branch_conditional_probs[ num_states + start_state ] = 1.0;
            
            // first calculate extinction likelihoods via a backward time pass
            double end_age = node.getParent().getAge();
            numericallyIntegrateProcess(branch_conditional_probs, 0, end_age, true, true);
            
            // now calculate conditional likelihoods along branch in forward time
            end_age        = node.getParent().getAge() - node.getAge();
            numericallyIntegrateProcess(branch_conditional_probs, 0, end_age, false, false);
            
            double total_prob = 0.0;
            for (size_t i = 0; i < num_states; ++i)
            {
                if ( char_state.isMissingState() == true || char_state.isGapState() == true || char_state.isStateSet(i) == true )
                {
                    total_prob += branch_conditional_probs[i];
                }
            }
            
            RandomNumberGenerator* rng = GLOBAL_RNG;
            size_t u = rng->uniform01() * total_prob;
            
            for (size_t i = 0; i < num_states; ++i)
            {
                
                if ( char_state.isMissingState() == true || char_state.isGapState() == true || char_state.isStateSet(i) == true )
                {
                    u -= branch_conditional_probs[i];
                    if ( u <= 0.0 )
                    {
                        endStates[node_index] = i;
                        break;
                    }
                    
                }
                
            }
            
        }
    }
    else
    {
        // sample characters by their probability conditioned on the branch's start state going to end states
        
        // initialize the conditional likelihoods for this branch
        state_type branch_conditional_probs = std::vector<double>(2 * num_states, 0);
        size_t start_state = startStates[node_index];
        branch_conditional_probs[ num_states + start_state ] = 1.0;
        
        // first calculate extinction likelihoods via a backward time pass
        double end_age = node.getParent().getAge();
        numericallyIntegrateProcess(branch_conditional_probs, 0, end_age, true, true);
        
        // now calculate conditional likelihoods along branch in forward time
        end_age        = node.getParent().getAge() - node.getAge();
        numericallyIntegrateProcess(branch_conditional_probs, 0, end_age, false, false);
        
        // TODO: if character mapping compute likelihoods for each time slice....
//        double current_time_slice = floor(begin_age / dt);
//        bool computed_at_least_one = false;
//        
//        // first iterate forward along the branch subtending this node to get the
//        // probabilities of the end states conditioned on the start state
//        while (current_time_slice * dt >= end_age || !computed_at_least_one)
//        {
//            double begin_age_slice = current_time_slice * dt;
//            double end_age_slice = (current_time_slice + 1) * dt;
//            numericallyIntegrateProcess(branch_conditional_probs, begin_age_slice, end_age_slice, false);
//            
//            computed_at_least_one = true;
//            current_time_slice--;
//        }
        
        std::map<std::vector<unsigned>, double> event_map;
        std::vector<double> speciation_rates;
        if ( use_cladogenetic_events == true )
        {
            // get cladogenesis event map (sparse speciation rate matrix)
//.........这里部分代码省略.........

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

void RateMap_Biogeography::calculateTransitionProbabilities(const TopologyNode& node, TransitionProbabilityMatrix &P, size_t charIdx) const
{

    double startAge = ( node.isRoot() ? 1e10 : node.getParent().getAge() );
    double endAge = node.getAge();
    double currAge = startAge;
    
    double r = ( branchHeterogeneousClockRates ? heterogeneousClockRates[node.getIndex()] : homogeneousClockRate );
    const std::vector<double>& glr = ( branchHeterogeneousGainLossRates ? heterogeneousGainLossRates[node.getIndex()] : homogeneousGainLossRates );
    
    // start at earliest epoch
    int epochIdx = getEpochIndex(startAge);
    

    // initialize P = I
    P[0][0] = 1.0;
    P[0][1] = 0.0;
    P[1][0] = 0.0;
    P[1][1] = 1.0;
    

    
    bool stop = false;
    while (!stop)
    {
        // get dispersal and extinction rates for site
        size_t idx = this->numCharacters * epochIdx + charIdx;
        
        double dispersalRate  = ( (availableAreaVector[ idx ] > 0) ? 1.0 : 0.0);
        double extinctionRate = ( (availableAreaVector[ idx ] > 0) ? 1.0 : 1e10);
        
        // get age of start of next oldest epoch
        double epochAge = epochs[ epochIdx ];
        
        // get next oldest age boundary (node or epoch)
        double incrAge = epochAge;
        
        // no more epochs, exit loop after computation
        if (endAge >= epochAge)
        {
            incrAge = endAge;
            stop = true;
        }
        
        // get branch length
        double diffAge = currAge - incrAge;
        
        // transition probabilities w/ sum-product
        double glr0 = glr[0] * extinctionRate;
        double glr1 = glr[1] * dispersalRate;
        double expPart = exp( -(glr0 + glr1) * r * diffAge);
        double p = glr0 / (glr0 + glr1);
        double q = 1.0 - p;
        
        double p00 = (p + q * expPart);
        double p01 = (q - q * expPart);
        double p10 = (p - p * expPart);
        double p11 = (q + p * expPart);

        double q00 = P[0][0];
        double q01 = P[0][1];
        double q10 = P[1][0];
        double q11 = P[1][1];
        
        P[0][0] = p00 * q00 + p01 * q10;
        P[0][1] = p00 * q01 + p01 * q11;
        P[1][0] = p10 * q00 + p11 * q10;
        P[1][1] = p10 * q01 + p11 * q11;
        
//        std::cout << P[0][0] << "," << P[0][1] << ";" << P[1][0] << "," << P[1][1] << "\n";
        
        
        
        if (!stop)
        {
            epochIdx += 1;
            currAge = epochAge;
        }
    }
}