19.3 The TopoPriorCalculator class

(Mac version)

< 19.2 | 19.3 | 19.4 >

With the introduction of polytomous trees, the computation of the tree topology prior becomes more interesting but also much more complicated. We now need to not only know the number of binary tree topologies (the highest resolution class), we need to know the number of tree topologies in all the other resolution classes.

There are also at least two interesting variations on the concept of a flat prior with respect to tree topology. The prior could be flat across all tree topologies, or it could be flat across resolution classes. In the first case, the star tree and any particular fully-resolved tree topology have the same prior probability. In the second case, the star tree (which is the only representative of its resolution class), gets much greater weight than any particular fully-resolved tree topology. For example, with just 10 taxa, the star tree would have (in the flat resolution class prior) a prior probabilility 2,027,025 times that of any resolved tree!

To help us compute the topology prior, create the following class TopoPriorCalculator in a header file named topo_prior_calculator.hpp as follows:

#pragma once    

#include "tree.hpp"
#include "lot.hpp"

namespace strom {

    class PolytomyTopoPriorCalculator {
    
        public:
        
            typedef std::shared_ptr&lt;PolytomyTopoPriorCalculator&gt; SharedPtr;
                                            
                                            PolytomyTopoPriorCalculator();
            virtual                         ~PolytomyTopoPriorCalculator();

            bool                            isResolutionClassPrior() const;
            bool                            isPolytomyPrior() const;

            bool                            isRooted() const;
            bool                            isUnrooted() const;

            void                            setNTax(unsigned n);
            unsigned                        getNTax() const;

            void                            chooseRooted();
            void                            chooseUnrooted();

            double                          getLogCount(unsigned n, unsigned m);
            double                          getLogSaturatedCount(unsigned n);
            double                          getLogTotalCount(unsigned n);
            std::vector&lt;double&gt;             getLogCounts();

            std::vector&lt;double&gt;             getCountsVect();
            std::vector&lt;int&gt;                getNFactorsVect();

            void                            chooseResolutionClassPrior();
            void                            choosePolytomyPrior();

            void                            setC(double c);
            double                          getC() const;

            void                            setLogScalingFactor(double lnf);
            double                          getLogScalingFactor() const;

            virtual double                  getLogTopoProb(Tree::SharedPtr t);

            double                          getLogTopologyPrior(unsigned m);
            double                          getLogNormalizedTopologyPrior(unsigned m);
            double                          getLogNormConstant();

            std::vector&lt;double&gt;             getTopoPriorVect();
            std::vector&lt;double&gt;             getRealizedResClassPriorsVect();

            unsigned                        sample(Lot::SharedPtr rng);

            void                            clear();
            void                            reset();

        private:

            void                            recalcCountsAndPriorsImpl(unsigned n);
            void                            recalcPriorsImpl();

            unsigned                        _ntax;
            bool                            _is_rooted;
            bool                            _is_resolution_class_prior;
            double                          _C;

            bool                            _topo_priors_dirty;
            bool                            _counts_dirty;

            double                          _log_scaling_factor;
            std::vector&lt;int&gt;                _nfactors;
            std::vector&lt;double&gt;             _counts;
            double                          _log_total_count;
            std::vector&lt;double&gt;             _topology_prior;
    };
    
    // Member function bodies go here
    
}   

The TopoPriorCalculator constructor, destructor, and clear function

The constructor calls the clear member function to perform initialization, and the destructor does nothing, as usual.

    PolytomyTopoPriorCalculator::PolytomyTopoPriorCalculator() {    
        //std::cout &lt;&lt; "PolytomyTopoPriorCalculator being created" &lt;&lt; std::endl;
        clear();
    }   

    PolytomyTopoPriorCalculator::~PolytomyTopoPriorCalculator() {   
        //std::cout &lt;&lt; "PolytomyTopoPriorCalculator being destroyed" &lt;&lt; std::endl;
    }   

    inline void PolytomyTopoPriorCalculator::clear() {
        _topo_priors_dirty            = true;
        _is_rooted                    = false;
        _is_resolution_class_prior    = true;
        _C                            = 1.0;
        _ntax                         = 4;
        _counts_dirty                 = true;
        _log_scaling_factor           = 10.0;
    }   

