17.4 Modify the TreeManip Class

(Mac version)

< 17.3 | 17.4 | 17.5 >

Three helper functions need to be added to the TreeManip class:

Here are the lines that must be added to the class declaration in tree_manip.hpp.

#pragma once    

#include <cassert>
#include <memory>
#include <stack>
#include <queue>
#include <set>
#include <regex>
#include <boost/range/adaptor/reversed.hpp>
#include <boost/format.hpp>
#include "tree.hpp"
#include "lot.hpp"
#include "xstrom.hpp"

namespace strom {

    class TreeManip {
        public:
                                        TreeManip();
                                        TreeManip(Tree::SharedPtr t);
                                        ~TreeManip();

            void                        setTree(Tree::SharedPtr t);
            Tree::SharedPtr             getTree();

            double                      calcTreeLength() const;
            unsigned                    countEdges() const;
            void                        scaleAllEdgeLengths(double scaler);

            void                        createTestTree();
            std::string                 makeNewick(unsigned precision, bool use_names = false) const;

            void                        buildFromNewick(const std::string newick, bool rooted, bool allow_polytomies);
            void                        storeSplits(std::set<Split> & splitset);
            void                        rerootAtNodeNumber(int node_number);
        
            Node *                      randomInternalEdge(Lot::SharedPtr lot);
            Node *                      randomChild(Lot::SharedPtr lot, Node * x, Node * avoid, bool parent_included);
            void                        LargetSimonSwap(Node * a, Node * b);
        
            void                        selectAll();
            void                        deselectAll();
            void                        selectAllPartials();
            void                        deselectAllPartials();
            void                        selectAllTMatrices();
            void                        deselectAllTMatrices();
            void                        selectPartialsHereToRoot(Node * a);
            void                        flipPartialsAndTMatrices();

            void                        clear();

        private:

            Node *                      findNextPreorder(Node * nd);
            void                        refreshPreorder();
            void                        refreshLevelorder();
            void                        renumberInternals();
            void                        rerootAtNode(Node * prospective_root);
            void                        extractNodeNumberFromName(Node * nd, std::set<unsigned> & used);
            void                        extractEdgeLen(Node * nd, std::string edge_length_string);
            unsigned                    countNewickLeaves(const std::string newick);
            void                        stripOutNexusComments(std::string & newick);
            bool                        canHaveSibling(Node * nd, bool rooted, bool allow_polytomies);

            Tree::SharedPtr             _tree;

        public:

            typedef std::shared_ptr< TreeManip > SharedPtr;
    };
    

Here is the function body for the randomInternalEdge member function.

