1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
|
#include "trees_impl.hpp"
#include "utils.hpp"
#include "validation.hpp"
#include <sstream>
#include <stack>
#include <terraces/rooting.hpp>
namespace terraces {
namespace {
bool is_rchild(const tree& t, index_t node) { return t[t[node].parent()].rchild() == node; }
void copy_subtree(const tree& src, tree& dst, index_t sub_root) {
std::stack<index_t> stack;
stack.push(sub_root);
while (!stack.empty()) {
auto idx = stack.top();
auto node = src[idx];
stack.pop();
if (is_leaf(node)) {
dst[idx].lchild() = none;
dst[idx].rchild() = none;
dst[idx].taxon() = node.taxon();
} else {
dst[idx].lchild() = node.lchild();
dst[idx].rchild() = node.rchild();
dst[idx].taxon() = none;
dst[node.lchild()].parent() = idx;
dst[node.rchild()].parent() = idx;
stack.push(node.lchild());
stack.push(node.rchild());
}
}
}
void copy_reversed(const tree& src, tree& dst, index_t child, index_t cur, index_t parent) {
auto coming_from_right = is_rchild(src, child);
auto opposite_idx = src[cur].child(!coming_from_right);
dst[cur].child(coming_from_right) = parent;
dst[cur].child(!coming_from_right) = opposite_idx;
dst[parent].parent() = cur;
dst[opposite_idx].parent() = cur;
copy_subtree(src, dst, opposite_idx);
}
void copy_reversed_final(const tree& src, tree& dst, index_t child, index_t cur,
index_t root_rchild) {
auto coming_from_right = is_rchild(src, child);
auto opposite_idx = src[cur].child(!coming_from_right);
dst[cur].child(coming_from_right) = root_rchild;
dst[cur].child(!coming_from_right) = opposite_idx;
dst[root_rchild].parent() = cur;
dst[opposite_idx].parent() = cur;
copy_subtree(src, dst, opposite_idx);
copy_subtree(src, dst, root_rchild);
}
void insert_root(tree& t, index_t lchild, index_t rchild, bool reverse) {
t[0].parent() = none;
t[0].taxon() = none;
t[0].child(reverse) = lchild;
t[0].child(!reverse) = rchild;
t[lchild].parent() = 0;
t[rchild].parent() = 0;
}
} // anonymous namespace
std::vector<bool> root_split(const tree& t, index_t num_leaves) {
std::vector<bool> split(num_leaves);
auto root = t[0];
foreach_preorder(
t,
[&](index_t i) {
auto node = t[i];
if (is_leaf(node)) {
split[node.taxon()] = true;
}
},
root.rchild());
return split;
}
tree reroot_at_node(const tree& t, index_t node_idx) {
utils::ensure<std::invalid_argument>(node_idx != 0, "can't reroot at the root");
terraces::check_rooted_tree(t);
// if the node is already below the root, we don't reroot.
if (t[node_idx].parent() == 0) {
return t;
}
tree out(t.size());
// everything below node_idx is unchanged
copy_subtree(t, out, node_idx);
// the root is placed between node_idx and its parent
// (keeping the orientation of node_idx wrt its parent)
bool is_node_rchild = is_rchild(t, node_idx);
insert_root(out, node_idx, t[node_idx].parent(), is_node_rchild);
auto prev = node_idx;
auto cur = t[node_idx].parent();
auto next = t[cur].parent();
while (next != 0) {
copy_reversed(t, out, prev, cur, next);
prev = utils::exchange(cur, utils::exchange(next, t[next].parent()));
}
// everything on the opposite side of the original root is unchanged
auto root_rchild = is_rchild(t, cur);
auto opposite_idx = t[0].child(!root_rchild);
copy_reversed_final(t, out, prev, cur, opposite_idx);
check_rooted_tree(out);
assert(is_isomorphic_unrooted(t, out));
return out;
}
void reroot_at_taxon_inplace(tree& t, index_t comp_taxon) {
// identify node corresponding to comp_taxon
index_t root_leaf = none;
for (index_t i = 0; i < t.size(); ++i) {
if (t[i].taxon() == comp_taxon) {
assert(root_leaf == none);
root_leaf = i;
}
}
assert(root_leaf != none && "The tree doesn't contain the given taxon");
terraces::check_rooted_tree(t);
// first: swap at inner nodes s.t. root_leaf is the rightmost leaf
index_t cur = root_leaf;
index_t p = t[cur].parent();
while (cur != 0) {
// if cur lies in the left subtree, swap it to the right
if (t[p].lchild() == cur) {
std::swap(t[p].lchild(), t[p].rchild());
}
cur = utils::exchange(p, t[p].parent());
}
// second: move the root down right until we meet root_leaf
auto& li = t[0].lchild();
auto& ri = t[0].rchild();
while (ri != root_leaf) {
auto& ln = t[li];
auto& rn = t[ri];
auto& rli = rn.lchild();
auto& rri = rn.rchild();
auto& rrn = t[rri];
std::swap(rrn.parent(), ln.parent());
std::tie(li, ri, rli, rri) = std::make_tuple(ri, rri, li, rli);
}
}
} // namespace terraces
|
