University of Connecticut University of UC Title Fallback Connecticut

The Chain class

The Chain class encapsulates a Markov chain. In this first draft, it mostly just maintains a TreeUpdater object that does essentially all of the work, but in future versions the Chain class will also maintain additional updaters that have the ability to update the GTR model parameters.

Create a file chain.hpp that contains the following code:

#pragma once

#include "data.hpp"
#include "likelihood.hpp"
#include "tree_manip.hpp"
#include "tree_updater.hpp"
#include "lot.hpp"

namespace strom
    {

    template<class T> class Likelihood;

    template <class T>
    class Chain
        {
        friend class Likelihood<T>;

        public:
                                                    Chain();
                                                    ~Chain();

            void                                    clear();

            void                                    start();
            void                                    nextStep();

            void                                    setTree(typename Tree<T>::SharedPtr tree);
            void                                    setLikelihood(typename Likelihood<T>::SharedPtr likelihood);
            void                                    setLot(typename Lot::SharedPtr lot);

            typedef boost::shared_ptr< Chain >      SharedPtr;

        private:

            typename TreeUpdater<T>::SharedPtr     _tree_updater;

        };

    template <class T>
    inline Chain<T>::Chain()
        {
        clear();
        }

    template <class T>
    inline Chain<T>::~Chain()
        {
        }

    template <class T>
    inline void Chain<T>::setTree(typename Tree<T>::SharedPtr tree)
        {
        if (!_tree_updater->_tree_manipulator)
            _tree_updater->_tree_manipulator.reset(new TreeManip<T>);
        _tree_updater->_tree_manipulator->setTree(tree);
        }

    template <class T>
    inline void Chain<T>::setLikelihood(typename Likelihood<T>::SharedPtr likelihood)
        {
        _tree_updater->setLikelihood(likelihood);
        }

    template <class T>
    inline void Chain<T>::setLot(typename Lot::SharedPtr lot)
        {
        _tree_updater->setLot(lot);
        }

    template <class T>
    inline void Chain<T>::clear()
        {
        _tree_updater.reset(new TreeUpdater<T>);
        _tree_updater->setLambda(1.0);
        _tree_updater->setTuning(true);
        _tree_updater->setTargetAcceptanceRate(0.1);
        }

    template <class T>
    inline void Chain<T>::start()
        {
        // temporary debugging output
        typename TreeUpdater<T>::kernel_t kernel = _tree_updater->calcLogPosteriorKernel();
        double TL = _tree_updater->_tree_manipulator->calcTreeLength();
        std::cout << boost::str(boost::format("%12s %12s %12s %12s %12s %12s %12s %12s %12s") % "topochg" % "accepted" % "% accept" % "lambda" % "logu" % "logdiff" % "loglike" % "logprior" % "TL") << std::endl;
        std::cout << boost::str(boost::format("%12s %12s %12s %12s %12s %12s %12.5f %12.5f %12.5f") % " " % "---" % "---" % "---" % "---" % "---" % kernel.first % kernel.second % TL) << std::endl;
        }

    template <class T>
    inline void Chain<T>::nextStep()
        {
        _tree_updater->update();
        }

    }

The start and nextStep functions

The only interesting member functions here are start() and nextStep(), and start() does nothing except print out column headers and nextStep simply calls the update function of TreeUpdater.