summaryrefslogtreecommitdiff
path: root/Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/Maeusle2012.lua
diff options
context:
space:
mode:
Diffstat (limited to 'Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/Maeusle2012.lua')
-rw-r--r--Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/Maeusle2012.lua778
1 files changed, 778 insertions, 0 deletions
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