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

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

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

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

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

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

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 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 AutocorrelatedBranchMatrixDistribution::recursiveLnProb( const TopologyNode& n ) {
    
    // get the index
    size_t nodeIndex = n.getIndex();
    
    double lnProb = 0.0;
    size_t numChildren = n.getNumberOfChildren();
    
    if ( numChildren > 0 ) {
        std::vector<double> parent = (*value)[nodeIndex].getStationaryFrequencies();

        std::vector<double>::iterator end = parent.end();
        for (std::vector<double>::iterator it = parent.begin(); it != end; ++it) {
            (*it) *= alpha->getValue();
        }

        for (size_t i = 0; i < numChildren; ++i) {
            const TopologyNode& child = n.getChild(i);
            lnProb += recursiveLnProb(child);
            
            size_t childIndex = child.getIndex();
            //        RateMatrix& rm = (*value)[childIndex];
            
            // compare if the child has a different matrix
            if ( matrixIndex[nodeIndex] == matrixIndex[childIndex] ) {
                // no change -> just the probability of no change
                lnProb += log( 1.0 - changeProbability->getValue() );
            } else {
                // change:
                
                // probability of change
                lnProb += log( changeProbability->getValue() );
                
                const std::vector<double>& descendant = (*value)[childIndex].getStationaryFrequencies();
        //            const std::vector<double>& descendant = uniqueMatrices[ matrixIndex[childIndex] ].getStationaryFrequencies();
                
                // probability of new descendant values
                double p = RbStatistics::Dirichlet::lnPdf(parent, descendant);
                lnProb += p;
            }
        } 
    }
    
    return lnProb;
    
}

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

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