Three helper functions need to be added to the TreeManip class:
randomInternalEdge function returns a pointer to the Node object managing a randomly selected internal edgerandomChild function chooses a child of the specified node at random, and other parameters allow including the parent as a possible “child” and specifying a node to avoid when making the choice; andLargetSimonSwap function effects a swap of nodes a and b.
In the LargetSimonSwap function, a and b are the ends of the 3-edge segment modified during a Larget-Simon move, with a the higher (furthest from the root) of the two in the tree. If b is not the bottom-most node in the path, node a and b are both removed from the tree, then a is placed back into the tree where b once was and vice versa. If b is bottom-most, then the other nodes along the path are swapped (which accomplishes the same NNI swap but is less trouble). This function works even if the nodes in the path are polytomies.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"
<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<Split> & 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<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();
}