summaryrefslogtreecommitdiff
path: root/Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/BalancedNearestNeighbourInterchange.lua
diff options
context:
space:
mode:
Diffstat (limited to 'Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/BalancedNearestNeighbourInterchange.lua')
-rw-r--r--Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/BalancedNearestNeighbourInterchange.lua372
1 files changed, 372 insertions, 0 deletions
diff --git a/Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/BalancedNearestNeighbourInterchange.lua b/Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/BalancedNearestNeighbourInterchange.lua
new file mode 100644
index 00000000000..a4fbc0cc3b9
--- /dev/null
+++ b/Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/BalancedNearestNeighbourInterchange.lua
@@ -0,0 +1,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