The getLogNormalizedTopologyPrior member function

This is the function that delivers what you need: the log of the prior probability of a tree topology having m internal nodes. This represents the resolution class prior if _is_resolution_class_prior is true, otherwise it represents the polytomy prior.

A flat resolution class prior places equal prior probability mass on every possible value of m. For example, for N=7 taxa, the resolution class variable m can range from 1 (star tree) to N-2 = 5 (fully-resolved tree) for the unrooted tree case. The prior probability of each resolution class is thus 0.2 and a particular tree in each class equals 0.2 divided by the number of distinct labeled tree topologies in that class:

m ntopo ptree ntopo*ptree
1 1 0.2/1 = 0.200000000 0.2
2 56 0.2/56 = 0.003571429 0.2
3 490 0.2/490 = 0.000408163 0.2
4 1260 0.2/1260 = 0.00015873 0.2
5 945 0.2/945 = 0.00021164 0.2
total 2752   1.0

where m is the resolution class (number of internal nodes), ntopo is the number of tree topologies in resolution class m (computed in the recalcCountsAndPriorsImpl function), and ptree is the probability of any given labeled tree topology in resolution class m.

The flat polytomy prior is simpler: it assumes the same prior probability for each tree, regardless of the resolution class. Because the total number of tree topologies over all resolution class is 1+56+490+1260+945 = 2752 in the 7-taxon example above, the prior for each tree would be 1/2752 = 0.000363372 under the flat polytomy prior.

The priors described above were flat priors. It is possible to make some resolution classes more probable in a resolution class prior, and trees with more or fewer internal nodes more probable in a polytomy prior. For example, suppose you wished resolution class m = 1 (the star tree) to be twice as probable as the class of trees having 2 internal nodes (m = 2), class m = 2 to be twice as probable as class m = 3, and so on. The data member _C determines the relative probability of resolution class m over resolution class m+1. Setting _C = 2 yields this resolution class prior distribution:

m ntopo ptree ntopo*ptree
1 1 (16/31)/1 = 0.516129032 0.516129032
2 56 (8/31)/56 = 0.004608295 0.258064516
3 490 (4/31)/490 = 0.000263331 0.129032258
4 1260 (2/31)/1260 = 0.000051203 0.064516129
5 945 (1/31)/945 = 0.000034136 0.032258065
total 2752   1.0

The data member _C in the case of a polytomy prior ensures that the probability of any particular tree in resolution class m is _C times more probable than any particular tree in resolution class m+1. To see how this prior is calculated, note that have these 5 constraints (again for the 7-taxon unrooted case):

n_1 p_1 + n_2 p_2 + n_3 p_3 + n_4 p_4 + n_5 p_5 = 1.0
p_1 = C p_2
p_2 = C p_3
p_3 = C p_4
p_4 = C p_5

These constraints leave only p_5 to be determined:

p_1 = C^4 p_5
p_2 = C^3 p_5
p_3 = C^2 p_5
p_4 = C^1 p_5
(n_1 C^4 + n_2 C^3 + n_3 C^2 + n_4 C^1 + n_5 C^0) p_5 = 1.0
p_5 = 1/[n_1 C^4 + n_2 C^3 + n_3 C^2 + n_4 C^1 + n_5 C^0]

For _C = 2 and a polytomy prior, we have

m ntopo qtree ptree ntopo*ptree
1 1 1*2^4 = 16 2^4/5889 = 0.00271693 1*2^4/5889 = 0.00271693
2 56 56*2^3 = 448 2^3/5889 = 0.001358465 56*2^3/5889 = 0.076074036
3 490 490*2^2 = 1960 2^2/5889 = 0.000679232 490*2^2/5889 = 0.332823909
4 1260 1260*2^1 = 2520 2^1/5889 = 0.000339616 1260*2^1/5889 = 0.427916454
5 945 945*2^0 = 945 2^0/5889 = 0.000169808 945*2^0/5889 = 0.16046867
total 2752 16+448+1960+2520+945 = 5889   1.0

where qtree is the unnormalized probability of a particular tree topology (i.e. ptree is qtree divided by the sum of qtree values).

