As you know, the Likelihood class uses BeagleLib to compute the likelihood of a tree. BeagleLib computes the likelihood by working down the tree in reverse level order, combining the conditional likelihood arrays for each pair of child lineages to compute the conditional likelihood vector for the parent node.
If trees contain polytomies, then a wrench is thrown in the works because a polytomous node has more than two children. We will take advantage of the fact that
_nodes array of the Tree object to represent a fully resolved tree, andBegin by adding the highlighted lines below to the Likelihood class declaration in likelihood.hpp:
#pragma once
#include <map>
#include <boost/algorithm/string.hpp>
#include <boost/format.hpp>
#include <boost/shared_ptr.hpp>
#include <boost/range/adaptor/reversed.hpp>
#include "libhmsbeagle/beagle.h"
#include "tree.hpp"
#include "tree_manip.hpp"
#include "data.hpp"
#include "model.hpp"
#include "xstrom.hpp"
namespace strom {
class Likelihood {
public:
Likelihood();
~Likelihood();
void setRooted(bool is_rooted);
void setPreferGPU(bool prefer_gpu);
void setAmbiguityEqualsMissing(bool ambig_equals_missing);
bool usingStoredData() const;
void useStoredData(bool using_data);
void useUnderflowScaling(bool do_scaling);
std::string beagleLibVersion() const;
std::string availableResources() const;
std::string usedResources() const;
void initBeagleLib();
void finalizeBeagleLib(bool use_exceptions);
double calcLogLikelihood(Tree::SharedPtr t);
Data::SharedPtr getData();
void setData(Data::SharedPtr d);
Model::SharedPtr getModel();
void setModel(Model::SharedPtr m);
void clear();
unsigned calcNumEdgesInFullyResolvedTree() const;
unsigned calcNumInternalsInFullyResolvedTree() const;
private:
struct InstanceInfo {
int handle;
int resourcenumber;
std::string resourcename;
unsigned nstates;
unsigned nratecateg;
unsigned npatterns;
unsigned partial_offset;
unsigned tmatrix_offset;
bool invarmodel;
std::vector<unsigned> subsets;
InstanceInfo() : handle(-1), resourcenumber(-1), resourcename(""), nstates(0), nratecateg(0), npatterns(0), partial_offset(0), tmatrix_offset(0), invarmodel(false) {}
};
typedef std::pair<unsigned, int> instance_pair_t;
unsigned getScalerIndex(Node * nd, InstanceInfo & info) const;
unsigned getPartialIndex(Node * nd, InstanceInfo & info) const;
unsigned getTMatrixIndex(Node * nd, InstanceInfo & info, unsigned subset_index) const;
void updateInstanceMap(instance_pair_t & p, unsigned subset);
void newInstance(unsigned nstates, int nrates, std::vector<unsigned> & subset_indices);
void setTipStates();
void setTipPartials();
void setPatternPartitionAssignments();
void setPatternWeights();
void setAmongSiteRateHeterogenetity();
void setModelRateMatrix();
void addOperation(InstanceInfo & info, Node * nd, Node * lchild, Node * rchild, unsigned subset_index);
void queuePartialsRecalculation(Node * nd, Node * lchild, Node * rchild, Node * polytomy = 0);
void queueTMatrixRecalculation(Node * nd);
void defineOperations(Tree::SharedPtr t);
void updateTransitionMatrices();
void calculatePartials();
double calcInstanceLogLikelihood(InstanceInfo & inst, Tree::SharedPtr t);
std::vector<InstanceInfo> _instances;
std::map<int, std::string> _beagle_error;
std::map<int, std::vector<int> > _operations;
std::map<int, std::vector<int> > _pmatrix_index;
std::map<int, std::vector<double> > _edge_lengths;
std::map<int, std::vector<int> > _eigen_indices;
std::map<int, std::vector<int> > _category_rate_indices;
double _relrate_normalizing_constant;
std::vector<int> _subset_indices;
std::vector<int> _parent_indices;
std::vector<int> _child_indices;
std::vector<int> _tmatrix_indices;
std::vector<int> _weights_indices;
std::vector<int> _freqs_indices;
std::vector<int> _scaling_indices;
Model::SharedPtr _model;
Data::SharedPtr _data;
unsigned _ntaxa;
bool _rooted;
bool _prefer_gpu;
bool _ambiguity_equals_missing;
bool _underflow_scaling;
bool _using_data;
std::vector<Node *> _polytomy_helpers;
std::map<int, std::vector<int> > _polytomy_map;
std::vector<double> _identity_matrix;
public:
typedef std::shared_ptr< Likelihood > SharedPtr;
};
Add the indicated line to the clear function.
inline void Likelihood::clear() {
finalizeBeagleLib(true);
_ntaxa = 0;
_rooted = false;
_prefer_gpu = false;
_ambiguity_equals_missing = true;
_underflow_scaling = false;
_using_data = true;
_data = nullptr;
_operations.clear();
_pmatrix_index.clear();
_edge_lengths.clear();
_eigen_indices.clear();
_category_rate_indices.clear();
_relrate_normalizing_constant = 1.0;
_subset_indices.assign(1, 0);
_parent_indices.assign(1, 0);
_child_indices.assign(1, 0);
_tmatrix_indices.assign(1, 0);
_weights_indices.assign(1, 0);
_freqs_indices.assign(1, 0);
_scaling_indices.assign(1, 0);
_identity_matrix.assign(1, 0.0);
_model = Model::SharedPtr(new Model());
// Store BeagleLib error codes so that useful
// error messages may be provided to the user
_beagle_error.clear();
_beagle_error[0] = std::string("success");
_beagle_error[-1] = std::string("unspecified error");
_beagle_error[-2] = std::string("not enough memory could be allocated");
_beagle_error[-3] = std::string("unspecified exception");
_beagle_error[-4] = std::string("the instance index is out of range, or the instance has not been created");
_beagle_error[-5] = std::string("one of the indices specified exceeded the range of the array");
_beagle_error[-6] = std::string("no resource matches requirements");
_beagle_error[-7] = std::string("no implementation matches requirements");
_beagle_error[-8] = std::string("floating-point range exceeded");
}
Initialize the _identity_matrix data member in the first few lines of newInstance:
inline void Likelihood::newInstance(unsigned nstates, int nrates, std::vector<unsigned> & subset_indices) {
unsigned num_subsets = (unsigned)subset_indices.size();
bool is_invar_model = (nrates < 0 ? true : false);
unsigned ngammacat = (unsigned)(is_invar_model ? -nrates : nrates);
// Create an identity matrix used for computing partials
// for polytomies (represents the transition matrix
// for the zero-length edges inserted to arbitrarily
// resolve each polytomy)
_identity_matrix.assign(nstates*nstates*ngammacat, 0.0);
for (unsigned k = 0; k < ngammacat; k++) {
unsigned offset = k*nstates*nstates;
_identity_matrix[0+offset] = 1.0;
_identity_matrix[5+offset] = 1.0;
_identity_matrix[10+offset] = 1.0;
_identity_matrix[15+offset] = 1.0;
}
//...
}
Add the polytomy parameter and the highlighted lines below in the queuePartialsRecalculation member function:
inline void Likelihood::queuePartialsRecalculation(Node * nd, Node * lchild, Node * rchild, Node * polytomy) {
for (auto & info : _instances) {
unsigned instance_specific_subset_index = 0;
for (unsigned s : info.subsets) {
if (polytomy) {
// nd has been pulled out of tree's _unused_nodes vector to break up the polytomy
// Note that the parameter "polytomy" is the polytomous node itself
// First get the transition matrix index
unsigned tindex = getTMatrixIndex(nd, info, instance_specific_subset_index);
// Set the transition matrix for nd to the identity matrix
// note: last argument 1 is the value used for ambiguous states (should be 1 for transition matrices)
int code = beagleSetTransitionMatrix(info.handle, tindex, &_identity_matrix[0], 1);
if (code != 0)
throw XStrom(boost::str(boost::format("Failed to set transition matrix for instance %d. BeagleLib error code was %d (%s)") % info.handle % code % _beagle_error[code]));
// Set the edgelength to 0.0 to maintain consistency with the transition matrix
nd->setEdgeLength(0.0);
// If employing underflow scaling, the scaling factors for these fake nodes need to be
// transferred to the polytomous node, as that will be the only node remaining after the
// likelihood has been calculated. Save the scaling factor index, associating it with
// the scaling factor index of the polytomy node.
if (_underflow_scaling) {
// Get the polytomy's scaling factor index
int spolytomy = getScalerIndex(polytomy, info);
// Get nd's scaling factor index
int snd = getScalerIndex(nd, info);
// Save nd's index in the vector associated with polytomy's index
_polytomy_map[spolytomy].push_back(snd);
}
}
addOperation(info, nd, lchild, rchild, instance_specific_subset_index);
++instance_specific_subset_index;
}
}
}
This is the main modification made to the Likelihood class to allow for polytomies in trees. The vector _polytomy_stack starts out empty. Currently unused nodes are assigned to each polytomy until the polytomy is resolved (the number of children is just 2). Each fake internal node is assigned a transition matrix equal to the identity matrix because its edge length is necessarily zero (so no calculation of the transition matrix is needed). Operations are added to compute partials for each fake internal node and also the actual (polytomous) node. Each unused node that is put into service is stored in _polytomy_stack so that, after the likelihood is computed, all these fake internal nodes can be returned to their former unused state.
Add or change the highlighted lines in the defineOperations member function:
inline void Likelihood::defineOperations(Tree::SharedPtr t) {
assert(_instances.size() > 0);
assert(t);
assert(t->isRooted() == _rooted);
assert(_polytomy_helpers.empty());
assert(_polytomy_map.empty());
_relrate_normalizing_constant = _model->calcNormalizingConstantForSubsetRelRates();
// Start with a clean slate
for (auto & info : _instances) {
_operations[info.handle].clear();
_pmatrix_index[info.handle].clear();
_edge_lengths[info.handle].clear();
_eigen_indices[info.handle].clear();
_category_rate_indices[info.handle].clear();
}
// Loop through all nodes in reverse level order
for (auto nd : boost::adaptors::reverse(t->_levelorder)) {
assert(nd->_number >= 0);
if (!nd->_left_child) {
// This is a leaf
if (nd->isSelTMatrix())
queueTMatrixRecalculation(nd);
}
else {
// This is an internal node
if (nd->isSelTMatrix())
queueTMatrixRecalculation(nd);
// Internal nodes have partials to be calculated, so define
// an operation to compute the partials for this node
if (nd->isSelPartial()) {
TreeManip tm(t);
if (tm.isPolytomy(nd)) {
// Internal node is a polytomy
unsigned nchildren = tm.countChildren(nd);
Node * a = nd->_left_child;
Node * b = a->_right_sib;
Node * c = 0;
for (unsigned k = 0; k < nchildren - 2; k++) {
c = tm.getUnusedNode();
c->_left_child = a;
_polytomy_helpers.push_back(c);
queuePartialsRecalculation(c, a, b, nd);
// Tackle next arm of the polytomy
b = b->_right_sib;
a = c;
}
// Now add operation to compute the partial for the real internal node
queuePartialsRecalculation(nd, a, b);
}
else {
// Internal node is not a polytomy
Node * lchild = nd->_left_child;
assert(lchild);
Node * rchild = lchild->_right_sib;
assert(rchild);
queuePartialsRecalculation(nd, lchild, rchild);
}
}
}
}
}
The polytomy helper nodes that were brought out of storage to help in computing the partials for polytomous nodes will be put back into storage after the likelihood calculation, so we must accumulate the scalers for these nodes now and transfer the sum of log scalers to the polytomous node itself if we don’t want to lose them.
inline void Likelihood::calculatePartials() {
assert(_instances.size() > 0);
if (_operations.size() == 0)
return;
int code = 0;
// Loop through all instances
for (auto & info : _instances) {
unsigned nsubsets = (unsigned)info.subsets.size();
if (nsubsets > 1) {
code = beagleUpdatePartialsByPartition(
info.handle, // Instance number
(BeagleOperationByPartition *) &_operations[info.handle][0], // BeagleOperation list specifying operations
(int)(_operations[info.handle].size()/9)); // Number of operations
if (code != 0)
throw XStrom(boost::format("failed to update partials. BeagleLib error code was %d (%s)") % code % _beagle_error[code]);
if (_underflow_scaling) {
// Accumulate scaling factors across polytomy helpers and assign them to their parent node
for (auto & m : _polytomy_map) {
for (unsigned subset = 0; subset < nsubsets; subset++) {
code = beagleAccumulateScaleFactorsByPartition(info.handle, &m.second[0], (int)m.second.size(), m.first, subset);
if (code != 0) {
throw XStrom(boost::format("failed to transfer scaling factors to polytomous node. BeagleLib error code was %d (%s)") % code % _beagle_error[code]);
}
}
}
}
}
else {
// no partitioning, just one data subset
code = beagleUpdatePartials(
info.handle, // Instance number
(BeagleOperation *) &_operations[info.handle][0], // BeagleOperation list specifying operations
(int)(_operations[info.handle].size()/7), // Number of operations
BEAGLE_OP_NONE); // Index number of scaleBuffer to store accumulated factors
if (code != 0)
throw XStrom(boost::format("failed to update partials. BeagleLib error code was %d (%s)") % code % _beagle_error[code]);
if (_underflow_scaling) {
// Accumulate scaling factors across polytomy helpers and assign them to their parent node
for (auto & m : _polytomy_map) {
code = beagleAccumulateScaleFactors(info.handle, &m.second[0], (int)m.second.size(), m.first);
if (code != 0) {
throw XStrom(boost::format("failed to transfer scaling factors to polytomous node. BeagleLib error code was %d (%s)") % code % _beagle_error[code]);
}
}
}
}
}
}
All we need to do to the calcLogLikelihood function is to return the nodes we pulled out of storage back to their previous, unused state.
inline double Likelihood::calcLogLikelihood(Tree::SharedPtr t) {
assert(_instances.size() > 0);
if (!_using_data)
return 0.0;
// Must call setData and setModel before calcLogLikelihood
assert(_data);
assert(_model);
if (t->_is_rooted)
throw XStrom("This version of the program can only compute likelihoods for unrooted trees");
// Assuming "root" is leaf 0
assert(t->_root->_number == 0 && t->_root->_left_child == t->_preorder[0] && !t->_preorder[0]->_right_sib);
setModelRateMatrix();
setAmongSiteRateHeterogenetity();
defineOperations(t);
updateTransitionMatrices();
calculatePartials();
double log_likelihood = 0.0;
for (auto & info : _instances) {
log_likelihood += calcInstanceLogLikelihood(info, t);
}
// We no longer need the internal nodes brought out of storage
// and used to compute partials for polytomies
TreeManip tm(t);
for (Node * h : _polytomy_helpers) {
tm.putUnusedNode(h);
}
_polytomy_helpers.clear();
_polytomy_map.clear();
return log_likelihood;
}