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