Additions to Node, Tree and TreeManip

We need to make a few small modifications to the Node, Tree and TreeManip classes before continuing.

Node class addition

In the Node class, add the following getter and setter member functions to the public part of the class declaration:

Node *        getLeftChild() {return _left_child;}
Node *        getRightSib() {return _right_sib;}
int           getNumber() {return _number;}
std::string   getName() {return _name;}
double        getEdgeLength() {return _edge_length;}
Split         getSplit() {return _split;}
void          setEdgeLength(double v) {_edge_length = v;}

Tree class addition

In the Tree class, add the following two accessor member functions to the public part of the class declaration:

bool      isRooted() const {return _is_rooted;}
unsigned  numLeaves() const {return _nleaves;}

Additions to the TreeManip class declaration

In the TreeManip class declaration, add these member function prototypes to the public part:

void    setTree(typename Tree<T>::SharedPtr tree);
T *     randomInternalEdge(double uniform01);
void    nniNodeSwap(T * a, T * b);
double  calcTreeLength() const;

TreeManip::setTree function body

This setter function allows us to assign a tree to a TreeManip object.

template <class T>
inline void TreeManip<T>::setTree(typename Tree<T>::SharedPtr tree)
    {
    assert(tree);
    _tree = tree;
    }

TreeManip::randomInternalEdge function body

This utility function makes it easy to obtain the node that maintains a randomly-determined internal edge. This function will be used in MCMC tree proposals to make modifications to the current tree topology.

template <class T>
inline T * TreeManip<T>::randomInternalEdge(double uniform_deviate)
    {
    assert(uniform_deviate >= 0.0);
    assert(uniform_deviate < 1.0);

    // Unrooted case:                        Rooted case:
    //
    // 2     3     4     5                   1     2     3     4
    //  \   /     /     /                     \   /     /     /
    //   \ /     /     /                       \ /     /     /
    //    8     /     /                         7     /     /
    //     \   /     /                           \   /     /
    //      \ /     /                             \ /     /
    //       7     /                               6     /
    //        \   /                                 \   /
    //         \ /                                   \ /
    //          6   nleaves = 5                       5   nleaves = 4
    //          |   num_internal_edges = 2            |   num_internal_edges = 2
    //          |   choose node 7 or node 8           |   choose node 6 or node 7
    //          1                                    root
    //
    // _preorder = [6, 7, 8, 2, 3, 4, 5]     _preorder = [5, 6, 7, 1, 2, 3, 4]
    //
    // Note: _preorder is actually a vector of T *, but is shown here as a
    // vector of integers solely to illustrate the algorithm below
    
    unsigned num_internal_edges = _tree->_nleaves - 2 - (_tree->_is_rooted ? 0 : 1);

    // Add one to skip first node in _preorder vector, which is an internal node whose edge
    // is either a terminal edge (if tree is unrooted) or invalid (if tree is rooted)
    unsigned index_of_chosen = 1 + (unsigned)std::floor(uniform_deviate*num_internal_edges);

    unsigned internal_nodes_visited = 0;
    T * chosen_node = 0;
    BOOST_FOREACH(T * nd, _tree->_preorder)
        {
        if (nd->_left_child)
            {
            if (internal_nodes_visited == index_of_chosen)
                {
                chosen_node = nd;
                break;
                }
            else
                ++internal_nodes_visited;
            }
        }
    assert(chosen_node);
    return chosen_node;
    }

TreeManip::nniNodeSwap function body

This utility function makes it easy to perform a Nearest Neighbor Interchange (NNI) modification to a tree topology. The two nodes provided are assumed to be related as indicated in the figure.

template <class T>
inline void TreeManip<T>::nniNodeSwap(T * a, T * b)
    {
    //     a                  b
    //      \   /              \   /
    //       \ /                \ /
    //        x     b            x     a
    //         \   /              \   /
    //          \ /      ==>       \ /
    //           |                  |
    //           |                  |
    //           |                  |
    //
    T * x = a->_parent;
    assert(x);
    
    T * y = b->_parent;
    assert(y);
    
    // Detach a from tree
    if (a->_right_sib)
        {
        x->_left_child = a->_right_sib;
        a->_parent = 0;
        a->_right_sib = 0;
        }
    else
        {
        x->_left_child->_right_sib = 0;
        a->_parent = 0;
        }

    // Detach b from tree
    if (b == x->_right_sib)
        {
        x->_right_sib = 0;
        b->_parent = 0;
        }
    else
        {
        y->_left_child = x;
        b->_right_sib = 0;
        b->_parent = 0;
        }
        
    // Reattach a to y
    assert(!x->_right_sib);
    x->_right_sib = a;
    a->_parent = y;
    
    // Reattach b to x
    assert(!x->_left_child->_right_sib);
    x->_left_child->_right_sib = b;
    b->_parent = x;
    
    refreshPreorder(_tree->_preorder[0]->_parent);
    }

TreeManip::calcTreeLength function body

This utility function simply sums the lengths of all edges in the tree.

template <class T>
inline double TreeManip<T>::calcTreeLength() const
    {
    double TL = 0.0;
    BOOST_FOREACH(T * nd, _tree->_preorder)
        {
        TL += nd->_edge_length;
        }
    return TL;
    }