19.4 Adding calcLogTopologyPrior to Updater

(Win version)

< 19.3 | 19.4 | 19.5 >

Updating the Updater class declaration

The following highlighted changes need to be made to the class declaration in updater.hpp:

#pragma once    

#include "tree.hpp"
#include "tree_manip.hpp"
#include "lot.hpp"
#include "xstrom.hpp"
#include "likelihood.hpp"
<span style="color:#0000ff"><strong>#include "topo_prior_calculator.hpp"</strong></span>

namespace strom {
    class Chain;

    class Updater {
    
        friend class Chain;

        public:
            typedef std::shared_ptr&lt;Updater&gt;        SharedPtr;
        
            TreeManip::SharedPtr                    getTreeManip() const;

                                                    Updater();
            virtual                                 ~Updater();

            void                                    setLikelihood(Likelihood::SharedPtr likelihood);
            void                                    setTreeManip(TreeManip::SharedPtr treemanip);
            void                                    setLot(Lot::SharedPtr lot);
            void                                    setLambda(double lambda);
            void                                    setHeatingPower(double p);
            void                                    setTuning(bool on);
            void                                    setTargetAcceptanceRate(double target);
            void                                    setPriorParameters(const std::vector&lt;double&gt; & c);
            <span style="color:#0000ff"><strong>void                                    setTopologyPriorOptions(bool resclass, double C);</strong></span>
            void                                    setWeight(double w);
            void                                    calcProb(double wsum);

            double                                  getLambda() const;
            double                                  getWeight() const;
            double                                  getProb() const;
            double                                  getAcceptPct() const;
            double                                  getNumUpdates() const;
            std::string                             getUpdaterName() const;

            virtual void                            clear();

            virtual double                          calcLogPrior() = 0;
            double                                  calcLogTopologyPrior() const;
            double                                  calcLogEdgeLengthPrior() const; 
            double                                  calcLogLikelihood() const;
            virtual double                          update(double prev_lnL);

            static double                           getLogZero();
            
        protected:

            virtual void                            reset();
            virtual void                            tune(bool accepted);

            virtual void                            revert() = 0;
            virtual void                            proposeNewState() = 0;

            Lot::SharedPtr                          _lot;
            Likelihood::SharedPtr                   _likelihood;
            TreeManip::SharedPtr                    _tree_manipulator;
            std::string                             _name;
            double                                  _weight;
            double                                  _prob;
            double                                  _lambda;
            double                                  _log_hastings_ratio;
            double                                  _log_jacobian;
            double                                  _target_acceptance;
            unsigned                                _naccepts;
            unsigned                                _nattempts;
            bool                                    _tuning;
            std::vector&lt;double&gt;                     _prior_parameters;

            double                                  _heating_power;
            <span style="color:#0000ff"><strong>mutable PolytomyTopoPriorCalculator     _topo_prior_calculator;</strong></span>
            
            static const double                     _log_zero;
    };
 
    // member function bodies go here
    

Modify the calcLogTopologyPrior function

This function must be modified to take account of polytomous trees. It uses the TopoPriorCalculator class to compute the prior for a tree topology with a specified number of internal nodes.

    inline double Updater::calcLogTopologyPrior() const {   
        <span style="color:#0000ff"><strong>Tree::SharedPtr tree = _tree_manipulator-&gt;getTree();</strong></span>
        <span style="color:#0000ff"><strong>assert(tree);</strong></span>
        <span style="color:#0000ff"><strong>if (tree-&gt;isRooted())</strong></span>
            <span style="color:#0000ff"><strong>_topo_prior_calculator.chooseRooted();</strong></span>
        <span style="color:#0000ff"><strong>else</strong></span>
            <span style="color:#0000ff"><strong>_topo_prior_calculator.chooseUnrooted();</strong></span>
        <span style="color:#0000ff"><strong>_topo_prior_calculator.setNTax(tree-&gt;numLeaves());</strong></span>
        <span style="color:#0000ff"><strong>unsigned m = tree-&gt;numInternals();</strong></span>
                
<span style="color:#0000ff"><strong></strong></span>
        <span style="color:#0000ff"><strong>double log_topology_prior = _topo_prior_calculator.getLogNormalizedTopologyPrior(m);</strong></span>

<span style="color:#0000ff"><strong></strong></span>
        <span style="color:#0000ff"><strong>return log_topology_prior;</strong></span>
    }   

Add the setTopologyPriorOptions function

This function calls the appropriate functions in the PolytomyTopoPriorCalculator object to set the C parameter for the topology prior as well as the type of polytomy prior (resolution class or standard polytomy).

    inline void Updater::setTopologyPriorOptions(bool resclass, double C) { 
        _topo_prior_calculator.setC(C);
        if (resclass)
            _topo_prior_calculator.chooseResolutionClassPrior();
        else
            _topo_prior_calculator.choosePolytomyPrior();
    }   

Modify the calcLogEdgeLengthPrior function

The calcLogEdgeLengthPrior function must be modified because, if polytomies may be present, we must now count the number of edges in the tree and can can no longer use a formula that only works for binary trees.

    inline double Updater::calcLogEdgeLengthPrior() const { 
        double log_prior = 0.0;
        Tree::SharedPtr tree = _tree_manipulator-&gt;getTree();
        assert(tree);

        double TL = _tree_manipulator-&gt;calcTreeLength();
        <span style="color:#0000ff"><strong>//double n = tree-&gt;numLeaves();</strong></span>
        <span style="color:#0000ff"><strong>double num_edges = _tree_manipulator-&gt;countEdges();</strong></span>

        double a = _prior_parameters[0];    // shape of Gamma prior on TL
        double b = _prior_parameters[1];    // scale of Gamma prior on TL
        double c = _prior_parameters[2];    // parameter of Dirichlet prior on edge length proportions

        // Calculate Gamma prior on tree length (TL)
        double log_gamma_prior_on_TL = (a - 1.0)*log(TL) - TL/b - a*log(b) - std::lgamma(a);

        // Calculate Dirichlet prior on edge length proportions
        //
        // Note that, for n edges, the Dirichlet prior density is
        //
        // p1^{c-1} p2^{c-1} ... pn^{c-1}
        // ------------------------------
        //    Gamma(c)^n / Gamma(n*c)
        //
        // where n = num_edges, pk = edge length k / TL and Gamma is the Gamma function.
        // If c == 1, then both numerator and denominator equal 1, so it is pointless
        // do loop over edge lengths.
        double log_edge_length_proportions_prior = std::lgamma(num_edges*c);
        if (c != 1.0) {
            for (auto nd : tree-&gt;_preorder) {
                double edge_length_proportion = nd-&gt;_edge_length/TL;
                log_edge_length_proportions_prior += (c - 1.0)*log(edge_length_proportion);
            }
            log_edge_length_proportions_prior -= std::lgamma(c)*num_edges;
        }

        log_prior = log_gamma_prior_on_TL + log_edge_length_proportions_prior;
        return log_prior;
    }   

< 19.3 | 19.4 | 19.5 >