summaryrefslogtreecommitdiff
path: root/Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics
diff options
context:
space:
mode:
authorKarl Berry <karl@freefriends.org>2014-02-27 00:01:50 +0000
committerKarl Berry <karl@freefriends.org>2014-02-27 00:01:50 +0000
commit90ae413a94da014dbf7a6a8ab584f7c668483a5a (patch)
tree4e0528a2b2fb2004988e7fc804ea2f1354f1f5d5 /Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics
parent800cc56380d1edfe567a6a66bef13e8950eb2b0c (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')
-rw-r--r--Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/AuthorDefinedPhylogeny.lua77
-rw-r--r--Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/BalancedMinimumEvolution.lua593
-rw-r--r--Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/BalancedNearestNeighbourInterchange.lua372
-rw-r--r--Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/DistanceMatrix.lua440
-rw-r--r--Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/Maeusle2012.lua778
-rw-r--r--Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/PhylogeneticTree.lua91
-rw-r--r--Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/SokalMichener1958.lua261
-rw-r--r--Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/library.lua35
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"
+
+