diff options
Diffstat (limited to 'graphics/pgf/base/tex/generic/graphdrawing/lua/pgf/gd/phylogenetics/BalancedNearestNeighbourInterchange.lua')
-rw-r--r-- | graphics/pgf/base/tex/generic/graphdrawing/lua/pgf/gd/phylogenetics/BalancedNearestNeighbourInterchange.lua | 372 |
1 files changed, 372 insertions, 0 deletions
diff --git a/graphics/pgf/base/tex/generic/graphdrawing/lua/pgf/gd/phylogenetics/BalancedNearestNeighbourInterchange.lua b/graphics/pgf/base/tex/generic/graphdrawing/lua/pgf/gd/phylogenetics/BalancedNearestNeighbourInterchange.lua new file mode 100644 index 0000000000..1980afe7ad --- /dev/null +++ b/graphics/pgf/base/tex/generic/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$ + + + +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 |