With this introduction, here is the deceptively simple getLogNormalizedTopologyPrior function. The hard part is computing the _topology_prior vector, which is described below in the explanation of the recalcPriorsImpl function. The log of the normalized topology prior is obtained as _topology_prior[m] minus _topology_prior[0] (the 0th element of _topology_prior holds the log of the normalization constant).

    double PolytomyTopoPriorCalculator::getLogNormalizedTopologyPrior(unsigned m) { 
        if (_topo_priors_dirty)
            reset();
        assert(m &lt; _topology_prior.size());
        return (_topology_prior[m] - _topology_prior[0]);
    }

The recalcPriorsImpl member function

This function is responsible for recomputing the _topology_prior vector. If _is_resolution_class_prior is true, this function requires knowledge of the number of tree topologies in each resolution class, so it is called inside the recalcCountsAndPriorsImpl after the counts are calculated.

    void PolytomyTopoPriorCalculator::recalcPriorsImpl() {  
        _topology_prior.clear();
        _topology_prior.push_back(0.0);    // This will hold the normalizing constant in the end

        // Figure out the maximum possible value for m, the number of internal nodes
        unsigned maxm = _ntax - (_is_rooted ? 1 : 2);

        if (_is_resolution_class_prior) {
            // _counts vector should have length equal to maxm if everything is ok
            assert(maxm == (unsigned)_counts.size());

            double logC = std::log(_C);
            double log_normalization_constant = 0.0;
            if (_C == 1.0)
                log_normalization_constant = std::log(maxm);
            else if (_C &gt; 1.0) {
                // factor out largest value to avoid overflow
                double term1  = logC*maxm;
                double term2a = std::exp(-logC*maxm);
                double term2b = std::log(1.0 - term2a);
                double term3  = std::log(_C - 1.0);
                log_normalization_constant = term1 + term2b + term3;
            }
            else {
                log_normalization_constant = std::log(1.0 - std::exp(logC*maxm)) - std::log(1.0 - _C);
            }
            _topology_prior[0] = log_normalization_constant;
            for (unsigned m = 1; m &lt;= maxm; ++m) {
                double logCterm = (double)(maxm - m)*logC;
                double log_count_m = std::log(_counts[m - 1]) + _log_scaling_factor*(double)_nfactors[m - 1];
                double log_v = logCterm - log_count_m;
                _topology_prior.push_back(log_v);
            }
        }
        else {
            double total = 0.0;
            double logC = std::log(_C);
            for (unsigned m = 1; m &lt;= maxm; ++m) {
                double logCterm = (double)(maxm - m)*logC;
                total += std::exp(logCterm);
                _topology_prior.push_back(logCterm);
            }
            _topology_prior[0] = std::log(total);
        }
    }   

The recalcCountsAndPriorsImpl member function

This function recomputes the _counts vector for the supplied number of internal nodes (n) using the method outlined by Joe Felsenstein in his 2004 book and also in Felsenstein (1978). Below I’ve reproduced the table from Felsenstein (2004) showing the number of possible rooted tree topologies for any number of internal nodes (rows) and up to 8 taxa (columns):

nodes 2 3 4 5 6 7 8 nfactors
1 1 1 1 1 1 1 1 0
2   3 10 25 56 119 246 0
3     15 105 490 1918 6825 0
4       0.0105 0.1260 0.945 5.6980 1
5         0.0945 1.7325 19.0575 1
6           1.0395 27.0270 1
7             13.5135 1

The recalcCountsAndPriorsImpl function works from left to right, calculating each column in turn. Within a column, it works down. Each cell except the first in a given column requires knowledge of two cells from the column to its left: the cell to its immediate left as well as the cell above the cell to its immediate left. This is because in order to know how many tree topologies there are for N taxa, one needs to know how many places a taxon could be added to a tree with N-1 taxa. The N-1 taxon tree could have the same number of internal nodes as the new tree (the new taxon was added to an existing node, either creating a new polytomy or enlarging an existing one), or the N-1 taxon tree could have one fewer internal nodes (in which case the new taxon inserts a new node).

