Testing the GTRModel class

To test the GTRModel class, you will need to tell your IDE about the Eigen library. The Eigen library comprises only header files, so there is nothing to build, and no libraries to which you need to link; you need only add its path to your project’s header search path.

Download and unpack Eigen

  • Download the latest stable release of Eigen 3 from http://eigen.tuxfamily.org/
  • Unpack the downloaded zip or tar.gz file to the directory of your choice (e.g. the same directory where you installed BeagleLib and Boost)
  • I will refer to the directory where Eigen is installed EIGEN_ROOT. Inside EIGEN_ROOT you should find a directory named Eigen (staring with a capital E).

Windows

In the Solution Explorer within the strom solution, right-click on the strom project and choose Properties. In the VC++ Directories section, add EIGEN_ROOT to Include Directories (remember that we are using EIGEN_ROOT to stand for the full path to the directory containing the Eigen code, so you should not type the word EIGEN_ROOT directly; instead, use the button to locate the Eigen directory).

Mac

  • Add an entry named EIGEN_ROOT to Xcode > Preferences… > Locations > Source Trees and assign it to the actual path to EIGEN_ROOT
  • Add $(EIGEN_ROOT) to the Header Search Paths in the Build Settings pane

Change the main function

Change main.cpp to the following and run your program in the debugger. The log-likelihood reported should now be -274.59466.

#include <iostream>
#include "node.hpp"
#include "tree_summary.hpp"
#include "likelihood.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("rbcL.nex");

        // Create a substitution model
        GTRModel::SharedPtr gtr = GTRModel::SharedPtr(new GTRModel());
        gtr->setExchangeabilitiesAndStateFreqs(
            {0.1307058, 0.1583282, 0.1598077, 0.2456609, 0.2202439, 0.08525337}, 
            {0.2249887, 0.2090694, 0.1375933, 0.4283486});
        gtr->setGammaShape(1.480126);
        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("rbcLjc.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 -274.59466)" << std::endl;
        }
    catch (XStrom x)
        {
        std::cerr << "Strom encountered a problem:\n  " << x.what() << std::endl;
        }

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

    return 0;
    }