19.1 Modifying TreeUpdater

(Win version)

< 19.0 | 19.1 | 19.2 >

We will need to make one modification to the TreeUpdater class. Specificaly, we need to provide a special case proposal for the case in which the current tree is the star tree. The star tree (only one internal node) is a special case because it does not contain a 3-edge segment to modify. In the case of the star tree, we will simply adjust two randomly-chosen edge proportions. Add the highlighted starTreeMove function prototype as well as the _star_tree_move data member to the class declaration

    class TreeUpdater : public Updater {    

        friend class Chain;

        public:

            typedef std::shared_ptr&lt; TreeUpdater &gt; SharedPtr;

                                                TreeUpdater();
                                                ~TreeUpdater();

            virtual double                      calcLogPrior();

        private:

            virtual void                        revert();
            virtual void                        proposeNewState();
        
            <span style="color:#0000ff"><strong>void                                starTreeMove();</strong></span>

            virtual void                        reset();

            double                              _orig_edgelen_top;
            double                              _orig_edgelen_middle;
            double                              _orig_edgelen_bottom;

            unsigned                            _case;
            bool                                _topology_changed;
            <span style="color:#0000ff"><strong>bool                                _star_tree_move;</strong></span>
            Node *                              _x;
            Node *                              _y;
            Node *                              _a;
            Node *                              _b;
    }; 

Now add the following function body to tree_updater.hpp.

    inline void TreeUpdater::starTreeMove() {    
        // Choose focal 2-edge segment to modify
        _orig_edgelen_middle = 0.0;
        
        // Choose the first edge
        _a = _tree_manipulator-&gt;randomChild(_lot, _x, 0, false);
        _orig_edgelen_top = _a-&gt;getEdgeLength();

        // Choose the second edge
        _b = _tree_manipulator-&gt;randomChild(_lot, _x, _a, true);
        if (!_b)
            _b = _x;
        _orig_edgelen_bottom = _b-&gt;getEdgeLength();

        // Note that _a must be a child of _x, but _b may either be a different child of _x or _x itself
        double u = _lot-&gt;uniform();
        double new_edgelen_top    = u*(_orig_edgelen_top + _orig_edgelen_bottom);
        double new_edgelen_bottom = (1.0 - u)*(_orig_edgelen_top + _orig_edgelen_bottom);

        // Hastings ratio and Jacobian are both 1 under Gamma-Dirichlet parameterization
        _log_hastings_ratio = 0.0;
        _log_jacobian = 0.0;

        // Change edge lengths and flag partials and transition matrices for recalculation
        _tree_manipulator-&gt;selectPartialsHereToRoot(_x);
        _a-&gt;setEdgeLength(new_edgelen_top);
        _a-&gt;selectTMatrix();
        _b-&gt;setEdgeLength(new_edgelen_bottom);
        _b-&gt;selectTMatrix();
    }   

Call the starTreeMove function near the beginning of the proposeNewState function if the current tree is the star tree:

    inline void TreeUpdater::proposeNewState() {    
        _case = 0;
        _topology_changed = false;
        assert(!_tree_manipulator-&gt;getTree()-&gt;isRooted());

        // Choose random internal node x that is not the root and has parent y that is also not the root.
        // After choosing x (and thus y), choose a and b randomly to create a contiguous 3-edge segment.
        //
        //        a
        //   \ | /
        //    \|/
        //     x
        //      \ | /
        //       \|/
        //        y
        //        |
        //        |
        //        b
        //
        <span style="color:#0000ff"><strong>// For the star tree, there is only one internal node. In this case, only choose</strong></span>
        <span style="color:#0000ff"><strong>// two edges and modify them (no change in tree topology is possible)</strong></span>
        <span style="color:#0000ff"><strong>//</strong></span>
        <span style="color:#0000ff"><strong>//           a</strong></span>
        <span style="color:#0000ff"><strong>//      \ | /</strong></span>
        <span style="color:#0000ff"><strong>//       \|/</strong></span>
        <span style="color:#0000ff"><strong>//        x</strong></span>
        <span style="color:#0000ff"><strong>//        |</strong></span>
        <span style="color:#0000ff"><strong>//        |</strong></span>
        <span style="color:#0000ff"><strong>//        b</strong></span>
        <span style="color:#0000ff"><strong>//</strong></span>
        
        _x = _tree_manipulator-&gt;randomInternalEdge(_lot);
        _orig_edgelen_middle = _x-&gt;getEdgeLength();
        
        <span style="color:#0000ff"><strong>// The only child of the root node will be chosen only if the tree equals the star tree</strong></span>
        <span style="color:#0000ff"><strong>// in which case we want to perform a starTreeMove rather than Larget-Simon</strong></span>
        <span style="color:#0000ff"><strong>_star_tree_move = false;</strong></span>
        <span style="color:#0000ff"><strong>if (_x-&gt;getParent() && !_x-&gt;getParent()-&gt;getParent()) {</strong></span>
            <span style="color:#0000ff"><strong>_star_tree_move = true;</strong></span>
            <span style="color:#0000ff"><strong>starTreeMove();</strong></span>
            <span style="color:#0000ff"><strong>return;</strong></span>
        <span style="color:#0000ff"><strong>}</strong></span>

        _y = _x-&gt;getParent();
        //...
        

Finally, add a section to the TreeUpdater::revert function to revert a star tree move, if the current tree is the star tree.

    inline void TreeUpdater::revert() { 
        <span style="color:#0000ff"><strong>if (_star_tree_move) {</strong></span>
            <span style="color:#0000ff"><strong>_a-&gt;setEdgeLength(_orig_edgelen_top);</strong></span>
            <span style="color:#0000ff"><strong>_b-&gt;setEdgeLength(_orig_edgelen_bottom);</strong></span>
        <span style="color:#0000ff"><strong>}</strong></span>
        <span style="color:#0000ff"><strong>else {</strong></span>
            assert(_case &gt; 0 && _case &lt; 9);
            if (_case == 2 || _case == 6)
                _tree_manipulator-&gt;LargetSimonSwap(_a, _b);
            else if (_case == 1 || _case == 5)
                _tree_manipulator-&gt;LargetSimonSwap(_b, _a);
            _a-&gt;setEdgeLength(_orig_edgelen_top);
            _x-&gt;setEdgeLength(_orig_edgelen_middle);
            if (_case == 1 || _case == 3 || _case == 5 || _case == 7)
                _b-&gt;setEdgeLength(_orig_edgelen_bottom);    // not actually necessary for case 3
            else
                _y-&gt;setEdgeLength(_orig_edgelen_bottom);    // not actually necessary for case 4
        <span style="color:#0000ff"><strong>}</strong></span>
    }   

< 19.0 | 19.1 | 19.2 >