17.4 Modify the TreeManip Class

(Linux 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 &lt;cassert&gt;
#include &lt;memory&gt;
#include &lt;stack&gt;
#include &lt;queue&gt;
#include &lt;set&gt;
#include &lt;regex&gt;
#include &lt;boost/range/adaptor/reversed.hpp&gt;
#include &lt;boost/format.hpp&gt;
#include "tree.hpp"
<span style="color:#0000ff"><strong>#include "lot.hpp"</strong></span>
#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&lt;Split&gt; & splitset);
            void                        rerootAtNodeNumber(int node_number);
        
            <span style="color:#0000ff"><strong>Node *                      randomInternalEdge(Lot::SharedPtr lot);</strong></span>
            <span style="color:#0000ff"><strong>Node *                      randomChild(Lot::SharedPtr lot, Node * x, Node * avoid, bool parent_included);</strong></span>
            <span style="color:#0000ff"><strong>void                        LargetSimonSwap(Node * a, Node * b);</strong></span>
        
            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&lt;unsigned&gt; & 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&lt; TreeManip &gt; 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-&gt;_preorder.size() - _tree-&gt;_nleaves - (_tree-&gt;_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-&gt;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-&gt;_preorder) {
            if (nd-&gt;_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-&gt;getLeftChild();
        while (child) {
            if (child != avoid)
                n++;
            child = child-&gt;getRightSib();
        }
        
        // Choose random child index
        unsigned upper = n + (parent_included ? 1 : 0);
        unsigned chosen = lot-&gt;randint(0,upper - 1);
        
        // If chosen &lt; n, then find and return that particular child
        if (chosen &lt; n) {
            child = x-&gt;getLeftChild();
            unsigned i = 0;
            while (child) {
                if (child != avoid && i == chosen)
                    return child;
                else if (child != avoid)
                    i++;
                child = child-&gt;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-&gt;_parent;
        assert(x);
        
        Node * y = x-&gt;_parent;
        assert(y);
        
        if (y == b-&gt;_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
            //          \|/     ==&gt;          \|/
            //          (y)                  (y)
            //           |                    |
            //           |                    |
            //           c                    c
            //

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

            // Reattach a to y
            a-&gt;_right_sib = y-&gt;_left_child;
            y-&gt;_left_child = a;
            a-&gt;_parent = y;
            
            // Reattach b to x
            b-&gt;_right_sib = x-&gt;_left_child;
            x-&gt;_left_child = b;
            b-&gt;_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)
            //          \|/     ==&gt;         \|/
            //          (y)                 (y)
            //           |                   |
            //           |                   |
            //          (b)                 (b)
            assert(b == y-&gt;_parent);
            
            // Remove x's children from tree and store in xchildren stack
            std::stack&lt;Node *&gt; xchildren;
            Node * child = x-&gt;_left_child;
            Node * prevchild = 0;
            while (child) {
                if (child == a) {
                    prevchild = child;
                    child = child-&gt;_right_sib;
                } else {
                    if (child == x-&gt;_left_child) {
                        x-&gt;_left_child = child-&gt;_right_sib;
                        child-&gt;_right_sib = 0;
                        child-&gt;_parent = 0;
                        xchildren.push(child);
                        child = x-&gt;_left_child;
                    } else if (child-&gt;_right_sib) {
                        prevchild-&gt;_right_sib = child-&gt;_right_sib;
                        child-&gt;_right_sib = 0;
                        child-&gt;_parent = 0;
                        xchildren.push(child);
                        child = prevchild-&gt;_right_sib;
                    } else {
                        assert(prevchild == a);
                        a-&gt;_right_sib = 0;
                        child-&gt;_parent = 0;
                        xchildren.push(child);
                        child = 0;
                        prevchild = 0;
                    }
                }
            }
            
            // Remove y's children from tree and store in ychildren stack
            std::stack&lt;Node *&gt; ychildren;
            child = y-&gt;_left_child;
            prevchild = 0;
            while (child) {
                if (child == x) {
                    prevchild = child;
                    child = child-&gt;_right_sib;
                } else {
                    if (child == y-&gt;_left_child) {
                        y-&gt;_left_child = child-&gt;_right_sib;
                        child-&gt;_right_sib = 0;
                        child-&gt;_parent = 0;
                        ychildren.push(child);
                        child = y-&gt;_left_child;
                    } else if (child-&gt;_right_sib) {
                        prevchild-&gt;_right_sib = child-&gt;_right_sib;
                        child-&gt;_right_sib = 0;
                        child-&gt;_parent = 0;
                        ychildren.push(child);
                        child = prevchild-&gt;_right_sib;
                    } else {
                        assert(prevchild == x);
                        x-&gt;_right_sib = 0;
                        child-&gt;_parent = 0;
                        ychildren.push(child);
                        child = 0;
                        prevchild = 0;
                    }
                }
            }
            
            // Reattach xchildren to y
            while (!xchildren.empty()) {
                Node * popped = xchildren.top();
                xchildren.pop();
                popped-&gt;_right_sib = y-&gt;_left_child;
                y-&gt;_left_child = popped;
                popped-&gt;_parent = y;
            }

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

< 17.3 | 17.4 | 17.5 >