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

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

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

size_t SpeciesTreeNodeSlideProposal::fillPreorderIndices(TopologyNode &node, size_t loc, std::vector<size_t> &indices)
{
    if ( node.isInternal() == true )
    {
        size_t l = fillPreorderIndices(node.getChild( 0 ), loc, indices);
        indices[node.getIndex()] = l;
        loc = fillPreorderIndices(node.getChild( 1 ), l + 1, indices);
    }
    
    return loc;
}

void MultivariateBrownianPhyloProcess::recursiveCorruptAll(const TopologyNode& from)    {
    
    dirtyNodes[from.getIndex()] = true;
    for (size_t i = 0; i < from.getNumberOfChildren(); ++i) {
        recursiveCorruptAll(from.getChild(i));
    }    
}

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

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

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 AutocorrelatedLognormalRateBranchwiseVarDistribution::recursiveSimulate(const TopologyNode& node, double parentRate) {
    
    // get the index
    size_t nodeIndex = node.getIndex();
    
    // compute the variance along the branch
	double scale = scaleValue->getValue();
    double variance = sigma->getValue()[nodeIndex] * node.getBranchLength() * scale;
	double mu = log(parentRate) - (variance * 0.5);
	double stDev = sqrt(variance);
    
    // simulate a new rate
    RandomNumberGenerator* rng = GLOBAL_RNG;
    double nodeRate = RbStatistics::Lognormal::rv( mu, stDev, *rng );
    
    // we store this rate here
    (*value)[nodeIndex] = nodeRate;
    
    // simulate the rate for each child (if any)
    size_t numChildren = node.getNumberOfChildren();
    for (size_t i = 0; i < numChildren; ++i) {
        const TopologyNode& child = node.getChild(i);
        recursiveSimulate(child,nodeRate);
    }
}

double AutocorrelatedLognormalRateBranchwiseVarDistribution::recursiveLnProb( const TopologyNode& n ) {
    
    // get the index
    size_t nodeIndex = n.getIndex();
    
    double lnProb = 0.0;
    size_t numChildren = n.getNumberOfChildren();
	double scale = scaleValue->getValue();
    
    if ( numChildren > 0 ) {
        double parentRate = log( (*value)[nodeIndex] );
        
        for (size_t i = 0; i < numChildren; ++i) {
            const TopologyNode& child = n.getChild(i);
            lnProb += recursiveLnProb(child);
            
            size_t childIndex = child.getIndex();
            // compute the variance
            double variance = sigma->getValue()[childIndex] * child.getBranchLength() * scale;
            
            double childRate = (*value)[childIndex];

			// the mean of the LN dist is parentRate = exp[mu + (variance / 2)],
			// where mu is the location param of the LN dist (see Kishino & Thorne 2001)
			double mu = parentRate - (variance * 0.5);
			double stDev = sqrt(variance);
            lnProb += RbStatistics::Lognormal::lnPdf(mu, stDev, childRate);
        } 
    }
    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));
    }
    
}

TopologyNode* SpeciesTreeNodeSlideProposal::mauReconstructSub(Tree &tree, size_t from, size_t to, std::vector<TopologyNode*> &order, std::vector<bool>&wasSwaped)
{
    if( from == to )
    {
        return order[2*from];
    }
    
    size_t node_index = -1;
    {
        double h = -1;
        
        for(size_t i = from; i < to; ++i)
        {
            double v = order[2 * i + 1]->getAge();
            if( h < v )
            {
                h = v;
                node_index = i;
            }
        }
    }
    
    TopologyNode* node = order[2 * node_index + 1];
    
    TopologyNode& lchild = node->getChild( 0 );
    TopologyNode& rchild = node->getChild( 1 );
    
    TopologyNode* lTargetChild = mauReconstructSub(tree, from, node_index, order, wasSwaped);
    TopologyNode* rTargetChild = mauReconstructSub(tree, node_index+1, to, order, wasSwaped);
    
    if( wasSwaped[node_index] )
    {
        TopologyNode* z = lTargetChild;
        lTargetChild = rTargetChild;
        rTargetChild = z;
    }
    
    if( &lchild != lTargetChild )
    {
        // replace the left child
        node->removeChild( &lchild );
        node->addChild( lTargetChild );
    }
    
    if( &rchild != rTargetChild )
    {
        // replace the right child
        node->removeChild( &rchild );
        node->addChild( rTargetChild );
    }
    
    return node;
}

/**
 * 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::recursiveGetStats(const TopologyNode& from, double& e1, double& e2, int& n) const {

    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) {
        recursiveGetStats(from.getChild(i),e1,e2,n);
    }
    
}

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