Note that this function calls recalcPriorsImpl before returning.

    void PolytomyTopoPriorCalculator::recalcCountsAndPriorsImpl(unsigned n) {   
        if (_is_resolution_class_prior)
            _counts_dirty = true;
        if (_counts_dirty) {
            double scaling_factor = exp(_log_scaling_factor);

            _counts.clear();
            _counts.push_back(1.0); // _counts are always 1 for m = 1

            int last_factor = 0;
            _nfactors.clear();
            _nfactors.push_back(last_factor); // never need to scale this one

            // temporary variables
            double a, b, c;
            double max_log_double = log(DBL_MAX);
            
            // The value of epsilon is arbitrary, but must be larger than
            // zero and less than scaling_factor
            double epsilon = scaling_factor/10.0;

            // Compute the vector of _counts for the number of internal nodes specified
            // This is the main loop over columns. z is the number of taxa minus 1
            // for rooted trees, and is the number of taxa minus 2 for unrooted trees
            for (unsigned z = 2; z &lt;= n; ++z) {
                // This column is one element longer than the column to its left
                _counts.push_back(0.0);
                _nfactors.push_back(last_factor);
                
                // _counts[0] is always 1.0 because there is only one star tree topology
                b = _counts[0];

                // This is the loop over rows within the current column.
                // m + 1 is the number of internal nodes.
                for (unsigned m = 1; m &lt; z; ++m) {
                    unsigned num_internal_nodes = m + 1;
                    double diff = (double)(_nfactors[m - 1] - _nfactors[m]);
                    double log_factor = diff*_log_scaling_factor;
                    assert(log_factor &lt; max_log_double);
                    a = b*exp(log_factor);
                    b = _counts[m];
                    c = a*((double)(z + num_internal_nodes - 1));
                    if (num_internal_nodes &lt; z) {
                        c += b*(double)num_internal_nodes;
                    }
                    if (c &gt; scaling_factor) {
                        double incr = floor(log(c - epsilon)/_log_scaling_factor);
                        _nfactors[m] += (unsigned)incr;
                        last_factor = _nfactors[m];
                        _counts[m] = exp(log(c) - incr*_log_scaling_factor);
                        b = exp(log(b) - incr*_log_scaling_factor);
                    }
                    else
                        _counts[m] = c;
                }
            }

            // Now compute the log of the total number of tree topologies
            // over all possible resolution classes (i.e. number of internal nodes)
            // Begin by creating a vector of log _counts and finding the
            // largest value (this will be factored out to avoid overflow)
            std::vector&lt;double&gt; v;
            unsigned sz = (unsigned)_nfactors.size();
            assert(sz == _counts.size());
            double max_log_count = 0.0;
            for (unsigned i = 0; i &lt; sz; ++i) {
                double num_factors = (double)_nfactors[i];
                double log_count = num_factors*_log_scaling_factor + log(_counts[i]);
                if (log_count &gt; max_log_count)
                    max_log_count = log_count;
                v.push_back(log_count);
            }

            // Compute log sum of _counts by factoring out the largest count.
            // Underflow will occur, but only for _counts that are so much
            // smaller than the dominant _counts that the underflow can be
            // ignored for all practical purposes
            double sum = 0.0;
            for (std::vector&lt;double&gt;::const_iterator it = v.begin(); it != v.end(); ++it) {
                double diff = (*it) - max_log_count;
                sum += exp(diff);
            }
            assert(sum &gt; 0.0);
            _log_total_count = log(sum) + max_log_count;
            _counts_dirty = false;
        }
        else {
            _nfactors.clear();
            _counts.clear();
            
            // _counts_dirty ensures that _counts will be calculated if requested,
            // for example by calling getCount
            _counts_dirty = true;
        }

        // Recalculate the _topology_prior vector too
        recalcPriorsImpl();

        _topo_priors_dirty = false;
    }   

The reset member function

This function forces recalculation of polytomy_prior if _is_resolution_class_prior is false, and both _counts and polytomy_prior if _is_resolution_class_prior is true (or if _counts_dirty is true).

    void PolytomyTopoPriorCalculator::reset() { 
        unsigned num_internal_nodes = (_is_rooted ? (_ntax - 1) : (_ntax - 2));
        recalcCountsAndPriorsImpl(num_internal_nodes);
    } 

The getLogCount member function

