Before we can test the new GTRModel class, we must integrate it into the Likelihood class.

Start by adding this include at the top of the likelihood.hpp file:

#include "gtr_model.hpp"

Add getModel and setModel member functions

Add these two declarations to the public section of the Likelihood class declaration:

void                setModel(GTRModel::SharedPtr model);
GTRModel::SharedPtr getModel();

Here are the bodies of these two functions. If a BeagleLib instance already exists, it must be destroyed (or finalized in the language of BeagleLib) and a new one created by the setModel function. This is because changing the model can change the dimension of the partials buffer.

template <class T>
inline void Likelihood<T>::setModel(GTRModel::SharedPtr model)
    {
    _model = model;
    if (_instance >= 0)
        {
        // init function was previously called, so set the model and create new BeagleLib instance
        int code = beagleFinalizeInstance(_instance);
        if (code != 0)
            throw XStrom(boost::str(boost::format("Likelihood setModel function failed to finalize BeagleLib instance. BeagleLib error code was %d (%s).") % code % _beagle_error[code]));
        _instance = -1;
        assert(_ntaxa > 0 && _nstates > 0 && _npatterns > 0);
        initBeagleLib();
        }
    }

template <class T>
inline GTRModel::SharedPtr Likelihood<T>::getModel()
    {
    return _model;
    }

Add a private data member _model

This data member should be added to the private section of the Likelihood class declaration.

GTRModel::SharedPtr  _model;

Initialize it in Likelihood constructor as follows:

_model = GTRModel::SharedPtr(new GTRModel());

Modify the initBeagleLib function

Specify the number of rate categories to be _model->_num_categ instead of 1.

_instance = beagleCreateInstance(
         _ntaxa,                    // tips
         num_internals,             // partials
         _ntaxa,                    // sequences
         _nstates,                  // states
         _npatterns,                // patterns
         1,                         // models
         num_transition_probs,      // transition matrices
         _model->_num_categ,        // rate categories
         0,                         // scale buffers
         NULL,                      // resource restrictions
         0,                         // length of resource list
         preferenceFlags,           // preferred flags
         requirementFlags,          // required flags
         &instance_details);        // pointer for details

Simplify the setDiscreteGammaShape function

The setDiscreteGammaShape member function can be simplified now because most of the work can be done by member functions of the GTRModel class. Replace the current version of setDiscreteGammaShape with this version.

template <class T>
inline void Likelihood<T>::setDiscreteGammaShape()
    {
    int code = _model->setBeagleAmongSiteRateVariationRates(_instance);
    if (code != 0)
        throw XStrom(boost::str(boost::format("failed to set category rates. BeagleLib error code was %d (%s)") % code % _beagle_error[code]));
    
    code = _model->setBeagleAmongSiteRateVariationProbs(_instance);
    if (code != 0)
        throw XStrom(boost::str(boost::format("failed to set category probabilities. BeagleLib error code was %d (%s)") % code % _beagle_error[code]));
    }

Simplify setModelRateMatrix function

We can also write setModelRateMatrix now using member functions of GTRModel.

template <class T>
inline void Likelihood<T>::setModelRateMatrix()
    {
    int code = _model->setBeagleStateFrequencies(_instance);
    if (code != 0)
        throw XStrom(boost::str(boost::format("failed to set state frequencies. BeagleLib error code was %d (%s)") % code % _beagle_error[code]));

    code = _model->setBeagleEigenDecomposition(_instance);
    if (code != 0)
        throw XStrom(boost::str(boost::format("failed to set eigen decomposition. BeagleLib error code was %d (%s)") % code % _beagle_error[code]));

    code = _model->setBeagleAmongSiteRateVariationRates(_instance);
    if (code != 0)
        throw XStrom(boost::str(boost::format("failed to set among-site rate variation rates. BeagleLib error code was %d (%s)") % code % _beagle_error[code]));

    code = _model->setBeagleAmongSiteRateVariationProbs(_instance);
    if (code != 0)
        throw XStrom(boost::str(boost::format("failed to set among-site rate variation weights. BeagleLib error code was %d (%s)") % code % _beagle_error[code]));
    }