The Chain class encapsulates a Markov chain. In this first draft, it contains a GammaRateVarUpdater object that modifies only one parameter, but in future versions the Chain class will include additional updaters that have the ability to update the other model parameters as well as the tree topology and edge lengths.
Create a file chain.hpp that contains the following code.
#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"
namespace strom {
class Chain {
friend class Likelihood;
public:
typedef std::vector<Updater::SharedPtr> updater_vect_t;
typedef std::shared_ptr<Chain> SharedPtr;
Chain();
~Chain();
void clear();
void startTuning();
void stopTuning();
void setTreeFromNewick(std::string & newick);
unsigned createUpdaters(Model::SharedPtr model, Lot::SharedPtr lot, Likelihood::SharedPtr likelihood);
TreeManip::SharedPtr getTreeManip();
Model::SharedPtr getModel();
double getLogLikelihood() const;
void setHeatingPower(double p);
double getHeatingPower() const;
void setChainIndex(unsigned idx);
double getChainIndex() const;
Updater::SharedPtr findUpdaterByName(std::string name);
std::vector<std::string> getUpdaterNames() const;
std::vector<double> getAcceptPercentages() const;
std::vector<unsigned> getNumUpdates() const;
std::vector<double> getLambdas() const;
void setLambdas(std::vector<double> & v);
double calcLogLikelihood() const;
double calcLogJointPrior() const;
void start();
void stop();
void nextStep(int iteration);
private:
Model::SharedPtr _model;
Lot::SharedPtr _lot;
TreeManip::SharedPtr _tree_manipulator;
updater_vect_t _updaters;
unsigned _chain_index;
double _heating_power;
double _log_likelihood;
};
// member function bodies go here
}
The constructor calls the clear function, which returns the Chain object to its just-constructed state. The destructor does nothing.
inline Chain::Chain() {
//std::cout << "Chain being created" << std::endl;
clear();
}
inline Chain::~Chain() {
//std::cout << "Chain being destroyed" << std::endl;
}
This function return the Chain object to its just-constructed state. It sets the chain index to 0 and the heating power to 1 (i.e. this is the “cold” chain by default) and calls startTuning to tell all updaters to begin automatic tuning.
inline void Chain::clear() {
_log_likelihood = 0.0;
_updaters.clear();
_chain_index = 0;
setHeatingPower(1.0);
startTuning();
}
Calling the Chain function startTuning calls, in turn, the setTuning function of all its component updaters, passing the value true to each. Similarly, the stopTuning function passes false to the setTuning function of all its component updaters.
inline void Chain::startTuning() {
for (auto u : _updaters)
u->setTuning(true);
}
inline void Chain::stopTuning() {
for (auto u : _updaters)
u->setTuning(false);
}
This function provides a copy of the tree in the form of a Newick tree description string to the Chain. If _tree_manipulator does not currently exist, then a TreeManip object is created to hold the tree built from the supplied string.
Our only current updater GammaRateVarUpdater does not modify tree topology or edge lengths, so it may be surprising that we call its setTreeFromNewick function here. The answer is that all updaters need to calculate the log likelihood (see Chain::calcLogLikelihood below), and a tree is required in order to calculate the likelihood.
inline void Chain::setTreeFromNewick(std::string & newick) {
assert(_updaters.size() > 0);
if (!_tree_manipulator)
_tree_manipulator.reset(new TreeManip);
_tree_manipulator->buildFromNewick(newick, false, false);
for (auto u : _updaters)
u->setTreeManip(_tree_manipulator);
}
This function creates an updater for each independent model parameter that has not been fixed by the user and adds the updater to the _updaters vector. Each updater is provided with a copy of the model, the pseudorandom number generator (lot), and the likelihood. The updater’s setWeight function is called to set the updater’s _weight data member. The weight is used to determine the probability that an updater’s update function will be called in any given iteration. The sum of weights of all updaters is used to normalize the weights (done by calling each updater’s calcProb function).
inline unsigned Chain::createUpdaters(Model::SharedPtr model, Lot::SharedPtr lot, Likelihood::SharedPtr likelihood) {
_model = model;
_lot = lot;
_updaters.clear();
double wstd = 1.0;
double sum_weights = 0.0;
// Add rate variance parameter updater 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);
}
for (auto u : _updaters) {
u->calcProb(sum_weights);
}
return (unsigned)_updaters.size();
}
This is an accessor for the _tree_manipulator data member.
inline TreeManip::SharedPtr Chain::getTreeManip() {
return _tree_manipulator;
}
Chain does not need to have a data member to store the model because the Likelihood object already provides housing for the model (and also provides a convenient getModel accessor function that we can use here).
inline Model::SharedPtr Chain::getModel() {
return _model;
}
This is an accessor for the _log_likelihood data member.
inline double Chain::getLogLikelihood() const {
return _log_likelihood;
}
The setHeatingPower function is responsible for setting the heating power for all updaters. Right now, there is only one updater (GammaRateVarUpdater) but more lines will be uncommented as other updaters are created. The getHeatingPower accessor function may be used to query a Chain object to find out its current heating power. When two chains swap, it is easier to simply exchange heating powers rather than swap states (which would involve copying lots of parameter values as well as the tree topology and edge lengths). As a result, the setHeatingPower function gets called often during an MCMC analysis involving multiple chains.
inline double Chain::getHeatingPower() const {
return _heating_power;
}
inline void Chain::setHeatingPower(double p) {
_heating_power = p;
for (auto u : _updaters)
u->setHeatingPower(p);
}
The setChainIndex function is responsible for setting the index of the chain (stored in the data member _chain_index). The chain index will become important when we begin to use heated chains; when two chains swap places, it is easier to simply swap their indices and chain powers rather than copy the entire state (tree topology, edge lengths, state frequencies, etc.).
inline double Chain::getChainIndex() const {
return _chain_index;
}
inline void Chain::setChainIndex(unsigned idx) {
_chain_index = idx;
}
Loops through all updaters in the _updaters vector, returning the first one whose getUpdaterName function returns name (or nullptr if an updater named name is not found).
inline Updater::SharedPtr Chain::findUpdaterByName(std::string name) {
Updater::SharedPtr retval = nullptr;
for (auto u : _updaters) {
if (u->getUpdaterName() == name) {
retval = u;
break;
}
}
assert(retval != nullptr);
return retval;
}
Part of the output generated will be a list of the tuning values for all updaters, so it is important to be able to get a list of the names of all updaters. This list may be obtained by calling the getUpdaterNames function.
inline std::vector<std::string> Chain::getUpdaterNames() const {
std::vector<std::string> v;
for (auto u : _updaters)
v.push_back(u->getUpdaterName());
return v;
}
It is important to provide feedback to the user regarding the performance of each updater. To that end, this function returns a vector of the acceptance percentages for all updaters.
inline std::vector<double> Chain::getAcceptPercentages() const {
std::vector<double> v;
for (auto u : _updaters)
v.push_back(u->getAcceptPct());
return v;
}
This function reports the number of attempts to update each parameter. Along with getAcceptPercentages, this provides feedback to the user regarding the performance of each updater.
inline std::vector<unsigned> Chain::getNumUpdates() const {
std::vector<unsigned> v;
for (auto u : _updaters)
v.push_back(u->getNumUpdates());
return v;
}
The getLambdas function returns a list of the values of the tuning parameters for all updaters. The values are returned in the same order as the names returned by getLambdaNames. The setLambdas function sets the tuning parameter values for all updaters. Right now, we have only one updater (and thus one tuning parameter value), but the setLambdas function accepts a vector argument so that all values may be specified at once. Currently, only the first value in the vector v is used.
inline std::vector<double> Chain::getLambdas() const {
std::vector<double> v;
for (auto u : _updaters)
v.push_back(u->getLambda());
return v;
}
inline void Chain::setLambdas(std::vector<double> & v) {
assert(v.size() == _updaters.size());
unsigned index = 0;
for (auto u : _updaters) {
u->setLambda(v[index++]);
}
}
This function calculates the current log likelihood. Any updater can be used for this, as all have the capability of computing the likelihood. Here we call upon the GammaRateVarUpdater to help us out because it is currently the only updater we have.
inline double Chain::calcLogLikelihood() const {
return _updaters[0]->calcLogLikelihood();
}
Each updater (right now there is just one) is called upon to calculate the log prior for the model parameter that it is responsible for updating. The log priors for each parameter are added together to yield the log of the joint prior, which is returned.
inline double Chain::calcLogJointPrior() const {
double lnP = 0.0;
for (auto u : _updaters) {
lnP += u->calcLogPrior();
}
return lnP;
}
The start and stop functions are called before the Chain object begins simulating a Markov chain and after it finishes, respectively. The stop function is simply a placeholder for now, but the start function is responsible for spitting out the starting state of the chain to the output.
inline void Chain::start() {
_tree_manipulator->selectAllPartials();
_tree_manipulator->selectAllTMatrices();
_log_likelihood = calcLogLikelihood();
}
inline void Chain::stop() {
}
Finally, we come to the nextStep function, whose job it is to take one step in the Markov chain. It calls the update function for each updater, then reports (every 100 iterations) the current state and saves parameter values (currently just the gamma rate variance) in the output every sampling_freq iterations.
inline void Chain::nextStep(int iteration) {
assert(_lot);
double u = _lot->uniform();
double cumprob = 0.0;
unsigned i = 0;
for (auto updater : _updaters) {
cumprob += updater->_prob;
if (u <= cumprob)
break;
i++;
}
assert(i < _updaters.size());
_log_likelihood = _updaters[i]->update(_log_likelihood);
}