This function returns the count of the number of distinct labeled tree topologies for n taxa and m internal nodes. It depends on the _counts vector being accurate, so it calls recalcCountsAndPriors function if n is not equal to _ntax. If _is_rooted is true, assumes m is less than _ntax. If _is_rooted is false, assumes m less than _ntax - 1.

    double PolytomyTopoPriorCalculator::getLogCount(unsigned n, unsigned m) {   
        assert((_is_rooted && (m &lt; n)) || (!_is_rooted && (m &lt; n - 1)));
        if (n != _ntax)
            setNTax(n);
        if (_counts_dirty)
            reset();
        double nf = (double)(_nfactors[m - 1]);
        double log_count = nf*_log_scaling_factor + log(_counts[m - 1]);
        return log_count;
    }   

The getLogSaturatedCount member function

Returns the number of saturated (i.e. fully-resolved and thus having as many internal nodes as possible) trees of n taxa. Calls recalcCountsAndPriors function if n is not equal to _ntax because that means that the _counts vector is no longer accurate and must be recalculated.

    double PolytomyTopoPriorCalculator::getLogSaturatedCount(unsigned n) {  
        if (n != _ntax)
            setNTax(n);
        if (_counts_dirty)
            reset();
        unsigned last = (unsigned)(_counts.size() - 1);
        double nf = (double)(_nfactors[last]);
        double log_count = nf*_log_scaling_factor + log(_counts[last]);
        return log_count;
    }   

The getLogCounts member function

This function constructs a vector in which the element having index m (i = 0, 1, …, maximum number of internal nodes) represents the natural logarithm of the number of tree topologies having m internal nodes. If _counts_dirty is true, it recomputes the _counts vectors first. The 0th element of the returned vector holds the natural log of the total number of tree topologies (log of the sum of all other elements).

    std::vector&lt;double&gt; PolytomyTopoPriorCalculator::getLogCounts() {   
        if (_is_resolution_class_prior)
            _counts_dirty = true;
        if (_counts_dirty)
            reset();

        std::vector&lt;double&gt; v;
        unsigned sz = _ntax - (_is_rooted ? 0 : 1);
        v.reserve(sz);
        v.push_back(_log_total_count);

        for (unsigned i = 1; i &lt; sz; ++i) {
            double log_Tnm = log(_counts[i - 1]) + _log_scaling_factor*(double)(_nfactors[i - 1]);
            v.push_back(log_Tnm);
        }

        return v;
    }   

The getLogTotalCount member function

This function returns the natural log of the total number of trees for n taxa, including all resolution classes from the star tree to fully resolved (saturated) trees. Calls recalcCountsAndPriors function if n is not equal to _ntax or if not using the resolution class prior (in which case _counts have not been calculated).

    double PolytomyTopoPriorCalculator::getLogTotalCount(unsigned n) {  
        if (n != _ntax)
            setNTax(n);
        if (_counts_dirty)
            reset();
        return _log_total_count;
    }   

The getRealizedResClassPriorsVect member function

Constructs a vector of realized resolution class priors from the values in the _topology_prior vector. If _topo_priors_dirty is true, it recomputes the _topology_prior vectors first. The mth element of the returned vector is set to T_{n,m}*_topology_prior[m] for m > 0. The 0th element of the returned vector holds the normalization constant (sum of all other elements). This function is not efficient because it is intended only to be used for providing information to the user on request. Table 2, p. 248, in the Lewis, Holder, and Holsinger (2005) presented (normalized) values from this vector, as do the columns labeled ntopo*ptree in the tables presented above in the documentation for the getLogNormalizedTopologyPrior member function.

    std::vector&lt;double&gt; PolytomyTopoPriorCalculator::getRealizedResClassPriorsVect() {  
        if (!_is_resolution_class_prior)
            _counts_dirty = true;
        if (_topo_priors_dirty || _counts_dirty)
            reset();

        std::vector&lt;double&gt; v;
        v.reserve(_topology_prior.size());
        v.push_back(0.0);

        unsigned sz = (unsigned)_topology_prior.size();

        // First loop will be to determine largest value, which will be factored out
        // the second time through so that the total does not overflow
        double log_factored_out = 0.0;
        for (unsigned i = 1; i &lt; sz; ++i) {
            double c = _counts[i - 1];
            double nf = (double)_nfactors[i - 1];
            double log_Tnm = log(c) + _log_scaling_factor*nf;
            double log_prior = log_Tnm + _topology_prior[i];
            v.push_back(log_prior);
            if (log_prior &gt; log_factored_out)
                log_factored_out = log_prior;
        }

        // Now we can compute the total
        double total = 0.0;
        std::vector&lt;double&gt;::const_iterator it = v.begin();
        for (++it; it != v.end(); ++it) {
            total += exp((*it) - log_factored_out);
        }
        v[0] = log(total) + log_factored_out;

        return v;
    }   

