weavingspace.topology

   1"""Classes for working with topology of tilings.
   2
   3Together the `Topology`, `Tile`, `Vertex`, `Edge`, and
   4`weavingspace.symmetry.Symmetries` classes enable extraction of the topological
   5structure of periodic `weavingspace.tileable.Tileable` objects so that
   6modification of equivalent tiles can be carried out while retaining
   7tileability.
   8
   9It is important to note that these are not fully generalised classes and
  10methods, that is, the Topology object that is supported is not a permanent
  11'backing' data structure for our Tileable objects. While it might become that
  12in time, it is not yet such a data structure. Instead usage is
  13
  14    tile = TileUnit(...)
  15    topology = Topology(tile)
  16    topology = topology.transform_*(...)
  17    new_tile = topology.tile_unit
  18
  19A Topology plot function is provided for a user to be able to see what they are
  20doing, because how edges and vertices in a tiling are labelled under tile
  21equivalences is an essential step in the process.
  22
  23Note that these classes do not precisely represent the distinctions made in the
  24mathematical literature between tiling vertices and tile corners, or between
  25tiling edges and tile sides.
  26"""
  27from __future__ import annotations
  28
  29import inspect
  30import itertools
  31import pickle
  32import string
  33from collections import defaultdict
  34from typing import TYPE_CHECKING
  35
  36import geopandas as gpd
  37import matplotlib.pyplot as plt
  38import networkx as nx
  39import numpy as np
  40import shapely.affinity as affine
  41import shapely.geometry as geom
  42from scipy import interpolate
  43
  44from weavingspace import (
  45  ShapeMatcher,
  46  Symmetries,
  47  Tileable,
  48  Transform,
  49  tiling_utils,
  50)
  51
  52if TYPE_CHECKING:
  53  from collections.abc import Callable, Iterable
  54
  55# all two letter combinations of the alphabet for labelling
  56LABELS = \
  57  list(string.ascii_letters.upper()) + \
  58  ["".join(x) for x in
  59   itertools.combinations(list(string.ascii_letters.upper()), 2)]
  60labels = \
  61  list(string.ascii_letters.lower()) + \
  62  ["".join(x) for x in
  63   itertools.combinations(list(string.ascii_letters.lower()), 2)]
  64
  65
  66class Topology:
  67  """Class to represent topology of a Tileable object.
  68
  69  NOTE: It is important that get_local_patch return the tileable elements and
  70  the translated copies in consistent sequence, i.e. if there are (say) four
  71  tiles in the unit, the local patch should be 1 2 3 4 1 2 3 4 1 2 3 4 ... and
  72  so on. This is because self.tiles[i % n_tiles] is frequently used to
  73  reference the base unit Tile which corresponds to self.tiles[i].
  74  """
  75
  76  tileable: Tileable
  77  """the Tileable on which the topology will be based."""
  78  tiles: list[Tile]
  79  """list of the Tiles in the topology. We use polygons returned by the
  80  tileable.get_local_patch method for these. That is the base tiles and 8
  81  adjacent copies (for a rectangular tiling), or 6 adjacent copies (for a
  82  hexagonal tiling)."""
  83  points: dict[int, Vertex]
  84  """dictionary of all points (vertices and corners) in the tiling, keyed by
  85  Vertex ID."""
  86  edges: dict[tuple[int,int], Edge]
  87  """dictionary of the tiling edges, keyed by Edge ID."""
  88  unique_tile_shapes: list[geom.Polygon]
  89  """a 'reference' tile shape one per tile shape (up to vertices, so two tiles
  90  might be the same shape, but one might have extra vertices induced by the
  91  tiling and hence is a different shape under this definition)."""
  92  dual_tiles: list[geom.Polygon]
  93  """list of geom.Polygons from which a dual tiling might be constructed."""
  94  n_tiles: int = 0
  95  """number of tiles in the base Tileable (retained for convenience)."""
  96  shape_groups: list[list[int]]
  97  """list of lists of tile IDs distinguished by shape and optionally tile_id"""
  98  tile_matching_transforms: list[tuple[float]]
  99  """shapely transform tuples that map tiles onto other tiles"""
 100  tile_transitivity_classes: list[tuple[int]]
 101  """list of lists of tile IDs in each transitivity class"""
 102  vertex_transitivity_classes: list[list[int]]
 103  """list of lists of vertex IDs in each transitivity class"""
 104  edge_transitivity_classes: list[list[tuple[int]]]
 105  """list of lists of edge IDs in each transitivity class"""
 106
 107  def __init__(
 108      self,
 109      unit: Tileable,
 110      ignore_tile_ids:bool = True,
 111    ) -> None:
 112    """Class constructor.
 113
 114    Args:
 115      unit (Tileable): the Tileable whose topology is required.
 116      ignore_tile_ids (bool): (EXPERIMENTAL) if True then only consider the tile
 117        shapes, not labels. If False consider any labels. Defaults to True.
 118
 119    """
 120    # Note that the order of these setup steps is fragile sometimes
 121    # not obviously so.
 122    self.tileable = unit # keep this for reference
 123    self.n_tiles = self.tileable.tiles.shape[0]
 124    self._initialise_points_into_tiles()
 125    self._setup_vertex_tile_relations()
 126    self._setup_edges()
 127    self._copy_base_tiles_to_patch()
 128    self._assign_vertex_and_edge_base_IDs()
 129    self._identify_distinct_tile_shapes(ignore_tile_ids)
 130    self._find_tile_transitivity_classes(ignore_tile_ids)
 131    self._find_vertex_transitivity_classes(ignore_tile_ids)
 132    self._find_edge_transitivity_classes(ignore_tile_ids)
 133    self.generate_dual()
 134
 135
 136  def __str__(self) -> str:
 137    """Return string representation of this Topology.
 138
 139    Returns:
 140      str: a message that recommends examining the tiles, points and edges
 141        attributes.
 142
 143    """
 144    return (f"""Topology of Tileable with {self.n_tiles} tiles.\n
 145            Examine .tiles, .points and .edges for more details.""")
 146
 147
 148  def __repr__(self) -> str:
 149    return str(self)
 150
 151
 152  def _initialise_points_into_tiles(self, debug:bool = False) -> None:
 153    """Set up dictionary of unique point locations and assign them to Tiles.
 154
 155    Args:
 156      debug (bool): if True prints useful debugging information.
 157
 158    """
 159    shapes = self.tileable.get_local_patch(r = 1, include_0 = True).geometry
 160    shapes = [tiling_utils.get_clean_polygon(s) for s in shapes]
 161    self.n_1_order_patch = len(shapes)
 162    labels = list(self.tileable.tiles.tile_id) * (len(shapes) // self.n_tiles)
 163    self.tiles = []
 164    self.points = {}
 165    for (i, shape), label in zip(enumerate(shapes), labels, strict = True):
 166      tile = Tile(i)
 167      tile.label = label
 168      tile.base_ID = tile.ID % self.n_tiles
 169      self.tiles.append(tile)
 170      tile.corners = []
 171      corners = tiling_utils.get_corners(shape, repeat_first = False)
 172      for c in corners:
 173        prev_vertex = None
 174        for p in reversed(self.points.values()):
 175          if c.distance(p.point) <= 2 * tiling_utils.RESOLUTION:
 176            # found an already existing vertex, so add to the tile and break
 177            prev_vertex = p
 178            tile.corners.append(prev_vertex)
 179            break
 180        if prev_vertex is None:
 181          # new vertex, add it to the topology dictionary and to the tile
 182          v = self.add_vertex(c)
 183          tile.corners.append(v)
 184          if debug:
 185            print(f"Added new Vertex {v} to Tile {i}")
 186
 187
 188  def _setup_vertex_tile_relations(self, debug:bool = False) -> None:
 189    """Determine relations between vertices and tiles.
 190
 191    In particular vertices along tile edges that are not yet included in their
 192    list of vertices are added. Meanwhile vertex lists of incident tiles and
 193    neighbours are set up.
 194
 195    Args:
 196      debug (bool): if True prints debugging information.
 197
 198    """
 199    # we do this for all tiles in the radius-1 local patch
 200    for tile in self.tiles:
 201      if debug:
 202        print(f"Checking for vertices incident on Tile {tile.ID}")
 203      corners = []
 204      # performance is much improved using the vertex IDs
 205      initial_corner_IDs = tile.get_corner_IDs()
 206      # we need the current shape (not yet set) to check for incident vertices
 207      shape = geom.Polygon([c.point for c in tile.corners])
 208      # get points incident on tile boundary, not already in the tile corners
 209      new_points = [v for v in self.points.values()
 210                    if v.ID not in initial_corner_IDs and
 211                    v.point.distance(shape) <= 2 * tiling_utils.RESOLUTION]
 212      # iterate over sides of the tile to see which the vertex is incident on
 213      for c1, c2 in zip(tile.corners, tile.corners[1:] + tile.corners[:1],
 214                        strict = True):
 215        all_points = [c1, c2]
 216        if len(new_points) > 0:
 217          if debug:
 218            print(f"{[v.ID for v in new_points]} incident on tile")
 219          ls = geom.LineString([c1.point, c2.point])
 220          to_insert = [v for v in new_points
 221                       if v.point.distance(ls) <= 2 * tiling_utils.RESOLUTION]
 222          if len(to_insert) > 0:
 223            # sort by distance along the side
 224            d_along = sorted([(ls.line_locate_point(v.point), v)
 225                              for v in to_insert], key = lambda x: x[0])
 226            to_insert = [v for d, v in d_along]
 227            all_points = all_points[:1] + to_insert + all_points[1:]
 228        corners.extend(all_points[:-1])
 229        for x1, x2 in zip(all_points[:-1], all_points[1:], strict = True):
 230          # x2 will add the tile and neigbour when we get to the next side
 231          # every vertex gets a turn!
 232          x1.add_tile(tile)
 233          x1.add_neighbour(x2.ID)
 234      tile.corners = corners
 235      tile.set_shape_from_corners()
 236
 237
 238  def _setup_edges(
 239      self,
 240      debug:bool = False,
 241    ) -> None:
 242    """Set up the tiling edges.
 243
 244    First vertices in the base tiles are classified as tiling vertices or not -
 245    only these can be classified reliably (e.g vertices on the perimeter are
 246    tricky).. Up to here all vertices have been considered tiling vertices.
 247
 248    Second edges are created by traversing tile corner lists. Edges are stored
 249    once only by checking for edges in the reverse direction already in the
 250    edges dictionary. Edge right and left tiles are initialised.
 251
 252    Third tile edge direction lists are initialised.
 253
 254    Args:
 255      debug (bool): if True print debug messages. Defaults to False.
 256
 257    """
 258    # classify vertices in the base tiles
 259    for tile in self.tiles[:self.n_tiles]:
 260      for v in tile.corners:
 261        v.is_tiling_vertex = len(v.neighbours) > 2
 262    if debug:
 263      print("Classified base tile vertices")
 264    self.edges = {}
 265    for tile in self.tiles:
 266      if debug:
 267        print(f"Adding edges from Tile {tile.ID}")
 268      tile.edges = []
 269      vertices = [v for v in tile.corners if v.is_tiling_vertex]
 270      # note that through here finding ints in lists is much faster than
 271      # finding Vertex objects, hence we use lists of IDs not Vertex objects
 272      if len(vertices) > 1:
 273        for v1, v2 in zip(vertices, vertices[1:] + vertices[:1], strict = True):
 274          corner_IDs = tile.get_corner_IDs()
 275          idx1 = corner_IDs.index(v1.ID)
 276          idx2 = corner_IDs.index(v2.ID)
 277          if idx1 < idx2:
 278            corners = corner_IDs[idx1:(idx2 + 1)]
 279          else:
 280            corners = corner_IDs[idx1:] + corner_IDs[:(idx2 + 1)]
 281          ID = (corners[0], corners[-1])
 282          if ID not in self.edges:
 283            # check that reverse direction edge is not present first
 284            r_ID = ID[::-1]
 285            if r_ID in self.edges:
 286              # if it is, then add it and set left_tile
 287              if debug:
 288                print(f"reverse edge {r_ID} found")
 289              e = self.edges[r_ID]
 290              e.left_tile = tile
 291              tile.edges.append(e)
 292            else:
 293              # we've found a new edge so make and add it
 294              if debug:
 295                print(f"adding new edge {corners}")
 296              e = self.add_edge([self.points[c] for c in corners])
 297              e.right_tile = tile
 298              tile.edges.append(e)
 299      # initialise the edge direction information in the tile
 300      tile.set_edge_directions()
 301
 302
 303  def _assign_vertex_and_edge_base_IDs(self) -> None:
 304    """Assign the base_ID attributes of vertices and edges.
 305
 306    These allow us to determine corrspondences between vertices and edges in
 307    the 'base' tiles in the Topology tileable, and those we have added at
 308    radius 1 for labelling and visualisation.
 309    """
 310    self._assign_vertex_base_IDs()
 311    self._assign_edge_base_IDs()
 312
 313
 314  def _assign_vertex_base_IDs(self) -> None:
 315    """Assign base_ID attribute of vertices."""
 316    # assign vertex base_ID frome the core tiles
 317    for tile0 in self.tiles[:self.n_tiles]:
 318      for v in tile0.get_corner_IDs():
 319        self.points[v].base_ID = v
 320    # assign others from their corresponding vertex in the core
 321    for t0 in self.tiles[:self.n_tiles]:
 322      for t1 in self.tiles[self.n_tiles:]:
 323        if t1.base_ID == t0.base_ID:
 324          for v0, v1 in zip(t0.corners, t1.corners, strict = True):
 325            v1.base_ID = v0.base_ID
 326
 327
 328  def _assign_edge_base_IDs(self) -> None:
 329    """Assign base_ID attribute of edges, based on base_IDs of vertices."""
 330    for tile0 in self.tiles[:self.n_tiles]:
 331      for e in tile0.get_edge_IDs():
 332        self.edges[e].base_ID = e
 333    for t0 in self.tiles[:self.n_tiles]:
 334      for t1 in self.tiles[self.n_tiles:]:
 335        if t1.base_ID == t0.base_ID:
 336          for e0, e1 in zip(t0.edges, t1.edges, strict = True):
 337            e1.base_ID = e0.base_ID
 338
 339
 340  def _copy_base_tiles_to_patch(self) -> None:
 341    """Copy attributes of base tiles to corresponding tiles in radius-1 patch.
 342
 343    This requires:
 344
 345    First, inserting any additional corners in the base tiles not found in the
 346    radius-1 tiles.
 347
 348    Second, any vertices in the base tiles that are NOT tiling vertices are
 349    applied to radius-1 tiles leading to merging of some edges.
 350    """
 351    # the number of tiles in the base + radius-1
 352    n_r1 = len(self.tiles)
 353    # first add any missing vertices to the non-base tiles
 354    # we add all the missing vertices before doing any merges
 355    for base in self.tiles[:self.n_tiles]:
 356      for other in self.tiles[base.ID:n_r1:self.n_tiles]:
 357        self._match_reference_tile_vertices(base, other)
 358    # then merge any edges that meet at a corner
 359    for base in self.tiles[:self.n_tiles]:
 360      for other in self.tiles[base.ID:n_r1:self.n_tiles]:
 361        self._match_reference_tile_corners(base, other)
 362
 363
 364  def _match_reference_tile_vertices(
 365      self,
 366      tile1:Tile,
 367      tile2:Tile,
 368    ) -> None:
 369    """Add vertices to tile2 so it matches tile1 adjusting edges as required.
 370
 371    This assumes the tiles are the same shape, but that tile2 may be missing
 372    some tiling vertices along some edges.
 373
 374    Args:
 375      tile1 (Tile): reference tile.
 376      tile2 (Tile): tile to change.
 377
 378    """
 379    to_add = len(tile1.corners) - len(tile2.corners)
 380    while to_add > 0:
 381      # find the reference x-y offset
 382      dxy = (tile2.centre.x - tile1.centre.x, tile2.centre.y - tile1.centre.y)
 383      for i, t1c in enumerate([c.point for c in tile1.corners]):
 384        t2c = tile2.corners[i % len(tile2.corners)].point
 385        if abs((t2c.x - t1c.x) - dxy[0]) > 10 * tiling_utils.RESOLUTION or \
 386           abs((t2c.y - t1c.y) - dxy[1]) > 10 * tiling_utils.RESOLUTION:
 387          # add vertex to t2 by copying the t1 vertex appropriately offset
 388          # note that this might alter the length of t2.corners
 389          v = self.add_vertex(geom.Point(t1c.x + dxy[0], t1c.y + dxy[1]))
 390          v.is_tiling_vertex = True
 391          old_edge, new_edges = tile2.insert_vertex_at(v, i)
 392          del self.edges[old_edge]
 393          for e in new_edges:
 394            self.edges[e.ID] = e
 395          to_add = to_add - 1
 396
 397
 398  def _match_reference_tile_corners(
 399      self,
 400      tile1:Tile,
 401      tile2:Tile,
 402    ) -> None:
 403    """Make vertices that are corners in tile1 corners in tile2.
 404
 405    Edges are merged as required.
 406
 407    Args:
 408        tile1 (Tile): reference tile.
 409        tile2 (Tile): tile to make match.
 410
 411    """
 412    vs_to_change = [
 413      vj for vi, vj in zip(tile1.corners, tile2.corners, strict = True)
 414      if not vi.is_tiling_vertex and vj.is_tiling_vertex]
 415    if len(vs_to_change) > 0:
 416      for v in vs_to_change:
 417        v.is_tiling_vertex = False
 418        # it's a corner not an edge so will have no more than 2 v.tiles
 419        old_edges, new_edge = v.tiles[0].merge_edges_at_vertex(v)
 420        for e in old_edges:
 421          del self.edges[e]
 422        self.edges[new_edge.ID] = new_edge
 423
 424  def _identify_distinct_tile_shapes(
 425      self,
 426      ignore_tile_id_labels:bool = True,
 427    ) -> None:
 428    """Identify unique tiles based on their symmetries and shapes.
 429
 430    At the same time assembles a list of the affine transforms under which
 431    matches occurs since these are potential symmetries of the tiling.
 432
 433    TODO: reimplement consideration of tile_id
 434
 435    Args:
 436      ignore_tile_id_labels (bool): if True only the shape of tiles matters; if
 437        False the tile_id label is also considered. Defaults to True.
 438
 439    """
 440    if ignore_tile_id_labels:
 441      matches = {}
 442      offsets = {}
 443      for tile in self.tiles[:self.n_tiles]:
 444        matches[tile.base_ID] = [tile.base_ID]
 445        matched = False
 446        s = Symmetries(tile.shape)
 447        for other in self.tiles[:self.n_tiles]:
 448          if other.ID > tile.ID:
 449            offset = s.get_corner_offset(other.shape)
 450            if offset is not None:
 451              offsets[tile.base_ID] = offset
 452              matches[tile.base_ID].append(other.base_ID)
 453              matched = True
 454        if not matched:
 455          offsets[tile.base_ID] = 0
 456      base_groups = list(
 457        nx.connected_components(nx.from_dict_of_lists(matches)))
 458      self.shape_groups = []
 459      for i, group in enumerate(base_groups):
 460        full_group = []
 461        for tile in self.tiles:
 462          if tile.base_ID in group:
 463            tile.shape_group = i
 464            tile.offset_corners(offsets[tile.base_ID])
 465            full_group.append(tile.ID)
 466        self.shape_groups.append(full_group)
 467    else:
 468      self.shape_groups = []
 469      for ti in self.tiles[:self.n_tiles]:
 470        self.shape_groups.append(
 471          [tj.ID for tj in self.tiles if tj.base_ID == ti.base_ID])
 472      for i, group in enumerate(self.shape_groups):
 473        for j in group:
 474          self.tiles[j].shape_group = i
 475
 476  def _find_tile_transitivity_classes(
 477      self,
 478      ignore_tile_id_labels:bool = True,
 479    ) -> None:
 480    """Find tiles equivalent under symmetries.
 481
 482    Also update the tile_matching_transforms attribute to contain only those
 483    transforms that pass this test.
 484
 485    Args:
 486      ignore_tile_id_labels (bool): if True then consider only shapes; if False
 487        also consider labels of tiles. Defaults to True.
 488
 489    """
 490    self.tile_matching_transforms = \
 491      self.get_potential_symmetries(ignore_tile_id_labels)
 492    if ignore_tile_id_labels:
 493      base_tiles = self.tiles[:self.n_tiles]
 494      # it is quicker (should be!) to only do within shape group tests
 495      # often there is only one when it will make no difference
 496      by_group_equivalent_tiles = []
 497      # maintain a set of transforms still potentially tiling symmetries
 498      poss_transforms = set(self.tile_matching_transforms.keys())
 499      # and a dictionary of booleans tracking which transforms are still valid
 500      eq_under_transform = dict.fromkeys(poss_transforms, True)
 501      for g, _ in enumerate(self.shape_groups):
 502        by_group_equivalent_tiles.append(set())
 503        source_tiles = [tile for tile in base_tiles if tile.shape_group == g]
 504        target_tiles = [tile for tile in self.tiles if tile.shape_group == g]
 505        for tr in poss_transforms:
 506          transform = self.tile_matching_transforms[tr].transform
 507          matched_tiles = {}
 508          eq_under_transform[tr] = True
 509          for source_tile in source_tiles:
 510            matched_tile_id = self._match_geoms_under_transform(
 511              source_tile, target_tiles, transform)
 512            if matched_tile_id == -1:
 513              eq_under_transform[tr] = False
 514              break
 515            matched_tiles[source_tile.ID] = matched_tile_id # actually a base_ID
 516          if eq_under_transform[tr]:
 517            for k, v in matched_tiles.items():
 518              # here we record the transform, in case it is later invalidated
 519              by_group_equivalent_tiles[g].add((tr, k, v))
 520        # remove valid transforms that didn't make it through this group
 521        poss_transforms = {t for t, x in eq_under_transform.items() if x}
 522      # compile equivalences from all groups made under still valid transforms
 523      # a dict of sets so singletons aren't lost in finding connected components
 524      equivalents = {i: set() for i in range(self.n_tiles)}
 525      for group_equivalents in by_group_equivalent_tiles:
 526        for (tr, tile_i, tile_j) in group_equivalents:
 527          if tr in poss_transforms:
 528            equivalents[tile_i].add(tile_j)
 529      self.tile_matching_transforms = {
 530        k: v for k, v in self.tile_matching_transforms.items()
 531        if k in poss_transforms}
 532      self.tile_transitivity_classes = []
 533      equivalents = nx.connected_components(nx.from_dict_of_lists(equivalents))
 534      for c, base_IDs in enumerate(equivalents):
 535        transitivity_class = []
 536        for tile in self.tiles:
 537          if tile.base_ID in base_IDs:
 538            transitivity_class.append(tile.ID)
 539            tile.transitivity_class = c
 540        self.tile_transitivity_classes.append(transitivity_class)
 541    else:
 542      # transitivity classes are just the individual tiles
 543      self.tile_transitivity_classes = []
 544      for i, tile in enumerate(self.tiles):
 545        tile.transitivity_class = tile.base_ID
 546        if i < self.n_tiles:
 547          self.tile_transitivity_classes.append([tile.ID])
 548        else:
 549          self.tile_transitivity_classes[tile.base_ID].append(tile.ID)
 550
 551
 552  def get_potential_symmetries(
 553      self,
 554      ignore_tile_id_labels:bool = True,
 555    ) -> dict[int, Transform]:
 556    """Assemble potential symmetries from symmetries of prototile and tiles.
 557
 558    Also remove any duplicates that result. The result is assigned to the
 559    tile_matching_transforms attribute.
 560
 561    TODO: consider retaining the Symmetry objects as these carry additional
 562    information that might facilitate labelling under a limited number of the
 563    symmetries not all of them.
 564
 565    Returns:
 566      dict[int, tuple[float]]: dictionary of the symmetries (transforms
 567        actually) in shapely affine transform 6-tuple format.
 568
 569    """
 570    self.tile_matching_transforms = {
 571      k: Transform("translation", 0, geom.Point(0, 0), v,
 572                   tiling_utils.get_translation_transform(v[0], v[1]))
 573      for k, v in enumerate(self.tileable.get_vectors())}
 574    if ignore_tile_id_labels:
 575      n_symmetries = len(self.tile_matching_transforms)
 576      ptile = self.tileable.prototile.loc[0, "geometry"]
 577      for tr in ShapeMatcher(ptile).get_polygon_matches(ptile):
 578        if tr.transform_type not in ["identity", "translation"]:
 579          self.tile_matching_transforms[n_symmetries] = tr
 580          n_symmetries = n_symmetries + 1
 581      for tile in self.tiles[:self.n_tiles]:
 582        for tr in ShapeMatcher(tile.shape).get_polygon_matches(tile.shape):
 583          if tr.transform_type not in ["identity", "translation"]:
 584            self.tile_matching_transforms[n_symmetries] = tr
 585            n_symmetries = n_symmetries + 1
 586      for tile in self.tiles[:self.n_tiles]:
 587        sm = ShapeMatcher(tile.shape)
 588        transforms = [sm.get_polygon_matches(self.tiles[i].shape)
 589          for i in self.shape_groups[tile.shape_group] if i < self.n_tiles]
 590        for tr in itertools.chain(*transforms):
 591          if tr.transform_type not in ["identity", "translation"]:
 592            self.tile_matching_transforms[n_symmetries] = tr
 593            n_symmetries = n_symmetries + 1
 594      self.tile_matching_transforms = self._remove_duplicate_symmetries(
 595        self.tile_matching_transforms)
 596    return self.tile_matching_transforms
 597
 598
 599  def _remove_duplicate_symmetries(
 600      self,
 601      transforms:dict[int,Transform],
 602    ) -> dict[int,Transform]:
 603    """Filter list of shapely affine transforms to remove duplicates.
 604
 605    Args:
 606      transforms (dict[int,Transform]): dictionary of Transforms to filter.
 607
 608    Returns:
 609      dict[int,Transform]: the filtered dictionary with duplicates removed.
 610
 611    """
 612    uniques = {}
 613    for k, v in transforms.items():
 614      already_exists = False
 615      for u in uniques.values():
 616        already_exists = np.allclose(
 617          v.transform, u.transform, atol = 1e-4, rtol = 1e-4)
 618        if already_exists:
 619          break
 620      if not already_exists:
 621        uniques[k] = v
 622    return uniques
 623
 624
 625  def _find_vertex_transitivity_classes(
 626      self,
 627      ignore_tile_id_labels:bool = True,
 628    ) -> None:
 629    """Find vertex transitivity classes.
 630
 631    This function checks which vertices align with which others under transforms
 632    in the tile_matching_transforms attribute. The process need only determine
 633    the classes for vertices in the core tileable.tiles, then assign those to
 634    all vertices by matched base_ID.
 635    """
 636    if ignore_tile_id_labels:
 637      equivalent_vertices = defaultdict(set)
 638      base_vertices = [v for v in
 639                      self.vertices_in_tiles(self.tiles[:self.n_tiles])
 640                      if v.is_tiling_vertex]
 641      for transform in self.tile_matching_transforms.values():
 642        for v in base_vertices:
 643          equivalent_vertices[v.ID].add(v.ID)
 644          match_ID = self._match_geoms_under_transform(
 645            v, base_vertices, transform.transform)
 646          if match_ID != -1:
 647            equivalent_vertices[v.ID].add(match_ID)
 648      equivalent_vertices = self._get_exclusive_supersets(
 649        [tuple(sorted(s)) for s in equivalent_vertices.values()])
 650      self.vertex_transitivity_classes = defaultdict(list)
 651      for c, vclass in enumerate(equivalent_vertices):
 652        for v in self.points.values():
 653          if v.base_ID in vclass:
 654            v.transitivity_class = c
 655            self.vertex_transitivity_classes[c].append(v.ID)
 656      self.vertex_transitivity_classes = list(
 657        self.vertex_transitivity_classes.values())
 658      # label vertices based on their transitivity class
 659      for v in self.points.values():
 660        if v.is_tiling_vertex:
 661          v.label = LABELS[v.transitivity_class]
 662    else:
 663      self.vertex_transitivity_classes = defaultdict(list)
 664      for v in self.points.values():
 665        if v.is_tiling_vertex:
 666          self.vertex_transitivity_classes[v.base_ID].append(v.ID)
 667          v.transitivity_class = v.base_ID
 668      self.vertex_transitivity_classes = list(
 669        self.vertex_transitivity_classes.values())
 670      for v in self.points.values():
 671        if v.is_tiling_vertex:
 672          v.label = LABELS[v.transitivity_class]
 673
 674
 675  def _find_edge_transitivity_classes(
 676      self,
 677      ignore_tile_id_labels:bool = True,
 678    ) -> None:
 679    """Find edge transitivity classes.
 680
 681    This function works by checking which edges align with which others under
 682    transforms in the tile_matching_transforms attribute. The process need only
 683    determine the classes for edges in the core tileable.tiles, then assign
 684    those to all edges by matched base_ID.
 685
 686    TODO: Note that this code is identical to the vertex transitivity code
 687    so it might make sense to merge.
 688    """
 689    if ignore_tile_id_labels:
 690      equivalent_edges = defaultdict(set)
 691      base_edges = self.edges_in_tiles(self.tiles[:self.n_tiles])
 692      for transform in self.tile_matching_transforms.values():
 693        for e in base_edges:
 694          equivalent_edges[e.ID].add(e.ID)
 695          match_id = self._match_geoms_under_transform(
 696            e, base_edges, transform.transform)
 697          if match_id != -1:
 698            equivalent_edges[e.ID].add(match_id)
 699      equivalent_edges = self._get_exclusive_supersets(
 700        [tuple(sorted(s)) for s in equivalent_edges.values()])
 701      self.edge_transitivity_classes = defaultdict(list)
 702      for c, eclass in enumerate(equivalent_edges):
 703        for e in self.edges.values():
 704          if e.base_ID in eclass:
 705            e.transitivity_class = c
 706            self.edge_transitivity_classes[c].append(e.ID)
 707      self.edge_transitivity_classes = list(
 708        self.edge_transitivity_classes.values())
 709      # label edges based on their transitivity class
 710      for e in self.edges.values():
 711        e.label = labels[e.transitivity_class]
 712    else:
 713      self.edge_transitivity_classes = defaultdict(list)
 714      for e in self.edges.values():
 715        self.edge_transitivity_classes[e.base_ID].append(e.ID)
 716      self.edge_transitivity_classes = list(
 717        self.edge_transitivity_classes.values())
 718      for i, eclass in enumerate(self.edge_transitivity_classes):
 719        for e in eclass:
 720          self.edges[e].transitivity_class = i
 721          self.edges[e].label = labels[i]
 722
 723
 724  def _match_geoms_under_transform(
 725      self,
 726      geom1:Tile|Vertex|Edge,
 727      geoms2:list[Tile|Vertex|Edge],
 728      transform:tuple[float,...],
 729    ) -> int|tuple[int]:
 730    """Determine if a geometry maps onto any in a patch under supplied symmetry.
 731
 732    Args:
 733      geom1 (Tile|Vertex|Edge): element whose geometry we want to match.
 734      geoms2 (list[Tile|Vertex|Edge]): set of elements among which a
 735        match is sought.
 736      transform (tuple[float]): shapely affine transform 6-tuple to apply.
 737
 738    Returns:
 739      int|tuple[int,int]: ID of the element in patch that matches the geom1
 740        element under the transform if one exists, otherwise returns -1. For
 741        edges note that the ID is a tuple.
 742
 743    """
 744    match_id = -1
 745    for geom2 in geoms2:
 746      if isinstance(geom1, Tile):
 747        # an area of intersection based test
 748        match = self.polygon_matches(
 749          affine.affine_transform(geom1.shape, transform), geom2.shape)
 750      elif isinstance(geom1, Vertex):
 751        # distance test
 752        match = affine.affine_transform(geom1.point, transform).distance(
 753          geom2.point) <= 10 * tiling_utils.RESOLUTION
 754      else: # must be an Edge
 755        # since edges _should not_ intersect this test should work in
 756        # lieu of a more complete point by point comparison
 757        c1 = geom1.get_geometry().centroid
 758        c2 = geom2.get_geometry().centroid
 759        match = affine.affine_transform(c1, transform) \
 760          .distance(c2) <= 10 *tiling_utils.RESOLUTION
 761      if match:
 762        return geom2.base_ID
 763    return match_id
 764
 765
 766  def _get_exclusive_supersets(
 767      self, sets:list[Iterable]) -> list[Iterable]:
 768    """Return sets of elements not found in the same set among those supplied.
 769
 770    The supplied sets share elements, i.e., they are non-exclusives sets. The
 771    returned sets are exclusive: each element will only appear in one of the
 772    sets in the returned list. This is accomplished using networkx's
 773    connected components applied to a graph where each intersection between two
 774    sets is an edge.
 775
 776    Args:
 777        sets (list[Iterable]): list of lists of possibly overlapping sets.
 778
 779    Returns:
 780      list[Iterable]: list of lists that include all the original
 781        elements without overlaps.
 782
 783    """
 784    overlaps = []
 785    for i, si in enumerate(sets):
 786      s1 = set(si)
 787      for j, sj in enumerate(sets):
 788        s2 = set(sj)
 789        if len(s1 & s2) > 0:
 790          overlaps.append((i, j))
 791    G = nx.from_edgelist(overlaps)
 792    result = []
 793    for component in nx.connected_components(G):
 794      s = set()
 795      for i in component:
 796        s = s.union(sets[i])
 797      result.append(tuple(s))
 798    return result
 799
 800
 801  def vertices_in_tiles(self, tiles:list[Tile]) -> list[Vertex]:
 802    """Get vertices incident on tiles in supplied list.
 803
 804    Args:
 805      tiles (list[Tile]): tiles whose vertices are required.
 806
 807    Returns:
 808      list[Vertex]: the required vertices.
 809
 810    """
 811    vs = set()
 812    for tile in tiles:
 813      vs = vs.union(tile.get_corner_IDs())
 814    return [self.points[v] for v in vs]
 815
 816
 817  def edges_in_tiles(self, tiles:list[Tile]) -> list[Edge]:
 818    """Get edges that are part of the boundary of tiles in supplied list.
 819
 820    Args:
 821      tiles (list[Tile]): tiles whose edges are required.
 822
 823    Returns:
 824      list[Edge]: the required edges.
 825
 826    """
 827    es = set()
 828    for tile in tiles:
 829      es = es.union(tile.get_edge_IDs())
 830    return [self.edges[e] for e in es]
 831
 832
 833  def generate_dual(self) -> list[geom.Polygon]:
 834    """Create the dual tiiing for the tiling of this Topology.
 835
 836    TODO: make this a viable replacement for the existing dual tiling
 837    generation.
 838
 839    TODO: also need to ensure that this finds a set of dual tiles that exhaust
 840    the plane...
 841
 842    Returns:
 843      list[geom.Polygon]: a list of polygon objects.
 844
 845    """
 846    for v in self.points.values():
 847      v.clockwise_order_incident_tiles()
 848    self.dual_tiles = {}
 849    base_id_sets = defaultdict(list)
 850    for v in self.points.values():
 851      base_id_sets[v.base_ID].append(v.ID)
 852    minimal_set = [self.points[min(s)] for s in base_id_sets.values()]
 853    for v in minimal_set:
 854    # for v in self.points.values():
 855      if v.is_interior() and len(v.tiles) > 2:
 856        self.dual_tiles[v.ID] = \
 857          geom.Polygon([t.centre for t in v.tiles])
 858
 859
 860  def get_dual_tiles(self) -> gpd.GeoDataFrame:
 861    """Return dual tiles as GeoDataFrame."""
 862    n = len(self.dual_tiles)
 863    return gpd.GeoDataFrame(
 864      data = {"tile_id": list(self.tileable.tiles.tile_id)[:n]},
 865      geometry = gpd.GeoSeries(self.dual_tiles.values()),
 866      crs = self.tileable.crs)
 867
 868
 869  def add_vertex(self, pt:geom.Point) -> Vertex:
 870    """Add and return Vertex at the specified point location.
 871
 872    No attempt is made to ensure Vertex IDs are an unbroken sequence: a new ID
 873    is generated one greater than the existing highest ID. IDs will usually be
 874    an unbroken sequence up to removals when geometry transformations are
 875    applied.
 876
 877    Args:
 878      pt (geom.Point): point location of the Vertex.
 879
 880    Returns:
 881      Vertex: the added Vertex object.
 882
 883    """
 884    n = 0 if len(self.points) == 0 else max(self.points.keys()) + 1
 885    v = Vertex(pt, n)
 886    self.points[n] = v
 887    return v
 888
 889
 890  def add_edge(self, vs:list[Vertex]) -> Edge:
 891    """Create an Edge from the suppled vertices and return it.
 892
 893    The new Edge is added to the edges dictionary. Edges are self indexing by
 894    the IDs of their end Vertices.
 895
 896    Args:
 897      vs (list[Vertex]): list of Vertices in the Edge to be created.
 898
 899    Returns:
 900        Edge: the added Edge.
 901
 902    """
 903    e = Edge(vs)
 904    self.edges[e.ID] = e
 905    return e
 906
 907
 908  def polygon_matches(
 909      self,
 910      geom1:geom.Polygon,
 911      geom2:geom.Polygon,
 912    ) -> bool:
 913    """Test if supplied polygons match geometrically.
 914
 915    Tests for equality of area, and equality of their area of overlap to their
 916    shared area, i.e. Area1 == Area2 == (Area 1 intersection 2).
 917
 918    Args:
 919      geom1 (geom.Polygon): first polygon.
 920      geom2 (geom.Polygon): second polygon.
 921
 922    Returns:
 923      bool: True if the polygons are the same, False otherwise.
 924
 925    """
 926    a, b = geom1.area, geom2.area
 927    return bool(
 928      np.isclose(a, b,
 929                 rtol = tiling_utils.RESOLUTION * 100,
 930                 atol = tiling_utils.RESOLUTION * 100) and
 931      np.isclose(a, geom1.intersection(geom2).area,
 932                 rtol = tiling_utils.RESOLUTION * 100,
 933                 atol = tiling_utils.RESOLUTION * 100))
 934
 935
 936  def _get_tile_geoms(self) -> gpd.GeoDataFrame:
 937    """Return GeoDataFrame of the tiles.
 938
 939    Returns:
 940      gpd.GeoDataFrame: the tiles as a GeoDataFrame.
 941
 942    """
 943    return gpd.GeoDataFrame(
 944      data = {"transitivity_class": [t.transitivity_class for t in self.tiles]},
 945      geometry = gpd.GeoSeries([t.shape for t in self.tiles]))
 946
 947
 948  def _get_tile_centre_geoms(self) -> gpd.GeoSeries:
 949    """Return the centre of tiles as a GeoDataFrame.
 950
 951    Returns:
 952      gpd.GeoDataFrame: the tiles centres.
 953
 954    """
 955    return gpd.GeoSeries([t.centre for t in self.tiles])
 956
 957
 958  def _get_point_geoms(self) -> gpd.GeoSeries:
 959    """Return tiling vertices and tile corners as a GeoDataFrame of Points.
 960
 961    Returns:
 962      gpd.GeoDataFrame: corners and vertices of the tiling.
 963
 964    """
 965    return gpd.GeoSeries([v.point for v in self.points.values()])
 966
 967
 968  def _get_vertex_geoms(self) -> gpd.GeoSeries:
 969    """Return tiling vertices (not corners) as a GeoDataFrame of Points.
 970
 971    Returns:
 972      gpd.GeoDataFrame: vertices of the tiling.
 973
 974    """
 975    return gpd.GeoSeries([v.point for v in self.points.values()
 976                          if v.is_tiling_vertex])
 977
 978
 979  def _get_edge_geoms(
 980      self,
 981      offset:float = 0.0,
 982    ) -> gpd.GeoSeries:
 983    """Return tiling edges a GeoSeries of LineStrings optionally offset.
 984
 985    Offsetting the edges allows their direction to be seen when visualised.
 986
 987    Returns:
 988      gpd.GeoDataFrame: edges of the tiling.
 989
 990    """
 991    return gpd.GeoSeries([e.get_topology().parallel_offset(offset)
 992                          for e in self.edges.values()])
 993
 994
 995  def plot(
 996      self,
 997      show_original_tiles: bool = True,
 998      show_tile_centres: bool = False,
 999      show_vertex_labels: bool = True,
1000      show_vertex_ids: bool = False,
1001      show_edges: bool = True,
1002      offset_edges: bool = True,
1003      show_edge_labels:bool = False,
1004      show_dual_tiles: bool = False,
1005    ) -> plt.Axes:
1006    """Delegate plotting of requested elements and return plt.Axes.
1007
1008    Args:
1009      show_original_tiles (bool, optional): if True show the tiles. Defaults to
1010        True.
1011      show_tile_centres (bool, optional): if True show tile centres with upper
1012        case alphabetical labels. Defaults to False.
1013      show_vertex_labels (bool, optional): if True show tiling vertices labelled
1014        to show their equivalence classes. Defaults to True.
1015      show_vertex_ids (bool, optional): if True show vertex IDs (i.e., sequence
1016        numbers) which is useful for debugging. Defaults to False.
1017      show_edges (bool, optional): if True show tiling edges (not tile sides).
1018        Defaults to True.
1019      offset_edges (bool, optional): if True offset edges a little from and
1020        parallel to their geometric position. Defaults to True.
1021      show_edge_labels (bool, optional): if true show lower case alphabetical
1022        labels identifying edge equivalence classes. Defaults to False.
1023      show_dual_tiles (bool, optional): if True show a candidate set of dual
1024        tiles as an overlay. Defaults to False.
1025
1026    Returns:
1027      plt.Axes: a plot of the Topology as requested.
1028
1029    """
1030    fig = plt.figure(figsize = (10, 10))
1031    ax = fig.add_axes(111)
1032    extent = gpd.GeoSeries([t.shape for t in self.tiles]).total_bounds
1033    dist = max([extent[2] - extent[0], extent[3] - extent[1]]) / 100
1034    if show_original_tiles:
1035      self._plot_tiles(ax)
1036    if show_tile_centres:
1037      self._plot_tile_centres(ax)
1038    if show_vertex_labels:
1039      self._plot_vertex_labels(ax, show_vertex_ids)
1040    if show_edge_labels or show_edges:
1041      self._plot_edges(ax, show_edges, show_edge_labels, dist, offset_edges)
1042    if show_dual_tiles:
1043      self._plot_dual_tiles(ax, dist)
1044    plt.axis("off")
1045    return ax
1046
1047
1048  def _plot_tiles(self, ax:plt.Axes) -> plt.Axes:
1049    """Plot Topology's tiles on supplied Axes.
1050
1051    Tiles are coloured by equivalence class.
1052
1053    Args:
1054      ax (plt.Axes): Axes on which to plot.
1055
1056    Returns:
1057      plt.Axes: the Axes.
1058
1059    """
1060    self._get_tile_geoms().plot(column = "transitivity_class",
1061      ax = ax, ec = "#444444", lw = 0.5, alpha = 0.25, cmap = "Greys")
1062    return ax
1063
1064
1065  def _plot_tile_centres(self, ax:plt.Axes) -> plt.Axes:
1066    """Print tile transitivity class at each tile centroid.
1067
1068    Args:
1069      ax (plt.Axes): Axes on which to plot.
1070
1071    Returns:
1072      plt.Axes: the Axes.
1073
1074    """
1075    for _, tile in enumerate(self.tiles):
1076      ax.annotate(tile.transitivity_class, xy = (tile.centre.x, tile.centre.y),
1077                  ha = "center", va = "center")
1078    return ax
1079
1080
1081  def _plot_vertex_labels(
1082      self,
1083      ax:plt.Axes,
1084      show_vertex_ids:bool = False,
1085    ) -> plt.Axes:
1086    """Plot either the Vertex transitivity class label or its sequential ID.
1087
1088    Args:
1089      ax (plt.Axes): Axes on which to plot.
1090      show_vertex_ids (bool, optional): If True plots the ID, else plots the
1091        transitivity class. Defaults to False.
1092
1093    Returns:
1094      plt.Axes: the Axes.
1095
1096    """
1097    for v in self.points.values():
1098      ax.annotate(v.ID if show_vertex_ids else v.label,
1099                  xy = (v.point.x, v.point.y), color = "k",
1100                  ha = "center", va = "center")
1101    return ax
1102
1103
1104  def _plot_edges(
1105      self,
1106      ax:plt.Axes,
1107      show_edges:bool = False,
1108      show_edge_labels:bool = False,
1109      dist:float = 0.0,
1110      offset_edges:bool = True,
1111    ) -> plt.Axes:
1112    """Plot edges, including an offset if specified and labels if requested.
1113
1114    Can also be called to only plot the labels.
1115
1116    Args:
1117      ax (plt.Axes): Axes on which to plot.
1118      show_edges (bool, optional): if True includes the edges as. a dotted blue
1119        line, optionally offset (for clarity) from the tile boundary. Defaults
1120        to False.
1121      show_edge_labels (bool, optional): if True shows an edge label. Defaults
1122        to False.
1123      dist (float, optional): a distance by which to offset the dotted line for
1124        the edge from the tile boundary. Defaults to 0.0.
1125      offset_edges (bool, optional): if True applies the edge drawing offset,
1126        if False the edge is drawn as a single line segment between its end
1127        vertices (and may not align with the sides of tiles. Defaults to True.
1128
1129    Returns:
1130      plt.Axes: the Axes.
1131
1132    """
1133    if show_edges:
1134      edges = self._get_edge_geoms(dist if offset_edges else 0)
1135      edges.plot(ax = ax, color = "dodgerblue", ls = ":")
1136    else:
1137      edges = [e.get_geometry() for e in self.edges.values()]
1138    if show_edge_labels:
1139      for ls, e in zip(edges, self.edges.values(), strict = True):
1140        c = ls.centroid
1141        ax.annotate(e.label, xy = (c.x, c.y), color = "k",
1142                    ha = "center", va = "center")
1143    return ax
1144
1145  def _plot_dual_tiles(self, ax:plt.Axes,
1146                       dist:float = 0.0) -> plt.Axes:
1147    gpd.GeoSeries(self.dual_tiles).buffer(
1148      -dist / 4, join_style = "mitre", cap_style = "square").plot(
1149        ax = ax, fc = "g", alpha = 0.25)
1150    return ax
1151
1152
1153  def plot_tiling_symmetries(
1154      self,
1155      **kwargs:float,
1156    ) -> None:
1157    """Plot the symmetries of Topology's tiling.
1158
1159    Most of the work here is delegated to `_plot_tiling_symmetry` which is run
1160    once per symmetry on a grid of plt.Axes built by this function.
1161
1162    Args:
1163      kwargs: passed through to `_plot_tiling_symmetry`
1164
1165    """
1166    n = len(self.tile_matching_transforms)
1167    nc = int(np.ceil(np.sqrt(n)))
1168    nr = int(np.ceil(n / nc))
1169    fig = plt.figure(figsize = (12, 12 * nr / nc))
1170    for i, tr in enumerate(self.tile_matching_transforms.values()):
1171      ax = fig.add_subplot(nr, nc, i + 1)
1172      self._plot_tiling_symmetry(tr, ax, **kwargs)
1173
1174
1175  def _plot_tiling_symmetry(
1176      self,
1177      transform:Transform,
1178      ax:plt.Axes,
1179      **kwargs:float,
1180    ) -> None:
1181    """Plot the supplied Transform on the supplied plt.Axes.
1182
1183    Args:
1184        transform (Transform): the Transform to plot.
1185        ax (plt.Axes): the Axes on which to plot.
1186        kwargs: passed through to plt.plot functions.
1187
1188    """
1189    tiles = gpd.GeoSeries([t.shape for t in self.tiles])
1190    base_tiles = tiles[:self.n_tiles]
1191    tiles.plot(ax = ax, fc = "k", alpha = .15, ec = "k", lw = .5)
1192    base_tiles.plot(ax = ax, fc = "#00000000", ec = "w", lw = 1, zorder = 2)
1193    transformed_base_tiles = gpd.GeoSeries(
1194      [transform.apply(g) for g in base_tiles])
1195    transformed_base_tiles.plot(ax = ax, fc = "k", alpha = .2, lw = 0, ec = "k")
1196    transform.draw(ax, **kwargs) # delegate to Transform.draw
1197    plt.axis("off")
1198
1199
1200  def transform_geometry(
1201      self,
1202      new_topology:bool,
1203      apply_to_tiles:bool,
1204      selector:str,
1205      transform_type:str,
1206      **kwargs:float) -> Topology:
1207    r"""Get a new Topology by applying specified transformation.
1208
1209    A transformation specified by `transform_type` and keyword arguments is
1210    applied to elements in the Topology whose labels match the selector
1211    parameter. The transform is optionally applied to update tiles and
1212    optionally requests a new Topology object.
1213
1214    Implemented in this way so that transformations can be applied one at a time
1215    without creating an intermediate set of new tiles, which may be invalid and
1216    fail. So, if you wish to apply (say) 3 transforms and generate a new
1217    Topology leaving the existing one intact:
1218
1219        new_topo = old_topo.transform_geometry(True,  False, "a", ...) \
1220                           .transform_geometry(False, False, "B", ...) \
1221                           .transform_geometry(False, True,  "C", ...)
1222
1223    The first transform requests a new Topology, subsequent steps do not, and it
1224    is only the last step which attempts to create the new tile polygons.
1225
1226    **kwargs supply named parameters for the requested transformation.
1227
1228    Args:
1229      new_topology (bool): if True returns a new Topology object, else returns
1230        the current Topology modified.
1231      apply_to_tiles (bool): if True attempts to create new Tiles after the
1232        transformation has been applied. Usually set to False, unless the last
1233        transformation in a pipeline, to avoid problems of topologically invalid
1234        tiles at intermediate steps.
1235      selector (str): label of elements to which to apply the transformation.
1236        Note that all letters in the supplied string are checked, so you can
1237        use e.g. "abc" to apply a transformation to edges labelled "a", "b" or
1238        "c", or "AB" for vertices labelled "A" or "B".
1239      transform_type (str): name of the type of transformation requested.
1240        Currently supported are `zigzag_edge`, `rotate_edge`, `scale_edge`,
1241        `push_vertex`, and `nudge_vertex`. Keyword arguments for each are
1242        documented in the corresponding methods.
1243      kwargs: contains any needed arguments of the requested transform_type.
1244
1245    Returns:
1246      Topology: if new_topology is True a new Topology based on this one with
1247        after transformation, if False this Topology is returned after the
1248        transformation.
1249
1250    """
1251    print("CAUTION: new Topology will probably not be correctly labelled. "
1252          "To build a correct Topology, extract the tileable attribute and "
1253          "rebuild Topology from that.")
1254    # pickle seems to be more reliable than copy.deepcopy here
1255    topo = (pickle.loads(pickle.dumps(self))   # noqa: S301
1256            if new_topology
1257            else self)
1258    transform_args = topo.get_kwargs(getattr(topo, transform_type), **kwargs)
1259    match transform_type:
1260      case "zigzag_edge":
1261        for e in topo.edges.values():
1262          if e.label in selector:
1263            topo.zigzag_edge(e, **transform_args)
1264      case "rotate_edge":
1265        for e in topo.edges.values():
1266          if e.label in selector:
1267            topo.rotate_edge(e, **transform_args)
1268      case "scale_edge":
1269        for e in topo.edges.values():
1270          if e.label in selector:
1271            topo.scale_edge(e, **transform_args)
1272      case "push_vertex":
1273        pushes = {}
1274        for v in topo.vertices_in_tiles(topo.tiles[:topo.n_tiles]):
1275          if v.label in selector:
1276            pushes[v.base_ID] = topo.push_vertex(v, **transform_args)
1277        for base_ID, (dx, dy) in pushes.items():
1278          for v in [v for v in topo.points.values() if v.base_ID == base_ID]:
1279            v.point = affine.translate(v.point, dx, dy)
1280      case "nudge_vertex":
1281         for v in topo.points.values():
1282          if v.label in selector:
1283            topo.nudge_vertex(v, **transform_args)
1284    if apply_to_tiles:
1285      for t in topo.tiles:
1286        t.set_corners_from_edges()
1287    topo.tileable.tiles.geometry = gpd.GeoSeries(
1288      [topo.tiles[i].shape for i in range(topo.n_tiles)])
1289    topo.tileable._setup_regularised_prototile()
1290    return topo
1291
1292
1293  def get_kwargs(
1294      self,
1295      fn:Callable,
1296      **kwargs:str|float,
1297    ) -> dict:
1298    """Filter the supplied kwargs to only contain arguments required by fn.
1299
1300    Args:
1301      fn (Callable): the function that is to be inspected.
1302      **kwargs (str|float): kwargs to be filtered.
1303
1304    Returns:
1305      str|float: filtered dictionary of kwargs.
1306
1307    """
1308    args = inspect.signature(fn).parameters
1309    return {k: kwargs.pop(k) for k in dict(kwargs) if k in args}
1310
1311
1312  def zigzag_edge(
1313      self,
1314      edge:Edge,
1315      start:str = "A",
1316      n:int = 2,
1317      h:float = 0.5,
1318      smoothness:int = 0,
1319    ) -> None:
1320    """Apply zigzag transformation to supplied Edge.
1321
1322    Currently this will only work correctly if h is even.
1323
1324    TODO: make it possible for odd numbers of 'peaks' to work (this may require
1325    allowing bidirectional Edges, i.e. storing Edges in both directions so that
1326    all Tile edges are drawn CW). The `start` parameter is a temporary hack for
1327    this.
1328
1329    Args:
1330      edge (Edge): Edge to transform
1331      start (str, optional): label at one end of edge which is used to determine
1332        the sense of h, enabling C-curves with an odd number n of zigs and zags
1333        to be applied. Defaults to 'A'.
1334      n (int, optional): number of zigs and zags in the edge. Defaults to 2.
1335      h (float, optional): width of the zig zags relative to edge length.
1336        Defaults to 0.5.
1337      smoothness (int, optional): spline smoothness. 0 gives a zig zag proper,
1338        higher values will produce a sinusoid. Defaults to 0.
1339
1340    """
1341    v0, v1 = edge.vertices[0], edge.vertices[1]
1342    if n % 2 == 1 and v0.label != start:
1343      h = -h
1344    ls = self.zigzag_between_points(v0.point, v1.point, n, h, smoothness)
1345    # remove current corners
1346    self.points = {k: v for k, v in self.points.items()
1347                   if k not in edge.get_corner_IDs()[1:-1]}
1348    # add the new ones
1349    new_corners = [self.add_vertex(geom.Point(xy)) for xy in ls.coords]
1350    edge.corners = edge.vertices[:1] + new_corners + edge.vertices[-1:]
1351    if edge.right_tile is not None:
1352      edge.right_tile.set_corners_from_edges(False)
1353    if edge.left_tile is not None:
1354      edge.left_tile.set_corners_from_edges(False)
1355
1356
1357  def zigzag_between_points(
1358      self,
1359      p0:geom.Point,
1360      p1:geom.Point,
1361      n:int,
1362      h:float = 1.0,
1363      smoothness:int = 0,
1364    ) -> geom.LineString:
1365    """Return a zig zag line optionally smoothed as a spline between two points.
1366
1367    Args:
1368      p0 (geom.Point): start point.
1369      p1 (geom.Point): end point.
1370      n (int): number of zig zags.
1371      h (float, optional): amplitude of zig zags relative to distance between
1372        points. Defaults to 1.0.
1373      smoothness (int, optional): number of spline smoothed points to add. If
1374        set to 0 a straight line zig zag is produced. Defaults to 0.
1375
1376    Returns:
1377      geom.LineString: the resulting zig zag line.
1378
1379    """
1380    r = p0.distance(p1)
1381    # make a sinusoidal template
1382    x = np.linspace(0, n * np.pi, n * 2 + 1, endpoint = True)
1383    y = [np.sin(x) for x in x]
1384    spline = interpolate.InterpolatedUnivariateSpline(x, y, k = 2)
1385
1386    spline_steps = (n + smoothness) * 2 + 1
1387    xs = np.linspace(0, n * np.pi, spline_steps, endpoint = True)
1388    ys = spline(xs)
1389
1390    sfx = 1 / max(x) * r
1391    sfy = h * r / 2
1392    theta = np.arctan2(p1.y - p0.y, p1.x - p0.x)
1393
1394    ls = geom.LineString(
1395      [geom.Point(x, y) for x, y in zip(xs, ys, strict = True)])
1396    ls = affine.translate(ls, 0, -(ls.bounds[1] + ls.bounds[3]) / 2)
1397    ls = affine.scale(ls, xfact = sfx, yfact = sfy, origin = (0, 0))
1398    ls = affine.rotate(ls, theta, (0, 0), use_radians = True)
1399    x0, y0 = next(iter(ls.coords))
1400    return affine.translate(ls, p0.x - x0, p0.y - y0)
1401
1402
1403  def rotate_edge(
1404      self,
1405      edge:Edge,
1406      centre:str = "",
1407      angle:float = 0,
1408    ) -> None:
1409    """Rotate edge.
1410
1411    Args:
1412      edge (Edge): the edge to rotate.
1413      centre (str): centre of rotation which should be label of one of its
1414        end point vertices or "" when rotation will be about the Edge centroid.
1415        Defaults to "".
1416      angle (float): angle of rotation.
1417
1418    """
1419    v0, v1 = edge.vertices
1420    ls = geom.LineString([v0.point, v1.point])
1421    if v0.label == centre:
1422      c = v0.point
1423    elif v1.label == centre:
1424      c = v1.point
1425    else:
1426      c = ls.centroid
1427    ls = affine.rotate(ls, angle, origin = c)
1428    v0.point, v1.point = [geom.Point(c) for c in ls.coords]
1429
1430
1431  def scale_edge(
1432      self,
1433      edge:Edge,
1434      sf:float = 1.0,
1435    ) -> None:
1436    """Scale edge.
1437
1438    Args:
1439      edge (Edge): the edge to scale.
1440      sf (float): amount by which to scale the edge.
1441
1442    """
1443    v0, v1 = edge.vertices
1444    ls = geom.LineString([v0.point, v1.point])
1445    ls = affine.scale(ls, xfact = sf, yfact = sf, origin = ls.centroid)
1446    v0.point, v1.point = [geom.Point(c) for c in ls.coords]
1447
1448
1449  def push_vertex(
1450      self,
1451      vertex:Vertex,
1452      push_d:float,
1453    ) -> tuple[float]:
1454    """Return displacement vector to push a vertex based on incident edges.
1455
1456    Args:
1457      vertex (Vertex): the vertex to push.
1458      push_d (float): the distance to push it.
1459
1460    """
1461    neighbours = [self.points[v] for v in vertex.neighbours]
1462    dists = [vertex.point.distance(v.point) for v in neighbours]
1463    x, y = vertex.point.x, vertex.point.y
1464    unit_vectors = [((x - v.point.x) / d, (y - v.point.y) / d)
1465                    for v, d in zip(neighbours, dists, strict = True)]
1466    return  (push_d * sum([xy[0] for xy in unit_vectors]),
1467             push_d * sum([xy[1] for xy in unit_vectors]))
1468
1469
1470  def nudge_vertex(
1471      self,
1472      vertex:Vertex,
1473      dx:float,
1474      dy:float,
1475    ) -> None:
1476    """Nudge vertex by specified displacement.
1477
1478    Args:
1479      vertex (Vertex): the vertext to nudge.
1480      dx (float): x displacement.
1481      dy (float): y displacement.
1482
1483    """
1484    vertex.point = affine.translate(vertex.point, dx, dy)
1485
1486
1487class Tile:
1488  """Class to represent essential features of polygons in a tiling."""
1489
1490  ID: int
1491  """integer ID number which indexes the Tile in the containing Topology tiles
1492  list."""
1493  base_ID: int
1494  """ID of corresponding Tile in the base tileable unit"""
1495  corners: list[Vertex]
1496  """list of Vertex objects. This includes all corners of the original polygon
1497  and any tiling vertices induced by (for example) a the corner of an adjacent
1498  tile lying halfway along an edge of the original polygon on which this tile
1499  is based. Vertex objects are stored in strictly clockwise sequence."""
1500  edges: list[Edge]
1501  """list of Edge objects that together compose the tile boundary."""
1502  edges_CW: list[bool]
1503  """list of Edge direction. Edges are stored only once in a Topology so some
1504  edges are in clockwise order and others  are in counter-clockwise order.
1505  These boolean flags are True if the corresponding Edge is clockwise, False if
1506  counter-clockwise."""
1507  label: str
1508  """tile_id label from the tileable source"""
1509  vertex_labels: list[str]
1510  """list of (upper case) letter labels of the tile corners (i.e. all corners,
1511  not only tiling vertices)."""
1512  edge_labels: list[str]
1513  """list of (lower case) letter labels of the tile edges (tiling edges, not
1514  tile sides)."""
1515  shape: geom.Polygon = None
1516  """the tile geometry (which may include some redundant points along sides
1517  where neighbouring tiles induce a tiling vertex). So for example a rectangle
1518  might have additional points along its sides:
1519
1520        +---+-------+
1521        |   |   2   |
1522        | 1 A---B---E---+
1523        |   |   |   4   |
1524        +---C 3 D-------+
1525            |   |
1526            +---+
1527
1528  In the above Tile 1 has additional point A, 2 has B and 3 has C and D induced
1529  by the corners of neighbouring tiles."""
1530  centre: geom.Point = None
1531  """a point centre for the Tile (determined by weavingspace.tiling_utils.
1532  incentre)."""
1533  shape_group: int = None
1534  """the tile shape group of this tile in its containing Topology."""
1535  transitivity_class: int = None
1536  """the tile transitivity class of this tile its containing Topology"""
1537
1538  def __init__(
1539      self,
1540      ID:int,
1541    ) -> None:
1542    """Class constructor.
1543
1544    Args:
1545      ID (int): sequence number ID of this vertex in the Topology.
1546
1547    """
1548    self.ID = ID
1549    self.corners = []
1550    self.edges = []
1551    self.edges_CW = []
1552    self.vertex_labels = []
1553    self.edge_labels = []
1554
1555
1556  def __str__(self) -> str:
1557    """Return string representation of the Tile.
1558
1559    Returns:
1560      str: string including Tile ID, list of corner vertex IDs and list of
1561        edge IDs.
1562
1563    """
1564    return (f"Tile {self.ID} Corners: {self.get_corner_IDs()} "
1565            f"Edges: {self.get_edge_IDs()}")
1566
1567
1568  def __repr__(self) -> str:
1569    return str(self)
1570
1571
1572  def get_corner_IDs(self) -> list[int]:
1573    """Return list of corner IDs (not Vertex objects).
1574
1575    Returns:
1576      list[int]: list of integer IDs of tile corners.
1577
1578    """
1579    return [c.ID for c in self.corners]
1580
1581
1582  def get_edge_IDs(self) -> list[tuple[int,int]]:
1583    """Return list of edge IDs (not Edge objects).
1584
1585    Returns:
1586      list[tuple[int]]: list of 2-element tuples of the start and end Vertex IDs
1587        of each edge.
1588
1589    """
1590    return [e.ID for e in self.edges]
1591
1592
1593  def set_shape_from_corners(self) -> None:
1594    """Set the shape attribute based on corners, and associated tile centre."""
1595    self.shape = geom.Polygon([c.point for c in self.corners])
1596    self.centre = tiling_utils.get_incentre(
1597      tiling_utils.get_clean_polygon(self.shape))
1598
1599
1600  def set_corners_from_edges(
1601      self,
1602      update_shape:bool = True,
1603    ) -> None:
1604    """Set corners attribute from the edges attribute.
1605
1606    Typically called after modification of topology edges. Optionally the shape
1607    attribute is NOT updated, which may save time when multiple changes to the
1608    edges of a tile are in process (i.e., only update the shape after all
1609    changes are complete).
1610
1611    Args:
1612      update_shape (bool, optional): if True the shape attribute will be
1613        updated, otherwise not. Defaults to True.
1614
1615    """
1616    self.corners = []
1617    for e, cw in zip(self.edges, self.edges_CW, strict = True):
1618      if cw: # clockwise to extend by all but the first corner
1619        self.corners.extend(e.corners[:-1])
1620      else: # counter-clockwise so extend in reverse
1621        self.corners.extend(e.corners[1:][::-1])
1622    if update_shape:
1623      self.set_shape_from_corners()
1624
1625
1626  def set_edge_directions(self) -> None:
1627    """Set up edges_CW attribute by inspection of the edges list.
1628
1629    It is (frankly!) hard to keep track of the correct sequence of CW/CCW order
1630    of edges as new ones are created or old ones merged. This method inspects
1631    the 'tail-head' relations between consecutive edges to set these flags
1632    correctly.
1633
1634    The test is simply to check if the 'tail' Vertex ID in each edge appears
1635    in the ID tuple of the following edge, i.e. if successive edge
1636    IDs are (0, 1) (2, 1) or (0, 1) (1, 2), then edge (0, 1) is in clockwise
1637    direction, but if we have (0, 1) (2, 3) then it is not.
1638    """
1639    edge_IDs = self.get_edge_IDs()
1640    self.edges_CW = [e1[-1] in e2 for e1, e2 in
1641                     zip(edge_IDs, edge_IDs[1:] + edge_IDs[:1], strict = True)]
1642
1643
1644  def insert_vertex_at(
1645      self,
1646      v:Vertex,
1647      i:int,
1648      update_shape:bool = False,
1649    ) -> tuple[tuple[int,int],tuple[tuple[int,int],...]]:
1650    """Insert the Vertex into tile at index position i.
1651
1652    Both corners and edges attributes are updated, and the old edge IDs for
1653    removal and the new edge itself are returned to the calling context (the
1654    containing Topology) for update of its edges collection. Optionally update
1655    the shape attribute.
1656
1657    This is NOT a generic vertex insertion method: it is only for use during
1658    Topology initialisation, and does not guarantee correct maintenance of
1659    all tile, edge and vertex relations in the general case---at any rate it
1660    has not been tested for this!
1661
1662    Args:
1663      v (Vertex): the Vertex to insert.
1664      i (int): index position in current corners after which to insert
1665        supplied Vertex.
1666      update_shape (bool, optional): if True shape attribute is updated.
1667        Defaults to False.
1668
1669    Returns:
1670      tuple: the (tuple) ID of the old edge which should be deleted, and
1671        the new Edges arising from insertion of this Vertex.
1672
1673    """
1674    self.corners = self.corners[:i] + [v] + self.corners[i:]
1675    old_edge = self.edges[i - 1]
1676    # store current ID of the affected edge for return to calling context
1677    old_edge_ID = old_edge.ID
1678    new_edges = old_edge.insert_vertex(v, self.corners[i - 1])
1679    self.edges = self.edges[:(i-1)] + new_edges + self.edges[i:]
1680    self.set_edge_directions()
1681    if update_shape:
1682      self.set_shape_from_corners()
1683    return old_edge_ID, new_edges
1684
1685
1686  def merge_edges_at_vertex(
1687      self,
1688      v:Vertex,
1689    ) -> tuple:
1690    """Merge edges that meet at the supplied Vertex.
1691
1692    It is assumed that only two tiles are impacted this one, and its neighbour
1693    across the Edge on which v lies. Both are updated. For this reason the work
1694    is delegated to `get_updated_edges_from_merge` which is run on both affected
1695    tiles, but only determines the edges to remove and the new edge to be added
1696    once. See that method for details.
1697
1698    Args:
1699      v (Vertex): Vertex at which to merge Edges. This should currently be an
1700        end
1701
1702    Returns:
1703      tuple: 2 item list of the edge IDs to be removed and a new Edge object to
1704        be added by the calling context (i.e. the containing Topology).
1705
1706    """
1707    to_remove, new_edge = self.get_updated_edges_from_merge(v)
1708    if len(v.tiles) > 1:
1709      v.tiles[1].get_updated_edges_from_merge(v, new_edge)
1710    return to_remove, new_edge
1711
1712
1713  def get_updated_edges_from_merge(
1714      self,
1715      v:Vertex,
1716      new_edge:Edge = None,
1717    ) -> None|tuple[tuple[tuple[int,int], tuple[int,int]], Edge]:
1718    """Update edges and edges_CW attributes based on insertion of Vertex.
1719
1720    If new_edge is supplied then the neighbour tile at v has already created
1721    the needed new Edge and this Edge is the one that will be 'slotted in' at
1722    the appropriate spot in the edges list.
1723
1724    The edges_CW is also updated to maintain correct directions of the edges.
1725    The corners attribute is unaffected by these changes.
1726
1727    Args:
1728      v (Vertex): existing Vertex at which to carry out the merge.
1729      new_edge (Edge, optional): if another Tile has already carried out this
1730        merge this should be the resulting new Edge for insertion into this
1731        Tile. Defaults to None (when the new Edge will be constructed).
1732
1733    Returns:
1734      None|tuple: either None (if a new edge was supplied) or a tuple
1735        of the two edge IDs to be removed and the new edge added for return to
1736        the calling context (i.e. the containing Topology).
1737
1738    """
1739    # get the two edge list index positions in which Vertex v is found
1740    i, j = self.get_edge_IDs_including_vertex(v)
1741    if new_edge is None: # then we must make a new one
1742      # also record existing edge IDs to be removed
1743      to_remove = [self.edges[i].ID, self.edges[j].ID]
1744      new_edge = self.get_merged_edge(i, j)
1745      return_edge_updates = True
1746    else:
1747      return_edge_updates = False
1748    if abs(i - j) != 1:
1749      # edge indices 'wrap' around from end of edge list to start so drop
1750      # first and last current edges and stick new one on at the end
1751      self.edges = self.edges[1:-1] + [new_edge]
1752    else:
1753      # insert new edge into list in place of the two old ones
1754      self.edges = self.edges[:i] + [new_edge] + self.edges[j+1:]
1755    # update the edge directions
1756    self.set_edge_directions()
1757    if return_edge_updates:
1758      return to_remove, new_edge
1759    return None
1760
1761
1762  def get_edge_IDs_including_vertex(
1763      self,
1764      v:Vertex,
1765    ) -> tuple[int]:
1766    """Get the (two) index positions of the edges that include supplied Vertex.
1767
1768    Args:
1769        v (Vertex): Vertex of interest.
1770
1771    Returns:
1772      tuple[int]: two index positions of Edges in edges list that contain v.
1773
1774    """
1775    return (i for i, e in enumerate(self.edges) if v.ID in e.ID)
1776
1777
1778  def get_merged_edge(self, i:int, j:int) -> Edge:
1779    """Return edge made by merging existing edges at i and j in the edges list.
1780
1781    For example, if the current list of edge IDs was
1782
1783        (0 1 2) (4 2) (4 5) (5 0)
1784
1785    and the merge requested is 0 and 1, the resulting new edge is constructed
1786    from vertices (0 1 2 4).
1787
1788    Returns:
1789      Edge: the requested new Edge.
1790
1791    """
1792    # if i and j are not consecutive, then j is predecessor edge
1793    if abs(i - j) != 1:
1794      i, j = j, i
1795    # get edges and their directions
1796    ei, ej = self.edges[i], self.edges[j]
1797    CWi, CWj = self.edges_CW[i], self.edges_CW[j]
1798    # DON'T MESS WITH THIS!!!
1799    # for predecessors (the head) we want everything including the Vertex
1800    # where the merge is occurring; for successors (the tail) we want all but
1801    # the first Vertex (which is the one where the merge is occurring). In both
1802    # cases contingent on whether existing Edges are CW or CCW we may need to
1803    # flip the Vertex sequence to ensure that the merge Vertex is in the middle
1804    # of the new edge that will be created
1805    head = ei.corners if CWi else ei.corners[::-1]
1806    tail = ej.corners[1:] if CWj else ej.corners[::-1][1:]
1807    v_sequence = (head if CWi else head[::-1]) + (tail if CWj else tail[::-1])
1808    return Edge(v_sequence)
1809
1810
1811  def offset_corners(
1812      self,
1813      offset:int,
1814    ) -> None:
1815    """Shift shape, corners, edges, and edges_CW by an offset amount.
1816
1817    This is used to align tiles that are similar, which is required for correct
1818    transfer of 'base' tile labelling on to 'radius 1' tiles during Topology
1819    construction.
1820
1821    Args:
1822      offset (int): the number of positions to shift the lists.
1823
1824    """
1825    if offset is not None or offset != 0:
1826      self.corners = self.corners[offset:] + self.corners[:offset]
1827      self.shape = geom.Polygon([c.point for c in self.corners])
1828      self.edges = self.edges[offset:] + self.edges[:offset]
1829      self.edges_CW = self.edges_CW[offset:] + self.edges_CW[:offset]
1830
1831
1832  def get_edge_label(self, e:Edge) -> str:
1833    """Return edge label of specified edge.
1834
1835    Args:
1836        e (Edge): Edge whose label is required.
1837
1838    Returns:
1839        str: requested edge label.
1840
1841    """
1842    return self.edge_labels[self.get_edge_IDs().index(e.ID)]
1843
1844
1845  def get_corner_label(self, v:Vertex) -> str:
1846    """Return corner label of specified corner.
1847
1848    Args:
1849        v (Vertex): corner whose label is required.
1850
1851    Returns:
1852        str: requested corner label.
1853
1854    """
1855    return self.edge_labels[self.get_corner_IDs().index(v.ID)]
1856
1857
1858  def angle_at(self, v:Vertex) -> float:
1859    """Return interior angle at the specified corner (in degrees).
1860
1861    Args:
1862        v (Vertex): corner where angle is requested.
1863
1864    Returns:
1865        float: angle at corner v in degrees.
1866
1867    """
1868    i = self.corners.index(v)
1869    n = len(self.corners)
1870    return tiling_utils.get_inner_angle(self.corners[i-1].point,
1871                                        self.corners[i].point,
1872                                        self.corners[(i + 1) % n].point)
1873
1874
1875class Vertex:
1876  """Class to store attributes of a vertex in a tiling."""
1877
1878  point: geom.Point
1879  """point (geom.Point): point location of the vertex."""
1880  ID: int
1881  """integer (mostly but not necessarily in sequence) of vertex keyed into the
1882  points dictionary of the containing Topology."""
1883  tiles: list[Tile]
1884  """list of Tiles incident on this vertex."""
1885  neighbours: list[int]
1886  """list of the immediately adjacent other corner IDs. Only required to
1887  determine if a point is a tiling vertex (when it will have) three or more
1888  neighbours, so only IDs are stored."""
1889  base_ID: int = 1_000_000
1890  """ID of corresponding Vertex in the tileable base_unit"""
1891  transitivity_class: int = None
1892  """transitivity class of the vertex under symmetries of the tiling"""
1893  label: str = ""
1894  """the (upper case letter) label of the vertex under the symmetries of the
1895  tiling."""
1896  is_tiling_vertex: bool = True
1897  """is_tiling_vertex (bool): True if this is a tiling vertex, rather than a
1898  tile corner. E.g., A below is a corner, not a tiling vertex. B is a tiling
1899  vertex:
1900
1901      +-------+
1902      | 1     |
1903      |   A---B---+
1904      |   | 2     |
1905      +---C   +---+
1906          |   |
1907          +---+"""
1908
1909
1910  def __init__(
1911      self,
1912      point:geom.Point,
1913      ID:int,
1914    ) -> None:
1915    """Class constructor.
1916
1917    Args:
1918      point (geom.Point): point location of the vertex.
1919      ID (int): a unique integer ID (which will be its key in the containing
1920        Topology points dictionary).
1921
1922    """
1923    self.point = point
1924    self.ID = ID
1925    self.base_ID = self.ID
1926    self.tiles = []
1927    self.neighbours = []
1928
1929
1930  def __str__(self) -> str:
1931    """Return string representation of Vertex.
1932
1933    Returns:
1934        str: string including ID, point and list of incident Tiles.
1935
1936    """
1937    return f"Vertex {self.ID} at {self.point} Tiles: {self.get_tile_IDs()}"
1938
1939
1940  def __repr__(self) -> str:
1941    return str(self)
1942
1943
1944  def get_tile_IDs(self) -> list[int]:
1945    """Return list of Tile IDs (not the Tiles themselves).
1946
1947    Returns:
1948        list[int]: list of Tile IDs
1949
1950    """
1951    return [t.ID for t in self.tiles]
1952
1953
1954  def add_tile(self, tile:Tile) -> None:
1955    """Add supplied Tile to the tiles list if it is not already present.
1956
1957    Args:
1958        tile (Tile): Tile to add.
1959
1960    """
1961    if tile not in self.tiles:
1962      self.tiles.append(tile)
1963
1964
1965  def add_neighbour(
1966      self,
1967      vertex_id:int,
1968    ) -> None:
1969    """Add supplied ID to the neighbours list if it is not already present.
1970
1971    Args:
1972      vertex_id (int): ID to add to the neighbours list.
1973
1974    """
1975    if vertex_id not in self.neighbours:
1976      self.neighbours.append(vertex_id)
1977
1978
1979  def clockwise_order_incident_tiles(self) -> None:
1980    """Reorder tiles list clockwise (this is for dual tiling construction)."""
1981    cw_order = self._order_of_pts_cw_around_centre(
1982      [t.centre for t in self.tiles], self.point)
1983    self.tiles = [self.tiles[i] for i in cw_order]
1984
1985
1986  def is_interior(self) -> bool:
1987    """Test if vertex is completely enclosed by its incident Tiles.
1988
1989    Based on summing the interior angles of the incident tiles at this vertex.
1990
1991    Returns:
1992        bool: True if vertex is completely enclosed by incident Tiles.
1993
1994    """
1995    return abs(360 - sum([t.angle_at(self) for t in self.tiles])) \
1996                     < tiling_utils.RESOLUTION
1997
1998
1999  def _order_of_pts_cw_around_centre(
2000      self,
2001      pts:list[geom.Point],
2002      centre:geom.Point,
2003    ) -> list[int]:
2004    """Return order of points clockwise relative to centre point.
2005
2006    Args:
2007      pts (list[geom.Point]): list of points to order.
2008      centre (geom.Point): centre relative to which CW order is determined.
2009
2010    Returns:
2011      list[int]: list of indices of reordered points.
2012
2013    """
2014    dx = [p.x - centre.x for p in pts]
2015    dy = [p.y - centre.y for p in pts]
2016    angles = [np.arctan2(dy, dx) for dx, dy in zip(dx, dy, strict = True)]
2017    d = dict(zip(angles, range(len(pts)), strict = True))
2018    return [i for angle, i in sorted(d.items(), reverse=True)]
2019
2020
2021class Edge:
2022  """Class to represent edges in a tiling (not tile sides)."""
2023
2024  ID: tuple[int]
2025  """IDs of the vertices at ends of the edge. Used as key in the containing
2026  Topology's edges dictionary."""
2027  vertices: list[Vertex]
2028  """two item list of the end vertices."""
2029  corners: list[Vertex]
2030  """list of all the vertices in the edge (including its end vertices). In a
2031  'normal' edge to edge tiling corners and vertices will be identical."""
2032  right_tile: Tile = None
2033  """the tile to the right of the edge traversed from its first to its last
2034  vertex. Given clockwise winding default, all edges will have a right_tile."""
2035  left_tile: Tile = None
2036  """the tile to the left of the edge traversed from its first to its last
2037  vertex. Exterior edges of the tiles in a Topology will not have a left_tile.
2038  """
2039  base_ID: tuple[int] = (1_000_000, 1_000_000)
2040  """ID of corresponding edge in the base tileable"""
2041  transitivity_class: int = None
2042  """transitivity class of the edge under symmetries of the tiling"""
2043  label: str = ""
2044  """the (lower case letter) label of the edge under the symmetries of the
2045  tiling."""
2046
2047  def __init__(
2048      self,
2049      corners:list[Vertex],
2050    ) -> None:
2051    """Class constructor.
2052
2053    Initialises the corners and vertices lists and sets ID to (vertices[0].ID,
2054    vertices[1].ID). The vertices list is all the corners with is_tiling_vertex
2055    property True -- Note that during initialisation the default of this
2056    property is True until after the relations between tiles and vertices have
2057    been determined.
2058
2059    Args:
2060      corners (list[Vertex]): list of all corners along the edge.
2061
2062    """
2063    self.corners = corners
2064    self.vertices = [v for v in self.corners if v.is_tiling_vertex]
2065    self.ID = tuple(v.ID for v in self.vertices)
2066
2067
2068  def __str__(self) -> str:
2069    """Return a string representation of the Edge.
2070
2071    Returns:
2072      str: include ID and a list of corner vertex IDs.
2073
2074    """
2075    return f"Edge {self.ID} Corners: {[c.ID for c in self.corners]}"
2076
2077
2078  def __repr__(self) -> str:
2079    return str(self)
2080
2081
2082  def get_corner_IDs(self) -> list[int]:
2083    """Return the IDs of edge corners.
2084
2085    Returns:
2086      list[int]: IDs of all corners.
2087
2088    """
2089    return [c.ID for c in self.corners]
2090
2091
2092  def get_vertex_IDs(self) -> list[int]:
2093    """Return IDs of edge vertices.
2094
2095    Returns:
2096      list[int]: list of IDs of the vertices.
2097
2098    """
2099    return [v.ID for v in self.vertices]
2100
2101
2102  def insert_vertex(self, v:Vertex, predecessor:Vertex) -> list[Edge]:
2103    """Insert vertex after predecessor and return modified edge and new edge.
2104
2105    If the initial edge was (say) (0 1 2 5) and the predecessor was set to 1
2106    the returned edges would be (0 1 v) and (v 2 5).
2107    """
2108    i = self.corners.index(predecessor)
2109    new_edge = Edge([v] + self.corners[(i+1):])
2110    if self.right_tile is not None:
2111      new_edge.right_tile = self.right_tile
2112    if self.left_tile is not None:
2113      new_edge.left_tile = self.left_tile
2114    self.corners = self.corners[:(i+1)] + [v]
2115    self.vertices = [self.vertices[0], v]
2116    self.ID = tuple(v.ID for v in self.vertices)
2117    return [self, new_edge]
2118
2119
2120  def get_geometry(
2121      self,
2122      forward:bool = True,
2123    ) -> geom.LineString:
2124    """Return geom.LineString representing geometry (including corners).
2125
2126    Args:
2127      forward (bool, optional): if True the returned LineString starts at
2128        corners[0], else at corners[-1]. Defaults to True.
2129
2130    Returns:
2131      geom.LineString: the required LineString.
2132
2133    """
2134    if forward:
2135      return geom.LineString([v.point for v in self.corners])
2136    return geom.LineString([v.point for v in self.corners[::-1]])
2137
2138
2139  def get_topology(
2140      self,
2141      forward:bool = True,
2142    ) -> geom.LineString:
2143    """Return LineString connecting vertices of this edge.
2144
2145    Args:
2146      forward (bool, optional): if True LineString starts at vertices[0],
2147        else at vertices[1]. Defaults to True.
2148
2149    Returns:
2150        geom.LineString: the required LineString.
2151
2152    """
2153    if forward:
2154      return geom.LineString([v.point for v in self.vertices])
2155    return geom.LineString([v.point for v in self.vertices[::-1]])