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