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
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
|
-- Copyright 2013 by Sarah Mäusle and Till Tantau
--
-- This file may be distributed an/or modified
--
-- 1. under the LaTeX Project Public License and/or
-- 2. under the GNU Public License
--
-- See the file doc/generic/pgf/licenses/LICENSE for more information
-- @release $Header: /cvsroot/pgf/pgf/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/BalancedNearestNeighbourInterchange.lua,v 1.2 2013/03/20 17:16:38 tantau Exp $
local BalancedNearestNeighbourInterchange = {}
-- Namespace
require("pgf.gd.phylogenetics").BalancedNearestNeighbourInterchange = BalancedNearestNeighbourInterchange
-- Imports
local InterfaceToAlgorithms = require("pgf.gd.interface.InterfaceToAlgorithms")
local DistanceMatrix = require("pgf.gd.phylogenetics.DistanceMatrix")
local lib = require("pgf.gd.lib")
-- Shorthand:
local declare = InterfaceToAlgorithms.declare
---
declare {
key = "balanced nearest neighbour interchange",
algorithm = BalancedNearestNeighbourInterchange,
phase = "phylogenetic tree optimization",
phase_default = true,
summary = [["
The BNNI (Balanced Nearest Neighbour Interchange) is a
postprocessing algorithm for phylogenetic trees. It swaps two
distant 3-subtrees if the total tree length is reduced by doing
so, until no such swaps are left.
"]],
documentation = [["
This algorithm is from Desper and Gascuel, \emph{Fast and
Accurate Phylogeny Reconstruction Algorithms Based on the
Minimum-Evolution Principle}, 2002.
"]]
}
---
declare {
key = "no phylogenetic tree optimization",
algorithm = { run = function(self) end },
phase = "phylogenetic tree optimization",
summary = [["
Switches off any phylogenetic tree optimization.
"]],
}
-- creates a binary heap, implementation as an array as described in
-- the respective wikipedia article
local function new_heap()
local heap = {}
function heap:insert(element, value)
local object = { element = element, value = value }
heap[#heap+1]= object
local i = #heap
local parent = math.floor(i/2)
-- sort the new object into its correct place
while heap[parent] and heap[parent].value < heap[i].value do
heap[i] = heap[parent]
heap[parent] = object
i = parent
parent = math.floor(i/2)
end
end
-- deletes the top element from the heap
function heap:remove_top_element()
-- replace first element with last and delete the last element
local element = heap[1].element
heap[1] = heap[#heap]
heap[#heap] = nil
local i = 1
local left_child = 2*i
local right_child = 2*i +1
-- sort the new top element into its correct place by swapping it
-- against its largest child
while heap[left_child] do
local largest_child = left_child
if heap[right_child] and heap[left_child].value < heap[right_child].value then
largest_child = right_child
end
if heap[largest_child].value > heap[i].value then
heap[largest_child], heap[i] = heap[i], heap[largest_child]
i = largest_child
left_child = 2*i
right_child = 2*i +1
else
return element
end
end
return element
end
return heap
end
-- BNNI (Balanced Nearest Neighbour Interchange)
-- [DESPER and GASCUEL: Fast and Accurate Phylogeny Reconstruction Algorithms Based on the Minimum-Evolution Principle, 2002]
-- swaps two distant-3 subtrees if the total tree length is reduced by doing so, until no such swaps are left
--
-- step 1: precomputation of all average distances between non-intersecting subtrees (already done by BME)
-- step 2: create heap of possible swaps
-- step 3: ( current tree with subtrees a,b,c,d: a--v-- {b, w -- {c, d}} )
-- (a): edge (v,w) is the best swap on the heap. Remove (v,c) and (w,b)
-- (b), (c), (d) : update the distance matrix
-- (e): remove the edge (v,w) from the heap; check the four edges adjacent to it for new possible swaps
-- (d): if the heap is non-empty, return to (a)
function BalancedNearestNeighbourInterchange:run()
local g = self.tree
-- create a heap of possible swaps
local possible_swaps = new_heap()
-- go over all arcs, look for possible swaps and add them to the heap [step 2]
for _, arc in ipairs (g.arcs) do
self:getBestSwap(arc, possible_swaps)
end
-- achieve best swap and update the distance matrix, until there is
-- no more swap to perform
while #possible_swaps > 0 do
-- get the best swap and delete it from the heap
local swap = possible_swaps:remove_top_element() --[part of step 3 (a)]
-- Check if the indicated swap is still possible. Another swap may
-- have interfered.
if g:arc(swap.v, swap.subtree1) and g:arc(swap.w, swap.subtree2) and g:arc(swap.v, swap.w) and g:arc(swap.a, swap.v) and g:arc(swap.d, swap.w) then
-- insert new arcs and delete the old ones to perform the swap [part of step 3 (a)]
-- disconnect old arcs
g:disconnect(swap.v, swap.subtree1)
g:disconnect(swap.subtree1, swap.v)
g:disconnect(swap.w, swap.subtree2)
g:disconnect(swap.subtree2, swap.w)
-- connect new arcs
g:connect(swap.v, swap.subtree2)
g:connect(swap.subtree2, swap.v)
g:connect(swap.w, swap.subtree1)
g:connect(swap.subtree1, swap.w)
--update distance matrix
self:updateBNNI(swap)
-- update heap: check neighbouring arcs for new possible swaps
-- [step 3 (e)]
self:getBestSwap(g:arc(swap.a,swap.v), possible_swaps)
self:getBestSwap(g:arc(swap.subtree2, swap.v), possible_swaps)
self:getBestSwap(g:arc(swap.d,swap.w), possible_swaps)
self:getBestSwap(g:arc(swap.subtree1, swap.w), possible_swaps)
end
end
end
--
-- Gets the distance between two nodes as specified in the distances
-- fields. Note: this function assumes that the distance from a to b
-- is the
-- same as the distance from b to a.
--
-- @param a,b The nodes
-- @return The distance between the two nodes
function BalancedNearestNeighbourInterchange:distance(a, b)
if a == b then
return 0
else
local distances = self.distances
return distances[a][b] or distances[b][a]
end
end
-- updates the distance matrix after a swap has been performed [step3(b),(c),(d)]
--
-- @param swap A table containing the information on the performed swap
-- subtree1, subtree2: the two subtrees, which
-- were swapped
-- a, d: The other two subtrees bordering the
-- swapping edge
-- v, w : the two nodes connecting the swapping edge
function BalancedNearestNeighbourInterchange:updateBNNI(swap)
local g = self.tree
local b = swap.subtree1
local c = swap.subtree2
local a = swap.a
local d = swap.d
local v = swap.v
local w = swap.w
local distances = self.distances
-- updates the distances in one of the four subtrees adjacent to the
-- swapping edge
function update_BNNI_subtree(swap, values)
local g = self.tree
local b = swap.farther
local c = swap.nearer
local a = swap.subtree
local v = swap.v
local d = swap.same
local w = swap.w
if not values then
values = {
visited = {[v] = v},
possible_ys = {v},
x = a,
y = v
}
-- if we're looking at subtrees in one of the swapped subtrees,
-- then need the old root (w) for the calculations
if swap.swapped_branch then values.possible_ys = {w} end
end
local visited = values.visited
local x = values.x
local y = values.y
local ys = values.possible_ys
local l = 0 -- number of edges between y and v
local dist_x_b = self:distance(x,b)
local dist_x_c = self:distance(x,c)
visited[x] = x --mark current x as visited
-- loop over possible y's:
for _, y in ipairs (ys) do
-- update distance [step 3(b)]
local distance = self:distance(x,y) - 2^(-l-2)*dist_x_b + 2^(-l-2)*dist_x_c
if y == w then y = v end -- the old distance w,x was used for the new distance calculation, but it needs to be
-- saved under its appropriate new name according to its new root. this case only arises when looking at x's
-- in one of the swapped subtrees (b or c)
distances[x][y] = distance
distances[y][x] = distance
l = l+1 -- length + 1, as the next y will be further away from v
end
-- update the distance between x and w (root of subtree c and d)
-- [step 3(c)]
local distance = 1/2 * (self:distance(x,b) + self:distance(x,d))
distances[x][w] = distance
distances[w][x] = distance
-- go to next possible x's
table.insert(ys, x) -- when we're at the next possible x, y can also be the current x
for _,arc in ipairs (g:outgoing(x)) do
if not visited[arc.head] then
values.x = arc.head
--go deeper
update_BNNI_subtree(swap, values)
end
end
end
-- name the nodes/subtrees in a general way that allows the use of the function update_BNNI_subtree
local update_a = {subtree = a, farther = b, nearer = c, v = v, same = d, w = w}
local update_b = {subtree = b, farther = a, nearer = d, v = w, same = c, w = v, swapped_branch = true}
local update_c = {subtree = c, farther = d, nearer = a, v = v, same = b, w = w, swapped_branch = true}
local update_d = {subtree = d, farther = c, nearer = b, v = w, same = a, w = v}
-- update the distances within the subtrees a,b,c,d respectively
update_BNNI_subtree(update_a)
update_BNNI_subtree(update_b)
update_BNNI_subtree(update_c)
update_BNNI_subtree(update_d)
-- update the distance between subtrees v and w [step 3 (d)]:
local distance = 1/4*( self:distance(a,b) + self:distance(a,d) + self:distance(c,b) + self:distance(c,d) )
distances[v][w] = distance
distances[w][v] = distance
end
-- finds the best swap across an arc and inserts it into the heap of
-- possible swaps
--
-- @param arc The arc, which is to be checked for possible swaps
-- @param heap_of_swaps The heap, containing all swaps, which
-- improve the total tree length
--
-- the following data of the swap are saved:
-- v,w = the nodes connecting the arc, across which the swap is
-- performed
-- subtree1,2 = the roots of the subtrees that are to be swapped
-- a,d = the roots of the two remaining subtrees adjacent to the arc
function BalancedNearestNeighbourInterchange:getBestSwap(arc, heap_of_swaps)
local g = self.tree
local possible_swaps = heap_of_swaps
local v = arc.tail
local w = arc.head
local is_leaf = self.is_leaf
-- only look at inner edges:
if not is_leaf[v] and not is_leaf[w] then
-- get the roots of the adjacent subtrees
local a, b, c, d
for _,outgoing in ipairs (g:outgoing(v)) do
local head = outgoing.head
if head ~= w then
a = a or head
b = head
end
end
for _,outgoing in ipairs (g:outgoing(w)) do
local head = outgoing.head
if head ~= v then
c = c or head
d = head
end
end
-- get the distances between the four subtrees
local a_b = self:distance(a,b)
local a_c = self:distance(a,c)
local a_d = self:distance(a,d)
local b_c = self:distance(b,c)
local b_d = self:distance(b,d)
local c_d = self:distance(c,d)
-- difference in total tree length between old tree (T) and new tree (T')
-- when nodes b and c are swapped
local swap1 = 1/4*(a_b + c_d - a_c - b_d )
-- difference in total tree length between old tree and new tree when nodes b and d are swapped
local swap2 = 1/4*(a_b + c_d - a_d - b_c)
-- choose the best swap that reduces the total tree length most (T-T' > 0)
if swap1 > swap2 and swap1 > 0 then
-- v,w = the nodes connecting the edge across which the swap is performed
-- subtree1 = one of the nodes to be swapped; connected to v
-- subtree2 = the other node to be swapped; connected to w
-- a = other node connected to v
-- d = other node connected to w
local swap = { v = v, w = w, subtree1 = b, subtree2 = c, a = a, d = d }
-- insert the swap into the heap
possible_swaps:insert(swap, swap1)
elseif swap2 > 0 then
local swap = { v = v, w = w, subtree1 = b, subtree2 = d, d = c, a = a }
possible_swaps:insert(swap, swap2)
end
end
end
return BalancedNearestNeighbourInterchange
|