summaryrefslogtreecommitdiff
path: root/Master/texmf-dist/tex/generic/pgf/graphdrawing/lua/pgf/gd/phylogenetics/SokalMichener1958.lua
blob: 8be5c1d51e93f952eeea404bc2462e67c468ed43 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
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$




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