University of Connecticut University of UC Title Fallback Connecticut

Testing the Chain class

Now all that is needed is to update the main function to set up Chain and run it. For this exercise we will use another variation on the rbcL data set used before.

Download rbcl10.zip, then unpack it and place the two files within (rbcl10.nex and rbcl10.tre) in your working directory with the other downloaded data and tree files. Then modify main.cpp as follows:

#include <iostream>
#include "node.hpp"
#include "tree_summary.hpp"
#include "likelihood.hpp"
#include "lot.hpp"
#include "chain.hpp"

using namespace strom;

int main(int argc, const char * argv[])
    {
    std::cout << "Starting..." << std::endl;

    try
        {
        // Read and store data
        Data::SharedPtr d(new Data());
        d->getDataFromFile("rbcl10.nex");

        // Create a substitution model
        GTRModel::SharedPtr gtr = GTRModel::SharedPtr(new GTRModel());
        gtr->setExchangeabilitiesAndStateFreqs(
            {0.06313,0.27216,0.06675,0.06353,0.48695,0.04748},
            {0.29311,0.18864,0.20845,0.30980});
        gtr->setGammaShape(0.257397);
        gtr->setGammaNCateg(4);
        std::cout << gtr->describeModel() << std::endl;

        // Create a likelihood object that will compute log-likelihoods
        Likelihood<Node>::SharedPtr likelihood(new Likelihood<Node>());
        likelihood->setData(d);
        likelihood->setModel(gtr);

        // Read in a tree
        TreeSummary<Node>::SharedPtr tree_summary(new TreeSummary<Node>());
        //tree_summary->readTreefile("rbcl738nj.tre", 0);
        tree_summary->readTreefile("rbcl10.tre", 0);
        Tree<Node>::SharedPtr tree = tree_summary->getTree(0);

        // Calculate the log-likelihood for the tree
        double lnL = likelihood->calcLogLikelihood(tree);
        std::cout << boost::str(boost::format("log likelihood = %.5f") % lnL) << std::endl;
        std::cout << "      (expecting -6723.248)" << std::endl;

        // Create a Chain object and take a step
        Lot::SharedPtr lot = Lot::SharedPtr(new Lot);
        lot->setSeed(12345);
        Chain<Node> chain;
        chain.setLot(lot);
        chain.setLikelihood(likelihood);
        chain.setTree(tree);
        chain.start();

        for (unsigned iteration = 0; iteration < 1000; ++iteration)
            chain.nextStep();

        }
    catch (XStrom x)
        {
        std::cerr << "Strom encountered a problem:\n  " << x.what() << std::endl;
        }

    std::cout << "\nFinished!" << std::endl;

    return 0;
    }

Upon running the program you should see output from 1000 calls to the nextStep function of the Chain object. A small portion of this output is shown below:

Starting...

----------------- GTR Model Info ------------------
State frequencies:   
  piA = 0.29311
  piC = 0.18864
  piG = 0.20845
  piT = 0.3098
Relative rates:    
  rAC = 0.06313
  rAG = 0.27216
  rAT = 0.06675
  rCG = 0.06353
  rCT = 0.48695
  rGT = 0.04748
Rate categories:   
  4
Gamma shape:       
  0.257397
Category boundaries and relative rate means:
    category        lower        upper    mean rate
           1      0.00000      0.01210      0.00247
           2      0.01210      0.18519      0.07208
           3      0.18519      1.05988      0.51683
           4      1.05988     infinity      3.40862
---------------------------------------------------

Constructing a TreeSummary
Sequence length:    1314
Number of taxa:     10
Number of patterns: 424
log likelihood = -6723.24766
      (expecting -6723.248)
     topochg     accepted     % accept       lambda         logu      logdiff      loglike     logprior           TL
                      ---          ---          ---          ---          ---  -6723.24766      7.32502      1.72969
                   reject          0.0      0.95000     -1.58689    -24.77552  -6723.24766      7.32502      1.72969
                   reject          0.0      0.90297     -0.77235    -98.83869  -6723.24766      7.32502      1.72969
                   reject          0.0      0.85871     -0.29073     -8.24978  -6723.24766      7.32502      1.72969
                   reject          0.0      0.81702     -0.91944    -13.53417  -6723.24766      7.32502      1.72969
                   accept         20.0      1.17054     -0.13676     -0.03321  -6723.46920      7.20182      1.74200
           *       reject         16.7      1.11480     -0.21538    -31.53573  -6723.46920      7.20182      1.74200
                   reject         14.3      1.06222     -0.82179     -4.61212  -6723.46920      7.20182      1.74200
                   reject         12.5      1.01258     -0.30489    -15.30611  -6723.46920      7.20182      1.74200
.
.
.
           *       reject         10.8      0.81917     -0.20811     -9.14952  -6732.68407      7.26181      1.73601
           *       reject         10.8      0.81543     -0.56954   -199.05239  -6732.68407      7.26181      1.73601
           *       reject         10.7      0.81171     -0.78569    -27.60402  -6732.68407      7.26181      1.73601
                   accept         10.8      0.84503     -1.56029     -0.73612  -6733.16510      7.35124      1.72706
                   reject         10.8      0.84118     -0.40550    -98.82861  -6733.16510      7.35124      1.72706
                   reject         10.8      0.83735     -3.30562    -13.58364  -6733.16510      7.35124      1.72706
           *       reject         10.8      0.83354     -1.64632    -37.10834  -6733.16510      7.35124      1.72706
Destroying a TreeSummary

Finished!