summaryrefslogtreecommitdiff
path: root/Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/BalancedNearestNeighbourInterchange.lua
blob: 1980afe7ad045c483e5cabb8ee24de0617105e94 (plain)
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$



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 Neighbor 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 Neighbor 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 neighboring 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