The sample member function

Samples a resolution class (i.e. number of internal nodes) from the realized resolution class distribution. This function is not particularly efficient because it calls PolytomyTopoPriorCalculator::getRealizedResClassPriorsVect, resulting in an unnecessary vector copy operation.

    unsigned PolytomyTopoPriorCalculator::sample(Lot::SharedPtr rng) {  
        std::vector&lt;double&gt; v = getRealizedResClassPriorsVect();
        double u = rng-&gt;uniform();
        double z = v[0];
        double cum = 0.0;
        for (unsigned i = 1; i &lt; v.size(); ++i) {
            cum += exp(v[i] - z);
            if (u &lt;= cum)
                return i;
        }
        assert(0);
        return (unsigned)(v.size() - 1);
    }   

The setNTax setter and getNTax accessor

The setNTax function returns immediately if _ntax equals new_ntax; however, if new_ntax differs from _ntax, setNTax sets _ntax to new_ntax and sets _topo_priors_dirty to true so that subsequent requests for prior probabilities or counts will trigger recalculation. The getNTax function returns the value of the data member _ntax.

    void PolytomyTopoPriorCalculator::setNTax(unsigned new_ntax) {  
        if (_ntax != new_ntax) {
            // Set _ntax to the new value
            assert(new_ntax &gt; (unsigned)(_is_rooted ? 1 : 2));
            _ntax = new_ntax;

            _counts_dirty = true;
            _topo_priors_dirty = true;
        }
    }   

    unsigned PolytomyTopoPriorCalculator::getNTax() const { 
        return _ntax;
    }  

The chooseRooted and chooseUnrooted member functions

The chooseRooted function sets the _is_rooted data member to true. There are more rooted than unrooted trees for the same value of _ntax, so this setting is important when asking questions that require knowledge of the numbers of possible trees. The chooseUnrooted function sets the _is_rooted data member to false.

    void PolytomyTopoPriorCalculator::chooseRooted() {  
        if (!_is_rooted) {
            _is_rooted = true;
            _topo_priors_dirty = true;
        }
    }   

    void PolytomyTopoPriorCalculator::chooseUnrooted() {    
        if (_is_rooted) {
            _is_rooted = false;
            _topo_priors_dirty = true;
        }
    }   

The chooseResolutionClassPrior and choosePolytomyPrior member functions

These two functions can be used to choose between a resolution class prior or the alternative (the polytomy prior).

    void PolytomyTopoPriorCalculator::chooseResolutionClassPrior() {    
        if (!_is_resolution_class_prior) {
            _is_resolution_class_prior = true;
            _topo_priors_dirty = true;
        }
    }   

    void PolytomyTopoPriorCalculator::choosePolytomyPrior() {   
        if (_is_resolution_class_prior) {
            _is_resolution_class_prior = false;
            _topo_priors_dirty = true;
        }
    }   

The setC and getC member functions

These two functions either set or return the value of the data member _C, which determines the nature of the resolution class or polytomy prior. For the resolution class prior, each resolution class m has probability _C times greater than resolution class m+1. For the polytomy prior, a particular tree with m internal nodes has probability _C times greater than a particular tree with m+1 internal nodes.

    void PolytomyTopoPriorCalculator::setC(double c) {  
        assert(c &gt; 0.0);
        if (c != _C) {
            _C = c;
            _topo_priors_dirty = true;
        }
    }   

    double PolytomyTopoPriorCalculator::getC() const {  
        return _C;
    }   

The setLogScalingFactor and getLogScalingFactor member functions

These two functions set or return the value of the data member _log_scaling_factor data member.

    void PolytomyTopoPriorCalculator::setLogScalingFactor(double lnf) { 
        assert(lnf &gt; 0.0);
        if (lnf != _log_scaling_factor) {
            _log_scaling_factor = lnf;
            _counts_dirty = true;
            _topo_priors_dirty = true;
        }
    }   

    double PolytomyTopoPriorCalculator::getLogScalingFactor() const {   
        return _log_scaling_factor;
    }   

