-- 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 BalancedMinimumEvolution = {} -- Namespace require("pgf.gd.phylogenetics").BalancedMinimumEvolution = BalancedMinimumEvolution -- Imports local InterfaceToAlgorithms = require("pgf.gd.interface.InterfaceToAlgorithms") local DistanceMatrix = require("pgf.gd.phylogenetics.DistanceMatrix") local Storage = require("pgf.gd.lib.Storage") local Digraph = require("pgf.gd.model.Digraph") local lib = require("pgf.gd.lib") -- Shorthand: local declare = InterfaceToAlgorithms.declare --- declare { key = "balanced minimum evolution", algorithm = BalancedMinimumEvolution, phase = "phylogenetic tree generation", summary = [[" The BME (Balanced Minimum Evolution) algorithm tries to minimize the total tree length. "]], documentation = [[" This algorithm is from Desper and Gascuel, \emph{Fast and Accurate Phylogeny Reconstruction Algorithms Based on the Minimum-Evolution Principle}, 2002. The tree is built in a way that minimizes the total tree length. The leaves are inserted into the tree one after another, creating new edges and new nodes. After every insertion the distance matrix has to be updated. "]], examples = [[" \tikz \graph [phylogenetic tree layout, balanced minimum evolution, grow'=right, sibling distance=0pt, distance matrix={ 0 4 9 9 9 9 9 4 0 9 9 9 9 9 9 9 0 2 7 7 7 9 9 2 0 7 7 7 9 9 7 7 0 3 5 9 9 7 7 3 0 5 9 9 7 7 5 5 0}] { a, b, c, d, e, f, g }; "]] } function BalancedMinimumEvolution:run() self.tree = Digraph.new(self.main_algorithm.digraph) self.distances = Storage.newTableStorage() local vertices = self.tree.vertices -- Sanity checks: if #vertices == 2 then self.tree:connect(vertices[1],vertices[2]) return self.tree elseif #vertices > 2 then -- Setup storages: self.is_leaf = Storage.new() -- First, build the initial distance matrix: local matrix = DistanceMatrix.graphDistanceMatrix(self.tree) -- Store distance information in the distance fields of the storages: for _,u in ipairs(vertices) do for _,v in ipairs(vertices) do self.distances[u][v] = matrix[u][v] end end -- Run BME self:runBME() -- Run postoptimizations local optimization_class = self.tree.options.algorithm_phases['phylogenetic tree optimization'] optimization_class.new { main_algorithm = self.main_algorithm, tree = self.tree, matrix = self.matrix, distances = self.distances, is_leaf = self.is_leaf, }:run() end -- Finish self:computeFinalLengths() self:createFinalEdges() return self.tree end -- the BME (Balanced Minimum Evolution) algorithm -- [DESPER and GASCUEL: Fast and Accurate Phylogeny Reconstruction -- Algorithms Based on the Minimum-Evolution Principle, 2002] -- -- The tree is built in a way that minimizes the total tree length. -- The leaves are inserted into the tree one after another, creating new edges and new nodes. -- After every insertion the distance matrix has to be updated. function BalancedMinimumEvolution:runBME() local g = self.tree local leaves = {} local is_leaf = self.is_leaf local distances = self.distances -- get user input for i, vertex in ipairs (g.vertices) do leaves[i] = vertex is_leaf[vertex] = true end -- create the new node which will be connected to the first three leaves local new_node = InterfaceToAlgorithms.createVertex( self.main_algorithm, { name = "BMEnode"..#g.vertices+1, generated_options = { { key = "phylogenetic inner node" } } } ) g:add {new_node} -- set the distances of new_node to subtrees local distance_1_2 = self:distance(leaves[1],leaves[2]) local distance_1_3 = self:distance(leaves[1],leaves[3]) local distance_2_3 = self:distance(leaves[2],leaves[3]) distances[new_node][leaves[1]] = 0.5*(distance_1_2 + distance_1_3) distances[new_node][leaves[2]] = 0.5*(distance_1_2 + distance_2_3) distances[new_node][leaves[3]] = 0.5*(distance_1_3 + distance_2_3) --connect the first three leaves to the new node for i = 1,3 do g:connect(new_node, leaves[i]) g:connect(leaves[i], new_node) end for k = 4,#leaves do -- compute distance from k to any subtree local k_dists = Storage.newTableStorage() for i = 1,k-1 do -- note that the function called stores the k_dists before they are overwritten self:computeAverageDistancesToAllSubtreesForK(g.vertices[i], { }, k,k_dists) end -- find the best insertion point local best_arc = self:findBestEdge(g.vertices[1],nil,k_dists) local head = best_arc.head local tail = best_arc.tail -- remove the old arc g:disconnect(tail, head) g:disconnect(head, tail) -- create the new node local new_node = InterfaceToAlgorithms.createVertex( self.main_algorithm, { name = "BMEnode"..#g.vertices+1, generated_options = { { key = "phylogenetic inner node" } } } ) g:add{new_node} -- gather the vertices that will be connected to the new node... local vertices_to_connect = { head, tail, leaves[k] } -- ...and connect them for _, vertex in pairs (vertices_to_connect) do g:connect(new_node, vertex) g:connect(vertex, new_node) end if not is_leaf[tail] then distances[leaves[k]][tail] = k_dists[head][tail] end if not is_leaf[head] then distances[leaves[k]][head] = k_dists[tail][head] end -- insert distances from k to subtrees into actual matrix... self:setAccurateDistancesForK(new_node,nil,k,k_dists,leaves) -- set the distance from k to the new node, which was created by inserting k into the graph distances[leaves[k]][new_node] = 0.5*( self:distance(leaves[k], head) + self:distance(leaves[k],tail)) -- update the average distances local values = {} values.s = head -- s--u is the arc into which k has been inserted values.u = tail values.new_node = new_node -- the new node created by inserting k self:updateAverageDistances(new_node, values,k,leaves) end end -- -- Updates the average distances from k to all subtrees -- -- @param vertex The starting point of the recursion -- @param values The values needed for the recursion -- - s, u The nodes which span the edge into which k has been -- inserted -- - new_node The new_node which has been created to insert k -- - l (l-1) is the number of edges between the -- new_node and the current subtree Y -- -- values.new_node, values.u and values.s must be set -- the depth first search must begin at the new node, thus vertex -- must be set to the newly created node function BalancedMinimumEvolution:updateAverageDistances(vertex, values, k, leaves) local g = self.tree local leaf_k = leaves[k] local y, z, x if not values.visited then values.visited = {} values.visited[leaf_k] = leaf_k -- we don't want to visit k! end -- there are (l-1) edges between new_node and y if not values.l then values.l = 1 end if not values.new_node then values.new_node = g:outgoing(leaf_k)[1].head end --values.s and values.u must be set -- the two nodes which connect the edge on which k was inserted: s,u local new_node = values.new_node local l = values.l local visited = values.visited visited[vertex] = vertex -- computes the distances to Y{k} for all subtrees X of Z function loop_over_x( x, y, values ) local l = values.l local y1= values.y1 -- calculate distance between Y{k} and X local old_distance -- the distance between Y{/k} and X needed for calculating the new distance if y == new_node then -- this y didn't exist in the former tree; so use y1 (see below) old_distance = self:distance(x,y1) else old_distance = self:distance(x,y) end local new_distance = old_distance + math.pow(2,-l) * ( self:distance(leaf_k,x) - self:distance(x,y1) ) self.distances[x][y] = new_distance self.distances[y][x] = new_distance -- symmetric matrix values.x_visited[x] = x --go deeper to next x for _, x_arc in ipairs (self.tree:outgoing(x)) do if not values.x_visited[x_arc.head] then local new_x = x_arc.head loop_over_x( new_x, y, values ) end end end --loop over Z's for _, arc in ipairs (self.tree:outgoing(vertex)) do if not visited[arc.head] then -- set y1, which is the node which was pushed further away from -- subtree Z by inserting k if arc.head == values.s then values.y1 = values.u elseif arc.head == values.u then values.y1 = values.s else assert(values.y1,"no y1 set!") end z = arc.head -- root of the subtree we're looking at y = arc.tail -- the root of the subtree-complement of Z x = z -- the first subtree of Z is Z itself values.x_visited = {} values.x_visited[y] = y -- we don't want to go there, as we want to stay within Z loop_over_x( z,y, values ) -- visit all possible subtrees of Z -- go to next Z values.l = values.l+1 -- moving further away from the new_node self:updateAverageDistances(z,values,k,leaves) values.l = values.l-1 -- moving back to the new_node end end end -- -- Computes the average distances of a node, which does not yet belong -- to the graph, to all subtrees. This is done using a depth first -- search -- -- @param vertex The starting point of the depth first search -- @param values The values for the recursion -- - distances The table in which the distances are to be -- stored -- - outgoing_arcs The table containing the outgoing arcs -- of the current vertex -- -- @return The average distance of the new node #k to any subtree -- The distances are stored as follows: -- example: distances[center][a] -- center is any vertex, thus if center is an inner vertex -- it has 3 neighbours a,b and c, which can all be seen as the -- roots of subtrees A,B,C. -- distances[center][a] gives us the distance of the new -- node k to the subtree A. -- if center is a leaf, it has only one neighbour, which -- can also be seen as the root of the subtree T\{center} -- function BalancedMinimumEvolution:computeAverageDistancesToAllSubtreesForK(vertex, values, k, k_dists) local is_leaf = self.is_leaf local arcs = self.tree.arcs local vertices = self.tree.vertices local center_vertex = vertex -- for every vertex a table is created, in which the distances to all -- its subtrees will be stored values.outgoing_arcs = values.outgoing_arcs or self.tree:outgoing(center_vertex) for _, arc in ipairs (values.outgoing_arcs) do local root = arc.head -- this vertex can be seen as the root of a subtree if is_leaf[root] then -- we know the distance of k to the leaf! k_dists[center_vertex][root] = self:distance(vertices[k], root) else -- to compute the distance we need the root's neighbouring vertices, which we can access by its outgoing arcs local arc1, arc2 local arc_back -- the arc we came from for _, next_arc in ipairs (self.tree:outgoing(root)) do if next_arc.head ~= center_vertex then arc1 = arc1 or next_arc arc2 = next_arc else arc_back = next_arc end end values.outgoing_arcs = { arc1, arc2, arc_back } -- go deeper, if the distances for the next center node haven't been set yet if not (k_dists[root][arc1.head] and k_dists[root][arc2.head]) then self:computeAverageDistancesToAllSubtreesForK(root, values, k,k_dists) end -- set the distance between k and subtree k_dists[center_vertex][root] = 1/2 * (k_dists[root][arc1.head] + k_dists[root][arc2.head]) end end end -- -- Sets the distances from k to subtrees -- In computeAverageDistancesToAllSubtreesForK the distances to ALL possbile -- subtrees are computed. Once k is inserted many of those subtrees don't -- exist for k, as k is now part of them. In this function all -- still accurate subtrees and their distances to k are -- extracted. -- -- @param center The vertex serving as the starting point of the depth-first search; -- should be the new_node function BalancedMinimumEvolution:setAccurateDistancesForK(center,visited,k,k_dists,leaves) local visited = visited or {} local distances = self.distances visited[center] = center local outgoings = self.tree:outgoing(center) for _,arc in ipairs (outgoings) do local vertex = arc.head if vertex ~= leaves[k] then local distance -- set the distance if not distances[leaves[k]][vertex] and k_dists[center] then distance = k_dists[center][vertex] -- use previously calculated distance distances[leaves[k]][vertex] = distance distances[vertex][leaves[k]] = distance end -- go deeper if not visited[vertex] then self:setAccurateDistancesForK(vertex,visited,k,k_dists,leaves) end end end end -- -- Find the best edge for the insertion of leaf #k, such that the -- total tree length is minimized. This function uses a depth first -- search. -- -- @param vertex The vertex where the depth first search is -- started; must be a leaf -- @param values The values needed for the recursion -- - visited: The vertices that already have been visited -- - tree_length: The current tree_length -- - best_arc: The current best_arc, such that the tree -- length is minimzed -- - min_length: The smallest tree_length found so far function BalancedMinimumEvolution:findBestEdge(vertex, values, k_dists) local arcs = self.tree.arcs local vertices = self.tree.vertices values = values or { visited = {} } values.visited[vertex] = vertex local c -- the arc we came from local unvisited_arcs = {} --unvisited arcs --identify arcs for _, arc in ipairs (self.tree:outgoing(vertex)) do if not values.visited[arc.head] then unvisited_arcs[#unvisited_arcs+1] = arc else c = arc.head --last visited arc end end for i, arc in ipairs (unvisited_arcs) do local change_in_tree_length = 0 -- set tree length to 0 for first insertion arc if not values.tree_length then values.tree_length = 0 values.best_arc = arc values.min_length = 0 else -- compute new tree length for the case that k is inserted into this arc local b = arc.head --current arc local a = unvisited_arcs[i%2+1].head -- the remaining arc local k_v = vertices[k] -- the leaf to be inserted change_in_tree_length = 1/4 * ( ( self:distance(a,c) + k_dists[vertex][b]) - (self:distance(a,b) + k_dists[vertex][c]) ) values.tree_length = values.tree_length + change_in_tree_length end -- if the tree length becomes shorter, this is the new best arc -- for the insertion of leaf k if values.tree_length < values.min_length then values.best_arc = arc values.min_length = values.tree_length end -- go deeper self:findBestEdge(arc.head, values, k_dists) values.tree_length = values.tree_length - change_in_tree_length end return values.best_arc end -- Calculates the total tree length -- This is done by adding up all the edge lengths -- -- @return the tree length function BalancedMinimumEvolution:calculateTreeLength() local vertices = self.tree.vertices local sum = 0 for index, v1 in ipairs(vertices) do for i = index+1,#vertices do local v2 = vertices[i] local dist = self.lengths[v1][v2] if dist then sum = sum + dist end end end return sum end -- generates edges for the final graph -- -- throughout the process of creating the tree, arcs have been -- disconnected and connected, without truly creating edges. this is -- done in this function function BalancedMinimumEvolution:createFinalEdges() local g = self.tree local o_arcs = {} -- copy arcs since createEdge is going to modify the arcs array... for _,arc in ipairs(g.arcs) do if arc.tail.event.index < arc.head.event.index then o_arcs[#o_arcs+1] = arc end end for _,arc in ipairs(o_arcs) do InterfaceToAlgorithms.createEdge( self.main_algorithm, arc.tail, arc.head, { generated_options = { { key = "phylogenetic edge", value = tostring(self.lengths[arc.tail][arc.head]) } }}) end end -- Gets the distance between two nodes as specified in their options -- or storage fields. -- Note: this function implies 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 BalancedMinimumEvolution:distance(a, b) if a == b then return 0 else local distances = self.distances return distances[a][b] or distances[b][a] end end -- -- computes the final branch lengths -- -- goes over all arcs and computes the final branch lengths, -- as neither the BME nor the BNNI main_algorithm does so. function BalancedMinimumEvolution:computeFinalLengths() local is_leaf = self.is_leaf local lengths = self.lengths local g = self.tree for _, arc in ipairs(g.arcs) do local head = arc.head local tail = arc.tail local distance local a,b,c,d -- assert, that the length hasn't already been computed for this arc if not lengths[head][tail] then if not is_leaf[head] then -- define subtrees a and b for _, arc in ipairs (g:outgoing(head)) do local subtree = arc.head if subtree ~= tail then a = a or subtree b = subtree end end end if not is_leaf[tail] then -- define subtrees c and d for _, arc in ipairs (g:outgoing(tail)) do local subtree = arc.head if subtree ~= head then c = c or subtree d = subtree end end end -- compute the distance using the formula for outer or inner edges, respectively if is_leaf[head] then distance = 1/2 * ( self:distance(head,c) + self:distance(head,d) - self:distance(c,d) ) elseif is_leaf[tail] then distance = 1/2 * ( self:distance(tail,a) + self:distance(tail,b) - self:distance(a,b) ) else --inner edge distance = self:distance(head, tail) -1/2 * ( self:distance(a,b) + self:distance(c,d) ) end lengths[head][tail] = distance lengths[tail][head] = distance end end end return BalancedMinimumEvolution