1# Copyright 2014, Sandia Corporation. Under the terms of Contract 2# DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain 3# rights in this software. 4 5"""Provides layout algorithms. 6""" 7 8 9import time 10 11import custom_inherit 12import numpy 13 14import toyplot.units 15 16def region( 17 xmin, 18 xmax, 19 ymin, 20 ymax, 21 bounds=None, 22 rect=None, 23 corner=None, 24 grid=None, 25 margin=None): 26 """Specify a rectangular target region relative to a parent region. 27 28 Parameters 29 ---------- 30 xmin: :class:`number<numbers.Number>`, required 31 Minimum X boundary of the parent region, specified in CSS pixel units. 32 xmax: :class:`number<numbers.Number>`, required 33 Maximum X boundary of the parent region, specified in CSS pixel units. 34 ymin: :class:`number<numbers.Number>`, required 35 Minimum Y boundary of the parent region, specified in CSS pixel units. 36 ymax: :class:`number<numbers.Number>`, required 37 Maximum Y boundary of the parent region, specified in CSS pixel units. 38 margin: :class:`number<numbers.Number>`, string, (:class:`number<numbers.Number>`, string) tuple, or tuple containing between one and four :class:`numbers<numbers.Number>`, strings, or (:class:`number<numbers.Number>`, string) tuples, optional 39 Padding around the target region, specified in real-world units. Defaults 40 to CSS pixel units. See :ref:`units` for details. Follows the same behavior as the CSS margin property. 41 42 Returns 43 ------- 44 xmin, xmax, ymin, ymax: :class:`number<numbers.Number>` 45 The boundaries of the target region, specified in CSS pixel units. 46 """ 47 48 if margin is None: 49 margin = 0 50 51 if isinstance(margin, tuple): 52 if len(margin) == 4: 53 margin_top = toyplot.units.convert(margin[0], "px", default="px") 54 margin_right = toyplot.units.convert(margin[1], "px", default="px") 55 margin_bottom = toyplot.units.convert(margin[2], "px", default="px") 56 margin_left = toyplot.units.convert(margin[3], "px", default="px") 57 elif len(margin) == 3: 58 margin_top = toyplot.units.convert(margin[0], "px", default="px") 59 margin_left = margin_right = toyplot.units.convert(margin[1], "px", default="px") 60 margin_bottom = toyplot.units.convert(margin[2], "px", default="px") 61 elif len(margin) == 2: 62 margin_top = margin_bottom = toyplot.units.convert(margin[0], "px", default="px") 63 margin_left = margin_right = toyplot.units.convert(margin[1], "px", default="px") 64 elif len(margin) == 1: 65 margin_top = margin_bottom = margin_left = margin_right = toyplot.units.convert(margin[0], "px", default="px") 66 else: 67 margin_top = margin_bottom = margin_left = margin_right = toyplot.units.convert(margin, "px", default="px") 68 69 def convert(vmin, vmax, value): 70 value = toyplot.units.convert( 71 value, "px", default="px", reference=vmax - vmin) 72 if value < 0: 73 return float(vmax + value) 74 return float(vmin + value) 75 76 # Specify explicit bounds for the region 77 if bounds is not None: 78 if isinstance(bounds, tuple) and len(bounds) == 4: 79 return ( 80 convert(xmin, xmax, bounds[0]), 81 convert(xmin, xmax, bounds[1]), 82 convert(ymin, ymax, bounds[2]), 83 convert(ymin, ymax, bounds[3]), 84 ) 85 raise TypeError("bounds parameter must be an (xmin, xmax, ymin, ymax) tuple.") # pragma: no cover 86 # Specify an explicit rectangle for the region 87 if rect is not None: 88 if isinstance(rect, tuple) and len(rect) == 4: 89 x = convert(xmin, xmax, rect[0]) 90 y = convert(ymin, ymax, rect[1]) 91 width = toyplot.units.convert( 92 rect[2], "px", default="px", reference=xmax - xmin) 93 height = toyplot.units.convert( 94 rect[3], "px", default="px", reference=ymax - ymin) 95 return (x, x + width, y, y + height) 96 raise ValueError("Unrecognized rect type") # pragma: no cover 97 # Offset a rectangle from a corner 98 if corner is not None: 99 if isinstance(corner, tuple) and len(corner) == 4: 100 position = corner[0] 101 inset = toyplot.units.convert(corner[1], "px", default="px") 102 width = toyplot.units.convert( 103 corner[2], "px", default="px", reference=xmax - xmin) 104 height = toyplot.units.convert( 105 corner[3], "px", default="px", reference=ymax - ymin) 106 else: 107 raise ValueError("corner parameter must be a (corner, inset, width, height) tuple.") # pragma: no cover 108 if position == "top": 109 return ((xmin + xmax - width) / 2, 110 (xmin + xmax + width) / 2, 111 ymin + inset, 112 ymin + inset + height) 113 elif position == "top-right": 114 return ( 115 xmax - 116 width - 117 inset, 118 xmax - 119 inset, 120 ymin + 121 inset, 122 ymin + 123 inset + 124 height) 125 elif position == "right": 126 return ( 127 xmax - width - inset, 128 xmax - inset, 129 (ymin + ymax - height) / 2, 130 (ymin + ymax + height) / 2) 131 elif position == "bottom-right": 132 return ( 133 xmax - 134 width - 135 inset, 136 xmax - 137 inset, 138 ymax - 139 inset - 140 height, 141 ymax - 142 inset) 143 elif position == "bottom": 144 return ((xmin + xmax - width) / 2, 145 (xmin + xmax + width) / 2, 146 ymax - inset - height, 147 ymax - inset) 148 elif position == "bottom-left": 149 return ( 150 xmin + 151 inset, 152 xmin + 153 inset + 154 width, 155 ymax - 156 inset - 157 height, 158 ymax - 159 inset) 160 elif position == "left": 161 return ( 162 xmin + inset, 163 xmin + inset + width, 164 (ymin + ymax - height) / 2, 165 (ymin + ymax + height) / 2) 166 elif position == "top-left": 167 return ( 168 xmin + 169 inset, 170 xmin + 171 inset + 172 width, 173 ymin + 174 inset, 175 ymin + 176 inset + 177 height) 178 else: 179 raise ValueError("Unrecognized corner") # pragma: no cover 180 181 # Choose a cell from an MxN grid, with optional column/row spanning. 182 if grid is not None: 183 # Cell n (in left-to-right, top-to-bottom order) of an M x N grid 184 if len(grid) == 3: 185 M, N, n = grid 186 i = n // N 187 j = n % N 188 rowspan, colspan = (1, 1) 189 elif len(grid) == 4: # Cell i,j of an M x N grid 190 M, N, i, j = grid 191 rowspan, colspan = (1, 1) 192 # Cells [i, i+rowspan), [j, j+colspan) of an M x N grid 193 elif len(grid) == 6: 194 M, N, i, rowspan, j, colspan = grid 195 196 cell_width = (xmax - xmin) / N 197 cell_height = (ymax - ymin) / M 198 199 return ( 200 (j * cell_width) + margin_left, 201 ((j + colspan) * cell_width) - margin_right, 202 (i * cell_height) + margin_top, 203 ((i + rowspan) * cell_height) - margin_bottom, 204 ) 205 # If nothing else fits, consume the entire region 206 return (xmin + margin_left, xmax - margin_right, ymin + margin_top, ymax - margin_bottom) 207 208 209class Graph(object): 210 """Stores graph layout information. 211 212 Typically used when sharing a layout among more than one graph. 213 """ 214 def __init__(self, vids, vcoordinates, edges, eshapes, ecoordinates): 215 self._vids = vids 216 self._vcoordinates = vcoordinates 217 self._edges = edges 218 self._eshapes = eshapes 219 self._ecoordinates = ecoordinates 220 221 @property 222 def vcount(self): 223 """Return the number of vertices :math:`V` in the graph.""" 224 return len(self._vids) 225 226 @property 227 def vids(self): 228 """Return :math:`V` graph vertex identifiers.""" 229 return self._vids 230 231 @property 232 def vcoordinates(self): 233 """Return graph vertex coordinates as a :math:`V \\times 2` matrix.""" 234 return self._vcoordinates 235 236 @property 237 def ecount(self): 238 """Return the number of edges :math:`E` in the graph.""" 239 return len(self._edges) 240 241 @property 242 def eshapes(self): 243 """Return a vector of :math:`E` string edge shapes.""" 244 return self._eshapes 245 246 @property 247 def ecoordinates(self): 248 """Return a matrix of edge coordinates. The number of coordinates is determined by the contents of the `eshapes` vector.""" 249 return self._ecoordinates 250 251 @property 252 def edges(self): 253 """Return the graph edges as a :math:`E \\times 2` matrix of source, target indices.""" 254 return self._edges 255 256def graph(a, b=None, c=None, olayout=None, layout=None, vcoordinates=None): 257 """Compute a graph layout.""" 258 259 stack = [c, b, a] 260 261 # Identify the graph's edges. 262 if isinstance(stack[-1], numpy.ndarray) and stack[-1].ndim == 2 and stack[-1].shape[1] == 2: 263 edges = stack.pop() 264 ecount = edges.shape[0] 265 else: 266 sources = toyplot.require.vector(stack.pop()) 267 ecount = len(sources) 268 targets = toyplot.require.vector(stack.pop(), ecount) 269 edges = numpy.column_stack((sources, targets)) 270 271 # Identify the vertex ids induced by the edges. 272 vids, edges = numpy.unique(edges, return_inverse=True) 273 edges = edges.reshape((ecount, 2), order="C") 274 275 # If the caller supplied extra vertex ids, merge them in. 276 if stack and stack[-1] is not None: 277 extra_vids = toyplot.require.vector(stack.pop()) 278 disconnected_vids = numpy.setdiff1d(extra_vids, vids) 279 vids = numpy.append(vids, disconnected_vids) 280 281 # Setup storage to receive vertex coordinates 282 vcount = len(vids) 283 ivcoordinates = numpy.ma.masked_all((vcount, 2)) 284 285 # If the caller supplied the layout for an external graph, merge those coordinates in. 286 if olayout is not None: 287 olayout = toyplot.require.instance(olayout, (toyplot.mark.Graph, toyplot.layout.Graph)) 288 oindices = numpy.in1d(olayout.vids, vids, assume_unique=True) 289 iindices = numpy.in1d(vids, olayout.vids, assume_unique=True) 290 ivcoordinates[iindices] = olayout.vcoordinates[oindices] 291 292 # If the caller supplied extra vertex coordinates, merge them in. 293 if vcoordinates is not None: 294 external_vcoordinates = toyplot.require.scalar_matrix(vcoordinates, rows=vcount, columns=2) 295 mask = numpy.ma.getmaskarray(external_vcoordinates) 296 ivcoordinates = numpy.ma.where(mask, ivcoordinates, external_vcoordinates) 297 298 # Apply the layout algorithm to whatever's left. 299 start = time.time() 300 if layout is None: 301 # If there are unspecified coordinates, use a force-directed layout. 302 if numpy.ma.is_masked(ivcoordinates): 303 layout = toyplot.layout.FruchtermanReingold() 304 else: 305 # Otherwise, we can ignore the vertices and just create edges. 306 layout = toyplot.layout.IgnoreVertices() 307 vcoordinates, eshapes, ecoordinates = layout.graph(ivcoordinates, edges) 308 toyplot.log.info("Graph layout time: %s ms", (time.time() - start) * 1000) 309 310 if numpy.ma.is_masked(vcoordinates): 311 raise RuntimeError("Graph layout cannot return masked vertex coordinates.") # pragma: no cover 312 313 return Graph(vids, vcoordinates, edges, eshapes, ecoordinates) 314 315def _add_at(target, target_indices, source): 316 """Add source values to the target and handle duplicate indices correctly. 317 318 With numpy, the expression `target[indices] += source` does not work intuitively 319 if there are duplicate indices. This function handles this case as you would 320 expect, by accumulating multiple values for a single target. 321 """ 322 if getattr(numpy.add, "at", None) is not None: 323 numpy.add.at(target, target_indices, source) 324 else: # Shim for numpy < 1.8 325 for source_index, target_index in enumerate(target_indices): # pragma: no cover 326 target[target_index] += source[source_index] 327 328 329#def _floyd_warshall_shortest_path(vcount, edges): 330# """Compute the (directed) shortest paths between every pair of vertices in a graph, using the Floyd-Warshall algorithm. 331# 332# Floyd-Warshall is a good choice for computing paths between all pairs of 333# vertices in dense graphs. 334# 335# Note that this algorithm is directed, i.e. a path from i -> j doesn't imply 336# a path from j -> i. The results will contain numpy.inf for vertices that 337# have no path. 338# 339# Returns 340# ------- 341# shortest_paths: :math:`V \\times V matrix` 342# A matrix where element :math:`E_ij` contains the shortest path distance 343# between vertex :math:`i` and vertex :math:`j`. 344# """ 345# distance = numpy.empty((vcount, vcount)) 346# distance[...] = numpy.inf 347# distance[numpy.diag_indices_from(distance)] = 0 348# distance[edges.T[0], edges.T[1]] = 1 349# for k in numpy.arange(vcount): 350# for i in numpy.arange(vcount): 351# for j in numpy.arange(vcount): 352# if distance[i, j] > distance[i, k] + distance[k, j]: 353# distance[i, j] = distance[i, k] + distance[k, j] 354# 355# return distance 356 357 358def _adjacency_list(vcount, edges): 359 """Return an adjacency list representation of a graph.""" 360 targets = [[] for i in numpy.arange(vcount)] 361 for source, target in edges: 362 targets[source].append(target) 363 return targets 364 365 366def _require_tree(children): 367 """Return the root vertex and maximum depth of a tree. 368 369 Parameters 370 ---------- 371 children: list of lists 372 Adjacency list representation of a graph. 373 374 Returns 375 ------- 376 root: integer 377 Index of the root vertex. 378 depth: integer 379 Maximum depth of the tree. 380 """ 381 roots = numpy.setdiff1d(numpy.arange(len(children)), numpy.concatenate(children)) 382 if len(roots) != 1: 383 raise ValueError("Not a tree.") # pragma: no cover 384 root = roots[0] 385 386 depth = [] 387 visited = numpy.zeros(len(children), dtype=bool) 388 def mark_visited(vertex, vdepth=0): 389 if visited[vertex]: 390 raise ValueError("Not a tree.") # pragma: no cover 391 depth.append(vdepth) 392 visited[vertex] = True 393 for child in children[vertex]: 394 mark_visited(child, vdepth + 1) 395 mark_visited(root) 396 397 return root, max(depth) 398 399class EdgeLayout(object, metaclass=custom_inherit.DocInheritMeta(style="numpy_napoleon")): 400 """Abstract interface for algorithms that compute graph edge coordinates.""" 401 def edges(self, vcoordinates, edges): 402 """Return edge coordinates for a graph. 403 404 Parameters 405 ---------- 406 vcoordinates : :math:`V \\times 2` matrix 407 Contains the coordinates for every graph vertex, in vertex order. 408 edges : :math:`E \\times 2` matrix 409 Contains the integer vertex indices for every graph edge in edge 410 order. The first and second matrix columns contain the source and 411 target vertices respectively. 412 413 Returns 414 ------- 415 eshapes : array of :math:`E` strings 416 Contains a shape string for each edge, in edge order. The shape 417 string contains drawing codes that define an arbitrary-complexity 418 path for the edge, using a set of current coordinates and a turtle 419 drawing model. The following codes are currently allowed: 420 421 * `M` - change the current coordinates without drawing (requires one set of coordinates). 422 * `L` - draw a straight line segment from the current coordinates (requires one set of coordinates). 423 * `Q` - draw a quadratic Bezier curve from the current coordinates (requires two sets of coordinates). 424 * `C` - draw a cubic Bezier curve from the current coordinates (requires three sets of coordinates). 425 ecoordinates : matrix containing two columns 426 Contains coordinates for each of the edge shape strings, in drawing-code order. 427 """ 428 raise NotImplementedError() # pragma: no cover 429 430 431class StraightEdges(EdgeLayout): 432 """Creates straight edges between graph vertices.""" 433 def edges(self, vcoordinates, edges): 434 loops = edges.T[0] == edges.T[1] 435 if numpy.any(loops): 436 toyplot.log.warning("Graph contains %s loop edges that will not be visible.", numpy.count_nonzero(loops)) 437 438 eshapes = numpy.tile("ML", len(edges)) 439 ecoordinates = numpy.empty((len(edges) * 2, 2)) 440 ecoordinates[0::2] = vcoordinates[edges.T[0]] 441 ecoordinates[1::2] = vcoordinates[edges.T[1]] 442 return eshapes, ecoordinates 443 444 445class CurvedEdges(EdgeLayout): 446 """Creates curved edges between graph vertices. 447 448 Parameters 449 ---------- 450 curvature : number 451 Controls the curvature of each edge's arc, as a percentage of the 452 distance between the edge's vertices. Use negative values to curve 453 edges in the opposite direction. 454 """ 455 def __init__(self, curvature=0.15): 456 self._curvature = curvature 457 458 def edges(self, vcoordinates, edges): 459 loops = edges.T[0] == edges.T[1] 460 if numpy.any(loops): 461 toyplot.log.warning("Graph contains %s loop edges that will not be visible.", numpy.count_nonzero(loops)) 462 463 eshapes = numpy.tile("MQ", len(edges)) 464 ecoordinates = numpy.empty((len(edges) * 3, 2)) 465 466 sources = vcoordinates[edges.T[0]] 467 targets = vcoordinates[edges.T[1]] 468 offsets = numpy.dot(targets - sources, [[0, 1], [-1, 0]]) * self._curvature 469 midpoints = ((sources + targets) * 0.5) + offsets 470 471 ecoordinates[0::3] = sources 472 ecoordinates[1::3] = midpoints 473 ecoordinates[2::3] = targets 474 475 return eshapes, ecoordinates 476 477 478class GraphLayout(object, metaclass=custom_inherit.DocInheritMeta(style="numpy_napoleon")): 479 """Abstract interface for algorithms that compute coordinates for graph vertices and edges.""" 480 def graph(self, vcoordinates, edges): 481 """Compute vertex and edge coordinates for a graph. 482 483 Parameters 484 ---------- 485 vcoordinates : :math:`V \\times 2` masked array 486 Coordinates for every graph vertex, in vertex order. Where 487 practical, only masked coordinates will have values assigned by the 488 underlying algorithm. 489 edges : :math:`E \\times 2` matrix 490 Contains the integer vertex indices for every graph edge in edge 491 order. The first and second matrix columns contain the source and 492 target vertices respectively. 493 494 Returns 495 ------- 496 vcoordinates : :math:`V \\times 2` matrix 497 Contains coordinates for every graph vertex, in vertex order. 498 eshapes : array of :math:`E` strings 499 Contains a shape string for each edge, in edge order. The shape 500 string contains drawing codes that define an arbitrary-complexity 501 path for the edge, using a set of current coordinates and a turtle 502 drawing model. The following codes are currently allowed: 503 504 * `M` - change the current coordinates without drawing (requires one set of coordinates). 505 * `L` - draw a straight line segment (requires one set of coordinates). 506 * `Q` - draw a quadratic Bezier curve (requires two sets of coordinates). 507 * `C` - draw a cubic Bezier curve (requires three sets of coordinates). 508 ecoordinates : matrix containing two columns 509 Contains coordinates for each of the edge shape strings, in drawing-code order. 510 """ 511 raise NotImplementedError() # pragma: no cover 512 513 514class IgnoreVertices(GraphLayout): 515 """Do-nothing graph layout for use when all vertices are already specified. 516 517 Parameters 518 ---------- 519 edges: :class:`toyplot.layout.EdgeLayout` instance, optional 520 The default will generate straight edges. 521 """ 522 def __init__(self, edges=None): 523 if edges is None: 524 edges = StraightEdges() 525 526 self._edges = edges 527 528 def graph(self, vcoordinates, edges): 529 eshapes, ecoordinates = self._edges.edges(vcoordinates, edges) 530 return vcoordinates, eshapes, ecoordinates 531 532 533class Random(GraphLayout): 534 """Compute a random graph layout. 535 536 Parameters 537 ---------- 538 edges: :class:`toyplot.layout.EdgeLayout` instance, optional 539 The default will generate straight edges. 540 seed: integer, optional 541 Random seed used to generate vertex coordinates. 542 """ 543 def __init__(self, edges=None, seed=1234): 544 if edges is None: 545 edges = StraightEdges() 546 547 self._edges = edges 548 self._seed = seed 549 550 def graph(self, vcoordinates, edges): 551 generator = numpy.random.RandomState(seed=self._seed) 552 mask = numpy.ma.getmaskarray(vcoordinates) 553 vcoordinates = numpy.ma.where(mask, generator.uniform(-1, 1, size=vcoordinates.shape), vcoordinates) 554 eshapes, ecoordinates = self._edges.edges(vcoordinates, edges) 555 return vcoordinates, eshapes, ecoordinates 556 557 558class Eades(GraphLayout): 559 """Compute a force directed graph layout using the 1984 algorithm of Eades. 560 561 Parameters 562 ---------- 563 edges: :class:`toyplot.layout.EdgeLayout` instance, optional 564 The default will generate straight edges. 565 c1, c2, c3, c4: numbers, optional 566 Constants defined in Eades' original paper. 567 M: integer, optional 568 Number of iterations to run the spring simulation. 569 seed: integer, optional 570 Random seed used to initialize vertex coordinates. 571 """ 572 def __init__(self, edges=None, c1=2, c2=1, c3=1, c4=0.1, M=100, seed=1234): 573 if edges is None: 574 edges = StraightEdges() 575 576 self._edges = edges 577 self._c1 = c1 578 self._c2 = c2 579 self._c3 = c3 580 self._c4 = c4 581 self._M = M 582 self._seed = seed 583 584 def graph(self, vcoordinates, edges): 585 generator = numpy.random.RandomState(seed=self._seed) 586 # Initialize coordinates 587 mask = numpy.ma.getmaskarray(vcoordinates) 588 vcoordinates = numpy.ma.where(mask, generator.uniform(-1, 1, size=vcoordinates.shape), vcoordinates) 589 590 # Repeatedly apply attract / repel forces to the vertices 591 vertices = numpy.column_stack(numpy.triu_indices(n=len(vcoordinates), k=1)) 592 for iteration in numpy.arange(self._M): 593 offsets = numpy.zeros_like(vcoordinates) 594 595 # Repel 596 a = vcoordinates[vertices.T[0]] 597 b = vcoordinates[vertices.T[1]] 598 delta = a - b 599 distance = numpy.linalg.norm(delta, axis=1)[:, None] 600 delta /= distance 601 force = self._c3 / numpy.square(distance) 602 delta *= force 603 _add_at(offsets, vertices.T[0], delta) 604 _add_at(offsets, vertices.T[1], -delta) 605 606 # Attract 607 a = vcoordinates[edges.T[0]] 608 b = vcoordinates[edges.T[1]] 609 delta = b - a 610 distance = numpy.linalg.norm(delta, axis=1)[:, None] 611 delta /= distance 612 force = self._c1 * numpy.log(distance / self._c2) 613 delta *= force 614 _add_at(offsets, edges.T[0], delta) 615 _add_at(offsets, edges.T[1], -delta) 616 617 # Sum offsets 618 vcoordinates = numpy.ma.where(mask, vcoordinates + self._c4 * offsets, vcoordinates) 619 620 eshapes, ecoordinates = self._edges.edges(vcoordinates, edges) 621 return vcoordinates, eshapes, ecoordinates 622 623 624class FruchtermanReingold(GraphLayout): 625 """Compute a force directed graph layout using the 1991 algorithm of Fruchterman and Reingold. 626 627 Parameters 628 ---------- 629 edges: :class:`toyplot.layout.EdgeLayout` instance, optional 630 The default will generate straight edges. 631 area, temperature: numbers, optional 632 Constants defined in the original paper. 633 M: integer, optional 634 Number of iterations to run the spring simulation. 635 seed: integer, optional 636 Random seed used to initialize vertex coordinates. 637 """ 638 def __init__(self, edges=None, area=1, temperature=0.1, M=50, seed=1234): 639 if edges is None: 640 edges = StraightEdges() 641 642 self._edges = edges 643 self._area = area 644 self._temperature = temperature 645 self._M = M 646 self._seed = seed 647 648 def graph(self, vcoordinates, edges): 649 generator = numpy.random.RandomState(seed=self._seed) 650 # Setup parameters 651 k = numpy.sqrt(self._area / len(vcoordinates)) 652 653 # Initialize coordinates 654 mask = numpy.ma.getmaskarray(vcoordinates) 655 vcoordinates = numpy.ma.where(mask, generator.uniform(-1, 1, size=vcoordinates.shape), vcoordinates) 656 657 # Repeatedly apply attract / repel forces to the vertices 658 vertices = numpy.column_stack(numpy.triu_indices(n=len(vcoordinates), k=1)) 659 for temperature in numpy.linspace(self._temperature, 0, self._M, endpoint=False): 660 offsets = numpy.zeros_like(vcoordinates) 661 662 # Repel 663 a = vcoordinates[vertices.T[0]] 664 b = vcoordinates[vertices.T[1]] 665 delta = a - b 666 distance = numpy.linalg.norm(delta, axis=1)[:, None] 667 delta /= distance 668 force = numpy.square(k) / distance 669 delta *= force 670 _add_at(offsets, vertices.T[0], +delta) 671 _add_at(offsets, vertices.T[1], -delta) 672 673 # Attract 674 a = vcoordinates[edges.T[0]] 675 b = vcoordinates[edges.T[1]] 676 delta = b - a 677 distance = numpy.linalg.norm(delta, axis=1)[:, None] 678 delta /= distance 679 force = numpy.square(distance) / k 680 delta *= force 681 _add_at(offsets, edges.T[0], +delta) 682 _add_at(offsets, edges.T[1], -delta) 683 684 # Limit offsets to the temperature 685 distance = numpy.linalg.norm(offsets, axis=1) 686 offsets /= distance[:, None] 687 offsets *= numpy.minimum(temperature, distance)[:, None] 688 689 # Sum offsets 690 vcoordinates = numpy.ma.where(mask, vcoordinates + offsets, vcoordinates) 691 692 eshapes, ecoordinates = self._edges.edges(vcoordinates, edges) 693 return vcoordinates, eshapes, ecoordinates 694 695 696class Buchheim(GraphLayout): 697 """Compute a tree layout using the 2002 algorithm of Buchheim, Junger, and Leipert. 698 699 Note: this layout currently ignores preexisting vertex coordinates. 700 """ 701 def __init__(self, edges=None, basis=None): 702 if edges is None: 703 edges = StraightEdges() 704 705 if basis is None: 706 basis = [[1, 0], [0, -1]] 707 708 self._edges = edges 709 self._basis = numpy.array(basis) 710 711 def graph(self, vcoordinates, edges): 712 # Convert the graph to an adjacency list 713 children = _adjacency_list(len(vcoordinates), edges) 714 # Ensure we actually have a tree 715 root, depth = _require_tree(children) 716 717 # Get rid of the mask, it complicates things. 718 vcoordinates = numpy.array(vcoordinates) 719 720 # Convert our flat adjacency list into a hierarchy, to make the implementation easier. 721 class Vertex(object): 722 def __init__(self, vertex, parent=None, number=0, depth=0): 723 self.vertex = vertex 724 self.parent = parent 725 self.number = number 726 self.depth = depth 727 self.children = [Vertex(child, self, number, depth+1) for number, child in enumerate(children[vertex])] 728 self.mod = 0 729 self.thread = None 730 self.ancestor = self 731 self.prelim = 0 732 self.change = 0 733 self.shift = 0 734 735 # We follow Appendix A of the original paper as closely as possible here. 736 distance = 1 737 def FirstWalk(v): 738 if not v.children: # v is a leaf 739 v.prelim = 0 740 if v.number: # v has a left sibling 741 v.prelim = v.parent.children[v.number-1].prelim + distance 742 else: # v is not a leaf 743 defaultAncestor = v.children[0] # leftmost child of v 744 for w in v.children: 745 FirstWalk(w) 746 defaultAncestor = Apportion(w, defaultAncestor) 747 ExecuteShifts(v) 748 midpoint = 0.5 * (v.children[0].prelim + v.children[-1].prelim) 749 if v.number: # v has a left sibling 750 v.prelim = v.parent.children[v.number-1].prelim + distance 751 v.mod = v.prelim - midpoint 752 else: 753 v.prelim = midpoint 754 755 def Apportion(v, defaultAncestor): 756 if v.number: # v has a left sibling 757 vip = vop = v 758 vim = v.parent.children[v.number-1] 759 vom = vip.parent.children[0] 760 sip = vip.mod 761 sop = vop.mod 762 sim = vim.mod 763 som = vom.mod 764 while NextRight(vim) and NextLeft(vip): 765 vim = NextRight(vim) 766 vip = NextLeft(vip) 767 vom = NextLeft(vom) 768 vop = NextRight(vop) 769 vop.ancestor = v 770 shift = (vim.prelim + sim) - (vip.prelim + sip) + distance 771 if shift > 0: 772 MoveSubtree(Ancestor(vim, v, defaultAncestor), v, shift) 773 sip += shift 774 sop += shift 775 sim += vim.mod 776 sip += vip.mod 777 som += vom.mod 778 sop += vop.mod 779 if NextRight(vim) and not NextRight(vop): 780 vop.thread = NextRight(vim) 781 vop.mod += sim - sop 782 if NextLeft(vip) and not NextLeft(vom): 783 vom.thread = NextLeft(vip) 784 vom.mod += sip - som 785 defaultAncestor = v 786 return defaultAncestor 787 788 def NextLeft(v): 789 if v.children: 790 return v.children[0] 791 else: 792 return v.thread 793 794 def NextRight(v): 795 if v.children: 796 return v.children[-1] 797 else: 798 return v.thread 799 800 def MoveSubtree(wm, wp, shift): 801 subtrees = wp.number - wm.number 802 wp.change -= shift / subtrees 803 wp.shift += shift 804 wm.change += shift / subtrees 805 wp.prelim += shift 806 wp.mod += shift 807 808 def ExecuteShifts(v): 809 shift = 0 810 change = 0 811 for w in v.children: 812 w.prelim += shift 813 w.mod += shift 814 change += w.change 815 shift += w.shift + change 816 817 def Ancestor(vim, v, defaultAncestor): 818 if vim.ancestor in v.parent.children: 819 return vim.ancestor 820 else: 821 return defaultAncestor 822 823 def SecondWalk(v, m): 824 vcoordinates[v.vertex][0] = v.prelim + m 825 vcoordinates[v.vertex][1] = v.depth 826 for w in v.children: 827 SecondWalk(w, m + v.mod) 828 829 r = Vertex(root) 830 FirstWalk(r) 831 SecondWalk(r, -r.prelim) 832 833 vcoordinates = numpy.dot(vcoordinates, self._basis) 834 835 eshapes, ecoordinates = self._edges.edges(vcoordinates, edges) 836 return vcoordinates, eshapes, ecoordinates 837