weavingspace.topology

   1#!/usr/bin/env python
   2# coding: utf-8
   3
   4"""Together the `Topology`, `weavingspace.topology_elements.Tile`, 
   5`weavingspace.topology_elements.Vertex`, `weavingspace.topology_elements.Edge`,
   6and `weavingspace.symmetry.Symmetries` classes enable extraction of the 
   7topological structure of periodic `weavingspace.tileable.Tileable` objects so 
   8that modification of equivalent tiles can be carried out while retaining 
   9tileability. It is important to note that these are not fully generalised 
  10classes and methods, that is, the Topology object that is supported is not a 
  11permanent 'backing' data structure for our Tileable objects. While it might 
  12become that in time, as at Feb 2024 it is not such a data structure. Instead 
  13usage is
  14
  15    tile = TileUnit(...)
  16    topology = Topology(tile)
  17    topology = topology.transform_*(...)
  18    new_tile = topology.tile_unit
  19
  20Topology plot function is necessary for a user to be able to see what they are
  21doing, because how edges and vertices in a tiling are labelled under tile
  22equivalences is an essential step in the process.
  23
  24Note also that these classes do not accurately represent the distinctions made
  25in the mathematical literature between tiling vertices and tile corners, or 
  26between tiling edges and tile sides.
  27"""
  28from collections import defaultdict
  29import inspect
  30import itertools
  31from typing import Callable, Union
  32from typing import Any
  33from typing import Iterable
  34import pickle
  35import string
  36
  37import numpy as np
  38import networkx as nx
  39import geopandas as gpd
  40import shapely.geometry as geom
  41import shapely.affinity as affine
  42import matplotlib.pyplot as pyplot
  43
  44from weavingspace import Tileable
  45from weavingspace import Symmetries
  46from weavingspace import Shape_Matcher
  47from weavingspace.topology_elements import *
  48import weavingspace.tiling_utils as tiling_utils
  49
  50ALPHABET = string.ascii_letters.upper()
  51alphabet = string.ascii_letters.lower()
  52
  53
  54class Topology:
  55  """Class to represent topology of a Tileable object.
  56  
  57  NOTE: It is important that get_local_patch return the tileable elements and 
  58  the translated copies in consistent sequence, i.e. if there are (say) four 
  59  tiles in the unit, the local patch should be 1 2 3 4 1 2 3 4 1 2 3 4 ... and
  60  so on. This is because self.tiles[i % n_tiles] is frequently used to reference
  61  the base unit Tile which corresponds to self.tiles[i].
  62  """
  63  tileable: Tileable
  64  """the Tileable on which the topology will be based."""
  65  tiles: list[Tile]
  66  """list of the Tiles in the topology. We use polygons returned by the
  67  tileable.get_local_patch method for these. That is the base tiles and 8 adjacent 
  68  copies (for a rectangular tiling), or 6 adjacent copies (for a hexagonal tiling)."""
  69  points: dict[int, Vertex]
  70  """dictionary of all points (vertices and corners) in the tiling, keyed by Vertex ID."""
  71  edges: dict[int, Edge]
  72  """dictionary of the tiling edges, keyed by Edge ID."""
  73  unique_tile_shapes: list[geom.Polygon]
  74  """a 'reference' tile shape one per tile shape (up to vertices, so two tiles
  75  might be the same shape, but one might have extra vertices induced by the
  76  tiling and hence is a different shape under this definition)."""
  77  dual_tiles: list[geom.Polygon]
  78  """list of geom.Polygons from which a dual tiling might be constructed."""
  79  n_tiles: int = 0
  80  """number of tiles in the base Tileable (retained for convenience)."""
  81  shape_groups: list[list[int]]
  82  """list of lists of tile IDs distinguished by shape and optionally tile_id"""
  83  tile_matching_transforms: list[tuple[float]] 
  84  """shapely transform tuples that map tiles onto other tiles"""
  85  tile_transitivity_classes: list[tuple[int]]
  86  """list of lists of tile IDs in each transitivity class"""
  87  vertex_transitivity_classes: list[list[int]]
  88  """list of lists of vertex IDs in each transitivity class"""
  89  edge_transitivity_classes: list[list[tuple[int]]]
  90  """list of lists of edge IDs in each transitivity class"""
  91
  92  def __init__(self, unit: Tileable, ignore_tile_ids:bool = True):
  93    """Class constructor.
  94
  95    Args:
  96      unit (Tileable): the Tileable whose topology is required.
  97    """
  98    # Note that the order of these setup steps is fragile sometimes
  99    # not obviously so.
 100    self.tileable = unit # keep this for reference
 101    self.n_tiles = self.tileable.tiles.shape[0]
 102    self._initialise_points_into_tiles()
 103    self._setup_vertex_tile_relations()
 104    self._setup_edges()
 105    self._copy_base_tiles_to_patch()
 106    self._assign_vertex_and_edge_base_IDs()
 107    self._identify_distinct_tile_shapes(ignore_tile_ids)
 108    self._find_tile_transitivity_classes()
 109    self._find_vertex_transitivity_classes()
 110    self._find_edge_transitivity_classes()
 111    self._label_tiles() # now no longer required...
 112    self.generate_dual()
 113
 114  def __str__(self) -> str:
 115    """Returns string representation of this Topology.
 116    
 117    Returns:
 118      str: a message that recommends examining the tiles, points and edges 
 119        attributes.
 120
 121    """
 122    return f"""Topology of Tileable with {self.n_tiles} tiles.\n
 123    Examine .tiles, .points and .edges for more details."""
 124
 125  def __repr__(self) -> str:
 126    return str(self)
 127
 128  def _initialise_points_into_tiles(self, debug:bool = False) -> None:
 129    """Sets up dictionary of unique point locations in the tiling and assigns
 130    them to Tiles.
 131    
 132    Args:
 133      debug (bool): if True prints useful debugging information.
 134    """
 135    shapes = self.tileable.get_local_patch(r = 1, include_0 = True).geometry
 136    shapes = [tiling_utils.get_clean_polygon(s) for s in shapes]
 137    self.n_1_order_patch = len(shapes)
 138    labels = list(self.tileable.tiles.tile_id) * (len(shapes) // self.n_tiles)
 139    self.tiles = []
 140    self.points = {}
 141    for (i, shape), label in zip(enumerate(shapes), labels):
 142      tile = Tile(i)
 143      tile.label = label
 144      tile.base_ID = tile.ID % self.n_tiles
 145      self.tiles.append(tile)
 146      tile.corners = []
 147      corners = tiling_utils.get_corners(shape, repeat_first = False)
 148      for c in corners:
 149        prev_vertex = None
 150        for p in reversed(self.points.values()):
 151          if c.distance(p.point) <= 2 * tiling_utils.RESOLUTION:
 152            # found an already existing vertex, so add to the tile and break
 153            prev_vertex = p
 154            tile.corners.append(prev_vertex)
 155            break
 156        if prev_vertex is None:
 157          # new vertex, add it to the topology dictionary and to the tile
 158          v = self.add_vertex(c)
 159          tile.corners.append(v)
 160          if debug: print(f"Added new Vertex {v} to Tile {i}")
 161
 162  def _setup_vertex_tile_relations(self, debug:bool = False):
 163    """Determines relations between vertices and tiles. In particular vertices
 164    along tile edges that are not yet included in their list of vertices are
 165    added. Meanwhile vertex lists of incident tiles and neighbours are set up.
 166    
 167    Args:
 168      debug (bool): if True prints debugging information.
 169    """
 170    # we do this for all tiles in the radius-1 local patch
 171    for tile in self.tiles:
 172      if debug: print(f"Checking for vertices incident on Tile {tile.ID}")
 173      corners = []
 174      # performance is much improved using the vertex IDs
 175      initial_corner_IDs = tile.get_corner_IDs()
 176      # we need the current shape (not yet set) to check for incident vertices 
 177      shape = geom.Polygon([c.point for c in tile.corners])
 178      # get points incident on tile boundary, not already in the tile corners
 179      new_points = [v for v in self.points.values()
 180                    if v.ID not in initial_corner_IDs and
 181                    v.point.distance(shape) <= 2 * tiling_utils.RESOLUTION]
 182      # iterate over sides of the tile to see which the vertex is incident on
 183      for c1, c2 in zip(tile.corners, tile.corners[1:] + tile.corners[:1]):
 184        all_points = [c1, c2]
 185        if len(new_points) > 0:
 186          if debug: print(f"{[v.ID for v in new_points]} incident on tile")
 187          ls = geom.LineString([c1.point, c2.point])
 188          to_insert = [v for v in new_points 
 189                       if v.point.distance(ls) <= 2 * tiling_utils.RESOLUTION]
 190          if len(to_insert) > 0:
 191            # sort by distance along the side
 192            d_along = sorted([(ls.line_locate_point(v.point), v)
 193                              for v in to_insert], key = lambda x: x[0])
 194            to_insert = [v for d, v in d_along]
 195            all_points = all_points[:1] + to_insert + all_points[1:]
 196        corners.extend(all_points[:-1])
 197        for x1, x2 in zip(all_points[:-1], all_points[1:]):
 198          # x2 will add the tile and neigbour when we get to the next side
 199          # every vertex gets a turn!
 200          x1.add_tile(tile)
 201          x1.add_neighbour(x2.ID)
 202      tile.corners = corners
 203      tile.set_shape_from_corners()
 204
 205  def _setup_edges(self, debug:bool = False):
 206    """Sets up the tiling edges. 
 207    
 208    First vertices in the base tiles are classified as tiling vertices or not -- 
 209    only these can be classified reliably (e.g vertices on the perimeter are 
 210    tricky).. Up to here all vertices have been considered tiling vertices.
 211    
 212    Second edges are created by traversing tile corner lists. Edges are stored
 213    once only by checking for edges in the reverse direction already in the 
 214    edges dictionary. Edge right and left tiles are initialised.
 215    
 216    Third tile edge direction lists are initialised.
 217    """
 218    # classify vertices in the base tiles
 219    for tile in self.tiles[:self.n_tiles]:
 220      for v in tile.corners:
 221        v.is_tiling_vertex = len(v.neighbours) > 2
 222    if debug: print(f"Classified base tile vertices")
 223    self.edges = {}
 224    for tile in self.tiles:
 225      if debug: print(f"Adding edges from Tile {tile.ID}")
 226      tile.edges = []
 227      vertices = [v for v in tile.corners if v.is_tiling_vertex]
 228      # note that through here finding ints in lists is much faster than
 229      # finding Vertex objects, hence we use lists of IDs not Vertex objects
 230      if len(vertices) > 1:
 231        for v1, v2 in zip(vertices, vertices[1:] + vertices[:1]):
 232          corner_IDs = tile.get_corner_IDs()
 233          idx1 = corner_IDs.index(v1.ID)
 234          idx2 = corner_IDs.index(v2.ID)
 235          if idx1 < idx2:
 236            corners = [c for c in corner_IDs[idx1:(idx2 + 1)]]
 237          else:
 238            corners = [c for c in corner_IDs[idx1:] + corner_IDs[:(idx2 + 1)]]
 239          ID = (corners[0], corners[-1])
 240          if not ID in self.edges:
 241            # check that reverse direction edge is not present first
 242            r_ID = ID[::-1]
 243            if r_ID in self.edges:
 244              # if it is, then add it and set left_tile
 245              if debug: print(f"reverse edge {r_ID} found")
 246              e = self.edges[r_ID]
 247              e.left_tile = tile
 248              tile.edges.append(e)
 249            else:
 250              # we've found a new edge so make and add it
 251              if debug: print(f"adding new edge {corners}")
 252              e = self.add_edge([self.points[c] for c in corners])
 253              e.right_tile = tile
 254              tile.edges.append(e)
 255      # initialise the edge direction information in the tile
 256      tile.set_edge_directions()
 257      
 258  def _assign_vertex_and_edge_base_IDs(self):
 259    """Assigns the base_ID attributes of vertices and edges. These allow us
 260    to determine corrspondences between vertices and edges in the 'base' tiles
 261    in the Topology tileable, and those we have added at radius 1 for labelling
 262    and visualisation.
 263    """
 264    self._assign_vertex_base_IDs(True, False)
 265    self._assign_vertex_base_IDs(False, True)
 266    self._assign_edge_base_IDs()
 267
 268  def _assign_vertex_base_IDs(self, first_pass:bool = True,
 269                              last_pass:bool = False):
 270    """Assigns base_ID attribute of vertices. Two passes are required.
 271    
 272    TODO: consider if there is a networkx.connected_communities option
 273    here, since this is in essence about finding connected groups of tiles,
 274    vertices and edges. At very least, the need to run the code twice suggests
 275    that a more elegant approach is out there...
 276
 277    Args:
 278      first_pass (bool, optional): it True the core vertices on the tiles in 
 279        the Topology tileable are assigned; if False this step is skipped and
 280        IDs are propagated to neighbours. Defaults to True.
 281      last_pass (bool, optional): if True, it is the last pass and IDs are 
 282        finally assigned based on the lowest ID assigned in previous rounds. 
 283        Defaults to False.
 284    """
 285    if first_pass:
 286      for i, v in enumerate(self.vertices_in_tiles(self.tiles[:self.n_tiles])):
 287        v.base_ID = set([i])
 288    for tile in self.tiles:
 289      for v, vb in zip(tile.corners, self.tiles[tile.base_ID].corners):
 290        if v.base_ID == v.ID: # 1_000_000:
 291          v.base_ID = set([x for x in vb.base_ID])
 292        else:
 293          v.base_ID = v.base_ID.union(vb.base_ID)
 294    if last_pass:
 295      for tile in self.tiles:
 296        for v in tile.corners:
 297          if isinstance(v.base_ID, set):
 298            v.base_ID = min(v.base_ID)
 299
 300  def _assign_edge_base_IDs(self):
 301    """Assigns base_ID attribute of edges, based on the base_ID attributes of
 302    vertices. This might be oversimplified...
 303    """
 304    for tile in self.tiles:
 305      for e, eb in zip(tile.edges, self.tiles[tile.base_ID].edges):
 306        v0 = self.points[eb.ID[0]].base_ID
 307        v1 = self.points[eb.ID[1]].base_ID
 308        e.base_ID = (v0, v1)
 309
 310  def _copy_base_tiles_to_patch(self):
 311    """Copies attributes of base tiles to their corresponding tiles in 
 312    the radius-1 patch. This requires:
 313     
 314    First, inserting any additional corners in the base tiles not found in the 
 315    radius-1 tiles. 
 316    
 317    Second, any vertices in the base tiles that are NOT tiling vertices are 
 318    applied to radius-1 tiles leading to merging of some edges.
 319    """
 320    # the number of tiles in the base + radius-1
 321    n_r1 = len(self.tiles)
 322    # first add any missing vertices to the non-base tiles
 323    # we add all the missing vertices before doing any merges
 324    for base in self.tiles[:self.n_tiles]:
 325      for other in self.tiles[base.ID:n_r1:self.n_tiles]:
 326        self._match_reference_tile_vertices(base, other)
 327    # then merge any edges that meet at a corner
 328    for base in self.tiles[:self.n_tiles]:
 329      for other in self.tiles[base.ID:n_r1:self.n_tiles]:
 330        self._match_reference_tile_corners(base, other)
 331
 332  def _match_reference_tile_vertices(self, tile1:Tile, tile2:Tile):
 333    """Adds vertices to tile2 so that it matches tile1 adjusting edges as
 334    required. This assumes the tiles are the same shape, but that tile2 may
 335    be missing some tiling vertices along some edges.
 336
 337    Args:
 338      tile1 (Tile): reference tile.
 339      tile2 (Tile): tile to change.
 340    """
 341    to_add = len(tile1.corners) - len(tile2.corners)
 342    while to_add > 0:
 343      # find the reference x-y offset 
 344      dxy = (tile2.centre.x - tile1.centre.x, tile2.centre.y - tile1.centre.y)
 345      for i, t1c in enumerate([c.point for c in tile1.corners]):
 346        t2c = tile2.corners[i % len(tile2.corners)].point
 347        if abs((t2c.x - t1c.x) - dxy[0]) > 10 * tiling_utils.RESOLUTION or \
 348           abs((t2c.y - t1c.y) - dxy[1]) > 10 * tiling_utils.RESOLUTION:
 349          # add vertex to t2 by copying the t1 vertex appropriately offset
 350          # note that this might alter the length of t2.corners
 351          v = self.add_vertex(geom.Point(t1c.x + dxy[0], t1c.y + dxy[1]))
 352          v.is_tiling_vertex = True
 353          old_edge, new_edges = tile2.insert_vertex_at(v, i)
 354          del self.edges[old_edge]
 355          for e in new_edges:
 356            self.edges[e.ID] = e
 357          to_add = to_add - 1
 358
 359  def _match_reference_tile_corners(self, tile1:Tile, tile2:Tile):
 360    """Finds vertices that are corners in tile2 but vertices in tile1 and 
 361    updates tile2 to match - merging edges as required.
 362
 363    Args:
 364        tile1 (Tile): reference tile.
 365        tile2 (Tile): tile to make match.
 366    """
 367    vs_to_change = [vj for vi, vj in zip(tile1.corners, tile2.corners)
 368                    if not vi.is_tiling_vertex and vj.is_tiling_vertex]
 369    if len(vs_to_change) > 0:
 370      for v in vs_to_change:
 371        v.is_tiling_vertex = False
 372        # it's a corner not an edge so will have no more than 2 v.tiles
 373        old_edges, new_edge = v.tiles[0].merge_edges_at_vertex(v)
 374        for e in old_edges:
 375          del self.edges[e]
 376        self.edges[new_edge.ID] = new_edge
 377
 378  def _identify_distinct_tile_shapes(self, ignore_tile_id_labels:bool = True):
 379    """Determines unique tiles based on their symmetries and shapes. At the
 380    same time assembles a list of the affine transforms under which matches
 381    occurs since these are potential symmetries of the tiling.
 382    
 383    TODO: reimplement consideration of tile_id
 384    
 385    Args:
 386      ignore_tile_id_labels (bool): if True only the shape of tiles matters; if 
 387        False the tile_id label is also considered. Defaults to True.
 388    """
 389    matches = {}
 390    offsets = {}
 391    for tile in self.tiles[:self.n_tiles]:
 392      matches[tile.base_ID] = [tile.base_ID]
 393      matched = False
 394      s = Symmetries(tile.shape)
 395      for other in self.tiles[:self.n_tiles]:
 396        if other.ID > tile.ID:
 397          offset = s.get_corner_offset(other.shape)
 398          if not offset is None:
 399            offsets[tile.base_ID] = offset
 400            matches[tile.base_ID].append(other.base_ID)
 401            matched = True
 402      if not matched:
 403        offsets[tile.base_ID] = 0
 404    base_groups = list(nx.connected_components(nx.from_dict_of_lists(matches)))
 405    self.shape_groups = []
 406    for i, group in enumerate(base_groups):
 407      full_group = []
 408      for tile in self.tiles:
 409        if tile.base_ID in group:
 410          tile.shape_group = i
 411          tile.offset_corners(offsets[tile.base_ID])
 412          full_group.append(tile.ID)
 413      self.shape_groups.append(full_group)
 414
 415  def _find_tile_transitivity_classes(self):
 416    """Finds tiles equivalent under symmetries, at the same time updating the 
 417    tile_matching_transforms attribute to contain only those transforms that 
 418    pass this test. 
 419    """
 420    self.tile_matching_transforms = self.get_potential_symmetries()
 421    base_tiles = [t for t in self.tiles[:self.n_tiles]]
 422    # it is quicker (should be!) to only do within shape group tests
 423    # often there is only one when it will make no difference
 424    by_group_equivalent_tiles = []
 425    # maintain a set of transforms still potentially tiling symmetries
 426    poss_transforms = set(self.tile_matching_transforms.keys())
 427    # and a dictionary of booleans tracking which transforms are still valid
 428    eq_under_transform = {tr: True for tr in poss_transforms}
 429    for g, group in enumerate(self.shape_groups):
 430      by_group_equivalent_tiles.append(set())
 431      source_tiles = [tile for tile in base_tiles if tile.shape_group == g]
 432      target_tiles = [tile for tile in self.tiles if tile.shape_group == g]
 433      for tr in poss_transforms:
 434        transform = self.tile_matching_transforms[tr].transform
 435        matched_tiles = {}
 436        eq_under_transform[tr] = True
 437        for source_tile in source_tiles:
 438          matched_tile_id = self._match_geoms_under_transform(
 439            source_tile, target_tiles, transform)
 440          if matched_tile_id == -1:
 441            eq_under_transform[tr] = False
 442            break
 443          else:
 444            matched_tiles[source_tile.ID] = matched_tile_id
 445        if eq_under_transform[tr]:
 446          for k, v in matched_tiles.items():
 447            # print(f"src {k} tgt {v} {tr=} {transform}")
 448            # here we record the transform, in case it is later invalidated
 449            by_group_equivalent_tiles[g].add((tr, k, v))
 450      # remove valid transforms that didn't make it through this group
 451      poss_transforms = set([t for t, x in eq_under_transform.items() if x])
 452    # compile equivalences from all groups made under still valid transforms
 453    # a dict of sets so singletons aren't lost in finding connected components
 454    equivalents = {i: set() for i in range(self.n_tiles)}
 455    for group_equivalents in by_group_equivalent_tiles:
 456      for (tr, tile_i, tile_j) in group_equivalents:
 457        if tr in poss_transforms:
 458          equivalents[tile_i].add(tile_j)
 459    self.tile_matching_transforms = {
 460      k: v for k, v in self.tile_matching_transforms.items() if k in 
 461      poss_transforms}
 462    self.tile_transitivity_classes = []
 463    equivalents = nx.connected_components(nx.from_dict_of_lists(equivalents))
 464    for c, base_IDs in enumerate(equivalents):
 465      transitivity_class = []
 466      for tile in self.tiles:
 467        if tile.base_ID in base_IDs:
 468          transitivity_class.append(tile.ID)
 469          tile.transitivity_class = c
 470      self.tile_transitivity_classes.append(transitivity_class)
 471
 472  def get_potential_symmetries(self) -> dict[int, tuple[float]]:
 473    """Assembles potential symmetries of the tiling from symmetries of the 
 474    tileable.prototile and of the tileable.tiles. Removes any duplicates that
 475    result. Result is assigned to the tile_matching_transforms attribute.
 476    
 477    TODO: consider retaining the Symmetry objects as these carry additional 
 478    information that might facilitate labelling under a limited number of the
 479    symmetries not all of them.
 480
 481    Returns:
 482      dict[int, tuple[float]]: dictionary of the symmetries (transforms 
 483        actually) in shapely affine transform 6-tuple format.
 484    """
 485    self.tile_matching_transforms = {}
 486    n_symmetries = 0
 487    pt = self.tileable.prototile.geometry[0]
 488    for tr in Shape_Matcher(pt).get_polygon_matches(pt):
 489      if not tr.transform_type in ["identity", "translation"]:
 490        self.tile_matching_transforms[n_symmetries] = tr
 491        n_symmetries = n_symmetries + 1
 492    for tile in self.tiles[:self.n_tiles]:
 493      for tr in Shape_Matcher(tile.shape).get_polygon_matches(tile.shape):
 494        if not tr.transform_type in ["identity", "translation"]:
 495          self.tile_matching_transforms[n_symmetries] = tr
 496          n_symmetries = n_symmetries + 1
 497    for tile in self.tiles[:self.n_tiles]:
 498      sm = Shape_Matcher(tile.shape)
 499      transforms = [sm.get_polygon_matches(self.tiles[i].shape) 
 500        for i in self.shape_groups[tile.shape_group] if i < self.n_tiles]
 501      for tr in itertools.chain(*transforms):
 502        if not tr.transform_type in ["identity", "translation"]:
 503          self.tile_matching_transforms[n_symmetries] = tr
 504          n_symmetries = n_symmetries + 1
 505    self.tile_matching_transforms = self._remove_duplicate_symmetries(
 506      self.tile_matching_transforms)
 507    return self.tile_matching_transforms
 508
 509  def _remove_duplicate_symmetries(self, transforms:dict[int,Any]):
 510    """Filters the supplied list of shapely affine transforms so that no 
 511    duplicates are retained.
 512    """
 513    uniques = {}
 514    for k, v in transforms.items():
 515      already_exists = False
 516      for i, u in uniques.items():
 517        already_exists = np.allclose(
 518          v.transform, u.transform, atol = 1e-4, rtol = 1e-4)
 519        if already_exists:
 520          break
 521      if not already_exists:
 522        uniques[k] = v
 523    return uniques
 524
 525  def _find_vertex_transitivity_classes(self):
 526    """Finds vertex transitivity classes by checking which vertices align
 527    with which others under transforms in the tile_matching_transforms 
 528    attribute. The process need only determine the classes for vertices
 529    in the core tileable.tiles, then assign those to all vertices by 
 530    matched base_ID.
 531    
 532    TODO: carefully consider here whether labelling only the tiling vertices
 533    is the right thing to do or not...
 534    """
 535    equivalent_vertices = defaultdict(set)
 536    base_vertices = [v for v in 
 537                     self.vertices_in_tiles(self.tiles[:self.n_tiles])
 538                     if v.is_tiling_vertex]
 539    for transform in self.tile_matching_transforms.values():
 540      for v in base_vertices:
 541        match_ID = self._match_geoms_under_transform(
 542          v, base_vertices, transform.transform)
 543        if match_ID != -1:
 544          equivalent_vertices[v.ID].add(match_ID)
 545    equivalent_vertices = self._get_exclusive_supersets(
 546      [tuple(sorted(s)) for s in equivalent_vertices.values()])
 547    self.vertex_transitivity_classes = defaultdict(list)
 548    for c, vclass in enumerate(equivalent_vertices):
 549      for v in self.points.values():
 550        if v.base_ID in vclass:
 551          v.transitivity_class = c
 552          self.vertex_transitivity_classes[c].append(v.ID)
 553    self.vertex_transitivity_classes = list(
 554      self.vertex_transitivity_classes.values())
 555    # label vertices based on their transitivity class
 556    for v in self.points.values():
 557      if v.is_tiling_vertex:
 558        v.label = ALPHABET[v.transitivity_class]
 559
 560  def _find_edge_transitivity_classes(self):
 561    """Finds edge transitivity classes by checking which edges align
 562    with which others under transforms in the tile_matching_transforms 
 563    attribute. The process need only determine the classes for edges
 564    in the core tileable.tiles, then assign those to all edges by 
 565    matched base_ID.
 566    
 567    TODO: Note that this code is identical to the vertex transitivity code
 568    so it might make sense to merge.
 569    """
 570    equivalent_edges = defaultdict(set)
 571    base_edges = self.edges_in_tiles(self.tiles[:self.n_tiles])
 572    for transform in self.tile_matching_transforms.values():
 573      for e in base_edges:
 574        match_id = self._match_geoms_under_transform(
 575          e, base_edges, transform.transform)
 576        if match_id != -1:
 577          equivalent_edges[e.ID].add(match_id)
 578    equivalent_edges = self._get_exclusive_supersets(
 579      [tuple(sorted(s)) for s in equivalent_edges.values()])
 580    self.edge_transitivity_classes = defaultdict(list)
 581    for c, eclass in enumerate(equivalent_edges):
 582      for e in self.edges.values():
 583        if e.base_ID in eclass:
 584          e.transitivity_class = c
 585          self.edge_transitivity_classes[c].append(e.ID)
 586    self.edge_transitivity_classes = list(
 587      self.edge_transitivity_classes.values())
 588    # label edges based on their transitivity class
 589    for e in self.edges.values():
 590      e.label = alphabet[e.transitivity_class]
 591    
 592  def _match_geoms_under_transform(
 593      self, geom1:Union[Tile,Vertex,Edge], 
 594      geoms2:list[Union[Tile,Vertex,Edge]], transform:tuple[float]) -> int:
 595    """Determines if there is a geometry in the supplied patch onto which the
 596    supplied geometry is mapped by the supplied symmetry.
 597
 598    Args:
 599      tile (Union[Tile,Vertex,Edge]): element whose geometry we want to match.
 600      patch (list[Union[Tile,Vertex,Edge]]): set of elements among which a 
 601        match is sought.
 602      transform (tuple[float]): shapely affine transform 6-tuple to apply.
 603
 604    Returns:
 605      Union[int,tuple[int]]: ID of the element in patch that matches the geom1
 606        element under the transform if one exists, otherwise returns -1.
 607    """
 608    match_id = -1
 609    for geom2 in geoms2:
 610      if isinstance(geom1, Tile):
 611        # an area of intersection based test
 612        match = tiling_utils.geometry_matches(
 613          affine.affine_transform(geom1.shape, transform), geom2.shape)
 614      elif isinstance(geom1, Vertex):
 615        # distance test
 616        match = affine.affine_transform(geom1.point, transform).distance(
 617          geom2.point) <= 10 * tiling_utils.RESOLUTION
 618      else: # must be an Edge
 619        # since edges _should not_ intersect this test should work in
 620        # lieu of a more complete point by point comparison
 621        c1 = geom1.get_geometry().centroid
 622        c2 = geom2.get_geometry().centroid
 623        match = affine.affine_transform(c1, transform) \
 624          .distance(c2) <= 10 *tiling_utils.RESOLUTION
 625      if match:
 626        return geom2.base_ID
 627    return match_id
 628
 629  def _get_exclusive_supersets(
 630      self, sets:list[Iterable[Any]]) -> list[Iterable[Any]]:
 631    """For the supplied list of lists, which may share elements, i.e. they are
 632    non-exclusives sets, returns a list of lists which are exclusive: each 
 633    element will only appear in one of the lists in the returned list. Uses
 634    networkx connected components function to achieve this based on a graph 
 635    where an intersection between two sets is an edge.
 636
 637    Args:
 638        sets (list[Iterable[Any]]): list of lists of possibly overlapping sets.
 639
 640    Returns:
 641      list[Iterable[Any]]: list of lists that include all the original 
 642        elements without overlaps.
 643    """
 644    overlaps = []
 645    for i, si in enumerate(sets):
 646      s1 = set(si)
 647      for j, sj in enumerate(sets):
 648        s2 = set(sj)
 649        if len(s1 & s2) > 0:
 650          overlaps.append((i, j))
 651    G = nx.from_edgelist(overlaps)
 652    result = []
 653    for component in nx.connected_components(G):
 654      s = set()
 655      for i in component:
 656        s = s.union(sets[i])
 657      result.append(tuple(s))
 658    return result
 659
 660  def vertices_in_tiles(self, tiles:list[Tile]) -> list[Vertex]:
 661    """Gets the vertices from self.points that are incident on the tiles in the 
 662    supplied list.
 663
 664    Args:
 665        tiles (list[Tile]): tiles whose vertices are required.
 666
 667    Returns:
 668        list[Vertex]: the required vertices.
 669    """
 670    vs = set()
 671    for tile in tiles:
 672      vs = vs.union(tile.get_corner_IDs())
 673    return [self.points[v] for v in vs]
 674
 675  def edges_in_tiles(self, tiles:list[Tile]) -> list[Edge]:
 676    """Gets the edges from self.edges that are part of the boundary of tiles in 
 677    the supplied list.
 678
 679    Args:
 680        tiles (list[Tile]): tiles whose edges are required.
 681
 682    Returns:
 683        list[Edge]: the required edges.
 684    """
 685    es = set()
 686    for tile in tiles:
 687      es = es.union(tile.get_edge_IDs())
 688    return [self.edges[e] for e in es]
 689
 690  def generate_dual(self) -> list[geom.Polygon]:
 691    """Create the dual tiiing for the tiling of this Topology.
 692    
 693    TODO: make this a viable replacement for the existing dual tiling 
 694    generation.
 695
 696    Returns:
 697      list[geom.Polygon]: a list of polygon objects.
 698    """
 699    for v in self.points.values():
 700      v.clockwise_order_incident_tiles()
 701    self.dual_tiles = {}
 702    base_id_sets = defaultdict(list)
 703    for v in self.points.values():
 704      base_id_sets[v.base_ID].append(v.ID)
 705    minimal_set = [self.points[min(s)] for s in base_id_sets.values()]
 706    for v in minimal_set:
 707    # for v in self.points.values():
 708      if v.is_interior() and len(v.tiles) > 2:
 709        self.dual_tiles[v.ID] = \
 710          geom.Polygon([t.centre for t in v.tiles])
 711
 712  def add_vertex(self, pt:geom.Point) -> Vertex:
 713    """Adds a Vertex at the specified point location, returning it to the 
 714    caller. No attempt is made to ensure Vertex IDs are an unbroken sequemce: a 
 715    new ID is generated one greater than the existing highest ID. IDs will 
 716    usually be an unbroken sequence up to removals when geometry transformations
 717    are applied.
 718
 719    Args:
 720        pt (geom.Point): point location of the Vertex.
 721
 722    Returns:
 723        Vertex: the added Vertex object.
 724    """
 725    n = 0 if len(self.points) == 0 else max(self.points.keys()) + 1
 726    v = Vertex(pt, n)
 727    self.points[n] = v
 728    return v
 729
 730  def add_edge(self, vs:list[Vertex]) -> Edge:
 731    """Adds an Edge to the edges dictionary. Edges are self indexing by the IDs
 732     of their end Vertices. Returns the new Edge to the caller.
 733
 734    Args:
 735      vs (list[Vertex]): list of Vertices in the Edge to be created.
 736
 737    Returns:
 738        Edge: the added Edge.
 739    """
 740    e = Edge(vs)
 741    self.edges[e.ID] = e
 742    return e
 743
 744  def _get_tile_geoms(self) -> gpd.GeoDataFrame:
 745    return gpd.GeoDataFrame(
 746      data = {"transitivity_class": [t.transitivity_class for t in self.tiles]},
 747      geometry = gpd.GeoSeries([t.shape for t in self.tiles]))
 748
 749  def _get_tile_centre_geoms(self) -> gpd.GeoSeries:
 750    return gpd.GeoSeries([t.centre for t in self.tiles])
 751
 752  def _get_point_geoms(self) -> gpd.GeoSeries:
 753    return gpd.GeoSeries([v.point for v in self.points.values()])
 754
 755  def _get_vertex_geoms(self) -> gpd.GeoSeries:
 756    return gpd.GeoSeries([v.point for v in self.points.values()
 757                          if v.is_tiling_vertex])
 758
 759  def _get_edge_geoms(self, offset:float = 0.0) -> gpd.GeoSeries:
 760    return gpd.GeoSeries([e.get_topology().parallel_offset(offset)
 761                          for e in self.edges.values()])
 762
 763  def plot(self, show_original_tiles: bool = True,  
 764           show_tile_centres: bool = False,
 765           show_tile_vertex_labels: bool = False, 
 766           show_tile_edge_labels: bool = False, 
 767           show_vertex_ids: bool = False,
 768           show_vertex_labels: bool = True,
 769           show_edges: bool = True,
 770           offset_edges: bool = True,
 771           show_edge_labels:bool = False,
 772           show_dual_tiles: bool = False) -> pyplot.Axes:
 773    fig = pyplot.figure(figsize = (10, 10))
 774    ax = fig.add_axes(111)
 775    extent = gpd.GeoSeries([t.shape for t in self.tiles]).total_bounds
 776    dist = max([extent[2] - extent[0], extent[3] - extent[1]])
 777    if show_original_tiles:
 778      self._get_tile_geoms().plot(column = "transitivity_class", cmap = "Set2",
 779        ax = ax, ec = "#444444", alpha = 0.2, lw = 0.75)
 780    if show_tile_centres:
 781      for i, tile in enumerate(self.tiles):
 782        ax.annotate(i, xy = (tile.centre.x, tile.centre.y), color = "b", 
 783                    fontsize = 10, ha = "center", va = "center")
 784    if show_vertex_labels:
 785      for v in self.points.values():
 786        ax.annotate(v.ID if show_vertex_ids else v.label, 
 787                    xy = (v.point.x, v.point.y), color = "r", fontsize = 10,
 788                    ha = "center", va = "center",
 789                    bbox = dict(boxstyle="circle", lw=0, fc="#ffffff40"))
 790    if show_tile_vertex_labels:
 791      for tile in self.tiles:
 792        for posn, label in zip(tile.get_vertex_label_positions(), 
 793                               tile.vertex_labels):
 794          ax.annotate(label, xy = (posn.x, posn.y), fontsize = 8, 
 795                      ha = "center", va = "center", color = "b")
 796    if show_tile_edge_labels:
 797      for tile in self.tiles:
 798        for posn, label in zip(tile.get_edge_label_positions(),
 799                               tile.edge_labels):
 800          ax.annotate(label, xy = (posn.x, posn.y), color = "b",
 801                      ha = "center", va = "center", fontsize = 8)
 802    if show_edge_labels:
 803      if show_edges:
 804        edges = self._get_edge_geoms(dist / 100 if offset_edges else 0)
 805        edges.plot(ax = ax, color = "forestgreen", ls = "dashed", lw = 1)
 806      else:
 807        edges = [e.get_geometry() for e in self.edges.values()]
 808      for l, e in zip(edges, self.edges.values()):
 809        c = l.centroid
 810        ax.annotate(e.label, xy = (c.x, c.y), ha = "center", va = "center",
 811          fontsize = 10 if show_edge_labels else 7, color = "k")
 812    if show_dual_tiles:
 813      gpd.GeoSeries(self.dual_tiles).buffer(
 814        -dist / 400, join_style = 2, cap_style = 3).plot(
 815          ax = ax, fc = "k", alpha = 0.2)
 816    pyplot.axis("off")
 817    return ax
 818
 819  def plot_tiling_symmetries(self, **kwargs):
 820    n = len(self.tile_matching_transforms)
 821    nc = int(np.ceil(np.sqrt(n)))
 822    nr = int(np.ceil(n / nc))
 823    gs = gpd.GeoSeries([t.shape for t in self.tiles])
 824    gsb = gs[:self.n_tiles]
 825    fig = pyplot.figure(figsize = (12, 12 * nr / nc))
 826    for i, tr in enumerate(self.tile_matching_transforms.values()):
 827      ax = fig.add_subplot(nr, nc, i + 1)
 828      gs.plot(ax = ax, fc = "b", alpha = 0.15, ec = "k", lw = 0.5)
 829      gsb.plot(ax = ax, fc = "#00000000", ec = "w", lw = 1, zorder = 2)
 830      gsm = gpd.GeoSeries([tr.apply(g) for g in gsb])
 831      gsm.plot(ax = ax, fc = "r", alpha = 0.2, lw = 0, ec = "r")
 832      tr.draw(ax, **kwargs)
 833      pyplot.axis("off")
 834
 835  def transform_geometry(self, new_topology:bool, apply_to_tiles:bool,
 836                         selector:str, type:str, **kwargs) -> "Topology":
 837    r"""Applies a specified transformation of elements in the Topology whose 
 838    labels match the selector parameter, optionally applying the transform to 
 839    update tile and optionally returning a new Topology object (or applying it 
 840    to this one). 
 841    
 842    Implemented in this way so that transformations can be applied one at a time
 843    without creating an intermediate set of new tiles, which may be invalid and
 844    fail. So, if you wish to apply (say) 3 transforms and generate a new 
 845    Topology leaving the existing one intact:
 846    
 847        new_topo = old_topo.transform_geometry(True,  False, "a", ...) \
 848                           .transform_geometry(False, False, "B", ...) \
 849                           .transform_geometry(False, True,  "C", ...)
 850    
 851    The first transform requests a new Topology, subsequent steps do not, and it
 852    is only the last step which attempts to create the new tile polygons.
 853    
 854    **kwargs supply named parameters for the requested transformation.
 855
 856    Args:
 857      new_topology (bool): if True returns a new Topology object, else returns 
 858        the current Topology modified.
 859      apply_to_tiles (bool): if True attempts to create new Tiles after the 
 860        transformation has been applied. Usually set to False, unless the last
 861        transformation in a pipeline, to avoid problems of topologically invalid
 862        tiles at intermediate steps.
 863      selector (str): label of elements to which to apply the transformation.   
 864        Note that all letters in the supplied string are checked, so you can 
 865        use e.g. "abc" to apply a transformation to edges labelled "a", "b" or 
 866        "c", or "AB" for vertices labelled "A" or "B".
 867      type (str): name of the type of transformation requested. Currently
 868        supported are `zigzag_edge`, `rotate_edge`, `push_vertex`, and 
 869        `nudge_vertex`. Keyword arguments for each are documented in the 
 870        corresponding methods.
 871
 872    Returns:
 873      Topology: if new_topology is True a new Topology based on this one with
 874        after transformation, if False this Topology is returned after the
 875        transformation.
 876    """
 877    print(
 878f"""CAUTION: new Topology will probably not be correctly labelled. To build a 
 879correct Topology, extract the tileable attribute and rebuild Topology from that.
 880""")
 881    topo = pickle.loads(pickle.dumps(self)) if new_topology else self
 882    transform_args = topo.get_kwargs(getattr(topo, type), **kwargs)
 883    if type == "zigzag_edge":
 884      for e in topo.edges.values():
 885        if e.label in selector:
 886          topo.zigzag_edge(e, **transform_args)
 887    # ---------------------------
 888    elif type == "rotate_edge":
 889      for e in topo.edges.values():
 890        if e.label in selector:
 891          topo.rotate_edge(e, **transform_args)
 892    # ---------------------------
 893    elif type == "push_vertex":
 894      pushes = {}
 895      for v in topo.vertices_in_tiles(topo.tiles[:topo.n_tiles]):
 896        if v.label in selector:
 897          pushes[v.base_ID] = topo.push_vertex(v, **transform_args)
 898      for base_ID, (dx, dy) in pushes.items():
 899        for v in [v for v in topo.points.values() if v.base_ID == base_ID]:
 900          v.point = affine.translate(v.point, dx, dy)
 901    # ---------------------------
 902    elif type == "nudge_vertex":
 903       for v in topo.points.values():
 904         if v.label in selector:
 905           topo.nudge_vertex(v, **transform_args)
 906    if apply_to_tiles:
 907      for t in topo.tiles:
 908        t.set_corners_from_edges()
 909    topo.tileable.tiles.geometry = gpd.GeoSeries(
 910      [topo.tiles[i].shape for i in range(topo.n_tiles)])
 911    topo.tileable.setup_regularised_prototile_from_tiles()
 912    return topo
 913
 914  def zigzag_edge(self, edge:Edge, start:str = "A",
 915                  n:int = 2, h:float = 0.5, smoothness:int = 0):
 916    """Applies a zigzag transformation to the supplied Edge. Currently this will
 917    only work correctly if h is even.
 918    
 919    TODO: make it possible for odd numbers of 'peaks' to work (this may require
 920    allowing bidirectional Edges, i.e. storing Edges in both directions so that
 921    all Tile edges are drawn CW). The `start` parameter is a temporary hack for 
 922    this
 923
 924    Args:
 925      edge (Edge): Edge to transform
 926      n (int, optional): number of zigs and zags in the edge. Defaults to 2.
 927      start (str, optional): label at one end of edge which is used to determine
 928        the sense of h, enabling C-curves with an odd number n of zigs and zags 
 929        to be applied. Defaults to 'A'.
 930      h (float, optional): width of the zig zags relative to edge length. 
 931        Defaults to 0.5.
 932      smoothness (int, optional): spline smoothness. 0 gives a zig zag proper,
 933        higher values will produce a sinusoid. Defaults to 0.
 934    """
 935    v0, v1 = edge.vertices[0], edge.vertices[1]
 936    if n % 2 == 1 and v0.label != start:
 937      h = -h
 938    ls = tiling_utils.zigzag_between_points(v0.point, v1.point, 
 939                                            n, h, smoothness)
 940    # remove current corners
 941    self.points = {k: v for k, v in self.points.items()
 942                   if not k in edge.get_corner_IDs()[1:-1]}
 943    # add the new ones
 944    new_corners = [self.add_vertex(geom.Point(xy)) for xy in ls.coords]
 945    edge.corners = edge.vertices[:1] + new_corners + edge.vertices[-1:]
 946    if not edge.right_tile is None:
 947      edge.right_tile.set_corners_from_edges(False)
 948    if not edge.left_tile is None:
 949      edge.left_tile.set_corners_from_edges(False)
 950
 951  def rotate_edge(self, edge:Edge, centre:str = "A", angle:float = 0) -> None:
 952    v0, v1 = edge.vertices[0], edge.vertices[1]
 953    if v0.label == centre:
 954      c = v0.point
 955    elif v1.label == centre:
 956      c = v1.point
 957    else:
 958      c = ls.centroid
 959    ls = affine.rotate(geom.LineString([v0.point, v1.point]), angle, origin = c)
 960    v0.point, v1.point = [geom.Point(c) for c in ls.coords]
 961
 962  def push_vertex(self, vertex:Vertex, push_d:float) -> tuple[float]:
 963    neighbours = [self.points[v] for v in vertex.neighbours]
 964    dists = [vertex.point.distance(v.point) for v in neighbours]
 965    x, y = vertex.point.x, vertex.point.y
 966    unit_vectors = [((x - v.point.x) / d, (y - v.point.y) / d)
 967                    for v, d in zip(neighbours, dists)]
 968    return  (push_d * sum([xy[0] for xy in unit_vectors]),
 969             push_d * sum([xy[1] for xy in unit_vectors]))
 970
 971  def nudge_vertex(self, vertex:Vertex, dx:float, dy:float):
 972    vertex.point = affine.translate(vertex.point, dx, dy)
 973
 974  def get_kwargs(self, fn:Callable, **kwargs) -> dict[Any]:
 975      args = list(inspect.signature(fn).parameters)
 976      return {k: kwargs.pop(k) for k in dict(kwargs) if k in args}
 977    
 978
 979  # ===========================================================================
 980  # potentially redundant code
 981  # ===========================================================================
 982  def _label_tiles(self):
 983    """Labels the base tile vertices and edges at a per tile level, then copies 
 984    to corresponding tiles in the radius-1 patch.
 985    
 986    NOTE: this is effectively legacy code, since its purpose originally was as 
 987    a (flawed) basis for determining vertex and edge labels.
 988    """
 989    first_letter = 0
 990    for i, group in enumerate(self.shape_groups):
 991      tiles = [t for t in self.tiles[:self.n_tiles] if t.ID in group]
 992      s = Symmetries(self.tiles[group[0]].shape)
 993      vlabels = list(s.get_unique_labels(offset = first_letter)["rotations"][0])
 994      elabels = self._get_edge_labels_from_vertex_labels(vlabels)
 995      for tile in tiles:
 996        tile.vertex_labels = vlabels
 997        tile.edge_labels = elabels
 998      first_letter = first_letter + len(set(vlabels))
 999    for tile in self.tiles: #[self.n_tiles:]:
1000      tile.vertex_labels = self.tiles[tile.base_ID].vertex_labels
1001      tile.edge_labels = self.tiles[tile.base_ID].edge_labels
1002
1003  def _get_edge_labels_from_vertex_labels(self, vlabels:list[str]) -> list[str]:
1004    r"""Given a list of vertex labels, returns a list of edge labels such that
1005    each distinct pair of vertex labels at the ends of an edge will be given a
1006    distinct lower case letter label. E.g., A B C B A will yield a+ b+ b- a- c 
1007    from AB -> a+ BC -> b+ CB -> b- BA -> a- AA -> c, where + designated CW
1008    traversal, - CCW, and no symbol means that edge appears only once.
1009    
1010    This may yet have use in more elaborate edge transformations. (See G&S 87 
1011    on tile incidence symbols.)
1012
1013    Args:
1014      vlabels (list[str]): ordered list of vertex labels.
1015
1016    Returns:
1017      list[str]: corresponding list of edge labels.
1018    """
1019    AB_labels = [a + b for a, b in zip(vlabels, vlabels[1:] + vlabels[:1])]
1020    letter = ALPHABET.index(min(list(vlabels)))
1021    edge_labels = {}
1022    for AB_label in AB_labels:
1023      if not AB_label in edge_labels:
1024        if AB_label[::-1] in edge_labels:
1025          label = edge_labels[AB_label[::-1]]
1026          edge_labels[AB_label] = label + "-"
1027          edge_labels[AB_label[::-1]] = label + "+"
1028        else:
1029          edge_labels[AB_label] = alphabet[letter]
1030          letter = letter + 1
1031    return [edge_labels[i] for i in AB_labels]
1032
1033  # ===========================================================================
1034  # redundant code
1035  # ===========================================================================
1036  def _label_vertices(self):
1037    """Labels vertices based on cycle of tile vertex labels around them.
1038    """
1039    uniques = set()
1040    # first label vertices in the core tiles
1041    vs = set()
1042    for t in self.tiles[:self.n_tiles]:
1043      vs = vs.union(t.get_corner_IDs())
1044    # Note: sorting is important for repeatable labelling!
1045    for vi in sorted(vs):
1046      v = self.points[vi]
1047      label = ""
1048      for tile in v.tiles:
1049        label = label + \
1050          tile.vertex_labels[tile.corners.index(v)]
1051      # TODO: resolve this question: cyclic sort seems more correct, but 
1052      # neither approach seems to work in all cases... see esp. cyclic sort 
1053      # applied to the cheese sandwich tiling.
1054      v.label = "".join(self._cyclic_sort_first(list(label)))
1055      # v.label = "".join(sorted(list(label)))
1056      # v.label = min(list(label))
1057      uniques.add(v.label)
1058    for vi in sorted(vs):
1059      v = self.points[vi]
1060      v.label = ALPHABET[sorted(uniques).index(v.label)]
1061    # now copy to corresponding vertices in the rest of tiling
1062    for ti, t in enumerate(self.tiles):
1063      if ti >= self.n_tiles:
1064        t0 = self.tiles[ti % self.n_tiles]
1065        for v0, v in zip(t0.corners, t.corners):
1066          # Vertices may appear in more than one tile, only label once!
1067          if self.points[v.ID].label == "":
1068            self.points[v.ID].label = self.points[v0.ID].label
1069
1070  def _label_edges(self):
1071    """Labels edges based on the tile edge label on each side.
1072    """
1073    uniques = set()
1074    labelled = set()
1075    # first label base tile edges from the central tile unit
1076    for t in self.tiles[:self.n_tiles]:
1077      for e in t.edges:
1078        rt, lt = e.right_tile, e.left_tile
1079        rlab, llab = rt.get_edge_label(e), lt.get_edge_label(e)
1080        e.label = "".join(sorted([rlab, llab]))
1081        uniques.add(e.label)
1082        labelled.add(e.ID)
1083    for ei in sorted(labelled):
1084      e = self.edges[ei]
1085      e.label = alphabet[sorted(uniques).index(e.label)]
1086    # now copy to corresponding edges in the rest of tiling
1087    for ti, t in enumerate(self.tiles):
1088      if ti >= self.n_tiles:
1089        t0 = self.tiles[ti % self.n_tiles]
1090        for e0, e in zip(t0.edges, t.edges):
1091          # edges may appear in more than one tile, only label once!
1092          if e.label == "":
1093            e.label = e0.label
1094
1095  def _cyclic_sort_first(self, lst:Iterable) -> Iterable:
1096    """Returns supplied Iterable in canonical cyclic sorted form. E.g. the 
1097    sequence ACABD, yields 5 possible cycles ACABD, CABDA, ABDAC, BDACA, and 
1098    DACAB. The lexically first of these is ABDAC, which would be returned. 
1099
1100    Args:
1101        lst (Iterable): an Iterable of elements (usu. strings).
1102
1103    Returns:
1104      Iterable: the lexically first of the possible cycles in the supplied 
1105        iterable.
1106    """
1107    cycle = lst + lst
1108    n = len(lst)
1109    cycles = [cycle[i:i + n] for i in range(n)]
1110    return sorted(cycles)[0]