The getCountsVect member function

Returns copy of the _counts vector, which contains in its (m-1)th element the number of tree topologies having exactly m internal nodes. Note that you will need to also call getNFactorsVect if there is any chance that some of the _counts are larger than exp(_log_scaling_factor). In such cases, the actual log count equals _nfactors[m]*_log_scaling_factor + log(count[m - 1].

    std::vector&lt;double&gt; PolytomyTopoPriorCalculator::getCountsVect() {  
        if (_counts_dirty)
            reset();
        return _counts;
    }   

The getNFactorsVect member function

Returns copy of the _nfactors vector, which contains in its (m-1)th element the number of times _counts[m] has been rescaled by dividing by the scaling factor (the log of which is _log_scaling_factor).

    std::vector&lt;int&gt; PolytomyTopoPriorCalculator::getNFactorsVect() {   
        if (_counts_dirty)
            reset();
        return _nfactors;
    }   

The getTopoPriorVect member function

Returns copy of the _topology_prior vector, which contains in its mth element the unnormalized prior for tree topologies having exactly m internal nodes. The 0th element of _topology_prior holds the normalizing constant.

    std::vector&lt;double&gt; PolytomyTopoPriorCalculator::getTopoPriorVect() {   
        if (_topo_priors_dirty)
            reset();
        return _topology_prior;
    }   

The getLogTopoProb member function

Returns result of call to getLogTopologyPrior(m), where m is the number of internal nodes in the supplied tree t.

    double PolytomyTopoPriorCalculator::getLogTopoProb(Tree::SharedPtr t) { 
        unsigned n = t-&gt;numLeaves();
        assert(n &gt; 3);
        if (n != getNTax())
            setNTax(n);
        unsigned m = t-&gt;numInternals();
        double topo_prior = getLogTopologyPrior(m);
        return topo_prior;
    }   

The getLogTopologyPrior member function

Returns the natural logarithm of the unnormalized topology prior. This represents the resolution class prior if _is_resolution_class_prior is true, otherwise it represents the polytomy prior.

    double PolytomyTopoPriorCalculator::getLogTopologyPrior(unsigned m) {   
        if (_topo_priors_dirty)
            reset();
        assert(m &lt; _topology_prior.size());
        return _topology_prior[m];
    }   

The getLogNormConstant member function

Returns the natural logarithm of the normalizing constant for the topology prior. This value is stored in _topology_prior[0].

    double PolytomyTopoPriorCalculator::getLogNormConstant() {  
        if (_topo_priors_dirty)
            reset();
        return _topology_prior[0];
    }   

The isResolutionClassPrior and isPolytomyPrior member functions

The isResolutionClassPrior function returns the current value of the data member _is_resolution_class_prior and isPolytomyPrior returns !isResolutionClassPrior().

    bool PolytomyTopoPriorCalculator::isResolutionClassPrior() const {  
        return _is_resolution_class_prior;
    }   

    bool PolytomyTopoPriorCalculator::isPolytomyPrior() const { 
        return !_is_resolution_class_prior;
    }   

The isRooted and isUnrooted member functions

The isRooted function returns the current value of the data member _is_rooted, and isUnrooted returns !isRooted().

    bool PolytomyTopoPriorCalculator::isRooted() const {    
        return _is_rooted;
    }   

    bool PolytomyTopoPriorCalculator::isUnrooted() const {  
        return !_is_rooted;
    }   

Literature Cited

Felsenstein, J. 1978. The number of evolutionary trees. Systematic Zoology 27:27-33 (see also the 1981 correction Systematic Zoology 30:122). https://doi.org/10.2307/2412810/

Felsenstein, J. 2004. Inferring Phylogeny. Sinauer, Sunderland, Massachusetts. ISBN: 9780878931774

Lewis, P. O., Holder, M. T., and Holsinger, K. E. 2005. Polytomies and bayesian phylogenetic inference. Systematic Biology 54(2):241–253. https://doi.org/10.1080/10635150590924208

< 19.2 | 19.3 | 19.4 >