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< TreeUpdater > SharedPtr;
TreeUpdater();
~TreeUpdater();
virtual double calcLogPrior();
private:
virtual void revert();
virtual void proposeNewState();
void starTreeMove();
virtual void reset();
double _orig_edgelen_top;
double _orig_edgelen_middle;
double _orig_edgelen_bottom;
unsigned _case;
bool _topology_changed;
bool _star_tree_move;
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->randomChild(_lot, _x, 0, false);
_orig_edgelen_top = _a->getEdgeLength();
// Choose the second edge
_b = _tree_manipulator->randomChild(_lot, _x, _a, true);
if (!_b)
_b = _x;
_orig_edgelen_bottom = _b->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->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->selectPartialsHereToRoot(_x);
_a->setEdgeLength(new_edgelen_top);
_a->selectTMatrix();
_b->setEdgeLength(new_edgelen_bottom);
_b->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->getTree()->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
//
// For the star tree, there is only one internal node. In this case, only choose
// two edges and modify them (no change in tree topology is possible)
//
// a
// \ | /
// \|/
// x
// |
// |
// b
//
_x = _tree_manipulator->randomInternalEdge(_lot);
_orig_edgelen_middle = _x->getEdgeLength();
// The only child of the root node will be chosen only if the tree equals the star tree
// in which case we want to perform a starTreeMove rather than Larget-Simon
_star_tree_move = false;
if (_x->getParent() && !_x->getParent()->getParent()) {
_star_tree_move = true;
starTreeMove();
return;
}
_y = _x->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() {
if (_star_tree_move) {
_a->setEdgeLength(_orig_edgelen_top);
_b->setEdgeLength(_orig_edgelen_bottom);
}
else {
assert(_case > 0 && _case < 9);
if (_case == 2 || _case == 6)
_tree_manipulator->LargetSimonSwap(_a, _b);
else if (_case == 1 || _case == 5)
_tree_manipulator->LargetSimonSwap(_b, _a);
_a->setEdgeLength(_orig_edgelen_top);
_x->setEdgeLength(_orig_edgelen_middle);
if (_case == 1 || _case == 3 || _case == 5 || _case == 7)
_b->setEdgeLength(_orig_edgelen_bottom); // not actually necessary for case 3
else
_y->setEdgeLength(_orig_edgelen_bottom); // not actually necessary for case 4
}
}