19.7 Adding PolytomyUpdater to the Chain class

(Mac version)

< 19.6 | 19.7 | 19.8 >

Now that we’ve created the PolytomyUpdater class, we need to insert it into the updater rotation. Start by including the polytomy_updater.hpp header near the top of the chain.hpp.

#pragma once    

#include &lt;memory&gt;
#include &lt;boost/format.hpp&gt;
#include "lot.hpp"
#include "data.hpp"
#include "tree.hpp"
#include "model.hpp"
#include "likelihood.hpp"
#include "tree_manip.hpp"
#include "updater.hpp"
#include "gamma_ratevar_updater.hpp"
#include "omega_updater.hpp"
#include "pinvar_updater.hpp"
#include "statefreq_updater.hpp"
#include "exchangeability_updater.hpp"
#include "subset_relrate_updater.hpp"
#include "tree_updater.hpp"
<span style="color:#0000ff"><strong>#include "polytomy_updater.hpp"</strong></span>
#include "tree_length_updater.hpp"

Modifying the setTreeFromNewick function

In the setTreeFromNewick function, make the last argument to buildFromNewick equal to true so that polytomies are allowed when creating (for example) the starting tree.

    inline void Chain::setTreeFromNewick(std::string & newick) { 
        assert(_updaters.size() &gt; 0);
        if (!_tree_manipulator)
            _tree_manipulator.reset(new TreeManip);
        <span style="color:#0000ff"><strong>_tree_manipulator-&gt;buildFromNewick(newick, false, true);</strong></span>
        for (auto u : _updaters)
            u-&gt;setTreeManip(_tree_manipulator);
    } 

Modifying the createUpdaters function

Add the highlighted lines to the createUpdaters function to add PolytomyUpdater to the list of updaters potentially called during an MCMC iteration.

    inline unsigned Chain::createUpdaters(Model::SharedPtr model, Lot::SharedPtr lot, Likelihood::SharedPtr likelihood) { 
        _model = model;
        _lot = lot;
        _updaters.clear();

        double wstd             = 1.0;
        double wtreelength      = 1.0;
        double wtreetopology    = 19.0;
        <span style="color:#0000ff"><strong>double wpolytomy        = 0.0;</strong></span>
        double sum_weights      = 0.0;
        
        <span style="color:#0000ff"><strong>if (_model-&gt;isAllowPolytomies()) {</strong></span>
            <span style="color:#0000ff"><strong>wstd             = 1.0;</strong></span>
            <span style="color:#0000ff"><strong>wtreelength      = 2.0;</strong></span>
            <span style="color:#0000ff"><strong>wtreetopology    = 9.0;</strong></span>
            <span style="color:#0000ff"><strong>wpolytomy        = 9.0;</strong></span>
        <span style="color:#0000ff"><strong>}</strong></span>
        
        //...
        
        // Add tree updater and tree length updater to _updaters  
        if (!_model-&gt;isFixedTree()) {
            double tree_length_shape = 1.0;
            double tree_length_scale = 10.0;
            double dirichlet_param   = 1.0;
                        
            Updater::SharedPtr u = TreeUpdater::SharedPtr(new TreeUpdater());
            u-&gt;setLikelihood(likelihood);
            u-&gt;setLot(lot);
            u-&gt;setLambda(0.5);
            u-&gt;setTargetAcceptanceRate(0.3);
            u-&gt;setPriorParameters({tree_length_shape, tree_length_scale, dirichlet_param});
            <span style="color:#0000ff"><strong>u-&gt;setTopologyPriorOptions(_model-&gt;isResolutionClassTopologyPrior(), _model-&gt;getTopologyPriorC());</strong></span>
            u-&gt;setWeight(wtreetopology); sum_weights += wtreetopology;
            _updaters.push_back(u);

            <span style="color:#0000ff"><strong>if (_model-&gt;isAllowPolytomies()) {</strong></span>
                <span style="color:#0000ff"><strong>Updater::SharedPtr u = PolytomyUpdater::SharedPtr(new PolytomyUpdater());</strong></span>
                <span style="color:#0000ff"><strong>u-&gt;setLikelihood(likelihood);</strong></span>
                <span style="color:#0000ff"><strong>u-&gt;setLot(lot);</strong></span>
                <span style="color:#0000ff"><strong>u-&gt;setLambda(0.05);</strong></span>
                <span style="color:#0000ff"><strong>u-&gt;setTargetAcceptanceRate(0.5);</strong></span>
                <span style="color:#0000ff"><strong>u-&gt;setPriorParameters({tree_length_shape, tree_length_scale, dirichlet_param});</strong></span>
                <span style="color:#0000ff"><strong>u-&gt;setTopologyPriorOptions(_model-&gt;isResolutionClassTopologyPrior(), _model-&gt;getTopologyPriorC());</strong></span>
                <span style="color:#0000ff"><strong>u-&gt;setWeight(wpolytomy); sum_weights += wpolytomy;</strong></span>
                <span style="color:#0000ff"><strong>_updaters.push_back(u);</strong></span>
            <span style="color:#0000ff"><strong>}</strong></span>

            u = TreeLengthUpdater::SharedPtr(new TreeLengthUpdater());
            u-&gt;setLikelihood(likelihood);
            u-&gt;setLot(lot);
            u-&gt;setLambda(0.2);
            u-&gt;setTargetAcceptanceRate(0.3);
            u-&gt;setPriorParameters({tree_length_shape, tree_length_scale, dirichlet_param});
            <span style="color:#0000ff"><strong>u-&gt;setTopologyPriorOptions(_model-&gt;isResolutionClassTopologyPrior(), _model-&gt;getTopologyPriorC());</strong></span>
            u-&gt;setWeight(wtreelength); sum_weights += wtreelength;
            _updaters.push_back(u);
        } 
        
        for (auto u : _updaters) {
            u-&gt;calcProb(sum_weights);
        }
        
        return (unsigned)_updaters.size();
    } 

The calcLogJointPrior function

We must not let either the TreeLengthUpdator nor the PolytomyUpdater contribute to the calculation of the joint prior because TreeUpdater is already handling the tree topology and edge length components of the joint prior.

    inline double Chain::calcLogJointPrior() const { 
        double lnP = 0.0;
        for (auto u : _updaters) {
            <span style="color:#0000ff"><strong>if (u-&gt;_name != "Tree Length" && u-&gt;_name != "Polytomies" )</strong></span>
                lnP += u-&gt;calcLogPrior();
        }
        return lnP;
    } 