double RateMap_Biogeography::getRate(const TopologyNode& node, std::vector<CharacterEvent*> from, CharacterEvent* to, unsigned* count, double age) const
{
    double rate = 0.0;
    int s = to->getState();
    
    if (from[ to->getIndex() ]->getState() == to->getState())
    {
        std::cout << count[0] << " " << count[1] << "\n";
        std::cout << node.getIndex() << " problem...\n";
        ;
    }
    
    // rate to extinction cfg is 0
    if (count[1] == 1 && s == 0 && forbidExtinction)
        return 0.0;
    
    // rate according to binary rate matrix Q(node)
    if (branchHeterogeneousGainLossRates)
        rate = heterogeneousGainLossRates[node.getIndex()][s];
    else
        rate = homogeneousGainLossRates[s];
    
    if (branchHeterogeneousClockRates)
        rate *= heterogeneousClockRates[node.getIndex()];
    else
        rate *= homogeneousClockRate;
    
    // apply rate modifiers
    if (useGeographyRateModifier) // want this to take in age as an argument...
        rate *= geographyRateModifier->computeRateModifier(node, from, to, age);
    
    return rate;

}

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 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 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 MultivariateBrownianPhyloProcess::recursiveCorruptAll(const TopologyNode& from)    {
    
    dirtyNodes[from.getIndex()] = true;
    for (size_t i = 0; i < from.getNumberOfChildren(); ++i) {
        recursiveCorruptAll(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 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;
        }
        
    }
    
}

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

double RateMap_Biogeography::getSiteRate(const TopologyNode& node, CharacterEvent* from, CharacterEvent* to, double age) const
{
    double rate = 0.0;
    int s = to->getState();
//    int charIdx = to->getIndex();
//    int epochIdx = getEpochIndex(age);
    
    // rate according to binary rate matrix Q(node)
    if (branchHeterogeneousGainLossRates)
        rate = heterogeneousGainLossRates[node.getIndex()][s];
    else
        rate = homogeneousGainLossRates[s];
    
    if (branchHeterogeneousClockRates)
        rate *= heterogeneousClockRates[node.getIndex()];
    else
        rate *= homogeneousClockRate;
    
    // area effects
    if (useGeographyRateModifier)
        rate *= geographyRateModifier->computeSiteRateModifier(node,from,to,age);

    
    return rate;
}

void SpeciesNarrowExchangeProposal::regraft( TopologyNode *node, TopologyNode *child )
{
    
    // regraft
    TopologyNode* newParent = &child->getParent();
    newParent->removeChild( child );
    newParent->addChild( node );
    node->setParent( newParent );
    node->addChild( child );
    child->setParent( node );
}

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 StitchTreeFunction::recursivelyStitchPatchClades(TopologyNode* node, size_t& index)
{
    
    // stich patch clade to matching tip taxon
    if (node->isTip())
    {
        for (size_t i = 0; i < patchTaxa.size(); i++)
        {
            if (node->getTaxon() == stitchTaxon[i])
            {
                // remove the node
                TopologyNode* parent = &node->getParent();
                parent->removeChild(node);
                node->setParent(NULL);
                delete node;
                
                // add the patch clade
                const Tree& t = patchClades->getValue()[i];
                
                // this memory is freed when the stitch tree is deleted in updateStitchTree()
                TopologyNode* patchRoot = new TopologyNode( t.getRoot() );
                
                // prune out non-patch taxa
                std::set<Taxon> remainingTaxa = prunedTaxa[i];
                recursivelyCleanPatchClade(patchRoot, patchRoot, remainingTaxa, index, i);
                recursivelyIndexPatchClade(patchRoot, index, i);

                // add patch clade to base tree
                parent->addChild(patchRoot);
                patchRoot->setParent(parent);
                return;
            }
        }
    }
    
    // recurse towards tips
    std::vector<TopologyNode*> children = node->getChildren();
    for (size_t i = 0; i < children.size(); i++)
    {
        recursivelyStitchPatchClades(children[i], index);
    }
    
    // assign index for first update only
    if (!haveIndex) {
        stitchTreeIndex[numPatches][node->getIndex()] = index++;
    }
    
    // set index as needed
    size_t old_index = node->getIndex();
    node->setIndex( stitchTreeIndex[numPatches][old_index] );
    
    return;
}

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

void SpeciesNarrowExchangeProposal::prune( TopologyNode *node, TopologyNode *child )
{
    // get the parent and brother node
    TopologyNode &parent = node->getParent();
    TopologyNode *brother = &node->getChild( 0 );
    if ( brother == child )
    {
        brother = &node->getChild( 1 );
    }
    // prune
    parent.removeChild( node );
    node->removeChild( brother );
    parent.addChild( brother );
    brother->setParent( &parent );
    
}

double RateMap_Biogeography::getUnnormalizedSumOfRates(const TopologyNode& node, std::vector<CharacterEvent*> from, unsigned* counts, double age) const
{
    size_t nodeIndex = node.getIndex();
    size_t epochIdx = getEpochIndex(age);
    
    // apply ctmc for branch
    const std::vector<double>& glr = ( branchHeterogeneousGainLossRates ? heterogeneousGainLossRates[nodeIndex] : homogeneousGainLossRates );
    
    // get sum of rates
    double sum = 0.0;
    for (size_t i = 0; i < from.size(); i++)
    {
        unsigned s = from[i]->getState();
        double v = availableAreaVector[ epochIdx * this->numCharacters + i ];
        
        if (forbidExtinction && s == 1 && counts[1] == 1)
            sum += 0.0;
        else if (s == 1 && v > 0)
            sum += glr[0];
        else if (s == 1 && v == 0)
            sum += 10e10;
        else  if (s == 0)
            sum += glr[1] * v;
    }
    
    // apply rate for branch
    if (branchHeterogeneousClockRates)
        sum *= heterogeneousClockRates[nodeIndex];
    else
        sum *= homogeneousClockRate;
    
    return sum;
}

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

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

/** Build random binary tree to size num_taxa. The result is a draw from the uniform distribution on histories. */
void HeterogeneousRateBirthDeath::buildRandomBinaryHistory(std::vector<TopologyNode*> &tips)
{
    
    if (tips.size() < num_taxa)
    {
        // Get the rng
        RandomNumberGenerator* rng = GLOBAL_RNG;
        
        // Randomly draw one node from the list of tips
        size_t index = static_cast<size_t>( floor(rng->uniform01()*tips.size()) );
        
        // Get the node from the list
        TopologyNode* parent = tips.at(index);
        
        // Remove the randomly drawn node from the list
        tips.erase(tips.begin()+long(index));
        
        // Add a left child
        TopologyNode* leftChild = new TopologyNode(0);
        parent->addChild(leftChild);
        leftChild->setParent(parent);
        tips.push_back(leftChild);
        
        // Add a right child
        TopologyNode* rightChild = new TopologyNode(0);
        parent->addChild(rightChild);
        rightChild->setParent(parent);
        tips.push_back(rightChild);
        
        // Recursive call to this function
        buildRandomBinaryHistory(tips);
    }
}