    inline Node * TreeManip::randomInternalEdge(Lot::SharedPtr lot) {   
        // Unrooted case:                        Rooted case:
        //
        // 2     3     4     5                   1     2     3     4
        //  \   /     /     /                     \   /     /     /
        //   \ /     /     /                       \ /     /     /
        //    8     /     /                         7     /     /
        //     \   /     /                           \   /     /
        //      \ /     /                             \ /     /
        //       7     /                               6     /
        //        \   /                                 \   /
        //         \ /                                   \ /
        //          6   nleaves = 5                       5   nleaves = 4
        //          |   preorder length = 7               |    preorder length = 7
        //          |   num_internal_edges = 7 - 5 = 2    |    num_internal_edges = 7 - 4 - 1 = 2
        //          1   choose node 7 or node 8          root  choose node 6 or node 7
        //
        // _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
        
        int num_internal_edges = (unsigned)_tree->_preorder.size() - _tree->_nleaves - (_tree->_is_rooted ? 1 : 0);

        // 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)
        double uniform_deviate = lot->uniform();
        unsigned index_of_chosen = 1 + (unsigned)std::floor(uniform_deviate*num_internal_edges);

        unsigned internal_nodes_visited = 0;
        Node * chosen_node = 0;
        for (auto 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;
    }   

Here is the function body for the randomChild member function.

    inline Node * TreeManip::randomChild(Lot::SharedPtr lot, Node * x, Node * avoid, bool parent_included) { 
        // Count number of children of x
        unsigned n = 0;
        Node * child = x->getLeftChild();
        while (child) {
            if (child != avoid)
                n++;
            child = child->getRightSib();
        }
        
        // Choose random child index
        unsigned upper = n + (parent_included ? 1 : 0);
        unsigned chosen = lot->randint(0,upper - 1);
        
        // If chosen < n, then find and return that particular child
        if (chosen < n) {
            child = x->getLeftChild();
            unsigned i = 0;
            while (child) {
                if (child != avoid && i == chosen)
                    return child;
                else if (child != avoid)
                    i++;
                child = child->getRightSib();
            }
        }

        // If chosen equals n, then the parent was chosen, indicated by returning NULL
        return NULL;
    } 

Here is the function body for the LargetSimonSwap member function.

    inline void TreeManip::LargetSimonSwap(Node * a, Node * b) {    
        // a and b are the ends of the selected 3-edge path in a Larget-Simon move
        // The 3-edge path is indicated by parentheses around the nodes involved.
        // x is always the parent of a
        // y can be the parent of b (case 1) or the child of b (case 2)
        
        Node * x = a->_parent;
        assert(x);
        
        Node * y = x->_parent;
        assert(y);
        
        if (y == b->_parent) {
            // Case 1: y is the parent of b
            //
            //    (a) d  e             (b) d  e
            //      \ | /                \ | /
            //       \|/                  \|/
            //       (x) f (b)            (x) f (a)    Swap a and b, leaving everything
            //         \ | /                \ | /      else as is
            //          \|/     ==>          \|/
            //          (y)                  (y)
            //           |                    |
            //           |                    |
            //           c                    c
            //

            // Detach a from tree
            if (a == x->_left_child) {
                x->_left_child = a->_right_sib;
            } else {
                Node * child = x->_left_child;
                while (child->_right_sib != a)
                    child = child->_right_sib;
                child->_right_sib = a->_right_sib;
            }
            a->_parent = 0;
            a->_right_sib = 0;
            
            // Detach b from tree
            if (b == y->_left_child) {
                y->_left_child = b->_right_sib;
            } else {
                Node * child = y->_left_child;
                while (child->_right_sib != b)
                    child = child->_right_sib;
                child->_right_sib = b->_right_sib;
            }
            b->_parent = 0;
            b->_right_sib = 0;

            // Reattach a to y
            a->_right_sib = y->_left_child;
            y->_left_child = a;
            a->_parent = y;
            
            // Reattach b to x
            b->_right_sib = x->_left_child;
            x->_left_child = b;
            b->_parent = x;
        }
        else {
            // Case 2: y is the child of b
            //
            //    (a) d  e             (a) f  c
            //      \ | /                \ | /
            //       \|/                  \|/
            //       (x) f  c            (x) d  e    swap x's children (except a)
            //         \ | /               \ | /     with y's children (except x)
            //          \|/     ==>         \|/
            //          (y)                 (y)
            //           |                   |
            //           |                   |
            //          (b)                 (b)
            assert(b == y->_parent);
            
            // Remove x's children from tree and store in xchildren stack
            std::stack<Node *> xchildren;
            Node * child = x->_left_child;
            Node * prevchild = 0;
            while (child) {
                if (child == a) {
                    prevchild = child;
                    child = child->_right_sib;
                } else {
                    if (child == x->_left_child) {
                        x->_left_child = child->_right_sib;
                        child->_right_sib = 0;
                        child->_parent = 0;
                        xchildren.push(child);
                        child = x->_left_child;
                    } else if (child->_right_sib) {
                        prevchild->_right_sib = child->_right_sib;
                        child->_right_sib = 0;
                        child->_parent = 0;
                        xchildren.push(child);
                        child = prevchild->_right_sib;
                    } else {
                        assert(prevchild == a);
                        a->_right_sib = 0;
                        child->_parent = 0;
                        xchildren.push(child);
                        child = 0;
                        prevchild = 0;
                    }
                }
            }
            
            // Remove y's children from tree and store in ychildren stack
            std::stack<Node *> ychildren;
            child = y->_left_child;
            prevchild = 0;
            while (child) {
                if (child == x) {
                    prevchild = child;
                    child = child->_right_sib;
                } else {
                    if (child == y->_left_child) {
                        y->_left_child = child->_right_sib;
                        child->_right_sib = 0;
                        child->_parent = 0;
                        ychildren.push(child);
                        child = y->_left_child;
                    } else if (child->_right_sib) {
                        prevchild->_right_sib = child->_right_sib;
                        child->_right_sib = 0;
                        child->_parent = 0;
                        ychildren.push(child);
                        child = prevchild->_right_sib;
                    } else {
                        assert(prevchild == x);
                        x->_right_sib = 0;
                        child->_parent = 0;
                        ychildren.push(child);
                        child = 0;
                        prevchild = 0;
                    }
                }
            }
            
            // Reattach xchildren to y
            while (!xchildren.empty()) {
                Node * popped = xchildren.top();
                xchildren.pop();
                popped->_right_sib = y->_left_child;
                y->_left_child = popped;
                popped->_parent = y;
            }

            // Reattach ychildren to x
            while (!ychildren.empty()) {
                Node * popped = ychildren.top();
                ychildren.pop();
                popped->_right_sib = x->_left_child;
                x->_left_child = popped;
                popped->_parent = x;
            }
        }
        
        refreshPreorder();
        refreshLevelorder();
    }   

< 17.3 | 17.4 | 17.5 >