14.3 The Chain Class

(Win version)

< 14.2 | 14.3 | 14.4 >

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 &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"

namespace strom {

    class Chain {
    
        friend class Likelihood;

        public:
        
            typedef std::vector&lt;Updater::SharedPtr&gt; updater_vect_t;
            typedef std::shared_ptr&lt;Chain&gt;          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&lt;std::string&gt;                getUpdaterNames() const;
            std::vector&lt;double&gt;                     getAcceptPercentages() const;
            std::vector&lt;unsigned&gt;                   getNumUpdates() const;
            std::vector&lt;double&gt;                     getLambdas() const;
            void                                    setLambdas(std::vector&lt;double&gt; & 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
    
}   

Constructor and destructor

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 &lt;&lt; "Chain being created" &lt;&lt; std::endl;
        clear();
    }

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

The clear member function

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();
    } 

The startTuning and stopTuning member functions

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-&gt;setTuning(true);
    }

    inline void Chain::stopTuning() {
        for (auto u : _updaters)
            u-&gt;setTuning(false);
    } 

The setTreeFromNewick member function

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() &gt; 0);
        if (!_tree_manipulator)
            _tree_manipulator.reset(new TreeManip);
        _tree_manipulator-&gt;buildFromNewick(newick, false, false);
        for (auto u : _updaters)
            u-&gt;setTreeManip(_tree_manipulator);
    } 

The createUpdaters member function

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-&gt;getRateVarParams();
        for (auto ratevar_shptr : ratevar_shptr_vect) {
            Updater::SharedPtr u = GammaRateVarUpdater::SharedPtr(new GammaRateVarUpdater(ratevar_shptr));
            u-&gt;setLikelihood(likelihood);
            u-&gt;setLot(lot);
            u-&gt;setLambda(1.0);
            u-&gt;setTargetAcceptanceRate(0.3);
            u-&gt;setPriorParameters({1.0, 1.0});
            u-&gt;setWeight(wstd); sum_weights += wstd;
            _updaters.push_back(u);
        }

        for (auto u : _updaters) {
            u-&gt;calcProb(sum_weights);
        }

        return (unsigned)_updaters.size();
    } 

The getTreeManip member function

This is an accessor for the _tree_manipulator data member.

    inline TreeManip::SharedPtr Chain::getTreeManip() { 
        return _tree_manipulator;
    } 

The getModel member function

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;
    } 

The getLogLikelihood member function

This is an accessor for the _log_likelihood data member.

    inline double Chain::getLogLikelihood() const {
        return _log_likelihood;
    } 

The getHeatingPower and setHeatingPower member functions

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-&gt;setHeatingPower(p);
    } 

The getChainIndex and setChainIndex member functions

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;
    } 

The findUpdaterByName member function

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-&gt;getUpdaterName() == name) {
                retval = u;
                break;
            }
        }
        assert(retval != nullptr);
        return retval;
    } 

The getUpdaterNames member function

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&lt;std::string&gt; Chain::getUpdaterNames() const { 
        std::vector&lt;std::string&gt; v;
        for (auto u : _updaters)
            v.push_back(u-&gt;getUpdaterName());
        return v;
    } 

The getAcceptPercentages member function

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&lt;double&gt; Chain::getAcceptPercentages() const { 
        std::vector&lt;double&gt; v;
        for (auto u : _updaters)
            v.push_back(u-&gt;getAcceptPct());
        return v;
    } 

The getNumUpdates member function

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&lt;unsigned&gt; Chain::getNumUpdates() const { 
        std::vector&lt;unsigned&gt; v;
        for (auto u : _updaters)
            v.push_back(u-&gt;getNumUpdates());
        return v;
    } 

The getLambdas and setLambdas member functions

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&lt;double&gt; Chain::getLambdas() const { 
        std::vector&lt;double&gt; v;
        for (auto u : _updaters)
            v.push_back(u-&gt;getLambda());
        return v;
    }

    inline void Chain::setLambdas(std::vector&lt;double&gt; & v) {
        assert(v.size() == _updaters.size());
        unsigned index = 0;
        for (auto u : _updaters) {
            u-&gt;setLambda(v[index++]);
        }
    } 

The calcLogLikelihood member function

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]-&gt;calcLogLikelihood();
    } 

The calcLogJointPrior member function

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-&gt;calcLogPrior();
        }
        return lnP;
    } 

The start and stop member functions

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-&gt;selectAllPartials();
        _tree_manipulator-&gt;selectAllTMatrices();
        _log_likelihood = calcLogLikelihood();
    }

    inline void Chain::stop() {
    } 

The nextStep member function

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-&gt;uniform();
        double cumprob = 0.0;
        unsigned i = 0;
        for (auto updater : _updaters) {
            cumprob += updater-&gt;_prob;
            if (u &lt;= cumprob)
                break;
            i++;
        }
        assert(i &lt; _updaters.size());
        _log_likelihood = _updaters[i]-&gt;update(_log_likelihood);
    } 

< 14.2 | 14.3 | 14.4 >