diff options
author | Karl Berry <karl@freefriends.org> | 2014-02-27 00:01:50 +0000 |
---|---|---|
committer | Karl Berry <karl@freefriends.org> | 2014-02-27 00:01:50 +0000 |
commit | 90ae413a94da014dbf7a6a8ab584f7c668483a5a (patch) | |
tree | 4e0528a2b2fb2004988e7fc804ea2f1354f1f5d5 /Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics | |
parent | 800cc56380d1edfe567a6a66bef13e8950eb2b0c (diff) |
pgf 3.0.0
git-svn-id: svn://tug.org/texlive/trunk@33057 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics')
8 files changed, 2647 insertions, 0 deletions
diff --git a/Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/AuthorDefinedPhylogeny.lua b/Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/AuthorDefinedPhylogeny.lua new file mode 100644 index 00000000000..dc5dea45b0f --- /dev/null +++ b/Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/AuthorDefinedPhylogeny.lua @@ -0,0 +1,77 @@ +-- Copyright 2013 by 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/AuthorDefinedPhylogeny.lua,v 1.2 2013/03/20 17:16:38 tantau Exp $ + + + + +local AuthorDefinedPhylogeny = {} + + +-- Namespace +require("pgf.gd.phylogenetics").AuthorDefinedPhylogeny = AuthorDefinedPhylogeny + +-- Imports +local InterfaceToAlgorithms = require "pgf.gd.interface.InterfaceToAlgorithms" +local Direct = require "pgf.gd.lib.Direct" + +-- Shorthand: +local declare = InterfaceToAlgorithms.declare + + +--- +declare { + key = "phylogenetic tree by author", + algorithm = AuthorDefinedPhylogeny, + phase = "phylogenetic tree generation", + phase_default = true, + + summary = [[" + When this key is used, the phylogenetic tree must be specified + by the author (rather than being generated algorithmically). + "]], + documentation = [[" + A spanning tree of the input graph will be computed first (it + must be connected, otherwise errors will result). + The evolutionary length of the edges must be specified through + the use of the |length| key for each edge. + "]], + examples = [[" + \tikz \graph [phylogenetic tree layout] { + a -- { + b [>length=2] --[length=1] { c, d }, + e [>length=3] + } + }; + "]] +} + + + +function AuthorDefinedPhylogeny:run() + + local spanning_tree = self.main_algorithm.digraph.options.algorithm_phases["spanning tree computation"].new { + ugraph = self.main_algorithm.ugraph, + events = {} -- no events + }:run() + + local phylogenetic_tree = Direct.ugraphFromDigraph(spanning_tree) + local lengths = self.lengths + + for _,a in ipairs(phylogenetic_tree.arcs) do + lengths[a.tail][a.head] = a:options('length') + end + + return phylogenetic_tree +end + + + +return AuthorDefinedPhylogeny diff --git a/Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/BalancedMinimumEvolution.lua b/Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/BalancedMinimumEvolution.lua new file mode 100644 index 00000000000..fd0fb62b3ad --- /dev/null +++ b/Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/BalancedMinimumEvolution.lua @@ -0,0 +1,593 @@ +-- 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/BalancedMinimumEvolution.lua,v 1.1 2013/02/08 17:14:05 tantau Exp $ + + + +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 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 diff --git a/Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/DistanceMatrix.lua b/Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/DistanceMatrix.lua new file mode 100644 index 00000000000..0bca7773212 --- /dev/null +++ b/Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/DistanceMatrix.lua @@ -0,0 +1,440 @@ +-- Copyright 2013 by 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/DistanceMatrix.lua,v 1.2 2013/12/20 14:44:47 tantau Exp $ + + + + +local DistanceMatrix = {} + + +-- Imports +local InterfaceToAlgorithms = require("pgf.gd.interface.InterfaceToAlgorithms") +local declare = InterfaceToAlgorithms.declare + + +--- + +declare { + key = "distance matrix vertices", + type = "string", + + summary = [[" + A list of vertices that are used in the parsing of the + |distance matrix| key. If this key is not used at all, all + vertices of the graph will be used for the computation of a + distance matrix. + "]], + + documentation = [[" + The vertices must be separated by spaces and/or + commas. For vertices containing spaces or commas, the vertex + names may be surrounded by single or double quotes (as in + Lua). Typical examples are |a, b, c| or |"hello world", 'foo'|. + "]] +} + + + +--- + +declare { + key = "distance matrix", + type = "string", + + summary = [[" + A distance matrix specifies ``desired distances'' between + vertices in a graph. These distances are used, in particular, in + algorithms for computing phylogenetic trees. + "]], + + documentation = [[" + When this key is parsed, the key |distance matrix vertices| is + considered first. It is used to determine a list of vertices + for which a distance matrix is computed, see that key for + details. Let $n$ be the number of vertices derived from that + key. + + The string passed to the |distance matrix| key is basically + a sequence of numbers that are used to fill an $n \times n$ + matrix. This works as follows: We keep track of a \emph{current + position $p$} in the matrix, starting at the upper left corner + of the matrix. We read the numbers in the string + one by one, write it to the current position of the matrix, and + advance the current position by going right one step; if we go + past the right end of the matrix, we ``wrap around'' by going + back to the left border of the matrix, but one line down. If we + go past the bottom of the matrix, we start at the beginning once + more. + + This basic behaviour can be modified in different ways. First, + when a number is followed by a semicolon instead of a comma or a + space (which are the ``usual'' ways of indicating the end of a + number), we immediately go down to the next line. Second, + instead of a number you can directly provide a \emph{position} + in the matrix and the current position will be set to this + position. Such a position information is detected by a + greater-than sign (|>|). It must be followed by + \begin{itemize} + \item a number or a vertex name or + \item a number or a vertex name, a comma, and another number or + vertex name or + \item a comma and a number and a vertex name. + \end{itemize} + Examples of the respective cases are |>1|, |>a,b|, and + |>,5|. The semantics is as follows: In all cases, if a vertex + name rather than a number is given, it is converted into a + number (namely the index of the vertex inside the matrix). Then, + in the first case, the column of the current position is set to + the given number; in the second case, the columns is set to the + first number and the column is set to the second number; and in + the third case only the row is set to the given number. (This + idea is that following the |>|-sign comes a ``coordinate pair'' + whose components are separated by a comma, but part of that pair + may be missing.) If a vertex name contains special symbols like + a space or a comma, you must surround it by single or double + quotation marks (as in Lua). + + Once the string has been parsed completely, the matrix may be + filled only partially. In this case, for each missing entry + $(x,y)$, we try to set it to the value of the entry $(y,x)$, + provided that entry is set. If neither are set, the entry is set + to $0$. + + Let us now have a look at several examples that all produce the + same matrix. The vertices are |a|, |b|, |c|. +\begin{codeexample}[code only, tikz syntax=false] +0, 1, 2 +1, 0, 3 +2, 3, 0 +\end{codeexample} +\begin{codeexample}[code only, tikz syntax=false] +0 1 2 1 0 3 2 3 0 +\end{codeexample} +\begin{codeexample}[code only, tikz syntax=false] +; +1; +2 3 +\end{codeexample} +\begin{codeexample}[code only, tikz syntax=false] +>,b 1; 2 3 +\end{codeexample} +\begin{codeexample}[code only, tikz syntax=false] +>b 1 2 >c 3 +\end{codeexample} + "]] +} + + +--- + +declare { + key = "distances", + type = "string", + + summary = [[" + This key is used to specify the ``desired distances'' between + a vertex and the other vertices in a graph. + "]], + + documentation = [[" + This key works similar to the |distance matrix| key, only it is + passed to a vertex instead of to a whole graph. The syntax is + the same, only the notion of different ``rows'' is not + used. Here are some examples that all have the same effect, + provided the nodes are |a|, |b|, and |c|. +\begin{codeexample}[code only, tikz syntax=false] +0, 1, 2 +\end{codeexample} +\begin{codeexample}[code only, tikz syntax=false] +0 1 2 +\end{codeexample} +\begin{codeexample}[code only, tikz syntax=false] +>b 1 2 +\end{codeexample} +\begin{codeexample}[code only, tikz syntax=false] +>c 2, >b 1 +\end{codeexample} + "]] +} + + + +local function to_index(s, indices) + if s and s ~= "" then + if s:sub(1,1) == '"' then + local _, _, m = s:find('"(.*)"') + return indices[InterfaceToAlgorithms.findVertexByName(m)] + elseif s:sub(1,1) == "'" then + local _, _, m = s:find("'(.*)'") + return indices[InterfaceToAlgorithms.findVertexByName(m)] + else + local num = tonumber(s) + if not num then + return indices[InterfaceToAlgorithms.findVertexByName(s)] + else + return num + end + end + end +end + +local function compute_indices(vertex_string, vertices) + local indices = {} + + if not vertex_string then + for i,v in ipairs(vertices) do + indices[i] = v + indices[v] = i + end + else + -- Ok, need to parse the vertex_string. Sigh. + local pos = 1 + while pos <= #vertex_string do + local start = vertex_string:sub(pos,pos) + if not start:find("[%s,]") then + local _, vertex + if start == '"' then + _, pos, vertex = vertex_string:find('"(.-)"', pos) + elseif start == "'" then + _, pos, vertex = vertex_string:find("'(.-)'", pos) + else + _, pos, vertex = vertex_string:find("([^,%s'\"]*)", pos) + end + local v = assert(InterfaceToAlgorithms.findVertexByName(vertex), "unknown vertex name '" .. vertex .. "'") + indices [#indices + 1] = v + indices [v] = #indices + end + pos = pos + 1 + end + end + + return indices +end + + +--- +-- Compute a distance matrix based on the values of a +-- |distance matrix| and a |distance matrix vertices|. +-- +-- @param matrix_string A distance matrix string +-- @param vertex_string A distance matrix vertex string +-- @param vertices An array of all vertices in the graph. +-- +-- @return A distance matrix. This matrix will contain both a +-- two-dimensional array (accessed through numbers) and also a +-- two-dimensional hash table (accessed through vertex indices). Thus, +-- you can write both |m[1][1]| and also |m[v][v]| to access the first +-- entry of this matrix, provided |v == vertices[1]|. +-- @return An index vector. This is an array of the vertices +-- identified for the |vertex_string| parameter. + +function DistanceMatrix.computeDistanceMatrix(matrix_string, vertex_string, vertices) + -- First, we create a table of the vertices we need to consider: + local indices = compute_indices(vertex_string, vertices) + + -- Second, build matrix. + local n = #indices + local m = {} + for i=1,n do + m[i] = {} + end + + local x = 1 + local y = 1 + local pos = 1 + -- Start scanning the matrix_string + while pos <= #matrix_string do + local start = matrix_string:sub(pos,pos) + if not start:find("[%s,]") then + if start == '>' then + local _, parse + _, pos, parse = matrix_string:find(">([^%s>;]*)", pos) + local a, b + if parse:find(",") then + _,_,a,b = parse:find("(.*),(.*)") + else + a = parse + end + x = to_index(a, indices) or x + y = to_index(b, indices) or y + elseif start == ';' then + x = 1 + y = y + 1 + elseif start == ',' then + x = x + 1 + else + local _, n + _, pos, n = matrix_string:find("([^,;%s>]*)", pos) + local num = assert(tonumber(n), "number expected in distance matrix") + m[x][y] = num + x = x + 1 + -- Skip everything up to first comma: + _, pos = matrix_string:find("(%s*,?)", pos+1) + end + end + pos = pos + 1 + if x > n then + x = 1 + y = y + 1 + end + if y > n then + y = 1 + end + end + + -- Fill up + for x=1,n do + for y=1,n do + if not m[x][y] then + m[x][y] = m[y][x] or 0 + end + end + end + + -- Copy to index version + for x=1,n do + local v = indices[x] + m[v] = {} + for y=1,n do + local u = indices[y] + m[v][u] = m[x][y] + end + end + + return m, indices +end + + + + +--- +-- Compute a distance vector. See the key |distances| for details. +-- +-- @param vector_string A distance vector string +-- @param vertex_string A distance matrix vertex string +-- @param vertices An array of all vertices in the graph. +-- +-- @return A distance vector. Like a distance matrix, this vector will +-- double indexed, once by numbers and once be vertex objects. +-- @return An index vector. This is an array of the vertices +-- identified for the |vertex_string| parameter. + +function DistanceMatrix.computeDistanceVector(vector_string, vertex_string, vertices) + -- First, we create a table of the vertices we need to consider: + local indices = compute_indices(vertex_string, vertices) + + -- Second, build matrix. + local n = #indices + local m = {} + local x = 1 + local pos = 1 + -- Start scanning the vector_string + while pos <= #vector_string do + local start = vector_string:sub(pos,pos) + if not start:find("[%s,]") then + if start == '>' then + local _, parse + _, pos, parse = vector_string:find(">([^%s>;]*)", pos) + x = to_index(parse, indices) or x + elseif start == ',' then + x = x + 1 + else + local _, n + _, pos, n = vector_string:find("([^,;%s>]*)", pos) + local num = assert(tonumber(n), "number expected in distance matrix") + m[x] = num + x = x + 1 + -- Skip everything up to first comma: + _, pos = vector_string:find("(%s*,?)", pos+1) + end + end + pos = pos + 1 + if x > n then + x = 1 + end + end + + -- Fill up + for x=1,n do + m[x] = m[x] or 0 + m[indices[x]] = m[x] + end + + return m, indices +end + + + +--- +-- Compute a distance matrix for a graph that incorporates all +-- information stored in the different options of the graph and the +-- vertices. +-- +-- @param graph A digraph object. +-- +-- @return A distance matrix for all vertices of the graph. + +function DistanceMatrix.graphDistanceMatrix(digraph) + local vertices = digraph.vertices + local n = #vertices + local m = {} + for i,v in ipairs(vertices) do + m[i] = {} + m[v] = {} + end + + local indices = {} + for i,v in ipairs(vertices) do + indices[i] = v + indices[v] = i + end + + if digraph.options['distance matrix'] then + local sub, vers = DistanceMatrix.computeDistanceMatrix( + digraph.options['distance matrix'], + digraph.options['distance matrix vertices'], + vertices + ) + + for x=1,#vers do + for y=1,#vers do + m[vers[x]][vers[y]] = sub[x][y] + end + end + end + + for i,v in ipairs(vertices) do + if v.options['distances'] then + local sub, vers = DistanceMatrix.computeDistanceVector( + v.options['distances'], + v.options['distance matrix vertices'], + vertices + ) + + for x=1,#vers do + m[vers[x]][v] = sub[x] + end + end + end + + -- Fill up number versions: + for x,vx in ipairs(vertices) do + for y,vy in ipairs(vertices) do + m[x][y] = m[vx][vy] + end + end + + return m +end + + + +return DistanceMatrix diff --git a/Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/Maeusle2012.lua b/Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/Maeusle2012.lua new file mode 100644 index 00000000000..ef0909a6d00 --- /dev/null +++ b/Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/Maeusle2012.lua @@ -0,0 +1,778 @@ +-- 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/Maeusle2012.lua,v 1.4 2013/05/23 20:01:29 tantau Exp $ + + +-- Imports +local Digraph = require 'pgf.gd.model.Digraph' +local Coordinate = require 'pgf.gd.model.Coordinate' +local Path = require 'pgf.gd.model.Path' + +local layered = require 'pgf.gd.layered' + +local lib = require 'pgf.gd.lib' + +local declare = require("pgf.gd.interface.InterfaceToAlgorithms").declare + + + +-- Main class of this file: + +local Maeusle2012 = lib.class {} + +-- Namespace +require("pgf.gd.phylogenetics").Maeusle2012 = Maeusle2012 + + + + +--- +declare { + key = "rooted rectangular phylogram", + algorithm = { + base_class = Maeusle2012, + run = function (self) + local root = self:getRoot() + self:setPosForRectangularLayout(root) + end + }, + phase = "phylogenetic tree layout", + phase_default = true, + + summary = [[" + A rooted rectangular phylogram is... + "]], + documentation = [[" + ... + "]], + examples = [[" + \tikz \graph [phylogenetic tree layout, + rooted rectangular phylogram, + balanced minimum evolution, + 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 }; + "]] +} + +--- +declare { + key = "rectangular phylogram", + use = { { key = "rooted rectangular phylogram" } }, + summary = "An alias for |rooted rectangular phylogram|" +} + +--- +declare { + key = "rooted straight phylogram", + algorithm = { + base_class = Maeusle2012, + run = function (self) + local root = self:getRoot() + self:setXPos(root) + self:setYPosForStraightLayout(root) + end + }, + phase = "phylogenetic tree layout", + + summary = [[" + A rooted straight phylogram is... + "]], + documentation = [[" + ... + "]], + examples = [[" + \tikz \graph [phylogenetic tree layout, + rooted straight phylogram, + balanced minimum evolution, grow=right, + 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 }; + "]]} + +--- +declare { + key = "straight phylogram", + use = { { key = "rooted straight phylogram" } }, + summary = "An alias for |rooted straight phylogram|" +} + +--- +declare { + key = "unrooted rectangular phylogram", + algorithm = { + base_class = Maeusle2012, + run = function (self) + local root1, root2 = self:getRoot() + self:setPosForUnrootedRectangular(root2, root1) + end + }, + phase = "phylogenetic tree layout", + + summary = [[" + A unrooted rectangular phylogram is... + "]], + documentation = [[" + ... + "]], + examples = [[" + \tikz \graph [phylogenetic tree layout, + unrooted rectangular phylogram, + balanced minimum evolution, grow=right, + 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 }; + "]] +} + +--- +declare { + key = "unrooted straight phylogram", + algorithm = { + base_class = Maeusle2012, + run = function (self) + local root1, root2 = self:getRoot() + self:setPosForUnrootedStraight(root2, root1) + end + }, + phase = "phylogenetic tree layout", + + summary = [[" + A unrooted straight phylogram is... + "]], + documentation = [[" + ... + "]], + examples = [[" + \tikz \graph [phylogenetic tree layout, + unrooted straight phylogram, + balanced minimum evolution, grow=right, + 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 }; + "]] +} + + +--- +declare { + key = "evolutionary unit length", + type = "length", + initial = "1cm", + + summary = [[" + Specifies how long a ``unit'' of evolutionary time should be on + paper. For instance, if two nodes in a phylogenetic tree have an + evolutionary distance of 3 and this length is set to |1cm|, then + they will be |3cm| apart in a straight-line phylogram. + "]], + documentation = [[" + (This key used to be called |distance scaling factor|.) + "]], +} + + + +-- +-- Gets the edge length between two nodes +-- +-- @param vertex1, vertex2 The two nodes +-- +-- @return The length of the edge between the two nodes +function Maeusle2012:edgeLength(vertex1, vertex2) + return self.lengths[vertex1][vertex2] +end + + +-- Sets the x and y coordinates for all nodes, using a depth first +-- search +-- +-- @param vertex The starting point; should usually be the root +-- @param values Values needed for the recursion +-- @param vertex2 A node that will not be visited; this parameter should only be set +-- for an unrooted layout to ensure that only the half of the tree is set. +function Maeusle2012:setPosForRectangularLayout(vertex, values, vertex2) + local arcs = self.tree.arcs + local vertices = self.tree.vertices + local adjusted_bb = self.main_algorithm.adjusted_bb + + values = values or { + length = 0, -- current path length + visited = {}, -- all nodes that have already been visited + leaves = {}, -- all leaves from left to right + } + + local vertex_is_leaf = true + values.visited[vertex] = true + + local children = {} -- a table containing all children of the + --current vertex (for the later determination of inner vertices + --x-positions) + + + for _, arc in ipairs (self.tree:outgoing(vertex)) do + if not values.visited[arc.head] and arc.head ~= vertex2 then + -- if arc.head hasn't been visited, the current vertex cannot be a leaf + vertex_is_leaf = false + local arc_length = self:edgeLength(vertex, arc.head) + + values.length = values.length + arc_length + + -- go deeper + self:setPosForRectangularLayout(arc.head, values, vertex2) + + -- get the children of the current vertex + children[#children+1] = arc.head + + values.length = values.length - arc_length + end + end + + if vertex_is_leaf then + -- subtract layer_pre, thus the leaf itself is NOT part of the + -- edge length + vertex.pos.y = - adjusted_bb[vertex].layer_pre + + values.leaves[#values.leaves+1] = vertex + + -- x coordinate: + -- the x coordinates of the leaves are the first to be set; the + -- first leave stays at x = 0, the x coordinates for the other + -- leaves is computed with help of the ideal_sibling_distance + -- function + if #values.leaves > 1 then + local left_sibling = values.leaves[#values.leaves-1] + local ideal_distance = layered.ideal_sibling_distance(adjusted_bb, self.tree, vertex, left_sibling ) + vertex.pos.x = left_sibling.pos.x + ideal_distance + end + + else -- the vertex is an inner node + -- the x position of an inner vertex is at the center of its children. + + -- determine the outer children + local left_child = children[1] + local right_child = left_child + for _, child in ipairs(children) do + if child.pos.x < left_child.pos.x then left_child = child end + if child.pos.x > right_child.pos.x then right_child = child end + end + + -- position between child with highest and child with lowest x-value, + -- if number of children is even + local index_of_middle_child = math.ceil(#children/2) + local even = #children/2 == index_of_middle_child + + if even then + vertex.pos.x = (left_child.pos.x + right_child.pos.x) / 2 + index_of_middle_child = 0 + else -- if number of children is odd, position above the middle child + vertex.pos.x = children[index_of_middle_child].pos.x + table.remove(children, index_of_middle_child) -- don't bend the edge to this node, as it it above it anyway + end + end + + -- set the node's y-coordinate, using the calculated length + -- and a scaling factor + vertex.pos.y = vertex.pos.y + (values.length * self.tree.options['evolutionary unit length']) + + -- if this is the second subtree to be set of an unrooted tree, have + -- it grow in the other direction + if values.second_subtree then + vertex.pos.y = -vertex.pos.y + end + + -- bend the edges for the rectangular layout + for i,child in ipairs(children) do + self:bendEdge90Degree(child, vertex) + end + + return values +end + + +-- Sets only the x-positions of all nodes using a depth-first search. +-- This is necessary for straight-edge layouts. +-- +-- @param vertex The starting point of the depth-first search; should usually be the root +-- @param values Values needed for the recursion +-- @param vertex2 A node that will not be visited; this parameter should only be set +-- for an unrooted layout to ensure that only the half of the tree is set. +function Maeusle2012:setXPos(vertex, values, vertex2) + local arcs = self.tree.arcs + local vertices = self.tree.vertices + if not values then + values = { + visited = {}, -- all nodes that have already been visited + leaves = {}, -- all leaves from left to right + } + end + + local vertex_is_leaf = true + values.visited[vertex] = true + local children = {} -- a table containing all children of the current vertex (for the later determination of inner vertices x-positions) + + for _, arc in ipairs (self.tree:outgoing(vertex)) do + if not values.visited[arc.head] and arc.head ~= vertex2 then + -- if arc.head hasn't been visited, the current vertex cannot be a leaf + vertex_is_leaf = false + + -- go deeper + self:setXPos(arc.head, values, vertex2) + + -- get the children of the current vertex + table.insert(children, arc.head) + end + end + + -- set the x-position of a leaf + if vertex_is_leaf then + + table.insert(values.leaves, vertex) + + if #values.leaves > 1 then + local left_sibling = values.leaves[#values.leaves-1] + local ideal_distance = layered.ideal_sibling_distance(self.main_algorithm.adjusted_bb, self.tree, vertex, left_sibling ) + vertex.pos.x = left_sibling.pos.x + ideal_distance + end + + -- set x position of an inner node, which is at the center of its + -- children + else + -- determine the outer children + local left_child = children[1] + local right_child = left_child + for _, child in ipairs(children) do + if child.pos.x < left_child.pos.x then left_child = child end + if child.pos.x > right_child.pos.x then right_child = child end + end + + -- position between child with highest and child with lowest x-value, + -- if number of children is even + local index_of_middle_child = math.ceil(#children/2) + local even = #children/2 == index_of_middle_child + + if even then + vertex.pos.x = (left_child.pos.x + right_child.pos.x) / 2 + else -- if number of children is odd, position above the middle child + vertex.pos.x = children[index_of_middle_child].pos.x + end + end + return values +end + + +-- +-- Sets only the y-positions of all nodes using a depth-first search. +-- This is needed for a straight-edge layout, as the x-positions have +-- to bet first so that the y-coordinates can be calculated correctly +-- here. +-- +-- @param vertex1 The starting point of the depth-first search +-- @param values Values needed for the recursion +-- @param vertex2 For unrooted layout only: The root of the second subtree. +-- This node and all its children will not be visited. +function Maeusle2012:setYPosForStraightLayout(vertex, values, vertex2) + local arcs = self.tree.arcs + local vertices = self.tree.vertices + local adjusted_bb = self.main_algorithm.adjusted_bb + + values = values or { + length = 0, -- current path length + visited = {}, -- all nodes that have already been visited + leaves = {}, -- all leaves from left to right + } + + local vertex_is_leaf = true + values.visited[vertex] = true + local children = {} -- a table containing all children of the current vertex (for the later determination of inner vertices x-positions) + + for _, arc in ipairs (self.tree:outgoing(vertex)) do + if not values.visited[arc.head] and arc.head ~= vertex2 then + -- if arc.head hasn't been visited, the current vertex cannot be a leaf + vertex_is_leaf = false + + -- calculate the arc length with the help of the Pythagorean + -- theorem + local a + local l = self:edgeLength(vertex, arc.head) * self.tree.options['evolutionary unit length'] + local b = math.abs(vertex.pos.x - arc.head.pos.x) + if b > l then + a = 0 + else + a = math.sqrt(l^2-b^2) + end + local arc_length = a + + + values.length = values.length + arc_length + + -- go deeper + self:setYPosForStraightLayout(arc.head, values, vertex2) + + -- get the children of the current vertex + table.insert(children, arc.head) + + values.length = values.length - arc_length + end + end + + if vertex_is_leaf then + -- subtract layer_pre, thus the leaf itself is NOT part of the + -- edge length + vertex.pos.y = - adjusted_bb[vertex].layer_pre + + table.insert(values.leaves, vertex) + end + + -- set the node's y-coordinate, using the calculated length + vertex.pos.y = vertex.pos.y + values.length + + -- if this is the second subtree to be set of an unrooted tree, have + -- it grow in the other direction + if values.second_subtree then vertex.pos.y = -vertex.pos.y end +end + +-- +-- Correct the x-positions in the unrooted layout for a more aesthetic result +-- +-- If the roots of the two subtrees have different x-positions, this is corrected +-- by shifting the x-positions of all nodes in one subtree by that difference. +-- +-- @param vertex1 The root of the first subtree +-- @param vertex2 The root of the second subtree. +function Maeusle2012:correctXPos(vertex1, vertex2, straight) + + -- correct the x-positions + -- + -- @param vertex Starting point of the depth-first search + -- @param values Values needed for the recursion + -- @param vertex2 The root of the subtree that will not be visited + local function x_correction(vertex, values, vertex2) + values.visited[vertex] = true + local children = {} + + for _, arc in ipairs (self.tree:outgoing(vertex)) do + if not values.visited[arc.head] and arc.head ~= vertex2 then + + table.insert(children, arc.head) + x_correction(arc.head, values, vertex2) + end + end + + vertex.pos.x = vertex.pos.x + values.diff + if not straight then + for i,child in ipairs(children) do + self:bendEdge90Degree(child, vertex) + end + end + + return values + end + + -- compute the difference of the x-positions of the two subtrees' + -- roots + local diff = vertex1.pos.x - vertex2.pos.x + local values = { visited = {} } + if diff < 0 then + values.diff = - diff + x_correction(vertex1, values, vertex2) + elseif diff > 0 then + values.diff = diff + x_correction(vertex2, values, vertex1) + end +end + + +-- +-- Sets the x- and y-positions of the vertices in an unrooted layout +-- +-- This is done using the function for setting the positions for a rooted layout: +-- Two neighboring vertices are chosen as roots; one half of the tree +-- is drawn in one direction, the other half 180° to the other +-- direction. +-- +-- @param vertex1, vertex2: The vertices functioning as roots +function Maeusle2012:setPosForUnrootedRectangular(vertex1, vertex2) + -- set positions for first half of the tree... + self:setPosForRectangularLayout(vertex2,false,vertex1) + local vals={ + length = self:edgeLength(vertex1, vertex2), -- the length between the two roots + visited={}, + leaves={}, + path={}, + second_subtree = true + } + -- ... and for the second half. + self:setPosForRectangularLayout(vertex1,vals,vertex2) + -- if the two roots have different x-values, correct the x-positions for nicer layout + self:correctXPos(vertex1, vertex2, false) +end + + +-- +-- Sets the x- and y-positions of the vertices in an unrooted straight layout +-- +-- This is done using the function for setting the positions for a rooted straight layout: +-- Two neighboring vertices are chosen as roots; one half of the tree +-- is drawn in one direction, the other half 180° to the other +-- direction. +-- +-- @param vertex1, vertex2: The vertices functioning as roots +function Maeusle2012:setPosForUnrootedStraight(vertex1, vertex2) + -- first set the x-positions of the two subtrees... + local vals = {visited = {}, leaves = {} } + self:setXPos(vertex2, vals, vertex1) + self:setXPos(vertex1, vals, vertex2) + + -- ... and then the y-positions + self:setYPosForStraightLayout(vertex2, false, vertex1) + local vals={ + length = self:edgeLength(vertex1, vertex2) * self.tree.options['evolutionary unit length'], + visited={}, + leaves={}, + path={}, + second_subtree = true + } + self:setYPosForStraightLayout(vertex1, vals, vertex2) + + -- if the two roots have different x-values, correct the x-positions for nicer layout + -- as the length between the roots of the two subtrees is set to the calculated value, + -- this step is mandatory for the unrooted, straight layout + self:correctXPos(vertex1, vertex2, true) +end + + + +-- Bends the arc between two nodes by 90 degree by updating the arc's +-- path +-- +-- @param head The head of the arc +-- @param tail The tail of the arc +function Maeusle2012:bendEdge90Degree(head, tail) + local arc = self.tree:arc(tail,head) + local syntactic_tail = arc:syntacticTailAndHead() + arc:setPolylinePath { Coordinate.new(head.pos.x, tail.pos.y) } +end + + + +-- Finds the longest path in a graph +-- +-- @ return A table containing the path (an array of nodes) and the +-- path length +function Maeusle2012:findLongestPath() + local starting_point = self.tree.vertices[1] -- begin at any vertex + -- get the path lengths from the starting point to all leaves: + local paths_to_leaves = self:getPathLengthsToLeaves(starting_point) + local path_lengths = paths_to_leaves.path_lengths + local paths = paths_to_leaves.paths + + -- looks for the longest path and identifies its end-point + local function find_head_of_longest_path(path_lengths, paths) + local longest_path + local node + -- to make sure that the same path is chosen every time, we go over all vertices with "ipairs"; if we would go over path_lengths directly, we could only use "pairs" + for _, vertex in ipairs(self.tree.vertices) do + local path_length = path_lengths[vertex] + if path_length then + -- choose longest path. if two paths have the same length, take the path with more nodes + if not longest_path or path_length > longest_path or (path_length == longest_path and #paths[vertex]>#paths[node]) then + longest_path = path_length + node = vertex + end + end + end + return node + end + + -- find the longest path leading away from the starting point and identify + -- the leaf it leads to. Use that leaf as the tail for the next path + -- search + local tail = find_head_of_longest_path(path_lengths, paths) + paths_to_leaves = self:getPathLengthsToLeaves(tail) -- gets new path information + -- paths_to leaves now has all paths starting at vertex "tail"; one of these paths is the + -- longest (globally) + path_lengths = paths_to_leaves.path_lengths + paths = paths_to_leaves.paths + local head = find_head_of_longest_path(path_lengths, paths) + + local path_information = + { path = paths_to_leaves.paths[head], -- longest path + length = path_lengths[head] } -- length of that path + + return path_information +end + + +-- a depth first search for getting all path lengths from a +-- starting point to all leaves +-- +-- @param vertex The vertex where the search is to start +-- @param values Table of values needed for the recursive computation +-- +-- @return A table containing: +-- a table of the leaves with corresponding path lengths +-- and a table containing the path to each leaf (an array of +-- nodes) +function Maeusle2012:getPathLengthsToLeaves(vertex, values) + local arcs = self.tree.arcs + local vertices = self.tree.vertices + if not values then + values = { + paths = {}, -- all paths we've found so far + path_lengths = {}, -- all path lengths that have so far been computed + length = 0, -- current path length + visited = {}, -- all nodes that have already been visited + path = {}, -- the current path we're on + leaves = {} -- all leaves from left to right + } + table.insert(values.path,vertex) + end + + local vertex_is_leaf = true + values.visited[vertex] = true + + for _, arc in ipairs (self.tree:outgoing(vertex)) do + if not values.visited[arc.head] then + -- the current vertex is not a leaf! note: if the starting vertex is a leaf, vertex_is_leaf + -- will be set to 'false' for it anyway. as we're not interested in the distance + -- of the starting vertex to itself, this is fine. + vertex_is_leaf = false + local arc_length = self.lengths[vertex][arc.head] + values.length = values.length + arc_length + + -- add arc.head to path... + table.insert(values.path,arc.head) + + -- ... and go down that path + self:getPathLengthsToLeaves(arc.head, values) + + -- remove arc.head again to go a different path + table.remove(values.path) + values.length = values.length - arc_length + end + end + + if vertex_is_leaf then -- we store the information gained on the path to this leaf + values.path_lengths[vertex] = values.length + values.paths[vertex] = {} + table.insert(values.leaves, vertex) + for i,k in pairs(values.path) do + values.paths[vertex][i] = k + end + end + -- the path_lengths and the paths are stored in one table and + -- returned together + local path_information = + { path_lengths = values.path_lengths, + paths = values.paths, + leaves = values.leaves } + return path_information +end + + +-- Gets the root of a tree +-- checks whether a tree is already rooted, if not, computeCenterOfPath() is +-- called, which defines a node in the center of the graph as the root +-- +-- @return The root +function Maeusle2012:getRoot() + -- check whether a root exists (vertex with degree 2) + local root = lib.find (self.tree.vertices, function(v) return #self.tree:outgoing(v) == 2 end) + if root then + return root, self.tree:outgoing(root)[1].head + else + return self:computeCenterOfPath() + end +end + + +-- +-- @return The newly computed root and its nearest neighbor +function Maeusle2012:computeCenterOfPath() + local longest_path = self:findLongestPath() + local path = longest_path.path + local root, neighbor_of_root + + local length = 0 --length between first vertex on the path and the current vertex we're looking at + for i = 1, #path-1 do + local node1 = path[i] + local node2 = path[i+1] + local node3 = path[i+2] + + local dist_node_1_2, dist_node_2_3 --distances between node1 and node2, and node2 and node3 + dist_node_1_2 = self:edgeLength(node1, node2) + if node3 then dist_node_2_3 = self:edgeLength(node2, node3) end + length = length + dist_node_1_2 -- length between first vertex on the path and current node2 + + if length == longest_path.length/2 then + root = node2 -- if there is a node exactly at the half of the path, use this node as root + + -- and find nearest neighbor of the root + if node3 == nil or dist_node_1_2 < dist_node_2_3 then -- neu 3.8 + neighbor_of_root = node1 + else + neighbor_of_root = node3 + end + break + + elseif length > longest_path.length/2 then + -- else find node closest to the center of the path and use it as the root; + local node2_length = math.abs(longest_path.length/2 - length) + local node1_length = math.abs(longest_path.length/2 - (length - dist_node_1_2)) + if node2_length < node1_length then + root = node2 + neighbor_of_root = node1 + -- if node3 is closer to node2 than node1 is, use node3 as neighbor! + if node3 and dist_node_2_3 < dist_node_1_2 then neighbor_of_root = node3 end + else + root = node1 + neighbor_of_root = node2 + --check if node i-1 is closer to node1 + local dist_node_0_1 + if i>1 then + node0 = path[i-1] + dist_node_0_1 = self:edgeLength(node0, node1) + if dist_node_0_1 < dist_node_1_2 then neighbor_of_root = node0 end + end + end + break + end + end + + return root, neighbor_of_root +end + + +return Maeusle2012 diff --git a/Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/PhylogeneticTree.lua b/Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/PhylogeneticTree.lua new file mode 100644 index 00000000000..16356aa7cfa --- /dev/null +++ b/Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/PhylogeneticTree.lua @@ -0,0 +1,91 @@ +-- 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/PhylogeneticTree.lua,v 1.2 2013/02/08 17:14:06 tantau Exp $ + + +local PhylogeneticTree = {} + + +-- Namespace +require("pgf.gd.phylogenetics").PhylogeneticTree = PhylogeneticTree + +-- Imports +local InterfaceToAlgorithms = require "pgf.gd.interface.InterfaceToAlgorithms" +local Storage = require "pgf.gd.lib.Storage" +local Direct = require "pgf.gd.lib.Direct" + +-- Shorthand: +local declare = InterfaceToAlgorithms.declare + + +--- +declare { + key = "phylogenetic tree layout", + algorithm = PhylogeneticTree, + + postconditions = { + upward_oriented = true + }, + + summary = [[" + Layout for drawing phylogenetic trees. + "]], + documentation = [[" + ... + "]], + examples = [[" + \tikz \graph [phylogenetic tree layout, upgma, + 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 }; + "]] +} + + +-- Computes a phylogenetic tree and/or visualizes it +-- - computes a phylogenetic tree according to what the "phylogenetic +-- algorithm" key is set to +-- - invokes a graph drawing algorithm according to what the +-- "phylogenetic layout" key is set to +function PhylogeneticTree:run() + + local options = self.digraph.options + + -- Two storages for some information computed by the phylogenetic + -- tree generation algorithm + local lengths = Storage.newTableStorage() + + -- First, compute the phylogenetic tree + local tree = options.algorithm_phases['phylogenetic tree generation'].new { + main_algorithm = self, + lengths = lengths + }:run() + + tree = Direct.ugraphFromDigraph(tree) + + -- Second, layout the tree + local layout_class = options.algorithm_phases['phylogenetic tree layout'] + layout_class.new { + main_algorithm = self, + distances = distances, + lengths = lengths, + tree = tree + }:run() + + tree:sync() +end + +return PhylogeneticTree diff --git a/Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/SokalMichener1958.lua b/Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/SokalMichener1958.lua new file mode 100644 index 00000000000..c55a7a7467c --- /dev/null +++ b/Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/SokalMichener1958.lua @@ -0,0 +1,261 @@ +-- 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/SokalMichener1958.lua,v 1.1 2013/02/08 17:14:06 tantau Exp $ + + + + +local SokalMichener1958 = {} + + +-- Namespace +require("pgf.gd.phylogenetics").SokalMichener1958 = SokalMichener1958 + +-- Imports +local InterfaceToAlgorithms = require("pgf.gd.interface.InterfaceToAlgorithms") +local DistanceMatrix = require("pgf.gd.phylogenetics.DistanceMatrix") +local lib = require("pgf.gd.lib") +local Storage = require("pgf.gd.lib.Storage") +local Digraph = require("pgf.gd.model.Digraph") + +-- Shorthand: +local declare = InterfaceToAlgorithms.declare + + +--- +declare { + key = "unweighted pair group method using arithmetic averages", + algorithm = SokalMichener1958, + phase = "phylogenetic tree generation", + + summary = [[" + The UPGMA (Unweighted Pair Group Method using arithmetic + Averages) algorithm of Sokal and Michener, 1958. It generates a + graph on the basis of such a distance matrix by generating nodes + and computing the edge lengths. + "]], + documentation = [[" + This algorithm uses a distance matrix, ideally an ultrametric + one, to compute the graph. + "]], + examples = [[" + \tikz \graph [phylogenetic tree layout, sibling distance=0pt, sibling sep=2pt, + unweighted pair group method using arithmetic averages, + 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 }; + "]] +} + + +--- +declare { + key = "upgma", + use = { { key = "unweighted pair group method using arithmetic averages" } }, + summary = "An shorthand for |unweighted pair group method using arithmetic averages|" +} + + + + +-- +-- The run function of the upgma algorithm. +-- +-- You must setup the following fields: The |main_algorithm| must +-- store the main algorithm object (for phase |main|). The |distances| +-- field must be a |Storage| object that will get filled with the +-- distances computed by this algorithm. The |lengths| field must also +-- be a |Storage| for the computed distances. +-- + +function SokalMichener1958:run() + self.distances = Storage.newTableStorage() + + self.tree = Digraph.new(self.main_algorithm.digraph) + + -- store the phylogentic tree object, containing all user-specified + -- graph information + self:runUPGMA() + self:createFinalEdges() + + return self.tree +end + + + +-- UPGMA (Unweighted Pair Group Method using arithmetic Averages) algorithm +-- (Sokal and Michener, 1958) +-- +-- this function generates a graph on the basis of such a distance +-- matrix by generating nodes and computing the edge lengths; the x- +-- and y-positions of the nodes must be set separately +-- +-- requirement: a distance matrix, ideally an ultrametric +function SokalMichener1958:runUPGMA() + local matrix = DistanceMatrix.graphDistanceMatrix(self.tree) + + local g = self.tree + local clusters = {} + + -- create the clusters + for _,v in ipairs(g.vertices) do + clusters[#clusters+1] = self:newCluster(v) + end + + -- Iniitialise the distances of these clusters: + for _,cx in ipairs(clusters) do + for _,cy in ipairs(clusters) do + cx.distances[cy] = matrix[cx.root][cy.root] + end + end + + -- search for clusters with smallest distance and merge them + while #clusters > 1 do + local minimum_distance = math.huge + local min_cluster1 + local min_cluster2 + for i, cluster in ipairs (clusters) do + for j = i+1,#clusters do + local cluster2 = clusters[j] + local cluster_distance = self:getClusterDistance(cluster, cluster2) + if cluster_distance < minimum_distance then + minimum_distance, min_cluster1, min_cluster2 = cluster_distance, i, j + end + end + end + self:mergeClusters(clusters, min_cluster1, min_cluster2, minimum_distance) + end +end + + +-- a new cluster is created +-- +-- @param vertex The vertex the cluster is initialized with +-- +-- @return The new cluster +function SokalMichener1958:newCluster(vertex) + return { + root = vertex, -- the root of the cluster + size = 1, -- the number of vertices in the cluster, + distances = {}, -- cached cluster distances to all other clusters + cluster_height = 0 -- this value is equivalent to half the distance of the last two clusters + -- that have been merged to form the current cluster; + -- necessary for determining the distances of newly generated nodes to their children. + } +end + + +-- gets the distance between two clusters +-- +-- @param cluster1, cluster2 The two clusters +-- +-- @return the distance between the clusters +function SokalMichener1958:getClusterDistance(c,d) + return c.distances[d] or d.distances[c] or 0 +end + + +-- merges two clusters by doing the following: +-- - deletes cluster2 from the clusters table +-- - adds all vertices from cluster2 to the vertices table of cluster1 +-- - updates the distances of the new cluster to all remaining clusters +-- - generates a new node, as the new root of the cluster +-- - computes the distance of the new node to the former roots (for +-- later computation of the y-positions) +-- - generates edges, connecting the new node to the former roots +-- - updates the cluster height +-- +-- @param clusters The array of clusters +-- @param index_of_first_cluster The index of the first cluster +-- @param index_of_second_cluster The index of the second cluster +-- @param distance The distance between the two clusters + +function SokalMichener1958:mergeClusters(clusters, index_of_first_cluster, index_of_second_cluster, distance) + + local g = self.tree + local cluster1 = clusters[index_of_first_cluster] + local cluster2 = clusters[index_of_second_cluster] + + --update cluster distances + for i,cluster in ipairs (clusters) do + if cluster ~= cluster1 and cluster ~= cluster2 then + local dist1 = self:getClusterDistance (cluster1, cluster) + local dist2 = self:getClusterDistance (cluster2, cluster) + local dist = (dist1*cluster1.size + dist2*cluster2.size)/ (cluster1.size+cluster2.size) + cluster1.distances[cluster] = dist + cluster.distances[cluster1] = dist + end + end + + -- delete cluster2 + table.remove(clusters, index_of_second_cluster) + + --add node and connect last vertex of each cluster with new node + local new_node = InterfaceToAlgorithms.createVertex( + self.main_algorithm, + { + name = "UPGMA-node ".. #self.tree.vertices+1, + generated_options = { { key = "phylogenetic inner node" } }, + } + ) + g:add{new_node} + -- the distance of the new node ( = the new root of the cluster) to its children (= the former roots) is + -- equivalent to half the distance between the two former clusters + -- minus the respective cluster height + local distance1 = distance/2-cluster1.cluster_height + self.distances[new_node][cluster1.root] = distance1 + local distance2 = distance/2-cluster2.cluster_height + self.distances[new_node][cluster2.root] = distance2 + + -- these distances are also the final edge lengths, thus: + self.lengths[new_node][cluster1.root] = distance1 + self.lengths[cluster1.root][new_node] = distance1 + + self.lengths[new_node][cluster2.root] = distance2 + self.lengths[cluster2.root][new_node] = distance2 + + g:connect(new_node, cluster1.root) + g:connect(new_node, cluster2.root) + + cluster1.root = new_node + cluster1.size = cluster1.size + cluster2.size + cluster1.cluster_height = distance/2 -- set new height of the cluster +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 SokalMichener1958: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 + o_arcs[#o_arcs+1] = arc + 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 + + +return SokalMichener1958 diff --git a/Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/library.lua b/Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/library.lua new file mode 100644 index 00000000000..194e809ef2e --- /dev/null +++ b/Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/library.lua @@ -0,0 +1,35 @@ +-- Copyright 2012 by 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/library.lua,v 1.2 2013/04/04 20:43:45 tantau Exp $ + + + +--- +-- A phylogenetic tree (or network) depicts the evolutionary histroy +-- of species or, more generally, so called taxa. The present library +-- includes a number of algorithms for drawing phylogenetic trees. +-- +-- @library + +local phylogenetics -- Library name + +-- Main layout: +require "pgf.gd.phylogenetics.PhylogeneticTree" + +-- Phylogenetic tree drawing: +require "pgf.gd.phylogenetics.Maeusle2012" + +-- Phylogenetic tree generation: +require "pgf.gd.phylogenetics.SokalMichener1958" +require "pgf.gd.phylogenetics.BalancedMinimumEvolution" +require "pgf.gd.phylogenetics.BalancedNearestNeighbourInterchange" +require "pgf.gd.phylogenetics.AuthorDefinedPhylogeny" + + |