Our new updater TreeUpdater and TreeLengthUpdater classes must be added to the Chain class in order to update the tree topology and edge lengths during the course of an MCMC analysis.
The changes that need to be made to chain.hpp are in bold, blue text. They are very similar to the changes you made to the Chain class to add the other updaters.
First, include the header files for TreeUpdater and TreeLengthUpdater at the top of chain.hpp.
#pragma once
#include <memory>
#include <boost/format.hpp>
#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"
<span style="color:#0000ff"><strong>#include "tree_updater.hpp"</strong></span>
<span style="color:#0000ff"><strong>#include "tree_length_updater.hpp"</strong></span>
Now modify the Chain::createUpdaters function to add TreeUpdater and TreeLengthUpdater to the mix.
inline unsigned Chain::createUpdaters(Model::SharedPtr model, Lot::SharedPtr lot, Likelihood::SharedPtr likelihood) {
_model = model;
_lot = lot;
_updaters.clear();
double wstd = 1.0;
<span style="color:#0000ff"><strong>double wtreelength = 1.0;</strong></span>
<span style="color:#0000ff"><strong>double wtreetopology = 19.0;</strong></span>
double sum_weights = 0.0;
// Add state frequency parameter updaters to _updaters
Model::state_freq_params_t & statefreq_shptr_vect = _model->getStateFreqParams();
for (auto statefreq_shptr : statefreq_shptr_vect) {
Updater::SharedPtr u = StateFreqUpdater::SharedPtr(new StateFreqUpdater(statefreq_shptr));
u->setLikelihood(likelihood);
u->setLot(lot);
u->setLambda(1.0);
u->setTargetAcceptanceRate(0.3);
u->setPriorParameters(std::vector<double>(statefreq_shptr->getStateFreqsSharedPtr()->size(), 1.0));
u->setWeight(wstd); sum_weights += wstd;
_updaters.push_back(u);
}
// Add exchangeability parameter updaters to _updaters
Model::exchangeability_params_t & exchangeability_shptr_vect = _model->getExchangeabilityParams();
for (auto exchangeability_shptr : exchangeability_shptr_vect) {
Updater::SharedPtr u = ExchangeabilityUpdater::SharedPtr(new ExchangeabilityUpdater(exchangeability_shptr));
u->setLikelihood(likelihood);
u->setLot(lot);
u->setLambda(1.0);
u->setTargetAcceptanceRate(0.3);
u->setPriorParameters({1.0, 1.0, 1.0, 1.0, 1.0, 1.0});
u->setWeight(wstd); sum_weights += wstd;
_updaters.push_back(u);
}
// Add rate variance parameter updaters to _updaters
Model::ratevar_params_t & ratevar_shptr_vect = _model->getRateVarParams();
for (auto ratevar_shptr : ratevar_shptr_vect) {
Updater::SharedPtr u = GammaRateVarUpdater::SharedPtr(new GammaRateVarUpdater(ratevar_shptr));
u->setLikelihood(likelihood);
u->setLot(lot);
u->setLambda(1.0);
u->setTargetAcceptanceRate(0.3);
u->setPriorParameters({1.0, 1.0});
u->setWeight(wstd); sum_weights += wstd;
_updaters.push_back(u);
}
// Add pinvar parameter updaters to _updaters
Model::pinvar_params_t & pinvar_shptr_vect = _model->getPinvarParams();
for (auto pinvar_shptr : pinvar_shptr_vect) {
Updater::SharedPtr u = PinvarUpdater::SharedPtr(new PinvarUpdater(pinvar_shptr));
u->setLikelihood(likelihood);
u->setLot(lot);
u->setLambda(0.5);
u->setTargetAcceptanceRate(0.3);
u->setPriorParameters({1.0, 1.0});
u->setWeight(wstd); sum_weights += wstd;
_updaters.push_back(u);
}
// Add omega parameter updaters to _updaters
Model::omega_params_t & omega_shptr_vect = _model->getOmegaParams();
for (auto omega_shptr : omega_shptr_vect) {
Updater::SharedPtr u = OmegaUpdater::SharedPtr(new OmegaUpdater(omega_shptr));
u->setLikelihood(likelihood);
u->setLot(lot);
u->setLambda(1.0);
u->setTargetAcceptanceRate(0.3);
u->setPriorParameters({1.0, 1.0});
u->setWeight(wstd); sum_weights += wstd;
_updaters.push_back(u);
}
// Add subset relative rate parameter updater to _updaters
if (!_model->isFixedSubsetRelRates()) {
Updater::SharedPtr u = SubsetRelRateUpdater::SharedPtr(new SubsetRelRateUpdater(_model));
u->setLikelihood(likelihood);
u->setLot(lot);
u->setLambda(1.0);
u->setTargetAcceptanceRate(0.3);
u->setPriorParameters(std::vector<double>(_model->getNumSubsets(), 1.0));
u->setWeight(wstd); sum_weights += wstd;
_updaters.push_back(u);
}
<span style="color:#0000ff"><strong>// Add tree updater and tree length updater to _updaters</strong></span>
<span style="color:#0000ff"><strong>if (!_model->isFixedTree()) {</strong></span>
<span style="color:#0000ff"><strong>double tree_length_shape = 1.0;</strong></span>
<span style="color:#0000ff"><strong>double tree_length_scale = 10.0;</strong></span>
<span style="color:#0000ff"><strong>double dirichlet_param = 1.0;</strong></span>
<span style="color:#0000ff"><strong>Updater::SharedPtr u = TreeUpdater::SharedPtr(new TreeUpdater());</strong></span>
<span style="color:#0000ff"><strong>u->setLikelihood(likelihood);</strong></span>
<span style="color:#0000ff"><strong>u->setLot(lot);</strong></span>
<span style="color:#0000ff"><strong>u->setLambda(0.5);</strong></span>
<span style="color:#0000ff"><strong>u->setTargetAcceptanceRate(0.3);</strong></span>
<span style="color:#0000ff"><strong>u->setPriorParameters({tree_length_shape, tree_length_scale, dirichlet_param});</strong></span>
<span style="color:#0000ff"><strong>u->setWeight(wtreetopology); sum_weights += wtreetopology;</strong></span>
<span style="color:#0000ff"><strong>_updaters.push_back(u);</strong></span>
<span style="color:#0000ff"><strong>u = TreeLengthUpdater::SharedPtr(new TreeLengthUpdater());</strong></span>
<span style="color:#0000ff"><strong>u->setLikelihood(likelihood);</strong></span>
<span style="color:#0000ff"><strong>u->setLot(lot);</strong></span>
<span style="color:#0000ff"><strong>u->setLambda(0.2);</strong></span>
<span style="color:#0000ff"><strong>u->setTargetAcceptanceRate(0.3);</strong></span>
<span style="color:#0000ff"><strong>u->setPriorParameters({tree_length_shape, tree_length_scale, dirichlet_param});</strong></span>
<span style="color:#0000ff"><strong>u->setWeight(wtreelength); sum_weights += wtreelength;</strong></span>
<span style="color:#0000ff"><strong>_updaters.push_back(u);</strong></span>
<span style="color:#0000ff"><strong>}</strong></span>
for (auto u : _updaters) {
u->calcProb(sum_weights);
}
return (unsigned)_updaters.size();
}
As with the other updaters, we are hard-coding the prior for tree topology and edge lengths. We will add options in Step 19 so that the user can change these values, but for now we are keeping things simple and baking in the priors. In this case, the tree length prior is an Exponential distribution with mean 10, and the prior on edge length proportions is a flat Dirichlet distribution. The prior on tree topology is Discrete Uniform over the space of all possible tree topologies (and tree topologies are assumed to be unrooted and binary).
I mentioned (when creating the TreeLengthUpdater class) that we would need to be careful to avoid counting the tree length prior and edge length proportions prior twice. We will now take care of that by adding a line to the calcLogJointPrior function that causes the TreeLengthUpdater to be skipped when computing the joint prior.
inline double Chain::calcLogJointPrior() const {
double lnP = 0.0;
for (auto u : _updaters) {
<span style="color:#0000ff"><strong>if (u->_name != "Tree Length")</strong></span>
lnP += u->calcLogPrior();
}
return lnP;
}