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

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 PhyloWhiteNoiseProcess::recursiveSimulate(const TopologyNode& from)
{
    
    if (! from.isRoot())
    {
        // get the index
        size_t index = from.getIndex();
    
        // compute the variance along the branch
        double mean = 1.0;
        double stdev = sigma->getValue() / sqrt(from.getBranchLength());
        double alpha = mean * mean / (stdev * stdev);
        double beta = mean / (stdev * stdev);
    
        // simulate a new Val
        RandomNumberGenerator* rng = GLOBAL_RNG;
        double v = RbStatistics::Gamma::rv( alpha,beta, *rng);
    
        // we store this val here
        (*value)[index] = v;
    
    }
    
    // simulate the val for each child (if any)
    size_t numChildren = from.getNumberOfChildren();
    for (size_t i = 0; i < numChildren; ++i)
    {
        const TopologyNode& child = from.getChild(i);
        recursiveSimulate(child);
    }
    
}

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

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

void RateMap_Biogeography::calculateTransitionProbabilities(const TopologyNode& node, TransitionProbabilityMatrix &P) const
{
   
    double branchLength = node.getBranchLength();
    
    double r = ( branchHeterogeneousClockRates ? heterogeneousClockRates[node.getIndex()] : homogeneousClockRate );
    const std::vector<double>& glr = ( branchHeterogeneousGainLossRates ? heterogeneousGainLossRates[node.getIndex()] : homogeneousGainLossRates );

    if (node.isRoot())
        branchLength = 1e10;

    double expPart = exp( -(glr[0] + glr[1]) * r * branchLength);
    double p = glr[0] / (glr[0] + glr[1]);
    double q = 1.0 - p;
    
    P[0][0] = p + q * expPart;
    P[0][1] = q - q * expPart;
    P[1][0] = p - p * expPart;
    P[1][1] = q + p * expPart;
}

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

}

double PhyloWhiteNoiseProcess::recursiveLnProb(const TopologyNode &from)   {

    double lnProb = 0.0;
    if (! from.isRoot())   {
        // compute the variance
        double mean = 1.0;
        double stdev = sigma->getValue() / sqrt(from.getBranchLength());
        double alpha = mean * mean / (stdev * stdev);
        double beta = mean / (stdev * stdev);
        double v = (*value)[from.getIndex()];
        lnProb += log( RbStatistics::Gamma::lnPdf(alpha,beta,v) );
    }
    
    size_t numChildren = from.getNumberOfChildren();
    for (size_t i = 0; i < numChildren; ++i) {
        const TopologyNode& child = from.getChild(i);
        lnProb += recursiveLnProb(child);
            
    }
    
    return lnProb;
}

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

/** Perform the move */
double GibbsPruneAndRegraft::performSimpleMove( void ) 
{
    
    // Get random number generator    
    RandomNumberGenerator* rng     = GLOBAL_RNG;
    
    TimeTree& tau = variable->getValue();
    
    // potential affected nodes for likelihood computation
    std::set<DagNode *> affected;
    variable->getAffectedNodes( affected );
    
    double backwardLikelihood = variable->getLnProbability();
    for (std::set<DagNode*>::const_iterator it = affected.begin(); it != affected.end(); ++it) 
    {
        backwardLikelihood += (*it)->getLnProbability();
    }
    int offset = (int) -backwardLikelihood;
    double backward = exp(backwardLikelihood + offset);
    
    // 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->getParent().isRoot() );
    
    TopologyNode* parent        = &node->getParent(); 
    TopologyNode& grandparent   = parent->getParent();
    TopologyNode& brother       = parent->getChild( 0 );
    // check if we got the correct child
    if ( &brother == node ) 
    {
        brother = parent->getChild( 1 );
    }
    
    // collect the possible reattachement points
    std::vector<TopologyNode*> new_brothers;
    findNewBrothers(new_brothers, *parent, &tau.getRoot());
    std::vector<double> weights = std::vector<double>(new_brothers.size(), 0.0);
    double sumOfWeights = 0.0;
    for (size_t i = 0; i<new_brothers.size(); ++i) 
    {
        // get the new brother
        TopologyNode* newBro = new_brothers[i];
        
        // do the proposal
        TopologyNode *newGrandparent = pruneAndRegraft(&brother, newBro, parent, grandparent);
        
        // flag for likelihood recomputation
        variable->touch();
        
        // compute the likelihood of the new value
        double priorRatio = variable->getLnProbability();
        double likelihoodRatio = 0.0;
        for (std::set<DagNode*>::const_iterator it = affected.begin(); it != affected.end(); ++it) 
        {
            likelihoodRatio += (*it)->getLnProbability();
        }
        weights[i] = exp(priorRatio + likelihoodRatio + offset);
        sumOfWeights += weights[i];
        
        // undo proposal
        pruneAndRegraft(newBro, &brother, parent, *newGrandparent);
        
        // restore the previous likelihoods;
        variable->restore();
    }
    
    if (sumOfWeights <= 1E-100) {
        // hack
        // the proposals have such a small likelihood that they can be neglected
//        throw new OperatorFailedException("Couldn't find another proposal with a decent likelihood.");
        return 0.0;
    }
    
    double ran = rng->uniform01() * sumOfWeights;
    size_t index = 0;
    while (ran > 0.0) {
        ran -= weights[index];
        index++;
    }
    index--;
    
    TopologyNode* newBro = new_brothers[index];
    
    // now we store all necessary values
    storedBrother       = &brother;
    storedNewBrother    = newBro;
    
    pruneAndRegraft(&brother, newBro, parent, grandparent);
    
    double forward = weights[index];
    
    double forwardProb = (forward / sumOfWeights);
    double backwardProb = (backward / (sumOfWeights - forward + backward));
    double hastingsRatio = log(backwardProb / forwardProb);
    
//.........这里部分代码省略.........